2
$\begingroup$

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 *) enter image description here

notebook

$\endgroup$
  • 5
    $\begingroup$ Try LeastSquares. $\endgroup$ – Henrik Schumacher Oct 14 '18 at 22:32
  • $\begingroup$ SVD ought to work here too, if you're concerned about saving factorization effort. $\endgroup$ – J. M. will be back soon Oct 15 '18 at 3:15
  • $\begingroup$ @J.M. I doubt that SingularValueDecomposition will be faster than LeastSquares... $\endgroup$ – Henrik Schumacher Oct 16 '18 at 12:34
  • $\begingroup$ @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... $\endgroup$ – J. M. will be back soon Oct 16 '18 at 12:42
3
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – J. M. will be back soon Oct 16 '18 at 13:08
  • 1
    $\begingroup$ For the normal equations approach, one can do LinearSolve[A.Transpose[A], Method -> "Cholesky"] to exploit the SPD structure of the problem. $\endgroup$ – J. M. will be back soon Oct 16 '18 at 13:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.