Problem Set 1

  1. Show that \(A^T A \neq A A^T\) in general. Please typeset your response using LaTeX mode in RStudio.

Proof

Let: \[\mathbf{A}=\left[\begin{array}{ll} a & b \\ c & d \end{array}\right]\]

Let: \[\mathbf{A}^{\mathbf{t}}=\left[\begin{array}{ll} w & x \\ y & z \end{array}\right]\]

\(A^T A\):

\[\mathbf{A}^{\mathbf{t}} \mathbf{A}=\left[\begin{array}{cc} w a+x c & w b+x d \\ y a+z c & y b+z d \end{array}\right]\]

\(A A^T\):

\[ \mathbf{A A}^{\mathbf{t}}=\left[\begin{array}{ll} a w+b y & a x+b z \\ c w+d y & c x+d z \end{array}\right] \]

Demonstration

# create matrix A
(A <- matrix(c(1,2,3,1,1,1,2,0,1),nrow=3))
##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    2    1    0
## [3,]    3    1    1
# transpose A
(A_t <- t(A))
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    1    1    1
## [3,]    2    0    1

Compare equations

(A_t%*%A == A%*%A_t)
##       [,1]  [,2]  [,3]
## [1,] FALSE FALSE FALSE
## [2,] FALSE FALSE FALSE
## [3,] FALSE FALSE FALSE
  1. For a special type of square matrix A, we get \(A^T A=A A^T\). Under what conditions could this be true? (Hint: The Identity matrix I is an example of such a matrix).

This would work for an Identity matrix:

# create identity matrix
(I <- diag(3))
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
# transpose identity matrix
(I_t <- t(diag(3)))
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

Compare

(I_t%*%I == I_t%*%I)
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE

Or any other symmetrical matrix

# create other symmetrical matrix
(O <- matrix(c(1,2,2,1), nrow = 2))
##      [,1] [,2]
## [1,]    1    2
## [2,]    2    1
# transpose
(O_t <- t(O))
##      [,1] [,2]
## [1,]    1    2
## [2,]    2    1
# compare
(O_t%*%O==O%*%O_t)
##      [,1] [,2]
## [1,] TRUE TRUE
## [2,] TRUE TRUE

Problem 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.

For this problem found a website explaining the Doolittle Algorithm: LU Decomposition Method which follows as such for a given square matrix \(A\):

Terms of U matrix are given by:

\[ \begin{aligned} & \forall j \\ & i=0 \rightarrow U_{i j}=A_{i j} \\ & i>0 \rightarrow U_{i j}=A_{i j}-\sum_{k=0}^{i-1} L_{i k} U_{k j} \end{aligned} \]

Terms of L matrix:

\[ \begin{aligned} & \forall i \\ & j=0 \rightarrow L_{i j}=\frac{A_{i j}}{U_{j j}} \\ & j>0 \rightarrow L_{i j}=\frac{A_{i j}-\sum_{k=0}^{j-1} L_{i k} U_{k j}}{U_{j j}} \end{aligned} \]

Write function to compute LU decomposition from square matrix.

matrix_LU <- function(a) {
  # n rows for matrix
  n <- nrow(a)
  # create upper and lower matrix
  U <- matrix(0, ncol = n, nrow = n) # same size as a all zeros
  # lower has to have 1's on diagonal
  L <- diag(x = 1, ncol = n, nrow = n) # identity matrix the size of a
  
  # loop through each row
  for (i in 1:n) {
    # create variables for i+1 and i-1
    i_p_1 <- i + 1
    i_m_1 <- i - 1
    for (j in 1:n) {
      # loop through matrix and copy a to U
      U[i,j] <- a[i,j]
      # if i-1 more than zero
      if (i_m_1 > 0) {
        # solve U matrix equation
        # iterate through U and L subtracting from L, U product second row
        # we know that row 1 will be the same in U 
        for (k in 1:i_m_1) {
          U[i,j] <- U[i,j] - L[i,k] * U[k,j]
        }
      }
    }
    if (i_p_1 <= n) {
      # solve L matrix equation
      # we know that upper portion is all zeros and diagonal are 1's
      # loop through j i+1 to n 
      for (j in i_p_1:n) {
        # copy lower portion of a to L
        L[j,i] <- a[j,i]
        # loop through rows starting at 2nd row
        if (i_m_1 > 0) {
          # solve lower half of matrix
          # for k in i-1 (starts on 2nd row)
          for (k in 1:i_m_1) {
            L[j,i] <- L[j,i] - L[j,k] * U[k,i]
          }
        }
        # divide portion of equation
        L[j,i] <- L[j,i]/ U[i,i]
      }
    } 
  }
  result <- list(L=L, U=U)
  return(result)
}

# try on matrix A from earlier
matrix_LU(A)
## $L
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    2    1    0
## [3,]    3    2    1
## 
## $U
##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    0   -1   -4
## [3,]    0    0    3