R version 2.9.0 (2009-04-17)
Copyright (C) 2009 The R Foundation for Statistical Computing
ISBN 3-900051-07-0

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

  Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

>  # svd1.r                 September 2007, 2009
>  #
>  # some examples of use of svd to do Singular Value Decomposition
>  #
>  # first example -- singular design matrix
>   X <- matrix(1,6,4)
>   X[,2] <- 1:6
>   X[,3] <- 19:24    # dependent on first two
>   X[,4] <- (1:6)*(1:6)
>   X                 # should look familiar
     [,1] [,2] [,3] [,4]
[1,]    1    1   19    1
[2,]    1    2   20    4
[3,]    1    3   21    9
[4,]    1    4   22   16
[5,]    1    5   23   25
[6,]    1    6   24   36
>  #                  first qr
>   qr(X)             # what does QR do in this case?
$qr
           [,1]       [,2]        [,3]          [,4]
[1,] -2.4494897 -8.5732141 -37.1505944 -5.266403e+01
[2,]  0.4082483  4.1833001  29.2831009  4.183300e+00
[3,]  0.4082483 -0.0537243   6.1101009 -1.134008e-15
[4,]  0.4082483 -0.2927700   0.6606006  3.087130e-16
[5,]  0.4082483 -0.5318157   0.3871728  1.182847e-01
[6,]  0.4082483 -0.7708615  -0.2135819 -8.588020e-01

$rank
[1] 3

$qraux
[1] 1.408248 1.185321 1.606702 1.498465

$pivot
[1] 1 2 4 3

attr(,"class")
[1] "qr"
>  #                  next svd
>   svd(X)
$d
[1] 6.872233e+01 2.096727e+01 7.845345e-01 2.236524e-16

$u
           [,1]       [,2]       [,3]        [,4]
[1,] -0.2170374 -0.5668386  0.6730656 -0.42017099
[2,] -0.2584600 -0.4899020 -0.0183671  0.79210474
[3,] -0.3189691 -0.3419531 -0.3908393 -0.07786869
[4,] -0.3985647 -0.1229921 -0.4443508 -0.31731245
[5,] -0.4972467  0.1669810 -0.1789019 -0.19933304
[6,] -0.6150152  0.5279663  0.4055077  0.22258042

$v
            [,1]        [,2]        [,3]          [,4]
[1,] -0.03354504 -0.03942996  0.05877906  9.969278e-01
[2,] -0.13767649  0.04474746 -0.98791465  5.538488e-02
[3,] -0.74148723 -0.66499174  0.07010847 -5.538488e-02
[4,] -0.65583276  0.74446554  0.12511778  5.204170e-17

>  #                  inner product matrix
>   xpx <- t(X) %*% X
>   print("X'X")
[1] "X'X"
>   xpx
     [,1] [,2] [,3] [,4]
[1,]    6   21  129   91
[2,]   21   91  469  441
[3,]  129  469 2791 2079
[4,]   91  441 2079 2275
>   eigen(xpx,symmetric=TRUE)            # eigenvectors and values of X'X
$values
[1] 4.722758e+03 4.396265e+02 6.154944e-01 3.039709e-15

$vectors
            [,1]        [,2]        [,3]          [,4]
[1,] -0.03354504 -0.03942996  0.05877906  9.969278e-01
[2,] -0.13767649  0.04474746 -0.98791465  5.538488e-02
[3,] -0.74148723 -0.66499174  0.07010847 -5.538488e-02
[4,] -0.65583276  0.74446554  0.12511778  1.425769e-14

>   xxp <- X %*% t(X)
>   print("XX'")
[1] "XX'"
>   xxp 
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]  364  387  412  439  468  499
[2,]  387  421  463  513  571  637
[3,]  412  463  532  619  724  847
[4,]  439  513  619  757  927 1129
[5,]  468  571  724  927 1180 1483
[6,]  499  637  847 1129 1483 1909
>   exxp <- eigen(xxp,symmetric=TRUE)    # vector and values of XX'
>   exxp$values                          # eigenvalues are same except 0's
[1]  4.722758e+03  4.396265e+02  6.154944e-01  4.188385e-14 -1.389260e-13
[6] -1.923924e-13
>   sqrt(exxp$values)                    # singular values
[1] 6.872233e+01 2.096727e+01 7.845345e-01 2.046554e-07          NaN
[6]          NaN
>   exxp$vectors                         # other singular vectors
           [,1]       [,2]       [,3]         [,4]          [,5]       [,6]
[1,] -0.2170374  0.5668386  0.6730656  0.000000000 -4.225771e-01  0.0000000
[2,] -0.2584600  0.4899020 -0.0183671 -0.286920311  7.606388e-01 -0.1787805
[3,] -0.3189691  0.3419531 -0.3908393  0.763822852 -2.231548e-13  0.2124761
[4,] -0.3985647  0.1229921 -0.4443508 -0.569946692 -3.380617e-01  0.4352545
[5,] -0.4972467 -0.1669810 -0.1789019 -0.003893929 -2.535463e-01 -0.7928154
[6,] -0.6150152 -0.5279663  0.4055077  0.096938080  2.535463e-01  0.3238653
>  # done
>   rm(list=ls())
>  q()