Is the upper triangular matrix a scipy.linalg.lu function always in series form?
I have an mxn matrix A, with n> m, and I am trying to identify independent rows using the row echelon form. Scipy.linalg.lu function returns the PLU factorization of my matrix, but the U-factor doesn't seem to be in echelon form, i.e. Anchor points are not in staircase view. As far as I know, the U-factor should always be in the form of a ladder.
Consider the following example:
from numpy import array
from scipy.linalg import lu
A = array([[1, 1, 1, 1, 0, 1, 1, 1, 1, 0],
[1, 1, 0, 0, 1, 0, 1, 0, 1, 1],
[1, 1, 0, 0, 0, 1, 1, 0, 0, 0],
[0, 1, 0, 1, 1, 0, 1, 0, 0, 1],
[1, 1, 0, 0, 1, 1, 1, 1, 1, 1]])
P, L, U = lu(A)
The U-factor is not in the lower case. For each row k, the pivot points must always be to the right of the axis of rotation in row k-1. See that in the fifth row, the pivot is not to the right of the pivot in the fourth row:
array([[ 1., 1., 1., 1., 0., 1., 1., 1., 1., 0.],
[ 0., 1., 0., 1., 1., 0., 1., 0., 0., 1.],
[ 0., 0., -1., -1., 0., 0., 0., -1., -1., 0.],
[ 0., 0., 0., 0., 1., -1., 0., 0., 1., 1.],
[ 0., 0., 0., 0., 1., 0., 0., 1., 1., 1.]])
source to share
I've never seen LUP decomposition before, but I suspect it just doesn't do what you think it does. I did the same decomposition in Matlab as you did in Python and I got the same thing as you did. So I don't think this is an issue with the LU Python function.
Update: I also implemented this in R (code below) and again U-matrix gives the same results as Matlab and Python. Unlike others, it gives the following warning:
Warning message:
In .local(x, ...) :
Exact singularity detected during LU decomposition: U[i,i]=0, i=4.
Matlab:
code:
A = ([[1, 1, 1, 1, 0, 1, 1, 1, 1, 0],
[1, 1, 0, 0, 1, 0, 1, 0, 1, 1],
[1, 1, 0, 0, 0, 1, 1, 0, 0, 0],
[0, 1, 0, 1, 1, 0, 1, 0, 0, 1],
[1, 1, 0, 0, 1, 1, 1, 1, 1, 1]]);
[L,U,P] = lu(A);
U
Matlab results:
U =
1 1 1 1 0 1 1 1 1 0
0 1 0 1 1 0 1 0 0 1
0 0 -1 -1 0 0 0 -1 -1 0
0 0 0 0 1 -1 0 0 1 1
0 0 0 0 1 0 0 1 1 1
Python results (from your code):
U =
array([[ 1., 1., 1., 1., 0., 1., 1., 1., 1., 0.],
[ 0., 1., 0., 1., 1., 0., 1., 0., 0., 1.],
[ 0., 0., -1., -1., 0., 0., 0., -1., -1., 0.],
[ 0., 0., 0., 0., 1., -1., 0., 0., 1., 1.],
[ 0., 0., 0., 0., 1., 0., 0., 1., 1., 1.]])
R:
code:
library(Matrix)
A <- Matrix( c( 1, 1, 1, 1, 0, 1, 1, 1, 1, 0 , 1, 1, 0, 0, 1, 0, 1, 0, 1, 1 , 1, 1, 0, 0, 0, 1, 1, 0, 0, 0 , 0, 1, 0, 1, 1, 0, 1, 0, 0, 1 , 1, 1, 0, 0, 1, 1, 1, 1, 1, 1 ),nrow=5,ncol=10,byrow=TRUE )
expand(mylu <-lu(A))
Results:
class "dtrMatrix" (unitriangular)
[,1] [,2] [,3] [,4] [,5]
[1,] 1 . . . .
[2,] 0 1 . . .
[3,] 1 0 1 . .
[4,] 1 0 1 1 .
[5,] 1 0 1 0 1
$U
5 x 10 Matrix of class "dgeMatrix"
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1 1 1 1 0 1 1 1 1 0
[2,] 0 1 0 1 1 0 1 0 0 1
[3,] 0 0 -1 -1 0 0 0 -1 -1 0
[4,] 0 0 0 0 1 -1 0 0 1 1
[5,] 0 0 0 0 1 0 0 1 1 1
$P
5 x 5 sparse Matrix of class "pMatrix"
[1,] | . . . .
[2,] . . . | .
[3,] . . | . .
[4,] . | . . .
[5,] . . . . |
source to share