VStoyanova_Assign4

Problem Set 1

A <- matrix(c(1, 2, 3, -1, 0, 4), nrow = 2, byrow = TRUE)
# compute X & Y
X <- A %*% t(A)
X

##      [,1] [,2]
## [1,]   14   11
## [2,]   11   17

Y <- t(A) %*% A
Y

##      [,1] [,2] [,3]
## [1,]    2    2   -1
## [2,]    2    4    6
## [3,]   -1    6   25

# eigenvalues
lambda_x <- eigen(X)$values
lambda_x

## [1] 26.601802  4.398198

lambda_y <- eigen(Y)$values
lambda_y

## [1] 2.660180e+01 4.398198e+00 1.058982e-16

# eigenvectors
s_x <- eigen(X)$vectors
s_x

##           [,1]       [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635  0.6576043

s_y <- eigen(Y)$vectors
s_y

##             [,1]       [,2]       [,3]
## [1,] -0.01856629 -0.6727903  0.7396003
## [2,]  0.25499937 -0.7184510 -0.6471502
## [3,]  0.96676296  0.1765824  0.1849001

# performing SVD of A
svd_a <- svd(A)
svd_a

## $d
## [1] 5.157693 2.097188
## 
## $u
##            [,1]       [,2]
## [1,] -0.6576043 -0.7533635
## [2,] -0.7533635  0.6576043
## 
## $v
##             [,1]       [,2]
## [1,]  0.01856629 -0.6727903
## [2,] -0.25499937 -0.7184510
## [3,] -0.96676296  0.1765824

u contains the eigenvalues for X and v contains the first two eigenvectors of Y. The first columns, corresponding to the first eigenvalues of X and Y are shown in the opposite direction in U and V. These vectors are equivalent, as they simply represent scalar multiplication and do not affect orthonormality.

Problem Set 2

myinverse <- function(A) {
  # check if matrix is full-rank and square
  if (nrow(A) == ncol(A) & Matrix::rankMatrix(A)[1] == nrow(A)) {
    fr_square <- TRUE
  } else {fr_square <- FALSE}
  
  if (fr_square == TRUE) {
    size <- nrow(A)
    C <- matrix(nrow = size, ncol = size)
    for (i in 1:size) {
      for (j in 1:size) {
        # M is A with row i & column j excluded
        M <- A[-i, -j]
        C[i, j] <- (-1)^(i + j) * det(M)
      }
    }
    # return inverse of T divided by det(A)
    B <- t(C) / det(A)
  } else {B <- "Matrix is not invertible"}
  B
}

Testing the function

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

##      [,1] [,2] [,3]
## [1,]    1    3    2
## [2,]    1    0    1
## [3,]    5    1    3

B<-myinverse(A)
B

##            [,1] [,2]       [,3]
## [1,] -0.1428571   -1  0.4285714
## [2,]  0.2857143   -1  0.1428571
## [3,]  0.1428571    2 -0.4285714

round(A %*% B, 14)

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

The returned results are tested to show that A×B=I and that the matrix calculated with myinverse is equal to the inverse returned by the solve function.