R Under development (unstable) (2023-11-28 r85645 ucrt) -- "Unsuffered Consequences" Copyright (C) 2023 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(expm) Loading required package: Matrix Attaching package: 'expm' The following object is masked from 'package:Matrix': expm > source(system.file("test-tools.R", package= "expm"), keep.source=FALSE)## -> assertError()... Loading required package: tools > > ## A matrix with 'Inf' > mI <- rbind(0, c(-Inf, Inf, 0, 0), 0, 0) > bal3 <- + list(dB = dgebal(mI, "B"), # = default + dP = dgebal(mI, "P"), + dN = dgebal(mI, "N")) Warning messages: 1: 'dgebal' is deprecated. Use 'balance' instead. See help("Deprecated") 2: 'dgebal' is deprecated. Use 'balance' instead. See help("Deprecated") 3: 'dgebal' is deprecated. Use 'balance' instead. See help("Deprecated") > str(bal3) List of 3 $ dB:List of 4 ..$ z : num [1:4, 1:4] Inf 0 0 0 -Inf ... ..$ scale: num [1:4] 1 1 3 4 ..$ i1 : int 1 ..$ i2 : int 1 $ dP:List of 4 ..$ z : num [1:4, 1:4] Inf 0 0 0 -Inf ... ..$ scale: num [1:4] 1 1 3 4 ..$ i1 : int 1 ..$ i2 : int 1 $ dN:List of 4 ..$ z : num [1:4, 1:4] 0 -Inf 0 0 0 ... ..$ scale: num [1:4] 1 1 1 1 ..$ i1 : int 1 ..$ i2 : int 4 > stopifnot(identical(mI, bal3$dN$z), + with(bal3, all.equal(dB, dP, tol=1e-14)), + all.equal(bal3$dB$z, rbind(c(Inf,-Inf,0,0), 0,0,0), tol=1e-14), + all.equal(bal3$dB$scale, c(1,1,3,4))) > assertError(dgebal(mI, "S"))# gave infinite loop > > > > ## Compare the two different "balance" pre-conditioning versions in Ward77: > set.seed(1) > mList <- lapply(integer(100), function(...) rSpMatrix(20, nnz=80)) > re20 <- sapply(mList, function(M) + relErr(expm(M, precond = "2bal"), + expm(M, precond = "1bal"))) > re20 ## ahh.: zero or ~ 1e-13 ... good [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 > table(re20 == 0) TRUE 100 > summary(re20[re20 != 0]) Min. 1st Qu. Median Mean 3rd Qu. Max. > ## Pentium M (ubuntu) > ## Min. 1st Qu. Median Mean 3rd Qu. Max. > ## 2.593e-14 8.703e-14 1.282e-13 2.434e-13 4.177e-13 6.295e-13 > > > demo(balanceTst) #-> the function definition and the first few examples demo(balanceTst) ---- ~~~~~~~~~~ > balanceTst <- function(A) { + + ## Purpose: Consistency checking of balance() {was "dgebal()"} + ## ---------------------------------------------------------------------- + ## Arguments: a square matrix + ## ---------------------------------------------------------------------- + ## Author: Martin Maechler, 20 Feb 2008 and on + + n <- dim(A)[1] + ## do *the* three calls and look at result + P <- balance(A, "P") + + doPerm <- function(A, pp, i1, i2) { + stopifnot(length(pp) == n, dim(A) == c(n,n), + 1 <= i1, i1 <= i2, i2 <= n) + A. <- A + if(i2 < n) { ## The upper part + for(i in n:(i2+1)) { # 'p2' in *reverse* order + ## swap i <-> pp[i] both rows and columns + tt <- A.[,i]; A.[,i] <- A.[,pp[i]]; A.[,pp[i]] <- tt + tt <- A.[i,]; A.[i,] <- A.[pp[i],]; A.[pp[i],] <- tt + } + } + if(i1 > 1) { ## The lower part + for(i in 1:(i1-1)) { # 'p1' in *forward* order + tt <- A.[,i]; A.[,i] <- A.[,pp[i]]; A.[,pp[i]] <- tt + tt <- A.[i,]; A.[i,] <- A.[pp[i],]; A.[pp[i],] <- tt + } + } + A. + } + + checkPerm <- function(P, orig.A) { + didPerm <- ((leftP <- (i1 <- P$i1) != 1L) | + (rightP <- (i2 <- P$i2) != n)) + if(didPerm) { ## *had* permutation -- now check my idea about it + pp <- as.integer(P$scale) + ## Permute A to become P$z : + A. <- doPerm(orig.A, pp = pp, i1=i1, i2=i2) + stopifnot(isTRUE(all.equal(A., P$z, tolerance = 1e-15))) + + ## Now the reverse: Use pp[] and permute A. "back to A": + if(leftP) { ## The lower part + for(i in (i1-1):1) { # 'p1' in *reverse* order + tt <- A.[,i]; A.[,i] <- A.[,pp[i]]; A.[,pp[i]] <- tt + tt <- A.[i,]; A.[i,] <- A.[pp[i],]; A.[pp[i],] <- tt + } + } + if(rightP) { ## The upper part + for(i in (i2+1):n) { # 'p2' in *forward* order + ## swap i <-> pp[i] both rows and columns + tt <- A.[,i]; A.[,i] <- A.[,pp[i]]; A.[,pp[i]] <- tt + tt <- A.[i,]; A.[i,] <- A.[pp[i],]; A.[pp[i],] <- tt + } + } + stopifnot(isTRUE(all.equal(A., orig.A, tolerance = 1e-15))) + } + } + checkPerm(P, orig.A = A) + + S <- balance(P$z, "S")# "S" starting from result of "P" + stopifnot(S$i1 == 1, S$i2 == n) + + ## Now check the scaling + checkScal <- function (d, A1, A2) { + stopifnot(length(d) == n, dim(A1) == dim(A2), dim(A2) == c(n,n)) + + ## A.scaled <- diag(1/d, n) \%*\% A1 \%*\% diag(d, n) + ## more efficiently: + A.scaled <- A1 * (rep(d, each = n) / d) + stopifnot(isTRUE(all.equal(A2, A.scaled, tolerance = 1e-15))) + ## Check the reverse: + S.rescaled <- A2 * (d * rep(1/d, each = n)) + stopifnot(isTRUE(all.equal(A1, S.rescaled, tolerance = 1e-15))) + } + checkScal(d = S$scale, A1 = P$z, A2 = S$z) + + B <- balance(A, "B")# "B" : B[oth] + stopifnot(P$i1 == B$i1, P$i2 == B$i2) + ## now check *both* permutation and scaling + + A.perm <- doPerm(A, pp = as.integer(B$scale), i1=B$i1, i2=B$i2) + ## checkPerm(B, orig.A = A) + + dB <- B$scale + dB[c(if(B$i1 > 1) 1:(B$i1-1), + if(B$i2 < n) (B$i2+1):n)] <- 1 + checkScal(d = dB, A1 = A.perm, A2 = B$z) + + ## return + list(P = P, S = S, B = B, Sz.eq.Bz = isTRUE(all.equal(S$z, B$z))) + } > m4. <- rbind(c(-1,-2, 0, 0), + c( 0, 0,10,11), + c( 0, 0,12, 0), + c( 0,13, 0, 0)) > str(b4. <- balanceTst(m4.)) List of 4 $ P :List of 4 ..$ z : num [1:4, 1:4] -1 0 0 0 -2 0 13 0 0 11 ... ..$ scale: num [1:4] 1 1 1 3 ..$ i1 : int 2 ..$ i2 : int 3 $ S :List of 4 ..$ z : num [1:4, 1:4] -1 0 0 0 -2 0 13 0 0 11 ... ..$ scale: num [1:4] 1 1 1 1 ..$ i1 : int 1 ..$ i2 : int 4 $ B :List of 4 ..$ z : num [1:4, 1:4] -1 0 0 0 -2 0 13 0 0 11 ... ..$ scale: num [1:4] 1 1 1 3 ..$ i1 : int 2 ..$ i2 : int 3 $ Sz.eq.Bz: logi TRUE > ## better (?) example > (m <- matrix(c(0,-1,0,-2,10, rep(0,11)), 4,4)) [,1] [,2] [,3] [,4] [1,] 0 10 0 0 [2,] -1 0 0 0 [3,] 0 0 0 0 [4,] -2 0 0 0 > str(ba <- balanceTst(m)) List of 4 $ P :List of 4 ..$ z : num [1:4, 1:4] 0 0 0 0 0 0 10 0 -2 -1 ... ..$ scale: num [1:4] 3 1 1 3 ..$ i1 : int 2 ..$ i2 : int 3 $ S :List of 4 ..$ z : num [1:4, 1:4] 0 0 0 0 0 0 2.5 0 -2 -4 ... ..$ scale: num [1:4] 1 0.25 1 1 ..$ i1 : int 1 ..$ i2 : int 4 $ B :List of 4 ..$ z : num [1:4, 1:4] 0 0 0 0 0 0 2.5 0 -2 -4 ... ..$ scale: num [1:4] 3 0.25 1 3 ..$ i1 : int 2 ..$ i2 : int 3 $ Sz.eq.Bz: logi TRUE > ## Hmm: here S$z *differs* from B$z > ## --- but at least, the scale[] and z[] returned seem ok > > > ## a non-empty ``less-balanced'' example --- > > m4 <- matrix(outer(2^(0:7),c(-1,1)), 4,4) > m4[lower.tri(m4)] <- 0 #--> upper triangular ==> will have many permutations > ## now permute it; so balance() will find the permutation > p <- c(4,2:1,3); m4 <- m4[p,p] > m4 [,1] [,2] [,3] [,4] [1,] 128 0 0 0 [2,] 32 -32 0 2 [3,] 16 -16 -1 1 [4,] 64 0 0 4 > str(dm4 <- balanceTst(m4)) # much permutation! i1 = i2 = 1 ! List of 4 $ P :List of 4 ..$ z : num [1:4, 1:4] -1 0 0 0 -16 -32 0 0 1 2 ... ..$ scale: num [1:4] 1 2 1 1 ..$ i1 : int 1 ..$ i2 : int 1 $ S :List of 4 ..$ z : num [1:4, 1:4] -1 0 0 0 -1 -32 0 0 0.25 8 ... ..$ scale: num [1:4] 16 1 4 1 ..$ i1 : int 1 ..$ i2 : int 4 $ B :List of 4 ..$ z : num [1:4, 1:4] -1 0 0 0 -16 -32 0 0 1 2 ... ..$ scale: num [1:4] 1 2 1 1 ..$ i1 : int 1 ..$ i2 : int 1 $ Sz.eq.Bz: logi FALSE > > dm4. <- dgebal(m4) Warning message: 'dgebal' is deprecated. Use 'balance' instead. See help("Deprecated") > storage.mode(m4) <- "integer" > stopifnot(identical(dm4., dgebal(m4))) Warning message: 'dgebal' is deprecated. Use 'balance' instead. See help("Deprecated") > > expm(m) [,1] [,2] [,3] [,4] [1,] -0.999786073 -0.06540707 0 0 [2,] 0.006540707 -0.99978607 0 0 [3,] 0.000000000 0.00000000 1 0 [4,] 0.013081414 -3.99957215 0 1 > expm(m,"Pade") ## are different indeed {when bug still existed} [,1] [,2] [,3] [,4] [1,] -0.999786073 -0.06540707 0 0 [2,] 0.006540707 -0.99978607 0 0 [3,] 0.000000000 0.00000000 1 0 [4,] 0.013081414 -3.99957215 0 1 attr(,"accuracy") [1] 7.918807e-14 > expm(m,"R_Pade")# same as Pade [,1] [,2] [,3] [,4] [1,] -0.999786073 -0.06540707 0 0 [2,] 0.006540707 -0.99978607 0 0 [3,] 0.000000000 0.00000000 1 0 [4,] 0.013081414 -3.99957215 0 1 > > > ## a non-empty ``non-balanced'' example --- > > expm.t.identity(m4, "Ward") [1] TRUE > > m6 <- zeroTrace(matrix(outer(2^(-8:9),c(-1,1)), 6,6)); m6 [,1] [,2] [,3] [,4] [,5] [,6] [1,] -75.2584635 -0.25000 -16.0000 0.00390625 0.25000 16.0000 [2,] -0.0078125 -75.75456 -32.0000 0.00781250 0.50000 32.0000 [3,] -0.0156250 -1.00000 -139.2546 0.01562500 1.00000 64.0000 [4,] -0.0312500 -2.00000 -128.0000 -75.22330729 2.00000 128.0000 [5,] -0.0625000 -4.00000 -256.0000 0.06250000 -71.25456 256.0000 [6,] -0.1250000 -8.00000 -512.0000 0.12500000 8.00000 436.7454 > m6[lower.tri(m6)] <- 0 ## plus one non-zero > m6[4,2] <- 77 > p <- c(6,4,5,2:1,3); m6 <- m6[p,p] > expm.t.identity(m6, "Ward") ## difference; indeed [1] TRUE > expm(m6) # is very different from [,1] [,2] [,3] [,4] [,5] [1,] 4.743903e+189 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [2,] 1.235301e+189 2.158311e-33 7.038631e-32 1.285266e-31 0.000000e+00 [3,] 2.390628e+189 0.000000e+00 1.133807e-31 0.000000e+00 0.000000e+00 [4,] 2.653829e+188 1.627057e-37 7.493870e-34 1.271561e-33 0.000000e+00 [5,] 1.326905e+188 8.096068e-38 3.743330e-34 -4.016897e-34 2.068543e-33 [6,] 5.312842e+188 5.266261e-37 1.683517e-33 3.098706e-35 0.000000e+00 [,6] [1,] 0.000000e+00 [2,] -6.346149e-32 [3,] 0.000000e+00 [4,] -6.406872e-34 [5,] -3.187638e-34 [6,] -1.529642e-35 > expm(m6,"R_Pade") [,1] [,2] [,3] [,4] [,5] [1,] 4.743903e+189 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [2,] 1.235301e+189 2.158311e-33 7.038631e-32 1.285266e-31 0.000000e+00 [3,] 2.390628e+189 0.000000e+00 1.133807e-31 0.000000e+00 0.000000e+00 [4,] 2.653829e+188 1.627057e-37 7.493870e-34 1.271561e-33 0.000000e+00 [5,] 1.326905e+188 8.096068e-38 3.743330e-34 -4.016897e-34 2.068543e-33 [6,] 5.312842e+188 5.266261e-37 1.683517e-33 3.098706e-35 0.000000e+00 [,6] [1,] 0.000000e+00 [2,] -6.346149e-32 [3,] 0.000000e+00 [4,] -6.406872e-34 [5,] -3.187638e-34 [6,] -1.529642e-35 > > str(dm6 <- balanceTst(m6)) List of 4 $ P :List of 4 ..$ z : num [1:6, 1:6] -75.3 0 0 0 0 ... ..$ scale: num [1:6] 3 1 1 1 3 1 ..$ i1 : int 2 ..$ i2 : int 4 $ S :List of 4 ..$ z : num [1:6, 1:6] -75.3 0 0 0 0 ... ..$ scale: num [1:6] 1 1 1 1 2 1 ..$ i1 : int 1 ..$ i2 : int 6 $ B :List of 4 ..$ z : num [1:6, 1:6] -75.3 0 0 0 0 ... ..$ scale: num [1:6] 3 1 1 1 3 1 ..$ i1 : int 2 ..$ i2 : int 4 $ Sz.eq.Bz: logi FALSE > ## Now, that's interesting: > ## > ## 1. 'S' scales *more* (2 .. 5) than just (2:4 == i1:i2) ! > ## > ## 2. 'B' has quite different scaling and it does (must!) obey rule > ## scale i1:i2 only > ## > ## 3. 'B'(oth) is better than "P" and "S" separately: > ## > kappa(eigen(m6)$vectors)# 597.5588 [1] 597.5588 > kappa(eigen(dm6$P$z)$vectors)# 597.5588 [1] 597.5588 > kappa(eigen(dm6$S$z)$vectors)# 42.58396 [1] 621.5746 > kappa(eigen(dm6$B$z)$vectors)# 22.20266 [1] 597.5588 > > > ## An n=17 example where octave's expm() is wrong too > m17 <- matrix(c(10,0, 0, 2, 3,-1, 0, 0, 0, 0, 0, 4, 0, 5, 0, 0,-2, + 0, 0, 0, 0,-3, 0, 0, 0, 0, 0, 0, 0, 0, 6, 0, 7, 0, + 0, 0,10, 0, 0,-4, 9, 0, 0, 0,-5, 0,-6, 0, 0, 0, 0, + 0, 0,-7, 0, 0, 0, 0, 0, 0,10, 0, 0, 0, 0, 0,11, 0, + 0, 0, 0, 0, 0, 0,12, 0, 0, 0, 0, 0,-8, 0, 0, 0, 0, + 0, 0,-9, 0, 0, 0, 0, 0, 0,-10,0,13,14,-11,-12,-13, 0, + 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0, 0,10, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0,-14,16,0,-10,0,17, 0, 0, 0, 0, 0, 0, + 0, 0,-16,0, 0,18,19, 0, 0, 0, 0, 0, 0, 0,20, 0, 21, + 22,0, 0, 0, 0, 0,-17,0, 0, 0,-10,-19,-20,0,0,0, 0, + 0,-21,0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0,23,24, 0,25,26, 0, 0,27,-22,0,28,-23,0,-24, + 0,-25,0,29, 0, 0, 0, 0, 0, 0, 0,30,31, 0, 0, 0, 0, + 0, 0,-26,32,0, 0, 0, 0, 0,-27,0,33,34, 0, 0, 0, 0, + 0,-28,-29,0,0, 0,35, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, 0,36,37, 0, 0, 0, 0, 0, 0, 0, 0,-10), + 17, 17) > str(dm17 <- balanceTst(m17)) List of 4 $ P :List of 4 ..$ z : num [1:17, 1:17] 0 0 0 0 0 0 0 0 0 0 ... ..$ scale: num [1:17] 7 8 9 1 1 1 1 1 1 1 ... ..$ i1 : int 4 ..$ i2 : int 16 $ S :List of 4 ..$ z : num [1:17, 1:17] 0 0 0 0 0 0 0 0 0 0 ... ..$ scale: num [1:17] 1 4 1 2 1 1 1 2 2 1 ... ..$ i1 : int 1 ..$ i2 : int 17 $ B :List of 4 ..$ z : num [1:17, 1:17] 0 0 0 0 0 0 0 0 0 0 ... ..$ scale: num [1:17] 7 8 9 2 2 1 1 2 2 1 ... ..$ i1 : int 4 ..$ i2 : int 16 $ Sz.eq.Bz: logi FALSE > sapply(dm17[1:3], `[[`, "scale") P S B [1,] 7 1.0 7 [2,] 8 4.0 8 [3,] 9 1.0 9 [4,] 1 2.0 2 [5,] 1 1.0 2 [6,] 1 1.0 1 [7,] 1 1.0 1 [8,] 1 2.0 2 [9,] 1 2.0 2 [10,] 1 1.0 1 [11,] 1 1.0 1 [12,] 1 2.0 2 [13,] 1 1.0 1 [14,] 1 1.0 1 [15,] 1 1.0 1 [16,] 1 1.0 1 [17,] 9 0.5 9 > > ## The balancing was really rather harmful -- cond(V) *not* improved: > condX <- function(x) kappa(x, exact=TRUE) > condX(eigen(m17)$vectors)# 8.9e16 [1] 3.485397e+16 > condX(eigen(dm17$P$z)$vectors)# 1.37e17 [1] 1.36865e+17 > condX(eigen(dm17$S$z)$vectors)# 1.44e17 [1] 6.451363e+16 > condX(eigen(dm17$B$z)$vectors)# 1.43e17 (very slightly smaller) [1] 1.437207e+17 > > proc.time() user system elapsed 1.67 0.34 2.00