Complete eigs implementation

[ViralBShah]

The eigs and svds implementation needs the following features to be completed:

• Support for the generalized eigenvalue problem
• Implement shifts (We now have dense and sparse direct solvers)
• Implement svds for the case where m < n
• Call eig and svd in cases where the user asks for all eigen/singular values
• Allow the user to specify a general linear operator instead of a matrix. Initial support for this is already in place
• Needs many more tests, and more comprehensive documentation

It has also been discussed in #1573 that this functionality should be merged into eig and svd. Once these are feature complete, we can unify the API.

[andrioni]

Allow the user to specify a general linear operator instead of a matrix. Initial support for this is already in place

What exactly does this mean? Are there any specific types other than matrices for this?

[StefanKarpinski]

Matlab's eigs and svds functions can work with arbitrary function handles that compute linear operators.

[stevengj]

@andrioni, what it means is that the user specifies a "black-box" function A(x) that computes A*x. This is useful in a variety of circumstances in which the matrix is stored (explicitly or implicitly) in some data structure that can be multiplied quickly by a vector, but isn't a standard matrix datastructure. For example, a circulant matrix can be stored implicitly by the FFT of its first row, and multipled in O(n log n) time by ifft(fft(x) .* fft(row)).

[StefanKarpinski]

Right – @stevengj's explanation is much better than mine :-)

[jtravs]

Did anyone have a go at implementing the shifts part? I dearly need this. I will try and have a look myself, but I'm afraid I am far from an expert on ARPACK.

[ViralBShah]

Not yet. This is not too difficult - and I can try. Do you need eigs or svds?

[jtravs]

What I really want is the sigma argument like in the SciPy or Matlab eigs functions, so that I can compute interior eigenvalues. I have been reading the SciPy source code (their eigs is faster than the Matlab one for the same problem which is weird as I though they both used ARPACK), but didn't get far yet:
https://github.com/scipy/scipy/blob/master/scipy/sparse/linalg/eigen/arpack/arpack.py

[ViralBShah]

The implementation goes into https://github.com/JuliaLang/julia/blob/master/base/linalg/arnoldi.jl

All the wrappers are already provided in arpack.jl. Basically, look at your deps/arpack-ng-3.1.3/EXAMPLES/COMPLEX/zndrv2.f for a driver example, and then replicate the flow in julia.

[jtravs]

OK, I'll have a go at it over the next few days. I'll may need some guidance later on.

[ViralBShah]

eigs is complete enough to close this. svds is discussed in JuliaLang/LinearAlgebra.jl#106

[mcbal]

What about the use of a general linear operator instead of a matrix? Basically, this requires passing a function Afun::Function to eigs, together with an integer n reflecting the size of the matrix representation A of the operator (similar to how it's done in MATLAB's eigs). The extension is straightforward, but I am not sure how it can be implemented in an elegant way, given that the properties of A cannot be deduced from the function A(x) alone and should be specified separately.

[stevengj]

@mcbal, this is accomplished by duck-typing (e.g. JuliaLinearAlgebra/IterativeSolvers.jl#2). The argument type of eigs is not specified, so you can pass any object that supports the requisite methods (I think *, size, and \ for inverse iteration size, issym, eltype, and *). To define your own linear operator, you just create a new type and give it these methods.

But this isn't really documented, and it is not as clean as it could be.

[mcbal]

Thank you very much for this Julia-esque solution!

[natschil]

Is there any documentation anywhere as to what functions exactly need to be implemented for this to work? The following example currently does not work, whereas I think it should:

import Base: *,size,transpose,issymmetric,eltype
type myLinearOperator
    m::Int
end
*(op::myLinearOperator, x::Array{Float64}) = x
size(op::myLinearOperator) = (op.m,op.m)
transpose(op::myLinearOperator) = op
ctranspose(op::myLinearOperator) = op
issymmetric(op::myLinearOperator) = false
eltype(op::myLinearOperator) = Float64
myoperator = myLinearOperator(10)
myArray = Array{Float64}(ones(10))
eigs(myoperator,v0=myArray)

The error I get from this is:

MethodError: no method matching A_mul_B!(::SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}, ::myLinearOperator, ::SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true})