Data Types

RegularizationProblem

This data type contains the cached matrices used in the inversion. The problem is initialized using the constructor setupRegularizationProblem with the design matrix A and the the Tikhonv matrix L as inputs. The hat quantities, e.g. Ā, is the calculated design matrix in standard form. ĀĀ, Āᵀ, F̄ are precomputed to speed up repeating inversions with different data. L⁺ₐ is cached to speed up the repeated conversion of data to_standard_form and to_general_form

Ā::Matrix{Float64}     # Standard form of design matrix
A::Matrix{Float64}     # General form of the design matrix (n×p)
L::Matrix{Float64}     # Smoothing matrix (n×p)
ĀĀ::Matrix{Float64}    # Cached value of Ā'Ā for performance
Āᵀ::Matrix{Float64}    # Cached value of Ā' for performance
F̄::SVD                 # Cached SVD decomposition of Ā 
Iₙ::Matrix{Float64}    # Cached identity matrix n×n
Iₚ::Matrix{Float64}    # Cached identity matrix p×p
K₀T⁻¹H₀ᵀ::Matrix{Float64} # Cached value to compute 2.42 and 2.44 in Hansen ch 2.

RegularizatedSolution

Data tpye to store the optimal solution x of the inversion. λ is the optimal λ used solution is the raw output from the Optim search.

x::AbstractVector
λ::AbstractFloat
solution::Optim.UnivariateOptimizationResults

Domain{T1<:Any,T2<:Number,T3<:Any}

Functor to map from a domain characterized by a list of setpoints [s], each associated with a list of numerical values [x] to a query value q.

s::AbstractVector{T1}
x::AbstractVector{T2}
q::T3

Γ(m::Int, order::Int)

Return the smoothing matrix L for Tikhonov regularization of a system of size m. Order can be 0, 1, ..., n. Code to generate matrix is based on the suggestion posted in Issue #7, which was inspired by Eilers, P. H. C. (2003). Analytical Chemistry, 75(14), 3631–3636. (Supporting Information).

L = Γ(m, 1)

setupRegularizationProblem(A::AbstractMatrix, order::Int)

Precompute matrices to initialize Reguluarization Problem based on design matrix A and zeroth, first, or second order difference operator. See Hanson (1998) and source code for details.

Example Usage

Ψ = setupRegularizationProblem(A, 0) # zeroth order problem
Ψ = setupRegularizationProblem(A, 2) # second order problem

setupRegularizationProblem(A::AbstractMatrix, L::AbstractMatrix)

Precompute matrices to initialize Reguluarization Problem based on design matrix and Tikhonov smoothing matrix. See Hansen (1998, Eq. 2.35)

Example Usage

Ψ = setupRegularizationProblem(A, L) 

to_standard_form(Ψ::RegularizationProblem, b::AbstractVector)

Converts vector b to standard form using (Hansen, 1998)

Example Usage (Regular Syntax)

b̄ = to_standard_form(Ψ, b)

Example Usage (Lazy Syntax)

b̄ = @>> b to_standard_form(Ψ)

to_standard_form(Ψ::RegularizationProblem, b::AbstractVector, x₀::AbstractVector)

Converts vector b and x₀ to standard form using (Hansen, 1998)

Example Usage (Regular Syntax)

b̄ = to_standard_form(Ψ, b, x₀)

to_general_form(Ψ::RegularizationProblem, b::AbstractVector, x̄::AbstractVector)

Converts solution $\bar {\rm x}$ computed in standard form back to general form ${\rm x}$ using (Hansen, 1998).

\[{\rm x}={\rm {\bf L^{+}_A}\bar{x} + x\_0}\]

where the matrices and vectors are defined in RegularizationProblem

Example Usage (Regular Syntax)

x = to_general_form(Ψ, b, x̄) 

Example Usage (Lazy Syntax)

x = @>> x̄ to_general_form(Ψ, b) 

solve(Ψ::RegularizationProblem, b̄::AbstractVector, λ::AbstractFloat)

Compute the Tikhonov solution for problem Ψ in standard form for regularization parameter λ and using zero as initial guess. Returns a vector $\rm {\bar x}_\lambda$.

\[{\rm x_{\lambda}}=\left({\rm {\bf \bar A}^{T}}{\rm {\bf \bar A}}+\lambda^{2}{\rm {\bf I}}\right)^{-1} {\rm {\bf {\bar A}}^{T}}{\rm {\bar b}} \]

Example Usage (Standard Syntax)

# A is a Matrix and b is a response vector. 
Ψ = setupRegularizationProblem(A, 1)     # Setup problem
b̄ = to_standard_form(Ψ, b)               # Convert to standard form
x̄ = solve(A, b̄, 0.5)                     # Solve the equation
x = to_general_form(Ψ, b, x̄)             # Convert back to general form

Example Usage (Lazy Syntax)

# A is a Matrix and b is a response vector. 
Ψ = setupRegularizationProblem(A, 1)     # Setup problem
b̄ = @>> b to_standard_form(Ψ)            # Convert to standard form
x̄ = solve(A, b̄, 0.5)                     # Solve the equation
x = @>> x̄ to_general_form(Ψ, b)          # Convert back to general form

solve(Ψ::RegularizationProblem, b̄::AbstractVector, x̄₀::AbstractVector, λ::AbstractFloat)

Compute the Tikhonov solution for problem Ψ in standard form for regularization parameter λ and using x̄₀ as initial guess.

\[{\rm x_{\lambda}}=\left({\rm {\bf \bar A}^{T}}{\rm {\bf \bar A}}+\lambda^{2}{\rm {\bf I}}\right)^{-1} \left({\rm {\bf {\bar A}}^{T}}{\rm {\bar b}} + \lambda^2 {\rm {\bar x}}_0 \right)\]

Example Usage (Standard Syntax)

# A is a Matrix and b is a response vector. 
Ψ = setupRegularizationProblem(A, 2)     # Setup problem
b̄, x̄₀ = to_standard_form(Ψ, b, x₀)       # Convert to standard form
x̄ = solve(A, b̄, x̄₀, 0.5)                 # Solve the equation
x = to_general_form(Ψ, b, x̄)             # Convert back to general form

function solve(
    Ψ::RegularizationProblem,
    b::AbstractVector;
    alg = :gcv_svd,
    method = Brent(),
    λ₁ = 0.0001,
    λ₂ = 1000.0,
    λ₀ = 1.0,
    kwargs...
)

Find the optimum regularization parameter λ using the algorithm alg and optimization method. Choices for algorithms are

    :gcv_tr - generalized cross validation using the trace formulation (slow)
    :gcv_svd - generalized cross validation using the SVD decomposition (fast)
    :L_curve - L-curve algorithm 

Optimization methods are those available in the Optim package. Methods that optimize on bounds (Brent() and GoldenSection()) will use [λ₁, λ₂] for bounds. Other methods will use λ₀ for the initial guess. The methods Brent() (default) and NelderMead() are recommended. Any additional keyword arguments will be passed on to the optimization method.

The solve function takes the original data, converts it to standard form, performs the search within the specified bounds and returns a RegularizatedSolution

Example Usage (Standard Syntax)

# A is a Matrix and b is a response vector. 
Ψ = setupRegularizationProblem(A, 2)     # Setup problem
sol = solve(Ψ, b)                        # Solve it

Example Usage (Lazy Syntax)

# A is a Matrix and b is a response vector. 
sol = @> setupRegularizationProblem(A, 1) solve(b)

function solve(
    Ψ::RegularizationProblem,
    b::AbstractVector,
    x₀::AbstractVector;
    alg = :gcv_svd,
    method = Brent(),
    λ₁ = 0.0001,
    λ₂ = 1000.0,
    λ₀ = 1.0,
    kwargs...
)

Same as above, but includes an initial guess x₀. Example Usage (Lazy Syntax)

# A is a Matrix and b is a response vector. 
sol = @> setupRegularizationProblem(A, 1) solve(b, x₀, alg = :L_curve, λ₂ = 10.0)

function solve(
    Ψ::RegularizationProblem, 
    b::AbstractVector,
    lower::AbstractVector, 
    upper::AbstractVector;
    kwargs...
)

Constraint minimization of RegularizationProblem Ψ, with observations b and upper and lower bounds for each xᵢ.

The function computes the algebraic solution using solve(Ψ, b; kwargs...), truncates the solution at the upper and lower bounds and uses this solution as initial condition for the minimization problem using a Least Squares numerical solver. The returned solution is using the regularization parameter λ obtained from the algebraic solution.

function solve(
    Ψ::RegularizationProblem, 
    b::AbstractVector,
    x₀::AbstractVector,
    lower::AbstractVector, 
    upper::AbstractVector;
    kwargs...
)

Constraint minimization of RegularizationProblem Ψ, with observations b, intial guess x₀ and upper and lower bounds for each xᵢ.

The function computes the algebraic solution using solve(Ψ, b; kwargs...), truncates the solution at the upper and lower bounds and uses this solution as initial condition for the minimization problem using a Least Squares numerical solver. The returned solution is using the regularization parameter λ obtained from the algebraic solution.

gcv_tr(Ψ::RegularizationProblem, b̄::AbstractVector, λ::AbstractFloat)

Compute the Generalized Cross Validation using the trace term. Requires that the vector b̄ is in standard form.

\[V(\lambda)=\frac{n\left\lVert ({\bf {\rm {\bf {\bf I}-}{\bf \bar {A}_{\lambda}}) {\bar{\rm b}}}}\right\rVert _{2}^{2}}{tr({\rm {\bf I}-{\rm {\bar {\bf A}_{\lambda}})}^{2}}}\]

Example Usage

using Underscores

Ψ = setupRegularizationProblem(A, 1)           # Setup problem
b̄ = to_standard_form(Ψ, b)                     # Convert to standard form
Vλ = gcv_tr(Ψ, b̄, 0.1)                         # V(λ) single λ value
Vλ = @_ map(gcv_tr(Ψ, b̄, _), [0.1, 1.0, 10.0]) # V(λ) for array of λ

gcv_tr(
    Ψ::RegularizationProblem,
    b̄::AbstractVector,
    x̄₀::AbstractVector,
    λ::AbstractFloat,
)

Compute the Generalized Cross Validation using the trace term and intial guess. Requires that the vectors b̄ and x̄₀ are in standard form.

\[V(\lambda)=\frac{n\left\lVert {\bf {\rm {\bf \bar{A}}{\rm \bar{x}{}_{\lambda}}- {\rm \bar{b}}}}\right\rVert _{2}^{2}}{tr({\rm {\bf I}-{\rm {\bar {\bf A}_{\lambda}})}^{2}}}\]

Example Usage

using Underscores

Ψ = setupRegularizationProblem(A, 1)               # Setup problem
b̄, x̄₀ = to_standard_form(Ψ, b, x₀)                 # Convert to standard form
Vλ = gcv_tr(Ψ, b̄, x̄₀, 0.1)                         # V(λ) single λ value
Vλ = @_ map(gcv_tr(Ψ, b̄, x̄₀, _), [0.1, 1.0, 10.0]) # V(λ) for array of λ

gcv_svd(Ψ::RegularizationProblem, b̄::AbstractVector, λ::AbstractFloat)

Compute the Generalized Cross Validation using the trace term using the SVD algorithm. Requires that the vector b̄ is in standard form.

Example Usage

using Underscores

Ψ = setupRegularizationProblem(A, 1)            # Setup problem
b̄, x̄₀ = to_standard_form(Ψ, b, x₀)              # Convert to standard form
Vλ = gcv_svd(Ψ, b̄, x̄₀, 0.1)                     # V(λ) single λ value
Vλ = @_ map(gcv_svd(Ψ, b̄, _), [0.1, 1.0, 10.0]) # V(λ) for array of λ

gcv_svd(
    Ψ::RegularizationProblem,
    b̄::AbstractVector,
    x̄₀::AbstractVector,
    λ::AbstractFloat,
)

Compute the Generalized Cross Validation using the SVD algorithm and intial guess. Requires that the vectors b̄ and x̄₀ are in standard form.

Example Usage

using Underscores

Ψ = setupRegularizationProblem(A, 1)               # Setup problem
b̄, x̄₀ = to_standard_form(Ψ, b, x₀)                 # Convert to standard form
Vλ = gcv_tr(Ψ, b̄, x̄₀, 0.1)                         # V(λ) single λ value
Vλ = @_ map(gcv_tr(Ψ, b̄, x̄₀, _), [0.1, 1.0, 10.0]) # V(λ) for array of λ

Lcurve_functions(Ψ::RegularizationProblem, b̄::AbstractVector)

Compute the L-curve functions to evaluate the norms L1, L2, and the curvature κ. Requires that the vectors b̄ is in standard form.

Example Usage

Ψ = setupRegularizationProblem(A, 1)
b̄ = to_standard_form(Ψ, b)
L1norm, L2norm, κ = Lcurve_functions(Ψ, b̄)

L1norm.([0.1, 1.0, 10.0])    # L1 norm for λ's
L2norm.([0.1, 1.0, 10.0])    # L2 norm for λ's
κ.([0.1, 1.0, 10.0])         # L-curve curvature for λ's

Lcurve_functions(Ψ::RegularizationProblem, b̄::AbstractVector, x̄₀::AbstractVector)

Compute the L-curve functions to evaluate the norms L1, L2, and the curvature κ. Requires that the vectors b̄ and x̄₀ are in standard form.

Example Usage

Ψ = setupRegularizationProblem(A, 1)
b̄, x̄₀ = to_standard_form(Ψ, b, x₀)                 
L1norm, L2norm, κ = Lcurve_functions(Ψ, b̄, x̄₀)

L1norm.([0.1, 1.0, 10.0])    # L1 norm for λ's
L2norm.([0.1, 1.0, 10.0])    # L2 norm for λ's
κ.([0.1, 1.0, 10.0])         # L-curve curvature for λ's

designmatrix(s::Any, q::Any, f::Function)::AbstractMatrix

• s is an array of setpoints
• q is an array of query points
• f is a function that maps the functor Domain to a solution y

The function to creates an array of domain nodes using the standard basis eᵢ of the vector space. The setpoint or query point can be any abstract notion of the input. For examples, a numerical value corresponding to a setting, a string label (["bin A", "bin B", ...), or a list of pixel coordinate [(1,1), (1,2), ...]. The function f must accept a single argument of type Domain and provide a numerical mapping between input $x$ and output $y$ at query point $q$. The function designmatrix then returns a design matrix $\mathbf{A}$ such that

$y = \mathbf{A}x$

where x is an array of numerical input values. In the case that [q] = [s], the shortcut designmatrix(s, f) can be used.

forwardmodel(
    s::AbstractVector{T1},
    x::AbstractVector{T2},
    f::Function,
)::AbstractArray where {T1<:Any,T2<:Number}

Forward model that maps numerical values [x] corresponding to setpoints [s] to output [y], given the function f. The function f must accept a single argument of type Domain and provide a numerical mapping between input $x$ and output $y$ at query point $q$.

Note that the designmatrix and forward model are related

A = designmatrix(s, q, f)
y1 = A*x

y2 = forwardmodel(s, x, q, f)
y1 == y2

function invert(A::Matrix, b::Vector, method::InverseMethod; kwargs...)

High-level API function to perform Tikhonov inversion. The function used the algebraic data type InverseMethod

    @data InverseMethod begin 
        Lₖ(Int)                            # Pure Tikhonov
        Lₖx₀(Int,Vector)                   # with initial guess
        LₖB(Int,Vector,Vector)             # with bounds
        Lₖx₀B(Int,Vector,Vector,Vector)    # with initial guess + bounds
        LₖDₓ(Int,Float64)                  # with filter 
        Lₖx₀Dₓ(Int,Vector,Float64)         # with initial guess + filter 
        LₖDₓB(Int,Float64,Vector,Vector)   # with filter + bound
        Lₖx₀DₓB(Int,Vector,Float64,Vector,Vector) # with initial guess + filter + bound
    end

to define the problem and then dispatches to the correct method. The kwargs... are passed to the solve function For example the standard way to perform second order regularization is

xλ = @> setupRegularizationProblem(A, 2) solve(b) getfield(:x)

this can alternativel written as

invert(A, b, Lₖ(2))

where Lₖ(2) denotes the second order method. The method nomenclature is Lₖ for regularization order, B for bounded search, x₀ for initial condition, and Dₓ for the Huckle and Sedlacek (2012) two-step data based regularization. The method data type takes hyper parameters to initialize the search. Examples of method initializations are

# Hyper parameters 
k: order, lb: low bound, ub: upper bound, ε: noise level, x₀: initial guess
k, lb, ub, ε, x₀ = 2, zeros(8), zeros(8) .+ 50.0, 0.02, 0.5*N

xλ = invert(A, b, Lₖ(k); alg = :gcv_tr, λ₁ = 0.1)
xλ = invert(A, b, Lₖ(k); alg = :gcv_svd, λ₁ = 0.1)
xλ = invert(A, b, LₖB(k, lb, ub); alg = :L_curve, λ₁ = 0.1)
xλ = invert(A, b, Lₖx₀(k, x₀); alg = :L_curve, λ₁ = 0.1)
xλ = invert(A, b, Lₖx₀(k, x₀); alg = :gcv_tr)
xλ = invert(A, b, Lₖx₀(k, x₀); alg = :gcv_svd)
xλ = invert(A, b, Lₖx₀B(k, x₀, lb, ub); alg = :gcv_svd)
xλ = invert(A, b, LₖDₓ(k, ε); alg = :gcv_svd)
xλ = invert(A, b, LₖDₓB(k, ε, lb, ub); alg = :gcv_svd)
xλ = invert(A, b, Lₖx₀Dₓ(k, x₀, ε); alg = :gcv_svd)
xλ = invert(A, b, Lₖx₀DₓB(k, x₀, ε, lb, ub); alg = :gcv_svd)