modifying matrix allocation functions for use with R

modifying matrix allocation functions for use with R

[Faheem Mitha]

Dear People,

Here is a relatively esoteric question. Some of you may be familiar with
R (sometimes called GNU S) which does statistical stuff. It is possible to
use C code compiled into a shared library as a backend for R. Now, I was
wanting to use some GSL routines in this code.

In particular, I need to use allocation functions like gsl_matrix_alloc.
Unfortunately, I'm not sure whether standard calloc/malloc that are used
work with R. I'm going to check this, but I think the answer is most
likely no. R uses a routine called R_alloc which takes storage from the R
heap and manages the cleanup itself (so you don't need to use free).

One option I was considering was rewriting gsl_matrix_alloc (leaving out
the error checking) in terms of R_alloc. However, I can't find the code
for FUNCTION(gsl_block, alloc) (I assume this means gsl_block_alloc)
used in the code for gsl_matrix_alloc. gsl_matrix_alloc itself seems to be
in matrix/init_source.c. Can someone point me to the code for
gsl_block_alloc, or suggest any other options?

I'm not subscribed to the list, so I would be obliged if you would cc me.

                            Sincerely, Faheem Mitha.

Re: modifying matrix allocation functions for use with R

[Timothy H. Keitt]

This code shows how to wrap a foreign pointer and have it garbage 
collected when the R object goes out of scope.

Tim

SEXP
PgSQLconnect(SEXP connInfo) {

  PGconn *conn;
  SEXP connPtr;

  char *connInfoStr = CHAR(STRING_ELT(connInfo, 0));

  PROTECT(connInfo = coerceVector(connInfo, STRSXP));
  conn = PQconnectdb(connInfoStr);
  UNPROTECT(1);

  PROTECT(connPtr = R_MakeExternalPtr((void *) conn,
                      mkChar("PgSQL connection"),
                      R_NilValue));

  R_RegisterCFinalizer(connPtr, (R_CFinalizer_t)PgSQLcloseConnection);

  UNPROTECT(1);

  return(connPtr);

}



Faheem Mitha wrote:

>Dear People,
>
>Here is a relatively esoteric question. Some of you may be familiar with
>R (sometimes called GNU S) which does statistical stuff. It is possible to
>use C code compiled into a shared library as a backend for R. Now, I was
>wanting to use some GSL routines in this code.
>
>In particular, I need to use allocation functions like gsl_matrix_alloc.
>Unfortunately, I'm not sure whether standard calloc/malloc that are used
>work with R. I'm going to check this, but I think the answer is most
>likely no. R uses a routine called R_alloc which takes storage from the R
>heap and manages the cleanup itself (so you don't need to use free).
>
>One option I was considering was rewriting gsl_matrix_alloc (leaving out
>the error checking) in terms of R_alloc. However, I can't find the code
>for FUNCTION(gsl_block, alloc) (I assume this means gsl_block_alloc)
>used in the code for gsl_matrix_alloc. gsl_matrix_alloc itself seems to be
>in matrix/init_source.c. Can someone point me to the code for
>gsl_block_alloc, or suggest any other options?
>
>I'm not subscribed to the list, so I would be obliged if you would cc me.
>
>                            Sincerely, Faheem Mitha.
>

-- 
Timothy H. Keitt
Department of Ecology and Evolution
State University of New York at Stony Brook
Stony Brook, New York 11794 USA
Phone: 631-632-1101, FAX: 631-632-7626
http://life.bio.sunysb.edu/ee/keitt/

Re: modifying matrix allocation functions for use with R

[Brian Gough]

Faheem Mitha writes:
 > One option I was considering was rewriting gsl_matrix_alloc (leaving out
 > the error checking) in terms of R_alloc. However, I can't find the code
 > for FUNCTION(gsl_block, alloc) (I assume this means gsl_block_alloc)
 > used in the code for gsl_matrix_alloc. gsl_matrix_alloc itself seems to be
 > in matrix/init_source.c. Can someone point me to the code for
 > gsl_block_alloc, or suggest any other options?

Based on my reading of the R documentation, I would suggest using the
standard gsl_matrix_alloc in this case, and handling any errors that
occur with the R libray function error().

e.g.,
         m = gsl_matrix_alloc(...)
         if (m == 0) { error ("could not allocate matrix"); } ;
         ....
         gsl_matrix_free(m);

This would only add a few lines to your code, and will work with the
standard GSL library so your extension will be more portable.

If you want to return an R object containing a gsl_matrix which can be
garbage collected then you could use a C++ wrapper, as the C++
interface in R allows the use of separate constructors and
destructors.

regards
Brian Gough

Re: modifying matrix allocation functions for use with R

[Faheem Mitha]

On Sat, 6 Oct 2001, Brian Gough wrote:

> Faheem Mitha writes:
>  > One option I was considering was rewriting gsl_matrix_alloc (leaving out
>  > the error checking) in terms of R_alloc. However, I can't find the code
>  > for FUNCTION(gsl_block, alloc) (I assume this means gsl_block_alloc)
>  > used in the code for gsl_matrix_alloc. gsl_matrix_alloc itself seems to be
>  > in matrix/init_source.c. Can someone point me to the code for
>  > gsl_block_alloc, or suggest any other options?
>
> Based on my reading of the R documentation, I would suggest using the
> standard gsl_matrix_alloc in this case, and handling any errors that
> occur with the R libray function error().
>
> e.g.,
>          m = gsl_matrix_alloc(...)
>          if (m == 0) { error ("could not allocate matrix"); } ;
>          ....
>          gsl_matrix_free(m);
>
> This would only add a few lines to your code, and will work with the
> standard GSL library so your extension will be more portable.

Thanks, this is a good sugestion. After considering all the alternatives,
I decided to go with something like this.

I think I spotted a typo. In the section LU Decomposition, it is written

 - Function: int gsl_linalg_LU_decomp (gsl_matrix * A, gsl_permutation
          * P, int *SIGNUM)
     This function factorizes the square matrix A into the LU
     decomposition PA = LU.  On output the diagonal and upper
     triangular part of the input matrix A contain the matrix R. The
     lower triangular part of the input matrix (excluding the diagonal)
     contains L.  The diagonal elements of L are unity, and are not
     stored.

Surely, this should be "...the input matrix A contain the matrix U"?

Incidentally, it is a little awkward to compute the inverse of a matrix,
since this involves calling both int gsl_linalg_LU_decomp and int
gsl_linalg_LU_invert. Would it not be useful to have a "wrapper" function
which uses both together?

I have one (slighly idle) question. I have noticed that macros are used
quite liberally in the gsl source code, for example you have a macro
called TYPE which seems to do nothing at all, there is another macro
called PUNCTION(foo, bar) which concatenates foo and bar etc. I'm just
wondering what purpose some of these serve, like the TYPE macro for
example. Do they just make the code easier to write?

                                    Sincerely, Faheem Mitha.

Re: modifying matrix allocation functions for use with R

[Brian Gough]

Faheem Mitha writes:

 > Surely, this should be "...the input matrix A contain the matrix U"?

Thanks. Now fixed.

 > Incidentally, it is a little awkward to compute the inverse of a matrix,
 > since this involves calling both int gsl_linalg_LU_decomp and int
 > gsl_linalg_LU_invert. Would it not be useful to have a "wrapper" function
 > which uses both together?

Maybe, but in this case computing the inverse matrix is not something
I want to encourage, as it's usually better to solve Ax=b instead.

 > I have one (slighly idle) question. I have noticed that macros are used
 > quite liberally in the gsl source code, for example you have a macro
 > called TYPE which seems to do nothing at all, there is another macro
 > called PUNCTION(foo, bar) which concatenates foo and bar etc. I'm just
 > wondering what purpose some of these serve, like the TYPE macro for
 > example. Do they just make the code easier to write?

These macros are used to make the functions for different types,
gsl_vector_float, gsl_vector_int, etc.  See the file templates_on.h
for details.

Brian

calculating matrix inverses (was Re: modifying matrix allocationfunctions for use with R)

[Faheem Mitha]

On Wed, 10 Oct 2001, Brian Gough wrote:

> Faheem Mitha writes:

>  > Incidentally, it is a little awkward to compute the inverse of a matrix,
>  > since this involves calling both int gsl_linalg_LU_decomp and int
>  > gsl_linalg_LU_invert. Would it not be useful to have a "wrapper" function
>  > which uses both together?
>
> Maybe, but in this case computing the inverse matrix is not something
> I want to encourage, as it's usually better to solve Ax=b instead.

This is an odd thing to say. I want to calculate the inverse of
(HH^t)^{-1} H, where H is a given matrix, as part of simulating from a
multivariate normal distribution (the given expression is the mean).  How
does this correspond to solving a linear system? I assume you mean here
that A is a matrix and x and b are vectors? I seem to have no option in
this case but to use gsl_linalg_LU_decomp aned gsl_linalg_LU_invert. If
there are any other options I would appreciate being made aware of them.

                                   Sincerely, Faheem Mitha.

Re: calculating matrix inverses (was Re: modifying matrix allocationfunctions for use with R)

[James Theiler]

On Wed, 10 Oct 2001, Faheem Mitha wrote:

]
]
] On Wed, 10 Oct 2001, Brian Gough wrote:
]
] > Faheem Mitha writes:
]
] >  > Incidentally, it is a little awkward to compute the inverse of a matrix,
] >  > since this involves calling both int gsl_linalg_LU_decomp and int
] >  > gsl_linalg_LU_invert. Would it not be useful to have a "wrapper" function
] >  > which uses both together?
] >
] > Maybe, but in this case computing the inverse matrix is not something
] > I want to encourage, as it's usually better to solve Ax=b instead.
]
] This is an odd thing to say. I want to calculate the inverse of
] (HH^t)^{-1} H, where H is a given matrix, as part of simulating from a
] multivariate normal distribution (the given expression is the mean).  How
] does this correspond to solving a linear system? I assume you mean here
] that A is a matrix and x and b are vectors? I seem to have no option in
] this case but to use gsl_linalg_LU_decomp aned gsl_linalg_LU_invert. If
] there are any other options I would appreciate being made aware of them.
]
]                                    Sincerely, Faheem Mitha.
]

This is getting away from software and API design and into numerical
methods, but for this kind of thing, singular value decomposition is
probably a better tool.  It's a little more expensive than LU, but you
have a lot more control over the stability of the estimated inverse.

jt

---------------------------------------------
James Theiler                     jt@lanl.gov
MS-D436, NIS-2, LANL        tel: 505/665-5682
Los Alamos, NM 87545        fax: 505/665-4414
----- Space and Remote Sensing Sciences -----

Re: calculating matrix inverses (was Re: modifying matrix allocationfunctions for use with R)

[Faheem Mitha]

On Wed, 10 Oct 2001, James Theiler wrote:

> On Wed, 10 Oct 2001, Faheem Mitha wrote:
>
> ]
> ]
> ] On Wed, 10 Oct 2001, Brian Gough wrote:
> ]
> ] > Faheem Mitha writes:
> ]
> ] >  > Incidentally, it is a little awkward to compute the inverse of a matrix,
> ] >  > since this involves calling both int gsl_linalg_LU_decomp and int
> ] >  > gsl_linalg_LU_invert. Would it not be useful to have a "wrapper" function
> ] >  > which uses both together?
> ] >
> ] > Maybe, but in this case computing the inverse matrix is not something
> ] > I want to encourage, as it's usually better to solve Ax=b instead.
> ]
> ] This is an odd thing to say. I want to calculate the inverse of
> ] (HH^t)^{-1} H, where H is a given matrix, as part of simulating from a
> ] multivariate normal distribution (the given expression is the mean).  How
> ] does this correspond to solving a linear system? I assume you mean here
> ] that A is a matrix and x and b are vectors? I seem to have no option in
> ] this case but to use gsl_linalg_LU_decomp aned gsl_linalg_LU_invert. If
> ] there are any other options I would appreciate being made aware of them.
> ]
> ]                                    Sincerely, Faheem Mitha.
> ]
>
> This is getting away from software and API design and into numerical
> methods, but for this kind of thing, singular value decomposition is
> probably a better tool.  It's a little more expensive than LU, but you
> have a lot more control over the stability of the estimated inverse.

Yes, but gsl does not currently seem to have inverses implemented using
singular value deccmposition, and I don't have the time, or the
appropriate algorithms to hand, to implement it myself. I am looking for
cut and dried, ready to use solutions. Thanks for the suggestion, anyway.

                   Sincerely, Faheem Mitha.

Re: calculating matrix inverses (was Re: modifying matrix allocationfunctions for use with R)

[Edwin Robert Tisdale]

Faheem Mitha wrote:

> I want to calculate the inverse of (HH^t)^{-1} H
> where H is a given matrix
> as part of simulating from a multivariate normal distribution
> (the given expression is the mean).

For a nonsingular square matrix

	A = (HH^T)

you need to solve

	AX = I

where I is the identity matrix and

	X = A^{-1}

is the inverse of A.

The problem with matrix inversion is that
there is no one right way to do it.
Perhaps, in this case, a Cholesky decomposition
might be the preferred solution method.
Computing the inverse in two steps --
decomposition followed by back substitution --
has the added advantage of allowing you
to evaluate the condition of the matrix
after the decomposition and
before you attempt to divide by zero
during the back substitution.
to divide by zero