Inconsistent results between LU decomposition in R

Question

I have the following matrix A in R:

           # [,1]       [,2]   [,3]       [,4]
# [1,] -1.1527778  0.4444444  0.375  0.3333333
# [2,]  0.5555556 -1.4888889  0.600  0.3333333
# [3,]  0.6250000  0.4000000 -1.825  0.8000000
# [4,]  0.6666667  0.6666667  0.200 -1.5333333

A <- structure(c(-1.15277777777778, 0.555555555555556, 0.625, 0.666666666666667, 
0.444444444444444, -1.48888888888889, 0.4, 0.666666666666667, 
0.375, 0.6, -1.825, 0.2, 0.333333333333333, 0.333333333333333, 
0.8, -1.53333333333333), .Dim = c(4L, 4L), .Dimnames = list(NULL, 
    NULL))

I compute its LU-decomposition as follows:

library(Matrix)
ex <- expand(lu(t(A)))
L <- ex$L
P <- ex$P
C <- U <- ex$U
C[lower.tri(U)] <- L[lower.tri(L)]

print(C)

# 4 x 4 Matrix of class "dgeMatrix"
           # [,1]       [,2]       [,3]          [,4]
# [1,] -1.1527778  0.5555556  0.6250000  6.666667e-01
# [2,] -0.3855422 -1.2746988  0.6409639  9.236948e-01
# [3,] -0.3253012 -0.6124764 -1.2291115  9.826087e-01
# [4,] -0.2891566 -0.3875236 -1.0000000 -2.220446e-16

On the other hand, this is the Python code for the same job:

lu, piv = scipy.linalg.lu_factor(A.T, check_finite=False)

print(lu)

# [[ -1.15277778e+00   5.55555556e-01   6.25000000e-01   6.66666667e-01]
 # [ -3.85542169e-01  -1.27469880e+00   6.40963855e-01   9.23694779e-01]
 # [ -2.89156627e-01  -3.87523629e-01   1.22911153e+00  -9.82608696e-01]
 # [ -3.25301205e-01  -6.12476371e-01  -1.00000000e+00   7.69432827e-16]]

I'm wondering why the two C and lu matrices in R and Python (respectively) are not the same. The point is that I have to get the same result as Python version (i.e. matrix lu). Do you have any idea what I am doing wrong?

Answer 1

It is embarrassing to realize 1.5 years later that my original answer is not fully correct. While it correctly pointed out that the rank-deficiency of A in the question is the cause, it is incorrect to attribute this as the root cause. The real issue is the non-unique choice of pivot, which could happen (although less likely) even if A is full-rank. Given that this post has been viewed 700+ times and has a score of 6, I might have misled many readers. SORRY!

I posted Write a trackable R function that mimics LAPACK's dgetrf for LU factorization and answered it just now. The question contains a LU function for LU factorization without pivoting, and the answer contains two functions LUP and LUP2 for a version with row pivoting that are consistent with LAPACK's dgetrf, which underlies the dense method of Matrix::lu and R base function solve. In particular, the LUP2 function allows one to track the factorization step by step. I will use this function here for investigation.

---

So you are factorization the transpose of A.

From the output of R and Python we see that they yield identical 1st pivot -1.1527778 and 2nd pivot -1.2746988, while the 3rd pivot differs. So there is definitely something interesting happening when the factorization has done the first two columns / rows and proceeds to the the 3rd column / row. Let's pause R's LU factorization at this point:

oo <- LUP2(t(A), to = 2)
#           [,1]       [,2]       [,3]       [,4]
#[1,] -1.1527778  0.5555556  0.6250000  0.6666667
#[2,] -0.3855422 -1.2746988  0.6409639  0.9236948
#[3,] -0.3253012 -0.6124764 -1.2291115  0.9826087
#[4,] -0.2891566 -0.3875236  1.2291115 -0.9826087
#attr(,"to")
#[1] 2
#attr(,"pivot")
#[1] 1 2 3 4

At this point, the Gaussian elimination has reduced t(A) to

getU(oo)
#           [,1]       [,2]       [,3]       [,4]
#[1,] -1.1527778  0.5555556  0.6250000  0.6666667
#[2,]  0.0000000 -1.2746988  0.6409639  0.9236948
#[3,]  0.0000000  0.0000000 -1.2291115  0.9826087
#[4,]  0.0000000  0.0000000  1.2291115 -0.9826087
#attr(,"to")
#[1] 2

Wow, we see something really interesting now: the 3rd and 4th rows only differ by a sign change. Then the Gaussian elimination is not unique, because either -1.2291115 or 1.2291115 can be a pivot as they have the same absolute value.

Clearly, R has chosen -1.2291115 as the pivot, but Python has chosen 1.2291115 as the pivot. A row exchange between 3rd and 4th rows will happen in Python. In your question, you didn't report what permutation index Python gives you, but it should 1, 2, 4, 3, instead of the 1, 2, 3, 4 in R.

---

Either way, the U factor ends up with a row of zeros in the bottom, so that t(A) or A is not full-rank. If you want to see a more consistent behavior between two software, you'd better try a full-rank matrix. In that case, it is less likely that you will have more than one choices for a pivot during LU factorization. You can generate a random full-rank matrix in R by

A <- matrix(runif(16), 4, 4)