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'\).

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.

See also

svd, the standard svd function

Eigen

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