3

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!

Improve this question
1
  • Can you compute A^T A ? If so, solve (A^T A + tiny I) x = A^T b with tiny say 1e-12 -- that helps when A^T A is nearly singular, see Tikhonov_regularization . Commented Nov 12, 2014 at 16:24

2 Answers 2

Reset to default
4

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

Improve this answer
Sign up to request clarification or add additional context in comments.

2 Comments

Thanks for your answer. I know a priori that the matrix is exactly determined. I should have included in my original post that the first method I tried was to use SciPy's sparse least-squares solvers! My worry, however, is that least squares is a very inefficient, iterative approach (involving gradually minimizing some cost function to some input tolerance) to solve something that could instead be solved exactly and much more efficiently using a direct linear solver. But I don't know how to e.g. take my matrix A and transform it somehow so I can feed it to scipy.sparse.linalg.spsolve...
I added the potential use of normal equations to the answer. This makes the system square and is equivalent to your problem if the matrix has full rank, as you describe. From this setting you can use spsolve or conjugate gradient to find a solution to the linear system
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.])
Improve this answer

3 Comments

Definitely correct, but there are tools to do this in scipy. scipy.linalg.lstsq does this decomposition if I am not mistaken and provides the corresponding solution (minimal l2 norm if underdetermined, otherwise closest solution in least squares sense). The most important reason not to decompose is that @JoshBurkart's matrix, according to him, is sparse (I conclude that it is potentially enormous). So you don't want to do O(nk^2) operations to get and SVD before solving for one single target variable. scip.sparse.lsqr solves iteratively for the specified target.
Update: scipy.linalg.lstsq calls a lapack function called gelss which "Computes the minimum norm least squares solution to an over- or under-determined system of linear equations AX=B, using the singular value decomposition of A." (source)
Thanks for your response. Unfortunately, as @eickenberg already commented, my matrices are huge and the SVD is a prohibitively expensive operation.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.