library(tinytex)
## Warning: package 'tinytex' was built under R version 4.3.2
library(matrixcalc)

##1. Problem set 1:

  1. Show that \(A^TA \neq AA^T\) in general. (Proof and demonstration)
A <- matrix(c(2,5,-1,1,0,3,-4,2,-2), nrow = 3, byrow = T)
At <- t(A)
At
##      [,1] [,2] [,3]
## [1,]    2    1   -4
## [2,]    5    0    2
## [3,]   -1    3   -2

\(A^TA\)

At %*% A
##      [,1] [,2] [,3]
## [1,]   21    2    9
## [2,]    2   29   -9
## [3,]    9   -9   14

\(AA^T\)

A %*% At
##      [,1] [,2] [,3]
## [1,]   30   -1    4
## [2,]   -1   10  -10
## [3,]    4  -10   24
  1. 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).
A2 <- matrix(c(1,0,0,0,1,0,0,0,1), nrow = 3, byrow = T)
A2
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

Here the transposed matrix A is the same as the regular A matrix

At2 <- t(A2)
At2
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

In this case \(A^TA = AA^T\) is true. A diagonally symmetrical matrix will result in the instance where \(A^TA = AA^T\)

At2 %*% A2
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
A2 %*% At2
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1

I can also try this with a different matrix instead of the identity matrix

B <- matrix(c(1,2,3,2,4,5,3,5,6), byrow = T, nrow = 3)

Check to see if my matrix is ‘Symmetrical’

B == t(B)
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE
Bt = t(B)
BtB <- Bt %*% B
BtB
##      [,1] [,2] [,3]
## [1,]   14   25   31
## [2,]   25   45   56
## [3,]   31   56   70
BBt <- B %*% Bt
BBt
##      [,1] [,2] [,3]
## [1,]   14   25   31
## [2,]   25   45   56
## [3,]   31   56   70

The same thing applies to other symmetrical non-identity matrices. \(B^TB = BB^T\)

BBt == BtB
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE

Problem Set 2:

Matrix 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 fights 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.

#LU Function

LU <- function(A){
    n <- nrow(A)
    U <- A
    L <- diag(n)
    
    for(i in 2:n){
        for(r in 1:(i-1)){ #Iteration pattern over the rows and columns
            mult <- -U[i,r] / U[r,r]
            U[i, ] <- mult * U[r, ] + U[i, ]
            L[i,r] <- -mult #adding the inverse multipliers into the lower triangular matrix (L)
        }
    }
    if (!all.equal(A, L %*% U)) {
        warning("L * U is not equal to A.") ##Fail safe if the function was not set up properly
    
}
    return(list(L,U)) #This will return the L and U matrices individually (in a list format) so they can be subset and tested to see if LU == A
}

Factorization application to A matrix

A <- matrix(c(1,3,-4,2,5,2,1,-1,6), nrow = 3, byrow = T)
LU(A)
## [[1]]
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    2    1    0
## [3,]    1    4    1
## 
## [[2]]
##      [,1] [,2] [,3]
## [1,]    1    3   -4
## [2,]    0   -1   10
## [3,]    0    0  -30

Setting L as the Lower Triangualar Matrix

L <- LU(A)[[1]]
L
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    2    1    0
## [3,]    1    4    1

Setting U as the Upper Triangualar Matrix

U <- LU(A)[[2]]
U
##      [,1] [,2] [,3]
## [1,]    1    3   -4
## [2,]    0   -1   10
## [3,]    0    0  -30

Does A = LU?

L %*% U
##      [,1] [,2] [,3]
## [1,]    1    3   -4
## [2,]    2    5    2
## [3,]    1   -1    6

While I have a checker built into the function I want to double check manually (the function returns an error message if they are not equal)

A == L%*%U
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE

Now that I know my function does decomposition I can try and try a LU package with my matrix and see if the values are the same. I will be using the MatrixCalc package

lu_decomp <- lu.decomposition(A)

I will follow the same process in vetting this information as I did above

A == lu_decomp$L %*% lu_decomp$U
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE