1. Write code in R to compute \(X = A{A^T}\) and \(Y = {A^T}A\).

#create a matrix
A = matrix(c(1,-1, 2,0, 3,4), nrow=2, byrow=FALSE)

#calculate transpose of matrix A
A.transpose = t(A)

#calculate X = matrix A dot product of transpose of matrix A
#results in 2X2 matrix
X = A%*%A.transpose

#calculate Y = transpose of matrix A dot product of matrix A
#results in 3X3 matrix
Y = A.transpose%*%A

Matrix \(A\) = \(\left[\begin{array}{cc}1 & 2 & 3 \\-1 & 0 & 4 \end{array}\right]\)

Transpose of matrix \({A^T}\) = \(\left[\begin{array}{cc}1 & -1 \\ 2 & 0 \\ 3 & 4\end{array}\right]\)

Dot product of matrix A and transpose of matrix A, \(X = A{A^T}\)

\(X\) = \(\left[\begin{array}{cc}14 & 11 \\11 & 17\end{array}\right]\)

Dot product of transpose of matrix A and matrix A, \(Y = {A^T}A\)

\(X\) = \(\left[\begin{array}{cc}2 & 2 & -1\\2 & 4 & 6 \\ -1 & 6 & 25\end{array}\right]\)

2. Compute the eigenvalues and eigenvectors of X and Y using the built-in commans in R.

The eigenvalues and eigenvectors will be calculated using R function eigen. Function eigen return two components values and vectors. The vectors are normalized to unit length.

#compute eigenvalues of X
X.eigenvalues = eigen(X)$values

#compute eigenvalues of X
X.eigenvectors = eigen(X)$vectors

Eigenvalues of \(X\) are \({\lambda}_1 = 26.6018017\) and \({\lambda}_2 = 4.3981983\)

Eigenvectors of \(X\) are

For \({\lambda}_1 = 26.6018017\), \({\cal E}_{X}(26.6018017) = \left[\begin{array}{cc}0.6576043 \\ 0.7533635\end{array}\right]\)

\({\lambda}_2 = 4.3981983\), \({\cal E}_{X}(4.3981983) = \left[\begin{array}{cc}-0.7533635 \\ 0.6576043\end{array}\right]\)

#compute eigenvalues of Y
Y.eigenvalues = eigen(Y)$values

#compute eigenvalues of Y
Y.eigenvectors = eigen(Y)$vectors

Eigenvalues of \(Y\) are \({\lambda}_1 = 26.6018017\), \({\lambda}_2 = 4.3981983\) and \({\lambda}_3 = -0.00000000000000060986\)

Eignevalue \({\lambda}_3\) has very small value, if rounded to 4 digits it is equal to zero. This value can be considered as mathematically significant, but practically non-significant.

Eigenvectors of \(Y\) are

For \({\lambda}_1 = 26.6018017\), \({\cal E}_{Y}(26.6018017) = \left[\begin{array}{cc}-0.0185663 \\ 0.2549994 \\ 0.966763 \end{array}\right]\)

\({\lambda}_2 = 4.3981983\), \({\cal E}_{Y}(4.3981983) = \left[\begin{array}{cc}0.6727903 \\ 0.718451 \\ -0.1765824\end{array}\right]\)

\({\lambda}_2 = -0.00000000000000060986\), \({\cal E}_{Y}(-0.00000000000000060986) = \left[\begin{array}{cc}0.7396003 \\ -0.6471502 \\ 0.1849001\end{array}\right]\)

3. Compute the left-singular, singular values, and right-singular vectors of A using the svd command

The SVD decomposition of the matrix is computed by: \(\mathbf{A = U \sum V^T}\).

Where \(\mathbf{U}\), \(\mathbf{V^T}\) are orthogonal matrices and \(\sum\) is diagonal matrix with the singular values.

\(\mathbf{U}\) = \(\left[\frac{1}{{\sigma}_1}AV_1 . \frac{1}{{\sigma}_2}AV_2 .... \frac{1}{{\sigma}_i}AV_i . \frac{{\cal N}(A^T)}{|{\cal N}(A^T)|}\right]\)

\({\sigma}_i = \sqrt{{\lambda}_i}\), eigenvectors \(\{{\lambda}_1, {\lambda}_2, ... {\lambda}_i\}\) of \({A^T}A\)

\(V_i =\) eigenvectors \(\{V_1, V_2, ... V_i\}\) of \({A^T}A\)

\({\cal N}(A^T)\) is NullSpace of transpose of matrix \(A\) and \(|{\cal N}(A^T)|\) is magnitude of NullSpace of transpose of matrix \(A\).

\(\mathbf{V^T}\) = transpose of eigenvectors of \({A^T}A\) = \({({A^T}A)}^T\)

\(\sum\) = diagonal matrix of square root of eigenvalues of \({A^T}A\)

Entire calculation can be computed by using R function svd.

#compute svd of matrix A
#dim(A)[1] is rows of matrix
#dim(A)[2] is cols of matrix

A.svd = svd(A, dim(A)[1], dim(A)[2])
#get u vectors
A.svd.u = A.svd$u
#get v vectors
A.svd.v = A.svd$v
#get diagonal matrix
A.svd.d = diag(A.svd$d)
A.svd.d.sq = diag(A.svd$d)%*%diag(A.svd$d)

Left-singular vectors of SVD are

\(\mathbf{U}\) = \(\left\{\left[\begin{array}{cc}-0.6576043 \\ -0.7533635 \end{array}\right] \left[\begin{array}{cc}-0.7533635 \\ 0.6576043 \end{array}\right] \right\}\)

When left-singular vectors of SVD \(\mathbf U\) are multiplied by \(-1\) scalar, they are same as eigenvectors of \(X\)

\(\left\{\left[\begin{array}{cc}0.6576043 \\ 0.7533635\end{array}\right] \left[\begin{array}{cc}-0.7533635 \\ 0.6576043\end{array}\right]\right\}\)

Right-singular vectors

\(\mathbf{V}\) = \(\left\{\left[\begin{array}{cc}0.0185663 \\ -0.2549994 \\ -0.966763\end{array}\right] \left[\begin{array}{cc}-0.6727903 \\ -0.718451\\0.1765824 \end{array}\right] \left[\begin{array}{cc}-0.7396003 \\ 0.6471502\\-0.1849001 \end{array}\right]\right\}\)

When right-singular vectors of SVD \(\mathbf V\) are multiplied by \(-1\) scalar, they are same as eigenvectors of \(Y\)

\(\left\{\left[\begin{array}{cc}-0.0185663 \\ 0.2549994 \\ 0.966763 \end{array}\right] \left[\begin{array}{cc}0.6727903 \\ 0.718451 \\ -0.1765824 \end{array}\right] \left[\begin{array}{cc}0.7396003 \\ -0.6471502 \\ 0.1849001 \end{array}\right]\right\}\)

Singular values

\(\sum\) = \(\left\{\left[\begin{array}{cc}5.1576934 \\ 0 \end{array}\right] \left[\begin{array}{cc}0 \\ 2.0971882 \end{array}\right] \right\}\)

4. Two non-zero eigenvalues (the 3rd value will be very close to zero, if not zero) of both \(X\) and \(Y\) are the same and are squares of the non-zero singular values of A

Square of singular values \(A\) are same as eigenvalues of \(X\) and \(Y\).

Eigenvalues of \(X\) are \(26.6018017\) and \(4.3981983\)

Eigenvalues of \(Y\) are \(26.6018017\) and \(4.3981983\), third one is left out because it is close to zero.

Square of singular values \(A\) are \(26.6018017\) and \(4.3981983\)

Using above \(\mathbf{U}\), \(\mathbf{V}\) and \(\sum\) we cannot arrive at orginal matrix \(A\), because \(\mathbf{V}\) is 3X3 matrix and other matrices are 2X2. In this case, we can ignore third vector of \(\mathbf{V}\) as it is near close to zero eigenvalue.

#if rows and cols dimensions are not passed svd will result in non-zero left-singular and right-singular vectors.
A.svd = svd(A)
#get u vectors
A.svd.u = A.svd$u
#get v vectors
A.svd.v = A.svd$v
#get diagonal matrix
A.svd.d = diag(A.svd$d)

#perform dot product of u, d and transpose of v to get orginal matrix
round(A.svd.u%*%A.svd.d %*%t(A.svd.v),0)
##      [,1] [,2] [,3]
## [1,]    1    2    3
## [2,]   -1    0    4
A==round(A.svd.u%*%A.svd.d %*%t(A.svd.v),0)
##      [,1] [,2] [,3]
## [1,] TRUE TRUE TRUE
## [2,] TRUE TRUE TRUE

Following is a custom function to calculate inverse of matrix. Function returns messages if given matrix is not a square matrix or its determinant is zero.

myinverse <- function(A){
  mRows <- nrow(A)
  mCols <- ncol(A)

  if (mRows != mCols) {
    #if rows != cols such matrix does not have inverse
    print ("Not a square matrix. Inverse cannot be calculated ...")
    return(-1)
  }
  
  if (round(det(A),0) == 0) {
    #if determinant is zero inverse cannot be calculated. Because 1 divide by determinant is infinity. 
    print ("This is a singular matrix as determinant is 0. Inverse cannot be calculated ...")
    return(-1)
  }
  
  #once matrix satisfies above conditions, initiate return matrix
  retM <- matrix(rep(0,length(A)), nrow = mRows, ncol = mCols)
  
  #loop columns and rows
  for(i in 1:mCols){
    for (j in 1:mRows){
      #calculate cofactors
      retM[j,i] = det(A[-j,-i])
      #change sign for alternative cell
      if ((j+i) %% 2 != 0){
        retM[j,i] = -1 * retM[j,i]
      }
    }
  }
  
  #apply tranpose and divide by determinant
  #resultant matrix is inverse matrix
  retM = t(retM)/det(A)

  return(retM)
}

Function test cases

#given matrix is A is non square matrix
A = matrix(c(1,-1, 2,0, 3,4), nrow=2, byrow=FALSE)
#calculate inverse using custom function
B = myinverse(A) 
## [1] "Not a square matrix. Inverse cannot be calculated ..."
#given matrix is A is a square matrix, with determinant as zero
A = matrix(c(1,3,8, 2,0,4, 1,4,10), nrow=3, byrow=FALSE)
det(A)
## [1] 0
#calculate inverse using custom function
B = myinverse(A) 
## [1] "This is a singular matrix as determinant is 0. Inverse cannot be calculated ..."
#given matrix is A is square matrix
A = matrix(c(1,3,2,1, 0,0,1,0, 2,0,4,5, -1,5,-3,0), nrow=4, byrow=FALSE)
#calculate inverse using custom function
B = myinverse(A)
#inverse matrix of A
B
##            [,1]        [,2] [,3]       [,4]
## [1,]  0.8333333  0.16666667    0 -0.3333333
## [2,] -2.5000000  0.10000000    1  0.2000000
## [3,] -0.1666667 -0.03333333    0  0.2666667
## [4,] -0.5000000  0.10000000    0  0.2000000
#dot product of matrix A and B should be equal to identity matrix
round(A%*%B,0)
##      [,1] [,2] [,3] [,4]
## [1,]    1    0    0    0
## [2,]    0    1    0    0
## [3,]    0    0    1    0
## [4,]    0    0    0    1
#given matrix is A is a square matrix
A = matrix(c(1,2,4, 2,3,9, 1,3,6), nrow=3, byrow=FALSE)
#calculate inverse using custom function
B = myinverse(A)
#inverse matrix of A
B
##      [,1]       [,2]       [,3]
## [1,]    3  1.0000000 -1.0000000
## [2,]    0 -0.6666667  0.3333333
## [3,]   -2  0.3333333  0.3333333
#dot product of matrix A and B should be equal to identity matrix
round(A%*%B,0)
##      [,1] [,2] [,3]
## [1,]    1    0    0
## [2,]    0    1    0
## [3,]    0    0    1