Statistical Component Separation for Targeted Signal Recovery in Noisy Mixtures
Authors: Bruno Régaldo-Saint Blancard
Michael Eickenberg
Venue: Accepted by TMLR ()

Abstract: Separating signals from an additive mixture may be an unnecessarily hard problem when one is only interested in specific properties of a given signal. In this work, we tackle simpler "statistical component separation" problems that focus on recovering a predefined set of statistical descriptors of a target signal from a noisy mixture. Assuming access to samples of the noise process, we investigate a method devised to match the statistics of the solution candidate corrupted by noise samples with those of the observed mixture. We first analyze the behavior of this method using simple examples with analytically tractable calculations. Then, we apply it in an image denoising context employing 1) wavelet-based descriptors, 2) ConvNet-based descriptors on astrophysics and ImageNet data. In the case of 1), we show that our method better recovers the descriptors of the target data than a standard denoising method in most situations. Additionally, despite not constructed for this purpose, it performs surprisingly well in terms of peak signal-to-noise ratio on full signal reconstruction. In comparison, representation 2) appears less suitable for image denoising. Finally, we extend this method by introducing a diffusive stepwise algorithm which gives a new perspective to the initial method and leads to promising results for image denoising under specific circumstances.

Published in Transactions on Machine Learning Research (02/2024)
Statistical Component Separation for Targeted Signal Recovery in Noisy Mixtures
Bruno Régaldo-Saint Blancard
bregoloofffcthonnnssttttt.oy
Center for Computational Mathematics
Flatiron Institute
New York, NY 10010
Michael Eickenberg
meichenbers@faatiiniiiitttttt.
Center for Computational Mothematics
Flatiron Institule
New York, NY 10010
Reviewed on OpenReview: https: // openreview. formm? id= y009604
Abstract
Separating signals from an additive mixture may be an unnecessarily hard problem when one is only interested in specific properties of a given signal. In this work, we tackle simpler "statistical component separation" problems that focus on recovering a predefined set of statistical descriptors of a target signal from noisy mixture. Assuming access to samples of the noise process, we investigate method devised to match the statisties of the solution candidate corrupted by noise samples with those of the observed mixture. We first analyze the behavior of this method using simple examples with analytically tractable calculations. Then, we apply it in an image denoising context employing 1) wavelet-based descriptors, 2) ConvNet-based descriptors on astrophysics and ImageNett data. In the case of 1), we show that our method better recovers the descriptors the target data than a standard denoising method in most situations. Additionally, despite not constructed for this purpose, it performs surprisingly well in terms of peak dignal-to-niis ratio on full signal reconstruction. In comparison, reprosentation 2) appears less suitable for image denoising. Finally, extend perspective
this method by introducing a diffusive stepwise algorithm which gives to the initial method and leads to promising results for image denoising under specific circumstances.
1 Introduction
We investigate the properties a new class of source separation algorithms known as statistical component separation methods that has recently emerged for the analysis of astrophysical data (Regaldo-Saing Blancard et al., 2021; Delouis al., 2022; Siahkoohi et al., 2023; Auclair al., 2024). Contrary to standard source separation algorithms, such as blind source separation techniques (Cardoso, 1998), these methods do not focus recovering the signal of interest, but on solely recovering certain statistics features derived from this signal. These have proven successful in scparating signals of distinct statistical natures in varicty of astrophysical contexts, such as the separation of interstellar dust emission and instrumental noise in data from the Planck satellite (Regaldo-Sain, Blancard et al., 2021; Delouis et al., 2022), the separation of interstellar dust emission from the cosmic infrared background (Auclair et al., 2024), or the removal of glitches in seismic data from the InSight Mars mission (Siahkoohi al., 2023). The methodology is not specific astrophysics, and could be of interest to other scientific fields.

---

Published in Transactions on Machine Learning Research (02/2024)
Problem. We observe noisy mixture y zo + co with zu the signal of interest and cp noise process. Signals can be viewed as vectors of TR... We refer to éo as the noise, but we make no assumption on its distribution p(Eo), so that it can include any form of contaminant of the signal z0. However, we assume that we have a way to sample p(co). Let ø be a function, called the representation, that maps 2 to a vector of Teatures or summary statistics in R& with lypically K < M. The ultimate goal of a statistical component. separation method is to recover 6(xo)) Regaldo-Saint Blancard et al. (2021) introduced a first algorithm to do that which consists in constructing to such that:
To E arg min [(x) with ((x) - Eemplea) [[lo(t +c) d(v))
(1)
and where the optimization of L is initialized with y. For suited o, previous works have demonstrated empirically that ((f)) can be a relevant estimate of Φ(zo), and that, while not expected, åp also seemed to
be a reasonable estimate of wo (see Regaldo-Saint Blancard et al. (2021); Delouis et al. (2022); Siahkoohi et al. (2023); Auclair et al. (2024)).² The goal of this paper is to give first formal elements to explain these results, establish performance basclines, and introduce new mothods for solving the problem. For numerical experiments, we will approximate the expected value involved in L by Monte Carlo estimates. Introducing Q independent noise samples 61,-..,60 .Q p(Eo); the corresponding empirical loss L reads:
o l0(x + a)) -
Σ [xx) = Q
Σ
(2)
Related work. Aside from the literature mentioned above, we are not aware of directly related work on this precise problem. However, we mention that adjacent problems of task-adapted reconstructions were explored in learning contexis. In particular, Mairal et al. (2012) investigated lask-adapted dictionary learning, for which sparse data representations can be tuned to specific tasks. More recently, Adler et al. (2022) established a framework for task-related solving of inverse problems and showed how deep neural networks can be used for it. Finally, we add that the approach investigated in this paper shares similarities with sparse regularization techniques proposed in the extensive denoising and component separation literature (e.g., Starck et al., 2005; Selesnick, 2012; Elad et al., 2023). The loss L can indeed be interpreted as a regularization term in a denoising or component separation context. However, we emphasize that the goals differ, since statistical component separation methods primarily focus on recovering @(zo) (and not on recovering zo).
Outline. In Sect. 2, we compute the analytical expressions of the global minimizers of L for different examples of representations ø and in the case of Gaussian noise cp- This will give us a sense of the constraints that & must respect for o(āp) to be a relevant estimate of Φ(zp). In Sect. 3, we describe a first algorithm to perform numerical experiments in typical image denoising set tings for two different representations: the first one based on wavelet phase harmonic statistical descriptors, and the second one based on ConvNet feature maps. Then, in Soct. 4, we discuss strategics to improve the results in the case of Gaussian noise using a new diffusive stepwise algorithm. We finally summarize our conclusions and perspectives in Sect. 5. Codes and data are provided on GitHub.?
Notations. We refer to the components of a vector x E CM as x4 or xfj.. The dot product between x and y is x x E JiWli, and the corresponding norm of x is || = (EM 1 The convolution of at y E CM is x y. We denote the matrix-vector product between A and æ by Ax, or A- x when there is ambiguity. We call sp(A) the set of eigenvalues of a matrix A. For A and B two matrices of same size, the Frobenius inner product of and B is (A, B)F = Tr(A' B) = Ee, B,y where A AT, and the corresponding norm of A is |Allp = E,, |A For P1.P2 two independent random processes, we write pr ~ P2 when P1 and pa follow the same distribution.
Note that alternatives losses C have been considered in Delouis al. (2022); Siahkoohi et al. (2023); Auclair et al. (2024), and have shown significant improvements for the estimation of o(xo), respectively. A formal investigation of these alternatives is left for future work. hhttts://gitthn http://gghhhhhcccc//rrattat///// cwt/bbrgggd/////taatttttommo _orppa____o__pp________
2

---

Published in Transactions on Machine Learning Research (02/2024)
2 A First Analytical Exploration
As first exploration, compute the set of global minimizers of L as defined in Eq. (1) for simple examples of 6, and choosing for 6n Gaussian white noise distribution of variance o², that is p(to) ~ of This assumption for 60 places us in typical denoising framework. We also provide in App. B a discussion of how our method relates to maximum likelihood estimation in this context.
2.1 Linear Representation (xx) AAz leads to mean subtraction
To start with, we consider the case where <(x) Az, with an injective matrix of size K M (with necessarily K 2 M). The vector <(x) then simply consists in a set of features that are linear combinations of the imput vector components. The following proposition establishes that the minimizer is necessarily E (co].
Proposition 2.1. For <(x) Ar with injective, C has unique globol minimizer equal y E [co].
This proposition proven in App. A.1.1. is a simple and informative result: our representation is linear, then the minimization of can only bring us back to the observation y diminished by the mean of the noise. In particular, when the mean of the noise is zero, this optimization has no effect, prompting us to turn to a nonlinear operator Φ
2.2 Quadratic Representation <(x) = =x² leads to aqrt-thresholding
A very simple nonlinear representation is the pointwise quadratic function. Without loss of generality, we only consider the case where the dimension of x M = The solution are given by the following proposition.
Proposition 2.2. For d(x) x² and p(co) ~ oN IM),), the global minimizers of C are 0 when y² s 3o and / =3d² when y² 302.
We refer the reader to App. A.1.2 for a proof. This solution introduces a threshold on 2². If y² not sufficiently large compared the variance of the noise o², then the solution is 0. Otherwise, we obtain - 3 This second case demonstrates an attempt to "subtract the noise magnitude from the signal magnitude".
40
20
The behavior, shown in Fig. 1, is similar to that of the piecewise linear soft- thresholding operator. One significant difference is that it asymptotically reaches the identity function.
20
4a
The intuition here is that, whenever the amplitude of the signal is too small compared to that of the noise, it is optimal to shrink it to zero. However, if the signal clearly stands out from the noise, then can be corrected accordingly through a small shrinkage.
Figure 1: sqrt-thresholding function involved in Sect. 2.2 (in blue), and asymptotic behavior (in black).
Note that related shrinkage functions are common in the context of sparse signal processing (see e.g. Selesnick, 2012; Al-Shabili et al., 2021).
2.3 Power Spectrum Representation
A standard summary statistic for stationary signals is the power spectrum its distribution of power over frequency components. We consider signals of arbitrary dimension M, and set of filters ,1...., d'x typically well localized in Fourier space (i.e. "bandpass filters"). We assume that these filters cover Fourier space, meaning that for any Fourier mode k, there exists in such that &. (k) # 0.3 We define the power spectrum representation
o(z) (1lv) xr| ..... 1 * x
(3)
3 here refers the discrete Fourier transform of
3

---

Published in Transactions on Machine Learning Research (02/2024)
It corresponds to vector of K coefficients measuring the power of the input signal on the passband of each of the filters ",,...,,...
If we assume that the frequency supports of these filters are disjoint, minimizing C is equivalent to minimizing independently for all I.....K:
(x,,,) = (1)) (|
(4)
where is the orthogonal projection of the input I' E RM in the subspace spanned by the Fourier modes of the passband of th and V, is the injective matrix representing the linear operation x - thh * x in that subspace. The following proposition (proof given in App. A.1.3), determines the minimizers of Lt.
Proposition 2.3. For ((x) = ||Ax|| with A injective and p(ea) (O, I²I,), introducing:
A sp(A A) such that ||Ay|l - E [[Ac||] 220
(5)
if A = 0, then the global minimizer of L is unique equal to 0, otherwise, the minimizers are the eigenvectors æ of A A associated with min such that ||Ax||²= || [|Ac||]]
20 min A.
Let us take a step back before breaking down this result. The power spectrum a statistic that is additive when computed on independent signals, so that for a and two independent processes, we have E [o(a + b)] = E [Ø(a]] + E [o(b)] In our setting, where zo is viewed as a deterministic quantity, this gives:
Bcorp(ee) [ø(v)] #(x0) + Begrptee) [e(E)]]
Therefore, an unbiased estimator of the power spectrum statistics of xo based on the observation is:
Φ(xo) <by) - Eeg Eeoop(c)) [o(Ea)]
Now, applied to (xxx)) - 1|, Prop. 2.3 tells us that there is a threshold below which the minimization of L; leads to a signal that is zero over the passband of vi, and above which this minimization leads to a 20² min A, = ox()) - 20 min Aş. For typical filters, we usually have 2 min A, < [] so that when the signal stands out from the noise, the global minimizers of L. have power spectra statistics that almost coincide with the unbiased estimator of the power spectrum coefficients of zo. In conclusion, in this setting, the power spectrum representation leads to minimizers ft such that ((+)) is an explicit estimate of ((xo).
signal z; such that |,,,,11== = ||,y|
We note that for filters with intersecting passbands, the previous analysis becomes significantly more technical. We prove in App. A.1.4 a generalization of Prop. 2.2 that addresses this case.
3 Statistical Component Separation and Image Denoising
Previous works have shown that statistical component separation methods can perform surprisingly well for image denoising provided that the representation & suited to the data (Regaldo-Saiee Saint Blancard et al., 2021; Delouis al., 2022; Auclair et al., 2024). Although the goal of these methods remains to estimate Φ(xx), in this section, we investigate numerically to which extent dg can also be relevant estimate of zi--
We focus typical denoising setting, where eo is colored Gaussian stationary noise. We employ the block-matching and 3D filtering (BM3D) algorithm (Dabov et al., 2007; Mäkinen et al., 2020) a benchmark. However, we emphasize that contrary to BM3D, statistical component separation methods can apply similarly to arbitrary noise processes, including non-Gaussian or non-stationary ones. This has already been illustrated in the previous literature on this subject, and we will consider an additional exotic noise in this section for this purpose.
We introduce in Sect. 3.1 the vanilla algorithm used for the experiments of this section. Then, we apply this algorithm for two distinct representations: Sect. 3.2 employs representation based on the wavelet phase harmonics (WPH) statistios, and Sect. 3.3 uses summary statistics defined from feature maps of ConvNet.
We can verify that for A matrix
1, we recover result. equivalent that of Prop. 2.2.
4

---

Published in Transactions on Machine Learning Research (02/2024)
3.1 Vanilla Algorithm
Algorithm Vanilla Statistical Component Separation
Inputs: y. >((E)), Q, T, gradient-bassee optimizer (e.g. LBFGS)
Initialize: 20 y
for i= 1..... ....Tor do
sample , ~ p(eo)
( = Ei= 100 + (x) - ~(v))111/ iQ
20 t LOOOON__TTO
f.. v((oo)
end for return 00
Analytically determining the global minimizer of L for arbitrary representations & quickly becomes intractable, in which case one has to solve this optimization problem numerically. A straightforward way to do that is to approximate L via Monte Carlo estimates. Introducing Q independent noise samples , p(eo), we define:
((x) - - Σ lo(x -
(8)
The minimization of C then lakes the form of a regular stochastic optimization described in Algorithm 1, and referred to as the vanilla algorithm in the following. This algorithm was used in Regaldo-Saint Blancard et al. (2021), where the authors had employed a L-BFGS optimizer (Byrd et al., 1995; Zhu et al., 1997) using y as the initial guess. We proceed similarly in the following, and fix the number of iterations to T 30 and the batch size to Q 100.
3.2 Wavelet Phase Harmonics Representation
Definition. Similarly to Regaldo-Saint Blancard et al. (2021); Auclair et al. (2024), we consider a repre- sentation based on wavelet phase harmonics (WPH) statistics (Mallat al., 2019; Zhang & Mallat, 2021; Allys et al., 2020). These statistics efficiently capture coherent structures in variety of non- Ganssansian stationary data. They rely on the wavelet transform, which locally decomposes the signal onto oriented scales. Formally, for a random image 2, the WPH statistics are estimates of covariances between pointwise nonlinear transformations of the wavelet transform of æ. Using a set of complex-valued wavelets ,i....,,. covering Fourier space with their respective passbands, we focus on covariances of the form:
8!1 = Cov (ar 1/, x s/),, 5 (a) = Cov [| whl, x v/l,
5 (a) = Cov [|x thal lar dall, C.) (a) = Cov [|x vhl, x x
(9)
(10)
We give in App. D the technical details related to the definition and computation of these statistics. In practice, ((x) is made K = 420 camplex-valued coefficients for x a M = 256 x 256 image.
Experimental Setting. We consider three different types of 256 x 256 images corresponding to a simulation of the emission of dust grains in the interstellar medium (the dusl image), a simulation of the large-scale structure of the Universe VVllaescusa-Navvrre et al., 2020) (the LSS image), and randomly selected images from the ImageNet dataset (Deng et al., 2009) (the ImageNet images). We give additional details on this data in App. C. We then consider four different noise processes: three colored Gaussian noises, namely pink, white, and blue noises, as well as a non-Gaussian noise made of small crosses (see Fig. 2, bottom right). We vary, for the colored noises, the amplitude a of the noise considering 10 different levels ranging from 0.1 to 2.14 (logarithmically spaced) in unit of the standard deviation of x, and for the "crosses" noises, the density of crosses p considering 10 different values ranging from 0.001 to 0.063 (logarithmically spaced) For each of
We have found empirically that the batch size should be kept sufficiently high for the L-BFGS optimizer to behave correctly. We also report that T 30 was sufficient achieve approximate convergence for experiments. 6The density p is defined the expected ratio of crosses the total number of pixels.
5

---

Published in Transactions on Machine Learning Research (02/2024)
Pink noise
White noise
Blue noise
"Cossee"" noise
Nonay 18.5 dEE)
Noiey 18.5 AT)
Noay 18.5 dP)
Noney 300. dB)
Original
Denaised (onz)) 25.8 aB
Demaaood (onus) 30.2 dT)
Demcotonn (cunre) 322.2 aT..
Demcttsc (onns) 27.7 ABB
50-
100-
150
200
250
Detoibett (BMND) [26.5 dB)
Deuiviodd (BMAD)) (a 22 dT)
Denoload (RMDDD 184. Is (TF)
100 200
Figure 2: Original dust image zo (center left), noisy realizations y for distinct colored Gaussian noise processes and a non- Gaussian noise (top row), denoised images èo using Algorithm with the WPH representation from Sect. 3.2 (middle row), and denoised images using BM3D) for colored noises (bottom row, except bottom right). Additionally, sample of the non-Gaussian "crosses" noise is shown (bottom right). The accompanying PSNR values are provided next to each noisy and denoised image.
theses cases, we apply Algorithm 1 for 30 different noise realizations. Each optimization takes 40 with a GPU-accelerated code on A100 GPU.
Results. In Fig. 2, compare the noise-free dust image xo with examples of its noisy y and denoised to versions for each noise process. In the case of the colored noises, we also include the BM3D-denoised images in the bottom row for reference. Our method effectively reduces noise while preserving the original image's structure. Even with the "crosses" noise, our algorithm successfully removes most the crosses, the remaining ones having been most likely confused with actual structures of zo- We evaluate the quality of the denoising in terms of peak signal-to-nois noise ratio (PSNR), significantly improving it all cases. However, BM3D outperforms our method for colored noises, which is not surprising as our approach was not explicitly designed for that.
We further evaluate our algorithm's performance for colored noises using PSNR and the relative error of ((tp). In Fig. 3, we present the PSNR and relative errors of WPH statistics for different coefficient classes as a function noise level a. Our method consistently improves PSNR, but BM3D performs better across all noise levels. We note that our method's performance degrades for very high o, potentially due to structure hallucination caused by extreme initial noise.? 7 In terms of WPH statistics relative error, our method effectively reduces noise impact and outperforms BM3D in most cases. Notably, excels in I and 5 coefficients, except for the high-noise regime in 811, However, for go0 and C coefficients, in the cases of blue and white noise, BM3D performs comparably or better than our method. Since a perfect PSNR implies perfect statistics recovery, we interpret this as the sign that a regular denoising algorithm, when adapted to the noise process, can also provide precise estimates of Φ(2)). However, we point out that the normalization of the WPH statistics may play a crucial role on these metrics (see App. D), and a fairer comparison should explore
WPH statistics may alao define generative models with sampling procedure sharing important similarities with Algorithm (see Allys et (2020); Zhang & Mallat (2021); Ragando-Saint Blancard et (2021): Régaldo-Saint Blancard al. (2023)).
6

---

Published in Transactions on Machine Learning Research (02/2024)
104
102-
Demoised Piak moise (ours)
Demoised Pink noise (BMMSD)
108
Rel. ceron 5t0
10"
Demoised
(ouns]
Rel. estor
Denoised White noise (BBMDD)
10
10-1
Demoised Blue noise (ours]
Demoised Bine noise (BM3D)
10
Higher betler
Lonee is better
10
109
PSVR (6D)
Rel. enter
10
Rel. orvor
10-14
10
10
10
10
108
Noise level
Noise level
Noke level
Figure 3: PSNR and relative errors of the WPH statistics per class of coefficients as function of the noise level for the denoised dust images described in Sect. 3.2 for each lype of colored noise.
10"
109
Rel. orror 500
Noiey Crosses
Rel. errer
Demoised Crosses (ours)
10-1
10-1
Higher better
10
Lower is better
35
10P
PSNR (1B)
[tel. eroor
Rel. error
25
20
15
10
10a2
10
10
Demsity of
Density of p
Density of
Figure 4: Same Fig. for the "crosses" noises as function of the density of crosses p.

---

Published in Transactions on Machine Learning Research (02/2024)
its role in this precise setting. We report similar results in Fig. for the "crosses" noise. These demonstrate clear mitigation of the noise in all regimes.
We report in Figs. F.1, F.2, F.4, F.5, and F.7 equivalent results for the LSS and ImageNet data. The overall conclusions are the same.
3.3 ConvNet-based Representation
Definition. To further explore the dependence on the representation & for the performance of statistical component separation for image denoising, we define representation based on feature maps of ConvNet. We make use of the VGG-19 BN network (Simonyan & Zisserman, 2015) which was trained for image classification on ImageNet. We only employ the first two convolutional blocks of the networks (each of these being made of a sequence of Conv2d, BatchNorm2d, and ReLU blocks), which, for each 224 x 224 RGB image 2, output a set of 64 224 x 224 feature maps f(z) = (fi(z),.... f6)(.) We define a representation d(x) of dimension = 64 from these leature maps by laking the squared Euclidean norm of each feature map, that.
<(x) -.../((")) (x))1)
(11)
This representation shares structural similarities with the WPH statistics and its parent the wavelet scattering transform statistics (Bruna & Mallat, 2013; Mallat, 2016). However, contrary to these previous transforms which employ generic wavelet filters, the filters of the VGG convolutional layers are specifically trained for the analysis of ImageNet data.
Results. We conduct the same experiment in Sect. 3.2 in this setting focusing on randomly selected ImageNet images. We show in Fig. F.8 examples of resulting images before and after denoising, and we show in Fig. 5 the root-mean-square error on the coefficients of the representation & as a function (mean and standard deviation computed across randomly selected test. images). While Fig. 5 indicates that our method significantly reduces the impact of the noise the representation &, regular denoiser in comparison to BM3D. This is interesting result, as it shows that the representation & we have built here is not as well suited for image denoising as the WPH representation was. Given the fact that VGG networks were trained for a classification problem, it likely that their feature maps are partly robust to the noise in the input. A spectral analysis of the denoised maps further shows that there remains significant amount of noise in the small scales, suggesting that the feature maps weakly impacted by the noise at small scales.
19
Denecned Pink noiee (ours)
Denotzod
Pink noiee (DMID)
White (ouns)
Doniasd White máse (BMBL)
Domiadd Blu miies (ouns)
Demilad Blue moise (BMR))
10-3
in Fig. F.8 that it performs poorly
10
Lomerr bester
10
Notse lcoel
Figure 5: Root-mean-square error on the coefficients the ConyNet- based representation & as a function of the noise level for the denoised ImageNet. images for each type of colored noise.
4 Diffusive Statistical Component Separation
Taking inspiration from the diffusion baatdd generative modeling literature (Ho et al., 2020; Song al., 2021), we investigate to which extent statistical component separation methods can benefit from the idea of breaking down the optimization into a sequence of optimization problems involving noises of smaller amplitude.
We introduce Sect. 4.1 a new algorithm that leverages this idea in the case of Gaussian noises and apply it the dust data introduced Sect. 3. We then study in Sect. 4.2 the limit regime where the amplitude of the noise tends zero. This gives an alternative way perform statistical component separation giving promising results in specific circumstances.
8

---

Published in Transactions on Machine Learning Research (02/2024)
4.1 A "Diffusive" Algorithm
Stable process. We assume that eo stable noise process, that is for any ,l,.,,,Ep i.d. p(Eo), we have Ei= 6, ~ 0ee + 3 for some scalar constants a and 3. A practical example is that of Gaussian processes, since for eo (0, E), we clearly have Σ 6 Peo (o, PE). Introducing E R such that aş 0 for all i, and ||a|| - Σ, až = 1, and provided that IE [ca] O, we can break down eo into "smaller" independent noise processes follows:
P
co Σ
a.g, with CI....,. add p(co);
(12)
where the variance of a,E; can be made arbitrarily small by taking sufficiently small value for ag.
Algorithm. In this setting, we leverage this decomposition by breaking down the minimization of C into "simpler" optimization problems involving noises of smaller variance, with the goal of finding a better optimum dg. We introduce Algorithm 2 for this purpose. This algorithm starts from &p = y and builds a sequence of signals fr-,,...,, to such that for all E (P- - 03:
i E arg min C(x; 0441, 1441) = Benplea) [lo( + as411) - of (B(+1))1²]
(13)
Algorithm 2 Diffusive Statistical Component Separation
Input: y, p(eo), Q, T, P, a E IR" with ||a|] = 1 and a; > 0, gradient based optimizer Initialize: åp = y
for i= P- do
4; = 8:+1
for j = 11.... T do
sample c EQ p(Eo)
((a.) -d (1o(E; + 0i4 - 6(1;+1))1 /Q
4; + ONE STEP OPTIM ii, vi(+:)
end for end for return to
A sufficient condition for perfect reconstruction to be achieved is that f; zo + é with é ~ - 1/2 We do not expect this strong condition to hold for arbitrary functions Ø but the design of Φ should be guided in that sense.
We also note that a stepwise approach had also been employed in Delouis et al. (2022), and we draw some connections between the two approaches in App. E. We show that these two approaches rely on related objective functions. A further exploration of the pros and cons of each of these two algorithms is left for future work.
Experiment. We apply Algorithm 2 to the dust data introduced in Sect. 3 in the case of the Gaussian white noise and the WPH representation with a; 1/VP and P= [10r]. We compare in Fig. the results of this algorithm to those of Algorithm 1. We see that this approach slightly improves the results for every metric except for the relative error on the Ta coefficients al an intermediate noise level where these are slightly deteriorated. Although these numerical experiments are not showing significant improvements, having at least consistent results suggest that, from a thcoretical perspective, component separation methods can be understood as the aggregation of these small optimization problems. We push this idea further in the rest of this section.

---

Published in Transactions on Machine Learning Research (02/2024)
10"
Denoised (ours Alg. 1]
Rel error gm
Hel. errie 540
10
Denolsed (ours Alg. 2)
10-1
Denoised (ours Alg. pert.)
Denoised (BM3D)
10-
10-2
Higher is better
Lower betser
45
40
Rel. error 50
Hel. (.01
PSNR (AB]
10
10
25
10
10
10P
10
10"
10"
Noise leseel
Naise level
Noise lovel
Figure 6: Same Fig. 3 for the white noise only, and comparing results obtained with Algorithm 1 (see Sect. 3.2), Algorithm 2 with non-perturbative loss (see Sect. 4.1), Algorithm with the perturbative loss introduced in Sect. 4.3, and BM3D.
4.2 Limit Regime for Infinitely Small Noises
We study the limit regime of infinitely small noises to give more formal perspective to Algorithm 2. We introduce f(x, a) = [| + ac) - +(y)11 and expands with respect to a E RR
Proposition 4.1. For d twice differentiable function, and to arbilrarily distribuled with E [co] 0, we have:
((x,a) -() [(J(())
E)r + (Ho(1), E)) (6(x) d(v))] + o(a?), (14)
where (x) is the Jucobian matria of (of size K x M), (x) is ils Hessian lensor (of rank 3, and size K M M), and E is the covariance matria of 6o.
The proof is given in App. A.2.1. The zeroth-order term of this expansion has a clear interpretation. It prevents the representation of æ from moving too far from that of y. However, the second order term is more intricate as combines first and second order derivatives of a, which are weighted by the covariance of the noise. As an example, let us consider ==o We then get:
(Jo(z) (z()). E/ = o² 12 (H&(r), E/ (6(z) o(v)) = 2 [H)())) (6(x) x(y)).
(15)
(16)
The first term is proportional to the squared norm of the Jacobian matrix, while the second term a dot product between the vector Hessian traces with the vector <(x) - d(y). The trace of Ha, (x) directly relates the mean curvature of the function d al point x, so that this dot product quantifies the alignment. between σ(x) Φ(y) and the vector of mean curvatures for each component of
We denote (H)( Fhe the vector K such that component, equals (Hay (s)), F.)e.
10

---

Published in Transactions on Machine Learning Research (02/2024)
4.3 A Perturbative Algorithm
We note that contrary to Eq. (1), Eq. (14) got rid of the expected value over e, and the second-order expansion of [(x, a:) now only depends on the covariance matrix E of the noise, which is often known easy to estimate from a set of noise samples. Then, provided L(x, a) can be correctly approximated by its second-order expansion in a, one can evaluate the loss in a much more straightforward way that does not rely on a noisy Monte Carlo estimate. In this limit regime, the computational challenge lies in efficiently computing the socond-order term of Eq. (14) in a differentiable way.
We compute in App. A.3 the relevant analytical expressions to do so for o the WPH representation introduced in Sect. 3.2. We then apply Algorithm 2 for the dust image in the Gaussian white noise case with a, 1/VP, P 100r], and T = 10, and by replacing the loss L with by its truncated sccond-order Taylor expansion as explicited in Eq- (14). We empirically found that the most stable setting is the one where the WPH statisties only include the St1 and S coefficients. We report in Fig. the quantitative results this setting for the PSNR and relative errors on these coofficients as function of a. Our method performs remarkably well in comparison to the previous results. It is very close to BM3D in terms of PSNR and outperforms in the high noise regime in terms of relative errors the S11 and ş01 coefficients. The fact that the relative error for the 5 coefficients is higher than for the noisy data in the low noise regime however suggests a form of instability for this range of a. Nevertheless, we find these results promising and clear demonstration of the relevance of this "diffusive" perspective on statistical component separation methods.
5 Conclusion
This paper has explored several aspects of statistical component separation methods. Section 2 has exhibited analytically the global minimizers of L for several examples of representations & in the case of Gaussian white noise. We have shown that a linear Φ cannot, extract any information representation leads to a form of sqrt-thresholding of the observation y. For 6 a power spectrum representation, which can be viowed as a more general quadratic representation, the minimizers of C lead to relevant estimations of (( Then, in Sect. 3, we have approached numerically the minimizers of L introducing Algorithm an image denoising setting for two representations where analytical calculations are intractable: 1) WPH statisties, 2) statistics derived from the feature maps of ConvNet. For 1), Algorithm acts as a regular image denoiser while not explicitly constructed for this. Although it does not outperform BM3D in terms of the PSNR metric, it better recovers the coefficients ((za) for most classes of coefficients and experimental settings. Additionally, it may extend to arbitrary noise processes as it was illustrated with an exotic noise made small crosses. For 2), the impact on the noise was clearly mitigated in o(fo). but the resulting images dg were still very noisy, showing that this representation is less suited for image denoising. Finally, in Sect. 4, we have introduced Algorithm 2, "diffusive" statistical component separation method that can be applied in contexts where the noise is a stable process. These ideas led in some cases to better results than Algorithm in denoising setting. And more importantly, it supported the idea that statistical component separation methods can be described as the sequence of optimization problems with noises of smaller amplitudes. This idea will be pushed further in future work.
No, while a simple quadratic
Acknowledgments
It is a pleasure to thank Erwan Allys, Jean-Mare Delouis, Stéphane Mallat, Loucas Pillaud- Pillaud-Vivien, and Léo Vacher for valuable discussions.
References
Jonas Adler, Sebastian Lunz, Olivier Verdier, Carola-Bibal Bibiane Schönlieb, and Ozan Öktem. Task adapted reconstruction for inverse problems. Inverse Problems, 38(7):075006, July 2022. doi: 10.1088/1361-642)/ ac28ec.
Abdullah H. Al-Shabili, Yining Feng, and Ivan Selesnick. Sharpening Sparse Regularizers via Smoothing. IEFE Open Journal of Signal Processing, 2:396-409 409, 2021. doi: O0.1109//JSP:2021..2 2000SP.2021.3111
11

---

Published in Transactions on Machine Learning Research (02/2024)
Erwan Allys, Tanguy Marchand, Jean-François Cardoso, Francisco Villasscusa-Navaroo, Shirley Ho, and Stéphane Mallat. New interpretable statistios for large-scale structure analysis and generation. Physical Review D, 102:103506, November 2020. doi:
Constant Auclair, Erwan Allys, François Boulanger, Matthieu Béthermin, Athanasia Gkogkou, Guilaine Lagache, Antoine Marchal, Marc-Antoine Miville-Deschènes, Bruno Régaldo-Saint Blancard, and Pablo Richard. Separation of dust emission from the cosmic infrared background in Herschel observations with wavelet phase harmonics. Astronomy & Astrophysics, 681:A1, January 2024. doi: 10.1051 000044600066611 202346814.
Joan Bruna and Stéphane Mallat. Invariant scattering convolution networks. IEEE Transactions on Pattern Analysis and Machine Intelligence, 35(8):1872-1) 1886, 2013. doi: 10.1109/TPAMI.
Joan Bruna and Stéphane Mallat. Multiscale sparse microcanonical models. Mathematical Statistics and Learning, 133))257-312 2019. doi: 10.4171/msl/7.
Richard Byrd, Peihuang Lu, Jorge Nocedal, and Ciyou Zhu. A Limited Memory Algorithm for Bound Constrained Optimization. STAM Journal on Scientific Computing, 16(5):1190-1208, 1995. doi: 10.1137/ 0916069.
Jean-François Cardoso. Blind signal separation: statistical principles. Proceedings of the IEEE, 86(10): 2009 2025, 1998. doi: 10.1109/5.720250.
Kostadin Dabov, Alessandro Foi, Vladimir Katkovnik, and Karen Egiazarian. Image denoising by sparse 3D monsform-domaite collaborative filtering. IEEE Transactions on Image Processing, 16(8):2080 2095, 2007. TSSN 22195491.
Jean-Marc Delouis, Erwan Allys, Edouard Gauvrit, and François Boulanger. Non-Gaussian modelling and statistical denoising of Planck dust polarisation full-sky maps using scattering transforms. Astronomy & Astroplegsics, A68:A122, 2022. doi: 10.1051 (0.101100004-661/2000 /20224456..
Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei. Imagenet: A large-scale hierarchical image database. In 2009 IEEE Conference on Computer Vision and Pattern Recognition, pp. 248 228, 2009. doi:
Michael Elad, Bahjat Kawar, and Gregory Vaksman. Image Denoising: The Deep Learning Revolution and Beyond-A Survey Paper. SIAM Journal on Imaging Sciences, 16(3):1594-1654, 2023. doi: 10.1137/ 23M1545859.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising Diffusion Probabilistic Models. In Adoances in Newral Information Processing Systems, 2020.
Julien Mairal, Francis Bach, and Jean Ponce. Task-Driven Dictionary Learning. IEEE Transactions on Pottern Analysis and Machine Intelligence, 34(4):791-800 2012. doi: 20.109//PAAM..01..5555
Stéphane Mallat. Understanding deep convolutional networks. Philosophical Transactions of the Royal Society of London Series A, (7((6)6)):2010000 April 2016. doi: 10.1098/rsta.2015.. 2atita..01500033
Stéphanc Mallat, Sixin Zhang, and Gaspar Rochette. Phase harmonic correlations and convolutional neural networks. Information and Inference: A Journal of the IMA, 933):721-747, 11 2019. doi: 10.1003//mmal iaz019.
Ymir Mäkinen, Lucio Awari, and Alessandro Foi. Collaborative Filtering of Correlated Noise: Exact. caansform-Dmmalle Variance for Improved Shrinkage and Patch Matching. IEEE Transactions Image Processing, 20:8339-8354, 2020. doi:
Bruno Regaldo-Saint Blancard, Erwan Allys, François Boulanger, François Levrier, and Niall Jeffrey. A new approach for the statistical denoising of Planck interstellar dust polarization data. Astronomy 8 Astrophysics, 649:L18, 2021. doi: 10.1051 /0.1051/0004-3361/ /200211400322
12

---

Published in Transactions on Machine Learning Research (02/2024)
Bruno Régaldo-Saint Blancard, Erwan Allys, Constant Auclair, François Boulanger, Michael Eickenberg, François Levrier, Téo Vacher, and Sixin Zhang. Generative Models of Multichannel Data from Single Example-Applicttion to Dust Emission. The Astrophysical Journal, 943(1):9, January 2023. doi: 10.3847/ 1538-4357 2588-435//cc5586 aca538.
Ivan Selosnick. Penalty and Shrinkage Functions for Sparse Signal Processing. Connerions, 11.22:23, 2012.
Ali Siahkoohi, Rudy Morel, Maarten V. De Hoop, Erwan Allys, Gregory Sainton, and Taichí Kawamura. Uncarthing InSights into Mars: Unsupervised Souree Separation with Limited Data. In International Conference Machine Learning, 2023.
Karen Simonyan and Andrew Zisserman. Very deep convolutional networks for large-scale image recognition. In International Conference on Learning Representations, 2015.
Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stelano Ermon, and Ben Poole. Score-Based Generative Modeling through Stochastic Differential Equations. In International Conference on Learning Represenlations, 2021.
Jcan-Luc Starck, Michael Elad, and David L. Donoho. Image decomposition via the combination of sparse representations and variational approach. IEEE Transaclions Image Processing, 14(10):1570 1582, 2005. doi: 10. N0.1109/TPP.005..1. 2022.8520..
Francisco Mataescussa-Navaryy ChangHoon Hahn, Elena Massara, Arka Banerjee, Ana Maria Delgado, Doogesh Kodi Ramanah, Tom Charnock, Elena Giusarma, Yin Li, Erwan Allys, Antoine Brochard, Cora Uhlemann, Chi-Ting Chiang, Siyu He, Alice Pisani, Andrej Obuljen, Yu Feng, Emanuele Castorina, Gabriella Contardo, Christina D. Kreisch, Andrina Nicola, Justin Alsing, Roman Scoccimarro, Licia Verde, Matteo Viel, Shirley Ho, Stephane Mallat, Benjamin Wandelt, and David N. Spergel. The Quijote Simulations. The Astrophysical Journal Supplement Series, 250(1):2, sep 2020. doi: (0.3847/1538-4455//
Sixin Zhang and Stéphane Mallat. Maximum entropy models from phase harmonic covariances. Applied and Computational Harmonic Analysis, 53:199-230, 2021. doi: https://doi.org/ http///i..oo/11111/////33111
10.00.0000......
Ciyou Zhu, Richard H. Byrd, Poihuang Lu, and Jorge Nocodal. Algorithm 778: L-BFGS-B: Fortran Subroutines for Large-Scale Bound- Conde-Conttr Optimization. ACM Transactions Matheratical Soflware, 23(4): 560--60, 1997. doi: 10. 279149969992..292236
13

---

Published in Transactions on Machine Learning Research (02/2024)
A Proofs
A.1 Proofs of Sect. 2
We give the proofs of the results presented in Sect. 2. Assuming that eo is a Gaussian white noise of variance o², that p(eo) ~ of we compute the set of global minimizers of C as defined in Eq. (1) in the following cases:
(xx) = Ax with A an injective matrix of size K x M,
2. (xx) x² for M 1,
((x) = ||rl with an injective matrix of size K' x M,
4. ((x) (I1,11) ..... ||Axl))) with A....., Ax K injective matrices of size K' M such that
A are co-diagonalizable.
For each case, note that L(x) is infinitely differentiable and obviously bounded from below by O, and will show below that lim C(x) too to demonstrate the existence of global minimum.
We will then determine the global minimizers by studying the zeros of the gradient of C, which, in its general form, reads:
v ((xx) = 2B.-p(eo) [J + +)" (ó(z + -(v))]
(17)
where Jø the Jacobian matrix of Φb
A.1.1 <(x) Ax with A injective
Proposition 2.1. For ((x) = Ax with injective, C. has a unique global minimizer equal to E [co].
Proof. We have:
=E [
(18)
If A injective of size K x M (in which case, we necessarily have K 2 M), A A is positive-defiite and we call Amin > 0 its smallest eigenvalue. We have ||Arr|| Amin ||~]| so that:
((x)= AminE [| + v|²] Amin ||]|| + O(I==])
(19)
This shows that lim C(x) too so that global minimum does exist.
Now, noticing that (o(z) A, have:
VL(x) = 2A" E [o(r + d(v)],
(20)
= 24 d(B [x + v)),
(21)
= 244". -0r - y + E [e]),
(22)
= 2AT A (r - + 1 [e]).
(23)
For injective A, AT A is invertible, and we have vL(x) = =0 x z= EE [c]. This shows that, in this case, the only global minimizer is y E [c].
Note that this result is independent of the noise distribution p(Eo).
This result can be extended to non-injective matrices A. In this case, there is an affine subspace solutions with the same loss value along shifts within the null space of A. By similar procodure as above, we obtain that = y [c] + v, where E AerA. Non-injectivity does not change the overall conclusion that, besides the subtraction of the mean of the noise, denoising does not take place.
14

---

Published in Transactions on Machine Learning Research (02/2024)
A.1.2 Φ(x) = = for M =
Proposition 2.2. For d(x) = x² and p(eo) ~ ~N(0 Ims), the global ménimizers of are 0 when 12 S 30";
and yV - when y² > 30²
Proof. For xxx) = x² and K = M = 1, we can show by expanding C that:
[(x) |-س|| 0(11=113)
x++00 to0.
(24)
This proves the existence of global minimum. Now, injecting ('(x) 2. in Eq. (17), get:
(x) 4E [ + ((x ,²)],
(25)
4E [ + c) (x (x+ c)y²
(26)
= 4E [I + é + 3 + 32E² - xy²² ey²]
(27)
= 4 (2 + xx30² - r²)),
(28)
where have used the fact that E [c] E [c]] 0 and E [c]] 0². The zeros of L are thus 0
x = + /1² 3 when y² 2 30 By looking at the values of ("(x)) = 12x² 2 + 12a² - 4y² at these points, we find that whenever 3² > 30" (0) < 0 so that 0 cannot be a local minimizer. In this case, we only have + / y² 30 left to be the global minimizers.
All of this shows that the global minimizers of L are when y² s 3o², and = y²² 30 when y² 30²
This denoising function sgn(y) / max y² 30 is similar to the soft-thresholding function stx stx(v) sgn(y) max(ly - Al, 0), sharing the flat region around the origin. The main difference that it contrary to soft-Uhresholding, this function approaches the identity function for high values. It also has infinite derivalives the border of its flat region.
A.1.3 <(x) ||x with A injective
Proposition 2.3. For ((x) ||xl|| with A injective and p(eo) o2 IMI, introducing:
A E A) such that ||Aull E ||Ac||]]
20 2 0},
(29)
if A = 0, then the global minimizer of L is unique equal to 0, otherwise, the minômizers are the eigenvectors x of AT A associated with min such that ||Az|| ||Ayll [||A||]] - -2 min A.9
Proof. With an injective matrix of size K' x M, we have:
[114C0+ + e)( [[[([=+)]--1AAy]1] -1|Ay11)
[(x)
(30)
[IA(x + e)1l +||Ay|* |1A( e)l ||Al|
(31)
Writing that ||A(x e)ll || + || 2 + 2Ax Ae, show that:
E
|A(z + +)]1?]
= || + E ||AE|]
(32)
E
|A(z e\I1
||ax||
+
||Ac||
+ 400
|| Az|
+ 2
[cc]
(33)
We verify that matrix
1,, we recover result equivalent, that of Prop.
15

---

Published in Transactions on Machine Learning Research (02/2024)
where we have used the facts that E [Ac] 0 and E [ec]] = We can then rewrite L as follows:
[(x)= = || + 2 || Axr 2 (R: [||4||]]] [lvll")
L()).
(34)
Same as Sect. A.1.1, with Amin min sp AA A) 0, where sp(A" A) the spectrum of ATA, we can show that:
( 2 X |
0(1==1)),
(35)
which leads to lim [(x) too and proves the existence of a global minimum. 200+00
Starting from Eq. (34) and using the fact that VØ(x) = 2.A Ax, we show that:
v((x) (||A1]]² + [|Ae]]]] - -1Ay)1 A A 80 (A
(36)
Since A is an invertible matrix, we then have:
VL(x) =0 = (a(x) + A² O.
(37)
with a(x) 2/2 20 || + ||Ae| |]] ||Avll").
Let us take a non-zero solution x of VC 0. Then, the matrix a<x) + A must be singular, which is equivalent to the fact that -a(x) X sp(A"AA A). Being zero of V. leads to A Ax Xx, that
| = x" AAArr Ar = A ||]| and then ||A" Axll A || By injecting this to Eq. (34), we finally get:
((x) .(0) - ||Ar||
(38)
[()) (l- - EE [|AA[|]] 20²2)²
(39)
Therefore, the global minimum is reached when is the smallest eigenvalue of AT A constrained by ||=l|| = ||yyll - E ||]]]]] - 20 2 0. Reciprocally, introducing
A E sp(A")) such that ||Ay|| E [|Ae||]] 20 2 0}.
(40)
we verify that, if A = 0, the global minimizer of L is O, and if A + 0, the global minimizers are cigenvectors æ of AT associated with min A such that ||0||² min [|yy]] - E ||Ae|| 20
A.1.4 <(z) (1|41111,... ..... 1A. with A Ax injective and AT ... AK Ak co-diagonalizable
For any matrix of size M and a subset Eq [k1 <... < lio) of (I....,M) we denote by A the Q x P submatrix of A such that [A = Ak... We also denote by A+ the Moore-Pemne inverse of the matrix A. Finally, we denote by span(e1...., en) the subspace generated by the vectors eI....,
Proposition 2.4. We consider d(x) = p(eo) ~ oT)) Introducing SA. = AT Ai, we assume that SAp,..,, ,AA are co-diagonalizable. We choose an orthonormal basis (e.....,Ems) to co-diagonalize them (there always exists one). We call = (N..))sssM.II/S5 the corvespouding matris of eigenvalues, where Ais is the eigenwalue of Sa, as- sociated with eq. We also assume that A is of rank K and that any submatriæ of A with size x K is invertūble.
with Ak injective and
MM))
In that case, introducing:
K
A = (Ai...,, AM) with A, = Σ 1=1
(A - E [IA)E|]]] (1 - 2022
(41)
T (=c (1, MD sisis
0).
(42)
T-= 0, then the global minimizer of is unique equal to otherwise, the minémizers are the r E span(e)) satisfying (1|4,1112.. ..... [|xx))?) = A Ae where E argmax
|
Ser
16

---

Published in Transactions on Machine Learning Research (02/2024)
Proof. We have:
K
((x) E
[| |1 + +]] 1
(43)
K
Σ
E
|| e)1l + +1AA - 2 11 + e) ||
(44)
Using the facts that E [Are] = O, IR Ase ||Ae|l]] = 0, and E [ec o²Ias, this loss reads:
((x) = Σ || + 40 ||A, ||A,A,01| + |A,x) (E [[AselP]] - [|A,v(²)] L()).
(45)
The existence of the global minimum is guaranteed for the same reasons as Prop. 2.3.
The expression of Vx(x) reads:
K
VL(x) Σ
4 (1A.2 +E [1IAce[]= ||A A7 Axx + (02 A.)
(46)
and we have:
K
VC(x)= =0 =
=0
Σ (.,()))A =Sy x = 0,
(47)
where we have introduced:
a. aA((F) (1A4:0112 E [|AAcc]]] - -|A /(202).
(48)
SA. A7 Ar.
(49)
We notice that:
VL(x) = 0 = x V[(x)= =0,
(50)
N
N
x.
T
S4 x =
Σ
a. (x)x²" SA.E.
(51)
i=1
N
|SA
aa (1)
(52)
Therefore, with global minimizer of L(x), we have:
K
((x) = Σ
||a1 - 400 -40--40-2000 + +2 |Aai|
(E [1A ||A L(0),
(53)
1=1
Σ
|
| 402a,,
2E [|Atel[]] |A:1) + C()),
(54)
k
C()) Σ
|
(55)
We take non-zero global minimizer -MM x.C,. Introducing Eq {i E (l,..., x, + 0), for all E Eq, Eq. (47) leads to:
Σ (I) + +X2,- =0,
(56)
17

---

Published in Transactions on Machine Learning Research (02/2024)
K
K Avs E [1a,c(P] - 2222220)
(57)
==
3=1
ANNal = -Ai
(58)
where we have introduced:
Az = (||Aaæ|12 ,...,, -IAxx11)))))))
(59)
K
A, Σ (Mes(llAgull [I[Asc(P]]
[1AAE(?] -00--211 [1AA)]|P)-2000 20022.))
(60)
In a more compact form, we have:
NAx A.O
(61)
By assumption, if Q 2 K then the K x K matrix AT Azo is invertible, and if Q < K, then the Q x Q matrix Asg AT is invertible. Using the fact that NA A Nz.50 where N₂ = (2î. 2M)", we show that:
(A Asa ) AA Asq if Q K,
NA
(62)
AT
(Aso AT Asqi if K,
that is:
NA = A AEo
(63)
We note that necessarily, all components of A A=a positive. By injecting this previous expression in C(x), we get:
((x) [()) ||N
(64)
C(0) A AEO
(65)
Reciprocally, introducing:
T= T= (BC (EC M) 14 sisiski (, 20)
(66)
we verify that, ifT = 0, the global minimizer of C is 0, otherwise, the global minimizers are the x E span(e))
satisfying NAr A As where E e argmax
A Ag
žer
For K = 1, we check that we recover Prop. 2.3.
A.2 Proofs of Sect.
A.2.1 Proof of Eq. (14)
Proposition 4.1. For à a twice differentiable function, and eo arbitrarily distributed with [co] 0, we
have:
C(x, a) = 11(()) -(y)) a² [J)(())][(()) E)p + (Ho(*), E)r (6(x) (())) + o(a?), (67)
where (x) is the Jacobian matriz of o (of size K x M), H& (x) is its Hession tensor (of rank 3, and size K M x M), and E is the covariance mairia of cn-
18

---

Published in Transactions on Machine Learning Research (02/2024)
Proof. We introduce R, and assume that is twice differentiable, and that cu is arbitrarily distributed with I [oa] = =0 In this context, we expand <(x + (6) as function of I at second order as follows:
ol + ac) = d(x) + aJo())e+ 1 a Ho()) o(a²)
(68)
where Jz(z)= = (bb (1)) (.)
is the Jacobian matrix of & (of size K x M), and Hz(z) = (0 8'
(x)
is a
02,0xx
ij
Hessian tensor of rank 3 and size K x M x M. The second order term must be understood as contraction on the second and third indices of the Hessian tensor.
Let us propagate this expansion in:
((I, a) -emp(ee) |106 + ae) -()))
(69)
=E -B -|1060)- - x(v)) aJo(1)e 2
d Ho()) +
(70)
= -l(((()) <(v)) a'E [1 2a (0(z) o(v)) E [Jo(z)d
(71)
+ - d()) E (o)H + o(a)))
(72)
= ((((z))((v))) + a² [E [ + (o(x) Av)) E [e" Ho(x)e)) + o(a²),
(73)
where we have used the fact that E [J)(z)d = (z) E [c] = 0. Introducing the covariance matrix E of e, we verify that:
E [ E. [[())
Jo(x) (L(()) E e Ha(x)e) = (x)
(74)
(75)
so that:
((z,a)= -())) az (]o(E))
E} + + (x),)))))) (6(x) <(x))] + o(a2).
(76)
A.3 Proofs of Sect. 4.3
We compute the explicit expressions of the Jacobian matrix (x) and the Hessian tensor Ha(x) for the operator giving the WPH statistics employed in Sect. 3.2. We break down these computations by first computing first and second-order derivatives for simpler & functions, namely 6(z) = v * and Φ(x) | + z/|- We assume periodic boundary conditions, so that for any e E KMM, we formally manipulate T E K defined by ö[i] uii mod M] for any i E Z, and the 'tilde' symbol is omitted for convenience. Below, the filters t, eh. and 1/2 are assumed to be complex-valued filters. For a filter v. we define the adjoint filter æे by [ = wi-i] The convolution operation corresponds to the periodic convolution. For z E C², we introduce sg(=) = =/|z|- The notation (x) refers to the average over the components of x, that is (a) = & El x4. We denote the element-wise product of x and by ry, or x ④ when there ambignity.
Derivatives of ((x) + * Since @ is linear, we simply have:
J사(x) i Hz(x) 0,
(77)
(78)
that is:
0
+ x((]) = bb12 - --],
(79)
dx,
(bb * x = 0.
(80)
0r Ox
Note that we always redefine 2g 30 + [c]] for this be the case.
19

---

Published in Transactions on Machine Learning Research (02/2024)
Derivatives of (xx) 11b æ|
Where o(x) is nonzero, we have:
(10 * =lii)) 0
V * x [[]
(81)
dxy
1
(Re(1 [i] Re (1) jl + Im( x>]((] Im())] - -)),
(82)
1/b z [i]
1
Re (v * x x|i|p] ji)
(83)
1/ * x [
Re (sg(v)
[d - []]
(84)
which leads to:
(z) = Re [(gg v 2) o (1b ((++-)]
(85)
The computation of Hx)x) now demands to derive x + sg( * æ)):
θ 0xz
+ * x [i]
(sg (v * xd))
vb * ar
1bb 2: ==((],
(86)
Dzy
1//
1/b
2
[1]
Oer
1 - 5
v * 20
Re (sg(1h * ( - 5l)
(87)
10b 20 [0][]]
J * [
(x
1b 51
[1]05
(88)
2/2b * x<[i]]
1 * xx|
[vblt - jl
sg(v * +x) " =
(89)
2 2/ x]|l
Therefore:
82
(10 x x(()) = Re
(sg(1) * xi]( vo(2 --]
(90)
1
Re pi
klbli jl
* [
AJ --]]
(91)
2 * x
Derivatives of 6(x) = (1bb 2013)
Using the above, we get:
0
a
(1pb * 2x = Or,
(1b (1
x 1b xb)
+ axz
(10 xx*x))= x 1o xx).
(92)
dry
M
2
Re
Σ 1=1
bbli - 5l x 2o * 2111
(93)
M
2 M
Re (wit * e * x) [i,
(94)
and:
82
2
(100 x10))
Ro (v) * 1b) [j kl.
(95)
0z
M
Derivatives of ((x) (1ub * x
We get:
1 M
(1p * xx1;2 = 2(1pb * x|) x
M Σ Re [sg( * x)(]ppii o]]
(96)
dx,
2
{11b * æ²)) Re [ast * sg(b x²] [
(97)
M
and:
8²
8²
(10 * xz)/²
= 201]==
(1 * xe))
(10 * æe)) + 2(1 * xx))
(l0 * x1),
(98)
0r,
"zz
0r
20

---

Published in Transactions on Machine Learning Research (02/2024)
2
Re [iat *gg((b * x)] [i] x Re [ast sg(( * x)] [k]]
(99)
M
+ M
1
(lv * r1) Σ Re
bb
ojbb [k- -1
bbli sg( * x 2 15 [k: 1]
(100)
1/b * x: [i]
16b * 2: [
2
Re [1bt * (gg((((a)))))
[i] x Re [ast sg( 1b * x]] [k]]
(101)
M2
1
1
sg(U )2 c * x
+ (lv= * x1> x Re
v [J. k] +0t
(102)
M
| x
where we have introduced the notation a a a Sobec [ū.1] ac* blj].
[i EM alj i]blū]c[k i]. Conveniently, note that
Derivatives of (xx) = (11b1 * a x 1/2 xx) We gel:
0 (1/1 w////22 * x²) ary
( Jo,
(1/2 * xxl)
x 4/2 * xt) + (115 xo| x dr. (th2 x)); >
(103)
Re [sg((h (1 =1] 2/2 * 2r + | * x z01i//bz2i =1]
(104)
M
M
Σ
|sg(1 * x - 5l002 * x11]] sg((21 * x - 5100 * x11]]
(105)
2M
i=1
M V/.2
1//2 * 20 [p]
(106)
|0] (sg(th) * xx) x 22 x) + vo * (sg(th
* xx) x 2/ x) b
(107)
2M
Vo * |1/1 * ++1]]. M
(108)
We break down the calculation of the second-order derivatives of & by first calculating:
M
(v2 * 1//1 * =1[s1)
Σ
Re Sg((12 l1l xi]]-
(109)
Ozx
1=1
M
Σ
|>>[1: i] sg(t Mis5 s]
(110)
+wi [k (3gg(h [WW - =1]
(111)
1
2 2 sg(gt) + 2 21 2 sg(U) xx) 2 2/ [.
(112)
We also have:
dax (sg(V1 * x)) x the xz)
2 x
1
t [k=
i}
sg(e) * 1]] s]
(113)
2/ x [
+ sg(U12
a].
(114)
(115)
so that:
0 w1 * (sg(t) * xx) x 2/2 * 1)] Oxy
Σ
M
:]
22
x
V2 +
[k
(116)
2
Jo x [i]
1 M
2 x x sg(1/) (2 * 2
Σ
L - 1)
w o]
(117)
2 1=1
1/ * x
21

---

Published in Transactions on Machine Learning Research (02/2024)
M
M i) sg(( 11
[d] 122 [1: a].
(118)
2/2 x x
422 * x x sg(
* x))
+
@ @
(119)
2
1/1 x2
1//1 x
+ +20
0
sg((
(120)
Therefore:
02
(101 * x
1
1/2 *
1/1
+ 1/1
x
& ms
(121)
0r
4M
2 +
1 * a
12 x x
sg(
1 *
sg 11 * x)2
@ vs
(122)
|
M * x
+ 201 @ sg(t) x 8 2/2
2011
X sg((1 1h * a 6 w
(123)
+ +201 sg(U) xx) 8 11 2 @ sg(
f1 R].
(124)
Summary for vpt = a and Σ = d For ((x) E C, we only focus on simplified expressions of:
M
OΦ
(1(())((o)))))) = Σ
(x)
(125)
1=1
M
(H)(2), E)r = Σ
x (z)o²
(126)
dzy
1=1 1
We get for (J)()))((x)) E)r:
4
M
)(x)= = 5TH (x)
Σ
Re (vb * xx) x 02(
(127)
M
M
((x) = 500 (x)
Re [bb * 16 x (1 x x xb sg( x)))) [1] x 2
(128)
M
i=1
((x) = 501(2) (x)
1
3x/ * 1 xr + 2 (sg( x 16 xr)
[i] c
(129)
4M2
((x) (OL(3)
1
1/2 (sg(v) * xx) Lb2 * x + Uh * (sg(v) xr) 4/22 * x))
(130)
4
+2xb | * x [t] x ??
(131)
And we get for (H²( E)e:
M Re (1b * 1/) [0] Σ>
2
x(x) = 511 (x)
0² [1]]
(132)
M
i=1
M
1
((x)= = (x)
Re
21 * (b(0)
Re (1) sg(b xx)]²
(133)
MM
M
(sg( * x²)² /// * xc|
]
+ (14 * xr/>
1/2
1/612
M [
(134)
1 *
M
=(x)= = (x)
6/ 1 * sg( 10 * x)) + +z2²
sg(bb (1b *
(135)
4M
22

---

Published in Transactions on Machine Learning Research (02/2024)
2 (sg( * x)
[t] x [²
(136)
M
6(x) (x)
2 2/1/1/ 12
W2 *
+ 4 * sg(t/1 x)) + 4401 1/ * sg((/1 * xx)
(137)
4M 1=1
W1
4/2 * x
x sg(11 * x)) 209
(pi* (
1/2 *
x sg((1 xx))
[d] x x
(138)
1//1 * x
B Relation to Maximum Likelihood Estimation
We investigate the links between statistical component separation methods and maximum likelihood estimation (MLE).
Estimation of zo- In a denoising context, the goal is to recover 2o given A naive application of MLE would be to maximize the likelihood p(yl zo). Assuming a Gaussian noise with covariance E, we have p(y (o) x exp (-2 11yly=lzl'' with lyy - zolls (y - =0)² E-1 (y zo); which leads to the trivial estimate do y. Since in this setting, naive MLE can only bring back to the noisy data, typically confuses MLE with the estimation of a maximum posteriori (MAP), which requires the definition of a prior distribution p(r) over the target signal (see e.g. Elad et al., 2023). With p(z) exp (-Ap(z)), still in the case of a Gaussian noise, a MAP estimate fn of zo reads:
2/ Ily - =1| + )
8o E arg min
(139)
This minimization problem aims to strike a balance between two constraints: one enforcing proximity the noisy data y, and the other imposing a prior constraint, typically in the form of regularity constraint over The parameter À controls the relative weight of these two constraints. A statistical component separation method can be remotely related to MLE in this picture. The estimate in defined in Eq. (1) can be interpreted as an MLE estimate by choosing p(z) = L(x) = Be-ples) [|0[ + e) du)l and picking
À + +00. Morcover, by initializing the optimization of ((x) with y, an implicit notion of proximity to y is implied.
Estimation of ((x)). When focusing on the estimation of (xx)) given y, which is the ultimate goal of a statistical component separation method, MLE methods require an explicit model connecting o(20) y. The approach taken in this paper share similarities with the sampling process of maximum entropy models defined in Bruna & Mallat (2019). For the purpose of this discussion, we model Zo realization of a macrocanonical model with density p(x) 2 exp (-0s Φ(z)), where (o(m)) is a vector of parameters fixed so that Eg<p(a) [d(x)] = @(xo). In this context, assuming a Gaussian noise with covariance Σ, the likelihood pyy ((xx)) can be written as:
ply o(zo)) & eep (-focen) o(x') -
1 11 2'11 dz'.
(140)
This likelihood remains generally intractable, which prevents further connections to the statistical component scparation method presented in this paper.
c Complementary Details on the Data
We give further details on the nature the data used for the experiments in Sect. 3 and 4.
The dust image This image corresponds to a simulated intensity map of the emission of dust grains the interstellar medium. It was taken from the work of Régaldo-Saint Blancard al. (2023), and we refer this paper for further details the way it generated.
23

---

Published in Transactions on Machine Learning Research (02/2024)
The LSS image This image was built from cosmological simulation taken from the Quijote suite (Villescusa-Navarro et al., 2020). This simulation describes the evolution of fluid of dark mat- ter in 1 (Gpe/h)a periodic box for a ACDM cosmology parameterized by (Am, Uz h, nay σa) = (0.3223, 0.04625, 0.7015, 0.9607, 0.8311) (high-resolution LH simulation). We take a snapshot at & = 0 of the dark malter density field at and project a 500 x 500 x 10 (Mpe/h)² slice along the thinnest. dimension. Our test image the logarithm of this projection.
The ImageNet image We have randomly picked a RGB image the ImageNet dataset (Deng et al., 2009) from the *digital clock" class. For the experiments of Sect. 3.2, to simplify, we have first preprocessed the image by averaging over the channel dimension.
D Definition and Computation of the WPH statistics
We use the same bump-steerable wavelets as in Regaldo-Saie Saint Blancard et al. (2021) with = and = 4, which leads a bank of are estimated follows:
28 wavelets. For a given image 2, the covariances introduced Sect. 3.2
(x) = (|x * thi
(x) = ((x * the 2hl (1x * 2 C (|x s * 2b
(141)
5 (|x the vt)
(142)
where (-) is spatial average operator. Moreoever, the numerical experiments of this paper, these coefficients are systematically normalized according the S11 coefficients of the noisy map y as follows:
(1) = (11(x)
(0) = 5:00 (x)
(143)
(1)
$i1 (y)
(x)
C (2)
(z) (7)
C9 (x) 5 ())
(144)
As it was reported in the related literature (Zhang & Mallat, 2021; Allys et al., 2020: Regaldo-Saint Blancard et al., 2021; Régaldo-Saint Blancard et al., 2023), this normalization acts preconditioning of the loss function L, and has then a direct impact on the optimum to one gets with a given optimizer. The quality of the results then a priori depends on this normalization, and it worth noticing than alternative approaches as in Delouis et al. (2022) could have been explored.
E Connection with Delouis et al. (2022)
We draw comnections between Algorithm 2 and the algorithm introduced in Delouis et al. (2022), which is formally described in Algorithm 3. These connections show that the two algorithms are very related. A further exploration of the pros and cons of each these two algorithms is left for future work.
The following lemma highlights the similarities between the loss function L defined in Eq. (1) and the one involved in Algorithm 3.
Lemma E. Introducing ax + c)² = E [lo(2 + e] [0(x +) that is the vector of variances of each component of Φ(x + E), we have:
((x) Eemp(ro) [110( +c) -(v)))?
(145)
= 10 + +)1l + ||E (0(x + e)] 2
(146)
24

---

Published in Transactions on Machine Learning Research (02/2024)
Algorithm Statistical Component Separation in Delouis et al. (2022)
Input: y, p(ea), Q, T, P
Initialize: y for i= 1. do
sample .Q p(eo)
m De 6(00 + <t/Q
Σ 1 (0(80 + Ek) - -m)?/// 1/2
B o(ăg)
for
j 1.... do
( = |1 + B - -v()) /012
2n f ONE_STEE STEP OPTMM ((00)
end for end for return 20
Proof.
[(x) = Bemp(ee)
|1 (x + <) ((v)
(147)
E |10( + e]]
+ ||( - 21E Re [6( + e) o(v)
(148)
E |1 + 11
+ ||o())11 Re (E [Ø(x + c)] o(v)
(149)
E || + +) )]1
||B [Ø(x )] + ||E [o(o + +)]1
(150)
+ ||o(v))1 Re
[ (16(=+))] + p(v)
(151)
= 1| + )) + |IE [d(x + c)] ((y) 2
(152)
Thanks to this lemma, now appears as the sum of two terms. The first one constrains the norm of the variance vector az(a + e) be minimal, while the second one constrains the mean vector [Ø(x + <)] to be close to d(y). In the light of this new expression of L, Algorithm aims to minimize the related following loss:
[(x) = (d(x + e)] -
0x +() | 2
(153)
which appears to be the second term of Eq. (146) normalized by the first term the same equation. However, instead of minimizing directly Σ, Algorithm makes the following approximations:
E [Ø(x + c) d(x) - B,
(154)
0사(x +c)
(155)
where B and a are repectively the bias and standard deviation terms explicited in Algorithm 2. Same as in Algorithm 2, Algorithm 3 adopts a stepwise approach, where and a are updated at each step i e {l...., P} using the f signal obtained from the previous step.
F Additional Results
We show Figs. F.1, F.2, and F.3 (respcctively, F.4, F.5, F.6, and F.7) equivalent results Figs. 2, 3, and 4, respectively, for the LSS (ImageNet) data. We also show in Fig. F.8 equivalent results Fig. for the ImageNet data and ConvNet-based representation.
25

---

Published in Transactions on Machine Learning Research (02/2024)
Pink
White
Blue "Crossss"
Noiey [18.77 dB
Noiey [19,2 an
Noiey 2000 dB
Noiey
Original
[17.4 dB
| eiB
Demonioon (oorns) 122.
Demoieet (ours)
Deaoionn (onus))
Denodooe
Denotedd (BMND) (17.6 d3)
Denoiedd (DM (DMD)) |21. dT))
Donosked (PM HDD 1:22 ATR)
100
200
100
200
100
200
100
200
Figure F.1: Same as Fig. 2 for the L.S image.
109
10P
Demoised Pink noiee (ours)
Demolsed Pink noise (BM3D)
crror 5T0
Demoised noise (ouns)
Rel. error
10
10-1
Damoised White noise (BM3D)
Demnieed Bhue noise (ouns)
10-2
Demoleed Blue noise (BM3D)
Higher better
Lower better
8
10
10a
30
eeror Ca
PSNR (aB)
Rel. arror
15-
10
10
109
10
10-1
10a
Noise level
Noiae level
Noise level
Figure F.2: Same as Fig. for the LSS image.
26

---

Published in Transactions on Machine Learning Research (02/2024)
10
Rel. ervor 59c
Noisy Crosens noise
Hel. ctror
Domoleed Crasses (ours)
10
Hligher
Lower is better
109
30
10"
25
FS PSNR (dB)
Itl. error
Rel. emot
15
10
10
10
2
Density of crosses
Density of crooses
Density of crosses
Figure F.3: Same as Fig. 4 for the LSS image.
Pink noise
White noise
Blue noise
"C Crossrss" noise
Noiey (6.2 dB)
Noiey 16.3 aTj
Noaey 1122 ABS
Nosey
Original
Decooieed (ears) 115.4 dtB)
Devasiond Jearss 118.2 dT.
Demonised Tonure 11:48 diB.
Decoised fours 117.4 dBB
50
150-
200
Detoked (BMAD) [i5. dB)
Denoiee (DM)))) HD) (22. 2 aD)
(DD len, i (TF)
Figure F.4: Same as Fig. for an ImageNet image.
27

---

Published in Transactions on Machine Learning Research (02/2024)
10a
10
Denoised Pink noiee (ours)
Denoised Pink nolee (BM3D)
10
Rel error 550
10"
Denoised Whitt (ours)
Rel. error
Denoised White voise (BMSD)
10
10-1
Denoised Blue voise (ours)
Denoised Blue noise (BM3D)
10
Highor botter
10aa
Lower is better
10P
30
orvor Co.
PSNR (dD)
Rel. croor
10-1
10-1
10
Noise Jevell
Noian level
Noise level
Figure F.5: Same as Fig. 3 for the ImageNet image shown in Fig. F.4.
10a
10"
109
Rel. ceror Sto
Nolgy Crosses
[tel. esror
Denoised
(ours)
10-1
Higher betler
Lonee is better
10
101
25
20
PSVR dB
Rel. error
109
orvor
10"
15
10
10-1
1
10
5
3
Density of
Density of
Density of
Figure F.6: Same as Fig. 4 for the ImageNet image.
28

---

Published in Transactions on Machine Learning Research (02/2024)
gu
Demoliani Pink molan (surs)
Denobred
Pank noise (DMGID)
10-1
Higher better
Lower better
30
PSNV (.88)
25-
5
orme
20-
10
10-1
5
10"1
10a
Notar level
Noian lend
Noiou Jevell
gu
10-1
Demonsoee White moiee (nurs)
10
orve
Denoised White zoise (BMBD)
10
10
Higher is betberr
Lenerr Twiteer
109
11P
PSWH (U1]]
200
erno
10-1
10
10
16
Notse level
Notse lovel
Noise Jevell
Bel gt
10-1
erne 200
Denoioni Blue (nurs)
Demoivee Rlos
10-2
10
10
Bligher is better
Lower better
100
108
05
30-
PSNT (18)
25
Bat. unee
10
20
10
Figure F.7: Same as Fig. for randomly selocted ImageNet images. We show the mean and standard deviation of the metrics computed across random batches of test ImageNet images.
29

---

Published in Transactions on Machine Learning Research (02/2024)
Pink noise
White noise
Blue noise
CCrosses" noise
Noiey [17. di3)
Notey 117.0
Nosey 117
Nojey
Original
Deturised (otzz)) [19. dE)
Deuoiso? (BMRD) 120.5 aT3) Demoieed (ouzs) 18.55 an)
Donoised (cuns) dP)
Doueieel (ouus) 139 dI
50
Deuoinst (BMSD) [22.8 dt8)
Donolsoo) (BMGD) [a ATB)
100
200
100
200
100
200
1001
2000
Figure F.8: Same as Fig. 2 for an example of ImageNet image and a ConvNet--ased representation (see Scct. 3.3).
30


Reviews:
[{"review_id": "iBQ710ccO0", "reviewer": "Reviewer_bxtv", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\nOverall, I think this paper is a good fit for TMLR. It analyzes an interesting recently-developed algorithm that could have a significant impact in astronomy and other scientific domains where data is hard to come by. While the performance of  extensions proposed in the paper is somewhat disappointing, I think this represents an interesting negative result that will inform the ML community.\n\nThe paper is relatively easy to follow.\n\n## Weaknesses\nThe paper might benefit from a high-level comparison with maximum likelihood estimation and other conventional estimation methods.\nThe performance of the proposed methods is underwhelming.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper analyzes and extends the \"statistical component separation\" algorithm developed by Regaldo-Saint Blancard et al. and others. The algorithm assumes one has a single noisy observation, y=x+\\epsilon, of a signal contaminated with noise and a method for sampling from the distribution of the noise. It then seeks a solution \\hat{x} that minimizes that minimizes the expected l2 difference between \\phi(x+noise) and \\phi(y), where \\phi is user-defined function and y is the noisy observation. Previous works have empirically demonstrated the resulting \\hat{x} and \\phi(\\hat{x}) are effective estimates of x and \\phi(x).\n\nThis paper...\n-Provides several examples of operators \\phi() where the above optimization problem has interesting analytic solutions.\n-Provides empirical evidence that with \\phi() defined as the covariances of wavelet coefficients of the images, the proposed method effectively denoises images (particular atmospheric images of dust) and can handle complex noise distributions (crosses randomly added to locations on the image).\n-It demonstrates that (somewhat counter-intuitively) defining \\phi() as the VGG feature mapping and matching the signal there results in poor overall denoising performance. I think this is an interesting negative result.\n-It develops (and partially analyzes) a diffusion-inspired denoising method which gradually denoises the signal in stages. This approach marginally improves performance over the single-shot approach.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "The optimizer was only run for 30 iterations. Please comment on the effect the number of iterations has on the performance of the algorithm. Is early stopping (à la deep image prior) responsible for some of the denoising performance?", "broader_impact_concerns": "No concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "vW13DrqLkA", "reviewer": "Reviewer_ShsH", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths: \n* The analysis aspects of the paper are clear, and a contribution to the community.\n* The overall language of the paper is very clear.\n* The experiments are fairly designed; both strengths and weaknesses of the proposed algorithm are portrayed.\n\nWeaknesses:\n1. Some language in this work seems to forget its position in the broader context of signal separation / denoising / optimization at large. Details are in 'requested changes' section.\n2. One specific example (with only one realization) of structured noise is provided (randomly placed crosses), i.e., the rightmost column of Fig 2. Everything else in the paper is _denoising_, *not* _signal component separation_. The paper presents itself as a \"component separation\" algorithm that works with \"mixtures\", but that's totally misleading when the only undesired component is Gaussian noise.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This work analyzes a recently proposed approach to extract signals from noisy mixtures. The signal to be extracted is specified in terms of a set of desired statistics (computed by representation transform $\\phi$), presented to the algorithm via a (possibly small) dataset of template signals. This paper analyzes the cost function of this problem for linear $\\phi$ and quadratic $\\phi$, leading up to analysis of a power-spectrum based summary statistic that has been shown to be effective in practice. The authors provide a statistical explanation of this problem's solution's relationship with the desired signal component (the ground truth used to generate the mixture signal under the assumed model).\n\nNumerical experiments quantitatively compare the algorithm in a denoising setting, and smaller experiments are shown for a \"structured noise\" (aka component separation with randomly placed crosses) setting.\n\nFinally, a diffusive extension of the algorithm is derived and also compared with BM3D in a white noise denoising setting.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Critical to securing recommendation:\n1. A recurring theme is that this work does an insufficient job presenting its position in the context of optimization, signal separation, & denoising at large. This is illustrated by your conclusion following the statement of Proposition 2.1, that a nonlinear $\\phi$ must be used to perform denoising, thereby discarding (literally) hundreds of years of denoising theory, not to mention highly successful, modern signal separation methods that use linear $\\phi$[1-3]. So I find this to be a fallacious (and a bit lazy?) justification for your choice to use nonlinear $\\phi$. The authors/ the paper/ the readers would be much better served if the authors owned up to the decision from a holistic perspective. For a journal paper, isn't it reasonable to ask: why do you choose to use the regularization function $R(x,y)=0$ in Cost Function (1)? How is the proposed approach to incorporating prior information better than the ``Bayesian'' way (via a choice of $R(x,y)\\neq0$)?\n\nI understand that Prop 2.1-2 are in the paper to lead up to Prop 2.3 (which is nicely done!), and I do not expect you to seriously numerically compare these. However they are brushing up against a whole world of techniques that are just barely referenced. In other words, no one would ever just \"shrink the signal to zero\" if the noise is just over said threshold. There are iterative methods and many continuous distributions of thresholding operators between soft thresholding and identity (see [4] and Fig 1 in [5]).\n\nProposition 2.1: Why not allow $E[\\epsilon]\\neq0$ show that you get a biased result? For component separation (as opposed to Gaussian denoising), of course $E[\\epsilon]$ would be nonzero.\n\n2. In my opinion, you should either (1) change all mention of \"component separation\" and \"mixtures\" to focus on what this is: denoising. You have a nice denoising analysis for different colors of Gaussian noise but only one single realization of the crosses noise. At the very least there should be some analysis on a series of corruption levels as in Figs 3-5. Or, (2), provide additive mixture examples where $\\epsilon$ is not just Gaussian noise. I think this paper would be so much stronger if it stood against other component separation approaches (MCA is the main one that comes to mind [1-3]). How, and in what sense, must $\\epsilon$ be distinct from $x_0$ for successful extraction? Your statistical analysis may be able to address this critical question that is also outstanding in MCA (see \"spark\" in [3]).\n\nOther suggestions:\n* Isn't $\\mathbb{K}=\\mathbb{C}$ ? Why don't you just use $\\mathbb{C}$?\n* A true novelty of statistical component separation compared to so much convex optimization literature is in the way a dataset (e.g. a way to sample \\epsilon) is used to shape the optimization landscape... and a novelty compared to other ML approaches is the possibility of using relatively small datasets. I think it will be more clear to the majority of readers (who just glance at the first page or two) if you state the empirical cost function ($\\hat{L}$ in Alg 1) near the beginning with (1)\n* The work performed in Lemma D is perhaps worth emphasizing \\textit{more}, e.g. providing a framework that joins the approaches.\n* Instead of picking a \"random\" image from imagenet (which I find hard to believe anyway...), why not perform the experiment on N random draws of imagenet and thus providing an interesting statistically significant experiment?\n* There are a number of places where the quotation marks do not match on either side of the word. Typically in LaTeX it looks good to have `` on the left and '' on the right (two apostrophes, not one right-quotation mark \")\n\n[1] Goehle, Geoff, et al. \"Approximate extraction of late-time returns via morphological component analysis.\" The Journal of the Acoustical Society of America 153.5 (2023): 2838-2838.\n\n[2] Parekh, Ankit, et al. \"Multichannel sleep spindle detection using sparse low-rank optimization.\" Journal of neuroscience methods 288 (2017): 1-16.\n\n[3] Starck, J-L., Michael Elad, and David L. Donoho. \"Image decomposition via the combination of sparse representations and a variational approach.\" IEEE transactions on image processing 14.10 (2005): 1570-1582.\n\n[4] Selesnick, Ivan. \"Penalty and shrinkage functions for sparse signal processing.\" Connexions 11.22 (2012): 23.\n\n[5] A. H. Al-Shabili, Y. Feng and I. Selesnick, \"Sharpening Sparse Regularizers via Smoothing,\" in IEEE Open Journal of Signal Processing, vol. 2, pp. 396-409, 2021, doi: 10.1109/OJSP.2021.3104497.", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "hwxY7NZ7Ya", "reviewer": "Reviewer_ZTZn", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strength**\n1. The paper is well-written and well-organized. It is easy to follow.\n2. The paper seems to be technically correct. Candid numerical experiments and comparison are provided to illustrate the performance of the proposed method.\n\n**Weakness**\n1. The significance of the paper is questionable. The proposed method has only been studied by the same set of authors in recent years, and I wonder what is the difference/relation between the proposed method and maximum-likelihood?\n2. I understand the method is not supposed to be a denoising algorithm, so one does not expect it to beat BM3D in denoising. However, I wonder what \"is\" the strong suit of the proposed method? In particular, after obtaining the WPH and VGG statistics (with a slightly better relative error), what can I use it for?\n3. The proposed diffusive algorithm only slightly improves the performance, but might introduce much more computational cost.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper aims to solve the problem of “statistical component separation” that focuses on recovering a predefined set of statistical descriptors of a target signal from a noisy mixture. Assuming prior knowledge of the noise distribution, the authors analyze a method devised to match the statistical discriptor of the solution. For simple discriptors (such as linear, quadratic, and power spectrum representations), analytic solutions are provided after simple calculations. The authors also empirically studied 1) wavelet-based descriptors and 2) ConvNet-based descriptors, and showcase the performance of their method in recovering (from noisy measurement) of the underlying noise-free discriptors. The authors also extend the method by introducing a diffusive stepwise algorithm, with slightly improved denoising performance.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "No", "audience_explanation": "", "requested_changes": "1. In proposition 2.3, is $A$ injective? Isn't $A$ the convolution with bandpass filters, and hence not injective?", "broader_impact_concerns": "no", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: Overall, all three reviewers agreed that the paper had interesting contributions, and hence were positive about publishing it. After looking at the paper myself, I agree that this is an interesting idea.

For minor revisions, please do a final check for typos, and fix a few capitalization issues in the bibliography (e.g., capitalize "herschel" since it is a proper noun).
Decision: Accept (poster)