1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
|
## tests of R functions based on the lapack module
## NB: the signs of singular and eigenvectors are arbitrary,
## so there may be differences from the reference ouptut,
## especially when alternative BLAS are used.
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)
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))
eigenok(sm, eigen(sm, sym=FALSE))
set.seed(123)
sm <- matrix(rnorm(25), 5, 5)
sm <- 0.5 * (sm + t(sm))
eigenok(sm, eigen(sm))
eigenok(sm, eigen(sm, sym=FALSE))
sm[] <- as.complex(sm)
Ceigenok(sm, eigen(sm))
Ceigenok(sm, eigen(sm, sym=FALSE))
sm[] <- sm + rnorm(25) * 1i
sm <- 0.5 * (sm + Conj(t(sm)))
Ceigenok(sm, eigen(sm))
Ceigenok(sm, eigen(sm, sym=FALSE))
## ------- tests of integer matrices -----------------
set.seed(123)
A <- matrix(rpois(25, 5), 5, 5)
A %*% A
crossprod(A)
tcrossprod(A)
solve(A)
qr(A)
determinant(A, log = FALSE)
rcond(A)
rcond(A, "I")
rcond(A, "1")
eigen(A)
svd(A)
La.svd(A)
As <- crossprod(A)
E <- eigen(As)
E$values
abs(E$vectors) # signs vary
chol(As)
backsolve(As, 1:5)
## ------- tests of logical matrices -----------------
set.seed(123)
A <- matrix(runif(25) > 0.5, 5, 5)
A %*% A
crossprod(A)
tcrossprod(A)
Q <- qr(A)
zapsmall(Q$qr)
zapsmall(Q$qraux)
determinant(A, log = FALSE) # 0
rcond(A)
rcond(A, "I")
rcond(A, "1")
E <- eigen(A)
zapsmall(E$values)
zapsmall(Mod(E$vectors))
S <- svd(A)
zapsmall(S$d)
S <- La.svd(A)
zapsmall(S$d)
As <- A
As[upper.tri(A)] <- t(A)[upper.tri(A)]
det(As)
E <- eigen(As)
E$values
zapsmall(E$vectors)
solve(As)
## quite hard to come up with an example where this might make sense.
Ac <- A; Ac[] <- as.logical(diag(5))
chol(Ac)
|
