Direct solving with Cholesky factors

[chaozg]

Hi Luiz,

Thank you very much for adding the Cholesky factorization feature – it’s fantastic!

I have a question regarding the use of the computed Cholesky factor for solving linear systems. Specifically, can the Cholesky factor be used to solve equations of the form Lx=b or L^Tx=b.

I have a minimal (but not working code) here and my goal is to use hchol.L as a direct solver:

using HMatrices, LinearAlgebra, StaticArrays, Random

# copied from README.md
const Point3D = SVector{3,Float64}
m = 800
X = Y = [Point3D(sin(θ)cos(ϕ),sin(θ)*sin(ϕ),cos(θ)) for (θ,ϕ) in zip(π*rand(m),2π*rand(m))]
function G(x,y) 
  d = norm(x-y) + 1e-8
  exp(-d)
end
K = KernelMatrix(G,X,Y)

# copied from test/cholesky_test.jl
splitter = CardinalitySplitter(; nmax = 50)
Xclt = ClusterTree(X, splitter)
Yclt = ClusterTree(Y, splitter)
adm = StrongAdmissibilityStd(10)
comp = PartialACA(; atol = 1e-6)
H = assemble_hmatrix(Hermitian(K), Xclt, Yclt; adm, comp, threads=false, distributed = false)

# solve Lx=b
b = randn(m)
hchol = cholesky(H; atol = 1e-6)
x = hchol.L\b

For now, the last line leads to the following error message:

ERROR: method `getindex(::HMatrix,args...)` has been disabled to
avoid performance pitfalls. Unless you made an explicit call to `getindex`,
this error usually means that a linear algebra routine involving an
`HMatrix` has fallen back to a generic implementation.
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] getindex(::Hermitian{Float64, HMatrix{ClusterTree{3, Float64}, Float64}}, i::Int64, j::Int64)
   @ HMatrices ~/.julia/packages/HMatrices/OBct7/src/hmatrix.jl:87
 [3] generic_trimatdiv!(C::Vector{…}, uploc::Char, isunitc::Char, tfun::typeof(identity), A::Hermitian{…}, B::Vector{…})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/triangular.jl:1199
 [4] _ldiv!
   @ ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/triangular.jl:752 [inlined]
 [5] ldiv!
   @ ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/triangular.jl:745 [inlined]
 [6] \(A::LowerTriangular{Float64, Hermitian{Float64, HMatrix{ClusterTree{3, Float64}, Float64}}}, B::Vector{Float64})
   @ LinearAlgebra ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/triangular.jl:1483
 [7] top-level scope
   @ REPL[30]:1
Some type information was truncated. Use `show(err)` to see complete types.

Could you provide any guidance on how to properly use the Cholesky factor in this context?

Thank you again for your help!

Charlie

[maltezfaria]

You can do

ldiv!(hchol.L,b)

to solve the system and write the solution inside b. It seems that hchol.L \ b actually dispatches to the three argument version of ldiv!, not implemented in HMatrices.jl. If you want the backslash to work, you could probably just add ldiv!(b,L,a) by calling the two argument version in the triangular.jl file. If you do add that, feel free to submit a PR with the changes so that the more friendly L \ b is supported.

[chaozg]

Hi @maltezfaria ,

Thanks very much for the hint; I was actually trying to play with ldiv!, but the solution seems wrong. Specifically, I first compute y=L*x1 and then compute x2=L\y. If x1==x2, I guess then it implies it's the correct way for computing L\y.

Do I need to permute the vectors into external/internal order somewhere?

x1 = randn(m)
y = hchol.L*x1

ldiv!(hchol.L, y) # y is actually x2 now since idiv! is in-place

plot(x1, marker=:x)
plot!(y, marker=:o)

[chaozg]

I was also trying to do permutation (by stealing some code from src/cholesky.jl) but unfortunately still haven't made it work. The minimal and full code is here

using HMatrices, LinearAlgebra, StaticArrays, Random
using HMatrices: coltree, rowtree, loc2glob, permute!, invpermute!

using Plots

# copied from README.md
const Point3D = SVector{3,Float64}
m = 800
X = Y = [Point3D(sin(θ)cos(ϕ),sin(θ)*sin(ϕ),cos(θ)) for (θ,ϕ) in zip(π*rand(m),2π*rand(m))]
function G(x,y) 
  d = norm(x-y) + 1e-8
  exp(-d)
end
K = KernelMatrix(G,X,Y)

# copied from test/cholesky_test.jl
splitter = CardinalitySplitter(; nmax = 50)
Xclt = ClusterTree(X, splitter)
Yclt = ClusterTree(Y, splitter)
adm = StrongAdmissibilityStd(10)
comp = PartialACA(; atol = 1e-6)
H = assemble_hmatrix(Hermitian(K), Xclt, Yclt; adm, comp, threads=false, distributed = false)

# Cholesky factorization
hchol = cholesky(H; atol = 1e-6)

# Experiment with ldiv! and (inv)permute!:
# 1) compute y=L*x1 and then 2) compute x2=L\y
# adapted from ldiv! of src/cholesky.jl
p = hchol.factors
ctree = coltree(p)
rtree = rowtree(p)

# 1. Compute y/y_copy=L*x with manual permutation
x1 = randn(m)
# permute x1 to internal order
permute!(x1, loc2glob(ctree))
# multiplication: L*x1
y = hchol.L*x1
# permute y to external order
invpermute!(y, loc2glob(ctree))

# 2. Compute x2= L\y with manual permutation
# permute y to internal order
permute!(y, loc2glob(ctree))
# division: L\y
ldiv!(hchol.L, y)
# permute y (actually x2 since it's in-place) back to external order
invpermute!(y, loc2glob(ctree))

plot(x1, marker=:x)
plot!(y, marker=:o)

[chaozg]

I just figured out the following code gives the same x1 and x2/y:

x = randn(m)
x1 = deepcopy(x)
y = hchol.L*x

# permute y to internal order
permute!(y, loc2glob(ctree))
ldiv!(hchol.L, y)
# permute y (actually x2 since it's in-place) back to external order
invpermute!(y, loc2glob(ctree))

plot(x1, marker=:x)
plot!(y, marker=:o)

Does the code make sense to you? So ldiv!(hchol.L, y) expect y to be in internal order while hchol.L*x expects x to be in external order?

[maltezfaria]

You are right: one has to be careful with the internal vs. external ordering here!

From a quick look at the code it does seem that ldiv! uses the internal ordering, which means this function should really be renamed to something like _ldiv! to prevent unintended external usage 🤔

[chaozg]

Thanks for all the discussion!

I'll have a bit more tests with permute and ldiv!. In case I see anything strange, I'll let you know.

BTW, is ctree retrieved from the Cholesky factor the same as Xclt? I mean

p = hchol.factors
ctree = coltree(p)
rtree = rowtree(p)

and

Xclt = ClusterTree(X, splitter)

[chaozg]

Hi Luiz, Just want to report to you that the result from (inv)permute and ldiv! matches that from gmres: (might be) an evidence that everything is working right.

[maltezfaria]

Thanks for testing and reporting the findings.

BTW, is ctree retrieved from the Cholesky factor the same as Xclt? I mean

I think so. AFAIR the cluster trees do not get mutated in the construction.

Hi Luiz, Just want to report to you that the result from (inv)permute and ldiv! matches that from gmres: (might be) an evidence that everything is working right.

I am still skeptical that this works as expected, so it is better to leave this issue open for now. I think using hchol \ x is fine, since that dispatches to an ldiv method in cholesky.jl that correctly handles the permutations, but once you start to do things like hchol.L \ x you have to be careful with the internal/external numbering.

I think the proper "fix" here is to clearly distinguish internal methods which always use the local numbering induced by the clustering algorithms and external methods which handle the permutations. If we do that, then things like ldiv!(hchol.L, x) will either error because the method does not exist, or accept the global_index keyword and internally handle the required permutations.

[chaozg]

Hi Luiz,

I wonder if it is possible to export/evalute the dense matrix behind the computed Cholesky factor. If that's possible, the verification of the correctness of chol.L\x might be easier as we can then directly compare with something like Matrix(chol.L)\x.

Charlie

[maltezfaria]

I wonder if it is possible to export/evalute the dense matrix behind the computed Cholesky factor. If that's possible, the verification of the correctness of chol.L\x might be easier as we can then directly compare with something like Matrix(chol.L)\x.

This starts to tweak with the internals a bit, but assuming you stored the data on the upper triangular part for the Cholesky factorization (the default AFAIU) you can do:

U = UpperTriangular(Matrix(hchol.factors))
L = adjoint(U)

Here is what it does:

1. Extract the underlying HMatrix behind the Cholesky factorization using hchol.factors
2. Convert that HMatrix into a Matrix. Here you can pass the global_index keyword argument (default to true) to handle the permutations. Note that only the upper triangular part will be meaningful (assuming hchol.uplo = 'U')
3. Wrap the matrix into an upper triangular
4. Transpose the meaningful bit to get the lower triangular

I will be busy for the next couple of weeks (summer conferences), but if you start implementing any of these changes, feel free to open a draft PR so that we can further discuss the implementation/issues.

[maltezfaria]

Hi @chaozg,

I am going over some old tickets in this repo, and I was wondering if you ever managed to do what you wanted with HMatrices. If not, can you remind me what the issue was? Perhaps I can give this another try...

[chaozg]

Hi @maltezfaria , thanks for the follow-up message. So far I'm happy with the computed Cholesky factor though I didn't spend much time on further experiments you suggested yet. I'm still keen on using HMatrices.jl in a project later and I'll definitely ask you other questions or raise an issue some time!

Feel free to close this issue if you need to.

I wonder if it is possible to export/evalute the dense matrix behind the computed Cholesky factor. If that's possible, the verification of the correctness of chol.L\x might be easier as we can then directly compare with something like Matrix(chol.L)\x.

This starts to tweak with the internals a bit, but assuming you stored the data on the upper triangular part for the Cholesky factorization (the default AFAIU) you can do:

U = UpperTriangular(Matrix(hchol.factors))
L = adjoint(U)

Here is what it does:

1. Extract the underlying HMatrix behind the Cholesky factorization using hchol.factors
2. Convert that HMatrix into a Matrix. Here you can pass the global_index keyword argument (default to true) to handle the permutations. Note that only the upper triangular part will be meaningful (assuming hchol.uplo = 'U')
3. Wrap the matrix into an upper triangular
4. Transpose the meaningful bit to get the lower triangular

I will be busy for the next couple of weeks (summer conferences), but if you start implementing any of these changes, feel free to open a draft PR so that we can further discuss the implementation/issues.

[maltezfaria]

Sounds good. At some point I will add better support for Cholesky, but I think it is OK to close this for now.