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
> 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.