bring LSMR implementation of jutho/krylovkit.jl#46 up to date #109
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
|
There was a problem hiding this comment. Probably it would be good to have an additional mention to the docstring that LSMR requires a different interface for the function handles, since it needs both the regular as well as the adjoint action. |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,214 @@ | ||
| # reference implementation https://github.com/JuliaLinearAlgebra/IterativeSolvers.jl/blob/master/src/lsmr.jl | ||
| function linsolve(operator, b, alg::LSMR) | ||
| return linsolve(operator, b, zerovector(apply_adjoint(operator, b)), alg) | ||
| end; | ||
| function linsolve(operator, b, x₀, alg::LSMR) | ||
| u = add!!(apply_normal(operator, x₀), b, 1, -1) | ||
| β = norm(u) | ||
|
|
||
| # initialize GKL factorization | ||
| iter = GKLIterator(operator, u, alg.orth) | ||
| fact = initialize(iter; verbosity=alg.verbosity - 2) | ||
| numops = 2 | ||
| sizehint!(fact, alg.krylovdim) | ||
|
|
||
| T = eltype(fact) | ||
| Tr = real(T) | ||
| alg.conlim > 0 ? ctol = convert(Tr, inv(alg.conlim)) : ctol = zero(Tr) | ||
| istop = 0 | ||
|
|
||
| # TODO: make this an explicit copy that works with the testing datatypes | ||
| x = x₀ | ||
|
|
||
| for topit in 1:(alg.maxiter)# the outermost restart loop | ||
| # Initialize variables for 1st iteration. | ||
| α = fact.αs[end] | ||
| ζbar = α * β | ||
| αbar = α | ||
| ρ = one(Tr) | ||
| ρbar = one(Tr) | ||
| cbar = one(Tr) | ||
| sbar = zero(Tr) | ||
|
|
||
| # Initialize variables for estimation of ||r||. | ||
| βdd = β | ||
| βd = zero(Tr) | ||
| ρdold = one(Tr) | ||
| τtildeold = zero(Tr) | ||
| θtilde = zero(Tr) | ||
| ζ = zero(Tr) | ||
| d = zero(Tr) | ||
|
|
||
| # Initialize variables for estimation of ||A|| and cond(A). | ||
| normA, condA, normx = -one(Tr), -one(Tr), -one(Tr) | ||
| normA2 = abs2(α) | ||
| maxrbar = zero(Tr) | ||
| minrbar = 1e100 | ||
|
|
||
| # Items for use in stopping rules. | ||
| normb = β | ||
| normr = β | ||
| normAr = α * β | ||
|
|
||
| hbar = scale(x, zero(T)) | ||
| h = scale(fact.V[end], one(T)) | ||
|
|
||
| while length(fact) < alg.krylovdim | ||
| β = normres(fact) | ||
| fact = expand!(iter, fact) | ||
| numops += 2 | ||
|
|
||
| v = fact.V[end] | ||
| α = fact.αs[end] | ||
|
|
||
| # Construct rotation Qhat_{k,2k+1}. | ||
| αhat = hypot(αbar, alg.λ) | ||
| chat = αbar / αhat | ||
| shat = alg.λ / αhat | ||
|
|
||
| # Use a plane rotation (Q_i) to turn B_i to R_i. | ||
| ρold = ρ | ||
| ρ = hypot(αhat, β) | ||
| c = αhat / ρ | ||
| s = β / ρ | ||
| θnew = s * α | ||
| αbar = c * α | ||
|
|
||
| # Use a plane rotation (Qbar_i) to turn R_i^T to R_i^bar. | ||
| ρbarold = ρbar | ||
| ζold = ζ | ||
| θbar = sbar * ρ | ||
| ρtemp = cbar * ρ | ||
| ρbar = hypot(cbar * ρ, θnew) | ||
| cbar = cbar * ρ / ρbar | ||
| sbar = θnew / ρbar | ||
| ζ = cbar * ζbar | ||
| ζbar = -sbar * ζbar | ||
|
|
||
| # Update h, h_hat, x. | ||
| hbar = add!!(hbar, h, 1, -θbar * ρ / (ρold * ρbarold)) | ||
| h = add!!(h, v, 1, -θnew / ρ) | ||
| x = add!!(x, hbar, ζ / (ρ * ρbar), 1) | ||
|
|
||
| ############################################################################## | ||
| ## | ||
| ## Estimate of ||r|| | ||
| ## | ||
| ############################################################################## | ||
|
|
||
| # Apply rotation Qhat_{k,2k+1}. | ||
| βacute = chat * βdd | ||
| βcheck = -shat * βdd | ||
|
|
||
| # Apply rotation Q_{k,k+1}. | ||
| βhat = c * βacute | ||
| βdd = -s * βacute | ||
|
|
||
| # Apply rotation Qtilde_{k-1}. | ||
| θtildeold = θtilde | ||
| ρtildeold = hypot(ρdold, θbar) | ||
| ctildeold = ρdold / ρtildeold | ||
| stildeold = θbar / ρtildeold | ||
| θtilde = stildeold * ρbar | ||
| ρdold = ctildeold * ρbar | ||
| βd = -stildeold * βd + ctildeold * βhat | ||
|
|
||
| τtildeold = (ζold - θtildeold * τtildeold) / ρtildeold | ||
| τd = (ζ - θtilde * τtildeold) / ρdold | ||
| d += abs2(βcheck) | ||
| normr = sqrt(d + abs2(βd - τd) + abs2(βdd)) | ||
|
|
||
| # Estimate ||A||. | ||
| normA2 += abs2(β) | ||
| normA = sqrt(normA2) | ||
| normA2 += abs2(α) | ||
|
|
||
| # Estimate cond(A). | ||
| maxrbar = max(maxrbar, ρbarold) | ||
| if length(fact) > 1 | ||
| minrbar = min(minrbar, ρbarold) | ||
| end | ||
| condA = max(maxrbar, ρtemp) / min(minrbar, ρtemp) | ||
|
|
||
| ############################################################################## | ||
| ## | ||
| ## Test for convergence | ||
| ## | ||
| ############################################################################## | ||
|
|
||
| # Compute norms for convergence testing. | ||
| normAr = abs(ζbar) | ||
| normx = norm(x) | ||
|
|
||
| # Now use these norms to estimate certain other quantities, | ||
| # some of which will be small near a solution. | ||
| test1 = normr / normb | ||
| test2 = normAr / (normA * normr) | ||
| test3 = inv(condA) | ||
|
|
||
|
Comment on lines
+143
to
+148
There was a problem hiding this comment. I'm a bit confused by the placement of these variables. Is there a reason to define these before the logger, and to then use There was a problem hiding this comment. Are these different convergence criteria currently used? If not, I would probably streamline them to be compatible with what we use in the other methods, even if this "reduces" the functionality or options to specify the convergence. I will take a stab at this myself. |
||
| t1 = test1 / (one(Tr) + normA * normx / normb) | ||
| rtol = alg.btol + alg.atol * normA * normx / normb | ||
| # The following tests guard against extremely small values of | ||
| # atol, btol or ctol. (The user may have set any or all of | ||
| # the parameters atol, btol, conlim to 0.) | ||
| # The effect is equivalent to the normAl tests using | ||
| # atol = eps, btol = eps, conlim = 1/eps. | ||
|
|
||
| if alg.verbosity > 2 | ||
| msg = "LSMR linsolve in iter $topit; step $(length(fact)-1): " | ||
| msg *= "normres = " | ||
| msg *= @sprintf("%.12e", normr) | ||
| @info msg | ||
| end | ||
|
|
||
| if 1 + test3 <= 1 | ||
|
There was a problem hiding this comment. Can this be re-expressed as |
||
| istop = 6 | ||
| break | ||
| end | ||
| if 1 + test2 <= 1 | ||
| istop = 5 | ||
| break | ||
| end | ||
| if 1 + t1 <= 1 | ||
| istop = 4 | ||
| break | ||
| end | ||
| # Allow for tolerances set by the user. | ||
| if test3 <= ctol | ||
| istop = 3 | ||
| break | ||
| end | ||
| if test2 <= alg.atol | ||
| istop = 2 | ||
| break | ||
| end | ||
| if test1 <= rtol | ||
| istop = 1 | ||
| break | ||
| end | ||
| end | ||
|
|
||
| u = add!!(apply_normal(operator, x), b, 1, -1) | ||
|
|
||
| istop != 0 && break | ||
|
|
||
| #restart | ||
| β = norm(u) | ||
| iter = GKLIterator(operator, u, alg.orth) | ||
| fact = initialize!(iter, fact) | ||
| end | ||
|
|
||
| isconv = istop ∉ (0, 3, 6) | ||
| if alg.verbosity > 0 && !isconv | ||
| @warn """LSMR linsolve finished without converging after $(alg.maxiter) iterations: | ||
| * norm of residual = $(norm(u)) | ||
| * number of operations = $numops""" | ||
| elseif alg.verbosity > 0 | ||
| if alg.verbosity > 0 | ||
| @info """LSMR linsolve converged due to istop $(istop): | ||
| * norm of residual = $(norm(u)) | ||
| * number of operations = $numops""" | ||
| end | ||
| end | ||
| return (x, ConvergenceInfo(Int(isconv), u, norm(u), alg.maxiter, numops)) | ||
| end | ||
|
There was a problem hiding this comment. I guess you should probably also test a non-square matrix here, since that is the primary use-case of the method? |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I think this docstring could use some formatting improvements and a bit of a clean up 😉