tests/log+sqrt.R

library(Matrix)
library(expm)

source(system.file("test-tools.R", package = "expm"), keep.source=FALSE)
doExtras

tst.sqrtm <- function(m, tol = 1e-12, zap.Im.tol = 1e-10) {
    r.m <- sqrtm(m)## should now work
    s <- r.m %*% r.m
    if(is.complex(s) && all(abs(Im(s)) < mean(abs(s)) * zap.Im.tol))
        s <- Re(s)
    all.equal(m, s, tolerance=tol)
}


options(verbose = TRUE) # -> get some messages from logm.Higham

### ---- Small exact : ----------
 L2 <- cbind(1, 0:1)
lL2 <- cbind(0:1, 0)
(L3 <- rbind(cbind(1,cbind(0:1,0)),1))
(lL3 <- cbind(rbind(0, cbind((2:1)/2,0:1)), 0))

assertError(logm(L2, method="Eigen"))
assertError(logm(L3, method="Eigen"))

logm.Higham08 <- expm:::logm.Higham08
l.L2 <- logm.Higham08(L2)
l.L3 <- logm.Higham08(L3)

all.equal(l.L2, lL2, tolerance=0)# 5.64 e-14 (32-bit *and* 64-bit)
all.equal(l.L3, lL3, tolerance=0)# 2.40 e-15   (ditto)
stopifnot(all.equal(l.L2, lL2, tolerance= 1000e-16),
          all.equal(l.L3, lL3, tolerance=   80e-16))
showProc.time()

### --------- More & larger randomly generated examples : -----------------
set.seed(101)

EA <- expm.Higham08(A <- matrix(round(rnorm(25),1), 5))
all.equal(EA, expm.Higham08(logm.Higham08(EA)), tolerance=0)
## "Mean relative difference: 1.020137e-13"
stopifnot(all.equal(EA, expm.Higham08(logm.Higham08(EA)), tolerance=1e-12))

S <- crossprod(A)
all.equal(S, sqrtm(S) %*% sqrtm(S), tolerance=0)
## "Mean relative difference: 2.26885e-15"
stopifnot(all.equal(S, sqrtm(S) %*% sqrtm(S), tolerance=1e-14))
showProc.time()

set.seed(3)

## n = 50 is already "too" slow (well: logm.Higham08(.) needs 2.2 sec
## --> CPU measurements below
for(n in c(2:5, 10:11, if(doExtras) 30)) {
    cat("n = ",n,": ")
    for(kk in seq_len(if(doExtras) 30 else 10)) {
        ## Testing  logm()
        EA <- expm.Higham08(A <- matrix(round(rnorm(n^2),2), n,n))
        stopifnot(all.equal(EA, expm.Higham08(logm.Higham08(EA)), tolerance=1e-12))
        cat(" ")
        ## Testing  sqrtm() --- for positive definite *and* arbitrary
        stopifnot(tst.sqrtm(A))# A is completely random
        S <- crossprod(A) + rnorm(n^2) / 1000
        stopifnot(tst.sqrtm(S))
        cat(".")
    }
    cat("\n")
}

showProc.time()


### CPU-measurements for  logm()
options(verbose = FALSE)# printing costs ..
set.seed(5)
if(doExtras) {
    n <- 50
    sim <- 32
} else {
    n <- 21
    sim <- 8
}

cpuT <- numeric(sim)
for(k in seq_len(sim)) {
    EA <- expm.Higham08(A <- matrix(rnorm(n^2), n,n))
    cat(".")
    cpuT[k] <- system.time(LEA <- logm.Higham08(EA))[1]
    stopifnot(all.equal(EA, expm.Higham08(LEA), tolerance=1e-12))
}; cat("\n")
summary(cpuT)
## cmath-5 {Feb.2009}; Michi's original code:
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
##   1.794   2.249   2.389   2.388   2.515   2.831

showProc.time()

Try the expm package in your browser

Any scripts or data that you put into this service are public.

expm documentation built on May 2, 2019, 5:25 p.m.