Edit on GitHub

RegularizationTools.jl

A Julia package to perform Tikhonov regularization for small to moderate size problems.

RegularizationTools.jl bundles a set routines to compute the regularized Tikhonov inverse using standard linear algebra techniques.

Package Features

  • Computes the Tikhonov inverse solution with optional boundary constraints.
  • Computes optimal regularization parameter using generalized cross validation or the L-curve.
  • Solves problems up to ~1000 equations.
  • Supports higher order regularization schemes out of the box.
  • Supports specifying an a-priori estimate of the solution.
  • Supports user specified smoothing matrices.
  • User friendly interface.
  • Extensive documentation.

About

Tikhonv regularization is also known as Phillips-Twomey-Tikhonov regularization or ridge regression (see Hansen, 2000 for a review). The Web-of-Sciences database lists more than 4500 peer-reviewed publications mentioning "Tikhonov regularization" in the title or abstract, with a current publication rate of ≈350 new papers/year.

Quick Start

The package computes the regularized Tikhonov inverse ${\rm x_{\lambda}}$ by solving the minimization problem

\[{\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\} \]

subject to the optional constraint ${\rm x}_{l}<{\rm x}<{\rm x}_{u}$. Here ${\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 solve function searches for the optimal $\lambda$ and returns the inverse. Optionally, ${\rm x}_{l}$ and ${\rm x}_{u}$ can be used to impose boundary constraints on the solution ${\rm x_{\lambda}}$.

The following script is a minimalist example how to use this package.

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

# This is a test problem for regularization methods
r = mdopen("shaw", 100, false)       # Load the "shaw" problem from MatrixDepot
A, x  = r.A, r.x                     # A is size(100,100), x is length(100)

y = A * x                            # y is the true response
b = y + 0.2y .* random(100)          # response with superimposed noise
x₀ = 0.4x                            # some a-priori estimate x₀

# Solve 2nd order Tikhonov inversion (L = uppertridiag(−1, 2, −1)) with intial guess x₀
xλ = invert(A, b, Lₖx₀(0, x₀))
n -150 -100 -50 0 50 100 150 200 250 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 -100 0 100 200 -100 -95 -90 -85 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 x xλ x0 h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 -2.5 0.0 2.5 5.0 -2.5 -2.4 -2.3 -2.2 -2.1 -2.0 -1.9 -1.8 -1.7 -1.6 -1.5 -1.4 -1.3 -1.2 -1.1 -1.0 -0.9 -0.8 -0.7 -0.6 -0.5 -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 2.0 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4.0 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 5.0 n -150 -100 -50 0 50 100 150 200 250 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 -100 0 100 200 -100 -95 -90 -85 -80 -75 -70 -65 -60 -55 -50 -45 -40 -35 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 180 185 190 195 200 y b h,j,k,l,arrows,drag to pan i,o,+,-,scroll,shift-drag to zoom r,dbl-click to reset c for coordinates ? for help ? -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 -5 0 5 10 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 7.2 7.4 7.6 7.8 8.0 8.2 8.4 8.6 8.8 9.0 9.2 9.4 9.6 9.8 10.0
Info

Per documentation, the way random numbers are generated in Julia is considered an implementation detail. Bug fixes and speed improvements may change the stream of numbers that are generated after a version change. To keep the examples in the documentation of this package consistent between Julia updates, the function random(n) is used to load n pre-generated normally distributed random numbers from a text file.

Installation

The package can be installed from the Julia package prompt with

julia> ]add RegularizationTools

The closing square bracket switches to the package manager interface and the add command installs the package and any missing dependencies. To return to the Julia REPL hit the delete key.

To load the package run

julia> using RegularizationTools

For optimal performance, also install the Intel MKL linear algebra library.

  • MultivariateStats: Implements ridge regression without a priori estimate and does not include tools to find the optimal regularization parameter.
  • RegularizedLeastSquares: Implements optimization techniques for large-scale scale linear systems.

Markus Petters, Department of Marine, Earth, and Atmospheric Sciences, NC State University .

Contribtions

Jonathan Stickel

  • Several bug fixes
  • Improvement of derivative matrix support

Citation and Acknowledgements

If you use this package and publish your work, please cite:

Petters, M. D.: Revisiting matrix-based inversion of scanning mobility particle sizer (SMPS) and humidified tandem differential mobility analyzer (HTDMA) data, Atmos. Meas. Tech., 14, 7909–7928, https://doi.org/10.5194/amt-14-7909-2021, 2021.

This work was supported by the United States Department of Energy, Office of Science, Biological and Environment Research, Grant number DE-SC0021074.