Question about LU decomposition for non-square and non-invertible matrices

[Tom H]

I came across the following issue when trying to use Sympy to compute an LU decomposition of a matrix. I'd like to determine the number of solutions a system of equations has, for example, if

A = sympy.Matrix([[1,2,3],[4,5,6],[7,8,9]])

b = sympy.Matrix([4, 13, 22])

then the system A*x=b has infinitely many solutions, all of which can be written as

x0 + t*n

where

x0 = sympy.Matrix([2,1,0])

n = sympy.Matrix([1,-2,1])

t = sympy.Symbol('t')

This solution can be computed from the LU decomposition of A, however the following code fails to compute the factors L and U.

>>> import sympy

>>> A = sympy.Matrix([[1,2,3],[4,5,6],[7,8,9]])

>>> LU, perm = A.LUdecomposition()

Traceback (most recent call last):

  File "<stdin>", line 1, in <module>

  File "/Library/Python/2.7/site-packages/sympy/matrices/matrices.py", line 1297, in LUdecomposition

    combined, p = self.LUdecomposition_Simple(iszerofunc=_iszero)

  File "/Library/Python/2.7/site-packages/sympy/matrices/matrices.py", line 1341, in LUdecomposition_Simple

    raise ValueError("No nonzero pivot found; inversion failed.")

ValueError: No nonzero pivot found; inversion failed.

A is square and noninvertible, which explains the message in the stack trace about the elimination algorithm encountering a zero subcolumn below the pivot position. Regardless, A can still be written as the product of lower triangular L with unit diagonal, and upper triangular U, where

L = sympy.Matrix([[1, 0, 0], [4, 1, 0], [7, 2, 1]])

U = sympy.Matrix([[1,  2,  3], [0, -3, -6], [0,  0,  0]])

Is A.LUdecomposition() the wrong function to call if I want to compute the LU decomposition of A in the case where A is noninvertible, or not square? I wrote my own LU decomposition routine to handle these cases, which I'd be happy to contribute, but I'd like to know if this functionality exists in Sympy.

Thanks in advance for your input,

Tom

[AMiT Kumar]

Hi Tom,

LUdecomposition() doesn't works for under determined systems, you can use , linsolve(), It supports all types of systems:

In [1]: from sympy import *

In [2]: from sympy.solvers.solveset import linsolve

In [3]: A = Matrix([[1,2,3],[4,5,6],[7,8,9]])

In [4]: b = Matrix([4, 13, 22])

In [5]: system = (A, b)

In [6]: r, s, t = symbols('r s t')

In [7]: linsolve((A, b), (r, s, t))
Out[8]: {(t + 2, -2⋅t + 1, t)}

If you are interested in Matrix form results, you can use gauss_jordan_solve() 

In [3]: A = Matrix([[1,2,3],[4,5,6],[7,8,9]])


In [4]: b = Matrix([4, 13, 22])


In [5]: sol, params = A.gauss_jordan_solve(b)

In [7]: init_printing()

In [8]: sol

Out[8]: 

⎡ τ₀ + 2  ⎤

⎢         ⎥

⎢-2⋅τ₀ + 1⎥

⎢         ⎥

⎣   τ₀    ⎦

Please note that, both these functions are not available in any release, you can use them from the development version: https://github.com/sympy/sympy

--

AMiT Kumar

[Tom H]

Hi Amit,

  Thanks for the speedy response! If I understand the code in gauss_jordan_solve correctly, solving A*x=b with multiple right hand sides requires a call to gauss_jordan_solve for each right hand side b, which means that the same Gaussian elimination steps are performed on all but the last column of the augmented matrix [A | b] for each b. Is this correct? I have attached what is probably the millionth implementation of Gaussian elimination for computing the LU decomposition of a matrix, in case someone finds it useful. Let me know if this is not a proper venue for sharing such code. Thanks again for your input.

Regards,

Tom

[AMiT Kumar]

Hi Tom,

Sorry for late reply, if you are interested in contributing

code, then you can send us a Pull Request on Github:

https://github.com/sympy/sympy/pulls/

BTW, Wolframalpha also doesn't calculates the LUDecomposition

of the Matrix you mentioned, since it's singular.

Though, your code works for singular matrices, but

I am not sure, If the LUDecomposition method in SymPy,

should compute LU Decomposition of a singular matrix.

Anyways, you can send us a pull Request someone would

surely help.

AMiT Kumar

[Tom H]

Hi Amit,

  I submitted a pull request. Thanks for the guidance!
Best,

Tom