[Rd] Adding a Matrix Exponentiation Operator

[Vincent Goulet]

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'

[snip]

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

> 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?

The rationale is: you seldom *really* need to inverse a matrix, so we  
won't help you go that route. If you *really* need the explicit  
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?

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

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

HTH

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