Charles-Alban Deledalle
email: charles-alban.deledalle@math.u-bordeaux.fr.

Loïc Denis
email: loic.denis@univ-st-etienne.fr

Sonia Tabti
email: sonia.tabti@telecom-paristech.fr

Florence Tupin
email: florence.tupin@telecom-paristech.fr

channel SAR speckle reduction?

Keywords: SAR, speckle, variance stabilization, ADMM, Wishart distribution

MuLoG, or How to apply Gaussian denoisers to multi-

I. INTRODUCTION

Synthetic aperture radar (SAR) imaging is a widely used technique for earth observation. It offers complementary information to the more common optical imaging. Among its distinctive features, one may cite its all-weather and dayand-night imaging capabilities [START_REF] Moreira | A tutorial on synthetic aperture radar[END_REF], the interferometric configurations that give access to 3-D reconstructions for digital elevation models [START_REF] Rabus | The shuttle radar topography mission -a new class of digital elevation models acquired by spaceborne radar[END_REF] and displacement estimation [START_REF] Massonnet | Radar interferometry and its application to changes in the earth's surface[END_REF], [START_REF] Ferretti | Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry[END_REF] or the polarimetric and tomographic modes that give access to estimates of the biomass in forested areas [START_REF] Le Toan | The BIOMASS mission: Mapping global forest biomass to better understand the terrestrial carbon cycle[END_REF].

Coherent combination of several radar echoes within each resolution cell results in interferences and the well-known speckle phenomenon [START_REF] Goodman | Some fundamental properties of speckle[END_REF]. Due to speckle, regions with homogeneous radar properties display strong fluctuations in SAR images. Direct estimation of the reflectivity, the interferometric phase or polarimetric properties is unusable given its prohibitively large variance. Speckle reduction is thus a longstanding topic in SAR imagery.

Since speckle statistics departs from the additive Gaussian noise model widely used for processing optical images, a whole line of denoising methods have been designed specifically for speckle reduction. Restoration of intensity images (i.e., single-channel SAR images) remains the first and most studied case. Many different schemes have been proposed in this context, see the recent reviews [START_REF] Argenti | A tutorial on speckle reduction in synthetic aperture radar images[END_REF] and [START_REF] Deledalle | Exploiting patch similarity for SAR image processing: the nonlocal paradigm[END_REF]. Among the different possible strategies, selection-based methods use various criteria to identify homogeneous collections of neighboring pixels. Pixels may be selected by locally choosing the best neighborhood within a fixed set of oriented windows [START_REF] Lee | Polarimetric SAR speckle filtering and its implication for classification[END_REF]. Another approach consists of selecting connected pixels by region growing [START_REF] Vasile | Intensity-driven adaptiveneighborhood technique for polarimetric and interferometric SAR parameters estimation[END_REF]. Patch-comparison has been shown to provide a robust way to get a (weighted) selection of similar pixels [START_REF] Deledalle | How to compare noisy patches? Patch similarity beyond Gaussian noise[END_REF], [START_REF] Deledalle | Exploiting patch similarity for SAR image processing: the nonlocal paradigm[END_REF]. Variational methods formulate the estimation problem as an optimization problem depending on the whole image. The objective function that is optimized is composed of two terms: a data-fitting term and a regularization term. Due to the non-Gaussian distribution of the intensity of SAR images, the data-fitting term differs from the usual sum of squared differences that arises from a Gaussian assumption. Several regularization terms have been considered in the literature, e.g., total variation (TV) [START_REF] Rudin | Multiplicative denoising and deblurring: Theory and algorithms[END_REF], [START_REF] Shi | A nonlinear inverse scale space method for a convex multiplicative noise model[END_REF], [START_REF] Denis | SAR image regularization with fast approximate discrete minimization[END_REF], [START_REF] Steidl | Removing multiplicative noise by Douglas-Rachford splitting methods[END_REF], [START_REF] Bioucas-Dias | Multiplicative noise removal using variable splitting and constrained optimization[END_REF], curvelets [START_REF] Durand | Multiplicative noise removal using L1 fidelity on frame coefficients[END_REF], Gaussian mixture models [START_REF] Tabti | Modeling the distribution of patches with shift-invariance: application to SAR image restoration[END_REF]. Yet another strategy is to transform the data so the noise becomes additive and approximately Gaussian and stationary by using a homomorphic transform. As further discussed in Sec. II, the log-transformed intensity is approximately Gaussian distributed and can thus be restored using a Gaussian denoiser, e.g., based on wavelet thresholding [START_REF] Xie | SAR speckle reduction using wavelet denoising and Markov random field modeling[END_REF], [START_REF] Achim | SAR image denoising via Bayesian wavelet shrinkage based on heavy-tailed modeling[END_REF], [START_REF] Bhuiyan | Spatially adaptive wavelet-based method using the Cauchy prior for denoising the SAR images[END_REF] or on patch redundancy [START_REF] Mäkitalo | Denoising of single-look SAR images based on variance stabilization and nonlocal filters[END_REF].

Although many applications are based on the intensity provided by a single SAR image (corresponding to the square of the modulus of the backscattered electromagnetic field), the SAR data can carry much more information. Indeed, the phase of a pixel carries information about the acquisition geometry. Combining two SAR acquisitions with a slightly different angle allows retrieving information on the elevation of the points and corresponds to the so-called interferometric mode. Processing this information can be done by considering the vector of the two complex SAR signals and their associated empirical covariance matrix. Another information is provided by the polarization of the emitted and backscattered wave. By emitting and receiving the signal with different polarizations (for instance vertical and horizontal polarizations), a vector of 4 complex values is recorded for each pixel. This vector gives valuable information on the backscattering mechanisms happening inside the pixel (double bounce scattering, volume scattering, etc.). This polarimetric mode can also be combined with the interferometric one, giving a D-dimensional complex vector for each pixel and allowing to retrieve the height of different scatterers inside the pixel. Finally, the tomographic mode exploits a set of interferometric acquisitions to recover the elevation profile inside the pixel. These multi-channel SAR data are of high interest for a wide range of applications for Earth observation. But as for the intensity data, the speckle phenomenon induces a high variability of the physical parameters (interferometric phase, coherence, inter-channel crosscorrelations, etc.).

Not all the aforementioned speckle reduction methods generalize well to multi-channel SAR images. Selection-based methods require extending the homogeneity or similarity criteria to multi-variate data. This can be done by using only part of the information, e.g., the span [START_REF] Vasile | Intensity-driven adaptiveneighborhood technique for polarimetric and interferometric SAR parameters estimation[END_REF] or the scattering properties [START_REF] Lee | Scattering-model-based speckle filtering of polarimetric SAR data[END_REF], or by exploiting the whole covariance matrices [START_REF] Chen | Nonlocal filtering for polarimetric SAR data: A pretest approach[END_REF], [START_REF] Deledalle | NL-SAR: A unified nonlocal framework for resolution-preserving (Pol)(In) SAR denoising[END_REF]. Once relevant pixels have been selected, estimation of polarimetric/interferometric information is straightforward, e.g., by using a linear mean square error approach [START_REF] Lee | Polarimetric SAR speckle filtering and the extended sigma filter[END_REF], [START_REF] Deledalle | NL-SAR: A unified nonlocal framework for resolution-preserving (Pol)(In) SAR denoising[END_REF] or a weighted maximum likelihood estimator [START_REF] Deledalle | Iterative weighted maximum likelihood denoising with probabilistic patch-based weights[END_REF], [START_REF] Deledalle | NL-InSAR: Nonlocal interferogram estimation[END_REF]. Extension of variational methods to multi-channel SAR data is more challenging. A simple remedy is to consider only restoring the diagonal elements [START_REF] Li | PolSAR image speckle reduction algorithm based on TV-PDE[END_REF] at a price of information loss about cross-channel correlations. Direct formulation of the objective function on the full covariance matrices raises several problems: (i) computational complexity due to the nonconvexity of the data-fitting term; (ii) difficulty to express regularity properties of the complex-valued terms of the covariance matrices; (iii) non-stationary variance of speckle that leads to over/under-smoothing in some areas. These difficulties explain that very few works were conducted in this direction, with the exception of recent works on multi-channel TV regularizations [START_REF] Nie | A variational model for PolSAR data speckle reduction based on the Wishart distribution[END_REF], [START_REF] Nie | A nonlocal TV-based variational method for PolSAR data speckle reduction[END_REF], [START_REF] Nie | A new algorithm for optimizing TV-based PolSAR despeckling model[END_REF]. Finally, while the homomorphic transform approach is well understood on intensity images, no variance stabilization transform is known for multi-channel images.

A generic methodology to apply denoising methods from the "additive Gaussian noise world" to multi-channel SAR data is definitely lacking. This paper attempts to fill this gap. The contributions of this paper are the following:

1) to provide a generic method, called MuLoG, to embed a Gaussian denoiser in a speckle reduction filter; 2) to apply as well to single-channel or to multi-channel SAR images; 3) to produce estimates of the complex covariance matrices that capture all polarimetric and/or interferometric information, 4) to better account for speckle statistics than other methods based on a homomorphic transform1 ; 5) to require no parameter tuning other than possibly within the Gaussian denoiser.

We introduce our generic methodology by first considering (single-channel) intensity images (Sec. II), then the extension to multi-channel SAR images (Sec. III). We discuss implementation issues (Sec. IV) before illustrating the proposed methodology with several Gaussian denoisers and different types of multi-channel SAR images (Sec. V).

II. SPECKLE REDUCTION FOR SAR INTENSITY IMAGES

A. Statistics of univariate SAR images a) Intensity: Univariate SAR images are by nature complex-valued and only the modulus (a.k .a., the amplitude) is informative. The square of the modulus (a.k .a., the intensity) is nevertheless easier to manipulate and, according to Goodman's model [START_REF] Goodman | Some fundamental properties of speckle[END_REF], it follows a gamma distribution G(R; L) with a probability density given by

p I (I|R) = L L I L-1 Γ(L)R L exp -L I R , (1)
where I ∈ R + is the observed intensity, R ∈ R + is the underlying reflectivity (related to the radar cross-section), L > 0 is the number of looks, and Γ the gamma function.

Note that we denote by R + the set of positive real values. The intensity I can be decomposed as a product of the reflectivity R and of a speckle component S distributed under a standard gamma distribution (S ∼ G(1; L)):

I = R × S, E[I] = R and Var[I] = R 2 L . (2)
As the variance depends on the expectation, fluctuations are said to be signal dependent. The top graph in Fig. 2(a) illustrates how a simple rectangle signal (solid curve) gets corrupted by speckle: the gray area shows values between the first and third quartiles, the dots represent the expectation and the dashed line a single noisy realization. The signaldependent nature of the fluctuations can be observed: the difference between first and third quartiles is larger when the underlying signal values are high. Last but not least, the gamma distribution has a heavier right-tail characterizing the typical bright outliers observed in SAR intensity images. b) Logarithm: The log-transform y = log I ∈ R is often employed to convert multiplicative fluctuations to additive ones. From the gamma distribution (1), y follows the Fisher-Tippett distribution defined as

p y (y|x) = L L Γ(L) e L(y-x) exp -Le y-x , (3)
where x = log R ∈ R. The Fisher-Tippett distribution, denoted by FT (x; L), models additive corruptions as [START_REF] Xie | Statistical properties of logarithmically transformed speckle[END_REF]

y = x + s , (4)
E[y] = x -log L + Ψ(L) , (5) and Var
[y] = Ψ(1, L) , (6)
where s ∼ FT (0; L), Ψ(•) is the digamma function and Ψ(•, L) is the polygamma function of order L. The log transform stabilizes the variance, i .e., the fluctuations are made signal independent. Equation [START_REF] Bellman | Introduction to matrix analysis[END_REF] shows that the noise has a non-zero mean. If noise is assumed zero-mean during a speckle-reduction step performed on log-transformed data, a subsequent debiasing step is then necessary. The bottom graph in Fig. 2(a) displays the intensities in decibels (i.e., 10 • log 10 of the intensity). By comparing with the top curve, variance stabilization can be observed, as well as the bias between the expectation of log-transformed intensities (black dots) and ground truth signal (solid curve).

B. Homomorphic approach

We now consider I and R as images, such that I ∈ R n + , and R ∈ R n + , where n is the number of pixels. Let y ∈ R n and x ∈ R n be the entry-wise logarithm of I and R respectively, i .e., such that y k = log I k and x k = log R k . A classical approach to estimate R is thus to approach the Fisher-Tippett distribution by a non-centered Gaussian distribution, which leads to the following estimation procedure

x = f Ψ(1,L) (y) + (log L -Ψ(L))1 n debiasing , (7)
where f σ 2 : R n → R n is a denoiser for images contaminated by zero-mean additive white Gaussian noise N (0; σ 2 Id). Typically, f can be a regularized least-square solver expressed as

f σ 2 (y) ∈ argmin x∈R n 1 2σ 2 ||y -x|| 2 + R(x) , (8)
where, within the Maximum A Posteriori (MAP) framework, R(x) =log p x (x) is a prior term enforcing some regularity on the solution. Finally, the estimate R is defined as Rk = exp xk for all pixel index k. This process is summarized on Fig. 1, top row. We illustrate the restored 1D signals obtained by applying L2+TV minimization [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] within a homomorphic procedure on Fig. 2(b).

The homomorphic approach has been extensively used, e.g., for wavelet prior in [START_REF] Xie | SAR speckle reduction using wavelet denoising and Markov random field modeling[END_REF], for patch-based priors as in KSVD [START_REF] Foucher | SAR image filtering via learned dictionaries and sparse representations[END_REF] or for non-local filtering [START_REF] Mäkitalo | Denoising of single-look SAR images based on variance stabilization and nonlocal filters[END_REF]. While for large values of L the Fisher-Tippett distribution approaches the Gaussian distribution, the two distributions differ significantly for low values of L. In particular, the Fisher-Tippett is asymmetrical with a heavier left-tail. Hence, this approximation often leads to remaining dark stains on the resulting images (see also Fig. 6(d)).

C. Variational approach

Rather than applying the homomorphic approach in order to reduce the original problem to the Gaussian case, an alternative proposed in [START_REF] Aubert | A variational approach to removing multiplicative noise[END_REF], [START_REF] Denis | SAR image regularization with fast approximate discrete minimization[END_REF] is to consider directly the gamma distribution of I ∈ R n + leading to a MAP solver of the form (Fig. 1 second row):

R ∈ argmin R∈R n + -log p I (I|R) + R(R) , (9)
where

-log p I (I|R) = L n k=1 log R k + I k R k + Cst.
Cst. denoting a constant for the optimization problem. Not only is this minimization constrained to positive-valued images, but the objective function is also nonconvex. This nonconvexity enforces some robustness to the bright outliers originating from the right tail of the gamma distribution. As a consequence, when standard iterative solvers are used (e.g., gradient descents [START_REF] Aubert | A variational approach to removing multiplicative noise[END_REF], or the forward-backward algorithm [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF]), the solution depends on the initialization and the internal parameters of the solver, even if R is chosen convex. Convexification by replacing the original objective function by its convex hull [START_REF] Kang | Two-level convex relaxed variational model for multiplicative denoising[END_REF] simplifies the minimization at the cost of a loss of accuracy of the statistical model. In particular, the convex hull does not capture the right tail of the gamma distribution, leading to remaining bright outliers on the resulting images. For some Markov prior regularization, namely convex pairwise regularizations, the global optimum can be obtained in finite time [START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF]. This optimization method has been applied in the 1D illustration of Fig. 2(c). It can be observed that, since large fluctuations in bright areas are more strongly penalized than the small fluctuations in lowlevel areas, speckle is reduced predominantly in the brightest areas. In the restored signal, the remaining noise variance is made uniform, i.e., signal-independent, in linear scale. Thus, small details with an identical signal to noise ratio will be more likely suppressed if placed in areas of large average value. This phenomenon is also observed on images: speckle noise is reduced more strongly in brighter areas [START_REF] Denis | SAR image regularization with fast approximate discrete minimization[END_REF], [START_REF] Steidl | Removing multiplicative noise by Douglas-Rachford splitting methods[END_REF], and bright punctual targets are found to be spread out.

D. Variational approach on log-transformed data

Another alternative proposed in [START_REF] Shi | A nonlinear inverse scale space method for a convex multiplicative noise model[END_REF], [START_REF] Steidl | Removing multiplicative noise by Douglas-Rachford splitting methods[END_REF], [START_REF] Bioucas-Dias | Multiplicative noise removal using variable splitting and constrained optimization[END_REF], [START_REF] Durand | Multiplicative noise removal using L1 fidelity on frame coefficients[END_REF] consists of formulating the estimation problem in the log domain. Rather than approximating the likelihood of log-transformed reflectivities by a Gaussian distribution (as done by homomorphic approaches discussed in Sec. II-B), the Fisher-Tippett distribution (3) is considered. The regularization is also expressed on the log-transformed reflectivities x ∈ R n , leading to a MAP optimization problem of the form:

x ∈ argmin x∈R n -log p y (y|x) + R(x) , (10)
wherelog p y (y|x) = L n k=1

x k + e y k -x k + Cst.

The final estimate is defined as Rk = exp xk for all pixel index k. Note that even though the solutions of problems ( 9) and ( 10) are different, their definitions only differ in terms of the prior regularization since log p y (y|x) = log p I (exp(y)| exp(x)) + Cst.. Nevertheless, this change of variable leads to several advantages compared to [START_REF] Chan | Plug-and-play ADMM for image restoration: Fixed point convergence and applications[END_REF]. First, the optimization is unconstrained as x can be any vector of R n . More importantly, the data fidelity becomes convex. As a consequence, if R is also chosen convex, solutions depend neither on the initialization nor the choice of the internal parameters of the solver. Such internal parameters only affect the speed of convergence and thus the number of iterations to perform in practice.

The multiplicative image denoising by augmented Lagrangian (MIDAL) algorithm [START_REF] Bioucas-Dias | Multiplicative noise removal using variable splitting and constrained optimization[END_REF] considers a convex regularization, TV, and minimizes (10) using the alternating direction method of multipliers (ADMM) algorithm [START_REF] Glowinski | Sur l'approximation, par éléments finis d'ordre un, et la résolution, par pénalisation-dualité d'une classe de problèmes de Dirichlet non linéaires[END_REF], [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite element approximation[END_REF], [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF] that repeats, for an internal parameter β > 0, the updates

ẑ ← argmin z β 2 ||z -x + d|| 2 + R(z) , (11)
d ← d + ẑ -x , (12)
x ← argmin

x β 2 ||x -ẑ -d|| 2 -log p y (y|x) . (13)
Clearly, the minimization for z in [START_REF] Combettes | Signal recovery by proximal forwardbackward splitting[END_REF] depends on the choice of R and can be solved by a dedicated regularized least square solver (corresponding to a Gaussian denoiser). This Fig. 1. Several approaches for speckle reduction discussed in this paper. The proposed scheme generalizes the variational approaches based on a log-transform to multi-channel SAR images. Similarly to homomorphic or variational approaches on log-transformed data, it embeds a Gaussian denoising step.

is schematized on Fig. 1, third row. The minimization for x in (13) amounts to solving n separable convex problems of the form

argmin x k β 2 (x k -a k ) 2 + L x k + e y k -x k , (14)
where a k = ẑk + dk . The explicit solution is given by the Lambert W functions [START_REF] Bioucas-Dias | Multiplicative noise removal using variable splitting and constrained optimization[END_REF], or can be computed more efficiently with a few iterations of Newton's method as

xk ← xk - β(x k -a k ) + L(1 -e y k -x k ) β + Le y k -x k . ( 15
)
Using ten iterations is usually enough to offer good performances within a reasonable computational time, see, e.g., [START_REF] Tabti | Modeling the distribution of patches with shift-invariance: application to SAR image restoration[END_REF].

As already mentioned, when R is convex, the parameter β only acts on the speed of convergence. Interestingly, as the noise variance is independent of the signal, the convergence is in practice uniform meaning that for a finite number of iterations the same amount of smoothing will be performed both in dark and bright regions. Small details with identical signal to noise ratio will be identically smoothed whatever the average value of the area.

When R is nonconvex but satisfies some weak conditions and β is chosen large enough, ADMM still converges to a local minimum [START_REF] Hong | Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems[END_REF]. In this case, the solution depends on both the initialization and the choice of β. We observe that choosing β as (1 + 2/L)Var[y] -1 , where Var[y] = Ψ(1, L), provides a similar smoothing whatever the number of iterations and the initial number of looks L.

This shows that taking the logarithm not only makes the data fidelity term convex, but also ensures a uniform speed of convergence of a non-constrained optimization problem. These are two key practical ingredients that challenge nonconvex strategies that directly deal with gamma distributed values.

We illustrate on Fig. 2(d) that minimization of the total variation of log-transformed intensities produces a smaller loss of contrast of bright structures and reduces speckle equally in all regions. Fig. 2. Illustration of speckle reduction by 1D total variation minimization: (a) speckle-corrupted rectangle signal (L = 2), (b) restored reflectivity using the homomorphic approach and "L2 + TV" minimization in log-domain, (c) restored reflectivity using a variational approach "gamma + TV" in the linear domain, (d) restored reflectivity using a variational approach "Fisher-Tippett + TV" in the log domain. Expectation and quartiles are computed from 50 000 noisy realizations. Restorations computed by exact discrete optimization [START_REF] Ishikawa | Exact optimization for Markov random fields with convex priors[END_REF] with 2 000 quantization levels.

III. EXTENSION TO MULTI-VARIATE SAR IMAGERY A. Statistics of covariance matrices

As mentioned in the introduction, multi-variate complex SAR images carry much more information than univariate SAR images since inter-channel cross-correlations capture geophysical features (e.g., height is related to the interferometric phase, geometrical configuration with the polarimetric properties). All relevant information of these complex-valued images can be gathered together at each pixel by forming a D × D complex covariance matrix C, with D = 2 in singlebaseline SAR interferometry, D = 3 in SAR polarimetry, D = 6 in PolInSAR (polarimetry + interferometry), and even larger values of D for multi-baseline interferometry or SAR tomography [START_REF] Moreira | A tutorial on synthetic aperture radar[END_REF]. Let us denote by k i the vector of D complex values recorded for a pixel i, then the empirical covariance matrix C is given by:

C = 1 L L t=1 k t . k * t
where * denotes the Hermitian transpose. Note that the L vectors k t are often extracted from a small window centered around the pixel of interest leading to a loss of resolution.

If spatial resolution must be preserved, a single-look setting should be used, i .e., L = 1, resulting to a matrix C with rank 1. Goodman's model [START_REF] Goodman | Some fundamental properties of speckle[END_REF] describes speckle in this matrix as being circular complex Wishart distributed, for L ≥ D:

p C (C|Σ) = L LD |C| L-D Γ D (L)|Σ| L exp -L tr(Σ -1 C) , ( 16
)
where Σ is the underlying covariance matrix encoding reflectivities and complex correlations, and L ≥ D is the number of looks. Both C and Σ belong to the open cone of complex Hermitian positive definite matrices. While fluctuations in univariate SAR images are multiplicative, the Wishart distribution, denoted by W(Σ; L) models fluctuations in multivariate SAR images as

C = Σ 1/2 SΣ 1/2 , (17)
E[C] = Σ , (18)
Var[C ij ] = 1 L Σ ii Σ jj , (19)
where S ∼ W(Id; L) (see, e.g., [START_REF] Maiwald | Calculation of moments of complex Wishart and complex inverse Wishart distributed matrices[END_REF]) 2 . Note that the variance for off-diagonal elements does not depend on Σ ij but on Σ ii and Σ jj , which indicates that the fluctuations are not only intra-channel signal dependent but inter-channel signal dependent. Interestingly, according to [START_REF] Goodman | The distribution of the determinant of a complex Wishart distributed matrix[END_REF], [START_REF] Maiwald | Calculation of moments of complex Wishart and complex inverse Wishart distributed matrices[END_REF], we have the following relations regarding the determinant and the trace

|C| = |Σ||S| , (20)
tr C = tr(ΣS) , (21)
E[tr C] = tr Σ , (22)
Var[tr C] = 1 L tr Σ 2 , (23)
and

E[tr(C 2 -Σ 2 )] = 1 L (tr Σ) 2 . (24)

B. Limit of a direct variational approach

As for the univariate case, a variational approach considering the Wishart distribution of C can be expressed as

Σ ∈ argmin Σ H ≻0 -log p Σ (Σ|C) + R(Σ) , (25)
where

-log p Σ (Σ|C) = L n k=1 log |Σ k | + tr(Σ -1 k C k ) + Cst.
where Σ H ≻ 0 reads as "Σ is Hermitian positive definite". Note that (25) boils down to (9) for D = 1. While estimating directly the reflectivity in a variational framework in the case of univariate data requires to enforce a positivity constraint, optimizing for Σ requires optimizing on the open cone of complex Hermitian positive definite matrices, which is much more challenging. Moreover, as noted for the univariate case, the neg-log-likelihood associated to the Wishart distribution is highly nonconvex, so that finding a good quality local optimum is very difficult. This approach is summarized on the second row of Fig. 1.

To circumvent the difficulty arising from the nonconvexity of the neg-log-likelihood, Nie et al. approximated it by its convex envelope [START_REF] Nie | A variational model for PolSAR data speckle reduction based on the Wishart distribution[END_REF], [START_REF] Nie | A nonlocal TV-based variational method for PolSAR data speckle reduction[END_REF]. In their first algorithm, WisTV [START_REF] Nie | A variational model for PolSAR data speckle reduction based on the Wishart distribution[END_REF], R was chosen as a matricial total-variation regularizer,

K :    α 1 . . . α D 2    →        α 1 (α D+1 + jα D+2 )/ √ 2 . . . (α D 2 -1 + jα D 2 )/ √ 2 (α D+1 -jα D+2 )/ √ 2 α 2 . . . . . . . . . (α D 2 -1 -jα D 2 )/ √ 2 . . . α D        (26)
and next in [START_REF] Nie | A nonlocal TV-based variational method for PolSAR data speckle reduction[END_REF] it was replaced by a matricial non-local totalvariation regularizer [START_REF] Gilboa | Nonlocal operators with applications to image processing[END_REF] to better preserve textures. Replacing the neg-log-likelihood by its convex envelope is a rather crude approximation since it leads to underestimating the right tail of the distribution, and thus prevents from being robust against bright outliers. The noise components in C are intensively signal dependent, so that speckle suppression is not as strong in all regions. Last but not least, as for the univariate case, it has been recently observed that the performance of WisTV can drop significantly when the dynamic range of the input SAR image becomes very high [START_REF] Nie | A new algorithm for optimizing TV-based PolSAR despeckling model[END_REF]. In this latter article, the authors solve this issue by normalizing each input matrix C k by a factor τ k defined as a non-linear sigmoid function of their span k . In addition, they reformulate (25) into an unconstrained optimization problem by writing

C k = Π k Π *
k + ǫId, for a fixed ǫ > 0, and then optimizing for Π k ∈ R D×D . The factorization for Π k being non unique, the resulting energy will necessarily present even more local minima. To cope with these different issues we have chosen a different path. We suggest mimicking the univariate case by making use of the matrix logarithm, thereby extending MIDAL approach [START_REF] Bioucas-Dias | Multiplicative noise removal using variable splitting and constrained optimization[END_REF] to multi-channel SAR images.

C. The Wishart-Fisher-Tippett distribution

Our objective is to generalize the use of the log transform in univariate SAR images to multi-variate SAR images. To that end, we resort to the matrix logarithm defined 3,4 as

Σ → Σ = log Σ = E diag( Λ)E -1 where Λi = log Λ i , (27)
E ∈ C D×D is the matrix whose column vectors are eigenvectors (with unit norm) of Σ, Λ ∈ R D + is the vector of corresponding eigenvalues, such that Σ = E diag(Λ)E -1 , and Λ ∈ R D . Its inverse transform is the matrix exponential defined similarly as

Σ → Σ = e Σ = E diag(Λ)E -1 where Λ i = exp Λi . ( 28
)
While Σ lies in the open cone of complex Hermitian positive definite matrices, Σ lies in the vector space of complex Hermitian matrices which is isomorphic to R D 2 . The change

3 the precise definition is Σ = ∞ n=1
Σn !n but is equivalent to [START_REF] Friedland | Convex spectral functions[END_REF] for Hermitian positive definite matrices. 4 the matrix logarithm is not to be confused with the log-determinant function Σ → log |Σ|, the latter equals to the trace of the former.

of variables C = log C and Σ = log Σ leads to the distribution of log-transformed matrices C:

p C ( C| Σ) = L LD f ( C) Γ D (L) e L tr( C-Σ) exp -L tr(e C e -Σ) (29
E[tr C] = tr Σ + D i=1 Ψ(0, L -i + 1) -D log L , (30)
and

Var[tr C] = D i=1 Ψ(1, L -i + 1) , (32)
where S ∼ WFT (0; L). This shows that the trace of the matrix logarithm suffers from additive non-zero-mean signalindependent corruptions. Note that (30) is a direct consequence of [START_REF] Deledalle | NL-SAR: A unified nonlocal framework for resolution-preserving (Pol)(In) SAR denoising[END_REF]. If follows that the D 2 channels of C can reasonably be assumed to have approximately a stabilized variance (see Sec. III-G for numerical evidence), which opens the door to regularization with iterative schemes. Note that we cannot use the matrix log transform directly as a variance stabilization procedure, as done in the univariate homomorphic case, because we do not have an inversion formula (i.e., a bias correction formula). We will thus adopt a variational strategy.

D. Log-channel decomposition

As mentioned in the previous paragraph, Σ lies in the vector space of complex Hermitian matrices that is isomorphic to R D 2 . In this section, we describe a re-parameterization of the log-transformed covariance matrix Σ as a vector of D 2 reals denoted x. We first define in (26) a unitary transform K that maps real-valued vectors α of R D 2 to Hermitian D × D matrices. We also introduce a whitening affine map (A, b) so that the D 2 channels of x are better decorrelated, and a scaling transform Φ (i .e., a diagonal matrix) to balance the variance of noise between channels. The transform Ω between (log-transformed) covariance matrices Σ and the vector of parameter x is:

Σ = Ω(x) = K(AΦx + b) . (33)
We denote its inverse

Ω -1 ( Σ) = Φ -1 A -1 (K -1 ( Σ) -b)
and introduce the linear operator Ω * = ΦA * K * , where * denotes the adjoint. The re-parameterized log-transformed data are noted y ∈ R D 2 and defined by y = Ω -1 ( C), see Fig. 1. The affine transform (A, b) is chosen such that the channels x i of the re-parameterized vector are decorrelated and can then be assumed close to independent:

p Σ( Σ) = p x (x) ≈ D 2 i=1 p x i (x i ) . (34)
The affine transform (A, b) is obtained by principal component analysis, as described in Sec. III-G. Moreover, the scaling matrix Φ is chosen such that the variance of y i is equal to 1 for all channels i, see Sec. III-G.

E. Proposed variational approach

We now extend x and y into images of n real-valued vectors each of dimension D 2 . Writinglog p x i (x i ) = R(x i ) and using the relations [START_REF] Hong | Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems[END_REF] within the MAP framework leads to the following minimization problem

x ∈ argmin x -log p y (y|x) + D 2 i=1 R(x i ) , (37)
where, from (29), we have

-log p y (y|x) = L n k=1
tr Ω(x k ) + e Ω(y k ) e -Ω(x k ) + Cst. ,

and the final estimate Σk at pixel index k is defined as Σk = exp Ω(x k ). Remark that this problem boils down to (10) when D = 1. Note also the difference in notations between x k a vector of D 2 coefficients at pixel k and x i a scalar image of n pixels corresponding to the i-th channel of x.

While the direct multivariate variational approach (25) requires optimizing on the cone of complex Hermitian matrices, a first major advantage is that Problem (37) is unconstrained on R D 2 . Even though the likelihood term [START_REF] Knaus | Dual-domain image denoising[END_REF] is still nonconvex, it appears to be much more suitable for optimization with less local minima. Figure 3 illustrates the evolution of this term (with D = 2, and n = 1) on three two-dimensional cross sections of R 4 showing that for each of such sections the energy appears to be convex. Nevertheless, note that convexity in low-dimensional spaces does not necessarily guarantee convexity in higher dimensional spaces. Unlike the direct multivariate variational approach ( 25), [START_REF] Knaus | Dual-domain image denoising[END_REF] appears to be convex in many scenarios, e.g., in the univariate case [START_REF] Chen | Nonlocal filtering for polarimetric SAR data: A pretest approach[END_REF], and more generally when xk is restricted to the vectorial subspace where Ω(x k ) commutes with Ω(y k ), hence satisfying e Ω(y k ) e -Ω(x k ) = e Ω(y k -x k ) . Convexity, in this case, follows from the fact that the trace is linear, hence convex, the function xk → tr e Ω(-x k +y k ) is a convex spectral function [START_REF] Friedland | Convex spectral functions[END_REF], [START_REF] Lewis | The convex analysis of unitarily invariant matrix functions[END_REF], and Ω is affine.

As for the mono-dimensional case, we will thus consider the ADMM algorithm which iteratively performs the updates

ẑ ← argmin z β 2 ||z -x + d|| 2 + D 2 i=1 R(z i ) , (39)
d ← d + ẑ -x , (40)
x ← argmin

x β 2 ||x -ẑ -d|| 2 -log p y (y|x) . (41)
Since the noise variance is approximately stabilized and close to 1 on each channel (see, III-G), a single value for β can be chosen for all channels. With such a choice, the same amount of smoothing is reached in practice in all regions after a finite number of iterations. Equation ( 39) is separable on the different channels, and can then be solved by applying D 2 times a regularized leastsquare solver. Equation [START_REF] Lee | Polarimetric SAR speckle filtering and its implication for classification[END_REF] amounts to solving n separable problems of the form

argmin x k β 2 ||x k -a k || 2 + L tr Ω(x k ) + e Ω(y k ) e -Ω(x k ) (42)
where a k = ẑk + dk . As for the univariate case, we will consider Newton's method to solve this optimization problem. Its gradient is given by

β(x k -a k ) + LΩ * Id - 1 0 e (u-1)Ω(x k ) e Ω(y k ) e -uΩ(x k ) du , ( 43
) ∆ k,i ∞ = β(x i k -a i k ) + L Ω * Id - 1 0 e (u-1)Ω(x k ) e Ω(y k ) e -uΩ(x k ) du i β + LΩ * 1 0 e (u-1)Ω(x k ) e Ω(y k ) e -uΩ(x k ) du i . ( 35
) ∆ k,i Q = β(x i k -a i k ) + L Ω * Id - 1 Q Q q=1 e (uq-1)Ω(x k ) e Ω(y k ) e -uqΩ(x k ) i β + LΩ * 1 Q Q q=1 e (uq-1)Ω(x k ) e Ω(y k ) e -uqΩ(x k ) i , u q = q -1 2 Q . ( 36
)
Dimension D see the development in Appendix A. By mimicking the univariate case, we consider the following approximation for the Hessian

∆ 0 ∆ 1 ∆ 2 ∆ 4 ∆ 8 ∆ 16
diag βId + LΩ * 1 0
e (u-1)Ω(x k ) e Ω(y k ) e -uΩ(x k ) du , [START_REF] Li | PolSAR image speckle reduction algorithm based on TV-PDE[END_REF] leading to the quasi-Newton iteration

xi k ← xi k -∆ k,i ∞ , (45)
where ∆ k,i ∞ is defined in eq. ( 35). Note that the expression becomes the exact Newton iteration when Ω(x k ) and Ω(y k ) are commuting matrices. As for the univariate case, we notice that ten iterations is enough to offer good performance with a reasonable computational time.

The update matrix ∆ k,i ∞ , defined in eq. ( 35), requires numerical integration. We use Riemann integral approximation with Q rectangles as defined in eq. [START_REF] Kang | Two-level convex relaxed variational model for multiplicative denoising[END_REF]. Of course, the smaller Q, the faster the iteration [START_REF] Liu | Noise estimation from a single image[END_REF]. Figure 4 shows the evolution of ||∆ 100 ||/||x|| with respect to D (with L = D and β = 10L) after 10 iterations of Newton's method where ∆ ∞ is substituted by ∆ Q , for Q = {1, 2, 4, 8, 16}. We furthermore define ∆ 0 by substituting the integral by exp(Ω(y kx k )) as a crude approximation of the gradient (and exact in the case where Ω(x k ) commutes with Ω(y k )). Ideally, the algorithm would reach the optimum solution such that ∆ ∞ ≈ ∆ 100 = 0 (as for the univariate case). This error is quite small whatever Q and D. For low values of D, the error is reduced by about a factor 10 each time Q is multiplied by 4. The increase of the error with D is due to an accumulation of numerical errors of the exponential matrix function, all the more important when Q gets larger. This experiment reveals that Q = 1 is a good choice in practice, reaching a relative error of about 2 • 10 -3 whatever D ∈ [2, 16] while requiring a reasonable computation time.

F. Adaptation to advanced filters

As in the univariate case, we notice that ADMM converges when the regularizer R is chosen as nonconvex. More remarkably, as observed in many contexts, e.g., [START_REF] Wen | Iterative algorithms based on decoupling of deblurring and denoising for image restoration[END_REF], [START_REF] Danielyan | BM3D frames and variational image deblurring[END_REF], [START_REF] Venkatakrishnan | Plug-and-play priors for model based reconstruction[END_REF], [START_REF] Dar | Postprocessing of compressed images via sequential denoising[END_REF], [START_REF] Rond | Poisson inverse problems by the plug-and-play scheme[END_REF], [START_REF] Brifman | Turning a denoiser into a superresolver using plug and play priors[END_REF], replacing the minimization problem in [START_REF] Le Toan | The BIOMASS mission: Mapping global forest biomass to better understand the terrestrial carbon cycle[END_REF] by the solution of a Gaussian denoiser -a scheme known as plug-and-play -leads to appealing results. The resulting algorithm, coined MUlti-channel LOgarithm with Gaussian denoising (MuLoG), is given as

ẑi ← f β -1/2 (x i -di ), for i = 1, . . . , D 2 , ( 46
)
d ← d + ẑ -x , (47)
x ← argmin

x β 2 ||x -ẑ -d|| 2 -log p y (y|x) , (48)
where f σ 2 is again a denoiser for images contaminated by zero-mean white Gaussian noise N (0; σ 2 Id), and ( 48) is solved as mentioned in the previous paragraph. This plug-andplay ADMM algorithm5 is proven to converge6 [START_REF] Chan | Plug-and-play ADMM for image restoration: Fixed point convergence and applications[END_REF] as soon as f σ 2 is a bounded denoiser, i.e., satisfying for all x

||f σ 2 (x) -x|| 2 ≤ nσ 2 C (49)
for some constant C independent of n and σ. As mentioned in [START_REF] Chan | Plug-and-play ADMM for image restoration: Fixed point convergence and applications[END_REF], this is a weak condition which can be expected from most denoisers. In particular, this condition implies that f σ 2 should tend to the identity when σ tends to 0. We depict this algorithm on the last row of Fig. 1.

G. Calibration

In practice, the operators A, b and Φ are obtained from the log matrix of the input image C with n pixels. Let {α 1 , . . . , α n } be the collection of n vectors extracted from the matrix logarithm of the input image as α k = K -1 ( Ck ). We thus define b = 1 n n k=1 α k and A the matrix whose columns are the eigenvectors (with unit norm) of

1 n n k=1 (α k - b)(α k -b) t .
As a result, the vector A -1 (α kb) is the representation of α k in the principal component analysis of {α 1 , . . . , α n } known to maximize inter-channel decorrelation, and thus leads to [START_REF] Hong | Convergence analysis of alternating direction method of multipliers for a family of nonconvex problems[END_REF]. By assuming the noise variance is stabilized on each channel of A -1 (αb), we define the diagonal matrix Φ as diag(σ 2 1 , . . . , σ2 D ), where σ2 1 , . . . , σ2 D are obtained by estimating separately the variance of the noise on each channel of A -1 (αb). This procedure ensures that the noise variance is about 1 on each channel of y. In practice, we resort to the median absolute deviation estimate (MAD) [START_REF] Hampel | The influence curve and its role in robust estimation[END_REF] for this task, but other strategies could be used, such as [START_REF] Liu | Noise estimation from a single image[END_REF], [START_REF] Liu | Single-image noise level estimation for blind denoising[END_REF], [START_REF] Sutour | Estimation of the noise level function based on a nonparametric detection of homogeneous image regions[END_REF].

Figure 5 displays the nine components of y for an image of 3 × 3 matrices C. While the noise variance in C is clearly signal-dependent, the amplitudes of the fluctuations in the different channels of y appear fairly constant whatever the underlying signal Σ. Nevertheless, the noise is not signalindependent, we can observe that some regions are a smidgen noisier than others. Remark also that the noise variance is also about the same for all channels. Then, we can reasonably claim from this experiment that our log channel decomposition approximately stabilizes the noise variance.

Thanks to this approximate stabilization, we notice that choosing 6 iterations with β = 1 + 2/L provides satisfying solutions in all tested situations, whatever the dimension D, the number of looks L and the embedded Gaussian denoiser. The initialization is chosen as x = y, ẑi = f 1 (y i ) and d = ẑx. We always choose Q = 1 according to Sec. III-E.

When L < D, the matrix C k is rank deficient, and the likelihood term [START_REF] Knaus | Dual-domain image denoising[END_REF], though convex when restricted to commutative matrices, is however no longer strictly convex and has an infinite number of minimizers. Worse, its logarithm Ck is undefined which prevents the computation of A, b and Φ used in our log-channel decomposition. To deal with these issues a practical solution is to perform a small diagonal loading of the input matrices C k to enforce their positive definite property. After trying several other alternatives, we found that satisfying results were reached by performing a spatially varying re-scaling of the off-diagonal elements as

(C k ) new i,j ← | l w k,l (C l ) i,j | l w k,l (C l ) ii l w k,l (C l ) jj (C k ) i,j , (50)
where w k,l are the weights of a Gaussian kernel with bandwidth 1 pixel. This procedure inevitably introduces some bias in the solution, but this bias does not seem to be significant.

IV. EFFICIENT MATRIX LOG AND EXP TRANSFORMS

Our algorithm requires to compute several times n matrix logarithms, matrix exponentials, and matrix-by-matrix products. It is thus important to make these operations as fast as

and the relation between both are given, for C = log C, as

   ã = [(a -b + δ)ℓ 1 -(a -b -δ)ℓ 2 ] /(2δ) , b = [(b -a + δ)ℓ 1 -(b -a -δ)ℓ 2 ] /(2δ) , c = c(ℓ 1 -ℓ 2 )/δ , (52)
and, for C = exp C, as

       a = (ã -b + δ)e 1 -(ã -b -δ)e 2 /(2 δ) , b = -(ã -b -δ)e 1 -(ã -b + δ)e 2 /(2 δ) , c = c(e 1 -e 2 )/ δ , (53)
where

δ = 4|c| 2 + (a -b) 2 , ℓ 1 = log (a + b + δ)/2 , ℓ 2 = log (a + b -δ)/2 , δ = 4|c| 2 + (ã -b) 2
, e 1 = exp (ã + b + δ)/2 and e 2 = exp (ã + bδ)/2 . With these established relations, we can compute the matrix logarithm and matrix exponential in a vectorial way, starting by computing for the n pixels the coefficient ã, and next b and finally c. We have implemented this vectorial procedure in Matlab (using the element-wise operators . * and ./). and observed an acceleration by a factor larger than 600 compared to calling n times the dedicated Matlab's function. We were also able to derive closed form formula for the case where D = 3, leading in this case to an acceleration by a factor 75. The formulas being too long to be inserted in this paper, we invite the interested reader to refer to [START_REF] Deledalle | Closed-form expressions of the eigen decomposition of 2 x 2 and 3 x 3 Hermitian matrices[END_REF]. The same workaround is applied for the matrix-by-matrix product.

V. NUMERICAL EXPERIMENTS a) Considered Gaussian denoisers: In this section, numerical experiments are conducted based on three Gaussian denoisers. First, we consider the classical isotropic totalvariation (TV) [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF], i.e., R(x i ) = λ n k=1 ||(∇x i ) k || which corresponds to the multivariate extension of [START_REF] Bioucas-Dias | Multiplicative noise removal using variable splitting and constrained optimization[END_REF]. The parameter λ has been tuned to 0.7 and kept the same for all data-sets, whatever the dimension D or the number of looks L. This can be achieved only because noise components in the different channels have all been stabilized to a variance of 1. Next, we consider DDID [START_REF] Knaus | Dual-domain image denoising[END_REF], a recent hybrid method based on bilateral filtering and shrinkage of Fourier coefficients. Finally, we consider BM3D [START_REF] Dabov | Image denoising by sparse 3-D transform-domain collaborative filtering[END_REF], a ten-years algorithm, reaching remarkable results with fast computation, based on patch similarity and three-dimensional wavelet shrinkage. The inner parameters of these two latter algorithms have been kept the same as the one provided in the authors' implementation.

b) Simulations in the univariate case: Figure 6 reports the results of a simulated experiment in the SAR amplitude context (D = 1). Speckle has been generated on top of a 409 × 409 "speckle-free" SAR image 7 of the CNES in Toulouse, France (obtained by ONERA's RAMSES sensor). These images having large dynamic ranges, they have been saturated and gamma-corrected for display purposes. Speckle with levels from L = 1 to 32 have been considered. The performance of MuLoG is compared to the homomorphic approach when both embed the aforementioned Gaussian de- noisers. As a baseline, the performance are also compared to a recent patch-based speckle reduction technique: NL-SAR [START_REF] Deledalle | NL-SAR: A unified nonlocal framework for resolution-preserving (Pol)(In) SAR denoising[END_REF].

In the single-look case (L = 1), close-ups of size 256 × 80 are given for visual inspection. PSNR is first used as a quantitative measure of amplitude reconstruction quality 8 . In this high dynamic range setting, PSNR values are extremely sensitive to a small discrepancy, all the more for tiny bright structures. As a consequence, it may appear irrelevant, for instance, by ranking MIDAL worst than the noisy image itself 9 . For this reason, we have also included SSIM [START_REF] Wang | Multiscale structural similarity for image quality assessment[END_REF] values. Visual inspection and quantitative criteria seem to indicate that MuLoG provides more relevant solutions than the homomorphic approach. In particular, the homomorphic approach tends to oversmooth bright targets and leaves residual dark stains, revealing its lack of robustness against the heavy left tail of the Fisher-Tippet distribution. Overall, results obtained with MuLoG are on a par with those of NL-SAR, while displaying small variations depending on the chosen Gaussian denoiser. Visual analysis of the output produced by each variant of MuLoG may be useful to discard structures appearing only with one Gaussian denoiser as artifacts. c) Results in the multivariate case: In the sub-figures (b) are displayed the results for the isotropic total-variation (TV) regularization. In this case, ADMM finds a local minimum of [START_REF] Karplus | A note on saturation in microwave spectroscopy[END_REF]. We can observe that the solution in each channel inherits from the well-known behavior of the univariate TV regularization: the solution appears piece-wise constant, small details and low contrasted features are lost, and a slight bias can be observed.

In the sub-figures (c) and (d) are displayed the results for DDID [START_REF] Knaus | Dual-domain image denoising[END_REF] and BM3D [START_REF] Dabov | Image denoising by sparse 3-D transform-domain collaborative filtering[END_REF]. In these cases, our proposed algorithm does not explicitly minimize a given energy but we can observe that it converges in practice to relevant solutions. As for TV, we can observe that the solution in each channel inherits from the behavior of the original method: in this case the small details and low contrasted features are well restored but some known typical artifacts (small oscillating features) of these two methods can be observed as well. Finally, we display in sub-figure 9(e) the result of NL-SAR [START_REF] Deledalle | NL-SAR: A unified nonlocal framework for resolution-preserving (Pol)(In) SAR denoising[END_REF] in the SAR polarimetric context. It provides a baseline for assessing the relative performance of MuLoG against a recent patch-based speckle reduction technique dedicated to SAR imagery. Compared to the different MuLoG versions, NL-SAR provides slightly over-smoothed results but is free of systematical oscillating artifacts inherent of DDID and BM3D. d) Computation time analysis: Table I reports the computation time of our proposed approach with different embedded Gaussian denoisers, different image sizes and different covariance matrix sizes. Our implementation uses parallelization obtained by running eq. ( 46) in parallel on the D 2 channels, and running eq. ( 48) on different subsets of the n pixels. Our experiments were conducted on a processor with 4 cores. In this experiment, we observe that the total computation time is always faster than running successively 6 × D 2 times the original algorithm. This paper only deals with the problem of applying Gaussian filters to multi-variate SAR data. A deeper study comparing their relative performance or comparing them to stateof-the-art multi-variate SAR filters is out of the scope of this paper. However, we believe that this study should be performed in a future work and is of main interest for the community. For this reason, we have released a Matlab script 10 , under CECILL license, that takes as input a multivariate SAR image, its number of looks and a Gaussian denoiser provided by the user, and outputs the filtered SAR image.

VI. CONCLUSION

While a large diversity of speckle reduction methods exist for intensity images, only few can be extended to process multi-channel SAR images. In particular, the extension of variational methods leads to challenging optimization problems due to the nonconvexity of the neg-log-likelihood and the positive definiteness constraint on the covariance matrices. Furthermore, signal and channel dependent variance lead to restoration results with an uneven suppression of speckle. This paper introduced a general scheme based on a matrix logarithm transform to approximately stabilize speckle variance and produce close to independent channels. Each channel can then be processed with a user-defined Gaussian denoiser. Upon re-iterating, a good fit of the restored multi-channel image with the Wishart distribution of input covariance matrices is enforced.

Special care is paid to ensure that the method requires no parameter tuning other than possibly within the Gaussian denoiser, and that these parameters, if any, can be tuned once for all. The resulting generic method can then include Gaussian denoisers selected by the user and tremendously extends the set of available speckle reduction methods applicable to multi-channel SAR images. We believe that this offers several for providing the RAMSES data, the Airbus Defence and Space for the TerraSAR-X data, the German Aerospace Center (DLR) for the F-SAR data. This work has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the GOTMI project (ANR-16-CE33-0010-01) the financial support from DGA (Direction Générale de l'Armement) and ANR in the frame of ALYS project (ANR-15-ASTR-0002).

S. Tabti also thanks the DGA and the Centre National de la Recherche Scientifique for funding her PhD program.

1 3

 1 tr C k . The output Σk are then rescaled back to their original range by multiplication by τ -1

  ) with f ( C) = |J exp ( C)|/ exp[D tr( C)] a normalization factor that involves the Jacobian of log transform |J exp ( C)| and that is equal to 1 when D = 1. We call such a distribution the Wishart-Fisher-Tippett distribution, denoted as WFT ( Σ; L), as it generalizes the Fisher-Tippett distribution to the case where D > 1. The expectation and variance of C do not seem to be known in closed form in the literature. Nevertheless, according to [2] its trace (which coincides with the logarithm determinant of C: tr C = log |C|) has the following statistics tr C = tr Σ + tr S ,

Fig. 3 .

 3 Fig.3. Three two-dimensional sections of the Wishart-Fisher-Tippett negative log likelihood with respect to (x 1 , x 2 ), (x 1 , x 3 ) and (x 3 , x 4 ).

Fig. 4 .

 4 Fig. 4. Evolution with respect to the matrix dimension D of the error ||∆ 100 ||/||x|| after 10 iterations of Newton's method for which the update ∆∞ is replaced by ∆ Q , for Q = {0, 1, 2, 4, 8, 16}. For each D, the results are averaged over 100 realizations of C (with L = D), each realization is obtained from a different version of Σ (also randomly generated).

  Fig. 5. (a) the nine log channels x 1 , . . . , x 9 , with approximate unit stabilized variance, of (c) the image C generated from (b) the image Σ displayed with a RGB representation based on its decomposition in Pauli's basis.

Fig. 6 .

 6 Fig. 6. (a) A "speckle-free" SAR image, (b) its simulated version with single-look speckle (L = 1), Results of speckle reduction with (c) NL-SAR [20],(d-f) MuLoG embedding respectively TV, DDID and BM3D, (g-i) the homomorphic approach embedding respectively TV, DDID and BM3D. (j-l) Measure of performance on this image in terms of PSNR (the larger the better) and SSIM (the larger the better) metrics as a function of the number of looks L.

Fig. 7 .

 7 Fig. 7. (a) A single-look interferometric image (D = 2, L = 1) of Toulouse (France) sensed by RAMSES ( c ONERA). Estimation with our proposed multivariate framework MuLoG for (b) Total-Variation (TV), (c) DDID and (d) BM3D. For each image, the four log channels are displayed on top x 1 , x 2 , x 3 , x 4 , below are the trace tr C, the phase arg C 12 and the coherence |C 12 |/ |C 11 C 22 | respectively.

  Figures 7,[START_REF] Brifman | Turning a denoiser into a superresolver using plug and play priors[END_REF] and 9 give three illustrations of MuLoG in real speckle reduction contexts. Figure 7 corresponds to a 256 × 256 airborne singlelook SAR interferometric image (D = 2, L = 1) of a building in Toulouse, France (sensed by RAMSES).

  Figure 8 corresponds to a 250 × 250 spaceborne single-look SAR interferometric image (D = 2, L = 1) of the dam of Serre-Ponc ¸on, France (sensed by TerraSAR-X).

Figure 9

 9 corresponds to a 512 × 512 airborne SAR polarimetric image (D = 3, L = 1) of the city of Kaufbeuren, Germany (sensed by F-SAR).

Fig. 8 .

 8 Fig. 8. (a) A single-look interferometric image (D = 2, L = 1) of Serre-Ponc ¸on (France) sensed by TerraSar-X (image courtesy to Airbus Defence and Space). Estimation with our proposed multi-variate framework MuLoG for (b) Total-Variation (TV), (c) DDID and (d) BM3D. For each image, the four log channels are displayed on top x 1 , x 2 , x 3 , x 4 , below are the trace tr C, the phase arg C 12 and the coherence |C 12 |/ |C 11 C 22 | respectively.

TABLE I COMPUTATION

 I TIME ON AN "INTEL(R) CORE(TM) I5-4590S CPU @ 3.00GHZ". TIME FOR GD IS THE RUNNING TIME OF THE GAUSSIAN DENOISER ON ONE SCALAR CHANNEL. TIME FOR NM IS FOR THE 10 ITERATIONS OF NEWTON'S METHOD WITH Q = 1. TOTAL TIME IS THE ACTUAL TIME OF OUR IMPLEMENTATION WITH 6 ITERATIONS OF THE ADMM SCHEME WHERE SOME PARTS OF THE CODE ARE RUN IN PARALLEL ON 4 DIFFERENT CORES.

	Image size 256 × 256 256 × 256 512 × 512	D 2 2 × 2 3 × 3 3 × 3	File size 2Mb 4.5Mb 18Mb	Gaussian denoiser BM3D TV DDID BM3D TV DDID BM3D TV DDID	Time for GD 0.87s 1.43s 38.4s 0.87s 1.43s 38.4s 3.6s 6.3s 2min 37s	Time for NM Total time 0.35s 8.97s 17.08s 4min 29s 2.20s 29.21s 48.18s 12min 6s 9.97s 2min 09s 4min 4s 48min 52s

note that these other methods are only applicable to intensity SAR images (i.e., single-channel)

for complex random variables Var[C ij ] = E[|C ij | 2 ] -|E[C ij ]| 2

this requires to add a small update of β during the iterations, according to the rule given in[START_REF] Chan | Plug-and-play ADMM for image restoration: Fixed point convergence and applications[END_REF].

convergence holds in our case by continuity of the gradient of[START_REF] Karplus | A note on saturation in microwave spectroscopy[END_REF], which implies that the gradient is Lipschitz in any compact subsets of R D 2 ×n .

Speckle in this 1 meter ground-resolution image has been reduced by a factor 100 by averaging a 10 centimeters ground-resolution image.

Since images are not 8bit, the standard PSNR has been modified with a peak value defined as the

99th quantile of the "speckle-free" SAR image.[START_REF] Chan | Plug-and-play ADMM for image restoration: Fixed point convergence and applications[END_REF] TV regularization produces a systematic loss of contrast (see,[START_REF] Strong | Edge-preserving and scale-dependent properties of total variation regularization[END_REF],[START_REF] Deledalle | CLEAR: Covariant LEAst-square Re-fitting with applications to image restoration[END_REF]), even stronger for punctual bright structures.

http://www.math.u-bordeaux.fr/ ∼ cdeledal/mulog notable advantages: (i) the SAR imaging community will

directly benefit from upcoming progress made in the field of image denoising; (ii) several images with reduced speckle can be produced by very different denoising algorithms, and these images can be compared to discard artifacts and confirm weak structures; (iii) this family of speckle reduction methods can serve as a reference to benchmark future specialized speckle reduction algorithms. This motivated the release of an opensource code implementing our method 10 .

APPENDIX A GRADIENT OF THE NEG LOG LIKELIHOOD

In this section, we establish the part of formula [START_REF] Lewis | The convex analysis of unitarily invariant matrix functions[END_REF] corresponding to the gradient of tr Ω(x) + e Ω(y) e -Ω(x) . To this end, recall that for a real valued differentiable function f ,

where g is a matrix valued function whose dimension depends on that of x, and * denotes the adjoint operator. Since Ω(x) = K(AΦx + b) is affine, this directly implies that

We now use that for any matrix valued function h, we have (proven in [START_REF] Karplus | A note on saturation in microwave spectroscopy[END_REF], [START_REF] Bellman | Introduction to matrix analysis[END_REF] according to [START_REF] Najfeld | Derivatives of the matrix exponential and their computation[END_REF])

It follows that, for any matrix L, we have

= tr L 1 0 e uh(x) (dh(x))e (1-u)h(x) du ,

= tr

Choosing h(x) = -Ω(x) and L = e Ω(y) leads to

the integral term being Hermitian since it reads as the integral of Hermitian matrices 1 2 0 e (u-1)Ω(x) e Ω(y) e -uΩ(x) + e -uΩ(x) e Ω(y) e (u-1)Ω(x) Hermitian du .

ACKNOWLEDGMENT

The authors would like to thank the anonymous reviewers as well as the associate editor for their comments, criticisms and encouragements. We also thank the Centre National d' Études Spatiales, the Office Nationale d' Études et de Recherches Aérospatiales and the Délégation Générale pour l'Armement


domain: info math

text:
Miguel Couceiro
email: miguel.couceiro@loria.fr

Pierre Mercuriali
email: pierre.mercuriali@loria.fr

Romain Péchoux
email: romain.pechoux@loria.fr

Abdallah Saffidine
email: abdallah.saffidine@gmail.com

Median based calculus for lattice polynomials and monotone Boolean functions

In this document, we consider a median-based calculus for efficiently representing polynomial functions over distributive lattices. We extend an equational specification of median forms from the domain of Boolean functions to the domain of lattice polynomials. We show that it is sound and complete, and we illustrate its usefulness when simplifying median formulas algebraically. Furthermore, we propose a definition of median normal forms (MNF), that are thought of as minimal median formulas with respect to a structural ordering of expressions. We also investigate related complexity issues and show that the problem of deciding whether a formula is in MNF is in Σ P 2 . Moreover, we explore polynomial approximations of solutions to this problem through a sound term rewriting system extracted from the proposed equational specification.

I. INTRODUCTION

Representing a function with various operators sometimes yields drastically different results in terms of the size of the formula used to represent it. In the case of Boolean functions, it has been shown that a median-based representation, that is, a representation based on the ternary median operator, yields asymptotically smaller formulas -in terms of number of connectives-than the classical DNF, CNF, and polynomial normal forms [START_REF] Couceiro | Composition of post classes and normal forms of boolean functions[END_REF]. Moreover, algorithmic procedures to obtain such median representations were given in [START_REF] Couceiro | An algorithm for producing median normal form representations for Boolean functions[END_REF]. However, these procedures may not produce median formulas of the lowest possible complexity (size), and consequently procedures for simplifying median formulas are required, in analogy with well known resolution procedures for DNF and CNF expressions.

Since median expressions can be translated into the language of lattices and conversely, lattice polynomials can also be represented by median expressions. Thus, in addition to applications in logic and circuit design, these simplifications become useful when efficiently representing noteworthy classes of aggregation functions such as the Sugeno integral ( [START_REF] Grabisch | Aggregation functions[END_REF]).

Motivated by these observations and following the work of [START_REF] Couceiro | Composition of post classes and normal forms of boolean functions[END_REF], in this paper we investigate the use of the ternary median connective to efficiently represent lattice polynomials, which in the case of 2-element lattices reduce to monotone Boolean functions. In Section II we recall basic background on lattice functions and median algebras that will be used throughout this paper. In Section III we give a finite equational specification (System 1) that allows the simplification of median formulas while preserving logical equivalence, and we show that this equational specification is both sound and complete (Theorem 1). As an immediate consequence, it then follows that we may rewrite any formula into any other that is logically equivalent to it, and that median formulas can be simplified algebraically according to this equational specification. We also propose a structural and lexicographic ordering of formulas in Section IV, and introduce median normal forms (MNFs) as being median formulas that are minimal with respect to this ordering.

In Section V we consider two decision problems related to the task of finding an MNF of a given formula, and we show that both are at most moderately intractable, by which we mean that it is likely to be intractable, but not beyond Σ P 2 . We also explore a different approach to the potential intractability issue by providing a term rewriting system based on the equational specification (System 1). Even though the resulting system is not complete, it runs in polynomial time, thus highlighting a trade-off in complexity: either we sacrifice completeness but preserve tractability, or we keep completeness at the cost of high complexity.

II. LATTICE POLYNOMIALS AND MEDIAN FORMULAS

A. Notation and Preliminaries

In this subsection we recall definitions and notations on lattice and lattice functions, while adopting the terminology of [START_REF] Couceiro | Characterizations of discrete sugeno integrals as polynomial functions over distributive lattices[END_REF]. A lattice is an algebraic structure L, ∧, ∨ where L is a nonempty set, called universe, and where ∧ and ∨ are two binary operations that satisfy the laws of commutativity, associativity, absorption, and idempotence; a lattice is said to be distributive if the two laws distribute over one another. With no danger of ambiguity, we will denote a lattice L, ∧, ∨ by its universe L.

In what follows, L will always denote an arbitrary bounded distributive lattice with least and greatest elements ⊥ and , respectively. For a, b ∈ L, a ≤ b means that a ∧ b = a or, equivalently, a ∨ b = b. For any integer n ≥ 1, we set [n] = {1, . . . , n}. For an arbitrary nonempty set A and a lattice L, the set L A of all functions from A to L also constitutes a lattice under the operations

(f ∧ g)(x) = f (x) ∧ g(x) and (f ∨ g)(x) = f (x) ∨ g(x)
for every f, g ∈ L A . In particular, any lattice L induces a lattice structure on the Cartesian product L n , n ≥ 1, by defining ∧ and ∨ componentwise, i.e.,

(a 1 , . . . , a n ) ∧ (b 1 , . . . , b n ) = (a 1 ∧ b 1 , . . . , a n ∧ b n ), (a 1 , . . . , a n ) ∨ (b 1 , . . . , b n ) = (a 1 ∨ b 1 , . . . , a n ∨ b n ).
We denote the elements of L by lower case letters a, b, c, . . . and the elements of L n , n > 1, by bold face letters a, b, c, . . . . We now recall the notion of lattice polynomial function. The class of lattice polynomial functions (or simply, polynomial functions) from L n to L is defined inductively as the set of functions represented by expressions constructed in the language of lattices: the projections x → x i , the constant functions x → c, c ∈ L, and if f and g are polynomial functions, then so are f ∧ g and f ∨ g.

For instance, the ternary median m, i.e., the function given by:

m(x 1 , x 2 , x 3 ) = (x 1 ∧ x 2 ) ∨ (x 2 ∧ x 3 ) ∨ (x 3 ∧ x 1 ) = (x 1 ∨ x 2 ) ∧ (x 2 ∨ x 3 ) ∧ (x 3 ∨ x 1 )
is an example of such a polynomial function.

B. Equational specification

We make use of the notations of [START_REF] Klop | Term rewriting systems[END_REF] and [START_REF] Baader | Term rewriting and all that[END_REF] to introduce equational specifications. An equational specification, (sometimes called equational system within this document) is a pair (Σ, E) of an alphabet or signature Σ and a set of equations E. The alphabet Σ consists of a countably infinite set of variables x 1 , x 2 , . . . and a nonempty set of function symbols or operator symbols. In the current setting, this set contains m and constants.

The set of terms (or expressions) over Σ is denoted by Ter(Σ) and it is defined inductively as follows:

(1) every variable and constant in Σ is in Ter(Σ), and (2) if f is an n-ary function symbol and t 1 , . . . , t n are terms, then f (t 1 , . . . , t n ) is in Ter(Σ). An equation is then an expression of the form s = t where s, t ∈ Ter(Σ). The set of all equations is denoted by E.

The set of median terms (also referred to as median expressions or median formulas in [START_REF] Couceiro | Composition of post classes and normal forms of boolean functions[END_REF]) M will denote the set of all formulas that are constructed using variables, constants and the median m. Even though [START_REF] Couceiro | Composition of post classes and normal forms of boolean functions[END_REF] focuses on median normal expressions for Boolean functions, this notion adapts rather naturally to lattice polynomial functions.

Given a formula φ in M, its depth is denoted by d(φ) and defined as follows:

(i) for every variable or constant a, d(a) = 0, (ii) for every formula φ = m(a, b, c) ∈ M,

d(φ) = max{d(a), d(b), d(c)} + 1.
The size |φ| of a median term φ is the number of medians in it.

Example 1. Let φ = m(m(x, x, y), m(x, y, z), v). Then d(φ) = 2 and |φ|= 3.

Two median terms φ and ψ are said to be equivalent, denoted by φ ≡ ψ, if they represent the same function.

Example 2. For instance, the median terms:

φ 1 = m(m(m(x, u, v), m(x, y, v), m(y, u, v)), x, v) φ 2 = m(m(u, y, v), x, v)
are equivalent. However, the median terms:

φ 3 = m(m(x, y, z), u, v) and φ 4 = m(x, m(y, z, u), v)
are not equivalent. In other words, m is not associative in the sense of [START_REF] Aczél | The associativity equation re-revisited[END_REF], [START_REF] Couceiro | Aczélian n-ary semigroups[END_REF], [START_REF] Post | Polyadic groups[END_REF].

A substitution is a map σ from Ter(Σ) to Ter(Σ) that satisfies

σ(F (t 1 , . . . , t n )) = F (σ(t 1 ), . . . , σ(t n ))
for every n-ary function symbol (here n ≥ 0).

Substitution together with the other rules recalled in Table I, give rise to the so-called derivable equations, i.e., equations obtained by applying a finite combination of these rules. If an equation s = t is derivable from the equations in E, than we write (Σ, E) s = t or s E t.

(Σ, E) t = t if t ∈ Ter(Σ) (Σ, E) s = t if s = t ∈ E (Σ, E) s = t (Σ, E) t = s (Σ, E) t 1 = t 2 , (Σ, E) t 2 = t 3 (Σ, E) t 1 = t 3 (Σ, E) s = t (Σ, E) σ(s) = σ(t)
for every substitution σ

(Σ, E) s 1 = t 1 , . . . , (Σ, E) sn = tn (Σ, E) F (s 1 , . . . , sn) = F (t 1 , . . . , tn) for every n-ary F ∈ Σ Table I EQUATIONAL INFERENCE SYSTEM.

C. Term Rewriting systems

A Term Rewriting System (TRS) is an equational specification with all its equations oriented. A pair (l, r) of terms in Ter(Σ), written as l -→ r, is a reduction rule if r is not a variable and all the variables in l are already contained in r.

Each term rewriting system yields a rewrite relation defined to be the closure by substitution and context of its reduction rules.

III. AXIOMATIZATION OF LATTICE POLYNOMIALS

In this section, we present an equational system for median calculus which is both sound and complete, and which we then use to manipulate median expressions.

First, recall that lattice polynomial functions f : L n → L are exactly the solutions of the median decomposition system [START_REF] Marichal | Weighted lattice polynomials[END_REF]:

f (x) = m(f (x ⊥ k ), x k , f (x k )) (1)
for all x = (x 1 , . . . , x n ), k ∈ [n] and c ∈ L, and where

x c k := (x 1 , . . . , x k-1 , c, x k+1 , . . . , x n ).
A direct consequence of this result is that every polynomial function can be represented by a median term. Indeed, recursive applications of (1) on each x 1 , . . . , x n of an n-ary polynomial function f produces a median formula representing f , see, e.g., [START_REF] Couceiro | An algorithm for producing median normal form representations for Boolean functions[END_REF]. However, this procedure may not produce optimal formulas (size wise), and this fact motivates the current study.

Example 3. Consider the 5-ary majority operator m 5 :

x → m 5 (x 1 , x 2 , x 3 , x 4 , x 5 ) over Boolean variables. Using the median decomposition algorithm mentioned above we can construct a median formula representation of this function using its values on every point of {⊥, } 5 . This representation, of size 31, is not minimal. There are smaller representations of size 4, as shown in [START_REF] Bandelt | The algebra of majority consensus[END_REF] and [START_REF] Miller | Majority-logic synthesis by geometric methods[END_REF].

Let us now recall an axiomatisation of the algebraic structure L, m, ⊥, , that is the set L with the ternary median function m and the 0-ary functions (constants) ⊥ and , by the following equational system.

System 1.

(M 1) m(x, y, z) = m(x, z, y) = m(z, x, y),

(M 2) m(x, x, y) = x, (M 3) m(m(x, u, v), m(y, u, v), z) = m(m(x, y, z), u, v), (M 4) m(⊥, , x) = x,
for all x, y, z, u, v in L.

Equations (M 1), (M 2), (M 3) are known as symmetry, majority, distributivity laws, respectively. Remark 1. An axiomatisation of the Boolean algebra {0, 1}, ∧, ∨,., 0, 1 was used in [START_REF] Amarú | Majority-inverter graph: A novel data-structure and algorithms for efficient logic optimization[END_REF] to simplify Boolean circuits. System 1 is an adaptation of these results to the case of lattice polynomials.

In order to manipulate median formulas, we need to know whether the axiomatisation given by System 1 is both sound and complete. Soundness means that every equation s = t that can be derived from System 1 is valid, i.e., that the formulas s and t are equivalent. Completeness means that every equation s = t that is valid can be derived from the axioms of the system. For our purposes, having a sound and complete system is interesting in order to rewrite median terms, and hopefully simplify them. Soundness indeed ensures that whatever simplification we do by applying an equation to a median term will preserve logical equivalence, while completeness ensures that we can infer a median term from an equivalent one by using System 1.

Theorem 1. The algebra L, m, ⊥, together with the axioms of System 1 is sound and complete. Sketch of Proof. Soundness is provable algebraically by deriving the axioms of the system using the properties of the lattice L. Completeness is an immediate consequence of the Birkhoff's Completeness Theorem for equational logic, see, e.g., [START_REF] Burris | A Course in universal algebra[END_REF], [START_REF] Klop | Term rewriting systems[END_REF].

IV. MEDIAN NORMAL FORMS

In this section, we propose a structural description of median formulas and introduce the notion of median normal forms (MNF) that, as we shall see, correspond to median formulas that are "minimal" with respect to the lexicographical ordering of their structural description.

A. Orders

A binary relation on a set S is a quasi-order, or preorder, if it is reflexive and transitive. A quasi-order is said to be a partial order if it is antisymmetric. If is a partial order, the structure S, is called a partially ordered set (or poset). A quasi-ordered set is well-founded if it satisfies the descending chain condition, i.e., there exists no infinite decreasing sequence • • • < t 2 < t 1 of elements of T . If for every pair (a, b), either a ≤ b or b ≤ a, then (T , ≤) is said to be totally ordered. A well-founded and totally ordered set is said to be a well-ordered set.

Let n be a positive integer, and T n the set of all ordered n-tuples over S. Let T = n≤1 T n . The lexicographical extension on T denoted by lex , is defined by: (x 1 , . . . , x m ) lex (y 1 , . . . , y n ) if

• m ≤ n and for all k ∈ [m], x k = y k 1 , or • there is k ∈ [min(m, n)] such that for allj ∈ [k -1],
we have x j = y j and x k ≺ y k (i.e., x k y k , x k = y k )

Example 4. The lexicographic order defined on a finite set of words is well-founded. On the other hand, for S = {a, b} with a b, the lexicographical extension on the infinite product S * is not a well-order:

. . . aaab aab ab b.

We now define a way to partially describe the structure of median terms that induces a well-order on the set of median formulas. We then show that for every polynomial function, there exists a set of minimal representations with respect to this order, which we call median normal forms (MNF). As it is the case for DNF and CNF representations, this representation is unique modulo some properties like commutativity or associativity. However, the general structure of these minimal representations still eludes us. Definition 1. Let φ be a median term of depth d. Let n 0 , . . . , n d be nonnegative integers, such that for all i ∈ {0} ∪ [d], n i is the number of medians at depth ≤ i. The structural representation of φ is the tuple S φ = (n d , . . . , n 0 ). Let ≤ S be the ordering of median formulas defined by:

φ 1 ≤ S φ 2 if S φ1 ≤ lex S φ2 .
Remark 2. Note that S φ is a decreasing sequence and that n d = |φ|. Also, the order ≤ S prioritizes the size of the formula over its depth. For instance, consider the following equivalent formulas

φ 1 = m(x 1 , x 2 , m(x 3 , x 4 , m(x 5 , x 6 , x 7 ))) φ 2 = m(m(x 1 , x 2 , x 3 ), m(x 1 , x 2 , x 4 ), m(x 5 , x 6 , x 7 )). Clearly, |φ 1 |= 3 < 4 = |φ 2 | while d(φ 1 ) = 3 > 2 = d(φ 1 ).
Looking at their structural representation, we have S φ1 = (3, 2, 1) whereas S φ1 = (4, 3), and hence φ 1 ≤ S φ 2 .

We now give a definition of a median normal form as a minimal median representation. Definition 2. We say that a median term φ is a median normal form (MNF) if for every median term φ ≡ φ, we have φ ≤ S φ .

Example 5. The formula φ = m(m(x, x, y), y, z) is not a median normal form since φ = m(x, y, z) is an equivalent formula, and S φ = (1) ≤ lex (2, 1) = S φ .

Remark 3. As it has been defined, the structural order cannot account for permutations of variables. For instance, the formulas m(x, y, z) and m(x, z, y) have the same structural tuple (1), but they are also equivalent, and both are median normal forms. Thus, a formula does not have a single median normal form, but rather a set of median normal forms: the formula φ from Example 5 has for set of normal forms {m(x, y, z), m(x, z, y), m(y, z, x)}.

V. CONSIDERATIONS ON COMPLEXITY

In this section we address the question of finding median normal forms. To this end, we formalize two decision problems that express the tasks of finding median formulas of smaller structural representation. We show that both problems are at most moderately intractable, and we propose a term rewriting system as a tool for approximating solutions to them.

A. Structurally smaller formulas

Although, the definition of the median normal form is expressed simply, a procedure to convert an input formula into an equivalent formula in median normal form is likely intractable. Indeed, we will show that the mere task of checking if a given formula is in MNF seems to be expensive. We formalize this decision problem in Definition 3.

Definition 3. Consider the decision problem SMALLMED:

Input: a median term φ and a decreasing sequence S Output: succeeds if there exists a formula ψ whose structural representation is strictly smaller than S. Fails if none exists.

Before proving that SMALLMED is intractable, we give a Σ P 2 upper-bound on its hardness. Recall that the Σ P 2 complexity class is on the second level of the polynomial hierarchy, between NP and PSPACE [START_REF] Savage | Models of computation: Exploring the power of computing[END_REF].

Theorem 2. SMALLMED is in the class Σ P 2 . Proof. A convenient characterization of Σ P 2 is that it contains decision problems such that the accepting instances can be expressed as a set of words {x : ∃c 1 ∀c 2 F (x, c 1 , c 2 )}, where c 1 and c 2 are certificates whose lengths are polynomial in |x| and F is computable in polynomial time. Consider Algorithm 1 which solves SMALLMED. The size of the first certificate ψ, is indeed polynomial in the size of the input, |φ|, because its structural representation is bounded by that of φ. The size of the second certificate σ is also polynomial in the input. Therefore, Algorithm 1 ensures SMALLMED is in Σ P 2 .

Algorithm 1 Finding a smaller equivalent median form.

Input: A formula φ, a decreasing sequence S.

1: Existentially guess a formula ψ such that S(ψ) < S S 2:

Universally guess an assignment σ 3:

Ensure that σ(φ) = σ(

If so, return SUCCESS 5: If none exist, then FAIL Note that Theorem 2 simply provides an asymptotic upper bound on the complexity of SMALLMED, but a corresponding lower-bound still eludes us.

Definition 3 assumed that the desired formula size was given as part of the input. If instead we assume a constant target size, then it is possible to obtain a better complexity bound. Definition 4. For any fixed decreasing sequence S, we define the decision problem SMALLMED S as: Input: a median term φ Output: succeeds if there is a formula ψ whose structural representation is smaller than S. Fails if none exists. Theorem 3. For any decreasing sequence S, SMALLMED S is in the class co-NP.

Proof. Let s be the first element of the sequence S, n be the number of variables occurring in φ, and V φ be the set of variables occurring in φ. Any formula ψ of structural representation smaller than S, has no more than s medians. Hence, such a formula ψ cannot involve more than N = 2s+1 variables. If φ is equivalent to a formula of structural representation smaller than S, then at most N variables among the ones that occur in φ are relevant.

For every subset of variables V ⊂ V φ of size at most N , there is a constant number of formulas smaller than S with variables drawn from V . For each such formula ψ, universally guess a variable assignment σ for the variables occurring in φ and check if φ and ψ agree on σ. If so, ψ is equivalent to φ and we can succeed. If no formula triggers a success, we can end the algorithm and fail.

The number of ways to specify the set V is bounded by n N . So the total number of universal guesses is bounded by O(n N ). Since N is constant, we conclude that we can determine if φ admits an equivalent formula of size smaller than S with a polynomial number of universal guesses. Therefore, SMALLMED S is in co-NP.

B. Determining the MNF of a formula

Algorithm 1 does not find an MNF for the input formula directly, but it can be used as a subroutine to an algorithm that does. Let φ be an input formula and let S φ be its structural representation, then a kind of binary search allows to identify the smallest structural representation such that Algorithm 1 succeeds. The output of the final call to Algorithm 1 is then an MNF of φ.

A binary search performs a number of comparison calls logarithmic in the size of the ordered domain. In our case, the domain is the set of sequences lexicographically smaller than S φ , and its cardinality is at most exponential in the size of φ. Therefore, the binary search performs a number of calls to Algorithm 1 that is polynomial in the input formula.

C. A term rewriting approach

A naive implementation of Algorithm 1 amounts to an exhaustive search for equivalent formulas among structurally smaller ones. A possible alternative, would be to search for a structurally smaller formula among equivalent ones.

Recall from Theorem 1 that System 1 is sound and complete. In other words, for any two formulas φ and ψ, the formulas are equivalent if and only if there exists a rewriting of φ into ψ by means of a sequence of equations from System 1. This idea is implemented in Algorithm 2. Update ψ ← ψ 5: return ψ This non-deterministic algorithm uses very little space, namely its space complexity is linear in the size of the input, and solves SMALLMED. It therefore constitutes a proof that SMALLMED is in the class NPSPACE. From Savitch's theorem (see, e.g., [START_REF] Savage | Models of computation: Exploring the power of computing[END_REF]) we know that NPSPACE = PSPACE, but Σ P 2 is contained in PSPACE so Theorem 2 is stronger.

We do not know whether Algorithm 2 is guaranteed to terminate in polynomial time, in the best case, when a structurally smaller formula exists. Indeed, equational reasoning to transform a formula into some smaller one may conceivably require an exponential number of rewriting steps. On the contrary, were Algorithm 2 to always terminate in polynomial time, it would constitute a proof of SMALLMED being in NP.

A solution to the problem of finding a derivation without unbounded searches is to orient the equations of the system using the order on the size of the formula: from bigger to smaller. As a result, applying any rule (except commutativity) to a formula will simplify it with regard to ≤ S : every derivation will be a simplification.

Let us then consider the following "term rewriting system" extracted from System 1 by orienting its equations according to the decreasing structural ordering.

System 2.

(R1) m(x, y, z) = m(x, z, y) = m(z, x, y),

(R2) m(x, x, y) -→ x, (R3) m(m(x, u, v), m(y, u, v), z) -→ m(m(x, y, z), u, v), (R4) m(⊥, , x) -→ x,
for all variables x, y, z, u, v in L.

The rewrite rule (R1), which is in fact (M 1), is kept the same, without orientation. Such systems are sometimes called rewriting systems modulo commutativity ([3]). A similar situation takes place when dealing with systems that involve commutative binary operations like ∨ or ∧, and in such cases commutativity and associativity may be kept as equational rules. As explained in [START_REF] Veroff | First-order proof of a median algebra problem[END_REF], in the case of median terms, commutativity (M 1) can be oriented by defining a total order on terms. Some effects of this orientation are the sorting of the terms by the application of the now-oriented commutativity rules, as well as the occasional blockage of derivation proofs (in these cases we thus lose completeness of the system).

By orientating the rules, we can more easily test a derivation between formulas. Lemma 1. Let φ and ψ be median formulas. If it exists, the derivation between φ and ψ is polynomial in the size of φ. Sketch of Proof. Every rewrite rule from System 2 save from the permutation (R1) removes a median m from any formula it is applied to. If such a derivation exists, necessarily φ has more medians than ψ. The derivation between φ and ψ thus contains n steps, with n being the difference between the number of medians in φ and the number of medians in ψ.

However, we may not be able to rewrite a formula into another equivalent one using System 2 (e.g., it is not possible to rewrite x into m(⊥, , x)), much less into a normal form. Example 6. Consider φ = m(m(m(x, y, z), u, v), y, z), which has for canonical form φ = m(m(x, u, v), y, z). It is not possible to simplify φ into φ using System 2, because it is not possible to apply any rule from System 2 to φ. Now, even though System 2 is no longer complete, the following result shows that it remains sound. Proposition 1. System 2 is sound but not complete. Sketch of Proof. Soundness follows from the fact that System 1 is sound (Theorem 1). Incompleteness follows from Example 6.

VI. CONCLUSION AND FUTURE WORK

In this paper, we discussed median-based formalism to efficiently represent monotone Boolean functions as well as polynomial functions over distributive lattices. This was achieved by proposing so-called median normal forms that were defined as median expressions that are minimal with respect to a structural ordering of formulas.

We also formalized the task of finding median formulas of smaller structural representation and investigated its computational complexity. This task turns out to be at most moderately intractable. In fact, we showed that the corresponding decision problem falls into Σ P 2 or co-NP according to whether the structural representation is given as part of the input.

However, the question of determining corresponding complexity lower bounds remains open. This and other complexity questions concerning decision problems that appear naturally in this median-based formalism are to be investigated in forthcoming collaborations.

Algorithm 2 3 :

 23 Bringing a formula closer to a median normal form. Input: A formula φ, a decreasing sequence S. Output: A formula ψ with structural representation S ψ < S, diverges if none exists. 1: Set ψ ← φ 2: While ψ ≥ S φ do Existentially guess a rewriting from System 1: ψ = ψ 4:

I.e., (x 1 , . . . , xm) is a prefix of (y 1 , . . . , yn).

* was supported by the Australian Research Council (project DE 150101351).