Let’s start with taking the transpose of \(A^{T}\)A:
\((A^{T}A)^{T}\)
Using the theorem of matrix multiplication, we see this equals:
\(A^{T}A\)
Which is not equal to \(AA^{T}\)
Thinking in terms of non-square matrices.. Suppose A is an mxn matrix. Then A^T is nxm matrix. If we multiply A * \(A^{T}\) then we’d be multiplying a mxn matrix by nxm matrix to return a mxm matrix. If we multiply \(A^{T}\) * A then we’d be multiplying a nxm matrix by mxn matrix to get a nxn matrix.
If the matrix A is not square meaning that m\(\neq\)n then the resulting matrices will be of different sizes.
#### Examples
Let A be a matrix of size 2x3:
A = matrix(c(2,1,4,3,5,4),nrow=2)
A## [,1] [,2] [,3]
## [1,] 2 4 5
## [2,] 1 3 4The transpose of A, \(A^{T}\), is a 3x2 matrix:
A = matrix(c(2,1,4,3,5,4),nrow=2)
t(A)## [,1] [,2]
## [1,] 2 1
## [2,] 4 3
## [3,] 5 4A * \(A^{T}\) returns 2x2 matrix:
A = matrix(c(2,1,4,3,5,4),nrow=2)
A%*%t(A)## [,1] [,2]
## [1,] 45 34
## [2,] 34 26\(A^{T}\) * A returns a 3x3 matrix:
A = matrix(c(2,1,4,3,5,4),nrow=2)
t(A) %*% A## [,1] [,2] [,3]
## [1,] 5 11 14
## [2,] 11 25 32
## [3,] 14 32 41It’s true if the matrix is symmetric because the transpose of a symmetric matrix is equal to itself. This is also true whenever the matrix is an upper or lower triangular matrix.
A = matrix(c(4,0,1,0,4,0,1,0,4),nrow=3)
A## [,1] [,2] [,3]
## [1,] 4 0 1
## [2,] 0 4 0
## [3,] 1 0 4A%*%t(A) ## [,1] [,2] [,3]
## [1,] 17 0 8
## [2,] 0 16 0
## [3,] 8 0 17t(A)%*%A## [,1] [,2] [,3]
## [1,] 17 0 8
## [2,] 0 16 0
## [3,] 8 0 17B = matrix(c(1,0,0,1,2,0,0,0,3),nrow=3)
B## [,1] [,2] [,3]
## [1,] 1 1 0
## [2,] 0 2 0
## [3,] 0 0 3B%*%t(B) ## [,1] [,2] [,3]
## [1,] 2 2 0
## [2,] 2 4 0
## [3,] 0 0 9t(B)%*%B## [,1] [,2] [,3]
## [1,] 1 1 0
## [2,] 1 5 0
## [3,] 0 0 9C = matrix(c(1,0,0,4,4,0,4,5,3),nrow=3)
C## [,1] [,2] [,3]
## [1,] 1 4 4
## [2,] 0 4 5
## [3,] 0 0 3C%*%t(C) ## [,1] [,2] [,3]
## [1,] 33 36 12
## [2,] 36 41 15
## [3,] 12 15 9t(C)%*%C## [,1] [,2] [,3]
## [1,] 1 4 4
## [2,] 4 32 36
## [3,] 4 36 50Matrix factorization is a very important problem. There are supercomputers built just to do matrix factorizations. Every second you are on an airplane, matrices are being factorized. Radars that track flights use a technique called Kalman filtering. At the heart of Kalman Filtering is a Matrix Factorization operation. Kalman Filters are solving linear systems of equations when they track your flight using radars. Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer. Please submit your response in an R Markdown document using our class naming convention, E.g. LFulton_Assignment2_PS2.png You don’t have to worry about permuting rows of A and you can assume that A is less than 5x5, if you need to hard-code any variables in your code. If you doing the entire assignment in R, then please submit only one markdown document for both the problems.
Here’s a function for only 2x2 matrices:
factor_matrix_2x2 <- function(matrix){
#Assume we're dealing with a square matrix
L = matrix(c(1,0,0,1),nrow=2)
#Make the position 2,1 a 0 in matrix
U = matrix
U[2,] = (matrix[2,1]/ matrix[1,1])*matrix[1,] - matrix[2,1]
L[2,1] = matrix[2,1]/ matrix[1,1]
print("Matrix L is:")
print(L)
print("Matrix U is:")
print(U)
}
A = matrix(c(1,-2,4,8),nrow=2)
factor_matrix_2x2(A)## [1] "Matrix L is:"
## [,1] [,2]
## [1,] 1 0
## [2,] -2 1
## [1] "Matrix U is:"
## [,1] [,2]
## [1,] 1 4
## [2,] 0 -6Here’s a function that can handles matrices up to size 5x5 assuming they are square matrices:
factor_matrix <- function(matrix){
#Assume we're dealing with a square matrix with 5 or less size but greater than 1
size <- dim(matrix)[1]
A <- matrix
L <- matrix(0, size, size)
diag(L) <- 1
U = A
#(2x2 matrix)
if (size > 1){
L[2,1] = (U[2,1] / U[1,1])
U[2,] = -(U[2,1] / U[1,1])*U[1,] + U[2,]
}
#(3x3 matrix)
if (size > 2){
L[3,1] = (U[3,1] / U[1,1])
U[3,] = -(U[3,1] / U[1,1])*U[1,] + U[3,]
L[3,2] = (U[3,2] / U[2,2])
U[3,] = -(U[3,2] / U[2,2])*U[2,] + U[3,]
}
#(4x4 matrix)
if (size > 3){
L[4,1] = (U[4,1] / U[1,1])
U[4,] = -(U[4,1] / U[1,1])*U[1,] + U[4,]
L[4,2] = (U[4,2] / U[2,2])
U[4,] = -(U[4,2] / U[2,2])*U[2,] + U[4,]
L[4,3] = (U[4,3] / U[3,3])
U[4,] = -(U[4,3] / U[3,3])*U[3,] + U[4,]
}
#(5x5 matrix)
if (size > 4){
L[5,1] = (U[5,1] / U[1,1])
U[5,] = -(U[5,1] / U[1,1])*U[1,] + U[5,]
L[5,2] = (U[5,2] / U[2,2])
U[5,] = -(U[5,2] / U[2,2])*U[2,] + U[5,]
L[5,3] = (U[5,3] / U[3,3])
U[5,] = -(U[5,3] / U[3,3])*U[3,] + U[5,]
}
print("Matrix L is:")
print(L)
print("Matrix U is:")
print(U)
}
B = matrix(c(1,-2,3,4,8,4,-3,5,7),nrow=3)
factor_matrix(B)## [1] "Matrix L is:"
## [,1] [,2] [,3]
## [1,] 1 0.0 0
## [2,] -2 1.0 0
## [3,] 3 -0.5 1
## [1] "Matrix U is:"
## [,1] [,2] [,3]
## [1,] 1 4 -3.0
## [2,] 0 16 -1.0
## [3,] 0 0 15.5C = matrix(c(2,1,-6,4,-4,-9,-4,3,5),nrow=3)
factor_matrix(C)## [1] "Matrix L is:"
## [,1] [,2] [,3]
## [1,] 1.0 0.0 0
## [2,] 0.5 1.0 0
## [3,] -3.0 -0.5 1
## [1] "Matrix U is:"
## [,1] [,2] [,3]
## [1,] 2 4 -4.0
## [2,] 0 -6 5.0
## [3,] 0 0 -4.5D = matrix(c(1,2,0,0,5,12,4,0,0,5,13,6,0,0,5,11),nrow=4)
factor_matrix(D)## [1] "Matrix L is:"
## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 2 1 0 0
## [3,] 0 2 1 0
## [4,] 0 0 2 1
## [1] "Matrix U is:"
## [,1] [,2] [,3] [,4]
## [1,] 1 5 0 0
## [2,] 0 2 5 0
## [3,] 0 0 3 5
## [4,] 0 0 0 1