Statistical approach to linear inverse problems

TL;DR

This paper discusses a classical statistical approach to linear inverse problems with additive noise, emphasizing stability, feasible solution regions, and a near-optimal nonlinear filtering method without Bayesian assumptions.

Contribution

It introduces a non-Bayesian statistical framework for stable inverse solutions, characterizes the feasible solution region, and presents a nonlinear filter close to the optimal Kolmogorov-Wiener filter.

Findings

Feasible solution region is an elongated ellipsoid shaped by Fisher information.

Spectrum of Fisher matrix determines inverse problem stability.

Proposes a nonlinear filter near the Kolmogorov-Wiener optimal filter.

Abstract

The main features of the statistical approach to inverse problems are described on the example of a linear model with additive noise. The approach does not use any Bayesian hypothesis regarding an unknown object; instead, the standard statistical requirements for the procedure for finding a desired object estimate are presented. In this way, it is possible to obtain stable and efficient inverse solutions in the framework of classical statistical theory. The exact representation is given for the feasible region of inverse solutions, i.e., the set of inverse estimates that are in agreement, in the statistical sense, with the data and available a priory information. The typical feasible region has the form of an extremely elongated hole ellipsoid, the orientation and shape of which are determined by the Fisher information matrix. It is the spectrum of the Fisher matrix that provides an exhaustive description of the stability of the inverse problem under consideration. The method of constructing a nonlinear filter close to the optimal Kolmogorov-Wiener filter is presented.

Equations

$$y_0=H{x}_0,$$

$$y_0=H{x}_0+\xi ,$$

$$x_{\gamma }=\arg \underset{x}{\min }\left[||y_0-Hx-a|{|}^2+\gamma \Phi (x)\right],$$

$$\Phi (x)=||x|{|}^2,$$

$$\left\{\begin{array}[]{ll}y_{0}=Hx_{0}+\xi,&\\ \langle\xi\rangle=a,\qquad{\rm cov}(\xi)\equiv\langle(\xi-a)(\xi-a)^{T}\rangle=C,\end{array}\right.$$

$$z_0=C^{-1/2}(y_0-a),\eta =C^{-1/2}(\xi -a),A=C^{-1/2}H,$$

$$\left\{\begin{array}[]{ll}z_{0}=Ax_{0}+\eta,&\\ \langle\eta\rangle=0,\qquad{\rm cov}(\eta)=E_{m},\end{array}\right.$$

$$A=U\Delta V^T,$$

$$U^{T}U=E_n,\Delta =\mathrm{diag}(\delta ),V^{-1}=V^T.$$

$$x_0=V{p}_0=\sum _{k=1}^n{p}_{0k}{v}_k,p_0=V^T{x}_0,$$

$$\varphi \equiv U^T{z}_0,\zeta \equiv U^T\eta ,$$

$$\left\{\begin{array}[]{ll}\phi=\Delta p_{0}+\zeta,&\\ \langle\zeta\rangle=0,\qquad{\rm cov}(\zeta)=E_{n}.&\\ \end{array}\right.$$

$$\Theta (y_0|x)\equiv \Vert \varphi -\Delta p{\Vert }^2,$$

$$t_{1-{\alpha }_2}^{(n)}\le \Theta (y_0|x)\le t_{1-{\alpha }_1}^{(n)}.$$

$$\Vert \varphi -\Delta p{\Vert }^2=t_{1-\alpha }^{(n)}.$$

$$x_{\ast }=(A^{T}A{)}^{-1}{A}^T{z}_0$$

$$x_{\ast }=V{p}_{\ast }=\sum _{k=1}^n{p}_{\ast k}{v}_k,p_{\ast }={\Delta }^{-1}\varphi .$$

$$⟨p_{\ast }⟩=p_0,\mathrm{cov}(p_{\ast })={\Lambda }^{-1},$$

$$\Lambda \equiv {\Delta }^2=\mathrm{diag}({\lambda }_1,\dots ,{\lambda }_n),{\lambda }_k={\delta }_k^2.$$

$$(x-x_{\ast }{)}^{T}I(x-x_{\ast })=t_{1-\alpha }^{(n)},I=A^{T}A=V\Lambda V^T.$$

$$I_{ik}(x_{0})\equiv\biggl{\langle}\frac{\partial}{\partial x_{0i}}\,\ln f(y_{0}|x_{0})\,\frac{\partial}{\partial x_{0k}}\,\ln f(y_{0}|x_{0})\biggr{\rangle},\,\quad i,k=1,\ldots,n,$$

$$x_w\equiv \sum _{k=1}^n{w}_k{p}_{\ast k}{v}_k=VW{p}_{\ast },W=\mathrm{diag}(w).$$

$${\epsilon }_w^2\equiv ⟨\Vert x_w-x_0{\Vert }^2⟩=\sum _{k=1}^n\left[w_k^2/{\lambda }_k+(1-w_k{)}^2{p}_{0k}^2\right].$$

$${\tilde{w}}_k={\lambda }_k{p}_{0k}^2/(1+{\lambda }_k{p}_{0k}^2),k=1,2,\dots ,n,$$

$${\tilde{x}}_w=\sum _{k=1}^n{\tilde{w}}_k{p}_{\ast k}{v}_k=V{\tilde{p}}_w,{\tilde{p}}_w=\tilde{W}(p_0)p_{\ast }.$$

$$\Vert (W-E_n)\varphi {\Vert }^2=t_{1-\alpha }^{(n)}.$$

$$⟨\Vert \tilde{W}(p_0)p_{\ast }-p_0{\Vert }^2⟩={\tilde{\epsilon }}_w^2(p_0).$$

$$\tilde{W}(p)=\mathrm{diag}\left[\tilde{w}(p)\right],{\tilde{w}}_k(p)={\lambda }_k{p}_k^2/(1+{\lambda }_k{p}_k^2),$$

$$⟨\Vert \tilde{W}(p)p_{\ast }-p{\Vert }^2⟩=\min .$$

$$\left\{\begin{array}[]{ll}\|[\tilde{W}(p)-E_{n}]\,\phi\,\|^{2}=t_{1-\alpha}^{(n)},&\\ &\\ \|\tilde{W}(p)p_{*}-p\,\|^{2}=\min.&\\ \end{array}\right.$$

Full text

Keywords: Inverse problems, image restoration

1. Introduction

In the usual sense, the inverse problem is to find the object $x_{0}$ from equation

[TABLE]

where the image $y_{0}$ describes the measured data, and $H$ is a known procedure (see, e.g., Tikhonov and Arsenin [1977], Bertero [1986]). We consider here the case when $H$ is a linear integral operator in the finite-dimensional space, so the unknown object $x_{0}$ and the observed image $y_{0}$ can be treated as the ${n\times 1}$ and ${m\times 1}$ vectors, respectively, while the point spread function (PSF) $H$ is the ${m\times n}$ matrix.

Inverse problems are especially characteristic for astronomy, which is still dealt predominantly with the interpretation of passive experiment (Feigelson and Babu [2003, 2012]). In recent years, classical problems of this kind have been supplemented by tasks associated with the creation of an early Universe model based on microwave background measurements and study of distant galaxies.

In practice, even if the unique solution of the problem (1) exists for any $y_{0}$, i.e., the problem is well-posed111For a linear problem, the continuity of the inverse mapping is a consequence of the stated requirements. in the sense of Hadamard [1923], the solution can be strongly unstable. The latter term means that relative error propagation from the image to the solution can be very large. Indeed, the data inevitably are randomly noised, so, for a model with an additive noise, one should rewrite equation (1) in the form

[TABLE]

where $\xi$ is an unknown random noise pattern. Only the mean value of noise, $a$, and its variance are known usually from the preliminary measurements. Thus, if we shall try to minimize some kind of the misfit, e.g. $||y_{0}-Hx-a||^{2}$, to find an appropriate inverse solution $x_{*}$, we obtain, as a rule, the function with huge oscillations, because the least square (in general, the maximum likelihood) solution $x_{*}$ is compelled to ‘explain’ sharp random noise fluctuations in the observed pattern $y_{0}$ purely by oscillations of the object’s profile. In view of smoothing nature of the operator $H$, the amplitude of these oscillations should be large.

One can find in the literature descriptions a number of methods aimed to reach the stable inverse solutions (see, e.g., Press et. al. [1992], Jansson [1997], Evans and Stark [2002]); the most widely applied now are the maximum entropy method (Janes [1957a,b], Burg [1967], Narayan and Nituananda [1986]) and the regularization method (Phillips [1962], Tikhonov [1963a,b], Tikhonov and Arsenin [1977]). Both approaches proceed from minimizing the sum of two functionals:

[TABLE]

where the first term, the misfit, measures the agreement of a trial inverse solution $x$ to the model (2), while $\Phi(x)$, the stabilizing or regularizing functional, describes some kind of “smoothness” of the desired solution. The regularization parameter $\gamma$ is introduced here to provide the trade-off between the accuracy and smoothness of the inverse solution.

It can be easily shown that the requirement (3) is equivalent to the Bayesian way of estimation given a priory information about the probabilistic ensemble of the allowable objects. The way is quite consistent and more efficient, comparing to the classic (“Simpsonian”, according to Eisenhart [1964]) estimation, but only if we really have the needed a priory information, i.e., if $\Phi(x)$ is known. Since this happens comparatively rarely, some intuitive forms of $\Phi(x)$ are usually applied, in particular, the quadratic norm

[TABLE]

or more general norm in the Sobolev’s space, or one of the (inequivalent to each other) ‘entropy’ presentations. From the statistical point of view, this way corresponds to the introduction of the Bayes’s hypothesis that was many times criticized due to unavoidable subjectivity and inherent contradictions (see, e.g., Feller [1957], Fisher [1959], Rao [1973], Cox and Hinkley [1974], Szekely [1986]). The frequently stated dissatisfaction in relation to the Bayesian hypothesis prompted Press et. al. [1992], p. 808, rather figuratively express the moods: “Courts have consistently held that academic license does not extend to shouting “Bayesian”  in a crowded hall.” We note only that just the wide variety of forms the regularizing functional $\Phi(x)$ that were proposed for a same problem clearly shows the absence of a natural form of a priory information about the searched object.

Meanwhile, it is possible to obtain the stable inverse solutions in a framework of the classic statistical theory. Some features of this way were described by Terebizh [1995a,b, 2003, 2004]; in full extent the approach is presented in the book Terebizh [2005], which also includes the statistical treatment of the other widely used inverse methods along with corresponding numerical algorithms. Both linear and non-linear models of data formation were considered; the first of them allows the quite general discussion, while the latter one is represented by the actual now phase problem and the long-standing problem of the time series spectral estimation.

The main purpose of the below consideration is to give a brief outline of the statistical approach within the linear data formation model with an additive noise. This model underlies more complicated cases and has very wide practical applications.

2. Statistical formulation of the inverse problem

Let us define the general linear model by equations

[TABLE]

where $n$-vector $x_{0}$ is an object, the ${m\times n}$ PSF matrix $H$ is assumed to be known, as well as the observed image $m$-vector $y_{0}$ (${m\geq n}$), and $\xi$ is the random noise pattern with the known mean level $a$ and the positive definite covariance ${m\times m}$ matrix $C$. It is assumed that the original vectors have the form of columns; the angle brackets mean averaging on the probabilistic ensemble. Evidently, the density distribution $f(y|x_{0})$ of the image is defined by the corresponding distribution of the noise.

The sought inverse solution $\tilde{x}$ is considered as a statistical estimate of the deterministic object $x_{0}$ given its image, a PSF, properties of the noise, and available a priori information about the object. Being a function of the stochastic image $y_{0}$, the inverse solution $\tilde{x}$ is also a random vector, as a rule, with the mutually dependent components. To define properly the notion of quality of a trial estimate $x$, we should carry out two preliminary procedures.

Firstly, correlations in the image $y_{0}$ should be eliminated, because mutual correlations of the noise components $\{\xi_{j}\}_{j=1}^{m}$ does not allow applying the direct definition of the misfit in the form $||y_{0}-H\tilde{x}-a||^{2}$.

Secondly, it is desirable to reduce the misfit vector dimension to the object’s length $n$. Indeed, it is possible to reach the quite good agreement of an image $y=Hx+\xi$ produced by some trial object $x$ with the observed image $y_{0}$ at the expense of good fit in the lengthy ‘wings’ of the images, especially when ${m\gg n}$. In fact, the wings represent mostly the noise patterns, while we have to fit primarily the part of the image that is caused by the smoothed object.

The transition to the independent data set is based on the known linear transform

[TABLE]

which converts the general model (5) to the standard model

[TABLE]

where $E_{m}$ is the unit ${m\times m}$ matrix. The matrix $C^{-1/2}$ in (6) is inverse to the square root of $C$; since the covariance matrix was assumed positive definite, its spectrum is positive and the square root $C^{1/2}$ exists. Thus, we have now the relatively more simple linear model with an additive white noise $\eta$ and the PSF matrix $A$.

The natural way to reach the second goal proceeds from the singular value decomposition (SVD) of the matrix $A$ (see, e.g., Golub and Van Loan [1989], Press et. al. [1992]). Assume that ${{\rm rank}(A)=n}$. Then

[TABLE]

where $U$ is an ${m\times n}$ column-orthogonal matrix, $\Delta$ is a diagonal ${n\times n}$ matrix with positive singular values ${\delta=[\delta_{1},\ldots,\delta_{n}]^{T}}$ of $A$, placed in the order of their decrease, and an ${n\times n}$ matrix ${V=[v_{1},\ldots,v_{n}]}$ is orthogonal:

[TABLE]

The corresponding decomposition of the object $x_{0}$ in the eigenvectors system $\{v_{k}\}$, namely

[TABLE]

defines the vector $p_{0}$ of the object’s principal components. Like the familiar Fourier coefficients, the principal components are often easier to recover than the object itself. The multiplication of (7) by $U^{T}$ is similar to the application of the Fourier transform. Designating $n$-vectors

[TABLE]

we obtain from (7) a final $n$-dimensional representation of the linear model:

[TABLE]

As was said above, the advantages of use of the ‘refined image’ $\phi$ of length $n$ are especially appreciable when ${m\gg n}$.

3. Feasible Region

Assume, for simplicity, that the noise $\xi$ is a Gaussian deviate. Then $\{\zeta_{k}\}$ in (12) are independent Gaussian deviates with zero mean value and unit variance, and the random variable $\|\phi-\Delta p_{0}\|^{2}=\sum_{k=1}^{n}\zeta_{k}^{2}$ has a $\chi^{2}$-distribution with $n$ degrees of freedom (Cramér [1946], Chapter 18). This result allows us to introduce a similar random variable, namely the misfit

[TABLE]

as a measure of the quality of a trial object’s estimate ${x=Vp}$.

Let ${t_{\gamma}^{(n)}\geq 0}$ be a quantile of the $\chi_{n}^{2}$ distribution $P_{n}(t)$, that is the root of equation ${P_{n}(t)=\gamma}$. Just as is usually done in mathematical statistics (Cramér [1946]), we shall choose the appropriate boundary significance levels for an inverse solution $\alpha_{1}$ and $\alpha_{2}$ ($0\leq\alpha_{1}\leq\alpha_{2}\leq 1$). By definition, a trial object’s estimate $x$ is called feasible, if

[TABLE]

We simply require of a feasible estimate $x$ that its image $y(x)$ should have moderate deviation, in the statistical sense, from the observed image $y_{0}$. Inequalities (14) define the feasible region (FR), consisting of all the object’s estimates $\{x\}$ that have feasible agreement with the data. It is convenient to call $x$ the estimate of significance level $\alpha$, if the misfit $\Theta(y_{0}|x)=t_{1-\alpha}^{(n)}$, that is

[TABLE]

According to the known Gauss–Markov theorem, the least squares estimate (LSE)

[TABLE]

has the smallest variance of all the unbiased object’s estimates (Lawson and Hanson [1974]). It follows from (8) and (16) that

[TABLE]

Equations (17) define the principal components of LSE $p_{*}$. Unlike the object’s principal components $\{p_{0k}\}$, the LSE components $\{p_{*k}\}$ are random variables. One can easily find the mean value and the covariance matrix of the LSE:

[TABLE]

where the matrix

[TABLE]

Thus, the LSE principal components $\{p_{*k}\}$ are the unbiased estimates of $p_{0k}$, and ${{\rm var}(p_{*k})=\lambda_{k}^{-1}}$. Usually, the ‘tail’ of the sequence $\{\lambda_{k}\}$ is very small, so the variance of corresponding $\{p_{*k}\}$ and consequently the variance of LSE are huge.

Let us remind the geometrical interpretation of this phenomenon. With the help of (8) and (17), it is easy to transform the definition (15) into the form

[TABLE]

Therefore, the feasible region consists of hollow ellipsoids, centered at the LSE, and the shape of ellipsoids is defined by the ${n\times n}$ matrix $I$. The latter is a representation of the Fisher matrix with the components

[TABLE]

for the particular inverse problem (5) under consideration (Terebizh [1995a,b]). The lengths of semi-axes of the FR ellipsoid are determined by expressions $\ell_{k}=\sqrt{t_{1-\alpha}^{(n)}/\lambda_{k}}$. Small values of the farthest eigenvalues $\{\lambda_{k}\}$ in the spectrum of matrix $I$ give rise to an extremely elongated shape for the FR. Just that phenomenon reveals itself in the well-known instability of inverse solutions. Indeed, a trial object’s estimate $x$ that situated very far from the true object $x_{0}$ can produce the image $y$ that is in feasible agreement, in a scale of natural noise fluctuations, with the really observed image $y_{0}$.

The feasible region usually does not include the LSE and the manifold in its vicinity. Specifying the said in the Introduction, the reason is that the object’s estimates close to LSE try to ‘explain’ all details of the observed image, irrespective of their statistical significance. Since the model (5) supposes essential smoothing of the object, one should admit large erroneous oscillations in the LSE in order to fit tiny random fluctuations in the image. The formal base of the corresponding requirement is given by defining two significance levels $(\alpha_{1},\alpha_{2})$, as it is usually done in mathematical statistics.

4. Optimal linear filter

It is possible to mitigate the harmful influence of the small eigenvalues of the Fisher matrix by introducing into (17) the appropriate set of weights ${w=[w_{1},\ldots,w_{n}]^{T}}$, so

[TABLE]

A number of known inverse solutions, in particular, Kolmogorov [1941] and Wiener [1942] optimal estimate, the regularized solution by Phillips [1962] and Tikhonov [1963], and the truncated estimate (Varah [1973], Hansen [1987, 1993], Press et. al. [1992]), belong to the class of linearly filtered estimates. It follows from (10) and (22) that the squared error of the filtered estimate

[TABLE]

As one can see, the error is minimized by the set of weights

[TABLE]

which constitutes the optimal Wiener filter $\tilde{W}(p_{0})={\rm diag}\left[\tilde{w}(p_{0})\right]$. Consequently, the best of linearly filtered estimates of the object is

[TABLE]

An important feature of the optimal filter is that the weights $\tilde{w}$ depend not only on the known properties of the PSF and the noise but also upon the object itself. For that reason, the filter can be applied only in the Bayesian approach to inverse problems. It is worth noting, in this connection, that the investigations of Kolmogorov [1941] and Wiener [1942] focused on time series analysis, where the Bayesian approach is well justified since the Gaussian nature of ensembles is ensured by the central limit theorem. For most other inverse problems, and in particular, image restoration, the availability of both object ensembles and prior probability distributions on those ensembles is unnatural.

We can simplify the general description of the FR for the linearly filtered estimates by substituting ${p_{w}=Wp_{*}}$ into (15) or (20). The result is:

[TABLE]

This condition imposes restrictions on the system of weights $w$. Then (22) enables the filtered estimate to be found.

One can expect that the requirement (26) is satisfied for the optimal filter ${W=\tilde{W}(p_{0})}$ at moderate values of the significance level $\alpha$. Indeed, extensive numerical simulations are in agreement with this assumption; the corresponding significance level usually is more than $0.70$.

5. Quasi-optimal filter

If it were possible to find a good approximation of the object’s principal components $\{p_{0k}\}$ in (24) with only the given and the observed quantities, the corresponding filter would doubtless have a practical value, but we have no a priori information for such immediate approximation. At the same time, and that is the key point of the quasi-optimal filtering, we have enough information about the structure of the optimal estimate $\tilde{x}_{w}$, in order to require similar properties for the estimate of the object searched for.

By substituting $x_{0}$ from (10) and $\tilde{x}_{w}$ from (25) into (23), and noting that the orthogonal transform does not change the vector norm, we obtain

[TABLE]

This equation simply gives another representation of the error of the optimal filter, which, by definition, is smallest in the class of linear filters.

Let us now consider a trial estimate $p$ close to $p_{0}$ (Fig. 1). Taking into account (27), we shall require that the filter

[TABLE]

which is based on such an estimate, had the minimal error:

[TABLE]

Note that the quasi-optimal filter (28) has the same structure as the optimal Wiener filter (24). Thus, we search for the estimate that most closely simulates behaviour of the best inverse solution.

If we depart from the averaging procedure, which is executable only in theory, and add the condition (26), which requires that the trial object’s estimate belongs to the feasible region, we obtain the simultaneous conditions

[TABLE]

The solution $p_{\min}$ of this system allows us to find the quasi-optimal estimates of the object and its principal components:

[TABLE]

Indeed, we are ultimately interested not in the $p_{\min}$ that is intended to replace $p_{0}$ only in argument of the filter (see Fig. 1), but in the filtered estimate of the principal components $\tilde{p}$, which is analogous to the optimal Wiener estimate $\tilde{p}_{w}$ in (25).

In the components of the corresponding vectors, equations (30) can be written as

[TABLE]

where $\tilde{w}_{k}(p)$ are given by (28) and the vector $\phi=\{\phi_{k}\}$ was defined by (11).

Unlike the Wiener filter, the quasi-optimal filter is nonlinear with respect to the LSE $p_{*}$, because a solution $p_{\min}$ of the system (30) is dependent upon $p_{*}$, and then we should apply filtering according to (31).

Since both functionals in equations (30) are positive definite, and the second functional is non-degenerate, the solution of the constrained minimization problem (30) is unique (Press et. al. [1992], § 18.4).

To understand the sense of the quasi-optimal filtering better, it is useful to bear in mind the following. The object and its least squares estimate were held fixed when searching for the optimal filter, whereas the filter structure has been optimized. On the contrary, equations (28) and (29) fix the previously determined structure of the filter (and the LSE, of course), concentrating attention on the search for an appropriate estimate of the object. Such an approach seems to be quite justified, because the simultaneous searches for both the best filter and the good inverse solution are possible only if complete information about the object is available. The efficiency of the optimal filtering should be high enough in the vicinity of the unknown object; so we do have reason to fix a form of the best filter for an estimate close to the object.

6. Model cases

Equations (28), (30) and (31) form the basis for an algorithm that can be programmed with a high-level programming language. To show the distinction between two filters under discussion more clearly, we deliberately consider here simple examples.

Figure 2 describes restoration of a low-frequency object that we have assumed to be the portion of a sinusoid having an amplitude $1000$. A space-invariant PSF

[TABLE]

was adopted, where ${{\rm sinc}(t)\equiv\sin(\pi t)/(\pi t)}$, and the characteristic radius $R$ was taken as $9$ pixels. Function (32) can be considered as the one-dimensional analogue of the Airy diffraction pattern. The mean level of the Gaussian white noise $a$ was taken as zero, its standard deviation $\sigma_{\xi}$ as $100$. The significance levels of the filters were equal each other.

As one can see from Figs. 2b and 2c, the quasi-optimal weight function and the principal components are practically coincide with the corresponding optimal values at low spatial frequencies. The same is true for the restored objects; both estimates are indistinguishable in the scale of Fig. 2d. Note the removal of the erroneous high-frequency oscillations in the object estimates, and the non-monotonic behaviour of both the optimal and the quasi-optimal weights, which is distinct from those for a truncated estimate. The latter leaves in the object’s estimate simply a few of the first principal components; the quasi-optimal filter leaves only those principal components that have the highest accuracy of restoration. The errors of the Wiener and the quasi-optimal filters were nearly the same for the considered example.

Figure 3 depicts a traditionally difficult model case that incorporates superposition of the sharp and smooth details. The Gaussian PSF has been applied this time with the standard deviation ${\sigma_{PSF}=3}$ pixels; the noise standard deviation has remained as above. As one can see from Figure 3, both the optimal and the quasi-optimal estimates have similar qualities.

We should not assert the claims to quality of restoration in the last case, because even the theoretically best filter, the Wiener’s one, shows rather inexpressive results in that case. Is is more important that the numerous model cases testify that the Wiener’s and the quasi-optimal filters provide very close results. So to speak, a rogue is able to cure any disease, whereas a true physician is only able to do what is possible under the circumstances.

The discussion of the non-negativity condition and the Poisson model is given elsewhere (Terebizh [2003, 2005]).

7. Concluding remarks

It is appropriate to emphasize the importance of the Fisher matrix (21) that plays a fundamental role not only in the linear model but also in the general inverse problem (Terebizh [1995a,b]). To simplify the discussion, we assumed above that the spectrum of the matrix $I$ can include the arbitrarily small, but strictly non-zero eigenvalues $\{\lambda_{k}\}$. This restriction is not essential for the final results; the case of some zero eigenvalues can be treated with the aid of the known additional procedure with the LSE (Press et. al. [1992]).

To avoid misunderstanding, let us repeat once more that the Bayesian way of estimation, in itself, is irreproachable. Moreover, since it incorporates an additional a priory information about the object, the quality of the Bayesian estimation is higher of those in the classical approach. Insurmountable contradictions arise only at the reference to the so called Bayes’s hypothesis, i.e., at substitution really available a priory information concerning the object by some speculative general principles.

Perhaps, the reasons are better visible at the analysis of the concrete procedures. Examples concerning the maximum entropy method were considered by Terebizh [2005], §3.3; we shall touch here the regularization procedure that is described more elaborately in Appendix 1.

In the Phillips–Tikhonov’s approach, the most widely used stabilizing requirement consists in condition of minimal ‘power’ of the inverse solution, namely

[TABLE]

given of course the acceptable misfit of the observed and the trial images. We have found out above an important role of the principal components ${p=\{p_{k}\}}$ that are associated with the estimate $x$ by the orthogonal transform ${p=V^{T}x}$. This role is based mainly on the mutual statistical independence of principal components, which allows to apply the theoretically most effective way of extracting information from an unstable least squares estimate: it is necessary to take principal components one after another according to descending order of the eigenvalues $\{\lambda_{k}\}$ of the Fisher matrix $I$. As the $I$ spectrum covers an extremely wide range, the estimates of the principal components have the essentially different accuracy. According to equations (18) and (19), ${{\rm var}(p_{*k})=\lambda_{k}^{-1}}$; for common in practice situations that means tens of orders in value. Meanwhile, the invariance of a vector norm at the orthogonal transformation entails the equivalence of the condition (33) and the requirement

[TABLE]

Evidently, the direct summation the variables of essentially different accuracy is not an optimal procedure. For example, it seems better to take $p_{k}^{2}$ with some weights $g(\lambda_{k})$ that are dependent on the corresponding variance of the least squares estimate. This way leads to the condition

[TABLE]

which, however, assumes some subjective choice of the weight function.

The objectively justified requirement is given in (30) by the close on sense, but more refined condition ${F(p)=\min}$, where

[TABLE]

Taking into account (28), we can rewrite the above expression as

[TABLE]

Usually, the production $\lambda_{k}p_{k}^{2}\gg 1$ for the small $k$ and then quickly decreases with the growth of $k$ to the values much less than 1. Therefore, the summation range in (37) can be approximately divided into two regions with the boundary value of $K$, such that ${\lambda_{K}p_{K}^{2}\simeq 1}$. Then

[TABLE]

We see that the sought estimates $p_{k}$ and the LSE principal components $p_{*k}$ should be close only for the large first eigenvalues $\lambda_{k}$ of the Fisher matrix. Such a requirement is quite reasonable in view of high accuracy of the first principal components of LSE. At the same time, just the summary ‘tail’ of $\{p_{k}^{2}\}$ of low accuracy is minimized when $\lambda_{k}p_{k}^{2}\ll 1$.

On the contrary, the condition (34) equally minimizes all principal components irrelatively of their accuracy. As the numerous model cases show, that entails too large systematic shift, the bias of the regularized inverse solution. In geometrical language it means that the point of a contact the ellipsoidal feasible region and the ‘stabilizing’ region noticeably depends on the shape of the latter region.

From the viewpoint of the regularization theory, it might seem that the functional (36) is a smoothing functional similar to $\|x\|^{2}$ or to one of the several forms of the ‘entropy’ ${\cal E}(x)$. Indeed, the condition ${F(p)=\min}$ promotes stabilization of the inverse solution, but the origin of this functional is of vital importance. The Bayesian hypothesis proposes to compensate the lack of a priori information by some general principle that directly concerns the properties of the sought object $x$ itself. Obviously, it is possible to offer an unlimited number of such principles. On the contrary, we rely on the intrinsic reserves of the inverse theory. The relatively much weaker assumption has been applied in the proposed above way, namely, that the optimal Wiener filter retains a high efficiency in the local vicinity of the unknown object. It appears that, instead of the prior information about the object, it is enough to lay down some reasonable statistical requirements only to the restoration procedure.

Appendix 1. Statistical treatment of the Phillips–Tikhonov

regularization

The statistical point of view on the Phillips–Tikhonov procedure is illustrated below by the example when regularizing functional $\Phi(x)$ is taken in the form (4). As it was stated in the main text, the corresponding feasible region (FR) is defined by equation (15). Combining this condition with (4) and taking into account that the norm of the vector ${x=Vp}$  remains the same under orthogonal transformation, we come to the simultaneous system

[TABLE]

where vector $\phi$ and the diagonal matrix $\Delta$ are defined by equations (11) and (9), respectively. The first of the above conditions describes an ellipsoidal FR which elements, by definition, provide the feasible misfit between the observed and trial images. Let us consider elements situated on the sphere ${\|p\,\|^{2}={\rm const}}$ of small radius centered at the origin. Gradually enlarging the radius of the sphere, we take into consideration elements of the larger power, so the point of contact of the spheres family and the fixed FR ellipsoid gives the Phillips–Tikhonov’s inverse estimate.

According to the method of Lagrange multipliers, the necessary extremum conditions follow from minimization of the auxiliary functional

[TABLE]

where the scalar $\gamma\geq 0$ as yet is free. If the vector $p_{\gamma}$ that provides minimum of the Lagrange function is found for any $\gamma$, then the desired value of the regularization parameter $\gamma$ is defined by substitution $p_{\gamma}$ into the first of conditions (A1.1).

In order to find an explicit form of $p_{\gamma}$, let us rewrite (A1.2) as

[TABLE]

where

[TABLE]

The diagonal matrix $\bar{\Delta}$ is of size ${n\times n}$, and $\bar{\phi}$ is the $n$-vector. The minimum of the functional (A1.3) gives the element

[TABLE]

where $\Lambda=\Delta^{2}$, and $p_{*}=\Delta^{-1}\phi$ is the least square estimate (LSE) of principal components (see the first of equations (A1.1)). By defining the diagonal ${n\times n}$ matrix

[TABLE]

we can write down the regularized vector of the principal components and the corresponding inverse solution as

[TABLE]

Comparing (A1.6) and (A1.7) with the general definition (22) of a filtered estimate, we come to conclusion that the regularized according to Phillips and Tikhonov inverse solution belongs to the class of estimates that were obtained by linear filtering of the least squares (maximum likelihood) estimate.

Remind that the definition of FR for a linear filtered estimate is given by (26), or, in an unfolded form,

[TABLE]

Substituting here the weights according to (A1.6), we come to the following equation for the inverse value of the regularization parameter ${\mu\equiv 1/\gamma}$:

[TABLE]

where the function

[TABLE]

According to (11) and (6), $\phi=U^{T}z_{0}=U^{T}C^{-1/2}(y_{0}-a)$, where $y_{0}$ is the observed image.

As one can see, function $f(\mu)$ monotonously descends from ${f(0)=\|\phi\|^{2}}$ down to zero when ${\mu\to\infty}$, so any of standard numerical methods can be easily applied to find the unique root of equation (A1.9), i.e., the regularization parameter $\gamma$ (see, e.g., Booth [1955], Press et al. [1992]). After that the formulae (A1.6) and (A1.7) allow us to find the statistical estimate of the object of given significance level.

The above algorithm is intended for the restoration of any objects, both positive and with some negative components. If a priory information about the object $x_{0}$ assumes non-negativity of the all its components, we have, instead of (A1.1), the following simultaneous system for estimating the vector ${x=Vp}$:

[TABLE]

It was shown by Terebizh [2005], §3.2 that the system can be reduced to the known constrained least squares problem (Lawson and Hanson [1974], Golub and Van Loan [1989], Kahaner et al. [1989]), which efficient solution is embedded to the powerful computing environment MatLab and some other modern systems.