Hi all I recently started to write a matrix exponentiation operator for R (by adding a new operator definition to names.c, and adding the following code to arrays.c). It is not finished yet, but I would like to solicit some comments, as there are a few areas of R's internals that I am still feeling my way around. Firstly: 1) Would there be interest in adding a new operator %^% that performs the matrix equivalent of the scalar ^ operator? I am implicitly assuming that the benefits of a native exponentiation routine for Markov chain evaluation or function generation would outstrip that of an R solution. Based on my tests so far (comparing it to a couple of different pure R versions) it does, but I still there is much room for optimization in my method. 2) Regarding the code below: Is there a better way to do the matrix multiplication? I am creating quite a few copies for this implementation of exponentiation by squaring. Is there a way to cut down on the number of copies I am making here (I am assuming that the lhs and rhs of matprod() must be different instances). Any feedback appreciated ! Thanks Rory <snip> /* Convenience function */ static void copyMatrixData(SEXP a, SEXP b, int nrows, int ncols, int mode) { for (int i=0; i < ncols; ++i) for (int j=0; j < nrows; ++j) REAL(b)[i * nrows + j] = REAL(a)[i * nrows + j]; } SEXP do_matexp(SEXP call, SEXP op, SEXP args, SEXP rho) { int nrows, ncols; SEXP matrix, tmp, dims, dims2; SEXP x, y, x_, x__; int i,j,e,mode; // Still need to fix full complex support mode = isComplex(CAR(args)) ? CPLXSXP : REALSXP; SETCAR(args, coerceVector(CAR(args), mode)); x = CAR(args); y = CADR(args); dims = getAttrib(x, R_DimSymbol); nrows = INTEGER(dims)[0]; ncols = INTEGER(dims)[1]; if (nrows != ncols) error(_("can only raise square matrix to power")); if (!isNumeric(y)) error(_("exponent must be a scalar integer")); e = asInteger(y); if (e < -1) error(_("exponent must be >= -1")); else if (e == 1) return x; else if (e == -1) { /* return matrix inverse via solve() */ SEXP p1, p2, inv; PROTECT(p1 = p2 = allocList(2)); SET_TYPEOF(p1, LANGSXP); CAR(p2) = install("solve.default"); p2 = CDR(p2); CAR(p2) = x; inv = eval(p1, rho); UNPROTECT(1); return inv; } PROTECT(matrix = allocVector(mode, nrows * ncols)); PROTECT(tmp = allocVector(mode, nrows * ncols)); PROTECT(x_ = allocVector(mode, nrows * ncols)); PROTECT(x__ = allocVector(mode, nrows * ncols)); copyMatrixData(x, x_, nrows, ncols, mode); // Initialize matrix to identity matrix // Set x[i * ncols + i] = 1 for (i = 0; i < ncols*nrows; i++) REAL(matrix)[i] = ((i % (ncols+1) == 0) ? 1 : 0); if (e == 0) { ; // return identity matrix } else while (e > 0) { if (e & 1) { if (mode == REALSXP) matprod(REAL(matrix), nrows, ncols, REAL(x_), nrows, ncols, REAL(tmp)); else cmatprod(COMPLEX(tmp), nrows, ncols, COMPLEX(x_), nrows, ncols, COMPLEX(matrix)); copyMatrixData(tmp, matrix, nrows, ncols, mode); e--; } if (mode == REALSXP) matprod(REAL(x_), nrows, ncols, REAL(x_), nrows, ncols, REAL(x__)); else cmatprod(COMPLEX(x_), nrows, ncols, COMPLEX(x_), nrows, ncols, COMPLEX(x__)); copyMatrixData(x__, x_, nrows, ncols, mode); e /= 2; } PROTECT(dims2 = allocVector(INTSXP, 2)); INTEGER(dims2)[0] = nrows; INTEGER(dims2)[1] = ncols; setAttrib(matrix, R_DimSymbol, dims2); UNPROTECT(5); return matrix; } [[alternative HTML version deleted]]
Antonio, Fabio Di Narzo
2008-Apr-05 16:51 UTC
[Rd] Adding a Matrix Exponentiation Operator
2008/4/5, Rory Winston <rory.winston at gmail.com>: <snip>> > /* Convenience function */ > static void copyMatrixData(SEXP a, SEXP b, int nrows, int ncols, int mode) { > for (int i=0; i < ncols; ++i) > for (int j=0; j < nrows; ++j) > REAL(b)[i * nrows + j] = REAL(a)[i * nrows + j]; > }I would use 'memcpy' here instead of the double for loop, i.e.: memcpy(REAL(b), REAL(a), length(b) * sizeof(double)) -- Antonio, Fabio Di Narzo Ph.D. student at Department of Statistical Sciences University of Bologna, Italy
>>>>> "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