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 907 - (view) (download)

1 : maechler 522 library(Matrix)
2 :    
3 : maechler 519 ## Matrix Exponential
4 :    
5 : maechler 907 source(system.file("test-tools.R", package = "Matrix"))
6 : maechler 519
7 : maechler 538 ## e ^ 0 = 1 - for matrices:
8 :     assert.EQ.mat(expm(Matrix(0, 3,3)), diag(3), tol = 0)# exactly
9 :     ## e ^ diag(.) = diag(e ^ .):
10 :     assert.EQ.mat(expm(as(diag(-1:4), "dgeMatrix")), diag(exp(-1:4)))
11 :     set.seed(1)
12 :     rE <- replicate(100,
13 :     { x <- rlnorm(12)
14 :     relErr(as(expm(as(diag(x), "dgeMatrix")),
15 :     "matrix"),
16 :     diag(exp(x))) })
17 :     stopifnot(mean(rE) < 1e-15,
18 :     max(rE) < 1e-14)
19 :     summary(rE)
20 :    
21 :     ## Some small matrices
22 :    
23 : maechler 519 m1 <- Matrix(c(1,0,1,1), nc = 2)
24 :     e1 <- expm(m1)
25 : maechler 521 assert.EQ.mat(e1, cbind(c(exp(1),0), exp(1)))
26 : maechler 519
27 :     m2 <- Matrix(c(-49, -64, 24, 31), nc = 2)
28 :     e2 <- expm(m2)
29 : maechler 532 ## The true matrix exponential is 'te2':
30 :     e_1 <- exp(-1)
31 :     e_17 <- exp(-17)
32 :     te2 <- rbind(c(3*e_17 - 2*e_1, -3/2*e_17 + 3/2*e_1),
33 :     c(4*e_17 - 4*e_1, -2 *e_17 + 3 *e_1))
34 :     assert.EQ.mat(e2, te2, tol = 1e-13)
35 :     ## See the (average relative) difference:
36 :     all.equal(as(e2,"matrix"), te2, tol = 0) # 1.48e-14 on "lynne"
37 :    
38 : maechler 520 ## The ``surprising identity'' det(exp(A)) == exp( tr(A) )
39 :     ## or log det(exp(A)) == tr(A) :
40 :     stopifnot(all.equal(determinant(e2)$modulus, sum(diag(m2))))
41 : maechler 519
42 :     m3 <- Matrix(cbind(0,rbind(6*diag(3),0)), nc = 4)
43 :     e3 <- expm(m3)
44 :     assert.EQ.mat(e3,
45 :     rbind(c(1,6,18,36),
46 :     c(0,1, 6,18),
47 :     c(0,0, 1, 6),
48 :     c(0,0, 0, 1)))
49 :    
50 : 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