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
L⁺ₐ::Matrix{Float64}   # Cached A-weighted generalized inverse of L(standard-form conversion)

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

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

Return the smoothing matrix L for zero, first and second order Tikhonov regularization based on the size of design matrix A. Order can be 0, 1 or 2.

L = Γ(A, 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). Solution is truncated to regularized space, given by the matrix L. If L is p × n and p < n, then only the solution 1:p is valid. The remaining parameters can be estiamted from the least-squares solution if needed.

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$.

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.

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,
    λ₁ = 0.0001,
    λ₂ = 1000.0,
)

Find the optimum regularization parameter λ between [λ₁, λ₂] using the algorithm alg. 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 

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,
    λ₁ = 0.0001,
    λ₂ = 1000.0,
)

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)

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

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

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.

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