R version 3.0.0 RC (2013-03-29 r62438) -- "Masked Marvel"
Copyright (C) 2013 The R Foundation for Statistical Computing
Platform: x86_64-unknown-linux-gnu (64-bit)
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.
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.
> ## tests of R functions based on the lapack module
>
> options(digits=4)
>
> ## ------- examples from ?svd ---------
>
> hilbert <- function(n) { i <- 1:n; 1 / outer(i - 1, i, "+") }
> Eps <- 100 * .Machine$double.eps
>
> X <- hilbert(9)[,1:6]
> (s <- svd(X)); D <- diag(s$d)
$d
[1] 1.668e+00 2.774e-01 2.224e-02 1.085e-03 3.244e-05 5.235e-07
$u
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] -0.7245 0.6266 0.27350 -0.08527 0.02074 -0.004025
[2,] -0.4282 -0.1299 -0.64294 0.55047 -0.27253 0.092816
[3,] -0.3122 -0.2804 -0.33633 -0.31418 0.61632 -0.440904
[4,] -0.2479 -0.3142 -0.06931 -0.44667 0.02945 0.530120
[5,] -0.2064 -0.3141 0.10786 -0.30242 -0.35567 0.237038
[6,] -0.1771 -0.3027 0.22106 -0.09042 -0.38879 -0.260449
[7,] -0.1553 -0.2877 0.29281 0.11551 -0.19286 -0.420945
[8,] -0.1384 -0.2722 0.33784 0.29313 0.11633 -0.160790
[9,] -0.1249 -0.2571 0.36543 0.43885 0.46497 0.434600
$v
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] -0.7365 0.6225 0.2550 -0.06976 0.01328 -0.001588
[2,] -0.4433 -0.1819 -0.6867 0.50860 -0.19627 0.041117
[3,] -0.3275 -0.3509 -0.2611 -0.50474 0.61606 -0.259216
[4,] -0.2626 -0.3922 0.1044 -0.43748 -0.40834 0.638902
[5,] -0.2204 -0.3946 0.3510 0.01612 -0.46428 -0.675827
[6,] -0.1904 -0.3832 0.5111 0.53856 0.44664 0.257249
> stopifnot(abs(X - s$u %*% D %*% t(s$v)) < Eps)# X = U D V'
> stopifnot(abs(D - t(s$u) %*% X %*% s$v) < Eps)# D = U' X V
>
> # The signs of the vectors are not determined here.
> X <- cbind(1, 1:7)
> s <- svd(X); D <- diag(s$d)
> stopifnot(abs(X - s$u %*% D %*% t(s$v)) < Eps)# X = U D V'
> stopifnot(abs(D - t(s$u) %*% X %*% s$v) < Eps)# D = U' X V
>
> # test nu and nv
> s <- svd(X, nu = 0)
> s <- svd(X, nu = 7) # the last 5 columns are not determined here
> stopifnot(dim(s$u) == c(7,7))
> s <- svd(X, nv = 0)
>
> # test of complex case
>
> X <- cbind(1, 1:7+(-3:3)*1i)
> s <- svd(X); D <- diag(s$d)
> stopifnot(abs(X - s$u %*% D %*% Conj(t(s$v))) < Eps)
> stopifnot(abs(D - Conj(t(s$u)) %*% X %*% s$v) < Eps)
>
>
>
> ## ------- tests of random real and complex matrices ------
> fixsign <- function(A) {
+ A[] <- apply(A, 2, function(x) x*sign(Re(x[1])))
+ A
+ }
> ## 100 may cause failures here.
> eigenok <- function(A, E, Eps=1000*.Machine$double.eps)
+ {
+ print(fixsign(E$vectors))
+ print(zapsmall(E$values))
+ V <- E$vectors; lam <- E$values
+ stopifnot(abs(A %*% V - V %*% diag(lam)) < Eps,
+ abs(lam[length(lam)]/lam[1]) < Eps || # this one not for singular A :
+ abs(A - V %*% diag(lam) %*% t(V)) < Eps)
+ }
>
> Ceigenok <- function(A, E, Eps=1000*.Machine$double.eps)
+ {
+ print(fixsign(E$vectors))
+ print(signif(E$values, 5))
+ V <- E$vectors; lam <- E$values
+ stopifnot(Mod(A %*% V - V %*% diag(lam)) < Eps,
+ Mod(A - V %*% diag(lam) %*% Conj(t(V))) < Eps)
+ }
>
> ## failed for some 64bit-Lapack-gcc combinations:
> sm <- cbind(1, 3:1, 1:3)
> eigenok(sm, eigen(sm))
[,1] [,2] [,3]
[1,] 0.5774 0.8452 0.9428
[2,] 0.5774 0.1690 -0.2357
[3,] 0.5774 -0.5071 -0.2357
[1] 5 1 0
> eigenok(sm, eigen(sm, sym=FALSE))
[,1] [,2] [,3]
[1,] 0.5774 0.8452 0.9428
[2,] 0.5774 0.1690 -0.2357
[3,] 0.5774 -0.5071 -0.2357
[1] 5 1 0
>
> set.seed(123)
> sm <- matrix(rnorm(25), 5, 5)
> sm <- 0.5 * (sm + t(sm))
> eigenok(sm, eigen(sm))
[,1] [,2] [,3] [,4] [,5]
[1,] 0.5899 0.1683 0.02315 0.471808 0.6329
[2,] 0.1936 0.2931 0.89217 -0.009784 -0.2838
[3,] 0.6627 -0.4812 -0.15825 0.082550 -0.5454
[4,] 0.1404 0.7985 -0.41848 0.094314 -0.3983
[5,] -0.3946 -0.1285 0.05768 0.872692 -0.2507
[1] 1.7814 1.5184 0.5833 -1.0148 -2.4908
> eigenok(sm, eigen(sm, sym=FALSE))
[,1] [,2] [,3] [,4] [,5]
[1,] 0.6329 0.5899 0.1683 0.471808 0.02315
[2,] -0.2838 0.1936 0.2931 -0.009784 0.89217
[3,] -0.5454 0.6627 -0.4812 0.082550 -0.15825
[4,] -0.3983 0.1404 0.7985 0.094314 -0.41848
[5,] -0.2507 -0.3946 -0.1285 0.872692 0.05768
[1] -2.4908 1.7814 1.5184 -1.0148 0.5833
>
> sm[] <- as.complex(sm)
> Ceigenok(sm, eigen(sm))
[,1] [,2] [,3] [,4] [,5]
[1,] 0.5899+0i 0.1683+0i 0.02315+0i 0.471808+0i 0.6329+0i
[2,] 0.1936+0i 0.2931+0i 0.89217+0i -0.009784+0i -0.2838+0i
[3,] 0.6627+0i -0.4812+0i -0.15825+0i 0.082550+0i -0.5454+0i
[4,] 0.1404+0i 0.7985+0i -0.41848+0i 0.094314+0i -0.3983+0i
[5,] -0.3946+0i -0.1285+0i 0.05768+0i 0.872692+0i -0.2507+0i
[1] 1.7814 1.5184 0.5833 -1.0148 -2.4908
> Ceigenok(sm, eigen(sm, sym=FALSE))
[,1] [,2] [,3] [,4] [,5]
[1,] 0.6329+0i 0.5899+0i 0.1683+0i 0.471808+0i 0.02315+0i
[2,] -0.2838+0i 0.1936+0i 0.2931+0i -0.009784+0i 0.89217+0i
[3,] -0.5454+0i 0.6627+0i -0.4812+0i 0.082550+0i -0.15825+0i
[4,] -0.3983+0i 0.1404+0i 0.7985+0i 0.094314+0i -0.41848+0i
[5,] -0.2507+0i -0.3946+0i -0.1285+0i 0.872692+0i 0.05768+0i
[1] -2.4908+0i 1.7814+0i 1.5184+0i -1.0148+0i 0.5833+0i
>
> sm[] <- sm + rnorm(25) * 1i
> sm <- 0.5 * (sm + Conj(t(sm)))
> Ceigenok(sm, eigen(sm))
[,1] [,2] [,3] [,4]
[1,] 0.5373+0.0000i 0.3338+0.0000i 0.02834+0.0000i 0.4378+0.0000i
[2,] 0.3051+0.0410i -0.0264-0.1175i -0.43963+0.7256i -0.0474+0.2975i
[3,] 0.3201-0.3756i 0.3379+0.4760i -0.09325-0.3281i 0.0536+0.2447i
[4,] 0.3394+0.2330i -0.1044-0.6839i 0.09966-0.3629i 0.1894+0.1979i
[5,] -0.2869+0.3483i -0.0766+0.2210i -0.14602+0.0132i 0.7449-0.1576i
[,5]
[1,] 0.6383+0.0000i
[2,] -0.1909-0.2093i
[3,] -0.4788-0.0861i
[4,] -0.3654+0.0418i
[5,] -0.2229-0.3012i
[1] 2.4043 1.3934 0.7854 -1.4050 -2.8006
> Ceigenok(sm, eigen(sm, sym=FALSE))
[,1] [,2] [,3] [,4]
[1,] 0.6383+0.0000i 0.5373+0.0000i 0.4283+0.0907i 0.0504-0.3300i
[2,] -0.1909-0.2093i 0.3051+0.0410i -0.1080+0.2813i -0.1201+0.0084i
[3,] -0.4788-0.0861i 0.3201-0.3756i 0.0018+0.2505i 0.5216-0.2622i
[4,] -0.3654+0.0418i 0.3394+0.2330i 0.1443+0.2329i -0.6918+0.0000i
[5,] -0.2229-0.3012i -0.2869+0.3483i 0.7614+0.0000i 0.2069+0.1091i
[,5]
[1,] 0.01468+0.02424i
[2,] -0.84836+0.00000i
[3,] 0.23232-0.24980i
[4,] 0.36200-0.10282i
[5,] -0.08698-0.11804i
[1] -2.8006+0i 2.4043+0i -1.4050+0i 1.3934+0i 0.7854+0i
>
>
> ## ------- tests of integer matrices -----------------
>
> set.seed(123)
> A <- matrix(rpois(25, 5), 5, 5)
>
> A %*% A
[,1] [,2] [,3] [,4] [,5]
[1,] 202 170 156 160 234
[2,] 161 124 145 147 185
[3,] 166 136 134 130 174
[4,] 218 156 169 204 234
[5,] 205 134 175 181 249
> crossprod(A)
[,1] [,2] [,3] [,4] [,5]
[1,] 226 160 153 174 240
[2,] 160 143 126 112 179
[3,] 153 126 171 137 205
[4,] 174 112 137 174 192
[5,] 240 179 205 192 293
> tcrossprod(A)
[,1] [,2] [,3] [,4] [,5]
[1,] 229 155 150 207 184
[2,] 155 144 140 184 161
[3,] 150 140 156 176 142
[4,] 207 184 176 251 209
[5,] 184 161 142 209 227
>
> solve(A)
[,1] [,2] [,3] [,4] [,5]
[1,] -0.048676 0.3390 -0.15756 -0.05892 -0.00854
[2,] -0.058711 -0.1262 0.19812 -0.03160 0.06426
[3,] 0.092656 0.2319 -0.02904 -0.12468 -0.09778
[4,] 0.062553 -0.1637 0.03800 -0.04270 0.12062
[5,] -0.002775 -0.2351 0.02391 0.22032 -0.02242
> qr(A)
$qr
[,1] [,2] [,3] [,4] [,5]
[1,] -15.0333 -10.64304 -10.1774 -11.5743 -15.965
[2,] 0.4656 -5.45212 -3.2430 2.0516 -1.667
[3,] 0.2661 0.97998 -7.5434 -3.4278 -4.920
[4,] 0.5322 -0.05761 -0.3972 4.9068 -1.269
[5,] 0.5987 -0.17944 -0.9104 -0.9086 3.088
$rank
[1] 5
$qraux
[1] 1.266 1.064 1.116 1.418 3.088
$pivot
[1] 1 2 3 4 5
attr(,"class")
[1] "qr"
> determinant(A, log = FALSE)
$modulus
[1] 9368
attr(,"logarithm")
[1] FALSE
$sign
[1] -1
attr(,"class")
[1] "det"
>
> rcond(A)
[1] 0.02466
> rcond(A, "I")
[1] 0.06007
> rcond(A, "1")
[1] 0.02466
>
> eigen(A)
$values
[1] 29.660+0.000i -4.631+0.000i 4.556+0.000i -2.292+3.117i -2.292-3.117i
$vectors
[,1] [,2] [,3] [,4] [,5]
[1,] 0.4698+0i 0.34581+0i 0.07933+0i 0.7246+0.0000i 0.7246+0.0000i
[2,] 0.3885+0i -0.30397+0i 0.31927+0i -0.2481-0.2591i -0.2481+0.2591i
[3,] 0.3760+0i -0.03566+0i 0.70330+0i 0.0496+0.3697i 0.0496-0.3697i
[4,] 0.5029+0i -0.74239+0i -0.32651+0i -0.2262+0.1063i -0.2262-0.1063i
[5,] 0.4839+0i 0.48539+0i -0.53901+0i -0.3375-0.1752i -0.3375+0.1752i
> svd(A)
$d
[1] 29.929 6.943 6.668 3.960 1.707
$u
[,1] [,2] [,3] [,4] [,5]
[1,] -0.4659 -0.04141 -0.8795 8.440e-02 0.02379
[2,] -0.3931 0.15031 0.2249 1.824e-05 0.87881
[3,] -0.3811 0.62205 0.2170 5.570e-01 -0.33239
[4,] -0.5166 0.14296 0.1852 -7.659e-01 -0.30290
[5,] -0.4652 -0.75386 0.3073 3.097e-01 -0.15780
$v
[,1] [,2] [,3] [,4] [,5]
[1,] -0.4831 -0.3264 0.47567 -0.1955 0.6289
[2,] -0.3627 0.3731 0.53450 0.5920 -0.3051
[3,] -0.3995 0.4779 -0.59210 0.2252 0.4590
[4,] -0.3983 -0.6985 -0.36298 0.3821 -0.2752
[5,] -0.5628 0.1948 -0.07549 -0.6439 -0.4743
> La.svd(A)
$d
[1] 29.929 6.943 6.668 3.960 1.707
$u
[,1] [,2] [,3] [,4] [,5]
[1,] -0.4659 -0.04141 -0.8795 8.440e-02 0.02379
[2,] -0.3931 0.15031 0.2249 1.824e-05 0.87881
[3,] -0.3811 0.62205 0.2170 5.570e-01 -0.33239
[4,] -0.5166 0.14296 0.1852 -7.659e-01 -0.30290
[5,] -0.4652 -0.75386 0.3073 3.097e-01 -0.15780
$vt
[,1] [,2] [,3] [,4] [,5]
[1,] -0.4831 -0.3627 -0.3995 -0.3983 -0.56284
[2,] -0.3264 0.3731 0.4779 -0.6985 0.19477
[3,] 0.4757 0.5345 -0.5921 -0.3630 -0.07549
[4,] -0.1955 0.5920 0.2252 0.3821 -0.64390
[5,] 0.6289 -0.3051 0.4590 -0.2752 -0.47430
>
> As <- crossprod(A)
> E <- eigen(As)
> E$values
[1] 895.737 48.201 44.468 15.678 2.915
> abs(E$vectors) # signs vary
[,1] [,2] [,3] [,4] [,5]
[1,] 0.4831 0.3264 0.47567 0.1955 0.6289
[2,] 0.3627 0.3731 0.53450 0.5920 0.3051
[3,] 0.3995 0.4779 0.59210 0.2252 0.4590
[4,] 0.3983 0.6985 0.36298 0.3821 0.2752
[5,] 0.5628 0.1948 0.07549 0.6439 0.4743
> chol(As)
[,1] [,2] [,3] [,4] [,5]
[1,] 15.03 10.643 10.177 11.574 15.965
[2,] 0.00 5.452 3.243 -2.052 1.667
[3,] 0.00 0.000 7.543 3.428 4.920
[4,] 0.00 0.000 0.000 4.907 -1.269
[5,] 0.00 0.000 0.000 0.000 3.088
> backsolve(As, 1:5)
[1] -0.009040 -0.005129 -0.006246 0.004158 0.017065
>
> ## ------- tests of logical matrices -----------------
>
> set.seed(123)
> A <- matrix(runif(25) > 0.5, 5, 5)
>
> A %*% A
[,1] [,2] [,3] [,4] [,5]
[1,] 2 2 2 1 3
[2,] 2 1 1 2 3
[3,] 2 2 1 1 3
[4,] 2 2 2 2 4
[5,] 2 1 2 2 3
> crossprod(A)
[,1] [,2] [,3] [,4] [,5]
[1,] 3 2 1 1 3
[2,] 2 3 2 0 3
[3,] 1 2 3 1 3
[4,] 1 0 1 2 2
[5,] 3 3 3 2 5
> tcrossprod(A)
[,1] [,2] [,3] [,4] [,5]
[1,] 3 1 2 2 2
[2,] 1 3 2 3 2
[3,] 2 2 3 3 1
[4,] 2 3 3 4 2
[5,] 2 2 1 2 3
>
> Q <- qr(A)
> zapsmall(Q$qr)
[,1] [,2] [,3] [,4] [,5]
[1,] -1.7321 -1.1547 -0.5774 -0.5774 -1.7321
[2,] 0.5774 -1.2910 -1.0328 0.5164 -0.7746
[3,] 0.0000 0.7746 -1.2649 -0.9487 -0.9487
[4,] 0.5774 0.2582 0.0508 0.7071 0.7071
[5,] 0.5774 -0.5164 -0.6803 -0.3136 0.0000
> zapsmall(Q$qraux)
[1] 1.000 1.258 1.731 1.950 0.000
> determinant(A, log = FALSE) # 0
$modulus
[1] 0
attr(,"logarithm")
[1] FALSE
$sign
[1] 1
attr(,"class")
[1] "det"
>
> rcond(A)
[1] 0
> rcond(A, "I")
[1] 0
> rcond(A, "1")
[1] 0
>
> E <- eigen(A)
> zapsmall(E$values)
[1] 3.163+0.000i 0.271+0.908i 0.271-0.908i -0.705+0.000i 0.000+0.000i
> zapsmall(Mod(E$vectors))
[,1] [,2] [,3] [,4] [,5]
[1,] 0.4358 0.3604 0.3604 0.6113 0.0000
[2,] 0.4087 0.4495 0.4495 0.3771 0.5774
[3,] 0.3962 0.5870 0.5870 0.2028 0.0000
[4,] 0.5340 0.2792 0.2792 0.6649 0.5774
[5,] 0.4483 0.4955 0.4955 0.0314 0.5774
> S <- svd(A)
> zapsmall(S$d)
[1] 3.379 1.536 1.414 0.472 0.000
> S <- La.svd(A)
> zapsmall(S$d)
[1] 3.379 1.536 1.414 0.472 0.000
>
> As <- A
> As[upper.tri(A)] <- t(A)[upper.tri(A)]
> det(As)
[1] 2
> E <- eigen(As)
> E$values
[1] 3.465 1.510 0.300 -1.000 -1.275
> zapsmall(E$vectors)
[,1] [,2] [,3] [,4] [,5]
[1,] -0.4023 0.2877 0.3638 0.5 0.6108
[2,] -0.5338 -0.3474 0.5742 -0.5 -0.1207
[3,] -0.4177 -0.5380 -0.6384 0.0 0.3585
[4,] -0.4959 0.0733 -0.1273 0.5 -0.6946
[5,] -0.3644 0.7084 -0.3378 -0.5 0.0369
> solve(As)
[,1] [,2] [,3] [,4] [,5]
[1,] 0 1 -1.0 0.0 0.0
[2,] 1 1 -1.0 0.0 -1.0
[3,] -1 -1 1.5 0.5 0.5
[4,] 0 0 0.5 -0.5 0.5
[5,] 0 -1 0.5 0.5 0.5
>
> ## quite hard to come up with an example where this might make sense.
> Ac <- A; Ac[] <- as.logical(diag(5))
> chol(Ac)
[,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
>
>
>