Pseudo-inverse without full SVD

[SebWouters]

I'm doing a dive into numerical linear algebra, since we need quite a few things in a memory safe language at our company.

One of the things we need is the pseudo-inverse solution of the system A * x = b, where A is an m x n matrix (m >=n and m < n are relevant) and b an m x 1 vector. For sand-boxing, I used the non-memory safe interface of mathdotnet numerics to one of the lapack implementations, which actually computes the full pseudo-inverse via SVD and then applies it to b.

However, when you only need the pseudo-inverse solution, a couple of things can be improved:

1. You don't need memory for U of A = U * s * VT, since you can apply the pieces which contribute to U directly into b.
2. You don't need to get the actual singular values, you just need to single out the null-space of A, and make sure x doesn't have components in this null space.

Consider m >= n, you could do A = Q * R, detect that all diagonal elements of R are numerically nonzero, and solve the system R * x = Q^T * b. Of course you do that without storing Q, but multiplying the pieces which contribute directly into b. This doesn't work when there's a numerical zero on the diagonal, even if you just make x[current row] = 0.0, since you then typically have components in your solution x which are in the nullspace of A.

You could work via bidiagonalization (for m>=n with upper diagonal):

A = Q1 * R2 * Q2

You apply the pieces of Q1 into b directly, but are required to keep Q2 either as such, or implicitly (Householder vecs...). If there's no numerical zero on the diagonal of R2, you can solve R2 * y = Q1^T * b easily, and obtain x = Q2^T * y. If there is a numerical zero, you can divide and conquer R2:

[ a b        ]
[   c d      ]
[     0 e    ]
[       f g  ]
[         h  ]

by flipping the off-diagonal elements to the other side, separating the bidiagonal matrix into disjoint blocks. This can be done in the example by post-multiplying with Givens rotations for column pairs (0, 1) and (1, 2), and pre-multiplying with Givens rotations for row pairs (3, 4) and (2, 3) [using counting from 0]. We then have A = Q3 * R3 * Q4, with R3:

[ p          ]
[ q r        ]
[     0      ]
[       s    ]
[       t u  ]

In the example, for the system R3 * y = Q3^T * b, y[2] is obviously 0. This correctly ensures there's no contribution of the null-space of A in x. The blocks before and after the zero can be dealt with separately.

Do you now of a routine in this lapack implementation, or any other lapack implementation, which performs pseudo-inverse application more efficient than full SVD building? If not, would this be something useful we could add to reference lapack?

[langou]

To solve overdetermined ("m >= n", least squares) or underdetermined ("m <= n", min "x" 2-norm) linear systems, you have to look at the GELS family.

Arbritrarily using D as the prefix, you have DGELS, DGELSY, DGELSS and DGELSD.

If the coefficient matrix is full, then have a look at DGELS. This is based on a straight up QR factorization. The method breaks if A is rank-deficient and the method does not perform well if A is "numerically rank-deficient".

If the coefficient matrix is rank-deficient or "numerically rank-deficient", then have a look at DGELSY, DGELSS or DGELSD.

DGELSY is based on the complete orthogonal decomposition of A. (So based on a QR factorization with pivoting, and then COD.)

The method you are describing is most related to DGELSD. DGELSD reduces to bidiagonal form, and then solves the bidiagonal linear least squares problem (or min x 2-norm problem or a combination of both).

DGELSS is based on the SVD of A.

[SebQonic]

Thanks!