From 59a3049c53fae95694b1eaa480b731b1ad84888a Mon Sep 17 00:00:00 2001 From: nHackel Date: Fri, 23 Aug 2024 10:52:56 +0200 Subject: [PATCH 1/2] Fix verbose printing in ADMM, FISTA, OptISTA, and POGM solvers --- src/ADMM.jl | 10 +++++----- src/FISTA.jl | 4 ++-- src/OptISTA.jl | 2 +- src/POGM.jl | 4 ++-- test/testSolvers.jl | 20 ++++++++++++++++++++ 5 files changed, 30 insertions(+), 10 deletions(-) diff --git a/src/ADMM.jl b/src/ADMM.jl index a837ade0..e4cd6722 100644 --- a/src/ADMM.jl +++ b/src/ADMM.jl @@ -204,7 +204,7 @@ performs one ADMM iteration. """ function iterate(solver::ADMM, state::ADMMState) done(solver, state) && return nothing - solver.verbose && println("Outer ADMM Iteration #$iteration") + solver.verbose && println("Outer ADMM Iteration #$(state.iteration)") # 1. solve arg min_x 1/2|| Ax-b ||² + ρ/2 Σ_i||Φi*x+ui-zi||² # <=> (A'A+ρ Σ_i Φi'Φi)*x = A'b+ρΣ_i Φi'(zi-ui) @@ -264,10 +264,10 @@ function iterate(solver::ADMM, state::ADMMState) end if solver.verbose - println("rᵏ[$i]/ɛᵖʳⁱ[$i] = $(solver.rᵏ[i]/solver.ɛᵖʳⁱ[i])") - println("sᵏ[$i]/ɛᵈᵘᵃ[$i] = $(solver.sᵏ[i]/solver.ɛᵈᵘᵃ[i])") - println("Δ[$i]/Δᵒˡᵈ[$i] = $(solver.Δ[i]/Δᵒˡᵈ)") - println("new ρ[$i] = $(solver.ρ[i])") + println("rᵏ[$i]/ɛᵖʳⁱ[$i] = $(state.rᵏ[i]/state.ɛᵖʳⁱ[i])") + println("sᵏ[$i]/ɛᵈᵘᵃ[$i] = $(state.sᵏ[i]/state.ɛᵈᵘᵃ[i])") + println("Δ[$i]/Δᵒˡᵈ[$i] = $(state.Δ[i]/Δᵒˡᵈ)") + println("new ρ[$i] = $(state.ρ[i])") flush(stdout) end end diff --git a/src/FISTA.jl b/src/FISTA.jl index db4d4cf4..cad3254d 100644 --- a/src/FISTA.jl +++ b/src/FISTA.jl @@ -154,7 +154,7 @@ function iterate(solver::FISTA, state::FISTAState) state.x .-= state.ρ .* state.res state.rel_res_norm = norm(state.res) / state.norm_x₀ - solver.verbose && println("Iteration $iteration; rel. residual = $(state.rel_res_norm)") + solver.verbose && println("Iteration $(state.iteration); rel. residual = $(state.rel_res_norm)") # the two lines below are equivalent to the ones above and non-allocating, but require the 5-argument mul! function to implemented for AHA, i.e. if AHA is LinearOperator, it requires LinearOperators.jl v2 # mul!(solver.x, solver.AHA, solver.xᵒˡᵈ, -solver.ρ, 1) @@ -170,7 +170,7 @@ function iterate(solver::FISTA, state::FISTAState) # gradient restart conditions if solver.restart == :gradient if real(state.res ⋅ (state.x .- state.xᵒˡᵈ) ) > 0 #if momentum is at an obtuse angle to the negative gradient - solver.verbose && println("Gradient restart at iter $iteration") + solver.verbose && println("Gradient restart at iter $(state.iteration)") state.theta = 1 end end diff --git a/src/OptISTA.jl b/src/OptISTA.jl index d835f0e3..33160191 100644 --- a/src/OptISTA.jl +++ b/src/OptISTA.jl @@ -184,7 +184,7 @@ function iterate(solver::OptISTA, state::OptISTAState) state.y .-= state.ρ * state.γ .* state.res state.rel_res_norm = norm(state.res) / state.norm_x₀ - solver.verbose && println("Iteration $iteration; rel. residual = $(state.rel_res_norm)") + solver.verbose && println("Iteration $(state.iteration); rel. residual = $(state.rel_res_norm)") # proximal map prox!(solver.reg, state.y, state.ρ * state.γ * λ(solver.reg)) diff --git a/src/POGM.jl b/src/POGM.jl index c5141816..ffc3ac74 100644 --- a/src/POGM.jl +++ b/src/POGM.jl @@ -183,7 +183,7 @@ function iterate(solver::POGM, state::POGMState) state.x .-= state.ρ .* state.res state.rel_res_norm = norm(state.res) / state.norm_x₀ - solver.verbose && println("Iteration $iteration; rel. residual = $(state.rel_res_norm)") + solver.verbose && println("Iteration $(state.iteration); rel. residual = $(state.rel_res_norm)") # inertial parameters state.thetaᵒˡᵈ = state.theta @@ -222,7 +222,7 @@ function iterate(solver::POGM, state::POGMState) if solver.restart == :gradient state.w .+= state.y .+ state.ρ ./ state.γ .* (state.x .- state.z) if real((state.w ⋅ state.x - state.w ⋅ state.z) / state.γ - state.w ⋅ state.res) < 0 - solver.verbose && println("Gradient restart at iter $iteration") + solver.verbose && println("Gradient restart at iter $(state.iteration)") state.σ = 1 state.theta = 1 else # decreasing γ diff --git a/test/testSolvers.jl b/test/testSolvers.jl index 09a9bb84..68eb41da 100644 --- a/test/testSolvers.jl +++ b/test/testSolvers.jl @@ -217,6 +217,25 @@ function testConvexLinearSolver(; arrayType = Array, elType = Float32) =# end +function testVerboseSolvers(; arrayType = Array, elType = Float32) + A = rand(elType, 3, 2) + x = rand(elType, 2) + b = A * x + + solvers = [ADMM, FISTA, POGM, OptISTA, SplitBregman] + + for solver in solvers + @test try + S = createLinearSolver(solver, arrayType(A), iterations = 3, verbose = true) + solve!(S, arrayType(b)) + true + catch e + @error e + false + end + end +end + @testset "Test Solvers" begin for arrayType in arrayTypes @@ -240,4 +259,5 @@ end end end end + testVerboseSolvers() end \ No newline at end of file From 12b8ef45b34cb3510a77752a4179a8c81369660c Mon Sep 17 00:00:00 2001 From: nHackel Date: Fri, 23 Aug 2024 11:20:46 +0200 Subject: [PATCH 2/2] Fix smaller formating issues in docu --- docs/src/literate/examples/compressed_sensing.jl | 2 +- docs/src/literate/examples/computed_tomography.jl | 4 ++-- docs/src/literate/examples/getting_started.jl | 6 +++--- docs/src/literate/howto/weighting.jl | 2 +- docs/src/solvers.md | 3 ++- 5 files changed, 9 insertions(+), 8 deletions(-) diff --git a/docs/src/literate/examples/compressed_sensing.jl b/docs/src/literate/examples/compressed_sensing.jl index 4e81d0fe..41ed15f5 100644 --- a/docs/src/literate/examples/compressed_sensing.jl +++ b/docs/src/literate/examples/compressed_sensing.jl @@ -46,7 +46,7 @@ fig # To recover the image from the measurement vector, we solve the TV-regularized least squares problem: # ```math # \begin{equation} -# \underset{\mathbf{x}}{argmin} \frac{1}{2}\vert\vert \mathbf{A}\mathbf{x}-\mathbf{b} \vert\vert_2^2 + \vert\vert\mathbf{x}\vert\vert_{\lambda\text{TV}} . +# \underset{\mathbf{x}}{argmin} \frac{1}{2}\vert\vert \mathbf{A}\mathbf{x}-\mathbf{b} \vert\vert_2^2 + \lambda\vert\vert\mathbf{x}\vert\vert_{\text{TV}} . # \end{equation} # ``` diff --git a/docs/src/literate/examples/computed_tomography.jl b/docs/src/literate/examples/computed_tomography.jl index 2d1fca28..1159dc8e 100644 --- a/docs/src/literate/examples/computed_tomography.jl +++ b/docs/src/literate/examples/computed_tomography.jl @@ -13,7 +13,7 @@ N = 256 image = shepp_logan(N, SheppLoganToft()) size(image) -# This produces a 64x64 image of a Shepp-Logan phantom. +# This produces a 256x256 image of a Shepp-Logan phantom. using RadonKA, LinearOperatorCollection angles = collect(range(0, π, 256)) @@ -43,7 +43,7 @@ fig # To recover the image from the measurement vector, we solve the $l^2_2$-regularized least squares problem # ```math # \begin{equation} -# \underset{\mathbf{x}}{argmin} \frac{1}{2}\vert\vert \mathbf{A}\mathbf{x}-\mathbf{b} \vert\vert_2^2 + \vert\vert\mathbf{x}\vert\vert^2_2 . +# \underset{\mathbf{x}}{argmin} \frac{1}{2}\vert\vert \mathbf{A}\mathbf{x}-\mathbf{b} \vert\vert_2^2 + \lambda\vert\vert\mathbf{x}\vert\vert^2_2 . # \end{equation} # ``` diff --git a/docs/src/literate/examples/getting_started.jl b/docs/src/literate/examples/getting_started.jl index 750d3341..d4050d96 100644 --- a/docs/src/literate/examples/getting_started.jl +++ b/docs/src/literate/examples/getting_started.jl @@ -23,7 +23,7 @@ using RegularizedLeastSquares # \underset{\mathbf{x}}{argmin} \frac{1}{2}\vert\vert \mathbf{A}\mathbf{x}-\mathbf{b} \vert\vert_2^2 + \mathbf{R(x)} . # \end{equation} # ``` -# where $\mathbf{A}$ is a linear operator, $\mathbf{y}$ is the measurement vector, and $\mathbf{R(x)}$ is an (optional) regularization term. +# where $\mathbf{A}$ is a linear operator, $\mathbf{b}$ is the measurement vector, and $\mathbf{R(x)}$ is an (optional) regularization term. # The goal is to retrieve an approximation of the unknown vector $\mathbf{x}$. In this first exampel we will just work with simple random arrays. For more advanced examples, please refer to the examples. A = rand(32, 16) @@ -41,11 +41,11 @@ isapprox(x, x_approx, rtol = 0.001) # The CGNR algorithm can solve optimzation problems of the form: # ```math # \begin{equation} -# \underset{\mathbf{x}}{argmin} \frac{1}{2}\vert\vert \mathbf{A}\mathbf{x}-\mathbf{b} \vert\vert_2^2 + \vert\vert\mathbf{x}\vert\vert^2_2 . +# \underset{\mathbf{x}}{argmin} \frac{1}{2}\vert\vert \mathbf{A}\mathbf{x}-\mathbf{b} \vert\vert_2^2 + \lambda\vert\vert\mathbf{x}\vert\vert^2_2 . # \end{equation} # ``` -# The corresponding solver can be built with the L2 regularization term: +# The corresponding solver can be built with the $l^2_2$-regularization term: solver = createLinearSolver(CGNR, A; reg = L2Regularization(0.0001), iterations=32); x_approx = solve!(solver, b) isapprox(x, x_approx, rtol = 0.001) diff --git a/docs/src/literate/howto/weighting.jl b/docs/src/literate/howto/weighting.jl index 54fc6ebc..d5682d48 100644 --- a/docs/src/literate/howto/weighting.jl +++ b/docs/src/literate/howto/weighting.jl @@ -11,7 +11,7 @@ # In the following, we will solve a weighted least squares problem of the form: # ```math # \begin{equation} -# \underset{\mathbf{x}}{argmin} \frac{1}{2}\vert\vert \mathbf{A}\mathbf{x}-\mathbf{b} \vert\vert_\mathbf{W}^2 + \vert\vert\mathbf{x}\vert\vert^2_2 . +# \underset{\mathbf{x}}{argmin} \frac{1}{2}\vert\vert \mathbf{A}\mathbf{x}-\mathbf{b} \vert\vert_\mathbf{W}^2 + \lambda\vert\vert\mathbf{x}\vert\vert^2_2 . # \end{equation} # ``` using RegularizedLeastSquares, LinearOperatorCollection, LinearAlgebra diff --git a/docs/src/solvers.md b/docs/src/solvers.md index 84a96dec..dc6294ac 100644 --- a/docs/src/solvers.md +++ b/docs/src/solvers.md @@ -82,4 +82,5 @@ SolverVariant(A; kwargs...) = Solver(A, VariantState(kwargs...)) function iterate(solver::Solver, state::VarianteState) # Custom iteration -end \ No newline at end of file +end +``` \ No newline at end of file