# Re: [Rd] Adding a Matrix Exponentiation Operator

Date: Fri, 11 Apr 2008 10:33:13 -0400

Dear Martin et al.,

Paul Gilbert forwarded me your discussions on this topic.

I have been interested in computing f(A), where A is any square matrix and f(.) is an analytic function whose Taylor-McLaurin series has a finite radius of convergence. There exist special algorithms when f(.) is square-root, exponential, logarithm (note: the TM series for this doesn't have a finite-radius of convergence), sine, and cosine. The easiest way to handle a general f(.) is via spectral decomposition, but that is limited to diagonalizable matrices, and is also ill-conditioned. A paper by Davies and Higham (SIAM 2003) provide an algorithm, the Schur-Parlett method, for doing this. I think this is the state-of-the-art for evaluating f(A). Higham also has a recent SIAM book (2008) on computing f(A). There are matlab codes provided with that, but it doesn't have the algorithm discussed in his (2003) paper. I asked him for the code, but he was unable to share it with me. I suspect this may be due to proprietary issues since this is implemented in the latest Matlab release. However, we should be able to do this in R based on the algorithmic description in is (2003) paper.

Would this be of interest to the R group?

Best,
Ravi.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

>>>>> "RW" == Rory Winston <rory.winston@gmail.com>
>>>>> on Sat, 5 Apr 2008 14:44:44 +0100 writes:

```    RW> Hi all I recently started to write a matrix
RW> exponentiation operator for R (by adding a new operator
RW> definition to names.c, and adding the following code to
RW> arrays.c). It is not finished yet, but I would like to
RW> solicit some comments, as there are a few areas of R's
RW> internals that I am still feeling my way around.

```

RW> Firstly:

```    RW> 1) Would there be interest in adding a new operator %^%
RW> that performs the matrix equivalent of the scalar ^
RW> operator?

```

Yes. A few weeks ago, I had investigated a bit about this and found several R-help/R-devel postings with code proposals and then about half dozen CRAN packages with diverse implementations of the matrix power (I say "power" very much on purpose, in order to not confuse it with the matrix exponential which is another much more interesting topic, also recently (+- two months?) talked about.

Consequently I made a few timing tests and found that indeed, the "smart matrix power" {computing m^2, m^4, ... and only those multiplications needed} as you find it in many good books about algorithms and e.g. also in *the* standard Golub & Van Loan "Matrix Computation" is better than "the eigen" method even for large powers.

matPower <- function(X,n)
## Function to calculate the n-th power of a matrix X {

if(n != round(n)) {

```        n <- round(n)
warning("rounding exponent `n' to", n)
```
}
if(n == 0)

return(diag(nrow = nrow(X)))
n <- n - 1
phi <- X
## pot <- X # the first power of the matrix.     while (n > 0)
{

```        if (n %% 2)
phi <- phi %*% X
if (n == 1) break
n <- n %/% 2
X <- X %*% X
```

}
return(phi)
}

"Simultaneously" people where looking at the matrix exponential expm() in the Matrix package,
and some of us had consequently started the 'expm' project on R-forge.
The main goal there has been to investigate several algorithms for the matrix exponential, notably the one buggy implementation

```(in the 'Matrix' package until a couple of weeks ago, the bug
stemming from 'octave' implementation).
The authors of 'actuar', Vincent and Christophe, notably also
```
had code for the matrix *power* in a C (building on BLAS) and I had added an R-interface '%^%' there as well.

Yes, with the goal to move that (not the matrix exponential yet) into standard R.
Even though it's not used so often (in percentage of all uses of R), it's simple to *right*, and I have seen very many versions of the matrix power that were much slower / inaccurate / ... such that a reference implementation seems to be called for.

install.packages("expm",repos="http://R-Forge.R-project.org")

• but only from tomorrow for Windows (which installs a pre-compiled package), since I found that we had accidentally broken the package trivially by small changes two weeks ago.

and then

library(expm)
?%^%

Best regards,
Martin Maechler, ETH Zurich

```    RW> operator? I am implicitly assuming that the benefits of
RW> a native exponentiation routine for Markov chain
RW> evaluation or function generation would outstrip that of
RW> an R solution. Based on my tests so far (comparing it to
RW> a couple of different pure R versions) it does, but I
RW> still there is much room for optimization in my method.
RW> 2) Regarding the code below: Is there a better way to do
RW> the matrix multiplication? I am creating quite a few
RW> copies for this implementation of exponentiation by
RW> squaring. Is there a way to cut down on the number of
RW> copies I am making here (I am assuming that the lhs and
RW> rhs of matprod() must be different instances).

```

RW> Any feedback appreciated ! Thanks Rory

RW> <snip>

```    RW> /* Convenience function */ static void
RW> copyMatrixData(SEXP a, SEXP b, int nrows, int ncols, int
RW> mode) { for (int i=0; i < ncols; ++i) for (int j=0; j <
RW> nrows; ++j) REAL(b)[i * nrows + j] = REAL(a)[i * nrows +
RW> j]; }

RW> SEXP do_matexp(SEXP call, SEXP op, SEXP args, SEXP rho)
```
RW> { int nrows, ncols; SEXP matrix, tmp, dims, dims2; SEXP     RW> x, y, x_, x__; int i,j,e,mode;

RW> // Still need to fix full complex support mode =     RW> isComplex(CAR(args)) ? CPLXSXP : REALSXP;

RW> SETCAR(args, coerceVector(CAR(args), mode)); x =     RW> CAR(args); y = CADR(args);

RW> dims = getAttrib(x, R_DimSymbol); nrows =     RW> INTEGER(dims)[0]; ncols = INTEGER(dims)[1];

RW> if (nrows != ncols) error(_("can only raise square     RW> matrix to power"));

RW> if (!isNumeric(y)) error(_("exponent must be a     RW> scalar integer"));

RW> e = asInteger(y);

RW> if (e < -1) error(_("exponent must be >= -1")); else     RW> if (e == 1) return x;

```    RW>     else if (e == -1) { /* return matrix inverse via
RW> solve() */ SEXP p1, p2, inv; PROTECT(p1 = p2 =
RW> allocList(2)); SET_TYPEOF(p1, LANGSXP); CAR(p2) =
RW> install("solve.default"); p2 = CDR(p2); CAR(p2) = x; inv
RW> = eval(p1, rho); UNPROTECT(1); return inv; }

RW>     PROTECT(matrix = allocVector(mode, nrows * ncols));
RW> PROTECT(tmp = allocVector(mode, nrows * ncols));
```
RW> PROTECT(x_ = allocVector(mode, nrows * ncols));     RW> PROTECT(x__ = allocVector(mode, nrows * ncols));

RW> copyMatrixData(x, x_, nrows, ncols, mode);

```    RW>     // Initialize matrix to identity matrix // Set x[i *
RW> ncols + i] = 1 for (i = 0; i < ncols*nrows; i++)
RW> REAL(matrix)[i] = ((i % (ncols+1) == 0) ? 1 : 0);

RW>     if (e == 0) { ; // return identity matrix } else
RW> while (e > 0) { if (e & 1) { if (mode == REALSXP)
RW> matprod(REAL(matrix), nrows, ncols, REAL(x_), nrows,
```
RW> ncols, REAL(tmp)); else cmatprod(COMPLEX(tmp), nrows,     RW> ncols, COMPLEX(x_), nrows, ncols, COMPLEX(matrix));

RW> copyMatrixData(tmp, matrix, nrows, ncols,     RW> mode); e--; }

```    RW>         if (mode == REALSXP) matprod(REAL(x_), nrows,
RW> ncols, REAL(x_), nrows, ncols, REAL(x__)); else
RW> cmatprod(COMPLEX(x_), nrows, ncols, COMPLEX(x_), nrows,
RW> ncols, COMPLEX(x__));

RW>         copyMatrixData(x__, x_, nrows, ncols, mode); e
```
RW> /= 2; }
```    RW>     PROTECT(dims2 = allocVector(INTSXP, 2));
RW> INTEGER(dims2)[0] = nrows; INTEGER(dims2)[1] = ncols;
RW> setAttrib(matrix, R_DimSymbol, dims2);

RW>     UNPROTECT(5); return matrix; }

RW> 	[[alternative HTML version deleted]]

RW> ______________________________________________
```
RW> R-devel_at_r-project.org mailing list     RW> https://stat.ethz.ch/mailman/listinfo/r-devel

R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Fri 11 Apr 2008 - 16:07:58 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Fri 11 Apr 2008 - 20:32:32 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.