`cholesky` specialized methods aren't triggering

[cgeoga]

Thanks again for the great package. Based on looking at the source, I would expect the following code to work and give me the hierarchical Cholesky approximation:

using HMatrices, StaticArrays

pts = rand(SVector{2,Float64}, 1000)
km  = KernelMatrix((x,y)->exp(-norm(x-y)), pts, pts)
hk  = assemble_hmatrix(km, atol=1e-10)
hkf = cholesky(hk; atol=1e-8)

But instead I get this error:

cg:hmatrices_chol> jp demo.jl 
ERROR: LoadError: MethodError: no method matching cholesky(::HMatrix{ClusterTree{2, Float64}, Float64}, ::NoPivot; atol::Float64)
This error has been manually thrown, explicitly, so the method may exist but be intentionally marked as unimplemented.

Closest candidates are:
  cholesky(::AbstractMatrix, ::NoPivot; check) got unsupported keyword argument "atol"
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/cholesky.jl:401
  cholesky(::SymTridiagonal, ::NoPivot; check) got unsupported keyword argument "atol"
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/tridiag.jl:909
  cholesky(::Union{Hermitian{var"#s5021", var"#s5020"}, Hermitian{Complex{var"#s5021"}, var"#s5020"}, Symmetric{var"#s5021", var"#s5020"}} where {var"#s5021"<:Real, var"#s5020"<:SymTridiagonal}, ::NoPivot; check) got unsupported keyword argument "atol"
   @ LinearAlgebra ~/.julia/juliaup/julia-1.11.5+0.x64.linux.gnu/share/julia/stdlib/v1.11/LinearAlgebra/src/special.jl:467
  ...

Stacktrace:
 [1] kwerr(::@NamedTuple{atol::Float64}, ::Function, ::HMatrix{ClusterTree{2, Float64}, Float64}, ::NoPivot)
   @ Base ./error.jl:165
 [2] top-level scope
   @ ~/Scratch/hmatrices_chol/demo.jl:8

And playing with @which on the call of just cholesky(hk) shows that it just falls back to the default indexed routines in LinearAlgebra.

Am I doing something wrong with types here? Or is this something that should work? It is honestly kind of confusing to me why your special methods aren't getting triggered, because from looking at the source I would expect that to work.

Thanks in advance!

[maltezfaria]

Hey @cgeoga ,

This should be fixed on #81:

using LinearAlgebra
using HMatrices, StaticArrays

pts = rand(SVector{2,Float64}, 1000)
km  = KernelMatrix((x, y) -> exp(-norm(x - y)), pts, pts) |> Hermitian
hk  = assemble_hmatrix(km; atol = 1e-10)
hkf = cholesky(hk; atol = 1e-8)

Note that you have to manually tell HMatrix that your kernel is Hermitian. This will be leveraged to save half the memory and flops when assembling hk, and more importantly will create a Hermitian HMatrix object that can be passed to cholesky.

Let me know if this works and I can merge the PR.

[cgeoga]

@maltezfaria, thanks so much! This Hermitian specialization on the KernelMatrix objects is great, and I can confirm that the cholesky call and everything works perfectly now. Thanks!

Also: what is your preferred way for users to cite this library in academic papers?

[cgeoga]

@maltezfaria, sorry to poke again. But I've been kicking the tires on your new bugfix branch and I'm hitting another issue where cholesky is failing for adm=StrongAdmissibilityStd(), but not WeakAdmissibilityStd(). Here is a MWE:

using BenchmarkTools, HMatrices, StaticArrays

struct GaussKernel{D}
  Mv::NTuple{D,Float64}
  prodmv::Float64
  perturbation::Float64
end

function GaussKernel(Mv::NTuple{D,Float64}; perturbation=1e-8) where{D}
  pmv = prod(Mv)
  GaussKernel(Mv, pmv, perturbation)
end

function (gk::GaussKernel{D})(x::SVector{D,Float64}, y::SVector{D,Float64}) where{D}
  inner = sum(j->abs2(gk.Mv[j]*(x[j] - y[j])), 1:D)
  out   = exp(-pi*inner)
  x == y ? out + gk.perturbation : out
end

pts = rand(SVector{2,Float64}, 1000)
ker = GaussKernel((22.36, 22.36), perturbation=1e-8)
km  = Hermitian(KernelMatrix(ker, pts, pts))
hk  = assemble_hmatrix(km, atol=1e-10) # if you put adm=WeakAdmissibilityStd() here, it works.
hkf = cholesky(hk; atol=1e-8)

which produces the following error:

ERROR: LoadError: BoundsError: attempt to access 136-element view(::Vector{Float64}, 1:136) with eltype Float64 at index [-1]
Stacktrace:
  [1] throw_boundserror(A::SubArray{Float64, 1, Vector{Float64}, Tuple{UnitRange{Int64}}, true}, I::Tuple{Int64})
    @ Base ./essentials.jl:14
  [2] checkbounds
    @ ./abstractarray.jl:699 [inlined]
  [3] getindex
    @ ./subarray.jl:342 [inlined]
  [4] _aca_partial(K::HMatrices.MulLinearOp{Float64, Tuple{Adjoint{Float64, HMatrix{ClusterTree{2, Float64}, Float64}}, HMatrix{ClusterTree{2, Float64}, Float64}}, Int64}, irange::Base.OneTo{Int64}, jrange::Base.OneTo{Int64}, atol::Float64, rmax::Int64, rtol::Float64, istart::Int64, buffer_::HMatrices.ACABuffer{Float64})
    @ HMatrices ~/.julia/packages/HMatrices/fzbbL/src/compressor.jl:133
[...]

But I have also seen this one before:

ERROR: LoadError: calling `getindex` of a `MulLinearOp` has been disabled
Stacktrace:
  [1] error(s::String)
    @ Base ./error.jl:35
  [2] getindex(::HMatrices.MulLinearOp{Float64, Tuple{HMatrix{…}, HMatrix{…}}, Int64}, ::Int64, ::Int64)
    @ HMatrices ~/.julia/packages/HMatrices/fzbbL/src/multiplication.jl:145
  [3] (::HMatrices.var"#38#41"{HMatrices.MulLinearOp{Float64, Tuple{…}, Int64}, Base.OneTo{Int64}})(j::Int64)
    @ HMatrices ~/.julia/packages/HMatrices/fzbbL/src/compressor.jl:140

Do you have any suggestions or ideas about what's going on here? I recognize that this is a kernel that decays very quickly, so maybe what is happening is that the Cholesky factor is much less amenable to the H-matrix representation. The only symmetric factorization algorithm for H-matrices I know is the Ambikasaran one for HODLR matrices, which probably doesn't have much tech in common.

Thanks again!

[maltezfaria]

Hey @cgeoga,

I won't have time to really look into this for a couple of weeks, but it seems that the issue is coming from assembling the HMatrix itself, and it is unrelated to cholesky. We may be reaching some corner cases of the PartialACA not properly compressing blocks which are essentially zero, but I can't say for sure without further inspection. Here is a possibly better starting point for this bug:

using BenchmarkTools, HMatrices, StaticArrays, LinearAlgebra

struct GaussKernel{D}
    Mv::NTuple{D,Float64}
    prodmv::Float64
    perturbation::Float64
end

function GaussKernel(Mv::NTuple{D,Float64}; perturbation = 1e-8) where {D}
    pmv = prod(Mv)
    return GaussKernel(Mv, pmv, perturbation)
end

function (gk::GaussKernel{D})(x::SVector{D,Float64}, y::SVector{D,Float64}) where {D}
    inner = sum(j -> abs2(gk.Mv[j] * (x[j] - y[j])), 1:D)
    out   = exp(-pi * inner)
    return x == y ? out + gk.perturbation : out
end

pts = rand(SVector{2,Float64}, 1000)
ker = GaussKernel((22.36, 22.36); perturbation = 1e-8)
km = KernelMatrix(ker, pts, pts)
hk = assemble_hmatrix(km; atol = 1e-2) # if you put adm=WeakAdmissibilityStd() here, it works.
any(isnan.(km)) # OK
any(isnan.(Matrix(hk))) # not OK, NaNs are present

[maltezfaria]

Also, as you pointed out, you may need to tweak the admissibility condition for these Gaussian kernels (if anything to exploit their fast decay). But I don't know much about it...

[cgeoga]

@maltezfaria, thanks for the note! Maybe I'll tinker a little on my end and see what progress I can make in figuring out what's going on and document what I find here. And then at a time when you can hook more in we could chat more and see how reasonable whatever my conclusions were? I really appreciate your time here, and if that isn't for a few weeks that is still great for me.

I was also wondering if fully zero blocks would be an issue here. Maybe the fix is as simple as forcing the blocks to always be at least rank one? But I'll see what I can do as I explore here.

Thanks again!

EDIT: I suspect what's going on is that the lu and cholesky routines don't know what to do with [object].data isa Nothing, which happens when a block has zero rank. So I'm going to try and chase down the method that gets hit in that case in those routines.

Here is a test of something that should probably work:

using StaticArrays, HMatrices

pts = rand(SVector{2,Float64}, 1000)
km  = KernelMatrix((x,y)->Float64(x==y), pts, pts)
hk  = assemble_hmatrix(km; rtol=1e-10)
hkf = cholesky(Hermitian(hk)) # fails for lu as well.

Certainly there is no rounding-type error here, but there are just a ton of blocks that are rank zero. If I can get to the bottom of this before you have time to look, I'll try to put together a PR and add a test like this.

EDIT 2: Following through the callstack a bit, the error seems to be coming up here. And it occurs when L.R and (or?) L.P are nothing. Because the types in RkMatrix are a union and not parameterized, catching this won't be possible in dispatch I don't think. But will experiment with an extra if statement in there to try and short-circuit in the nothing case.

EDIT 3: PR with fix created.

[maltezfaria]

@cgeoga sorry for the super long delay! I will take a look at the PR right now and get back to you.