balance.RdBalance a square matrix via LAPACK's DGEBAL.
This is an R interface, mainly used for experimentation.
This LAPACK routine is used internally for Eigenvalue decompositions, but also, in Ward(1977)'s algorithm for the matrix exponential.
The name balance() is preferred nowadays, and “dgebal()”
has been deprecated (finally, after 9 years ...).
balance(A, job = c("B", "N", "P", "S"))
## Deprecated now:
## dgebal(A, job = c("B", "N", "P", "S"))A list with components
the transformation of matrix A, after permutation and
or scaling.
numeric vector of length \(n\), containing the permutation and/or scale vectors applied.
integers (length 1) in \(\{1,2,\dots,n\}\), denoted by
ILO and IHI respectively in the LAPACK
documentation. Only relevant for "P" or "B", they
describe where permutations and where scaling took place; see the
‘Details’ section.
An excerpt of the LAPACK documentation about DGEBAL() or
ZGEBAL(), respectively, describing the result
(output) integer
(output) integeri1 and i2 are set to integers such that on exit
z[i,j] = 0 if \(i > j\) and \(j = 1,...,i1-1\) or \(i = i2+1,...,n\).
If job = 'N' or 'S', i1 = 1 and i2 = n.
(output) numeric vector of length n.
Details of the permutations and scaling factors applied to
A. If P[j] is the index of the row and column interchanged
with row and column j and D[j] is the scaling factor
applied to row and column j, then
scale[j] = P[j] for \(j = 1,...,i1-1\) = D[j] for \(j = i1,...,i2\), = P[j] for \(j = i2+1,...,n\).
The order in which the interchanges are made is n to i2+1,
then 1 to i1-1.
Look at the LAPACK documentation for more details.
LAPACK Reference Manual, https://netlib.org/lapack/, balancing ‘gebal’, currently at https://www.netlib.org/lapack/explore-html/df/df3/group__gebal.html.
m4 <- rbind(c(-1,-1, 0, 0),
c( 0, 0,10,10),
c( 0, 0,10, 0),
c( 0,10, 0, 0))
(b4 <- balance(m4))
#> $z
#> [,1] [,2] [,3] [,4]
#> [1,] -1 -1 0 0
#> [2,] 0 0 10 10
#> [3,] 0 10 0 0
#> [4,] 0 0 0 10
#>
#> $scale
#> [1] 1 1 1 3
#>
#> $i1
#> [1] 2
#>
#> $i2
#> [1] 3
#>
## --- for testing and didactical reasons : ----
if(expm:::doExtras()) withAutoprint({
sessionInfo()
packageDescription("Matrix")
"expm installed at"
dirname(attr(packageDescription("expm"), "file"))
})
demo(balanceTst) # also defines the balanceTst() function
#>
#>
#> 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))
#>
#> > op <- options(str = strOptions(vec.len = 12))
#>
#> > 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 0 0 0 10 0 12
#> ..$ 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 0 0 0 10 0 12
#> ..$ 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 0 0 0 10 0 12
#> ..$ scale: num [1:4] 1 1 1 3
#> ..$ i1 : int 2
#> ..$ i2 : int 3
#> $ Sz.eq.Bz: logi TRUE
#>
#> > with(b4., all.equal(P, B)) # TRUE (everywhere?)
#> [1] 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 0 0 0 0 0 0
#> ..$ 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 0 0 0 0 0 0
#> ..$ 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 0 0 0 0 0 0
#> ..$ scale: num [1:4] 3 0.25 1 3
#> ..$ i1 : int 2
#> ..$ i2 : int 3
#> $ Sz.eq.Bz: logi TRUE
#>
#> > (eq <- with(ba, all.equal(S$z, B$z))) # TRUE now (everywhere?)
#> [1] TRUE
#>
#> > ba$Sz.eq.Bz # ditto
#> [1] TRUE
#>
#> > ## 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 4 0 16 32 64 128
#> ..$ 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 4 0 1 32 16 128
#> ..$ 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 4 0 16 32 64 128
#> ..$ scale: num [1:4] 1 2 1 1
#> ..$ i1 : int 1
#> ..$ i2 : int 1
#> $ Sz.eq.Bz: logi FALSE
#>
#> > ##----------- Complex examples
#> > zba4 <- balanceTst(m4 + 3i * m4)
#>
#> > str(zba4)
#> List of 4
#> $ P :List of 4
#> ..$ z : cplx [1:4, 1:4] -1-3i 0+0i 0+0i 0+0i -16-48i -32-96i 0+0i 0+0i 1+3i ...
#> ..$ scale: num [1:4] 1 2 1 1
#> ..$ i1 : int 1
#> ..$ i2 : int 1
#> $ S :List of 4
#> ..$ z : cplx [1:4, 1:4] -1-3i 0+0i 0+0i 0+0i -1-3i -32-96i 0+0i 0+0i 0.25+0.75i ...
#> ..$ scale: num [1:4] 16 1 4 1
#> ..$ i1 : int 1
#> ..$ i2 : int 4
#> $ B :List of 4
#> ..$ z : cplx [1:4, 1:4] -1-3i 0+0i 0+0i 0+0i -16-48i -32-96i 0+0i 0+0i 1+3i ...
#> ..$ scale: num [1:4] 1 2 1 1
#> ..$ i1 : int 1
#> ..$ i2 : int 1
#> $ Sz.eq.Bz: logi FALSE
#>
#> > zba <- balanceTst(m*(1 + 1i))
#>
#> > str(zba)
#> List of 4
#> $ P :List of 4
#> ..$ z : cplx [1:4, 1:4] 0+0i 0+0i 0+0i 0+0i 0+0i 0+0i 10+10i 0+0i -2-2i ...
#> ..$ scale: num [1:4] 3 1 1 3
#> ..$ i1 : int 2
#> ..$ i2 : int 3
#> $ S :List of 4
#> ..$ z : cplx [1:4, 1:4] 0+0i 0+0i 0+0i 0+0i 0+0i 0+0i 2.5+2.5i 0+0i -2-2i ...
#> ..$ scale: num [1:4] 1 0.25 1 1
#> ..$ i1 : int 1
#> ..$ i2 : int 4
#> $ B :List of 4
#> ..$ z : cplx [1:4, 1:4] 0+0i 0+0i 0+0i 0+0i 0+0i 0+0i 2.5+2.5i 0+0i -2-2i ...
#> ..$ scale: num [1:4] 3 0.25 1 3
#> ..$ i1 : int 2
#> ..$ i2 : int 3
#> $ Sz.eq.Bz: logi TRUE
#>
#> > stopifnot(exprs = {
#> + all.equal(ba$ S$z, Re(zba$ S$z))
#> + all.equal(ba$ S$z, Im(zba$ S$z))
#> + all.equal(dm4$ S$z, Re(zba4$ S$z))
#> + all.equal(dm4$ S$z * 3, Im(zba4$ S$z))
#> + })
#>
#> > options(op) # revert
# which in its tests ``defines'' what
# the return value means, notably (i1,i2,scale)