Computational Math and Matrix Factorization
Problem Set 1
- Show that A^T A≠AA^T in general. (Proof and demonstration.)
- For a special type of square matrix A, we get A^T A=AA^T. Under what conditions could this be true? (Hint: The Identity matrix I is an example of such a matrix).
Load Packages
library(knitr)The proofs were done using Word’s own Equation editor. After writing the proof, the screen shot of the word doc was written on to a ProblemSet_1.png file, which is inserted below:
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 fight 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.
Answer:
I have taken the following 5 x 5 matrix, for processing.
[,1] [,2] [,3] [,4] [,5]
[1,] 1 7 4 1 2
[2,] 5 4 1 2 3
[3,] 4 2 5 3 2
[4,] 2 5 2 3 1
[5,] 3 5 4 4 2 Please find below the solution to the above problem. First I wrote a function my_gauss_jordan_transformation(), which accepts two parameters–a matrix and the constraints. The comments explain the code. Right below the function, I made a call to the function, using the matrix and constraints of the given equation (1).
matrix_triangulation <- function(m)
{
#
# Setting up variables.
#
cols <- c(1:nrow(m))
rows <- c(1:ncol(m))
dynamic_rows <- rows
num_of_rows <- nrow(m)
L = diag(1, num_of_rows, num_of_rows)
diag_flag = FALSE
#
# Processing starts here.
#
for(col in cols) {
for(row in dynamic_rows) {
if( (m[row, col]) == 0 ) {
if(row == col) {
diag_flag = TRUE
diag = row
}
next
} else {
k = row
#
if(diag_flag) {
m[c(diag, k),] <- m[c(k, diag),]
temp = k
k = diag
diag = temp
diag_flag = FALSE
}
dynamic_rows <- dynamic_rows [!dynamic_rows %in% c(k)]
#
l <- k + 1
for(row in k:num_of_rows) {
if((row != k) & (m[row, col] != 0)) {
multiplier = m[row, col] / m[k, col]
m[row,] <- m[row,] + (-1) * m[k,] * multiplier
L[row, col] <- multiplier
}
}
#
break
}
}
}
U = m
result <- list(Lower = L, Upper = U)
return(result)
}In the following, I made call to the above above function matrix_triangulation, with the 5 x 5 matrix.
m <- matrix(c(1, 5, 4, 2, 3, 7, 4, 2, 5, 5, 4, 1, 5, 2, 4, 1, 2, 3, 3, 4, 2, 3, 2, 1, 2), nrow = 5)
## m <- matrix(c(0, 5, 4, 2, 3, 0, 4, 2, 5, 5, 4, 1, 5, 2, 4, 1, 2, 3, 3, 4, 2, 3, 2, 1, 2), nrow = 5)
#
print('Original matrix:')## [1] "Original matrix:"print(m)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 7 4 1 2
## [2,] 5 4 1 2 3
## [3,] 4 2 5 3 2
## [4,] 2 5 2 3 1
## [5,] 3 5 4 4 2#
processed_m <- matrix_triangulation(m)
#
print('Processed matrices:')## [1] "Processed matrices:"print('L-matrix:')## [1] "L-matrix:"print(processed_m$Lower)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 0.0000000 0.00000000 0.000000 0
## [2,] 5 1.0000000 0.00000000 0.000000 0
## [3,] 4 0.8387097 1.00000000 0.000000 0
## [4,] 2 0.2903226 -0.09803922 1.000000 0
## [5,] 3 0.5161290 0.36601307 0.987055 1print('U-matrix:')## [1] "U-matrix:"print(processed_m$Upper)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 7 4.000000 1.000000 2.0000000
## [2,] 0 -31 -19.000000 -3.000000 -7.0000000
## [3,] 0 0 4.935484 1.516129 -0.1290323
## [4,] 0 0 0.000000 2.019608 -0.9803922
## [5,] 0 0 0.000000 0.000000 0.6278317Now, I’ll test whether the factors indeed yield the original matrix, as product.
print('Test whether product of Lower and Upper is the original given matrix:')## [1] "Test whether product of Lower and Upper is the original given matrix:"product <- processed_m$Lower %*% processed_m$Upper
print(product)## [,1] [,2] [,3] [,4] [,5]
## [1,] 1 7 4 1 2
## [2,] 5 4 1 2 3
## [3,] 4 2 5 3 2
## [4,] 2 5 2 3 1
## [5,] 3 5 4 4 2Marker: 605-02