# balance: Balance a Square Matrix via LAPACK's DGEBAL In expm: Matrix Exponential, Log, 'etc'

## Description

Balance 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, where “dgebal” will probably become deprecated.

## Usage

 ```1 2``` ```balance(A, job = c("B", "N", "P", "S")) dgebal(A, job = c("B", "N", "P", "S")) ```

## Arguments

 `A` a square (n x n) numeric matrix. `job` a one-letter string specifying the ‘job’ for DGEBAL. PPermutation SScaling BBoth permutation and scaling NNone

## Details

An excerpt of the LAPACK documentation about DGEBAL(), describing the result

i1 ("ILO")

(output) integer

i2 ("IHI")

(output) integer
`i1` 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`.

scale

(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.

## Value

A list with components

 `z` the transformation of matrix `A`, after permutation and or scaling. `scale` numeric vector of length n, containing the permutation and/or scale vectors applied. `i1,i2` integers (length 1) in \{1,2,…,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.

Martin Maechler

## References

LAPACK Reference Manual

`eigen`, `expm`.

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11``` ```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)) ## --- for testing and didactical reasons : ---- demo(balanceTst) # also defines the balanceTst() function # which in its tests ``defines'' what # the return value means, notably (i1,i2,scale) ```

### Example output

```Loading required package: Matrix

Attaching package: 'expm'

The following object is masked from 'package:Matrix':

expm

\$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

demo(balanceTst)
---- ~~~~~~~~~~

> dgebalTst <- function(A) {
+
+     ## Purpose: Consistency checking of	 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 <- dgebal(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 <- dgebal(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 <- dgebal(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. <- dgebalTst(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 <- dgebalTst(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 -1 -4 ...
..\$ scale: num [1:4] 1 0.125 0.5 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 FALSE

> ## 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 dgebal() 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 <- dgebalTst(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 -32 -32 0 0 4 4 ...
..\$ scale: num [1:4] 1 2 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
```

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