library(knitr)

knitr::opts_chunk$set(echo = FALSE)
knitr::opts_knit$set(root.dir= normalizePath('..'))
knitr::opts_chunk$set(error = FALSE)

Problem Set 1

1. Show that A^T.A NOT EQUAL TO A.A^T in general.

Solution:

Let \(M = \begin{bmatrix}a & b\\ c & d\end{bmatrix}.\) Then \(M^T = \begin{bmatrix}a & c\\ b & d\end{bmatrix}\).

The product of \(MM^T\) will be: \(\begin{bmatrix}a\times a + b\times b & a\times c + b\times d\\ c\times a + d\times b & c\times c + d \times d\end{bmatrix}\).

The product of \(M^TM\) will be: \(\begin{bmatrix}a\times a + c\times c & a\times b + c\times d\\ a\times b + c\times d & b\times b + d \times d\end{bmatrix}\).

From above example, we can see that matrix multiplication does not have the commutative property. The product of \(MM^T\) is different from the product of \(M^TM\)

2. For a special type of square matrix \(A\), we get \(A^TA = AA^T\). Under what conditions could this be true? (Hint: The identity matrix \(I\) is an example of such a matrix.)

Solution

Theorem: For any matrix \(A\), \(AA^T\) and \(A^TA\) are symmetric matrices.
A symmetric matrix is a square matrix that is equal to its transpose The product of two symmetric matrices is symmetric. i.e. \(AA^T = A^TA\).

Problem Set 2

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

Solution:

Helper Functions

# get dimensions of given matrix
getdim = function(matrix){
  m = nrow(matrix)
  n = ncol(matrix)
  size = c(m, n)
  size
}


# set identity matirx with given dimension
identity = function(size)
  I = diag(size[1])

function to solve A = LU

getLU = function(M){
  
  #find dimensions
  dimension = getdim(M)
  
  if(dimension[1] != dimension[2]){
    return("Not a square matrix")
  }
  
  #initialize L as identity matrix
  L = identity(dimension)
  
  #initialize U as a copy of M
  U = M
  print(M)
  
  for (i in 2:dimension[1]){
    for(j in 1:(i - 1)){
        # initialize a new E for each iteration
        E <-  identity(dimension)
        # Find E based on state of U
        E[i, j] <- -(U[i,j]/U[j,j])
        
        # update U using E
        U <- E %*% U
        # update L based on E^-1
        L <- L %*% solve(E)
    }
  }

  return(list(L,U))
}

Function test

a = matrix(c(1,2,4,9), 2, byrow=T)
res = getLU(a)
##      [,1] [,2]
## [1,]    1    2
## [2,]    4    9
print("L:")
## [1] "L:"
print(res[[1]])
##      [,1] [,2]
## [1,]    1    0
## [2,]    4    1
print("U:")
## [1] "U:"
print(res[[2]])
##      [,1] [,2]
## [1,]    1    2
## [2,]    0    1
print("LU:")
## [1] "LU:"
LU <- res[[1]] %*% res[[2]] 
print(LU)
##      [,1] [,2]
## [1,]    1    2
## [2,]    4    9
print("A:")
## [1] "A:"
print(a)
##      [,1] [,2]
## [1,]    1    2
## [2,]    4    9