Singular Value Decomposition of a Matrix

Compute the singular-value decomposition of a matrix \(X\) either by Jacobi rotations (the default) or from the eigenstructure of \(X'X\) using Eigen. Both methods are iterative. The result consists of two orthonormal matrices, \(U\), and \(V\) and the vector \(d\) of singular values, such that \(X = U diag(d) V'\).

Usage

SVD(
  X,
  method = c("Jacobi", "eigen"),
  tol = sqrt(.Machine$double.eps),
  max.iter = 100
)

Arguments

X: a square symmetric matrix
method: either "Jacobi" (the default) or "eigen"
tol: zero and convergence tolerance
max.iter: maximum number of iterations

Value

a list of three elements: d– singular values, U– left singular vectors, V– right singular vectors

Details

The default method is more numerically stable, but the eigenstructure method is much simpler. Singular values of zero are not retained in the solution.

Author

John Fox and Georges Monette

Examples

C <- matrix(c(1,2,3,2,5,6,3,6,10), 3, 3) # nonsingular, symmetric
C
#>      [,1] [,2] [,3]
#> [1,]    1    2    3
#> [2,]    2    5    6
#> [3,]    3    6   10
SVD(C)
#> $d
#> [1] 14.93303437  1.00000000  0.06696563
#> 
#> $U
#>           [,1]          [,2]       [,3]
#> [1,] 0.2505248  4.308756e-14  0.9681102
#> [2,] 0.5370109  8.320503e-01 -0.1389662
#> [3,] 0.8055164 -5.547002e-01 -0.2084492
#> 
#> $V
#>           [,1]          [,2]       [,3]
#> [1,] 0.2505248  2.960306e-15  0.9681102
#> [2,] 0.5370109  8.320503e-01 -0.1389662
#> [3,] 0.8055164 -5.547002e-01 -0.2084492
#> 

# least squares by the SVD
data("workers")
X <- cbind(1, as.matrix(workers[, c("Experience", "Skill")]))
head(X)
#>           Experience Skill
#> Abby    1          0     2
#> Betty   1          5     5
#> Charles 1          5     8
#> Doreen  1         10     6
#> Ethan   1         10    10
#> Francie 1         15     7
y <- workers$Income
head(y)
#> [1] 20 35 40 30 50 50
(svd <- SVD(X))
#> $d
#> [1] 66.864276  8.885751  1.100924
#> 
#> $U
#>         Experience       Skill            
#> Abby    0.01366778  0.21621629  0.62512040
#> Betty   0.10046906  0.26924846  0.34867115
#> Charles 0.12007689  0.56866664 -0.06013463
#> Doreen  0.17419846  0.12266851  0.34475908
#> Ethan   0.20034223  0.52189275 -0.20031530
#> Francie 0.24792785 -0.02391144  0.34084702
#> Georges 0.34780102  0.22873284 -0.20813943
#> Harry   0.40192259 -0.21726529  0.19675429
#> Isaac   0.50179576  0.03537900 -0.35223216
#> Juan    0.55591732 -0.41061914  0.05266156
#> 
#> $V
#>            [,1]       [,2]        [,3]
#> [1,] 0.03984368  0.1475406  0.98825314
#> [2,] 0.89856830 -0.4378649  0.02914291
#> [3,] 0.43702116  0.8868518 -0.15002142
#> 
VdU <- svd$V %*% diag(1/svd$d) %*%t(svd$U)
(b <- VdU %*% y)
#>            [,1]
#> [1,] 14.7808811
#> [2,]  0.1130069
#> [3,]  3.4053991
coef(lm(Income ~ Experience + Skill, data=workers))
#> (Intercept)  Experience       Skill 
#>  14.7808811   0.1130069   3.4053991