Sixin Zhang

Stéphane Mallat

Maximum Entropy Models from Phase Harmonic Covariances

Keywords: covariance, stationary process, phase, Fourier, wavelets, turbulence

The covariance of a stationary process X is diagonalized by a Fourier transform. It does not take into account the complex Fourier phase and defines Gaussian maximum entropy models. We introduce a general family of phase harmonic covariance moments, which rely on complex phases to capture non-Gaussian properties. They are defined as the covariance of H(LX), where L is a complex linear operator and H is a non-linear phase harmonic operator which multiplies the phase of each complex coefficient by integers. The operator H can also be calculated from rectifiers, which relates H(LX) to neural network coefficients. If L is a Fourier transform then the covariance is a sparse matrix whose non-zero off-diagonal coefficients capture dependencies between frequencies. These coefficients have similarities with high order moment, but smaller statistical variabilities because H is Lipschitz. If L is a complex wavelet transform then off-diagonal coefficients reveal dependencies across scales, which specify the geometry of local coherent structures. We introduce maximum entropy models conditioned by these wavelet phase harmonic covariances. The precision of these models is numerically evaluated to synthesize images of turbulent flows and other stationary processes.

Introduction

Many phenomena in physics, finance, signal processing and image analysis can be modeled as realizations of a non-Gaussian stationary process X. We often observe a single realization of dimension d from which the model must be estimated. Models of stochastic processes may be defined as maximum entropy distributions conditioned by a predefined family of moments [START_REF] Jaynes | Information theory and statistical mechanics[END_REF][START_REF] Cover | Elements of Information Theory[END_REF]. To estimate accurately these moments from a able image texture syntheses have been obtained in [START_REF] Ustyuzhaninov | What does it take to generate natural textures[END_REF], by computing the covariance of output coefficients of an one-layer convolutional neural network with a rectifier. Section 3.2 proves that a non-linear phase harmonic operator H can be calculated with a rectifier followed by a Fourier transform along a phase variable. These neural network covariances are therefore equivalent to phase harmonic covariances, but for a convolutional operator L which is not a wavelet transform. The analysis of phase harmonic operators and the properties of covariance matrices gives a general mathematical framework to understand these algorithms. It guides the choice of the linear transform L, of the diagonal and off-diagonal covariance moments, to model specific classes of random processes.

Section 4 reviews the properties of maximum entropy models conditioned by empirical estimators of covariance moments. Sampling a maximum entropy distribution requires to use expensive algorithms that iterate over Gibbs samplers [START_REF] Lustig | Microcanonical monte carlo simulation of thermodynamic properties[END_REF], which is not feasible over large size images. We thus rather use microcanonical models studied in [START_REF] Bruna | Multiscale sparse microcanonical models[END_REF]. Samples of microcanonical models are estimated by transporting an initial maximum entropy measure with a gradient descent which adjusts its moments. Section 5 defines low-dimensional foveal models conditioned by a small subset of wavelet phase harmonics covariances. These moments capture local dependencies within and across scales. We evaluate the numerical precision of foveal wavelet phase harmonic models, to approximate two-dimensional non-Gaussian processes including images of turbulent flows and multifractals. All calculations can be reproduced by a Python software available at https://github.com/kymatio/phaseharmonics.

Notations: We write z * the complex conjugate of z ∈ C, Y * is the complex transpose of a matrix Y and A * is the adjoint of an operator A. The covariance of two random variables A and B is written Cov(A, B) = E(A B * ) -E(A) E(B) * . An inner product is written x, y = x(u) y(u) du in L 2 (R d ) and x, y = u x(u) y(u) in R d . The cardinal of a set S is |S|. The norm of x in L 2 (R d ) or R d is written x = x, x 1/2 .

Fourier Phase Harmonic Representation

Section 2.1 introduces general properties of maximum entropy models conditioned by covariance coefficients of a non-linear representation R(X) over a graph. Section 2.2 then shows that the Fourier coefficients of a stationary process are uncorrelated because of phase fluctuations. It motivates the use of higher order moments to define a representation which can cancel random phase fluctuations and capture dependence across frequencies. Section 2.3 reviews the properties of phase harmonic operators H introduced in [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF]. If L is a Fourier transform, Section 2.4 shows that R(X) = H(LX) is a bi-Lipschitz representation, whose covariance captures non-Gaussian properties through dependencies across frequencies, similarly to high order moments.

Maximum Entropy Covariance Graph Models

We introduce maximum entropy models of a stationary random vector X, conditioned by covariance coefficients of a representation R(X) which is also a random vector:

K R = Cov(R(X), R(X)) = E (R(X) -M R )(R(X) -M R ) *
with M R = E(R(X)). We suppose that X(u) is defined for u in a cube Λ d ⊂ Z r with d grid points, for example in [1, d 1/r ] r . If r = 2 then each realization is an image of d pixels. This section reviews general properties of maximum entropy models defined on a graph, which take advantage of known symmetries of the distribution of X.

Maximum entropy on a graph with symmetries Let us write R(X) = {R v (X)} v∈V , where R v (X) ∈ C and V is a finite set of vertices in a graph that we shall define. The representation R(X) is a complex-valued vector in C |V | .

We introduce a maximum entropy model conditioned by the covariance between vertices

K R (v, v ) = Cov(R v (X), R v (X)),
where v is in a neighborhood N v ⊂ V of v. We suppose that v ∈ N v . If v ∈ N v then v ∈ N v so it defines a reflexive and symmetric undirected graph (V, E), where the set of edges relates all neighbors

E = {(v, v ) : v ∈ V , v ∈ N v }. The edge weights are the covariances K R (v, v ) over E. In statistical physics, (R v (x)-M R (v))(R v (x)-M R (v ))
* is called the interaction potential of (v, v ) ∈ E, where the mean M R is considered here as a predefined constant vector.

Let p be the probability density of X. We can reduce the constraints of the maximum entropy model if we also know a finite group G of symmetries of p. We consider linear unitary symmetries g from R d to R d , which satisfy p(g.x) = p(x) for all x ∈ R d . The stationarity of X in Λ d means that p is invariant to periodic translations, that we write g.x(u) = x(u -g), so G includes all translations. If p is invariant to rotations then G also includes rotations. Since p is invariant to the action of G, its covariance is also invariant to the action of any g ∈ G

Cov(R v (g.X), R v (g.X)) = Cov(R v (X), R v (X)).
Let |G| be the total number of symmetries. These |G| conditions will be included in the maximum entropy model.

The entropy of a density p on R d is H(p) = -p(x) log p(x) dx.

A maximum entropy macrocanonical model X conditioned by the covariance K R over E and by the symmetry group G has a probability density p which maximizes H(p) and satisfies covariance moment conditions for all (v, v , g)

∈ E × G (R v (g.x) -M R (v)) (R v (g.x) -M R (v )) * p(x) dx = K R (v, v ). (1)
If there exist Lagrange multipliers

β v,v ∈ C such that p(x) = Z -1 exp - 1 |G| (v,v ,g)∈E×G β v,v (R v (g.x)-M R (v)) (R v (g.x)-M R (v )) * (2)
is a probability distribution which satisfies each equality constraint [START_REF] Jaynes | Information theory and statistical mechanics[END_REF], then it is the unique solution to the maximum-entropy problem [START_REF] Cover | Elements of Information Theory[END_REF]. The sum in the exponential is the Gibbs energy, and Z is the partition function.

Since K R (v, v ) = K * R (v , v) we have β v,v = β * v ,v .
We verify that for all g ∈ G, p(g.x) = p(x) so G is also a group of symmetries of p. If G includes all translations then p is stationary.

In the particular case where all R v (x) are linear operators then the Gibbs energy is bilinear and p(x) is thus a Gaussian distribution. Appendix A shows that the Lagrange multipliers can then be computed efficiently [START_REF] Cover | Elements of Information Theory[END_REF]. If some of the R v (x) are non-linear then computing the Lagrange multipliers is computationally very expensive when the dimension d of X and the total number |E| of moments is large.

Covariance estimation from symmetries

The covariance estimation from a single realization x of X is calculated with an average over the symmetries of G. The orbit of the action of G on x is the set of {g.x} g∈G . The empirical estimation of

M R (v) = E(R v (X)
) is computed as an empirical average over this orbit

M R = 1 |G| g∈G R(g.x). (3)
Similarly, K R is estimated with an average on the same orbit

K Rx = 1 |G| g∈G R(g.x) -M R R(g.x) -M R * . (4)
Summing over all transformations of x by g ∈ G is called "data augmentation" in machine learning. The larger |G| the more accurate the estimation. The maximum entropy model is estimated from a single realization x by replacing the mean M R and covariance K R in (1) by their estimation M R and K Rx .

The existence of symmetries also reduces the number of covariances that must to be estimated. We define E G as a minimum set of edges (v, v ) such that for any

(v, v ) ∈ E there exists (v 1 , v 1 , g) ∈ E G × G such that R v (x) = R v 1 (g.x) and R v (x) = R v 1 (g.x). Since p is invariant to the action of G, its covariance K R is also invariant: K R (v, v ) = Cov(R v 1 (g.X), R v 1 (g.X)) = K R (v 1 , v 1 ).
It is therefore sufficient to estimate covariance coefficients indexed by E G to specify all covariances indexed by E. The set E G is interpreted as a set of sufficient statistics. Given a realization x of dimension d, to control the overall covariance estimation errors we must insure that |E G | d.

Bi-Lipschitz continuity

The estimation of covariance coefficients may have a large variance in the presence of rare outliers. These outliers induce a large variability in the empirical sum (4) depending upon the realization x of X. To avoid amplifying these outliers, we impose that R is bi-Lipschitz, so that the variations of R(X) are of the same order as the variations of X. It also insures that R is an invertible operator. The representation R is bi-Lipschitz if there exists A R > 0 and B R such that for all (x, x

) ∈ R 2d A R x -x 2 ≤ R(x) -R(x ) 2 ≤ B R x -x 2 .
(5)

For any random vector Z, we write σ 2 (Z) = E( Z -E(Z) 2 ), which is the trace of its covariance. The following proposition proves that the variance of R(X) and X have the same order of magnitude.

Proposition 2.1. If R is bi-Lipschitz then A R σ 2 (X) ≤ σ 2 (R(X)) ≤ B R σ 2 (X). ( 6
)
Proof: If X and X are two independent random vectors having same probability distribution then the bi-Lipschitz bounds (5) of R imply that

A R E( X -X 2 ) ≤ E( R(X) -R(X ) 2 ) ≤ B H E( X -X 2 ). (7)
If Y and Y are two i.i.d random variables then we verify that

E(|Y -Y | 2 ) = 2 E(|Y -E(Y )| 2 ).
where the first expected value is taken relatively to the joint distribution of Y and

Y . It results that E( X -E(X) 2 ) = 2 -1 E( X -X 2 ) and E( R(X) -E(R(X)) 2 ) = 2 -1 E( R(X) -R(X )) 2 )
. Inserting these equalities in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] proves [START_REF] Grossmann | Reading and Understanding Continuous Wavelet Transforms[END_REF].

Maximum entropy reduction with sparsity Zhu, Wu and Mumford [START_REF] Zhu | Minimax entropy principle and its application to texture modeling[END_REF] have proposed to optimize maximum entropy parameterized models by minimizing the resulting maximum entropy. Indeed, if X has a probability density p then model error can be evaluated with a Kullback-Leibler divergence

D KL (p||p) = p(x) log p(x) p(x) dx.
We verify that p(x) log p(x) dx = p(x) log p(x) dx by inserting (2) and by using the equalities [START_REF] Jaynes | Information theory and statistical mechanics[END_REF]. Since H(p) = p(x) log p(x) dx, it results that the Kullback-Leibler divergence is equal to the excess of entropy of the maximum entropy model:

D KL (p||p) = H(p) -H(p) ≥ 0. ( 8
)
We thus reduce the model error by reducing the maximum entropy H(p) to obtain a lower entropy model. The minimum H(p) = H(p) is reached if and only if p = p almost-surely.

In our case, the model depends upon the choice of R and of the edges E ⊂ V 2 of the covariance graph. Optimizing R by calculating H(p) is computationally too expensive. However, we explain below that we can minimize an upper bound of H(p) by finding a representation such that R(X) is as sparse as possible, which gives a partial control on the maximum entropy. Let X be a maximum entropy model of density p. If R has a stable inverse, which we guarantee with the bi-Lipschitz condition, then an upper bound of H(p) can be computed from the entropy of each marginal R v ( X). Such an entropy is small if R v ( X) is sparse, because it then has a narrow probability density centered at 0. To impose this sparsity we must find a representation R such that R v (X) is sparse. This sparsity must also be captured by diagonal covariance coefficients K R (v, v) so that the maximum entropy model X, conditioned by these moments retains this sparsity. This is the case for Fourier and wavelet phase harmonic representations studied in Sections 2.4 and 3.2. For non-Gaussian processes, it amounts to represent "coherent structures" with as few non-zero coefficients as possible.

Increasing the number |E| of moment conditions can further reduce H(p) but it also increases the statistical error when estimating these moments from a single realization of X. The choice of E is thus a trade-off between the model error (bias) and the estimation error (variance).

Fourier Phase and High Order Moments

The covariance of a stationary random vector is diagonalized by the discrete Fourier transform, because of random phase fluctuations. This suggests defining a covariance representation in a Fourier basis. For non-Gaussian processes, the dependence of Fourier coefficients is partly captured by high order moments which cancel random phase fluctuations.

For R = Id, since X(u) for u ∈ Λ d is stationary, K R (u, u ) = Cov(X(u), X(u )) only depends on u -u . We write |u| the norm of u ∈ Λ d and u .u the inner product. A low-dimensional maximum entropy model is constructed on V = Λ d by restricting covariances to neighborhoods of fixed radius:

N v = {v : |v -v | ≤ c}.
The radius c is chosen to be large enough to capture long range correlations. Since X is stationary, the symmetry group includes translations and we can thus define E G by setting v = 0. The maximum entropy model then defines a Gauss-Markov random vector [START_REF] Cover | Elements of Information Theory[END_REF] specified by covariances in a small neighborhood. This model does not capture non-Gaussian properties. If Λ d is an r-dimensional grid, then the model size is |E G | = O(c r ). It may be large if the process has long-range spatial correlations.

To better understand how to capture non-Gaussian properties, we study these covariance coefficients in a Fourier basis. We write x = F u x the discrete Fourier transform of x over Λ d :

x(ω) = u∈Λ d x(u) e -iω.u for ω = 2πd -1/r m with m ∈ Λ d . (9)
The Fourier representation

R = F u is indexed by ω ∈ V = 2πd -1/r Λ d . The covariance for (ω, ω ) ∈ V 2 is K R (ω, ω ) = Cov( X(ω), X(ω )).
If ω = 0 then E( X(ω)) = 0. If ω = ω then K R (ω, ω ) = 0 because of random phase fluctuations. Indeed, translating X(u) by any τ ∈ G q multiplies X(ω) by e -iτ.ω . Since X is stationary, this translation is a symmetry which does not modify Cov( X(ω), X(ω )). It results that Cov( X(ω), X(ω )) = e i(ω-ω ).τ Cov( X(ω), X(ω )). [START_REF] Bruna | Multiscale sparse microcanonical models[END_REF] Since this is true for any τ ∈ Λ d it implies that Cov( X(ω), X(ω

)) = 0 if ω = ω . (11)
Since the discrete Fourier transform is periodic beyond [0, 2π] r , any equality or nonequality between frequencies must be understood modulo 2π along the r directions. If X is Gaussian then X is also Gaussian so this non-correlation implies that X(ω) and X(ω ) are independent. However, if X is non-Gaussian then X(ω) and X(ω ) are typically not independent.

Phase cancellation with higher order moments To capture the dependencies of Fourier coefficients across frequencies, one can use high order moments [START_REF] Cramér | Mathematical Methods of Statistics (PMS-9)[END_REF]. For any k ∈ Z, similarly to [START_REF] Bruna | Multiscale sparse microcanonical models[END_REF], translating

X by τ ∈ Λ d yields Cov( X(ω) k , X(ω ) k ) = e i(kω-k ω ).τ Cov( X(ω) k , X(ω ) k ). (12)
A high-order Fourier representation

R v (X) = X(ω) k is indexed by v = (ω, k) ∈ V with 0 ≤ k ≤ k max . It results from (12) that K R (v, v ) = Cov( X(ω) k , X(ω ) k ) = 0 if kω = k ω . (13)
If kω = k ω and X is not Gaussian then this covariance is typically non-zero because the phase variations of X(ω ) k * cancel the phase variations of X(ω) k . This is also the key idea behind the use of bi-spectrum moments [START_REF] Rao | An introduction to bispectral analysis and bilinear time series models[END_REF]. By adjusting (k, k ) these moments provide some dependency information between X(ω) and X(ω ) for ω = ω . For example, let X(u) = x(u -S) be a random shift vector, where x(u) is a fixed signal supported in Λ d and S is a random periodic shift which is uniformly distributed in Λ d . It is a stationary process whose Fourier coefficients have a random phase:

X(ω) = x(ω) e -iSω . In this case, if kω = k ω then Cov( X(ω) k , X(ω ) k ) = x(ω) k x(ω ) * k .
It is non-zero at frequencies where x does not vanish. This shows that covariances of the high order Fourier exponents R(X) can capture the dependence of Fourier coefficients at different frequencies. However, this high order Fourier representation is not Lipschitz. Indeed, when k > 1, the exponent k amplifies the variability of each random variables X(ω), so estimators of covariance coefficients have a large variance.

Phase Harmonic Operator

We review the properties of phase harmonic operators H introduced in [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF], which replace high order exponents by an exponent on the phase only. These operators are bi-Lipschitz and can be computed by applying a rectifier followed by a Fourier transform on a phase variable. If L is a Fourier tranform, next section shows that H(LX) can have non-zero covariance coefficients across frequencies because of relative phase cancellations.

If z = |z|e iϕ(z) ∈ C then z k = |z| k e ikϕ(z) . If k > 1 then |z| k amplifies large values of |z| and hence the variance of a random variable z. A phase harmonic computes a power k ∈ Z of the phase only:

[z] k = |z| e ikϕ(z) .
It preserves the modulus:

|[z] k | = |z|. A phase harmonic operator computes all such phase harmonics H(z) = { h(k) [z] k } k∈Z , (14)
with finite energy weights

h 2 = k∈Z | h(k)| 2 < ∞.
The harmonic weights h(k) amplify or eliminate different phase harmonics.

Phase windows and phase Fourier transform

We show that H is the Fourier transform along a phase variable of a non-linear phase windowing operator, which can be computed with rectifiers used in neural networks.

Let us write h = F α (h) the Fourier transform along phases:

h(k) = 1 2π 2π 0 h(α) e -ikα dα. (15)
Observe that F α h(ϕ(z) + α) = { ĥ(k) e ikϕ(z) } k∈Z . It results that

H(z) = { h(k) [z] k } k∈Z = F α H(z) with H(z) = {|z| h(ϕ(z) + α)} α∈[0,2π] . (16)
The operator H translates the phase ϕ(z) of z by α ∈ [0, 2π] and applies a 2π periodic phase windowing h(α).

A rectifier ρ(a) = max(a, 0) is an important example of non-linearity which acts as a phase windowing when applied on the real part of z. Indeed

ρ(Real(z)) = |z| ρ(cos ϕ(z)) , so {ρ(Real(e iα z))} α∈[0,2π] = H(z) with h(α) = ρ(cos α). (17)
The rectifier phase window ρ(cos α) is positive and supported in [-π/2, π/2]. The corresponding harmonic weights are computed in [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF] with the Fourier integral [START_REF] Unser | A unifying parametric framework for 2d steerable wavelet transforms[END_REF]:

h(k) =      (-1) k/2+1 π(k 2 -1)
if k is even

1 4 if k = ±1 0 if |k| > 1 is odd . ( 18
)
Their decay is slow because h(α) has discontinuous derivatives at ±π/2.

Bi-Lipschitz continuity

We mentioned that polynomial exponents amplify the variability of random variables around their mean. Indeed if (z, z

) ∈ C 2 then |z k - z k |/|z -z | may be arbitrarily large if k > 1.
On the contrary, a phase harmonic preserves the modulus which provides a bound on such amplification. It is proved in [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF] that it is Lipschitz continuous

∀(z, z ) ∈ C 2 , |[z] k -[z ] k | ≤ max(|k|, 1) |z -z |. (19)
The distance |z -z | is therefore amplified by at most |k|. This Lipschitz continuity is extended to phase windowed Fourier transforms, which are shown to be bi-Lipschitz. The Fourier transform preserves norms and distances. Calculating the norm of [START_REF] Simoncelli | The steerable pyramid: a flexible architecture for multi-scale derivative computation[END_REF] gives

H(z) = H(z) = h 2 |z| 2 , with h 2 = k | h(k)| 2 .
Following [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF], we also derive from [START_REF] Bortoli | Macrocanonical models for texture synthesis[END_REF] that a phase windowed Fourier transform is bi-Lipschitz over complex numbers (z, z

) ∈ C 2 A H |z -z | 2 ≤ H(z) -H(z) 2 ≤ B H |z -z | 2 , (20)
with

A H = 2 | h(1)| 2 and B H = | h(0)| 2 + k∈Z k 2 | h(k)| 2 .

Phase Harmonics of Fourier Coefficients

If L = F u computes the Fourier transform of X(u) then the resulting phase harmonic representation is

R(X) = H(F u X) = h(k) [ X(ω)] k k∈Z,ω∈R
.

It is indexed by v = (ω, k). We show that the covariance of H(F u X) is sparse but not necessarily diagonal. Section 2.2 explains that if X is stationary then X(ω) and X(ω

)
are not correlated if ω = ω but X(ω) k and X(ω ) k may be correlated if X is not Gaussian. Similarly we show that the covariance of [ X(ω)] k and [ X(ω )] k may become non-zero for ω = ω , which reveals non-Gaussian properties. Since H is bi-Lipschitz and the Fourier transform is unitary up to a factor d, it

results that R = HF u is bi-Lipschitz with Lipschitz constants A R = dA H and B R = dB H . Proposition 2.1 implies that dA H σ 2 (X) ≤ σ 2 ( H(X)) ≤ dB H σ 2 (X). ( 21
)
Sparse phase harmonic covariance Phase harmonic covariance coefficients are

K HFu (ω, k, ω , k ) = h(k) h(k ) * Cov([ X(ω)] k , [ X(ω )] k ). If h is compactly supported in [0, k max ], since X is of dimension d, the covariance K HF has (k max + 1) 2 d 2 coefficients.
The following proposition proves that coefficients of K HF are mostly zero and it is diagonal if X is Gaussian. Diagonal values are specified.

Theorem 2.2. If X is real stationary then Cov([ X(ω)] k , [ X(ω )] k ) = 0 if kω = k ω . (22)
Along the diagonal, if ω = 0 and k = 0 or if ω = 0 and k is odd then

Cov([ X(ω)] k , [ X(ω)] k ) = Cov( X(ω), X(ω)) . (23)
If ω = 0 and k is even then

Cov([ X(0)] k , [ X(0)] k ) = Cov(| X(0)|, | X(0)|) . ( 24
)
For all ω = 0 Cov(| X(ω)|, | X(ω)|) Cov( X(ω), X(ω)) = 1 - E(| X(ω)|) 2 E(| X(ω)| 2 ) . ( 25
)
If X is Gaussian then K HF is diagonal and E(| X(ω)|) 2 /E(| X(ω)| 2 ) = π/4 if ω = 0.
The proof is in Appendix C. All equalities or non-equalities on frequencies must be understood modulo 2π in dimension r. Property [START_REF] Georgii | Gibbs measures and phase transitions[END_REF] proves that K HFu is highly sparse. If X is not Gaussian then off-diagonal coefficients are typically not zero when kω = k ω . For k = k = 0, we get d 2 modulus covariances Cov(| X(ω)|, | X(ω )|) which are a priori non-zero for all (ω, ω ). It provides no information on the phase of X(ω) and X(ω ). Phase correlations are captured when k = 0. If ω = 0, it results from [START_REF] Georgii | Gibbs measures and phase transitions[END_REF] that Cov([

X(ω)] k , [ X(ω )] k ) = 0 only if ω is colinear with ω. For a grid Λ d of size d, there are O(d 1/r ) such frequencies ω . The total number of non-zero off-diagonal covariance coefficients is thus at most d 2 + O(d 1+1/r ).
For a non-Gaussian process, off-diagonal coefficients with kω = k ω are typically non-zero. This is illustrated with a random shift process

X(u) = x(u -S) where S is uniformly distributed in Λ d . If kω = k ω then one can verify that Cov( X(ω) k , X(ω ) k ) = [ x(ω)] k [ x(ω )] -k if ω = 0, ω = 0.
If x(ω) does not vanish then all these coefficients are non-zero.

Along the diagonal, ( 23) and [START_REF] Nocedal | Numerical Optimization[END_REF] prove that all coefficients are either equal to Cov( X(ω), X(ω)) or to Cov(| X(ω)|, | X(ω)|). Property (25) also proves that the ratio between these covariance values depend upon the ratio

E(| X(ω)|) 2 /E(| X(ω)| 2
). This ratio measures the sparsity of multiple realizations of the random variable X(ω). Indeed, expected values are estimated by normalized sums, so ratio of expected values can be interpreted as ratios of 1 over 2 norms, which is a measure or sparsity. For all Gaussian random variables, this ratio is equal to π/4. If the ratio is smaller than π/4 then X(ω) has a high probability to be relatively small and it has large amplitude outliers of low probability. The probability distribution of X(ω) is then highly peaked in zero with a long tail. The smaller this ratio the lower the entropy of this marginal probability distribution.

We saw in ( 8) that a maximum entropy model X introduces errors if its entropy is much larger than the entropy of X. The entropy H(p) is bounded by the sum of the entropy of each random variable X(ω), which is conditioned by Cov( X(ω), X(ω) and Cov(| X(ω)|, | X(ω)|). It gives an upper bound on the entropy H(p), which decreases if the sparsity of X(ω) increases. However, if X(u) has sharp localized transitions then its Fourier coefficients have typically a large amplitude across most frequencies and are not sparse. On the contrary, wavelet coefficients may be sparse. In this case, we shall capture this sparsity and thus compute maximum entropy models of lower entropy, by replacing the Fourier transform by a wavelet transform.

Wavelet Phase Harmonics

To model random vectors X whose realizations include singularities and sharp transitions, we define a phase harmonic representation R = H L where L is a wavelet transform as opposed to a Fourier transform. Indeed, wavelet transform provides sparse representations of such signals and Section 2.1 explains that imposing sparsity constraints typically improves the accuracy of maximum entropy models. We concentrate on two-dimensional image applications. The same approach applies to signals of any dimension, with an appropriate wavelet transform. Section 3.1 reviews the construction of complex steerable wavelet frames and the properties of their covariance matrices. Section 3.2 studies the covariance of wavelet phase harmonics.

Steerable Wavelet Frame Covariance

Complex steerable wavelet frames were introduced in [START_REF] Simoncelli | The steerable pyramid: a flexible architecture for multi-scale derivative computation[END_REF] and are further studied in [START_REF] Unser | A unifying parametric framework for 2d steerable wavelet transforms[END_REF], to easily compute wavelet coefficients of rotated images. For simplicity, we begin by introducing wavelets as a localized functions ψ(u) for u ∈ R 2 , with ψ(u) du = 0. Complex steerable wavelets have a Fourier ψ(ω) concentrated over one-half of the Fourier domain. We impose that ψ(-u) = ψ * (u) so that ψ(ω) is real. This Fourier transform is centered at a frequency ξ ∈ R 2 and is non negligible for ω ∈ R 2 such that |ω -ξ| ≤ C |ξ| for some C > 0. Figure 1 gives an example of such a wavelet, which is specified in Appendix D.

Let r be a rotation by an angle 2 π/Q. Multiscale steerable wavelets are derived from ψ with dilations by 2 j for j ∈ Z, and rotations over

Q angles θ = 2 π/Q for 0 ≤ < Q ψ λ (u) = 2 -j ψ(2 -j r -u) ⇒ ψ λ (ω) = 2 j ψ(2 j r ω) with λ = 2 -j r -ξ. ( 26
)
Since ψ(ω) is non negligible for |ω -ξ| ≤ C|ξ| it results that ψ λ (ω) is centered at λ and non negligible for |ω-λ| ≤ C|λ|. In space, ψ λ (u) is non-negligible for |u| ≤ C |λ| -1 . We limit the scale 2 j to a maximum 2 J . The lowest frequencies are captured by a wavelet centered at λ = 0. It is computed by dilating a function φ(u) such that φ(u) du = 1:

ψ 0 (u) = 2 -J φ(2 -J u) ⇒ ψ 0 (ω) = 2 J φ(2 J ω). (27)
A wavelet frame is constructed by translating each ψ λ for λ = 0 by u = 2 j-1 n and ψ 0 by u = 2 J-1 n for all n ∈ Z 2 . It introduces a factor 2 oversampling relatively to a wavelet orthonormal basis [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF], which creates some redundancy. The wavelet transform of x ∈ L 2 (R 2 ) is defined by

Wx = {x ψ λ (u)} (λ,u)∈Γ
where Γ is a frequency-space index set with (λ, u)

= (2 -j r -ξ, 2 j-1 n) for 1 ≤ j ≤ J, 0 ≤ < Q, n ∈ Z, or (λ, u) = (0, 2 J-1 n).
Under appropriate conditions on ψ, the wavelet family [START_REF] Unser | A unifying parametric framework for 2d steerable wavelet transforms[END_REF]. This means that there exists 0

{ψ λ (•-u)} (λ,u)∈Γ is a frame of L 2 (R 2 )
< A W ≤ B W such that for any x ∈ L 2 (R 2 ) A W x 2 ≤ Wx 2 ≤ B W x 2 (28)
with

Wx 2 = (λ,u)∈Γ |x ψ λ (u)| 2 .
The wavelet transform can be redefined over discrete images x of d pixels supported in a two-dimensional square grid Λ d , uniformly sampled at intervals 1 in [1, d 1/2 ] 2 . It requires to discretize and modify "boundary wavelets" whose supports intersect image boundaries. This can be done over steerable wavelets [START_REF] Simoncelli | The steerable pyramid: a flexible architecture for multi-scale derivative computation[END_REF][START_REF] Unser | A unifying parametric framework for 2d steerable wavelet transforms[END_REF], while preserving the frame constants A W and B W . The resulting wavelet

ψ λ (• -u) are supported in Λ d . They are still indexed by (λ, u) ∈ Γ. If λ = r -ξ = 0 for 1 ≤ j ≤ J, 0 ≤ < Q then u∈Λ d ψ λ (u) = 0. If λ = 0 then u∈Λ d ψ 0 (u) = 2 J . Since J ≤ (log 2 d)/2, there are at most Q(log 2 d)/2 + 1 different frequency channels λ. For λ = r -ξ = 0, ψ λ is translated by u = 2 j-1 n ∈ Λ d which yields 2 -j+1 d wavelet coefficients. The total number of wavelets coefficients is about 4Qd/3 if J = (log 2 d)/2.
The mother wavelet ψ is chosen in order to obtain a sparse wavelet representation of realizations of X, with few large amplitude wavelet coefficients. This sparsity highlights non-Gaussian properties. Figure 2 displays the modulus and phase of wavelet coefficients of the vorticity field of a turbulent flow. This flow is obtained by running the 2D Navier Stokes equation with periodic boundary conditions, initialized with a random Gaussian field [START_REF] Schneider | Coherent vortex extraction and simulation of 2d isotropic turbulence[END_REF]. After a fixed time, it defines a stationary but non-Gaussian random process. For each scale and orientation, large amplitude modulus coefficients are located at positions where the image has sharp transitions, and the phase depends upon the position of these sharp transitions.

Wavelet covariance A wavelet transform defines a linear representation R = W indexed by v = (λ, u). Similarly to Fourier coefficients, we show that wavelet coefficients have a covariance which nearly vanish at different frequencies.

Similarly to Fourier coefficients, wavelet coefficients have a zero mean at non-zero frequencies.

If λ = 0 then E(X ψ λ (u)) = 0 because u ψ λ (u) = 0. If λ = 0 then u ψ 0 (u) = 2 J so E(X ψ 0 (u)) = 2 J E(X(u)). The covariance at v = (λ, u) and v = (λ , u ) is K W (v, v ) = Cov(X ψ λ (u) , X ψ λ (u )) .
It depends on u -u because X is stationary. Let ψ λ (ω) be the discrete Fourier transform of ψ λ (u) defined in [START_REF] Lustig | Microcanonical monte carlo simulation of thermodynamic properties[END_REF]. Since wavelet coefficients are convolutions, covariance values can be rewritten from the power spectrum

K(ω) = 1 d Cov( X(ω), X(ω)) of X K W (v, v ) = ω∈Λ d K(ω) ψ λ (ω) ψ * λ (ω) e i(u-u ).ω . ( 29
) It results that K W (v, v ) = 0 if ψ λ (ω) ψ λ (ω) = 0 for all ω. Since ψ λ (ω) is non-negligible only if |ω -λ| ≤ C|λ|, the covariance K W (v, v ) is non-negligible for λ = λ only if |λ -λ | |λ| + |λ | ≤ C . ( 30
)
It shows that similarly to Fourier coefficients, wavelet covariances are negligible across frequencies which are sufficiently far apart.

Maximum entropy wavelet graph model A maximum entropy model conditioned by wavelet covariances is Gaussian because the wavelet transform is linear. Figure 3(a) gives a realization of a stationary turbulent flow X. A low-dimensional maximum entropy model is defined on a graph of covariance coefficients K W (v, v ) for v = (λ , u ) in a small neighborhood of v = (λ, u). Since wavelet coefficients are nearly decorrelated across frequencies, at each scale 2 j , the neighborhood of v = (λ, u) is 14 defined as the set of v = (u , λ ) such that λ = λ and |u -u | ≤ 2 j-1 ∆ for a fixed ∆. Covariances are thus specified over a spatial range proportional to the scale. This is a foveal neighborhood which is sufficient to approximate the covariances of large classes of random processes such as fractional Brownian motions [START_REF] Andreux | Foveal Autoregressive Neural Time-Series Modeling[END_REF]. Since (u, u ) = 2 j-1 (n, n ), the neighborhoods of all v have the same size, which is smaller than (2∆ + 1) 2 . Since X is stationary, its probability distribution is invariant to the group G of translations. The number of covariance coefficients that must be estimated is equal to the number |E G | of edges in the graph modulo translations. It is equal to the number of wavelet frequencies λ multiplied by the size of each neighborhood, and thus bounded by (JQ + 1)(2∆ + 1) 2 = O(log 2 d). This model size is much smaller than the image size d when d is large. 2) is explained in Appendix A. To measure the accuracy of this model we compare the power spectrum of X and the Gaussian model X. Both processes are isotropic so Figure 3(c) gives the radial log power spectrum of X and X as a function of |ω|. These power spectrum are nearly the same, which means that the wavelet covariance graph gives an accurate estimation of the second order moments of X from a single realization.

The realization of the Gaussian model in Figure 3(b) has a geometry which is very different from the turbulence flow X in Figure 3(a). It shows that X is highly non-Gaussian, which also appears in Figure 2. The wavelet coefficients of a stationary process X are not correlated at different scales and angles, which would imply that they are independent if X was Gaussian. On the contrary, Figure 2 shows that the the modulus and phases of wavelet coefficients of X are strongly dependent across scales and angles. High amplitude modulus coefficients are located in the same spatial neighborhoods because they are produced by the same sharp transitions of the flow. Next section explains how to capture this dependence with phase harmonics.

Wavelet Phase Harmonic Covariance

A wavelet phase harmonic representation Φ = H L is computed with a linear wavelet transform L = W. As in the Fourier case, the covariance of Φ(X) = H(WX) has nonzero covariance coefficients across different frequency bands. We study the properties of the resulting covariance matrix, and the role of sparsity.

Wavelet phase harmonics To specify the dependence across frequencies, we apply a phase harmonic operator to wavelet coefficients:

H(Wx) = h(k) [x ψ λ (u)] k
K HW (v, v ) = h(k) h(k ) * Cov([X ψ λ (u)] k , [X ψ λ (u )] k ).
Since X is stationary, it only depends on u -u . Wavelet harmonic covariances only need to be calculated for k ≥ 0 and k ≥ 0. Indeed, since X is real and ψ(-u) = ψ * (u) one can verify that K HW (v, v ) does not change its value if (λ, k) becomes (-λ, -k) or if (k, k ) becomes (-k, -k ). Such wavelet harmonic covariances have first been computed by Portilla and Simoncelli [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] to characterize the statistics of image textures. Their representation correspond to (k, k ) equal to (0, 0), (1, 1), [START_REF] Jaynes | Information theory and statistical mechanics[END_REF][START_REF] Cover | Elements of Information Theory[END_REF], which amounts to choosing ĥ(k) = 1 [0,2] (k).

Since H and W are bi-Lipschitz the operator HW is also bi-Lipschitz, with lower and upper bounds A H A W and B H B W . Proposition 2.1 implies that

A H A W σ 2 (X) ≤ σ 2 ( H(WX)) ≤ B H B W σ 2 (X), (31)
which controls the variance of wavelet harmonic coefficients.

Rectified neural network coefficients Ustyuzhaninov et. al. in [START_REF] Ustyuzhaninov | What does it take to generate natural textures[END_REF] have shown that one can get good texture synthesis from the covariance of a one-layer convolutional neural network, computed with a rectifier. In the following we show that these statistics are equivalent to phase harmonic covariances, computed with a rectifier phase window h(α). Section 2.3 proves that H = F α H, where H computes a phase windowing of wavelet coefficients

H(Wx) = |x ψ λ (u)| h(ϕ(x ψ λ (u)) + α)} (λ,u)∈Γ,α∈[0,2π] .
The covariance of H(WX) and H(WX) thus satisfy K HW = F α K HW F -1 α . The following proposition proves that K HW gives the covariance of rectified wavelet coefficients if h is a rectifier phase window. Proposition 3.1. Let v = (λ, α, u) and v = (λ , α , u ). For a rectifier phase window h(α) = ρ(cos α)

K HW (v, v ) = Cov ρ(X ψ λ,α (u)) , ρ(X ψ λ ,α (u )) ( 32
)
with ψ λ,α (u) = Real(e -iα ψ λ (u)).

Proof: We proved in ( 17) that if h(α) = ρ(cos α) then

H(z) = {ρ(Real(e iα z))} α∈[0,2π] .
It results that

K HW (v, v ) = Cov(ρ(Real(e iα X ψ λ (u))), ρ(Real(e iα X ψ λ (u ))) ,
which proves (32). Rectified wavelet coefficients ρ(X ψ λ,α (u)) can be interpreted as one-layer convolutional network coefficients, computed with wavelet filters ψ λ,α of different frequencies λ and different phases α. This result applies to any linear operator L, not necessarily a wavelet transform. The statistics used by Ustyuzhaninov et. al. in [START_REF] Ustyuzhaninov | What does it take to generate natural textures[END_REF] use local cosine or random filters in their network as opposed to steerable wavelets. It corresponds to convolutional operators L that are not wavelet transforms.

Sparse harmonic covariance Proposition 2.2 specifies the properties of a Fourier phase harmonic covariance. We qualitatively explain, without proof, why Cov([X ψ λ (u)] k , [X ψ λ (u )] k ) has similar sparsity properties.

For a fixed λ and k, [X ψ λ (u)] k is a stationary random vector in u. The covariance of [X ψ λ (u)] k and [X ψ λ (u )] k is non-zero only if their power-spectrum have a support which overlap. We give a necessary condition by approximating these spectrum supports.

For k = 1, the spectrum of X ψ λ (u) has an energy concentrated at frequencies where the Fourier transform of ψ λ is concentrated, which is included in a ball centered at λ of radius C|λ|. For |k| > 1, [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF] explains that the Fourier transform of [X ψ λ (u)] k and hence its power spectrum is concentrated in a ball centered at kλ of radius |k|C|λ|. If k = 0 then the spectrum of |X ψ λ (u)| is concentrated in a ball centered at the 0 frequency, of radius C|λ|. This is illustrated numerically in Figure 4 which displays the power spectrum of [X ψ λ (u)] k for k = 0, 1, 2. In this case, X is a stationary vorticity field shown at the top of Figure 3(a). These power spectrum are estimated for a fixed λ from 100 independent realizations of X. For k = 1, the power spectrum is supported over the Fourier support of ψ λ . For k = 2 its support is approximately dilated by 2, whereas for k = 0 it is centered at the zero frequency. ) occur for kλ ≈ k λ . This is similar to the vanishing property [START_REF] Georgii | Gibbs measures and phase transitions[END_REF] of Fourier harmonic coefficients, at frequencies which are colinear. Since λ = 2 -j r -ξ and λ = 2 -j r -ξ, it requires that 2 j-j ≈ |k |/|k| and that ≈ .

Diagonal covariance coefficients have similar properties as Fourier coefficients. They are specified by Cov(|X ψ λ (u)|, |X ψ λ (u)|) and Cov(X ψ λ (u), X ψ λ (u)), which do not depend upon u. Moreover when λ = 0,

Cov(|X ψ λ (u)|, |X ψ λ (u)|) Cov(X ψ λ (u), X ψ λ (u)) = 1 - E(|X ψ λ (u)|) 2 E(|X ψ λ (u)| 2 )
.

The ratio

E(|X ψ λ (u)|) 2 /E(|X ψ λ (u)| 2
) measures the sparsity of wavelet coefficients X ψ λ (u). Indeed, the sparsity of a time series is usually quantified by its 1 norm for a fixed 2 norm [START_REF] Mallat | A Wavelet Tour of Signal Processing: The Sparse Way[END_REF], which amounts to compute their ratios. Since X ψ λ (u) is stationary, spatial sums over u converge to expected values. It results that a ratio of an 1 over an

2 norm ( u |X ψ λ (u)|) 2 /( u |X ψ λ (u)| 2
) is an estimator of the ratio of expected values. If X is Gaussian then X ψ λ (u) is a complex Gaussian random variable and one can verify that this ratio is π 4 for all λ. A ratio smaller than π 4 implies that wavelet coefficients X ψ λ (u) are sparse and hence that X is non-Gaussian.

To increase the accuracy of a maximum entropy model X, the Kullback-Leibler divergence [START_REF] Ustyuzhaninov | What does it take to generate natural textures[END_REF] shows that we must minimize its entropy. As in the Fourier case, an upper bound of the entropy is obtained as a sum of the entropy of the marginals of wavelet coefficients. The marginal entropies get smaller by reducing the sparsity ratio

E(|X ψ λ (u)|) 2 /E(|X ψ λ (u)| 2 ).
To minimize this upper bound of the model entropy, this suggests by finding a mother wavelet ψ which yields sparse coefficients.

Gaussianity test Even though ψ λ and ψ λ may have disjoint Fourier supports, the previous analysis showed that one can find k and k such that the spectrum of [X ψ λ (u)] k and [X ψ λ (u )] k overlap, for example with k = k = 0. A priori, the covariance K HW (v, v ) is then non-zero, unless X is Gaussian in which case these coefficients vanish, as proved by the following theorem. Proposition 3.2. Let (λ, λ ) be such that ψ λ ψ λ = 0. If X is Gaussian and stationary then

K HW (v, v ) = 0 if v = (λ, k, u) and v = (λ , k , u ), for any (k, u, k , u ).
Proof: If X is Gaussian then X ψ λ (u) and X ψ λ (u ) are jointly Gaussian random variables. Equation [START_REF] Leonarduzzi | Maximum-entropy Scattering Models for Financial Time Series[END_REF] proves that if ψ λ ψ λ = 0 then Cov(X ψ λ (u), X ψ λ (u )) = 0. These Gaussian random variables are uncorrelated and therefore independent. It results that [X ψ λ (u)] k and [X ψ λ (u )] k are also independent and thus have a covariance which is zero.

This proposition gives a test of Gaussianity by considering two frequencies λ and λ where |λ -λ | is sufficiently large so that the support of ψ λ and ψ λ do not overlap. If kλ ≈ k λ then wavelet harmonic covariances are zero if X is Gaussian but it is typically non-zero if X is not Gaussian.

Microcanonical Models

Section 2.1 introduces macrocanonical maximum entropy models conditioned by covariance coefficients of a representation R(X). The resulting maximum entropy distribution p depend upon Lagrange coefficients which are computationally very expensive to calculate, despite the development of efficient algorithms [START_REF] Bortoli | Macrocanonical models for texture synthesis[END_REF]. Section 4.1 reviews maximum entropy microcanonical models which avoid computing these Lagrange multipliers. Section 4.2 gives an alternative microcanonical model which is calculated with a faster algorithm, which transports a maximum entropy measure by gradient descent.

Maximum Entropy Microcanonical Models

Microcanonical models avoid the calculation of Lagrange multipliers and are guaranteed to exist as opposed to macrocanonical models. We first specify covariance estimators and then briefly review the properties of these microcanonical models.

Microcanonical models rely on an ergodicity property which insures that covariance estimations concentrate near the true covariance K R when d is sufficiently large. Let G be a known group of linear unitary symmetries of the density p of X. It includes translations because X is stationary. An estimation of the covariance K R over the edge set E G is computed in (4) from a single realization x of X, by transforming x with all g ∈ G. To differentiate realizations of X from other x ∈ R d , we associate a similar covariance estimation to any x ∈ R d . This covariance is centered on the empirical mean M R computed in (3) and defined by

K Rx (v, v ) = 1 |G| g∈G R v (g.x) -M R (v) R v (g.x) -M R (v ) * . ( 34
)
It is invariant to the action of any g ∈ G on x.

We denote K R 2

E G = (v,v )∈E G |K R (v, v )| 2 .
We shall suppose that X satisfies the following covariance ergodicity property over E G :

∀ > 0, lim d→∞ Prob( K RX -K R E G ≤ ) = 1. ( 35
)
To satisfy this property the number |E G | of covariance moments must be small compared to the dimension d of X. The ergodicity hypothesis implies that when d is sufficiently large, with high probability, the empirical covariance K Rx of a realization x of X is close to the true covariance K R over E G . Given a realization x of X, a microcanonical set of width is the set of all x ∈ R d which have nearly the same covariance estimations as x:

Ω = {x ∈ R d : K Rx -K Rx E G ≤ }. (36)
A maximum entropy microcanonical model has a probability distribution of maximum entropy supported in Ω . The set Ω is bounded. Indeed, if x ∈ Ω then

K Rx E G ≤ K Rx E G + . Since E G includes all diagonal coefficients (v, v) for v ∈ V ,
one can derive an upper bound for Rx . Since R is bi-Lipschitz, we also obtain an upper bound for x . Since > 0, the Lebesgue measure of Ω is non-zero, therefore the maximum entropy distribution is uniform in Ω relatively to the Lebesgue measure. A major issue in statistical physics is to find sufficient conditions to prove the Boltzmann equivalence principle which guarantees the convergence of microcanonical and macrocanonical models towards the same Gibbs measures when d goes to ∞ [START_REF] Dembo | Large Deviations Techniques and Applications[END_REF]. This involves the proof of a large deviation principle which expresses concentration properties of the covariance of R(X). If R is continuous and bounded so that the interaction potential (R -M R )(R -M R ) * is also continuous and bounded, and if there is no phase transition, which means that the limit is a unique Gibbs measure, then one can prove that microcanonical and macrocanonical measures converge to the same limit for an appropriate topology [START_REF] Deuschel | Microcanonical distributions for lattice gases[END_REF][START_REF] Georgii | Gibbs measures and phase transitions[END_REF]. The bounded hypothesis is not necessary and may not be satisfied. The ergodicity property (35) guarantees that X concentrates in Ω with a high probability when d is sufficiently large. The microcanonical set Ω may however be much larger than the set where X concentrates, which means that the maximum entropy microcanonical distribution may have a much larger entropy than the entropy of X. As in the macrocanonical case, the representation R must be optimized in order to reduce the maximum entropy, which motives the use of sparse representations.

Gradient-Descent Microcanonical Models

Sampling a maximum entropy microcanonical set requires to use Monte Carlo algorithms. They are computationally very expensive when the number of moments |E G | and the dimension d are large, because their mixing time become prohibitive [START_REF] Lustig | Microcanonical monte carlo simulation of thermodynamic properties[END_REF]. Following the approach in [START_REF] Bruna | Multiscale sparse microcanonical models[END_REF], we approximate these microcanonical models with a gradient descent algorithm.

Gradient-descent microcanonical models are computed by transporting an initial Gaussian white noise measure into the microcanonical set Ω . This transport is calculated with a gradient descent which progressively minimizes

f (x) = K Rx -K Rx 2 E G , (37)
with a sufficiently large number of iterations so that f (x) < and hence x ∈ Ω . The initial Gaussian white noise is a maximum entropy distribution conditioned by a variance σ 2 . This variance is chosen to be an upper bound of the empirical variance of x(u) along u, calculated from diagonal coefficients of K Rx . It guarantees that the resulting white noise has an entropy larger than the microcanonical model. The gradient descent progressively reduces this entropy while transporting the measure towards Ω . Entropy reduction bounds are computed in [START_REF] Bruna | Multiscale sparse microcanonical models[END_REF].

The initial x 0 is Gaussian white noise. For all t ≥ 0, the gradient descent iteratively computes

x t+1 = x t -η ∇f (x t ) .
This operation transports the probability measure µ t of x t into a measure µ t+1 of x t+1 . We stop the algorithm at a time t = T which is large enough so that f (x T ) < with a high probability. Since f (x) is not convex, there is no guarantee that f (x T ) < even for large T . In numerical calculations in Section 5, we may use several initializations, typically 10, to evaluate the model. The differentiability of f relies on the smoothness of the function z → [z] k , which is differentiable at z = 0. The implementation uses an analytic formula of its derivative at z = 0 and applies the same formula at z = 0. The empirical covariance K Rx is computed in (34) with an average over all symmetries of a known group G of linear unitary operators in R d . The following theorem proves that G is also a group of symmetries of the probability measures obtained by gradient descent. Theorem 4.1. For any t ≥ 0, the probability measure µ t of x t is invariant to the action of G.

The proof is in Appendix E. It is a minor modification of the proof in [START_REF] Bruna | Multiscale sparse microcanonical models[END_REF]. We replace the gradient descent method by L-BFGS algorithm with line search [START_REF] Nocedal | Numerical Optimization[END_REF]. It has a faster and more accurate numerical convergence. Since G includes translations, this theorem implies that each x t is stationary. It is shown in [START_REF] Bruna | Multiscale sparse microcanonical models[END_REF] that the transported density µ T supported in Ω may be different from the maximum entropy density in Ω because the gradient descent reduces too much the entropy. In general, the precision of these microcanonical gradient descent models are not well understood. However, this theorem proves that the gradient descent preserves all known symmetries of the distribution of X, which is also true for a maximum entropy measure. The sampling algorithm in [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] is different since it relies on iterative gradient projections for x to satisfy constraints of various nature, which is not guaranteed to preserve symmetries.

The gradient descent may converge faster by preconditioning f (x). This is done by replacing estimated covariances K Rx (v, v ) by normalized correlation coefficients

K Rx (v, v ) K Rx (v, v) 1/2 K Rx (v , v ) 1/2 . ( 38
)
Any g ∈ G is a symmetry of K Rx and is therefore a symmetry of normalized correlation coefficients, so Theorem 4.1 remains valid with this preconditioning.

Foveal Wavelet Harmonic Covariance Models

We study microcanonical models conditioned by wavelet harmonic covariances, computed with different symmetry groups G, and different sufficient statistics set E G . Section 5.2 introduces foveal models which limit the multiscale spatial range of coefficients in E G . Section 5.3 gives a methodology to evaluate the precision of different foveal models, with numerical results.

Rotation and Reflection Symmetries

The covariance K R is estimated from a single realization of X, by taking advantage of a known group of symmetries G of the probability distribution of X. For a wavelet phase harmonic representation R = HW, we specify the properties of K R when G includes sign changes, reflections and rotations.

The following proposition proves that some covariance coefficients vanish when symmetries include a sign change or a central reflection. These covariances thus do not need to be included in the sufficient statistics set E G . The proof is in Appendix F. u) is a symmetry of the probability distribution of X and h is real then K HW (v, v ), and K HWx (v, v ) are real.

Proposition 5.1. Let v = (λ, k, u) and v = (λ , k , u ). (i) If the sign change g.x = -x is a symmetry of the probability distribution of X then K HW (v, v ) = K HWx (v, v ) = 0 if k + k is odd. (ii) If the central reflection g.x(u) = x(-
Isotropic models An isotropic random process X has a probability distribution which is invariant by rotations. The group G then includes all translations and rotations. We show that K HW becomes sparse after applying a Fourier transform on rotation angles.

If g = r η is a rotation by 2πη/Q for 0 ≤ η < Q then (r η .x) ψ λ (u) = x ψ rηλ (r -η u) where r η λ = 2 -j r η-ξ is the rotation of λ = 2 -j r -ξ. To eliminate the effect of the change of position r -η u on covariance coefficients, we only keep wavelet harmonic covariances at a same spatial position. This means that v = (λ , k , u ) is a neighbor of v = (λ, k, u) in the covariance graph model only if u = u . It implies that (v, v ) ∈ E G only if v and v have the same spatial position.

Isotropy is a form of stationarity along rotations angles. To diagonalize angular covariance matrices, we use a discrete Fourier transform along rotations written F . The discrete Fourier transform of y( ) for 0

≤ < Q at a frequency 0 ≤ m < Q is (F y)(m) = Q-1 =0 y( ) e -i2mπ /Q . The representation R(X) = F H(WX) computes F ([X ψ λ (u)] k ) for λ = 2 -j r -ξ with (j, u, k) fixed and varying. It is indexed by v = (j, m, k, u). The covariance matrix of R(X) is K R = F K HW F -1 .
We consider the restriction of K R to E G . The next theorem proves that if X is isotropic then K R has diagonal angular Fourier matrices. Isotropic processes X may also have a probability distribution which are invariant to line reflections. A line reflection of orientation η computes g.x(u) = x(u η ), where u η is symmetric to u relatively to a line going through the origin in R 2 , with an orientation η ∈ [0, 2π]. If X is isotropic and invariant to a line reflection for an angle η then it is invariant to line reflections for any η ∈ [0, 2π]. The following theorem applies to the bump steerable wavelets used in numerical experiments. Theorem 5.1. Let v = (j, m, k, u), v = (j , m , k , u) and u = (u 1 , u 2 ). If the probability distribution of X stationary and isotropic and R = F HW then

K R (v, v ) = K Rx (v, v ) = 0 if m = m . (39)
Furthermore, if the distribution of X is invariant to line reflections and if the wavelet satisfies

ψ(u 1 , -u 2 ) = ψ(u 1 , u 2 ) and φ(u 1 , -u 2 ) = φ(u 1 , u 2 ) then K R (v, v ) and K Rx (v, v ) are real if m = m .
The proof is in Appendix G. This theorem proves that the sufficient statistics set can be reduced to diagonal angular coefficients (v, v ) ∈ E G with m = m . It reduces its size by a factor Q. Invariance to line reflections implies that it is sufficient to keep the real part of these diagonal values. The rotation invariance strictly applies to the process defined on a continuous variable u ∈ R 2 . On discrete images it is not valid at the finest scale because of the square sampling grid, and it is not valid at the largest scale because of their square support. For preconditioning, the covariance of the isotropic model is also normalized with (38) and the Fourier transform F is then applied on the normalized coefficients.

Foveal Wavelet Harmonic Covariance

A wavelet harmonic covariance model is defined by the choice of the harmonic weights ĥ, by the symmetry group G and neighborhood relations which defines the edge set E. In the following we shall impose that ĥ(k) = 1 [k min ,kmax] (k), which limits harmonic exponents in the range [k min , k max ]. We define several foveal models which capture different properties.

Foveal models Wavelet harmonic coefficients are indexed by v = (λ, k, u), with λ = 2 -j r -ξ and u = 2 j-1 n. The covariance graph model is specified by the neighborhoods N v . A foveal model defines neighborhoods whose size doe not depend upon v. The range of spatial, scale and angular parameters is limited by three parameters ∆ n , ∆ j and ∆ .

A vertex v = (λ , k , u ) with λ = 2 -j r -ξ and u = 2 j -1 n is a neighbor of v = (λ, k, u) only if |n -n | ≤ ∆ n , |j -j | ≤ ∆ j , | -| ≤ ∆ , (k, k ) ∈ [k min , k max ] 2 .
A foveal model has a spatial range proportional to the scale. Long range spatial correlations are partly captured because |u -u | = |n2 j-1 -n 2 j -1 | become large at large scales. It provides high frequency correlations between close points and lowfrequency correlations between far away points. It is similar to a visual fovea [START_REF] Mallat | Foveal detection and approximation for singularities[END_REF]. Such foveal models have been used by Portilla and Simoncelli [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] to synthesize image textures.

Because of translation invariance, the sufficient statistics E G can be defined by setting n = 0. Since there are Q angles and at most (log 2 d)/2 scales j, the size of E G is at most

|E G | = O(Q∆ (k max -k min + 1) 2 (2∆ n + 1) 2 ∆ j log 2 d).
Increasing the values of k max , ∆ j , ∆ , ∆ n decreases the model bias but it also increases the size |E G | and hence the variance of the estimation. To ensure the bi-Lipschitz continuity of R = HW, we impose that k min ≤ 1 ≤ k max . In the following we describe several models of different sizes, which capture different properties of X. Each realization x is an image of d = 256 2 pixels. We use the bump steerable wavelets of Appendix D, computed on J = 5 scales and Q = 16 angles. The maximum scale 2 J depends on the integral scale of X, which is the distance beyond which all coefficients are nearly independent. A large number of angles Q gives a finer angular resolution. We specify 4 models corresponding to different choices of neighborhood parameters. For the first three models, the symmetry group G is reduced to translations where as the last model is also invariant to rotations.

• Model A with k min = k max = 1. It corresponds to a Gaussian maximum entropy wavelet model of Section 3.1. The covariance of wavelet coefficients is neglected across scales and angles by setting ∆ j = 0, ∆ = 0 and ∆ n = 2. The relative dimension of this model is G |/d = 1.2 10 -2 . Theorem 4.1 proves that the resulting gradient descent microcanonical distribution is invariant to these Q rotations. The Python implementation uses GPU (graphics processing unit) to synthesize non-Gaussian samples. Each synthesis sample (at d = 256 2 ) takes 30min for Model B, 50min for Model C, and 8min for Model D, by running a maximal 5000 number of L-BFGS iterations. The GPU memory cost is between 3GB to 6GB.

|E A G |/d = 3.

Visual evaluation

The quality of different microcanonical models is first evaluated visually. Among 10 synthesis of each model, we retain the one which yields the smallest loss in (37). The top row of Figure 5 shows a realization x of different stationary processes X. The first and second columns display isotropic and non-istotropic vorticity fields of two-dimensional turbulent flows, with a factor 2 zoom on the central part of each image. The third and fourth columns show an image of bubbles and a realization of Multifractal Random Walk [START_REF] Bacry | Multifractal random walk[END_REF]. An MRW is a self-similar random process with long range dependencies. We model the increments of MRW, which is a stationary process and limit its maximum correlation scale to 2 5 . The next rows give realizations of the microcanonical models A, B, C and D computed from the same realization x shown at the top.

As in Figure 3, the foveal Gaussian model A looses most geometric structures. On the contrary the models B, C, D recover most of this geometry. The model B captures the correlations of modulus coefficients across angles whereas the model C also impose covariance conditions on phases. There are little visual differences on turbulent flows but it is more apparent on bubble images. Model B does not reproduce closed bubbles whereas model C recovers bubbles having a better geometry. The model D is an isotropic version of model C. Its realizations are therefore isotropic. When X is isotropic then it is visually as precise as C but the variance reduction is not visible. If X is not isotropic, as in the turbulence of the second column, then model D does not reproduce the angular anisotropy.

Evaluations of Foveal Covariance Models

A microcanonical model X is conditioned by values of the covariance of R v (X) on a graph whose edges E relate neighbor vertices (v, v ). We differentiate two types of model errors. Type I errors result from the use of graph neighborhoods which are too small. In this case, the covariance of R(X) and R( X) may be different for (v, v ) ∈ V 2 -E. Type II errors are due to the choice of the representation R. Such errors are evaluated by comparing high order moments of X and X, which are not calculated by R. These comparisons are performed for R = HW.

Wavelet harmonic covariance error Comparisons are performed over correlation matrices, which normalizes covariance values. Type I errors are calculated by comparing the correlation values of X and X over a neighborhood which is much larger than the ones used by the different models. This paragraph considers the case where this neighborhood includes all scales and angles, but over a limited spatial range.

We compute correlation coefficients (38) by normalizing K R with its diagonal:

C R = D -1/2 R K R D -1/2 R where D R = diag(K R ). ( 40
)
We estimate K R with an empirical average over 100 realizations of X. This gives an accurate estimation of both D R and C R . An empirical estimator C Rx of C R is computed from a single realization x of X, with the same normalization

C Rx = D -1/2 R K Rx D -1/2 R
.

It is calculated on V 2 0 where V 0 ⊂ V is a foveal subset of all vertices. It incorporates correlations between all scales and angles, across a limited spatial range.

Let C R op,V 2 0 be the operator norm and hence the largest eigenvalue of the symmetric matrix C R restricted to V 2 0 . The following empirical error measures the estimator error of C R from one realization x

emp = C R -C Rx op,V 0 C R op,V 0 . (41)
It is a variance term due to variabilities of realizations x of X. This empirical error is compared to the estimation error of C R computed from the different microcanonical models. Let K model Rx be the covariance matrix of wavelet harmonic coefficients of a microcanonical model X instead of X. This matrix is estimated from 10 realizations of X. We compare K model

= C R -C model Rx op,V 0 C R op,V 0 . ( 42
)
This error has a variance term because the microcanonical model depends upon a particular realization x of X. It has also a bias term because the model is calculated

E G ⊂ V 2 0 is too small and is unable to reproduce the correlations C R in V 2 0 -E G then model is larger than emp . If E G = V 2
0 then the variance term dominates and model = emp . We define V 0 with k min = 0, k max = 4, a maximum range of scales ∆ j = J and a maximum range of angles ∆ = Q/2, but a small translation range ∆ n = 2. Its size is |V 0 | = 8025, so errors are evaluated over 8025 × 8025 correlation values. Table 1 compares the empirical error emp in (41) with the model error model in (42), for models A, B, C, D. We report the mean of these estimated errors by averaging the operator norms over 10 realizations x of X. The standard deviation is given in brackets. These error values are consistent with visual evaluations.

Table 1 show that the Gaussian model A gives a larger error than emp , because it captures no dependence across scales and angles. It introduces a large model bias. For models B and C, model < emp . The increase of E G dramatically reduces the model bias, which clearly appears visually. These models are only conditioned by wavelet harmonic covariances over neighbor scales, but the microcanonical model propagates these constraints across all scales. It provides a good approximation of covariance values between far away scales, evaluated in V 2 0 . The error of model C is smaller than model B, suggesting that scale correlations through phase also play an important role. This also appears in the visual quality of synthesized bubble images. The model D has a similar error as model C for the isotropic turbulence although it has a lower variance because estimators are averaged across Q angular directions. This indicates that the error is dominated by the model bias term and hence that a better model would be obtained by further increasing E G . This is verified on model D by increasing ∆ j = J and k max = 4, which yields a smaller covariance error 0.20 (0.02) on the isotropic turbulence. On the contrary, the model D has an error which is much larger than model C on the non-isotropic turbulence, which also appears visually.

Long-range spatial correlations Our microcanonical models are computed with foveal sufficient statistics sets, which are conditioned over relatively small spatial neighborhood at each scale 2 j . Such models can therefore introduce a large error if there exists long range spatial correlations between wavelet harmonic coefficients. This para- that in this case the covariance terms across scales are sufficient to capture high order moments. On the contrary, the model A has significant errors for q = 2 because it includes no covariance terms across scales and angles. Table 3 reports the value of these errors for the MRW. The errors of model D remain below the errors of model A but are significant in this case. This is probably due to the fact that the model D does not capture well the long range spatial correlations of the MRW, as shown in Figure 6(b). It introduces a large model error which appears in higher order moments. Table 2: The structure function error st (j, q) of model A,D for the isotropic turbulent vorticity field as a function of q, for j = 1, 2.

This analysis shows that microcanonical models computed from phase harmonic covariances in foveal neighborhoods capture important non-Gaussian properties. It gives accurate models of large classes of processes such as the turbulence examples. The eval-q 1 2 3 4 5

A st (1, q) 0.06 (0.02) 0.04 (0.03) 0.16 (0.06) 0.38 (0.06) 0.60 (0.04) D st (1, q) 0.02 (0.01) 0.04 (0.03) 0.05 (0.05) 0.13 (0.08) 0.28 (0.09) Table 3: The structure function error st (j, q) of model A,D for the increment of MRW process as a function of q, for j = 1.

uation methodology with moments is also able to detect model errors. It shows that long range correlations of wavelet harmonic coefficients are not well captured by foveal models. To circumvent this issue without increasing too much the model size, one can capture these long range correlations with another wavelet transform at each scale, which defines a scattering transform [START_REF] Leonarduzzi | Maximum-entropy Scattering Models for Financial Time Series[END_REF]. It amounts to cascade several wavelet transforms and non-linearities to define R as in a multilayer convolutional neural network. Very good visual texture syntheses have been obtained with covariances of multi-layer convolutional neural network coefficients [START_REF] Gatys | Texture synthesis using convolutional neural networks[END_REF], but their mathematical properties are not understood.

A Maximum Entropy Wavelet Covariance model

We explain how to compute the Lagrangian multiples of a maximum entropy stationary Gaussian random vector specified by wavelet covariance values, which is used for the model A.

To solve the maximum entropy problem, we follow [START_REF] Cover | Elements of Information Theory[END_REF] and derive a dual problem to minimize the entropy with respect to the Lagrangian multiples. To simplify the deviation, we assume that X is a zero-mean stationary process, and we consider only the translation group G. The wavelet transform is indexed by v = (λ, u). The maximumentropy wavelet covariance model can be written p

(x) = 1 Z exp - 1 2 (v,v )∈E β v,v g∈G x ψ λ (u -g) x ψ * λ (u -g) .
The sum over E is real-valued because β v,v = β * v ,v . With the Parseval formula, the sum over translations in G can be written in the Fourier domain as a sum over frequencies p

(x) = 1 Z exp - 1 2d (v,v )∈E β v,v ω | x(ω)| 2 ψ λ (ω) ψ * λ (ω) e -iω.(u-u ) .
The power spectrum P (ω) of p is

P (ω) = 1 d |x(ω)| 2 p(x)dx = (v,v )∈E β v,v ψ λ (ω) ψ * λ (ω)e -iω.(u-u ) -1 .
To find the optimal Lagrange multipliers, we minimize the entropy of p with respect to all the dual variables β v,v . The entropy has a closed form,

H(p) = 1 2 (v,v )∈E β v,v ω P (ω) ψ λ (ω) ψ λ (ω) * e -iω.(u-u ) + log(Z). (43)
The partition function Z also has a closed form,

log(Z) = 1 2 ω log P (ω) + d 2 log(2π).
Note that the β v,v should be constrained so that P (ω) > 0, ∀ω. This constrained optimization problem is addressed by setting the entropy loss to infinity if any condition is violated. The derivative of H(p) with respect to β v,v is calculated from (43). We can thus use an unconstrained optimization solver L-BFGS to solve this problem from the real and imaginary parts of β v,v and β v ,v . The optimal solution gives the power spectrum P (ω) of the Gaussian model A. Samples of model A are obtained by sampling a stationary Gaussian random vector whose power spectrum is P (ω).

B Proof of Theorem 2.1

We show that for the rectifier ρ(a) = max(a, 0),

H(z) -H(z ) 2 ≤ 1 4 |z -z | 2 , ∀(z, z ) ∈ C 2 . ( 44
and for any a ∈ R

1 2 |a -a | 2 ≤ |ρ(a) -ρ(a )| 2 + |ρ(-a) -ρ(-a )| 2 ≤ |a -a | 2 . ( 46
)
For any z ∈ C, Real(e i(α+π) z) = Real(-e iα z) and Real(e i(α+π/2) z) = Imag(e iα z), so (45) and (46) imply that C Proof of Theorem 2.2

Property ( 22) is proved similarly to (13) by translating X by τ ∈ Λ d and by using its stationarity to verify that Cov([ X(ω)] k , [ X(ω )] k ) = e i(kω-k ω ).τ Cov([ X(ω)] k , [ X(ω )] k ) so both terms vanish if kω = k ω . Property ( 25) is verified for ω = 0 by decomposing nominator and denominator with Cov(A, B) = E(AB * ) -E(A) E(B * ). Property [START_REF] Lustig | Microcanonical monte carlo simulation of thermodynamic properties[END_REF] uses the same decomposition and the fact that E( X(ω)) = 0 and E([ X(ω)] k ) = 0 if kω = 0 . This is proved again by translating X by τ ∈ Λ d , which transforms X(ω) into e -iτ.ω X(ω). Since X is stationary, it does not modify this expected value and hence E([ X(ω)] k ) = e -ikτ.ω E([ X(ω)] k ). Since this is valid for any τ ∈ Λ d , E([ X(ω)] k ) = 0 if kω = 0 modulo 2π in dimension r.

Property [START_REF] Lustig | Microcanonical monte carlo simulation of thermodynamic properties[END_REF] for ω = 0 and property [START_REF] Nocedal | Numerical Optimization[END_REF] are also proved by decomposing Cov(A, B) = E(AB * ) -E(A) E(B * ) and showing that

E([ X(0)] 2k ) = E(| X(0)|) , E([ X(0)] 2k+1 ) = E( X(0)) = d E(X(u)) . ( 47
)
This results from the fact that X(0) = u∈Λ d X(u) is real so [ X(0)] 2k = | X(0)| and [ X(0)] 2k+1 = X(0). If X is Gaussian then X(ω) and X(ω ) are independent if ω = ω because they are Gaussian random variables with a zero covariance. It results that [ X(ω)] k and [ X(ω )] k are also independent for any (k, k ) ∈ Z 2 and hence have a zero covariance if ω = ω . If ω = ω and k = k then [START_REF] Georgii | Gibbs measures and phase transitions[END_REF] proves that their covariance remains zero. It results that Cov([ X(ω)] k , [ X(ω )] k ) = 0 only if (ω, k) = (ω , k ) and hence that the covariance is diagonalized. If ω = 0 then X(ω) is a zero-mean complex Gaussian random variable whose real and imaginary parts are not correlated. As a result E(| X(ω)|) 2 /E(| X(ω)| 2 ) = π/4.

D Bump Steerable wavelet

We review the bump steerable wavelet introduced in [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF]. It is illustrated in Figure 1 and provides a sparse representation of images with oriented structures. A general way of constructing steerable wavelets is to use the polar coordinate in the Fourier frequency domain [START_REF] Unser | A unifying parametric framework for 2d steerable wavelet transforms[END_REF].

We specify the bump steerable wavelet along the radius |ω| and angle arg(ω) for ω = |ω|e i•arg (ω) . We assume that the number of angles Q is even and its central frequency is λ = (ξ 0 , 0). Its formula is

ψ(ω) = c • exp -(|ω| -ξ 0 ) 2 ξ 2 0 -(|ω| -ξ 0 ) 2 1 [0,2ξ 0 ] (|ω|) • cos Q/2-1 (arg(ω))1 arg(ω)< π 2 ,
where c is a normalization constant. The radial function along |ω| is a bump function which is a compact-support approximation of a Gaussian window. The same angular window function along arg(ω) is used in [START_REF] Simoncelli | The steerable pyramid: a flexible architecture for multi-scale derivative computation[END_REF].

As in [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF], the father wavelet φ is an isotropic Gaussian window function,

φ(ω) = exp - |ω| 2 2σ 2 .
We choose ξ 0 = 1.7π, σ = 0.702 × √ 2 × 2 -0.55 ξ 0 and c = 1.29

-1 2 Q 2 -1 ( Q 2 -1)! ( Q 2 )(Q-2)!
. These hyper-parameters are also used in [START_REF] Mallat | Phase harmonic correlations and convolutional neural networks[END_REF]. For the wavelet transform Wx with oversampling, we obtain numerically the frame constants A W = 2.0 and B W = 4.6 for d = 128 2 , J = 5 and Q = 16. It is therefore a complete and stable representation. These mother and father wavelets also satisfy ψ(u 1 , -u 2 ) = ψ(u 1 , u 2 ) and φ(u 1 , -u 2 ) = φ(u 1 , u 2 ).

E Proof of Theorem 4.1

Observe first that µ 0 is invariant to the action of any g ∈ G. Indeed, the probability measure µ 0 of a Gaussian white noise is uniform on an ball of R d centered in 0. As a result µ 0 is invariant to the action of any linear unitary operator and hence invariant to the action of any g ∈ G which is linear and unitary. Since µ 0 is invariant to the action of g, if the gradient descent is covariant to the action of g then the derivations of Theorem 3.4 in [START_REF] Bruna | Multiscale sparse microcanonical models[END_REF] prove that the probability measure µ t of x t is invariant to g for t ≥ 0. We thus need to show that the L-BFGS gradient-descent algorithm is covariant to the action g ∈ G at each iteration (before stopped by line-search conditions). Let xt and x t be two sequences generated by the gradient descent from the initial conditions x0 and x 0 . The covariance to the action of g means that if x0 = g.x 0 then xt = g.x t for any t ≥ 0. Although invariant properties of a gradient-based algorithm are well studied in the optimization literature [START_REF] Nocedal | Numerical Optimization[END_REF][START_REF] Fletcher | Practical Methods of Optimization, 2nd Edition[END_REF], we give below a prove for completeness.

L-BFGS is a variant of the quasi-Newton method BFGS. It estimates the inverse Hessian matrix of the objective function f based on s t = x t+1 -x t and y t = ∇f (x t+1 ) -∇f (x t ).

Assume there is a limited memory size m to compute an approximation H m t of the Hessian, from s t-a and y t-a for 1 ≤ a ≤ m. At each iteration t, the algorithm chooses an appropriate step-size α t satisfying strong Wolfe conditions [START_REF] Nocedal | Numerical Optimization[END_REF] and updates x t with x t+1 = x t -α t H m t ∇f (x t ).

Since f (x) = K Rx -K Rx 2 E G and K HWg.x = K HWx , it results that f (g.x) = f (x). It implies that ∇f (g.x) = g.∇f (x).

(48)

We prove the covariance property by induction. We suppose that xt = g.x t for t ≤ t in order to prove that xt+1 = g.x t+1 . By definition, xt+1 = xt -αt Hm t ∇f (x t ).

The induction assumption implies that st = g.s t and ỹt = g.y t for t ≤ t. It is therefore sufficient to verify that Hm t = g.H m t .g T , αt = α t , ∀t ≥ 0 Each H m t is defined recursively by

H a t = V T t+a-m-1 H a-1 t V t+a-m-1 + ρ t+a-m-1 s t+a-m-1 s T t+a-m-1 ,
for 1 ≤ a ≤ m and for all t ≥ 0 V t = Id -ρ t s t y T t and ρ t = (s T t y t ) -1 .

We thus have Ṽt = g.V t .g T . By convention, H 0 0 = Id and H 0 t = (s T t-1 y t-1 ) -1 Id for t ≥ 1. We thus conclude that Hm t = gH m t g T . Note that if t ≤ m, the index t+a-m-1 is negative, and in this case V t+a-m-1 = Id and ρ t+a-m-1 = 0.

To verify that αt = α t , we observe that the strong Wolfe-line search conditions are invariant to any linear unitary g ∈ G. Indeed these conditions are invariant to any affine operator [START_REF] Nocedal | Numerical Optimization[END_REF].

F Proof of Proposition 5.1

If p is invariant to an action of G then M HW and K HW are not modified if X is transformed into g.X. Similarly, we verify from (3) and (4) that M HW and K HW are not modified if x is transformed into g.

Figure 1 :

 1 Figure 1: (a): Real-part of steerable bump wavelet ψ(u). (b): Fourier transform ψ(ω).

Figure 2 :

 2 Figure 2: Top: turbulent vorticity field. Bottom: Each image gives the modulus (above) or the phase (below) of wavelet coefficients x ψ λ (u) for different frequency channels λ = 2 -j r -ξ. Large modulus coefficients are shown in black. The columns correspond to different scales and angles (j, ). (a): (j, ) = (1, 0). (b): (j, ) = (2, 0). (c): (j, ) = (1, Q/4). (d): (j, ) = (2, Q/4).

Figure 3 (

 3 b) shows a realization of the maximum entropy Gaussian model X. It is conditioned by wavelet covariances on a foveal graph with ∆ = 2. The wavelet transform is computed with a bump wavelet specified in Appendix D, with J = 5, Q = 16 and d = 256 2 . In this case |E G |/d = 3.6 10 -2 . The covariances are estimated from a single realization of X. The calculation of the Lagrange multipliers in (

Figure 3 :

 3 Figure 3: (a): Realization of a stationary turbulent vorticity field X. (b): Realization of a Gaussian maximum entropy model X calculated from wavelet covariances estimated on a foveal graph. (c): The full line and dashed lines are the power spectrum of X and the power spectrum of X respectively, as a function of the radial frequency |ω| (shown in log-log scale).

Figure 4 :

 4 Figure 4: Power spectrum of [X ψ λ (u)] k for a fixed λ, for the turbulent vorticity flow at the top of Figure 5(a), and different k. The power spectrum is shown in log base 10. (a): k = 0. (b): k = 1. (c): k = 2.

6 10 - 2 .

 2 • Model B with k min = 0 and k max = 1. It includes the covariance of the modulus of wavelet coefficients across angles. Covariance across angles at most distant by π/4 are computed for k, k ∈ {0, 1} by setting ∆ = Q/4. We set ∆ j = 0 and ∆ n = 2, and thus only incorporate covariance across spatial positions. Compared to Model A, spatial correlations are included for k = k = 1 but also for k = k = 0. The relative model size is |E B G |/d = 1.1 10 -1 . • Model C with k min = 0 and k max = 2. It incorporates covariances of the modulus and phase of wavelet coefficients across scales and angles. Neighborhoods are limited by ∆ j = 1, ∆ = Q/4 and ∆ n = 2. Compared to Model B, it incorporates j = j + 1 to capture scale interactions. The phases of wavelet coefficients at a scale 2 j and 2 j = 2 j+1 are correlated with the harmonic exponents (k, k ) = (1, 2). The set of (k, k ) are restricted to k = 0, 1, 2 when k = 0, and k = 1, 2 when k = 1. We use the same spatial correlation range as Model B. The relative model size is |E C G |/d = 1.7 10 -1 . • Model D with k min = 0 and k max = 2 is a rotation invariant version of Model C, with the same ∆ j , ∆ and (k, k ). This model sets ∆ n = 0 and thus does not capture spatial correlations explicitly. The symmetry group G is composed of translations and Q rotations by 2πη/Q. This rotation invariance is represented by computing a Fourier transform along angles and setting to zero the covariance coefficients according to Theorem 5.1. The model size is therefore much smaller with |E D

Figure 5 :

 5 Figure 5: The first row shows a realization x of X. The second to fifth row give realizations of microcanonical models A, B, C, D computed from x. Each column corresponds to a different X. First column: Isotropic turbulent vorticity field. Second column: Non-isotropic turbulent vorticity field. Third column: Bubbles. Fourth column: Multifractal Random Walk.

RxwithK

  Rx by using the same normalization. We define C model Rx by normalizing the K model Rx with D R as in (40). The relative error of the model is defined by model

Figure 6 :

 6 Figure 6: Each graph shows correlation values C R (k, j, a) in full line, as a function of the normalized distance a for a fixed scale 2 j and k. The values of C model R (k, j, a) are shown with dash-dot lines for model A and dashed lines for model D. The value of k and j is given at the top of each graph. The graphs in (a) are computed for isotropic turbulences and in (b) for increments of a Multifractal Random Walk.

  q) 0.04 (0.01) 0.03 (0.02) 0.05 (0.04) 0.14 (0.05) 0.24 (0.05) D st (1, q) 0.01 (0.01) 0.02 (0.03) 0.04 (0.05) 0.06 (0.07) 0.10 (0.10) A st (2, q) 0.03 (0.02) 0.03 (0.02) 0.05 (0.03) 0.08 (0.05) 0.14 (0.07) D st (2, q) 0.01 (0.01) 0.03 (0.03) 0.04 (0.05) 0.06 (0.07) 0.08 (0.09)

)

  By definition, H(z) = {ρ(Real(e iα z))} α∈[0,2π] . Since H(z) -H(z ) 2 = H(z) -H(z ) 2 , we do the proof on H(z) -H(z ) 2 . We prove the theorem by first showing that for a rectifier ρ(a) = max(a, 0) and any (a, a ) ∈ R 2 |ρ(a) -ρ(a )| + |ρ(-a) -ρ(-a )| = |a -a |. Indeed ρ(a) = a and ρ(-a) = 0 or ρ(a) = 0 and ρ(-a) = -a . The equality is verified by considering separately the cases where a and a have same or different signs. For a = 0 we get |ρ(a)| 2 + |ρ(-a)| 2 = |a| 2

3 n=0 2 . 4 |z -z | 2 ,

 3242 |ρ(Real(e i(α+nπ/2) z))| 2 = |z| 2 . (e i(α+nπ/2) z)) -ρ(Real(e i(α+nπ/2) z ))| 2 ≤ |z -z | Integrating both of these double inequalities over α ∈ [0, π/2] (e iα z)) -ρ(Real(e iα z ))| 2 dα ≤1which proves (44).

  x. The proposition results are obtained by relating [g.x ψ λ ] k with [x ψ λ (u)] k and applying this relation to x = X and x = x.(i). If g.x = -x then[g.x ψ λ (u)] k = (-1) k [x ψ λ (u)] k . (49)For v = (λ, k, u) we thus derive that M HW (v) = (-1) k M HW (v) and M HW (v) = (-1) k M HW (v) and thus vanish if k is odd. Similarly, for v = (λ , k , u ) we verify that

Table 1 :

 1 Covariance errors (42) models A,B,C,D compared to the empirical error (41), for the isotropic and non-isotropic vorticity fields.from covariances over a limited set E G ⊂ V 2 0 . Optimizing E G is a trade-off between the variance and bias terms. If the sufficient statistic set

		Isotropic Non Isotropic
	emp	0.58 (0.05)	0.68 (0.09)
	A	0.82 (0.01)	0.80 (0.01)
	B	0.29 (0.02)	0.30 (0.02)
	C	0.24 (0.02)	0.25 (0.03)
	D	0.25 (0.02)	1.6 (0.09)

* This work is supported by the ERC InvariantClass 320959. This work was supported by the PRAIRIE 3IA Institute of the French ANR-19-P3IA-0001 program.

1

(k, j, a) as a function of the normalized distance a, for k = 0, 1 and j = 1, 2, 3. The notion of short and long range correlations is defined relatively to each scale 2 j through the parameter a. For k = 1, it corresponds to correlations of wavelet coefficients, which have a fast decay when the distance a increases. At all scales 2 j the long range correlations of wavelet coefficients for models A and D are close to the correlations obtained with the original turbulence and MRW processes.

For k = 0, Figure 6 gives the long range correlation of the modulus of wavelet coefficients. It still has a fast decay for the turbulence field but a slow decay for the MRW. The model D gives a much better approximation of long range spatial correlations of turbulence than model A. However there is a residual error at the finest scale j = 1 because these foveal models do not capture long range correlations at fine scales. The error is more dramatic for the MRW process where the modulus of wavelet coefficients are correlated over much longer ranges at fine scales. As expected, this evaluation shows that long range correlations at different scales are not fully captured by foveal models which are constrained on a correlation range of the order of 2 j ∆ n at each scale 2 j . One could incorporate these long range spatial covariances by increasing ∆ n but it would increase considerably the model size which is proportional to (2∆ n + 1) 2 .

High order moments Moments of multiple order q can be measured with structure functions at different scale 2 j [START_REF] Van'yan | Structure function of the velocity field in turbulent flows[END_REF][START_REF] Farge | Turbulence analysis, modelling and computing using wavelets[END_REF]. We compute the maximum moment of order q of increments between points (u, u

It is compared with the structure function S model Rx (j, q) calculated for each microcanonical model X, with an error st (j, q) = |S(j, q) -S model Rx (j, q)| |S(j, q)| .

Table 2 reports the mean and standard deviation (in brackets) of st (j, q) for the model A and model D, at two different scales 2 j , for the isotropic turbulence. They are estimated over 10 independent realizations x. As expected, the errors increase with the order q because foveal covariance model only imposes moments of order q = 1 and q = 2. The model D creates no significant error on these high-order moments, since the estimated average error is comparable to the estimator standard deviation. It indicates

The same property applies to M HW (v). Applying (50) to K HW (v, v ) and K HW x(v, v ) also proves that they are real.

G Proof of Theorem 5.1

Since λ = 2 -j r -ξ the index of R = HW can be written v = (j, , k, u). We compute the covariance K HW (v, v ) for v = (j , , k , u ) with u = u. Since the distribution is stationary the covariance value remains the same if we set u = u = 0.

To prove the theorem we first observe that if g.x = r η x then

If X is isotropic then K HW (v, v ) and K HW x(v, v ) are invariant by rotations. The covariance property (51) thus implies that they are a function of the difference of angles -. It is therefore a periodic convolution kernel along angles. It results that the angular dependence is diagonalized by a Fourier transform along the angle variable, which proves (39) for R = F HW.

A central reflection is a rotation by π. Since the distribution of X is invariant to all rotations it is invariant to a central reflection. Property (ii) of Proposition (5.1) implies that K HW is real.

Let us consider the line reflection g relatively to the horizontal axis. Since ψ(u 1 , -u 2 ) = ψ(u 1 , u 2 ) and φ(u 1 , -u 2 ) = φ(u 1 , u 2 ) we verify that for any λ = 2 -j r -ξ and λ = 2 -j r ξ we have

As a result, if λ = 2 -j r -ξ then [(g.x) ψ λ (0)] k = [x ψ λ (0)] k with λ = 2 -j r ξ.

If K HW is invariant to line reflections then (53) implies that K HW is not modified when ( , ) is transformed into (-, -). Since K HW is real and invariant to a change of sign of ( , ) it results that the Fourier coefficients of K F HW are also real.