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
|
#' Header for spatstat.sparse/tests/*R
#'
require(spatstat.sparse)
ALWAYS <- FULLTEST <- TRUE
##
## tests/linalgeb.R
##
## checks validity of linear algebra code
##
## $Revision: 1.13 $ $Date: 2022/04/18 03:16:26 $
##
local({
p <- 3
n <- 4
k <- 2
## correctness of 'quadform'
x <- matrix(1:(n*p), n, p)
v <- matrix(runif(p^2), p, p)
zW <- zU <- numeric(n)
for(i in 1:n) {
xi <- x[i,,drop=FALSE]
zW[i] <- xi %*% v %*% t(xi)
zU[i] <- xi %*% t(xi)
}
if(!all(zU == quadform(x)))
stop("quadform gives incorrect result in Unweighted case")
if(!all(zW == quadform(x,v)))
stop("quadform gives incorrect result in Weighted case")
## correctness of 'sumouter'
w <- runif(n)
y <- matrix(1:(2*n), n, k)
zUS <- zWS <- matrix(0, p, p)
zUA <- zWA <- matrix(0, p, k)
for(i in 1:n) {
zUS <- zUS + outer(x[i,],x[i,])
zWS <- zWS + w[i] * outer(x[i,],x[i,])
zUA <- zUA + outer(x[i,],y[i,])
zWA <- zWA + w[i] * outer(x[i,],y[i,])
}
if(!identical(zUS, sumouter(x)))
stop("sumouter gives incorrect result in Unweighted Symmetric case")
if(!identical(zWS, sumouter(x,w)))
stop("sumouter gives incorrect result in Weighted Symmetric case")
if(!identical(zUA, sumouter(x, y=y)))
stop("sumouter gives incorrect result in Unweighted Asymmetric case")
if(!identical(zWA, sumouter(x, w, y)))
stop("sumouter gives incorrect result in Weighted Asymmetric case")
#' complex quadratic forms - execute only
dimnames(x) <- list(letters[1:nrow(x)], LETTERS[1:ncol(x)])
a <- quadform(x + 1i)
b <- quadform(x + 1i, v+1i)
d <- quadform(x , v+1i)
a <- sumouter(x + 1i)
b <- sumouter(x + 1i, w + 1i)
d <- sumouter(x + 1i, w + 1i, x - 1i)
#' NA values
xna <- x; xna[1,1] <- NA
wna <- w; w[2] <- NA
vna <- v; v[1,2] <- NA
o <- quadform(xna)
o <- quadform(xna, vna)
o <- sumouter(xna)
o <- sumouter(xna, w)
o <- sumouter(xna, wna)
#' sumsymouter
x <- array(as.numeric(1:(p * n * n)), dim=c(p, n, n))
w <- matrix(1:(n*n), n, n)
y <- matrix(numeric(p * p), p, p)
#' check correctness
for(i in 1:n)
for(j in (1:n)[-i])
y <- y + w[i,j] * outer(x[,i,j], x[,j,i])
z <- sumsymouter(x, w)
if(!identical(y,z))
stop("sumsymouter gives incorrect result")
#' cover code blocks
o <- sumsymouter(x, distinct=FALSE)
a <- sumsymouter(x + 1i)
b <- sumsymouter(x + 1i, w + 1i)
if(require(Matrix)) o <- sumsymouter(x, as(w, "sparseMatrix"))
#' matrixpower, matrixsqrt, matrixinvsqrt
checkit <- function(A, B=diag(nrow(A)), what) {
discrep <- max(abs(A-B))
if(discrep > sqrt(.Machine$double.eps))
stop(paste("large discrepancy", discrep, "in", what), call.=FALSE)
return(discrep)
}
#' (a) power of matrix is complex
M <- diag(c(4,-4))
dimnames(M) <- list(letters[1:2], letters[1:2])
V <- matrixsqrt(M)
U <- matrixinvsqrt(M)
V2 <- matrixpower(M, 1/2)
checkit(V %*% V, M, "square of matrixsqrt")
checkit(V %*% U, what="matrixsqrt * matrixinvsqrt")
checkit(V2 %*% V2, M, "square of matrixpower(1/2)")
W <- matrixsqrt(abs(M), complexOK=FALSE)
#' (b) power of asymmetric complex matrix
Z <- matrix(c(1+1i, 2+1i, 2+3i, 5+5i), 2, 2)
V <- matrixsqrt(Z)
U <- matrixinvsqrt(Z)
V2 <- matrixpower(Z, 1/2)
checkit(V %*% V, Z, "square of matrixsqrt (complex)")
checkit(V %*% U, what="matrixsqrt * matrixinvsqrt (complex)")
checkit(V2 %*% V2, Z, "square of matrixpower(1/2) (complex)")
#' infrastructure
A <- matrix(1:12, 3, 4)
B <- matrix(1:8, 4, 2)
check.mat.mul(A, B)
check.mat.mul(A, B[,1])
check.mat.mul(A, A, fatal=FALSE)
D <- diag(c(1,4,9))
checksolve(D)
D[1,1] <- 0
checksolve(D, "silent")
})
|
