Re: [Gretl-devel] psdroot vs. M^0.5
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