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).
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.
Example with 100 point discretization and zero initial guess.
using RegularizationTools, MatrixDepot, Lazy
r = mdopen("phillips", 100, false)
A, x = r.A, r.x
y = A * x
b = y + 0.1y .* randn(100)
xλ = @> setupRegularizationProblem(A, 2) solve(b) getfield(:x)
Note
The random perturbation b = y + 0.1y .* randn(100) in each of the examples uses a fixed random seed to ensure reproducibility. The random seed and plot commands are hidden for clarity.
Note
The example system is a test problem for regularization methods is taken from MatrixDepot.jl and is the same system used in Hansen (2000).
Zeroth order example with 500 point discretization and moderate initial guess.
using RegularizationTools, MatrixDepot, Lazy, Random
r = mdopen("shaw", 500, false)
A, x = r.A, r.x
Random.seed!(850) #hide
y = A * x
b = y + 0.1y .* randn(500)
x₀ = 0.6x
xλ = @> setupRegularizationProblem(A, 0) solve(b, x₀) getfield(:x)
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$.
The solve function searches for the optimum regularization parameter $\lambda$ between $[\lambda_1, \lambda_2]$. The default search range is [0.001, 1000.0] and the interval range can be modified through keyword parameters. The optimality criterion is either the minimum of the Generalized Cross Validation function, or the the maximum curvature of the L-curve (see L-Curve Algorithm). The algorithm can be specified through the alg keyword. Valid algorithms are :L_curve, :gcv_svd, and :gcv_tr (see Solve).
using RegularizationTools, MatrixDepot, Lazy, Random
r = mdopen("baart", 100, false)
A, x = r.A, r.x
y = A * x
b = y + 0.1y .* randn(100)
xλ1 = @> setupRegularizationProblem(A,1) solve(b, alg=:L_curve, λ₁=100.0, λ₂=1e6) getfield(:x)
xλ2 = @> setupRegularizationProblem(A,1) solve(b, alg=:gcv_svd) getfield(:x)
Note that the output from the L-curve and GCV algorithm are nearly identical.
Note
The L-curve algorithm is more sensitive to the bounds and slower than the gcv_svd algorithm. There may, however, be cases where the L-curve approach is preferable.
The solution is obtained by first transforming the problem to standard form (see Transformation to Standard Form). The following example can be used to extract the GCV function.
using RegularizationTools, MatrixDepot, Lazy, Random, Underscores, Printf
r = mdopen("shaw", 100, false)
A, x = r.A, r.x
y = A * x
b = y + 0.1y .* randn(100)
Ψ = setupRegularizationProblem(A,1)
λopt = @> solve(Ψ, b, alg=:gcv_svd) getfield(:λ)
b̄ = to_standard_form(Ψ, b)
λs = exp10.(range(log10(1e-1), stop = log10(10), length = 100))
Vλ = @_ map(gcv_svd(Ψ, b̄, _), λs)
The calculated λopt from
λopt = @> solve(Ψ, b, alg=:gcv_svd) getfield(:λ)
is 1.1 and corresponds to the minimum of the GCV curve.
Alternatively, the L-curve is retrieved through the L-curve Functions
using RegularizationTools, MatrixDepot, Lazy, Random, Underscores, Printf
r = mdopen("shaw", 100, false)
A, x = r.A, r.x
y = A * x
b = y + 0.1y .* randn(100)
Ψ = setupRegularizationProblem(A,1)
λopt = @> solve(Ψ, b, alg=:L_curve, λ₂ = 10.0) getfield(:λ)
b̄ = to_standard_form(Ψ, b)
λs = exp10.(range(log10(1e-3), stop = log10(100), length = 200))
L1norm, L2norm, κ = Lcurve_functions(Ψ, b̄)
L1, L2 = L1norm.(λs), L2norm.(λs)
Systems up to a 1000 equations are unproblematic. The setup for much larger system slows down due to the $\approx O(n^2)$ (or worse) time complexity of the SVD and generalized SVD factorization of the design matrix. Larger systems require switching to SVD free algorithms, which are currently not supported by this package.
using RegularizationTools, TimerOutputs, MatrixDepot
to = TimerOutput()
function benchmark(n)
r = mdopen("shaw", n, false)
A, x = r.A, r.x
y = A * x
b = y + 0.05y .* randn(n)
Ψ = setupRegularizationProblem(A, 2)
for i = 1:1
@timeit to "Setup (n = $n)" setupRegularizationProblem(A, 2)
@timeit to "Invert (n = $n)" solve(Ψ, b)
end
end
map(benchmark, [10, 100, 1000])
show(to)