1##------OVERVIEW----------------------------------------------------------------
2
3## Input:  A; nxn Matrix, no eigenvalues <=0, not singular
4## Output: log(A); Matrixlogarithm; nxn Matrix
5
6
7## Function for Calculation of log(A) with the Inverse Scaling&Squaring Method
8
9## Step 0:    Schur Decompostion Tr
10## Step 1:    Scaling (root of Tr)
11## Step 2:    Padé-Approximation
12## Step 3:    Squaring
13## Step 4:    Reverse Schur Decomposition
14
15## R-Implementation of Higham's Algorithm from the Book
16## "Functions of Matrices - Theory and Computation", Chapter 11, Algorithm 11.9
17
18##-------CODE-------------------------------------------------------------------
19
20## The coefficients for the Padé-approximation can be computed at install time:
21
22## r: exponents are in  (-51):(-56)
23## p: exponents are in  c((-47):(-53), -56)
24
25
26logm.H08.r <-
27    rbind(c(5003999585967230*2^(-54), 8006399337547537*2^(-54), 5/18, 0,0,0,0),
28          c(5640779706068081*2^(-51), 8899746432686114*2^(-53),
29            8767290225458872*2^(-54), 6733946100265013*2^(-55), 0,0,0),
30          c(5686538473148996*2^(-51), 4670441098084653*2^(-52),
31            5124095576030447*2^(-53), 5604406634440294*2^(-54),
32            8956332917077493*2^(-56), 0,0),
33          c(5712804453675980*2^(-51), 4795663223967718*2^(-52),
34            5535461316768070*2^(-53), 6805310445892841*2^(-54),
35            7824302940658783*2^(-55), 6388318485698934*2^(-56), 0),
36          c(5729264333934497*2^(-51), 4873628951352824*2^(-52),
37            5788422587681293*2^(-53), 7529283295392226*2^(-54),
38            4892742764696865*2^(-54), 5786545115272933*2^(-55),
39            4786997716777457*2^(-56)))
40
41logm.H08.p <-
42    - rbind(c(7992072898328873*2^(-53), 1/2, 8121010851296995*2^(-56), 0,0,0,0),
43            c(8107950463991866*2^(-49), 6823439817291852*2^(-51),
44              6721885580294475*2^(-52), 4839623620596807*2^(-52), 0,0,0),
45            c(6000309411699298*2^(-48),  4878981751356277*2^(-50), 2,
46              5854649940415304*2^(-52), 4725262033344781*2^(-52),0,0),
47            c(8336234321115872*2^(-48), 6646582649377394*2^(-50),
48              5915042177386279*2^(-51), 7271968136730531*2^(-52),
49              5422073417188307*2^(-52), 4660978705505908*2^(-52), 0),
50            c(5530820008925390*2^(-47), 8712075454469181*2^(-50),
51              7579841581383744*2^(-51), 4503599627370617*2^(-51),
52              6406963985981958*2^(-52), 5171999978649488*2^(-52),
53              4621190647118544*2^(-52)))
54
55
56logm.Higham08 <- function(x) {
57    ## work with "Matrix" too: x<-as.matrix(x)
58
59    ##MM: No need to really check here; we get correct error msg later anyway
60    ##    and don't need to compute det() here, in the good cases !
61    ##    if (det(x) == 0) stop("'x' is singular")
62
63    ##-------Step 0: Schur Decomposition-----------------------------------------
64
65    ## Schur() checks for square matrix also:
66    Sch.x <- Schur(Matrix(x, sparse=FALSE))
67    ## FIXME 'sparse=FALSE' is workaround - good as long Matrix has no sparse Schur()
68    ev <- Sch.x@EValues
69    if(getOption("verbose") && any(abs(Arg(ev) - pi) < 1e-7))
70## Let's see what works: temporarily *NOT* stop()ping :
71	message(sprintf("'x' has negative real eigenvalues; maybe ok for %s", "logm()"))
72    n <- Sch.x@Dim[1]
73    Tr <- as.matrix(Sch.x@T)
74    Q  <- as.matrix(Sch.x@Q)
75
76    ##----- Step 1: [Inverse] Scaling -------------------------------------------
77    I <- diag(n)
78    thMax <- 0.264
79    theta <- c(0.0162, 0.0539, 0.114, 0.187, thMax)
80    p <- k <- 0 ; t.o <- -1
81    ## NB: The following could loop forever, e.g., for  logm(Diagonal(x=1:0))
82    repeat{
83        t <- norm(Tr - I, "1") # norm(x, .) : currently x is coerced to dgeMatrix
84	if(is.na(t)) {
85	    warning(sprintf(ngettext(k,
86				     "NA/NaN from %s after %d step.\n",
87				     "NA/NaN from %s after %d steps.\n"),
88			    " || Tr - I || ", k),
89		    "The matrix logarithm may not exist for this matrix.")
90	    return(array(t, dim=dim(Tr)))
91	}
92        if (t < thMax) {
93            ## FIXME: use findInterval()
94            j2 <- which.max( t <= theta)
95            j1 <- which.max( (t/2) <= theta)
96            if ((j2-j1 <= 1) || ((p <- p+1) == 2)) {
97                m <- j2 ## m := order of the Padé-approximation
98                break
99            }
100        } else if(k > 20 && abs(t.o - t) < 1e-7*t) {
101            ##
102	    warning(sprintf("Inverse scaling did not work (t = %g).\n", t),
103		    "The matrix logarithm may not exist for this matrix.",
104		    "Setting m = 3 arbitrarily.")
105            m <- 3
106            break
107        }
108        Tr <- rootS(Tr)##--> Matrix Square root of Jordan T
109        ##    -----    [see below;  compare with ./sqrtm.R
110        t.o <- t
111        k <- k+1
112    }
113    if(getOption("verbose"))
114	message(sprintf("logm.Higham08() -> (k, m) = (%d, %d)", k,m))
115
116    ##------ Step 2: Padé-Approximation -----------------------------------------
117
118    ## of order m :
119    r.m <- logm.H08.r[m,]
120    p.m <- logm.H08.p[m,]
121    X <- 0
122    Tr <- Tr-I
123    for (s in 1:(m+2)) {
124        X <- X + r.m[s]*solve(Tr - p.m[s]*I, Tr)
125    }
126
127    ##--- Step 3 & 4: Squaring & reverse Schur Decomposition -----------------
128
129    2^k* Q %*% X %*% solve(Q)
130}
131
132### --- was  rootS.r -----------------------------------------------------------
133###          ~~~~~~~
134
135##------OVERVIEW----------------------------------------------------------------
136
137## Input:  UT; nxn upper triangular block matrix (real Schur decomposition)
138## Output: root of matrix UT, nxn upper triangular Matrix
139
140## Function for calculation of UT^(1/2), which is used for the logarithm function
141
142## Step 0:    Analyse block structure
143## Step 1:    Calculate diagonal elements/blocks
144## Step 2:    Calculate superdiagonal elements/blocks
145
146## R-Implementation of Higham's Algorithm from the Book
147## "Functions of Matrices - Theory and Computation", Chapter 6, Algorithm 6.7
148
149## NB: Much in parallel with sqrtm() in ./sqrtm.R  <<< keep in sync
150##                           ~~~~~        ~~~~~~~
151rootS <- function(x) {
152    ## Generate Basic informations of Matrix x
153    stopifnot(length(d <- dim(x)) == 2, is.numeric(d),
154              (n <- d[1]) == d[2], n >= 1)
155    ## FIXME : should work for "Matrix" too: not     S <- as.matrix(x)
156    S <- x
157
158    ##------- STEP 0: Analyse block structure ----------------------------------
159    if(n > 1L) {
160        ## Count 2x2 blocks (as Schur(x) is the real Schur Decompostion)
161        J.has.2 <- S[cbind(2:n, 1:(n-1))] != 0
162        k <- sum(J.has.2) ## := number of non-zero SUB-diagonals
163    } else k <- 0L
164
165    ## Generate Blockstructure and save it as R.index
166    R.index <- vector("list",n-k)
167    l <- 1L
168    i <- 1L
169    while(i < n) { ## i advances by 1 or 2, depending on 1- or 2- Jordan Block
170	if (S[i+1L,i] == 0) {
171	    R.index[[l]] <- i
172	}
173	else {
174            i1 <- i+1L
175	    R.index[[l]] <- c(i,i1) # = i:(i+1)
176	    i <- i1
177	}
178	i <- i+1L
179	l <- l+1L
180    }
181    if (is.null(R.index[[n-k]])) { # needed; FIXME: should be able to "know"
182        ##message(sprintf("R.index[n-k = %d]] is NULL, set to n=%d", n-k,n))
183	R.index[[n-k]] <- n
184    }
185
186    ##---------STEP 1: Calculate diagonal elements/blocks------------------------
187    ## Calculate the root of the diagonal blocks of the Schur Decompostion S
188    I <- diag(2)
189    X <- matrix(0,n,n)
190    for (j in seq_len(n-k)) {
191        ij <- R.index[[j]]
192	if (length(ij) == 1L) {
193	    ## Sij <- S[ij,ij]
194	    ## if(Sij < 0)
195	    ##	   ## FIXME(?) : in sqrtm(), we take *complex* sqrt() if needed :
196	    ##	   ## -----  but afterwards  norm(Tr - I, "1") fails with complex
197	    ##	   ## Sij <- complex(real = Sij, imaginary = 0)
198	    ##	   stop("negative diagonal entry -- matrix square does not exist")
199	    ## X[ij,ij] <- sqrt(Sij)
200	    X[ij,ij] <- sqrt(S[ij,ij])
201	}
202        else {
203            ## "FIXME"(better algorithm): only need largest eigen value
204            ev1 <- eigen(S[ij,ij], only.values=TRUE)$values[1]
205            r1 <- Re(sqrt(ev1)) ## sqrt(<complex>) ...
206            X[ij,ij] <-
207                r1*I + 1/(2*r1)*(S[ij,ij] - Re(ev1)*I)
208        }
209    }
210
211### ___ FIXME __ code re-use: All the following is identical to 'STEP 3' in sqrtm()
212###     -----    and almost all of STEP 1 above is == 'STEP 2' of sqrtm()
213
214    ##---------STEP 2: Calculate superdiagonal elements/blocks-------------------
215
216    ## Calculate the remaining, not-diagonal blocks
217    if (n-k > 1L) for (j in 2L:(n-k)) {
218        ij <- R.index[[j]]
219        for (i in (j-1L):1L) {
220            ii <- R.index[[i]]
221            sumU <- 0
222
223            ## Calculation for 1x1 Blocks
224	    if (length(ij) == 1L & length(ii) == 1L ) {
225		if (j-i > 1L) for (l in (i+1L):(j-1L)) {
226		    il <- R.index[[l]]
227		    sumU <- sumU + {
228			if (length(il) == 2 ) X[ii,il]%*%X[il,ij]
229			else		      X[ii,il] * X[il,ij]
230		    }
231		}
232		X[ii,ij] <- solve(X[ii,ii]+X[ij,ij],S[ii,ij]-sumU)
233	    }
234
235	    ## Calculation for	1x2 Blocks
236	    else if (length(ij) == 2 & length(ii) == 1L ) {
237		if (j-i > 1L) for (l in(i+1L):(j-1L)) {
238		    il <- R.index[[l]]
239		    sumU <- sumU + {
240			if (length(il) == 2) X[ii,il]%*%X[il,ij]
241			else		     X[ii,il] * X[il,ij]
242		    }
243		}
244		X[ii,ij] <- solve(t(X[ii,ii]*I + X[ij,ij]),
245				  as.vector(S[ii,ij] - sumU))
246	    }
247
248	    ## Calculation for	2x1 Blocks
249	    else if (length(ij) == 1L & length(ii) == 2 ) {
250		if (j-i > 1L) for (l in(i+1L):(j-1L)) {
251		    il <- R.index[[l]]
252		    sumU <- sumU + {
253			if (length(il) == 2 ) X[ii,il]%*%X[il,ij]
254			else		      X[ii,il] * X[il,ij]
255		    }
256		}
257		X[ii,ij] <- solve(X[ii,ii]+X[ij,ij]*I,S[ii,ij]-sumU)
258	    }
259	    ## Calculation for	2x2 Blocks with special equation for solver
260	    else if (length(ij) == 2 & length(ii) == 2 ) {
261		if (j-i > 1L) for (l in(i+1L):(j-1L)) {
262		    il <- R.index[[l]]
263		    sumU <- sumU + {
264			if (length(il) == 2 ) X[ii,il] %*%   X[il,ij]
265			else		      X[ii,il] %*% t(X[il,ij])
266		    }
267		}
268		tUii <- matrix(0,4,4)
269		tUii[1:2,1:2] <- X[ii,ii]
270		tUii[3:4,3:4] <- X[ii,ii]
271		tUjj <- matrix(0,4,4)
272		tUjj[1:2,1:2] <- t(X[ij,ij])[1L,1L]*I
273		tUjj[3:4,3:4] <- t(X[ij,ij])[2L,2L]*I
274		tUjj[1:2,3:4] <- t(X[ij,ij])[1L,2L]*I
275		tUjj[3:4,1:2] <- t(X[ij,ij])[2L,1L]*I
276		X[ii,ij] <- solve(tUii+tUjj,as.vector(S[ii,ij]-sumU))
277	    }
278	}
279    }
280    X
281}
282
283