Crout LU decomposition in R

Solution 1:

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?

Solution 2:

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