Iterated Tikhonov Filters

Why Iterate Tikhonov Regularization?

In an inverse problem, recovering a component along a small singular value requires dividing by a small number. Ordinary Tikhonov regularization controls that instability by reducing the response of such components. Iterated Tikhonov applies a regularized correction repeatedly: each step works on the residual left by the previous step. To understand what repetition changes, it is essential to distinguish the residual filter from the multiplier that maps observed data to the estimated solution.

We work with a compact linear operator $A:\mathcal H\to\mathcal G$ between Hilbert spaces, and data $y\in\mathcal G$. Let $(\sigma_j,u_j,v_j)$ be a singular system for the positive singular values: $Av_j=\sigma_j u_j$, $A^*u_j=\sigma_j v_j$, and $\sigma_j>0$. Write $y_j=\langle y,u_j\rangle$. Components of $y$ orthogonal to the closure of the range of $A$ do not contribute to $A^*y$, and components of a solution in the null space of $A$ are not determined by the data.

Definition: Fix a regularization parameter $\alpha>0$ and a positive integer $m$. The $m$-fold iterated Tikhonov estimate starts at $x_0=0$ and is defined recursively by $$ x_k=x_{k-1}+(A^*A+\alpha I)^{-1}A^*(y-Ax_{k-1}),\qquad k=1,\ldots,m. $$ Each step applies an ordinary Tikhonov correction to the current residual $y-Ax_{k-1}$, using the same $\alpha$.

The inverse in this definition is well-defined and bounded: $A^*A$ is nonnegative, so $A^*A+\alpha I$ has spectrum bounded below by $\alpha$. The iteration is therefore well-defined even when $A$ does not have a bounded inverse.

Deriving the Iterated Spectral Filter

Consider one singular coordinate. If $x_{k-1,j}=\langle x_{k-1},v_j\rangle$, then the data residual in direction $u_j$ is $y_j-\sigma_jx_{k-1,j}$. The correction adds $\sigma_j(y_j-\sigma_jx_{k-1,j})/(\sigma_j^2+\alpha)$ to that solution coordinate. Equivalently, the residual is multiplied by $\alpha/(\sigma_j^2+\alpha)$ at each step.

Theorem: For every positive singular value $\sigma_j$, the $m$-fold iterated Tikhonov estimate has coordinate $$ x_{m,j}=\frac{1-\left(\frac{\alpha}{\sigma_j^2+\alpha}\right)^m}{\sigma_j}\,y_j. $$ Writing $t=\sigma_j^2$, its residual filter is $r_{m,\alpha}(t)=\left(\frac{\alpha}{t+\alpha}\right)^m$, its data filter is $q_{m,\alpha}(t)=1-r_{m,\alpha}(t)$, and its solution filter is $$ g_{m,\alpha}(t)=\frac{q_{m,\alpha}(t)}{t} =\frac{1}{t+\alpha}\sum_{k=0}^{m-1}\left(\frac{\alpha}{t+\alpha}\right)^k,\qquad t>0. $$ Thus $x_m=g_{m,\alpha}(A^*A)A^*y$, or, in singular coordinates, $x_m=\sum_j \sigma_jg_{m,\alpha}(\sigma_j^2)y_jv_j$.

Proof. Let $e_{k,j}=y_j-\sigma_jx_{k,j}$ denote the residual coordinate after step $k$. The update gives $$ x_{k,j}=x_{k-1,j}+\frac{\sigma_j}{\sigma_j^2+\alpha}e_{k-1,j}. $$ Consequently, $$ e_{k,j}=y_j-\sigma_jx_{k,j} =e_{k-1,j}-\frac{\sigma_j^2}{\sigma_j^2+\alpha}e_{k-1,j} =\frac{\alpha}{\sigma_j^2+\alpha}e_{k-1,j}. $$ Since $x_{0,j}=0$, we have $e_{0,j}=y_j$, and repeated application yields $e_{m,j}=(\alpha/(\sigma_j^2+\alpha))^m y_j$. For $\sigma_j>0$, the identity $e_{m,j}=y_j-\sigma_jx_{m,j}$ now gives the asserted coordinate formula. Dividing its numerator by $\sigma_j$ gives the data-to-solution multiplier. Finally, the finite geometric-sum identity $1-a^m=(1-a)\sum_{k=0}^{m-1}a^k$, with $a=\alpha/(t+\alpha)$ and $1-a=t/(t+\alpha)$, gives the stated expression for $g_{m,\alpha}$. This also shows that $g_{m,\alpha}(t)$ extends continuously to $t=0$ with value $m/\alpha$. The spectral formula is well-defined as an operator: for $t\leq\alpha$, $\sqrt t\,g_{m,\alpha}(t)\leq m/\sqrt\alpha$; for $t\geq\alpha$, $q_{m,\alpha}(t)\leq1$ gives $\sqrt t\,g_{m,\alpha}(t)=q_{m,\alpha}(t)/\sqrt t\leq1/\sqrt\alpha$. Thus the multiplier on $A^*y$ is bounded. $\square$

Ordinary Tikhonov is the special case $m=1$: $q_{1,\alpha}(t)=t/(t+\alpha)$, $g_{1,\alpha}(t)=1/(t+\alpha)$, and the data-to-solution multiplier is $\sigma/(\sigma^2+\alpha)$. Iteration changes these filters in a specific way; it does not simply replace $\alpha$ by a smaller regularization parameter.

Worked Example: A Two-Step Correction on One Coordinate

Take one singular value $\sigma=2$, data coordinate $y_j=3$, and $\alpha=1$. The first correction is $$ x_{1,j}=0+\frac{2}{2^2+1}(3)=\frac65. $$ The residual is $3-2(6/5)=3/5$. The second correction is $2(3/5)/5=6/25$, so $x_{2,j}=6/5+6/25=36/25=1.44$.

The filter formula gives $q_{2,1}(4)=1-(1/5)^2=24/25$, hence $g_{2,1}(4)=(24/25)/4=6/25$. The data-to-solution multiplier is $\sigma g=12/25$, and $(12/25)(3)=36/25$, matching the direct iteration.

What Happens at Small Singular Values?

The filter comparison is easiest when $t=\sigma^2$ is small relative to $\alpha$. For ordinary Tikhonov, $q_{1,\alpha}(t)=t/(t+\alpha)$ is approximately $t/\alpha$. For fixed $m$, the iterated data filter is approximately $mt/\alpha$, while its residual filter is approximately $1-mt/\alpha$. These statements follow directly from the exact finite-sum relation $$ \frac{q_{m,\alpha}(t)}{q_{1,\alpha}(t)} =\sum_{k=0}^{m-1}\left(\frac{\alpha}{t+\alpha}\right)^k,\qquad t>0. $$ As $t$ decreases to zero, each of the $m$ terms tends to 1, so this ratio tends to $m$. In particular, the solution filter $g=q/t$ tends to $m/\alpha$, and the data-to-solution multiplier $\sigma g$ behaves like $m\sigma/\alpha$.

Comparison: For fixed $\alpha>0$ and positive integer $m$, $0\leq r_{m,\alpha}(t)\leq r_{1,\alpha}(t)$ and $q_{m,\alpha}(t)\geq q_{1,\alpha}(t)$ for $t\geq0$. When $m>1$ and $t>0$, the second inequality is strict. Iteration leaves less residual in each positive singular direction, so its data filter passes more of that direction than ordinary Tikhonov does. Its solution multiplier is correspondingly larger; it is not more suppressive of small-singular-value data components.

To verify the comparison, put $a=\alpha/(t+\alpha)$, which lies in $[0,1]$. Then $r_{m,\alpha}(t)=a^m\leq a=r_{1,\alpha}(t)$, so $q_{m,\alpha}(t)=1-a^m\geq1-a=q_{1,\alpha}(t)$. If $t>0$, then $a<1$, and for $m>1$, $a^m<a$, giving strict inequality. Dividing the data-filter comparison by $\sigma>0$ also compares the solution multipliers.

Worked Example: Small Singular Value, Exact Comparison

Let $\sigma=0.1$, $\alpha=0.25$, and $m=3$, with data coordinate $y_j=1$. Then $t=0.01$ and $\alpha/(t+\alpha)=25/26$. Ordinary Tikhonov has $$ q_{1,\alpha}(t)=1-\frac{25}{26}=\frac1{26},\qquad \frac{q_{1,\alpha}(t)}{\sigma}=\frac5{13}\approx0.3846. $$ For three iterations, $$ q_{3,\alpha}(t)=1-\left(\frac{25}{26}\right)^3 =\frac{17576-15625}{17576} =\frac{1951}{17576}\approx0.1110, $$ so the data-to-solution multiplier is $q_{3,\alpha}(t)/\sigma\approx1.110$. The iterated estimate is larger because the filter passes more of this coordinate. The ratio of the two data filters is $1+25/26+(25/26)^2$, about $2.89$, approaching $m=3$ as $t/\alpha$ gets smaller.

At the other end of the spectrum, if $t$ is large relative to $\alpha$, then $\alpha/(t+\alpha)$ is small. Its $m$th power is smaller still, so $q_{m,\alpha}(t)$ is close to 1. The data coordinate is then nearly passed through the inverse multiplier $1/\sigma$. Iteration primarily changes the transition and small-$t$ behavior; it does not change the limiting inverse response at large $t$.

Worked Example: Applying the Filter to Two Coordinates

Suppose a two-coordinate problem has singular values $0.5$ and $2$, data coordinates $2$ and $-1$, regularization parameter $\alpha=0.25$, and $m=2$. For $\sigma=0.5$, $t=0.25$, so $\alpha/(t+\alpha)=1/2$. Thus $q=1-(1/2)^2=3/4$, and the solution multiplier is $q/\sigma=3/2$. This coordinate of the estimate is $(3/2)(2)=3$.

For $\sigma=2$, $t=4$, so $\alpha/(t+\alpha)=1/17$. Thus $q=1-(1/17)^2=288/289$, and the solution multiplier is $q/\sigma=144/289$. This coordinate is $-144/289$. Therefore $x_2$ has coordinates $(3,-144/289)$. As a check, ordinary Tikhonov has multipliers $0.5/(0.25+0.25)=1$ and $2/(4+0.25)=8/17$, yielding coordinates $(2,-8/17)$. Both iterated coordinates have larger magnitude than their ordinary-Tikhonov counterparts, consistent with the filter comparison.

Interpretation and a Common Pitfall

Iteration can reduce the residual substantially: after $m$ steps, its factor in direction $u_j$ is $(\alpha/(\sigma_j^2+\alpha))^m$. But “more residual reduction” and “more attenuation of the corresponding solution coordinate” are not interchangeable claims. The data filter is $q_{m,\alpha}$, and the solution multiplier is $q_{m,\alpha}(\sigma^2)/\sigma$; both are larger than their ordinary-Tikhonov counterparts for $m>1$ and $\sigma>0$. Small singular values are still suppressed relative to exact inversion, because $q_{m,\alpha}(\sigma^2)$ tends to zero with $\sigma$. They are simply suppressed less, for fixed $\alpha$, than under ordinary Tikhonov.

Common mistake: “The residual factor is raised to the $m$th power, so the iterated method must suppress small singular-value solution coordinates more strongly.” This confuses the residual factor $r_{m,\alpha}$ with the data filter $q_{m,\alpha}=1-r_{m,\alpha}$.
Correction: Raising the residual factor to a power makes it smaller, so the data filter becomes larger. The exact ratio $q_{m,\alpha}/q_{1,\alpha}=\sum_{k=0}^{m-1}(\alpha/(t+\alpha))^k$ shows how much larger: it tends to $m$ as $t$ tends to zero.
Another caution: At $t=0$, $g_{m,\alpha}(t)$ has the finite continuous extension $m/\alpha$, but a zero singular value has no data-to-solution inverse coordinate. The actual update uses $A^*y$, which has no component along $\ker(A)$.

Key takeaway: The $m$-fold iterated Tikhonov residual filter is $(\alpha/(t+\alpha))^m$, and the data filter is its complement. Repetition reduces residuals more and passes more signal in every positive singular direction than ordinary Tikhonov; for small singular values, the solution multiplier is approximately $m\sigma/\alpha$, rather than $\sigma/\alpha$.

Check Your Understanding

Use the filter definitions and comparisons above to answer the following.

For a positive singular value $\sigma$, what factor multiplies its data residual after $m$ iterations?
With $\alpha=1$, $t=1$, and $m=2$, calculate the residual filter and the data filter.
Why does $q_{m,\alpha}(t)\geq q_{1,\alpha}(t)$ imply a larger solution multiplier for each positive singular value?
As $t$ tends to zero with $\alpha$ and $m$ fixed, what are the limits of $g_{m,\alpha}(t)$ and $q_{m,\alpha}(t)/q_{1,\alpha}(t)$?
Explain why a zero singular value does not give a data-to-solution coordinate, even though $g_{m,\alpha}(t)$ has a finite extension at $t=0$.

What You'll Learn