Exercise 9: Singular Value Decomposition 
Before you actually start working on this Exercise, it is essential that you have read Chapter 2 of Numerical Recipes (particularly the section Singular Value Decomposition), and that you understand the various concepts mentioned there (rank, range, orthogonality, etc.). That will make it much easier to understand what is going on when the various control parameters are changed in the experiments below.

Subsections  

Preparation section  15 min 
SVD experimentarium  
Improving SVD  60 min 
Setting up the matrix system from scratch  25 min 
Home Work 
CVS 
[about 1 minute] 
To extract the exercise files for this weeks exercise, do
cd ~/ComputerPhysics cvs update d
In case of problems, see the CVS update help page.
SVD preparatory section 
[about 15 minutes] 
Before you start playing with the SVD experimentarium, let's check that you picked up enough from the lecture and the Numerical Recipes to make the exercise meaningful.
The following is a short summary, combined with questions. If you have prepared enough, you should have no problem with the questions. Conversely, if you find the summary too condensed and incomprehensible, you probably have not prepared enough.
The basic theorem on which the SVD method is based says that an arbitrary MxN (M rows, N columns) matrix A may be decomposed into a product of three matrices, denoted U, W, and V^{t} in the Lecture Notes and in Numerical Recipes:
A = U W V^{t}



The linear equation system
A x = b
may be interpreted as a mapping from an Ndimensional xspace to an Mdimensional bspace, where we assume that M is larger than N. Using the SVD decomposition of A into U W V^{t}, we may write A x as U ( W ( V^{t} x)), where the innermost operation is a generalized rotation in Nspace, the next one is a scaling of each of the N components, and the last one is a linear combination of the N column vectors in U, each with M components.
Since, according to this formula, b is a linear combination of the column vectors of U, it can never lie outside the Ndimensional space spanned by these. But b has M (> N) components; hence it formally belongs to an Mdimensional space.
Which of the following statements are true?

If M>N, an Mdimensional vector has more components than a vector expressed as components of the N orthogonal column vectors in U. Hence it is possible to have b vectors for which the equation system has no solution; there is no x that will produce a b outside the range of A. However, the expression for the inversion of the equation,
x = A^{1} b = (V (W^{1} (U^{t }b)))
may still be formally applied, and indeed still produces a useful solution; the solution for which the remainder (A x  b) is "smallest", in the sense that it has the smallest sum of squares of its components. To see that this is indeed the case, consider the following:
The three successive steps indicated by the parentheses may be interpreted as finding the components of b along the N column vectors of U, then scaling the components, and finally rotating the answer back to the original coordinate system.
Which of the following statements are true?




Each column vector in U is multiplied with a singular value w_{i} when going from x to b. It may happen, that some changes of x have practically no effect on b (think about b as a measurement, and x the "true" function  a limited resolution in the measurement means that details of x can hardly be "seen" in b).


Thus, if b is obtained from measurements, and is influenced by noise in the measurements, then the components of the noise along those particular U column vectors get multiplied with very large coefficients, which leads to huge uncertainties in the answer.
At this point, it should be easy to understand that one is better off ignoring these components completely than by blowing up their accidental values with large numbers. Hence the rule to put 1/w_{i} = 0 if w_{i} is too small.
In the experimentarium, synthetic measurements b are produced by first smearing the input profile x with a (slightly noisy) Gaussian profile, and then adding some noise to the "measurements".
SVD experimentarium 
[about 30 minutes] 
To run the SVD experimentarium:





Improving the SVD solution 
[about 60 minutes] 
The SVD solution is much better than the LUdecomposition (less noise sensitive), but it shows some ringing near the sharp corners of the original function.
solve.pro
),
which tapers off the amplitude
smoothly instead of abruptly, with the sharpness contolled by a
parameter p
for j=0,n1 do wp[j,j] = 1./(w[j]^data.p+(max(w)*data.eps/(w[j]+1e30))^data.p)^(1./data.p)
which has a sharper dropoff for larger p. To try this out:
In order to pass the test below your xsvd.pro needs to have:


Setting up the equations from scratch 
[about 25 minutes] 
In the previous part of this exercise you became familiar with SVD using a prebuilt setup. In this part you perform the basic calculations yourself, starting from a series of data files and information about how these relate to the final result.
where r_{j} are values for the Earth's reflectivity. This may be written as a linear system of equations:
where s is given by the data in the file response.dat and the matrix A_{i,j} = a_{ij} may be constructed from the profile given in signal.dat.
The reflectivity as a function of depth, r, can now be determined using SVD:
IDL> journal,'svd.jou'
IDL> openr,1,'signal.dat' & Na=25 & a=fltarr(Na) & readf,1,a & close,1
From the formula above you can convince yourself that A becomes a "lower band diagonal" square matrix (A_{i,j} is nonzero only for band of values with j≤i). Assign values to it:
IDL> aa=fltarr(Ns,Ns) ; make the 2D array IDL> for i=0,Ns1 do for j=(iNa+1)>0,i do aa[i,j]=a[ij] ; set nonzero values
IDL> ? svdc
The use of SVDC (as well as other routines that have been imported to IDL from the Numerical Recipes library) can be confusing, in that the matrix storage order is opposite to the normal mathematical one, unless the /COLUMN keyword is used. To follow standard math [row,column] matrix notation one should thus set the /COLUMN keyword. For a more detailed discussion of matrix storage order, search for "Columns, Rows" in IDLDE Help. 
IDL> ww=fltarr(Ns,Ns) & for i=0,Ns1 do if w[i] gt max(w)*threshold then ww[i,i]=1./w[i]
IDL> r1=v#ww#transpose(u)#s
IDL matrix multiplication with #operators works just like in math; it takes rows from the matrix on the left (running over the 2nd matrix index), multiplies with columns from the matrix to the right (running over the first matrix index), and puts the result at the 'intersection' of the row and column in question. So, the line above corresponds exactly to what you can find in Numerical Recipes and in the Lecture Notes. (The ##notation used in the IDL help file for SVDC works with the opposite matrix notation [column,row] that is default when the /COLUMN keyword is NOT used.) For a more detailed discussion of IDL matrix operations search for "matrix operators" in IDLDE Help. 
Compute the rms value of the difference between your derived solution and the one from the data file.
IDL> openr,1,'reflectivity.dat' & r2=fltarr(Ns) & readf,1,r2 & close,1 IDL> rms=sqrt(total((r1r2)^2)/Nr) IDL> print,'RMS error: ',rms
IDL> journalOpen the journal file in an editor (e.g. the IDLDE editor) and remove lines related to mistakes and failed attempts, such that the edited journal file can be used to reexecute the whole thing. Try it out with
IDL> @svd.jou
In order to pass the test below your svd.jou journal file must assign the RMS value to an IDL variable 'rms'. 

Use any remaining time for example to play with the Xwidget interface, either with this weeks interface, or the one from last week.
Home Work 
[about 2 hours] 
Spend at least one hour this week reading the chapter about Linear Algebra in Numerical Recipes; understanding what goes on in this exercise is nontrivial! Details about what to read are given in the Lecture Notes.
If you did not have enough time to play with IDL last week, or during the exercise hours, then do it in your home work time. The Xwidget interface this week is similar to the one last week; in fact if you would like to experiment with modifications of the Xwidgets, you can pick either one of the two exercises.
If you want to prepare for the next Lecture, then check out the Chapter in Numerical Recipes about Ordinary Differential Equations.