# Getting least norm solution

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

• Try LeastSquares. – Henrik Schumacher Oct 14 '18 at 22:32
• SVD ought to work here too, if you're concerned about saving factorization effort. – J. M. will be back soon Oct 15 '18 at 3:15
• @J.M. I doubt that SingularValueDecomposition will be faster than LeastSquares... – Henrik Schumacher Oct 16 '18 at 12:34
• @Henrik: a test would be needed to be sure, but computing the SVD once, and taking the pseudoinverse of the diagonal factor would be the only precomputation needed. Unless LeastSquares[] actually stores a factorization internally instead of computing it for each right-hand-side... – J. M. will be back soon Oct 16 '18 at 12:42

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.

• Nicely done. Just as a note for everybody else: you can manually get the Moore-Penrose inverse of a singular diagonal matrix by inverting all the nonzero diagonal elements and leaving the zero ones untouched. In the inexact arithmetic case, one has to decide on how big a diagonal element must be before one takes its reciprocal, which is efffectively how one decides numerical rank. – J. M. will be back soon Oct 16 '18 at 13:08
• For the normal equations approach, one can do LinearSolve[A.Transpose[A], Method -> "Cholesky"] to exploit the SPD structure of the problem. – J. M. will be back soon Oct 16 '18 at 13:10