Getting least norm solution

Question

Noob question -- how should I get least norm solution in Mathematica for an under-constrained problem? Matrix is not full rank.

I could use pseudo-inverse, but inverting a matrix to get a single solution seems like an overkill, is there a more efficient approach?

ones = {{1, 1}}; H = Transpose[ones].ones; ones.PseudoInverse[H] (* gives 1/2, 1/2 *) LinearSolve[H, {1, 1}] (* gives 1, 0 *)

notebook

Answer 1

Since you have a low rank factorization for your actual matrix, you can speed this up by using PseudoInverse, SingularValueDecomposition, or LeastSquares on the low rank factor A. Here a model problem:

n = 3000;
m = 300;
A = RandomReal[{-1, 1}, {m, n}];
b = RandomReal[{-1, 1}, {n}];
H = Transpose[A].A;

Applying all proposed solutions to the matrix H:

x1 = LeastSquares[H, b]; // AbsoluteTiming // First
x2 = PseudoInverse[H].b; // AbsoluteTiming // First
x3 = Module[{U, Σ, V},
     {U, Σ, V} = SingularValueDecomposition[H];
     U.LeastSquares[SparseArray[Σ], b.V]
     ]; // AbsoluteTiming // First
Max[Abs[x1 - x2]]
Max[Abs[x1 - x3]]

3.79227

5.75832

6.10518

6.72036*10^-18

6.31378*10^-18

So, for a single solve, LeastSquares performs best. For multiple solves, PseudoInverse may be the better choice because it can be reused multiple times.

Now, exploiting the low rank factorization:

y1 = LeastSquares[A, LeastSquares[Transpose[A], b]]; // AbsoluteTiming // First
y2 = With[{P = PseudoInverse[A]}, P.(b.P)]; // AbsoluteTiming // First
y3 = Module[{U, Σ, V},
     {U, Σ, V} = SingularValueDecomposition[A];
     V.LeastSquares[Transpose[#].# &@SparseArray[Σ], b.V]
     ]; // AbsoluteTiming // First
Max[Abs[x1 - y1]]
Max[Abs[x1 - y2]]
Max[Abs[x1 - y3]]

0.069174

0.057212

0.37002

1.04626*10^-17

6.23416*10^-18

7.01852*10^-18

Bang! Same accuracy but almost two orders of magnitude faster (best in the second run compared to best in first run). Needless to say that the speed-up get better the smaller the the row count of A is.

If you know a priorily that A has full rank (and accuracy is not your major concern), then you can also use the following:

z1 = With[{S = LinearSolve[A.Transpose[A], Method -> "Cholesky"]},
    S[S[A.b]].A
    ]; // AbsoluteTiming // First
Max[Abs[z1 - x1]]

0.013487

7.94178*10^-18

However, this employs the infamous normal equations and if A badly conditioned, this can lead to severe errors.

PS.: No, LeastSquares is not the bottleneck in the applications of SingularValueDecomposition though it can be optimized away by ``inverting'' the diagonal of Σ by hand.