DATA 605: Assignment 4

Problem Set 1

(1) Given a 3 × 2 matrix A

\[ A = \begin{pmatrix} 1 & 2 & 3 \\ -1 & 0 & 4 \\ \end{pmatrix} \]

Write code in R to compute \(X = A A^T\) and \(Y = A^TA\). Then, compute the eigenvalues and eigenvectors of X and Y using the built-in commands in R.

A = matrix(c(1, 2, 3, -1, 0, 4), nrow=2, ncol=3, byrow=T)

\(X = A A^T\) is calculated in R by:

X = A %*% t(A)

\(Y = A^TA\) is calculated in R by:

Y = t(A) %*% A

The eigenvalues and eigenvalues of \(X\) are calculated using the eigen function:

x_eigen <- eigen(X)

x_eigenvalues  <- x_eigen$values
x_eigenvectors <- x_eigen$vectors

The eigenvalues of \(X\) are:

\(\lambda_1\) = 26.6018017
\(\lambda_2\) = 4.3981983

and the eigenvalues of \(X\) are:

\[ v_1 = \begin{pmatrix} 0.6576043 \\ 0.7533635 \\ \end{pmatrix} \]

\[ v_2 = \begin{pmatrix} -0.7533635 \\ 0.6576043 \\ \end{pmatrix} \]

Similarly the The eigenvalues and eigenvalues of \(Y\) are calculated as:

y_eigen <- eigen(Y)

y_eigenvalues <- y_eigen$values
y_eigenvectors <- y_eigen$vectors

The eigenvalues of \(Y\) are:

\(\lambda_1\) = 26.6018017
\(\lambda_2\) = 4.3981983
\(\lambda_3\) = 1.05898210^{-16}

and the eigenvalues of \(Y\) are:

\[ v_1 = \begin{pmatrix} -0.0185663 \\ 0.2549994 \\ 0.966763 \\ \end{pmatrix} \]

\[ v_2 = \begin{pmatrix} -0.6727903 \\ -0.718451 \\ 0.1765824 \\ \end{pmatrix} \]

\[ v_3 = \begin{pmatrix} 0.7396003 \\ -0.6471502 \\ 0.1849001 \\ \end{pmatrix} \]

Then, compute the left-singular, singular values, and right-singular vectors of A using the svd command. Examine the two sets of singular vectors and show that they are indeed eigenvectors of X and Y. In addition, the 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. Your code should compute all these vectors and scalars and store them in variables.

The SVD decomposition of the matrix is computed by:

\[\mathbf{X = U D V} \]

where \(\mathbf{U}\) and \(\mathbf{V}\) are orthogonal, \(\mathbf{V}\) means \(V\) transposed, and \(\mathbf{D}\) is a diagonal matrix with the singular values \(D[i,i]\).

The returned value is a list with components

d = a vector containing the singular values of x, of length min(n, p).
u = a matrix whose columns contain the left singular vectors of x, present if nu > 0. Dimension c(n, nu).
v = a matrix whose columns contain the right singular vectors of x, present if nv > 0. Dimension c(p, nv).

Calculate the left-singular, singular values, and right-singular vectors of \(A\) using the svd command in R:

svdA <- svd(A, dim(A)[1], dim(A)[2])

The left-singular vectors (u) of \(A\) =

svdA$u

\[ u_1 = \begin{pmatrix} -0.6576043 \\ -0.7533635 \\ \end{pmatrix} \]

\[ u_2 = \begin{pmatrix} -0.7533635 \\ 0.6576043 \\ \end{pmatrix} \]

The singular values (d) of \(A\) = 5.1576934 and 2.0971882

The right-singular vectors (v) of \(A\) =

svdA$v

\[ v_1 = \begin{pmatrix} 0.0185663 \\ -0.2549994 \\ -0.966763 \\ \end{pmatrix} \]

\[ v_2 = \begin{pmatrix} -0.6727903 \\ -0.718451 \\ 0.1765824 \\ \end{pmatrix} \]

\[ v_3 = \begin{pmatrix} -0.7396003 \\ 0.6471502 \\ -0.1849001 \\ \end{pmatrix} \]

Examine the two sets of singular vectors and show that they are indeed eigenvectors of X and Y

Looking at the eigenvectors of \(X\) and the left-singular values, we see that the only difference between them is that the left-singular \(v_1\) appears to be multipled by a scalar value of -1. We know that a eigenvector multiplied by a scalar is in the same place so we can conclude that these are indeed equal.

Eigenvectors of \(X\):

\[ v_1 = \begin{pmatrix} 0.6576043 \\ 0.7533635 \\ \end{pmatrix} \]

\[ v_2 = \begin{pmatrix} -0.7533635 \\ 0.6576043 \\ \end{pmatrix} \]

Left-singular values

\[ u_1 = \begin{pmatrix} -0.6576043 \\ -0.7533635 \\ \end{pmatrix} \]

\[ u_2 = \begin{pmatrix} -0.7533635 \\ 0.6576043 \\ \end{pmatrix} \]

Looking at the eigenvectors of \(Y\) and the right-singluar values, we see that the only difference again is a scalar multiple.

The eigenvalues of \(Y\) are:

\[ v_1 = \begin{pmatrix} -0.0185663 \\ 0.2549994 \\ 0.966763 \\ \end{pmatrix} \]

\[ v_2 = \begin{pmatrix} -0.6727903 \\ -0.718451 \\ 0.1765824 \\ \end{pmatrix} \]

\[ v_3 = \begin{pmatrix} 0.7396003 \\ -0.6471502 \\ 0.1849001 \\ \end{pmatrix} \]

The right-singular vectors (v) of \(A\) =

\[ v_1 = \begin{pmatrix} 0.0185663 \\ -0.2549994 \\ -0.966763 \\ \end{pmatrix} \]

\[ v_2 = \begin{pmatrix} -0.6727903 \\ -0.718451 \\ 0.1765824 \\ \end{pmatrix} \]

\[ v_3 = \begin{pmatrix} -0.7396003 \\ 0.6471502 \\ -0.1849001 \\ \end{pmatrix} \]

Lastly, the 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.

This can be seen inspecting the eigenvalues of \(X\), \(Y\), and comparing the the squared singular values:

Eignenvalues of \(X\): 26.6018017 and 4.3981983
Eignenvalues of \(Y\): 26.6018017 and 4.3981983
Squared singular values: 26.6018017 and 4.3981983

Problem Set 2

Using the procedure outlined in section 1 of the weekly handout, write a function to compute the inverse of a well-conditioned full-rank square matrix using co-factors. In order to compute the co-factors, you may use built-in commands to compute the determinant. Your function should have the following signature:

B = myinverse(A)

where \(A\) is a matrix and \(B\) is its inverse and \(A×B = I\). The off-diagonal elements of \(I\) should be close to zero, if not zero. Likewise, the diagonal elements should be close to 1, if not 1. Small numerical precision errors are acceptable but the function myinverse should be correct and must use co-factors and determinant of \(A\) to compute the inverse.

myinverse <- function(A) {
    
    
  if (det(A) == 0) {stop("The determinant of the matrix provided is 0. Processing stopped...")}
    
  rows <- nrow(A)
  cols <- ncol(A)
  

  # identify if this is a 2xm matrix
  # if so, this will be used later in the matrix of minors calculation
  
  m2x2 = FALSE
  if ((rows == 2) & (cols == 2))  { m2x2 = TRUE }
 
  #initialize matrices for the matrix of minors and cofactors
  matrixMinors <- matrix(0L, rows, cols)
  matrixCofactors <- matrix(0L, rows, cols)
  
  ## Step 1:
  # --------------------------------------------
  #        A - calculate matrix of minors
  #        B - caluate the matrix of cofactors
  #
  #  2x2 matrices are handled specially since you cannot take the determinant

  
  for (i in 1:nrow(A)) {
      
      for (j in 1:ncol(A)) {
          
          if (m2x2) {
              matrixMinors[i, j] <- A[-i, -j]
          }
          else {
             matrixMinors[i, j] <- det(A[-i, -j])
          }
          
          matrixCofactors[i, j] <- matrixMinors[i, j] * (-1)^(i+j)
          
      }
  }
  
  ## Step 2: 
  ## ----------------------------
  ## create the adjugate of the provided matrix A
  ## This is the transpose of the cofactor matrix
  
  adj_A <- t(matrixCofactors)
  
  ## Step 3: calculate the determinant of the provided matrix A
  ## 
  ## This is simplified now that we have the matrix of minors calculated
  
  
  # This code calculates the determinant using the matrix A and the calculated matrix of minors
      # det <- 0
      # 
      # for (i in 1:ncol(A)) {
      #     
      #     det <- det + A[1, i] * (matrixMinors[1, i] * (-1)^(i-1))
      #     
      # }
      # 
  
  ## Final Step:
  ## ----------------------------
  ## Calculate the inverse of A using
  ## A inverse = 1/det(A) * adjoint(A)
  
  invA <- 1/det(A) * adj_A
  
}

Test the myinverse function

2x2 Matrix Test:

# 2x2 Matrix
A = matrix(c(2, 1, 1, 1), nrow=2, ncol=2, byrow=T)

B <- myinverse(A)

The inverse of \(A\):

##      [,1] [,2]
## [1,]    1   -1
## [2,]   -1    2

Confirm that \(AB = I\). The identity matrix looks good for the 2x2 matrix test.

B%*%A

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

3x3 Matrix Test:

# 3x3 Matrix
A = matrix(c(2, 1, 1, 1, 2, 3, 1, -1, 1), nrow=3, ncol=3, byrow=T)

B <- myinverse(A)

The inverse of \(A\):

##            [,1]       [,2]       [,3]
## [1,]  0.5555556 -0.2222222  0.1111111
## [2,]  0.2222222  0.1111111 -0.5555556
## [3,] -0.3333333  0.3333333  0.3333333

Confirm that \(AB = I\).

B%*%A

##               [,1]          [,2]          [,3]
## [1,]  1.000000e+00  1.387779e-16  1.942890e-16
## [2,] -1.110223e-16  1.000000e+00 -2.220446e-16
## [3,] -1.110223e-16 -5.551115e-17  1.000000e+00

We see some small numerical precision errors which can be accounted for:

round(B%*%A, 6)

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

4x4 Matrix Test:

# 4x4 Matrix
A = matrix(c(1, 2, 1, 0, 2, 1, 1, 1, -1, 2, 1, -1, 1, 1, 1, 2), nrow=4, ncol=4, byrow=T)

B <- myinverse(A)

The inverse of \(A\):

##            [,1]       [,2]       [,3]       [,4]
## [1,]  0.3333333  0.3333333 -0.3333333 -0.3333333
## [2,]  1.6666667 -1.3333333 -0.6666667  0.3333333
## [3,] -2.6666667  2.3333333  1.6666667 -0.3333333
## [4,]  0.3333333 -0.6666667 -0.3333333  0.6666667

Confirm that \(AB = I\).

round(B%*%A, 6)

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

5x5 Matrix Test:

# 5x5 Matrix
A = matrix(c(1, 2, 1, 0, 2, 1, 1, 1, -1, 2, 1, -1, 1, 1, 1, 2, 1, 3, 2, 4, 3, 2, 4, 3, 2), nrow=5, ncol=5, byrow=T)

B <- myinverse(A)

The inverse of \(A\):

##            [,1]       [,2]       [,3]        [,4]        [,5]
## [1,]  1.4285714 -0.2142857  1.7142857 -1.07142857  0.07142857
## [2,]  0.4285714 -0.2142857 -0.2857143 -0.07142857  0.07142857
## [3,] -1.8571429  0.9285714 -1.4285714  0.64285714  0.35714286
## [4,]  0.5714286 -0.7857143  0.2857143  0.07142857 -0.07142857
## [5,]  0.2857143 -0.1428571  0.1428571  0.28571429 -0.28571429

Confirm that \(AB = I\).

round(B%*%A, 6)

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

DATA 605: Assignment 4

Keith Folsom

September 16, 2016

Problem Set 1

(1) Given a 3 × 2 matrix A

Examine the two sets of singular vectors and show that they are indeed eigenvectors of X and Y

Lastly, the 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.

Problem Set 2

Test the myinverse function