[R] quadratic form

[Prof Brian Ripley]

On Thu, 3 Nov 2005, Alvarez Pedro wrote:

> On page 22 of the R-introduction guide it's written:
>
> the quadratic form x^{'} A^{-1} x which is used in
> multivariate computations, should be computed by
> something like x%*%solve(A,x), rather than computing
> the inverse of A.
>
> Why isn't it good to compute t(x) %*% solve(A) %*% x?

The answer is only two lines above;

   Numerically, it is both inefficient and potentially unstable to compute
   @code{x <- solve(A) %*% b} instead of @code{solve(A,b)}.

See also the footnote.

-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595