Solve seemingly (but not actually!) overdetermined sparse linear system in Python

Question

I have a sparse matrix A (using scipy.sparse) and a vector b, and want to solve Ax = b for x. A has more rows than columns, so it appears to be overdetermined; however, the rows of A are linearly dependent, so that in actuality the row rank of A is equal to the number of columns. For example, A could be

A = np.array([[1., 1.], [-1., -1.], [1., 0.]])

while b is

b = np.array([0., 0., 1.])

The solution is then x = [1., -1.]. I'm wondering how to solve this system in Python, using the functions available in scipy.sparse.linalg. Thanks!

Answer 1

Is your system possibly underdetermined? If it is not, and there is actually a solution, then the least squares solution will be that solution, so you can try

from scipy.sparse.linalg import lsqr
return_values = lsqr(A, b)
x = return_values[0]

If your system is actually underdetermined, this should find you the minimum L2 norm solution. If it doesn't work, set the parameter damp to something very small (e.g. 1e-5).

If your system is exactly determined (i.e. A is of full rank) and has a solution, and your matrix A is tall, as you describe it, then you can find an equivalent system in the normal equations:

A.T.dot(A).dot(x) == A.T.dot(b)

has a unique solution in x. This is a square linear system and is thus solvable using linear system solvers such as scipy.sparse.linalg.spsolve

Answer 2

The formally correct way of solving your problem is to use SVD. You have a system of the form

A [MxN] * x [Nx1] = b [Mx1]

The SVD decomposes the matrix A into three others, so you get:

U [MxM] * S[MxN] * V[N*N] * x[Nx1] = b[Mx1]

The matrices U and V are both orthogonal (their inverse is their transpose), and S is a diagonal matrix. If we rewrite the above we get:

S[MxN] * V [N * N] * x[Nx1] = U.T [MxM] * b [Mx1]

If M > N then the matrix S will have its last M - N rows full of zeros, and if your system is truly determined, then U.T b should also have the last M - N rows zero. That means that you can solve your system as:

>>> a = np.array([[1., 1.], [-1., -1.], [1., 0.]])
>>> b = np.array([0., 0., 1.])
>>> u, s, v = np.linalg.svd(a)
>>> np.allclose(u.T.dot(b)[-m+n:], 0) #check system is not overdetermined
True
>>> np.linalg.solve(s[:, None] * v, u.T.dot(b)[:n])
array([ 1., -1.])