Merge pull request #6053 from JuliaLang/pr/5776

Default to shift-and-invert mode for eigs for numeric sigma

NEWS.md

@@ -158,6 +158,8 @@ Library improvements
 
       * `sparse(A) \ B` now supports a matrix `B` of right-hand sides ([#5196]).
 
+      * `eigs(A, sigma)` now uses shift-and-invert for nonzero shifts `sigma` and inverse iteration for `which="SM"`. If `sigma==nothing` (the new default), computes ordinary (forward) iterations. ([#5776])
+
     * Dense linear algebra for special matrix types
 
       * Interconversions between the special matrix types `Diagonal`, `Bidiagonal`,
@@ -290,6 +292,7 @@ Deprecated or removed
 [#5427]: https://github.com/JuliaLang/julia/pull/5427
 [#5748]: https://github.com/JuliaLang/julia/issues/5748
 [#5511]: https://github.com/JuliaLang/julia/issues/5511
+[#5776]: https://github.com/JuliaLang/julia/issues/5776
 [#5819]: https://github.com/JuliaLang/julia/issues/5819
 [#4871]: https://github.com/JuliaLang/julia/issues/4871
 [#4996]: https://github.com/JuliaLang/julia/issues/4996

base/linalg/arnoldi.jl

@@ -3,7 +3,7 @@ using .ARPACK
 ## eigs
 
 function eigs(A;nev::Integer=6, ncv::Integer=20, which::ASCIIString="LM",
-          tol=0.0, maxiter::Integer=1000, sigma=0,v0::Vector=zeros((0,)),
+          tol=0.0, maxiter::Integer=1000, sigma=nothing, v0::Vector=zeros((0,)),
           ritzvec::Bool=true, complexOP::Bool=false)
     n = chksquare(A)
     (n <= 6) && (nev = min(n-1, nev))
@@ -13,18 +13,20 @@ function eigs(A;nev::Integer=6, ncv::Integer=20, which::ASCIIString="LM",
     T     = eltype(A)
     cmplx = T<:Complex
     bmat  = "I"
+
     if !isempty(v0)
         length(v0)==n || throw(DimensionMismatch(""))
         eltype(v0)==T || error("Starting vector must have eltype $T")
     else
         v0=zeros(T,(0,))
     end
 
-    if sigma == 0
+    if sigma===nothing # normal iteration
         mode = 1
+        sigma = 0
         linop(x) = A * x
-    else
-        C = lufact(A - sigma*eye(A))
+    else # inverse iteration with level shift
+        C = lufact(sigma==0 ? A : A - sigma*eye(A))
         if cmplx
             mode = 3
             linop(x) = C\x

doc/stdlib/base.rst

@@ -3561,6 +3561,10 @@ Constructors
 
    m-by-n identity matrix
 
+.. function:: eye(A)
+
+   Constructs an identity matrix of the same dimensions and type as ``A``.
+
 .. function:: linspace(start, stop, n)
 
    Construct a vector of ``n`` linearly-spaced elements from ``start`` to ``stop``.

doc/stdlib/linalg.rst

@@ -383,19 +383,43 @@ Linear algebra functions in Julia are largely implemented by calling functions f
 
    The conjugate transposition operator (``'``).
 
-.. function:: eigs(A; nev=6, which="LM", tol=0.0, maxiter=1000, sigma=0, ritzvec=true, op_part=:real,v0=zeros((0,))) -> (d,[v,],nconv,niter,nmult,resid)
+.. function:: eigs(A; nev=6, which="LM", tol=0.0, maxiter=1000, sigma=nothing, ritzvec=true, op_part=:real,v0=zeros((0,))) -> (d,[v,],nconv,niter,nmult,resid)
 
-   ``eigs`` computes eigenvalues ``d`` of A using Arnoldi factorization. The following keyword arguments are supported:
+   ``eigs`` computes eigenvalues ``d`` of ``A`` using Lanczos or Arnoldi iterations for real symmetric or general nonsymmetric matrices respectively. The following keyword arguments are supported:
     * ``nev``: Number of eigenvalues
-    * ``which``: type of eigenvalues ("LM", "SM")
+    * ``which``: type of eigenvalues to compute. See the note below.
+
+      ========= ======================================================================================================================
+      ``which`` type of eigenvalues
+      --------- ----------------------------------------------------------------------------------------------------------------------
+      ``"LM"``  eigenvalues of largest magnitude
+      ``"SM"``  eigenvalues of smallest magnitude
+      ``"LA"``  largest algebraic eigenvalues (real symmetric ``A`` only)
+      ``"SA"``  smallest algebraic eigenvalues (real symmetric ``A`` only)
+      ``"BE"``  compute half of the eigenvalues from each end of the spectrum, biased in favor of the high end. (symmetric ``A`` only)
+      ``"LR"``  eigenvalues of largest real part (nonsymmetric ``A`` only)
+      ``"SR"``  eigenvalues of smallest real part (nonsymmetric ``A`` only)
+      ``"LI"``  eigenvalues of largest imaginary part (nonsymmetric ``A`` only)
+      ``"SI"``  eigenvalues of smallest imaginary part (nonsymmetric ``A`` only)
+      ========= ======================================================================================================================
+
     * ``tol``: tolerance (:math:`tol \le 0.0` defaults to ``DLAMCH('EPS')``)
     * ``maxiter``: Maximum number of iterations
-    * ``sigma``: find eigenvalues close to ``sigma`` using shift and invert
+    * ``sigma``: Specifies the level shift used in inverse iteration. If ``nothing`` (default), defaults to ordinary (forward) iterations. Otherwise, find eigenvalues close to ``sigma`` using shift and invert iterations.
     * ``ritzvec``: Returns the Ritz vectors ``v`` (eigenvectors) if ``true``
-    * ``op_part``: which part of linear operator to use for real A (:real, :imag)
-    * ``v0``: starting vector from which to start the Arnoldi iteration
+    * ``op_part``: which part of linear operator to use for real ``A`` (``:real``, ``:imag``)
+    * ``v0``: starting vector from which to start the iterations
+
    ``eigs`` returns the ``nev`` requested eigenvalues in ``d``, the corresponding Ritz vectors ``v`` (only if ``ritzvec=true``), the number of converged eigenvalues ``nconv``, the number of iterations ``niter`` and the number of matrix vector multiplications ``nmult``, as well as the final residual vector ``resid``.
-
+
+   .. note:: The ``sigma`` and ``which`` keywords interact: the description of eigenvalues searched for by ``which`` do _not_ necessarily refer to the eigenvalues of ``A``, but rather the linear operator constructed by the specification of the iteration mode implied by ``sigma``. 
+
+      =============== ================================== ==================================
+      ``sigma``       iteration mode                     ``which`` refers to eigenvalues of
+      --------------- ---------------------------------- ----------------------------------
+      ``nothing``     ordinary (forward)                 :math:`A`
+      real or complex inverse with level shift ``sigma`` :math:`(A - \sigma I )^{-1}`
+      =============== ================================== ==================================
 
 .. function:: svds(A; nev=6, which="LA", tol=0.0, maxiter=1000, ritzvec=true)
 

test/arpack.jl

@@ -1,4 +1,4 @@
-begin
+    begin
     local n,a,asym,d,v
     n = 10
     areal  = randn(n,n)
@@ -24,6 +24,11 @@ begin
 	(d,v) = eigs(apd, nev=3)
 	@test_approx_eq apd*v[:,3] d[3]*v[:,3]
 #        @test_approx_eq eigs(apd; nev=1, sigma=d[3])[1][1] d[3]
+
+    # test (shift-and-)invert mode
+    (d,v) = eigs(apd, nev=3, sigma=0)
+    @test_approx_eq apd*v[:,3] d[3]*v[:,3]
+
     end
 end