Crout LU decomposition in R-CodePudding

Just wondering if anybody would share an implementation of the Crout algorithm for a LU decomposition (A = L * U) in R. There's a lu function in pracma library but uses Doolite instead.

> A
     [,1] [,2] [,3]
[1,]  -10   30   50
[2,]   -6   16   22
[3,]   -2    1   -5
> lu(A)
$L
     [,1] [,2] [,3]
[1,]  1.0  0.0    0
[2,]  0.6  1.0    0
[3,]  0.2  2.5    1

$U
     [,1] [,2] [,3]
[1,]  -10   30   50
[2,]    0   -2   -8
[3,]    0    0    5

Whereas for the Crout algorithm you'd get something like this instead:

$L
     [,1] [,2] [,3]
[1,]  -10    0    0
[2,]   -6   -2    0
[3,]   -2   -5    5

$U
     [,1] [,2] [,3]
[1,]    1   -3   -5
[2,]    0    1    4
[3,]    0    0    1

I've been googling for something like that for a while but didn't found any working implementation in R.

Thanks!

CodePudding user response：

It is not difficult to convert some MATLAB code to an R implementation:

LUcrout <- function(A) {
    n <- nrow(A)
    L <- matrix(0, n, n); U <- matrix(0, n, n)
    for (i in 1:n) {
        L[i, 1] <- A[i, 1]
        U[i, i] <- 1
    }
    for (j in 2:n) {
        U[1, j] <- A[1, j] / L[1, 1]
    }
    for (i in 2:n) {
        for (j in 2:i) {
            L[i, j] <- A[i, j] - L[i, 1:(j-1)] %*% U[1:(j-1), j]
        }
        if (i < n) {
            for (j in ((i 1):n)) {
                U[i, j] = (A[i, j] - L[i, 1:(i-1)] %*% U[1:(i-1), j]) / L[i, i]
            }
        }
    }
    return(list(L = L, U = U))
}

Applied to your matrix A example it returns

A = matrix(c(-10, 30, 50,
              -6, 16, 22,
              -2,  1, -5), 3, 3, byrow = TRUE)

> LUcrout(A)
$L
     [,1] [,2] [,3]
[1,]  -10    0    0
[2,]   -6   -2    0
[3,]   -2   -5    5

$U
     [,1] [,2] [,3]
[1,]    1   -3   -5
[2,]    0    1    4
[3,]    0    0    1

which is not the same as what you suggested, but is identical to what MATLAB returns (see the 'Crout-matrix-decomposition' page on Wikipedia).

LU decompositions are not unique. Which advantages does the Crout algorithm have in your opinion?

CodePudding user response：

Take the Doolittle decomposition of A transpose, swap L and U and take their transposes.