# SCM Repository

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

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

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: