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 -----