Re: [Gretl-devel] psdroot vs. M^0.5
On Sun, 23 Dec 2018, Allin Cottrell wrote:
On Sun, 23 Dec 2018, Sven Schreiber wrote:
>>>> On Sat, 22 Dec 2018, Sven Schreiber wrote:
>>>
>>>>> I'm wondering what's then the role of the psdroot() function
now. Is it
>>>>> covering semi-definite matrices (the "s" in psd) for which
the new
>>>>> syntax would fail?
>> OK, that's now in git: A^0.5 calls gretl_matrix_psd_root().
>
> I guess this should be reverted.
OK, if it's potentially misleading. Done in git.
In fact, that made me think that the present cholesky() function is just a
second-class psdroot() function: for pd matrices they produce the same
result, but cholesky() fails for singular matrices. However, the
cholesky() function calls a LAPACK function which is much more efficient
for larger matrices, since it's automatically parallelised (psdroot()
appears to be quicker for smaller matrices). I've experimented via the
following script:
<hansl>
set verbose off
n = 50
H = 4000
X = mnormal(n+1,n)
A = X'X
tc = 0
tp = 0
loop H -q
set stopwatch
cholesky(A)
tc += $stopwatch
endloop
loop H -q
set stopwatch
psdroot(A)
tp += $stopwatch
endloop
print tc tp
</hansl>
We may consolidate the two functions and make one an alias of the other,
doing something like this:
(a) check if the matrix is "small"; if so -> psdroot()
(b) try cholesky(); if error -> psdroot()
unless, of course, we find an inexpensive way to check for singularity
before (b). Or perhaps, use the LU decomposition? (BTW: we don't have the
LU decomposition as per LAPACK's dgetrf() as a user-level function,
although we use it internally: maybe we should)
-------------------------------------------------------
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
-------------------------------------------------------