library('matrixcalc')
library('pracma')

1

A = cbind(c(1,-1), c(2,0), c(3,4))
A
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]   -1    0    4
# AAt and AtA
X = A %*% t(A)
print('X')
## [1] "X"
X
##      [,1] [,2]
## [1,]   14   11
## [2,]   11   17
Y = t(A) %*% A
print('Y')
## [1] "Y"
Y
##      [,1] [,2] [,3]
## [1,]    2    2   -1
## [2,]    2    4    6
## [3,]   -1    6   25
# eigens of X and Y
X_evals_evecs = eigen(X)
print('X eigenvals and eigenvecs')
## [1] "X eigenvals and eigenvecs"
X_evals_evecs
## eigen() decomposition
## $values
## [1] 26.601802  4.398198
## 
## $vectors
##           [,1]       [,2]
## [1,] 0.6576043 -0.7533635
## [2,] 0.7533635  0.6576043
Y_evals_evecs = eigen(Y)
print('Y eigenvals and eigenvecs')
## [1] "Y eigenvals and eigenvecs"
Y_evals_evecs
## eigen() decomposition
## $values
## [1] 2.660180e+01 4.398198e+00 1.058982e-16
## 
## $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
# left singular vectors, singular values, right singular vectors
svd_A = svd(A)
print('SVD of A')
## [1] "SVD of 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
# check left singular vectors. Check out (first eigenvector matches, just flip sign)
print('Checking left singular vectors. First eigenvector scale by -1 thus match')
## [1] "Checking left singular vectors. First eigenvector scale by -1 thus match"
print(round(svd_A$u, 3))
##        [,1]   [,2]
## [1,] -0.658 -0.753
## [2,] -0.753  0.658
print(round(X_evals_evecs$vectors,3))
##       [,1]   [,2]
## [1,] 0.658 -0.753
## [2,] 0.753  0.658
# check singular values
# singular value is defined as sqrt of eigenvalues
# first verify that two non-zero eigenvalues of AAt and AtA match
print('Verifying that nonzero eigenvalues from X and Y match')
## [1] "Verifying that nonzero eigenvalues from X and Y match"
print(round(X_evals_evecs$values, 3))
## [1] 26.602  4.398
print(round(Y_evals_evecs$values, 3))
## [1] 26.602  4.398  0.000
# next verify the sqrt matches the diagonal d
print('Verifying that singular values are sqrt of eigenvalues')
## [1] "Verifying that singular values are sqrt of eigenvalues"
print(round(sqrt(X_evals_evecs$values), 3))
## [1] 5.158 2.097
print(round(svd_A$d, 3))
## [1] 5.158 2.097
print('Checking right singular vectors')
## [1] "Checking right singular vectors"
print(round(svd_A$v, 3))
##        [,1]   [,2]
## [1,]  0.019 -0.673
## [2,] -0.255 -0.718
## [3,] -0.967  0.177
print(round(Y_evals_evecs$vectors,3))
##        [,1]   [,2]   [,3]
## [1,] -0.019 -0.673  0.740
## [2,]  0.255 -0.718 -0.647
## [3,]  0.967  0.177  0.185
print('verify the last eigenvector corresponding to the 0 eigenvalue. Can see it is in null space of Y')
## [1] "verify the last eigenvector corresponding to the 0 eigenvalue. Can see it is in null space of Y"
# make augemented matrix and rref
rref(cbind(A, c(0,0)))
##      [,1] [,2] [,3] [,4]
## [1,]    1    0 -4.0    0
## [2,]    0    1  3.5    0
print(-3.5 * .185)
## [1] -0.6475
print(-4 * .185)
## [1] -0.74
print('verify svd function')
## [1] "verify svd function"
round(svd_A$u %*% diag(svd_A$d) %*% t(svd_A$v),3)
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]   -1    0    4
print('verify svd with parts of X/Y eigenvals and eigenvecs')
## [1] "verify svd with parts of X/Y eigenvals and eigenvecs"
round(X_evals_evecs$vectors %*% cbind(diag(svd_A$d),c(0,0)) %*% t(Y_evals_evecs$vectors), 3)
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]   -1    0    4

2

myinverse = function(A) {
  # find inverse using cofactors
  # cofactor = (sign) * det(submatrix)
  # A-inv = Ct/det(A)
  # Need Ct and det(A)
  det_A = det(A)
  # initialize C
  C = matrix(data=NA, nrow = nrow(A), ncol=ncol(A))
  for (i in 1:nrow(A)) {
    for (j in 1:ncol(A)) {
      # always use first row
      sign = (-1)**(i+j)
      sub_m = A[-i, -j]
      C[i,j] = sign * det(sub_m)
    }
  }
  return(t(C)/det_A)
}

Test some examples

A = cbind(c(1,2,-1), c(1,-1,-2), c(3,5,4))
A
##      [,1] [,2] [,3]
## [1,]    1    1    3
## [2,]    2   -1    5
## [3,]   -1   -2    4
round(A %*% myinverse(A), 3)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1
B = cbind(c(1,0, 2,1), c(1,3,3, 0), c(1,1,1,2), c(0,2,0,1))
B
##      [,1] [,2] [,3] [,4]
## [1,]    1    1    1    0
## [2,]    0    3    1    2
## [3,]    2    3    1    0
## [4,]    1    0    2    1
round(B %*% myinverse(B), 3)
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    0    1    0    0
## [3,]    0    0    1    0
## [4,]    0    0    0    1