# Adding a Matrix Exponentiation Operator Classic List Threaded 11 messages Open this post in threaded view
|

## Adding a Matrix Exponentiation Operator

 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 /* 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);     ncols = INTEGER(dims);     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) = nrows;     INTEGER(dims2) = ncols;     setAttrib(matrix, R_DimSymbol, dims2);     UNPROTECT(5);     return matrix; }         [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Open this post in threaded view
|

## Re: Adding a Matrix Exponentiation Operator

 2008/4/5, Rory Winston <[hidden email]>: > >  /* 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 ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Open this post in threaded view
|

## Re: Adding a Matrix Exponentiation Operator

Open this post in threaded view
|

## Re: Adding a Matrix Exponentiation Operator

 In reply to this post by Antonio, Fabio Di Narzo Thanks Antonio, thats a good suggestion. On Sat, Apr 5, 2008 at 5:51 PM, Antonio, Fabio Di Narzo < [hidden email]> wrote: > 2008/4/5, Rory Winston <[hidden email]>: > > > > >  /* 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 >         [[alternative HTML version deleted]] ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Open this post in threaded view
|

## Re: Adding a Matrix Exponentiation Operator

Open this post in threaded view
|

## Re: Adding a Matrix Exponentiation Operator

Open this post in threaded view
|

## Re: Adding a Matrix Exponentiation Operator

Open this post in threaded view
|

## Re: Adding a Matrix Exponentiation Operator

Open this post in threaded view
|