12:26 p.m.
Hi,
in the brand new 2018d there is the facility to raise positive-definite
matrices to fractional powers, which is great. (See the changelog.)
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?
(I could just check and test, but first I'm lazy and secondly I'm asking
about the intended coverage, which may differ from what is the case now.)
thanks,
sven
2:08 p.m.
On Sat, 22 Dec 2018, Sven Schreiber wrote:
Hi,
in the brand new 2018d there is the facility to raise positive-definite
matrices to fractional powers, which is great. (See the changelog.)
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?
(I could just check and test, but first I'm lazy and secondly I'm asking
about the intended coverage, which may differ from what is the case now.)
Good question. Actually the new code also works for psd matrices
(this was part of a late-breaking change from Jack). But it turns
out that the specialized code for psdroot is slightly more accurate
for the pow = 0.5 case.
<hansl>
# make a psd matrix
n = 10
A = zeros(n, n)
loop i=1..8 -q
x = muniform(n, 1)
A += x*x'
endloop
# expose non-positive eigenvalues
eval eigensym(A)
# apply psdroot()
r = psdroot(A)
eval A - r*r'
# apply new fractional power code
r = A^0.5
eval A - r*r'
</hansl>
I'm seeing mostly zeros reported on the first check, and mostly very
small but non-zero values on the second.
Allin
2:44 p.m.
Am 22.12.2018 um 20:08 schrieb Allin Cottrell:
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?
Good question. Actually the new code also works for psd matrices
(this
was part of a late-breaking change from Jack). But it turns out that the
specialized code for psdroot is slightly more accurate for the pow = 0.5
case.
Interesting, thanks. In that case it would seem that internally an
expression like A^0.5 should be remapped to psdroot(A) or the respective
underlying C code?
Also, I just ran a small speed comparison and it seems that psdroot() is
over 30 times faster than A^0.5 (on a pd matrix, no semi).
thanks,
sven
3:01 p.m.
On Sat, 22 Dec 2018, Sven Schreiber wrote:
Am 22.12.2018 um 20:08 schrieb Allin Cottrell:
> 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?
> Good question. Actually the new code also works for psd matrices (this was
> part of a late-breaking change from Jack). But it turns out that the
> specialized code for psdroot is slightly more accurate for the pow = 0.5
> case.
Interesting, thanks. In that case it would seem that internally an expression
like A^0.5 should be remapped to psdroot(A) or the respective underlying C
code?
Also, I just ran a small speed comparison and it seems that psdroot() is over
30 times faster than A^0.5 (on a pd matrix, no semi).
OK, that's now in git: A^0.5 calls gretl_matrix_psd_root().
Allin
8:49 a.m.
Am 22.12.2018 um 21:01 schrieb Allin Cottrell:
On Sat, 22 Dec 2018, Sven Schreiber wrote:
> Am 22.12.2018 um 20:08 schrieb Allin Cottrell:
>> 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?
>
>> Good question. Actually the new code also works for psd matrices
>> (this was part of a late-breaking change from Jack). But it turns out
>> that the specialized code for psdroot is slightly more accurate for
>> the pow = 0.5 case.
OK, hold on. As clearly documented psdroot() is not actually calculating
what is "typically" called the square root of a psd matrix.
With "typically" I mean the --apparently quite universal-- convention
that M = RR should hold in order to call R the square root of M. R will
also be psd (and thus symmetric).
In contrast, psdroot generalizes the Choleski decomposition and returns
a triangular matrix. So we have M = PP', but not = PP or = P'P.
(M = RR' also holds because, again, R is symmetric anyway. That's why
your comparison PP' vs. RR' was OK.)
> Also, I just ran a small speed comparison and it seems that
psdroot()
> is over 30 times faster than A^0.5 (on a pd matrix, no semi).
OK, that's now in git: A^0.5 calls gretl_matrix_psd_root().
I guess this should be reverted.
Perhaps the naming of psdroot() is actually a little misleading, at
least it misled us. But it has long standing of course.
cheers,
sven
11:18 a.m.
On Sun, 23 Dec 2018, Sven Schreiber wrote:
Am 22.12.2018 um 21:01 schrieb Allin Cottrell:
> On Sat, 22 Dec 2018, Sven Schreiber wrote:
>
>> Am 22.12.2018 um 20:08 schrieb Allin Cottrell:
>>> 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?
>>
>>> Good question. Actually the new code also works for psd matrices (this
>>> was part of a late-breaking change from Jack). But it turns out that the
>>> specialized code for psdroot is slightly more accurate for the pow = 0.5
>>> case.
OK, hold on. As clearly documented psdroot() is not actually calculating what
is "typically" called the square root of a psd matrix.
With "typically" I mean the --apparently quite universal-- convention that M
= RR should hold in order to call R the square root of M. R will also be psd
(and thus symmetric).
In contrast, psdroot generalizes the Choleski decomposition and returns a
triangular matrix. So we have M = PP', but not = PP or = P'P.
(M = RR' also holds because, again, R is symmetric anyway. That's why your
comparison PP' vs. RR' was OK.)
>> Also, I just ran a small speed comparison and it seems that psdroot() is
>> over 30 times faster than A^0.5 (on a pd matrix, no semi).
>
> 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.
Allin
12:16 p.m.
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
-------------------------------------------------------
12:47 p.m.
Am 23.12.2018 um 18:16 schrieb Riccardo (Jack) Lucchetti:
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()
Sounds good. Would there be any problem just leaving the name cholesky()
and explaining that it became actually more general?
cheers,
sven
2:37 p.m.
On Sun, 23 Dec 2018, Sven Schreiber wrote:
Am 23.12.2018 um 18:16 schrieb Riccardo (Jack) Lucchetti:
>
> 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()
Sounds good. Would there be any problem just leaving the name cholesky() and
explaining that it became actually more general?
Hmm, cholesky and psdroot are not exactly nested at present: cholesky
checks for both symmetry and pd-ness, while psdroot just assumes
symmetry and does not check rigorously for psd-ness unless the
@psdcheck argument is given a non-zero value.
Allin
2156
days inactive
2157
days old
8 comments
3 participants
participants (3)
-
Allin Cottrell -
Riccardo (Jack) Lucchetti -
Sven Schreiber