Problem Set 1

(1) Show that \(A^TA \ne AA^T\) in general.

If \(A = \left[ \begin{array}{ccc} 1 & 1 & 2 \\ 2 & 1 & 0\\3 & 1 & 1\end{array}\right]\), then \(A^T = \left[ \begin{array}{ccc} 1 & 2 & 3 \\ 1 & 1 & 1\\2 & 0 & 1\end{array}\right]\) (which is kind of general)

\(AA^T = \left[ \begin{array}{ccc} 1 & 1 & 2 \\ 2 & 1 & 0\\3 & 1 & 1\end{array}\right] \left[ \begin{array}{ccc} 1 & 2 & 3 \\ 1 & 1 & 1\\2 & 0 & 1\end{array}\right] = \left[ \begin{array}{ccc} 6 & 3 & 6 \\ 3 & 5 & 7 \\ 6 & 7 & 11 \end{array}\right]\)

\(A^TA = \left[ \begin{array}{ccc} 1 & 2 & 3 \\ 1 & 1 & 1\\2 & 0 & 1\end{array}\right] \left[ \begin{array}{ccc} 1 & 1 & 2 \\ 2 & 1 & 0\\3 & 1 & 1\end{array}\right] = \left[ \begin{array}{ccc} 14 & 6 & 5 \\ 6 & 3 & 3 \\ 5 & 3 & 5 \end{array}\right]\)

\(\left[ \begin{array}{ccc} 14 & 6 & 5 \\ 6 & 3 & 3 \\ 5 & 3 & 5 \end{array}\right] \ne \left[ \begin{array}{ccc} 6 & 3 & 6 \\ 3 & 5 & 7 \\ 6 & 7 & 11 \end{array}\right]\), so \(A^TA \ne AA^T\)

Since we can only multiply two matrices together when the row and column count are correct, \(A^TA\) is only going to work with a square matrix. We tried an abritary square 3 X 3 matrix and found \(A^TA \ne AA^T\). That leaves trying to find a case when it is equal.

If \(B = \left[ \begin{array}{ccc} 3 & 0 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 3\end{array}\right]\), then \(B^T = \left[ \begin{array}{ccc} 3 & 0 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 3\end{array}\right]\)

\(BB^T = \left[ \begin{array}{ccc} 3 & 0 & 0 \\ 0 & 3 & 0\\ 0 & 0 & 3\end{array}\right] \left[ \begin{array}{ccc} 3 & 0 & 0 \\ 0 & 3 & 0 \\ 0 & 0 & 3\end{array}\right] = \left[ \begin{array}{ccc} 9 & 0 & 0 \\ 0 & 9 & 0 \\ 0 & 0 & 9 \end{array}\right] = B^TB\)

So, \(A^TA = AA^T\) only when A is some multiple of a square identity matrix, with any non-zero values along the diagonal. This means \(A^TA \ne AA^T\) in general.

A <- matrix(c(1, 2, 3, 1, 1, 1, 2, 0, 1), nrow=3)
B <- matrix(c(3, 0, 0, 0, 3, 0, 0, 0, 3), nrow=3)
C <- matrix(c(1, 0, 0, 0, 2, 0, 0, 0, 3), nrow=3)
D <- matrix(c(21, 0, 0, 0, 5, 0, 0, 0, 3.14), nrow=3)
A
##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    2    1    0
## [3,]    3    1    1
t(A)
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    1    1    1
## [3,]    2    0    1
A%*%t(A)
##      [,1] [,2] [,3]
## [1,]    6    3    6
## [2,]    3    5    7
## [3,]    6    7   11
t(A)%*%A
##      [,1] [,2] [,3]
## [1,]   14    6    5
## [2,]    6    3    3
## [3,]    5    3    5
B
##      [,1] [,2] [,3]
## [1,]    3    0    0
## [2,]    0    3    0
## [3,]    0    0    3
t(B)
##      [,1] [,2] [,3]
## [1,]    3    0    0
## [2,]    0    3    0
## [3,]    0    0    3
B%*%t(B)
##      [,1] [,2] [,3]
## [1,]    9    0    0
## [2,]    0    9    0
## [3,]    0    0    9
t(B)%*%B
##      [,1] [,2] [,3]
## [1,]    9    0    0
## [2,]    0    9    0
## [3,]    0    0    9
C
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    2    0
## [3,]    0    0    3
t(C)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    2    0
## [3,]    0    0    3
C%*%t(C)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    4    0
## [3,]    0    0    9
t(C)%*%C
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    4    0
## [3,]    0    0    9
D
##      [,1] [,2] [,3]
## [1,]   21    0 0.00
## [2,]    0    5 0.00
## [3,]    0    0 3.14
t(D)
##      [,1] [,2] [,3]
## [1,]   21    0 0.00
## [2,]    0    5 0.00
## [3,]    0    0 3.14
D%*%t(D)
##      [,1] [,2]   [,3]
## [1,]  441    0 0.0000
## [2,]    0   25 0.0000
## [3,]    0    0 9.8596
t(D)%*%D
##      [,1] [,2]   [,3]
## [1,]  441    0 0.0000
## [2,]    0   25 0.0000
## [3,]    0    0 9.8596

Problem Set 2

(2) Write an R function to factorize a square matrix A into LU or LDU, whichever you prefer.

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.

ldu <- function(A) {
  # Factorize a square matrix A into LU, where L is a Lower Triangular matrix and U is a Upper Triangular matrix

  arows <- dim(A)[1]
  acols <- dim(A)[2]
  adim <- dim(A)[1]
  # Final version of funtion should check for square matrix
  # if arows <> acols Print "Not a Square matrix. Cannot run." Return
  
  # Final version of funtion should check for zeros on diagonals and swap rows when we have a zero
  # This swaps row 1 to bottom row 7 for matrix A, A[ c(2:7,1), ]
  
  # Need a counter for the eliminations and to initialize our variables
  k <- 1
  U <- A
  L <- diag(adim)
  id <- diag(adim)
  
  # Make a list to store our inverses for L and our results
  ls <- rep(list(diag(adim)), adim)
  results <- rep(list(diag(adim)), 2)
  
  # Go through all the columns
  for (j in 1:(acols-1)) {
    
    # Go through all the rows
    for (i in j:(arows-1)) {
      M <- id
      M[i+1, j] <- -U[i+1, j]
      U <- M %*% U
      ls[[k]] <- solve(M)
      k <- k + 1
    }
  }
  # Why do I have to use a for-loop? I could not find a better way.
  for (m in 1:adim) {
    L <- L %*% ls[[m]]
  }
  # Need a prettier way to return these results
  results[[1]] <- L
  results[[2]] <- U
  print("L %*% U")
  print(L %*% U)
  print("L")
  print(L)
  print("U")
  print(U)
}

mat <- matrix(c(1, 2, 3, 1, 1, 1, 2, 0, 1),nrow=3)
ldu(mat)
## [1] "L %*% U"
##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    2    1    0
## [3,]    3    1    1
## [1] "L"
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    2    1    0
## [3,]    3   -2    1
## [1] "U"
##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    0   -1   -4
## [3,]    0   -4  -13
mat
##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    2    1    0
## [3,]    3    1    1