integralequations / hmatrices.jl Goto Github PK
View Code? Open in Web Editor NEWA Julia library for hierarchical matrices
License: MIT License
A Julia library for hierarchical matrices
License: MIT License
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.
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:
norm
of the entry instead of its abs
.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).
SparseArrays
is loaded, the axpy!
or axpby!
function is defined to add a sparse matrix in-place to an HMatrix
HMatrix
is low-rank, decide what to doThis 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 thatf(X[i],Y[j])
make sense.
However, I run into problems constructing ClusterTree
s 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
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
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 GC
ed, 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.
A declarative, efficient, and flexible JavaScript library for building user interfaces.
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
An Open Source Machine Learning Framework for Everyone
The Web framework for perfectionists with deadlines.
A PHP framework for web artisans
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
Some thing interesting about web. New door for the world.
A server is a program made to process requests and deliver data to clients.
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Some thing interesting about visualization, use data art
Some thing interesting about game, make everyone happy.
We are working to build community through open source technology. NB: members must have two-factor auth.
Open source projects and samples from Microsoft.
Google ❤️ Open Source for everyone.
Alibaba Open Source for everyone
Data-Driven Documents codes.
China tencent open source team.