From maechler at stat.math.ethz.ch  Tue Apr  8 15:57:18 2008
From: maechler at stat.math.ethz.ch (Martin Maechler)
Date: Tue, 8 Apr 2008 15:57:18 +0200
Subject: [Rd] Adding a Matrix Exponentiation Operator
In-Reply-To: <2BA1FDBF-4075-42E7-AED9-AD56F091C45D@act.ulaval.ca>
References: <3f446aa30804050644y2fa4cc77l5c22dabc80fe4cec@mail.gmail.com>
	<18423.48201.828727.639538@cmath-5.math.ethz.ch>
	<3f446aa30804060401u34871311h575edf6861cd9285@mail.gmail.com>
	<2BA1FDBF-4075-42E7-AED9-AD56F091C45D@act.ulaval.ca>
Message-ID: <18427.31166.634180.543962@stat.math.ethz.ch>

>>>>> "VG" == Vincent Goulet <vincent.goulet at act.ulaval.ca>
>>>>>     on Tue, 08 Apr 2008 09:28:00 -0400 writes:

    VG> Le dim. 6 avr. ? 07:01, Rory Winston a ?crit :
    >> Hi Martin
    >> 
    >> Thanks for the detailed reply. I had a look at the matrix power
    >> implementation in the actuar package and the modified version in the  
    >> expm
    >> package. I have a couple of questions/comments:
    >> 
    >> 1. Firstly, I seem to have trouble loading expm.
    >> 
    >>> install.packages("expm",repos="http://R-Forge.R-project.org")
    >> trying URL '
    >> http://R-Forge.R-project.org/bin/macosx/universal/contrib/2.6/expm_0.9-1.tgz
    >> '
    >> Content type 'application/x-gzip' length 149679 bytes (146 Kb)
    >> opened URL
    >> ==================================================
    >> downloaded 146 Kb
    >> ...
    >>> library("expm")
    >> Error in namespaceExport(ns, exports) : undefined exports :matpow
    >> Error: package/namespace load failed for 'expm'

    VG> [snip]

    VG> The current version of the package on R-Forge is 0.9-2 and the  
    VG> NAMESPACE issue should be fixed there.

Yes, and I told Rory explicitly that he should wait a day before
downloading the windows version of 'expm' ...
{Siiiggh, it's not only the documentation that people are not reading ...}

    >> 2. On to the package implementation, I see we actually have very  
    >> similar
    >> implementations. The main differences are:
    >> 
    >> i) For an exponent equal to -1, I call back into R for the solve()  
    >> function
    >> using eval() and CAR/CDR'ing the arguments into place. The actuar  
    >> package
    >> calls dgesv() directly. I suspect that the direct route is more  
    >> efficient
    >> and thus the more preferable one. I see that your implementation  
    >> doesnt
    >> calculate the inverse for an exponent of -1,is there a specific  
    >> reason for
    >> doing that?

    VG> The rationale is: you seldom *really* need to inverse a matrix, so we  
    VG> won't help you go that route. If you *really* need the explicit  
    VG> inverse, then use solve() directly (as the error message says).

    >> ii) Regarding complex matrices: I guess we should have support for  
    >> this, as
    >> its not unreasonable that someone may do this, and it should be  
    >> pretty easy
    >> to implement. My function doesnt have full support yet.
    >> 
    >> iii) A philosophical question - where the the "right" place for the  
    >> %^%
    >> operator? Is it in the operator list at a C level along with %*% and  
    >> the
    >> like? Or is it in an R file as a function definition?

    VG> Well... both. The operator %^% is defined at the R level but the  
    VG> computations are done at the C level by function matpow(). Or perhaps  
    VG> I didn't understand what you mean, here.

Thank you, Vincent.
I think Rory was asking about how the integration into "R base"
should happen.  In that case there's much more choice than just
the .Call() or .External() one: There's also .Internal() and ".Primitive",
and for instance %*% is primitive.

One current advantage of having %^% be primitive would be that
it can automatically also be an S4 and S3 generic function.
However there are plans to make this easier for R 2.8.0++ and we
are talking about that version of R anyway.

    >> I dont really have a
    >> preference either way...have you an opinion on this?
    >> 
    >> Thanks
    >> Rory

    VG> HTH

    VG> Vincent


    >> 
    >> 
    >> 
    >> On Sat, Apr 5, 2008 at 6:52 PM, Martin Maechler <maechler at stat.math.ethz.ch 
    >> >
    >> wrote:
    >> 
    >>>>>>>> "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
    >>> 
    >> 
    >> [[alternative HTML version deleted]]
    >> 
    >> ______________________________________________
    >> R-devel at r-project.org mailing list
    >> https://stat.ethz.ch/mailman/listinfo/r-devel

    VG> ______________________________________________
    VG> R-devel at r-project.org mailing list
    VG> https://stat.ethz.ch/mailman/listinfo/r-devel



[Rd] Adding a Matrix Exponentiation Operator
Martin Maechler