Problem Set 1 1. Show that \(A^TA \neq AA^T\) in general. (Proof and demonstration)
Proof: For an \(m \times n\) matrix \(A\), its transpose \(A^T\) is an \(n \times m\) matrix.
Therefore \(A^TA\) is an \(n \times n\) matrix, while \(AA^T\) is an \(m \times m\) matrix (according to the definition of matrix multiplication.) Therefore \(A^TA \neq AA^T\), in general, as long as \(m \neq n.\)
Demonstration:
A <- matrix(c(1,2,3,4,5,6), c(2,3), byrow = TRUE)
A_T <- t(A)
A_T %*% A## [,1] [,2] [,3]
## [1,] 17 22 27
## [2,] 22 29 36
## [3,] 27 36 45A %*% A_T## [,1] [,2]
## [1,] 14 32
## [2,] 32 77As expected, when we started with a \(2 \times 3\) matrix, the matrix multiplication of \(A\) and its transpose led to both a \(3 \times 3\) matrix and a \(2 \times 2\) matrix.
\(A^TA = AA^T \iff A\) is a symmetric matrix. A symmetric matrix is the same when it is “flipped” along its main diagonal. By definition, a symmetric matrix is a matrix which is equal to its transpose (\(A^T = A\),) so \(A^TA = AA^T\) is true because \(AA = AA\). Therefore \(A\) being symmetric \(\implies A^TA = AA^T\).
Conversely, \(A^TA = AA^T \implies A\) is symmetric. Since multiplication is not commutative in general for matrices, this would only hold if \(A^T = A\), meaning \(A\) is symmetric.
Example with a symmetric matrix:
S <- matrix(c(1,-2,9,-2,5,16,9,16,7), c(3,3), byrow = TRUE)
S_T <- t(S)
S_T %*% S## [,1] [,2] [,3]
## [1,] 86 132 40
## [2,] 132 285 174
## [3,] 40 174 386S %*% S_T## [,1] [,2] [,3]
## [1,] 86 132 40
## [2,] 132 285 174
## [3,] 40 174 386Problem Set 2
“Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer… 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.”
Factorization into LU:
The function below uses an if/else structure to find the lower triangular matrix and the upper triangular matrix of an \(n \times n\) matrix \(A\) where \(2 \leq n \leq 4\). Setting the value of \(n\) as well as the elements of the matrix \(A\) and running the code outputs both L and U, as well as evidence that \(LU = A\).
When creating each elimination matrix, the element to be eliminated is negated and divided by the element in the same column that is in the row to be added with the current row. Multiplying these elimination matrices by \(A\) gives the upper triangular matrix, since the product of the inverses of the elimination matrices is the lower triangular matrix.
#Define n as the size of the square matrix:
n <- 3
#Enter the matrix to be factorized by row:
A <- matrix(c(1,2,3,4,5,6,7,8,9),c(n,n), byrow= TRUE)
if (n == 2) {
E21 <- matrix(c(1,0,-A[2,1]/A[1,1],1),c(2,2), byrow= TRUE)
U <- E21 %*% A
L <- solve(E21)
print("The Lower Triangular Matrix:")
print(L)
print("The Upper Triangular Matrix:")
print(U)
print("Does LU equal A?")
print(L%*%U == A)
} else if (n == 3) {
E21 <- matrix(c(1,0,0,-A[2,1]/A[1,1],1,0,0,0,1),c(3,3), byrow= TRUE)
E31 <- matrix(c(1,0,0,0,1,0,-A[3,1]/A[1,1],0,1),c(3,3), byrow= TRUE)
E32 <- matrix(c(1,0,0,0,1,0,0,-A[3,2]/A[2,2],1),c(3,3), byrow= TRUE)
U <- E32 %*% E31 %*% E21 %*% A
L <- solve(E21) %*% solve(E31) %*% solve(E32)
print("The Lower Triangular Matrix:")
print(L)
print("The Upper Triangular Matrix:")
print(U)
print("Does LU equal A?")
print(L%*%U == A)
} else if (n == 4) {
E21 <- matrix(c(1,0,0,0,-A[2,1]/A[1,1],1,0,0,0,0,1,0,0,0,0,1),c(4,4), byrow= TRUE)
E31 <- matrix(c(1,0,0,0,0,1,0,0,-A[3,1]/A[1,1],0,1,0,0,0,0,1),c(4,4), byrow= TRUE)
E41 <- matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,-A[4,1]/A[1,1],0,0,1),c(4,4), byrow= TRUE)
E32 <- matrix(c(1,0,0,0,0,1,0,0,0,-A[3,2]/A[2,2],1,0,0,0,0,1),c(4,4), byrow= TRUE)
E42 <- matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,0,-A[4,2]/A[2,2],0,1),c(4,4), byrow= TRUE)
E43 <- matrix(c(1,0,0,0,0,1,0,0,0,0,1,0,0,0,-A[4,3]/A[3,3],1),c(4,4), byrow= TRUE)
U <- E43 %*% E42 %*% E32 %*% E41 %*% E31 %*% E21 %*% A
L <- solve(E21) %*% solve(E31) %*% solve(E41) %*% solve(E32) %*% solve(E42) %*% solve(E43)
print("The Lower Triangular Matrix:")
print(L)
print("The Upper Triangular Matrix:")
print(U)
print("Does LU equal A?")
print(L%*%U == A)
} else {
print("The given matrix is outside the bounds of this assignment.")
}## [1] "The Lower Triangular Matrix:"
## [,1] [,2] [,3]
## [1,] 1 0.0 0
## [2,] 4 1.0 0
## [3,] 7 1.6 1
## [1] "The Upper Triangular Matrix:"
## [,1] [,2] [,3]
## [1,] 1 2.0 3.0
## [2,] 0 -3.0 -6.0
## [3,] 0 -1.2 -2.4
## [1] "Does LU equal A?"
## [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUEAlgebraic Approach:
Initially, I tried to derive L and U algebraically. This was clunky but proved effective for both a \(2 \times 2\) matrix and a \(3 \times 3\) matrix, but gave a dividing by zero error for a \(4 \times 4\) matrix, probably due to human error.
For a \(2 \times 2\) matrix \[ \begin{bmatrix} a & b \\ c & d \\ \end{bmatrix} \]
a <- 1
b <- 2
c <- 3
d <- 4
A <- matrix(c(a,b,c,d), c(2,2), byrow= TRUE)
L <- matrix(c(1,0,c/a,1), c(2,2), byrow= TRUE)
U <- matrix(c(a,b,0,-b*c/a+d), c(2,2), byrow= TRUE)
print("The Lower Triangular Matrix:")## [1] "The Lower Triangular Matrix:"print(L)## [,1] [,2]
## [1,] 1 0
## [2,] 3 1print("The Upper Triangular Matrix:")## [1] "The Upper Triangular Matrix:"print(U)## [,1] [,2]
## [1,] 1 2
## [2,] 0 -2print("Does LU equal A?")## [1] "Does LU equal A?"print(L%*%U == A)## [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUEFor a \(3 \times 3\) matrix \[ \begin{bmatrix} a & b & c \\ d & e & f \\ g & h & i \\ \end{bmatrix} \]
a <- 1
b <- 2
c <- 3
d <- 4
e <- 5
f <- 6
g <- 7
h <- 8
i <- 9
A <- matrix(c(a,b,c,d,e,f,g,h,i),nrow=3, byrow= TRUE)
L <- matrix(c(1,0,0,d/a,1,0,g/a,-(-h+g*b/a)/(-d*b/a+e),1),c(3,3), byrow= TRUE)
U <- matrix(c(a,b,c,0,-d*b/a+e,-d*c/a+f,0,0,(-h+g*b/a)/(-d*b/a+e)*(-d*c/a+f)+-g*c/a+i), c(3,3), byrow= TRUE)
print("The Lower Triangular Matrix:")## [1] "The Lower Triangular Matrix:"print(L)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 4 1 0
## [3,] 7 2 1print("The Upper Triangular Matrix:")## [1] "The Upper Triangular Matrix:"print(U)## [,1] [,2] [,3]
## [1,] 1 2 3
## [2,] 0 -3 -6
## [3,] 0 0 0print("Does LU equal A?")## [1] "Does LU equal A?"print(L%*%U == A)## [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUEFor a \(4 \times 4\) matrix \[ \begin{bmatrix} a & b & c & d \\ e & f & g & h \\ i & j & k & l \\ m & n & o & p \\ \end{bmatrix} \]
a <- 1
b <- 2
c <- 3
d <- 4
e <- 5
f <- 6
g <- 7
h <- 8
i <- 9
j <- 10
k <- 11
l <- 12
m <- 13
n <- 14
o <- 15
p <- 16
A <- matrix(c(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p),c(4,4), byrow= TRUE)
L <- matrix(c(1,0,0,0,e/a,1,0,0,i/a,-(b*i/a-j)/(-b*e/a+f),1,0,m/a,-(b*m/a-n)/(-b*e/a+f),-(-((b*m/a-n)/(-b*e/a+f)*(-c*e/a+g)+-c*m/a+o))/((b*i/a-j)/(-b*e/a+f)*(-c*e/a+g)+-c*i/a+k),1),c(4,4), byrow= TRUE)
U <- matrix(c(a,b,c,d,0,-b*e/a+f,-c*e/a+g,-d*e/a+h,0,0,(b*i/a-j)/(-b*e/a+f)*(-c*e/a+g)+-c*i/a+k,(b*i/a-j)/(-b*e/a+f)*(-d*e/a+h)+-d*i/a+l,0,0,0,((-((b*m/a-n)/(-b*e/a+f)*(-c*e/a+g)+-c*m/a+o))/((b*i/a-j)/(-b*e/a+f)*(-c*e/a+g)+-c*i/a+k))*((b*i/a-j)/(-b*e/a+f)*(-d*e/a+h)+-d*i/a+l)+((b*m/a-n)/(-b*e/a+f)*(-d*e/a+h)+-d*m/a+p)), c(4,4), byrow= TRUE)
print("The Lower Triangular Matrix:")## [1] "The Lower Triangular Matrix:"print(L)## [,1] [,2] [,3] [,4]
## [1,] 1 0 0 0
## [2,] 5 1 0 0
## [3,] 9 2 1 0
## [4,] 13 3 NaN 1print("The Upper Triangular Matrix:")## [1] "The Upper Triangular Matrix:"print(U)## [,1] [,2] [,3] [,4]
## [1,] 1 2 3 4
## [2,] 0 -4 -8 -12
## [3,] 0 0 0 0
## [4,] 0 0 0 NaNprint("Does LU equal A?")## [1] "Does LU equal A?"print(L%*%U == A)## [,1] [,2] [,3] [,4]
## [1,] TRUE TRUE TRUE NA
## [2,] TRUE TRUE TRUE NA
## [3,] TRUE TRUE TRUE NA
## [4,] NA NA NA NA