Madelaine, I'd be interested in hearing more about how do you plan to solve your problem in Matlab if B is singular. The backtraces that Matlab's warnings and errors produce suggest that Matlab also factors B. I'm not really familiar with the methods for the generalized eigenvalue problem when B is only semidefinite and the ARPACK manual seems to have lost the math. The "note that [blank] thus [blank]" leaves much to the reader. See
http://www.caam.rice.edu/software/ARPACK/UG/node53.html#SECTION00851000000000000000 On Saturday, August 6, 2016 at 5:32:55 PM UTC-4, Madeleine Udell wrote: > > Oh gosh! It's terrible that we're trying to factor the second matrix. It's > quite common that the second matrix in a generalized eigenvalue problem is > not positive definite. > > If we can't get this fixed soon, I'm going to have to resort to Matlab for > my current project...! > > On Sat, Aug 6, 2016 at 1:53 PM, Ralph Smith <[email protected] > <javascript:>> 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 >>> >> > > > -- > Madeleine Udell > Assistant Professor, Operations Research and Information Engineering > Cornell University > https://people.orie.cornell.edu/mru8/ > (415) 729-4115 >
