Re: [Gretl-devel] slight qform() generalization
On Wed, 1 Oct 2014, Allin Cottrell wrote:
Let me just add a footnote to this discussion. My claim, infra, is
more or less, "Don't try to second-guess gretl in respect of matrix
operations, just write the 'natural' code (e.g. "X'X" for X'X)
and
chances are you'll get optimal performance."
I'm reasonably confident about this in general, but there may well
be some counter-examples. If you find one, please let us know and
we'll try to fix it.
As shown by Jack's timing comparison script (also infra) there's no
need to speculate about whether hansl code variant A is or is not
faster than hansl code variant B. Just wrap the calculation in a big
fat loop and apply $stopwatch to both versions -- the empirical
method.
To exemplify both points, consider the following script:
<hansl>
set echo off
set messages off
scalar dim = 1000
scalar rep = 600
scalar H = 10
matrix times = zeros(H,3)
loop h=1..H --quiet
a = mnormal(dim, h)
b = a'
c = a
matrix x = zeros(dim, dim)
set stopwatch
loop rep --quiet
matrix x = a'a
endloop
times[h,1] = $stopwatch
set stopwatch
loop rep --quiet
matrix x = b*a
endloop
times[h,2] = $stopwatch
set stopwatch
loop rep --quiet
matrix x = c'a
endloop
times[h,3] = $stopwatch
endloop
print times
</hansl>
Running the above will give you all sorts of interesting points to ponder
on.
On Thu, 2 Oct 2014, Sven Schreiber wrote:
Well, the whole thing started because of the help text of qform():
"Using this function instead of ordinary matrix multiplication
guarantees more speed and better accuracy." So this would seem like an
exception to your rule?
(If so, I guess it would come from leveraging the extra information that
the matrix A is symmetric.)
Uh, I guess you're right. I edited the doc file to that effect.
-------------------------------------------------------
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
-------------------------------------------------------