Least-square inversion with inexact adjoints. Method of conjugate directions: A tutorial

PROGRAM

The program in Table 1 implements one iteration of the conjugate-direction method. It is based upon Jon Claerbout's cgstep() program (, ) and uses an analogous naming convention. Vectors in the data space are denoted by double letters.

void sf_cdstep(bool forget     /* restart flag */, 
	       int nx          /* model size */, 
	       int ny          /* data size */, 
	       float* x        /* current model [nx] */, 
	       const float* g  /* gradient [nx] */, 
	       float* rr       /* data residual [ny] */, 
	       const float* gg /* conjugate gradient [ny] */) 
/*< Step of conjugate-direction iteration. 
  The data residual is rr = A x - dat
>*/
{
    float *s, *si, *ss;
    double alpha, beta;
    int i, n, ix, iy;

    s = sf_floatalloc(nx+ny);
    ss = s+nx;

    for (ix=0; ix < nx; ix++) {  s[ix] = g[ix]; }
    for (iy=0; iy < ny; iy++) { ss[iy] = gg[iy]; }

    sf_llist_rewind (steps);
    n = sf_llist_depth (steps);

    for (i=0; i < n; i++) {
	sf_llist_down (steps, &si, &beta);
	alpha = - cblas_dsdot(ny, gg, 1, si+nx, 1) / beta;
	cblas_saxpy(nx+ny,alpha,si,1,s,1);
    }
    
    beta = cblas_dsdot(ny, s+nx, 1, s+nx, 1);
    if (beta < DBL_EPSILON) return;

    sf_llist_add (steps, s, beta);
    if (forget) sf_llist_chop (steps);
    alpha = -  cblas_dsdot(ny, rr, 1, ss, 1) / beta;
    
    cblas_saxpy(nx,alpha,s, 1,x, 1);
    cblas_saxpy(ny,alpha,ss,1,rr,1);
}

Table 1. The source of this program is RSF/api/c/cdstep.c

In addition to the previous steps ${\bf s}_j$ and their conjugate counterparts ${\bf A s}_j$ (array s), the program stores the squared norms $\Vert{\bf A s}_j\Vert^2$ (variable beta) to avoid recomputation. For practical reasons, the number of remembered iterations can actually be smaller than the total number of iterations.

Least-square inversion with inexact adjoints. Method of conjugate directions: A tutorial