SCM

SCM Repository

[matrix] Annotation of /pkg/tests/matr-exp.R
ViewVC logotype

Annotation of /pkg/tests/matr-exp.R

Parent Directory Parent Directory | Revision Log Revision Log


Revision 561 - (view) (download)

1 : maechler 522 library(Matrix)
2 :    
3 : maechler 519 ## Matrix Exponential
4 :    
5 :    
6 : maechler 538 ## checking; 'show' is for convenience of the developer
7 :     assert.EQ.mat <- function(M, m, tol = if(show) 0 else 1e-15, show=FALSE) {
8 : bates 561 ## temporary fix for R-2.0.1
9 :     MM <- as(M, "matrix")
10 :     attr(MM, "dimnames") <- NULL
11 :     if(show) all.equal(MM, m, tol = tol)
12 :     else stopifnot(all.equal(MM, m, tol = tol))
13 : maechler 538 }
14 :     ## The relative error typically returned by all.equal:
15 :     relErr <- function(target, current)
16 :     mean(abs(target - current)) / mean(abs(target))
17 : maechler 519
18 : maechler 538 ## e ^ 0 = 1 - for matrices:
19 :     assert.EQ.mat(expm(Matrix(0, 3,3)), diag(3), tol = 0)# exactly
20 :     ## e ^ diag(.) = diag(e ^ .):
21 :     assert.EQ.mat(expm(as(diag(-1:4), "dgeMatrix")), diag(exp(-1:4)))
22 :     set.seed(1)
23 :     rE <- replicate(100,
24 :     { x <- rlnorm(12)
25 :     relErr(as(expm(as(diag(x), "dgeMatrix")),
26 :     "matrix"),
27 :     diag(exp(x))) })
28 :     stopifnot(mean(rE) < 1e-15,
29 :     max(rE) < 1e-14)
30 :     summary(rE)
31 :    
32 :     ## Some small matrices
33 :    
34 : maechler 519 m1 <- Matrix(c(1,0,1,1), nc = 2)
35 :     e1 <- expm(m1)
36 : maechler 521 assert.EQ.mat(e1, cbind(c(exp(1),0), exp(1)))
37 : maechler 519
38 :     m2 <- Matrix(c(-49, -64, 24, 31), nc = 2)
39 :     e2 <- expm(m2)
40 : maechler 532 ## The true matrix exponential is 'te2':
41 :     e_1 <- exp(-1)
42 :     e_17 <- exp(-17)
43 :     te2 <- rbind(c(3*e_17 - 2*e_1, -3/2*e_17 + 3/2*e_1),
44 :     c(4*e_17 - 4*e_1, -2 *e_17 + 3 *e_1))
45 :     assert.EQ.mat(e2, te2, tol = 1e-13)
46 :     ## See the (average relative) difference:
47 :     all.equal(as(e2,"matrix"), te2, tol = 0) # 1.48e-14 on "lynne"
48 :    
49 : maechler 520 ## The ``surprising identity'' det(exp(A)) == exp( tr(A) )
50 :     ## or log det(exp(A)) == tr(A) :
51 :     stopifnot(all.equal(determinant(e2)$modulus, sum(diag(m2))))
52 : maechler 519
53 :     m3 <- Matrix(cbind(0,rbind(6*diag(3),0)), nc = 4)
54 :     e3 <- expm(m3)
55 :     assert.EQ.mat(e3,
56 :     rbind(c(1,6,18,36),
57 :     c(0,1, 6,18),
58 :     c(0,0, 1, 6),
59 :     c(0,0, 0, 1)))
60 :    
61 : maechler 538 proc.time() # for ``statistical reasons''

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge