integralequations / hmatrices.jl Goto Github PK

A question please: how do I go about calculating the log-determinant of an HMatrix? Thanks.

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

It would be nice to use DataFlowTasks for the parallel LU factorisation. I am opening this issue to keep track of the changes and PRs needed to make that happen.

• Have a more robust set of benchmarks to judge progress/regression (#38)
• #39
• Add @dspawn for the LU factorisation

Hi there,

It seems the plot(H) gives an error for even a small problem. To reproduce the error, I used the code from README.md and changed m to a very small integer, like 600.

using HMatrices, LinearAlgebra, StaticArrays
const Point3D = SVector{3,Float64}
# sample some points on a sphere
m = 600 #500 is fine but it's dense so nothing is interesting
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
  1/(4π*d)
end
K = KernelMatrix(G,X,Y)

H = assemble_hmatrix(K;atol=1e-6)

using Plots
plot(H)

And the error is as follows:

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(::HMatrix{ClusterTree{3, Float64}, Float64}, i::Int64, j::Int64)
    @ HMatrices ~/.julia/packages/HMatrices/OBct7/src/hmatrix.jl:87
  [3] _getindex
    @ ./abstractarray.jl:1341 [inlined]
  [4] getindex
    @ ./abstractarray.jl:1291 [inlined]
  [5] iterate
    @ ./abstractarray.jl:1217 [inlined]
  [6] iterate
    @ ./abstractarray.jl:1215 [inlined]
  [7] copyto_unaliased!
    @ ./abstractarray.jl:1093 [inlined]
  [8] copyto!(dest::Matrix{Float64}, src::HMatrix{ClusterTree{3, Float64}, Float64})
    @ Base ./abstractarray.jl:1068
  [9] copy_similar
    @ ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/LinearAlgebra.jl:416 [inlined]
 [10] eigencopy_oftype
    @ ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/eigen.jl:538 [inlined]
 [11] svdvals(A::HMatrix{ClusterTree{3, Float64}, Float64})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/svd.jl:242
 [12] rank(A::HMatrix{ClusterTree{3, Float64}, Float64}; atol::Float64, rtol::Float64)
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/generic.jl:991
 [13] rank(A::HMatrix{ClusterTree{3, Float64}, Float64})
    @ LinearAlgebra ~/.julia/juliaup/julia-1.10.4+0.aarch64.apple.darwin14/share/julia/stdlib/v1.10/LinearAlgebra/src/generic.jl:989
 [14] macro expansion
    @ ~/.julia/packages/HMatrices/OBct7/src/hmatrix.jl:485 [inlined]
 [15] macro expansion
    @ ~/.julia/packages/RecipesBase/BRe07/src/RecipesBase.jl:340 [inlined]
 [16] macro expansion
    @ ~/.julia/packages/HMatrices/OBct7/src/hmatrix.jl:482 [inlined]
 [17] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, hmat::HMatrix)
    @ HMatrices ~/.julia/packages/RecipesBase/BRe07/src/RecipesBase.jl:300
 [18] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/user_recipe.jl:38
 [19] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
    @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/RecipesPipeline.jl:72
 [20] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
    @ Plots ~/.julia/packages/Plots/sxUvK/src/plot.jl:223
 [21] #plot#188
    @ ~/.julia/packages/Plots/sxUvK/src/plot.jl:102 [inlined]
 [22] plot(args::Any)
    @ Plots ~/.julia/packages/Plots/sxUvK/src/plot.jl:93
 [23] top-level scope
    @ REPL[21]:1

evaluate_greens_func evaluates to complex number. I get error that it is not supported."MethodError: no method matching assemble_hmat(::Vector{ComplexF64}; threads=false)"
`using HMatrices, LinearAlgebra, StaticArrays
const Point3D = SVector{3,Float64}
X = rand(Point3D,10_000)
Y = rand(Point3D,10_000)
GF,GF1,GF2,GF3 = KernelMatrix(evaluate_greens_func,X,Y)

Xclt = Yclt = ClusterTree(X)
adm = StrongAdmissibilityStd()
comp = PartialACA(;atol=1e-6)
H = assemble_hmat(GF,threads = false)`

It would be nice to have a convenient wrapper so that we easily construct an HMatrix from an IntegralOperator from WavePropBase. Something along the lines of the following syntax:

iop = WavePropBase.IntegralOperator(K,X,Y)
H = HMatrix(iop;kwargs...)

where kwargs... are passed down to the assemble_hmatrix constructor.

It should be possible to make HMatrix work with kernels which return entries that are e.g. SMatrix. This would be useful for compressing Nystrom operators for tensorial PDEs such as elastostatic and Maxwell. Some things which need to be changed include:

• The ACA compressor needs to be modified so that the choice of pivot is based on some norm of the entry instead of its abs.
• The RkMatrix needs to support matrices of SMatrix.

It would be useful to be able to use hcat and vcat on hierarchical matrices, and there are some problems where one naturally has a "block matrix" structure. The concatenated object could just be a light wrapper, in the form of an AbstractMatrix, around the HMatrices being concatenated.

As Makie is becoming more mature, I think it makes sense to provide recipes which make use of Makie instead of Plots.

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 Luiz,

I am playing around with this package for Helmholtz problems. Do you recommend using the assemble_hmat from this package of the assemble_ifgf from your IFGF.jl package?

Cheers,
Mikkel

In some boundary integral equation applications, it makes sense to add a sparse correction matrix to an HMatrix approximation. In those cases, as it turns out, the sparse matrix only has non-zero entries precisely where the HMatrix has dense blocks, so adding the two should be rather simple and cost-effective (i.e. no increase in the rank).

• Add an extension where, if SparseArrays is loaded, the axpy! or axpby! function is defined to add a sparse matrix in-place to an HMatrix
• If the sparse matrix is nonzero in blocks where the HMatrix is low-rank, decide what to do

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Hi, running the following block of code return an error.

using HMatrices, LinearAlgebra, StaticArrays
const Point3D = SVector{3,Float64}
# sample some points on a sphere
m = 10_000
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
  1/(4π*d)
end
K = KernelMatrix(G,X,Y)
H = assemble_hmat(K;atol=1e-6)
using Plots
plot(H)

The returned error is :

UndefVarError: Leaves not defined

In Kernel matrices documentation it's stated that elements can be of any type:

The row and columns elements may be as simple as points in $\mathbb{R}^d$ (as is the case for Nyström methods), but they can also be more complex objects such as triangles or basis functions –- the only thing required is that f(X[i],Y[j]) make sense.

However, I run into problems constructing ClusterTrees if I try to exercise this ability. Some assorted observations:

• as suggested in WavePropBase, I've implemented HMatrices.center(::MyType)::SVector{3,Float64}

• in the documentation, it's stated that ClusterTree comes from WavePropBase, but actually it seems this is no longer true.

• whilst WavePropBase's ClusterTree supports arbitrary element types, the current version in HMatrices only supports SVectors.

• Trying to use WavePropBase.ClusterTree explicitly results in a MethodError for isleaf.

Am I approaching this case incorrectly? Or is this a bug? In any case, I'm looking forward to using HMatrices, as it should be directly applicable to my problem if it can handle this case.
Thanks!

Hi there,

I am currently using this nice package for various matrices in relation to the acoustical boundary element method. An example of such matrix is

$$ B = \begin{bmatrix} G(k,t_1,y_1) & G(k,t_1,y_2) & \dots & G(k,t_1,y_{NQ}) \\ G(k,t_2,y_1) & G(k,t_2,y_2) & \dots & G(k,t_2,y_{NQ}) \\ \vdots & \vdots & \ddots & \vdots\\ G(k,t_n,y_1) & G(k,t_n,y_2) & \dots & G(k,t_n,y_{NQ}) \\ \end{bmatrix}, $$

where G is the Green's function for the scalar Helmholtz equation. Notice that the system is not square, meaning that i need to solve the resulting linear system using e.g. lsqr. Doing so requires products with the conjugate transpose, which fails for the HMatrix format (due to getindex not being defined). Was there another reason than time why the conjugate transpose was left out? Would it not be "easily" by implemented by conjugating each of the blocks?

My workaround so far is to simply define another HMatrix as

$$ B^\text{H} = \begin{bmatrix} G(-k, t_1, y_1) & G(-k, t_2, y_1) & \dots & G(-k, t_n, y_{1}) \\ G(-k, t_1, y_2) & G(-k, t_2, y_2) & \dots & G(-k, t_{n}, y_{2}) \\ \vdots & \vdots & \ddots & \vdots\\ G(-k, t_1, y_{NQ}) & G(-k, t_2, y_{NQ}) & \dots & G(-k, t_n, y_{NQ}) \\ \end{bmatrix}, $$

but this of course doubles both assembly time and storage, so I would rather not have to do that 😄

Cheers,
Mikkel

Turns out WeakKeyDict does not work as I though: only the keys are weak references.

This mean that an HMatrix for which a partition has been created and cached in the CACHED_PARTITION global variable will not be GCed, even when it seems to have gone out of scope.

The simplest way to fix this is probably to store WeakRef(partition) instead of the partition itself.

The command

pip install hmatrices

produces

JuliaError: Dev path `/Users/lfaria/Research/WaveProp/HMatrices` does not exist.