DATA 605 FUNDAMENTALS OF COMPUTATIONAL MATHEMATICS

1. Problem set 1

1. Show that \(A^TA \neq AA^T\) in general. (Proof and demonstration.)

Let A be a \(A^{1x1}\),

\(A = I_n\) (identity)

\(A = A^T\) (symmetric),

\(A^{-1} = A^T\) (orthogonal)

or \(AA^T = A^TA = diagonal_n\) (Both products are the same diagonal matrix).

Then \(A^TA = AA^T\).

Let A be a \(A^{1x1}\),

\(A \neq I_n\) (identity)

\(A \neq A^T\) (symmetric),

\(A^{-1} \neq A^T\) (orthogonal)

or \(AA^T \neq A^TA \neq diagonal_n\).

Then \(A^TA \neq AA^T\)

True, in the random sampling of small matrices, the ~86% of the 2 matrix products were not equal.

max_dim <- 3
max_range <- 2
iterations <- 20000

dim_range <- c(1:max_dim)
samp_space <- c(-max_range:max_range)

results_df <- data.frame(mat = character(), m = integer(), n = integer(), equality_test = logical(), 
    matrix_type = character(), stringsAsFactors = F)

matrix_type_test <- function(mat, tmat, m_dim, n_dim, A_result, B_result, rand_samp) {
    if (all(m_dim == 1, n_dim == 1)) {
        return("1 by 1 matrix")
    } else if (isTRUE(all.equal(mat, diag(x = 1, m_dim)))) {
        return("Is identity matrix")
    } else if (identical(mat, tmat)) {
        return("Is symmetric matrix")
    } else if (m_dim == n_dim && identical(mat %*% tmat, diag(m_dim))) {
        return("Is orthogonal matrix")
    } else if (identical(A_result, B_result) && identical(A_result, diag(A_result[1, 
        1], nrow(A_result)))) {
        return("Both products are the same diagonal matrix")
    } else {
        return("N/A")
    }
}


for (i in 1:iterations) {
    
    # select matrix dims
    m_dim <- sample(dim_range, 1)
    n_dim <- sample(dim_range, 1)
    samp_size <- m_dim * n_dim
    
    # select matrix elements
    rand_samp <- sample(samp_space, samp_size, replace = T)
    
    # create matrices
    current_matrix <- matrix(rand_samp, m_dim)
    t_current_matrix <- t(current_matrix)
    
    # multiply
    A_result <- t_current_matrix %*% current_matrix
    B_result <- current_matrix %*% t_current_matrix
    
    
    # test & store the results
    match_test <- identical(A_result, B_result)
    matrix_type <- matrix_type_test(current_matrix, t_current_matrix, m_dim, 
        n_dim, A_result, B_result, rand_samp)
    
    new_row <- list(c(paste(rand_samp, collapse = ",")), m_dim, n_dim, match_test, 
        matrix_type)
    
    results_df[i, ] <- new_row
    
}

prop.table(table(results_df$equality_test))

## 
##  FALSE   TRUE 
## 0.8628 0.1372

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

I identified 5 conditions for when the 2 matrix products are equal to each other.

• 1 by 1 matrices, which are also symmetric matrices

• Identity matrices

• Symmetric matrices

• Othogonal matrices

• When both matrix products produce the same diagonal matrix

t(prop.table(table(results_df[, 4:5])))

##                                             equality_test
## matrix_type                                    FALSE    TRUE
##   1 by 1 matrix                              0.00000 0.11300
##   Both products are the same diagonal matrix 0.00000 0.00355
##   Is identity matrix                         0.00000 0.00020
##   Is orthogonal matrix                       0.00000 0.00030
##   Is symmetric matrix                        0.00000 0.02380
##   N/A                                        0.85900 0.00015

2. 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.

LU_decomp <- function(A) {
    size <- nrow(A)
    U <- A
    I <- diag(size)
    myvector <- c()
    
    for (i in 2:size) {
        j_max <- i - 1
        for (j in 1:j_max) {
            elim_mat <- I
            pivot <- ifelse(U[i - 1, j] != 0, U[i - 1, j], U[i - 2, j])
            elim_mat[i, j] <- -U[i, j]/pivot
            (U <- elim_mat %*% U)
            var_name <- paste0("E", i, j)
            assign(var_name, elim_mat)
            myvector <- c(myvector, var_name)
        }
    }
    
    L <- diag(size)
    for (m in myvector) {
        L <- L %*% solve(get(m))
    }
    
    print(U)
    print(L)
    ifelse(identical(L %*% U, A), "The decomposition works", "The decomposition doesn't work")
}


A <- matrix(c(1, 2, 3, 1, 1, 1, 2, 0, 1), nrow = 3)
LU_decomp(A)

##      [,1] [,2] [,3]
## [1,]    1    1    2
## [2,]    0   -1   -4
## [3,]    0    0    3
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    2    1    0
## [3,]    3    2    1

## [1] "The decomposition works"

A <- matrix(c(1, -2, 3, 4, 8, 4, -3, 5, 7), 3)
LU_decomp(A)

##      [,1] [,2] [,3]
## [1,]    1    4 -3.0
## [2,]    0   16 -1.0
## [3,]    0    0 15.5
##      [,1] [,2] [,3]
## [1,]    1  0.0    0
## [2,]   -2  1.0    0
## [3,]    3 -0.5    1

## [1] "The decomposition works"