Adding a Matrix Exponentiation Operator

classic Classic list List threaded Threaded
11 messages Options
Reply | Threaded
Open this post in threaded view
|

Adding a Matrix Exponentiation Operator

Rory Winston
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]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Adding a Matrix Exponentiation Operator

Antonio, Fabio Di Narzo
2008/4/5, Rory Winston <[hidden email]>:
<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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Adding a Matrix Exponentiation Operator

Martin Maechler
In reply to this post by Rory Winston
>>>>> "RW" == Rory Winston <[hidden email]>
>>>>>     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> [hidden email] mailing list
    RW> https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Adding a Matrix Exponentiation Operator

Rory Winston
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]>:
> <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
>

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Adding a Matrix Exponentiation Operator

Rory Winston
In reply to this post by Martin Maechler
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 <[hidden email]>
wrote:

> >>>>> "RW" == Rory Winston <[hidden email]>
> >>>>>     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> [hidden email] mailing list
>    RW> https://stat.ethz.ch/mailman/listinfo/r-devel
>

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Adding a Matrix Exponentiation Operator

Rory Winston
Hi Martin

I have changed the implementation slightly so that it now handles complex
matrices as well. Take a look and see what you think. I realise that a lot
of the if/else mode checking could be merged.

Cheers
Rory

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;

    // necessary?
    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));

    if (mode == REALSXP)
        Memcpy(REAL(x_), REAL(x), (size_t)nrows*ncols);
    else
        Memcpy(COMPLEX(x_), COMPLEX(x), (size_t)nrows*ncols);

    // Initialize matrix to identity matrix
    // Set x[i * ncols + i] = 1
    if (mode == REALSXP)
        for (i = 0; i < ncols*nrows; i++)
            REAL(matrix)[i] = ((i % (ncols+1) == 0) ? 1 : 0);
    else
        for (i = 0; i < ncols*nrows;i++) {
            COMPLEX(matrix)[i].i = 0.0;
            COMPLEX(matrix)[i].r = ((i % (ncols+1) == 0) ? 1.0 : 0.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(matrix), nrows, ncols,
                        COMPLEX(x_), nrows, ncols, COMPLEX(tmp));

            //copyMatrixData(tmp, matrix, nrows, ncols, mode);
            if (mode == REALSXP)
                Memcpy(REAL(matrix), REAL(tmp), (size_t)nrows*ncols);
            else
                Memcpy(COMPLEX(matrix), COMPLEX(tmp), (size_t)nrows*ncols);
            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);
        if (mode == REALSXP)
            Memcpy(REAL(x_), REAL(x__), (size_t)nrows*ncols);
        else
            Memcpy(COMPLEX(x_), COMPLEX(x__), (size_t)nrows*ncols);

        e >>= 1;
    }

    PROTECT(dims2 = allocVector(INTSXP, 2));
    INTEGER(dims2)[0] = nrows;
    INTEGER(dims2)[1] = ncols;
    setAttrib(matrix, R_DimSymbol, dims2);

    UNPROTECT(5);
    return matrix;
}



On Sun, Apr 6, 2008 at 12:01 PM, Rory Winston <[hidden email]>
wrote:

> 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 <
> [hidden email]> wrote:
>
> > >>>>> "RW" == Rory Winston <[hidden email]>
> > >>>>>     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> [hidden email] mailing list
> >    RW> https://stat.ethz.ch/mailman/listinfo/r-devel
> >
>
>

        [[alternative HTML version deleted]]

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Adding a Matrix Exponentiation Operator

Vincent Goulet
In reply to this post by Rory Winston

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 <[hidden email]
> >
> wrote:
>
>>>>>>> "RW" == Rory Winston <[hidden email]>
>>>>>>>    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> [hidden email] mailing list
>>   RW> https://stat.ethz.ch/mailman/listinfo/r-devel
>>
>
> [[alternative HTML version deleted]]
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Adding a Matrix Exponentiation Operator

Martin Maechler
>>>>> "VG" == Vincent Goulet <[hidden email]>
>>>>>     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 <[hidden email]
    >> >
    >> wrote:
    >>
    >>>>>>>> "RW" == Rory Winston <[hidden email]>
    >>>>>>>> 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> [hidden email] mailing list
    RW> https://stat.ethz.ch/mailman/listinfo/r-devel
    >>>
    >>
    >> [[alternative HTML version deleted]]
    >>
    >> ______________________________________________
    >> [hidden email] mailing list
    >> https://stat.ethz.ch/mailman/listinfo/r-devel

    VG> ______________________________________________
    VG> [hidden email] mailing list
    VG> https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Adding a Matrix Exponentiation Operator

Rory Winston
Hi martin
Actually I did read your email. If you had actually read *my* email you would have seen that I am not using the windows version, which is what you explicitly referred to earlier. Lack of observation aside, that was indeed what I was asking.

Thanks
Rory
Sent using BlackBerry® from Orange

-----Original Message-----
From: Martin Maechler <[hidden email]>

Date: Tue, 8 Apr 2008 15:57:18
To:Vincent Goulet <[hidden email]>
Cc:Rory Winston <[hidden email]>, [hidden email]
Subject: Re: [Rd] Adding a Matrix Exponentiation Operator


>>>>> "VG" == Vincent Goulet <[hidden email]>
>>>>>     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 <[hidden email]
    >> >
    >> wrote:
    >>
    >>>>>>>> "RW" == Rory Winston <[hidden email]>
    >>>>>>>> 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> [hidden email] mailing list
    RW> https://stat.ethz.ch/mailman/listinfo/r-devel
    >>>
    >>
    >> [[alternative HTML version deleted]]
    >>
    >> ______________________________________________
    >> [hidden email] mailing list
    >> https://stat.ethz.ch/mailman/listinfo/r-devel

    VG> ______________________________________________
    VG> [hidden email] mailing list
    VG> https://stat.ethz.ch/mailman/listinfo/r-devel
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Adding a Matrix Exponentiation Operator

Martin Maechler
>>>>> "rw" == rory winston <[hidden email]>
>>>>>     on Tue, 8 Apr 2008 14:06:26 +0000 writes:

    rw> Hi martin
    rw> Actually I did read your email. If you had actually read *my* email you would have seen that I am not using the windows version, which is what you explicitly referred to earlier. Lack of observation aside, that was indeed what I was asking.

Ah, ok,
please excuse me.

I only quickly scanned that you had installed the outdate binary
version of the package, and wrongly concluded you were using Windows.

Indeed, you used the MacOS one which for this case suffers from the
identical problem: It installs a *binary* package, the binaries
are only rebuilt from the source once a day (or so), and so you got the
old version instead of the - then current - version.

    rw> Thanks
    rw> Rory

Martin

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Adding a Matrix Exponentiation Operator

Ravi Varadhan
In reply to this post by Martin Maechler
Dear Martin et al.,

Paul Gilbert forwarded me your discussions on this topic.

I have been interested in computing f(A), where A is any square matrix and
f(.) is an analytic function whose Taylor-McLaurin series has a finite
radius of convergence.  There exist special algorithms when f(.) is
square-root, exponential, logarithm (note: the TM series for this doesn't
have a finite-radius of convergence), sine, and cosine.  The easiest way to
handle a general f(.) is via spectral decomposition, but that is limited to
diagonalizable matrices, and is also ill-conditioned.  A paper by Davies and
Higham (SIAM 2003) provide an algorithm, the Schur-Parlett method, for doing
this.  I think this is the state-of-the-art for evaluating f(A).  Higham
also has a recent SIAM book (2008) on computing f(A).  There are matlab
codes provided with that, but it doesn't have the algorithm discussed in his
(2003) paper.  I asked him for the code, but he was unable to share it with
me.  I suspect this may be due to proprietary issues since this is
implemented in the latest Matlab release.  However, we should be able to do
this in R based on the algorithmic description in is (2003) paper.  

Would this be of interest to the R group?

Best,
Ravi.

----------------------------------------------------------------------------
-------

Ravi Varadhan, Ph.D.

Assistant Professor, The Center on Aging and Health

Division of Geriatric Medicine and Gerontology

Johns Hopkins University

Ph: (410) 502-2619

Fax: (410) 614-9625

Email: [hidden email]

Webpage:  http://www.jhsph.edu/agingandhealth/People/Faculty/Varadhan.html

 

----------------------------------------------------------------------------
--------


-----Original Message-----
From: [hidden email] [mailto:[hidden email]]
On Behalf Of Martin Maechler
Sent: Saturday, April 05, 2008 1:52 PM
To: Rory Winston
Cc: [hidden email]
Subject: Re: [Rd] Adding a Matrix Exponentiation Operator

>>>>> "RW" == Rory Winston <[hidden email]>
>>>>>     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> [hidden email] mailing list
    RW> https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel