diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index cc7aa29e..5dbf31dd 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-06-28T18:03:13","documenter_version":"1.5.0"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.4","generation_timestamp":"2024-07-02T15:08:31","documenter_version":"1.5.0"}} \ No newline at end of file diff --git a/dev/API/regularization/index.html b/dev/API/regularization/index.html index 84d4db61..4d3ffad0 100644 --- a/dev/API/regularization/index.html +++ b/dev/API/regularization/index.html @@ -1,5 +1,5 @@ -Regularization Terms · RegularizedLeastSquares.jl

API for Regularizers

This page contains documentation of the public API of the RegularizedLeastSquares. In the Julia REPL one can access this documentation by entering the help mode with ?

RegularizedLeastSquares.L21Regularization — Type
L21Regularization

Regularization term implementing the proximal map for group-soft-thresholding.

Arguments

  • λ - regularization paramter

Keywords

  • slices=1 - number of elements per group
source
RegularizedLeastSquares.LLRRegularization — Type
LLRRegularization

Regularization term implementing the proximal map for locally low rank (LLR) regularization using singular-value-thresholding.

Arguments

  • λ - regularization paramter

Keywords

  • shape::Tuple{Int} - dimensions of the image
  • blockSize::Tuple{Int}=(2,2) - size of patches to perform singular value thresholding on
  • randshift::Bool=true - randomly shifts the patches to ensure translation invariance
  • fullyOverlapping::Bool=false - choose between fully overlapping block or non-overlapping blocks
source
RegularizedLeastSquares.NuclearRegularization — Type
NuclearRegularization

Regularization term implementing the proximal map for singular value soft-thresholding.

Arguments:

  • λ - regularization paramter

Keywords

  • svtShape::NTuple - size of the underlying matrix
source
RegularizedLeastSquares.TVRegularization — Type
TVRegularization

Regularization term implementing the proximal map for TV regularization. Calculated with the Condat algorithm if the TV is calculated only along one real-valued dimension and with the Fast Gradient Projection algorithm otherwise.

Reference for the Condat algorithm: https://lcondat.github.io/publis/Condat-fast_TV-SPL-2013.pdf

Reference for the FGP algorithm: A. Beck and T. Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems", IEEE Trans. Image Process. 18(11), 2009

Arguments

  • λ::T - regularization parameter

Keywords

  • shape::NTuple - size of the underlying image
  • dims - Dimension to perform the TV along. If Integer, the Condat algorithm is called, and the FDG algorithm otherwise.
  • iterationsTV=20 - number of FGP iterations
source

Projection Regularization

Nested Regularization

RegularizedLeastSquares.innerreg — Method
innerreg(reg::AbstractNestedRegularization)

return the inner regularization term of reg. Nested regularization terms also implement the iteration interface.

source

Scaled Regularization

Misc. Nested Regularization

RegularizedLeastSquares.MaskedRegularization — Type
MaskedRegularization

Nested regularization term that only applies prox! and norm to elements of x for which the mask is true.

Examples

julia> positive = PositiveRegularization();
+Regularization Terms · RegularizedLeastSquares.jl        
         
 
API for Regularizers
This page contains documentation of the public API of the RegularizedLeastSquares. In the Julia REPL one can access this documentation by entering the help mode with ?
RegularizedLeastSquares.L21Regularization — Type
L21Regularization
Regularization term implementing the proximal map for group-soft-thresholding.
Arguments
  • λ                  - regularization paramter
Keywords
  • slices=1           - number of elements per group
source
RegularizedLeastSquares.LLRRegularization — Type
LLRRegularization
Regularization term implementing the proximal map for locally low rank (LLR) regularization using singular-value-thresholding.
Arguments
  • λ                  - regularization paramter
Keywords
  • shape::Tuple{Int}            - dimensions of the image
  • blockSize::Tuple{Int}=(2,2)  - size of patches to perform singular value thresholding on
  • randshift::Bool=true         - randomly shifts the patches to ensure translation invariance
  • fullyOverlapping::Bool=false - choose between fully overlapping block or non-overlapping blocks
source
RegularizedLeastSquares.NuclearRegularization — Type
NuclearRegularization
Regularization term implementing the proximal map for singular value soft-thresholding.
Arguments:
  • λ           - regularization paramter
Keywords
  • svtShape::NTuple  - size of the underlying matrix
source
RegularizedLeastSquares.TVRegularization — Type
TVRegularization
Regularization term implementing the proximal map for TV regularization. Calculated with the Condat algorithm if the TV is calculated only along one real-valued dimension and with the Fast Gradient Projection algorithm otherwise.
Reference for the Condat algorithm: https://lcondat.github.io/publis/Condat-fast_TV-SPL-2013.pdf
Reference for the FGP algorithm: A. Beck and T. Teboulle, "Fast Gradient-Based Algorithms for Constrained Total Variation Image Denoising and Deblurring Problems", IEEE Trans. Image Process. 18(11), 2009
Arguments
  • λ::T                    - regularization parameter
Keywords
  • shape::NTuple           - size of the underlying image
  • dims                    - Dimension to perform the TV along. If Integer, the Condat algorithm is called, and the FDG algorithm otherwise.
  • iterationsTV=20         - number of FGP iterations
source
Projection Regularization
Nested Regularization
RegularizedLeastSquares.innerreg — Method
innerreg(reg::AbstractNestedRegularization)
return the inner regularization term of reg. Nested regularization terms also implement the iteration interface.
source
Scaled Regularization
Misc. Nested Regularization
RegularizedLeastSquares.MaskedRegularization — Type
MaskedRegularization
Nested regularization term that only applies prox! and norm to elements of x for which the mask is true.
Examples
julia> positive = PositiveRegularization();
 
 julia> masked = MaskedRegularization(reg, [true, false, true, false]);
 
@@ -8,11 +8,11 @@
   0.0
  -1.0
   0.0
- -1.0
source
RegularizedLeastSquares.TransformedRegularization — Type
TransformedRegularization(reg, trafo)
Nested regularization term that applies prox! or norm on z = trafo * x and returns (inplace) x = adjoint(trafo) * z.
Example
julia> core = L1Regularization(0.8)
 L1Regularization{Float64}(0.8)
 
 julia> wop = WaveletOp(Float32, shape = (32,32));
 
 julia> reg = TransformedRegularization(core, wop);
 
-julia> prox!(reg, randn(32*32)); # Apply soft-thresholding in Wavelet domain
source
RegularizedLeastSquares.PlugAndPlayRegularization — Type
    PlugAndPlayRegularization
Regularization term implementing a given plug-and-play proximal mapping. The actual regularization term is indirectly defined by the learned proximal mapping and as such there is no norm implemented.
Arguments
  • λ                  - regularization paramter
Keywords
  • model       - model applied to the image
  • shape       - dimensions of the image
  • input_transform - transform of image before model
source
Miscellaneous Functions
RegularizedLeastSquares.prox! — Method
prox!(reg::AbstractParameterizedRegularization, x)
perform the proximal mapping defined by reg on x. Uses the regularization parameter defined for reg.
source
RegularizedLeastSquares.prox! — Method
prox!(regType::Type{<:AbstractParameterizedRegularization}, x, λ; kwargs...)
construct a regularization term of type regType with given λ and kwargs and apply its prox! on x
source
LinearAlgebra.norm — Method
norm(reg::AbstractParameterizedRegularization, x)
returns the value of the reg regularization term on x. Uses the regularization parameter defined for reg.
source
LinearAlgebra.norm — Method
norm(regType::Type{<:AbstractParameterizedRegularization}, x, λ; kwargs...)
construct a regularization term of type regType with given λ and kwargs and apply its norm on x
source

+julia> prox!(reg, randn(32*32)); # Apply soft-thresholding in Wavelet domain
source
RegularizedLeastSquares.PlugAndPlayRegularization — Type
    PlugAndPlayRegularization

Regularization term implementing a given plug-and-play proximal mapping. The actual regularization term is indirectly defined by the learned proximal mapping and as such there is no norm implemented.

Arguments

  • λ - regularization paramter

Keywords

  • model - model applied to the image
  • shape - dimensions of the image
  • input_transform - transform of image before model
source

Miscellaneous Functions

RegularizedLeastSquares.prox! — Method
prox!(reg::AbstractParameterizedRegularization, x)

perform the proximal mapping defined by reg on x. Uses the regularization parameter defined for reg.

source
RegularizedLeastSquares.prox! — Method
prox!(regType::Type{<:AbstractParameterizedRegularization}, x, λ; kwargs...)

construct a regularization term of type regType with given λ and kwargs and apply its prox! on x

source
LinearAlgebra.norm — Method
norm(reg::AbstractParameterizedRegularization, x)

returns the value of the reg regularization term on x. Uses the regularization parameter defined for reg.

source
LinearAlgebra.norm — Method
norm(regType::Type{<:AbstractParameterizedRegularization}, x, λ; kwargs...)

construct a regularization term of type regType with given λ and kwargs and apply its norm on x

source
diff --git a/dev/API/solvers/index.html b/dev/API/solvers/index.html index 6bed59ad..681af621 100644 --- a/dev/API/solvers/index.html +++ b/dev/API/solvers/index.html @@ -40,10 +40,10 @@ end plot_trace (generic function with 1 method) -julia> x_approx = solve!(S, b; callbacks = [conv, plot_trace]);

The keyword callbacks allows you to pass a (vector of) callable objects that takes the arguments solver and iteration and prints, stores, or plots intermediate result.

See also StoreSolutionCallback, StoreConvergenceCallback, CompareSolutionCallback for a number of provided callback options.

source

ADMM

RegularizedLeastSquares.ADMM — Type
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.

Required Arguments

  • A - forward operator

OR

  • AHA - normal operator (as a keyword argument)

Optional Keyword Arguments

  • 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
  • iterations::Int - maximum number of (outer) ADMM iterations
  • iterationsCG::Int - maximum number of (inner) CG iterations
  • absTol::Real - absolute tolerance for stopping criterion
  • relTol::Real - relative tolerance for stopping criterion
  • 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, solve!.

source

CGNR

RegularizedLeastSquares.CGNR — Type
CGNR(A; AHA = A' * A, reg = L2Regularization(zero(real(eltype(AHA)))), normalizeReg = NoNormalization(), iterations = 10, relTol = eps(real(eltype(AHA))))
-CGNR( ; AHA = ,       reg = L2Regularization(zero(real(eltype(AHA)))), normalizeReg = NoNormalization(), iterations = 10, relTol = eps(real(eltype(AHA))))

creates an CGNR object for the forward operator A or normal operator AHA.

Required Arguments

  • A - forward operator

OR

  • AHA - normal operator (as a keyword argument)

Optional Keyword Arguments

  • AHA - normal operator is optional if A is supplied
  • reg::AbstractParameterizedRegularization - regularization term; can also be a vector of regularization terms
  • normalizeReg::AbstractRegularizationNormalization - regularization normalization scheme; options are NoNormalization(), MeasurementBasedNormalization(), SystemMatrixBasedNormalization()
  • iterations::Int - maximum number of iterations
  • relTol::Real - tolerance for stopping criterion

See also createLinearSolver, solve!.

source

Kaczmarz

RegularizedLeastSquares.Kaczmarz — Type
Kaczmarz(A; reg = L2Regularization(0), normalizeReg = NoNormalization(), randomized=false, subMatrixFraction=0.15, shuffleRows=false, seed=1234, iterations=10)

Creates a Kaczmarz object for the forward operator A.

Required Arguments

  • A - forward operator

Optional Keyword Arguments

  • reg::AbstractParameterizedRegularization - regularization term
  • normalizeReg::AbstractRegularizationNormalization - regularization normalization scheme; options are NoNormalization(), MeasurementBasedNormalization(), SystemMatrixBasedNormalization()
  • randomized::Bool - randomize Kacmarz algorithm
  • subMatrixFraction::Real - fraction of rows used in randomized Kaczmarz algorithm
  • shuffleRows::Bool - randomize Kacmarz algorithm
  • seed::Int - seed for randomized algorithm
  • iterations::Int - number of iterations

See also createLinearSolver, solve!.

source

FISTA

RegularizedLeastSquares.FISTA — Type
FISTA(A; AHA=A'*A, reg=L1Regularization(zero(real(eltype(AHA)))), normalizeReg=NoNormalization(), iterations=50, verbose = false, rho = 0.95 / power_iterations(AHA), theta=1, relTol=eps(real(eltype(AHA))), restart = :none)
-FISTA( ; AHA=,     reg=L1Regularization(zero(real(eltype(AHA)))), normalizeReg=NoNormalization(), iterations=50, verbose = false, rho = 0.95 / power_iterations(AHA), theta=1, relTol=eps(real(eltype(AHA))), restart = :none)

creates a FISTA object for the forward operator A or normal operator AHA.

Required Arguments

  • A - forward operator

OR

  • AHA - normal operator (as a keyword argument)

Optional Keyword Arguments

  • 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
  • normalizeReg::AbstractRegularizationNormalization - regularization normalization scheme; options are NoNormalization(), MeasurementBasedNormalization(), SystemMatrixBasedNormalization()
  • rho::Real - step size for gradient step; the default is 0.95 / max_eigenvalue as determined with power iterations.
  • theta::Real - parameter for predictor-corrector step
  • relTol::Real - tolerance for stopping criterion
  • iterations::Int - maximum number of iterations
  • restart::Symbol - :none, :gradient options for restarting
  • verbose::Bool - print residual in each iteration

See also createLinearSolver, solve!.

source

OptISTA

RegularizedLeastSquares.OptISTA — Type
OptISTA(A; AHA=A'*A, reg=L1Regularization(zero(real(eltype(AHA)))), normalizeReg=NoNormalization(), iterations=50, verbose = false, rho=0.95 / power_iterations(AHA), theta=1, relTol=eps(real(eltype(AHA))))
-OptISTA( ; AHA=,     reg=L1Regularization(zero(real(eltype(AHA)))), normalizeReg=NoNormalization(), iterations=50, verbose = false, rho=0.95 / power_iterations(AHA), theta=1, relTol=eps(real(eltype(AHA))))

creates a OptISTA object for the forward operator A or normal operator AHA. OptISTA has a 2x better worst-case bound than FISTA, but actual performance varies by application. It stores 2 extra intermediate variables the size of the image compared to FISTA.

Reference:

  • Uijeong Jang, Shuvomoy Das Gupta, Ernest K. Ryu, "Computer-Assisted Design of Accelerated Composite Optimization Methods: OptISTA," arXiv:2305.15704, 2023, [https://arxiv.org/abs/2305.15704]

Required Arguments

  • A - forward operator

OR

  • AHA - normal operator (as a keyword argument)

Optional Keyword Arguments

  • AHA - normal operator is optional if A is supplied
  • reg::AbstractParameterizedRegularization - regularization term
  • normalizeReg::AbstractRegularizationNormalization - regularization normalization scheme; options are NoNormalization(), MeasurementBasedNormalization(), SystemMatrixBasedNormalization()
  • rho::Real - step size for gradient step; the default is 0.95 / max_eigenvalue as determined with power iterations.
  • theta::Real - parameter for predictor-corrector step
  • relTol::Real - tolerance for stopping criterion
  • iterations::Int - maximum number of iterations
  • verbose::Bool - print residual in each iteration

See also createLinearSolver, solve!.

source

POGM

RegularizedLeastSquares.POGM — Type
POGM(A; AHA = A'*A, reg = L1Regularization(zero(real(eltype(AHA)))), normalizeReg = NoNormalization(), iterations = 50, verbose = false, rho = 0.95 / power_iterations(AHA), theta = 1, sigma_fac = 1, relTol = eps(real(eltype(AHA))), restart = :none)
-POGM( ; AHA = ,     reg = L1Regularization(zero(real(eltype(AHA)))), normalizeReg = NoNormalization(), iterations = 50, verbose = false, rho = 0.95 / power_iterations(AHA), theta = 1, sigma_fac = 1, relTol = eps(real(eltype(AHA))), restart = :none)

Creates a POGM object for the forward operator A or normal operator AHA. POGM has a 2x better worst-case bound than FISTA, but actual performance varies by application. It stores 3 extra intermediate variables the size of the image compared to FISTA. Only gradient restart scheme is implemented for now.

References:

  • A.B. Taylor, J.M. Hendrickx, F. Glineur, "Exact worst-case performance of first-order algorithms for composite convex optimization," Arxiv:1512.07516, 2015, SIAM J. Opt. 2017 [http://doi.org/10.1137/16m108104x]

  • Kim, D., & Fessler, J. A. (2018). Adaptive Restart of the Optimized Gradient Method for Convex Optimization. Journal of Optimization Theory and Applications, 178(1), 240–263. [https://doi.org/10.1007/s10957-018-1287-4]

    Required Arguments

    • A - forward operator

    OR

    • AHA - normal operator (as a keyword argument)

    Optional Keyword Arguments

    • AHA - normal operator is optional if A is supplied
    • reg::AbstractParameterizedRegularization - regularization term
    • normalizeReg::AbstractRegularizationNormalization - regularization normalization scheme; options are NoNormalization(), MeasurementBasedNormalization(), SystemMatrixBasedNormalization()
    • rho::Real - step size for gradient step; the default is 0.95 / max_eigenvalue as determined with power iterations.
    • theta::Real - parameter for predictor-corrector step
    • sigma_fac::Real - parameter for decreasing γ-momentum ∈ [0,1]
    • relTol::Real - tolerance for stopping criterion
    • iterations::Int - maximum number of iterations
    • restart::Symbol - :none, :gradient options for restarting
    • verbose::Bool - print residual in each iteration

See also createLinearSolver, solve!.

source

SplitBregman

RegularizedLeastSquares.SplitBregman — Type
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.

Required Arguments

  • A - forward operator

OR

  • AHA - normal operator (as a keyword argument)

Optional Keyword Arguments

  • 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
  • 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
  • relTol::Real - relative tolerance for stopping criterion
  • 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), i.e. ||R(x)||₁ such that ||Ax -b||₂² < σ². In order to solve the unconstraint problem (Eq. (4.8) in Tom Goldstein and Stanley Osher), 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)).

See also createLinearSolver, solve!.

source

Miscellaneous Functions

RegularizedLeastSquares.StoreSolutionCallback — Type
StoreSolutionCallback(T)

Callback that accumlates the solvers solution per iteration. Results are stored in the solutions field.

source
RegularizedLeastSquares.StoreConvergenceCallback — Type
StoreConvergenceCallback()

Callback that accumlates the solvers convergence metrics per iteration. Results are stored in the convMeas field.

source
RegularizedLeastSquares.CompareSolutionCallback — Type
CompareSolutionCallback(ref, cmp)

Callback that compares the solvers current solution with the given reference via cmp(ref, solution) per iteration. Results are stored in the results field.

source
RegularizedLeastSquares.linearSolverList — Function

Return a list of all available linear solvers

source
RegularizedLeastSquares.createLinearSolver — Function
createLinearSolver(solver::AbstractLinearSolver, A; kargs...)

This method creates a solver. The supported solvers are methods typically used for solving regularized linear systems. All solvers return an approximate solution to Ax = b.

TODO: give a hint what solvers are available

source
RegularizedLeastSquares.applicableSolverList — Function
applicable(args...)

list all solvers that are applicable to the given arguments. Arguments are the same as for isapplicable without the solver type.

See also isapplicable, linearSolverList.

source
RegularizedLeastSquares.isapplicable — Function
isapplicable(solverType::Type{<:AbstractLinearSolver}, A, x, reg)

return true if a solver of type solverType is applicable to system matrix A, data x and regularization terms reg.

source
+julia> x_approx = solve!(S, b; callbacks = [conv, plot_trace]);

The keyword callbacks allows you to pass a (vector of) callable objects that takes the arguments solver and iteration and prints, stores, or plots intermediate result.

See also StoreSolutionCallback, StoreConvergenceCallback, CompareSolutionCallback for a number of provided callback options.

source

ADMM

RegularizedLeastSquares.ADMM — Type
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.

Required Arguments

  • A - forward operator

OR

  • AHA - normal operator (as a keyword argument)

Optional Keyword Arguments

  • 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
  • iterations::Int - maximum number of (outer) ADMM iterations
  • iterationsCG::Int - maximum number of (inner) CG iterations
  • absTol::Real - absolute tolerance for stopping criterion
  • relTol::Real - relative tolerance for stopping criterion
  • 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, solve!.

source

CGNR

RegularizedLeastSquares.CGNR — Type
CGNR(A; AHA = A' * A, reg = L2Regularization(zero(real(eltype(AHA)))), normalizeReg = NoNormalization(), iterations = 10, relTol = eps(real(eltype(AHA))))
+CGNR( ; AHA = ,       reg = L2Regularization(zero(real(eltype(AHA)))), normalizeReg = NoNormalization(), iterations = 10, relTol = eps(real(eltype(AHA))))

creates an CGNR object for the forward operator A or normal operator AHA.

Required Arguments

  • A - forward operator

OR

  • AHA - normal operator (as a keyword argument)

Optional Keyword Arguments

  • AHA - normal operator is optional if A is supplied
  • reg::AbstractParameterizedRegularization - regularization term; can also be a vector of regularization terms
  • normalizeReg::AbstractRegularizationNormalization - regularization normalization scheme; options are NoNormalization(), MeasurementBasedNormalization(), SystemMatrixBasedNormalization()
  • iterations::Int - maximum number of iterations
  • relTol::Real - tolerance for stopping criterion

See also createLinearSolver, solve!.

source

Kaczmarz

RegularizedLeastSquares.Kaczmarz — Type
Kaczmarz(A; reg = L2Regularization(0), normalizeReg = NoNormalization(), randomized=false, subMatrixFraction=0.15, shuffleRows=false, seed=1234, iterations=10)

Creates a Kaczmarz object for the forward operator A.

Required Arguments

  • A - forward operator

Optional Keyword Arguments

  • reg::AbstractParameterizedRegularization - regularization term
  • normalizeReg::AbstractRegularizationNormalization - regularization normalization scheme; options are NoNormalization(), MeasurementBasedNormalization(), SystemMatrixBasedNormalization()
  • randomized::Bool - randomize Kacmarz algorithm
  • subMatrixFraction::Real - fraction of rows used in randomized Kaczmarz algorithm
  • shuffleRows::Bool - randomize Kacmarz algorithm
  • seed::Int - seed for randomized algorithm
  • iterations::Int - number of iterations

See also createLinearSolver, solve!.

source

FISTA

RegularizedLeastSquares.FISTA — Type
FISTA(A; AHA=A'*A, reg=L1Regularization(zero(real(eltype(AHA)))), normalizeReg=NoNormalization(), iterations=50, verbose = false, rho = 0.95 / power_iterations(AHA), theta=1, relTol=eps(real(eltype(AHA))), restart = :none)
+FISTA( ; AHA=,     reg=L1Regularization(zero(real(eltype(AHA)))), normalizeReg=NoNormalization(), iterations=50, verbose = false, rho = 0.95 / power_iterations(AHA), theta=1, relTol=eps(real(eltype(AHA))), restart = :none)

creates a FISTA object for the forward operator A or normal operator AHA.

Required Arguments

  • A - forward operator

OR

  • AHA - normal operator (as a keyword argument)

Optional Keyword Arguments

  • 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
  • normalizeReg::AbstractRegularizationNormalization - regularization normalization scheme; options are NoNormalization(), MeasurementBasedNormalization(), SystemMatrixBasedNormalization()
  • rho::Real - step size for gradient step; the default is 0.95 / max_eigenvalue as determined with power iterations.
  • theta::Real - parameter for predictor-corrector step
  • relTol::Real - tolerance for stopping criterion
  • iterations::Int - maximum number of iterations
  • restart::Symbol - :none, :gradient options for restarting
  • verbose::Bool - print residual in each iteration

See also createLinearSolver, solve!.

source

OptISTA

RegularizedLeastSquares.OptISTA — Type
OptISTA(A; AHA=A'*A, reg=L1Regularization(zero(real(eltype(AHA)))), normalizeReg=NoNormalization(), iterations=50, verbose = false, rho=0.95 / power_iterations(AHA), theta=1, relTol=eps(real(eltype(AHA))))
+OptISTA( ; AHA=,     reg=L1Regularization(zero(real(eltype(AHA)))), normalizeReg=NoNormalization(), iterations=50, verbose = false, rho=0.95 / power_iterations(AHA), theta=1, relTol=eps(real(eltype(AHA))))

creates a OptISTA object for the forward operator A or normal operator AHA. OptISTA has a 2x better worst-case bound than FISTA, but actual performance varies by application. It stores 2 extra intermediate variables the size of the image compared to FISTA.

Reference:

  • Uijeong Jang, Shuvomoy Das Gupta, Ernest K. Ryu, "Computer-Assisted Design of Accelerated Composite Optimization Methods: OptISTA," arXiv:2305.15704, 2023, [https://arxiv.org/abs/2305.15704]

Required Arguments

  • A - forward operator

OR

  • AHA - normal operator (as a keyword argument)

Optional Keyword Arguments

  • AHA - normal operator is optional if A is supplied
  • reg::AbstractParameterizedRegularization - regularization term
  • normalizeReg::AbstractRegularizationNormalization - regularization normalization scheme; options are NoNormalization(), MeasurementBasedNormalization(), SystemMatrixBasedNormalization()
  • rho::Real - step size for gradient step; the default is 0.95 / max_eigenvalue as determined with power iterations.
  • theta::Real - parameter for predictor-corrector step
  • relTol::Real - tolerance for stopping criterion
  • iterations::Int - maximum number of iterations
  • verbose::Bool - print residual in each iteration

See also createLinearSolver, solve!.

source

POGM

RegularizedLeastSquares.POGM — Type
POGM(A; AHA = A'*A, reg = L1Regularization(zero(real(eltype(AHA)))), normalizeReg = NoNormalization(), iterations = 50, verbose = false, rho = 0.95 / power_iterations(AHA), theta = 1, sigma_fac = 1, relTol = eps(real(eltype(AHA))), restart = :none)
+POGM( ; AHA = ,     reg = L1Regularization(zero(real(eltype(AHA)))), normalizeReg = NoNormalization(), iterations = 50, verbose = false, rho = 0.95 / power_iterations(AHA), theta = 1, sigma_fac = 1, relTol = eps(real(eltype(AHA))), restart = :none)

Creates a POGM object for the forward operator A or normal operator AHA. POGM has a 2x better worst-case bound than FISTA, but actual performance varies by application. It stores 3 extra intermediate variables the size of the image compared to FISTA. Only gradient restart scheme is implemented for now.

References:

  • A.B. Taylor, J.M. Hendrickx, F. Glineur, "Exact worst-case performance of first-order algorithms for composite convex optimization," Arxiv:1512.07516, 2015, SIAM J. Opt. 2017 [http://doi.org/10.1137/16m108104x]

  • Kim, D., & Fessler, J. A. (2018). Adaptive Restart of the Optimized Gradient Method for Convex Optimization. Journal of Optimization Theory and Applications, 178(1), 240–263. [https://doi.org/10.1007/s10957-018-1287-4]

    Required Arguments

    • A - forward operator

    OR

    • AHA - normal operator (as a keyword argument)

    Optional Keyword Arguments

    • AHA - normal operator is optional if A is supplied
    • reg::AbstractParameterizedRegularization - regularization term
    • normalizeReg::AbstractRegularizationNormalization - regularization normalization scheme; options are NoNormalization(), MeasurementBasedNormalization(), SystemMatrixBasedNormalization()
    • rho::Real - step size for gradient step; the default is 0.95 / max_eigenvalue as determined with power iterations.
    • theta::Real - parameter for predictor-corrector step
    • sigma_fac::Real - parameter for decreasing γ-momentum ∈ [0,1]
    • relTol::Real - tolerance for stopping criterion
    • iterations::Int - maximum number of iterations
    • restart::Symbol - :none, :gradient options for restarting
    • verbose::Bool - print residual in each iteration

See also createLinearSolver, solve!.

source

SplitBregman

RegularizedLeastSquares.SplitBregman — Type
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.

Required Arguments

  • A - forward operator

OR

  • AHA - normal operator (as a keyword argument)

Optional Keyword Arguments

  • 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
  • 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
  • relTol::Real - relative tolerance for stopping criterion
  • 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), i.e. ||R(x)||₁ such that ||Ax -b||₂² < σ². In order to solve the unconstraint problem (Eq. (4.8) in Tom Goldstein and Stanley Osher), 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)).

See also createLinearSolver, solve!.

source

Miscellaneous Functions

RegularizedLeastSquares.StoreSolutionCallback — Type
StoreSolutionCallback(T)

Callback that accumlates the solvers solution per iteration. Results are stored in the solutions field.

source
RegularizedLeastSquares.StoreConvergenceCallback — Type
StoreConvergenceCallback()

Callback that accumlates the solvers convergence metrics per iteration. Results are stored in the convMeas field.

source
RegularizedLeastSquares.CompareSolutionCallback — Type
CompareSolutionCallback(ref, cmp)

Callback that compares the solvers current solution with the given reference via cmp(ref, solution) per iteration. Results are stored in the results field.

source
RegularizedLeastSquares.linearSolverList — Function

Return a list of all available linear solvers

source
RegularizedLeastSquares.createLinearSolver — Function
createLinearSolver(solver::AbstractLinearSolver, A; kargs...)

This method creates a solver. The supported solvers are methods typically used for solving regularized linear systems. All solvers return an approximate solution to Ax = b.

TODO: give a hint what solvers are available

source
RegularizedLeastSquares.applicableSolverList — Function
applicable(args...)

list all solvers that are applicable to the given arguments. Arguments are the same as for isapplicable without the solver type.

See also isapplicable, linearSolverList.

source
RegularizedLeastSquares.isapplicable — Function
isapplicable(solverType::Type{<:AbstractLinearSolver}, A, x, reg)

return true if a solver of type solverType is applicable to system matrix A, data x and regularization terms reg.

source
diff --git a/dev/gettingStarted/index.html b/dev/gettingStarted/index.html index 8d7f4c1a..60989928 100644 --- a/dev/gettingStarted/index.html +++ b/dev/gettingStarted/index.html @@ -8,4 +8,4 @@ y = A*vec(I)

To recover the image, we solve the TV-regularized least squares problem

\[\begin{equation} \underset{\mathbf{x}}{argmin} \frac{1}{2}\vert\vert \mathbf{A}\mathbf{x}-\mathbf{y} \vert\vert_2^2 + λTV(\mathbf{x}) . \end{equation}\]

For this purpose we build a TV regularizer with regularization parameter $λ=0.01$

reg = TVRegularization(0.01; shape=(N,N))

To solve the CS problem, the Alternating Direction Method of Multipliers can be used. Thus, we build the corresponding solver

solver = createLinearSolver(ADMM, A; reg=reg, ρ=0.1, iterations=20)

and apply it to our measurement

Ireco = solve!(solver,y)
-Ireco = reshape(Ireco,N,N)

The original phantom and the reconstructed image are shown below

Phantom Reconstruction

+Ireco = reshape(Ireco,N,N)

The original phantom and the reconstructed image are shown below

Phantom Reconstruction

diff --git a/dev/index.html b/dev/index.html index d260eb2e..1a24ba6b 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,3 +1,3 @@ Home · RegularizedLeastSquares.jl

RegularizedLeastSquares.jl

Solvers for Linear Inverse Problems using Regularization Techniques

Introduction

RegularizedLeastSquares.jl is a Julia package for solving large scale linear systems using different types of algorithms. Ill-conditioned problems arise in many areas of practical interest. To solve these problems, one often resorts to regularization techniques and non-linear problem formulations. This packages provides implementations for a variety of solvers, which are used in fields such as MPI and MRI.

The implemented methods range from the $l_2$-regularized CGNR method to more general optimizers such as the Alternating Direction of Multipliers Method (ADMM) or the Split-Bregman method.

For convenience, implementations of popular regularizers, such as $l_1$-regularization and TV regularization, are provided. On the other hand, hand-crafted regularizers can be used quite easily. For this purpose, a Regularization object needs to be build. The latter mainly contains the regularization parameter and a function to calculate the proximal map of a given input.

Depending on the problem, it becomes unfeasible to store the full system matrix at hand. For this purpose, RegularizedLeastSquares.jl allows for the use of matrix-free operators. Such operators can be realized using the interface provided by the package LinearOperators.jl. Other interfaces can be used as well, as long as the product *(A,x) and the adjoint adjoint(A) are provided. A number of common matrix-free operators are provided by the package LinearOperatorColection.jl.

Installation

Within Julia, use the package manager:

using Pkg
-Pkg.add("RegularizedLeastSquares")

This adds the latest release of the package is added. To install a different version, please consult the Pkg documentation.

Usage

+Pkg.add("RegularizedLeastSquares")

This adds the latest release of the package is added. To install a different version, please consult the Pkg documentation.

Usage

diff --git a/dev/regularization/index.html b/dev/regularization/index.html index ab4d7186..208cf547 100644 --- a/dev/regularization/index.html +++ b/dev/regularization/index.html @@ -36,4 +36,4 @@ julia> foreach(r -> println(nameof(typeof(r))), reg) TransformedRegularization -L1Regularization +L1Regularization diff --git a/dev/solvers/index.html b/dev/solvers/index.html index 181d1bee..ab6123ad 100644 --- a/dev/solvers/index.html +++ b/dev/solvers/index.html @@ -4,4 +4,4 @@ ... solver = createLinearSolver(ADMM, A; params...)

This notation can be convenient when a large number of parameters are set manually.

It is possible to check if a given solver is applicable to the wanted arguments, as not all solvers are applicable to all system matrix and data (element) types or regularization terms combinations. This is achieved with the isapplicable function:

isapplicable(Kaczmarz, A, x, [L21Regularization(0.4f0)])
-false

For a given set of arguments the list of applicable solvers can be retrieved with applicableSolverList.

+false

For a given set of arguments the list of applicable solvers can be retrieved with applicableSolverList.