text:
# On Optimal Detection of Point Sources in CMB Maps

## 1 Introduction

The separation of different physical components is an important issue in the analysis of Cosmic Microwave Background (CMB) data. Among the foreground components, point-sources deserve especial attention given their extreme non-Gaussianity and highly variable spectral index. A brief summary of different methods can be found in Sanz et al. (2001) (henceforth SHM).

Methods used to detect point sources should act as high-pass filters to detect the high frequency structure introduced by the sources in the CMB data while masking other foreground contamination (like dust, synchrotron and free-free emission) characterized by lower frequencies. In principle, the good space/frequency characteristics of wavelets should make these functions an attractive choice for such tasks (e.g., Cayón et al. 2000; Vielva et al. 2001a, b, and references therein). Indeed, wavelets have been proved optimal for detection of point like singularities – at least for one-dimensional signals. Point sources in CMB maps, however, are not point singularities because the signal is smoothed by the beam of the antenna; sources are expected to have the shape of the antenna’s pattern. The question is then how to include this beam profile information in the wavelet analysis. SHM considered this question and, were lead to define optimal pseudo-filters for source detection in CMB maps. However, we show that nothing seems to be gained with these wavelet filters and that other simpler techniques lead to better source detection methods.

We first set up the basic framework. Although the detection of point-like sources in CMB maps is a two-dimensional problem, we present our arguments in $`R^n`$, as in SHM, because the same methods may be used in other applications.

The sources are assumed to be point-like signals convolved with the beam of the measuring instrument and are thus assumed to have a known profile $`\tau (𝒙)`$. The signal $`y(𝒙)`$, $`𝒙R^n`$, is modeled as

$$y(𝒙)=\underset{j}{}s_j(𝒙)+z(𝒙)$$
(1)
where

$$s_j(𝒙)=A_j\tau (𝒙𝒙_j),$$
(2)
$`A_j`$ and $`𝒙_j`$ are, respectively, unknown source amplitudes and locations, and $`z(𝒙)`$ is a zero-mean background with power-spectrum $`P(𝒒)`$

$$\mathrm{E}[z(𝒒)z^{}(𝒒^{})]=P(𝒒)\delta ^n(𝒒𝒒^{}).$$
(3)
Henceforth $`\mathrm{E}[]`$ and “$`^{}`$” will denote the expectation and complex conjugate operators, respectively, $`\delta ^n(𝒒𝒒^{})`$ the $`n`$-dimensional Dirac distribution, and $`z(𝒒)`$ the Fourier transform of $`z(𝒙)`$

$$z(𝒒)=_{\mathrm{}}^+\mathrm{}z(𝒙)\mathrm{e}^{i𝒒𝒙}𝒅𝒙.$$
(4)
To properly remove the point sources from the signal we need to estimate the locations $`\{𝒙_j\}`$ and amplitudes (fluxes) $`\{A_j\}`$ of the sources.

A classical method used to estimate source locations is based on identifying peaks in the cross-correlation function

$$c(𝒙)=_{\mathrm{}}^+\mathrm{}y(𝒙+𝒃)\tau (𝒃)𝒅𝒃.$$
(5)
The rationale is that $`c(𝒙_𝒐)`$ measures the similarity between the source profile with a section of $`y(𝒙)`$ centered at $`𝒙_𝒐`$; a peak in $`c(𝒙_𝒐)`$ is an indication of a source signal at $`𝒙_𝒐`$. Eq. (5) is a filtering of the signal $`y(𝒙)`$ with a filter $`\tau (𝒃)`$ that amplifies the characteristic frequencies of the source. Once the sources have been located, their amplitudes can be estimated by means of classical fitting procedures like least squares.

The cross-correlation function (5) does not take into account the background characteristics. This is a great disadvantage in cases where the power spectrum $`P(𝒒)`$ is known or a good estimate is available. In Sects. 2 and 3 we consider other methods that take into account this information and that may be considered extensions of the cross-correlation method.

The basic procedure we consider is as follows. The signal is first filtered to enhance the sources with respect to the background. This is done by cross-correlating the signal $`y(𝒙)`$ with a filter $`\mathrm{\Phi }`$ as in (5) (with $`\mathrm{\Phi }`$ in place of $`\tau `$). The source locations are then determined by selecting the peaks in the filtered signal that are above a selected threshold. Finally, the source amplitudes are estimated with the values of the filtered signal at the estimated locations. The question we consider first is the selection of an optimal filter $`\mathrm{\Phi }`$ for such procedure.

## 2 Designing an optimal filter

The optimality criteria we use are based on the following assumptions (for futher justification of these assumptions see SHM). The source profile and background spectrum are known. The profile is spherically symmetric, characterized by a scale $`R_s`$, and the background is assumed to be isotropic. In this case, we can write $`s(𝒙)s(x)`$, where $`x=𝒙`$, and $`P(𝒒)P(q)`$ for $`q=𝒒`$. In addition, source overlap is assumed negligible.

We consider the general family of spherically symmetric filters $`\mathrm{\Phi }(𝒙;𝒃)`$ of the form

$$\mathrm{\Phi }(𝒙;𝒃)=\varphi (𝒙𝒃)$$
(6)
with Fourier transform $`\varphi (q)`$. The filtered field is

$$\begin{array}{cc}\hfill w(𝒃;\varphi )& =_{\mathrm{}}^+\mathrm{}y(𝒙)\mathrm{\Phi }(𝒙;𝒃)𝒅𝒙\hfill \\ & =_{\mathrm{}}^+\mathrm{}y(𝒒)\varphi (q)\mathrm{e}^{i𝒒𝒃}𝒅𝒒,\hfill \end{array}$$
(7)
where $`y(𝒒)`$ and $`\varphi (q)`$ are, respectively, the Fourier transforms of $`y(𝒙)`$ and $`\varphi (𝒙)`$. The mean and variance of $`w(𝒃;\varphi )`$ are

$$\begin{array}{cc}\hfill \mu (𝒃;\varphi )& =\mathrm{E}w(𝒃;\varphi )\hfill \\ & =\alpha _0^+\mathrm{}q^{n1}s(q)\varphi (q)\mathrm{e}^{i𝒒𝒃}𝑑q;\hfill \end{array}$$
(8)
$$\begin{array}{cc}\hfill \sigma ^2(\varphi )& =\mathrm{E}w^2(𝒃,\varphi )\mu (𝒃;\varphi )^2\hfill \\ & =\alpha _0^+\mathrm{}q^{n1}P(q)\varphi (q)𝑑q,\hfill \end{array}$$
(9)
where $`\alpha =2\pi ^{n/2}\mathrm{\Gamma }^1(n/2)`$.

The first constraint on the filter concerns the second stage of the procedure; the source locations are assumed known and the objective is to estimate the amplitudes. Given the assumed distance between sources, it is enough to consider a field $`y(𝒙)`$ as in (1) with a single source at the origin, $`s(𝒙)=A\tau (𝒙)`$. To estimate its amplitude we ask that $`w(\mathrm{𝟎};\varphi )`$ be an unbiased estimator of $`A`$ – i.e., $`\mu (\mathrm{𝟎};\varphi )=A`$ – so that $`\varphi `$ is required to satisfy the equation

$$_0^+\mathrm{}q^{n1}\tau (q)\varphi (q)𝑑q=\frac{1}{\alpha }.$$
(10)
To enhance the magnitude of the source relative to the background we determine the filter $`\mathrm{\Phi }`$ that minimizes the variance $`\sigma ^2(\varphi )`$. This has the effect of maximizing, among unbiased estimators, the detection level

$$𝒟(\varphi )=\frac{\mu (\mathrm{𝟎};\varphi )}{\sigma (\varphi )},$$
(11)
which measures the capability of the filter to correctly detect a source at the prescribed location (see SHM).

Since $`\mathrm{\Phi }`$ is chosen so that $`w(\mathrm{𝟎};\varphi )`$ is a minimum variance linear – in $`y(𝒙)`$ – unbiased estimator of $`A`$, it follows that (Gauss-Markov theorem) $`w(\mathrm{𝟎};\varphi )`$ is the (generalized) least squares estimate of $`A`$ achieved by the filter

$$\varphi (q)=\frac{1}{\alpha a}\frac{\tau (q)}{P(q)},a_0^+\mathrm{}q^{n1}\frac{\tau ^2}{P}𝑑q,$$
(12)
with minimum variance

$$\sigma ^2(\varphi )=\frac{1}{\alpha a}.$$
(13)
Filter (12) is a particular case of a well known class of filters, known as matched filters in the engineering literature, that are designed to optimize signal-to-noise ratio (e.g., Kozma & Kelley 1965; Pratt 1991). The arguments in this section also show that the filter introduced by Tegmark & Oliveira-Costa (1998) is a particular case of the matched filter (12) with $`P(q)`$ representing the background power spectrum before the smoothing of the antenna.

For white noise, $`P(q)=\mathrm{const}=D`$, filter (12) simplifies to

$$\varphi (q)=\frac{1}{\alpha aD}\tau (q)$$
(14)
which, up to a constant factor, is identical to the filter used in the classical cross-correlation function (5). This provides a justification for the use of the cross-correlation for source detection in white noise.

Having recognized $`w(\mathrm{𝟎};\varphi )`$ as a least squares estimator of $`A`$, we close this section with some remarks from least squares methodology that we consider relevant. First note that, regardless of the spectrum $`P(q)`$, the source profile $`\tau (q)`$ properly normalized, that is $`\varphi (q)=\tau (q)/K`$ for $`K=\alpha q^{n1}\tau (q)\varphi (q)\mathrm{e}^{i𝒒𝒃}𝑑q`$, also provides an unbiased estimator $`w(\mathrm{𝟎};\varphi )`$ of $`A`$. This is the (ordinary) least squares estimate that does not take into account the covariance of the background; it is unbiased but not minimum variance. However, it is well known that when the covariance is actually estimated from the data, the ordinary least squares estimate may be better than the generalized one (e.g., Draper & Smith 1998). In other words, uncertainities in the spectrum estimates may lead to worse amplitude estimates than those obtained with the simpler cross-correlation filter. Uncertainties in the spectrum will also affect the selection of a detection threshold.

Note also that the unbiasesness of $`w(\mathrm{𝟎};\varphi )`$ as an estimator of $`A`$ depends on knowing the correct source location, it is not necessarily unbiased once the source locations are estimated from the data. This is shown in Sect. 4.2.

## 3 Pseudo-filters

In the pseudo-filter approach of SHM the filters are of the form (6) with an additional scale dependence

$$\mathrm{\Psi }(𝒙;R,𝒃)=\frac{1}{R^N}\psi \left(\frac{𝒙𝒃}{R}\right),$$
(15)
for some spherically symmetric function $`\psi `$. From (15) we get a set of wavelets of different scales $`R`$ and shifts $`𝒃`$ provided the regularity conditions $`\tau (q)/P(q)0`$, $`(1/P)(d\tau /d\mathrm{ln}q)0`$ for $`q0`$ are satisfied. The filtered field $`w(𝒃,R;\psi )`$ at scale $`R`$ is defined as in (7) but with $`\psi (Rq)`$ in place of $`\varphi (q)`$.

To determine an optimal filter $`\psi `$, SHM minimize the variance of the filtered field subject to two constraints: First, $`w(\mathrm{𝟎},R_0;\psi )`$ is required to be, as in the previous section, an unbiased estimator of $`A`$ for some known $`R_0R_s`$. For the second constraint $`\psi `$ is selected so that $`\mu (\mathrm{𝟎},R;\psi )`$ has a local maximum at scale $`R_0`$. This constraint translates to

$$_0^+\mathrm{}q^{n1}\tau (q)\psi (R_0q)\left(n+\frac{d\mathrm{ln}\tau }{d\mathrm{ln}q}\right)𝑑q=0.$$
(16)
Minimizing $`\sigma ^2(R_0;\psi )`$ with the two constraints yields the filter (SHM)

$$\psi (R_0q)=\frac{1}{\alpha \mathrm{\Delta }}\frac{\tau (q)}{P(q)}\left[nb+c(na+b)\frac{d\mathrm{ln}\tau (q)}{d\mathrm{ln}q}\right],$$
(17)
where $`\mathrm{\Delta }=acb^2`$,

$`b`$ $`{\displaystyle _0^+\mathrm{}}q^{n1}{\displaystyle \frac{\tau }{P}}{\displaystyle \frac{d\tau }{d\mathrm{ln}q}},`$ (18)
$`c`$ $`{\displaystyle _0^+\mathrm{}}q^{n1}{\displaystyle \frac{1}{P}}\left({\displaystyle \frac{d\tau }{d\mathrm{ln}q}}\right)^2𝑑q,`$
and $`a`$ is as in (12). This filter provides a field of variance

$$\sigma ^2(R_0;\psi )=\frac{n^2a+2nb+c}{\alpha \mathrm{\Delta }}.$$
(19)
The estimator of the amplitude $`A`$ obtained with this filter is again linear and unbiased and therefore, by the optimality of least squares, $`\sigma ^2(R_0;\psi )\sigma ^2(\varphi )`$ regardless of the source profile and background spectrum. This means that the detection level of $`\mathrm{\Phi }`$ is at least as high, or higher, than that achieved with $`\mathrm{\Psi }`$.

### 3.1 Is $`\mathrm{\Psi }`$ optimal for source detection?

We have seen that the second constraint (16) does not increase the detection level when the source locations are known. But this is not surprising since the constraint is defined to take advantage of the known source scale to help determine source locations. We will show that even when source location uncertainty is taken into account, (16) does not improve on the simpler filter $`\mathrm{\Phi }`$ based on the single contraint (10). In other words, enough information about the scale of the source is already included in the derivation of the matched filter. Therefore, it seems that nothing is gained with wavelet methods.

In principle, as explained by SHM, constraint (16) can be used to test if a detection corresponds to a true source by checking for its maximum at scale $`R_0`$. But finding a spike at the correct scale is not enough to make sure it is not a spurious noise artifact, we have to follow it across scales to make sure that it scales appropriately. In any case, a similar test can also be designed for the matched filter: given that the source profile is assumed known, the matched filter leads to an amplitude versus scale dependence (see example in Sect. 4.1) that can be determined and used for such a test. This functional dependence provides more information than the existence of a local maxima and should, in the presence of noise, lead to better location estimates.

## 4 Examples

To compare the theoretical performances of the filters $`\mathrm{\Phi }`$ and $`\mathrm{\Psi }`$ we use the gain as defined by SHM

$$\begin{array}{cc}\hfill g(\psi ,\varphi )& =𝒟(\varphi )/𝒟(\psi ,R_0)\hfill \\ & =\sigma (R_0;\psi )/\sigma (\varphi ),\hfill \end{array}$$
(20)
where $`𝒟(\psi ,R_0)`$ is as in (11) but defined for $`\mathrm{\Psi }(𝒙;R_0,\mathrm{𝟎})`$. We start with a specific example with Gaussian profiles and power-law spectra.

### 4.1 Gaussian sources and $`P(q)=Dq^\gamma `$

Gaussian sources

$$\tau (q)=\theta ^n\mathrm{e}^{(q\theta )^2/\mathrm{\hspace{0.17em}2}},$$
(21)
where $`\theta `$ is the “standard deviation” defining the scale, provide an important family of source profiles. Indeed, in many practical applications the instrument’s profile, and thus the point sources, are effectively characterized by Gaussian profiles. For the noise process we take the power-law spectrum $`P(q)=Dq^\gamma `$. This family of spectra can be used to locally approximate spectra of other homogeneous processes and has as special cases white ($`\gamma =0`$) and $`1/f`$ ($`\gamma =1`$) noise processes.

For a Gaussian profile and a power-law spectrum, Eq. (12) leads to

$$\varphi (q)=\frac{\mathrm{\Gamma }(n/2)}{\mathrm{\Gamma }(m)}\frac{(q\theta )^\gamma }{\pi ^{n/2}}\mathrm{e}^{(q\theta )^2/\mathrm{\hspace{0.17em}2}},$$
(22)
where $`m=(n+\gamma )/2`$. ¿From (13) and (19) we obtain the following variance and gain

$`\sigma ^2(\varphi )`$ $`={\displaystyle \frac{D}{\theta ^{n\gamma }\pi ^{n/2}}}{\displaystyle \frac{\mathrm{\Gamma }(n/2)}{\mathrm{\Gamma }(m)}},`$ (23)
$`g(\psi ,\varphi )`$ $`=\left[\mathrm{\hspace{0.17em}1}+{\displaystyle \frac{(n\gamma )^2}{4m}}\right]^{1/2}.`$
We see that $`g1`$, which shows that $`\varphi `$ has a higher detection level than $`\psi `$. For $`n=\gamma `$ —for example, one-dimensional process with $`1/f`$ noise— the two filters lead to the same detection levels. This is expected since $`\varphi (q)=\psi (Rq)`$ for $`n=\gamma `$ (compare Eq. (22) with Eq. (24) in SHM). In other words, the second constraint is just redundant in this case. Note also that $`\varphi (q)`$ is the Mexican hat wavelet for $`n=2`$ and $`\gamma =2`$. This results show that, contrary to what has been claimed before (e.g., Cayón et al. 2000), the optimality of the Mexican hat wavelet does depend on the background spectrum.

### 4.2 A numerical experiment

We have presented theoretical arguments showing that $`\varphi (q)`$ has a better source detection capability than $`\psi (Rq)`$ when the source location is known. We now confirm that $`\varphi `$ is still better when the uncertainity of source location estimates is taken into account. To compare with the results in SHM, it is enough to consider a simple example with one-dimensional Gaussian sources and white noise ($`\gamma =0`$) (we already know that $`\varphi =\psi `$ for $`n=1`$ and $`1/f`$ noise). In this case (14) and (23) become, respectively,

$`\varphi (q)`$ $`={\displaystyle \frac{1}{\sqrt{\pi }}}\mathrm{e}^{(q\theta )^2/\mathrm{\hspace{0.17em}2}},`$ (24)
$`g(\psi ,\varphi )`$ $`=\left({\displaystyle \frac{3}{2}}\right)^{1/2}>1.`$
That is, the detection level of $`\varphi `$ is about 20% larger than that of $`\psi `$. The amplitude dependence on $`\theta `$ of a Gaussian source of scale $`R_s`$ filtered with $`\varphi `$ is

$$A(\theta )=\frac{\sqrt{2}R_s}{(R_s^2+\theta ^2)^{1/2}}A.$$
(25)
A fit to this dependence can be used to determine if a detection corresponds to a source of the correct scale, just as large wavelet coefficients would be tracked across different wavelet scales.

Fig. 1 shows the filters $`\mathrm{\Phi }(𝒙)`$ and $`\mathrm{\Psi }(𝒙;R_0,\mathrm{𝟎})`$ and their corresponding Fourier transforms. It shows that the two filters are quite different. For example, to provide filtered sources with the scale $`R_0`$, $`\psi `$ has to pass higher frequencies than $`\varphi `$. This can be a problem for signals contaminated by high frequency noise.

Table 1 shows the average number of correct and incorrect detections obtained with $`\mathrm{\Psi }`$ and $`\mathrm{\Phi }`$ and a fixed $`3\sigma `$ threshold for signal-to-noise (S/N) ratios equal to 1, 2 and 3. The two filters give equivalent results for higher S/N sources. We see that $`\varphi `$ leads to a higher number of correct detections and a lower number of incorrect ones. But a proper comparison should take into account that the filters require different thresholds. The second row for $`\mathrm{\Psi }`$ shows the corresponding results when the threshold is chosen to lead to the same average number of rejections as $`\mathrm{\Phi }`$. For low S/N sources $`\mathrm{\Phi }`$ leads again to a higher number of correct detections while for larger S/N they give similar results.

To compare amplitude estimates that include location uncertainty, we take the average of the amplitudes (since all the generated sources have the same amplitude) of all detections. The results are shown in Table 1. The errors in the amplitudes are of the order of 0.5% or less. We see that amplitude estimates are biased when the source locations are estimated, and that the bias is larger for $`\mathrm{\Psi }`$. For low S/N the amplitude is overestimated because high peaks are easier to detect and noise peaks are incorrectly classified as sources. For high S/N sources we also have centering problems but this time the smaller amplitude in the noise peaks leads to underestimated source amplitudes. We can draw similar conclusions from the results of two-dimensional simulations shown in Tables 2 and 3;

To conclude, note that the spatial support of the filters is an important factor when the assumption of nonoverlapping filters is invalid. The support must be small compared to the distance between the sources. Figs. 2-3 show (for $`n=1,2`$) the cumulative energies

$$\mathrm{E}_\psi (x)=\frac{_0^x\mathrm{\Psi }^2(b;R_0,0)𝑑b}{_0^{\mathrm{}}\mathrm{\Psi }^2(b;R_0,0)𝑑b}$$
(26)
and

$$\mathrm{E}_\mathrm{\Phi }(x)=\frac{_0^x\mathrm{\Phi }^2(b)𝑑b}{_0^{\mathrm{}}\mathrm{\Phi }^2(b)𝑑b},$$
(27)
as a function of $`\gamma `$ for a background spectrum $`P(q)=Dq^\gamma `$. These functions measure the “energy concentration” of the filters and provide information about their “spatial” support. The figures show that also in this respect $`\mathrm{\Phi }(x)`$ has similar, if not better, characteristics than $`\mathrm{\Psi }(x;R_0,0)`$. In particular, the filter $`\mathrm{\Phi }`$ has a tighter support for faster decaying noise spectra.

#### 4.2.1 Selecting the detection level

We briefly justify our selection of $`3\sigma `$ level. This threshold should be chosen large enough to reduce the number of false detections but small enough not to miss too many sources. To properly choose a threshold we have to understand the statistics of local maxima above a chosen level. For a general homogeneous (Gaussian) random field this is a difficult question (some asymptotic results can be found in Adler 1981) but simulations can be carried out when the background power spectrum is known. Our simulations showed that the traditional $`5\sigma `$ level is too conservative for the signal lengths used in the examples. If $`P_M`$ is the proportion of local maxima above $`3\sigma `$ for a field without sources. We found that the probability that $`P_M`$ is higher than 0.001 is about $`9\%`$ for the signal filtered with $`\mathrm{\Psi }`$ and less than $`10^5`$ for the signal filtered with $`\mathrm{\Phi }`$.

## 5 Summary and Conclusions

We have revisited the problem of estimating point sources of known profile in an isotropic background. The methods we considered are based on two basic interrelated procedures: source detection by thresholding of local maxima, and amplitude estimation by linear filtering. We have compared the effects of different constraints on the selection of an optimal filter.

The first constraint is typical in matched filter methodology where S/N is maximized. The optimal filter provides unbiased least squares estimates of source amplitudes at known source locations. By the optimality of least squares, these amplitude estimates can not be improved with any other unbiased linear filter. However, uncertainities in source locations introduce a bias in amplitude estimates. Amplitudes are overestimated at low S/N and underestimated at high S/N. A second constraint introduced by SHM is designed to improve estimates of source locations by maximizing amplitudes at the correct scale. We found that this constraint does not lead to better estimates as compared to those obtained with a simpler matched filter, especially for lower S/N sources. For high S/N sources, as it was observed by Cayón et al. (2000), the results of the two methods are the same. These results contradict previous optimality studies of wavelet-based filters for detection of point-sources of known profile in an isotropic CMB background of known spectrum.


format: text

text:
# A theorem for the beam splitter entangler

## Abstract

It is conjectured that the an entanglement output states from a beam splitter requires the non-classicality in the input state(M.S. Kim, W. Son, V. Buzek and P. L. Knight, Phys. Rev. A, 65, 032323(2002)). Here we give a proof for this conjecture.

The beam splitter is one of the few quantum devices that may act as the entangler. The entangler properties of a beam splitter have been studied in the past. In particular, Kim et al studied the entangler property with many different input states such as the Fock number state, the coherent state, the squeezed state and the mixed states in Gaussian form. It has been conjectured there that, to obtain an entangled output state, a necessary condition is that the input state should be non-classical. Unfortunately, there is no proof for this conjecture. In this paper, we give a very simple proof for this conjecture.
Consider a loseless beam splitter(see figure 1). We can distinguish the field mode $`a`$ and mode $`b`$ by the different propagating direction. Most generally, the property of a beam splitter operator $`\widehat{B}`$ in Schrodinger picture can be summarized by the following equations(see e.g., ref)

$`\rho _{out}=\widehat{B}\rho _{in}\widehat{B}^1,`$ (1)
$`\widehat{B}^{}=\widehat{B}^1,`$ (2)
$`\widehat{B}\left(\begin{array}{c}\widehat{a}\\ \widehat{b}\end{array}\right)\widehat{B}^1=M_B\left(\begin{array}{c}\widehat{a}\\ \widehat{b}\end{array}\right),`$ (7)
$`M_B=\left(\begin{array}{cc}\mathrm{cos}\theta e^{i\varphi _0}& \mathrm{sin}\theta e^{i\varphi _1}\\ \mathrm{sin}\theta e^{i\varphi _1}& \mathrm{cos}\theta e^{i\varphi _0}\end{array}\right)`$ (10)
$`\widehat{B}|00=|00.`$ (11)
Here $`\rho _{in}`$ and $`\rho _{out}`$ are the density operator for the input and output states respectively. Both of them are two mode states including mode $`a`$ and mode $`b`$. The elements in the matrix $`M_B`$ are determined by the beam splitter itself, $`\widehat{a},\widehat{b}`$ are the annihilation operators for mode $`a`$ and mode $`b`$ respective, $`|00`$ is the vacuum state for both mode. Equation (11) is due to the simple fact of no input no output. Without any loss of generality, we can express $`\rho _{in}`$ by the $`P`$ representation in the following form:

$`\rho _{in}={\displaystyle _{\mathrm{}}^{\mathrm{}}}P(\alpha _a,\alpha _b,\alpha _a^{},\alpha _b^{})|\alpha _a,\alpha _b\alpha _a,\alpha _b|\mathrm{d}^2\alpha _a\mathrm{d}^2\alpha _b,`$ (12)
where $`|\alpha _a,\alpha _b`$ is the coherent state in two mode Fock space, i.e.

$`|\alpha _a,\alpha _b=\widehat{D}_{ab}(\alpha _a,\alpha _b)|00`$ (13)
and

$`\widehat{D}_{ab}(\alpha _a,\alpha _b)=e^{\widehat{a}^{}\alpha _a\widehat{a}\alpha _a^{}+\widehat{b}^{}\alpha _b\widehat{b}\alpha _b^{}}.`$ (14)
If $`\rho _{in}`$ is a classical state, the distribution function $`P(\alpha _a,\alpha _b,\alpha _a^{},\alpha _b^{})`$ must be non-negative definite in the whole complex planes. In such a case, the ouput state is

$`\rho _{out}={\displaystyle _{\mathrm{}}^{\mathrm{}}}P(\alpha _a,\alpha _b,\alpha _a^{},\alpha _b^{})\widehat{B}|\alpha _a,\alpha _b\alpha _a,\alpha _b|\widehat{B}^1\mathrm{d}^2\alpha _a\mathrm{d}^2\alpha _b`$ (15)
which is equivalent to

$`\rho _{out}={\displaystyle _{\mathrm{}}^{\mathrm{}}}P(\alpha _a,\alpha _b,\alpha _a^{},\alpha _b^{})\widehat{B}\widehat{D}_{ab}(\alpha _a,\alpha _b)\widehat{B}^1\widehat{B}|0000|\widehat{B}^1\widehat{B}\widehat{D}_{ab}\widehat{B}^{}.`$ (16)
From equation(11) we know that $`\widehat{B}|0000|\widehat{B}^1=|0000|`$. By equation(7) we can see that

$`\widehat{B}\widehat{D}_{ab}(\alpha _a,\alpha _b)\widehat{B}^1=\widehat{D}_{ab}(\alpha _a^{},\alpha _b^{})`$ (17)
and

$`(\alpha _a^{},\alpha _b^{})=(\alpha _a,\alpha _b)M_B.`$ (18)
In short, the following equatoin can be easily obtained by equation(7,10,11 )

$`\widehat{B}|\alpha _a,\alpha _b\alpha _a,\alpha _b|\widehat{B}^1=|\alpha _a^{},\alpha _b^{}\alpha _a^{},\alpha _b^{}|`$ (19)
Since $`\mathrm{det}M_B=1`$, we have the following formula for the output state

$`\rho _{out}={\displaystyle _{\mathrm{}}^{\mathrm{}}}P(\alpha _a,\alpha _b,\alpha _a^{},\alpha _b^{})|\alpha _a^{},\alpha _b^{}\alpha _a^{},\alpha _b^{}|\mathrm{d}^2\alpha _a^{}\mathrm{d}^2\alpha _b^{}.`$ (20)
This is equivalent to

$`\rho _{out}={\displaystyle _{\mathrm{}}^{\mathrm{}}}P^{}(\alpha _a,\alpha _b,\alpha _a^{},\alpha _b^{})|\alpha _a,\alpha _b\alpha _a,\alpha _b|\mathrm{d}^2\alpha _a\mathrm{d}^2\alpha _b`$ (21)
and

$`P^{}(\alpha _a,\alpha _b,\alpha _a^{},\alpha _b^{})=P(\alpha _a^{\prime \prime },\alpha _b^{\prime \prime },\alpha _a^{\prime \prime },\alpha _b^{\prime \prime })`$ (22)
and

$`(\alpha _a^{\prime \prime },\alpha _b^{\prime \prime })=(\alpha _a^{\prime \prime },\alpha _b^{\prime \prime })M_B^1.`$ (23)
Since $`P(\alpha _a,\alpha _b,\alpha _a^{},\alpha _b^{})`$ is non-negative, functional $`P^{}(\alpha _a,\alpha _b,\alpha _a^{},\alpha _b^{})`$ must be also non-negative. By the definition of separability, state $`\rho _{out}`$ defined by equation(21) must be separable. Therefore we have the following theorem:
Theorem: If the input state is a classical state, the output state of a beam splitter must be a separable state.
This is equivalent to say, in order to obtain an entangled output state, the non-classicality of the input state is a necessary condition. This theorem can be extended to a more general situation in the multimode Fock space. Let’s consider the rotation operator $`\widehat{R}`$ in n-mode Fock space. We have

$`\widehat{R}\mathrm{\Lambda }\widehat{R}^1=M_R\mathrm{\Lambda },`$ (24)
where $`\mathrm{\Lambda }=(\widehat{c}_1,\widehat{c}_2\mathrm{}\widehat{c}_n)^T`$, $`\widehat{c}_i`$ are the annihilation operators of the $`i^{}`$th mode and $`M_R`$ is a $`n`$dimensional unitary matrix. By using the BCH formula

$`e^\mu \nu e^\mu =\nu +[\mu ,\nu ]+{\displaystyle \frac{1}{2!}}[\mu ,[\mu ,\nu ]]+\mathrm{}`$ (25)
we have the following explicitly formula for the operator $`\widehat{R}`$:

$`\widehat{R}=\mathrm{exp}\left(\mathrm{\Lambda }^{}\mathrm{ln}M_R\mathrm{\Lambda }\right).`$ (26)
Therefore we know that

$`\widehat{R}|00\mathrm{}0=|00\mathrm{}0.`$ (27)
Any classical multimode state in Fock space can be written in the following probabilistic distribution

$`\rho ={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝐏(\alpha ,\alpha ^{})|\alpha \alpha |\mathrm{d}^2\alpha ,`$ (28)
where $`|\alpha =|\alpha _\mathrm{𝟏}\alpha _\mathrm{𝟐}\mathrm{}\alpha _𝐧`$ and $`𝐏(\alpha ,\alpha ^{})`$ is a non-negative functional provided that $`\rho `$ is a classical state. Similar to the two mode case, we can shown that,

$`\rho ={\displaystyle _{\mathrm{}}^{\mathrm{}}}𝐏^{}(\alpha ,\alpha ^{})|\alpha \alpha |\mathrm{d}^2\alpha ,`$ (29)
and

$`𝐏^{}(\alpha ,\alpha ^{})=𝐏(\alpha ^{\prime \prime },\alpha ^{{}_{}{}^{\prime \prime }})`$ (30)
and

$`(\alpha ^{\prime \prime })=(\alpha )M_R^1.`$ (31)
Obviously, the funtional $`𝐏^{}(\alpha ,\alpha ^{})`$ is non-negative. Thus we draw the following conclusion in the multimode Fock space: A classical density operator in the multimode Fock space is separable under arbitrary rotation.
Although the nonclassicality in the input state is a necessary condition, it is obviously not a sufficient condition for the entanglement in the output state of a beam splitter. Since a beam splitter operater is unitary, it is reversible. It has been shown in ref that nonclassical separable input state can be changed to an entangled state in the output. The inverse of such a process makes examples that even though the input state is nonclassical, the output could be still separable. Some specific examples are given in .
I thank Prof Imai for support. I thank Dr Huang WY, Dr Winter, Dr Yura H, Dr Matsumoto K,and Dr Tomita A for useful discussions.


text:
# Naive mean field approximation for image restoration

## I Introduction

In this paper, we investigate the image restoration problem. This problem involves synthesis of an image from a corrupted image with a model of the information available (or assumed) on the source image and the corruption process Geman84  Pryce95  Geiger91 . The problem lends itself naturally to a Bayesian formulation, which estimates a probability (posterior) for the original image on the basis of the model probability (priors) of the assumed model of the source image and corruption process.

One strategy in image restoration is to use the Bayesian inference by adopting the image that maximizes the posterior probability. This method is called the maximum a posteriori (posterior) probability (MAP) inference. Given a corrupted input image, the MAP inference accepts the image that maximizes the posterior probability as the restored result. The logarithm of posterior probability can be regarded as the energy, so we can consider the MAP inference as an energy minimization problem. Therefore, the image restoration problem can be regarded as an optimization problem. Geman & Geman demonstrated that the MAP inference result can be applied to the image restoration problem by using simulated annealing, which is a tool for searching the ground state Geman84 .

Another strategy is the inference in which the expectation value with respect to the maximized marginal posterior probability at each site in thermal equilibrium is regarded as the original image. This method is called maximizer of the posterior marginals (MPM) inference Marroquin87  Rujan93  Sourlas94 . In the MAP inference, the posterior probability is given for each set of pixels. In contrast, in the MPM inference, the posterior marginal probability controlled by the temperature $`T`$ is given for each pixel value. To find a restoration image by the MPM inference, we should calculate thermal average for each pixel value. The MPM inference includes the MAP inference Nishimori00  NishimoriBook01  because , at the limit of the temperature $`T0`$, the MPM inference becomes equivalent to the MAP inference. In general, however, the MPM inference is evaluated at $`T>0`$, so this method is also called “finite temperature restoration”, and the MPM inference has an advantage over the MAP inference in the restoration ability of each pixel Nishimori00  NishimoriBook01 .

Recently, the MPM inference has been discussed in the field of error correcting code Rujan93  Sourlas94  Nishimori00 . The problem of error correcting code is similar to the image restoration problem in the sense that the received bit sequence, which corresponds to the image represented by a set of pixels, is corrupted by noise, and the receiver tries to retrieve the original bit sequence/image from the noisy one. The major difference between error correcting code and image restoration is that the image restoration problem usually gives only the corrupted image and not other additional redundant information. Therefore, in image restoration, we usually assume alternative information for the original image such as the prior probability. In the field of error correcting code, Rujan has pointed out that carrying out the decoding procedure not at the ground state but at a finite temperature is effective Rujan93 . Sourlas used the Bayes formula to re-derive the finite-temperature decoding of Rujan’s result under more general conditions Sourlas94 .

Assuming the number of pixels as $`N`$ and each pixel as a binary unit that can take $`\{1,+1\}`$, the number of feasible combinations of pixel values becomes $`2^N`$. Because of the existence of local minima, finding the ground state task is very difficult with a method such as the gradient decent. Thus, to apply the MAP or MPM inference, we usually use a relaxation method as typified by annealing technique. Moreover, the MPM inference involves a more important problem as described below. In the MAP inference, once the system reaches the macroscopic equilibrium state, each pixel value can be properly determined with its probability as $`1`$. On the other hand, in the MPM inference, when the system reaches the macroscopic equilibrium state, each pixel value cannot be determined uniquely because the probabilities of each binary state have some finite values in the finite temperature decoding. Therefore, we should calculate thermal averages for each pixel, and this requires many samplings.

In this study, we discuss an approximation that replaces the stochastic dynamics of the MPM inference with deterministic dynamics. In statistical mechanics, this approximation is called the “naive mean field equation (NMFE)” Bray86 . The NMFE has usually been applied for emulating the behavior of the stochastic unit. By Applying the NMFE, each restoring unit, which is defined as a stochastic binary unit, is replaced by a deterministic analog unit that takes $`[1,+1]`$. One important advantage of applying NMFE is the ability to reduce calculation cost by eliminating the need to calculate thermal average , which requires many samplings. The NMFE has been applied to several combinatorial optimization problems such as the traveling salesman problem (TSP) Hopfield86 . However, almost all researchers who have used NMFE have not pursued quantitative issues such as the accuracy of the approximation.

In §II, we formulate the image restoration problem using the Bayes inference in the manner of Nishimori & Wong’s formulationNishimori00 . Recently, from the statistical-mechanical point of view, Nishimori & Wong have formulated the image restoration problem by introducing a mean field model for binary image restoration. They analyzed the model theoretically by the replica method. We applied NMFE to the formulation and analyzed it by the replica method in the manner of Bray et al. Bray86 .

In §III, we compare the results between our analysis and computer simulations. Within the limits of the mean field approximation, our analysis showed agreement with the computer simulations. However, the prior probability derived from the mean field approximation, which is called an infinite range model, is usually not practical for a real image prior. Therefore, we discuss the difference between infinite range model and its origin, called the nearest 0 interaction model.

## II Model and Analysis

### II.1 Formulation of the Image restoration

In this section, we apply the mean field theory to the image restoration problem in the manner of Nishimori & Wong. Figure 1 shows a schematic diagram of image restoration. The original image is represented by $`\{\xi _i\}`$, and each unit is a binary unit that takes two states $`\{1,+1\}`$. The number of units corresponding to the image size is denoted by $`N`$. In addition, the original image $`\{\xi _i\}`$ is assumed to have the following prior probabilities:

$`P_s\left(\left\{\xi _i\right\}\right)`$ $`={\displaystyle \frac{1}{Z(\beta _s)}}\mathrm{exp}({\displaystyle \frac{\beta _s}{N}}{\displaystyle \underset{i<j}{}}\xi _i\xi _j),`$ (1)
$`Z(\beta _s)`$ $`=\mathrm{Tr}_\xi \mathrm{exp}({\displaystyle \frac{\beta _s}{N}}{\displaystyle \underset{i<j}{}}\xi _i\xi _j),`$ (2)
where $`\beta _s>0`$. The operator $`\mathrm{Tr}_\xi `$ means trace, a sum over all possible $`2^N`$ states of $`\xi `$. $`Z(\beta _s)`$ indicates the partition function of this prior probability. Assuming $`\beta _s>0`$, this probability suggests that pixel values $`\xi _i`$ and $`\xi _j`$ have a tendency to take the same value.

In image restoration, the interaction between each pixel should be described by a local rule, such as a sum of the nearest neighbors. For analysis by the mean field theory, however, we replace this local interaction rule with a global one, that means each pixel has interactions to all other pixels.

In the conventional image restoration framework, the transmitting signal is only the raw code $`\xi `$ (upper path in Fig.1). However, considering the similarity of the formulation of the error-correcting code problem using a Bayesian inference Rujan93  Sourlas94  Nishimori00 , it is natural to introduce transmission of a redundant code for better image restoration. Nishimori & Wong proposed the introduction of a redundant code called the “parity check code” (lower path in Fig. 1) as well as the raw code. For comparison with their results, in this study we also adopted transmission of the parity check code. For the parity check code, we adopted a 2-body interaction term denoted by $`\xi _i\xi _j/N`$. The reason for dividing $`\xi _i\xi _j`$ by $`N`$ is to apply the mean field theory. At the receiver side, the transmitting signal $`\xi _i`$ corresponds to $`\tau _i`$, and the parity check code $`\xi _i\xi _j/N`$ corresponds to $`J_{ij}`$. The posterior probability based on observation signal $`P(\{\xi _i\}|\{J_{ij}\},\{\tau _i\})`$ can be represented by the Bayes formula:

$$P(\{\xi _i\}|\{J_{ij}\},\{\tau _i\})P_{\mathrm{out}}(\{J_{ij}\},\{\tau _i\}|\{\xi _i\})P_s(\{\xi _i\}).$$
(3)
When the receiver can guess the form of the original image posterior probability, the restoring posterior probability can be assumed by replacing the original pixel value $`\xi _i`$ with the estimated pixel value $`\sigma _i`$. In general, however, the receiver could not guess the original image prior probability accurately. In order to evaluate the performance of the restoration ability by the difference between the original image prior and the restoring prior, we assume that the restoring image prior probability $`P_m()`$ has a different parameter $`\beta _m`$ from the original image prior $`P_s()`$.

$`P_m\left(\left\{\sigma _i\right\}\right)`$ $`={\displaystyle \frac{1}{Z(\beta _m)}}\mathrm{exp}({\displaystyle \frac{\beta _m}{N}}{\displaystyle \underset{i<j}{}}\sigma _i\sigma _j),`$ (4)
$`Z(\beta _m)`$ $`=\mathrm{Tr}_\sigma \mathrm{exp}({\displaystyle \frac{\beta _m}{N}}{\displaystyle \underset{i<j}{}}\sigma _i\sigma _j),`$ (5)
where $`\beta _m>0`$. Each form of prior is identical mathematically except for the parameters of interaction strength. We primarily discuss the $`\beta _s\beta _m`$ case.

$`P_{\mathrm{out}}(\{J_{ij}\},\{\tau _i\}|\{\xi _i\})`$ is a conditional probability of the observation signal, and it is a probability expression of the corrupting process in the transmission channel. In this study, we assume that the noise added in the corrupting process is similar to Gauss distribution as follows:

$`P_{\mathrm{out}}(\{J_{ij}\},\{\tau _i\}|\{\xi _i\})`$
$`\mathrm{exp}\left({\displaystyle \frac{N\left(J_{ij}\frac{J_0}{N}\xi _i\xi _j\right)^2}{2J^2}}{\displaystyle \frac{\left(\tau _i\tau _0\xi _i\right)^2}{2\tau ^2}}\right)`$ (6)
The first term in the exponential indicates that the noise corresponds to the parity check code, and the second term corresponds to the raw code. For the first term in the exponential, the random variables $`J_{ij}`$ follow the normal distribution, whose average is $`J_0\xi _i\xi _j/N`$, and the variance is $`J^2/N`$. In the second term in the exponential, the random variables $`\tau _i`$ follow the normal distribution whose average is $`\tau _0\xi _i`$, and the variance is $`\tau ^2`$.

Nishimori & Wong discussed the macroscopic characteristic of this system by using statistical mechanics Nishimori00 . They introduced a system that has the following Hamiltonian:

$$H=\beta \underset{i<j}{}J_{ij}\sigma _i\sigma _j\frac{\beta _m}{N}\underset{i<j}{}\sigma _i\sigma _jh\underset{i}{}\tau _i\sigma _i.$$
(7)
They calculated the free energy by the replica method for averaging probabilities $`P_s\left(\left\{\xi _i\right\}\right)`$ and $`P_{\mathrm{out}}(\{J_{ij}\},\{\tau _i\}|\{\xi _i\})`$. Ignoring the 2-body interaction term, that is $`\beta =0`$, this Hamiltonian consists of the prior probability term and the observed data fitting term. This form is equivalent to the conventional cost function form.

In the MAP inference, the estimation value of restored pixels $`\sigma `$ results in the minimum of the Hamiltonian defined in eq. (7). On the other hand, the MPM inference, which corresponds to the finite temperature estimation, estimates pixel values by the $`\mathrm{sgn}(\sigma _i)`$, where $`\sigma _i`$ means the thermal average. Considering the limit of the temperature at $`0`$, the MPM inference becomes equivalent to the MAP inference.

We introduced the quantity overlap $`M_o`$ to evaluate the restoration ability measure for the MPM inference.

$$M_o=\frac{1}{N}\underset{i}{}\xi _i\mathrm{sgn}(\sigma _i).$$
(8)
Nishimori & Wong suggested that the minimization of the Hamiltonian (7) and the maximization of the overlap are not equivalent Nishimori00 . Moreover they also indicated that the MPM inference (finite temperature estimation) has better restoration ability than the MAP inference in the sense of the overlap:

$$\begin{array}{cc}\hfill M_o& =\mathrm{Tr}_{\xi ,J,\tau }P_s\left(\left\{\xi _i\right\}\right)P_{\mathrm{out}}(\{J_{ij}\},\{\tau _i\}|\{\xi _i\})\hfill \\ & \times \xi _i\mathrm{sgn}(\sigma _i).\hfill \end{array}$$
(9)
This overlap can be evaluated by averaging probability variables with the statistical-mechanical technique.

### II.2 Naive mean field approximation

In the finite temperature decoding which corresponds to the MPM inference, we should calculate the thermal average $`\sigma _i`$, which means the estimation value of $`\sigma _i`$. We evaluated the calculation cost of MPM inference by computer simulation. As we will show that calculating this average requires about 50 times longer computational time than that for achieving to the macroscopic equilibrium state. This result will be given in §III.

In this study, to find the ground state, we introduce the analog neural network model for approximation proposed by Hopfield & Tank Hopfield86 . By this approximation, each Ising unit is replaced by an analog unit that can take continuous value $`[1,+1]`$, and the output of each analog unit can be regarded as the thermal average of the corresponding Ising spin unit, which can take two states $`\sigma _i=\pm 1`$. For replacing Ising unit to analog one, we introduced a Hamiltonian as substitution for eq. (7) :

$$=M\left[\beta \underset{i<j}{}J_{ij}\widehat{s}_i\widehat{s}_j\frac{\beta _m}{N}\underset{i<j}{}\widehat{s}_i\widehat{s}_jh\underset{i}{}\tau _i\widehat{s}_i\right],$$
(10)
where the unit $`\widehat{s}_i`$ is an analog unit. $`M`$ is a scaling factor described as below. Following the manner of Bray et alBray86 , we assumed each $`i`$th site consists of $`M`$ Ising units, and the analog unit output $`\widehat{s}_i`$ is calculated by the average of $`M`$ Ising units $`\widehat{\sigma }_{ia}`$ (see Fig. 2).

$$\widehat{s}_i=\frac{1}{M}\underset{a}{\overset{M}{}}\widehat{\sigma }_{ia}.$$
(11)
Considering the limit $`M\mathrm{}`$, each output $`\widehat{s}_i`$ can take an analog value $`[1,+1]`$. When $`M`$ is a finite value, the analog unit output $`\widehat{s}_i`$ is called ’binominal spins’ which can take $`1,1+\frac{2}{M},\mathrm{},1\frac{2}{M},1`$ with binominal distribution. So that, we can define ‘spin weight function’ as

$`W(\widehat{s}_i)`$ $`=M\mathrm{Tr}\delta \left(M\widehat{s}_i{\displaystyle \underset{a=1}{\overset{M}{}}}\widehat{\sigma }_{ia}\right)`$
$`={\displaystyle \frac{M}{2\pi \mathrm{j}}}{\displaystyle _\mathrm{j}\mathrm{}^{+\mathrm{j}\mathrm{}}}𝑑u_i\mathrm{exp}(M(u_i\widehat{s}_i+\mathrm{ln}2\mathrm{cosh}(u_i))).`$ (12)
The partition function $`Z`$ can be described as

$`Z`$ $`={\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle _1^{+1}}𝑑\widehat{s}_iW(\widehat{s}_i)\mathrm{exp}()`$
$`=M{\displaystyle _\mathrm{j}\mathrm{}^{+\mathrm{j}\mathrm{}}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left({\displaystyle \frac{du_i}{2\pi j}}\right){\displaystyle _1^{+1}}{\displaystyle \underset{i=1}{}}d\widehat{s}_i`$
$`\mathrm{exp}(M(\beta {\displaystyle \underset{i<j}{}}J_{ij}\widehat{s}_i\widehat{s}_j+{\displaystyle \frac{\beta _m}{N}}\widehat{s}_i\widehat{s}_j+h{\displaystyle \underset{i}{}}\tau _i\widehat{s}_i`$
$`u_i\widehat{s}_i+\mathrm{ln}2\mathrm{cosh}(u_i)))`$ (13)
In the limit $`M\mathrm{}`$, the integrals over $`\{u_i\}`$ and $`\{\widehat{s}_i\}`$ in eq. (13) could be evaluated by the saddle-point method. The stationary equations are

$`0`$ $`=\beta {\displaystyle \underset{i<j}{}}J_{ij}\widehat{s}_i\widehat{s}_j+{\displaystyle \frac{\beta _m}{N}}{\displaystyle \underset{i<j}{}}\widehat{s}_i\widehat{s}_j+h{\displaystyle \underset{i}{}}\tau _i\widehat{s}_iu_i,`$ (14)
$`0`$ $`=\widehat{s}_i+\mathrm{tanh}u_i.`$ (15)
Eliminating $`u_i`$, we can obtain

$$\widehat{s}_i=\mathrm{tanh}(\beta \underset{j}{}J_{ij}\widehat{s}_j+\frac{\beta _m}{N}\underset{j}{}\widehat{s}_j+h\tau _i).$$
(16)
From the manner of Hopfield & Tank Hopfield86 , we can also consider a discrete synchronous updating rule:

$$s_i^{t+1}=\mathrm{tanh}(\beta \underset{j}{}J_{ij}s_j^t+\frac{\beta _m}{N}\underset{j}{}s_j^t+h\tau _i),$$
(17)
where $`s_i^t`$ denotes the analog unit output at time $`t`$. When the system described by eq. (17) reached to the equilibrium state $`s_i^{\mathrm{}}`$, whole units should satisfy eq. (16). Therefore, to investigate the equilibrium state of dynamics (17), we should take the analog Hamiltonian described by eq. (10) correctly. In the equilibrium state, each analog unit state expressed by $`s_i^{\mathrm{}}`$ corresponds to $`\sigma _i`$, i.e. $`s_i^{\mathrm{}}`$ can be regarded as the thermal average of $`\sigma _i`$.

This analog unit replacement is sometimes called the “naive mean field equation (NMFE) approximation”. From eq. (17), this system follows deterministic dynamics, so that it is no need to calculate the thermal average of stochastic unit; therefore, the calculation cost is lower than the stochastic Ising model.

However this replacement from stochastic Ising unit $`\sigma _i`$ to the deterministic analog unit $`s_i`$ is an merely substitution. Thus, we should analytically investigate the result of applying NMFE by analog model to compare it with the estimated thermal average of the Ising model. Hence, to investigate the equilibrium state of the system described by eq. (16), we should analyze the system that has the Hamiltonian described by eq. (10).

#### II.2.1 Replica method

We now analyze the system that has the Hamiltonian described by eq. (10) by the “replica method”, which is a standard statistical analysis tool. The $`n`$ replicated partition function can be expressed as:

$$\begin{array}{cc}\hfill [Z^n]& =Tr\left[\underset{i<j}{}dJ_{ij}\sqrt{\frac{N}{2\pi J^2}}\right]\left[\underset{i}{}\frac{d\tau _i}{\sqrt{2\pi }\tau }\right]\hfill \\ & \times P_s\left(\left\{\xi _i\right\}\right)P_{\mathrm{out}}(\{J_{ij}\},\{\tau _i\}|\{\xi _i\})\hfill \\ & \times \mathrm{exp}(\frac{\beta }{M}_{i<j,\alpha ,a,b}J_{ij}\widehat{\sigma }_{ia}^\alpha \widehat{\sigma }_{jb}^\alpha \hfill \\ & +\frac{\beta _m}{MN}\underset{i<j,\alpha ,a,b}{}\widehat{\sigma }_{ia}^\alpha \widehat{\sigma }_{jb}^\alpha +\underset{i,\alpha ,a}{}h\tau _i\widehat{\sigma }_{ia}^\alpha ),\hfill \end{array}$$
(18)
where operator $`\mathrm{Tr}`$ means the sum over all states about $`\{\widehat{\sigma }_{ia}^\alpha \}`$ and $`\{\xi _i\}`$. We analyzed this replicated partition function by the standard replica method. The replica symmetry solution can be described as:

$`m_0`$ $`=\mathrm{tanh}(\beta _sm_0),`$ (19)
$`m`$ $`={\displaystyle \frac{Tr_\xi e^{\beta _sm_0\xi }Dx\widehat{F}(U(x))}{2\mathrm{cosh}(\beta _sm_0)}},`$ (20)
$`t`$ $`={\displaystyle \frac{Tr_\xi e^{\beta _sm_0\xi }\xi Dx\widehat{F}(U(x))}{2\mathrm{cosh}(\beta _sm_0)}},`$ (21)
$`q`$ $`={\displaystyle \frac{Tr_\xi e^{\beta _sm_0\xi }Dx\widehat{F}(U(x))^2}{2\mathrm{cosh}(\beta _sm_0)}},`$ (22)
$`\chi `$ $`={\displaystyle \frac{1}{\sqrt{h^2\tau ^2+\beta ^2J^2q}}}{\displaystyle \frac{Tr_\xi e^{\beta _sm_0\xi }Dxx\widehat{F}(U(x))}{2\mathrm{cosh}(\beta _sm_0)}},`$ (23)
where $`U()`$ means

$$U(x)=\sqrt{h^2\tau ^2+\beta ^2J^2q}x+\beta _mm+(h\tau _0+\beta J_0t)\xi ,$$
(24)
and the function $`\widehat{F}()`$ is a solution of the self consistent equation:

$`F(x)`$ $`=`$ $`\mathrm{tanh}(U(x)+\beta ^2J^2\chi F(x)).`$ (25)
Using these solutions, we could obtain the overlap $`M_o`$:

$$\begin{array}{cc}\hfill M_o& =\frac{1}{N}\underset{i}{}\xi _i\text{sgn}(s_i)\hfill \\ & =\frac{Tr_\xi e^{\beta _sm_0\xi }Dx\xi \mathrm{sgn}(U(x))}{2\mathrm{cosh}(\beta _sm_0)}.\hfill \end{array}$$
(26)
## III Results

In this section, we compare the analysis results between applying the NMFE approximation (analog model) and estimated thermal average of the Ising model. We also compare these theoretical analyses with the computer simulation results.

In the following subsection from III.1 to III.3, we treated infinite range Ising/analog model. In the infinite range model, the original image prior probability is defined by eq. (1), and the image is corrupted by the noisy channel according to eq. (6). In the infinite range Ising model, we assumed restoring prior probability as eq. (4), and the Hamiltonian could be derived as eq. (7). On the other hand, in the infinite range analog model, we defined the Hamiltonian as eq. (10).

However, for image prior, the infinite range interaction is just a little strange assumption. Thus, we introduced models whose prior have only nearest neighbor interaction in III.4. We call them as Ising model with nearest neighbor interactions and analog model with nearest neighbor interactions respectively. Unfortunately, it is difficult to treat the nearest neighbor interaction models analytically. Therefore, we calculated the result of nearest neighbor interaction models by computer simulation and compared them with the result of the infinite range models. In the nearest neighbor interaction models, we assume that the corruption process in transmitting is same as infinite range model, that is defined by eq. (6).

In the Ising model Monte Carlo simulation, we use the asynchronous Glauber dynamics, i.e. we selected one site $`\sigma _i`$ randomly, and calculate probability as

$$P(\sigma _i=\pm 1)=\frac{(1\pm \mathrm{tanh}(\beta _jJ_{ij}\sigma _j+\frac{\beta _m}{N}_j\sigma _j+h\tau _i))}{2}.$$

According to the probability the selected site $`\sigma _i`$ is decided.

On the other hand, for analog simulation, we used discrete synchronous update rule described as eq. (17). i.e. whole units are updated simultaneously. We consider each unit value in the equilibrium state as its thermal average.

### III.1 Without parity check code

First, we discuss the case where no parity check code is used ($`\beta =0`$). This case corresponds to the conventional image restoration problem. In the analysis, when $`\beta =0`$, the order parameter equations obtained in the previous section are identical to the result of the Ising analysis suggested by Nishimori & Wong Nishimori00 . Thus, the MPM inference by the stochastic infinite range Ising model and deterministic infinite range analog model are equivalent. Figure 4 illustrates the restoring ability of finite temperature decoding using the infinite range Ising stochastic units. Figure 4 shows the restoring ability using infinite range analog model. In each figure, the vertical axis means the overlap $`M_o`$, and the horizontal axis means the restoring temperature $`T_m`$, which is corresponds to the restoring image prior coefficients $`\beta _m`$ ($`T_m=\beta _m^1`$). The solid line means the analysis result by the replica method, and each error bar shows the quartile deviation obtained by the computer simulation.

We set the parameter $`T_s`$, which is the reciprocal of the original image prior coefficients $`\beta _s`$ ($`T_s=\beta _s^1`$), at $`T_s=0.9`$. At the limit of $`T_m0`$, the MPM inference becomes equivalent to the MAP inference. In both figures 4 and 4, the restoration ability around $`T_m=(T_s)=0.9`$ is better than that of $`T_m0`$. This means that the MPM inference has better image restoration ability rather than the MAP inference as far as the overlap.

### III.2 Convergence speed

The biggest advantage of the analog model over the Ising model is its lower calculation cost. In this section, we discuss the calculation time that should be set in the computer simulation for each method.

We observed the convergence time of the overlap in each simulation. Figure 5 shows the convergence speed. In the infinite range Ising model, we adopt the synchronous update and define $`N_m`$ updates as one Monte-Carlo step (1 MCS). $`N_m`$ means the number of units in the simulation. A synchronous analog update corresponds to 1 MCS in the sense of the number of update units. The horizontal axis of Fig. 5 indicates the calculation time measured by MCS. The vertical axis means the average of overlap at that time. Each model achieved the macroscopic equilibrium state in 20 MCS. The analog model adopts a deterministic algorithm, we can calculate $`s_i^{\mathrm{}}`$ at that time, and the $`M_o`$ converges at 20 MCS. However, assuming an Ising stochastic algorithm, we should calculate the thermal average of $`s_i`$ by sampling. Figure 5 shows that needs about 1000 MCS to converge. In this result, the deterministic infinite range analog model converged 50 times faster than the stochastic infinite range Ising model in finite temperature decoding.

### III.3 Using parity check code

In the replica method analysis, the infinite range analog model and the infinite range Ising model are equivalent in the case using $`\beta =0`$. This is derived from the fact that the susceptibility $`\chi `$ remained at $`0`$ under $`\beta =0`$. Using the parity check code for better restoration, the analysis results differ. In the case of $`\beta >0`$, $`\chi `$ becomes greater than $`0`$, and the term corresponding to $`\chi `$ in eq. (25) is effective. This term is called the Onsager reaction term in statistical mechanics. As a result, a difference in analysis between the Ising model and the infinite range analog model appears.

When the noise variance of the parity check code is small (e.g., $`J0.6`$), the infinite range analog model and the infinite range Ising model are almost the same in the sense of overlap. In this case, little noise exists in the parity check code, and better image restoration is possible. Figure 7 shows a comparison of overlaps obtained by each replica analysis. The horizontal axis indicates the parameter $`\beta `$, and the vertical axis indicates the overlap $`M_o`$. We assumed that the decoding temperature is optimal ($`T_m=T_s(=0.9)`$) and that the parameter is $`h=1`$. Comparing $`\beta `$ with $`h`$, $`\beta `$ becomes large, the restoring system tends to rely on the parity check rather than the raw code. Each restoration result appears the same, and better restoration is achieved than in the case transmitting the raw code only.

Figure 7 shows the case of large noise variance in the parity check code channel ($`J=1.0`$). The horizontal axis signifies $`\beta `$. In this case, each method restores better than when only the raw code channel is used. However, relying on $`\beta `$ too much results in poorer restoration result than when only the raw code channel is used. A difference between the infinite range analog model and the infinite range Ising model appears as restoration worsens. When the noise variance of the parity check code becomes large, the restoration ability of infinite range analog model is not as good as that of the infinite range Ising model.

Figure 9 shows the computer simulation result of the infinite range Ising model under the noise variance parameter $`J=1.0`$, and Fig. 9 shows the computer simulation result of the infinite range analog model. Each simulation matches the analysis curve, and the performance of infinite range analog model simulation is not much better than that of the infinite range Ising simulation. However, when we properly set the ratio of the coefficients $`\beta ,\beta _m`$, and $`h`$, the performance of each is optimized and the peak of ability for each is almost the same. Thus when we can estimate the hyper-parameters properly, there is no difference between these two methods.

### III.4 Ising/analog model with nearest neighbor interaction

In each infinite range model we compared, we assumed that all units have interactions with each other in prior probabilities to use the mean field approximation. However, it is more natural that these probabilities are defined by their nearest-neighbor interactions. We call the model whose encoding/decoding prior probability is defined by the nearest neighbor interaction on the pixel lattice the Ising/analog model with nearest neighbor interaction. It is difficult to analyze the nearest neighbor interaction model with a statistical tool because of its local interactions. Nishimori & Wong compared the analysis result of the Ising model with nearest neighbor interaction solved by computer simulation Nishimori00 . Figure 10 shows the results of our reproduction. The image size is $`400\times 400`$ pixels, and the original image’s prior probability is denoted:

$`P_s\left(\left\{\xi _i\right\}\right)`$ $`=`$ $`{\displaystyle \frac{1}{Z(\beta _s)}}\mathrm{exp}(\beta _s{\displaystyle \underset{ij}{}}\xi _i\xi _j),`$ (27)
$`Z(\beta _s)`$ $`=`$ $`\mathrm{Tr}_𝝃\mathrm{exp}(\beta _s{\displaystyle \underset{ij}{}}\xi _i\xi _j),`$ (28)
where the summation $`_{ij}`$ denotes a sum extending over all pairs of neighboring sites on the pixel lattice.

In this simulation, we adopt $`T_s=1/\beta _s=2.15`$, generate the image $`\{\xi _i\}`$ by the prior probability (27), and corrupt by flipping each pixel by the probability $`p=0.1`$. We use only the raw code channel, not the parity check code channel in this simulation. By using the mean field theory, Nishimori & Wong proved that the overlap $`M_o`$ is maximized despite of the dimensions of the system (image size), when the decoding temperature $`T_m`$ equals the image generating temperature $`T_s`$. In the simulation, we set $`T_s=T_m=2.15`$. Figure 10 shows the reproducing result of the Ising/analog model with nearest neighbor interaction obtained by the computer simulation. With the nearest-neighbor interaction model, it is clear that the overlap of the analog model is smaller than that of the Ising model. We believe that this phenomenon occurred because the analog model does not gather in the effect of fluctuation by other pixel values.

Thus, we should consider the nearest neighbor interaction prior effect in the computer simulation. Figure 11 shows an example of the computer simulation result. Each image has an area of $`400\times 400`$ pixels. Figure 11(a) shows the original image. In this image, noise-like patterns are inherently contained. Figure 11(b) shows the corrupted image. Figure 11(c) shows the restoration result by the Ising model, and Figure 11(d) shows the result by the analog model. The overlap of the Ising model (Figure 11(c)) becomes $`M_o=0.879`$, and the analog model (Figure 11(c)) becomes $`M_o=0.837`$. In addition, when we applied the image restoration by MAP inference, the noise-like patterns contained in the original image (Figure 11(a)) were eliminated by the prior probability term. Thus, to maintain the fine structure of the original image, the MPM inference is an effective method. This characteristics was confirmed in the analog model by the computer simulation.

## IV Discussion

In this research, we discussed the image restoration problem by finite temperature decoding. In the MPM inference proposed by Nishimori & Wong, many trials are required to calculate the thermal average for each unit. We tried to replace this operation with deterministic dynamics called naive mean field approximation. At first, we applied mean field approximation to the infinite range prior model and analyzed this approximated system by the replica method. We also quantitatively compared the restoration ability of the Ising stochastic model with that of the analog deterministic model. As a result, we proved that their abilities are equivalent when we used only the raw code signal, that is, equivalent to conventional image restoration. Assuming a redundant signal, called the parity check code, the result of each method differed slightly. However, when we adopted the ratio of the coefficients properly, that is, the coefficients $`\beta `$, $`h`$, and $`\beta _m`$, the restoring ability of each method showed no difference. This means that if we can estimate these coefficients, called hyper-parameters, we can apply the analog deterministic method without worrying about the restoration ability. Moreover, the analog deterministic method is $`50`$ times faster than the Ising stochastic method.

In the future problem, In analog model, Horiguchi pointed out that the inclination of output function around $`\pm 1`$ cause constitutive difference Horiguchi91 . Thus, we should consider this effect.

Moreover, we will attempt to propose a more effective analog approximation method for the nearest neighbor interaction model.