More Stable invquad(A::PDMat, z::AbstractVector)

[bayesthm]

For A a PDMat and z a Vector, would it not be more numerically stable to compute invquad(A, z) as

zhalf  = A.chol.L \ z
q = dot(zhalf, zhalf)

rather than the current:

PDMats.jl/src/pdmat.jl

Line 72 in 3c74b56

invquad(a::PDMat, x::AbstractVector) = dot(x, a \ x)

In some experiments my suggested way of computing invquad(A,z) more reliably gives positive values (as it always should!) when A and z contain very large values and A is close to singular. This is important when computing logpdf() of a multivariate normal distribution at fairly extreme values using the Distributions.jl package, which appears to use the PDMats invquad() function.

If I understand the source code correctly this is already how the computation is handled in the case of z being a matrix:

PDMats.jl/src/pdmat.jl

Line 110 in 3c74b56

function Xt_invA_X(a::PDMat, x::StridedMatrix)

An analogous suggestion could also pertain to the computing of quad(A, z), though the numerical stability properties there are probably less of a problem. I'm running Julia 1.1.0 and PDMats 0.9.7.

[andreasnoack]

It would be more stable and faster as well. Would you be able to prepare a PR?

[bayesthm]

Yes, I will do so.

[devmotion]

Fixed by #95.