Spectral regularization and parameter choice

Since $A$ is compact and $A^{\ast }$ is bounded, [F3] makes $A^{\ast }A$ compact. The same lemma makes $A{A}^{\ast }$ compact because it is the composition of the compact map $A$ after the bounded map $A^{\ast }$. By [F4], $(A^{\ast }A{)}^{\ast }=A^{\ast }A$ and $(A{A}^{\ast }{)}^{\ast }=A{A}^{\ast }$, so both are self-adjoint.

For every $x\in X$, [F1] and [F7] give $⟨A^{\ast }Ax,x{⟩}_X=⟨Ax,Ax{⟩}_Y=\mathrm{\parallel }Ax{\mathrm{\parallel }}^2\ge 0$. For every $y\in Y$, the same identity gives $⟨A{A}^{\ast }y,y{⟩}_Y=⟨A^{\ast }y,A^{\ast }y{⟩}_X=\mathrm{\parallel }{A}^{\ast }y{\mathrm{\parallel }}^2\ge 0$. Thus both operators are positive.

If $Ax=0$, then $A^{\ast }Ax=0$. Conversely, if $A^{\ast }Ax=0$, step 1.2 implies $\mathrm{\parallel }Ax{\mathrm{\parallel }}^2=⟨A^{\ast }Ax,x⟩=0$, so positive definiteness [F6] gives $Ax=0$. Hence $\ker (A^{\ast }A)=\ker A$. The same argument with $A^{\ast }$ in place of $A$ proves $\ker (A{A}^{\ast })=\ker A^{\ast }$.

Write $x=P_{0}x+m$ with $m\in (\ker S{)}^{\perp }$ using [F1]–[F2]. For a finite $F\subseteq \mathrm{\Sigma }$, the vectors $P_{t}x$ are mutually orthogonal by [F1] and [F3]; their sum is the orthogonal projection of $m$ onto the finite orthogonal sum of the spaces $E_t$. Pythagoras therefore gives ${\sum }_{t\in F}\mathrm{\parallel }{P}_{t}x{\mathrm{\parallel }}^2\le \mathrm{\parallel }m{\mathrm{\parallel }}^2\le \mathrm{\parallel }x{\mathrm{\parallel }}^2$.

The set $\{\mathrm{\mid }q(t)\mathrm{\mid }:t\in \{0\}\cup \mathrm{\Sigma }\}$ is nonempty and bounded, so [F8] gives its finite supremum $M_q$. For finite $F\subseteq \mathrm{\Sigma }$, orthogonality and step 1.1 give $\mathrm{\parallel }{\sum }_{t\in F}q(t)P_{t}x{\mathrm{\parallel }}^2={\sum }_{t\in F}\mathrm{\mid }q(t){\mathrm{\mid }}^2\mathrm{\parallel }{P}_{t}x{\mathrm{\parallel }}^2\le M_q^2\mathrm{\parallel }x{\mathrm{\parallel }}^2$. Put $a_t=\mathrm{\parallel }{P}_{t}x{\mathrm{\parallel }}^2$. Step 1.1 bounds every finite subsum of $(a_t)$ by $\mathrm{\parallel }x{\mathrm{\parallel }}^2$, so [F7] gives a finite sum with tails arbitrarily small. If $M_q=0$, every weighted spectral term vanishes. Otherwise, for any $\epsilon >0$, choose a finite $F_0\subseteq \mathrm{\Sigma }$ with ${\sum }_{t\mathrm{\notin }{F}_0}{a}_t<(\epsilon \mathrm{/}{M}_q{)}^2$. For finite $F,H\supseteq F_0$, orthogonality gives $\mathrm{\parallel }{\sum }_{t\in F}q(t)P_{t}x-{\sum }_{t\in H}q(t)P_{t}x{\mathrm{\parallel }}^2\le M_q^2{\sum }_{t\in F\mathrm{△}H}{a}_t<{\epsilon }^2$. Thus the finite-subset net is Cauchy; completeness of $H$ makes it norm-convergent.

Finite partial sums are linear in $x$, so taking their norm limits shows that $x\mapsto q(0)P_{0}x+{\sum }_{t\in \mathrm{\Sigma }}q(t)P_{t}x$ is linear. Orthogonality and convergence of the finite-subset net give $\mathrm{\parallel }{Q}_{q}x{\mathrm{\parallel }}^2=\mathrm{\mid }q(0){\mathrm{\mid }}^2\mathrm{\parallel }{P}_{0}x{\mathrm{\parallel }}^2+{\sum }_{t\in \mathrm{\Sigma }}\mathrm{\mid }q(t){\mathrm{\mid }}^2\mathrm{\parallel }{P}_{t}x{\mathrm{\parallel }}^2\le M_q^2\mathrm{\parallel }x{\mathrm{\parallel }}^2$, with the nonnegative sum understood as in [F7]. Thus $Q_q$ is bounded and [F6] gives $\mathrm{\parallel }{Q}_q\mathrm{\parallel }\le M_q$.

Put $G=A^{T}A+\lambda I$. For every nonzero $z$, $z^{T}Gz=\mathrm{\parallel }Az{\mathrm{\parallel }}^2+\lambda \mathrm{\parallel }z{\mathrm{\parallel }}^2>0$, so $G$ is positive definite and invertible; [F1]–[F2] give $x_{\lambda }=G^{-1}{A}^{T}b$. The SVD [F3] shows $H_{\lambda }$ has eigenvalues ${\sigma }_j^2\mathrm{/}({\sigma }_j^2+\lambda )$ on the left singular vectors and zero on their orthogonal complement. Every eigenvalue is in $[0,1)$, so $I-H_{\lambda }$ is positive definite and $1-(H_{\lambda }{)}_{ii}>0$ for every coordinate vector.

Fix $i$, write $a=a_i$, $h=a^T{G}^{-1}a=(H_{\lambda }{)}_{ii}$, and $G_{-}=G-a{a}^T=A_{-i}^T{A}_{-i}+\lambda I$. The same positive-definiteness argument makes $G_{-}$ invertible. Direct multiplication verifies $G_{-}^{-1}=G^{-1}+G^{-1}a{a}^T{G}^{-1}\mathrm{/}(1-h)$; the denominator is positive by step 1.1. Substituting $A^{T}b-a{b}_i$ for the reduced right-hand side and simplifying gives $x_{\lambda }^{(-i)}=x_{\lambda }-G^{-1}a((I-H_{\lambda })b{)}_i\mathrm{/}(1-h)$. Multiplying by $a^T$ and using $a^T{G}^{-1}a=h$ yields the stated deleted residual identity.

Squaring the deleted residual identity and summing over $i$ proves the exact leave-one-out sum; dividing by $m$ gives its mean. Replacing all $h_i=(H_{\lambda }{)}_{ii}$ by $\bar{h}=m^{-1}{\sum }_i{h}_i=\mathrm{tr}(H_{\lambda })\mathrm{/}m$ in that mean gives $m^{-1}{\sum }_i((I-H_{\lambda })b{)}_i^2\mathrm{/}(1-\bar{h}{)}^2=\mathrm{\parallel }(I-H_{\lambda })b{\mathrm{\parallel }}^2\mathrm{/}[m(1-\bar{h}{)}^2]$, the displayed design convention for GCV. If all $h_i=\bar{h}$, the exact mean has the same formula. Conversely, if the exact mean and GCV agree for every $b$, then $z=(I-H_{\lambda })b$ ranges over all of ${\mathbb{R}}^m$ because step 1.1 makes $I-H_{\lambda }$ positive definite and hence invertible. Equality of the two diagonal quadratic forms for all $z$ forces $(1-h_i{)}^{-2}=(1-\bar{h}{)}^{-2}$ for every $i$; all terms $1-h_i$ are positive by step 1.1, so $h_i=\bar{h}$ for every $i$. This is equality as formulas in the data; individual data vectors may yield accidental equality.

For each fixed $n$, $\mathrm{\mid }{x}_n^{(k)}-x_n^{(j)}\mathrm{\mid }\le \mathrm{\parallel }{x}^{(k)}-x^{(j)}{\mathrm{\parallel }}_2$, so $(x_n^{(k)}{)}_k$ is Cauchy in $\mathbb{R}$. By [F2] let $x_n={\lim }_k{x}_n^{(k)}$; each coordinate limit is unique, so these values define a sequence $x$. To see $(x^{(k)})$ is bounded, choose $K$ so that $\mathrm{\parallel }{x}^{(k)}-x^{(j)}{\mathrm{\parallel }}_2<1$ for all $k,j\ge K$; then $\mathrm{\parallel }{x}^{(k)}{\mathrm{\parallel }}_2\le \mathrm{\parallel }{x}^{(K)}{\mathrm{\parallel }}_2+1$ for $k\ge K$, and the finitely many earlier norms have a finite maximum. Thus $\mathrm{\parallel }{x}^{(k)}{\mathrm{\parallel }}_2\le M$ for some finite $M$. For any finite $F\subset \mathbb{N}$, ${\sum }_{n\in F}\mathrm{\mid }{x}_n{\mathrm{\mid }}^2={\lim }_k{\sum }_{n\in F}\mathrm{\mid }{x}_n^{(k)}{\mathrm{\mid }}^2\le M^2$. Taking the supremum over finite $F$ gives ${\sum }_n\mathrm{\mid }{x}_n{\mathrm{\mid }}^2\le M^2$, hence $x\in {\mathrm{\ell }}^2$.

Let $\epsilon >0$. Choose $K$ such that $\mathrm{\parallel }{x}^{(k)}-x^{(j)}{\mathrm{\parallel }}_2<\epsilon \mathrm{/}2$ for all $k,j\ge K$. Fix $k\ge K$. For every finite $F\subset \mathbb{N}$, coordinate convergence and continuity of finite sums give ${\sum }_{n\in F}\mathrm{\mid }{x}_n^{(k)}-x_n{\mathrm{\mid }}^2={\lim }_{j\to \mathrm{\infty }}{\sum }_{n\in F}\mathrm{\mid }{x}_n^{(k)}-x_n^{(j)}{\mathrm{\mid }}^2\le {\epsilon }^2\mathrm{/}4$. Taking the supremum over finite $F$ yields $\mathrm{\parallel }{x}^{(k)}-x{\mathrm{\parallel }}_2\le \epsilon \mathrm{/}2<\epsilon$. Thus $x^{(k)}\to x$ in ${\mathrm{\ell }}^2$. Every Cauchy sequence converges, proving completeness and, by [F3], the Hilbert-space claim.

The input grid has finite length $q$. Each positive parameter gives an invertible matrix $A^{T}A+{\lambda }_{j}I$ because $z^T(A^{T}A+{\lambda }_{j}I)z=\mathrm{\parallel }Az{\mathrm{\parallel }}^2+{\lambda }_j\mathrm{\parallel }z{\mathrm{\parallel }}^2>0$ for $z\mathrm{\ne }0$. Thus the solution and residual calculations in [F1] are defined; a finite number of scalar norm, positivity, logarithm, determinant, and arithmetic operations determine the retained points and curvatures. Every loop terminates after at most $q$ iterations.

If fewer than three points remain, [F1] returns the prescribed failure status. Otherwise the finite nonempty set of assigned curvature values has a maximum; the finite set of indices attaining it has a least original grid index. The tie rule in [F1] therefore returns a maximizer and terminates.

By [F1], $\ker S=\ker A$. Thus $x^{†}\perp \ker S$, and the spectral decomposition in [F2] gives $x^{†}={\sum }_{t\in \mathrm{\Sigma }}{P}_t{x}^{†}$ as a norm-convergent finite-subset net.

Since $y=A{x}^{†}$, $A^{\ast }y=S{x}^{†}$. The filter in [F3] therefore gives $x_{\lambda }=g_{\lambda }(S)S{x}^{†}$, whose coefficient on $E_t$ is $t\mathrm{/}(t+\lambda )$. The coefficient on $\ker S$ is zero. Subtracting from the expansion in step 1.1 proves the displayed bias identity.

Given $\epsilon >0$, choose a finite $F\subseteq \mathrm{\Sigma }$ so the norm of the tail ${\sum }_{t\mathrm{\notin }F}{P}_t{x}^{†}$ is less than $\epsilon \mathrm{/}2$; this is possible because the finite-subset net in [F2] converges. Each scalar $\lambda \mathrm{/}(t+\lambda )$ for $t\in F$ tends to zero, so for all sufficiently small $\lambda >0$ the norm of the finite part of the bias sum is less than $\epsilon \mathrm{/}2$. The tail has norm at most its unfiltered norm because $0\le \lambda \mathrm{/}(t+\lambda )\le 1$. The triangle inequality then gives $\mathrm{\parallel }{x}^{†}-x_{\lambda }\mathrm{\parallel }<\epsilon$.

Put $T=A{A}^{\ast }$. The normal equation $(A^{\ast }A+\lambda I)x_{\lambda }(g)=A^{\ast }g$ follows from [F1] and [F5]. Applying $A$ gives $(T+\lambda I)A{x}_{\lambda }(g)=Tg$, whence $g-A{x}_{\lambda }(g)=\lambda (T+\lambda I{)}^{-1}g\mathrm{.}$ The inverse exists by [F5] applied to $A^{\ast }:Y\to X$. The multiplier $q(s)=\lambda \mathrm{/}(s+\lambda )$, including $q(0)=1$, has supremum at most one; [F4] gives the operator norm bound.

Let $P_s$ be the eigenspace projections of $T$ for its nonzero eigenvalues $s$, and $P_0$ the projection onto $\ker T$. By [F3] and step 1.1, $r_g(\lambda {)}^2=\mathrm{\parallel }{P}_{0}g{\mathrm{\parallel }}^2+{\sum }_{s>0}{\left(\frac{\lambda }{s+\lambda }\right)}^2\mathrm{\parallel }{P}_{s}g{\mathrm{\parallel }}^2\mathrm{.}$ Each coefficient $\lambda \mathrm{/}(s+\lambda )$ is nondecreasing in $\lambda$; therefore the sum and its nonnegative square root are nondecreasing.

On any compact interval $[a,b]\subset (0,\mathrm{\infty })$, $\mid \frac{\lambda }{s+\lambda }-\frac{\mu }{s+\mu }\mid =\frac{s\mathrm{\mid }\lambda -\mu \mathrm{\mid }}{(s+\lambda )(s+\mu )}\le \frac{\mathrm{\parallel }T\mathrm{\parallel }}{a^2}\mathrm{\mid }\lambda -\mu \mathrm{\mid }$ for $0\le s\le \mathrm{\parallel }T\mathrm{\parallel }$ and $\lambda ,\mu \in [a,b]$. Hence the residual multipliers vary continuously in operator norm, which makes $\lambda \mapsto \lambda (T+\lambda I{)}^{-1}g$ and its norm continuous.

As $\lambda \downarrow 0$, the coefficient on $P_{0}g$ stays $1$ and each coefficient on $P_{s}g$ for $s>0$ tends to zero. For a given $\epsilon >0$, choose a finite spectral set whose orthogonal tail of $g$ has norm below $\epsilon \mathrm{/}2$; the multipliers have norm at most one, and the finite positive-eigenvalue sum tends to zero. Thus the residual tends to $P_{0}g$. As $\lambda \to \mathrm{\infty }$, the bound ${\sup }_{0\le s\le \mathrm{\parallel }T\mathrm{\parallel }}\mathrm{\mid }1-\lambda \mathrm{/}(s+\lambda )\mathrm{\mid }\le \mathrm{\parallel }T\mathrm{\parallel }\mathrm{/}\lambda$ shows that the residual tends to $g$.

The adjoint identity [F8] gives $(\mathrm{Ran}A{)}^{\perp }=\ker A^{\ast }$: orthogonality to every $Ax$ is equivalent to $⟨x,A^{\ast }z⟩=0$ for every $x$, which is equivalent to $A^{\ast }z=0$. Taking closures leaves the orthogonal complement unchanged. By [F2], $\ker T=\ker A^{\ast }$, and [F7] identifies $\mathrm{dist}(g,\overline{\mathrm{Ran}A})=\mathrm{\parallel }{P}_{0}g\mathrm{\parallel }$. Step 2.3 gives the claimed lower endpoint, while its upper endpoint is $\mathrm{\parallel }g\mathrm{\parallel }$.

On the eigenspace $E_t$ with $t>0$, the exact-data regularization error $x^{†}-x_{\lambda }$ multiplies the source coefficient by $\lambda t^{\nu }\mathrm{/}(t+\lambda )$; on $\ker S$ the source multiplier $t^{\nu }$ is zero. Orthogonality in [F2] and the definition in [F3] therefore give $x^{†}-x_{\lambda }={\sum }_{t\in \mathrm{\Sigma }}\frac{\lambda t^{\nu }}{t+\lambda }{P}_{t}w\mathrm{.}$

Put $s=t\mathrm{/}\lambda$. If $0<s\le 1$, then $s^{\nu }\mathrm{/}(1+s)\le s^{\nu }\le 1$. If $s\ge 1$, then $s^{\nu }\mathrm{/}(1+s)\le s^{\nu -1}\le 1$ because $\nu \le 1$. Thus $\lambda t^{\nu }\mathrm{/}(t+\lambda )={\lambda }^{\nu }{s}^{\nu }\mathrm{/}(1+s)\le {\lambda }^{\nu }$ for every $t>0$.

Applying the orthogonal norm formula to step 1.1 and using step 1.2 gives $\mathrm{\parallel }{x}^{†}-x_{\lambda }{\mathrm{\parallel }}^2={\sum }_{t\in \mathrm{\Sigma }}{\mid \frac{\lambda t^{\nu }}{t+\lambda }\mid }^2\mathrm{\parallel }{P}_{t}w{\mathrm{\parallel }}^2\le {\lambda }^{2\nu }{\sum }_{t\in \mathrm{\Sigma }}\mathrm{\parallel }{P}_{t}w{\mathrm{\parallel }}^2\le {\lambda }^{2\nu }\mathrm{\parallel }w{\mathrm{\parallel }}^2\mathrm{.}$ Taking nonnegative square roots proves the bound. The scalar inequality in step 1.2 is the qualification condition from [F3] with $C_{\nu }=1$.

Each ${\lambda }_j>0$ makes $A^{T}A+{\lambda }_{j}I$ positive definite, so all matrix inverses in [F1] exist. By [F2], the denominator of each score is strictly positive; the numerator is a finite squared norm. Thus each score is a finite real number. The procedure performs finitely many operations for finitely many grid entries and terminates.

A nonempty finite list of real scores has a minimum, and a nonempty finite set of indices attaining the minimum has a least element. By the selection rule in [F1], the algorithm returns that index and hence a least-score grid parameter. No part of this finite comparison identifies the score with the unknown solution error.