library(stats)For Matrix A
A B
C D
And Matrix At
A C
B D
The product A*At is equal to
AA + BB AC + BD
CA + DB CC + DD
And At*A is equal to
AA+BB AB+CD
BA+DC CC+DD
For unique values of B and C, these are not equivalent.
A <- matrix(c(1,2,3,4),nrow = 2,byrow = 2)
At <- t(A)
AAt <- A%*%At
AtA <- At%*%A
identical(AAt,AtA)## [1] FALSEIn cases where B=C, we see that both products are equal. Using the above example we see they are both equal to:
AA + BB AB + BD
BA + DB BB + DD
A more general form of this relationship, say for larger matrices, is the observation that the matrix must be equal to its transposition (At == A). This is to say, the matrix must be symmetrical.
x <- round(runif(6,min=0,max=9))
A <- matrix(data = c(x[1],x[2],x[3],x[2],x[4],x[5],x[3],x[5],x[6]),nrow=3)
At <- t(A)
identical(A,At)## [1] TRUEAAt <- A%*%At
AtA <- At%*%A
identical(AAt,AtA)## [1] TRUELUDecomp <- function(M){
#Function provides original matrix, TRUE/FALSE value indicating successful factorization, upper and lower triangle matrix, in order.
#Create L and U Matrix (Easy way to keep row count the same)
L <- M
U <- M
#Solve First Row in each Matrix (Row 1 of U is copied from M)
L[1,] <- rep(0,ncol(M))
L[1,1] <- 1
#Solve Second Rows
L[2,] <- rep(0,ncol(M))
L[2,2]<- 1
L[2,1] <- M[2,1]/M[1,1]
U[2,] <- U[2,]-L[2,1]*U[1,]
#Below can be rewritten in general form using looping and relative reference (to solve any matrix size). For purposes of this exercise, this is cleaner as is:
#Solve Third Row (If needed)
if(nrow(M) > 2){
L[3,] <- rep(0,ncol(M))
L[3,3] <- 1
L[3,1] <- M[3,1]/M[1,1]
U[3,] <- U[3,]-L[3,1]*U[1,]
L[3,2] <- U[3,2]/U[2,2]
U[3,] <- U[3,]-L[3,2]*U[2,]
}
#Solve Fourth Row (If needed)
if(nrow(M) > 3){
L[4,] <- rep(0,ncol(M))
L[4,4] <- 1
L[4,1] <- M[4,1]/M[1,1]
U[4,] <- U[4,]-L[4,1]*U[1,]
L[4,2] <- U[4,2]/U[2,2]
U[4,] <- U[4,]-L[4,2]*U[2,]
L[4,3] <- U[4,3]/U[3,3]
U[4,] <- U[4,]-L[4,3]*U[3,]
}
#Return:
MTest <- L%*%U
print(M)
print(identical(M,MTest))
print(L)
print(U)
#print(L%*%U)
}#Test 3X3
M <- matrix(data = c(1,4,-3,-2,8,5,3,4,7),nrow = 3, byrow = TRUE)
LUDecomp(M)## [,1] [,2] [,3]
## [1,] 1 4 -3
## [2,] -2 8 5
## [3,] 3 4 7
## [1] TRUE
## [,1] [,2] [,3]
## [1,] 1 0.0 0
## [2,] -2 1.0 0
## [3,] 3 -0.5 1
## [,1] [,2] [,3]
## [1,] 1 4 -3.0
## [2,] 0 16 -1.0
## [3,] 0 0 15.5#Test 4x4
set.seed(321)
M <- matrix(data = round(runif(16,min=0,max = 9)),nrow=4)
LUDecomp(M)## [,1] [,2] [,3] [,4]
## [1,] 9 4 4 7
## [2,] 8 3 7 0
## [3,] 2 4 5 5
## [4,] 2 3 3 2
## [1] TRUE
## [,1] [,2] [,3] [,4]
## [1,] 1.0000000 0.0 0.0000000 0
## [2,] 0.8888889 1.0 0.0000000 0
## [3,] 0.2222222 -5.6 1.0000000 0
## [4,] 0.2222222 -3.8 0.6495726 1
## [,1] [,2] [,3] [,4]
## [1,] 9 4.0000000 4.000000 7.000000
## [2,] 0 -0.5555556 3.444444 -6.222222
## [3,] 0 0.0000000 23.400000 -31.400000
## [4,] 0 0.0000000 0.000000 -2.803419#Test 2x2
set.seed(321)
M <- matrix(data = round(runif(4,min=0,max = 9)),nrow=2)
LUDecomp(M)## [,1] [,2]
## [1,] 9 2
## [2,] 8 2
## [1] TRUE
## [,1] [,2]
## [1,] 1.0000000 0
## [2,] 0.8888889 1
## [,1] [,2]
## [1,] 9 2.0000000
## [2,] 0 0.2222222#First draft of 3x3 test (Just as a reference)
LUDecomp3 <- function(M){
#Decomposes 3x3 Matrix into LU
#Set Knowns
L = c(0,0)
U = c(0,0)
L[1] = 1
L[2] = 0
L[3] = 0
L[5] = 1
L[6] = 0
L[9] = 1
U[4] = 0
U[7] = 0
U[8] = 0
U[1] = M[1,1]
U[2] = M[1,2]
U[3] = M[1,3]
#Solve others
#Second Row, Col 1
L[4] <- M[2,1]/M[1,1]
U[5] <- M[2,2] - M[1,2]*L[4]
U[6] <- M[2,3] - M[1,3]*L[4]
#Third Row, Col 1
L[7] <- M[3,1]/M[1,1]
U[8] <- M[3,2] - M[1,2]*L[7]
U[9] <- M[3,3] - M[1,3]*L[7]
#Third Row, Col 2
L[8] <- U[8]/U[5]
U[8] <- 0
U[9] <- U[9] - L[8]*U[6]
LMatrix <- matrix(data = L,nrow = 3, byrow = TRUE)
UMatrix <- matrix(data = U,nrow=3,byrow = TRUE)
print(LMatrix)
print(UMatrix)
}