Re: [Help-gsl] Best way to Calculating inv(A) . B

[John D Lamb]

On Fri, 2010-04-02 at 14:27 +1100, Srimal Jayawardena wrote: 
> One way I can think of is to use,
> 
>  gsl_linalg_LU_solve (&m.matrix, p, &b.vector, x);
> 
> as explained in
> 
> http://www.gnu.org/software/gsl/manual/html_node/Linear-Algebra-Examples.html
> 
> I would need to find the solution vector x for each col vector b (of
> matrix B) and and concatenate the x vectors to obtain my inv(A).B
> matrix.
> 
> Is there a better / more direct approach that is more numerically
> stable and accurate ?

Two possibilities occur to me. Both use gsl_linalg_LU_refine(). This
uses some form of iterative improvement (perhaps Jacobi — I haven’t
checked) and computes a vector of residuals r = Ax – b, which should be
zero if x solves Ax = b. Something like the following ought to give more
numerically accurate solutions.

double d;
gsl_vector* r;

for( d = 10; d > TOLERANCE; ){
  gsl_linalg_LU_refine( A, LU, p, b, x, r );
  d = gsl_blas_dnrm2( r );
}

This should improve the value of x. It would make sense to check that d
is actually decreasing because (i) iterative improvements are not
guaranteed to converge (though they typically do), (ii) the improvement
will be limited by the numerical accuracy of double-precision
calculations.

The first possibility is to improve each estimate x of a column of
inv(A).B in turn. The second is to estimate inv(A) directly and use
iterative improvement to improve the estimates of the columns of inv(A).
I think I would use the first method.

This is not more direct, but it should typically be more numerically
accurate.

-- 
JDL