在 Julia 中计算复厄尔米特稀疏矩阵的特征值
Compute eigenvalues of complex-hermitian sparsematrix in Julia
我正在使用大约 100000x100000 hermitian complex 稀疏矩阵,大约有 5填充的条目百分比,并希望计算 eigenvalues/eigenvectors.
到目前为止我一直在使用 Arpack.jl eigs(A)
。
但是,一旦我将尺寸调到高于 5000,这就不能正常工作了。
对于基准测试,我一直在使用以下代码生成一些 TestMatrices:
using Arpack
using SparseArrays
using ProgressMeter
pop = 0.05
n = 3000 # for example
A = spzeros(Complex{Float64}, n, n)
@showprogress for _ in 1:round(Int,pop * (n^2))
A[rand(1:n), rand(1:n)] = rand(Complex{Float64})
end
# make A hermite
A = A + conj(A)
t = @elapsed eigs(A,maxiter=1500) # ends up being ~ 13 seconds
对于 n ~ 3000,eigs()
调用在我的机器上已经花费了 13 秒,对于更大的 n,它不会在任何 'reasonable' 时间内完成或直接退出。
有专门的 package/method 吗?
感谢任何帮助
https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl 好像是你想要的:
julia> let pop=0.05, n=3000
A = sprand(Complex{Float64},n,n, 0.05)
A = A + conj(A)
@time eigs(A; maxiter=1500)
@time decomp, history = partialschur(A, nev=10, tol=1e-6, which=LM());
end;
10.521786 seconds (50.73 k allocations: 3.458 MiB, 0.04% gc time)
2.244129 seconds (19 allocations: 1.892 MiB)
完整性检查:
julia> a,(b,c) = let pop=0.05, n=300
A = sprand(Complex{Float64},n,n, 0.05)
A = A + conj(A)
eigs(A; maxiter=2500), partialschur(A, nev=6, tol=1e-6, which=LM());
end;
julia> a[1]
6-element Vector{ComplexF64}:
14.5707071003175 + 8.218901803015509e-16im
4.493079744504954 - 0.8390429567118733im
4.493079744504933 + 0.8390429567118641im
-0.3415176925293196 + 4.254184281244591im
-0.3415176925293088 - 4.25418428124452im
0.49406553681541177 - 4.229680489599233im
julia> b
PartialSchur decomposition (ComplexF64) of dimension 6
eigenvalues:
6-element Vector{ComplexF64}:
14.570707100307926 + 7.10698633463049e-12im
4.493079906269516 + 0.8390429076809746im
4.493079701528448 - 0.8390430155670777im
-0.3415174262177961 + 4.254183175902487im
-0.34151626930774975 - 4.25418321627979im
0.49406543866702 + 4.229680079205066im
我正在使用大约 100000x100000 hermitian complex 稀疏矩阵,大约有 5填充的条目百分比,并希望计算 eigenvalues/eigenvectors.
到目前为止我一直在使用 Arpack.jl eigs(A)
。
但是,一旦我将尺寸调到高于 5000,这就不能正常工作了。
对于基准测试,我一直在使用以下代码生成一些 TestMatrices:
using Arpack
using SparseArrays
using ProgressMeter
pop = 0.05
n = 3000 # for example
A = spzeros(Complex{Float64}, n, n)
@showprogress for _ in 1:round(Int,pop * (n^2))
A[rand(1:n), rand(1:n)] = rand(Complex{Float64})
end
# make A hermite
A = A + conj(A)
t = @elapsed eigs(A,maxiter=1500) # ends up being ~ 13 seconds
对于 n ~ 3000,eigs()
调用在我的机器上已经花费了 13 秒,对于更大的 n,它不会在任何 'reasonable' 时间内完成或直接退出。
有专门的 package/method 吗?
感谢任何帮助
https://github.com/JuliaLinearAlgebra/ArnoldiMethod.jl 好像是你想要的:
julia> let pop=0.05, n=3000
A = sprand(Complex{Float64},n,n, 0.05)
A = A + conj(A)
@time eigs(A; maxiter=1500)
@time decomp, history = partialschur(A, nev=10, tol=1e-6, which=LM());
end;
10.521786 seconds (50.73 k allocations: 3.458 MiB, 0.04% gc time)
2.244129 seconds (19 allocations: 1.892 MiB)
完整性检查:
julia> a,(b,c) = let pop=0.05, n=300
A = sprand(Complex{Float64},n,n, 0.05)
A = A + conj(A)
eigs(A; maxiter=2500), partialschur(A, nev=6, tol=1e-6, which=LM());
end;
julia> a[1]
6-element Vector{ComplexF64}:
14.5707071003175 + 8.218901803015509e-16im
4.493079744504954 - 0.8390429567118733im
4.493079744504933 + 0.8390429567118641im
-0.3415176925293196 + 4.254184281244591im
-0.3415176925293088 - 4.25418428124452im
0.49406553681541177 - 4.229680489599233im
julia> b
PartialSchur decomposition (ComplexF64) of dimension 6
eigenvalues:
6-element Vector{ComplexF64}:
14.570707100307926 + 7.10698633463049e-12im
4.493079906269516 + 0.8390429076809746im
4.493079701528448 - 0.8390430155670777im
-0.3415174262177961 + 4.254183175902487im
-0.34151626930774975 - 4.25418321627979im
0.49406543866702 + 4.229680079205066im