Manual

A theoretical background of Tikhonov regularization is provided in the The Inverse Problem section.

Solving a Problem

A problem consists of a design matrix ${\rm {\bf A}}$ and a vector ${\rm b}$ such that

\[{\rm b} = {\bf {\rm {\bf A}{\rm x}}} + \epsilon\]

where $\epsilon$ is noise and the objective is to reconstruct the original parameters ${\rm x}$. The solution to the problem is to minimize

\[{\rm {\rm x_{\lambda}}}=\arg\min\left\{ \left\lVert {\bf {\rm {\bf A}{\rm x}-{\rm b}}}\right\rVert _{2}^{2}+\lambda^{2}\left\lVert {\rm {\bf L}({\rm x}-{\rm x_{0}})}\right\rVert _{2}^{2}\right\} \]

where ${\rm x_{\lambda}}$ is the regularized estimate of ${\rm x}$, $\left\lVert \cdot\right\rVert _{2}$ is the Euclidean norm, ${\rm {\bf L}}$ is the Tikhonov filter matrix, $\lambda$ is the regularization parameter, and ${\rm x_{0}}$ is a vector of an a priori guess of the solution. The initial guess can be taken to be ${\rm x_{0}}=0$ if no a priori information is known.

The basic steps to solve the problem are

Ψ = setupRegularizationProblem(A, 2)   # Setup the problem 
solution = solve(Ψ, b)                 # Compute the solution 
xλ = solution.x                        # Extract the x

The solve function finds Optimal Regularization Parameter, by default using Generalized Cross Validation. It applies the optimal $\lambda$ value to compute ${\rm x}_\lambda$. The solution is of type RegularizatedSolution, which contains the optimal solution (solution.x), the optimal $\lambda$ (solution.λ) and output from the otimization routine (solution.solution).

An convenient way to find x is using Lazy pipes:

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

Specifying the Order

Common choices for the ${\bf{\rm{L}}}$ matrix are finite difference approximations of a derivative. There are termed zeroth, first, and second order inversion matrices.

\[{\rm {\bf L}_{0}=\left(\begin{array}{ccccc} 1 & & & & 0\\ & 1\\ & & \ddots\\ & & & 1\\ 0 & & & & 1 \end{array}\right)}\]

\[{\rm {\bf L}_{1}=\left(\begin{array}{cccccc} 1 & -1 & & & & 0\\ & 1 & -1\\ & & \ddots & \ddots\\ & & & 1 & -1\\ 0 & & & & 1 & -1 \end{array}\right)}\]

\[{\rm {\bf L}_{2}=\left(\begin{array}{ccccccc} -1 & 2 & -1 & & & & 0\\ & -1 & 2 & -1\\ & & \ddots & \ddots & \ddots\\ & & & -1 & 2 & -1\\ 0 & & & & -1 & 2 & -1 \end{array}\right)}\]

You can specify which of these matrices to use in

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

where order = 0, 1, 2 corresponds to ${\bf{\rm{L}}_0}$, ${\bf{\rm{L}}_1}$, and ${\bf{\rm{L}}_2}$. Higher orders can be specified if desired.

Example : Phillips Problem

Example with 100 point discretization and zero initial guess.

using RegularizationTools, MatrixDepot, Lazy, DelimitedFiles
random(n) = @> readdlm("random.txt") vec x -> x[1:n]

r = mdopen("phillips", 100, false)
A, x = r.A, r.x

y = A * x
b = y + 0.1y .* random(100)
xλ = @> setupRegularizationProblem(A, 2) solve(b) getfield(:x)

Example 2: Shaw Problem

Zeroth order example with 500 point discretization and moderate initial guess.

using RegularizationTools, MatrixDepot, Lazy, DelimitedFiles
random(n) = @> readdlm("random.txt") vec x -> x[1:n]

r = mdopen("shaw", 500, false)
A, x = r.A, r.x

y = A * x
b = y + 0.1y .* random(500)
x₀ = 0.6x

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

Using a Custom L Matrix

You can specify custom ${\rm {\bf L}}$ matrices when setting up problems.

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

For example, Huckle and Sedlacek (2012) propose a two-step data based regularization

\[{\rm {\bf L}} = {\rm {\bf L}}_k {\rm {\bf D}}_{\hat{x}}^{-1} \]

where ${\rm {\bf L}}_k$ is one of the finite difference approximations of a derivative, ${\rm {\bf D}}_{\hat{x}}=diag(|\hat{x_{1}}|,\ldots|\hat{x_{n}}|)$, $\hat{x}$ is the reconstruction of $x$ using ${\rm {\bf L}}_k$, and $({\rm {\bf D}}_{\hat{x}})_{ii}=\epsilon\;\forall\;|\hat{x_{i}}|<\epsilon$, with $\epsilon << 1$.

Example 3: Heat Problem

This examples illustrates how to implement the Huckle and Sedlacek (2012) matrix. Note that Γ(n, 2) returns the Tikhonov Matrix of order 2.

using RegularizationTools, MatrixDepot, Lazy, LinearAlgebra, DelimitedFiles
random(n) = @> readdlm("random.txt") vec x -> x[1:n]

r = mdopen("heat", 100, false)
A, x = r.A, r.x

y = A * x
b = y + 0.04y .* random(100)
x₀ = zeros(length(b))

L₂ = Γ(length(b), 2)
xλ1 = @> setupRegularizationProblem(A, L₂) solve(b) getfield(:x)
x̂ = deepcopy(abs.(xλ1))
x̂[abs.(x̂) .< 0.1] .= 0.1
L = L₂*Diagonal(x̂)^(-1)
xλ2 = @> setupRegularizationProblem(A, L) solve(b) getfield(:x)

The solution xλ2 is improved over the regular L₂ solution.

Adding Boundary Constraints

To add boundary constraints (e.g. enforce that all solutions are positive), the following procedure is implemented. Compute the algebraic solution without constraints, truncate the solution at the upper and lower bounds, and use the result as initial condition for solving the minimization problem with a least squares numerical solver LeastSquaresOptim. The regularization parameter λ obtained from the algebraic solution is used for a single pass optimization. See solve for a complete list of methods.

using RegularizationTools, MatrixDepot, Lazy, LinearAlgebra, DelimitedFiles
random(n) = @> readdlm("random.txt") vec x -> x[1:n]

r = mdopen("heat", 100, false)
A, x = r.A, r.x

y = A * x
b = y + 0.04y .* random(100)
x₀ = zeros(length(b))

L₂ = Γ(length(b),2)
xλ1 = @> setupRegularizationProblem(A, L₂) solve(b) getfield(:x)
x̂ = deepcopy(abs.(xλ1))
x̂[abs.(x̂) .< 0.1] .= 0.1
ψ = @>> L₂*Diagonal(x̂)^(-1) setupRegularizationProblem(A)
lower = zeros(100)
upper = ones(100)
xλ2 = @> solve(ψ, b, lower, upper) getfield(:x)