[Ether]

Originally Posted by Ether

...the following two computational methods yield virtually identical results for min L2 norm of b-Ax:


Method 1

1a) [U,S,V] = svd(A)

1b) Replace the diagonal elements of S with their reciprocals, except when abs(Sjj)<threshold (I used 1e-4 for threshold), in which case make Sjj zero.

1c) compute x = V*S*(U'*b)


Method 2

2a) N = A'*A

2b) d= A'*b

2c) compute x = N\d ..... (Octave mldivide notation)

2d) compute m = mean(x)

2e) subtract m from each element of x


Notice Method 1 factors A, not A'A, resulting in less rounding error.

There's a simpler way to do Method#1 above if you are using Matlab or Octave (hat tip to AGPapa):

x = pinv(A,tol)*b;

pinv() is explained in detail here:

http://www.mathworks.com/help/matlab/ref/pinv.html

(well worth reading; explains the interesting difference between x₁=pinv(A) and x₂=A\b)