From mboxrd@z Thu Jan 1 00:00:00 1970 From: Faheem Mitha To: James Theiler Cc: Subject: Re: calculating matrix inverses (was Re: modifying matrix allocationfunctions for use with R) Date: Wed, 19 Dec 2001 13:20:00 -0000 Message-id: References: X-SW-Source: 2001/msg00592.html 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.