IntegralEquations
Cholesky factorization
Hi there,
I just found HMatrices.jl and everything looks great. Meanwhile, I'm curious if there is any plan to add methods to perform Cholesky factorization of the approximated matrix.
Regards, Charlie
Hi Charlie,
This should be relatively easy to do by mirroring the first 80 lines or so of the lu.jl file, but I have not done it since SPD matrices don't show up that frequently in my research (and I can always do an LU in those cases at 2x performance hit). What is your use case? Is LU not enough?
Regards, Luiz
p.s. it would probably also be good to indicate whether an HMatrix is symmetric by e.g. passing a symmetric keyword argument to the constructor and wrapping the returned HMatrix object in a (lazy) Symmetric type. That way you can also save a 2x time and memory at assembly.
Hi Luiz,
Thanks for the swift reply!
Yes, Cholesky factorization will be particularly useful in my case/plan for generating Gaussian random signals. Let's say we have a covariance/kernel matrix K = LL^T, then Lx, with x being an iid draw from a standard normal distribution, will be one sample of the desired distribution.
It will be great if we have Cholesky in HMatrices.jl : )
Regards, Charlie
I see, thanks for the explanation. If I find some free time in the next couple of days I will try to add it.
Depending on how much Julia you know, you can also take a crack at it by replicating the lu.jl file (I would be happy to review/contribute to a PR if you put a draft together).
Hi Luiz,
I'll try to give it a go and let you know if I have any progress or give up at some point : )
Charlie
Did you get a chance to look at it? I can probably give this a shot in the next couple of days. Let me know.
Hi Luiz,
That's GREAT! please go ahead and I'm looking forward to your implementation. I haven't started my try yet...sorry
Charlie
I haven't forgotten about this, just did not find the time to get it yet!
In the meantime, a quick question: at twice the computational cost, you can use the lu factorization right? If the matrix is SPD, lu will essentially compute a chol doing a bunch of redundant operations. You just have to be careful with the normalization of the diagonal entries in L, essentially taking its square root to get your Cholesky correct at the main diagonal.
No worries at all. I was also busy with other work recently and didn't have time to explore HMatrices.jl. Please take your time. And please don't feel pressured: it's definitely not an urgent task or something : )
No pressure felt, don't worry 😉
I have a half-working version on a local branch and was wondering if you could provide a MWE using e.g. a dense matrix representation so I can test against it.
Ok, I tested this, and it seems things are working on the dev branch (see the cholesky_test.jl file). There may still be a few functionalities missing for what you need, but the core should be there. Let me know.