Editing Guide to madagascar programs (section)

==sfconjgrad==
{| class="wikitable" align="center" cellspacing="0" border="1"
! colspan="4" style="background:#ffdead;" | Generic conjugate-gradient solver for linear inversion 
|-
! colspan="4" | sfconjgrad < dat.rsf mod=mod.rsf mwt=mwt.rsf known=known.rsf x0=x0.rsf > to.rsf < from.rsf > out.rsf niter=1
|-
| ''string '' || '''known=''' ||   || 	auxiliary input file name
|-
| ''file   '' || '''mod=''' ||   || 	auxiliary input file name
|-
| ''string '' || '''mwt=''' ||   || 	auxiliary input file name
|-
| ''int    '' || '''niter=1''' ||   || 	number of iterations
|-
| ''string '' || '''x0=''' ||   || 	auxiliary input file name
|}

<tt>sfconjgrad</tt> is a generic program for least-squares linear
inversion with the [http://en.wikipedia.org/wiki/Conjugate_gradient_method conjugate-gradient method]. Suppose you have an
executable program <tt><prog></tt> that takes an RSF file from the
standard input and produces an RSF file in the standard output. It may
take any number of additional parameters but one of them must be
<tt>adj=</tt> that sets the forward (<tt>adj=0</tt>) or adjoint
(<tt>adj=1</tt>) operations.  The program <tt><prog></tt> is typically
an RSF program but it could be anything (a script, a multiprocessor
MPI program, etc.) as long as it implements a linear operator
<math>\mathbf{L}</math> and its adjoint. There are no restrictions on the data
size or shape. You can easily test the adjointness with
<tt>sfdottest</tt>. The <tt>sfconjgrad</tt> program searches for a
vector <math>\mathbf{m}</math> that minimizes the least-square misfit 
<math>\|\mathbf{d - L\,m}\|^2</math> for the given input data vector <math>\mathbf{d}</math>. 

The pseudocode for <tt>sfconjgrad</tt> is given at the end of the [https://ahay.org/RSF/book/gee/lsq/paper.pdf "Model fitting with least squares" chapter] of ''Imaging Estimation by Example'' by Jon Claerbout, with the earliest form published in [http://sepwww.stanford.edu/data/media/public/oldreports/sep48/48_25.pdf "Conjugate Gradient Tutorial"] (SEP-48, 1986, same author). A simple toy implementation with a small matrix shows that this algorithm produces the same steps as the algorithm described in equations 45-49 of [http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf "An introduction to the Conjugate Gradient Method Without the Agonizing Pain"] by J.R. Shewchuk, 1994, when the equation <math>A^T A x = A^T b</math> (in Shewchuk's notation) is solved. Multiplying with the transpose ensures a correct solution even when matrix A is square but not symmetric at all. The program [https://ahay.org/RSF/sfcconjgrad.html sfcconjgrad] implements this algorithm for the case when inputs are complex.

Here is an example. The <tt>sfhelicon</tt> program implements Claerbout's multidimensional helical filtering
(Claerbout, 1998<ref>Claerbout, J.,  1998, Multidimensional recursive filters via a helix:  Geophysics, '''63''', 1532--1541.</ref>). It requires a filter to be specified in addition to the input and output vectors. We create a helical  2-D filter using the Unix <tt>echo</tt> command.
<pre>
bash$ echo 1 19 20 n1=3 n=20,20 data_format=ascii_int in=lag.rsf > lag.rsf
bash$ echo 1 1 1 a0=-3 n1=3 data_format=ascii_float in=flt.rsf > flt.rsf
</pre>
Next, we create an example 2-D model and data vector with <tt>sfspike</tt>.
<pre>
bash$ sfspike n1=50 n2=50 > vec.rsf
</pre>
The <tt>sfdottest</tt> program can perform the dot product test to
check that the adjoint mode works correctly.
<pre>
bash$ sfdottest sfhelicon filt=flt.rsf lag=lag.rsf mod=vec.rsf dat=vec.rsf
sfdottest:  L[m]*d=5.28394
sfdottest: L'[d]*m=5.28394
</pre>
Your numbers may be different because <tt>sfdottest</tt> generates new
random input on each run.
Next, let us make some random data with <tt>sfnoise</tt>.
<pre>
bash$ sfnoise seed=2005 rep=y < vec.rsf > dat.rsf
</pre>
and try to invert the filtering operation using <tt>sfconjgrad</tt>:
<pre>
bash$ sfconjgrad sfhelicon filt=flt.rsf lag=lag.rsf mod=vec.rsf < dat.rsf > mod.rsf niter=10
sfconjgrad: iter 1 of 10
sfconjgrad: grad=3253.65
sfconjgrad: iter 2 of 10
sfconjgrad: grad=289.421
sfconjgrad: iter 3 of 10
sfconjgrad: grad=92.3481
sfconjgrad: iter 4 of 10
sfconjgrad: grad=36.9417
sfconjgrad: iter 5 of 10
sfconjgrad: grad=18.7228
sfconjgrad: iter 6 of 10
sfconjgrad: grad=11.1794
sfconjgrad: iter 7 of 10
sfconjgrad: grad=7.26941
sfconjgrad: iter 8 of 10
sfconjgrad: grad=5.15945
sfconjgrad: iter 9 of 10
sfconjgrad: grad=4.23055
sfconjgrad: iter 10 of 10
sfconjgrad: grad=3.57495
</pre>
The output shows that, in 10 iterations, the norm of the gradient vector decreases by almost 1000. 
We can check the residual misfit before
<pre>
bash$ < dat.rsf sfattr want=norm
norm value = 49.7801
</pre>
and after
<pre>
bash$ sfhelicon filt=flt.rsf lag=lag.rsf < mod.rsf | sfadd scale=1,-1 dat.rsf | sfattr want=norm
norm value = 5.73563
</pre>
The misfit decreased by an order of magnitude in 10 iterations. Running the program for more iterations can improve the result.

An equivalent implementation for complex-valued inputs is [https://ahay.org/RSF/sfcconjgrad.html sfcconjgrad]. A lightweight Python implementation can be found in [https://github.com/ahay/src/blob/master/user/fomels/conjgrad.py $PYTHONPATH/rsf/conjgrad.py].