Bug detected in Arpack

General Usage
1

I can say with certainty that the Arpack.eigs has a very serious bug. I am running the following simple code to find the 4 smallest eigenvalues and eigenvectors of a 16X16 symmetric matrix. The file xxx.txt is uploaded here.

using Serialization;
using Arpack;

K = deserialize("xxx.txt");
display(K);
k=4;
(vals, vecs) = Arpack.eigs(K; nev = k, which=:SM);

This is the output I get. You will notice that the matrix is pretty harmless.

16×16 Array{Float64,2}:
  0.659233   -0.195475   -0.0260898     0.0        -0.0917547  -0.234595      0.0       -0.14918     0.0        -0.02051     0.0         0.0        -0.031589    0.0         0.0         0.0
 -0.195475    0.61303    -0.141618     -0.0249161  -0.140206   -0.0322586     0.0        0.0         0.0        -0.11133     0.0         0.0        -0.014688    0.0        -0.0718231   0.0
 -0.0260898  -0.141618    0.853604      0.0        -0.130495   -0.0147248     0.0        0.0         0.0        -0.14347    -0.0275495   0.0         0.0        -0.0106849   0.0        -0.000746843
  0.0        -0.0249161   0.0           0.789886    0.0         0.0          -0.223093   0.0         0.0         0.0         0.0        -0.0826823   0.0        -0.232067   -0.0259092   0.0
 -0.0917547  -0.140206   -0.130495      0.0         0.666155   -0.0479611     0.0       -0.182288    0.0        -0.102586    0.0         0.0         0.0         0.0         0.0         0.0
 -0.234595   -0.0322586  -0.0147248     0.0        -0.0479611   0.714565      0.0       -0.155445    0.0        -0.025811   -0.0136642   0.0        -0.0280463   0.0        -0.0307591  -0.000370425
  0.0         0.0         0.0          -0.223093    0.0         0.0           0.75254    0.0         0.0         0.0         0.0        -0.091713    0.0        -0.257414    0.0         0.0
 -0.14918     0.0         0.0           0.0        -0.182288   -0.155445      0.0        0.687753   -0.124747   -0.0093745  -0.0215335  -0.0798659  -0.106118    0.0        -0.0608073  -0.0740178
  0.0         0.0         0.0           0.0         0.0         0.0           0.0       -0.124747    0.731166   -0.062372   -0.0831172  -0.107406   -0.159367    0.0        -0.0885275  -0.118728
 -0.02051    -0.11133    -0.14347       0.0        -0.102586   -0.025811      0.0       -0.0093745  -0.062372    0.748047   -0.185609    0.0        -0.0142867   0.0         0.0        -0.0705392
  0.0         0.0        -0.0275495     0.0         0.0        -0.0136642     0.0       -0.0215335  -0.0831172  -0.185609    0.687992   -0.0446934  -0.0328169   0.0         0.0        -0.210717
  0.0         0.0         0.0          -0.0826823   0.0         0.0          -0.091713  -0.0798659  -0.107406    0.0        -0.0446934   0.797492   -0.0933103  -0.0954023  -0.0550552  -0.168394
 -0.031589   -0.014688    0.0           0.0         0.0        -0.0280463     0.0       -0.106118   -0.159367   -0.0142867  -0.0328169  -0.0933103   0.847823    0.0        -0.0769094  -0.0866605
  0.0         0.0        -0.0106849    -0.232067    0.0         0.0          -0.257414   0.0         0.0         0.0         0.0        -0.0954023   0.0         0.449437    0.0         0.0
  0.0        -0.0718231   0.0          -0.0259092   0.0        -0.0307591     0.0       -0.0608073  -0.0885275   0.0         0.0        -0.0550552  -0.0769094   0.0         0.696568   -0.037958
  0.0         0.0        -0.000746843   0.0         0.0        -0.000370425   0.0       -0.0740178  -0.118728   -0.0705392  -0.210717   -0.168394   -0.0866605   0.0        -0.037958    0.660025
 1ERROR: LoadError: SingularException(16)
Stacktrace:
 [1] checknonsingular at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/factorization.jl:19 [inlined]
 [2] checknonsingular at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/factorization.jl:21 [inlined]
 [3] lu!(::Array{Float64,2}, ::Val{true}; check::Bool) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lu.jl:85
 [4] #lu#136 at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lu.jl:273 [inlined]
 [5] lu at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/lu.jl:272 [inlined] (repeats 2 times)
 [6] factorize(::Array{Float64,2}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.5/LinearAlgebra/src/dense.jl:1279
 [7] _eigs(::Array{Float64,2}, ::UniformScaling{Bool}; nev::Int64, ncv::Int64, which::Symbol, tol::Float64, maxiter::Int64, sigma::Nothing, v0::Array{Float64,1}, ritzvec::Bool) at /home/nodi/.julia/packages/Arpack/o35I5/src/Arpack.jl:164
 [8] #eigs#10 at /home/nodi/.julia/packages/Arpack/o35I5/src/Arpack.jl:46 [inlined]
 [9] #eigs#9 at /home/nodi/.julia/packages/Arpack/o35I5/src/Arpack.jl:45 [inlined]
 [10] top-level scope at /home/nodi/Documents/Julia/t.jl:10
 [11] include(::String) at ./client.jl:457
 [12] top-level scope at REPL[13]:1
in expression starting at /home/nodi/Documents/Julia/t.jl:8

I urge the Julia developers to give this problem some serious consideration. A language meant to do high performance computing should be able to find 4 eigenvalues of a 16X16 symmetric matrix. I could do it myself on pen and paper.

2

Looks like people are already aware of this bug: eigs with singular matrix · Issue #87 · JuliaLinearAlgebra/Arpack.jl · GitHub

and if you’re not willing to use pen and paper, you can also use eigen or GitHub - JuliaLinearAlgebra/ArnoldiMethod.jl: Implicitly Restarted Arnoldi Method, natively in Julia

4 Likes
3

Dear @fgerick thank you very much for the links. I will pursue them.

4

I am wondering if you are aware of

using LinearAlgebra
(vals, vecs) = eigen(Symmetric(K))

That’s very clever of you! Note, however, that this a bug in a package that is not part of the language per se. Bugs happen, and that does not reflect on what Julia is or isn’t able to do. Also, while I am as much of a fan as drama as the next person, it is unlikely to lead to the timelier resolution of bugs.

4 Likes
5

On the contrary it is very much a part of the language, and

keep it up

ভাই আরেকটা কথা বলি আপনাকে। ভুল স্বীকার করতে না শিখলে উন্নতির পথ বন্ধ হয়ে যায়।

6

I can assure you. ARPACK is a package.
The list of standard libraries is here, even for ones in seperate repos (they are listed at the bottom).

All software has bugs. Bugs are intenvitable.
Obviously everyone tried to write bug free code, and has tests trying to prove it.
but it is ineviatable that things slip though.
You can look at any nontrivial project in any language.
Even CPUs are regularly found to have bugs (see e.g. a log of the recent vulnerabilities, but even simpler things like calculators have been found to give incorrect results)
This is why there are issue trackers.

On that note: please report bugs on the issue trackers; not Discourse.
Maintainers don’t nescecarily check Discourse.
and when you maintain dozens of packages trying to track what is going on in Discourse is just not reasonable.

9 Likes
7

Thank you for the link.