On Sat, 5 Nov 2016, oleg_komashko@ukr.net wrote:

> Dear all,
> Let X be a matrix and A a symmetric matrix
> (dimentions are compatible)
> What is the best way to compute
>
> x'A^-1x

Well, it looks like a case of inversion + qform() versus
"right division". So let the horses race:

<hansl>
/* x' A^{-1} x : best method? */

set seed 775632

scalar m = 50
scalar n = 10
scalar p = 3

matrix x = mnormal(n, p)
matrix Z = mnormal(m, n)
matrix A = Z'Z

# check equivalence
eval qform(x', inv(A))
eval x'(A\x)

set stopwatch

loop 5000 -q
   matrix R = qform(x', inv(A))
endloop
printf "method 1: %gs\n", $stopwatch

loop 5000 -q
   matrix R = x'(A\x)
endloop
printf "method 2: %gs\n", $stopwatch
</hansl>

<output>
method 1: 0.0922s
method 2: 0.0410s
</output>

And the spotted pony wins.

Allin
_______________________________________________
Gretl-users mailing list
Gretl-users@lists.wfu.edu
http://lists.wfu.edu/mailman/listinfo/gretl-users