The purpose of this project is to create a repository of functions for commonly used algorithms in Linear Algebra. The following functions are given:
LU factorization.
LDM^T factorization.
Cholesky factorization by LDM^T factorization.
Cholesky factorization by outer product method.
QR factorization by Householder reflections.
QR factorization by Givens rotations.
Inverse of Upper traingular Matrix.
Least Square solution to Linear System of equations by QR factorization.
We use the algorithms mentioned in the following reference:
Reference : Golub, G.H., Van Loan, C.F. Matrix Computations
For a given matrix A, find lower triangular matrix L and upper triangular U such that A = LU
LUfactorization <- function(mat){
num_rows = nrow(mat)
num_cols = ncol(mat)
U = mat
L = diag(num_rows)
num_iter = min(num_rows, num_cols) - 1
sapply(1:num_iter, function(curr_col){
M <- diag(num_rows)
Minv <- diag(num_rows)
curr_pivot_col <- U[((curr_col+1):num_rows),
curr_col]/U[curr_col, curr_col]
M[((curr_col+1):num_rows),curr_col] <- curr_pivot_col * (-1)
Minv[((curr_col+1):num_rows),curr_col] <- curr_pivot_col
U <<- M %*% U
L <<- L %*% Minv
})
return(list("L" = L, "U" = U))
}
A <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 10), ncol = 3)
lu_A <- LUfactorization(A)
A## [,1] [,2] [,3]
## [1,] 1 4 7
## [2,] 2 5 8
## [3,] 3 6 10lu_A## $L
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 2 1 0
## [3,] 3 2 1
##
## $U
## [,1] [,2] [,3]
## [1,] 1 4 7
## [2,] 0 -3 -6
## [3,] 0 0 1norm(A - lu_A$L %*% lu_A$U, type = "F")## [1] 0For a given matrix A, find lower triangular matrices L and M, and diagonal matrix D such that A = LDM^T
LDMfactorization <- function(mat){
num_rows <- nrow(mat)
num_cols <- ncol(mat)
LUfactors <- LUfactorization(mat)
U <- LUfactors$U
D <- matrix(data = 0, num_rows, num_rows)
Mtran <- U
sapply(1:num_rows, function(curr_row){
D[curr_row, curr_row] <<- Mtran[curr_row, min(curr_row, num_cols)]
Mtran[curr_row,] <<- Mtran[curr_row,]/D[curr_row, curr_row]
})
M <- t(Mtran)
return(list("L" = LUfactors$L, "D" = D, "M" = M))
}
A <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 10), ncol = 3)
ldm_A <- LDMfactorization(A)
A## [,1] [,2] [,3]
## [1,] 1 4 7
## [2,] 2 5 8
## [3,] 3 6 10ldm_A## $L
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 2 1 0
## [3,] 3 2 1
##
## $D
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 -3 0
## [3,] 0 0 1
##
## $M
## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 4 1 0
## [3,] 7 2 1norm(A - ldm_A$L %*% ldm_A$D %*% t(ldm_A$M), type = "F")## [1] 0For a given symmetric positive definite matrix A, find lower triangular matrix L such that A = LL^T using LDM^T factorization. If A is symmetric positve definite, then chol(A) = L %*% sqrt(D), where L and D is given by LDM^T factorization
CholeskyFactorizationByLDM <- function(mat){
num_rows <- nrow(mat)
LDLTfactors <- LDMfactorization(mat)
L <- LDLTfactors$L
D <- LDLTfactors$D
sqrtD <- matrix(data = 0, num_rows, num_rows)
sapply(1:num_rows, function(curr_row){
sqrtD[curr_row, curr_row] <<- sqrt(D[curr_row, curr_row])
})
L <- L %*% sqrtD
return(L)
}
Asym <- matrix(c(2, -1, 0, -1, 2, -1, 0, -1, 2), ncol = 3)
cholAsym <- CholeskyFactorizationByLDM(Asym)
Asym## [,1] [,2] [,3]
## [1,] 2 -1 0
## [2,] -1 2 -1
## [3,] 0 -1 2cholAsym## [,1] [,2] [,3]
## [1,] 1.4142136 0.0000000 0.000000
## [2,] -0.7071068 1.2247449 0.000000
## [3,] 0.0000000 -0.8164966 1.154701norm(Asym - cholAsym %*% t(cholAsym), type = "F")## [1] 6.28037e-16Cholesky lower triangular factor can be found by an alternative method which uses outer product method. More details about this method can be found in the Reference book.
CholeskyByOuter <- function(mat){
num_rows <- nrow(mat)
L <- mat
sapply(1:num_rows, function(curr_row){
alpha <- L[curr_row, curr_row]
beta <- sqrt(alpha)
if(curr_row < num_rows){
B <- L[(curr_row+1):num_rows, (curr_row+1):num_rows]
v <- L[(curr_row+1):num_rows,curr_row]
vbeta <- v/beta
v_outer <- (v %*% t(v)) / alpha
G <- B - v_outer
L[curr_row, curr_row] <<- beta
L[(curr_row+1):num_rows, curr_row] <<- vbeta
L[curr_row, (curr_row+1):num_rows] <<- 0
L[(curr_row+1):num_rows, (curr_row+1):num_rows] <<- G
}else{
L[curr_row, curr_row] <<- beta
}
})
return(L)
}
Asym <- matrix(c(2, -1, 0, -1, 2, -1, 0, -1, 2), ncol = 3)
cholOutAsym <- CholeskyByOuter(Asym)
Asym## [,1] [,2] [,3]
## [1,] 2 -1 0
## [2,] -1 2 -1
## [3,] 0 -1 2cholOutAsym## [,1] [,2] [,3]
## [1,] 1.4142136 0.0000000 0.000000
## [2,] -0.7071068 1.2247449 0.000000
## [3,] 0.0000000 -0.8164966 1.154701norm(Asym - cholOutAsym %*% t(cholOutAsym), type = "F")## [1] 7.691851e-16For a given matrix A, find orthogonal matrix Q and upper triangular matrix R such that A = QR using Householder reflections.
QRfactorizationHouseholder <- function(mat) {
num_cols <- ncol(mat)
num_rows <- nrow(mat)
R <- mat
num_iter <- min(num_cols, num_rows)
H <- lapply(1:num_iter, function(i){
x <- R[i:num_rows, i]
len_x <- length(x)
x_ext <- matrix(c(rep(0, (i-1)),
x), ncol = 1)
e <- matrix(c(rep(0, (i-1)),
1,
rep(0, (len_x - 1))), ncol = 1)
u <- sqrt(sum(x^2))
v <- x_ext + sign(x[1]) * u * e
v_inner <- (t(v) %*% v)[1]
Hi <- diag(num_rows) - 2 * (v %*% t(v)) / v_inner
R <<- Hi %*% R
return(Hi)
})
Q <- Reduce("%*%", H)
res <- list('Q'=Q,'R'=R)
return(res)
}
A <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 10), ncol = 3)
qrhh_A <- QRfactorizationHouseholder(A)
A## [,1] [,2] [,3]
## [1,] 1 4 7
## [2,] 2 5 8
## [3,] 3 6 10qrhh_A## $Q
## [,1] [,2] [,3]
## [1,] -0.2672612 0.8728716 -0.4082483
## [2,] -0.5345225 0.2182179 0.8164966
## [3,] -0.8017837 -0.4364358 -0.4082483
##
## $R
## [,1] [,2] [,3]
## [1,] -3.741657e+00 -8.552360 -14.1648458
## [2,] 7.251946e-16 1.963961 3.4914862
## [3,] -3.391787e-16 0.000000 -0.4082483norm(A - qrhh_A$Q %*% qrhh_A$R, type = "F")## [1] 7.809533e-15QR factorization can be found by an alternative method which uses Givens Rotations. More details about this method can be found in the Reference book.
QRfactorizationGivens <- function(mat) {
givens <- function(a, b){
c <- 0
s <- 0
if (b == 0){
c <- 1
s <- 0
}else{
if(abs(b) > abs(a)){
tau <- -1 * a/b
s <- 1 / sqrt(1 + tau^2)
c <- s * tau
}else{
tau <- -1 * b/a
c <- 1 / sqrt(1 + tau^2)
s <- c * tau
}
}
return(list("c" = c, "s" = s))
}
num_cols <- ncol(mat)
num_rows <- nrow(mat)
R <- mat
num_iter <- min(num_cols, num_rows) - 1
G <- diag(num_rows)
sapply(1:num_iter, function(j){
sapply(num_rows:(j+1), function(i){
givens_curr <- givens(R[(i-1), j], R[i, j])
rotation_mat <- matrix(c(givens_curr$c, -givens_curr$s,
givens_curr$s, givens_curr$c), ncol = 2)
R[(i-1):i, j:num_cols] <<- t(rotation_mat) %*% R[(i-1):i, j:num_cols]
G_curr <- diag(num_rows)
G_curr[(i-1):i, (i-1):i] <- rotation_mat
G <<- G %*% G_curr
})
})
Q <- G
res <- list('Q'=Q,'R'=R)
return(res)
}
A <- matrix(c(1, 2, 3, 4, 5, 6, 7, 8, 10), ncol = 3)
qrgiv_A <- QRfactorizationGivens(A)
A## [,1] [,2] [,3]
## [1,] 1 4 7
## [2,] 2 5 8
## [3,] 3 6 10qrgiv_A## $Q
## [,1] [,2] [,3]
## [1,] 0.2672612 0.8728716 -0.4082483
## [2,] 0.5345225 0.2182179 0.8164966
## [3,] 0.8017837 -0.4364358 -0.4082483
##
## $R
## [,1] [,2] [,3]
## [1,] 3.741657e+00 8.552360e+00 14.1648458
## [2,] 1.110223e-16 1.963961e+00 3.4914862
## [3,] 0.000000e+00 -1.110223e-16 -0.4082483norm(A - qrgiv_A$Q %*% qrgiv_A$R, type = "F")## [1] 4.167405e-15Given an upper triangular matrix U, find Uinv such that Uinv %*% U = I
InverseUpperTriangular <- function(U){
if(any(diag(U) == 0)){
stop("Upper triangular matrix not invertible, 0 found in diagonal")
}
num_rows <- nrow(U)
UI <- cbind(U, diag(num_rows))
sapply(num_rows:1, function(curr_row){
curr_pivot <- UI[curr_row, curr_row]
UI[curr_row, ] <<- UI[curr_row, ]/curr_pivot
if(curr_row > 1){
sapply((curr_row-1):1, function(i){
curr_row_mult <- UI[i, curr_row]
UI[i, ] <<- UI[i, ] - UI[curr_row, ] * curr_row_mult
})
}
})
Uinv <- UI[, (num_rows+1):(2*num_rows)]
return(Uinv)
}
U <- matrix(c(1, 0, 0, 4, 5, 0, 7, 8, 10), ncol = 3)
Uinv <- InverseUpperTriangular(U)
U## [,1] [,2] [,3]
## [1,] 1 4 7
## [2,] 0 5 8
## [3,] 0 0 10Uinv## [,1] [,2] [,3]
## [1,] 1 -0.8 -0.06
## [2,] 0 0.2 -0.16
## [3,] 0 0.0 0.10Uinv %*% U## [,1] [,2] [,3]
## [1,] 1 0 -8.881784e-16
## [2,] 0 1 0.000000e+00
## [3,] 0 0 1.000000e+00For a given matrix A and vector b, find least square solution x for Ax = b using QR decomposition. For an overdetermined system, the solution minimizes L2 norm ||Ax - b||.
LeastSquareSolutionByQR <- function(A, b){
y <- matrix(data = b, ncol=1)
num_rows <- nrow(A)
num_cols <- ncol(A)
if(num_cols > num_rows){
stop("Underdetermined system")
}
if(nrow(y) != num_rows){
stop("Incompatible dimensions of A and b")
}
QRfactor <- QRfactorizationHouseholder(A)
Q <- QRfactor$Q
R1 <- QRfactor$R[1:num_cols,]
Rinv <- InverseUpperTriangular(R1)
Qty <- t(Q) %*% y
c <- Qty[1:num_cols]
x_LS <- Rinv %*% c
return(x_LS)
}
Atri <- matrix(c(1,3,5,2,4,6), ncol=2)
b <- matrix(c(3,2,5), ncol = 1)
x_LS <- LeastSquareSolutionByQR(Atri, b)
Atri## [,1] [,2]
## [1,] 1 2
## [2,] 3 4
## [3,] 5 6b## [,1]
## [1,] 3
## [2,] 2
## [3,] 5x_LS## [,1]
## [1,] -1.333333
## [2,] 1.833333norm(b - Atri %*% x_LS, type = "2")## [1] 1.632993