Solving square singular system?

Question

I have a square singular system A and vector b which is in range space of A. Because b is in range space of A and A is singular there are infinitely many solutions. Now what I want is some solution to Ax=b, not necessarily the minimum norm. As my system A is large I want to avoid svd based solutions(finding pseudo inverse). I was thinking of simply doing LU factorization and setting free variables to 0 but all the methods I tried failed to give me solution when A is singular.
I have tried scipy.linalg.solve but that requires system to be non-singular. I have also tried scipy lu_factor and lu_solve but while doing lu_factor it gave me runtime warning saying "Diagonal number %d is exactly zero. Singular matrix.".
So my question is this - Is there a way(using scipy) to find a solution to a singular system using LU factorization given b is in range space of A. Any suggestions are greatly appreciated. Thanks a lot.

Hmm, in https://stackoverflow.com/a/43426991/500207 I used the SVD to find the nullspace of a matrix (via http://scipy-cookbook.readthedocs.io/items/RankNullspace.html), I’ll see if you can use LU to find the nullspace instead of SVD. This answer: https://stackoverflow.com/a/2219160/500207 endorses using QR. Is your problem too big for QR? — Ahmed Fasih, Sep 26 '17 at 10:20
@AhmedFasih Actually I don't need multiple solutions, just need single solution to Ax=b when A is singular (avoiding svd/qr as they are heavy for my application). I have tried LU routines of scipy but they fails sometime give runtime warning added in the question. Thanks. — user1131274, Sep 26 '17 at 11:53
you just need to use [scipy.linalg.lstsq](https://docs.scipy.org/doc/scipy-0.19.1/reference/generated/scipy.linalg.lstsq.html#scipy.linalg.lstsq) — percusse, Sep 27 '17 at 02:36
@percusse they don’t want use the SVD (which `lstsq` uses) because SVD slow. — Ahmed Fasih, Sep 27 '17 at 03:17
@AhmedFasih You can select the lapack driver `gelsy` which uses QR factorization. — percusse, Sep 27 '17 at 09:03
@percusse oh cool, good to know! But they don’t want to use QR either. — Ahmed Fasih, Sep 27 '17 at 10:36
@AhmedFasih I didn't see where OP mentions it but QR is quite cheaper than SVD. — percusse, Sep 27 '17 at 10:54
@percusse Right, it was in a [comment](https://stackoverflow.com/questions/46421252/solving-square-singular-system?noredirect=1#comment79808482_46421252), “avoiding svd/qr as they are heavy for my application”. It’s reasonable to wonder how to do this with a dirt-cheap LU factorization… I’m trying to figure it out too. — Ahmed Fasih, Sep 27 '17 at 12:56
@AhmedFasih I mean its hard to believe why LU fails for singular matrices? Is it because matrix is singular so its hard to distinguish between very small value and 0. In theory shouldn't by setting free variables to something(let say 0) we should get some solution. Am I missing something? — user1131274, Sep 27 '17 at 13:06
It's because then you have `PLU x = b` and after having inverses you have `U x = b2` with at the bottom right an upper triangular part zeros on the diagonal. So you basically accomplished nothing yet about free variables. You can then truncate and invert the nonsingular part but still you need to solve the new smaller least squares problem. — percusse, Sep 27 '17 at 14:35
@percusse But shouldn't the backward(or forward) substitution work because U is upper triangular. For example, after setting free variables to zero we start from the bottom of U finding one variable each time(textbook procedure). Isn't this how it works? If I am not wrong inverses are also found out in similar way. Thanks — user1131274, Sep 27 '17 at 14:52
Try to solve `[0 1 2;0 0 3; 0 0 0]` with arbitrary RHS. How would you proceed? It is still possible but not with the tools what you are trying to use. — percusse, Sep 27 '17 at 15:05
@precusse Since b is in range space of A in my question lets take b=[1,1,0] for example. Then we have 3x_3 = 1 giving x_3=1/3 and x_2+2x_3=1 giving x_2=1/3. x_1 can be set arbitrarily. Am I missing something? Are u saying LU factorization can not be done for singular matrices? Sorry I am little confused. — user1131274, Sep 27 '17 at 16:28

score 1 · Answer 1 · answered Sep 26 '17 at 08:36

Try scipy.linalg.lu(). It Computes pivoted LU decompostion of a matrix.

    import pprint
    import scipy
    import scipy.linalg   # SciPy Linear Algebra Library

    A = scipy.array([ [7, 3, -1, 2], [3, 8, 1, -4], [-1, 1, 4, -1], [2, -4, -1, 6] ])
    P, L, U = scipy.linalg.lu(A)

    print "A:"
    pprint.pprint(A)

    print "P:"
    pprint.pprint(P)

    print "L:"
    pprint.pprint(L)

    print "U:"
    pprint.pprint(U)

The output from the code is given below:

   A:
   array([[ 7,  3, -1,  2],
          [ 3,  8,  1, -4],
          [-1,  1,  4, -1],
          [ 2, -4, -1,  6]])
   P:
   array([[ 1.,  0.,  0.,  0.],
          [ 0.,  1.,  0.,  0.],
          [ 0.,  0.,  1.,  0.],
          [ 0.,  0.,  0.,  1.]])
   L:
   array([[ 1.        ,  0.        ,  0.        ,  0.        ],
          [ 0.42857143,  1.        ,  0.        ,  0.        ],
          [-0.14285714,  0.21276596,  1.        ,  0.        ],
          [ 0.28571429, -0.72340426,  0.08982036,  1.        ]])
   U:
   array([[ 7.        ,  3.        , -1.        ,  2.        ],
          [ 0.        ,  6.71428571,  1.42857143, -4.85714286],
          [ 0.        ,  0.        ,  3.55319149,  0.31914894],
          [ 0.        ,  0.        ,  0.        ,  1.88622754]])

Um, could you please elaborate a bit on this? I assume op knows how to use LU decomposition in Scipy. How can LU be used to obtain the nullspace, and thus multiple solutions? — Ahmed Fasih, Sep 26 '17 at 10:18
@N.K my question is regarding singular matrices. I know of LU routine but it sometimes give runtime warning for singular matrices. I have added some explanation in the question itself. — user1131274, Sep 26 '17 at 11:44

sascha · Answer 2 · 2017-09-26T12:11:45.657

You might try out all these iterative-solvers for your task. As i'm not an expert in regards to these singular-systems, i just mention some links and show a small example:

Code using lsqr (tuned for sparse applications though):

import numpy as np
from scipy.sparse.linalg import lsqr

A = np.array([[0,0],[0,1]])
b = [1, 2]
x = lsqr(A, b)[0]
print(x)
# [ 0.  2.]

A whole dissertation about the topic:
- ITERATIVE METHODS FOR SINGULAR LINEAR EQUATIONS AND LEAST-SQUARES PROBLEMS (Sou-Cheng (Terrya) Choi )
SOLVING ILL-CONDITIONED AND SINGULAR LINEAR SYSTEMS: A TUTORIAL ON REGULARIZATION
Some discussion at math.stackexchange, recommending to solve: AtAx = Atb

Usually i would recommend linalg.lstsq, but it's probably one of those SVD-based solutions you don't want, although this one is highly optimized for the dense case.

Thanks a lot. But my system is ill-conditioned and convergence of iterative procedure depends upon condition numbers. Adding regularization would help but will change the problem slightly. I just can't understand why simple LU dosen't work. I mean in theory they teach us that when system are singular you can set free variables to zero and find a solution but all these procedures fail on singular matrices. Am I missing something. Thanks again. — user1131274, Sep 26 '17 at 12:41

Solving square singular system?

2 Answers2