Re: [Gretl-users] matrix operation efficiency question
On Sat, 5 Nov 2016, oleg_komashko(a)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