Commit 57ac07ef authored by Martin Uecker's avatar Martin Uecker

update iterative algorithms and nonlinear inversion

parent de4c1bfb
......@@ -106,47 +106,10 @@ void landweber_sym(unsigned int maxiter, float epsilon, float alpha, long N,
/**
* Store information about iterative algorithm.
* Used to flexibly modify behavior, e.g. continuation
*
* @param rsnew current residual
* @param rsnot initial residual
* @param iter current iteration
* @param maxiter maximum iteration
*/
struct iter_data {
double rsnew;
double rsnot;
unsigned int iter;
const unsigned int maxiter;
};
/**
* Continuation for regularization. Returns fraction to scale regularization parameter
*
* @param itrdata state of iterative algorithm
* @param delta scaling of regularization in the final iteration (1. means don't scale, 0. means scale to zero)
*
*/
static float ist_continuation(struct iter_data* itrdata, const float delta)
{
/*
// for now, just divide into evenly spaced bins
const float num_steps = itrdata->maxiter - 1;
int step = (int)(itrdata->iter * num_steps / (itrdata->maxiter - 1));
float scale = 1. - (1. - delta) * step / num_steps;
return scale;
*/
float a = logf( delta ) / (float) itrdata->maxiter;
return expf( a * itrdata->iter );
}
......@@ -156,8 +119,6 @@ static float ist_continuation(struct iter_data* itrdata, const float delta)
* @param maxiter maximum number of iterations
* @param epsilon stop criterion
* @param tau (step size) weighting on the residual term, A^H (b - Ax)
* @param lambda_start initial regularization weighting
* @param lambda_end final regularization weighting (for continuation)
* @param N size of input, x
* @param vops vector ops definition
* @param op linear operator, e.g. A
......@@ -166,45 +127,37 @@ static float ist_continuation(struct iter_data* itrdata, const float delta)
* @param b observations
* @param monitor compute objective value, errors, etc.
*/
void ist(unsigned int maxiter, float epsilon, float tau,
float continuation, bool hogwild, long N,
void ist(unsigned int maxiter, float epsilon, float tau, long N,
const struct vec_iter_s* vops,
ist_continuation_t ist_continuation,
struct iter_op_s op,
struct iter_op_p_s thresh,
float* x, const float* b,
struct iter_monitor_s* monitor)
{
struct iter_data itrdata = {
struct ist_data itrdata = {
.rsnew = 1.,
.rsnot = 1.,
.iter = 0,
.maxiter = maxiter,
.tau = tau,
.scale = 1.,
};
float* r = vops->allocate(N);
itrdata.rsnot = vops->norm(N, b);
float ls_old = 1.;
float lambda_scale = 1.;
int hogwild_k = 0;
int hogwild_K = 10;
for (itrdata.iter = 0; itrdata.iter < maxiter; itrdata.iter++) {
iter_monitor(monitor, vops, x);
ls_old = lambda_scale;
lambda_scale = ist_continuation(&itrdata, continuation);
if (NULL != ist_continuation)
ist_continuation(&itrdata);
if (lambda_scale != ls_old)
debug_printf(DP_DEBUG3, "##lambda_scale = %f\n", lambda_scale);
iter_op_p_call(thresh, tau, x, x);
iter_op_p_call(thresh, itrdata.scale * itrdata.tau, x, x);
iter_op_call(op, r, x); // r = A x
......@@ -217,19 +170,7 @@ void ist(unsigned int maxiter, float epsilon, float tau,
if (itrdata.rsnew < epsilon)
break;
vops->axpy(N, x, tau * lambda_scale, r);
if (hogwild)
hogwild_k++;
if (hogwild_k == hogwild_K) {
hogwild_K *= 2;
hogwild_k = 0;
tau /= 2;
}
vops->axpy(N, x, itrdata.tau, r);
}
debug_printf(DP_DEBUG3, "\n");
......@@ -245,8 +186,6 @@ void ist(unsigned int maxiter, float epsilon, float tau,
* @param maxiter maximum number of iterations
* @param epsilon stop criterion
* @param tau (step size) weighting on the residual term, A^H (b - Ax)
* @param lambda_start initial regularization weighting
* @param lambda_end final regularization weighting (for continuation)
* @param N size of input, x
* @param vops vector ops definition
* @param op linear operator, e.g. A
......@@ -254,22 +193,23 @@ void ist(unsigned int maxiter, float epsilon, float tau,
* @param x initial estimate
* @param b observations
*/
void fista(unsigned int maxiter, float epsilon, float tau,
float continuation, bool hogwild,
void fista(unsigned int maxiter, float epsilon, float tau,
long N,
const struct vec_iter_s* vops,
ist_continuation_t ist_continuation,
struct iter_op_s op,
struct iter_op_p_s thresh,
float* x, const float* b,
struct iter_monitor_s* monitor)
{
struct iter_data itrdata = {
struct ist_data itrdata = {
.rsnew = 1.,
.rsnot = 1.,
.iter = 0,
.maxiter = maxiter,
.tau = tau,
.scale = 1.,
};
float* r = vops->allocate(N);
......@@ -280,24 +220,15 @@ void fista(unsigned int maxiter, float epsilon, float tau,
itrdata.rsnot = vops->norm(N, b);
float ls_old = 1.;
float lambda_scale = 1.;
int hogwild_k = 0;
int hogwild_K = 10;
for (itrdata.iter = 0; itrdata.iter < maxiter; itrdata.iter++) {
iter_monitor(monitor, vops, x);
ls_old = lambda_scale;
lambda_scale = ist_continuation(&itrdata, continuation);
if (lambda_scale != ls_old)
debug_printf(DP_DEBUG3, "##lambda_scale = %f\n", lambda_scale);
if (NULL != ist_continuation)
ist_continuation(&itrdata);
iter_op_p_call(thresh, lambda_scale * tau, x, x);
iter_op_p_call(thresh, itrdata.scale * itrdata.tau, x, x);
ravine(vops, N, &ra, x, o); // FISTA
iter_op_call(op, r, x); // r = A x
......@@ -310,18 +241,7 @@ void fista(unsigned int maxiter, float epsilon, float tau,
if (itrdata.rsnew < epsilon)
break;
vops->axpy(N, x, tau, r);
if (hogwild)
hogwild_k++;
if (hogwild_k == hogwild_K) {
hogwild_K *= 2;
hogwild_k = 0;
tau /= 2;
}
vops->axpy(N, x, itrdata.tau, r);
}
debug_printf(DP_DEBUG3, "\n");
......@@ -342,6 +262,7 @@ void landweber(unsigned int maxiter, float epsilon, float alpha, long N, long M,
struct iter_op_s op,
struct iter_op_s adj,
float* x, const float* b,
struct iter_op_s callback,
struct iter_monitor_s* monitor)
{
float* r = vops->allocate(M);
......@@ -365,6 +286,9 @@ void landweber(unsigned int maxiter, float epsilon, float alpha, long N, long M,
iter_op_call(adj, p, r);
vops->axpy(N, x, alpha, p);
if (NULL != callback.fun)
iter_op_call(callback, x, x);
}
vops->del(r);
......@@ -470,10 +394,14 @@ cleanup:
* Iteratively Regularized Gauss-Newton Method
* (Bakushinsky 1993)
*
* y = F(x) = F x0 + DF dx + ...
* y = F(x) = F xn + DF dx + ...
*
* IRGNM: DF^H ((y - F x_0) + DF (xn - x0)) = ( DF^H DF + alpha ) (dx + xn - x0)
* DF^H ((y - F x_0)) - alpha (xn - x0) = ( DF^H DF + alpha) dx
* IRGNM: DF^H ((y - F xn) + DF (xn - x0)) = ( DF^H DF + alpha ) (dx + xn - x0)
* DF^H ((y - F xn)) - alpha (xn - x0) = ( DF^H DF + alpha) dx
*
* This version only solves the second equation for the update 'dx'. This corresponds
* to a least-squares problem where the quadratic regularization applies to the difference
* to 'x0'.
*/
void irgnm(unsigned int iter, float alpha, float alpha_min, float redu, long N, long M,
const struct vec_iter_s* vops,
......@@ -481,7 +409,8 @@ void irgnm(unsigned int iter, float alpha, float alpha_min, float redu, long N,
struct iter_op_s adj,
struct iter_op_p_s inv,
float* x, const float* xref, const float* y,
struct iter_op_s callback)
struct iter_op_s callback,
struct iter_monitor_s* monitor)
{
float* r = vops->allocate(M);
float* p = vops->allocate(N);
......@@ -489,7 +418,7 @@ void irgnm(unsigned int iter, float alpha, float alpha_min, float redu, long N,
for (unsigned int i = 0; i < iter; i++) {
// printf("#--------\n");
iter_monitor(monitor, vops, x);
iter_op_call(op, r, x); // r = F x
......@@ -520,6 +449,67 @@ void irgnm(unsigned int iter, float alpha, float alpha_min, float redu, long N,
}
/**
* Iteratively Regularized Gauss-Newton Method
* (Bakushinsky 1993)
*
* y = F(x) = F xn + DF dx + ...
*
* IRGNM: R(DF^H, DF^H DF, alpha) ((y - F xn) + DF (xn - x0)) = (dx + xn - x0)
*
* This version has an extra call to DF, but we can use a generic regularized
* least-squares solver.
*/
void irgnm2(unsigned int iter, float alpha, float alpha_min, float alpha_min0, float redu, long N, long M,
const struct vec_iter_s* vops,
struct iter_op_s op,
struct iter_op_s der,
struct iter_op_p_s lsqr,
float* x, const float* xref, const float* y,
struct iter_op_s callback,
struct iter_monitor_s* monitor)
{
float* r = vops->allocate(M);
float* q = vops->allocate(M);
for (unsigned int i = 0; i < iter; i++) {
iter_monitor(monitor, vops, x);
iter_op_call(op, r, x); // r = F x
vops->xpay(M, -1., r, y); // r = y - F x
debug_printf(DP_DEBUG2, "Step: %u, Res: %f\n", i, vops->norm(M, r));
if (NULL != xref)
vops->axpy(N, x, -1., xref);
iter_op_call(der, q, x);
vops->axpy(M, r, +1., q);
iter_op_p_call(lsqr, alpha, x, r);
if (NULL != xref)
vops->axpy(N, x, +1., xref);
alpha = (alpha - alpha_min) / redu + alpha_min;
if (alpha < alpha_min0)
alpha = alpha_min0;
if (NULL != callback.fun)
iter_op_call(callback, x, x);
}
vops->del(q);
vops->del(r);
}
/**
* Alternating Minimzation
*
......
......@@ -14,6 +14,7 @@
#endif
#include "misc/types.h"
#include "misc/nested.h"
struct vec_iter_s;
......@@ -67,7 +68,8 @@ float conjgrad(unsigned int maxiter, float l2lambda, float epsilon,
long N,
const struct vec_iter_s* vops,
struct iter_op_s linop,
float* x, const float* b, struct iter_monitor_s* monitor);
float* x, const float* b,
struct iter_monitor_s* monitor);
void landweber(unsigned int maxiter, float epsilon, float alpha,
......@@ -76,6 +78,7 @@ void landweber(unsigned int maxiter, float epsilon, float alpha,
struct iter_op_s op,
struct iter_op_s adj,
float* x, const float* b,
struct iter_op_s callback,
struct iter_monitor_s* monitor);
void landweber_sym(unsigned int maxiter, float epsilon, float alpha,
......@@ -85,19 +88,44 @@ void landweber_sym(unsigned int maxiter, float epsilon, float alpha,
float* x, const float* b,
struct iter_monitor_s* monitor);
/**
* Store information about iterative algorithm.
* Used to flexibly modify behavior, e.g. continuation
*
* @param rsnew current residual
* @param rsnot initial residual
* @param iter current iteration
* @param maxiter maximum iteration
* @param tau tau
* @param scale scaling of regularization
*/
struct ist_data {
double rsnew;
double rsnot;
unsigned int iter;
const unsigned int maxiter;
float tau;
float scale;
};
typedef void CLOSURE_TYPE(ist_continuation_t)(struct ist_data* itrdata);
void ist(unsigned int maxiter, float epsilon, float tau,
float continuation, _Bool hogwild,
long N,
const struct vec_iter_s* vops,
ist_continuation_t ist_continuation,
struct iter_op_s op,
struct iter_op_p_s thresh,
float* x, const float* b,
struct iter_monitor_s* monitor);
void fista(unsigned int maxiter, float epsilon, float tau,
float continuation, _Bool hogwild,
long N,
const struct vec_iter_s* vops,
ist_continuation_t ist_continuation,
struct iter_op_s op,
struct iter_op_p_s thresh,
float* x, const float* b,
......@@ -111,7 +139,18 @@ void irgnm(unsigned int iter, float alpha, float alpha_min, float redu,
struct iter_op_s adj,
struct iter_op_p_s inv,
float* x, const float* x0, const float* y,
struct iter_op_s callback);
struct iter_op_s callback,
struct iter_monitor_s* monitor);
void irgnm2(unsigned int iter, float alpha, float alpha_min, float alpha0, float redu,
long N, long M,
const struct vec_iter_s* vops,
struct iter_op_s op,
struct iter_op_s der,
struct iter_op_p_s lsqr,
float* x, const float* xref, const float* y,
struct iter_op_s callback,
struct iter_monitor_s* monitor);
void altmin(unsigned int iter, float alpha, float redu,
long N,
......
......@@ -76,6 +76,54 @@ static bool checkeps(float eps)
}
static bool check_ops(long size,
const struct operator_s* normaleq_op,
unsigned int D,
const struct operator_p_s* prox_ops[D],
const struct linop_s* ops[D])
{
if (NULL != normaleq_op) {
auto dom = operator_domain(normaleq_op);
if (size != 2 * md_calc_size(dom->N, dom->dims)) // FIXME: too weak
return false;
auto cod = operator_codomain(normaleq_op);
if (size != 2 * md_calc_size(cod->N, cod->dims)) // FIXME: too weak
return false;
}
for (unsigned int i = 0; i < D; i++) {
long cosize = size;
if ((NULL != ops) && (NULL != ops[i])) {
auto dom = linop_domain(ops[i]);
if (size != 2 * md_calc_size(dom->N, dom->dims)) // FIXME: too weak
return false;
auto cod = linop_codomain(ops[i]);
cosize = 2 * md_calc_size(cod->N, cod->dims);
}
if ((NULL != prox_ops) && (NULL != prox_ops[i])) {
auto dom2 = operator_p_domain(prox_ops[i]);
if (cosize != 2 * md_calc_size(dom2->N, dom2->dims))
return false; // FIXME: too weak
}
}
return true;
}
void iter2_conjgrad(iter_conf* _conf,
const struct operator_s* normaleq_op,
unsigned int D,
......@@ -92,6 +140,8 @@ void iter2_conjgrad(iter_conf* _conf,
assert(NULL == biases);
UNUSED(xupdate_op);
assert(check_ops(size, normaleq_op, D, prox_ops, ops));
auto conf = CAST_DOWN(iter_conjgrad_conf, _conf);
float eps = md_norm(1, MD_DIMS(size), image_adj);
......@@ -127,6 +177,8 @@ void iter2_ist(iter_conf* _conf,
#endif
UNUSED(xupdate_op);
assert(check_ops(size, normaleq_op, D, prox_ops, ops));
auto conf = CAST_DOWN(iter_ist_conf, _conf);
float eps = md_norm(1, MD_DIMS(size), image_adj);
......@@ -134,10 +186,15 @@ void iter2_ist(iter_conf* _conf,
if (checkeps(eps))
goto cleanup;
assert((conf->continuation >= 0.) && (conf->continuation <= 1.));
// This was probably broken for IST until v0.4.04
// better turn of it off with an error
assert(1 == conf->continuation);
ist(conf->maxiter, eps * conf->tol, conf->INTERFACE.alpha * conf->step, conf->continuation, conf->hogwild, size, select_vecops(image_adj),
OPERATOR2ITOP(normaleq_op), OPERATOR_P2ITOP(prox_ops[0]), image, image_adj, monitor);
// Let's see whether somebody uses it...
assert(!conf->hogwild);
ist(conf->maxiter, eps * conf->tol, conf->INTERFACE.alpha * conf->step, size, select_vecops(image_adj),
NULL, OPERATOR2ITOP(normaleq_op), OPERATOR_P2ITOP(prox_ops[0]), image, image_adj, monitor);
cleanup:
;
......@@ -167,15 +224,43 @@ void iter2_fista(iter_conf* _conf,
float eps = md_norm(1, MD_DIMS(size), image_adj);
assert(check_ops(size, normaleq_op, D, prox_ops, ops));
if (checkeps(eps))
goto cleanup;
return; // clang limitation
// goto cleanup;
assert((conf->continuation >= 0.) && (conf->continuation <= 1.));
fista(conf->maxiter, eps * conf->tol, conf->INTERFACE.alpha * conf->step, conf->continuation, conf->hogwild, size, select_vecops(image_adj),
OPERATOR2ITOP(normaleq_op), OPERATOR_P2ITOP(prox_ops[0]), image, image_adj, monitor);
__block int hogwild_k = 0;
__block int hogwild_K = 10;
cleanup:
NESTED(void, continuation, (struct ist_data* itrdata))
{
float a = logf(conf->continuation) / (float)itrdata->maxiter;
itrdata->scale = expf(a * itrdata->iter);
if (conf->hogwild) {
/* this is not exactly identical to what was implemented
* before as tau is now reduced at the beginning. But this
* seems more correct. */
hogwild_k++;
if (hogwild_k == hogwild_K) {
hogwild_k = 0;
hogwild_K *= 2;
itrdata->tau /= 2;
}
}
};
fista(conf->maxiter, eps * conf->tol, conf->INTERFACE.alpha * conf->step, size, select_vecops(image_adj),
continuation, OPERATOR2ITOP(normaleq_op), OPERATOR_P2ITOP(prox_ops[0]), image, image_adj, monitor);
// cleanup:
;
}
......@@ -207,6 +292,8 @@ void iter2_chambolle_pock(iter_conf* _conf,
assert((long)md_calc_size(iv->N, iv->dims) * 2 == size);
assert(check_ops(size, normaleq_op, D, prox_ops, ops));
// FIXME: sensible way to check for corrupt data?
#if 0
float eps = md_norm(1, MD_DIMS(size), image_adj);
......@@ -240,6 +327,8 @@ void iter2_admm(iter_conf* _conf,
{
auto conf = CAST_DOWN(iter_admm_conf, _conf);
assert(check_ops(size, normaleq_op, D, prox_ops, ops));
struct admm_plan_s admm_plan = {
.maxiter = conf->maxiter,
......@@ -322,6 +411,8 @@ void iter2_pocs(iter_conf* _conf,
UNUSED(xupdate_op);
UNUSED(image_adj);
assert(check_ops(size, normaleq_op, D, prox_ops, ops));
struct iter_op_p_s proj_ops[D];
for (unsigned int i = 0; i < D; i++)
......
......@@ -41,94 +41,15 @@ const struct iter3_irgnm_conf iter3_irgnm_defaults = {
.nlinv_legacy = false,
};
const struct iter3_landweber_conf iter3_landweber_defaults = {
struct irgnm_s {
.INTERFACE.TYPEID = &TYPEID2(iter3_landweber_conf),
INTERFACE(iter_op_data);
struct iter_op_s der;
struct iter_op_s adj;
float* tmp;
long size;
int cgiter;
float cgtol;
bool nlinv_legacy;
.iter = 8,
.alpha = 1.,
.epsilon = 0.1,
};
DEF_TYPEID(irgnm_s);
static void normal(iter_op_data* _data, float* dst, const float* src)
{
auto data = CAST_DOWN(irgnm_s, _data);
iter_op_call(data->der, data->tmp, src);
iter_op_call(data->adj, dst, data->tmp);
}
static void inverse(iter_op_data* _data, float alpha, float* dst, const float* src)
{
auto data = CAST_DOWN(irgnm_s, _data);
md_clear(1, MD_DIMS(data->size), dst, FL_SIZE);
float eps = data->cgtol * md_norm(1, MD_DIMS(data->size), src);
/* The original (Matlab) nlinv implementation uses
* "sqrt(rsnew) < 0.01 * rsnot" as termination condition.
*/
if (data->nlinv_legacy)
eps = powf(eps, 2.);
conjgrad(data->cgiter, alpha, eps, data->size, select_vecops(src),
(struct iter_op_s){ normal, CAST_UP(data) }, dst, src, NULL);
}
void iter3_irgnm(iter3_conf* _conf,
struct iter_op_s frw,
struct iter_op_s der,
struct iter_op_s adj,
long N, float* dst, const float* ref,
long M, const float* src,
struct iter_op_s cb)
{
auto conf = CAST_DOWN(iter3_irgnm_conf, _conf);
float* tmp = md_alloc_sameplace(1, MD_DIMS(M), FL_SIZE, src);
struct irgnm_s data = { { &TYPEID(irgnm_s) }, der, adj, tmp, N, conf->cgiter, conf->cgtol, conf->nlinv_legacy };
irgnm(conf->iter, conf->alpha, conf->alpha_min, conf->redu, N, M, select_vecops(src),
frw, adj,
(struct iter_op_p_s){ inverse, CAST_UP(&data) },
dst, ref, src, cb);
md_free(tmp);
}
void iter3_landweber(iter3_conf* _conf,
struct iter_op_s frw,
struct iter_op_s der,
struct iter_op_s adj,
long N, float* dst, const float* ref,
long M, const float* src)
{
auto conf = CAST_DOWN(iter3_landweber_conf, _conf);
assert(NULL == der.fun);
assert(NULL == ref);
float* tmp = md_alloc_sameplace(1, MD_DIMS(N), FL_SIZE, src);
landweber(conf->iter, conf->epsilon, conf->alpha, N, M,
select_vecops(src), frw, adj, dst, src, NULL);
md_free(tmp);
}
......@@ -5,24 +5,11 @@
*/
#include "misc/types.h"
typedef struct iter3_conf_s { TYPEID* TYPEID; } iter3_conf;
struct iter_op_s;
typedef void iter3_fun_f(iter3_conf* _conf,
struct iter_op_s frw,
struct iter_op_s der,
struct iter_op_s adj,
long N, float* dst, const float* ref,
long M, const float* src);
typedef void iter3_irgnm_f(iter3_conf* _conf,
struct iter_op_s frw,
struct iter_op_s der,
struct iter_op_s adj,
long N, float* dst, const float* ref,
long M, const float* src,
struct iter_op_s cb);
struct iter3_irgnm_conf {
......@@ -41,7 +28,6 @@ struct iter3_irgnm_conf {
};
iter3_irgnm_f iter3_irgnm;
......@@ -55,9 +41,8 @@ struct iter3_landweber_conf {
};
iter3_fun_f iter3_landweber;
extern const struct iter3_irgnm_conf iter3_irgnm_defaults;
// extern const struct iter3_landweber_conf iter3_landweber_defaults;
extern const struct iter3_landweber_conf iter3_landweber_defaults;
......@@ -3,6 +3,9 @@
* a BSD-style license which can be found in the LICENSE file.
*/
#include <assert.h>
#include <math.h>
#include "num/ops.h"
#include "num/multind.h"
#include "num/flpmath.h"
......@@ -18,8 +21,9 @@
#include "iter/italgos.h"
#include "iter/vec.h"
#include "iter/iter3.h"
#include "iter/iter2.h"