Assignment 2

Fundamentals of Computational Mathematics

CUNY MSDS DATA 605, Fall 2018

Rose Koh

09/05/2018

# Let's assume that t(A) %*% A == A %*% t(A)
# Here is the testing matrix
A <- matrix(c(2,1,2,6), 2, byrow=T)
print(A)

##      [,1] [,2]
## [1,]    2    1
## [2,]    2    6

# t(A) %*% A
p1 <- t(A) %*% A
print(p1)

##      [,1] [,2]
## [1,]    8   14
## [2,]   14   37

# A %*% t(A)
p2 <- A %*% t(A)
print(p2)

##      [,1] [,2]
## [1,]    5   10
## [2,]   10   40

# The order matters
p1 == p2

##       [,1]  [,2]
## [1,] FALSE FALSE
## [2,] FALSE FALSE

# We can conclude that the assumption t(A) %*% A == A %*% t(A) is false and generalize.

# Commutative law: Attempt to multiply a matrix by its transposed matrix doesn't work unless for symmetric matrices.
# A %*% t(A) or t(A) %*% A to get symmetric matrix.
# Since t(S) == S, there is no difference in the dot prroduct as a result.

# create sample matrix
A <- matrix(sample(1:100, 16),4)
A

##      [,1] [,2] [,3] [,4]
## [1,]   29    7   70   61
## [2,]   64   56   57   14
## [3,]   98   24   30   87
## [4,]   55   74   50   49

# create symmetric matrix by multiplying its own transpose
S <- t(A) %*% A
S

##       [,1]  [,2]  [,3]  [,4]
## [1,] 17566 10209 11368 13886
## [2,] 10209  9237  8102  6925
## [3,] 11368  8102 11549 10128
## [4,] 13886  6925 10128 13887

# symmetric matrix multiply by its own transpose
p1 <- S %*% t(S)
p1

##           [,1]      [,2]      [,3]      [,4]
## [1,] 734840457 461895913 554330046 622588787
## [2,] 461895913 303143879 354600484 383952930
## [3,] 554330046 354600484 430829613 471578206
## [4,] 622588787 383952930 471578206 536201774

# transposed symmetric matrix multiplied by symmetric matrix
p2 <- t(S) %*% S
p2

##           [,1]      [,2]      [,3]      [,4]
## [1,] 734840457 461895913 554330046 622588787
## [2,] 461895913 303143879 354600484 383952930
## [3,] 554330046 354600484 430829613 471578206
## [4,] 622588787 383952930 471578206 536201774

# check whether the two are the same
p1 == p2

##      [,1] [,2] [,3] [,4]
## [1,] TRUE TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE TRUE
## [3,] TRUE TRUE TRUE TRUE
## [4,] TRUE TRUE TRUE TRUE

ALU <- function(A){
  # Factorize A, a square matrix, into Lower Triangular(L) and Upper Triangular(U) matrix
  
  # This function only works for square matrix, thus check whether the matrix is square first.
  if(nrow(A) != ncol(A)){
    stop("This is not a square matrix. Try again.")
  }
  
  # initialize variables
  counter <- as.integer(nrow(A))
  U <- A
  L <- diag(counter)
  # m for rows, n for columns
  
  # iterate through columns
  for (n in 1:(counter-1)){
    
    # iterate through the rows 
    m <- n+1
    
    for (m in (n+1):counter){
      x <- -U[m,n]/U[n,n] # U[m,n] <- U[m,n]+x*U[n,n]    
      U[m,n] <- 0
      
      # multiply the remiainings of current row by x
      n2 <- n+1
      
      for (n2 in (n+1):counter){
        U[m,n2] <- U[m,n2] + x*U[n,n2]
      }
      
    # assign the x to L
    L[m,n] = -x
    } 
  }
  
  print("A")
  print(A)
  print("==============================")
  print("U")
  print(U)
  print("==============================")
  print("L")
  print(L)
  print("==============================")
  print("A = LU")
  print(A == (L %*% U))
  print(L %*% U)
}

# Test
A <- matrix(sample(as.integer(1:100), 16),4)
ALU(A)

## [1] "A"
##      [,1] [,2] [,3] [,4]
## [1,]   43   86   52   35
## [2,]   65    9    5    8
## [3,]   51    3   41   56
## [4,]   97   27   58   96
## [1] "=============================="
## [1] "U"
##      [,1] [,2]      [,3]      [,4]
## [1,]   43   86  52.00000  35.00000
## [2,]    0 -121 -73.60465 -44.90698
## [3,]    0    0  39.54757  51.23044
## [4,]    0    0   0.00000  24.24997
## [1] "=============================="
## [1] "L"
##          [,1]      [,2]   [,3] [,4]
## [1,] 1.000000 0.0000000 0.0000    0
## [2,] 1.511628 1.0000000 0.0000    0
## [3,] 1.186047 0.8181818 1.0000    0
## [4,] 2.255814 1.3801653 1.0692    1
## [1] "=============================="
## [1] "A = LU"
##       [,1] [,2] [,3] [,4]
## [1,]  TRUE TRUE TRUE TRUE
## [2,]  TRUE TRUE TRUE TRUE
## [3,] FALSE TRUE TRUE TRUE
## [4,]  TRUE TRUE TRUE TRUE
##      [,1] [,2] [,3] [,4]
## [1,]   43   86   52   35
## [2,]   65    9    5    8
## [3,]   51    3   41   56
## [4,]   97   27   58   96