Re: [Gretl-devel] slight qform() generalization
On Wed, 1 Oct 2014, Sven Schreiber wrote:
Hi,
when I'm feeling really nerdy I try to use qform() instead of directly
writing the matrix product Z'Z. But since qform() is of course actually
for x*A*x', I have to insert an identity matrix as the second argument,
like this:
qform(x', I(cols(x)))
So two questions:
1) Is it actually really faster to use qform() for this kind of inner
product where A is the identity matrix? (I know the speed difference
won't be noticeable anyway, but in theory.)
Absolutely not :) Try this:
<hansl>
set echo off
set messages off
H = 20
dim = 100
times = zeros(H, 2)
loop h=1..H --quiet
x = mnormal(dim,h)
set stopwatch
loop 100
matrix q = x'x
endloop
times[h,1] = $stopwatch
loop 100
q = qform(x',I(dim))
endloop
times[h,2] = $stopwatch
endloop
print times
</hansl>
2) If the answer is yes, could qform() be generalized such that the
second argument is optional and defaults to a conformable identity
matrix? Then we could write qform(x') for x'x.
2b) Actually it now occurs to me: with this mysterious but indoubtedly
great new syntax parser that gretl has, could gretl be clever enough to
spot occurrences of x'x and internally execute that as qform(x')? (Still
provided the answer to 1 is yes.)
I had a look at the source and it occurred to me that we may try to
intercept expressions like x'x and apply a specialised algorithm (which
should be somewhat faster) for those: Allin, correct me if I'm wrong, but
we could check for "l" and "r" for being the same guy in the B_TRMUL
branch of the big "switch" inside eval (line 10501 of geneval.c) and then,
if so, simply call matrix_multiply_self_transpose(). No?
-------------------------------------------------------
Riccardo (Jack) Lucchetti
Dipartimento di Scienze Economiche e Sociali (DiSES)
Università Politecnica delle Marche
(formerly known as Università di Ancona)
r.lucchetti(a)univpm.it
http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------