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'
Possibly a namespace file issue? My version is:
platform i386-apple-darwin8.10.1
arch i386
os darwin8.10.1
system i386, darwin8.10.1
status
major 2
minor 6.1
year 2007
month 11
day 26
svn rev 43537
language R
version.string R version 2.6.1 (2007-11-26)
>
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?
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? I dont really have a
preference either way...have you an opinion on this?
Thanks
Rory
On Sat, Apr 5, 2008 at 6:52 PM, Martin Maechler <[EMAIL PROTECTED]>
wrote:
> >>>>> "RW" == Rory Winston <[EMAIL PROTECTED]>
> >>>>> 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@r-project.org mailing list
> RW> https://stat.ethz.ch/mailman/listinfo/r-devel
>
[[alternative HTML version deleted]]
______________________________________________
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
