Ralph, ARPACK seems to be okay with the singularity (of A) in the
non-generalized case. In the problem here B=I so it's only a generalized
problem by algorithm, not by math and
julia> eigs(A)[1]
6-element Array{Float64,1}:
-100.0
100.0
4.44392e-15
-3.0281e-18
3.5842e-33
1.69212e-33
so here ARPACK is doing something different. Maybe expanding with vectors
from the null space which I guess is fine as long the αs and βs are zeroed
and you reorthogonalize.
On Saturday, August 6, 2016 at 5:02:42 PM UTC-4, Ralph Smith wrote:
>
> The 0.5 behavior is the same as Fortran, and is technically correct.
> ARPACK cannot create more (orthogonal) Krylov vectors than the rank of the
> matrix A,
> so your example needs the additional keyword ncv=2. This doesn't seem to
> be adequately documented in arpack-ng either. The 0.4 behavior looks
> unreliable.
>
> On Saturday, August 6, 2016 at 2:03:21 AM UTC-4, Madeleine Udell wrote:
>>
>> Actually: I'm not sure we should chalk the error up to ARPACK. Julia 0.4
>> is able to produce a (correct, I think) answer to
>>
>> eigs(A, C)
>>
>> but 0.5 gives an ARPACK error. I don't suppose ARPACK changed between
>> Julia versions...!?
>>
>> On Sat, Aug 6, 2016 at 1:54 AM, Madeleine Udell <[email protected]>
>> wrote:
>>
>>> Andreas, thanks for the investigation. I'll use 0.5 for now, and hope
>>> the real problems I encounter are within the capabilities of ARPACK.
>>>
>>> It's embarrassing to be bested by Matlab...
>>>
>>> On Fri, Aug 5, 2016 at 9:23 PM, Andreas Noack <[email protected]>
>>> wrote:
>>>
>>>> I've looked a bit into this. I believe there is a bug in the Julia
>>>> wrappers on 0.4. The good news is that the bug appears to be fixed on 0.5.
>>>> The bad news is the ARPACK seems to have a hard time with the problem. I
>>>> get
>>>>
>>>> julia> eigs(A,C,nev = 1, which = :LR)[1]
>>>> ERROR: Base.LinAlg.ARPACKException("unspecified ARPACK error: -9999")
>>>> in aupd_wrapper(::Type{T},
>>>> ::Base.LinAlg.#matvecA!#67{Array{Float64,2}}, ::Base.LinAlg.##64#71{Array
>>>> {Float64,2}}, ::Base.LinAlg.##65#72, ::Int64, ::Bool, ::Bool, ::String,
>>>> ::Int64, ::Int64, ::String, :
>>>> :Float64, ::Int64, ::Int64, ::Array{Float64,1}) at ./linalg/arpack.jl:53
>>>> in #_eigs#60(::Int64, ::Int64, ::Symbol, ::Float64, ::Int64, ::Void,
>>>> ::Array{Float64,1}, ::Bool, ::B
>>>> ase.LinAlg.#_eigs, ::Array{Float64,2}, ::Array{Float64,2}) at
>>>> ./linalg/arnoldi.jl:271
>>>> in (::Base.LinAlg.#kw##_eigs)(::Array{Any,1}, ::Base.LinAlg.#_eigs,
>>>> ::Array{Float64,2}, ::Array{Floa
>>>> t64,2}) at ./<missing>:0
>>>> in #eigs#54(::Array{Any,1}, ::Function, ::Array{Float64,2},
>>>> ::Array{Float64,2}) at ./linalg/arnoldi.
>>>> jl:80
>>>> in (::Base.LinAlg.#kw##eigs)(::Array{Any,1}, ::Base.LinAlg.#eigs,
>>>> ::Array{Float64,2}, ::Array{Float6
>>>> 4,2}) at ./<missing>:0
>>>>
>>>> and since SciPy ends up with the same conclusion I conclude that the
>>>> issue is ARPACK. Matlab is doing something else because they are able to
>>>> handle this problem.
>>>>
>>>> Given that 0.5 is almost release, I'll not spend more time on the issue
>>>> on 0.4. Thought, if anybody is able to figure out what is going on, please
>>>> let us know.
>>>>
>>>> On Friday, August 5, 2016 at 8:47:26 AM UTC-4, Madeleine Udell wrote:
>>>>>
>>>>> Setting `which=:LR, nev=1` does not return the generalized eigenvalue
>>>>> with the largest real parts, and does not give a warning or error:
>>>>>
>>>>> n = 10
>>>>> C = eye(n)
>>>>> A = zeros(n,n)
>>>>> A[1] = 100
>>>>> A[end] = -100
>>>>> @assert eigs(A, C, nev=1, which=:LR)[1][1] == maximum(eigs(A, C)[1])
>>>>>
>>>>> Am I expected to set nev greater than the number of eigenvalues I
>>>>> truly desire, based on my intuition as a numerical analyst? Or has eigs
>>>>> broken its implicit guarantee?
>>>>>
>>>>>
>>>
>>>
>>> --
>>> Madeleine Udell
>>> Assistant Professor, Operations Research and Information Engineering
>>> Cornell University
>>> https://people.orie.cornell.edu/mru8/
>>> (415) 729-4115
>>>
>>
>>
>>
>> --
>> Madeleine Udell
>> Assistant Professor, Operations Research and Information Engineering
>> Cornell University
>> https://people.orie.cornell.edu/mru8/
>> (415) 729-4115
>>
>
