ASSIGNMENT 2 IS 605 FUNDAMENTALS OF COMPUTATIONAL MATHEMATICS- 2023

  1. Problem set 1 (1) Show that AT A = AAT in general. (Proof and demonstration.)

(2) For a special type of square matrix A, we get AT A = AAT . Under what conditions could this be true? (Hint: The Identity matrix I is an example of such a matrix). Please typeset your response using LaTeX mode in RStudio. If you do it in paper, please either scan or take a picture of the work and submit it. Please ensure that your image is legible and that your submissions are named using your first initial, last name, assignment and problem set within the assignment. E.g. LFulton_Assignment2_PS1.png

The statement “A^T * A = A * A^T” is not true in general. This statement only holds true if the matrix A is square, meaning that it has the same number of rows and columns.

Let A be an m x n matrix. The transpose of A, denoted as A^T, is an n x m matrix where the elements are obtained by swapping the rows and columns of A. The product of A^T and A is given by A^T * A, which is an n x n matrix. The product of A and A^T is given by A * A^T, which is an m x m matrix.

a <- matrix(1:12, ncol = 3, nrow = 4)

at<-t(a)

aat<-a%*%at

aat
##      [,1] [,2] [,3] [,4]
## [1,]  107  122  137  152
## [2,]  122  140  158  176
## [3,]  137  158  179  200
## [4,]  152  176  200  224
ata<-at%*%a

ata
##      [,1] [,2] [,3]
## [1,]   30   70  110
## [2,]   70  174  278
## [3,]  110  278  446

Like stated above, for a square matrix, the ATA = AAT only if a matrix is symmetrical

a <- matrix(c(1,2,3,2,4,5,3,5,6), ncol = 3, nrow = 3)

a
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]    2    4    5
## [3,]    3    5    6
at<-t(a)

aat<-a%*%at

aat
##      [,1] [,2] [,3]
## [1,]   14   25   31
## [2,]   25   45   56
## [3,]   31   56   70
ata<-at%*%a

ata
##      [,1] [,2] [,3]
## [1,]   14   25   31
## [2,]   25   45   56
## [3,]   31   56   70

If a square matrix is not symetrical, then equality won’t hold.

a <- matrix(1:9, ncol = 3, nrow = 3)

at<-t(a)

aat<-a%*%at

aat
##      [,1] [,2] [,3]
## [1,]   66   78   90
## [2,]   78   93  108
## [3,]   90  108  126
ata<-at%*%a

ata
##      [,1] [,2] [,3]
## [1,]   14   32   50
## [2,]   32   77  122
## [3,]   50  122  194
  1. 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 flights 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.

The code below is an implementation of the LU decomposition method, which is a way of factorizing a square matrix into the product of two matrices - a lower triangular matrix (L) and an upper triangular matrix (U).

The code first checks if the input matrix (A) is a square matrix, and returns NULL if it’s not. Then it initializes the two matrices L and U and the dimensions of the input matrix.

If the input matrix is a 1x1 matrix, the function returns the diagonal matrix and the input matrix as the output.

The code then uses a nested for-loop to perform the LU decomposition algorithm. The loop calculates a multiplier mx for each iteration and updates the corresponding row of the matrix M (which will be used to represent U). The diagA matrix is updated to store the multipliers.

Finally, the function returns a list of two matrices - the lower triangular matrix L (diagA) and the upper triangular matrix U (M).

luDecomposer <- function(A) 
{
  # if matrix A isn't a square matrix then it won't work; return null
  if (dim(A)[1]!=dim(A)[2]) 
  {
    return(NULL)
  }
  
  M <- A
  dimOfA <- dim(A)[1]
  diagA <- diag(dimOfA)
  #if matrix A is a 1x1 matrix, return the diagonal matrix and the input matrix as the output.
  if (dimOfA == 1) 
  {
    return(list(diagA, M))
  }
  
  for(a in 2:dimOfA)
  {
    for(b in 1:(a-1)) 
    {
      mx <- (-M[a,b] / M[b,b]) 
      M[a, ] <- mx * M[b, ] + M[a, ]
      diagA[a,b] <- -mx
    }
  }
  return(list(diagA,M))
}

Let’s check if it works:

Example 1:

(M <- matrix(runif(12), nrow = 4, ncol = 3))
##           [,1]      [,2]      [,3]
## [1,] 0.8738831 0.1624580 0.4033208
## [2,] 0.7613575 0.9201777 0.4671804
## [3,] 0.2365978 0.3061699 0.3526402
## [4,] 0.1016645 0.1006549 0.4521525
luDecomposer(M)
## NULL

Perfect! It wasn’t a square matrix and therefore returned NULL.

Example 2:

(A <- matrix(runif(9), nrow = 3, ncol = 3))
##           [,1]      [,2]       [,3]
## [1,] 0.7246337 0.8414531 0.07579553
## [2,] 0.7204475 0.8836921 0.31064713
## [3,] 0.5447012 0.1069917 0.36831274
luDecomposed <- luDecomposer(A)
L<-luDecomposed[[1]]
U<-luDecomposed[[2]]
A == L %*% U
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE