[Rd] Adding a Matrix Exponentiation Operator

[Martin Maechler]

>>>>> "RW" == Rory Winston <rory.winston at 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.

So, please consider
    
    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