# Learning from Incomplete Features by Simultaneous Training of Neural
Networks and Sparse Coding

Cesar F. Caiafa12, Ziyao Wang3, Jordi Solé-Casals4, and Qibin Zhao1  Tensor
Learning Team, Center for Advanced Intelligence Project, RIKEN, JAPAN1
Instituto Argentino de Radioastronomía, CONICET-CCT La Plata / CIC-PBA / UNLP,
ARGENTINA2
School of Automation, Southeast University, CHINA23
University of Vic - Central University of Catalonia, SPAIN4
Corresponding author: C. F. Caiafa (Email<EMAIL_ADDRESS>
###### Abstract

In this paper, the problem of training a classifier on a dataset with
incomplete features is addressed. We assume that different subsets of features
(random or structured) are available at each data instance. This situation
typically occurs in the applications when not all the features are collected
for every data sample. A new supervised learning method is developed to train
a general classifier, such as a logistic regression or a deep neural network,
using only a subset of features per sample, while assuming sparse
representations of data vectors on an unknown dictionary. Sufficient
conditions are identified, such that, if it is possible to train a classifier
on incomplete observations so that their reconstructions are well separated by
a hyperplane, then the same classifier also correctly separates the original
(unobserved) data samples. Extensive simulation results on synthetic and well-
known datasets are presented that validate our theoretical findings and
demonstrate the effectiveness of the proposed method compared to traditional
data imputation approaches and one state-of-the-art algorithm.

##  I Introduction

Learning methods from limited or imperfect data has attracted great attention
in the literature recently. Datasets with limited, weak, noisy labels or
incomplete features represent an important and still open problem. In this
paper, we address the problem of training a classifier on a dataset with
incomplete features, which arises in many machine learning applications where
sometimes the measurements are incomplete, noisy or affected by artifacts.
Examples of this situation include: recommendation systems built upon the
information gathered by different users where not all the users have fully
completed their forms; medical datasets where typically not all tests can be
performed on every patient; or a self-driving vehicle or robot where objects
in the view field can be partially occluded.

Handling correctly the incomplete-features problem is a classical challenge in
machine learning. Skipping missing features by setting them to zero values
damages the classification accuracy [1]. Most previous studies addressed this
problem by using an imputation approach, which consists of performing data
completion followed by training the classifier with those reconstructions
(referred here as the sequential method). However, this strategy cannot ensure
the statistical consistence of the classifier, as data completion is usually
fully unsupervised or label information is partially or inefficiently
exploited.

In this work, a new supervised learning method is developed to train a general
classifier, such as a logistic regression or a deep neural network, using only
a subset of features per sample, while assuming sparse representations of data
vectors on an unknown dictionary. The proposed method simultaneously learns
the classifier, the dictionary and the corresponding sparse representation of
each input data sample. In this way, we combine the approximation power and
simplicity of sparse coding with the extraordinary ability of neural networks
(NNs) to model complex decision functions (classifiers) with the goal to
successfully train a classifier based on incomplete features.

We analyze the limitations of the sequential approach (section I-C), i.e.
imputation followed by training, and introduce the simultaneous classification
and coding approach in section II. Our method consists of incorporating a
sparse data representation model into a single cost function that is optimized
for training the classifier and, at the same time, finding the best
representation of the observed data. A learning algorithm is presented in
section II-A to train a classifier on incomplete features and sufficient
conditions under which such a classifier performs as good as the ideal
classifier, i.e. the one that can be obtained from complete observations, is
identified (section II-B). Extensive experimental results are presented in
section III, using synthetic and well known benchmark datasets that validate
our theoretical findings and illustrate the effectiveness of the proposed
method.

###  I-A Related work

Practical sequential methods based on statistical imputation, such as
computing the “mean”, “regression” and “multiple” imputation techniques are
common practice [2]. Remarkably, it was shown that imputing with a constant,
e.g. the mean, is Bayes-risk consistent only when missing features are not
informative [3]. More elaborated completion methods were also explored, such
as $K$-nearest neighbor estimators, multilayer or recurrent NNs and others,
see [4] and references therein.

However, sequential methods do not fully exploit label information. Data
labels can provide valuable information about missing features that could
potentially improve the classifier learning process. Recent advances on
probabilistic generative models have allowed for a formulation of supervised
learning with incomplete features as a statistical inference problem arriving
at algorithms that significantly outperformed sequential methods. In the
seminal work [5], a framework for maximum likelihood density estimation based
on mixtures models was proposed and successfully applied to small incomplete-
features problems. In particular, a Gaussian Mixture Model (GMM) was fitted to
incomplete-data through an Expectation-Maximization (EM) algorithm. Building
upon this generative model strategy, some approaches have considered
integrating out the missing values based on a simple logistic regression
function [6, 7]. Other versions of this approach proposed an explicit
simultaneous learning of the model and the decision function [8, 9]. While
probabilistic generative models provided a nice and elegant approach to the
incomplete-features problem showing good results on small datasets, they are
not suitable for many modern machine learning applications because: (1)
despite some acceleration techniques were explored, e.g. [10, 11], those
algorithms are computationally expensive becoming prohibitive for moderate to
large datasets; (2) GMM is impractical for modeling high-dimensional datasets
because the number of parameters to achieve good approximations becomes
unmanageable; and (3) they do not consider complex classification functions as
the ones provided by deep NN architectures.

Recently, some approaches based on the low-rank property of the features data
matrix were investigated and algorithms for data completion were proposed
incorporating the label information [12, 13, 14]. Since the rank estimation of
a matrix is a computationally expensive task, usually based on the Singular
Value Decomposition (SVD), the obtained algorithms are prohibitive to solve
modern machine learning problems with large datasets. Additionally, as in the
case of the probabilistic generative models, none of these methods considered
complex classification functions. To overcome this drawback, more recently, a
framework based on various NN architectures such as autoencoders, multilayer
perceptrons and Radial Basis Function Networks (RBFNs), was proposed for
handling missing input data by setting a probabilistic model, e.g. a GMM, for
every missing feature, which is trained together with the NN weights [15].
This method combined the great capability of NNs to approximate complex
decision functions with the nice formulation of the GMM to model missing data.
However, it inherited the drawbacks of GMMs, i.e. they are not well suited to
higher-dimensional datasets.

On the other hand, during the last few years in the signal processing
community, there has been a rapid development of theory and algorithms for
sparse coding approximations which, by exploiting the redundancy of natural
signals, are able to provide simple and accurate models of complex data
distributions, see [16, 17, 18, 19, 20] and references therein. Sparse coding
is nearly ubiquitous in Nature, for example, it is found in the way that
neurons encode sensory information [21, 22]. Sparse representations of data
showed to be useful also in classification problems. In [23], a Linear
Discriminant Analysis (LDA) classifier was trained on corrupted data providing
a robust classification method. In [24], algorithms for learning
discriminative sparse models, instead of purely reconstructive ones, were
proposed based on simple linear and bilinear classifiers. Similar methods were
also studied by either using class-specific dictionaries [25, 26] or using a
single one for all classes [27]. However, these proposed methods neither were
applied to the incomplete-features problem nor considered deep NN classifiers.

###  I-B Problem formulation

We assume a supervised learning scenario with vector samples and labels
$\\{{\mathbf{x}}_{i},y_{i}\\}$, $i=1,2,\dots,I$,
${\mathbf{x}}_{i}\in\mathds{R}^{N}$ and $y_{i}\in\\{0,1,\dots,C-1\\}$ ($C$
classes). However, we are constrained to observe only subsets of features and
their labels: $\\{{\mathbf{x}}^{o}_{i},y_{i}\\}$,
${\mathbf{x}}^{o}_{i}\in\mathds{R}^{M_{i}}$ with $M_{i}<N$. Unobserved
(missing) features are denoted by
${\mathbf{x}}^{m}_{i}\in\mathds{R}^{N-M_{i}}$. We consider arbitrary patterns
of missing features, which are allowed to be different for each data instance
$i$. The set of indices of missing features at sample $i$ is denoted by
${\cal{M}}_{i}$, i.e. ${\mathbf{x}}^{m}_{i}={\mathbf{x}}_{i}({\cal{M}}_{i})$
and
${\mathbf{x}}^{o}_{i}={\mathbf{x}}_{i}\left({\overline{\mathcal{M}}}_{i}\right)$.

We define the set of all $K$-sparse vectors
$\Sigma_{K}^{P}=\\{{\mathbf{s}}\in\mathds{R}^{P}\text{ s.t.
}\|{\mathbf{s}}\|_{0}\leq K\\}$ (containing at most $K$ non-zero entries) and
assume that data vectors ${\mathbf{x}}_{i}$ admit $K$-sparse representations
over an unknown dictionary $\mathbf{D}\in\mathds{R}^{N\times P}$ ($P\geq N$):

${\mathbf{x}}_{i}=\mathbf{D}{\mathbf{s}}_{i},\mbox{ with
}{\mathbf{s}}_{i}\in\Sigma_{K}^{P}.$ (1)

The columns of a dictionary are called “atoms” because every data vector can
be written as a linear combination of at most $K$ elementary components.
Sometimes dictionaries are orthogonal such as the ones derived from the
Discrete Cosine or Wavelet [19] transforms. However, overcomplete ($P\geq N$)
nonorthogonal dictionaries have demonstrated to play an important role in
image processing tasks such as denoising, inpainting, etc [17, 28].

By partitioning $\mathbf{D}$ according to the pattern of missing features at
sample $i$, we obtain
$\mathbf{D}^{o}_{i}=\mathbf{D}\left({\overline{\mathcal{M}}}_{i},:\right)\in\mathds{R}^{M_{i}\times
P}$ and
$\mathbf{D}^{m}_{i}=\mathbf{D}({\cal{M}}_{i},:)\in\mathds{R}^{{(N-M_{i})}\times
P}$, which according to equation (1) implies:

${\mathbf{x}}^{o}_{i}=\mathbf{D}^{o}_{i}{\mathbf{s}}_{i},\mbox{ and
}{\mathbf{x}}^{m}_{i}=\mathbf{D}^{m}_{i}{\mathbf{s}}_{i}.$ (2)

Let us assume that a perfect classifier, e.g. a logistic regression or deep
NN, that assigns probability $p_{\Theta}(\hat{y}|{\mathbf{x}})$ to predicted
label $\hat{y}$ given data ${\mathbf{x}}$ can be trained on the complete
dataset $\\{{\mathbf{x}}_{i},y_{i}\\}$, such that, in a two-classes scenario
($C=2$), $p_{\Theta}(\hat{y}=y_{i}|{\mathbf{x}}_{i})>p_{\Theta}(\hat{y}\neq
y_{i}|{\mathbf{x}}_{i})$, $\forall i=1,2,\dots,I$, where $\Theta$ is the set
of trained parameters. Our goal is to develop a method to obtain an estimate
$\hat{\Theta}$ of parameters using only the incomplete dataset
$\\{{\mathbf{x}}^{o}_{i},y_{i}\\}$ and to identify conditions under which such
a classifier is compatible with the ideal one.

###  I-C Why training after imputation is difficult?

If the $K$-sparse representations of the observations ${\mathbf{x}}^{o}_{i}$
were unique, then ${\mathbf{x}}_{i}$ can be perfectly reconstructed from the
incomplete observations and the classifier can be successfully trained using
these reconstructions. In the particular case where the dictionary is known in
advance, there exist conditions on the sampling patterns based on the
coherence, spark or RIP (Restricted Isometry Property) of matrix
$\mathbf{D}^{o}_{i}$ that can guarantee uniqueness [20]. However, these
conditions are difficult to meet in practice and determining RIP/Spark
properties are NP-hard in general [29]. Moreover, in the general case where
the dictionary $\mathbf{D}$ is unknown and needs to be learned from data, it
is even more difficult to obtain well separated reconstructions which
certainly leads to suboptimal or wrong classifiers.

Next, we provide some intuition about the limitation of the sequential
approach through a toy example. Let us consider the classification of hand-
written digit images belonging to two classes: “3s” and “8s” and assume that
they admit $2$-sparse representations over a dictionary. Fig. 1 (a-b) shows
the representations of two example vectors ${\mathbf{x}}_{i}$ and
${\mathbf{x}}_{j}$ belonging to classes “3” and “8”, respectively. If only the
right halves of the images are observed and no label information is provided,
we are clearly faced with a problem because our observed samples from two
different classes are identical, i.e.
${\mathbf{x}}^{o}_{i}={\mathbf{x}}^{o}_{j}$. It is obvious that at least two
possible $2$-sparse representations for the observed data exist as illustrated
in Fig. 1(c). When the sparse solution is not unique, we may end up
reconstructing wrong vectors that could not be even well separated as
illustrated in Fig. 1 (d-e). In general, sequential methods using only the
information of observed features are prone to fail because the non-uniqueness
of solutions can make the training of a good classifier an impossible task.
However, we could solve this problem by incorporating the labelling
information from the very beginning as it is proposed in the following
section.

##  II Simultaneous learning and coding approach

We propose to train the classifier and find the proper representation, not
only as sparse as possible but also providing the best separation of classes.
We want to combine the training of the classifier together with the learning
of a dictionary and optimal sparse representations such that the reconstructed
data vectors are compatible with observations and well separated. To do that
we propose to minimize the following global cost function:

$\displaystyle J(\Theta,\mathbf{D},{\mathbf{s}}_{i})=$
$\displaystyle\textstyle\underbrace{\frac{1}{I}\sum_{i=1}^{I}\big{\\{}J_{0}(\Theta,\hat{{\mathbf{x}}}_{i},y_{i})+\lambda_{1}J_{1}(\mathbf{D},{\mathbf{s}}_{i})\big{\\}}}_{F(\Theta,\mathbf{D},{\mathbf{s}}_{i})}+\underbrace{\frac{1}{I}\sum_{i=1}^{I}\big{\\{}\lambda_{2}J_{2}({\mathbf{s}}_{i})\big{\\}}}_{G({\mathbf{s}}_{i})},$
(3)

with respect to $\Theta$, $\mathbf{D}$ and ${\mathbf{s}}_{i}$
($i=1,2,\dots,I$), where $\Theta$ contains the classifier parameters, i.e. the
vector of weights in a deep NN classifier architecture;
$\mathbf{D}\in\mathds{R}^{N\times P}$ ($P\geq N$) is a dictionary and
${\mathbf{s}}_{i}\in\Sigma_{K}^{P}$ are the representation coefficients such
that the reconstructed data vectors are
$\hat{{\mathbf{x}}}_{i}=\mathbf{D}{\mathbf{s}}_{i}$.

$J_{0}(\Theta,\hat{{\mathbf{x}}}_{i},y_{i})$ is a measure of the
classification error for the reconstructed sample vector
$\hat{{\mathbf{x}}}_{i}$. Typically, we use the crossentropy measure, i.e.
$J_{0}(\Theta,\hat{{\mathbf{x}}}_{i},y_{i})=-\log[p_{\Theta}(y_{i}|\hat{{\mathbf{x}}}_{i})]$,
where $p_{\Theta}(y_{i}|\hat{{\mathbf{x}}}_{i})$ is the probability assigned
by the classifier to sample $\hat{{\mathbf{x}}}_{i}$ as belonging to class
$y_{i}$. $J_{1}(\mathbf{D},{\mathbf{s}}_{i})$ is a measure of the
approximation error of the reconstruction when it is restricted to observed
features, which is defined as follows:
$J_{1}(\mathbf{D},{\mathbf{s}}_{i})=\frac{M_{i}}{N}\|{\mathbf{m}}_{i}\odot({\mathbf{x}}_{i}-\mathbf{D}{\mathbf{s}}_{i})\|^{2}$,
where $\odot$ stands for the entry-wise product,
${\mathbf{m}}_{i}\in\mathds{R}^{N}$ is the observation mask for sample $i$,
i.e. $m_{i}(n)=0$ ($1$) if data entry ${\mathbf{x}}_{i}(n)$ is missing
(available); and $J_{2}({\mathbf{s}}_{i})=\frac{1}{N}\|{\mathbf{s}}_{i}\|_{1}$
is proportional to the $\ell_{1}$-norm whose minimization promotes the
sparsity of the representation since $\ell_{1}$-norm is a convenient proxy for
$\ell_{0}$-norm [30]. Finally, the hyper-parameters $\lambda_{1}$ and
$\lambda_{2}$ allow us to give more or less importance to the representation
accuracy and its sparsity, with respect to the classification error.
Intuitively, minimizing equation (II) favors solutions that not only have
sparse representations compatible with observed features, but also providing
reconstructions that are best separated in the given classes.

Figure 1: Toy example: (a) 4 out $P$ dictionary elements $\mathbf{d}_{i}$
(atoms). (b) Digits “3” and “8” can be represented by combining only two atoms
in the dictionary ($2$-sparse representations). (c) A left-half occluded digit
“3” or “8” admits more than one $2$-sparse representation (sum of
$\mathbf{d}^{o}_{1}$ and $\mathbf{d}^{o}_{2}$, or $\mathbf{d}^{o}_{3}$ and
$\mathbf{d}^{o}_{4}$). (d) Linearly separable samples from two classes (A and
B) having two features: $\mathbf{x}=[x_{1},x_{2}]\in\mathds{R}^{2}$ where
incomplete observations are taken by observing only one feature. Note that
$\mathbf{x}_{A}$ and $\mathbf{x}_{B}$ belong to different classes but their
observations are identical. (e) Without using label information, the
sequential method could lead to wrong reconstructions of data vectors, i.e.
$\hat{\mathbf{x}}_{A}\neq\mathbf{x}_{A}$ and
$\hat{\mathbf{x}}_{B}\neq\mathbf{x}_{B}$ making the set of reconstructed
vectors not linearly separable.

###  II-A A sparsity-promoting sub-gradient optimization algorithm

To minimize the cost function in equation (II) we propose to alternate between
the optimization over ${\mathbf{s}}_{i}$ ($i=1,2,\dots,I$) and
$\\{\Theta,\mathbf{D}\\}$ using the training dataset (incomplete).

For fixed $\\{\Theta,\mathbf{D}\\}$, the optimization with respect to
${\mathbf{s}}_{i}$ is a non-smooth separable minimization sub-problem, which
was extensively studied in the literature [31, 32]. In this sub-problem, the
objective function is written as the sum of
$F(\Theta,\mathbf{D},{\mathbf{s}}_{i})$ and a non-smooth separable function
$G({\mathbf{s}}_{i})$, for which highly specialized, efficient and provable
convergent solvers, namely the Coordinate Gradient Descent (CGD), already
exists. However, the following key differences in our setting makes it not
suitable for the CGD approach: first, our function
$F(\Theta,\mathbf{D},{\mathbf{s}}_{i})$ involves evaluation of a multi-layer
NN classifier, which can be non-smooth due to involved activation functions
like ReLU or others; second, and more importantly, the computation of its
second derivatives (Hessian) becomes prohibitive. Therefore, we choose a
simpler and standard first order (stochastic sub-gradient based) search of
local minima with back-propagation. We take the strategy similar to the
heuristics used in [33]. To update ${\mathbf{s}}_{i}$, we need to subtract
$\sigma_{{\mathbf{s}}}\frac{\partial J}{\partial{\mathbf{s}}_{i}}(j)$ from
each coordinate $j$ provided that we do not cross zero in the process in order
to avoid escaping from a region where $G({\mathbf{s}}_{i})$ is differentiable.
In such a case, we let the new value of ${\mathbf{s}}_{i}(j)$ be exactly zero.
More specifically, we define
$\Delta_{i}(j)=-\sigma_{{\mathbf{s}}}\frac{\partial
J}{\partial{\mathbf{s}}_{i}}(j)$ and, if
${\mathbf{s}}_{i}(j)[{\mathbf{s}}_{i}(j)+\Delta_{i}(j)]<0$ (zero crossing
condition), we re-define $\Delta_{i}(j)=-{\mathbf{s}}_{i}(j)$; finally we
update ${\mathbf{s}}_{i}\leftarrow{\mathbf{s}}_{i}+\Delta_{i}$. It is noted
that, once a coefficient ${\mathbf{s}}_{i}(j)$ reaches zero at a coordinate
$j$, it becomes fixed, in other words, sparsity of solution ${\mathbf{s}}_{i}$
is monotonically increasing with iterations.

When ${\mathbf{s}}_{i}$ is fixed, our problem is reduced to minimize
$F(\Theta,\mathbf{D},{\mathbf{s}}_{i})$ with respect to $\Theta$ and
$\mathbf{D}$, which is easily done by standard first order (stochastic
gradient based) search of local minima. The algorithm proposed for the
training phase is presented as Algorithm 1.

In addition, for the testing phase, if the test dataset is incomplete, we need
to find first the sparsest representation for the given observations, compute
the reconstructions $\hat{{\mathbf{x}}}_{i}=\mathbf{D}{\mathbf{s}}_{i}$ and
then apply the previously learned classifier to them as presented in Supp.
material, Algorithm 2.

Algorithm 1 : Simultaneous classification and coding

0: $\\{{\mathbf{x}}^{o}_{i},y_{i}\\}$, $i=1,2,\dots,I$, hyper-parameters
$\lambda_{1}$ and $\lambda_{2}$, $N_{iter}$ and update rates
$\sigma_{\Theta}$, $\sigma_{\mathbf{D}}$ and $\sigma_{{\mathbf{s}}}$

0: Weights $\Theta$ and reconstructions
$\hat{{\mathbf{x}}}_{i}=\mathbf{D}{\mathbf{s}}_{i},\forall i$

1: Randomly initialize $\Theta,\mathbf{D},{\mathbf{s}}_{i},\forall i$

2: for $n\leq N_{iter}$ do

3: Fix ${\mathbf{s}}_{i}$, update $\Theta$ and $\mathbf{D}$:

4: $\Theta=\Theta-\sigma_{\Theta}\frac{\partial J}{\partial\Theta}$

5: $\mathbf{D}=\mathbf{D}-\sigma_{\mathbf{D}}\frac{\partial
J}{\partial\mathbf{D}}$

6: Normalize columns of matrix $\mathbf{D}$

7: Fix $\Theta$ and $\mathbf{D}$, update ${\mathbf{s}}_{i}$, $\forall i$:

8: $\Delta_{i}=-\sigma_{{\mathbf{s}}}\frac{\partial
J}{\partial{\mathbf{s}}_{i}}$, $\forall i$

9: if ${\mathbf{s}}_{i}(j)[{\mathbf{s}}_{i}(j)+\Delta_{i}(j)]<0$ then

10: $\Delta_{i}(j)=-{\mathbf{s}}_{i}(j),\forall(i,j)$;

11: end if

12: ${\mathbf{s}}_{i}={\mathbf{s}}_{i}+\Delta_{i}$, $\forall i$

13: end for

14: return
$\Theta,\mathbf{D},{\mathbf{s}}_{i},\hat{{\mathbf{x}}}_{i}=\mathbf{D}{\mathbf{s}}_{i},\forall
i$

###  II-B Theoretical analysis

Here, we investigate about conditions under which a perfect classifier of the
complete data can be obtained from incomplete data samples.

####  II-B1 Logistic regression

Let us first consider a logistic regression classifier [34] where the set of
parameters $\Theta=\\{{\mathbf{w}},b\\}$ are a vector
${\mathbf{w}}\in\mathds{R}^{N}$ and a scalar $b$ (bias). A perfect classifier
exists if there is a hyperplane that separates both classes, i.e., for each
data vector ${\mathbf{x}}_{i}$:
$f({\mathbf{x}}_{i})=\langle{\mathbf{w}},{\mathbf{x}}_{i}\rangle+b>0$ if
$y_{i}=1$, and $f({\mathbf{x}}_{i})\leq 0$ if $y_{i}=0$. We consider data
samples admitting a $K$-sparse representations
${\mathbf{x}}_{i}=\mathbf{D}{\mathbf{s}}_{i}$ with dictionary
$\mathbf{D}\in\mathds{R}^{N\times P}$ having unit-norm columns. We also assume
an arbitrary pattern of missing features $\mathcal{M}_{i}$ such that, data
samples and dictionary are partitioned as
$\\{{\mathbf{x}}^{m}_{i},{\mathbf{x}}^{o}_{i}\\}$ and
$\\{\mathbf{D}^{m}_{i},\mathbf{D}^{o}_{i}\\}$, respectively. The following
lemma identifies a sufficient condition under which, if we are able to train a
classifier on incomplete observations such that the reconstructed data points
are well separated by a hyperplane, then the same classifier correctly
separates the original (unobserved) data vectors.

######  Lemma II.1 (Sufficient condition type I).

Suppose that we have obtained an alternative dictionary
$\mathbf{D}^{\prime}\neq\mathbf{D}\in\mathds{R}^{N\times P}$ such that, for
the incomplete observations ${\mathbf{x}}^{o}_{i}\in\mathds{R}^{M_{i}}$, the
$K$-sparse representation solutions are non-unique, i.e.
$\exists{\mathbf{s}}_{i},{\mathbf{s}}_{i}^{\prime}\in\Sigma_{K}^{P}$ such that
${\mathbf{x}}^{o}_{i}=\mathbf{D}^{o}_{i}{\mathbf{s}}_{i}=\mathbf{D}^{\prime
o}_{i}{\mathbf{s}}_{i}^{\prime}$, where ${\mathbf{s}}_{i}\in\mathds{R}^{P}$
are the vectors of coefficients of the true data and
${\mathbf{s}}_{i}^{\prime}$ provides reconstructions
$\hat{{\mathbf{x}}}_{i}=\mathbf{D}^{\prime}{\mathbf{s}}_{i}^{\prime}$. If a
perfect classifier $\\{{\mathbf{w}},b\\}$ of the reconstructions
$\hat{{\mathbf{x}}}_{i}$ exists s.t.
$|f(\hat{{\mathbf{x}}}_{i})|>\epsilon_{i}>0$ and

$\displaystyle\epsilon_{i}>|\langle{\mathbf{w}}^{m}_{i},{\mathbf{e}}^{m}_{i}\rangle|$
(4)

with ${\mathbf{e}}^{m}_{i}={\mathbf{x}}^{m}_{i}-\hat{{\mathbf{x}}}^{m}_{i}$,
then the full data vectors ${\mathbf{x}}_{i}$ are also perfectly separated
with this classifier, in other words:
$f({\mathbf{x}}_{i})=\langle{\mathbf{w}}_{i},{\mathbf{x}}_{i}\rangle+b>0$
($\leq 0$) if $y_{i}=1$ ($y_{i}=0$).

###### Proof.

By using the missing/observed partition and omitting the sample index $i$, we
can write:
$f({\mathbf{x}})=\langle{\mathbf{w}},{\mathbf{x}}\rangle+b=\langle{\mathbf{w}}^{o},{\mathbf{x}}^{o}\rangle+\langle{\mathbf{w}}^{m},{\mathbf{x}}^{m}\rangle+b$.
If we add and subtract the term
$\langle{\mathbf{w}}^{m},\hat{{\mathbf{x}}}^{m}\rangle$ on the right left
hand, arrange terms and use the fact that
$\hat{{\mathbf{x}}}^{o}=\mathbf{D}^{\prime
o}{\mathbf{s}}^{\prime}={\mathbf{x}}^{o}$, we get:

$f({\mathbf{x}})=f(\hat{{\mathbf{x}}})+\langle{\mathbf{w}}^{m},{\mathbf{e}}^{m}\rangle.$
(5)

Since we assumed that $f(\hat{{\mathbf{x}}})>\epsilon>0$ (for $y_{i}=1$) and
$|\langle{\mathbf{w}}^{m},{\mathbf{e}}^{m}\rangle|<\epsilon$, it implies that
$f({\mathbf{x}})>0$. ∎

Basically, condition (4) means requiring that reconstruction vector
$\hat{{\mathbf{x}}}$ has a distance $\epsilon$ to the separating hyperplane
larger than the absolute dot product between ${\mathbf{w}}$ and the residual
${\mathbf{e}}={\mathbf{x}}-\hat{{\mathbf{x}}}$, which of course is true when
the reconstruction is accurate, i.e. ${\mathbf{x}}\approx\hat{{\mathbf{x}}}$.
However, in practice, reconstructions are not accurate so we are interested in
conditions under which Lemma II.1 can still holds. Below, we derive a more
restrictive but useful sufficient condition:

######  Proposition II.1 (Sufficient condition type II).

Under the same hypothesis of Lemma II.1, the following condition is enough to
guarantee a proper classifier trained on incomplete data:

$\epsilon>|\langle{\mathbf{w}}^{m},{\mathbf{x}}^{m}\rangle|+|\langle{\mathbf{w}}^{m},\hat{{\mathbf{x}}}^{m}\rangle|.$
(6)

###### Proof.

By using the fact that
$|\langle{\mathbf{w}}^{m},{\mathbf{e}}^{m}\rangle|=|\langle{\mathbf{w}}^{m},{\mathbf{x}}^{m}-\hat{{\mathbf{x}}}^{m}\rangle|\leq|\langle{\mathbf{w}}^{m},{\mathbf{x}}^{m}\rangle|+|\langle{\mathbf{w}}^{m},\hat{{\mathbf{x}}}^{m}\rangle|$,
and applying Lemma II.1 the proof is completed. ∎

We highlight that, in our experiments, we were able to verify that Sufficient
Condition type II is met in practice (see section III, Fig. 3).

In Supp. material section V-B, we derive an additional sufficient condition
based on the Restricted Isometry Property (RIP) of the dictionary $\mathbf{D}$
and sparsity level $K$, showing that sufficient condition (6) is easier to
hold for datasets admitting highly sparse representations on dictionaries as
close to orthogonal ones as possible.

####  II-B2 Multilayer-perceptron

Lemma II.1 can be straightforwardly generalized to multilayer-perceptron NNs
where, if a softmax function is used at the output of the last layer then, as
before, the prediction is based on the sign of the linear function:

$\small f({\mathbf{x}})=\langle{\mathbf{w}},{\mathbf{x}}^{(L)}\rangle+b,\mbox{
with
}{\mathbf{x}}^{(l)}=h\left(\mathbf{W}^{T}_{l}{\mathbf{x}}^{(l-1)}+\mathbf{B}_{l}\right),$
(7)

$l=1,2,\dots,L$, where $L+1$ is the total number of layers, $N_{l}$ is the
number of neurons in layer $l$, ${\mathbf{w}}\in\mathds{R}^{N_{L+1}}$ contains
the weights in the last layer, $h(\cdot)$ is an activation function, e.g.
ReLU, $\mathbf{W}_{l}\in\mathds{R}^{N_{l-1}\times N_{l}}$ and
$\mathbf{B}_{l}\in\mathds{R}^{N}_{l}$ contain the weights and biases
associated to neurons at layer $l$; and ${\mathbf{x}}^{(0)}={\mathbf{x}}$ is
the input data vector. In this case, the first layer matrix
$\mathbf{W}_{1}\in\mathds{R}^{N\times N_{1}}$ can be partitioned into
submatrices $\mathbf{W}^{o}_{1i}\in\mathds{R}^{M\times N_{1}}$ and
$\mathbf{W}^{m}_{1i}\in\mathds{R}^{(N-M)\times N_{1}}$ according to the
observed and missing input features, respectively.

######  Proposition II.2.

Under the same conditions of Lemma II.1, if a NN-classifier
$\\{{\mathbf{W}_{l},\mathbf{B}_{l}}(l=1,2,\dots L),{\mathbf{w}},b,h=ReLU\\}$
of the reconstruction $\hat{{\mathbf{x}}}_{i}$ exists such that

$\epsilon_{i}>A\max_{j}|\langle\mathbf{W}^{m}_{1i}(:,j),{\mathbf{e}}^{m}_{i}\rangle|,$
(8)

where $A=\|{\mathbf{w}}\|\prod_{l=2}^{L}\|\mathbf{W}_{l}\|_{2}$ and
${\mathbf{e}}^{m}_{i}={\mathbf{x}}^{m}_{i}-\hat{{\mathbf{x}}}^{m}_{i}$, then
the full data vector ${\mathbf{x}}_{i}$ is also perfectly separated, in other
words: $f({\mathbf{x}}_{i})>0$ ($<0$) if $y_{i}=1$ ($y_{i}=0$).

In the proof of Lemma II.1, we were interested in finding a bound of the
output error when the input ${\mathbf{x}}$ of a classifier is perturbed, i.e.
we found conditions such that
$|f({\mathbf{x}})-f({\mathbf{x}}+\mathbf{\delta})|<\epsilon$. By generalizing
the classifier to the case of a multilayer perceptron we can derive the proof
as follows:

###### Proof.

Given a perturbation $\mathbf{\delta}^{(l-1)}\in\mathds{R}^{N_{l}}$ at the
input of layer $l-1$, i.e.
$\hat{{\mathbf{x}}}^{(l-1)}={\mathbf{x}}^{(l-1)}+\mathbf{\delta}^{(l-1)}$, it
is propagated to the output of layer $l$. By writing the error at the output
we obtain:

$\scriptstyle\mathbf{\delta}^{(l)}=h\left(\mathbf{W}_{l}^{T}{\mathbf{x}}^{(l-1)}+\mathbf{B}_{l}+\mathbf{W}_{l}^{T}\mathbf{\delta}^{(l-1)}\right)-h\left(\mathbf{W}_{l}^{T}{\mathbf{x}}^{(l-1)}+\mathbf{B}_{l}\right),$
(9)

and, by using the sub-additivity of ReLU function $h(\cdot)$, i.e. $h(a+b)\leq
h(a)+h(b)$, we derive the following entry-wise inequality:

$\small\mathbf{\delta}^{(l)}\leq
h\left(\mathbf{W}_{l}^{T}\mathbf{\delta}^{(l-1)}\right),$ (10)

and, by considering the property of ReLU activation function
$\|h({\mathbf{x}})\|\leq\|{\mathbf{x}}\|$, it turns out:

$\small\|\mathbf{\delta}^{(l)}\|\leq\|\mathbf{W}_{l}^{T}\mathbf{\delta}^{(l-1)}\|,$
(11)

Since the last layer of the NN is a linear classifier as in the case of Lemma
II.1, we can ask that
$\langle{\mathbf{w}},\mathbf{\delta}^{(L)}\rangle<\epsilon$. Thus, by
recursively using equation (11), we write

$\langle{\mathbf{w}},\mathbf{\delta}^{(L)}\rangle\leq\|{\mathbf{w}}\|\|\mathbf{\delta}^{(L)}\|\leq\|{\mathbf{w}}\|\|\mathbf{W}_{L}\|_{2}\|\mathbf{W}_{l-1}\|_{2}\cdots\|\mathbf{W}_{2}\|_{2}\|\mathbf{\delta}^{(1)}\|.$
(12)

By defining $A=\|{\mathbf{w}}\|\prod_{l=2}^{L}\|\mathbf{W}_{l}\|_{2}$,
evaluating equation (11) with $l=1$ and taking into account that perturbation
at the input of first layer is $\mathbf{\delta}^{(0)}={\mathbf{e}}$ with
${\mathbf{e}}^{o}=\mathbf{0}$, we arrive at:

$\small\langle{\mathbf{w}},\mathbf{\delta}^{(L)}\rangle\leq
A\|\mathbf{W}_{1}^{mT}{\mathbf{e}}^{m}\|\leq
A\max_{j}|\langle\mathbf{W}^{m}_{1}(:,j),{\mathbf{e}}^{m}\rangle|<\epsilon,$
(13)

which completes the proof.

∎

It is interesting to note that $A=1$ is attained when unit-norm filters
(columns of $\mathbf{W}_{l}$) are orthogonal, which can be imposed by using
orthogonality regularization [35].

##  III Experimental results

We implemented all the algorithms in Pytorch 1.0.0 on a single GPU.
Implementation details are reported in Supplemental material, sections V-C1
and V-C2. The code is available at
111https://github.com/ccaiafa/SimultCodClass.

Synthetic datasets: We synthetically generated $I=11,000$ ($10,000$ training
$+$ $1,000$ test) $K$-sparse data vectors
${\mathbf{x}}_{i}\in\mathds{R}^{100}$ using a dictionary
$\mathbf{D}\in\mathds{R}^{100\times 200}$ obtained from a Gaussian
distribution with normalized atoms, i.e. $\|\mathbf{D}(:,j)\|=1,\forall j$. A
random hyperplane $\\{{\mathbf{w}},b\\}$ with ${\mathbf{w}}\in\mathds{R}^{N}$,
$b\in\mathds{R}$ was randomly chosen dividing data vectors into two classes
according to the sign of the expression
$\langle{\mathbf{w}},{\mathbf{x}}_{i}\rangle+b$, which defined the label
$y_{i}$. We also controlled the degree of separation between classes by
discarding all data vectors with distances to the hyperplane lower than a pre-
specified threshold, i.e. $|\langle{\mathbf{w}},{\mathbf{x}}_{i}\rangle+b|<d$.
We used $n=10$ repetitions of each experiment with different masks and input
data in order to compute statistics.

We applied our simultaneous method (Simult.) with hyperparameters
$\lambda_{1}$ and $\lambda_{2}$ in the cost function (II) tuned via cross-
validation to train a logistic regression classifier on incomplete datasets
with randomly distributed missing features. Then, we computed the
classification accuracy on the complete test dataset and compared the results
against the following standard sequential methods:
Sequential Sparsity based (Seq. Sp.): reconstructions are obtained by finding
the sparsest representation compatible with the observations solving a LASSO
problem. We used Algorithm 3 as shown in the Supp. material;
Zero Fill (ZF): missing features are filled with zeros, which is equivalent to
ignore unknown values;
Mean Unsupervised (MU): missing features are filled with the mean computed on
the available values;
Mean Supervised (MS): as in the previous case but the mean is computed on the
same class vectors only;
K-Nearest Neighbor (KNN): as in the previous case but the mean is computed on
the K-nearest neighbors of the same class vectors only.

To compare the performance of classifiers, we computed the mean accuracy $\pm$
standard error of the mean (s.e.m.), with $n=10$, on complete test datasets
using all the methods for two levels of separation between classes
($d=0.0,0.2$), two levels of sparsity ($K=4,32$) and missing features in the
training dataset ranging from $25\%$ to $95\%$ as shown in Fig. 2. Our results
show that the simultaneous algorithm clearly outperforms all the sequential
methods. A t-test was performed to evaluate the statistical significance with
$p<0.05$ of the difference between our algorithm and MS. It is interesting to
note that, when classes has some degree of separation ($d=0.2$), using the
simple MS method, can give good results but not better than our algorithm.

Figure 2: Experimental results on synthetic dataset with random masks using
our algorithm (red) and compared to various sequential methods. Test accuracy
(mean $\pm$ s.e.m with $n=10$) is shown as a function of the percentage of
missing features for separation of classes $d=0.0,0.2$ and levels of sparsity
$K=4,32$. Statistical significance for the difference between Simult. and MS
is shown ($p<0.05$).

In the second experiment, we generated $I=10,000$ $K$-sparse data vectors
${\mathbf{x}}_{i}\in\mathds{R}^{100}$ using
$\mathbf{D}\in\mathds{R}^{100\times 100}$ and we evaluated the sufficient
condition of equation (6) on $n=10$ repetitions of the experiment with $95\%$
missing features and separation $d=0.0$. Fig. 3 clearly shows that the
sufficient condition is mostly met in practice, especially for highly sparse
representations of input data (small $K$). This means that in practice it is
not necessary to accurately reconstruct the input vectors, it is enough to
capture the intrinsic characteristics of the classes such that the distances
of reconstructions to the separating hyperplane satisfy the sufficient
condition (3).

Figure 3: Verification of the sufficient condition (6) for various levels of
sparsity $K$: 2D-histogram of $\epsilon$ versus
$g=|\langle{\mathbf{w}}^{m},{\mathbf{x}}^{m}\rangle|+|\langle{\mathbf{w}}^{m},\hat{{\mathbf{x}}}^{m}\rangle|$.
Mean + s.e.m ($n=10$) percentage of correctly classified data samples are
shown for $\epsilon>g$ and $\epsilon<g$.

Benchmark datasets: We also considered three popular computer vision datasets:
MNIST [36] and Fashion [37] consisting of 70,000 images (60,000 train + 10,000
test) each; and CIFAR10 [38] having 60,000 images (50,000 train + 10,000
test). MNIST/Fashion datasets contains $28\times 28$ gray scale images while
CIFAR10 dataset is built upon $32\times 32\times 3$ color images of different
objects. The corresponding data sample size is $N=28\times 28=784$ for
MNIST/Fashion and $N=32\times 32\times 3=3,072$ for CIFAR10. We considered a
dictionary of size $784\times 784$ (MNIST/Fashion) and $1,024\times 1,024$
(CIFAR10) and applied our simultaneous algorithm to learn the classifier on
incomplete data using uniform random missing masks with several levels of
missing data (25%, 50% and 75%) and 50% for random partial occlusions with
MNIST/Fashion.

We used a logistic regression classifier (single layer NN) and a 4-layer
convolutional neural network [39] (CNN4) for the MNIST/Fashion dataset using
batch normalization (BN) [40] in the Fashion dataset. For CIFAR10 dataset, an
18-layer residual neural network, Resnet-18 [41] was implemented. We did not
use any data augmentation strategy. The hyper-parameters $\lambda_{1}$ and
$\lambda_{2}$ in cost function (II) were adjusted by cross-validation through
a grid-search, as shown in Supp. material (Table IV and Fig. 5). We compared
our proposed algorithm with the following standard sequential methods: ZF, MS,
KNN-10, KNN-20, KNN-50 and KNN-100; and against the recently proposed method
from [15], referred here as NN-GMM, which uses the same NN classifier as in
our method and models missing features through
GMM222https://github.com/lstruski/Processing-of-missing-data-by-neural-
networks. We trained the classifiers on incomplete data with random masks and
tested them on complete data for MNIST and CIFAR10 datasets. The obtained mean
Test Accuracy $\pm$ s.e.m ($n=10$) are reported in Table I. It is noted that
NN-GMM provided good results with MNIST dataset compared to sequential
methods, however, our simultaneous method outperformed all the methods.
Interestingly, NN-GMM performed worst than any other method with CIFAR10
dataset. It seems that NN-GMM is not robust to large amount of missing data
because, when we reduced the missing entries to $10\%$, the test accuracy
sensibly increased to $52.57\%$. Additionally, our method showed to have
little variability (small s.e.m) compared to the second best method (NN-GMM
for MNIST and ZF for CIFAR10).

TABLE I: Test accuracy (mean $\pm$ s.e.m with $n=10$) of various methods
trained on incomplete data and tested on complete ones for MNIST and CIFAR10.

MNIST (CNN4)
---
Miss. | ZF | MS | KNN10 | KNN20 | KNN50 | KNN100 | NN-GMM | Simult.
75% | $84.86\pm 0.02$ | $83.79\pm 0.01$ | $88.16\pm 0.01$ | $87.94\pm 0.01$ | $87.03\pm 0.002$ | $86.52\pm 0.01$ | $96.36\pm 0.12$ | $\bf 98.09\pm 0.04$
50% | $90.13\pm 0.06$ | $88.55\pm 0.01$ | $91.36\pm 0.02$ | $91.11\pm 0.02$ | $90.87\pm 0.01$ | $90.82\pm 0.01$ | $97.57\pm 0.37$ | $\bf 98.23\pm 0.10$
CIFAR10 (Resnet18)
Miss. | ZF | MS | KNN10 | KNN20 | KNN50 | KNN100 | NN-GMM | Simult.
75% | 32.22 $\pm$ 2.09 | 21.30 $\pm$ 0.40 | 22.84 $\pm$ 0.87 | 25.67 $\pm$ 0.80 | 26.52 $\pm$ 0.70 | 26.01 $\pm$ 0.52 | $12.10\pm 0.61$ | $\bf 54.81\pm 0.47$
50% | 46.37 $\pm$ 1.93 | 17.90 $\pm$ 0.94 | 30.94 $\pm$ 0.54 | 29.68 $\pm$ 0.46 | 30.01 $\pm$ 0.51 | 26.23 $\pm$ 1.01 | $14.02\pm 0.75$ | $\bf 62.50\pm 0.95$

In Table II, test accuracies obtained when the learned model is applied to
incomplete and complete test data, are shown. The right-most column shows the
baseline results obtained by training the model on complete datasets using a
CNN4 [39] and a Resnet-18 [41], whose implementations can be found at
333https://github.com/pytorch/examples/tree/master/mnist and
444https://github.com/kuangliu/pytorch-cifar. It is interesting to note that
for the logistic regression classifier, we obtained better results when
training with incomplete data rather than using complete data. Also, it is
highlighted that training on incomplete data with 50% or fewer random missing
features, provides similar test accuracy as training on complete data for
MNIST dataset. This could be explained by noting two facts: (1) random missing
features is similar to applying dropout, with the exception that missing data
do not change during training; and (2) our model has more parameters
(Dictionary + sparse coefficients + Linear layer) compared to the baseline
logistic regression classifier. To provide a deeper understanding of this
effect, we ran the baseline with Dropout at the input and we obtained: 91.95%,
91.97% and 92.01% for $p=0.0$, $0.1$ and $0.25$, respectively, which shows
that the improvement we obtained with our method is not solely caused by a
dropout alike behavior.

TABLE II: Test Accuracies obtained with our method on MNIST, Fashion and
CIFAR10 datasets training with incomplete data and testing on
incomplete/complete data. Baseline results obtained by training the models on
complete data are shown for reference in the right-most column.

Dataset | Classifier | Random missing features | Occlusion | Baseline
---|---|---|---|---
%Train / %Test | %Train / %Test | %Train / %Test
75/75 | 50/50 | 25/25 | 75/0 | 50/0 | 25/0 | 50/50 | 50/0 | 0/0
MNIST | Log. Reg. | $90.45$ | $93.68$ | $94.14$ | $91.94$ | $93.44$ | $94.43$ | - | - | $91.95$
CNN4 | $94.62$ | $98.34$ | $98.94$ | $98.09$ | $98.23$ | $98.95$ | $88.55$ | $91.37$ | $98.95$
Fashion | CNN4+BN | $83.71$ | $86.09$ | $86.38$ | $86.39$ | $87.11$ | $87.04$ | $81.73$ | $82.47$ | $90.76$
CIFAR10 | Resnet18 | $53.82$ | $61.08$ | $63.73$ | $54.81$ | $62.50$ | $63.87$ | - | - | $80.13$

In Fig. 4, we present some randomly selected visual examples comparing the
original images in the MNIST/Fashion test dataset, their observations using
random masks and partial occlusions, and the reconstructions using the
dictionary learned from the incomplete training data. It is clear that,
despite the reconstructions may be not very similar to the original images
(see “5” digit example), they clearly own the properties of the class to which
they belong to. Additional examples are provided in Supp. material (Figs. 6
and 7).

Figure 4: Original (top), observed (middle) and reconstructed (bottom) MNIST
and Fashion test images.

##  IV Discussion

It is well known that sparse coding has the ability to accurately model
complex distributions of data, such as natural signals (images, audio, EEG,
etc). In this work, we demonstrated that assuming a sparse representation for
input data allows for the successful training of a general NN when incomplete
data is given outperforming traditional sequential approaches and other start-
of-the-art methods. It is highlighted that our method can be used with
potentially any deep NN architecture, thus relying on their extraordinary
capability to accommodate complex decision boundaries as usually needed in
modern machine learning.

Our method overcomes well known issues of previous approaches: (1) compared to
imputation methods, our algorithm successfully incorporates the labelling
information into the modeling of missing features; (2) sparse coding allows
for a simple way to train dictionaries through linear methods such as
stochastic gradient descent with back-propagation compared to the very
expensive EM estimators for GMM used in probabilistic generative models, or
SVD based algorithms for matrix rank minimization in matrix completion; (3)
sparse coding can be more accurate modeling missing values in natural signals
compared to GMM, especially for high dimensional data where GMM may require a
huge number of parameters making it computationally prohibitive.

We analyzed the limitations of the classical imputation approach and
demonstrated through experiments with synthetical and real-world datasets that
our simultaneous algorithm always outperforms them for various cases such as
LASSO, zero-filling, supervised/unsupervised mean and KNN based methods as
well as the state-of-the-art method based on NNs and GGM recently proposed in
[15]. Nevertheless, our experimental results on synthetic and real-world
dataset showed that, even though we only constrained dictionaries to have
unit-norm columns but not enforcing any other kind of constraint like maximum
coherence, the obtained results seem to be satisfactory enough. However,
further analysis on the required properties of dictionaries could provide
deeper insights and alternative ways to improve the algorithm, which we aim to
address in a future work.

While current simple sub-gradient based optimization approach provided
satisfactory results in terms of performance, it is remarked that observed
convergence is slow requiring a thousand of iterations sometimes. We believe,
it could be improved by trying to incorporate some second-order derivatives
information for computing the updates. Although, full Hessian computation
becomes prohibitive with multi-layer NNs a diagonal approximation approach
could be explored. Also, a rigorous convergence analysis in the line of the
analysis in [31, 32] and taking special properties of multi-layer NN
classifier functions can be conducted in a future work.

Finally, we provided theoretical insights of the problem by providing
sufficient conditions under which, if it is possible to train a classifier on
incomplete observations so that its reconstructions are well separated by a
hyperplane, then the same classifier also correctly separates the original
(unobserved) data samples.

Acknowledgments: We are thankful for the RAIDEN computing system and its
support team at RIKEN AIP, Tokyo. This work was supported by the JSPS KAKENHI
(Grant No. 20H04249, 20H04208).

## References

  * [1] G. Chechik, G. Heitz, G. Elidan, P. Abbeel, and D. Koller, “Max-margin Classification of Data with Absent Features.” _Journal of Machine Learning Research (JMLR)_ , vol. 9, pp. 1–21, 2008.
  * [2] R. J. A. Little and D. B. Rubin, _Statistical Analysis with Missing Data_. John Wiley & Sons, Aug. 2014\.
  * [3] J. Josse, N. Prost, E. Scornet, and G. Varoquaux, “On the consistency of supervised learning with missing values,” _arXiv.org_ , p. arXiv:1902.06931, Feb. 2019.
  * [4] P. J. García-Laencina, J.-L. Sancho-Gómez, and A. R. Figueiras-Vidal, “Pattern classification with missing data: a review,” _Neural Computing and Applications_ , vol. 19, no. 2, pp. 263–282, 2009.
  * [5] Z. Ghahramani and M. I. Jordan, “Supervised learning from incomplete data via an EM approach.” in _NIPS_ , 1993.
  * [6] D. Williams, X. Liao, Y. Xue, and L. Carin, “Incomplete-data classification using logistic regression.” _ICML_ , 2005.
  * [7] C. Bhattacharyya, P. K. Shivaswamy, and A. J. Smola, “A Second Order Cone programming Formulation for Classifying Missing Data.” in _NIPS_ , 2004\.
  * [8] X. Liao, H. Li, and L. Carin, “Quadratically gated mixture of experts for incomplete data classification.” in _ICML_ , 2007.
  * [9] U. Dick, P. Haider, and T. Scheffer, “Learning from incomplete data with infinite imputations.” in _ICML_ , 2008.
  * [10] T. I. Lin, J. C. Lee, and H. J. Ho, “On fast supervised learning for normal mixture models with missing information,” _Pattern Recognition_ , vol. 39, no. 6, pp. 1177–1187, 2006.
  * [11] O. Delalleau, A. C. Courville, and Y. Bengio, “Efficient EM Training of Gaussian Mixtures with Missing Data,” _arXiv_ , vol. cs.LG, p. 1209.0521, 2012.
  * [12] A. B. Goldberg, X. Zhu, B. Recht, J.-M. Xu, and R. D. Nowak, “Transduction with Matrix Completion - Three Birds with One Stone.” in _NIPS_ , 2010.
  * [13] E. Hazan, R. Livni, and Y. Mansour, “Classification with Low Rank and Missing Data.” in _ICML_ , 2015.
  * [14] S.-J. Huang, M. Xu, M.-K. Xie, M. Sugiyama, G. Niu, and S. Chen, “Active Feature Acquisition with Supervised Matrix Completion,” _arXiv_ , vol. cs.LG, p. 1802.05380, 2018.
  * [15] M. Smieja, L. Struski, J. Tabor, B. Zielinski, and P. Spurek, “Processing of missing data by neural networks.” in _NeurIPS_ , 2018.
  * [16] H. Lee, A. Battle, R. Raina, and A. Y. Ng, “Efficient sparse coding algorithms.” in _NIPS_ , 2006.
  * [17] J. Mairal, F. R. Bach, J. Ponce, and G. Sapiro, “Online dictionary learning for sparse coding.” in _ICML_ , 2009.
  * [18] K. Gregor and Y. LeCun, “Learning Fast Approximations of Sparse Coding.” in _ICML_ , 2010.
  * [19] S. Mallat, _A Wavelet Tour of Signal Processing_ , ser. The Sparse Way. Academic Press, 2009.
  * [20] M. Davenport, M. F. Duarte, Y. C. Eldar, and G. Kutyniok, _Introduction to compressed sensing_ , ser. Theory and Applications. Cambridge: Cambridge University Press, 2012.
  * [21] B. A. Olshausen and D. J. Field, “Emergence of simple-cell receptive field properties by learning a sparse code for natural images.” _Nature_ , vol. 381, no. 6583, pp. 607–609, 1996.
  * [22] B. Olshausen and D. Field, “Sparse coding with an overcomplete basis set: A strategy employed by V1?” _Vision research_ , vol. 37, no. 23, pp. 3311–3325, 1997.
  * [23] K. Huang and S. Aviyente, “Sparse Representation for Signal Classification.” in _NIPS_ , 2006.
  * [24] J. Mairal, J. Ponce, G. Sapiro, A. Z. A. i. neural, and 2009, “Supervised dictionary learning,” in _NIPS_ , 2008.
  * [25] I. Ramirez, P. Sprechmann, and G. Sapiro, “Classification and clustering via dictionary learning with structured incoherence and shared features.” in _CVPR_ , 2010.
  * [26] P. Sprechmann and G. Sapiro, “Dictionary learning and sparse coding for unsupervised clustering.” in _ICASSP_ , 2010.
  * [27] I. Tosic and P. Frossard, “Dictionary Learning,” _IEEE Signal Processing Magazine_ , vol. 28, no. 2, pp. 27–38, 2011.
  * [28] R. Jenatton, J. Mairal, G. Obozinski, and F. R. Bach, “Proximal Methods for Sparse Hierarchical Dictionary Learning.” in _ICML_ , 2010.
  * [29] J. Weed, “Approximately Certifying the Restricted Isometry Property is Hard,” _IEEE Transactions on Information Theory_ , vol. 64, no. 8, pp. 5488–5497, 2017.
  * [30] E. J. Candès and T. Tao, “Decoding by Linear Programming,” _IEEE Transactions on Information Theory_ , vol. 51, no. 12, pp. 4203–4215, 2005.
  * [31] P. Tseng and S. Yun, “A coordinate gradient descent method for nonsmooth separable minimization,” _Mathematical Programming_ , vol. 117, no. 1-2, pp. 387–423, 2007.
  * [32] E. T. Hale, W. Yin, and Y. Zhang, “Fixed-Point Continuation for $\ell_{1}$-Minimization: Methodology and Convergence,” _SIAM J. OPTIM._ , vol. 19, no. 3, pp. 1107–1130, 2008.
  * [33] S. Shalev-Shwartz and A. Tewari, “Stochastic Methods for l1-regularized Loss Minimization.” _Journal of Machine Learning Research (JMLR)_ , 2011.
  * [34] T. Hastie, R. Tibshirani, and J. H. Friedman, _The elements of statistical learning - data mining, inference, and prediction, 2nd Edition._ , ser. Springer Series in Statistics. New York, NY: Springer New York, 2009.
  * [35] N. Bansal, X. Chen, and Z. Wang, “Can We Gain More from Orthogonality Regularizations in Training Deep CNNs?” in _NeurIPS_ , 2018.
  * [36] Y. LeCun, B. E. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. E. Hubbard, and L. D. Jackel, “Handwritten Digit Recognition with a Back-Propagation Network.” in _NIPS_ , 1989.
  * [37] H. Xiao, K. Rasul, and R. Vollgraf, “Fashion-MNIST: a Novel Image Dataset for Benchmarking Machine Learning Algorithms,” _arXiv_ , vol. cs.LG, p. arXiv:1708.07747, 2017.
  * [38] A. Krizhevsky, “Learning multiple layers of features from tiny images,” Ph.D. dissertation, Toronto University, Toronto, 2009.
  * [39] Y. LeCun, P. Haffner, L. Bottou, and Y. Bengio, “Object Recognition with Gradient-Based Learning.” in _Shape, Contour and Grouping in Computer Vision_ , 1999.
  * [40] S. Ioffe and C. Szegedy, “Batch Normalization - Accelerating Deep Network Training by Reducing Internal Covariate Shift.” in _ICML_ , 2015.
  * [41] K. He, X. Zhang, S. Ren, and J. Sun, “Deep Residual Learning for Image Recognition,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition (CVPR_. IEEE, 2016, pp. 770–778.
  * [42] V. Naumova and K. Schnass, “Dictionary learning from incomplete data for efficient image restoration.” _EUSIPCO_ , 2017.
  * [43] J. Mairal, M. Elad, and G. Sapiro, “Sparse Representation for Color Image Restoration,” _IEEE Transactions on Image Processing_ , vol. 17, no. 1, pp. 53–69, Jan. 2008.

##  V Supplemental material

###  V-A Additional pseudocodes

Here, additional pseudocode of the algorithms discussed in the paper are
provided. Once the classifier is trained by using Algorithm 1, we are able to
apply it to incomplete test data by using Algorithm 2, where for fixed
$\Theta$ and $\mathbf{D}$, we need to find the corresponding sparse
coefficients ${\mathbf{s}}_{i}$, compute the full data vector estimations and,
finally, apply the classifier.

A sparsity-based sequential method is presented in Algorithm 3 (sequential
approach), which consists on learning first the optimal dictionary
$\mathbf{D}$ and sparse coefficients ${\mathbf{s}}_{i}$ compatible with the
incomplete observations (dictionary learning and coding phase), followed by
the training phase, where the classifier weights are tuned in order to
minimize the classification error of the reconstructed input data vectors
$\hat{{\mathbf{x}}}_{i}=\mathbf{D}{\mathbf{s}}_{i}$. It is noted that for the
imputation stage (lines 2-12) other and more specialized dictionary learning
algorithms with missing data can be applied, such as the ones proposed in [42]
for high-dimensional data or [43] for color image data.

Algorithm 2 : Testing on incomplete data

0: Incomplete data vectors $\\{{\mathbf{x}}^{o}_{i}\\}$, $i=1,2,\dots,I$,
classifier parameters $\Theta$, dictionary $\mathbf{D}$, hyper-parameters
$\lambda_{1}$ and $\lambda_{2}$, number of iterations $N_{iter}$ and update
rate $\sigma_{{\mathbf{s}}}$

0: ${\hat{y}_{i}}$ and reconstructions
$\hat{{\mathbf{x}}}_{i}=\mathbf{D}{\mathbf{s}}_{i},\forall i$

1: Sparse coding stage: for fixed dictionary $\mathbf{D}$ find sparse
representations of observations ${\mathbf{x}}^{o}_{i}$

2: Initialize ${\mathbf{s}}_{i},\forall i$ randomly

3: for $n\leq N_{iter}$ do

4: $\Delta_{i}=-\sigma_{{\mathbf{s}}}\big{[}\lambda_{1}\frac{\partial
J_{1}}{\partial{\mathbf{s}}_{i}}+\lambda_{2}\frac{\partial
J_{2}}{\partial{\mathbf{s}}_{i}}\big{]}$, $\forall i$

5: if ${\mathbf{s}}_{i}(j)[{\mathbf{s}}_{i}(j)+\Delta_{i}(j)]<0$ then

6: $\Delta_{i}(j)=-{\mathbf{s}}_{i}(j)$; avoid zero crossing

7: end if

8: ${\mathbf{s}}_{i}={\mathbf{s}}_{i}+\Delta_{i}$, $\forall i$

9: end for

10: $\hat{{\mathbf{x}}}_{i}=\mathbf{D}{\mathbf{s}}_{i},\forall i$; Compute
reconstructions

11: Classification stage: apply classifier to reconstructions
$\hat{{\mathbf{x}}}_{i}$

12:
$\hat{y}_{i}=\operatorname*{arg\,max}_{y}(p^{y}_{\Theta}(\hat{{\mathbf{x}}}_{i}))$

13: return $\Theta,\hat{y}_{i},{\mathbf{s}}_{i},\hat{{\mathbf{x}}}_{i},\forall
i$

Algorithm 3 : Sequential sparsity based approach

0: Incomplete data vectors and their labels
$\\{{\mathbf{x}}^{o}_{i},y_{i}\\}$, $i=1,2,\dots,I$, hyper-parameters
$\lambda_{1}$ and $\lambda_{2}$, number of iterations $N_{iter}$ and update
rate $\sigma_{\Theta}$, $\sigma_{\mathbf{D}}$ and $\sigma_{{\mathbf{s}}}$

0: Classifier weights $\Theta$ and reconstructions
$\hat{{\mathbf{x}}}_{i}=\mathbf{D}{\mathbf{s}}_{i},\forall i$

1: Randomly initialize $\mathbf{D},{\mathbf{s}}_{i},\forall i$

2: Imputation stage: learning of $\mathbf{D}$ and ${\mathbf{s}}_{i}$

3: for $n\leq N_{iter}$ do

4: $\mathbf{D}=\mathbf{D}-\sigma_{\mathbf{D}}\frac{\partial
J_{1}}{\partial\mathbf{D}}$

5: Normalize columns of matrix $\mathbf{D}$

6: $\Delta_{i}=-\sigma_{{\mathbf{s}}}\big{[}\lambda_{1}\frac{\partial
J_{1}}{\partial{\mathbf{s}}_{i}}+\lambda_{2}\frac{\partial
J_{2}}{\partial{\mathbf{s}}_{i}}\big{]}$, $\forall i$

7: if ${\mathbf{s}}_{i}(j)[{\mathbf{s}}_{i}(j)+\Delta_{i}(j)]<0$ then

8: $\Delta_{i}(j)=-{\mathbf{s}}_{i}(j)$; avoid zero crossing

9: end if

10: ${\mathbf{s}}_{i}={\mathbf{s}}_{i}+\Delta_{i}$, $\forall i$

11: end for

12: $\hat{{\mathbf{x}}}_{i}=\mathbf{D}{\mathbf{s}}_{i},\forall i$; Compute
reconstructions

13: Training stage: update $\Theta$

14: for $n\leq N$ do

15: $\Theta=\Theta-\sigma_{\Theta}\frac{\partial J_{0}}{\partial\Theta}$;

16: end for

17: return
$\Theta,\mathbf{D},{\mathbf{s}}_{i},\hat{{\mathbf{x}}}_{i}=\mathbf{D}{\mathbf{s}}_{i},\forall
i$

###  V-B A condition based on RIP and sparsity

The Restricted Isometry Property (RIP): An overcomplete dictionary
$\mathbf{D}$ satisfies the RIP of order $K$ if there exists
$\delta_{K}\in[0,1)$ s.t.

$(1-\delta_{K})\|{\mathbf{s}}\|^{2}_{2}\leq\|\mathbf{D}{\mathbf{s}}\|^{2}_{2}\leq(1+\delta_{K})\|{\mathbf{s}}\|^{2}_{2},$
(14)

holds for all ${\mathbf{s}}\in\Sigma_{K}^{P}$. RIP was introduced in [30] and
characterizes matrices which are nearly orthonormal when operating on sparse
vectors.

In the following theorem, we show that, by imposing conditions on the sparsity
level of the representation and the RIP constant of a sub-matrix of the
dictionary, we can guarantee to meet the sufficient condition (6).

######  Theorem V.1.

Given a dataset $\\{{\mathbf{x}}_{i},y_{i}\\}$, $i=1,2,\dots,I$ with
normalized data vectors ($\|{\mathbf{x}}_{i}\|\leq 1$) admitting a $K$-sparse
representation over a dictionary $\mathbf{D}\in\mathds{R}^{N\times P}$ with
unit-norm columns, whose sub-matrices $\mathbf{D}^{m}_{i}$ satisfy the RIP of
order $K$ with constant $\delta^{i}_{K}$, and suppose that, we have obtained
an alternative dictionary $\mathbf{D}^{\prime}\in\mathds{R}^{N\times P}$,
whose sub-matrices $\mathbf{D}^{\prime m}_{i}$ also satisfy the RIP of order
$K$ with constant $\delta^{i}_{K}$ such that, for the incomplete observation
${\mathbf{x}}^{o}_{i}\in\mathds{R}^{M_{i}}$, the $K$-sparse representation
solution is non-unique, i.e.
$\exists{\mathbf{s}}_{i},{\mathbf{s}}_{i}^{\prime}\in\Sigma_{K}^{P}$ such that
${\mathbf{x}}^{o}_{i}=\mathbf{D}^{o}_{i}{\mathbf{s}}_{i}=\mathbf{D}^{\prime
o}_{i}{\mathbf{s}}_{i}^{\prime}$, where ${\mathbf{s}}_{i}\in\mathds{R}^{P}$ is
the vector of coefficients of the true data, i.e.
${\mathbf{x}}_{i}=\mathbf{D}{\mathbf{s}}_{i}$ and ${\mathbf{s}}_{i}^{\prime}$
provides a plausible reconstruction through
$\hat{{\mathbf{x}}}_{i}=\mathbf{D}^{\prime}{\mathbf{s}}_{i}^{\prime}$ with
$\|\hat{{\mathbf{x}}}_{i}\|\leq 1$. If a perfect classifier
$\\{{\mathbf{w}},b\\}$ of the reconstruction $\hat{{\mathbf{x}}}_{i}$ exists
such that
$|f(\hat{{\mathbf{x}}})|=|\langle{\mathbf{w}},\hat{{\mathbf{x}}}\rangle+b|>\epsilon_{i}>0$
and

$\epsilon_{i}>2\|{\mathbf{w}}^{m}_{i}\|_{1}\sqrt{\frac{K}{1-\delta^{i}_{K}}},$
(15)

then the full data vector ${\mathbf{x}}_{i}$ is also perfectly separated with
this classifier, in other words:
$f({\mathbf{x}}_{i})=\langle{\mathbf{w}}_{i},{\mathbf{x}}_{i}\rangle+b>0$
($\leq 0$) if $y_{i}=1$ ($y_{i}=0$).

###### Proof.

Let us prove that the sufficient condition (6) is met under the hypothesis of
Theorem V.1. Taking into account that
${\mathbf{x}}_{i}^{m}=\mathbf{D}_{i}^{m}{\mathbf{s}}_{i}$, we can write

$\displaystyle\small|\langle{\mathbf{w}}^{m},\mathbf{D}_{i}^{m}{\mathbf{s}}_{i}\rangle|$
$\displaystyle=\Big{|}\sum_{j=1}^{N-M_{i}}{\mathbf{w}}^{m}(j)\sum_{n=1}^{N}\mathbf{D}_{i}^{m}(j,n){\mathbf{s}}_{i}(n)\Big{|}$
$\displaystyle\leq\sum_{j=1}^{N-M_{i}}|{\mathbf{w}}^{m}(j)|\sum_{n=1}^{N}|\mathbf{D}_{i}^{m}(j,n)||{\mathbf{s}}_{i}(n)|.$
(16)

Since we assumed normalized vectors $\|{\mathbf{x}}_{i}\|\leq 1$, by applying
the left-hand side of the RIP we obtain: $\|{\mathbf{s}}_{i}\|\leq
1/\sqrt{1-\delta^{i}_{K}}$, and, taking into account that
$\|{\mathbf{s}}_{i}\|_{1}\leq\sqrt{K}\|{\mathbf{s}}_{i}\|$ and
$|\mathbf{D}_{i}^{m}(j,n)|\leq 1$ (columns of $\mathbf{D}$ are unit-norm), we
obtain:

$\small|\langle{\mathbf{w}}^{m},\mathbf{D}_{i}^{m}{\mathbf{s}}_{i}\rangle|\leq\sqrt{\frac{K}{1-\delta^{i}_{K}}}\sum_{j=1}^{N-M_{i}}|{\mathbf{w}}^{m}(j)|=\sqrt{\frac{K}{1-\delta^{i}_{K}}}\|{\mathbf{w}}^{m}_{i}\|_{1}.$
(17)

Similarly, using $\hat{{\mathbf{x}}_{i}}^{m}=\mathbf{D}_{i}^{\prime
m}{\mathbf{s}}_{i}^{\prime}$, we can obtain that

$|\langle{\mathbf{w}}^{m},\mathbf{D}_{i}^{\prime
m}{\mathbf{s}}_{i}^{\prime}\rangle|\leq\sqrt{\frac{K}{1-\delta^{i}_{K}}}\|{\mathbf{w}}^{m}_{i}\|_{1}.$
(18)

Putting equations (17) and (18) together with equation (15) complete the proof
of the sufficient condition (6). ∎

###  V-C Experimental results details

####  V-C1 Implementation

We implemented all the algorithms in Pytorch 1.0.0 on a single GPU. The code
is available at555 https://github.com/ccaiafa/SimultCodClass..

Initializations of dictionary $\mathbf{D}$ and coefficients ${\mathbf{s}}_{i}$
were made at random. However, we think some improvements in convergence could
be achieved by using some dedicated dictionaries such as the case of Wavelet
or Cosine Transform matrices.

To update NN weights ($\Theta$), we used standard Stochastic Gradient Descent
(SGD) with learning rate $\sigma_{\Theta}$ and momentum $m$, while for
updating dictionary $\mathbf{D}$ and vector coefficients ${\mathbf{s}}_{i}$,
we used fixed update rate $\sigma=\sigma_{\mathbf{D}}=\sigma_{{\mathbf{s}}}$.
It is noted that we used different update rates for training and testing
stages. In Table III, we report the settings used for experiments for MNIST
and CIFAR10 datasets, which includes Number of iterations $N_{iter}$, batch
size $bs$, learning rate $\sigma_{\Theta}$, momentum $m$, update rate $\sigma$
(training and test).

####  V-C2 Hyperparameter tuning

In Table IV we present the results of the grid search for hyper-parameter
tuning on MNIST and CIFAR10 datasets. We fit our model to the training dataset
for a range of values of parameters $\lambda_{1}$ and $\lambda_{2}$ and apply
it to a validation data set. Figure 5 shows the validation accuracy obtained
with different classifiers and levels of missing entries for MNIST dataset.

TABLE III: Experimental settings for MNIST and CIFAR10 datasets: Number of
iterations $N_{iter}$, batch size $bs$, learning rate $\sigma_{\Theta}$,
momentum $m$, update rate $\sigma$ (training and test)

Dataset | Classifier | $N_{iter}$ | $bs$ | $\sigma_{\Theta}$ | $m$ |  $\sigma$ (train) |  $\sigma$ (test)
---|---|---|---|---|---|---|---
MNIST | Log. Reg. | 3000 | 250 | 0.1 | 0.5 | 1.0 | 5.0
CNN4 | 3500 | 250 | 05 | 0.5 | 0.4 | 0.5
CIFAR10 | Resnet18 | 1000 | 64 | 0.01 | 0.5 | 1.0 | 2.5

####  V-C3 Additional visual results

To visually evaluate our results, additional randomly selected examples of
original (complete) images of the test dataset in MNIST and Fashion, together
with their given incomplete observations and obtained reconstructions, are
shown in Figure 6 and Figure 7

TABLE IV: Hyper-parameter tuning: crossvalidated hyperparameters $\lambda_{1}$
and $\lambda_{2}$ obtained for MNIST and CIFAR10 datasets with the classifiers
used in our experiments.

Dataset | Classifier | Random missing entries | Occlusion
---|---|---|---
75% | 50% | 25% | 50%
$\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{1}$ | $\lambda_{2}$ | $\lambda_{1}$ | $\lambda_{2}$
MNIST | Log. Reg. | 0.32 | 1.28 | 0.64 | 1.28 | 0.64 | 1.28 | - | -
CNN4 | 1.28 | 1.28 | 2.56 | 1.28 | 5.12 | 1.28 | 10.24 | 10.24
CIFAR10 | Resnet18 | 0.024 | 0.008 | 0.032 | 0.004 | 0.032 | 0.01 | - | -

Figure 5: Test accuracy in the grid search for hyper-parameter tuning in MNIST
dataset: $\lambda_{1}$ and $\lambda_{2}$ were tuned by cross-validation for
various levels of missing entries: 25%, 50% and 75%. Figure 6: Reconstructions
of incomplete test MNIST dataset images by applying our simultaneous
classification and coding algorithm with the CNN4 architecture. Figure 7:
Reconstructions of incomplete test Fashion dataset images by applying our
simultaneous classification and coding algorithm with the CNN4 architecture
with Batch Normalization (BN).


# Direct Evidence of a Dual Cascade in Gravitational Wave Turbulence

Sébastien Galtier  Laboratoire de Physique des Plasmas, École polytechnique,
F-91128 Palaiseau Cedex, France  Université Paris-Saclay, IPP, CNRS,
Observatoire Paris-Meudon, France  Institut universitaire de France
<EMAIL_ADDRESS>Sergey V. Nazarenko  Institut de
Physique de Nice, Université Nice-Sophia Antipolis, Parc Valrose, 06108 Nice,
France

###### Abstract

We present the first direct numerical simulation of gravitational wave
turbulence. General relativity equations are solved numerically in a periodic
box with a diagonal metric tensor depending on two space coordinates only,
$g_{ij}\equiv g_{ii}(x,y,t)\delta_{ij}$, and with an additional small-scale
dissipative term. We limit ourselves to weak gravitational waves and to a
freely decaying turbulence. We find that an initial metric excitation at
intermediate wavenumber leads to a dual cascade of energy and wave action.
When the direct energy cascade reaches the dissipative scales, a transition is
observed in the temporal evolution of energy from a plateau to a power-law
decay, while the inverse cascade front continues to propagate toward low
wavenumbers. The wavenumber and frequency-wavenumber spectra are found to be
compatible with the theory of weak wave turbulence and the characteristic
time-scale of the dual cascade is that expected for four-wave resonant
interactions. The simulation reveals that an initially weak gravitational wave
turbulence tends to become strong as the inverse cascade of wave action
progresses with a selective amplification of the fluctuations $g_{11}$ and
$g_{22}$.

††preprint: 1

## Introduction.

Wave turbulence (WT) is a state of a continuous medium with random mutually
interacting waves of weak amplitude excited over a broad range of wavenumbers.
The long-time statistical properties of such a medium have a natural
asymptotic closure induced by the large separation of linear and nonlinear
time scales (Benney and Saffman, 1966; Benney, 1967; Benney and Newell,
1967a). The dynamics of WT is driven by kinetic equations which describe the
redistribution of spectral densities via mainly three- or four-wave resonant
interactions. The kinetic equations have two types of exact stationary power-
law solutions: the zero-flux equilibrium thermodynamic spectra and the finite
flux non-equilibrium Kolmogorov-Zakharov spectra (Zakharov  _et al._ , 1992).
The latter solutions are much more interesting because they describe the
spectral transfer of conserved quantities, such as energy or wave action,
generally between a source and a sink (Nazarenko, 2011; Newell and Rumpf,
2011). The direction of the cascade, direct or inverse, can be found by a
numerical evaluation of the sign of the associated flux. The theory also
offers the possibility to predict the Kolmogorov constant. All these
properties makes WT a very interesting regime to understand in depth the
mechanisms underlying turbulence.

WT is of interest to many physical systems for which theoretical predictions
have been made and verified numerically or experimentally. We have, among
others, capillary waves (Benney and Newell, 1967b; Zakharov and Filonenko,
1967; Lommer and Levinsen, 2002; Pushkarev and Zakharov, 2000; Brazhnikov  _et
al._ , 2001; Deike  _et al._ , 2014; Pan and Yue, 2014; Berhanu  _et al._ ,
2019) and gravity waves (Zakharov and Filonenko, 1966; Falcon  _et al._ ,
2007; Clark di Leoni  _et al._ , 2014; Aubourg and Mordant, 2015) on fluid
surfaces, inertial waves in rotating hydrodynamics (Galtier, 2003; Bellet  _et
al._ , 2006; Yarom and Sharon, 2014; Clark di Leoni and Mininni, 2016; Le Reun
_et al._ , 2017; Monsalve  _et al._ , 2020; Le Reun  _et al._ , 2020; Yokoyama
and Takaoka, 2021), elastic waves on thin vibrating plates (Düring  _et al._ ,
2006; Boudaoud  _et al._ , 2008; Mordant, 2008; Cobelli  _et al._ , 2009;
Chibbaro and Josserand, 2016; Hassaini  _et al._ , 2019), optical waves in
optical fibers (Dyachenko  _et al._ , 1992; Laurie  _et al._ , 2012), waves in
Bose-Einstein condensate (Nazarenko, 2006; Proment  _et al._ , 2012), Kelvin
waves on quantum vortex filaments (Kivotides  _et al._ , 2001; Kozik and
Svistunov, 2004; Nazarenko and Onorato, 2006), magnetostrophic waves in
geodynamo (Galtier, 2014; Bell and Nazarenko, 2019) and magnetohydrodynamic
waves in space plasmas (Galtier  _et al._ , 2000; Galtier, 2006; Meyrand  _et
al._ , 2015, 2016, 2018; David and Galtier, 2019). Recently, a theory of WT
has been developed for gravitational waves (GW) (Galtier and Nazarenko, 2017),
a few years after their first direct detection (Abbott and et al., 2016). A
promising application concerns the primordial universe shortly after the
hypothetical initial singularity. During this period, GW can be produced by
different mechanisms like e.g. first order phase transition Witten (1984);
Krauss (1992) or the merger of primary black holes which can be formed from
the primordial space-time fluctuations (Carr  _et al._ , 2016). A typical
length scale of GW excitation can be $10^{-30}$m. Following this idea, a
scenario of cosmological inflation has been proposed relying on presence of
weak or strong GW turbulence and rapid formation of a condensate via an
inverse cascade (Galtier  _et al._ , 2020). In this scenario, the initial
strong GW bursts are quickly diluted as they propagate through the surrounding
space, resulting in a statistically quasi-homogeneous GW field that is weakly
or strongly non-linear depending on the strength and density of forcing
events. For the weak WT case, a kinetic equation that describes the dynamics
of energy and wave action via four-wave resonant interactions was derived. It
has exact stationary scaling solutions for the one-dimensional (1D) isotropic
spectrum of wave action: $k^{-1}$ for the direct energy cascade and $k^{-2/3}$
for the inverse wave action cascade. Further, an explosive front propagation
in the inverse cascade is predicted and observed numerically with a
phenomenological nonlinear diffusion model where strongly local interactions
are retained (Galtier  _et al._ , 2019). With this model, it is also shown
that the non-stationary isotropic spectrum of wave action is slightly
different from the Kolmogorov-Zakharov solution with a power-law index
$\sim-0.6517$ instead of $-2/3$ for the stationary solution. However, this
regime of WT is realized under assumptions that the initial condition consists
of weak waves with random phases, and that the nonlinearity remains weak and
the phases remain random during the evolution. Study of validity of such
assumptions and robustness of the results with respect to more general (not
necessarily wave dominated) initial conditions requires numerical simulations.
Besides, the underlying nonlinear dynamics has never been explored in the
physical space and this requires the use of numerical simulations.

In this Letter, we report the first direct numerical simulations (DNS) of weak
GW turbulence. This regime is studied in both physical and spectral spaces.
Our study reveals that, despite having its own distinct properties, space-time
turbulence behaves in a classical way from the point of view of the general
turbulence theory.

## Einstein’s equations.

Figure 1: Time evolution (rms values) of the basic fields $\tilde{\alpha}$,
$\tilde{\beta}$ and $\lambda$ (insert).

Following (Galtier and Nazarenko, 2017; Hadad and Zakharov, 2014), we shall
consider the vacuum general relativity equations (Einstein, 1915, 1916) with
the Hadad–Zakharov ($2.5+1$) diagonal metric:

$\displaystyle g_{00}$ $\displaystyle=$
$\displaystyle-(1+\tilde{\gamma})^{2}e^{-2\lambda},\quad
g_{11}=(1+\tilde{\beta})^{2}e^{-2\lambda},$ (1) $\displaystyle g_{22}$
$\displaystyle=$ $\displaystyle(1+\tilde{\alpha})^{2}e^{-2\lambda},\quad
g_{33}=e^{2\lambda}.$ (2)

The metric depends on the space-time variables $x$, $y$ and $t$ but is
independent of $z$. It can describe GW of any amplitude, however, we limit
ourselves to weak GW amplitude, i.e.
$\tilde{\alpha},\tilde{\beta},\tilde{\gamma},\lambda\ll 1$. Within this limit,
Einstein’s equations in the leading order are (Galtier and Nazarenko, 2017):

$\displaystyle\partial_{x}\dot{\tilde{\alpha}}=-2\dot{\lambda}(\partial_{x}\lambda),$
$\displaystyle\quad\partial_{y}\dot{\tilde{\beta}}=-2\dot{\lambda}(\partial_{y}\lambda),$
(3) $\displaystyle\partial_{x}\partial_{y}\tilde{\gamma}$ $\displaystyle=$
$\displaystyle-2(\partial_{x}\lambda)(\partial_{y}\lambda),$ (4)
$\displaystyle\partial_{t}{\left[(1+\tilde{\alpha}+\tilde{\beta}-\tilde{\gamma})\dot{\lambda}\right]}\,$
$\displaystyle=$
$\displaystyle\partial_{x}\left[(1+\tilde{\alpha}-\tilde{\beta}+\tilde{\gamma})\partial_{x}\lambda\right]$
$\displaystyle+$
$\displaystyle\partial_{y}\left[(1-\tilde{\alpha}+\tilde{\beta}+\tilde{\gamma})\partial_{y}\lambda\right],$

where we define $\dot{f}\equiv\partial_{t}f$. The linear solution of this
system (see Eq. (Einstein’s equations.)) is a plus-polarized GW with a
dispersion relation $\omega=k$ (with the speed of light $c=1$ and $k=|{\bf
k}|$) (Maggiore, 2008). For the analysis below, it is useful to recall the
link between the canonical variable $a_{{\bf k}}$ and the primary variables,
$a_{{\bf k}}=\sqrt{\frac{k}{2}}\lambda_{{\bf
k}}+i\sqrt{\frac{1}{2k}}{\dot{\lambda}}_{{\bf k}}$, where $\lambda_{{\bf k}}$
is the Fourier transform of $\lambda$ (Galtier and Nazarenko, 2017). The two-
dimensional (2D) wave action spectrum will be computed numerically using the
relation $N({\bf k})=|a_{{\bf k}}|^{2}$ (homogeneity will be used as well as
plane waves), and the 2D energy spectrum is $E({\bf k})=\omega N({\bf k})$
(Galtier and Nazarenko, 2017).

Figure 2: Time evolution of the normalized energy. Insert: the same variation
in log-log reveals a decay close to $t^{-1/4}$. The colored vertical lines at
the top correspond to the times chosen to plot the spectra in Fig. 3.

DNS of Eqs. (3)–(Einstein’s equations.) are performed with an additional
dissipative term acting at small scales to avoid numerical instabilities
(Pretorius, 2005). This term, added in the RHS of Eq. (Einstein’s equations.),
takes the form $-\nu k^{4}{\dot{\lambda}}_{{\bf k}}$ (for $k\geq k_{diss}$) in
Fourier space (see eg. (Düring  _et al._ , 2006)). Physically, GW dissipation
by matter is expected through e.g. Landau damping (Baym  _et al._ , 2017).
(Note that the non-dissipative Einstein equations become invalid at the Planck
length.) However, in our case, the dissipation has no precise origin and must
be considered as a mechanism mimicking existence of a positive energy flux
which would otherwise be blocked due to presence of the maximum wavenumber in
our numerical method. We developed a pseudo-spectral code using FFTW3
routines, with periodic boundary conditions and dealiasing (Shukla  _et al._ ,
2019). Since we deal with real fields only, the Fourier space is restricted to
$k_{y}\geq 0$. An Adams-Bashforth numerical scheme is used to integrate the
nonlinear terms. For solving the double time derivative in Eq. (Einstein’s
equations.), an intermediate variable $\Lambda=\partial_{t}\lambda$ is
introduced. We also introduce intermediate variables
$A=\partial_{t}\tilde{\alpha}$, $B=\partial_{t}\tilde{\beta}$ and
$G=\partial_{t}\tilde{\gamma}$ for solving Eqs. (3) and (Einstein’s
equations.) where several time-derivative terms appear. Therefore, in practice
$8$ equations are numerically solved at each timestep $\Delta T$. The
simulation shown is made with a spatial resolution of $512$ points in each
direction, $\nu=4\times 10^{-11}$ and $k_{diss}=140$. Initially (at $t=0$)
only $\lambda$ is excited around wavenumber $k_{i}=89$ (in order to see the
dual cascade), with random phases, and such that $|\lambda_{{\bf
k}}|^{2}=C(k^{2}-88^{2})(90^{2}-k^{2})$ with $k_{x,y}\in[88,90]/\sqrt{2}$. The
constant $C$ is chosen so that the total wave action is initially equal to
$50$. The time will be normalized in the linear GW time unit
$t_{GW}=1/\omega_{i}=1/k_{i}$ (characteristic time of the initial excitation).
We take $\Delta T=10^{-5}t_{GW}$, which is relatively small but is necessary
to avoid numerical instabilities and to ensure the conservation of invariants
during the primary phase. The present simulation parameters appear to give the
most representative illustration of the processes we study (we have performed
a large number of simulations with various values of these parameters).

## Results.

Fig. 1 shows the global evolution in time, i.e. the root mean square (rms)
values, of the basic fields $\tilde{\alpha}$, $\tilde{\beta}$ and $\lambda$.
We see that ${\tilde{\alpha}}_{rms}$ and ${\tilde{\beta}}_{rms}$ are strongly
correlated with two distinct stages. At the first stage, until $\sim
3500t_{GW}$, we see two almost identical signals in which a slight growth is
superimposed with oscillations with a typical period of $100t_{GW}$. After
that the signals become different from each other: the high-frequency
oscillations remain but in addition there appear long-scale oscillation (with
a period about ten times longer) in which ${\tilde{\alpha}}_{rms}$ and
${\tilde{\beta}}_{rms}$ oscillate in counter-phase. At the large scale, both
signals increase over time reflecting local amplifications of
${\tilde{\alpha}}(x,y)$ and ${\tilde{\beta}}(x,y)$. Overall, $\tilde{\alpha}$
and $\tilde{\beta}$ behave like twin variables with similar dynamics. The
behavior of ${\lambda}_{rms}$ (shown in the figure insert) and
${\tilde{\gamma}}_{rms}$ (not shown) is different with no significant increase
in amplitude and faster oscillations with a period close to $t_{GW}$.

In Fig. 2, we show the temporal evolution of the normalized energy $E$ (with
initially $E(t=0)=28185$). Two phases are clearly present. First, the energy
is conserved until $\sim 3500t_{GW}$. This phase corresponds to the time
interval needed to develop a direct energy cascade and to reach the small
dissipative scales (see Fig. 3). This observation can be seen as an evidence
of the accuracy of the numerical code (the wave action is also conserved
during this phase). Then, the energy slowly decreases. To appreciate this
decay law, a log-log plot is given in the insert: this reveals a power law
decay in time. This type of behavior is classical in freely decaying (strong
or weak) turbulence in presence of a direct cascade (see eg. (Galtier  _et
al._ , 1997; Bigot  _et al._ , 2008)).

Figure 3: Time evolution of $k^{2/3}N_{1d}(k)$ (solid lines) with the initial
spectrum in black dashed line. Corresponding times are given in Fig. 2. The
same spectrum, renormalized in scale and amplitude, is shown for an inviscid
simulation at a resolution of $1024$ points (blue dashed line). The 1D spectra
$2\pi k\int_{0}^{2\pi}|\tilde{\alpha}_{\bf k}|^{2}d\phi$ and $2\pi
k\int_{0}^{2\pi}|\tilde{\beta}_{\bf k}|^{2}d\phi$ are also shown (dotted
lines) at the final time of the simulation (at a resolution of $512$ points)
and compared with the power law $k^{-3.5}$.

The evolution in time of the 1D wave action spectrum $N_{1d}(k)=2\pi
k\int_{0}^{2\pi}N({\bf k})d\phi$, where $\phi$ is the polar angle in the
k-space, is displayed in Fig. 3. The spectra are compensated by the
theoretical prediction $k^{-2/3}$ (Galtier and Nazarenko, 2017). The
propagation towards small scales is interpreted as a signature of a direct
energy cascade. If we come back to Fig. 2 we see that the spectrum shown in
light green is close to the end of the plateau: this is the moment when the
dissipative scales are reached. Afterwards, the small-scale propagation of the
spectrum is stopped. On the other hand, the development of the front is
interpreted as an inverse cascade of wave action, which could be viewed as a
strongly non-equilibrium Bose-Einstein condensation process. We see that
beyond $3500t_{GW}$ (from red curve onwards) the inverse cascade is preserved
with a further expansion of the inertial range characterized by a bump in the
front propagation (see also Fig. 4). We can already see the presence of a
plateau and conclude that our result is qualitatively in agreement with the
theoretical prediction. The 1D spectra $2\pi
k\int_{0}^{2\pi}|\tilde{\alpha}_{\bf k}|^{2}d\phi$ and $2\pi
k\int_{0}^{2\pi}|\tilde{\beta}_{\bf k}|^{2}d\phi$ are shown in dotted lines:
as expected, $\alpha$ and $\beta$ are much smaller than $\lambda$ (since
$N_{1d}(k)\sim k^{2}|\lambda_{\bf k}|^{2}$) in the inertial range of WT.

Figure 4: $N({\bf k})$ around the final time of the simulation (mean over 4
times). The center of the domain of initial excitation is indicated by the
symbol ’o’.

However, we observe a significant selective amplification of $\alpha$ and
$\beta$ at large scales which is not described by the weak WT theory.
Respectively, it means an amplification of the metric components $g_{11}$ and
$g_{22}$ at the large scales while $g_{00}$ and $g_{33}$ remain in fast
oscillations only (see also Fig. 6). For the GW, such large-scale variations
of $g_{11}$ and $g_{22}$ are perceived as slow variation of the scale factor
of the underlying space, i.e. its expansions and contractions. Superimposed to
these plots, we also show the spectrum of an inviscid (i.e. $\nu=0$)
simulation at $1024$ points resolution with an initial excitation at
$k_{i}=128$ (with same timestep and type of initial condition). The simulation
is stopped before reaching the smallest scales. The spectrum is renormalized
in scale such that the initial excitation is moved to $89$. The purpose of
this plot is to confirm the dual cascade while we are still in the
conservative phase, with a plateau at large scales and, at small scales, a
wider inertial range compared to the $512^{2}$ simulation.

In Fig. 4 the normalized 2D wave action spectrum is shown around the final
time of the simulation to appreciate the number of modes excited during these
cascades. The presence of a bow at large wavenumbers can be seen as the
signature of the dissipation which starts abruptly at $k=140$. This plot
reveals that a large domain in Fourier space is affected by the cascades which
start around the symbol ’o’. Note that we are still far from the axes and no
condensation is found. Note also a dark spot at the smallest excited
wavenumbers corresponding to a spectrum amplification. It points out an
overshoot nature of the propagating condensation front with a significant
localized lump of wave action moving toward the smaller $k$’s. Also the 2D
spectrum looks like a wedge/angle, which means that the propagation toward the
large and the small scales takes place without significant spreading in the
angular distribution of the wave propagation, i.e. without a full
isotropization.

WT is a state dominated by waves of weak amplitude; all non-wave initial
disturbances eventually die out. Such a system is characterized by a very
specific $\omega$–$k$ spectrum that concentrates near the dispersion relation
curve of the wave in question (Nazarenko and Onorato, 2006; Cobelli  _et al._
, 2009; Clark di Leoni  _et al._ , 2014; Yarom and Sharon, 2014; Aubourg and
Mordant, 2015; Meyrand  _et al._ , 2015; Le Reun  _et al._ , 2017; Bell and
Nazarenko, 2019; Hassaini  _et al._ , 2019). It is also the case for GW
turbulence as we can see in Fig. 5 where a $\omega$–$k$ spectrum is plotted.
This plot is obtained by taking the time evolution of the canonical variable
$a_{{\bf k}}$ (real part) for $t/t_{GW}\in[17000,20000]$. Then, we analyze
signals corresponding to $k\in[1,140]$, such that $k_{x}=k_{y}$. A Fourier
transform in time is then applied to each signal weighted with a Hamming
function. The modulus squared of each signal normalized by its maximum is then
plotted. We see that a signal is obtained along the dispersion relation
(dotted line) confirming the wave-like character of this turbulence.

Figure 5: $\omega$–$k$ spectrum of wave action. The dotted line corresponds to
the dispersion relation of a GW.

In Fig. 6, the four metric components are plotted at the final time of the
simulation. We see that the dominant components are $g_{11}$ and $g_{22}$
whose large-scale oscillations are approximately in anti-phase. The two other
metric components $g_{00}$ and $g_{33}$ are characterized by relatively small
small-scale fluctuations around $-1$ and $+1$, respectively. These
fluctuations remain relatively small during the entire simulation. We have
found that $g_{11}$ and $g_{22}$ behave similarly to each other during the
simulation with mainly an anti-correlation and a gradual increase of the
fluctuations. This observation is consistent with the behavior observed in
Fig. 1 where an increase of $\tilde{\alpha}$ and $\tilde{\beta}$ is also
reported. The final time reported here is a reasonable limit to stop the
simulation because the original equations used (3)–(Einstein’s equations.) are
only valid for weak GW (which correspond to small fluctuations around $\pm 1$
of the metric components). Interestingly this final time, $t\sim
10^{4}t_{GW}$, corresponds to the expected time required to develop the weak
turbulence regime when the small parameter $\epsilon$ used for the expansion
is $\sim 0.1$: indeed, a phenomenological evaluation gives a typical cascade
time $\tau_{\rm cascade}\sim\epsilon^{4}t_{GW}$ for four-wave interactions
(whereas it would be $\sim\epsilon^{2}t_{GW}$ for three-wave interactions)
(Nazarenko, 2011). Note that this is an extremely slow time-scale from a
numerical point of view that limits us to a relatively low spatial resolution
(which however proved here to be sufficient to obtain physical results) 111A
run at resolution $1024^{2}$ takes about $6$ months..

Figure 6: Metric components $g_{00}$ (blue), $g_{11}$ (orange), $g_{22}$
(green) and $g_{33}$ (red) at the final time of the simulation. $g_{00}$ and
$g_{33}$ have been vertically shifted by $+1.9$ and $+0.15$, respectively.

## Conclusion.

In this Letter, we have reported the first DNS of GW turbulence. Specifically,
the weak regime was studied for which analytical predictions exist. By using
both physical and Fourier spaces we have been able to show that WT can emerge
from an initial excitation of the space-time metric with a dual cascade of
energy and wave action. This behavior is understood as the result of four-wave
resonant interactions of the $2\leftrightarrow 2$ type for which the wave
action is an invariant. Further, we have observed a new effect which is beyond
the weak WT predictions— emergence and continuous amplification of strong
large-scale fluctuations of metric components $g_{11}$ and $g_{22}$ while the
other components exhibit only weak small-scale oscillations. In particular the
temporal component $g_{00}$ remains close to one which provides a natural
cosmic time for weak GW turbulence. The amplification of $g_{11}$ and $g_{22}$
limits our study since initially weak WT tends to become strong at large
scales (Galtier  _et al._ , 2020). In principle, the regime of strong WT can
then be studied numerically with the metric (1)–(2) (whose form is preserved
at all times as proved by (Hadad and Zakharov, 2014)) by including all
nonlinear terms of Einstein’s equations.

The main conclusion of this work is that it is possible to produce turbulence
in general relativity. Unlike the classical hydrodynamic turbulence, it does
not consist of randomly interacting vortices but, rather, it takes a form of
random interacting waves—the wave turbulence. Further we show that GW
turbulence is a dual cascade system. Namely, in addition to the direct energy
cascade, there is an inverse cascade of wave action. The latter is important
as it may shed light on the processes in early universe (Galtier  _et al._ ,
2020). We can also mention a strong similarity to elastic wave turbulence in
the high tension limit. Indeed, both problems involve four-wave interactions
with an inverse cascade of wave action (Hassaini  _et al._ , 2019). This
similarity could be a motivation to pursue the comparison in an analog
laboratory experiment to better understand strong GW turbulence in cosmology
and the formation of a metric condensate.

## References

  * Benney and Saffman (1966) D. Benney and P. Saffman, Proc. R. Soc. Lond. A 289, 301 (1966).
  * Benney (1967) D. Benney, J. Math. Phys. 46, 115 (1967).
  * Benney and Newell (1967a) D. Benney and A. Newell, J. Math. Phys. 46, 363 (1967a).
  * Zakharov  _et al._ (1992) V. Zakharov, V. L’Vov, and G. Falkovich,  _Kolmogorov spectra of turbulence I: Wave turbulence_ (Springer Series in Nonlinear Dynamics, Berlin: Springer, 1992).
  * Nazarenko (2011) S. Nazarenko,  _Wave Turbulence_ , Lecture Notes in Physics, Vol. 825 (Berlin Springer Verlag, 2011).
  * Newell and Rumpf (2011) A. C. Newell and B. Rumpf, Ann. Rev. Fluid Mech. 43, 59 (2011).
  * Benney and Newell (1967b) D. J. Benney and A. C. Newell, Phys. Fluids 10, S281 (1967b).
  * Zakharov and Filonenko (1967) V. Zakharov and N. Filonenko, J. Applied Mech. Tech. Phys. 8, 37 (1967).
  * Lommer and Levinsen (2002) M. Lommer and M. Levinsen, J. Florescence 12, 45 (2002).
  * Pushkarev and Zakharov (2000) A. Pushkarev and V. Zakharov, Physica D Nonlinear Phenomena 135, 98 (2000).
  * Brazhnikov  _et al._ (2001) M. Y. Brazhnikov, G. Kolmakov, A. Levchenko, and L. Mezhov-Deglin, Sov. J. Exp. Theo. Phys. Lett. 73, 398 (2001).
  * Deike  _et al._ (2014) L. Deike, D. Fuster, M. Berhanu, and E. Falcon, Phys. Rev. Lett. 112, 234501 (2014).
  * Pan and Yue (2014) Y. Pan and D. K. P. Yue, Phys. Rev. Lett. 113, 094501 (2014).
  * Berhanu  _et al._ (2019) M. Berhanu, E. Falcon, G. Michel, C. Gissinger, and S. Fauve, EPL (Europhysics Letters) 128, 34001 (2019).
  * Zakharov and Filonenko (1966) V. Zakharov and N. Filonenko, Doclady Akad. Nauk. SSSR 170, 1292 (1966).
  * Falcon  _et al._ (2007) É. Falcon, C. Laroche, and S. Fauve, Phys. Rev. Lett. 98, 094503 (2007).
  * Clark di Leoni  _et al._ (2014) P. Clark di Leoni, P. J. Cobelli, and P. D. Mininni, Phys. Rev. E 89, 063025 (2014).
  * Aubourg and Mordant (2015) Q. Aubourg and N. Mordant, Phys. Rev. Lett. 114, 144501 (2015).
  * Galtier (2003) S. Galtier, Phys. Rev. E 68, 015301(R) (2003).
  * Bellet  _et al._ (2006) F. Bellet, F. Godeferd, J. Scott, and C. Cambon, J. Fluid Mech. 562, 83 (2006).
  * Yarom and Sharon (2014) E. Yarom and E. Sharon, Nature Physics 10, 510 (2014).
  * Clark di Leoni and Mininni (2016) P. Clark di Leoni and P. Mininni, J. Fluid Mech. 809, 821 (2016).
  * Le Reun  _et al._ (2017) T. Le Reun, B. Favier, A. J. Barker, and M. Le Bars, Phys. Rev. Lett. 119, 034502 (2017).
  * Monsalve  _et al._ (2020) E. Monsalve, M. Brunet, B. Gallet, and P.-P. Cortet, Phys. Rev. Lett. 125, 254502 (2020).
  * Le Reun  _et al._ (2020) T. Le Reun, B. Favier, and M. Le Bars, EPL (Europhysics Letters) 132, 64002 (2020).
  * Yokoyama and Takaoka (2021) N. Yokoyama and M. Takaoka, J. Fluid Mech. 908, A17 (2021).
  * Düring  _et al._ (2006) G. Düring, C. Josserand, and S. Rica, Phys. Rev. Lett. 97, 025503 (2006).
  * Boudaoud  _et al._ (2008) A. Boudaoud, O. Cadot, B. Odille, and C. Touzé, Phys. Rev. Lett. 100, 234504 (2008).
  * Mordant (2008) N. Mordant, Phys. Rev. Lett. 100, 234505 (2008).
  * Cobelli  _et al._ (2009) P. Cobelli, P. Petitjeans, A. Maurel, V. Pagneux, and N. Mordant, Phys. Rev. Lett. 103, 204301 (2009).
  * Chibbaro and Josserand (2016) S. Chibbaro and C. Josserand, Physical Review E 94, 011101(R) (2016).
  * Hassaini  _et al._ (2019) R. Hassaini, N. Mordant, B. Miquel, G. Krstulovic, and G. Düring, Phys. Rev. E 99, 033002 (2019).
  * Dyachenko  _et al._ (1992) S. Dyachenko, A. Newell, A. Pushkarev, and V. Zakharov, Physica D 57, 96 (1992).
  * Laurie  _et al._ (2012) J. Laurie, U. Bortolozzo, S. Nazarenko, and S. Residori, Phys. Rep. 514, 121 (2012).
  * Nazarenko (2006) S. Nazarenko, JETP Lett. 83, 198 (2006).
  * Proment  _et al._ (2012) D. Proment, S. Nazarenko, and M. Onorato, Physica D 241, 304 (2012).
  * Kivotides  _et al._ (2001) D. Kivotides, J. C. Vassilicos, D. C. Samuels, and C. F. Barenghi, Phys. Rev. Lett. 86, 3080 (2001).
  * Kozik and Svistunov (2004) E. Kozik and B. Svistunov, Phys. Rev. Lett. 92, 035301 (2004).
  * Nazarenko and Onorato (2006) S. Nazarenko and M. Onorato, Physica D 219, 1 (2006).
  * Galtier (2014) S. Galtier, J. Fluid Mech. 757, 114 (2014).
  * Bell and Nazarenko (2019) N. Bell and S. Nazarenko, J. Phys. A: Math. Theo. 52, 445501 (2019).
  * Galtier  _et al._ (2000) S. Galtier, S. Nazarenko, A. Newell, and A. Pouquet, J. Plasma Phys. 63, 447 (2000).
  * Galtier (2006) S. Galtier, J. Plasma Phys. 72, 721 (2006).
  * Meyrand  _et al._ (2015) R. Meyrand, K. Kiyani, and S. Galtier, J. Fluid Mech. 770, R1 (2015).
  * Meyrand  _et al._ (2016) R. Meyrand, S. Galtier, and K. H. Kiyani, Phys. Rev. Lett. 116, 105002 (2016).
  * Meyrand  _et al._ (2018) R. Meyrand, K. H. Kiyani, Ö. D. Gürcan, and S. Galtier, Phys. Rev. X 8, 031066 (2018).
  * David and Galtier (2019) V. David and S. Galtier, Astrophys. J. Lett. 880, L10 (2019).
  * Galtier and Nazarenko (2017) S. Galtier and S. V. Nazarenko, Phys. Rev. Lett. 119, 221101 (2017).
  * Abbott and et al. (2016) B. Abbott and et al., Phys. Rev. Lett. 116, 061102 (2016).
  * Witten (1984) E. Witten, Phys. Rev. D 30, 272 (1984).
  * Krauss (1992) L. M. Krauss, Physics Lett. B 284, 229 (1992).
  * Carr  _et al._ (2016) B. Carr, F. Kühnel, and M. Sandstad, Phys. Rev. D 94, 083504 (2016).
  * Galtier  _et al._ (2020) S. Galtier, J. Laurie, and S. Nazarenko, Universe 6, 98 (2020).
  * Galtier  _et al._ (2019) S. Galtier, S. V. Nazarenko, E. Buchlin, and S. Thalabard, Physica D 390, 84 (2019).
  * Hadad and Zakharov (2014) Y. Hadad and V. Zakharov, J. Geometry and Physics 80, 37 (2014).
  * Einstein (1915) A. Einstein, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin), Seite 844-847. (1915).
  * Einstein (1916) A. Einstein, Sitzungsberichte der Königlich Preußischen Akademie der Wissenschaften (Berlin), Seite 688-696. (1916).
  * Maggiore (2008) M. Maggiore,  _Gravitational Waves, Volume 1_ (Oxford Univ. Press, 2008).
  * Pretorius (2005) F. Pretorius, Phys. Rev. Lett. 95, 121101 (2005).
  * Baym  _et al._ (2017) G. Baym, S. P. Patil, and C. J. Pethick, Phys. Rev. D 96, 084033 (2017).
  * Shukla  _et al._ (2019) V. Shukla, P. D. Mininni, G. Krstulovic, P. C. di Leoni, and M. E. Brachet, Phys. Rev. A 99, 043605 (2019).
  * Galtier  _et al._ (1997) S. Galtier, H. Politano, and A. Pouquet, Phys. Rev. Lett. 79, 2807 (1997).
  * Bigot  _et al._ (2008) B. Bigot, S. Galtier, and H. Politano, Phys. Rev. Lett. 100, 074502 (2008).
  * Note (1) A run at resolution $1024^{2}$ takes about $6$ months.