Problem 1

Show that ATA != AAT in general. (Proof and demonstration.)

Here is an example to demonstrate that ATA != AAT.
The below 3 x 3 matrix A will be used to demonstrate ATA != AAT.

library(matrixcalc)

A<- matrix(c(2,-3, 4, 1, 0,1,2,1,2), nrow = 3, byrow = TRUE)
A
##      [,1] [,2] [,3]
## [1,]    2   -3    4
## [2,]    1    0    1
## [3,]    2    1    2

Below, matrix A is transposed.

At <- t(A)
At
##      [,1] [,2] [,3]
## [1,]    2    1    2
## [2,]   -3    0    1
## [3,]    4    1    2

As demonstrated below, AAt and AtA generally do not produce identical matrix. This because transposing matrix results in different coordinates in the transposed matrix. For instance Matrix A \(\begin{smallmatrix}a & b \\c & d\end{smallmatrix}\) transposed equals \(\begin{smallmatrix}a & c \\b & d\end{smallmatrix}\). The order of operations in matrix multiplication will result in different products since the tranposed matrix has different row and column values.

A%*%At
##      [,1] [,2] [,3]
## [1,]   29    6    9
## [2,]    6    2    4
## [3,]    9    4    9
At%*%A
##      [,1] [,2] [,3]
## [1,]    9   -4   13
## [2,]   -4   10  -10
## [3,]   13  -10   21

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

ATA = AAT when the matrix is symmetrical along it’s diagonal. The transpose matrix is identical to the original matrix, so the multiplying the matrix by it’s transpose is equivalent of the matrix multiplied by itself.

Below ATA = AAT is demonstrated using the following matrix:

A<- matrix(c(2,1, 0, 1, 2, 1 ,0 , 1, 2), nrow = 3, byrow = FALSE)
A
##      [,1] [,2] [,3]
## [1,]    2    1    0
## [2,]    1    2    1
## [3,]    0    1    2

Here is the matrix transposed (notice it is identitcal to matrix A):

At <- t(A)
At
##      [,1] [,2] [,3]
## [1,]    2    1    0
## [2,]    1    2    1
## [3,]    0    1    2

Below AtA=AAt because the matrix A is identical along it’s diagonal.

At%*%A
##      [,1] [,2] [,3]
## [1,]    5    4    1
## [2,]    4    6    4
## [3,]    1    4    5
A%*%At
##      [,1] [,2] [,3]
## [1,]    5    4    1
## [2,]    4    6    4
## [3,]    1    4    5

Problem 3

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

The below matrix will be used for decomposition:

A<- matrix(c(3, 1, 2, 1, 3, 1, 1, 1, 2), nrow = 3, byrow = TRUE)
A
##      [,1] [,2] [,3]
## [1,]    3    1    2
## [2,]    1    3    1
## [3,]    1    1    2

Location [2,1] needs to equal zero for decomposition. The below matrix is used to accomplish this. Location [2,1] in the below matrix E21 is set to -1/3 so location [2,1] in matrix A results in zero after matrix multiplication.

     E21=matrix(c(1,-1/3,0,0,1,0,0,0,1),nrow=3)
     E21 %*% A
##      [,1]     [,2]      [,3]
## [1,]    3 1.000000 2.0000000
## [2,]    0 2.666667 0.3333333
## [3,]    1 1.000000 2.0000000

Next, loction [3,1] needs to equal zero to accomplish decomposition. Location [3,1] in the below matrix E31 is set to -1/3 accordingly. 

     E31=matrix(c(1,0,-1/3,0,1,0,0,0,1),nrow=3)
     E31 %*% E21 %*% A
##      [,1]      [,2]      [,3]
## [1,]    3 1.0000000 2.0000000
## [2,]    0 2.6666667 0.3333333
## [3,]    0 0.6666667 1.3333333

Finally, location [3,2] needs to equal zero to accomplish decomposition. Location [3,2] in the below matrix E32 is set to -0.6666667/2.6666667 accordingly.

     E32=matrix(c(1,0,0,0,1,-0.6666667/2.6666667,0,0,1),nrow=3)
     E32
##      [,1]  [,2] [,3]
## [1,]    1  0.00    0
## [2,]    0  1.00    0
## [3,]    0 -0.25    1
     E32 %*%E31 %*% E21 %*% A
##      [,1]         [,2]      [,3]
## [1,]    3  1.000000000 2.0000000
## [2,]    0  2.666666667 0.3333333
## [3,]    0 -0.000000025 1.2500000

Making a function to accomplish decomposition into LU.

LU<-function(x) {
  b<-E21 %*% x
  c<-E31 %*% E21 %*% x
  U=E32%*%E31 %*% E21 %*% x
  U
  L<-solve(E21) %*% solve(E31) %*% solve(E32)
  L%*%U
  (L%*%U==x)
  result <- list( L=L, U=U )
     return( result )
  }

LU(A)
## $L
##           [,1] [,2] [,3]
## [1,] 1.0000000 0.00    0
## [2,] 0.3333333 1.00    0
## [3,] 0.3333333 0.25    1
## 
## $U
##      [,1]         [,2]      [,3]
## [1,]    3  1.000000000 2.0000000
## [2,]    0  2.666666667 0.3333333
## [3,]    0 -0.000000025 1.2500000

Below, the decomposition function is used to validate that LU above are correct.
The results are for the most part identical, deviating at location [3,2] by -0.000000025.

lu.decomposition(A)
## $L
##           [,1] [,2] [,3]
## [1,] 1.0000000 0.00    0
## [2,] 0.3333333 1.00    0
## [3,] 0.3333333 0.25    1
## 
## $U
##      [,1]     [,2]      [,3]
## [1,]    3 1.000000 2.0000000
## [2,]    0 2.666667 0.3333333
## [3,]    0 0.000000 1.2500000