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

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

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.