Solving square singular system?

user_1131274 0
Viewed 1604

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.

scipy linear-algebra python
2 Answers

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]])
user_81702850

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.]

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.

user_23200350
Related