Problem Set 1

\[A = \left[ \begin{array}{ccc}1 & 2 & 3 \\-1 & 0 & 4 \end{array} \right]\]

A <- matrix(c(1, 2, 3, -1, 0, 4), nrow = 2, byrow = TRUE)
# compute X & Y
(X <- A %*% t(A))
     [,1] [,2]
[1,]   14   11
[2,]   11   17
(Y <- t(A) %*% A)
     [,1] [,2] [,3]
[1,]    2    2   -1
[2,]    2    4    6
[3,]   -1    6   25
# eigenvalues
(lambda_x <- eigen(X)$values)
[1] 26.601802  4.398198
(lambda_y <- eigen(Y)$values)
[1] 2.660180e+01 4.398198e+00 1.058982e-16
# eigenvectors
(s_x <- eigen(X)$vectors)
          [,1]       [,2]
[1,] 0.6576043 -0.7533635
[2,] 0.7533635  0.6576043
(s_y <- eigen(Y)$vectors)
            [,1]       [,2]       [,3]
[1,] -0.01856629 -0.6727903  0.7396003
[2,]  0.25499937 -0.7184510 -0.6471502
[3,]  0.96676296  0.1765824  0.1849001
# perform SVD of 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

Comparing the results, it can be seen that \(U\) contains the eigenvalues of \(X\), that is \[U = S_X = \left[ \begin{array}{cc}-0.6576 & -0.7534 \\-0.7534 & 0.6576 \end{array} \right]\]

Additionally, \(V\) contains the first two eigenvectors of \(Y\).
\[V = S_{Y_{1:2}} = \left[ \begin{array}{ccc}0.0186 & -0.6728 \\-0.2550 & -0.7185 \\-0.9668 & 0.1766 \end{array} \right]\]

In both cases above, 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.

Finally, the non-zero eigenvalues of \(V\) are equivalent to the eigenvalues of \(U\) – these are equivalent to the square of the singular values \(d\): \[\Lambda_X = \Lambda_{Y_{1:2}} = \left[ \begin{array}{cc}26.6018 & 0 \\0 & 4.3982 \end{array} \right] = \Sigma^2\] Where \[\Sigma = \left[ \begin{array}{cc}5.1577 & 0 \\0 & 2.0972 \end{array} \right]\]

Problem Set 2

The function below returns the inverse of a full-rank square matrix using co-factors.

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

Testing with the example matrix \(A\) from Assignment 2

\[A = \left[ \begin{array}{cc}1 & 1 & 2 \\2 & 1 & 0 \\3 & 1 & 1 \end{array} \right]\]

A <- matrix(c(1, 2, 3, 1, 1, 1, 2, 0, 1), nrow=3)
(B <- myinverse(A))
           [,1]       [,2]       [,3]
[1,] -0.3333333 -0.3333333  0.6666667
[2,]  0.6666667  1.6666667 -1.3333333
[3,]  0.3333333 -0.6666667  0.3333333

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

round(A %*% B, 14)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    1    0
[3,]    0    0    1
identical(round(B, 14), round(solve(A), 14))
[1] TRUE