Hi,
I came across the situation where a relatively large symmetric but
singular matrix was given (quite a few eigenvalues very close to zero)
which I had to "invert". Understandably, gretl's invpd didn't like the
input. So I used ginv instead, but it then turned out the result was not
"sufficiently" symmetric; for example gretl's own qform function
complained about the alleged non-symmetry.
So I'm wondering whether hansl/gretl might benefit from a specialized
ginv function for symmetric input. Below is a proposal in hansl based on
standard algebra. Perhaps we even discussed something similar before?
thanks
sven
<hansl>
function matrix ginvsym(const matrix m, scalar thresh[0::1e-10])
# symmetry is only implicitly checked by eigensym
matrix evecs
matrix evals = eigensym(m, &evecs)
if min(evals) > thresh
print "Warning: generalized inv not required"
matrix out = invpd(m)
# (or alternatively:)
# out = (evecs .* (1 / evals')) * evecs'
else
r = rows(m)
matrix ev = selifc(evecs, evals' .> thresh)
matrix out = ev .* (1 / evals[r - cols(ev) + 1: r]') * ev'
endif
return out
end function
M = mnormal(40,20)
MM = M*M'
matrix GMM = ginvsym(MM)
# check symmetry:
eval max(abs(GMM - GMM'))
# check g-inverse property:
eval max(abs(MM - MM*GMM*MM))
</hansl>