Re: [Gretl-devel] psdroot vs. M^0.5

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.