From 043c0bfd0002e33326abe8e828f31c5df1198fd3 Mon Sep 17 00:00:00 2001 From: Jakob Asslaender Date: Tue, 2 Jan 2024 18:33:36 -0500 Subject: [PATCH 1/2] implement changes to ADMM and SplitBregman that allows removal of this type of regularization --- docs/src/API/regularization.md | 1 - src/ADMM.jl | 43 ++++++++----------- src/CGNR.jl | 2 +- src/FISTA.jl | 2 +- src/Kaczmarz.jl | 2 +- src/OptISTA.jl | 2 +- src/POGM.jl | 2 +- src/PrimalDualSolver.jl | 10 +++-- .../ConstraintTransformedRegularization.jl | 24 ----------- src/Regularization/Regularization.jl | 4 -- src/RegularizedLeastSquares.jl | 8 ++-- src/SplitBregman.jl | 42 +++++++++--------- src/deprecated.jl | 16 +++++++ 13 files changed, 70 insertions(+), 88 deletions(-) delete mode 100644 src/Regularization/ConstraintTransformedRegularization.jl create mode 100644 src/deprecated.jl diff --git a/docs/src/API/regularization.md b/docs/src/API/regularization.md index c18e147a..3887043b 100644 --- a/docs/src/API/regularization.md +++ b/docs/src/API/regularization.md @@ -39,7 +39,6 @@ RegularizedLeastSquares.FixedParameterRegularization ```@docs RegularizedLeastSquares.MaskedRegularization RegularizedLeastSquares.TransformedRegularization -RegularizedLeastSquares.ConstraintTransformedRegularization RegularizedLeastSquares.PlugAndPlayRegularization ``` diff --git a/src/ADMM.jl b/src/ADMM.jl index 4bffb0bc..b6d19379 100644 --- a/src/ADMM.jl +++ b/src/ADMM.jl @@ -40,8 +40,8 @@ mutable struct ADMM{matT,opT,R,ropT,P,vecT,rvecT,preconT,rT} <: AbstractPrimalDu end """ - ADMM(A; AHA = A'*A, precon = Identity(), reg = L1Regularization(zero(real(eltype(AHA)))), normalizeReg = NoNormalization(), rho = 1e-1, vary_rho = :none, iterations = 10, iterationsCG = 10, absTol = eps(real(eltype(AHA))), relTol = eps(real(eltype(AHA))), tolInner = 1e-5, verbose = false) - ADMM( ; AHA = , precon = Identity(), reg = L1Regularization(zero(real(eltype(AHA)))), normalizeReg = NoNormalization(), rho = 1e-1, vary_rho = :none, iterations = 10, iterationsCG = 10, absTol = eps(real(eltype(AHA))), relTol = eps(real(eltype(AHA))), tolInner = 1e-5, verbose = false) + ADMM(A; AHA = A'*A, precon = Identity(), reg = L1Regularization(zero(real(eltype(AHA)))), regTrafo = opEye(eltype(AHA), size(AHA,1)), normalizeReg = NoNormalization(), rho = 1e-1, vary_rho = :none, iterations = 10, iterationsCG = 10, absTol = eps(real(eltype(AHA))), relTol = eps(real(eltype(AHA))), tolInner = 1e-5, verbose = false) + ADMM( ; AHA = , precon = Identity(), reg = L1Regularization(zero(real(eltype(AHA)))), regTrafo = opEye(eltype(AHA), size(AHA,1)), normalizeReg = NoNormalization(), rho = 1e-1, vary_rho = :none, iterations = 10, iterationsCG = 10, absTol = eps(real(eltype(AHA))), relTol = eps(real(eltype(AHA))), tolInner = 1e-5, verbose = false) Creates an `ADMM` object for the forward operator `A` or normal operator `AHA`. @@ -54,6 +54,7 @@ Creates an `ADMM` object for the forward operator `A` or normal operator `AHA`. * `AHA` - normal operator is optional if `A` is supplied * `precon` - preconditionner for the internal CG algorithm * `reg::AbstractParameterizedRegularization` - regularization term; can also be a vector of regularization terms + * `regTrafo` - transformation to a space in which `reg` is applied; if `reg` is a vector, `regTrafo` has to be a vector of the same length. Use `opEye(eltype(AHA), size(AHA,1))` if no transformation is desired. * `normalizeReg::AbstractRegularizationNormalization` - regularization normalization scheme; options are `NoNormalization()`, `MeasurementBasedNormalization()`, `SystemMatrixBasedNormalization()` * `rho::Real` - penalty of the augmented Lagrangian * `vary_rho::Symbol` - vary rho to balance primal and dual feasibility; options `:none`, `:balance`, `:PnP` @@ -64,6 +65,8 @@ Creates an `ADMM` object for the forward operator `A` or normal operator `AHA`. * `tolInner::Real` - relative tolerance for CG stopping criterion * `verbose::Bool` - print residual in each iteration +ADMM differs from ISTA-type algorithms in the sense that the proximal operation is applied separately from the transformation to the space in which the penalty is applied. This is reflected by the interface which has `reg` and `regTrafo` as separate arguments. E.g., for a TV penalty, you should NOT set `reg=TVRegularization`, but instead use `reg=L1Regularization(λ), regTrafo=RegularizedLeastSquares.GradientOp(Float64; shape=(Nx,Ny,Nz))`. + See also [`createLinearSolver`](@ref), [`solve!`](@ref). """ ADMM(; AHA, kwargs...) = ADMM(nothing; kwargs..., AHA = AHA) @@ -72,6 +75,7 @@ function ADMM(A ; AHA = A'*A , precon = Identity() , reg = L1Regularization(zero(real(eltype(AHA)))) + , regTrafo = opEye(eltype(AHA), size(AHA,1)) , normalizeReg::AbstractRegularizationNormalization = NoNormalization() , rho = 1e-1 , vary_rho::Symbol = :none @@ -86,23 +90,15 @@ function ADMM(A T = eltype(AHA) rT = real(T) - reg = vec(reg) + reg = isa(reg, AbstractVector) ? reg : [reg] + regTrafo = isa(regTrafo, AbstractVector) ? regTrafo : [regTrafo] + @assert length(reg) == length(regTrafo) "reg and regTrafo must have the same length" - regTrafo = [] indices = findsinks(AbstractProjectionRegularization, reg) proj = [reg[i] for i in indices] proj = identity.(proj) deleteat!(reg, indices) - # Retrieve constraint trafos - for r in reg - trafoReg = findfirst(ConstraintTransformedRegularization, r) - if isnothing(trafoReg) - push!(regTrafo, opEye(T,size(AHA,2))) - else - push!(regTrafo, trafoReg.trafo) - end - end - regTrafo = identity.(regTrafo) + deleteat!(regTrafo, indices) if typeof(rho) <: Number rho = [rT.(rho) for _ ∈ eachindex(reg)] @@ -116,10 +112,10 @@ function ADMM(A β_y = similar(x) # fields for primal & dual variables - z = [similar(x, size(regTrafo[i],1)) for i ∈ eachindex(vec(reg))] - zᵒˡᵈ = [similar(z[i]) for i ∈ eachindex(vec(reg))] - u = [similar(z[i]) for i ∈ eachindex(vec(reg))] - uᵒˡᵈ = [similar(u[i]) for i ∈ eachindex(vec(reg))] + z = [similar(x, size(regTrafo[i],1)) for i ∈ eachindex(reg)] + zᵒˡᵈ = [similar(z[i]) for i ∈ eachindex(reg)] + u = [similar(z[i]) for i ∈ eachindex(reg)] + uᵒˡᵈ = [similar(u[i]) for i ∈ eachindex(reg)] # statevariables for CG # we store them here to prevent CG from allocating new fields at each call @@ -135,8 +131,7 @@ function ADMM(A # normalization parameters reg = normalize(ADMM, normalizeReg, reg, A, nothing) - return ADMM(A,reg,regTrafo,proj,AHA,β,β_y,x,xᵒˡᵈ,z,zᵒˡᵈ,u,uᵒˡᵈ,precon,rho,iterations - ,iterationsCG,cgStateVars,rᵏ,sᵏ,ɛᵖʳⁱ,ɛᵈᵘᵃ,rT(0),Δ,rT(absTol),rT(relTol),rT(tolInner),normalizeReg,vary_rho,verbose) + return ADMM(A, reg, regTrafo, proj, AHA, β, β_y, x, xᵒˡᵈ, z, zᵒˡᵈ, u, uᵒˡᵈ, precon, rho, iterations, iterationsCG, cgStateVars, rᵏ, sᵏ, ɛᵖʳⁱ, ɛᵈᵘᵃ, rT(0), Δ, rT(absTol), rT(relTol), rT(tolInner), normalizeReg, vary_rho, verbose) end """ @@ -144,7 +139,7 @@ end (re-) initializes the ADMM iterator """ -function init!(solver::ADMM, b; x0 = 0) +function init!(solver::ADMM, b; x0=0) solver.x .= x0 # right hand side for the x-update @@ -156,7 +151,7 @@ function init!(solver::ADMM, b; x0 = 0) # primal and dual variables for i ∈ eachindex(solver.reg) - solver.z[i] .= solver.regTrafo[i]*solver.x + solver.z[i] .= solver.regTrafo[i] * solver.x solver.u[i] .= 0 end @@ -181,7 +176,7 @@ solverconvergence(solver::ADMM) = (; :primal => solver.rᵏ, :dual => solver.s performs one ADMM iteration. """ function iterate(solver::ADMM, iteration=1) - if done(solver, iteration) return nothing end + done(solver, iteration) && return nothing solver.verbose && println("Outer ADMM Iteration #$iteration") # 1. solve arg min_x 1/2|| Ax-b ||² + ρ/2 Σ_i||Φi*x+ui-zi||² @@ -195,7 +190,7 @@ function iterate(solver::ADMM, iteration=1) end solver.verbose && println("conjugated gradients: ") solver.xᵒˡᵈ .= solver.x - cg!(solver.x, AHA, solver.β, Pl = solver.precon, maxiter = solver.iterationsCG, reltol = solver.tolInner, statevars = solver.cgStateVars, verbose = solver.verbose) + cg!(solver.x, AHA, solver.β, Pl=solver.precon, maxiter=solver.iterationsCG, reltol=solver.tolInner, statevars=solver.cgStateVars, verbose=solver.verbose) for proj in solver.proj prox!(proj, solver.x) diff --git a/src/CGNR.jl b/src/CGNR.jl index 14866b55..69c2ebae 100644 --- a/src/CGNR.jl +++ b/src/CGNR.jl @@ -62,7 +62,7 @@ function CGNR(A ζl = zero(T) #temporary scalar # Prepare regularization terms - reg = vec(reg) + reg = isa(reg, AbstractVector) ? reg : [reg] reg = normalize(CGNR, normalizeReg, reg, A, nothing) idx = findsink(L2Regularization, reg) if isnothing(idx) diff --git a/src/FISTA.jl b/src/FISTA.jl index f2648c8a..cb7586dc 100644 --- a/src/FISTA.jl +++ b/src/FISTA.jl @@ -76,7 +76,7 @@ function FISTA(A end # Prepare regularization terms - reg = vec(reg) + reg = isa(reg, AbstractVector) ? reg : [reg] indices = findsinks(AbstractProjectionRegularization, reg) other = [reg[i] for i in indices] deleteat!(reg, indices) diff --git a/src/Kaczmarz.jl b/src/Kaczmarz.jl index c4265b62..c751e41f 100644 --- a/src/Kaczmarz.jl +++ b/src/Kaczmarz.jl @@ -66,7 +66,7 @@ function Kaczmarz(A end # Prepare regularization terms - reg = vec(reg) + reg = isa(reg, AbstractVector) ? reg : [reg] reg = normalize(Kaczmarz, normalizeReg, reg, A, nothing) idx = findsink(L2Regularization, reg) if isnothing(idx) diff --git a/src/OptISTA.jl b/src/OptISTA.jl index bc916ced..ea0dbd65 100644 --- a/src/OptISTA.jl +++ b/src/OptISTA.jl @@ -88,7 +88,7 @@ function OptISTA(A θn = (1 + sqrt(1 + 8 * θn^2)) / 2 # Prepare regularization terms - reg = vec(reg) + reg = isa(reg, AbstractVector) ? reg : [reg] indices = findsinks(AbstractProjectionRegularization, reg) other = [reg[i] for i in indices] deleteat!(reg, indices) diff --git a/src/POGM.jl b/src/POGM.jl index 86c44d48..5b86054e 100644 --- a/src/POGM.jl +++ b/src/POGM.jl @@ -99,7 +99,7 @@ function POGM(A rho /= abs(power_iterations(AHA)) end - reg = vec(reg) + reg = isa(reg, AbstractVector) ? reg : [reg] indices = findsinks(AbstractProjectionRegularization, reg) other = [reg[i] for i in indices] deleteat!(reg, indices) diff --git a/src/PrimalDualSolver.jl b/src/PrimalDualSolver.jl index 43ebc2fa..8906165d 100644 --- a/src/PrimalDualSolver.jl +++ b/src/PrimalDualSolver.jl @@ -50,9 +50,11 @@ function PrimalDualSolver(A::Matrix{T} ) where T M,N = size(A) - if (reg isa Vector && reg[1] isa L1Regularization) || reg isa L1Regularization + reg = isa(reg, AbstractVector) ? reg : [reg] + + if reg[1] isa L1Regularization gradientOp = opEye(T,N) #UniformScaling(one(T)) - elseif (reg isa Vector && reg[1] isa TVRegularization) || reg isa TVRegularization + elseif reg[1] isa TVRegularization gradientOp = gradientOperator(T,shape) end @@ -63,9 +65,9 @@ function PrimalDualSolver(A::Matrix{T} y2 = zeros(T,size(gradientOp*x,1)) # normalization parameters - reg = normalize(PrimalDualSolver, normalizeReg, vec(reg), A, nothing) + reg = normalize(PrimalDualSolver, normalizeReg, reg, A, nothing) - return PrimalDualSolver(A,vec(reg),gradientOp,u,x,cO,y1,y2,T(σ),T(τ),T(ϵ),T(PrimalDualGap),enforceReal,enforcePositive,iterations,shape, + return PrimalDualSolver(A,reg,gradientOp,u,x,cO,y1,y2,T(σ),T(τ),T(ϵ),T(PrimalDualGap),enforceReal,enforcePositive,iterations,shape, normalizeReg) end diff --git a/src/Regularization/ConstraintTransformedRegularization.jl b/src/Regularization/ConstraintTransformedRegularization.jl deleted file mode 100644 index 37a0dfc3..00000000 --- a/src/Regularization/ConstraintTransformedRegularization.jl +++ /dev/null @@ -1,24 +0,0 @@ -export ConstraintTransformedRegularization, transform - -""" - ConstraintTransformedRegularization - -Nested regularization term that associates the `nested` regularization term with a transform. - -# Arguments -* `reg` - inner regularization term -* `trafo` - transform associated with `reg` -""" -struct ConstraintTransformedRegularization{S, R<:AbstractRegularization, TR} <: AbstractNestedRegularization{S} - reg::R - trafo::TR - ConstraintTransformedRegularization(reg::R, trafo::TR) where {R<:AbstractRegularization, TR} = new{R, R, TR}(reg, trafo) - ConstraintTransformedRegularization(reg::R, trafo::TR) where {S, R<:AbstractNestedRegularization{S}, TR} = new{S,R, TR}(reg, trafo) -end -innerreg(reg::ConstraintTransformedRegularization) = reg.reg -""" - transform(reg::ConstraintTransformedRegularization) - -return the transform associated with `innerreg(reg)`. -""" -transform(reg::ConstraintTransformedRegularization) = reg.trafo \ No newline at end of file diff --git a/src/Regularization/Regularization.jl b/src/Regularization/Regularization.jl index 2e82112b..760114b5 100644 --- a/src/Regularization/Regularization.jl +++ b/src/Regularization/Regularization.jl @@ -65,7 +65,6 @@ include("ScaledRegularization.jl") include("NormalizedRegularization.jl") include("TransformedRegularization.jl") include("MaskedRegularization.jl") -include("ConstraintTransformedRegularization.jl") include("PlugAndPlayRegularization.jl") @@ -88,9 +87,6 @@ end findsinks(::Type{S}, reg::Vector{<:AbstractRegularization}) where S <: AbstractRegularization = findall(x -> sinktype(x) <: S, reg) -Base.vec(reg::AbstractRegularization) = AbstractRegularization[reg] -Base.vec(reg::AbstractVector{AbstractRegularization}) = reg - """ RegularizationList() diff --git a/src/RegularizedLeastSquares.jl b/src/RegularizedLeastSquares.jl index 6f2a328b..3236ac8e 100644 --- a/src/RegularizedLeastSquares.jl +++ b/src/RegularizedLeastSquares.jl @@ -120,7 +120,7 @@ end Pass `cb` as the callback to `solve!` # Examples -```julia +```julia julia> x_approx = solve!(solver, b) do solver, iteration println(iteration) end @@ -181,6 +181,8 @@ include("PrimalDualSolver.jl") include("Callbacks.jl") +include("deprecated.jl") + """ Return a list of all available linear solvers """ @@ -259,6 +261,4 @@ function createLinearSolver(solver::Type{T}; AHA, kargs...) where {T<:AbstractLi return solver(; [key=>kargs[key] for key in filtered]..., AHA = AHA) end -@deprecate createLinearSolver(solver, A, x; kargs...) createLinearSolver(solver, A; kargs...) - -end +end \ No newline at end of file diff --git a/src/SplitBregman.jl b/src/SplitBregman.jl index 1d7a8091..86faba81 100644 --- a/src/SplitBregman.jl +++ b/src/SplitBregman.jl @@ -19,7 +19,7 @@ mutable struct SplitBregman{matT,opT,R,ropT,P,vecT,rvecT,preconT,rT} <: Abstract # other parameters precon::preconT ρ::rvecT - iterations::Int64 + iterationsOuter::Int64 iterationsInner::Int64 iterationsCG::Int64 # state variables for CG @@ -40,8 +40,8 @@ mutable struct SplitBregman{matT,opT,R,ropT,P,vecT,rvecT,preconT,rT} <: Abstract end """ - SplitBregman(A; AHA = A'*A, precon = Identity(), reg = L1Regularization(zero(real(eltype(AHA)))), normalizeReg = NoNormalization(), rho = 1e-1, iterations = 1, iterationsInner = 10, iterationsCG = 10, absTol = eps(real(eltype(AHA))), relTol = eps(real(eltype(AHA))), tolInner = 1e-5, verbose = false) - SplitBregman( ; AHA = , precon = Identity(), reg = L1Regularization(zero(real(eltype(AHA)))), normalizeReg = NoNormalization(), rho = 1e-1, iterations = 1, iterationsInner = 10, iterationsCG = 10, absTol = eps(real(eltype(AHA))), relTol = eps(real(eltype(AHA))), tolInner = 1e-5, verbose = false) + SplitBregman(A; AHA = A'*A, precon = Identity(), reg = L1Regularization(zero(real(eltype(AHA)))), regTrafo = opEye(eltype(AHA), size(AHA,1)), normalizeReg = NoNormalization(), rho = 1e-1, iterationsOuter = 10, iterationsInner = 10, iterationsCG = 10, absTol = eps(real(eltype(AHA))), relTol = eps(real(eltype(AHA))), tolInner = 1e-5, verbose = false) + SplitBregman( ; AHA = , precon = Identity(), reg = L1Regularization(zero(real(eltype(AHA)))), regTrafo = opEye(eltype(AHA), size(AHA,1)), normalizeReg = NoNormalization(), rho = 1e-1, iterationsOuter = 10, iterationsInner = 10, iterationsCG = 10, absTol = eps(real(eltype(AHA))), relTol = eps(real(eltype(AHA))), tolInner = 1e-5, verbose = false) Creates a `SplitBregman` object for the forward operator `A` or normal operator `AHA`. @@ -54,9 +54,10 @@ Creates a `SplitBregman` object for the forward operator `A` or normal operator * `AHA` - normal operator is optional if `A` is supplied * `precon` - preconditionner for the internal CG algorithm * `reg::AbstractParameterizedRegularization` - regularization term; can also be a vector of regularization terms + * `regTrafo` - transformation to a space in which `reg` is applied; if `reg` is a vector, `regTrafo` has to be a vector of the same length. Use `opEye(eltype(AHA), size(AHA,1))` if no transformation is desired. * `normalizeReg::AbstractRegularizationNormalization` - regularization normalization scheme; options are `NoNormalization()`, `MeasurementBasedNormalization()`, `SystemMatrixBasedNormalization()` * `rho::Real` - weights for condition on regularized variables; can also be a vector for multiple regularization terms - * `iterations::Int` - maximum number of outer iterations. Set to 1 for unconstraint split Bregman + * `iterationsOuter::Int` - maximum number of outer iterations. Set to 1 for unconstraint split Bregman (equivalent to ADMM) * `iterationsInner::Int` - maximum number of inner iterations * `iterationsCG::Int` - maximum number of (inner) CG iterations * `absTol::Real` - absolute tolerance for stopping criterion @@ -64,6 +65,10 @@ Creates a `SplitBregman` object for the forward operator `A` or normal operator * `tolInner::Real` - relative tolerance for CG stopping criterion * `verbose::Bool` - print residual in each iteration +This algorithm solves the constraint problem (Eq. (4.7) in [Tom Goldstein and Stanley Osher](https://doi.org/10.1137/080725891)), i.e. ||R(x)||_1 such that ||Ax -b||_2^2 < σ^2. In order to solve the unconstraint problem (Eq. (4.8) in [Tom Goldstein and Stanley Osher](https://doi.org/10.1137/080725891)), i.e. ||Ax -b||_2^2 + λ ||R(x)||_1, you can either set `iterationsOuter=1` or use ADMM instead, which is equivalent (`iterationsOuter=1` in SplitBregman in implied in ADMM and the SplitBregman variable `iterationsInner` is simply called `iterations` in ADMM) + +Like ADMM, SplitBregman differs from ISTA-type algorithms in the sense that the proximal operation is applied separately from the transformation to the space in which the penalty is applied. This is reflected by the interface which has `reg` and `regTrafo` as separate arguments. E.g., for a TV penalty, you should NOT set `reg=TVRegularization`, but instead use `reg=L1Regularization(λ), regTrafo=RegularizedLeastSquares.GradientOp(Float64; shape=(Nx,Ny,Nz))`. + See also [`createLinearSolver`](@ref), [`solve!`](@ref). """ SplitBregman(; AHA, kwargs...) = SplitBregman(nothing; kwargs..., AHA = AHA) @@ -72,9 +77,10 @@ function SplitBregman(A ; AHA = A'*A , precon = Identity() , reg = L1Regularization(zero(real(eltype(AHA)))) + , regTrafo = opEye(eltype(AHA), size(AHA,1)) , normalizeReg::AbstractRegularizationNormalization = NoNormalization() , rho = 1e-1 - , iterations::Int = 1 + , iterationsOuter::Int = 10 , iterationsInner::Int = 10 , iterationsCG::Int = 10 , absTol::Real = eps(real(eltype(AHA))) @@ -86,23 +92,15 @@ function SplitBregman(A T = eltype(AHA) rT = real(T) - reg = vec(reg) + reg = isa(reg, AbstractVector) ? reg : [reg] + regTrafo = isa(regTrafo, AbstractVector) ? regTrafo : [regTrafo] + @assert length(reg) == length(regTrafo) "reg and regTrafo must have the same length" - regTrafo = [] indices = findsinks(AbstractProjectionRegularization, reg) proj = [reg[i] for i in indices] proj = identity.(proj) deleteat!(reg, indices) - # Retrieve constraint trafos - for r in reg - trafoReg = findfirst(ConstraintTransformedRegularization, r) - if isnothing(trafoReg) - push!(regTrafo, opEye(T,size(AHA,2))) - else - push!(regTrafo, trafoReg.trafo) - end - end - regTrafo = identity.(regTrafo) + deleteat!(regTrafo, indices) if typeof(rho) <: Number rho = [rT.(rho) for _ ∈ eachindex(reg)] @@ -116,9 +114,9 @@ function SplitBregman(A β_y = similar(x) # fields for primal & dual variables - z = [similar(x, size(regTrafo[i],1)) for i ∈ eachindex(vec(reg))] - zᵒˡᵈ = [similar(z[i]) for i ∈ eachindex(vec(reg))] - u = [similar(z[i]) for i ∈ eachindex(vec(reg))] + z = [similar(x, size(regTrafo[i],1)) for i ∈ eachindex(reg)] + zᵒˡᵈ = [similar(z[i]) for i ∈ eachindex(reg)] + u = [similar(z[i]) for i ∈ eachindex(reg)] # statevariables for CG # we store them here to prevent CG from allocating new fields at each call @@ -136,7 +134,7 @@ function SplitBregman(A # normalization parameters reg = normalize(SplitBregman, normalizeReg, reg, A, nothing) - return SplitBregman(A,reg,regTrafo,proj,y,AHA,β,β_y,x,z,zᵒˡᵈ,u,precon,rho,iterations,iterationsInner,iterationsCG,cgStateVars,rᵏ,sᵏ,ɛᵖʳⁱ,ɛᵈᵘᵃ,rT(0),rT(absTol),rT(relTol),rT(tolInner),iter_cnt,normalizeReg,verbose) + return SplitBregman(A,reg,regTrafo,proj,y,AHA,β,β_y,x,z,zᵒˡᵈ,u,precon,rho,iterationsOuter,iterationsInner,iterationsCG,cgStateVars,rᵏ,sᵏ,ɛᵖʳⁱ,ɛᵈᵘᵃ,rT(0),rT(absTol),rT(relTol),rT(tolInner),iter_cnt,normalizeReg,verbose) end """ @@ -252,4 +250,4 @@ function converged(solver::SplitBregman) return true end -@inline done(solver::SplitBregman,iteration::Int) = converged(solver) || (iteration == 1 && solver.iter_cnt > solver.iterations) \ No newline at end of file +@inline done(solver::SplitBregman,iteration::Int) = converged(solver) || (iteration == 1 && solver.iter_cnt > solver.iterationsOuter) \ No newline at end of file diff --git a/src/deprecated.jl b/src/deprecated.jl new file mode 100644 index 00000000..a64cc23e --- /dev/null +++ b/src/deprecated.jl @@ -0,0 +1,16 @@ +@deprecate createLinearSolver(solver, A, x; kargs...) createLinearSolver(solver, A; kargs...) + +function Base.vec(reg::AbstractRegularization) + Base.depwarn("vec(reg::AbstractRegularization) will be removed in a future release. Use `reg = isa(reg, AbstractVector) ? reg : [reg]` instead.", reg; force=true) + return AbstractRegularization[reg] +end + +function Base.vec(reg::AbstractVector{AbstractRegularization}) + Base.depwarn("vec(reg::AbstractRegularization) will be removed in a future release. Use reg = `isa(reg, AbstractVector) ? reg : [reg]` instead.", reg; force=true) + return reg +end + +export ConstraintTransformedRegularization +function ConstraintTransformedRegularization(args...) + error("ConstraintTransformedRegularization has been removed. ADMM and SplitBregman now take the regularizer and the transform as separat inputs.") +end \ No newline at end of file From 0ec3f2b251823ed068bac2d455d151c9e54f0e6a Mon Sep 17 00:00:00 2001 From: Jakob Asslaender Date: Tue, 2 Jan 2024 18:45:52 -0500 Subject: [PATCH 2/2] docstring beautifying --- src/SplitBregman.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/SplitBregman.jl b/src/SplitBregman.jl index 86faba81..f665e240 100644 --- a/src/SplitBregman.jl +++ b/src/SplitBregman.jl @@ -65,7 +65,7 @@ Creates a `SplitBregman` object for the forward operator `A` or normal operator * `tolInner::Real` - relative tolerance for CG stopping criterion * `verbose::Bool` - print residual in each iteration -This algorithm solves the constraint problem (Eq. (4.7) in [Tom Goldstein and Stanley Osher](https://doi.org/10.1137/080725891)), i.e. ||R(x)||_1 such that ||Ax -b||_2^2 < σ^2. In order to solve the unconstraint problem (Eq. (4.8) in [Tom Goldstein and Stanley Osher](https://doi.org/10.1137/080725891)), i.e. ||Ax -b||_2^2 + λ ||R(x)||_1, you can either set `iterationsOuter=1` or use ADMM instead, which is equivalent (`iterationsOuter=1` in SplitBregman in implied in ADMM and the SplitBregman variable `iterationsInner` is simply called `iterations` in ADMM) +This algorithm solves the constraint problem (Eq. (4.7) in [Tom Goldstein and Stanley Osher](https://doi.org/10.1137/080725891)), i.e. `||R(x)||₁` such that `||Ax -b||₂² < σ²`. In order to solve the unconstraint problem (Eq. (4.8) in [Tom Goldstein and Stanley Osher](https://doi.org/10.1137/080725891)), i.e. `||Ax -b||₂² + λ ||R(x)||₁`, you can either set `iterationsOuter=1` or use ADMM instead, which is equivalent (`iterationsOuter=1` in SplitBregman in implied in ADMM and the SplitBregman variable `iterationsInner` is simply called `iterations` in ADMM) Like ADMM, SplitBregman differs from ISTA-type algorithms in the sense that the proximal operation is applied separately from the transformation to the space in which the penalty is applied. This is reflected by the interface which has `reg` and `regTrafo` as separate arguments. E.g., for a TV penalty, you should NOT set `reg=TVRegularization`, but instead use `reg=L1Regularization(λ), regTrafo=RegularizedLeastSquares.GradientOp(Float64; shape=(Nx,Ny,Nz))`.