Chebyshev as iterator (#145)

* Improve tests: with preconditioner and a more sane eigenspectrum for the SPD matrix (no eigenvalues close to 0 to improve convergence)

* Use the iterator idea

Author: haampie

src/chebyshev.jl

@@ -1,67 +1,129 @@
+import Base: next, start, done
+
 export chebyshev, chebyshev!
 
+type ChebyshevIterable{precT, matT, vecT, realT <: Real}
+    Pl::precT
+    A::matT
+    b::vecT
+
+    x::vecT
+    r::vecT
+    u::vecT
+    c::vecT
+
+    α::realT
+
+    λ_avg::realT
+    λ_diff::realT
+
+    resnorm::realT
+    reltol::realT
+    maxiter::Int
+    mv_products::Int
+end
+
+converged(c::ChebyshevIterable) = c.resnorm ≤ c.reltol
+start(::ChebyshevIterable) = 0
+done(c::ChebyshevIterable, iteration::Int) = iteration ≥ c.maxiter || converged(c)
+
+function next(cheb::ChebyshevIterable, iteration::Int)
+    T = eltype(cheb.u)
+
+    solve!(cheb.c, cheb.Pl, cheb.r)
+
+    if iteration == 1
+        cheb.α = T(2) / cheb.λ_avg
+        copy!(cheb.u, cheb.c)
+    else
+        β = (cheb.λ_diff * cheb.α / 2) ^ 2
+        cheb.α = inv(cheb.λ_avg - β)
+
+        # Almost an axpy u = c + βu
+        scale!(cheb.u, β)
+        @blas! cheb.u += one(T) * cheb.c
+    end
+
+    A_mul_B!(cheb.c, cheb.A, cheb.u)
+    cheb.mv_products += 1
+
+    @blas! cheb.x += cheb.α * cheb.u
+    @blas! cheb.r -= cheb.α * cheb.c
+
+    cheb.resnorm = norm(cheb.r)
+
+    cheb.resnorm, iteration + 1
+end
+
+chebyshev_iterable(A, b, λmin::Real, λmax::Real; kwargs...) =
+    chebyshev_iterable!(zerox(A, b), A, b, λmin, λmax; kwargs...)
+
+function chebyshev_iterable!(x, A, b, λmin::Real, λmax::Real;
+    tol = sqrt(eps(real(eltype(b)))), maxiter = size(A, 1), Pl = Identity(), initially_zero = false)
+
+    λ_avg = (λmax + λmin) / 2
+    λ_diff = (λmax - λmin) / 2
+
+    T = eltype(b)
+    r = copy(b)
+    u = zeros(x)
+    c = similar(x)
+
+    # One MV product less
+    if initially_zero
+        resnorm = norm(r)
+        reltol = tol * resnorm
+        mv_products = 0
+    else
+        A_mul_B!(c, A, x)
+        @blas! r -= one(T) * c
+        resnorm = norm(r)
+        reltol = tol * norm(b)
+        mv_products = 1
+    end
+
+    ChebyshevIterable(Pl, A, b,
+        x, r, u, c,
+        zero(real(T)),
+        λ_avg, λ_diff,
+        resnorm, reltol, maxiter, mv_products
+    )
+end
+
 ####################
 # API method calls #
 ####################
 
 chebyshev(A, b, λmin::Real, λmax::Real; kwargs...) =
-    chebyshev!(zerox(A, b), A, b, λmin, λmax; kwargs...)
+    chebyshev!(zerox(A, b), A, b, λmin, λmax; initially_zero = true, kwargs...)
 
 function chebyshev!(x, A, b, λmin::Real, λmax::Real;
-    n::Int=size(A,2), tol::Real = sqrt(eps(typeof(real(b[1])))),
-    maxiter::Int = n^3, log::Bool=false, kwargs...
-    )
-	K = KrylovSubspace(A, n, 1, Adivtype(A, b))
-	init!(K, x)
+    Pl = Identity(),
+    tol::Real=sqrt(eps(real(eltype(b)))),
+    maxiter::Int=size(A, 1),
+    log::Bool=false,
+    verbose::Bool=false,
+    initially_zero::Bool=false
+)
     history = ConvergenceHistory(partial=!log)
     history[:tol] = tol
-    reserve!(history,:resnorm,maxiter)
-    chebyshev_method!(history, x, K, b, λmin, λmax; tol=tol, maxiter=maxiter, kwargs...)
-    log && shrink!(history)
-    log ? (x, history) : x
-end
+    reserve!(history, :resnorm, maxiter)
 
-#########################
-# Method Implementation #
-#########################
+    verbose && @printf("=== chebyshev ===\n%4s\t%7s\n","iter","resnorm")
 
-function chebyshev_method!(
-    log::ConvergenceHistory, x, K::KrylovSubspace, b, λmin::Real, λmax::Real;
-    Pr = 1, tol::Real = sqrt(eps(typeof(real(b[1])))), maxiter::Int = K.n^3,
-    verbose::Bool=false
-    )
+    iterable = chebyshev_iterable!(x, A, b, λmin, λmax; tol=tol, maxiter=maxiter, Pl=Pl, initially_zero=initially_zero)
+    history.mvps = iterable.mv_products
+    for (iteration, resnorm) = enumerate(iterable)
+        nextiter!(history)
+        history.mvps = iterable.mv_products
+        push!(history, :resnorm, resnorm)
+        verbose && @printf("%3d\t%1.2e\n", iteration, resnorm)
+    end
+    verbose && println()
+    setconv(history, converged(iterable))
+    log && shrink!(history)
 
-    verbose && @printf("=== chebyshev ===\n%4s\t%7s\n","iter","resnorm")
-	local α, p
-	K.order = 1
-    tol = tol*norm(b)
-    log.mvps=1
-	r = b - nextvec(K)
-	d::eltype(b) = (λmax + λmin)/2
-	c::eltype(b) = (λmax - λmin)/2
-	for iter = 1:maxiter
-        nextiter!(log, mvps=1)
-		z = solve(Pr,r)
-		if iter == 1
-			p = z
-			α = 2/d
-		else
-			β = (c*α/2)^2
-			α = 1/(d - β)
-			p = z + β*p
-		end
-		append!(K, p)
-		@blas! x += α*p
-		@blas! r -= α*nextvec(K)
-		#Check convergence
-		resnorm = norm(r)
-        push!(log, :resnorm, resnorm)
-        verbose && @printf("%3d\t%1.2e\n",iter,resnorm)
-        resnorm < tol && break
-	end
-    verbose && @printf("\n")
-    setconv(log, 0<=norm(r)<tol)
-	x
+    log ? (x, history) : x
 end
 
 #################
@@ -108,7 +170,7 @@ $arg
 
 ## Keywords
 
-`Pr = 1`: right preconditioner of the method.
+`Pl = 1`: left preconditioner of the method.
 
 `tol::Real = sqrt(eps())`: stopping tolerance.
 

test/chebyshev.jl

@@ -5,24 +5,41 @@ using LinearMaps
 
 srand(1234321)
 
+function randSPD(T, n)
+    A = rand(T, n, n) + n * eye(T, n)
+    A = A' + A
+    A = A' * A
+end
+
 #Chebyshev
 facts("Chebyshev") 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'
-    A=A'*A #Construct SPD matrix
-    b=convert(Vector{T}, randn(n))
-    T<:Complex && (b+=convert(Vector{T}, im*randn(n)))
-    b=b/norm(b)
-    tol = 0.1 #For some reason Chebyshev is very slow
-    v = eigvals(A)
-    mxv = maximum(v)
-    mnv = minimum(v)
-    x_cheby, c_cheby= chebyshev(A, b, mxv+(mxv-mnv)/100, mnv-(mxv-mnv)/100, tol=tol, maxiter=10^5, log=true)
-    @fact c_cheby.isconverged --> true
-    @fact norm(A*x_cheby-b) --> less_than(tol)
+        A = randSPD(T, n)
+        b = rand(T, n)
+        b /= norm(b)
+
+        tol = sqrt(eps(real(T)))
+        mnv, mxv = eigmin(A), eigmax(A)
+        Δ = (mxv - mnv) / 100
+
+        x, history = chebyshev(A, b, mnv - Δ, mxv + Δ, tol=tol, maxiter=10n, log=true)
+        @fact history.isconverged --> true
+        @fact norm(A * x - b) --> less_than(tol)
+
+        context("Preconditioned") do
+            B = randSPD(T, n)
+            B_fact = cholfact!(B)
+            BA = B_fact \ A
+            λs = eigvals(BA)
+            mnv, mxv = minimum(real(λs)), maximum(real(λs))
+            Δ = (mxv - mnv) / 100
+
+            x, history = chebyshev(A, b, mnv - Δ, mxv + Δ, Pl = B_fact, tol=tol, maxiter=10n, log=true)
+            @fact history.isconverged --> true
+            @fact norm(A * x - b) --> less_than(tol)
+        end
     end
 end
 end