I have singular matrix A(10*10) which is rank deficient(rank=9) and I have vector b which is in range space of A. Now I am interested in some solution to Ax=b. For concreteness here is my A
array([[ 0. , 0. , 0. , 0.86826141, 0. ,
0. , 0.88788426, 0. , 0.4089203 , 0.88134901],
[ 0. , 0. , 0.46416372, 0. , 0. ,
0. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0.31303966,
0. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0.3155742 , 0. , 0.64059294, 0. ],
[ 0. , 0. , 0. , 0. , 0.51349938,
0. , 0. , 0. , 0.53593509, 0. ],
[ 0. , 0.01252787, 0. , 0.6870415 , 0. ,
0. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0. ],
[ 0. , 0. , 0. , 0.16643105, 0. ,
0. , 0. , 0. , 0. , 0. ],
[ 0.08626592, 0. , 0. , 0. , 0. ,
0. , 0. , 0. , 0. , 0.66939531],
[ 0.43694586, 0. , 0. , 0. , 0. ,
0.95941661, 0. , 0.52936733, 0.79687149, 0.81463887]])
b is generated using A.dot(np.ones(10)). Now I wanted to solve this using lu factorization and for that I did this following
lu_fac=scipy.linalg.lu_factor(X)
scipy.linalg.lu_solve(lu_fac,b)
Which gives
array([ nan, nan, nan, nan, nan, nan, nan, nan, nan, nan])
Also lu_factor seems to be working fine in this case(some time it does give run time warning saying "Diagonal number %d is exactly zero. Singular matrix"). For completeness here is the code for verifying PLU from lu_factor is same as A :
L=np.tril(lu_fac[0])
np.fill_diagonal(L,1)
U=np.triu(lu_fac[0])
perm=np.arange(10)
ipiv=lu_factor[1]
for i in range(10):
temp=perm[i]
perm[i]=perm[ipiv[i]]
perm[ipiv[i]]=temp
np.allclose(X[perm,:],L.dot(U))
Now I know my matrix is singular and there are infinitely many solutions to my problem. But I am interested in any solution and I am just confused why lu factorization fails, can't it set free variables to 0 and find some solution as we are taught? Also what is the deal with the run time warning "Diagonal number %d is exactly zero. Singular matrix". Note I am not interested in svd/qr approach to solve this, I am just curious to know why lu fails for singular matrices. Any suggestions are greatly appreciated. Thanks.
0 / lu_fac[0][9, 9]
returns nan because that entry - the last diagonal entry of U, is zero. So this nan becomes the value of the 9-th variable. Then it's substituted in the equations above, and naturally, the rest comes out as nan, too. SciPy's LU code, or rather the Fortran code it wraps, is not designed for rank-deficient matrices, so it is not going to make up values for the variables that can't be determined.
Also what is the deal with the run time warning "Diagonal number %d is exactly zero. Singular matrix".
The warning is clear: the algorithm detected a singular matrix, which is not expected. It also tells you that the implementation is not intended for use with singular matrices.
have vector b which is in range space of A
That's theoretically. In practice, one can't be sure about anything being in the range space of a rank-deficient matrix because of the errors inherent in the floating point arithmetics. You can compute b = A.dot(...) and then try to solve Ax=b, and there won't be a solution because of the errors introduced when manipulating floating point numbers.
By the way: you mentioned that PLU factorization exists for every square matrix, but SciPy is not necessarily designed to compute it. For example,
scipy.linalg.lu_factor(np.array([[0, 1], [0, 0]]))
returns a matrix with NaNs. In your case, NaN appear later, when attempting to find a solution and encountering a zero diagonal element of factor U.