The Inverse Problem

Consider the following linear system of equations

\[{\bf {\rm {\bf A}{\rm x}={\rm y}}}\]

where ${\bf {\rm {\bf A}}}$ is a square design matrix, ${\rm x}$ is a vector of input parameters and ${\rm y}$ is a vector of responses. To estimate unknown inputs from response, the matrix inverse can be used

\[{\rm x={\rm {\bf A}}^{-1}y}\]

However, if a random measurement error $\epsilon$ is superimposed on ${\rm y}$, i.e. $b_{i}=y_{i}+\epsilon_{i}$, the estimate ${\rm \hat{x}}$ from the matrix inverse

\[{\rm \hat{x}={\rm {\bf A}}^{-1}b}\]

becomes dominated by contributions from data error for large systems.

Example

Consider the following example system of 100 equations. The matrix $\rm{\bf{A}}$ is 100x100.

A

100×100 Matrix{Float64}:
 4.71979e-13  4.7207e-11   7.1812e-10   …  0.000123953  3.10037e-5
 4.7207e-11   3.44407e-10  2.18592e-9      0.00027885   0.000123953
 7.1812e-10   2.18592e-9   7.39566e-9      0.000495319  0.000278567
 4.72173e-9   9.94389e-9   2.33509e-8      0.000772739  0.000494292
 1.98422e-8   3.4489e-8    6.52639e-8      0.00111017   0.000770268
 6.32328e-8   9.79521e-8   1.61791e-7   …  0.00150629   0.00110528
 1.67027e-7   2.39805e-7   3.61585e-7      0.00195933   0.00149769
 3.85205e-7   5.24189e-7   7.41078e-7      0.00246701   0.00194536
 8.01584e-7   1.04834e-6   1.41335e-6      0.00302642   0.00244557
 1.53875e-6   1.95195e-6   2.53782e-6      0.00363393   0.00299494
 ⋮                                      ⋱               
 0.00244557   0.00302642   0.00366901      1.04834e-6   8.01584e-7
 0.00194536   0.00246701   0.00305023      5.24189e-7   3.85205e-7
 0.00149769   0.00195933   0.00248243      2.39805e-7   1.67027e-7
 0.00110528   0.00150629   0.00196869      9.79521e-8   6.32328e-8
 0.000770268  0.00111017   0.00151147   …  3.4489e-8    1.98422e-8
 0.000494292  0.000772739  0.00111262      9.94389e-9   4.72173e-9
 0.000278567  0.000495319  0.000773564     2.18592e-9   7.1812e-10
 0.000123953  0.00027885   0.000495319     3.44407e-10  4.7207e-11
 3.10037e-5   0.000123953  0.000278567     4.7207e-11   4.71979e-13

The vector $\rm{x}$ of input variables has 100 elements.

x

100-element Vector{Float64}:
 0.1079137578052813
 0.12297053360334534
 0.13957600811265403
 0.15779962367244207
 0.17769967237809117
 0.19932086589081324
 0.22269187777914254
 0.2478229164774714
 0.2747033932325347
 0.3032997543151104
 ⋮
 0.43775867610162883
 0.35986269701602813
 0.2923532021212666
 0.23471831462134501
 0.18623171118659587
 0.14602510434716662
 0.11315366183414295
 0.08665170515974066
 0.06557729627915371

Computing $\rm{y}$ and $\rm{b}$ using the pseudo-inverse pinv shows that the error in $\rm{b}$ makes the inversion unusable.

random(n) = @> readdlm("random.txt") vec x -> x[1:n] # hide


y = A * x
b = y + 0.1y .* random(100)
x = pinv(A) * y
x̂ = pinv(A) * b

Tikhonov Regularization

Tikhonov regularization is a means to filter this noise 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\} \]

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 matrix ${\rm {\bf A}}$ does not need to be square.

For $\lambda=0$ the Tikhonov problem reverts to the ordinary least squares solution. If ${\rm {\bf A}}$ is square and $\lambda=0$, the least-squares solution is ${\rm \hat{x}={\rm {\bf A}}^{-1}b}$. For large $\lambda$ the solution reverts to the initial guess., i.e. $\lim_{\lambda\rightarrow\infty}{\rm x_{\lambda}}={\rm x_{0}}$. Therefore, the regularization parameter $\lambda$ interpolates between the initial guess and the noisy ordinary least squares solution. The filter matrix ${\rm {\bf L}}$ provides additional smoothness constraints on the solution. The simplest form is to use the identity matrix, ${\rm {\bf L}}={\rm {\bf I}}$.

The formal solution to the Tikhonov problem is given by

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

The equation is readily derived by writing $f=\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}$, take $\frac{df}{d{\rm {\rm x}}}=0$, and solve for ${\rm x}$. Use http://www.matrixcalculus.org/ to validate symbolic matrix derivatives.

Example

Here is a simple regularized inversion for the same system using ${\rm {\bf L}}={\rm {\bf I}}$ and ${\rm x_{0}}=0$. The regularized solution is

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

The regularized inverse can be trivially computed assuming a value for $\lambda = 0.11$.

random(n) = @> readdlm("random.txt") vec x -> x[1:n] #hide

y = A * x
b = y + 0.1y .* random(100)
Iₙ = Matrix{Float64}(I, 100, 100)
λ = 0.11
xλ = inv(A'A + λ^2.0 * Iₙ) * A' * b

The solution is not perfect, but it is free of random error and a reasonable approximation of the true $\rm{x}$.

Optimal Regularization Parameter

The choice of the optimal regularization parameter is not obvious. If we pick $\lambda$ too small, the solution is dominated by noise. If we pick $\lambda$ too large the solution will not approximate the correct solution.

random(n) = @> readdlm("random.txt") vec x -> x[1:n] # hide


y = A * x
b = y + 0.1y .* random(100)
Iₙ = Matrix{Float64}(I, 100, 100)
f(λ) = inv(A'A + λ^2.0 * Iₙ) * A' * b
xλ1 = f(0.01)
xλ2 = f(0.1)
xλ3 = f(10.0)