Power method iterators (#139)

* Power method as iterator

* Some doc fixes and getting around the lufact bug in base

* Improve presentation & docs

Author: haampie

src/simple.jl

@@ -1,96 +1,115 @@
+import Base: start, next, done
+
 #Simple methods
 export powm, invpowm
 
-####################
-# API method calls #
-####################
+type PowerMethodIterable{matT, vecT <: AbstractVector, numT <: Number, eigvalT <: Number}
+    A::matT
+    x::vecT
+    tol::numT
+    maxiter::Int
+    θ::eigvalT
+    r::vecT
+    Ax::vecT
+    residual::numT
+end
 
-function powm(A;
-    x=nothing, tol::Real=eps(real(eltype(A)))*size(A,2)^3, maxiter::Int=size(A,2),
-    log::Bool=false, kwargs...
-    )
-    K = KrylovSubspace(A, 1)
-    x==nothing ? initrand!(K) : init!(K, x/norm(x))
 
-    history = ConvergenceHistory(partial=!log)
-    history[:tol] = tol
-    reserve!(history,:resnorm, maxiter)
-    eig, v = powm_method!(history, K; tol=tol, maxiter=maxiter, kwargs...)
-    log && shrink!(history)
-    log ? (eig, v, history) : (eig, v)
+##
+## Iterators
+##
+
+@inline converged(p::PowerMethodIterable) = p.residual ≤ p.tol
+
+@inline start(p::PowerMethodIterable) = 0
+
+@inline done(p::PowerMethodIterable, iteration::Int) = iteration > p.maxiter || converged(p)
+
+function next(p::PowerMethodIterable, iteration::Int)
+
+    A_mul_B!(p.Ax, p.A, p.x)
+
+    # Rayleigh quotient θ = x'Ax
+    p.θ = dot(p.x, p.Ax)
+
+    # (Previous) residual vector r = Ax - λx
+    copy!(p.r, p.Ax)
+    BLAS.axpy!(-p.θ, p.x, p.r)
+
+    # Normed residual
+    p.residual = norm(p.r)
+
+    # Normalize the next approximation
+    copy!(p.x, p.Ax)
+    scale!(p.x, one(eltype(p.x)) / norm(p.x))
+
+    p.residual, iteration + 1
 end
 
-#########################
-# Method Implementation #
-#########################
-
-function powm_method!{T}(log::ConvergenceHistory, K::KrylovSubspace{T};
-    tol::Real=eps(real(T))*K.n^3, maxiter::Int=K.n, verbose::Bool=false
-    )
-    verbose && @printf("=== powm ===\n%4s\t%7s\n","iter","resnorm")
-    θ = zero(T)
-    v = Array(T, K.n)
-    for iter=1:maxiter
-        nextiter!(log,mvps=1)
-        v = lastvec(K)
-        y = nextvec(K)
-        θ = dot(v, y)
-        resnorm = real(norm(y - θ*v))
-        push!(log, :resnorm, resnorm)
-        verbose && @printf("%3d\t%1.2e\n",iter,resnorm)
-        resnorm <= tol*abs(θ) && (setconv(log, resnorm >= 0); break)
-        appendunit!(K, y)
-    end
-    verbose && @printf("\n")
-    θ, v
+# Transforms the eigenvalue back whether shifted or inversed
+@inline transform_eigenvalue(θ, inverse::Bool, σ) = σ + (inverse ? inv(θ) : θ)
+
+function powm_iterable!(A, x; tol = eps(real(eltype(A))) * size(A, 2) ^ 3, maxiter = size(A, 1))
+    T = eltype(x)
+    PowerMethodIterable(A, x, tol, maxiter, zero(T), similar(x), similar(x), realmax(real(T)))
+end
+
+function powm_iterable(A; kwargs...)
+    x0 = rand(Complex{real(eltype(A))}, size(A, 1))
+    scale!(x0, one(eltype(A)) / norm(x0))
+    powm_iterable!(A, x0; kwargs...)
 end
 
 ####################
 # API method calls #
 ####################
 
-function invpowm(A;
-    x=nothing, shift::Number=0, tol::Real=eps(real(eltype(A)))*size(A,2)^3,
-    maxiter::Int=size(A,2), log::Bool=false, kwargs...
-    )
-    K = KrylovSubspace(A, 1)
-    x==nothing ? initrand!(K) : init!(K, x/norm(x))
+function powm(A; kwargs...)
+    x0 = rand(Complex{real(eltype(A))}, size(A, 1))
+    scale!(x0, one(eltype(A)) / norm(x0))
+    powm!(A, x0; kwargs...)
+end
 
-    history = ConvergenceHistory(partial=!log)
+function powm!(A, x;
+    tol = eps(real(eltype(A))) * size(A, 2) ^ 3,
+    maxiter = size(A, 1),
+    shift = zero(eltype(A)),
+    inverse::Bool = false,
+    log::Bool = false,
+    verbose::Bool = false
+)
+    history = ConvergenceHistory(partial = !log)
     history[:tol] = tol
-    reserve!(history,:resnorm, maxiter)
-    eig, v = invpowm_method!(history, K, shift; tol=tol, maxiter=maxiter, kwargs...)
+    reserve!(history, :resnorm, maxiter)
+    verbose && @printf("=== powm ===\n%4s\t%7s\n", "iter", "resnorm")
+
+    iterable = powm_iterable!(A, x, tol = tol, maxiter = maxiter)
+
+    for (iteration, residual) = enumerate(iterable)
+        nextiter!(history, mvps = 1)
+        verbose && @printf("%3d\t%1.2e\n", iteration, residual)
+    end
+
+    setconv(history, converged(iterable))
+
+    verbose && println()
+
     log && shrink!(history)
-    log ? (eig, v, history) : (eig, v)
+
+    λ = transform_eigenvalue(iterable.θ, inverse, shift)
+    x = iterable.x
+
+    log ? (λ, x, history) : (λ, x)
 end
 
-#########################
-# Method Implementation #
-#########################
-
-function invpowm_method!{T}(log::ConvergenceHistory, K::KrylovSubspace{T}, σ::Number=zero(T);
-    tol::Real=eps(real(T))*K.n^3, maxiter::Int=K.n, verbose::Bool=false
-    )
-    verbose && @printf("=== invpowm ===\n%4s\t%7s\n","iter","resnorm")
-    θ = zero(T)
-    v = Array(T, K.n)
-    y = Array(T, K.n)
-    σ = convert(T, σ)
-    for iter=1:maxiter
-        nextiter!(log,mvps=1)
-        v = lastvec(K)
-        y = (K.A-σ*eye(K))\v
-        θ = dot(v, y)
-        resnorm = norm(y - θ*v)
-        push!(log, :resnorm, resnorm)
-        verbose && @printf("%3d\t%1.2e\n",iter,resnorm)
-        resnorm <= tol*abs(θ) && (setconv(log, resnorm >= 0); break)
-        appendunit!(K, y)
-    end
-    verbose && @printf("\n")
-    σ+1/θ, y/θ
+function invpowm(B; kwargs...)
+    x0 = rand(Complex{real(eltype(B))}, size(B, 1))
+    scale!(x0, one(eltype(B)) / norm(x0))
+    invpowm!(B, x0; kwargs...)
 end
 
+invpowm!(B, x0; kwargs...) = powm!(B, x0; inverse = true, kwargs...)
+
 #################
 # Documentation #
 #################
@@ -99,13 +118,25 @@ let
 #Initialize parameters
 doc1_call = """    powm(A)
 """
-doc2_call = """    invpowm(A)
+doc2_call = """    invpowm(B)
 """
-doc1_msg = """Find biggest eigenvalue of `A` and its associated eigenvector
+doc1_msg = """Find the largest eigenvalue `λ` (in absolute value) of `A` and its associated eigenvector `x`
 using the power method.
 """
-doc2_msg = """Find closest eigenvalue of `A` to `shift` and its associated eigenvector
-using the inverse power iteration method.
+doc2_msg = """For an eigenvalue problem Ax = λx, find the closest eigenvalue in the complex plane to `shift`
+together with its associated eigenvector `x`. The first argument `B` should be a mapping that has
+the effect of B * x = inv(A - shift * I) * x and should support the `A_mul_B!` operation.
+
+# Examples
+
+```julia
+using LinearMaps
+σ = 1.0 + 1.3im
+A = rand(Complex128, 50, 50)
+F = lufact(A - UniformScaling(σ))
+Fmap = LinearMap((y, x) -> A_ldiv_B!(y, F, x), 50, Complex128, ismutating = true)
+λ, x = invpowm(Fmap, shift = σ, tol = 1e-4, maxiter = 200)
+```
 """
 doc1_karg = ""
 doc2_karg = "`shift::Number=0`: shift to be applied to matrix A."
@@ -124,13 +155,11 @@ $call
 
 $msg
 
-If `log` is set to `true` is given, method will output a tuple `eig, v, ch`. Where
-`ch` is a `ConvergenceHistory` object. Otherwise it will only return `eig, v`.
+If `log` is set to `true` is given, method will output a tuple `λ, x, ch`. Where
+`ch` is a `ConvergenceHistory` object. Otherwise it will only return `λ, x`.
 
 # Arguments
 
-`K::KrylovSubspace`: krylov subspace.
-
 `A`: linear operator.
 
 ## Keywords
@@ -152,15 +181,15 @@ containing extra information of the method execution.
 
 **if `log` is `false`**
 
-`eig::Real`: eigen value
+`λ::Number`: eigenvalue
 
-`v::Vector`: eigen vector
+`x::Vector`: eigenvector
 
 **if `log` is `true`**
 
-`eig::Real`: eigen value
+`eig::Real`: eigenvalue
 
-`v::Vector`: eigen vector
+`x::Vector`: eigenvector
 
 `ch`: convergence history.
 

test/simple_eigensolvers.jl

@@ -9,32 +9,44 @@ srand(1234321)
 # Eigensystem solvers #
 #######################
 
-facts("simple eigensolvers") do
+facts("Simple eigensolvers") do
+
+n = 10
+
 for T in (Float32, Float64, Complex64, Complex128)
+
     context("Matrix{$T}") do
-    A=convert(Matrix{T}, randn(n,n))
-    T<:Complex && (A+=convert(Matrix{T}, im*randn(n,n)))
-    A=A+A' #Symmetric/Hermitian
 
-    tol = (eltype(T) <: Complex ?2:1)*n^2*cond(A)*eps(real(one(T)))
-    v = eigvals(A)
+    A = rand(T, n, n)
+    A = A' * A
+    λs = eigvals(A)
+
+    tol = (eltype(T) <: Complex ? 2 : 1) * n^2 * cond(A) * eps(real(one(T)))
 
     ## Simple methods
 
     context("Power iteration") do
-    eval_big = maximum(v) > abs(minimum(v)) ? maximum(v) : minimum(v)
-    eval_pow = powm(A; tol=sqrt(eps(real(one(T)))), maxiter=2000)[1]
-    @fact norm(eval_big-eval_pow) --> less_than(tol)
+        λ, x = powm(A; tol = tol, maxiter = 10n)
+        @fact λs[end] --> roughly(λ)
+        @fact norm(A * x - λ * x) --> less_than(tol)
     end
 
     context("Inverse iteration") do
-    irnd = ceil(Int, rand()*(n-2))
-    eval_rand = v[1+irnd] #Pick random eigenvalue
-    # Perturb the eigenvalue by < 1/4 of the distance to the nearest eigenvalue
-    eval_diff = min(abs(v[irnd]-eval_rand), abs(v[irnd+2]-eval_rand))
-    σ = eval_rand + eval_diff/2*(rand()-.5)
-    eval_ii = invpowm(A; shift=σ, tol=sqrt(eps(real(one(T)))), maxiter=2000)[1]
-    @fact norm(eval_rand-eval_ii) --> less_than(tol)
+        # Set a target near the middle eigenvalue
+        idx = div(n, 2)
+        σ = T(0.75 * λs[idx] + 0.25 * λs[idx + 1])
+
+        # Construct F = inv(A - σI) "matrix free"
+        # Make sure we use complex arithmetic everywhere,
+        # because of the follow bug in base: https://github.com/JuliaLang/julia/issues/22683
+        F = lufact(complex(A) - UniformScaling(σ))
+        Fmap = LinearMap((y, x) -> A_ldiv_B!(y, F, x), size(A, 1), complex(T), ismutating = true)
+
+        λ, x = invpowm(Fmap; shift = σ, tol = tol, maxiter = 10n)
+
+        @fact norm(A * x - λ * x) --> less_than(tol)
+        @fact λ --> roughly(λs[idx])
+
     end
 
     end