text:
# Universal amplitude-exponent relation for the Ising model on sphere-like lattices

## 1 Introduction

The observation that the symmetry of systems at a critical point goes beyond scale invariance emerging from the divergence of correlation lengths, and additionally comprises rotational and translational as well as inversional invariance has greatly enlarged the scope of exact results that can be derived for such systems on field-theory grounds. While customary scaling theory and the real-space renormalisation ansatz  give rise to general relations between the scaling exponents and the concept of universality  of exponents and certain amplitude ratios, conformal field theory (CFT)  allows for a classification of models according to their operator contents, thus providing the values of exponents and amplitude ratios, and, in special cases, even the amplitudes themselves. By interpreting the finite size of a system as an additional field in the scaling ansatz, these considerations apply not only to thermal scaling but as well to the scaling of observables of finite systems at a critical point in the limit of diverging system sizes, i.e. to finite-size scaling (FSS) . Concerning universality in FSS, based on renormalisation-group considerations it has been argued that for semi-finite systems below their upper critical dimension the correlation lengths (measured in units of the lattice spacing) themselves are universal quantities . Exploiting conformal invariance this statement can be strengthened and the amplitude of this universal correlation-length scaling can be calculated at least for a special geometry. To see this, consider the logarithmic map

$$w=\frac{L}{2\pi }\mathrm{ln}z,z,$$
(1)
in two dimensions, which wraps the complex plane $``$ around an infinite-length cylinder of circumference $`L`$. Exploiting the facts that the critical, connected two-point correlation function in the plane is fixed by the postulate of conformal invariance and that the given map is conformal, the critical two-point function of any conformally transforming (primary) operator on the cylinder $`S^1\times `$ follows immediately. In particular, correlations of the considered operator along the $``$-axis decline exponentially with a correlation length of

$$\xi _{}=\frac{L}{2\pi x},$$
(2)
where $`x`$ denotes the scaling dimension of the operator. Thus, one arrives at a FSS relation including the amplitude $`A=1/2\pi x`$, implying universality of the correlation length itself.

Since there is an essential difference between the conformal groups in two dimensions, where it happens to be infinite dimensional, and higher dimensions, where it is reduced to a finite-dimensional Lie group, very few CFT results exceeding those already known from renormalisation-group arguments survive a transition to higher dimensions. However, considering the transformation (1) as a simple change of coordinates for a moment, rewriting it in polar coordinates allows for a generalisation to higher dimensions, now mapping $`^d`$ to $`S^{d1}\times `$. Noting that such a map only acts on the radial part of the coordinates, leaving the angular part invariant, Cardy  conjectured another linear amplitude-exponent relation for the correlation lengths in the geometry $`S^{d1}\times `$, namely:

$$\xi _{}=\frac{R}{x},$$
(3)
where $`R`$ denotes the radius of $`S^{d1}`$ and $`x=x(d)`$ is again the scaling dimension. Note that this includes the $`d=2`$ result (2) with $`L=2\pi R`$. This relation, if it can be confirmed, constitutes one of the few exact results for non-trivial three-dimensional systems. Since the notion of primarity of a scaling operator is a priori not well defined for $`d>2`$, however, it is not entirely clear to which operators eq. (3) should apply. To find the analogues of primary operators for $`d>2`$ one might start with the operator product expansion (OPE), cp. the explorative analysis of the operator content of the O($`n`$) models in three dimensions that can be found in ref. . Note that this generalised mapping now connects different geometries, whereas the two-dimensional mapping can be understood as a meromorphism of the Riemann sphere onto itself.

Considering alternative approaches, attempts of numerical investigation of the FSS of systems on such spherical geometries are severely hampered by the fact that there is only a finite number of triangulations, i.e. regular polyhedra, for each sphere $`S^{d1}`$. A numerical transfer matrix study for the case of $`d=3`$ using Platonic solids as a regular discretisation of the sphere was found to be inconclusive due to the restriction in system sizes . Here we follow a different line and use slightly irregular lattices as sphere discretisations for $`d=3`$, exploiting and explicitly probing for the universality of the considered observables. Considering universality, it might be favourable to study ratios of correlation lengths which are known to be quite insensitive to lattice distortions. The attempt to escape such complications by considering the flat geometry $`S^1\times S^1\times \mathrm{}\times `$ instead surprisingly yields similar results, which are, however, up to now not connected to CFT calculations .

## 2 Lattice discretisation

The most obvious model lattice with spherical topology is given by a cube covered by a rectangular mesh of lattice points, a lattice type that will be denoted by (C) in the following. This inevitably introduces irregularities as compared to an idealized regular (large volume) triangulation of the sphere, which are given by a concentration of curvature around the cube corners and the defective coordination number of the corner points (three neighbours instead of four). In view of this problem several refinements have been proposed such as the replacement of the cube corners by triangular plaquettes, a construction that hides the coordination number defects away into the dual lattice, or the projection of the cube geometry (C) on the sphere $`S^2`$ by the application of appropriate weight factors to the links between lattice sites (S) . This latter lattice can be shown to arrive at the “right” continuum field theory in the thermodynamic limit . On the other hand, even simpler approximations could be devised such as for example a “pillow” resulting from the glueing together of rectangles along their four sides. Since, at least for bulk quantities, differences between those different choices of model lattices appear to be small , especially between (C) and its projected version (S), we here concentrate on the cubic discretisation (C) which is computationally much more convenient than (S). As far as universality is concerned, the ratios of asymptotic correlation lengths are found universal even with respect to marginal perturbations , so that we do not expect to see artefacts of discretisation defects here. Amplitude universality, however, might be sensitive to such distortions. Other types of lattice discretisations will be considered in a forthcoming publication .

The cube geometry (C) is being realised by six $`L\times L`$ square lattices appropriately glued together along their edges. Since the model lattices including type (C) are only topologically equivalent to spheres one has to conceive a definition of effective radii to enable a proper FSS analysis in the sense of eq. (3). Assigning a unit volume to each lattice site, to each pair of bonds, or to each individual square of the model lattice, respectively, on counting these entities one arrives at areas of

$$A=\{\begin{array}{ccc}6L& (L2)+8\hfill & \mathrm{`}\mathrm{`}\mathrm{sites}\mathrm{"},\hfill \\ 6L& (L2)+6\hfill & \mathrm{`}\mathrm{`}\mathrm{bonds}\mathrm{"},\mathrm{`}\mathrm{`}\mathrm{squares}\mathrm{"},\hfill \end{array}$$
(4)
from which the effective radii result via the relation $`R=\sqrt{A/4\pi }`$. Since for a “quadrangulation” one obviously has

$$4\mathrm{\#}(\mathrm{squares})=2\mathrm{\#}(\mathrm{bonds}),$$
(5)
the “bonds” and “squares” definitions coincide, but differ from the “sites” definition by a constant shift, leading to a slightly different approach to the leading FSS behaviour.

## 3 Model and simulation

We consider a classical, ferromagnetic, nearest-neighbour Ising model with Hamiltonian

$$=J\underset{i,j}{}\sigma _i\sigma _j,\sigma _i=\pm 1.$$
(6)
The spins reside on lattices compound of the cubical sphere discretisation (C) with edge lengths $`L`$ times a further linear lattice direction of length $`L_z`$. We apply periodic boundary conditions in $`z`$-direction to eliminate surface effects in modelling the infinite $``$ part. To control the effect of finite $`L_z`$ still present we enforce the condition $`L_z/\xi _{}1`$ in a systematic way, i.e. we scale $`L_z`$ proportionally to the radii $`R`$ and thus (to leading order) proportionally to $`L`$. We use a ratio of $`L_z/\xi _{}15`$ which empirically proves sufficient to let the effect of finite $`L_z`$ drop far below the threshold of statistical fluctuations. Simulations were done at a quite precise recent estimate for the Ising model critical coupling in three dimensions , $`\beta _c=\mathrm{0.221\hspace{0.17em}654\hspace{0.17em}4}(3)`$; the effect of critical coupling uncertainty was checked by a temperature reweighting technique and found negligible compared to statistical errors. For the FSS analysis we performed single-cluster Monte Carlo simulations of nine systems with cube edge lengths ranging from $`L=4`$ to $`L=12`$, corresponding to effective radii of $`R2`$ to $`R8`$, and collected up to about $`810^6`$ approximately independent measurements after an initial equilibration phase. Taking the scaled longitudinal sizes $`L_z`$ into account, system sizes thus range up to about $`310^5`$ spins.

Within the framework of CFT, the densities of magnetisation and energy turn out to be primary operators of spin models in two dimensions. Taking this as a hint for higher dimensions, we thus determine the longitudinal correlation lengths of these two observables to check whether Cardy’s conjecture eq. (3) holds in the discretised case. We start measuring the corresponding connected, longitudinal correlation functions $`G_{\sigma /ϵ}^{c,}(z)`$, maximally enhancing the statistics by appropriate averaging and projection techniques, for details see . We find that for the extraction of the correlation lengths from this information at a level of high precision the customary way of fitting the correlation function data to the functional form $`G(z)^{c,}=a\mathrm{exp}(z/\xi _{})+b`$ for large distances $`z`$ has severe drawbacks. Due to the statistical nature of the data one is not allowed to assume $`b=0`$, i.e. the theoretical limit for infinite-length time series, a priori for finite-length measurements and is thus forced to use intrinsically unstable non-linear fits. Therefore, we resort to a difference-ratio type estimator of the form

$$\widehat{\xi }_{}(z)=\mathrm{\Delta }\left[\mathrm{ln}\frac{G^{c,}(z)G^{c,}(z\mathrm{\Delta })}{G^{c,}(z+\mathrm{\Delta })G^{c,}(z)}\right]^1,$$
(7)
which eliminates the inconvenient constants $`a`$ and $`b`$ above. The parameter $`\mathrm{\Delta }1`$ controls the signal-noise ratio of the estimator and is being adapted for a minimum statistical error. To find final estimates for the correlation lengths $`\xi _{}(R)`$ of each system, we average the estimates $`\widehat{\xi _{}}(z)`$ over a regime of distances $`z`$ delimited by short distance deviations from the purely exponential behaviour and exploding fluctuations in the limit of very large distances. The averaging is done in a statistically optimized way entailing an expensive analysis of variances and covariances as well as autocorrelations, based on the “jackknife” resampling technique .

## 4 Results

Collecting the final estimates for the correlation lengths $`\xi _{,\sigma /ϵ}`$ of the densities of magnetisation and energy we arrive at the finite-size results shown in fig. 1(a), which scale linearly to leading order in the effective radii $`R`$ defined above. Differences between the two radii definitions can hardly be distinguished at this scale. Plotting the effective amplitudes $`\xi _{}/R`$, however, reveals the presence of corrections to the leading scaling behaviour that are clearly resolvable at the given level of accuracy, cp. fig. 1(b). We thus fit our data to a FSS ansatz including a leading correction term,

$$\xi _{}(R)=AR+BR^\alpha ,$$
(8)
treating $`\alpha `$ as an additional fit parameter that thus constitutes an effective correction exponent including parts of the higher-order corrections as well. This higher-order effect is kept small by successively removing points from the small $`R`$ end while monitoring the quality-of-fit parameter $`Q`$. The range of points used is indicated by the range of the fit lines in fig. 1. Combining the fits to the data for the two different radii definitions we arrive at final estimates for the leading correlation length scaling amplitudes and their ratio of

$$\begin{array}{c}A_\sigma =1.996(20),\hfill \\ A_ϵ=0.710(38),\hfill \\ A_\sigma /A_ϵ=2.81(15).\hfill \end{array}$$
(9)
Using weighted averages of recent estimates for the critical exponents of the three-dimensional Ising model , namely $`\nu =0.63005(18)`$ and $`\gamma =1.23717(28)`$, the conjectured amplitudes for comparison are

$$\begin{array}{c}A_\sigma ^{\mathrm{conj}}=1/x_\sigma =2/(d\gamma /\nu )=1.9298(13),\hfill \\ A_ϵ^{\mathrm{conj}}=1/x_ϵ=1/(d1/\nu )=0.70780(23),\hfill \\ A_\sigma ^{\mathrm{conj}}/A_ϵ^{\mathrm{conj}}=x_ϵ/x_\sigma =2.7264(13),\hfill \end{array}$$
(10)
where $`d=3`$ is understood, cp. eq. (3). Comparing the amplitude ratios, which are known to be more robust with respect to lattice inhomogeneities , we find very good agreement between the simulation and Cardy’s conjecture. The same holds true for the amplitude of the energy-energy correlation length whose estimate is naturally less precise than that of the spin-spin correlation length. For the latter we still obtain reasonable agreement, but with a deviation of the numerical estimate towards larger values which is on the edge of statistical significance. To check whether this result indicates a systematic deviation due to the irregularity of the used lattice discretisation (C), simulations for the other discretisations mentioned above in the introduction should be performed .

## 5 Conclusions

Approximating the sphere $`S^2`$ by a cube covered by a rectangular mesh of lattice points, we find that Cardy’s conjecture of a linear amplitude-exponent relation in the FSS of correlation lengths on the sphere-like geometry $`S^{d1}\times `$ holds true for the case of the three-dimensional Ising model and the densities of magnetisation and energy. In terms of the analogy to CFT in two dimensions the considered operators are thus analogues of “primary” operators in three dimensions. The ratio of leading scaling amplitudes does, at the given level of accuracy, not show any effect of the irregularity of lattice discretisation; for the amplitude of the spin-spin correlation length further simulations using different discretisation schemes are necessary to judge the extent of universality that is being obeyed. Our results indicate that the geometry $`S^2\times `$ considered here is closer to the original cylinder $`S^1\times `$ (as far as the prevalence of universal features of correlation lengths scaling is concerned) than the toroidal geometry $`T^2\times `$ which is not conformally flat, but nevertheless exhibits a scaling law quite similar to that found in two dimensions .

\***

MW gratefully acknowledges support by the “Deutsche Forschungsgemeinschaft” through the Graduiertenkolleg “Quantenfeldtheorie”.


format: text

text:
# Why Use a Hamilton Approach in QCD?

(1 August 2000)

## Abstract

We discuss $`QCD`$ in the Hamiltonian frame work. We treat finite density $`QCD`$ in the strong coupling regime. We present a parton-model inspired regularisation scheme to treat the spectrum ($`\theta `$-angles) and distribution functions in $`QED_{1+1}`$. We suggest a Monte Carlo method to construct low-dimensionasl effective Hamiltonians. Finally, we discuss improvement in Hamiltonian $`QCD`$.

1. Introduction
There are presently the following problems in $`QCD`$: How to compute non-perturbatively (on the lattice): (a) $`QCD`$ at non-zero chemical potential - finite quark density. (b) Hadron structure functions, at small $`x_B`$ and small $`Q`$. (c) $`S`$-matrix, scattering cross sections, phase shifts, (d) Decay amplitudes. Here, we want to discuss recent progress in the Hamiltonian formulation, $`QCD`$ at non-zero chemical potential in the strong coupling regime, distribution function in $`QED_{1+1}`$, the Monte Carlo Hamiltonian, and improvement of the continuum behavior of the Hamiltonian.

2. Finite density QCD
According to the big bang model of cosmology creation of quarks and gluons occured at about $`10^{35}`$ sec, quark-gluon plasma ($`QGP`$) was formed at about $`10^6`$ sec, then there was a phase transition and the formation of protons and neutrons began at $`10^4`$ sec. Today the universe is in a low temperature and low density phase: quarks are confined. Quark-gluon plasma may exist in cores of very dense stars like neutron stars. Experimental search for $`QGP`$ is carried out at the Relativistic Heavy Ion Collider (RHIC) at BNL and at the Large Hadron Collider (LHC) at CERN. Theorists search for a precise determination of phase structure of $`QCD`$ at finite temperature $`T`$ and finite chemical potential $`\mu `$ (finite density). Shuryak  has predicted that $`QCD`$ at high temperature and density displays beyond the $`QGP`$ pase transition a much richer phase structure, due to strong instanton effects.

Here we suggest to consider finite temperature and finite density $`QCD`$ in the Hamiltonian formulation. Consider the grand canonical partition function at temperature $`T`$ and chemical potential $`\mu `$

$`Z`$ $`=`$ $`Tre^{\beta (H\mu N)},N={\displaystyle d^3x\psi ^{}(x)\psi (x)},\beta =(k_BT)^1.`$ (1)
$`H`$ is the Hamiltonian and $`N`$ is the fermion (quark) particle number operator, No complex action problem occurs! Let us consider as example the energy density for free quarks

$$ϵ=\frac{1}{V}\frac{1}{Z}TrHe^{\beta (H\mu N)}=\frac{1}{V}\frac{lnZ}{\beta }|_{\mu ,\beta }.$$
(2)
Taking the zero-temperature limit ($`T0`$) and the chiral limit ($`m0`$) the energy density density (with the contribution $`\mu =0`$ subtracted) becomes $`ϵ_{sub}=\frac{\mu ^2}{4\pi ^2}`$. The same result is obtained by starting from a lattice Hamiltonian with Wilson fermions.

Strong coupling QCD at non-zero chemical potential.
Latttice $`QCD`$ at $`\mu =0`$ confines quarks and spontaneously breaks chiral symmetry. Note: The strong coupling regime $`1/g^2<<1`$ is not compatible with continuum limit $`a0`$. The goal of the following calculation is to get a better understanding of mechanism of chiral phase transition. We make an ansatz for the lattice Hamiltonian in the strong coupling regime following Ref.. Then the following results are obtained . The vacuum energy becomes

$`E_\mathrm{\Omega }`$ $`=`$ $`2N_fN_cN_s\left[(m_{dyn}^{(0)}\mu )\mathrm{\Theta }(\mu m_{dyn}^{(0)})m_{dyn}^{(0)}\mu \right]`$
$`m_{dyn}^{(0)}`$ $`=`$ $`{\displaystyle \frac{d}{ag^2C_N}},C_N=(N_c^21)/(2N_c),`$ (3)
where $`m_{dyn}^{(0)}`$ is the dynamical quark mass at $`\mu =0`$. In order to compute the chiral condensate and the critical chemical potential one can use the Feynman-Hellmann theorem. The chiral condensate at $`\mu =0`$ becomes

$$<\overline{\psi }\psi >^{(0)}=2N_C[14d/(g^2C_N^2)],$$
(4)
and the chiral condensate at $`\mu 0`$ becomes

$$<\overline{\psi }\psi >=<\overline{\psi }\psi >^{(0)}\left[1\mathrm{\Theta }(\mu m_{dyn}^{(0)})\right].$$
(5)
As a result one has

$`<\overline{\psi }\psi >0\text{for}\mu <m_{dyn}^{(0)}:\text{chiral symmetry broken},`$
$`<\overline{\psi }\psi >=0\text{for}\mu >m_{dyn}^{(0)}:\text{chiral symmetry restored}.`$ (6)
Thus we make the prediction of a first order chiral phase transition. The critical chemical potential is $`\mu _{crit}=m_{dyn}^{(0)}`$. Our result for $`\mu _{crit}`$ is in agreement with the result obtained by Le Yaouanac et al. However, Le Yaouanac et al. find a phase transition of 2nd order. There have been predictions of the nucleon mass  $`M_{nucl}^{(0)}3m_{dyn}^{(0)}`$. Thus our result is $`\mu _{crit}M_{nucl}^{(0)}/3`$, being consistent with other theoretical predictions . We have moreover computed the quark number density, the corresponding susceptibility and some thermal masses .

3. Hamiltonian in a fast moving frame
At hand of the massive Schwinger model, $`QED_{1+1}`$, we want to discuss the following topics: (i) A physically motivated many-body method to solve a Hamiltonian many-body problem and compute structure functions. (ii) The treatment of the $`\theta `$-angle in Hamiltonian gauge theory. The motivation for the Hamiltonian many-body method comes from the parton model of deep inelastic lepton-hadron scattering. Experimentally, the Breit frame is a convenient choice for the interpretation of structure functions. We suggest to use a fast moving frame also in a computational study . The Hamiltonian in axial gauge $`A^3=0`$ is given by

$$H=_L^L𝑑x^3(\overline{\psi }\gamma ^3i_3\psi +m\overline{\psi }\psi )+\frac{g^2}{2}_L^L𝑑x^3(\psi ^{}\psi )\frac{1}{_3^2}(\psi ^{}\psi ).$$
(7)
We start from a space-time lattice of spacing $`a`$ and length $`[L,L]`$. Then we go over to a momentum lattice with cut-off $`\mathrm{\Lambda }=\pi /a`$ and resolution $`\mathrm{\Delta }p=\pi /L`$. We use the assumption (parton model): Left moving particles are not dynamically important, when physical particles are considered from a reference frame with velocity $`v=P/Ec`$. On the momentum lattice the following bounds hold: $`0p^0,p^3P`$. A fast moving particle is Lorentz contracted. It fits into a small lattice volume (compared to rest frame). A small lattice size implies a small number of virtual particle pairs created from vacuum: Number of vacuum pairs $``$ vacuum density $`\times `$ lattice size. For the scaling window holds: $`a\text{Compton wavelength}L`$ and $`a\xi aL`$, where $`\xi `$ is the correlation length in dimensionless units and $`M=\frac{1}{\xi a}`$ is the physical mass of particle in the ground state. In the strong relativistic regime holds $`M<<P<\mathrm{\Lambda }`$. In the limit $`a0`$ the scaling window is replaced by: $`1/P1/ML`$. In Ref. we have computed the mass spectrum: Vector boson and scalar boson, and compared it with predictions of chiral perturbation theory. We also determined the dependence on the $`\theta `$ angle, defined by a global gauge transformation $`\mathrm{\Psi }[A]e^{in\theta }\mathrm{\Psi }[A],n=0,\pm 1,\pm 2,\mathrm{}`$. The $`\theta `$-action $`S`$ and topological charge $`q`$ are related by $`S[A,\theta ]=S[A]\theta q[A]`$. In Ref.  we have also studied the massive Schwinger model and its dependence on the $`\theta `$-angle in the chiral limit. For the massless Schwinger model the chiral condensate is analytically known,

$$\frac{<\overline{\psi }\psi >}{g}=\frac{e^\gamma }{2\pi \sqrt{\pi }}\mathrm{cos}(\theta )0.1599\mathrm{cos}(\theta ).$$
(8)
Using the Feynman-Hellmann theorem to relate the vector mass and the chiral condensate, we determine the slope from numerical data and obtain $`\frac{<\overline{\psi }\psi >}{g}0.16\mathrm{cos}(\theta )`$, which is close to the analytical result.

The Hamiltonian approach also allows to compute wave functions, giving information on the parton structure.

$$f(x_B)=<\mathrm{\Psi }_V(P)|\frac{1}{2}[b_k^{}b_k+d_k^{}d_k|\mathrm{\Psi }_V(P)>,$$
(9)
where $`x_B=\frac{k}{P}`$ is the fraction of momentum of vector boson carried by parton. One observes that a large lattice (N=200) is needed to resolve distribution functions precisely (see Ref.).

4. Stochastic Hamiltonian
Lagrangian lattice field theory and its non-perturbative numerical solution by Monte Carlo methods is very successful. The typical number of configurations (depending on the statistical error) is in the order of $`10^2`$. In Ref., we have propopsed to construct an effective Hamiltonian, by treating the many-body problem by applying as far as possible Monte Carlo with importance sampling (e.g. Metropolis algorithm) in analogy to Lagrangian lattice simulations. Let us consider as example 1-D Q.M. We compute transition matrix elements in imaginary time

$$M_{ij}(T)=<x_i|e^{HT/\mathrm{}}|x_j>,i,j1,\mathrm{},N,$$
(10)
where $`x_i`$ are discrete position eigen states. We define an effective Hamiltonian via eigenfunctions and eigenvalues:

$$M(T)=U^{}D(T)U,U_{ik}^{}=<e_i|E_k^{eff}>,D_i(T)=e^{E_i^{eff}T/\mathrm{}}.$$
(11)
We compute matrix elements using Monte Carlo with importance sampling (e.g. Metropolis algorithm). Splitting the action $`S=S_0+S_V`$, we treat $`O=exp[S_V/\mathrm{}]`$ as observable and compute the matrix element by Monte Carlo

$`M_{ij}(T)=M_{ij}^{(0)}(T)`$ $`\times `$ $`{\displaystyle [dx]exp[S_V[x]/\mathrm{}]exp[S_0[x]/\mathrm{}]}|_{x_j,0}^{x_i,T}`$ (12)
$`/`$ $`{\displaystyle [dx]exp[S_0[x]/\mathrm{}]}|_{x_j,0}^{x_i,T}.`$
The construction of a stochastic basis has been proposed in Ref.. The result of a numerical study is shown in Fig..

5. Improvement in Hamiltonian QCD
Szymanzik  suggested improvement program for lattice field theory: Construction of ”improved” action, such that error due to finite lattice spacing $`a`$ becomes smaller. Progress in improvement for lattice pure gauge theory was introduced by Lepage et al. , suggesting tadpole improvement. We have suggested how to construct an improved Hamiltonian on the lattice. For pure gauge theory the construction of an improved Hamiltonian for $`SU(3)`$ gauge theory via the transfer matrix has been reported by by Luo et al. in Ref.. Improvement of the fermionic sector of a lattice Hamiltonian has been carried out for $`QCD_{1+1}`$ by Luo et al. in Ref.. A numerical study was carried out by Jiang et al. . An example is shown in Fig., showing that improvement reduces the dependence on Wilson’s r-parameter.

6. Conclusion
Hamiltonian formulation of QCD is attractive, for the following reasons: (a) There is no sign or complex action problem! (b) It is suitable for the computation of wave functions and structure functions. (c) It is suitable for the computation of the S-matrix. The problem is: How to solve the Hamiltonian many-body system? (d) Proposal: Construct an effective Hamiltonian via Monte Carlo. (e) Improvement program: Also possible in Hamiltonian formulation.

Acknowledgement. H. Kröger and K.J.M. Moriarty are grateful for support by NSERC(Canada). X.Q. Luo has been supported by NNSF(China) and NSF(China) for Distinguished Young Scholars.


text:
# A Bayesian approach to source separation Presented at the 19th Int. worskhop on Bayesian and Maximum Entropy methods (MaxEnt 1999), Aug. 2-6, 1999, Boise, Idaho, USA

## 1. Introduction

The simplest model for a source separation is

$$𝒙(t)=𝑨𝒔(t),$$
(1)
where $`𝑨`$ is a mixing matrix, $`𝒔(t)`$ is a vector of sources and $`𝒙(t)`$ a vector of independent measurements. The main task is then to recover $`𝒔(t)`$, but one may instead be interested in recovering a separating matrix $`𝑩`$ such that $`\widehat{𝒔}(t)=𝑩𝒙(t)`$. When $`𝑨`$ is invertible, it is natural to assume that $`𝑩=𝑨^1`$ or $`𝑩=𝚺𝚲𝑨^1`$ where $`𝚺`$ is a permutation matrix and $`𝚲`$ a diagonal scaling matrix.

Many source separation algorithms have been recently proposed based on likelihood , contrast function , estimating function , information theory , and more generally on principle component analysis (PCA) , Independent factor analysis (IFA)  and independent component analysis (ICA) . All these methods assume that the mixing matrix $`𝑨`$ is invertible and mainly search for a separating matrix $`𝑩`$ such that the components of $`𝒚(t)=𝑩𝒙(t)`$ be independent. This means that all these methods implicitly assume that the sources $`𝒔(t)`$ are independent. This may not be the case in some applications. However, the main differences between these methods are in the way they try to insure this independence.

* Maximum likelihood (ML) techniques use directly the independence property by assuming

$$p(𝒔)=\underset{i}{}p_i(s_i)$$
(2)
and as a result

$$p(𝒙|𝑨)=\frac{1}{\text{det}(𝑨)}\underset{i}{}p_i\left([𝑨^1𝒙]_i\right),$$
(3)
or equivalently

$$p(𝒙|𝑩)=|\text{det}(𝑩)|\underset{i}{}p_i\left([𝑩𝒙]_i\right)=|\text{det}(𝑩)|\underset{i}{}p_i(y_i),$$
(4)
where $`y_i=[𝑨^1𝒙]_i=[𝑩𝒙]_i`$ and where $`p(𝒔)`$ is the probability density function of the source vector $`𝒔`$. The ML estimate of the separating matrix is defined as

$`\widehat{𝑩}`$ $`=`$ $`\underset{𝑩}{\text{arg}\mathrm{max}}\left\{\mathrm{log}p(𝒙|𝑩)\right\}`$ (5)
$`=`$ $`\underset{𝑩}{\text{arg}\mathrm{max}}\left\{{\displaystyle \underset{i}{}}r_i(y_i)+\mathrm{log}|\text{det}(𝑩)|\right\}\text{with}r_i(y_i)=\mathrm{log}p_i(y_i).`$
A great number of algorithms have been proposed to perform this optimization .
* Infomax techniques use the entropy of $`𝒚=𝑩𝒙`$ as a measure of independence :

$$S=\underset{i}{}p_i(y_i)\mathrm{log}p_i(y_i).$$
(6)
Thus $`S`$ is a function of the separating matrix $`𝑩`$ and one tries to optimize $`S`$ with respect to $`𝑩`$.
* M-estimation techniques define an estimate for the separating matrix $`𝑩`$ such that

$$\frac{1}{T}\underset{t}{}𝑯\left(𝒚(t)\right)=\left[\begin{array}{c}0\end{array}\right]\text{with}𝒚(t)=𝑩𝒙(t),$$
(7)
where it is assumed to have $`T`$ independent observations $`\{𝒙(1),\mathrm{},𝒙(T)\}`$ and where $`𝑯`$ is an appropriately defined matrix valued function. $`\left[\begin{array}{c}0\end{array}\right]`$ represents a matrix whose elements are all equal to zero. We can note that M-estimate methods generalize the ML estimation method since the latter can be obtained by taking

$$𝑯\left(𝒚\right)=\frac{\mathrm{log}p(𝒚)}{𝑩}.$$
(8)
* Contrast function minimization techniques are based on the optimization of a contrast function $`c(𝒚)=c\left(𝑩𝒙\right)`$ which takes its extremal value when $`𝑩`$ is a separating matrix . Typical examples are the contrast functions measuring, in some way, the independence of the components of $`𝒚`$, sometimes subject to the constraint that $`𝒚(t)`$ be spatially white

$$\frac{1}{T}\underset{t}{}𝒚(t)𝒚^{}(t)=𝑰.$$
(9)
* Higher order statistics (HOS) techniques try to insure the independence of the components of $`𝒚`$ by minimizing, under the whiteness constraint, a contrast function related to the statistics of the order greater that two such as the cumulants .

The main limitations of these techniques are the following:

* None of these techniques consider the possible errors on the model or the measurement (sensor) noises;
* All these methods assume that the mixing matrix $`𝑨`$ is invertible and cannot account for the cases in which $`𝑨`$ is rectangular (number of sensors different from the number of sources).
* All these methods assume that the sources are independent. Some assume the sources to be also temporally white.

Recently, a few works using the Bayesian approach have been presented to push farther the limits of these methods .

In the following, we first present the basics of the Bayesian approach, then we show how some of the preceeding techniques can be obtained as special cases, and finally, we propose new ideas to account for spatial correlation between neighbor sources or time correlation of the sources.

## 2. Bayesian approach

The main idea in the Bayesian approach is to use not only the likelihood
$`p(𝒙(1),\mathrm{},𝒙(T)|𝑨)`$ but also some prior knowledge about the sources $`𝒔`$ and the mixing matrix $`𝑨`$ through the assignment of prior probabilies $`p(𝒔)`$ and $`p(𝑨)`$. Then, noting $`𝒙_{1..T}=\{𝒙(1),\mathrm{},𝒙(T)\}`$ and $`𝒔_{1..T}=\{𝒔(1),\mathrm{},𝒔(T)\}`$ and using these direct probability laws we determine the posterior law

$$\mathrm{log}p(𝑨,𝒔_{1..T}|𝒙_{1..T})=\mathrm{log}p\left(𝒙_{1..T}|𝑨,𝒔_{1..T}\right)+\mathrm{log}p(𝑨)+T\mathrm{log}p(𝒔)+cte,$$
(10)
where we assumed the independence of the sources $`𝒔`$ and the mixing matrix $`𝑨`$.

From this posterior probability law we can deduce any inference about $`𝑨`$ and $`𝒔`$. For example, we can estimate both $`𝑨`$ and $`𝒔`$ by a joint maximum a posteriori  (JMAP) criterion using an alternate maximization algorithm. We can also focus on the estimation of the mixing matrix $`𝑨`$ by marginalizing this posterior law with respect to $`𝒔`$ to obtain $`p(𝑨|𝒙_{1..T})`$ and use the resulting MAP criterion to estimate $`𝑨`$. Finally, we can integrate $`𝑨`$ from the joint law to obtain $`p(𝒔_{1..T}|𝒙_{1..T})`$ and estimate the sources from this marginal posterior law.

Now, before going further in details of these three methods, we are going to illustrate some special cases which result in some classical techniques.

### 2.1. Exact invertible model and independent sources

If we assume that the model $`𝒙=𝑨𝒔`$ is exact and that there is not any measurement noise and that the mixing matrix $`𝑨`$ is invertible and well conditioned, then we can only look for a separating matrix $`𝑩=𝑨^1`$. Indeed, as in conventional methods, if we assume that the sources $`𝒔`$ are independent, we have the following relations:

$$p(𝒔)=\underset{i}{}p_i(s_i)$$
(11)
and so

$$p(𝒙|𝑩)=|\text{det}(𝑩)|\underset{i}{}p_i\left([𝑩𝒙]_i\right),$$
(12)
where $`p_i(s_i)`$ is the probability density function of the source component $`i`$. Using these relations, and noting by $`𝒚(t)=𝑩𝒙(t)`$ we have

$$\mathrm{log}p\left(𝑩|𝒙_{1..T}\right)=\mathrm{log}p\left(𝒙_{1..T}|𝑩\right)+\mathrm{log}p(𝑩)+cte,$$
(13)
where

$$\mathrm{log}p\left(𝒙_{1..T}|𝑩\right)=T\mathrm{log}|\text{det}(𝑩)|+\underset{t}{}\underset{i}{}\mathrm{log}p_i\left(y_i(t)\right)$$

and $`p(𝑩)`$ is a probability distribution on the separating matrix $`𝑩`$. Here we assume that we can assign a probability law $`p(𝑨)`$ to the mixing matrix $`𝑨`$ or equivalently $`p(𝑩)`$ to the separating matrix $`𝑩`$ to translate any prior knowledge we have about (or we wish to impose to) them. For example, we may know (or assume) that the mixing matrix is such that

𝑨2

=
kl|ak,l|2ϵ,superscriptnorm𝑨2

=
subscript𝑘subscript𝑙superscriptsubscript𝑎𝑘𝑙2italic-ϵ\|{\mbox{\boldmath$A$}}\|^{2}\,\shortstack{$\triangle$\\
 =}\,\sum_{k}\sum_{l}|a_{k,l}|^{2}\leq\epsilon, (14)
for some $`ϵ`$; or we may wish that the separating matrix $`𝑩`$ be such that its determinant $`|\text{det}(𝑩)|0`$ and not very far from one. In the first case we can choose

$$p(𝑨)\mathrm{exp}\left[\frac{1}{2\sigma _a^2}𝑨^2\right]=\mathrm{exp}\left[\frac{1}{2\sigma _a^2}\underset{k}{}\underset{l}{}a_{k,l}^2\right]$$
(15)
and in the second case

$$p(𝑩)|\text{det}(𝑩)|.$$
(16)
Some other possibilities are:

$$p(𝑨)\mathrm{exp}\left[\frac{1}{2\sigma _a^2}𝑰𝑨^2\right]=\mathrm{exp}\left[\frac{1}{2\sigma _a^2}\left[\underset{k}{}(1a_{k,k})^2+\underset{lk}{}a_{k,l}^2\right]\right]$$
(17)
which tries to impose $`|a_{k,l}|1,k=l`$ and $`|a_{k,l}|0,kl`$;

$$p(𝑨)\mathrm{exp}\left[\frac{1}{2\sigma _a^2}𝑰𝑨𝑨^t^2\right]=\mathrm{exp}\left[\frac{1}{2\sigma _a^2}\left[\underset{k}{}(1a_{k,}^2)^2\underset{lk}{}[𝑨𝑨^t]_{k,l}^2\right]\right]$$
(18)
when the number of sources is less than the number of the sensors; and

$$p(𝑨)\mathrm{exp}\left[\frac{1}{2\sigma _a^2}𝑰𝑨^t𝑨^2\right]=\mathrm{exp}\left[\frac{1}{2\sigma _a^2}\left[\underset{l}{}(1a_{,l}^2)^2\underset{kl}{}[𝑨^t𝑨]_{k,l}^2\right]\right]$$
(19)
when the number of sources is greater than the number of the sensores. These two last expressions have been proposed and used by Knuth . Other choices based on prior knowledge of the geometrical positions of the sources and receivers and knowledge of the signal propagation law for an acoustical application have been used by .

Now, if we consider the MAP estimation, the MAP criterion to optimize becomes

$$J(𝑩)=\mathrm{log}p\left(𝑩|𝒙_{1..T}\right)=T\mathrm{log}|\text{det}(𝑩)|+\underset{t}{}\underset{i}{}\mathrm{log}p_i\left(y_i(t)\right)+\mathrm{log}p(𝑩)+cte.$$
(20)
Searching now for the MAP solution, the necessary condition is

$$\frac{J(𝑩)}{𝑩}=\left[\begin{array}{c}0\end{array}\right]\underset{t}{}𝑯\left(𝒚(t)\right)=\left[\begin{array}{c}0\end{array}\right],$$
(21)
where $`𝑯`$ is a matrix valued function given by

$$H(𝒚)=\frac{}{𝑩}\left[\underset{i}{}\mathrm{log}p_i\left(y_i\right)+\mathrm{log}|\text{det}(𝑩)|+\frac{1}{T}\mathrm{log}p(𝑩)\right].$$
(22)
As an example, consider a uniform a priori law for $`𝑩`$. Then we obtain the classical ML estimate which satisfies

$$\underset{t}{}𝑯\left(𝒚(t)\right)=\left[\begin{array}{c}0\end{array}\right]\text{ with }𝑯\left(𝒚\right)=\mathit{\varphi }(𝒚)𝒚^t𝑰,$$
(23)
where $`\mathit{\varphi }(𝒚)=[\varphi _1(y_1),\mathrm{},\varphi _n(y_n)]^t`$ with

$$\varphi _i(z)=\frac{p_i^{}(z)}{p_i(z)}.$$
(24)
One can add some extra constraints to this optimization. For example, we can optimize the MAP criterion subject to the constraint $`\frac{1}{T}_t𝒚(t)𝒚^t(t)=𝑰`$ which leads again to

$$\underset{t}{}𝑯\left(𝒚(t)\right)=\left[\begin{array}{c}0\end{array}\right]\text{ with }𝑯\left(𝒚\right)=\alpha (𝒚𝒚^t𝑰)+\beta \left(\mathit{\varphi }(𝒚)𝒚^t+𝒚\mathit{\varphi }^t(𝒚)\right).$$
(25)
Note that in all these relations, $`\varphi _i(z)`$ is related to the probability distribution of the source number $`i`$. The following table gives the expression of this function for a few known cases.

$$\begin{array}{ccc}& & \\ & & \\ \text{Gauss}\hfill & p(z)\mathrm{exp}\left[\alpha z^2\right]\hfill & \varphi (z)=2\alpha z\hfill \\ & & \\ \text{Laplace}\hfill & p(z)\mathrm{exp}\left[\alpha |z|\right]\hfill & \varphi (z)=\alpha \text{sign}(z)\hfill \\ & & \\ \text{Cauchy}\hfill & p(z)\frac{1}{1+(z/\alpha )^2}\hfill & \varphi (z)=\frac{2z/\alpha ^2}{1+(z/\alpha )^2}\hfill \\ & & \\ \text{Gamma}\hfill & p(z)z^\alpha \mathrm{exp}\left[\beta z\right]\hfill & \varphi (z)=\alpha /z+\beta \hfill \\ & & \\ \text{sub-Gaussian law}\hfill & p(z)\mathrm{exp}\left[\frac{1}{2}z^2\right]\text{sech}^2(z)\hfill & \varphi (z)=z+\mathrm{tanh}(z)\hfill \\ & & \\ \text{Mixture of Gaussians}\hfill & \begin{array}{cc}p(z)\hfill & \mathrm{exp}\left[\frac{1}{2}(z\alpha )^2\right]\hfill \\ & +\mathrm{exp}\left[\frac{1}{2}(z+\alpha )^2\right]\hfill \end{array}\hfill & \varphi (z)=\alpha z\alpha \mathrm{tanh}(\alpha z)\hfill \\ & & \end{array}$$

Remark:
$`H(𝒚)`$ in equations (21) and (23) corresponds to the gradient of MAP and ML criteria. A common technique to obtain the MAP or the ML solutions is then to use a gradient based algorithm such as

$$𝑩^{(k+1)}=𝑩^{(k)}\gamma 𝑯(𝒚)$$
(26)
where $`(k)`$ and $`(k+1)`$ stand for two successive iterations in static case or two successive time instants for dynamic case. This equation forms the main body of a great number of neural network (NN) based algorithms for source separation.

### 2.2. Accounting for errors

Here we relax the inversibility of the matrix $`𝑨`$ and take also account of the errors on the data. As an example, we consider the case where the errors can be modelled by an additive term $`\mathit{ϵ}(t)`$:

$$𝒙(t)=𝑨𝒔(t)+\mathit{ϵ}(t),t=1,\mathrm{},T.$$
(27)
We assume also that we can assign a probability law $`p(\mathit{ϵ})`$ to $`\mathit{ϵ}`$. In general, it is natural to assume that $`\mathit{ϵ}(t)`$ has independent components and is centered and white, i.e.

$$\mathrm{log}p(ϵ(1),\mathrm{},ϵ(T))=\underset{t}{}\underset{i}{}\mathrm{log}p_i\left(ϵ_i(t)\right).$$
(28)
From this assumption, we obtain

$$\mathrm{log}p\left(𝒙_{1..T}|𝑨,𝒔_{1..T}\right)=\underset{t}{}\underset{i}{}q_i\left(x_i(t)[𝑨𝒔]_i(t)\right)$$
(29)
with $`q_i(.)=\mathrm{log}p_i(.)`$.

Now, we can give the expression of the posterior law which is

$`\mathrm{log}p(𝑨,𝒔_{1..T}|𝒙_{1..T})`$ $`=`$ $`\mathrm{log}p\left(𝒙_{1..T}|𝑨,𝒔_{1..T}\right)+\mathrm{log}p(𝒔_{1..T})+\mathrm{log}p(𝑨)+cte`$
$`=`$ $`{\displaystyle \underset{t}{}}{\displaystyle \underset{i}{}}q_i\left(x_i(t)[𝑨𝒔]_i(t)\right)+\mathrm{log}p(𝒔_{1..T})+\mathrm{log}p(𝑨)+cte.`$
As mentioned before, from here, we can go in at least three directions:

* First integrate $`p(𝑨,𝒔_{1..T}|𝒙_{1..T})`$ with respect to $`𝑨`$ to obtain $`p(𝒔_{1..T}|𝒙_{1..T})`$ and estimate $`𝒔_{1..T}`$ by

$$\widehat{𝒔}_{1..T}=\underset{𝒔_{1..T}}{\text{arg}\mathrm{max}}\left\{p(𝒔_{1..T}|𝒙_{1..T})\right\}.$$
(31)
* Second integrate $`p(𝑨,𝒔_{1..T}|𝒙_{1..T})`$ with respect to $`𝒔_{1..T}`$ to obtain $`p(𝑨|𝒙_{1..T})`$ and estimate $`𝑨`$ by

$$\widehat{𝑨}=\underset{𝑨}{\text{arg}\mathrm{max}}\left\{p(𝑨|𝒙_{1..T})\right\}.$$
(32)
But, here, when $`\widehat{𝑨}`$ is obtained, we still have to obtain $`\widehat{𝑩}=\widehat{𝑨}^1`$ and $`\widehat{𝑨}`$ may not be invertible.
* Third, optimize $`p(𝑨,𝒔_{1..T}|𝒙_{1..T})`$ simultaneously with respect to both $`𝒔_{1..T}`$ and $`𝑨`$ by using an alternating optimization procedure such as

$$\{\begin{array}{ccc}\widehat{𝒔}_{1..T}^{(k)}\hfill & =& \underset{𝒔_{1..T}}{\text{arg}\mathrm{max}}\left\{p(\widehat{𝑨}^{(k1)},𝒔_{1..T}|𝒙_{1..T})\right\}\hfill \\ \widehat{𝑨}^{(k)}\hfill & =& \underset{𝑨}{\text{arg}\mathrm{max}}\left\{p(𝑨,\widehat{𝒔}_{1..T}^{(k1)}|𝒙_{1..T})\right\}\hfill \end{array}$$
(33)

In the two first cases, the integrations can be done analytically only in the Gaussian case. We then obtain closed form expressions for the solutions.

In any case, before applying any optimization, we have to ensure that the criterion to be optimized has at least an optimum and that this optimum is unique.

### 2.3. Spatially independent and white sources

The case where we can assume that the sources are independent and white is the simplest one. We have:

$$\mathrm{log}p(𝒔_{1..T})=\underset{t}{}\underset{j}{}r_j(s_j(t))$$

and

$$\mathrm{log}p(𝑨,𝒔_{1..T}|𝒙_{1..T})=\underset{t}{}\underset{i}{}q_i\left(x_i(t)z_i(t)\right)+\underset{t}{}\underset{j}{}r_j(s_j(t))+\mathrm{ln}p(𝑨)+cte.$$
(34)
with $`z_i=[𝑨𝒔]_i`$. Then, we can omit the time summation. To simplify the details of the derivations, let first assume

$$p(𝑨)\mathrm{exp}\left[\frac{1}{2\sigma _a^2}𝑨^2\right]=\mathrm{exp}\left[\frac{1}{2\sigma _a^2}\underset{k}{}\underset{l}{}a_{k,l}^2\right].$$
(35)
Later, we will also consider other possibilities such as (17), (18) or (19).

#### Joint MAP estimation

First we consider the joint estimation of $`𝑨`$ and $`𝒔`$ where the alternating optimization algorithm becomes

$$\{\begin{array}{ccc}\widehat{𝒔}^{(k)}\hfill & =& \underset{𝒔}{\text{arg}\mathrm{max}}\left\{\underset{i}{}q_i\left(x_iz_i\right)+\underset{j}{}r_j(s_j)\right\}\hfill \\ \widehat{𝑨}^{(k)}\hfill & =& \underset{𝑨}{\text{arg}\mathrm{max}}\left\{\underset{i}{}q_i\left(x_iz_i\right)\frac{1}{2\sigma _a^2}\underset{k}{}\underset{l}{}a_{kl}^2\right\}\hfill \end{array}$$
(36)
The solution at each iteration has to satisfy

$$\{\begin{array}{c}\frac{}{s_j}=\underset{i}{}a_{ij}q_i^{}\left(x_iz_i\right)+r_j^{}(s_j)=0\hfill \\ \frac{}{a_{ij}}=s_jq_i^{}\left(x_iz_i\right)\frac{1}{\sigma _a^2}a_{ij}=0\hfill \end{array}$$
(37)
These equations are in general nonlinear and depend on the expressions of $`q`$ and $`r`$. One exception is the particular case of Gaussian laws

$$p_i(n)𝒩(0,\sigma _ϵ^2)q_i(n)=\frac{1}{2\sigma _ϵ^2}n^2q_i^{}(n)=\frac{1}{\sigma _ϵ^2}n$$

and

$$p_j(s)𝒩(0,\sigma _s^2)r_j(s)=\frac{1}{2\sigma _s^2}s^2r_i^{}(s)=\frac{1}{\sigma _s^2}s$$

where we obtain two sets of linear equations to solve for $`s_j`$ and $`a_{ij}`$ :

$$\{\begin{array}{c}\frac{1}{\sigma _ϵ^2}\underset{i}{}a_{ij}\left(x_i[𝑨𝒔]_i\right)\frac{1}{\sigma _s^2}s_j=0\hfill \\ \frac{1}{\sigma _ϵ^2}s_j\left(x_i[𝑨𝒔]_i\right)\frac{1}{\sigma _a^2}a_{ij}=0\hfill \end{array}\{\begin{array}{c}\underset{i}{}a_{ij}\left(x_i[𝑨𝒔]_i\right)\lambda s_j=0\hfill \\ s_j\left(x_i[𝑨𝒔]_i\right)\mu a_{ij}=0\hfill \end{array}$$
(38)
with $`\lambda =\sigma _ϵ^2/\sigma _s^2`$ and $`\mu =\sigma _ϵ^2/\sigma _a^2`$.

These two equations have to be solved in each iteration of alternating optimization procedure. Two strategies can be used :

* Solve these equations for each $`s_j`$ and then for each $`a_{ij}`$ at each iteration:

$$\{\begin{array}{ccc}s_j=\frac{_ia_{ij}(x_i\widehat{x}_i)}{\lambda +𝒂_j^2}\hfill & & \\ a_{ij}=\frac{s_j(x_i\widehat{x}_i)}{s_j^2+\mu }\hfill & & \end{array}$$
(39)
with $`\widehat{x}_i=_{kj}a_{ik}s_k`$, and $`𝒂_i^2=_ja_{i,j}^2`$. This is a single coordinate-wise gradient descent based algorithm.
* Solve these equations for all $`s_j`$ and then for all $`a_{ij}`$ at each iteration:

$$\{\begin{array}{ccc}𝑨^t(𝒙𝑨𝒔)\lambda 𝒔\hfill & =& \mathrm{𝟎}\hfill \\ (𝒙𝑨𝒔)𝒔^t\mu 𝑨\hfill & =& \mathrm{𝟎}\hfill \end{array}\{\begin{array}{ccc}𝒔\hfill & =& (𝑨^t𝑨+\lambda 𝑰)^1𝑨^t𝒙\hfill \\ 𝑨\hfill & =& 𝒙𝒔^t(𝒔𝒔^t+\mu 𝑰)^1\hfill \\ & =& \frac{𝒙𝒔^t}{\mu }\left[𝑰\frac{𝒔𝒔^t}{𝒔^t𝒔+\mu }\right]\hfill \end{array}$$
(40)
This is a bloc coordinate-wise gradient descent based algorithm.

Remark 1:
These two last closed form expressions give us the possibility to discuss the convergency of the joint MAP algorithm in the considered Gaussian case. We may immediately note that $`𝑨`$ obtained by this algorithm is not invertible. This means that, in the Gaussian hypothesis, this algorithm does not really separate the signals. Actually, we could remark this from the expression of the joint criterion in this case which is

$$J(𝑨,𝒔)=\mathrm{log}p(𝑨,𝒔|𝒙)=𝒙𝑨𝒔^2+\lambda 𝒔^2+\mu 𝑨^2+cte.$$
(41)
As we can see, in this case, $`J(𝑨,𝒔)`$ is a quadratic function of $`𝒔`$ for given $`𝑨`$ and a quadratic function of $`𝑨`$ for given $`𝒔`$, but it is a biquadratic function of both $`𝑨`$ and $`𝒔`$. This symetry property means that the joint MAP solution in this case is not unique. This criterion may even have an infinite equivalent optima. The proposed iterative algorithm may then converge to any of these solutions depending on the initialization.

Unfortunately, with any non Gaussian hypothesis, we can not obtain any closed form solution and the existance and the uniqueness of a global optimum is very hard to study. However, we can always propose either a fixed point or a gradient descent based algorithm to compute them numerically. For example, if we assume a Gaussian law for the noise, but non Gaussian prior laws for $`𝒔`$ and for $`𝑨`$, we have

$$J(𝑨,𝒔)=\mathrm{log}p(𝑨,𝒔|𝒙)=𝒙𝑨𝒔^2+\lambda \varphi (𝒔)+\mu \psi (𝑨)+cte.$$
(42)
where $`\varphi (𝒔)\mathrm{log}p(𝒔)`$ and $`\psi (𝒔)\mathrm{log}p(𝑨)`$. Note that the choice of these prior laws is then important if we want to eliminate the above mentionned symetry property and to be able to find a unique solution to the problem. Then, a gradient based algorithm writes:

$$\{\begin{array}{ccc}\widehat{𝒔}^{(k+1)}\hfill & =& \widehat{𝒔}^{(k)}\alpha ^{(k)}\frac{J}{𝒔}(\widehat{𝑨}^{(k)},\widehat{𝒔}^{(k)})\hfill \\ \widehat{𝑨}^{(k+1)}\hfill & =& \widehat{𝑨}^{(k)}\beta ^{(k)}\frac{J}{𝑨}(\widehat{𝑨}^{(k)},\widehat{𝒔}^{(k)})\hfill \end{array}$$
(43)
where $`\alpha ^{(k)}`$ and $`\beta ^{(k)}`$ are two step parameters which can be either constant (fixed step gradient) or adaptive during the iterations $`(k)`$. Replacing for the gradient expressions we obtain:

$$\{\begin{array}{ccc}\widehat{𝒔}^{(k+1)}\hfill & =& \widehat{𝒔}^{(k)}+\alpha ^{(k)}\left[2\widehat{𝑨}^{t^{(k)}}(𝒙\widehat{𝒙}^{(k)})+\lambda \frac{\varphi }{𝒔}(\widehat{𝒔}^{(k)})\right]\hfill \\ \widehat{𝑨}^{(k+1)}\hfill & =& \widehat{𝑨}^{(k)}+\beta ^{(k)}\left[2\widehat{𝒔}^{(k)}(𝒙\widehat{𝒙}^{(k)})^t+\mu \frac{\psi }{𝑨}(\widehat{𝑨}^{(k)})\right]\hfill \end{array}$$
(44)
with $`\widehat{𝒙}^{(k)}=\widehat{𝑨}^{(k)}\widehat{𝒔}^{(k)}`$.

A fixed point based algorithm writes:

$$\{\begin{array}{ccc}\frac{\varphi }{𝒔}(\widehat{𝒔}^{(k)})\hfill & =& \frac{1}{\lambda }\left[\widehat{𝑨}^{t^{(k)}}(𝒙\widehat{𝒙}^{(k)})\right]\hfill \\ \frac{\psi }{𝑨}(\widehat{𝑨}^{(k)})\hfill & =& \frac{1}{\mu }\left[\widehat{𝒔}^{(k)}(𝒙\widehat{𝒙}^{(k)})^t\right]\hfill \end{array}$$
(45)
One can make the comparison with different neural network based algorithm.

Remark 2:
In the Gaussian hypothesis case, if we use the prior law (17), we obtain similar closed form expressions equivalent to (39)

$$\{\begin{array}{ccc}s_j=\frac{_ia_{ij}(x_i\widehat{x}_i)}{\lambda +𝒂_j^2}\hfill & & \\ a_{ij}=\frac{s_j(x_i\widehat{x}_i)}{s_j^2+\mu },\text{for }i=j\text{and}a_{ij}=\frac{s_j(x_i\widehat{x}_i)}{s_j^2\mu },\text{for }ij.\hfill & & \end{array}$$
(46)
Remark 3:
The main interest of this approach is that we can, at least in theory, account for the existance of any correlation between $`s_j`$ and $`s_k`$ or to model the temporal behavior of any source $`s_j(t)`$, for example, via a markov model. We can also account for any prior information we may have about the mixing matrix $`𝑨`$ or impose any desired structure for the separating matrix $`𝑩`$. For example, if we know that the sources are labelled in such a way that the sensor $`x_i`$ is closer to the sources $`s_i`$, $`s_{i1}`$ and $`s_{i+1}`$ than to any others, we can use it by choosing a prior probability law

$$p(𝑨)\mathrm{exp}\left[\frac{1}{\sigma _a^2}\underset{i}{}\underset{j}{}(w_{ij}a_{ij})^2\right]$$
(47)
with

$$w_{ij}=\{\begin{array}{ccc}1\hfill & if& i=j\hfill \\ 1/(2|ij+1])\hfill & if& ij\hfill \end{array}$$
(48)
or $`w_{ii}=1`$, $`w_{i,i1}=w_{i1,i}=\alpha `$ and $`(1\alpha )`$ for all the other coefficients $`w_{ij}`$ for some $`0.5<\alpha <1`$. Then the equations (39) and (40) become

$$\{\begin{array}{ccc}s_j=\frac{_ia_{ij}(x_i\widehat{x}_i)}{\lambda +𝒂_i^2}\hfill & & \\ a_{ij}=\frac{s_j(x_i\widehat{x}_i)}{w_{ij}^2(s_j^2+\mu )}\hfill & & \end{array}$$
(49)
and

$$\{\begin{array}{ccc}𝒔\hfill & =& (𝑨^t𝑨+\lambda 𝑰)^1𝑨^t𝒙\hfill \\ 𝑨\hfill & =& 𝒙𝒔^t(𝒔𝒔^t+\mu 𝑾𝑾^t)^1\hfill \end{array}$$
(50)
#### Marginal MAP estimations

Now, we consider the two other approaches of marginal MAP estimations. First we note that we can rewrite $`𝒙=𝑨𝒔`$ with $`𝑨`$ is a $`(m\times n)`$ matrix as

$$𝒙=𝑨𝒔=𝑺𝒂$$
(51)
where $`𝑺`$ a $`(m\times mn)`$ bloc Toeplitz matrix and $`𝒂`$ a vector of dimension $`mn`$ obtained by pilling up all the rows of the matrix $`𝑨`$:

$$𝑺=\left(\begin{array}{ccccc}𝒔^t& \mathrm{𝟎}& \mathrm{}& & \mathrm{𝟎}\\ \mathrm{𝟎}& 𝒔^t& \mathrm{}& & \mathrm{}\\ \mathrm{}& & & \\ \mathrm{𝟎}& \mathrm{}& & 𝒔^t& \mathrm{𝟎}\\ \mathrm{𝟎}& \mathrm{𝟎}& \mathrm{}& & 𝒔^t\end{array}\right)\text{and}𝒂=\left(\begin{array}{cccc}𝒂_1& 𝒂_2& \mathrm{}& 𝒂_m\end{array}\right)^t$$
(52)
To be able to obtain closed form expression, in the following we consider only the Gaussian case:

$`p(𝒙|𝑨,𝒔)`$ $``$ $`\mathrm{exp}\left[{\displaystyle \frac{1}{2\sigma _ϵ^2}}𝒙𝑨𝒔^2\right]`$ (53)
$`p(𝑨)`$ $``$ $`\mathrm{exp}\left[{\displaystyle \frac{1}{2\sigma _a^2}}\psi (𝑨)\right]\text{ with }\psi (𝑨)=𝑨^2`$ (54)
$`p(𝒔)`$ $``$ $`\mathrm{exp}\left[{\displaystyle \frac{1}{2\sigma _s^2}}\varphi (𝒔)\right]\text{ with }\varphi (𝒔)=𝒔^2`$ (55)
which gives

$$p(𝑨,𝒔|𝒙)\mathrm{exp}\left[\frac{1}{2\sigma _ϵ^2}J(𝑨,𝒔)\right]$$
(56)
with

$`J(𝑨,𝒔)`$ $`=`$ $`𝒙𝑨𝒔^2+\lambda \varphi (𝒔)+\mu \psi (𝑨)=𝒙𝑨𝒔^2+\lambda 𝒔^2+\mu 𝑨^2`$ (57)
$`=`$ $`𝒙𝑺𝒂^2+\lambda \varphi (𝑺)+\mu \psi (𝒂)=𝒙𝑺𝒂^2+{\displaystyle \frac{\lambda }{m}}𝑺^2+\mu 𝒂^2`$ (58)
Note that, in this case, $`J(𝑨,𝒔)`$ is a quadratic function of $`𝒔`$ for fixed $`𝑨`$ and a quadratic function of $`𝑨`$ for fixed $`𝒔`$, but it is a bilinear function of both $`𝑨`$ and $`𝒔`$. This remark means that the joint MAP solution in this case is not unique.

Note also that $`J(𝑨,𝒔)`$ can be rewritten as

$$J(𝑨,𝒔)=(𝒔\widehat{𝒔})^t\widehat{𝑷}_s^1(𝒔\widehat{𝒔})\widehat{𝒔}^t\widehat{𝑷}_s^1\widehat{𝒔}+𝒙^t𝒙+\mu 𝑨^2$$
(59)
with

$$\widehat{𝒔}=(𝑨^t𝑨+\lambda 𝑰)^1𝑨^t𝒙\text{and}\widehat{𝑷}_s=(𝑨^t𝑨+\lambda 𝑰)^1;$$

or as

$$J(𝒂,𝑺)=(𝒂\widehat{𝒂})^t\widehat{𝑷}_a^1(𝒂\widehat{𝒂})\widehat{𝒂}^t\widehat{𝑷}_a^1\widehat{𝒂}+𝒙^t𝒙+\lambda 𝑺^2$$
(60)
with

$$\widehat{𝒂}=(𝑺^t𝑺+\mu 𝑰)^1𝑺^t𝒙\text{and}\widehat{𝑷}_a=(𝑺^t𝑺+\mu 𝑰)^1.$$

With these notations, it is then easy to obtain the expression of the marginals laws $`p(𝒔|𝒙)`$ and $`p(𝑨|𝒙)`$:

$`\mathrm{ln}p(𝑨|𝒙)`$ $``$ $`\mathrm{ln}\left|\text{det}(\widehat{𝑷}_s^1)\right|J(𝑨,\widehat{𝒔})`$ (61)
$`=`$ $`\mathrm{ln}\left|\text{det}(𝑨^t𝑨+\lambda 𝑰)\right|𝒙^t(𝒙𝑨\widehat{𝒔})+\mu 𝑨^2`$
$`\mathrm{ln}p(𝒔|𝒙)`$ $``$ $`\mathrm{ln}\left|\text{det}(\widehat{𝑷}_a^1)\right|J(\widehat{𝒂},𝒔)`$ (62)
$`=`$ $`\mathrm{ln}\left|\text{det}(𝑺^t𝑺+\mu 𝑰)\right|𝒙^t(𝒙𝑺\widehat{𝒂})+\lambda 𝒔^t𝒔`$
In non Gaussian hypothesis for $`𝒔`$ and $`𝒂`$ where $`\mathrm{ln}p(𝒔)=\mathit{\varphi }(𝒔)+cte`$ and $`\mathrm{ln}p(𝒂)=𝝍(𝒂)+cte`$, we can always use the Laplace approximation of the posterior laws. Then, we can use the first lines of these two last equations by replacing $`\text{det}(\widehat{𝑷}_s^1)`$ and $`\text{det}(\widehat{𝑷}_a^1)`$ by, the jacobians of their respective log probability densities: $`_𝒔^2J(𝑨,𝒔)`$ and $`_𝒂^2J(𝒂,𝑺)`$ computed for the MAP estimates $`\widehat{𝒔}=\text{arg}\mathrm{min}_𝒔\left\{J(𝑨,𝒔)\right\}`$ and $`\widehat{𝒂}=\text{arg}\mathrm{min}_𝒂\left\{J(𝒂,𝑺)\right\}`$.

The marginal MAP solutions of $`𝑨`$ and $`𝒔`$ respectively satisfy $`\frac{\mathrm{ln}p(𝑨|𝒙)}{𝑨}=0`$ and $`\frac{\mathrm{ln}p(𝒔|𝒙)}{𝒔}=0`$. Note that, excepted the Gaussian case where these equations have analytical solutions, we need a numerical optimisation algorithm to compute the solutions. For example, the marginal MAP estimate of $`𝑨`$ can be computed by the following iterative algorithm:

$$\begin{array}{ccc}\widehat{𝑨}^{(k)}\hfill & =& \underset{𝑨}{\text{arg}\mathrm{max}}\left\{\mathrm{ln}\left|\text{det}(\widehat{𝑷}_s^{(k1)})\right|+J(𝑨,\widehat{𝒔}^{(k1)})\right\}\hfill \\ & =& \underset{𝑨}{\text{arg}\mathrm{max}}\left\{\mathrm{ln}\left|\text{det}(𝑨^t𝑨+\lambda 𝑰)\right|+𝒙^t𝑨\widehat{𝒔}^{(k1)}+\mu \psi (𝑨)\right\}\hfill \end{array}$$
(63)
where

$$\widehat{𝒔}^{(k)}=\underset{𝒔}{\text{arg}\mathrm{min}}\left\{J(\widehat{𝑨}^{(k1)},𝒔)\right\}=\underset{𝒔}{\text{arg}\mathrm{min}}\left\{𝒙\widehat{𝑨}^{(k1)}𝒔^2+\lambda \varphi (𝒔)\right\}.$$
(64)
In general, neither of these equations have explicite solutions and we have to do the optimization numerically. For example, a gradient based algorithm to compute the marginal MAP estimate of $`𝑨`$ writes:

$$\{\begin{array}{ccc}\mathrm{\Delta }\widehat{𝒔}^{(k)}\hfill & & \widehat{𝑨}^{t^{(k1)}}(𝒙\widehat{𝑨}^{(k1)}𝒔)+\lambda \mathit{\varphi }^{}(𝒔)\hfill \\ \mathrm{\Delta }\widehat{𝑨}^{(k)}\hfill & & 𝑨^t(𝑨^t𝑨+\lambda 𝑰)^1+𝒙\widehat{𝒔}^{(k1)}+\mu 𝝍^{}(𝑨)\hfill \end{array}$$
(65)
### 2.4. Spatially correlated sources

As a first extension we still assume that the sources are white, but they are spatially correlated, i.e.

$$\mathrm{log}p(𝒔_{1..T})=\underset{t}{}r(s_1(t),\mathrm{},s_N(t))$$

where $`r(s_1,\mathrm{},s_N)`$ represents the joint probability law of the sources. Then we obtain

$$\mathrm{log}p(𝑨,𝒔_{1..T}|𝒙_{1..T})=\underset{t}{}\underset{i}{}\left[q_i\left(x_i(t)y_i(t)\right)+r(s_1(t),\mathrm{},s_N(t))\right]+\mathrm{log}p(𝑨)+cte.$$
(66)
But the main difficulty here will be the modeling of these dependencies and simplification of the expression of the joint probability law $`r(s_1,\mathrm{},s_N)`$. For example, if the sources are labelled in such a way that only the immediate neighbor sources are correlated, then we can use a first order Markov model and write

$$r(s_1,\mathrm{},s_N)=\underset{j}{}r(s_j|s_{j1},s_{j+1})$$
(67)
However, the complexity of the optimization algorithms depend on the expression of the potential function $`r(s_j|s_{j1},s_{j+1})`$. A simple case is a Gaussian model where

$$\underset{j}{}r(s_j|s_{j1},s_{j+1})=\frac{1}{2\sigma _s^2}\left|2s_j(s_{j1}+s_{j+1})\right|^2=\frac{1}{2\sigma _s^2}𝑫𝒔^2$$

with $`𝑫`$ a tri-diagonal Toeplitz matrix with diagonal elements equal to 2 and off-diagonal elements equal to -1. In this case, the equations (40) become

$$\{\begin{array}{ccc}𝑨^t(𝒙𝑨𝒔)\lambda 𝑫^t𝑫𝒔\hfill & =& \mathrm{𝟎}\hfill \\ (𝒙𝑨𝒔)𝒔^t\mu 𝑨\hfill & =& \mathrm{𝟎}\hfill \end{array}\{\begin{array}{ccc}𝒔\hfill & =& (𝑨^t𝑨+\lambda 𝑫^t𝑫)^1𝑨^t𝒙\hfill \\ 𝑨\hfill & =& 𝒙𝒔^t(𝒔𝒔^t+\mu 𝑰)^1\hfill \end{array}$$
(68)
### 2.5. Spatially independent but colored sources

Here we assume that the sources are mutually independent but that they are temporally colored, i.e.

$$\mathrm{log}p(𝒔_{1..T})=\underset{j}{}r_j(s_j(1),\mathrm{},s_j(T))$$

where $`r_j(s_j(1),\mathrm{},s_j(T))`$ represents the joint probability law of the different samples of source number $`j`$. Then we obtain

$`\mathrm{log}p(𝑨,𝒔_{1..T}|𝒙_{1..T})`$ $`=`$ $`{\displaystyle \underset{t}{}}{\displaystyle \underset{i}{}}q_i\left(x_i(t)y_i(t)\right)+{\displaystyle \underset{j}{}}r_j(s_1(1),\mathrm{},s_j(T))`$ (69)
$`+{\displaystyle \underset{k}{}}{\displaystyle \underset{l}{}}a_{kl}^2+cte.`$
Here again, the main difficulty is the modelization and simplification of the expression of the joint probability laws $`r_j(s_j(1),\mathrm{},s_j(T))`$. For example, we can use a first order markov chain model and write

$$r_j(s_j(1),\mathrm{},s_j(T))=\underset{t}{}r_j(s_j(t)|s_j(t1))$$
(70)
As an example here we consider the Gaussian case (or equivalently first order AR models) :

$$r_j(s_j(1),\mathrm{},s_j(T))=\underset{t}{}\alpha _j(s_j(t)s_j(t1))^2$$
(71)
With this assumption we have

$$\{\begin{array}{ccc}\widehat{𝒔}_{1..T}^{(k)}\hfill & =& \underset{𝒔}{\text{arg}\mathrm{max}}\left\{\underset{t}{}\underset{i}{}q_i\left(x_i(t)y_i(t)\right)\underset{t}{}\underset{j}{}\alpha _j(s_j(t)s_j(t1))^2\right\}\hfill \\ \widehat{𝑨}^{(k)}\hfill & =& \underset{𝑨}{\text{arg}\mathrm{max}}\left\{\underset{t}{}\underset{i}{}q_i\left(x_i(t)y_i(t)\right)+\frac{1}{\sigma _a^2}\underset{k}{}\underset{l}{}a_{kl}^2\right\}\hfill \end{array}$$
(72)
The solution at each iteration has to satisfy

$$\{\begin{array}{ccc}\underset{t}{}\underset{i}{}a_{ij}q_i^{}\left(x_i(t)y_i(t)\right)\underset{t}{}\underset{j}{}2\alpha _j(s_j(t)s_j(t1))=0\hfill & & \\ \underset{t}{}s_jq_i^{}\left(x_i(t)y_i(t)\right)+\frac{2}{\sigma _a^2}a_{ij}=0\hfill & & \end{array}$$
(73)
For the special case of Gaussian noise we obtain

$$\{\begin{array}{ccc}\underset{t}{}\underset{i}{}a_{ij}\left(x_i(t)[𝑨𝒔]_i(t)\right)\underset{t}{}\underset{j}{}2\lambda _j(s_j(t)s_j(t1))\hfill & =& 0\hfill \\ \underset{t}{}s_j(t)\left(x_i(t)[𝑨𝒔]_i(t)\right)+\mu a_{ij}\hfill & =& 0\hfill \end{array}$$
(74)
with $`\lambda _j=\alpha _j\sigma _ϵ^2`$ and $`\mu =\frac{2\sigma _ϵ^2}{\sigma _a^2}`$.

The two algorithms of (39) and (40) in this case become:

$$\{\begin{array}{ccc}s_j(t)=\frac{\lambda _js_j(t1)+_ia_{ij}(x_i\widehat{x}_i)}{\lambda +𝒂_i^2}\hfill & & \\ a_{ij}=\frac{s_j(x_i\widehat{x}_i)}{s_j^2+\mu }\hfill & & \end{array}$$
(75)
and

$$\{\begin{array}{ccc}𝒔(t)\hfill & =& (𝑨^t𝑨+\lambda 𝑰)^1\left[\text{diag}\{\lambda _1,\mathrm{},\lambda _n\}𝒔(t1)+𝑨^t𝒙(t)\right]\hfill \\ 𝑨\hfill & =& 𝒙𝒔^t(𝒔𝒔^t+\mu 𝑰)^1\hfill \end{array}$$
(76)
Here, we conclude the presentation of the Bayesian approach to source separation. Beside some of the details of implementation, we showed that the Bayesian approach gives us the possibility to push further some of the limitations of the classical techniques in source separation. In the next section we give a few preliminary numerical results to show the performances of the proposed algorithems.

## 3. Simulation results

In the following we give a few preliminary examples of simple source separation problem to show some performances of the proposed methods. In all these examples, we used the following algorithm:

$$\{\begin{array}{ccc}𝒚\hfill & =& (𝑨^t𝑨+\lambda 𝑰)^1𝑨^t𝒙\hfill \\ 𝒔\hfill & =& 𝒈(𝒚)\hfill \\ \mathrm{\Delta }𝑨\hfill & & 𝑨^t(𝑨^t𝑨+\lambda 𝑰)^1+𝒙𝒔+\mu 𝝍^{}(𝑨)\hfill \end{array}$$
(77)
with $`\lambda =\mu =.1`$, $`N=100`$ and appropriate $`𝒈`$.

### 3.1. Example 1

Hier, we considered two sources

$$\{\begin{array}{ccc}s_1(t)\hfill & =& \mathrm{sin}(500t+10\mathrm{cos}(50t))\hfill \\ s_2(t)\hfill & =& \mathrm{sin}(300t)\hfill \end{array},t=[0:.001:.499].$$
(78)
and used the mixing matrice $`𝑨=\left(\begin{array}{cc}1& .4\\ .6& 1\end{array}\right)`$ to obtain the two set of data $`x_1(t)`$ and $`x_2(t)`$. Then we applied the algorithm given in (77). The following figures show the obtained results $`\widehat{s}_1(t)`$ and $`\widehat{s}_2(t)`$.

Fig. 1: Results of the source separation in Example 1.

Fig. 2: Results of the source separation in Example 1: phase space distribution of sources, mixed signals and separated sources

Fig. 3: Results of the source separation in Example 1: histograms of sources, mixed signals and separated sources

The two sources are well separated.

### 3.2. Example 2

Hier, we considered three sources

$$\{\begin{array}{ccc}s_1(t)\hfill & =& \mathrm{sin}(500t+10\mathrm{cos}(50t))\hfill \\ s_2(t)\hfill & =& \mathrm{sin}(300t)\hfill \\ s_3(t)\hfill & =& \text{sign}(\mathrm{cos}(120t5cos(50t)))\hfill \end{array},t=[0:.001:.499].$$
(79)
and used the mixing matrice

$$𝑨=.3\left(\begin{array}{ccc}1.& .5& .2\\ .5& 1.& .5\\ .5& .5& 1.\end{array}\right)$$

to obtain the three set of data $`x_1(t)`$, $`x_2(t)`$ and $`x_3(t)`$. The following figures show the obtained results.

Fig. 4: Results of the source separation in Example 2.

Fig. 5: Results of the source separation in Example 2: phase space distribution of sources, mixed signals and separated sources

Fig. 6: Results of the source separation in Example 2: histograms of sources, mixed signals and separated sources

Here also the three sources are well separated.

### 3.3. Example 3

Hier, we considered the two sources of the first example, but we simulated the case where there are three receivers using the mixing matrice

$$𝑨=\left(\begin{array}{cc}1.& .5\\ .5& 1.\\ .2& .5\end{array}\right)$$

to obtain the three set of data $`x_1(t)`$, $`x_2(t)`$ and $`x_3(t)`$. Then we applied again the algorithm given in (77). The following figures show the obtained results.

Fig. 7: Results of the source separation in Example 3.

Fig. 8: Results of the source separation in Example 3: phase space distribution of sources, mixed signals and separated sources

Fig. 9: Results of the source separation in Example 3: histograms of sources, mixed signals and separated sources

### 3.4. Example 4

Hier, we considered the three sources of the Exammple 2 and simulated the case where there are only two receivers using the mixing matrice

$$𝑨=\left(\begin{array}{ccc}1.& .2& 1\\ .5& 1.& .2\end{array}\right)$$

to obtain the two set of data $`x_1(t)`$ and $`x_2(t)`$. Then we applied again the algorithm given in (77). The following figures show the obtained results.

Fig. 7: Results of the source separation in Example 4.

Fig. 8: Results of the source separation in Example 4: phase space distribution of sources, mixed signals and separated sources

Fig. 9: Results of the source separation in Example 4: histograms of sources, mixed signals and separated sources

## 4. Conclusions

We investigated the use of the Bayesian estimation theory to source separation and showed that this approach has the potential to push farther the limits of the classical methods. This work is not really yet finished. We are going now to compare the performances of the proposed methods to other conventional ones on simulated and real data.


text:
# Anomalous aging phenomena caused by drift velocities

11institutetext:  Service de Physique Théorique, CEA Saclay, 91191 Gif-sur-Yvette cedex, France
ICTP, Strada Costiera 11, 34100 Trieste, Italy, and S N Bose National Centre for Basic Sciences, Block JD Sector 3, Salt Lake, Calcutta 700098, India

Nonequilibrium thermodynamics, irreversible processes

## Abstract

We demonstrate via several examples that a uniform drift velocity gives rise to anomalous aging, characterized by a specific form for the two-time correlation functions, in a variety of statistical-mechanical systems far from equilibrium. Our first example concerns the oscillatory phase observed recently in a model of competitive learning. Further examples, where the proposed theory is exact, include the voter model and the Ohta-Jasnow-Kawasaki theory for domain growth in any dimension, and a theory for the smoothing of sandpile surfaces.

Aging is a characteristic feature of systems such as glasses and spin glasses in their low-temperature phase. More generally, aging is a consequence of slow dynamics, which can be observed in a broad range of systems far from equilibrium . Aging phenomena are commonly described in terms of the two-time correlation and response functions of the relevant variables, which often exhibit scaling behavior in the two-time plane .

In this Letter we demonstrate that this aging scenario is strongly affected by the presence of a uniform drift velocity, either intrinsic to the microscopic dynamics, or imposed by external conditions. On general grounds, a uniform velocity field is expected to be relevant for long-time dynamics. Consider for simplicity a diffusive relaxational dynamics, characterized by a diffusion constant $`D`$. To a non-zero velocity $`V`$ naturally corresponds a finite relaxation time, $`D/V^2`$, beyond which the dynamical evolution can be qualitatively different from that with $`V=0`$.

Non-equilibrium phenomena in the presence of velocity fields, either uniform (drift) or non-uniform (shear), and more generally in a range of situations with a preferred direction, have been described in recent works. For instance, it has been shown that a shear flow induces a strong anisotropy in coarsening dynamics . The influence of drift on phase-ordering dynamics has also been investigated in Ref. . The peculiar scaling properties of the flocking model of Ref.  are yet another consequence of directionality in non-equilibrium systems.

Here we show that a variety of models exhibit a novel type of non-equilibrium dynamics in the presence of a uniform drift velocity $`V`$, characterized by the specific scaling form (3) for the two-time correlation function. We refer to this phenomenon as anomalous aging. Uniform drift velocities cannot here be transformed away naively, by invoking Galilean invariance.

We first investigate in detail the oscillatory phase described recently in the context of a model of competitive learning . The anomalous aging behavior of this specific model will be explained at a quantitative level by the combined effects of drift and diffusion. Anomalous aging will then be shown to be a common feature of a range of non-equilibrium models in the presence of a drift, including the voter model, the OJK theory for domain growth, and the smoothing of sandpile surfaces .

We begin with a brief description of our model of competitive learning . Two types of individuals, considered as quick and slow learners, are respectively classified by their efficiencies (coded as Ising variables $`\eta _𝐱=\pm 1`$), and occupy the sites of a regular lattice. The updating dynamics of the efficiencies involves two steps. First, a majority rule favors the conversion of every individual to the dominant type among its neighbors. In the next step, there is a majority rule involving the outcomes of individuals in some stochastic game, favoring the conversion of every individual to the type which has had the most favorable outcomes among its neighbors at the last iteration. Our model has two variants, the ‘interfacial’ and the ‘cooperative’, which differ in the precise formulation of the second step.

We consider the situation where both types of individuals, while distinguishable, are equally slow learners, so that no type wins a priori. We refer to this case as coexistence, in analogy with ferromagnetic models. Several kinds of ordered and disordered phases, i.e., types of non-equilibrium steady states, have been found in the phase diagram of both variants of the model ; most phase transitions were found to belong either to the Ising or to the voter universality classes. However, an unusual type of coarsening was found in the cooperative version of the model, where the order parameter underwent more or less regular oscillations in time, hence the name oscillatory phase.

As the model of competitive learning is entirely defined by dynamical rules, without reference to a Hamiltonian (and therefore without the concept of detailed balance), details of the dynamics can influence its steady states. It turns out that the updating scheme, including the order in which sites are updated (sequential or parallel, ordered or random), does affect the phase diagram. From now on, we focus on the cooperative model at coexistence on the square lattice, characterized by a single parameter, $`0p1`$, which represents an individual’s probability of obtaining a successful outcome. In Ref.  ordered sequential dynamics were used, and an oscillatory phase was found for $`p<p_o`$ $`(p_o0.136)`$. On the other hand, parallel dynamics leads to a ferromagnetically ordered phase for $`p<p_c`$ $`(p_c0.259)`$, while random sequential dynamics leads to a disordered phase with very weak correlations between efficiencies.

The special feature of ordered sequential dynamics that is responsible for the oscillatory phase turns out to be its directionality. Consider a square sample with sites $`𝐱=(m,n)`$ with $`m,n=1,\mathrm{},L`$. For ordered sequential dynamics, neglecting boundary effects, the next site updated after $`𝐱_t=(m,n)`$ is $`𝐱_{t+\delta t}=(m,n+1)`$, with $`\delta t=1/L^2`$. We define the directionality of this updating scheme as $`𝐝=𝐱_{t+\delta t}𝐱_t`$. We have $`𝐝=(0,1)`$ for ordered sequential dynamics, while $`𝐝=\mathrm{𝟎}`$ for parallel dynamics or random sequential dynamics. Other updating schemes have been considered, including alternating sequential dynamics: the sample is first swept in a regular way, such as e.g. $`𝐱_{t+\delta t}=(m,n+1)`$ if $`𝐱_t=(m,n)`$, and then in a different regular way, such as $`𝐱_{t+\delta t}=(m,n1)`$, or $`𝐱_{t+\delta t}=(m+1,n)`$. We have $`𝐝=\mathrm{𝟎}`$ in the first example, and $`𝐝=\frac{1}{2}(1,1)\mathrm{𝟎}`$ in the second. The model exhibits an oscillatory phase at low $`p`$ if, and only if, the updating scheme has a net directionality $`𝐝\mathrm{𝟎}`$.

In the following, we restrict the analysis to the prototypical example of ordered sequential dynamics. Spin (efficiency) patterns are advected with a constant drift velocity $`𝐕`$, due to the directionality $`𝐝`$ of the dynamics. The direction of $`𝐕`$ is given by the arrow of $`𝐝`$, while its strength $`V`$ depends e.g. on the parameter $`p`$. If periodic boundary conditions are used, spin clusters wrap around the sample, so that the magnetization oscillates more or less regularly, with period $`T_1=2L/V`$, which defines a time scale associated with the drift velocity. This phenomenon was depicted in Figure 11 of Ref. . A correlation analysis of a long magnetization record allows for the accurate determination of $`V`$. We find that $`V`$ decreases roughly linearly with $`p`$, from $`V1.14`$ at $`p=0`$ to $`V0.97`$ at $`p=p_o0.136`$. The drift velocity thus remains of order unity throughout the oscillatory phase, including the critical point $`p_o`$, where it does not seem to exhibit any kind of critical behavior.

Next, a coarsening phenomenon takes place if the system starts from a random initial state with uncorrelated efficiencies. The commonly accepted picture of phase-ordering dynamics with no conservation law  involves diffusive behavior, hence the emergence of a second time scale, $`T_2=L^2/D`$, with the diffusion coefficient $`D`$ being a microscopic constant. The competition between these two time scales will turn out to be central to anomalous aging \[see paragraph below Eq. (5)\].

Along the lines of investigations of aging properties in glassy systems, we next measured the two-time correlation function of the efficiency of a given individual, $`C(t,s)=\eta _𝐱(s)\eta _𝐱(t)`$, with $`0s`$ (waiting time) $`t=s+\tau `$ (observation time). We have chosen to work at $`p=0`$ for definiteness. For each value of the waiting time $`s`$, data were averaged over more than 1,000 samples of size $`300\times 300`$. The particular case $`s=0`$ corresponds to the correlation of the efficiencies at time $`t`$ with their random initial values. Figures 1 and 2 respectively show plots of $`C(t,0)`$ and $`C(t,s)`$ for $`s>0`$. Figure 1 demonstrates that $`C(t,0)`$ decays exponentially, just as in any equilibrium process, with a microscopic relaxation time $`\tau _01/0.342.9`$, while Figure 2 clearly shows non-stationarity or aging effects: correlations decay more and more slowly when the waiting time is increased.

In order to elucidate these observations, and especially to understand better the role of the drift velocity $`V`$ on the non-equilibrium dynamics of the model, we propose the following phenomenological linear theory meant to incorporate the combined effects of drift and diffusion, in a natural fashion and at a minimal level of sophistication. The efficiencies $`\eta _𝐱`$ are replaced by a continuous field $`\varphi (𝐱,t)`$. The equal-time correlation function $`C(𝐱,t)=\varphi (𝐱,t)\varphi (\mathrm{𝟎},t)`$ is assumed to obey the diffusion equation

$$C(𝐱,t)/t=D_1^2C(𝐱,t),$$
(1)
as the drift velocity is expected to have no effect on one-time observables. A disordered initial state corresponds to $`C(𝐱,0)=\delta ^{(d)}(𝐱)`$ at $`t=0`$, in dimension $`d`$. The two-time correlation function $`C(𝐱,t,s)=\varphi (𝐱,t)\varphi (\mathrm{𝟎},s)`$ is assumed to be affected by both the drift and the diffusion mechanisms:

$$C(𝐱,t,s)/t=D_2^2C(𝐱,t,s)𝐕C(𝐱,t,s)$$
(2)
$`(ts)`$, with initial condition $`C(𝐱,s,s)=C(𝐱,s)`$ at $`t=s`$. The diffusion coefficients $`D_1`$ and $`D_2`$ are possibly different. Convoluting the Green’s functions of both equations, and neglecting pre-exponential factors, we obtain

$$C(t=s+\tau ,s)\mathrm{exp}\left(\frac{V^2\tau ^2}{4(D_1s+D_2\tau )}\right).$$
(3)
This analytical result is in excellent agreement with our results above. First, Eq. (3) predicts that $`C(t,0)`$ decays exponentially, with a relaxation time $`\tau _0=4D_2/V^2`$, as seen in Figure 1, where the values of $`V`$ and $`\tau _0`$ yield $`D_20.94`$. Second, in order to check the specific form of Eq. (3), the data for $`\mathrm{ln}C(t,s)`$ have been plotted as a function of $`\tau ^2/(\tau +\alpha s)`$. A convincing data collapse and linear behavior are observed (Figure 3) for the optimal value $`\alpha 0.55`$, yielding $`D_1=\alpha D_20.52`$. Moreover, Eq. (3) predicts that the slopes of Figures 1 and 3 are equal; the measured slopes of the least-square fits, $`0.33`$ and $`0.34`$, coincide within errors, in agreement with the above.

Our result (3) is the centerpiece of anomalous aging. The correlation function therein is clearly aging, in the general sense that it is not stationary, and decays more and more slowly as the waiting time $`s`$ is increased. The expression (3) does not, however, exhibit any simple scaling form as a function of the time ratio $`\tau /s`$, since it involves the intrinsic time scale $`\tau _0=4D_2/V^2`$, emblematic of the competition between drift and diffusion, as expected from the discussion in the beginning of this Letter.

For $`\tau s`$, the correlation function in Eq. (3) has a Gaussian scaling form

$$C(t,s)\mathrm{exp}\left(\tau ^2/(\tau _1s)\right)(\tau =tss),$$
(4)
with $`\tau _1=4D_1/V^2=\tau _0D_1/D_2`$. The correlation is characterized by an effective relaxation time

$$\tau _{\mathrm{eff}}(s)(\tau _1s)^{1/2},$$
(5)
which diverges, albeit anomalously slowly, with the waiting time $`s`$.

Eq. (5) can be recast as $`(V\tau _{\mathrm{eff}})^24D_1s`$, to be put in perspective with the similar quadratic relationship $`(VT_1)^2=4DT_2=4L^2`$ which holds between both time scales introduced above. The stationary or time-translationally-invariant (TTI) regime usually observed for $`\tau s`$, where the correlation function only depends on $`\tau `$, is absent in the present case, having already died away before our continuum description applied. The TTI regime is also too small to be visible on Figure 2.

For $`\tau s`$, the result (3) only depends on $`\tau `$, and falls off exponentially, as

$$C(t,s)\mathrm{exp}(\tau /\tau _0)(\tau =tss),$$
(6)
just as $`C(t,0)`$. This paradoxical emergence of a translationally-invariant regime for long time separations is one of the reasons we refer to this aging as anomalous. The crossover between both kinds of behavior (4) and (6) takes place for $`\tau s`$, where the correlation is already exponentially small, of order $`\mathrm{exp}(s/\tau _0)`$.

In agreement with general expectations on aging phenomena , the result (3) can be recast as $`C(t,s)F[h(t)/h(s)]`$, with $`h(t)=\mathrm{exp}\left(2(t/\tau _1)^{1/2}\right)`$ and $`F(x)=\mathrm{exp}\left((\mathrm{ln}x)^2\right)`$. Anomalous aging thus appears as a special case of sub-aging, using the terminology of Ref. , as the effective age function $`h(t)`$ grows much more rapidly than real time $`t`$. The effective relaxation time (5) is recovered as $`\tau _{\mathrm{eff}}(t)h(t)/(\mathrm{d}h(t)/\mathrm{d}t)`$. A similar behavior, with the same age function $`h(t)`$, has been found  in the spherical model with random asymmetric couplings at zero temperature. The ferromagnetic spherical model with a conserved order parameter  provides another interesting instance of anomalous aging, involving a scaling function similar to the above $`F(x)`$.

Since Eq. (3) is a simple consequence of the combined effects of drift and diffusion, it might be expected to apply to a variety of situations. Several examples are discussed below, where anomalous aging applies exactly in the regime of long times.

First, our linear theory is exact for the voter model. The dynamics of this model  amounts to saying that a given site, say $`𝐱`$, takes the opinion of any of its neighbors, say $`𝐲=𝐱+𝐞`$, with the lattice vector $`𝐞`$ being chosen at random with some weight $`w(𝐞)`$. If all vectors $`𝐞`$ have equal weights ($`w(𝐞)=1/z`$ for all $`𝐞`$, with $`z`$ the coordination number), the usual voter model is recovered, with its well-known diffusive dynamics . In the general (asymmetric) case where some neighbors are preferred over others, patterns of opinions are advected with a velocity $`𝐕=_𝐞𝐞w(𝐞)`$. We have shown that the equal-time and two-time correlation functions obey Eqs. (1) and (2) in the continuum limit, with the identity

$$D_1=2D_2,$$
(7)
so that anomalous aging applies exactly to the voter model.

Second, the Ohta-Jasnow-Kawasaki (OJK) theory , a linear approximation commonly used in the study of phase-ordering dynamics , also yields, in the presence of drift, the anomalous aging result (3). The main difference between the OJK theory and our linear theory is that linear equations are assumed to be obeyed by the field $`\varphi (𝐱,t)`$ itself in the OJK theory, and by its correlations in our approach. We have checked that, despite this difference, anomalous aging (Eq. (3)) applies exactly to the OJK theory with drift, again with the simultaneous validity of Eq. (7).

Our last example relates to the smoothing of sandpile surfaces  represented in earlier work by noisy coupled equations for the fields $`h(𝐱,t)`$ and $`\rho (𝐱,t)`$. Here $`h(𝐱,t)`$ represents the local height fluctuations of the sandpile surface ( caused by clusters of ‘stuck’ grains) defined with respect to a mean surface, while $`\rho (𝐱,t)`$ represents the local density of mobile grains rolling down the clusters . The simplest example of these coupled equations is fully linear:

$$\{\begin{array}{c}h/t=D_h^2h+𝐜h+\eta (𝐱,t),\hfill \\ \rho /t=D_\rho ^2\rho 𝐜h,\hfill \end{array}$$
(8)
where $`\eta (𝐱,t)`$ is Gaussian white noise, such that $`\eta (𝐱,t)\eta (𝐱^{},t^{})=\mathrm{\Delta }^2\delta ^{(d)}(𝐱𝐱^{})\delta (tt^{})`$. The equation for $`h(𝐱,t)`$ is nothing but the Edwards-Wilkinson (EW) equation  with a flow (drift) term added; in the context of sandpiles, this represents the effect of tilt .

In Ref.  it was found that while $`S_h(𝐤,t)`$, a one-time quantity, shows only the smoothing exponents of the EW equation (i.e., the flow term can indeed be Galilean transformed away in this case), $`S_h(𝐤,\omega )`$, a quantity which involves an integration over many times, shows a clear crossover from EW exponents to a smoothing state with $`\alpha =\beta =0`$, when the drift velocity is large enough. Evaluating the two-time correlation function of the $`h`$-field, we find

$$C_h(t,s)=h(\mathrm{𝟎},t)h(\mathrm{𝟎},s)=\mathrm{\Delta }^2_0^sdu\left(4\pi D_h(t+s2u)\right)^{d/2}\mathrm{exp}\left(\frac{c^2(ts)^2}{4D_h(t+s2u)}\right).$$
(9)
In the scaling regime, where both $`s`$ and $`ts=\tau `$ are large compared to the intrinsic time scale $`\tau _0=4D_h/c^2`$, the integral is dominated by microscopic times $`u\tau _0`$. We thus recover the anomalous aging result (3), with $`D_1=2D_h`$, $`D_2=D_h`$, so that Eq. (7) again holds.

A corollary of this result is that anomalous aging (or, equivalently, the anomalous smoothing referred to in Ref. ) will manifest itself at experimentally observable times only for large $`c`$ – this, in the sandpile context, would correspond to unstable packings of grains susceptible to rolling under tilt . It would be interesting to see if measurements on such a surface indicate that it is composed of large domains characterized by EW roughening exponents stacked end on end. Physically, this is one possible realization of the anomalous smoothing scenario – inter-domain roughness measurements would here show no significant variation in roughness, while intra-domain measurements would show EW roughening exponents.

To sum up, we have demonstrated that a uniform drift velocity is relevant for the long-time dynamics of systems far from equilibrium. We have shown that the prediction (3) of a simple linear theory, referred to as anomalous aging, provides an accurate or even exact description of the non-equilibrium dynamics of a wide variety of systems, including the oscillatory phase observed recently in a model of competitive learning, the voter model, the Ohta-Jasnow-Kawasaki theory for domain growth, and a theory for the smoothing of sandpile surfaces.

###### Acknowledgements.

It is a pleasure for us to thank L. Berthier, L. Cugliandolo, S. Franz, and M. Laessig for fruitful discussions.