Transposes and LDU Factorization
Show that
\[A^T \cdot A \ \neq \ A \cdot A^T\]
The dot product of 2 matrices is an operation on the row of operand 1 and the column of operand 2
It will be shown that by flipping the operands the values within the respective row and column can change
This function decomposes a dot product position row 1, column 1, shows the values that effect position 1,1,
This will support the posit that it is not a commutative operation.
decompose_dot_product<-function(A,B) {
print("The decomposition of position (1,1) is : ")
print(paste("(", A[1,1], " * ", B[1,1], ") + (", A[1,2], "* ", B[2,1], ") + (", A[1,3], " * ", B[3,1], ")"))
print("The dot product matrix value in position (1,1) is ")
print((A %*% B)[1,1])
# explicitly
# print((A[1,1] * B[1,1]) + (A[1,2] * B[2,1]) + (A[1,3] * B[3,1]))
}Here we create matrix A
\[A \ = \ \begin{bmatrix} 1&2&3\\ 4&5&6\\ 7&8&9 \end{bmatrix}\]
a<-matrix(seq(1,9), byrow = TRUE, ncol = 3) Here we create the transpose of T(A).
\[T(A) \ = \ \begin{bmatrix} 1&4&7\\ 2&5&8\\ 3&6&9 \end{bmatrix}\]
t_a<-t(a) For any 2 Matrices, the dot product result will have a row/column corresponding the row of operand 1 and the column of operand 2.
So in this case, row 1 of A (1,2,3) times column 1 of T(A) (1,2,3) = 14
decompose_dot_product(a,t_a)## [1] "The decomposition of position (1,1) is : "
## [1] "( 1 * 1 ) + ( 2 * 2 ) + ( 3 * 3 )"
## [1] "The dot product matrix value in position (1,1) is "
## [1] 14When the operands are exchanged, the operands are sourced from different values.
Now row 1 of T(A) (1,4,7) times column 1 of A (1,4,7) = 14 = 66
decompose_dot_product(t_a,a)## [1] "The decomposition of position (1,1) is : "
## [1] "( 1 * 1 ) + ( 4 * 4 ) + ( 7 * 7 )"
## [1] "The dot product matrix value in position (1,1) is "
## [1] 66We previously saw commutativity of A and T(A) is only assured if the columns of operand 2 equal the rows of operand 1
By definition, a symmetric matrix is identical to its transpose and thus in that case, the dot product is commutative
Here we create the following matrix.
\[\begin{bmatrix} 4&6&2\\ 6&0&5\\ 2&5&0 \end{bmatrix}\]
a<-matrix(c(4,6,2,6,0,5,2,5,0), byrow = TRUE, ncol = 3)
t_a<-t(a) Its clear the rows match the columns. Now pass A and its transpose to the decomposition function. Do it again, flipping the operands.
decompose_dot_product(a, t_a)## [1] "The decomposition of position (1,1) is : "
## [1] "( 4 * 4 ) + ( 6 * 6 ) + ( 2 * 2 )"
## [1] "The dot product matrix value in position (1,1) is "
## [1] 56decompose_dot_product(t_a,a)## [1] "The decomposition of position (1,1) is : "
## [1] "( 4 * 4 ) + ( 6 * 6 ) + ( 2 * 2 )"
## [1] "The dot product matrix value in position (1,1) is "
## [1] 56An LU Matrix refers to the factorization of a matrix into an upper triangular and lower triangular matrix.
\[A \ = \ L \ \cdot \ U\]
All numbers below the diagonal are zero
All numbers below the diagonal are zero
There are several algorithms to perform LU or LDU ( seperating the diagonal component ) factorization.
And there are several functions in R.
This homework translated the VBA submission to the rosettacode challenge to code different algorithms in different languages.
The original code is found here
My translated function is below.
# returns a permutation function to flip rows
pivotize<-function(A) {
l<-length(A)
n = dim(A)[1]
# create an identity matrix
im<-diag(n)
for (i in (1:n)) {
mx<-abs(A[i,i]) # the diagonal value
row_<-i
for (j in (i:n)) { # loop from diag point to the end
if (abs(A[j,i]) > mx) {
mx<-abs(A[j,i]) # save off the new max
row_<-j
}
} # end of j loop 1
for (j in 1:n) {
tmp<-im[i,j]
im[i,j]<-im[row_,j]
im[row_,j]<-tmp
} # end of j loop 2
} # end of i loop
return (im)
}
lu<-function(A) {
l<-length(A)
n = dim(A)[1]
L<-matrix(rep(0,l), byrow = TRUE, ncol = n)
U<-matrix(rep(0,l), byrow = TRUE, ncol = n)
# this returns a permutation matrix that flips rows
P<-pivotize(A)
A2<-P %*% A
for ( j in 1:n) {
L[j, j] = 1
for ( i in 1:j) {
sum1<-0
for (k in 1:i) {
sum1<-sum1 + U[k,j] * L[i,k]
}
U[i,j]<-A2[i,j] - sum1
} # end of i loop
# we only want to iterate betwen j+1 up to n-1
# VBA is different than R, so this could be cleaned up
jplus<-j+1
if (jplus > n) {
jplus<-n
}
for ( i in jplus:n) {
sum2<-0
for (k in 1:j) {
sum2<-sum2 + U[k,j] * L[i,k]
}
if(i != j) {
L[i,j]<- (A2[i,j]-sum2) / U[j,j]
}
} # end of second i loop
} # end of j loop
cat("Lower ")
prmatrix(round(L,2), rowlab=rep("",3), collab=rep("",3))
cat("\n\nUpper \n")
prmatrix(round(U,2), rowlab=rep("",3), collab=rep("",3))
cat("\n\nPermutation \n")
prmatrix(P, rowlab=rep("",3), collab=rep("",3))
cat("\n\nPermutation %*% Lower %*% Upper \n")
prmatrix(P %*% L %*% U, rowlab=rep("",3), collab=rep("",3))
}Invoke the function for the following :
\[\begin{bmatrix} 1&2&3\\ 7&8&9\\ 4&5&6 \end{bmatrix}\]
A<-matrix(c(1,2,3, 7,8,9, 4,5,6), byrow = TRUE, ncol = 3)
B<-lu(A)## Lower
## 1.00 0.0 0
## 0.14 1.0 0
## 0.57 0.5 1
##
##
## Upper
##
## 7 8.00 9.00
## 0 0.86 1.71
## 0 0.00 0.00
##
##
## Permutation
##
## 0 1 0
## 1 0 0
## 0 0 1
##
##
## Permutation %*% Lower %*% Upper
##
## 1 2 3
## 7 8 9
## 4 5 6The Matrix and matrixcalc packages offer lu.decomposition(), qr(), and chol()
Another way is to manually employ Gausian Elimination to “zero out the elements” below the diagonal, creating an Upper Triangular Matrix.
Those row operations, which zero out the elements, combine into a single matrix that is now the inverse of the Lower Triangular Matrix, i.e.
\[L^{-1} \ \cdot \ A \ = \ U\] where your row operations represent \(L^{-1}\)
An excellent tutorial on the manual method can be found here