tion network. Amortization gap can then theoretically be reduced by taking more samples for the expectation approximation, its precision increasing with L. However, Monte-Carlo approximation performs poorly, and the precision gained by increasing L does not worth the computational cost of decoding L more values.

**Importance-Weighted Auto-Encoders** A well-known alternative to VAEs are *Importance Weighting Auto-Encoders* (IWAE), that replace the single Monte-Carlo approximation described in previous section by a importance-weighted estimate of the ELBO, taken over several samples from the recognition distribution. The new gradient estimate of the lower-bound $\mathcal{L}_k(q)$ then corresponds to the k-sampled importance weighting of the log-likelihood

$$ 
\begin{aligned}
\nabla_{\theta,\phi} \mathcal{L}_k[q] &= \mathbb{E}_{z_1, \dots, z_n \sim q(z|x)} \left[ \log \frac{1}{k} \sum_{i=1}^{k} \frac{p(x, z)}{q(z|x)} \right] \\
&= \nabla_{\theta,\phi} \mathbb{E}_{\epsilon_1, \dots, \epsilon_k} \left[ \sum_{i=1}^{k} \tilde{w}_k \nabla_{\theta,\phi} \log w(x, z(x, \epsilon)) \right] \\
w_i &= \frac{p(x, z(x, \epsilon_i))}{q(z(x, \epsilon_i)|x)}
\end{aligned}
$$

where we applied the reparametrization trick to obtain k stochastic components $[\epsilon_1, \dots, \epsilon_k]$. In the latter expression, $w_i$ denotes the density ratios, and the weight $\tilde{w}_i = w_i / \sum_{j=1}^k w_j$ is the importance weight of the corresponding sample. Importance-weighting approximation of the ELBO improves the quality of the reconstructions, and at the same time can help to prevent latent over-pruning (see section 2.3.2.3). Cremer & al. argue that the true recognition model recovered by IWAE is slightly different the one recovered by the VAE, as it takes into account the real posterior $p(z|x)$ [177]

$$ q_{IW} = \mathbb{E}_{z_2, \dots, z_k \sim q(z|x)} \left[ \frac{p(x, z)}{\frac{1}{k} \left( \frac{p(x, z)}{q(z|x)} + \sum_{j=2}^{k} \frac{p(x, z_j)}{q(z_j|x)} \right)} \right] $$

where the additional term brought by IWAE is $\sum_{j=2}^{k} p(x, z_j) / q(z_j|x)$. Such methods have then be generalized under the name of *Monte-Carlo objectives* for global log-likelihood estimators [178], giving a series of high-precision estimators such as Filtering Variational Objectives [179] for sequential variational data.

### 2.3.2.3 Latent over-pruning and posterior collapse

One of the main reasons of the AEVB success in the machine learning literature is the emerging structural properties of the extracted representations [22]. According to manifold hypothesis, this tendency would point out that VAEs are able to retrieve the underlying manifold of a given dataset, and recovering the true factors of variation structuring the data. Furthermore, a natural *latent pruning* effect has been noticed during training, pointing out that VAEs automatically prunes latent dimensions that are considered unnecessary for data

<figure><img src="image_33.png" /><figcaption>Figure 2.31: Importance weighting samples several samples from the distribution $q(z|x)$, that are then weighted by their probability in the final loss.</figcaption></figure>
representation. This would indicate that the VAE retrieves the intrinsic dimensionality of the data manifold, and would then structurally protected against over-fitting. However, further investigations of AEVB training dynamics tempered these observations, exhibiting undesirable effects occurring during the training process and mitigating the natural *disentanglement* nature of the extracted manifold.

**Latent over-pruning** While dimensional pruning can be taken as a desirable property of the system, VAEs has been showed to deactivate latent dimensions even when the reconstruction quality was insufficient, *over-pruning* the latent representation. This deficiency can be intuited recalling the formulation (2.18) of the ELBO

$$p(\mathbf{x}) \geq \mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x})}[p_{\theta}(\mathbf{x}|\mathbf{z})] - D_{KL}[q_{\phi}(\mathbf{z}|\mathbf{x})||p(\mathbf{z})]$$

observing that the distribution $q_{\phi}(\mathbf{z}|\mathbf{x})$ is pushed to match the overall prior $p(\mathbf{z}) \sim \mathcal{N}(0, I)$ for each sample $\mathbf{x}$ of the data by the regularization term. The two terms are then antagonistic, as the reconstruction term will shrink the latent projections to atomic points (similarly to regular auto-encoders), while the regularization term enforces each data point $\mathbf{x}$ to cover the whole support of the representation (see section 2.1.4.2). This behavior is due to the amortization of variational inference, that from one hand confers strong generalization properties to the system, and from the other hand can also *over-regularize* some dimensions, pushing some axis of the posterior $q(\mathbf{z}|\mathbf{x})$ to match the prior and then being deactivated during the rest of the training process. This effect, called *posterior collapse*, has been clearly highlighted, as can be shown whether by studying the decoder weights [180] or by analyzing statistics from latent projections [164, 181]. As shown by Yeung & al., *active* dimensions have typically tiny variances in average and a strong variance in the distributions' means, while *inactive* dimensions are almost perfectly fitted to the prior, projecting means near to zero and setting the variances to 1. Such regime is then called *polarized*, as the distinction between active and inactive dimensions is significant [182]. This observation clearly demonstrates the double bind of the ELBO, where the inactive units are the dimensions matching the prior.

**Bits-back coding** Latent over-pruning is particularly sensible to decoder's capacity. Bowman indeed noticed that, when using an auto-regressive decoder, the information given by the latent space was totally bypassed, such that almost every dimension was matching the prior [161]. Chen & al. proposed to explain this strange effect using a *bits-back coding* interpretation of the ELBO (2.14). Remembering the formulation of Bit-Back code length (see section 2.2.3)

$$ \begin{aligned} \mathcal{C}_{\text{BitsBack}} &= \mathbb{E}_{X \sim p(\mathbf{x}), Z \sim q(\mathbf{z}|\mathbf{x})}[\log q(\mathbf{z}|\mathbf{x}) - \log p(\mathbf{z}) - \log p(\mathbf{z}|\mathbf{x})] \\ &= \mathbb{H}[p(\mathbf{x})] + \mathbb{E}_{X \sim p(\mathbf{x})}[D_{KL}[q(\mathbf{z}|\mathbf{x})||p(\mathbf{z}|\mathbf{x})]] \end{aligned} $$

the extra-cost is then $D_{KL}[q(\mathbf{z}|\mathbf{x})||p(\mathbf{z}|\mathbf{x})]$, that can be easily avoided by the system by taking $p(\mathbf{z}|\mathbf{x}) = p(\mathbf{z})$, such that $q(\mathbf{z}|\mathbf{x}) = p(\mathbf{z})$. This

<figure><img src="image_34.png" /><figcaption>**Figure 2.32:** (a) an auto-encoder associates to a code $\mathbf{z}$ to a given example $\mathbf{x}$, but is not trained between the examples. (b) a well-balances VAE enforces the distributions to overlap, being able to find mixtures between the two distributions (c) if the recognition model $q(\mathbf{z}|\mathbf{x})$ matches the prior, the variables $\mathbf{z}$ and $\mathbf{x}$ are independent.</figcaption></figure>
principle, called *information preference*, is particularly strong at the beginning of the training, as the encoder can fit much faster the prior $p(\mathbf{z})$ than the decoder can fit the data $p(\mathbf{x})$, especially if the decoder has high capacity and thus train more slowly. This *information preference* property of the VAE thus weakens the representation learning process if the decoder has too much capacity, and then somehow weakens the adaptability of the model.

**$D_{KL}$ scheduling and data selection.** Different approaches have been developed to counter over-pruning in VAEs. The most common method to prevent this effect is to schedule the weight of regularization term $D_{KL}[q(\mathbf{z}|\mathbf{x})\|p(\mathbf{z}|\mathbf{x})]$, initiated at 0 and raising to 1 during training [183]. This *warm-up* procedure allows to prevent hasty posterior collapses, especially in hierarchical models where higher latent layers can converge really quickly to degenerated solutions [184]. Warm-up is a simple and efficient way of avoiding latent over-pruning at the beginning of the training, but unfortunately does not prevent it once the regularization term is weighted as usual. Another solution if to threshold the regularization term, such as the *free bits* constraint proposed by Kingma & al. [173] that replaces the original $D_{KL}$ by a clamped regularization term $\max\{D_{KL}[q(\mathbf{z}|\mathbf{x})\|p(\mathbf{z})], \beta\}$, such that the encoding model is forced to not exactly match the posterior.

Alternatively, over-pruning can be alleviated by skewing the amount of information given to the decoder, such that the decoder cannot recover the full data $\mathbf{x}$ without $\mathbf{z}$. This idea, called by Chen & al. *explicit information placement*, can be done by limiting the target scope of the generative distribution [185], or to train the encoder on the future of a sequence rather than its past in the case of sequential data [186]. Alternatively, Yeung & al. propose a latent gating mechanism that automatically selects the latent dimensions used for each example [164], enforcing the system to use all of its latent dimensions.

**Aggregated posterior matching.** Over-pruning occurs mainly because of the term $D_{KL}[q(\mathbf{z}|\mathbf{x})\|p(\mathbf{z})]$, that tries to fit the approximated posterior distribution to match the prior *for each data sample x*, then over-regularizing the representation *z*. This can be formally expressed by another decomposition of the ELBO, aiming to exhibit the information loss caused by regularizing $q(\mathbf{z}|\mathbf{x})$ instead of $q(\mathbf{z})$ [187]. In order to express the aggregated posterior $q(\mathbf{z})$ from $q(\mathbf{z}|\mathbf{x})$, Hoffman & al. propose to express the index $n$ as a random variable, that according to training procedures are drawn from a uniform distribution $q(n) = p(n) = 1/N$
$$ \mathcal{L}[q] = \underbrace{\left[ \frac{1}{N} \sum_{n=1}^{N} \mathbb{E}_{q(z_n|x_n)}[\log p(x_n|z_n)] \right]}_{\text{average reconstruction}} - \underbrace{(\log N - \mathbb{E}_{q(z)}[\mathbb{H}[q(n|z)]]))}_{\text{index-code mutual information}} \qquad (2.57) $$

$$ - \underbrace{D_{KL}[(\mathbf{z}||p(\mathbf{z}))]}_{\text{marginal divergence}} \qquad (2.58) $$

where the term $(\log N - \mathbb{E}_{q(z)}[\mathbb{H}[q(n|z)]]$ acts as a regularizer, enforcing the variational distribution $q(\mathbf{z}|\mathbf{x})$ to overlap between $n$. While this term regularizes the representation, it also prevents the full model to maximize the ELBO. This decomposition shows an inner trade-off of amortized variational inference: while conditioning the approximating distribution $q$ on the data $\mathbf{x}$ allows an explicit variational posterior $q(\mathbf{z}|\mathbf{x})$, it also prevents the aggregated posterior $q(\mathbf{z})$ to match the prior $p(\mathbf{z})$. Some authors have then been investigating methods to match the aggregated posterior $q(\mathbf{x})$ to the prior $q(\mathbf{x})$ (see fig. 2.33). Aggregated posterior regularization can also be justified from an optimal transport perspective, as overlapping projections $q(\mathbf{z}|\mathbf{x})$ inevitably results in producing blurry samples (see section 2.1.3.3).

<figure><img src="image_35.png" /></figure>

$D_{KL}[q(\mathbf{z}|\mathbf{x})||p(\mathbf{z})]$

<figure><img src="image_36.png" /></figure>

$D_{KL}[q(\mathbf{z})||p(\mathbf{z})]$

<caption>**Figure 2.33:** Matching the aggregated prior $q(\mathbf{z})$ instead of $q(\mathbf{z}|\mathbf{x})$ allows to prevent over-pruning, as its matching the overall distribution rather than the posterior.</caption>

As recovering the aggregated posterior $p(\mathbf{x})$ is generally intractable, we have to use implicit distribution matching $\mathcal{D}[p, q]$ to replace the $D_{KL}$ term in the ELBO 2.18. Moreover, implicit distribution matching approaches also allow arbitrary priors $p(\mathbf{z})$, providing that we can sample from it. As we saw section 2.1.3.2, MMD provides an efficient way to estimate the divergence of implicit distributions [75], as proposed by the Wasserstein Auto-Encoder (WAE) [84], or leveraging Sinkhorn distances, computing the transport cost between the two distributions [188]. Adversarial Auto-Encoders, proposed by Makhzani & al., rather propose to use density ratio estimation (see sec. 2.19) to address the classification problem

$$ r(\mathbf{x}) = \frac{p(\mathbf{x})}{q(\mathbf{x})} = \frac{p(\mathbf{x}|y=1)}{q(\mathbf{x}|y=0)} = \frac{\mathcal{D}(\mathbf{x})}{1 - \mathcal{D}(\mathbf{x})} $$

where $\mathcal{D}(\mathbf{x})$ is the output of a discriminator. Based on a similar intuition, Mescheder & al. proposed an unified Adversarial Variational Bayes (AVB) framework, using adversarial methods to rather match the full joint distributions $p(\mathbf{x}, \mathbf{z})$ and $q(\mathbf{x}, \mathbf{z})$; however, as underlined by Bousquet & al., the actual formulation of AVB does not implicitly perform aggregated posterior. While using implicit distribution methods allows us to match the marginal distribution $q(\mathbf{z})$ to the prior, and are well-defined from an optimal transport perspective, they can yet provide degenerated representation as it encourages disjoint encoding distributions $q(\mathbf{z}|\mathbf{x})$, as shown by (2.57). Thus, some approaches propose to optimize both $D_{KL}[q(\mathbf{z}|\mathbf{x})||p(\mathbf{z})]$ and $D_{KL}[q(\mathbf{z})||p(\mathbf{z})]$, hence maximizing the mutual information $\mathbb{I}[\mathbf{x}, \mathbf{z}]$ [189, 190].
**Prior learning.** Alternative approaches rather suggest to provide richer priors, assuming that the prior $p(z) \sim N(0, I)$ unavoidably yield to over-regularized representations. Tomczak & al. propose *variational mixtures of posterior priors* (VAMP), formulating the prior as a mixture of normal distributions inferred from pseudo-inputs learned with an auxiliary distribution. Interestingly, the pseudo-inputs inferred by VAMP priors implicitly represent existing patterns in the dataset, helping to perform unsupervised clustering in the latent representation. Multi-layer latent approaches, such as Deep Latent Gaussian Models [160], an also be thought as learning priors on the lower levels [32], while keeping the regularization term on the top-level. Similarly, Klushyn & al. propose to use a hierarchical prior model learning to fit the aggregated posterior $q(z)$. Some methods also propose to automatically adjust the number of latent dimensions, keeping the optimal dimensionality that is needed to maximize the reconstructions. Finally, the infinite VAEs proposed by Abbasnejad & al. propose to dynamically create single dimension VAEs once a sufficient number of data instances is not represented, using random affectation variables drawn from a Dirichlet process [191].

### 2.3.2.4 Code efficiency and disentanglement

**Disentangling and $D_{KL}$ weighting** The natural disentangling properties of VAE raised an important interest in the machine-learning community, explaining the popularity of this model. The concept of *disentangling* relates to the amount of correlation between the axis of the representation, assuming that having independent axis means that we retrieve independent factors of variation. Disentangling can thus be related to *interpretability*, as an interpretable representation would have dimensions corresponding to independent human-understandable factors of variation [192]. However, the actual origin of the witnessed tendency is hard to identify precisely. From a geometric perspective, if we want two axis of variation to be independent, these axis should be *orthogonal* in the representation, and even correspond to its axis (see fig. 2.34).

The disentangling properties of existing methods such as probabilistic PCA, robust PCA [180] or non-linear ICA [193] have been investigated, all of these methods enforcing the orthogonality of representation. Rowlin & al. recently demonstrated the role of the stochastic part of the reconstruction loss, promoting local orthogonality of the decoder. More specifically, Higgins & al., underlined the role of the regularization term $D_{KL}[q(z|x)\|p(z)]$ in disentangling, investigating the alternative objective [194]

$$ \mathcal{L}_{\beta} \mathcal{E}_{q(z)}[\log p(x|z)] - \beta D_{KL}[q(z|x)\|p(z)] \quad (2.59) $$

and showing that increasing the $\beta$ parameter was enhancing the *disentangling* properties of the representation, recovering factors of variation from image such as translation, rotation, and shape. Indeed, in VAEs, the prior is defined as a mean-field normal distribution $p(z)$, such that each dimension $z_i$ of the distribution is independent. Increasing the weight of $D_{KL}[q(z|x)\|p(z)]$ then enforces the variational distributions

<figure><img src="image_37.png" /><figcaption>**Figure 2.34:** Disentangling: given two factors of variation (running and color), here neither axis are *covariant* with axis of latent space, and are not *orthogonal* between them.</figcaption></figure>
to match the prior, and then to be factorized among these dimensions. However, as shown in previous sections, increasing $\beta$ also favors posterior-collapse and brings more distortion to the data, then blurring the reconstructions. Indeed, the objective 2.59 is a reformulation of objective (2.36), such that increasing $\beta$ lowers the representation rate, but increases the distortion. Alternative methods then propose to rather regularize the total covariance of the aggregated posterior [150]

$$TC(z) = KL[q(z) \parallel \prod_{i} q_i(z_i)] \qquad (2.60)$$

that enforces the factorization of the aggregated posterior $q(z)$ without increasing over-pruning. We can also increase mutual information between $x$ and $z$ by reducing the pair-wise distances between samples in the encoding space, as proposed in [195], or directly approximating the mutual information with optimal transport methods [196]. However, retrieving independent axis does not ensure their interpretability, as independence is rotation-invariant and may fail to extract accurate axis of variation from the data.

**Disentangling and Riemannian geometry.** Further insights on the disentangling properties of AEVB can be given under the scope of *Riemannian geometry*, bridging Bayesian learning and differential geometry. As summarized previously, *disentanglement* is closely linked to the idea of orthogonality, providing that the axis of the representation correspond to independent factor of variations underlying the data. In other words, an *interpolation* in the latent space should then correspond to a *flat* trajectory on the underlying data manifold. The quality of interpolations in the latent space are thus related to how much the representation has caught the inner manifold metrics, that in practice is rather unsatisfying for complex datasets.

Riemannian geometry can explain why the orthogonality of latent axis not a sufficient criterion for extracting disentangled representations. Remembering the definition of a manifold (see sec. 2.2.1), each point lying on the manifold defines a tangent space, with an infinity of possible orthogonal coordinates. According to the manifold a given metrics, the optimal or *parallel* translation between these two points is given by its geodesics, that is a property of the manifold. These inner metrics can be obtained with the intrinsic riemannian metrics, that describe the manifold curvature. Then, retrieving the ideal translation between two points in the latent space by making a linear interpolation imply to have extracted *flat* coordinates of manifold, providing an additional criterion for latenet disentanglement (see fig. 2.35). Formulating the data output $x = f(x) + \sigma(z) \odot \epsilon$, the corresponding Riemannian metrics tensor can be obtained by [197]

$$M_z = (J_z^{\mu})^{\top} (J_z^{\mu}) + (J_z^{\sigma})^{\top} (J_z^{\sigma})$$

that allows us to get an estimation of the Riemannian metrics defined by the latent representation. Geodesic interpolation between two

<figure><img src="image_38.png" /><figcaption>**Figure 2.35:** Evaluating the true underlying factors of a given dataset implies not only to extract the exact sub-manifold from the original space, but also to infer its underlying metrics.</figcaption></figure>
points can then be obtained by minimizing the curve length in the manifold, as proposed by Shao & al. [198] [199] using finite-time differences and arc length optimization. Geodesic interpolations have been shown to qualitatively improving translations between two points, indicating that the latent space coordinates are not flat respectively to the Riemannian metrics of the underlying manifold. To our knowledge, adding a flatness criterion to the ELBO has not been yet investigated in the literature, even if this idea is closely related to contrastive auto-encoders 2.45.

Riemannian approaches to variational auto-encoders then give consequent conclusions on the inner behavior of these algorithms, underlining the role of differential geometry on the emerging properties of the latent representation. Recent approaches then investigate this topic further to define non-Euclidean latent spaces such as spheres or toruses [200], or investigating hyperbolic geometries [201, 202]. Riemannian metrics have also been investigated to regularize the latent space in a supervised way by Hadjeres & al., enforcing latent geometrical properties to match with user-defined attribute functions [203].

**Multi-object approaches.** Disentangling in VAEs has been mostly studied in multi-object recognition, where several objects are depicted in a single image. Multi-object recognition is a seminal machine-learning task, challenging the model abilities to understand *compositionality* of images. While very efficient methods have been developed using deep convolutional networks architectures [204] [205], these methods still do not provide interpretable representations such as the one provided by AEVB. Chazan & al. proposed *deep clustering* algorithm for object recognition by pre-training multiple variational auto-encoders on every object, and then use these models as a mixture of experts based on the reconstruction score [206]. However, pre-training specific auto-encoders for each object can be tedious, and for many objects the reconstruction procedure can be too slow. Alternatively, Greff & al. propose a multi-object recognition model using VAEs by splitting the latent spaces in subsets that are fed to multiple decoders, recomposing the original image by generating a mask in addition to the image [207]. Multi-object recognition can also be processed in a sequential way, generating the salient objects one after the other on a canvas, as done by the DRAW model proposed by Gregor & al. [208]. Similarly, the *Attend, Infer, Repeat* procedure proposed by Eslami & al is based on a single structured decoder, that is used sequentially to recompose all the objects present in a scene until the inference model stops the procedure. Latent space is then explicitly disentangled to infer separately spatial and object properties, then enforcing its interpretability [209]. Explicit latent disentanglement can then be enforced by designing the generation process, and is certainly an interesting perspective for variational auto-encoding systems.
<figure><img src="image_39.png" /><figcaption>**Figure 2.36:** Relaxation of discrete probability function: (a) discrete probability function $p(x)$ (b) obtain cumulative distribution function $\int p(x)dx$ (c) relaxation by approximating the CDF with a continuous distribution. Provided its tractable reciprocal, this distribution can be sampled.</figcaption></figure>

### 2.3.3 Discrete latent spaces

Finally, we will address the problem of discrete latent spaces, that is a challenging objective as back-propagation through discrete variables is generally not possible [210]. Indeed, discrete variational distributions $q_{\phi}(z|x)$ are not re-parametrizable, and approximating the gradient with a Monte-Carlo estimator is often too noisy to ensure convergence. Yet, extracting discrete representations could be useful for several applications such as unsupervised classification or topic detection [211]. Discrete representations can also be a very efficient method to reduce over-pruning when using auto-regressive decoders, as shown by van den Oord & al [212], or to achieve latent input selection by inferring a mask [213].

**Continuous relaxation.** Discrete latent variables are usually modeled using Categorical distributions, that can be thought as a generalization of Bernoulli distributions with $k > 2$ possible issues. Categorical distributions are typically used at the end of a classification network, jointly with a final Softmax non-linearity layer, and directly back-propagated with cross-entropy estimators. However, categorical distribution do not allow re-parametrization trick and estimation of its gradient with REINFORCE converge slowly. A solution to ease the back-propagation through discrete random variables is continuous relaxation, that consists in projecting the discrete representation in a continuous space using a smoothing transformation. Categorical distribution can be smoothed using Gumbel distribution, such that sampling from Categorical distribution can be done using samples from $g_i \sim \text{Gumbel}(0, 1)$ [214, 215]

$$ 
\begin{aligned}
\mathbf{z} &= \text{onehot}(\operatorname{argmax}_i[g_i + \pi_i]) \\
&= \frac{\exp(\log \pi_i + g_i)/\tau}{\sum_{i=1}^k \exp(\log \pi_i + g_i)/\tau}
\end{aligned}
$$

by approximating the argmax function with a Softmax. Such variables are distributed as Gumbel-Softmax distribution (also called Concrete), whose probability density function is tractable and can be used to approximate Categorical variables. The $\tau$ parameter is called the temperature of the distribution, such as it tends to be atomic when $\tau \to 0$, but equals a uniform distribution as $\tau \to +\infty$. As the Gumbel-Softmax distribution is obtained from a tractable inverse cumulative distribution, it can be applied the reparametrization trick, and then be efficiently integrated into the AEVB framework. The gradient of a Categorical distribution can then be approximated by using the one from a Gumbel-Softmax, a method called straight-through Gumbel estimator, that tends to the real gradient when $\tau \to 0$. While the Gumbel-Softmax trick can provide efficient approximations to Categorical distributions, they can be too smooth to accurately perform discrete sampling. Rolfe & al., followed by Vahdat & al., rather propose spike-and-exponential smoothing transformations, that only activates the unit if the real threshold of the corresponding latent unit is passed. Xu & al. also propose to define the simplicial probabilities by taking planar projections
of samples obtained from a Von-Mises distribution, whose support is defined on an hyper-sphere $S^K$ [216].

**Unsupervised and semi-supervised learning** Some approaches rather extract discrete information from the continuous latent space by quantifying it, then labelling different zones of the latent spaces in a similar way to clustering. By example, the GM-VAE proposed by Dilokthanakul & al. [217] uses a Gaussian-Mixture prior, where the index of the corresponding mixture is inferred from the latent space. Similarly, van den Oord & al. propose to use vector quantization on the latent space to extract a discrete dictionary of latent centroids, using straight-through gradient estimation to back-propagate in the encoder. In that case the regularization constraint then become a $\ell^2$ norm between latent projections and the closest embedding vector, then preventing over-pruning but enforcing the sparsity of the representation.

Extracting discrete information from continuous latent-space can also be performed in a supervised manner, adding a discriminator on the top of the latent space and trained on corresponding symbols y. Classification can then be processed afterwards, such that the latent space provide a compressed representation to the discriminator. More interestingly, the auto-encoder and the discriminator can be jointly processed, such that the representation is also influenced by the gradient information coming from the discriminator. An advantage over traditional supervised approach is that, in case of incomplete label information, the auto-encoder can still be trained by inferring the corresponding label. Such methods have been proposed by Kingma & al. as semi-supervised, and provide a flexible method for task-oriented representation learning (see fig. 2.37). Probability densities of data x, labels y and latent variables z can be described

$$p(y) = \text{Cat}(y|\pi(\mathbf{z})); \quad p(\mathbf{z}) = \mathcal{N}(\mathbf{z}|\mathbf{0}, \mathbf{I}); \quad p(\mathbf{x}|y, \mathbf{z}, \theta) = f(\mathbf{x}; y, \mathbf{z}, \theta)$$

where $p(y)$ is a prior over discrete variables. Two different bounds are then derived, depending on if the availability of the label information. If the label is available, the lower-bound is the direct application of (2.14) to the generative model:

$$\mathcal{L}(\mathbf{x}, y) = \mathbb{E}_{q_{\phi}(\mathbf{z}|\mathbf{x}, y)}[\log p_{\theta}(\mathbf{x}|y, \mathbf{z}) + \log p_{\theta}(y) + \log p(\mathbf{z}) - \log q_{\phi}(\mathbf{z}|\mathbf{x}, y)]$$

that is generally added an explicit classification cost $\alpha \cdot \mathbb{E}_{p_{\theta}(\mathbf{x}|y)}[-\log q_{\phi}(y, \mathbf{x})]$. If the label information is missing, we compute the expectation over the full inference distribution $q_{\phi}(y, \mathbf{z}|\mathbf{x}, y)$

$$\mathcal{U}(\mathbf{x}) = \mathbb{E}_{q_{\phi}(y, \mathbf{z}|\mathbf{x}, y)}[\log p_{\theta}(\mathbf{x}|y, \mathbf{z}) + \log p_{\theta}(y) + \log p(\mathbf{z}) - \log q_{\phi}(y, \mathbf{z}|\mathbf{x})]$$

where the expectation respectively to y can be used by summing the possible values of y, y being discrete. As the influence of the discriminator can be important and remove the unsupervised features of the latent space, the discrimination step can be proceeded on the second layer of a multi-layer VAE as proposed by Kingma & al., or using auxiliary variables [218]. Bouchacourt also propose a disentangling

<figure><img src="image_40.png" /><figcaption>**Figure 2.37: Semi-supervised learning : a** supervised task is perform in the latent space, and the extracted label is given to the decoder for class-dependant generation.</figcaption></figure>
method by grouping observations, enforcing the same support for data
items with common label properties [219].
Learning spectral representations and audio regularization strategies

In the previous chapter, we introduced two different frameworks for the extraction of invertible representations: Bayesian latent models, based on modeling the data with a generative distribution $p_{\theta}(x|z)$ that is conditioned on inferred latent variables $z$, and methods coming from representation learning, leveraging high-capacity function approximators to extract features from the data according to a suitable loss function. We saw that these two frameworks could be beneficially combined to obtain a hybrid method for generative models, called Auto-Encoding Variational Bayes (AEVB), that is based on both Bayesian statistical inference and neural networks to model both the generative model $p_{\theta}(x|z)$ and the approximated posterior $q_{\theta}(z|x)$.

In this chapter, we will invest AEVB as an generative model for audio signals. We focus on AEVB because it combines two different processes: the **generation process** $p_{\theta}(x|z)$, that generates new samples from the analysed dataset conditioned on a latent representation $z$, and an **inference process** $q_{\theta}(z|x)$, that extracts this latent representation in an unsupervised way. This PhD thesis focuses on this idea of using this framework to mirror audio *analysis-synthesis* methods, that are seminal DSP techniques based on two similar steps: a *synthesis* process, generating a signal with respect to a set of parameters, and an *analysis* process, that extracts the generative parameters from an incoming audio signal (see 3.1). As an introduction, we will shortly summarize some analysis and synthesis methods in the DSP domain, linking them to machine-learning methods, and then reviewing the actual analysis-synthesis techniques and see how our work takes place in the current state of the art.

**Analysis and inference**

Analysis methods, in the domain of audio DSP, are commonly based on the extraction of higher-order features from the signal, aiming to reflect structural, acoustical or musical properties of the signal. These high-order features, that can also be called *descriptors* in the audio DSP literature, can be based on the extraction of features from intermediate transforms, that are generally derived from *spectral* transformations (see 3.10), to reflect perceptual properties of the signal: statistical or perceptual descriptors such that fundamental frequency, spectral centroid, and many others. Contrary to the features extracted with standard machine-learning, that are not domain-dependent, audio descriptors are generally crafted by researchers to reflect these qualities, linked to perceptual attributes that can be confirmed with psycho acoustical experiments (see sec. 3.4.1). While these descriptors are generally short-term, some analysis methods rather target to model the *dynamics* of the signal, whether by tracking the evolution of these descriptors through time [220] or to model its temporal structure (for speech analysis [221]), energy profiles

<figure><img src="image_41.png" /><figcaption>**Figure 3.1:** The inference/generation processes of AEVB can be thought as mirroring a DSP analysis/synthesis system, extracting invertible features from a given corpus of audio signals.</figcaption></figure>
<figure><img src="image_42.png" /><figcaption>**Figure 3.2:** Most MIR approaches extract acoustical features from the signal, and/or from signal representations to extract higher-level features, as could be used here from chord detection.</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>**Figure 3.3:** Classical audio synthesis methods are based on a physical (analog) or algorithmic (digital) structure, some parameters being offered to the user to shape the sound.</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>**Figure 3.4:** Even if a distribution $p(x)$ perfectly represents the modelled data, the lack of parameters makes it unusable.</figcaption></figure>

[222] or even audio watermarking [223]). The domain of music information retrieval (MIR) [224–226] focuses on the extraction of musical features from an audio signal, that are typically high-order descriptors mixing temporal and spectral dimensions such as pitch and chord tracking [227, 228], beat extraction and tempo estimation [229–232], structural segmentation using similarity matrices [233] or spectral clustering [234–236] and more high-level tasks like emotion recognition [237, 238], genre recognition [239, 240] and music recommendation [241–243].

As MIR is focused on high-level tasks that often involve several descriptors on multiple temporal scopes, the resort to machine-learning techniques for feature extraction can be relevantly introduced to support audio analysis frameworks. The recent raise of machine learning techniques based on convolutional neural networks motivated their direct application on spectrograms, by example on classification/regression tasks [244–247]. However, the lack of interpretability and the data-centred nature of these algorithms aroused a mixed reception of the MIR community, despite a substantial increase in the results. Therefore, hybrid methods between MIR and machine-learning have since been investigated, combining best of both worlds to address complex tasks such as source separation [248–251], transcription [252] or voice recognition [253].

**Synthesis and generation** Audio synthesis is one of the most seminal domain of audio DSP since its advent [254, 255], partly due to the early use of analogical and then numerical devices in music creation. Digital synthesis appeared almost simultaneously with computers, and is still a flourishing domain of research led by many different approaches: additive and subtractive synthesis [256, 257], FM synthesis [3], sampling, granular synthesis [258], physical modeling [Avanzini01controllingmaterial, 259–261] and many other, that are now widely integrated in modern practices of music production and performance. Digital synthesis can be loosely defined as a generating system whose state is defined by a set of parameters, that conditions the structure of the generated signal (see fig 3.3). Parameters of digital synthesis techniques are often tied to structural properties of the generation model, such that the effect of each generative parameter is explicitly known. Some parameters are then offered to the performer or producer for direct interaction, providing controls for the generation model. Therefore, the identity of synthesis method is not only defined by the range of sounds it can produce, but also by the controls it gives to the human user.

Machine-learning based generative models are slightly more complicated, as their development come from different origins and initially bypassed this interaction aspect. As we summarized last chapters, generation in the machine-learning domain was rather formulated as the modeling an underlying stochastic process $p(x)$ under a finite set
of observations {$x_1x_2...x_n$}$_{i=1}^N$. However, this distribution can be unsatisfactory for musical purposes as it only allows *sampling*, and do not allow explicit control (see fig. 3.4). Therefore, we saw in the last section a way to define an extended model $p(x, z)$ with parameters or additional inputs $z$, such that $z$ can shape the density of $x$ and then control the generation. However, the number of machine-learning generative method specifically developed for generation is not as flourishing than in digital audio synthesis, and the generative models can be gathered in four families : *auto-regressive* [262, 263], *auto-encoder* approaches (see sec. 2.2.3), *generative adversarial networks* [23, 264, 265]), and recently flow-based models [171, 266, 267]. However, only auto-encoding and flow-based approaches are invertible, and can be thus used for our purpose.

Therefore, a seminal difference between digital audio synthesis and machine-learning generative models is that, in the first case, *the model shapes the generated sound*, while in the second *data is shaping the sound*. We think that these solutions are equally good, and that this complementarity can significantly help human creativity. However, we think machine-learning based approaches have still to catch up digital audio synthesis on various aspects, especially in terms of perceptual qualities, usability and interaction. While we investigate the second chapter 5, the two following sections will propose diverse ameliorations of AEVB generative models for sound synthesis, tackling two different aspects : latent space regularizations, and temporal modeling.

**Analysis-synthesis** The complementarity between analysis and synthesis also exists in the signal processing domain, such that powerful analysis methods do not generally allow direct regeneration of the signal, and that conversely digital synthesis algorithms do not provide any way to infer corresponding parameters from a given sound. This duality then motivated the development of frameworks allowing both analysis and synthesis, providing bijective methods between audio and parameters domains. These models, called *analysis-synthesis* methods, often rely on *invertible* transformations that both extract features from the signal and allow the regeneration of original data. The obtained features can then be used whether for analysis purposes (such as spectral transforms, that ground most of perceptual descriptors) or sound transformations, modifying the incoming sound in the obtained representation. *Phase-Vocoder* algorithm is a seminal example of such method, using the invertibility of Short-Term Fourier Transform (STFT, see sec. 3.10) to transform existing sounds (such as transposition, time stretching, or partial tracking), or for direct generation by designing sounds directly in the spectral domain. Other analysis-synthesis methods, such as ESPRIT [268] or cepstrum-based representations [269], also provides bijective representations that are based on different modeling properties. However, a drawback of these methods is that their general invertibility enforces the bijectivity of the transform, such that the dimensionality of data and representation domains are equivalent. Hence, these methods are often used as intermediary signal
representations for higher-level analysis / generation methods.

Hence, we propose to use the AEVB framework as an analysis-synthesis method for digital audio synthesis, combining modern probabilistic generative and latent models to extract high-level invertible audio features. This model gains from other generative models such as GANs by also learning an inference process on the latent representation, thus allowing the system to structure the space without supervision and to retrieve the parameters corresponding to an incoming data. We also think that this framework provides a complementary approach to classical DSP analysis-synthesis methods, targeting to extract meaningful features from a restricted set of data rather than requiring the global invertibility of the representation. The relaxation of this requirement thus allows us to reduce the input dimensionality, rather focusing to model a dataset with a finite amount of samples.

This section is divided in three parts. First, we will first provide some baselines of performances obtained with diverse AEVB systems, that we train with the magnitude of audio spectra. These baselines, developed during the early phase of the PhD, were at our knowledge the first application of AEVB to audio signals. Then, we will present two latent regularization strategies that we developed specifically for audio processing and generation : a first method based on external symbolic information (such as pitch, octave and dynamics) using the latent space as a *translation space* between signal and symbol domains, and a second method based on perceptual constraints, obtained from psycho-acoustical measurements.

**3.1 Spectral representations of audio signals**

In this section, we evaluate the performances obtained with diverse AEVB systems when trained on audio signals. However, learning directly audio waveforms is not straightforward, containing short and long-term temporal dependencies that would involve higher-order time series models (see section 4). Hence, we rather perform training on *magnitude spectra*, that allows us to train on stationary and parsimonious data that more easily reflects structural perceptual attributes of audio signals. As such transforms are widely used in the DSP community and implemented in most audio software, it also allows to easily develop a real-time audio engine to generate from the obtained models (see section 5), hence offering a good compromise between usability and complexity.

**3.1.1 Invertible spectral representations of audio
signals**

**Raw representations.** The *raw* representation of an audio signal is the direct description of the corresponding physical pressure wave, that can be represented using one or several signals $s(t)$, or by an acoustical field $s(t, x)$, $x$ being the position in space. This signal can be
<figure><img src="image_45.png" /><figcaption>**Figure 3.5:** Diagrams showing the three models used in this section. (left) simple AEVB system, trained on audio reconstruction (middle) AEVB with translating latent space, using audio x and symbol data y 3.3 (right) perceptual regularization of the latent space, using dissimilarity experiments of instruments sounds from audio perception studies 3.4</figcaption></figure>

recorded at one or several positions in space using a transducer, or di-
rectly generated with a synthesis engine. In the numerical domain, this
wave is sampled both in time and amplitude in order to be encoded
numerically. The sampling frequency, corresponding to the number of
samples recorded by seconds, thresholds the maximum frequency that
we are able to render by the Shannon-Nyquist equation $f_s = 2f_{max}$
(see fig 3.6). The standard value $f_s$ for CD audio recordings, that is
44100 Hz, thus means that we can restore the original frequency con-
tent of the signal up to 22050 Hz, that is the standard upper perceptive
threshold of human audition. The number of bits used to encode the
numerical value of each sample, called the *bit rate*, defines the range
of *dynamics* that the recording is able to render. This bit rate is usually
expressed in decibels (dB), such that a 16 bits resolution allows a 90dB
range.

Individual raw waveforms are then 1-dimensional discrete time se-
ries arbitrarily, that are defined onto a range [-1; 1]. These signals
can be modelled through time with several analysis methods such as
auto-regressive, integrated, or moving-average models (see sec 4.1).
However, the high temporal density and non-stationary nature of these
signals make most of these methods cumbersome on the long-term,
such that their direct use for audio feature extraction is difficult. More-
over, the temporal scope needed to observe some high-level audio
features (from beat detection to genre recognition) can be very large,
and then make their extraction from time series dauntingly hard.
Modern machine learning techniques, based on the modeling of com-
plex auto-regressive distributions p(x<sub>t</sub>|x<sub>t-1</sub>, h) such as Hidden Markov
Models, are also trained with difficulties, then requiring more subtle
approaches (see chapter 4). As most time series analysis are defined

<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /><figcaption>**Figure 3.6:** (a) analog signal x(t) (b) sampled series x[n] (c) aliasing: the discretization of the waveform introduces an infinite amount of possible frequencies $f_{alias} = f_s - nf - \pi$</figcaption></figure>
<figure><img src="image_49.png" /><figcaption>**Figure 3.7:** Decomposition of a square signal over the 15 first spectral components.</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>**Figure 3.8:** A complex number can be whether expressed in Cartesian coordinates (a, b) or polar coordinates (r, θ).</figcaption></figure>

in the real domain, some machine-learning algorithms rather use the underlying binarization of numerical signals to define a classification problem, sparsifying the data and speeding up the convergence [262]. μ-law is a common non-linear encoding scheme that encodes more precisely high amplitude values, thus reducing the signal/noise ratio

$$F(x) = \operatorname{sgn}(x) \frac{\log 1 + \mu|x|}{\log(1 + \mu)}$$

However, this methods also significantly increase the dimensionality of the input space, transforming each sample into a $\mathbb{R}^d$ vector whose dimensionality scales exponentially with the number of bits (16 bits amounts to 65 536). Therefore, some downsampling is required, altering the quality of the audio signal. Raw waveform learning is then one of the most difficult machine-learning challenge, is cannot be tackled on sufficiently large time scales with standard AEVB frameworks. Hence, we rather resort to alternative representations with smoother properties, defined in the spectral domain, that significantly ease the training of the proposed models on audio information.

**Fourier transform and time-frequency trade-off** While raw audio waveforms is the closest representation of the corresponding acoustical phenomenon, it is difficult to model with the proposed techniques. Furthermore, this representation does not provide explicit information about the frequency content of the signal, that provides a more direct interpretation of several properties of human perception such as pitch, harmonicity, or timbre-related attributes. Given a continuous signal $x(t)$, its frequency content can be extracted using the Fourier transform, a representation defined in the complex domain $X(f) \in \mathbb{C}$

$$X(f) = \int_{-\infty}^{\infty} x(t)e^{-2\pi x f} dx \qquad (3.1)$$

that can be thought as the projection of the signal onto an infinite function space, whose basis are the complex exponentials $\exp ix = \cos x + i \sin x$ (see fig. 3.7). The Fourier transform then describes the signal as a function of frequencies $f$ instead of time $t$, and can be intuited as the convolution of the signal with a complex exponential of frequency $f$. As $X(f)$ is defined in the complex plane $X(f) \in \mathbb{C}$, it is composed by a real part $a$ and an imaginary part $b$, such that $X(f) = a + bi$. The complex frequencies $X(f)$ can then also be expressed using polar coordinates $r = \sqrt{x^2 + y^2}$ and $\theta = \operatorname{atan}\frac{y}{x}$ (see fig. 3.8)

$$X(f) = re^{i\theta}$$

where $r$ is the radius and $\theta$ is the angle of $X(f)$. Expressing spectral components in polar coordinates provides a more explicit interpretation of $X(f)$: the radius is the magnitude of the spectral component $X(f)$, while the angle $\theta$ corresponds to its phase, the relative position of the component relatively to its period (see fig. 3.9). This decomposition is widely used in audio signal processing, as it allows us to detach the magnitude of each spectral component from its temporal information. Furthermore, this transformation is entirely invertible, such that we
can recover the original temporal $x(t)$ from $X(f)$ by performing inverse Fourier transform

$$x(t) = \int_{-\infty}^{\infty} X(f)e^{2\pi x f} df \qquad (3.2)$$

who is defined on the real domain $\mathbb{R}$, as the Fourier transform of a real signal is an even function.

The expression of Fourier transform defined in 3.1 allows us to get an invertible and decomposable representation of a signal $x(t)$ in the frequency space $f \in \mathbb{C}$, and is then a seminal transformation in signal processing. However, in the numerical domain, the signal is a discrete time series $x[n]$ sampled at rate $f_s$, such that transforms 3.1 and 3.2 are not applicable. The discrete Fourier transform (DFT), that is an extension of continuous Fourier decomposition to discrete signals $x[n] = \{x_0, x_1...x_{N-1}\}$, can be expressed as

<figure><img src="image_51.png" /><figcaption>Figure 3.9: (top) signal $x[n]$ (middle) magnitude spectrum $|x[f]|$, with the aliasing $k > f_s/2$ and (down) the phase $\angle X[f]$.</figcaption></figure>

$$X[k] = \sum_{i=0}^{N-1} x[n] e^{-\frac{i2\pi}{N}kn} x[n] = \sum_{k=0}^{N-1} X[k] e^{-\frac{i2\pi}{N}kn} \qquad (3.3)$$

that defines a finite complex vector $X \in \mathbb{C}^N$, that is periodic with period $N$ because of the aliasing caused by sampling. The discrete frequencies $k$ can be converted back in Hertz using the identity $f_k = f_s/k$, as the period $N$ of the transform is the sampling frequency $f_s$ of the signal.

An seminal property of spectral representations is their inherent duality between time and frequency. Indeed, the exact estimation of the frequency content requires the observation of the signal over an infinite time range, while conversely the exact invertibility of the spectral representation needs to integrate the spectrum over an infinite frequency domain. Therefore, the higher we want the estimated spectrum to be precise, the larger we will have to integrate the signal, hence discarding its dynamics. This time-frequency trade-off, also called *Gabor limit*, is seminal in spectral representation, implying that the spectral resolution of a signal is complementary to the time range of the observation. The hypothesis underlying this duality is the *stationarity* of the signal, assuming that its spectral content stays the same over an infinite time. In addition to being difficult in practical applications, the key point of audio signal analysis is also to retrieve the dynamics of these spectral features, such that this trade-off is problematic for signal analysis.

A trade-off between time and frequency resolution can be found by slicing the signal in short temporal windows, that are then used for spectral analysis. This can be obtained by the Short-Term Fourier Transform (STFT), that is then dependent of both time and frequency

$$\text{STFT}[x(t)](\tau, f) = \int_{-\infty}^{+\infty} x(t)w(t-\tau)e^{-j\omega t}dt$$

where $w(t)$ is a window function with finite time support, such that the corresponding spectral content is extracted around $t = \tau$. If the sum of
<figure><img src="image_52.png" /><figcaption>**Figure 3.10:** Windowing the signal and performing FFT on each segment allows to represent the frequency evolving through time.</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>**Figure 3.11:** Examples of wavelets. The input wave is decomposed in *wavelets*, allowing to perform multi-resolution analysis.</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>**Figure 3.12:** Discrete Cosine Transform (DCT) encodes the phase in the negative components of the spectrum. It can be understood as a DFT in a different base.</figcaption></figure>

windowing functions sums to one over time, such that $\forall t \sum_{\tau} w(t - \tau) = 1$, the original signal can be exactly recovered by inverting each frequency slice and performing the same windowing function (then called *overlap-add* reconstruction). STFT allows us to disentangle frequencies and time, and then to describe the evolution of spectral content over time (see fig 3.10). STFT then provides the baseline of most digital audio analysis techniques developed so far, and allowed the development of analysis-synthesis techniques such as the phase vocoder, providing an efficient method to both represent and generating signals [270].

**Alternative time-frequency representations.** Several others time-frequency representations based on the same principle have been developed. *Discrete Cosine Transform*, that is an alternative formulation of STFT defined on the real domain [271]. DCT allows to represent the phase by introducing the idea of *negative magnitude*, then entangling phase and magnitude. However, we found that machine-learning were struggling to learn this representation, as this disentanglement of phase and magnitude was preventing the algorithm to catch salient features.

A well-known alternative to Fourier transforms are *wavelet transforms*, that project the signal on arbitrary families of multi-resolution frequency bands. Wavelets allow to avoid the time-frequency trade-off of STFT, as the different frequency bands are defined on multiple time scales (see fig. 3.11). Contrary to DFT, that assumes a linear frequency scale, wavelet transforms allow to model an arbitrary range of frequencies, allowing to be arbitrarily precise over a specific spectral span. *Constant-Q transform* is an example of wavelet-based transform widely used in musical signal processing, defining a logarithmic frequency range centred on octave fractions, providing a precise representation for the analysis of pitched audio signals [272]. Unfortunately, this multi-resolution property makes wavelet transforms non-invertible, and then not suitable for generation. However, the recent *Non-Stationary Gabor Transform* (NSGT) alleviate this problem by proposing an invertible wavelet representation, using multi-resolution alignment of windows for accurate offline inversion [273]. This transform is used for the proposed perceptual regularization method of latent spaces (see sec. 3.28).

**Learning complex representations.** Time-frequency representations are then an efficient solution for training machine-learning methods on audio signals, providing a parsimonious representation that extract relevant features in terms of structure and human perception. However, most time-frequency are defined on a complex domain, such that the usual gradient-based methods based on derivability are not applicable. Neural networks developed specifically for complex-domain learning have been investigated [274–276], that could be combined with complex normal distributions over spectra [250, 277]; however,
we leave this investigation for future works.

We take a more straightforward solution by rather resorting to the *polar decomposition* described above, decomposing a spectral component $X(f)$ into a positive real *magnitude* vector $|X(f)| \in \mathbb{R}^+$ and a *phase* vector, bounded between $\theta(f) \in [0; 2\pi(\cdot)$. In the audio domain, the information provided by the spectrum magnitude contains most of perceptual properties, corresponding to the amplitude of all the frequency component. Furthermore, magnitude spectrum of simple signals are generally sparse, then providing a suitable input space for machine learning algorithms.

Conversely, the phase information is difficult to learn because of its high variance and its invariance through time. Indeed, phase represents the temporal part of the spectral transform, and is seminal to preserve the temporal consistency of a given sequence of spectral frames. However, as in this chapter we will perform training on individual frequency frames, direct interaction from the model's latent space will not preserve the temporal coherency of initial samples, so it seems reasonable to drop the information contained in the phase. However, as phase is important when inverting a sequence of several spectral frames, we have to leverage *phase reconstruction* algorithms to reconstruct the correct shape of the obtained generation. Phase reconstruction is a difficult task, as it implies recovering half of the information lost from the original spectrum. A standard method for phase reconstruction is the Griffin-Lim algorithm, based on iterative STFT and ISTFT passes to retrieve the phase distribution ensuring the consistency of the transformation [278, 279]. Several methods have then been proposed, such as phase gradient-heap integration proposed by Pruša and Søndergaard [280]. In this work, we use Griffin-Lim reconstruction with 30 iterations for offline examples.

## 3.2 Audio analysis-synthesis with variational auto-encoders

In this section, we will then provide a first baseline of AEVB performances when trained on spectrum magnitudes. As we propose to use AEVB as an analysis-synthesis method, our objective is twofold. First, we have to evaluate the generation capacities of the model, i.e. evaluating the reconstruction ability of the model over data samples. As usual in the machine-learning domain, this evaluation is performed on another set, called the *test set*, that is different from the train set, to measure the generalization ability of the model. We also have to study the emerging properties of the extracted representation $z$, and if the obtained space is explaining the underlying structure of the data. As we summarized last section, the flexibility of AEVB techniques provide numerous possibilities in terms of architecture, regularizations, and training criteria; we then propose to evaluate the influence of these choices on the emerging properties of the obtained latent spaces.
Furthermore, as machine-learning methods are highly dependent on
the involved data, we also evaluate these methods on three datasets,
with separated generative factors and different complexity.

**3.2.1 Datasets.**

The data-centred approach proposed by machine learning then re-
quires a particular attention on the choice of evaluation datasets, as
their design will condition the properties we will observe. Moreover,
the comparison of different models involves to be applied on the same
dataset, thus requiring reference datasets that are designed to reflect the
main challenges of the targeted task. A typical example of such dataset
is the MNIST corpus used in image processing, containing over 80,000
examples of hand-written digits [159]. This dataset still provides a very
efficient baseline for evaluation purposes, as it encompasses most of
the image processing challenges : classification of digits, style recogni-
tion, and limited variability (limited number of lines, 10 classes, and
black & white). The comparison of different machine-learning models
is then straight-forward, providing a direct feedback on how the con-
cerned models are behaving on these different factors of variations.

While some reference datasets for machine learning are established
in the image processing domain for decades, such datasets are unfor-
tunately missing in the audio processing domain. Valuable reference
datasets exist for main MIR domains such as pitch extraction, structural
segmentation or style recognition, but datasets focused on evaluating
the generative models performance are, to our knowledge, missing.
Hence, we built different datasets for the purpose of this thesis, focus-
ing on different factors of variation. The proposed datasets are split in
two types : first, we designed toy datasets using sound synthesis, allow-
ing us to focus on seminal generative factors such as pitch, harmonics,
and spectral content. Then, we evaluate the proposed models on audio
data coming from real audio recordings, then having more variance
and allowing us to evaluate the credibility of generated samples.

3.2.1.1 Toy datasets.

As we aim to study the performances of AEVB on spectral magnitudes,
we have to design our toy datasets on generative factors that are trans-
parent under this representation. The proposed toy datasets, that can
be easily generated using simple synthesizers, are very convenient to
rapidly evaluate the VAE performances over a short range of audio
factors with an arbitrary number of data samples. These toy datasets
are available here, and can be synthesized with the vschaos tool-
box.

**Additive dataset.** The toy_additive dataset has been obtained by synthesizing harmonic sounds over a range of 100 different fundamental frequencies, then focusing on the *pitch* and *harmonic* content
of the sound. Additive synthesis is based on the following generative process

$$x[n] = \sum_{l=1}^{L} \exp\left(\frac{-ld}{L}\right) \cos\left(2\pi f_0 l \frac{n}{f_s} + \phi\right) s \quad (3.4)$$

where $f_s$ is the sampling frequency of the target signal. This generator has 3 parameters: the fundamental frequency $f_0$, the number of partials $L$, and the harmonic decay $d$. The fundamental frequencies are chosen between 64 and 1024 Hz, the number of partials between 1 and 7, and the harmonic decay between 0.2 and 4. The phase is sampled individually from a normal distribution, and the sound is then enveloped with simple linear ramps of 100 milliseconds to avoid discontinuities. The total number of examples then amounts to 3500, and provides a simple dataset to grasp the VAE performances on fundamental properties of sounds.

**FM dataset.** We also provide a more complex toy dataset, toy fm, based on FM synthesis, that is a widely used technique in signal processing and digital sound synthesis. This method allows us generate sounds covering a wider space of frequencies and inharmonic relations with a limited number of parameters, based on the generative process

$$x_{\text{carrier}}[n] = \cos\left(2\pi m * f_0 \frac{n}{f_s} + \phi_c\right) \quad (3.5)$$

$$x[n] = \cos\left(\alpha x_{\text{carrier}}[n] * f_0 \frac{n}{f_s} + \phi\right) \quad (3.6)$$

where $f_0$ is the carrier frequency of the signal, $m$ the ratio of the modulating frequency, and $\alpha$ the modulation amplitude. Indeed, if $m \in \mathbb{N}$, the obtained will be harmonic, producing an amount of frequency components that is controlled by the modulation amplitude. However, for non integer values of $m$, the generator will produce complex spectra involving non harmonic partials, whose strength can be controlled by $\alpha$. We then generated a total amount of 21 250 examples, with 25 center frequencies from 60 to 600Hz, 25 FM ratios from 0 to 3, and 34 modulation amplitudes from 0 to 30.

### 3.2.1.2 Ground-truth datasets.

To evaluate the capabilities of the models on real-life examples, in this PhD we selected two datasets: a first one, acidsInstruments-ordinario consisting of orchestral instruments samples, and a second one, diva - dataset, consisting of sounds coming from presets of a commercial synthesizer, Diva*.

#### Orchestral ordinario sounds.

<footnote>* https://u-he.com/products/diva/</footnote>
<figure><img src="image_55.png" /><figcaption>**Figure 3.13:** Examples of spectra generated by the toy dataset. (top) : examples from *toy_additive*, generated from a fundamental frequency and a given number of partials, whose decay are paramatrized with the exponential decay parameter. (down) : examples from *toy FM*, that are typically much more complex than the ones from *toy_additive*, where the partials are not spread according to harmonic relationships.</figcaption></figure>
The first dataset, *acidsInstruments-ordinario*, is built from orchestral instruments samples coming from the *Studio On Line* database [281]. This sample bank offers a wide variety of instruments playing in various modes, and is thus a precious source for building datasets. However, to restrict this bank to a tractable database for our VAE, we only retrieved the samples taken from *ordinario* playing modes of the following instruments: English-Horn, French-Horn, Tenor-Trombone, Trumpet-C, Piano, Violin, Violoncello, Alto-Sax, Bassoon, Clarinet-Bb, Flute, and Oboe. The dataset contains samples from the full span of the tessitura for each instrument, and is available in four different dynamics: *pp*, *p*, *mf*, *ff*, amounting to 1885 files. The size of this dataset is reasonable, quite small compared to most image processing datasets (80000 for MNIST), seems realistic for practical applications, reflecting how the system would be have if the user wanted to train it on a custom set of sounds (such as sample banks).

**Diva dataset** The second dataset, *diva_dataset*, is built from the audio signals generated with the synthesizer *Diva* from U-He, that is a cutting-edge polyphonic virtual synthesizer widely used in music production. *Diva* is based on the simulation of a analogue synthesizers, offering many features and a generous set of parameters allowing to generate a wide diversity of different sounds. This dataset was chosen because most of its parameters are MIDI CC controllable, allowing us to automatically generate samples according to a sampling scheme, and because of the amounts of presets available online.

We manually established the correspondence between the synth and
MIDI parameters, as well as their values range and distributions. We
only kept the continuous parameters, normalizing all their values to [0,
1], the other being set to their fixed default value. We used RenderMan*
to batch-generate all the audio files by playing the note for 3 sec. and
recording for 4 sec. to capture the note release. The files are saved
in 22050Hz and 16bit floating point format. This procedure allowed
us to gather a total amount of 11,000 presets, providing a precious
ground-truth database for synthesized sounds. We use this dataset in
the chapter 4, as we more interested into the temporal evolution of
these sounds than into their spectral content.

<figure><img src="image_56.png" /><figcaption>**Figure 3.14:** Examples of spectra from
the *acidsInstruments-ordinario*
dataset.</figcaption></figure>

**3.2.2 Evaluation of generative and representational**
**properties.**

In this section, we provide preliminary results on the performances
achieved by VAE when trained on magnitude spectra. The evaluation
of generative models is a difficult task, especially in music where the
creative and perceptive dimensions of the generated sounds are sub-
ject dependant. While we let the *the creative* evaluation to a further
work (yet initiated, see chapter 5), we follow the *function-behavior on-
tology of Gero & Kannegiesser [282, 283]* as a first evaluation of our
model, implying first to define an expected behaviour and then how
to measure it.

<footnote>* https://github.com/fedden/RenderMan</footnote>
As our method targets both in generation and inference, we propose a twofold evaluation procedure. The first expected behaviour of the system is being able to reconstruct the samples from the original dataset. Hence, we evaluate the reconstruction performances performed by the models with different spectral losses, on a test set that has not been used for training.

The second aims to evaluate the emerging properties of the obtained latent representation. This is a rather hard task, this multi-dimensionality of this space preventing to obtain accurate visualizations. In addition to measure how much the obtained representation differs from the prior p(z), that is a training criterion, we also perform independence measures, disentanglement measures, and descriptor-based topology maps, to get additional insights on the organization of the latent space. Regarding baselines, the originality of inference and generation process of AEVB systems make them difficult to compare with other models. Therefore, we use invertible dimension reduction methods (PCA and ICA) to mirror the encoding / decoding processes, the latent space corresponding to the obtained projection.

**Models.** The flexibility of AEVB framework allows the full parametrization of both inference and generation processes, granting the right to choose almost every architecture parameters : encoding/decoding distributions, relative loss weights, and regularization loss. Moreover, we reviewed during last chapter the numerous evolutions proposed in the literature to tackle different aspects of the system, providing a plenteous amount of design choices.

As we can not obviously test all the models presented in the last chapter, we first chose to study the impact of latent space dimensionality (hence the compression performance of the model), the $D_{KL}$ weighting $\beta$ (hence, the impact of the representation rate), and the regularization strategy. To isolate the respective effects of these design choices we start from a standard VAE, and vary these factors one by one. We obtain the following experiments

* **latent dimensions**: four different VAEs of respectively 2, 4, 8 or 16 latent dimensions

* **regularization weighting**: three different $\beta$-VAEs, comparing $\beta \in [1, 4, 10]$

* **regularization divergence**: three different explicit divergences ($D_{KL}$, Renyi with $\alpha = 2$, Jensen-Shannon divergence) and two implicit divergences (MMD, Adversarial).

all of these experiments are applied to the three datasets toy_additive, toy_fm, and acidsInstruments-ordinario. The encoders and decoders are symmetric 4-layers convolutional auto-encoders, with respectively 3, 7, 9 and 11 kernels, 64, 32, 16 and 8 channels, ReLU non-linearities [284] and batch-norm normalization [285]. The models were all trained using the ADAM optimizer [123], with an initial learning rate of $1e-4$ and a batch size of 64. For the adversarial loss, the classifier is defined as a fully-connected two-layer discriminator with 500 linear units,
trained separately with the same optimization parameters.

### 3.2.2.1 Reconstruction evaluation.

We evaluate the performances our model in reconstruction on both *train* and *test* partitions of our model. The table containing the full results can be consulted in appendix A.1. As evaluating the reconstruction loss on the same criterion than the one used for training, we provide several spectral losses :

* *log-density*, the likelihood of the data under the generative model $p_{\theta}(x|\mu_p(z), \sigma_p^2(z))$ (training loss)

* $\ell^2$ loss, the mean-squared error between $\|\cdot\|^2$ the generated spectrum and the target

* *spectral convergence* (SC) : $\|x - \hat{x}\|^F / \|x\|^F$, where $\|\cdot\|^F$ is the Frobenius norm of the magnitude spectra $x$. This loss is mainly focused on the large frequency components [286]

* *Itakura-Saito divergence* (ISD) : a divergence specifically used for audio signals, related to perceptual divergence [287]

* *log-$\ell^1$*, the absolute $\|\cdot\|^1$ loss between the generated spectrum and the target

* *log-$\ell^2$*, the Euclidean $\|\cdot\|^2$ loss between the generated spectrum and the target

As our systems model the generative distribution $p(x|z)$ and that most of these criteria are taking definite values, we always choose the MAP estimation, that is the mean $\mu_p$ of $N(\mu_p(z); \sigma_p^2(z))$.

**Influence of dimensionality.** To study the influence of latent dimensionality on reconstruction results, we train models on each dataset with different latent shapes. A short results table is given 3.1 for the *acidsInstruments-ordinario* dataset, showing the log-density (training criterion), mean-squared error, and spectral convergence obtained with the test dataset.

The influence of latent dimensionality on reconstruction can be clearly observed, decreasing sharply from 2 to stagnate around 8 dimensions. This indicates that the model requires enough space to organize the whole dataset, emphasizing the role of VAE as a compression algorithm. However, even if the performances of two dimensions are clearly lower, reconstructions shown fig. 3.15 indicate that the VAE still honourably generates the given examples, showing that it is able to compress data with a very high rate (2 vs. 1024 dimensions). We also indicate the performances obtained with our baselines (PCA and ICA provide similar results in terms of reconstructions), that are clearly outperformed by the proposed model.

**Impact of regularization.** The divergence used to regularize the latent space has also an impact on the reconstruction, as shown by table 3.2. As suggested by the theoretical analysis of AEVB exposed in the

<table><caption>**Table 3.1:** Results obtained on test data with the *acidsInstruments-ordinario* dataset, for different latent space dimensionalities.</caption><thead><tr><th></th><th>-log p(x|z)</th><th>ℓ<sup>2</sup></th><th>SC</th></tr></thead><tbody><tr><td><strong>2d VAE</strong></td><td>334.34</td><td>68.12</td><td>21.5</td></tr><tr><td>Baseline</td><td>-</td><td>104.76</td><td>16.47</td></tr><tr><td><strong>4d VAE</strong></td><td>287.0</td><td>31.93</td><td>18.21</td></tr><tr><td>Baseline</td><td>-</td><td>92.7</td><td>17.02</td></tr><tr><td><strong>8d VAE</strong></td><td>269.9</td><td>20.11</td><td>12.28</td></tr><tr><td>Baseline</td><td>-</td><td>71.54</td><td>16.79</td></tr><tr><td><strong>16d VAE</strong></td><td>276.69</td><td>23.67</td><td>20.13</td></tr><tr><td>Baseline</td><td>-</td><td>48.65</td><td>15.16</td></tr></tbody></table>

<figure><img src="image_57.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>**Figure 3.15:** Examples or reconstructions obtained from a latent space with (a) 2 dimensions (b) 8 dimensions. While the system performs honourably with two dimensions, it cannot generate some examples.</figcaption></figure>
last section, using implicit regularizations strategies (MMD, Adversarial) indeed achieve better reconstruction results, outperforming methods using explicit divergences $D_{KL}$, Renyi, JSD. However, this gap is lower on test examples, pointing out that implicit methods may be more subject to over-fitting. Indeed, adversarial regularization performs very well with datasets of low variability (toy_additive_mini, acidsInstruments-ordinario) but performs badly with toy_fm, indicating that this regularization method prevents to efficiently compress the dataset. This may be explained by the inherent discriminating of this regularization process, isolating each individual example. Contrastly, MMD regularization (similar to WAE [84], see sec. 2.3.2.3) seems to avoid this over-fitting tendency, achieving good performances in any case.

Regarding explicit regularization methods, Rényi divergence with $\alpha = 2$ seems to perform slightly better than $D_{KL}$, suggesting that enforcing a little the zero-avoiding behaviour of the regularization may help the system to reconstruct the examples. However, Jenson-Shannon divergence seems to performs very badly, having the worst reconstruction scores for all the datasets. This may be due to the two-side construction of the JSD, that may be unsuitable as an optimization criterion.

**Impact of $D_{KL}$ weighting.** We also evaluated the effect of the $\beta$ term, that weighs the regularization term in the ELBO formulation (2.59), that can be used to control the rate of the representation (see sec. 2.16). In figure ., we can verify that increasing the rate of the representation, that is here the divergence $D_{KL}[q_{\phi}(z|x)\|p_{\theta}(z))]$, reduces the quality of the reconstruction. A short table is displayed fig. 3.3, showing the results obtained with the test set. We can thus verify that rate and distortion are complementary, the quality of reconstructions decreasing as $\beta$ is increased. However, as shown by the work of Higgins & al., increasing $\beta$ should increase the significance of the representation, enforcing it to represent factors of variation underlying the data. We will investigate that in the following section, studying the emerging properties of the latent space.

### 3.2.2.2 Representation evaluation

While evaluating the quality of reconstructions is a rather easy task, evaluating the emerging properties of the obtained latent spaces is more difficult. Indeed, the relatively high amount of dimensions prevents us to get a direct visualization of the latent space, and the projections applied to PCA or ICA dimensionality reduction methods may flatten the intrinsic latent organization. Therefore, we rely on different evaluators to describe the organization of latent projections, aiming to provide additional descriptors of its topology.

**Evaluation strategies.** As we said, visualizing the space is an hard task, as the number of dimensions is relatively high and that common dimensionality reduction techniques flatten the inner organization of

<table><caption>**Table 3.2:** Reconstruction scores obtained for different regularization divergences, using test partition of the datasets.</caption><thead><tr><th></th><th>-log p(x|z)</th><th>ℓ<sup>2</sup></th><th>SC</th></tr></thead><tbody><tr><td><strong>additive</strong></td><td></td><td></td><td></td></tr><tr><td>D<sub>KL</sub></td><td>488.37</td><td>13.21</td><td>11.07</td></tr><tr><td>MMD</td><td>474.75</td><td>6.06</td><td>6.67</td></tr><tr><td>D<sub>α=2</sub></td><td>487.84</td><td>13.01</td><td>9.03</td></tr><tr><td>JSD</td><td>519.1379</td><td>24.36</td><td>7.60</td></tr><tr><td>Adv</td><td>485.88</td><td>10.18</td><td>10.0</td></tr><tr><td>Baseline</td><td>-</td><td>137.29</td><td>160.42</td></tr><tr><td><strong>fm</strong></td><td></td><td></td><td></td></tr><tr><td>D<sub>KL</sub></td><td>910.10</td><td>340.40</td><td>15.81</td></tr><tr><td>MMD</td><td>766.8</td><td>221.78</td><td>62.19</td></tr><tr><td>D<sub>α=2</sub></td><td>210.86</td><td>3.40e16</td><td>218.86</td></tr><tr><td>JSD</td><td>259.07</td><td>242.90</td><td>27.60</td></tr><tr><td>Adv</td><td>789.56</td><td>525.16</td><td>6.46</td></tr><tr><td>Baseline</td><td>-</td><td>859.12</td><td>527.28</td></tr><tr><td><strong>aIns-o</strong></td><td></td><td></td><td></td></tr><tr><td>D<sub>KL</sub></td><td>276.69</td><td>23.67</td><td>20.13</td></tr><tr><td>MMD</td><td>248.19</td><td>8.43</td><td>5.31</td></tr><tr><td>D<sub>α=2</sub></td><td>255.74</td><td>12.38</td><td>10.53</td></tr><tr><td>JSD</td><td>259.07</td><td>14.22</td><td>7.60</td></tr><tr><td>Adv</td><td>246.7</td><td>7.71</td><td>10.0</td></tr><tr><td>Baseline</td><td>-</td><td>48.65</td><td>15.16</td></tr></tbody></table>

<table><caption>**Table 3.3:** Reconstruction scores obtained for different $\beta$, using test partition of the datasets.</caption><thead><tr><th></th><th>-log p(x|z)</th><th>ℓ<sup>2</sup></th><th>SC</th></tr></thead><tbody><tr><td><strong>toy_additive</strong></td><td></td><td></td><td></td></tr><tr><td>β = 1</td><td>276.69</td><td>23.67</td><td>20.13</td></tr><tr><td>β = 4</td><td>518.12</td><td>31.43</td><td>19.46</td></tr><tr><td>β = 10</td><td>553.48</td><td>57.91</td><td>37.85</td></tr><tr><td>PCA</td><td>-</td><td>137.29</td><td>160.42</td></tr><tr><td><strong>toy fm</strong></td><td></td><td></td><td></td></tr><tr><td>β = 1</td><td>825.75</td><td>262.16</td><td>15.817</td></tr><tr><td>β = 4</td><td>897.32</td><td>316.24</td><td>15.61</td></tr><tr><td>β = 10</td><td>1099.0</td><td>453.92</td><td>15.97</td></tr><tr><td>PCA</td><td>-</td><td>859.12</td><td>527.28</td></tr><tr><td><strong>aIns-o</strong></td><td></td><td></td><td></td></tr><tr><td>β = 1</td><td>261.37</td><td>15.34</td><td>13.74</td></tr><tr><td>β = 4</td><td>279.83688</td><td>27.77</td><td>16.56</td></tr><tr><td>β = 10</td><td>304.09</td><td>45.21</td><td>44.54</td></tr><tr><td>PCA</td><td>-</td><td>859.12</td><td>527.28</td></tr></tbody></table>
the latent space (see fig. 3.16). A first possible evaluation is to measure how much the encoding distribution fits the prior using four different divergences: $D_{KL}[q||p]$, $D_{KL}[p||q]$, $D_{\alpha=2}[q||p]$, $JSD[q||p]$, and MMD. Remembering sec. 2.1.3.2, note that the four first divergences fit individual projections $q(z|x)$, while the last one fit the aggregated posterior $q(z)$, and hence evaluates the global projection.

Another desired property of the representation is disentanglement, quantifying how much the axis of the obtained representation is reflecting the generative factors of variation of the dataset. However, we do not always have the required information to track the corresponding organization in the latent space. A way of measuring the disentanglement of a representation is to compute the total covariance (TC) of the representation, that is the $D_{KL}$ between the full posterior distribution and the product of its components [150]

$$TC(z) = D_{KL}[q(z) \parallel \prod_{i} q_i(z_i)] \quad (3.7)$$

However, the computation of the aggregate posterior $q(z)$ is impossible, as its expression is intractable [69, 288]). A way to estimate it is to marginalize on $x$, $q(x)$ being the empirical distribution of the data:

$$q(\mathbf{z}) = \int q(\mathbf{z}|\mathbf{x})q(\mathbf{x})d\mathbf{x} = \int q(\mathbf{z}|\mathbf{x})\frac{1}{N}\sum_{n=1}^{N}\delta[\mathbf{x} = x_n]d\mathbf{x} = \frac{1}{N}\sum_{n=1}^{N}q(\mathbf{z}|\mathbf{x}_n) \quad (3.8)$$

However, the $D_{KL}$ of total correlation is still intractable, as $q(z)$ is a mixture of normal distributions. Another way to estimate $TC[q(z)]$ is to use DRE (see 2.13), that is the density ratio between the estimated aggregated posterior and the product of its components $\hat{q}(z) = \prod_{n=1}^{N} q_i(z_i)$ [289]. In that case, the contrasting distribution $\hat{q}(z)$ is obtained by randomly shuffling batches across the latent dimensions. Total correlation can thus be estimated with

$$T\hat{C}[q(\mathbf{z})] = \mathbb{E}_{q(\mathbf{z})} \frac{D(\mathbf{z}^{\text{true}})}{1 - D(\mathbf{z}^{\text{false}})} \quad (3.9)$$

where the discriminator $D$ is also learned during learning. An additional whitening procedure is also proposed in [290] using PCA decomposition, such that the total covariance is rather processed in the PCA space.

While Total Covariance estimate the independence between representation axis, they do not explicitly link them to existing factors of variation of the data. In their β-VAE paper, Higgins et al. propose a validation procedure using external label information [194]. Based on the intuition that variance of latent projections should be lower when a factor is fixed, this procedure trains a linear classifier to predict the targeted labels when given the pairwise distances of their corresponding projections. While the original paper uses a data generator to perform this validation procedure, this evaluation can still be performed with the corresponding metadata. This procedure is based on the following steps:

<figure><img src="image_59.png" /><figcaption>Figure 3.16: PCA representation of a VAE trained in acidsInstruments-ordinario, each point being a data point projection coloured by its pitch class. We can see that the extracted representation has been successfully shaped as an isotropic normal distribution, but that this visualization does not provide any information on the classwise repartition.</figcaption></figure>
<figure><img src="image_60.png" /><figcaption>**Figure 3.17:** Statistics of latent projections obtained with two VAEs, one trained with $D_{KL}$ and one trained with adversarial criterion. (a) variances of latent means, obtained with $D_{KL}$ regularization (b) means of variances obtained, with $D_{KL}$ regularization (c) variance of means, obtained with adversarial regularization (d) means of variances, obtained with adversarial regularization</figcaption></figure>

<table><caption>**Table 3.4:** Latent organization descriptors obtained from the three datasets.</caption><thead><tr><th></th><th>JSD<sub>0.5</sub>[p||q]</th><th>MMD[p, q]</th><th>TC[p, q]</th></tr></thead><tbody><tr><td><strong>additive</strong></td><td></td><td></td><td></td></tr><tr><td>D<sub>KL</sub></td><td>104.05</td><td>159.52</td><td>22.89</td></tr><tr><td>MMD</td><td>521.5e9</td><td>588.83</td><td>22.67</td></tr><tr><td>Adv</td><td>8.81e+31</td><td>117231</td><td>211.75</td></tr><tr><td>Baseline</td><td>-</td><td>759.39</td><td>400.13</td></tr><tr><td><strong>fm</strong></td><td></td><td></td><td></td></tr><tr><td>D<sub>KL</sub></td><td>2176</td><td>2406</td><td>65.32</td></tr><tr><td>MMD</td><td>9.67e27</td><td>405.44</td><td>35.98</td></tr><tr><td>Adv</td><td>1.13e+29</td><td>143165</td><td>409.4</td></tr><tr><td>Baseline</td><td>-</td><td>207162</td><td>406.33</td></tr><tr><td><strong>aIns-0</strong></td><td></td><td></td><td></td></tr><tr><td>D<sub>KL</sub></td><td>59.63</td><td>35.72</td><td>23.36</td></tr><tr><td>MMD</td><td>65.98e12</td><td>588.83</td><td>19.35</td></tr><tr><td>Adv</td><td>4.37e30</td><td>471229</td><td>535</td></tr><tr><td>Baseline</td><td>-</td><td>39409</td><td>204.82</td></tr></tbody></table>

* create K classifiers for each of the K independent factors of variations

* randomly choose a factor of variation $k \in [0, K]$

* obtain two sets $\{x_1\}$ and $\{x_2\}$, with a factor K fixed to given value k

* obtain corresponding latent vectors $z^1$ and $z^2$

* compute the average L1 norms $z_{diff} = |z^1 - z^2|$

* train classifier K on $z_{diff}$ to predict the corresponding factor k

Concomitantly, a related procedure was proposed by Eastwood & Williams, where all the input $z$ is directly used to feed each factor-wise classifier. This way, we can directly obtain the responsibilities of each axis, and identify the corresponding permutation of the axis to get the desired representation [291]. A disentanglement score can then be computed using the entropy of each axis responsibility $D_i = (1 + \sum_{k=0}^{K} P_{ik} \log_K P_{ik})$, where $P_{ij} = R_{ij} / \sum_{k=0}^{K} R_{ik}$ is the responsibility of $z_i$ for the corresponding factor. This approach allows to evaluate the completeness of the representation, evaluating simultaneously the different tasks.

**Impact of divergence strategies.** The full results table is displayed in appendix A.2. A summarized table is provided 3.4, where we show three different latent descriptors

* the Jensen-Shannon divergence, that is a symmetrized version of $D_{KL}$ that measures the average distance between $q_\phi(z|x)$ and $(z)$,

* the MMD, that computes the divergence between the aggregated posterior $q(z)$ and $p(z)$

* the total covariance $TC[q(z)]$, obtained with the method presented in (3.9)

We can see that the divergence obtained with JSD on models trained with implicit distributions is very high. This can be observed by comparing the statistics of the obtained latent projections with two different regularizations, one with $D_{KL}$ and the other with adversarial (see fig 3.17). In the case of $D_{KL}$, the variances of the variational distributions' means are very close to 1, while their variance varies between 0 and 1, showing that the space is well regularized. Conversely, the latent space obtained with an adversarial criterion is *atomic*, such that no information is shared between examples. However, the total covariance shows that the MMD regularization efficiently provides orthogonal axis, showing that it successfully matched the aggregated posterior with the prior.

**Impact of $\beta$ strategies.** Table 3.18 show the divergence obtained when varying the $\beta$ factor with three different estimators: the training criterion $D_{KL}[q||p]$, plus two divergences $MMD[p, q]$ and $TC[p, q]$, that describe how the aggregate posterior fits the prior $N(0, I)$. We can see that, as expected, increasing the $\beta$ factor will reduce the $D_{KL}$; however, more surprisingly, we can see that increasing $\beta$ does not decrease the TC, showing that it does not necessarily enforce the independence of
latent axis (the TC is even increasing with *toy-fm*). A further analysis of the latent statistics indeed points out that increasing *β* does not seem to *orthogonalize* the latent representation but rather to *polarize* it, at least for magnitude spectra. Indeed, as can be seen fig. 3.18, all the dimensions of models trained with *β* = 1 seem to bring an equal contribution of the latent code, while with *β* = 10 seems to clearly deactivate some dimensions, where the mean variances of the corresponding latent dimension *z*<sub>j</sub> equals 1 (see section 2.3.2.3). Therefore, increasing the regularization weight does not spread the latent representation, but rather concentrates the information on a fewer number of axis.

**Disentangling.** While the estimators used so far allow to measure the performance of the regularization and the latent axis independence, it does not quantify if the obtained dimensions reflect the generative factors of each dataset. We show table 3.6 the disentanglement results obtained with the dataset *acidsInstruments-ordinario*, using the Higgins' disentangling measure. To ensure the consistency of these measurements, we also perform prediction from the trained classifier on random labels, and compare the scores with the ones obtained with true labels.

While the best results are obtained, for octave and pitch, with the model with 8 latent dimensions, we can see that the difference with the results taken with random labels is little, showing that the classifier could have learn the distribution of the point-wise distances. We also see that increasing the $\beta$ do not especially lead to better disentanglement measures, even worsening it in the case of the *instrument* class, that is surprisingly best disentangled by ICA. We obtain similar behaviours with other dataset, showing that the latent spaces do not automatically disentangle the generative factors of audio datasets, at least with the Higgins evaluation method.

**Perceptive maps.** Finally, another method to analyse the properties of the latent space is rather to explore its topology in the data domain, performing grid sampling of the latent space and computing acoustical descriptors of the corresponding generation. However, grid-sampling scales exponentially with the number of dimensions, and can then be expensive to compute. Hence, we resort to invertible dimensionality reduction algorithms (PCA and ICA) to project the latent space in three dimensions, and then obtain the *perceptual maps* for the corresponding model. These projections are computed on the latent variables

<table><thead><tr><th></th><th>octave</th><th>random</th><th>pitch</th><th>random</th><th>instrument</th><th>random</th></tr></thead><tbody><tr><td>2d</td><td>2.52</td><td>2.46</td><td>2.72</td><td>2.74</td><td>2.74</td><td>2.72</td></tr><tr><td>4d</td><td>2.55</td><td>2.42</td><td><strong>2.41</strong></td><td>2.67</td><td>2.79</td><td>2.76</td></tr><tr><td>8d</td><td>2.51</td><td>2.47</td><td>2.60</td><td>2.65</td><td>2.75</td><td>2.65</td></tr><tr><td>D<sub>KL</sub></td><td>2.77</td><td>2.47</td><td>2.76</td><td>2.66</td><td>2.46</td><td>2.69</td></tr><tr><td>MMD</td><td>2.4</td><td>2.35</td><td>2.50</td><td>2.64</td><td>1.68</td><td>1.75</td></tr><tr><td>Adv</td><td>6.94</td><td>7.52</td><td>10.112</td><td>10.03</td><td>8.74</td><td>8.25</td></tr><tr><td>β = 4</td><td>2.64</td><td>2.45</td><td>2.75</td><td>2.62</td><td>2.38</td><td>2.61</td></tr><tr><td>β = 10</td><td><strong>2.33</strong></td><td>2.33</td><td>2.65</td><td>2.58</td><td>2.78</td><td>2.57</td></tr><tr><td>PCA</td><td>3.00</td><td>3.04</td><td>3.35</td><td>3.27</td><td>3.19</td><td></td></tr><tr><td>ICA</td><td>2.41</td><td>2.265</td><td>2.53</td><td>2.53</td><td><strong>1.61</strong></td><td>2.53</td></tr></tbody></table>

<figure><img src="image_61.png" /><figcaption>**Figure 3.18:** Statistics obtained from a VAE trained with $\beta$ = 10. (a) variances of means obtained with $D_{KL}$ (b) means of variances obtained with $D_{KL}$ (c) variances of means obtained with adversarial (d) means of variances with adversarial</figcaption></figure>

<table><caption>**Table 3.5:** Table of results for latent descriptors obtained with models with various $\beta$. As we can see, increasing $\beta$ does not lead to reduce the total covariance TC</caption><thead><tr><th></th><th>D<sub>KL</sub>[q||p]</th><th>MMD[p, q]</th><th>TC[p, q]</th></tr></thead><tbody><tr><td>additive</td><td></td><td></td><td></td></tr><tr><td>β = 1</td><td>35.20</td><td>159.52</td><td>22.89</td></tr><tr><td>β = 4</td><td>18.28</td><td>143.10</td><td>19.83</td></tr><tr><td>β = 10</td><td>12.47</td><td>140.6</td><td>20.10</td></tr><tr><td>PCA</td><td>-</td><td>759.39</td><td>202.13</td></tr><tr><td>ICA</td><td>-</td><td>8875.41</td><td>6.77e-4</td></tr><tr><td>toy-fm</td><td></td><td></td><td></td></tr><tr><td>β = 1</td><td>49.20</td><td>2406</td><td>65.33</td></tr><tr><td>β = 4</td><td>34.77</td><td>2534</td><td>89.76</td></tr><tr><td>β = 10</td><td>24.87</td><td>2419.8</td><td>113.10</td></tr><tr><td>PCA</td><td>-</td><td>205598</td><td>405.70</td></tr><tr><td>ICA</td><td>-</td><td>7594.93</td><td>6.10e-4</td></tr><tr><td>aIns-o</td><td></td><td></td><td></td></tr><tr><td>D<sub>KL</sub></td><td>19.28</td><td>2840.25</td><td>50.27</td></tr><tr><td>β = 4</td><td>9.36</td><td>2527.10</td><td>50.90</td></tr><tr><td>β = 10</td><td>5.423</td><td>1807.7</td><td>45.52</td></tr><tr><td>PCA</td><td>-</td><td>39409</td><td>204.82</td></tr><tr><td>ICA</td><td>-</td><td>7881.27</td><td>3.92e-4</td></tr></tbody></table>

<caption>**Table 3.6:** Disentanglement scores for tasks *octave*, *pitch*, *instrument* obtained with the dataset *acidsInstruments-ordinario*.</caption>
obtained from examples of the dataset, and provide an interesting analysis of the *local* behaviour of the latent space. These perceptive maps provide intuitive feedback of the latent organization, and could be advantageously used for interactive exploration.

First, we use the *spectral centroid*, that is the weighted-mean of frequencies according to their corresponding amplitude [292]

$$ \text{Centroid}[X(f)] = \frac{\sum_{i=0}^{N-1} f(i)X_i(f)}{\sum_{i=0}^{N-1} f(i)} \quad (3.10) $$

that has a connection with the perception of *brightness* of the corresponding sound. The second descriptor is *bandwidth*, that is defined [293]

$$ \text{Bandwidth}[X(f)] = \left[ \sum_{i=0}^{N-1} X(f)(f(i) - \text{Centroid}[X(f)])^p \right]^{1/p} \quad (3.11) $$

that estimates the width of a the spectra over its frequency range, and is related to the perception of *narrowness* of a sound. Finally, we also compute the *flatness* descriptor [294]

$$ \text{Flatness}[X(f)] = \frac{\exp\left(\frac{1}{2\pi}\int_{-\pi}^{\pi}\log X(f)df\right)}{\frac{1}{2\pi}\int_{-\pi}^{\pi}X(f)df} \quad (3.12) $$

that can be roughly described as the difference between the marginal entropy $E[X(f)]$ of the spectrum distribution $X(f)$ and the Kolmogorov-Sinai entropy, that estimates the entropy of the conditional distribution $E[X(f)_k|X(f)_{k-1}]$. Flatness is then related to the *noisiness* of the spectrum, where signals close to noise will verify $\text{Flatness}[X(f)] \to 0$, while signal with a more *deterministic* structure such that pure sounds will rather have high values.

An example of such maps are given fig. 3.19 for *toy_additive* and *acidsInstruments-ordinario* using PCA and ICA. These maps allow us to draw very interesting conclusions on the *perceptual* organization of this space. First, we can see that perceptual maps obtained with PCA, that lies close to the projected points, provides perceptually smooth maps, showing the regularization effect of the VAE. However, this contrasts with the maps obtained with ICA, showing rather complex latent structures. While the maps obtained with ICA are sufficiently smooth when trained with *toy_additive* with $D_{KL}$ regularization, the ones obtained on *acidsInstruments-ordinario*, with adversarial divergence are more irregular. We can explain this by the fact that ICA, contrary to PCA, enforces the orthogonality of its axis and thus can represent untrained zone of the latent space.

**Overall conclusions.** These experiments then allow to draw important conclusions on how vanilla auto-encoders behave with magnitude spectrum. First, we can see that VAEs can compress the dataset very efficiently, even providing convincing reconstruction results with only
<figure><img src="image_62.png" /><figcaption>**Figure 3.19:** Descriptor plots obtained from (left) toy_additive with $D_{KL}$, with PCA (above) and ICA (below), and from (right) acidsInstruments-ordinario with adversarial, with PCA.</figcaption></figure>

two latent dimensions. Regarding the impact of regularization methods, MMD and adversarial divergences allow to reduce even more the reconstruction error, at the cost of developing atomic representations that discourage information sharing between examples. We also showed that the $\beta$ parameter, at least for magnitude spectra, does not encourage the orthogonality of latent axis but rather the *sparsity* of the latent dimensions, encouraging fewer active latent units at the cost of lower reconstruction abilities. However, regardless of the chosen divergence, the system manages to regularize the representation on the chosen prior, showing the efficiency of the regularization process. However, we saw that VAE *did not* disentangle inner factors of variation for audio data, at least for this standard configuration*. Further investigations of disentanglement abilities for audio signals could be pursued by taking even simpler datasets and black-box analysis approaches, such as the Riemannian approaches such as the ones proposed in [197]; however, we leave that to future work. In this PhD, this disentanglement issues rather motivated us to develop regularization strategies specific to audio signals, relying on the flexibility of variational process. Such strategies would allow us to specify externally latent organization policies for the system, that would be beneficial from a design point of view. We thus propose two different regulariza-

<footnote>* we did similar experiments with fully-connected instead of convolutional networks, and obtained similar conclusions</footnote>
tion strategies: one based on *symbol-signal translation spaces*, allowing the user to define custom vocabularies that would be reflected in the latent space, and another one based on metrics matching, based on the perceptual regularization of the latent space.

## 3.3 Bijective signal-symbol regularization

In the previous section, we studied the emerging properties of the latent representations extracted with AEVB trained on magnitude spectra. We concluded that, while the obtained representations were able to drastically reduce the dimensionality of data, the interpretability of latent features was still difficult and not representing inner factors of variation. However, in some applications, we can access complete or partial metadata on the dataset, such that we may want to enforce the representation to reflect these properties. Furthermore, some methods take inspiration of *domain transfer* to propose multi-modal latent spaces, such that we can infer the symbolic information from the signal and reversely, using the latent space as a *translation space*. Such models can be very beneficial for the proposed approach, as this could allow us to explicitly model generation parameters without interfering too much with the unsupervised learning process of AEVB. In this section, we will summarize the most common methods to influence the latent representation with external information, from full supervised approaches such as *conditioning* to loose semi-supervised methods such as *domain translation* (see fig. 3.20). We will then propose an application of such methods to musical signals, using the latent space as signal-symbol translation spaces between some orchestral instruments and symbolic information (pitch, octave, and dynamics).

### 3.3.1 (Semi-)supervised methods for latent regularization

Several methods can be used to introduce external symbolic information in the learning process of AEVB models. Some approaches are based on *conditioning*, where the symbols are directly provided to the encoding/decoding modules; however, these methods are entirely supervised, such that the model is untrainable if we miss some label information. Other methods propose to alleviate this constraint with *semi-supervised* learning, where an symbolic inference module is added to AEVB such that the system is able to recover the missing information. Finally, some methods take inspiration from *domain transfer* to provide shared representation between the data and corresponding symbols.

**Conditioning.** We can influence the latent representation by conditioning whether the variational or generative distributions on the label information. Given a dataset {$x, y$}_{i=0...N} composed by data x and corresponding symbols y, provided as class labels, we can *condition* the encoder by providing the label information as an additional input,
modeling a variational distribution $q(z|x, y)$. Hence, the symbol information can directly influence the latent space through the encoder. Similarly, we can also condition the decoder to model the generative distribution $p(x|z, y)$, where we also add the label information to the module's input. As the variables $z$ and $y$ are assumed independent, conditioning the generation model on symbols enforces the *invariance* of the latent representation to the corresponding label, then disentangling the space from the class information. The obtained latent space is thus shared between classes, and we can then *conditionally generate* from the decoder to obtain samples belonging to the target label. These models can then be used for between-class translation, encoding a sample and decoding the corresponding latent vector with another class label. This translation have been showed to also transfer the stylistic properties of the encoded data to the decoded sample, thus proving the consistency of the obtained representation [295]. Furthermore, Yan & al. showed that model was able to perform *class interpolation*, interpolating between the two one-hot vectors $y_1$ and $y_2$ [296].

This method, that we call *concatenative conditioning*, allows us to efficiently obtain class-invariant latent representations. However, concatenating the symbolic information to the data may be not sufficient to accurately condition the decoder to the class information. Conditional Normalization methods rather propose to condition the normalization layers of neural networks, then to directly condition some of the network parameters to the class information. The *Conditional Instance Normalization* proposed by Dumoulin, Shlens & al., have been successfully applied to style transfer in the image processing domain, providing interesting translation results [297]. This method have then been extended by Perez, Strub & al. with Feature-wise Linear Modulation (FiLM), that uses embeddings extracted from class information to directly modulate the features provided by each layer of the network [298]. While these conditioning methods are widely used in image processing, their use in the audio domain is still marginal. However, FiLM modulators have been successfully applied to many-to-many instrument translation, showing the efficiency of these methods for audio [299].

**Semi-supervised learning.** Despite the efficiency of advanced conditioning methods described above, the availability of symbolic information for the whole dataset can sometimes be burdensome. We can circumvent this requirement by modeling the label information $y$ as a random variable, and adding an inference model $q(y|x, z)$ to the system that is able to predict the class information in the case of incomplete data. The symbolic variables $y$ can also be modelled as additional latent variables, that are extracted with a second recognition model $q(y|x)$. The latent representation is then disentangled between class-related and class-insensitive features [300].

Alternatively, we can infer the symbol information using the latent space, running a discrimination task on a high-level representation of the data. This framework, called *semi-supervised learning*, combines the benefits of both supervised and unsupervised approaches, and allow us to train the model on data even if the corresponding label

<figure><img src="image_63.png" /><figcaption>**Figure 3.20:** Four different approaches for conditioned data generation. (a) un-conditioned system (b) *direct conditioning*: the label information is directly given to the decoder, making the latent representation independent to the corresponding task. (c) *semi-supervised learning*: an additional classification system is added on the top of latent space, predicting the label it is missing (d) *translation*: the latent space is used as a *translation space*, that is shared between labels and data</figcaption></figure>
information is missing. This framework was first proposed by Kingma & al. [295], proposing a 2-layer variational model where label information was extracted from the first latent layer, and the given to the decoder for concatenative conditioning. Provided the low capacity of the discriminator (such as a simple linear classifier), the variational distribution will then be enforced to linearly separate the latent projections of the different classes. Semi-supervised learning is then also an efficient way to regularize the latent space of the variational model on external symbolic information, influencing the representation to discriminate examples belonging to different classes. However, the original model from Kingma & al. does not explicitly model the distribution $p(\mathbf{z}|\mathbf{y})$, limiting the interpretability of the process. We then take inspiration from domain-translation, a more general framework that can be applied for multi-modal inference.

**Multi-modal domain-translation.** The different approaches described in the last paragraphs can be considered as a particular case of *multi-modal learning*. Multi-modal learning addresses the problem of modeling a dataset composed by $K$ different views of a same example $\mathbf{X} = \{\mathbf{x}^1, \mathbf{x}^2...\mathbf{x}^K\}$, that can possibly be incomplete when some views are missing. A way to extract useful features from multi-modal datasets is to infer a representation $q(\mathbf{z}|\mathbf{X}) = q(\mathbf{z}|\mathbf{x}_1) = ... = q(\mathbf{z}|\mathbf{x}^K)$ shared by the different views, such that all views $\mathbf{x}^1, \mathbf{x}^2...\mathbf{x}^K$ are represented by the same latent vector $\mathbf{z}$. Provided the corresponding generation models $p(\mathbf{x}_k|\mathbf{z})$, this shared representation can be then used as a translation space by modeling the distribution $q(\mathbf{x}_k|\mathbf{x}_m) = p(\mathbf{x}_k|\mathbf{z})q(\mathbf{z}|\mathbf{x}_m)$, using the same vector $\mathbf{z}$. Sharing latent representations then allow to perform both domain translation and direct generation, by directly sampling all the views from a latent vector $\mathbf{z}$.

Hence, a variational formulation of this setting would require a set of encoders $q(\mathbf{z}|\mathbf{x}_1), q(\mathbf{z}|\mathbf{x}_2), ..., q(\mathbf{z}|\mathbf{x}_K)$ and a set of decoders $p(\mathbf{x}_K|\mathbf{z}_1), p(\mathbf{x}_K|\mathbf{z}_2), ..., p(\mathbf{x}_K|\mathbf{z}_K)$, and formulating the overall loss by accumulating the individual ELBOs of each domain

$$ \mathcal{L}_{\text{total}} = \mathcal{L}_1(q^1, p^1) + \mathcal{L}_2(q^2, p^2) + \dots + \mathcal{L}_K(q^K, p^K) $$

however, this criterion does not provide any explicit transfer objective between different views. Taking two domains $\mathbf{x}$ and $\mathbf{y}$, Vedentam proposes to model both single-domain encoders $q(\mathbf{z}|\mathbf{x})$ and $q(\mathbf{z}|\mathbf{y})$ and joint encoder $q(\mathbf{z}|\mathbf{x}, \mathbf{y})$, and learn on the triple ELBO objective $\mathcal{L}[q(\mathbf{z}|\mathbf{x}, \mathbf{y})] + \mathcal{L}[q(\mathbf{z}|\mathbf{x})] + \mathcal{L}[q(\mathbf{z}|\mathbf{x}, \mathbf{y})]$ with shared decoders $p(\mathbf{x}|\mathbf{z})$ and $p(\mathbf{z}|\mathbf{x})$ [301]. This objective ensures the symmetry of the inference process between $\mathbf{x}$ and $\mathbf{z}$, and enforces the single-domain encoders to match the joint encoder. Alternatively, the SCAN method proposed by Higgins & al. bypass the addition of a supplementary by adding an explicit transfer penalty to the ELBO, giving [302]

$$ \mathcal{L}_{\text{SCAN}} = \mathcal{L}_1(q^1, p^1) + \mathcal{L}_2(q^2, p^2) + D_{KL}[q^1 \| q^2] \quad (3.13) $$

where the sense of the additional $D_{KL}$ between the two posterior
approximations has been chosen to encourage the mass-covering behaviour of $q(\mathbf{z}|\mathbf{y})$ with respect to $q(\mathbf{z}|\mathbf{x})$, such that the symbolic posterior $q(\mathbf{z}|\mathbf{y})$ does not under-match any example $\mathbf{x}$ of the class. Adding this transfer term then allows us to directly match the latent representations of both models without any additional module, thus speeding up the process. Suzuki & al. propose to match both approaches, modeling the joint distribution and adding an explicit transfer term [303]. Liu & al. also propose to enforce even more this transfer by sharing the weights of upper layers of each encoder-decoder pairs, and use an additional adversarial criterion to evaluate the transferred data [304].

### 3.3.2 Signal-symbol translation with latent space matching

**Motivations.** Our motivation for investigating latent regularization strategies based on external symbolic information is twofold. From an analysis perspective, recovering symbols from latent spaces obtain from audio signals encompass most of MIR techniques such as instrument recognition, pitch extraction, and many others. However, most machine-learning based MIR methods are entirely *supervised*, such that the system is untrainable if the data-symbol pairs are incomplete. Furthermore, models trained from an supervised manner often provide representations that are too task-specific, such that extracted features can represent poorly the underlying structure of data and cannot be used for other tasks. Semi-supervised methods then alleviate this drawbacks, being able to both extract robust representations from data and efficiently extract the corresponding symbol information. These methods can then allow us to perform *multi-task learning*, where the latent space provide efficient representations that is enforced to represent diverse modalities of the data.

Moreover, using latent representations as a bijective *translation space* between symbols audio signals is interesting from a creative point of view. Indeed, considering symbolic information as alternative modalities of the data then allows to learn on arbitrary data/symbol pairs, and can be customized by the user during the training process. The method can thus be used for *constrained generation*, where the signal generator is conditioned on complete or partial symbol information, but also reversely for inferring the corresponding symbols from a given audio signal. Furthermore, this framework still allows free navigation in the latent space, such that both data and symbols are generated from $\mathbf{z}$, then providing useful symbolic feedback for generation.

**Formulation.** We then extend the audio data $\mathbf{x}$ with symbol labels $\mathbf{y}_1, \dots, \mathbf{y}_K$ to obtain our multi-modal dataset $\{\mathbf{x}, \mathbf{y}_1, \dots, \mathbf{y}_K\}_{i=1}^N$.

<figure><img src="image_64.png" /><figcaption>Figure 3.21: Two VAEs, one in the signal domain $\mathbf{x}$, one in the signal domain $\mathbf{y}$, that share the same representation.</figcaption></figure>

$$q(\mathbf{z}|\mathbf{y}) = \mathcal{N}(\mu_q(\mathbf{y}), \sigma_q^2(\mathbf{y}))$$

$$p(\mathbf{y}|\mathbf{z}) = \prod_{i=1}^{L} \mathcal{B}(\mu_{p,i}(\mathbf{z})) \quad \text{or} \quad \prod_{i=1}^{L} \text{Cat}(\mu_{p,i}(\mathbf{z}))$$
where we model the label density with a Bernoulli distribution $\mathcal{B}(\mu_{p,i}(\mathbf{z}))$ if the label is a binary vector, or a Categorical distribution $\text{Cat}(\mu_{p,i}(\mathbf{z}))$ if the label correspond to a given class. Note that, in the Bernoulli case, each dimension of a label information $[\mu_{k,1}, \mu_{k,2}, \dots, \mu_{k,i}]$ is independent while in the Categorical case class probabilities $[\mu_{k,1}, \mu_{k,2}, \dots, \mu_{k,i}]$ sum to one. Taking inspiration from the SCAN model proposed by Higgins & al. [302], we train our auto-encoders on the multi-modal data dataset $\{\mathbf{x}, \mathbf{y}_1, \dots, \mathbf{y}_K\}_{i=1}^N$ using a VAE for the signal domain $\mathbf{x}$, and another VAE on the joint label $\{\mathbf{y}_1, \dots, \mathbf{y}_K\}$ space, on that we train on individual ELBO plus an additional transfer term $D_{KL}[q(\mathbf{z}|\mathbf{x})||q(\mathbf{z}|\mathbf{y})]$ (see 3.13).

We propose to evaluate the performance of this framework on bijective audio-signal translation by addressing an *audio transcription* challenge, where label information is a triplet [*pitch*, *octave*, *dynamics*]. The corresponding generative distribution is then

$$p(\mathbf{y}|\mathbf{z}) = \prod_{i=1}^{K} p(\mathbf{y}_i^p|\mathbf{z})p(\mathbf{y}_i^o|\mathbf{z})p(\mathbf{y}_i^d|\mathbf{z}) \quad (3.14)$$

where the index $i$ denotes the number of target instruments. In this case, we use a single symbolic VAE to output the concatenated symbols $[\mathbf{y}_1^p, \mathbf{y}_1^o, \mathbf{y}_1^d, \dots, \mathbf{y}_K^p, \mathbf{y}_K^o, \mathbf{y}_K^d]$ to lighten the process, and such that sharing the encoding/decoding functions allows the VAE to internally model correlations between the different symbols. The total loss of the joint model is then, reformulating 3.13

$$ \mathcal{L}_{\text{SCAN}} = \mathbb{E}_{z \sim q(\mathbf{z}|\mathbf{x})} \left[ \log p(\mathbf{x}|\mathbf{z}) - \beta D_{KL}[q(\mathbf{z}|\mathbf{x})||p(\mathbf{z})] \right] \quad \text{data-domain ELBO} \quad (3.15) $$

$$ + \mathbb{E}_{z \sim q(\mathbf{z}|\mathbf{y})} \left[ \log p(\mathbf{y}|\mathbf{z}) - \beta D_{KL}[q(\mathbf{z}|\mathbf{y})||p(\mathbf{z})] \right] \quad \text{symbol-domain ELBO} \quad (3.16) $$

$$ -\alpha D_{KL}q(\mathbf{z}|\mathbf{x})||q(\mathbf{z}|\mathbf{y})] \quad \text{weighted transfer penalty} \quad (3.17) $$

**Experiment details.** We evaluate the proposed approach on two datasets. We first perform training on the proposed *acidsInstruments-ordinario* dataset, composed by individual samples from classical instruments, where the corresponding symbolic information is entirely available. In the multi-instrument detection case, we obtain the mixtures of audio signals by randomly sampling data samples over different instruments and sum the obtained spectra, provided the linearity of the spectral transform. In this study, we selected five different instruments (violin, alto-sax, flute, piano, and trumpet-C) and evaluate the model on each individual domain, and on instruments mixtures up to three different sources (saxophone + violin, and saxophone + violin + trumpet-C). We obtained the spectral representation from audio signals using an NSGT transform, based on a Constant-Q scale with 48 bins per octave. We only take the stationary part of each sound for training, that we obtain by cropping each spectra from the 20th to the 50th frame. We also evaluate our models on a validation dataset,
ensuring the robustness of obtained representation. However, datasets
of single instruments recordings with corresponding annotations are
uncommon, such that we could not provide validation measures for
every instruments. Fortunately, Cantos & al. proposed a dataset of
annotated flute recordings * gathering single notes, arpeggios, and im-
provision fragments, that we can use to validate our models trained
on flute sounds [305].

Regarding the models, the models used for both data and symbolic
domains are simple VAEs with 32 dimensions, with MLPs modules
as encoding/decoding functions. For the data domain, the MLPs are
given 2 layers of 2000 units when learning single instruments, and
5000 units when learning mixtures of two instruments or above, using
in both cases ReLU non-linearity and batch-normalization. For the
symbol domain, networks' architectures was set to 2 layers of 1000
units, regardless of the instruments number. Both systems were jointly
optimized using ADAM gradient descent with an initial learning rate
of 1e-3, progressively shrunk during training depending on the loss
derivative. The β weighting ot the regularization term is set to 1, with
a warm-up procedure during the first 100 epochs. The α weight of the
transfer loss is set to 10, as recommanded in the original SCAN article.
We also compared the obtained results in symbolical inference with a
baseline model that mimicking the VAE procedure from audio to data,
first applying a PCA with 32 dimensions and then using a two-layer
MLP for classification.

3.3.3 Evaluation of domain transfer and further
applications.

Our evaluation strategy is split in three parts : we first obtain the
performances achieved in data domain, both on reconstruction and
transfer, and comparing the two to catch the consistency of translation
space. We then do the same analysis in the symbolic domain, focusing
more on transfer as the reconstruction quality of the symbolic part
is quite secondary in this application. Then, we evaluate the perfor-
mances of the proposed model on a validation dataset made of flute
samples, and provide examples of the diverse creative applications us-
ing the proposed system. This article has been submitted and accepted
at DaFX 2019 conference, and will be released soon through the arXiv
plateform. The corresponding code and support page can found on
github.

**Audio reconstruction results.** Reconstruction results are shown ta-
ble 3.7, where we compare the log-likelihood and Itakura-Saito Diver-
gence (ISD) obtained of reconstructions – log p(x|z) and transferred
data from symbols – log p(x|y) (see 3.2). Both scores are presented for
signal-to-signal reconstruction (left) and symbol-to-signal inference

<footnote>* The dataset is available at: https://zenodo.org/record/1408985</footnote>
<figure><img src="image_65.png" /><figcaption>**Figure 3.22:** Reconstruction and transfer performed by the system. We can see that both transfer and reconstructions are perfect with one instrument. While the performances are still satisfactory with two instruments, they decrease with three instruments, showing that the system struggles to disentangle correctly the sources. This can be observed by observing the reconstructions and transfers in the symbolic domain: while the most symbols are accurately found, they are affected to the wrong instrument, showing the limitations of the followed approach.</figcaption></figure>

<table><caption>**Table 3.7:** Signal reconstruction and transfer performances</caption><thead><tr><th></th><th>-log p(x|z)</th><th>ISD</th><th>-log p(x|y)</th><th>ISD</th></tr></thead><tbody><tr><td>Alto-Sax (Sax)</td><td>-694.1</td><td>0.093</td><td>-416.6</td><td>0.177</td></tr><tr><td>Violin (Vn)</td><td>-671.4</td><td>0.104</td><td>-551.1</td><td>0.151</td></tr><tr><td>Trumpet-C (TpC)</td><td>-706.9</td><td>0.073</td><td>276.71</td><td>0.35</td></tr><tr><td>Flute (Fl)</td><td>-706.2</td><td>0.076</td><td>-379.2</td><td>0.147</td></tr><tr><td>Piano (Pn)</td><td>-813.5</td><td>0.044</td><td>-361.13</td><td>0.112</td></tr><tr><td>Sax + Vn</td><td>-358.71</td><td>0.364</td><td>-27.37</td><td>0.852</td></tr><tr><td>Sax + Vn + Fl</td><td>-268.7</td><td>0.624</td><td>692.4</td><td>3.813</td></tr></tbody></table>

(right). We can see that performances in both signal reconstruction and transfer decrease with the number of instruments, as the complexity of the incoming signal increases. Both reconstruction and signal-to-transfer scores are almost perfect in the case of solo instruments, providing convincing and high-quality sound samples generation. In the case of mixtures of two or more instruments, reconstruction scores maintain an acceptable performance, but symbol-to-signal transfer scores clearly decrease. This observation correlates with the decrease of performance observed in the symbolic domain, as discussed in the following sub-section.

**Symbolical inference.** Results obtained in the symbolic domain are displayed table 3.8. We provide three different evaluations: the log-density of the generated symbol, the classification ratio performed by the symbolical generator, and then the corresponding classification scores obtained from our baseline. The results within parenthesis are the scores achieved by transfer, when available. In the multi-instrument case we also provide a *loose classification ratio*, that consider a given label correct regardless of the corresponding instrument.

In the single-instrument case, both generation and transfer scores are near to perfection, proving the efficiency of the proposed method. Our model also clearly outperforms the baseline model, especially with dynamics. However, the performances degrade with the number of
<table><caption>**Table 3.8:** Symbolic inference reconstruction and classification results (succesively pitch, octave and dynamics). Scores without parenthesis are reconstruction scores obtained within the symbolic domain, while scores in parenthesis are obtained when performing transfer from the signal domain</caption><thead><tr><th></th><th></th><th>-log p(y|z)</th><th>Success Ratio (%)</th><th>loose (%)</th><th>Baseline (loose)</th></tr></thead><tbody><tr><td rowspan="3">Sax</td><td>p</td><td>-1.0 (-1.0)</td><td>100% (100%)</td><td>-</td><td>94%</td></tr><tr><td>o</td><td>-1.0 (-1.0)</td><td>100% (100%)</td><td>-</td><td>97%</td></tr><tr><td>d</td><td>-1.0 (-1.0)</td><td>100% (100%)</td><td>-</td><td>46%</td></tr><tr><td rowspan="3">Vn</td><td>p</td><td>-1.0 (-1.0)</td><td>100% (100%)</td><td>-</td><td>89.8%</td></tr><tr><td>o</td><td>-1.0 (-1.0)</td><td>100% (100%)</td><td>-</td><td>99.0%</td></tr><tr><td>d</td><td>-1.0 (-1.0)</td><td>100% (100%)</td><td>-</td><td>35.3%</td></tr><tr><td rowspan="3">TpC</td><td>p</td><td>-1.0 (-1.0)</td><td>99.9% (100%)</td><td>-</td><td>76.1%</td></tr><tr><td>o</td><td>-1.0 (-1.0)</td><td>100% (100%)</td><td>-</td><td>99.8%</td></tr><tr><td>d</td><td>-0.998 (-1.0)</td><td>99.7% (100%)</td><td>-</td><td>47.8%</td></tr><tr><td rowspan="3">Fl</td><td>p</td><td>-1.0 (-1.0)</td><td>100% (100%)</td><td>-</td><td>52.4%</td></tr><tr><td>o</td><td>-1.0 (-1.0)</td><td>100% (100%)</td><td>-</td><td>81.8%</td></tr><tr><td>d</td><td>-1.0 (-1.0)</td><td>100% (100%)</td><td>-</td><td>41.4%</td></tr><tr><td rowspan="3">Pn</td><td>p</td><td>-1.0 (-1.0)</td><td>100% (100%)</td><td>-</td><td>51.6%</td></tr><tr><td>o</td><td>-1.0 (-1.0)</td><td>100% (100%)</td><td>-</td><td>63.9%</td></tr><tr><td>d</td><td>-1.0 (-0.999)</td><td>99.9% (100.0%)</td><td>-</td><td>40.0%</td></tr><tr><td rowspan="3">Sax + Vn</td><td>p</td><td>-0.534 (-0.871)</td><td>54.0% (87.9%)</td><td>62.6% (81.6%)</td><td>65.3%</td></tr><tr><td>o</td><td>-0.782 (-0.980)</td><td>84.6% (99.2%)</td><td>94.9% (88.7%)</td><td>79.1%</td></tr><tr><td>d</td><td>-0.712 (-0.939)</td><td>74.4% (95.9%)</td><td>82.4% (66.3%)</td><td>52.0%</td></tr><tr><td rowspan="3">Sax + Vn + TpC</td><td>p</td><td>-0.381 (-0.725)</td><td>38.6% (75.0%)</td><td>62.6% (84.5%)</td><td>56.6%</td></tr><tr><td>o</td><td>-0.377 (-0.641)</td><td>42.4% (67.8%)</td><td>79.3% (88.7%)</td><td>62.3%</td></tr><tr><td>d</td><td>-0.347 (-0.616)</td><td>34.6% (62.4%)</td><td>66.9% (69.5%)</td><td>41.2%</td></tr></tbody></table>

involved instruments, being still acceptable in the two-instruments mixtures case but clearly insufficient with more. This can be explained by observing the loose ratio classification scores: indeed, the results are significantly increased if we compare the obtained symbols regardless of the instrument, even outperforming the baseline. This shows that the VAE is facing a *permutation* problem, where the accurate symbols are retrieved from the mixture but given to the wrong instrument. This issue then prevent the actual formulation of the problem to be applied on instrument mixture; we leave that to a future work, restraining ourselves for the moment to monophonic signals.

<table><caption>**Table 3.9:** Results obtained with the validation flute dataset.</caption><thead><tr><th>Likelihood</th><th>ISD</th><th>Class. Ratio</th><th>Baseline</th></tr></thead><tbody><tr><td rowspan="2">2648 (1057)</td><td rowspan="2">1.065 (0.632)</td><td>65.4%</td><td>63.8%</td></tr><tr><td>81.9%</td><td>76.8%</td></tr></tbody></table>

**Validation results.** We then evaluate the performance of our approach on a validation dataset, to challenge the robustness of the translation processed by our method when facing unknown data. This validation step is important, as checking if the system still performs accurately on different instruments and/or recordings ensures that the extracted features are meaningful. To this end, we take the annotated flute recordings dataset given by Cantos & al. [305], and evaluate the performance of our model trained on monophonic flute signals from the *acidsInstruments-ordinario*. As the annotations provided by the dataset are given in a MIDI format, we then extracted a piano-roll representation that we converted into a symbolic activation matrix, where we sliced the time scale accordingly to the time positions of the corresponding audio spectrum. Unfortunately, dynamics annotations were not provided, so we could not evaluate the robustness of the corresponding inference. Results are shown table 3.9, among with those obtained with the corresponding baseline.

We can observe the achieved performances are still decent, proving that our model is able to generalize on new datasets. More interestingly, the performances obtained in the signal domain are better in transfer
<figure><img src="image_66.png" /><figcaption>**Figure 3.24:** Functional graph of the symbol-signal translation system.</figcaption></figure>

<figure><img src="image_67.png" /></figure>

than reconstruction, showing that our system successfully managed to "abstract" the learned symbolic information. Also, the obtained results have to be contextualized on two points. First, the reconstruction results in table 3.9 have been obtained by forwarding full sounds into the model, while the system was only trained on the stationary parts of the recording. Then, an analysis of spectrograms showed that the flute used in the validation dataset presented a stronger octave-lower harmonic that is not present in the *acidsInstruments-ordinario* dataset, and may explain the drop in octave classification. However, these problems have to be alleviated in order to develop robust signal-symbol translation with this method, with more complete datasets involving several instruments, and training on the full range of the signal. Furthermore, a more complete approach would be to model full sounds rather than single spectral bins, and then model more accurately temporal features such as dynamics; this precise point motivates the models developed in the chapter 4 of this PhD.

**Creative applications.** In addition to allow efficient multi-modal modeling, this system allows additional creative application than with uni-modal variational auto-encoders. Indeed, this system then allows
the user to use custom symbolic information to regularize the latent
space, and can be then provided explicit control parameters for the
generation of feedback during free exploration. The smooth regular-
ization method performed by domain translation then opens up the
system to novel possible interactions that direct supervision methods
would not allow. We summarize below the new features allowed by
using such framework:

* **free trajectory:** we can directly sample the latent spaces to generate the corresponding spectrum. The symbolic decoder can thus been used to provide interesting feedback to the user, relaxing the lack of interpretability of the latent representation. Examples are provided

* **sequence generation:** the navigation in the latent space can be constrained by an external sequence of symbols, such that the user can navigate the distribution $q(z|y)$ while specifying some of the attributes of the target sounds (see fig. 3.23)

* **spectral morphing:** we can retrieve smooth interpolations between two sounds or two symbols in the latent space, then achieving a *spectral morphing* between these two sources.

Fig. 3.24 provides a graphical summary of the interactions allowed by the system, with example taken from the validation flute datasets. Audio examples of these different generation methods are available on the support page.

## 3.4 Perceptual regularization of orchestral instruments

As we saw in section 3.2, the latent representation provided by the AEVB when trained on *ordinario* orchestral instruments is not immediately interpretable, as pointed out by the obtained perceptual maps 3.19. This difficult interpretability can not only make the direct exploration of the latent space unsatisfactory, but also limit the use of the obtained representation as a high-order feature space for analysis. Hence, we may need an additional regularization constraint to enforce the space to match properties of the human perception, making this exploration more intuitive and the representation more likely to encode relevant features about the audio content of incoming signals. To this end, we propose in this section a regularization based on *timbre*, a seminal notion in audio perception.

### 3.4.1 Audio perception and timbre spaces

The notion of timbre is central in numerous music-related domains such as composition, orchestration, musical analysis, instrument making, or digital synthesis. However, giving a formal definition of timbre is a hard task, as this notion is intrinsically linked to human perception and then vary significantly among individuals and musical practices.

<figure><img src="image_68.png" /><figcaption>**Figure 3.25:** Example of random trajectory navigating the latent space obtained from a model trained on violin.</figcaption></figure>
However, the discovery of digital audio synthesis in the early 60s enabled an experimental approach to timbre understanding, underlining the connection of between the perceived timbre properties of a signal and its spectral envelopes, as investigated by Risset & al. [2] and Grey & al. [306]. The link between spectral content and timbre thus motivated the development of quantitative descriptors based on spectral representations, hand-crafted to correspond given timbre properties such as *brilliance*, that can be described using *spectral centroid*, or the *attack quality*, that can be evaluated using slope time estimation [307]. Alternative methods rather derive descriptions for the *cepstrum* of a signal, that takes the inverse Fourier transform of the spectrum logarithm of the signal, estimating the variation rate [308]. Cepstral analysis thus provide an interesting representation of the spectrum envelope and led to very compressed timbre representation, such as Mel-Frequency Cepstrum Coefficients (MFCC) that are still widely used in speech processing and timbre analysis [309, 310]. Timbre analysis being an interesting bridge between digital signal processing and audio perception, this field has then been extensively investigated and provided an impressive number of different features, up to 300 as MIR toolbox [311].

While the development of these audio descriptors is helpful to evaluate timbre properties of a given signal, we also have to ensure that the proposed descriptors correspond to actual perceptual properties of the human hearing [312, 313]. Such psychological approaches then study the consistency of these descriptors on perception by organizing perceptive tests, aiming to find a correlation between the proposed descriptor and the target timbre property [292, 314]. However, an issue raised by such approaches is the problem of *vocabulary*, that can vary a lot depending on the social and psychological context of the human subject. A way to remove this bias is to rather obtain measures of global dissimilarity of sounds, and then analyse the respective pairwise similarity distances between the different types of sounds [315]. This approach is very interesting has it subsumes the idea of an underlying *perceptive measure* behind similarity judgement, opening the idea of a *continuous space* behind human perception.

### 3.4.2 Defining perceptual regularization and similarity space matching

**Timbre spaces from similarity matrices.** Dissimilarity measurements then led to the idea of *timbre space* equipped from a similarity measure, such that two near points correspond to perceptually similar sounds while two distant ones would be very contrasting. Such spaces can be obtained using *Multi-Dimensional Scaling* (MDS), a method that project a set of points in an Euclidean space provided the pairwise distances between the points [110]. Remembering section 2.2.1, given *N* points and the corresponding pairwise dissimilarity measures *d*<sub>*n*,m*</sub> between the points *x*<sub>*n*</sub> and *x*<sub>*m*</sub>, performing MDS comes back to the optimization

<figure><img src="image_69.png" /><figcaption>**Figure 3.26:** MDS can be leveraged to obtain an Euclidean space, whose distances are enforced to matching a given set of dissimilarity measures.</figcaption></figure>
problem

$$ \min_{x_1, \dots, x_N} \sum_{i \le j} (\|x_i - x_j\| - d_{i,j})^2 \qquad (3.18) $$

to obtain a space $x_1...x_N$ that minimizes the mean-squared error between the dissimilarity measures and the Euclidean measures of each pair of points (see fig.3.26). MDS then allows to extract continuous similarity spaces from human dissimilarity measurements, then providing valuable analysis spaces where the relative properties of acoustical instruments could be easily visualized [316, 317]. Timbre spaces obtained from dissimilarity measurements have been first proposed by Grey, using 16 samples from classical instruments rated by 22 subjects over a scale from 0 (similar) to 1 (dissimilar). Consequently, similar experiments were conducted by Krumhasnl [318] using 21 instruments sounds rated by 9 subjects using a discrete scale, Iverson & al. [319] with 16 samples and 10 subjects that was specifically targeting the role of *dynamics* in timbre perception, McAdams & al. with 18 instruments and 24 instruments [320], and Lakatos & al. with 17 subjects and several sets of instruments [321].

In these studies, the dimensions obtained from MDS were generally semantically annotated after visualization and qualitative analysis of the space, as the *brightness-roughness* two-dimensional timbre space obtained by Grey or the spectral centroid / spectral flux / rise-time space obtained by McAdams & al [320]. A strength of this approach is then to extract continuous timbre spaces directly from similarity measurements, without further structural assumptions on the signals nor using specific vocabularies that could bias the experiments. However, while these spaces are very interesting for analytic purposes, they are also limited because of their static construction, implying the re-computation of the MDS for each additional input data, and does not allow *perceptual inference* of external sounds nor direct generation from the inferred timbre space.

**Latent perceptual regularization.** As we saw in the previous sections, perception of timbre is strongly correlated to the spectral content of the perceived signal. We could then beneficially use the latent spaces extracted with AEVB from classical instrument sounds to infer their perceptual properties. AEVB framework would then allow not only to infer the perceptual properties of instruments unknown from the dissimilarity measurements, but also to generate new spectra directly from the timbre space, that would be regularized with perceptive criteria. To this end, we thus have to find a regularization strategy to enforce our latent space to represent whether the pairwise perceptual dissimilarity between two encoded latent vectors, whether more global properties of the timbre space obtained with MDS. To this end, we propose different methods to enforce the geometry of the latent space z to match some properties of the similarity space T obtained with MDS. To our knowledge, we proposed the first regularization strategies based on external metric constraints.

**Prior regularization.** The idea underlying prior regularization is to encode the perceptual information given by the MDS directly in the prior
$p(\mathbf{z}$ of the latent space, then replacing the uninformative isotropic normal $\mathcal{N}(0, \mathbf{1})$ by an informative prior $\mathcal{N}(\mathbf{z}|\boldsymbol{\mu}_i, \boldsymbol{\sigma}^2)$ where the mean $\boldsymbol{\mu}_i$ and variances $\boldsymbol{\sigma}^2$ are conditioned by some information relating the target instrument. Therefore, we model the distribution of each instrument in the dissimilarity MDS with a normal distribution $\mathcal{N}(\boldsymbol{\mu}(\mathcal{C}_i), \boldsymbol{\sigma}^2(\mathcal{C}_i))$, where the mean $\boldsymbol{\mu}(\mathcal{C}_i)$ is the average and $\boldsymbol{\sigma}(\mathcal{C}_i)$ is the estimated variance of the same corresponding projections. We then condition the prior $p(\mathbf{z})$ on the instrument class during training, then regularizing the inference distribution $q(\mathbf{z}_i|\mathbf{x}_i)$ on the corresponding distribution $p(\mathbf{z}_i)$. This regularization allows to directly embed the information given by the perceptive space $\mathcal{T}$ in the ELBO, and then naturally integrates into the Bayesian formulation of the AEVB.

**$\ell^2$ regularization.** In a similar way to the MDS criterion, the $\ell^2$ strategy regularizes the distance of two inferred latent positions $\mathcal{D}(\mathbf{z}_i, \mathbf{z}_j)$ on their distance in the perceptive space $\mathcal{D}(\mathcal{T}_i, \mathcal{T}_j)$. We thus add another criterion $\mathcal{C}(\mathbf{z})$ to the ELBO $\mathcal{L}(q)$ (2.18)

$$ \mathcal{L}_{\ell^2}(q) = \mathbb{E}_{q(\mathbf{z}|\mathbf{x})}[p(\mathbf{x}|\mathbf{z})] + \beta D_{KL}[q(\mathbf{z}|\mathbf{x}), p(\mathbf{z})] + \alpha \mathcal{C}_{\ell^2}(\mathbf{z}) \quad (3.19) $$

where $\mathcal{C}_{\ell^2}(\mathbf{z}) = \|\mathcal{D}(\mathcal{T}_i, \mathcal{T}_j) - \mathcal{D}(\mathbf{z}_i, \mathbf{z}_j)\|^2$, such that this regularization enforces distances in the latent space to match those of the perceptual space using mean-squared error. The $\alpha$ hyper-parameter is the weight controlling the perceptive regularization weight. However, this criterion does not model the uncertainty over the perceptual distances $\mathcal{D}(\mathcal{T}_i, \mathcal{T}_j)$, and thus may be subject to over-fitting.

**Gaussian regularization.** Extending the idea of $\ell^2$ regularization, we rather assume the distances $\mathcal{D}(\mathcal{T}_i, \mathcal{T}_j)$ as drawn from a normal distribution $\{\mathcal{D}_{\mathcal{T}_i, \mathcal{T}_j}^\mathcal{T}\}_{inst} \sim \mathcal{N}(\boldsymbol{\mu}_{i,j}, \boldsymbol{\sigma}_{i,j})$, where the parameters $\boldsymbol{\mu}_{i,j}$ and $\boldsymbol{\sigma}_{i,j}$ are estimated on the mean and variances of inter-class perceptual distances. This regularization method thus models the uncertainty of the distances $\mathcal{D}(\mathcal{T}_i, \mathcal{T}_j)$. We similarly add this criterion to the ELBO 2.18

$$ \mathcal{L}_{\mathcal{N}}(q) = \mathbb{E}_{q(\mathbf{z}|\mathbf{x})}[p(\mathbf{x}|\mathbf{z})] + \beta D_{KL}[q(\mathbf{z}|\mathbf{x}), p(\mathbf{z})] + \alpha \sum_{i,j} p(\mathcal{D}(\mathbf{z}_i, \mathbf{z}_j)|\boldsymbol{\mu}_t, \boldsymbol{\sigma}_t^2) \quad (3.20) $$

**Student-t regularization.** Pushing further the idea of Gaussian regularization, with the Student-t regularization we also target to model the uncertainty of both mean and variances of the normal distribution $\mathcal{N}(\boldsymbol{\mu}_{i,j}, \boldsymbol{\sigma}_{i,j})$. To this end, we then replace this distribution by a Student-t distribution and, taking inspiration from t-SNE, (see sec. 2.2.1) we rewrite the density probability of the perceptual distance as

$$ \mathcal{D}(\mathcal{T}_i, \mathcal{T}_j) = \frac{(1 + \|\mathbf{z}_i - \mathbf{z}_j\|^2)^{-1}}{\sum_{k \neq l} (1 + \|\mathbf{z}_i - \mathbf{z}_k\|^2)^{-1}} \quad (3.21) $$

and the density probability in the latent space as

$$ \mathcal{D}(\mathbf{z}_i, \mathbf{z}_j) = \frac{\exp(-\|\mathbf{z}_i - \mathbf{z}_j\|^2 / 2\sigma^2)}{\sum_{k \neq l} \exp(-\|\mathbf{z}_i - \mathbf{z}_k\|^2 / 2\sigma_i^2)} \quad (3.22) $$
we then formulate the additional objective $C_{tsne}$ as the $D_{KL}$ between these two distributions, then providing the ELBO

$$ \mathcal{L}_N(q) = \mathbb{E}_{q(\mathbf{z}|\mathbf{x})}[p(\mathbf{x}|\mathbf{z})] + \beta D_{KL}[q(\mathbf{z}|\mathbf{x}), p(\mathbf{z})] + \alpha \sum_{i,j} \mathcal{D}(\mathbf{z}_i, \mathbf{z}_j) \frac{\mathcal{D}(\mathbf{z}_i, \mathbf{z}_j)}{\mathcal{D}(\mathcal{T}_i, \mathcal{T}_j)} \quad (3.23) $$

The three first regularization strategies were accepted to ISMIR2018 conference, where it gained the best presentation award. The further results regarding the student-t regularization, that obtained the better results, have been also accepted in DaFX2018. The corresponding paper is available on arXiv* [322]. The code and corresponding page can be found here: https://github.com/acids-ircam/variational-timbre/.

### 3.4.3 Evaluation, perceptual inference and descriptor-based generation

#### 3.4.3.1 Experiment details.

For this experiment, we then had not only to get a suitable set of audio data, but also the corresponding similarity ratings. In this section we describe precisely the chosen experimental protocol, and evaluate the performance of both the AEVB and the effect of the perceptual regularization, comparing the different methods proposed in the last section.

**Dissimilarity ratings.** For this study, we gathered a collection of perceptual dissimilarity measurements from the different sources cited in the previous sections. We ensured the statistical consistency of the overall ratings by only selecting the instruments that were fairly represented across the different studies, and normalized their respective similarity evaluation scales to get uniform ratings between 0 and 1. Despite these constraints, we managed to obtain a consistent set of 12 instruments (Piano, Cello, Violin, Flute, Clarinet, Trombone, French Horn, English Horn, Oboe, Saxophone and Trumpet), for a total amount of 11845 ratings coming from an overall number of 1217 subjects. We then perform a MSD over *d* dimensions on these normalized perceptual ratings, whose PCA can be visualized fig. 3.27. We can see that the obtained representation is perceptually consistent, gathering family of instruments together (brass, woodbrass, strings, and piano as an outlier), and resembling to the timbre spaces obtained from previous studies.)

<footnote>* 322.</footnote>
<figure><img src="image_70.png" /><figcaption>**Figure 3.27:** Obtained perceptual spaces from 11845 dissimilarity measures, representing averages human perceptive distances between orchestral instruments.</figcaption></figure>

**Audio datasets.** As the dissimilarity measurements have been obtained on classical instruments, we then used the *acidsInstruments-ordinario* dataset collected from the Studio On-Line database (see sec. ) as audio dataset. Indeed, this database contains samples of all the 12 required instruments over the full range of their respective tessiture, and additionally provides for each three dynamical ranges (*pp*, *mf*, and *ff*). Furthermore, as within the scope of this study we only train on individual spectral bins, the amount of data available for each pitch and instruments is comfortably sufficient, ensuring enough variability to overcome over-fitting. The total number of sample then amounts to 2,200 samples, that we randomly split across notes to obtain two different datasets for training and testing, with a respective balance of 90%-10%. Moreover, the additional instruments contained in this database can be projected in the extracted timbre space, allowing us to evaluate the *perceptual inference* performance of our models. Regarding the transformation, we trained our models on four different transforms, computed on signal resampled at 22050Hz. First, we used STFT and DCT time-frequency representations, using a Hamming window of 40ms and an hop size of 10ms. We then used NSGT on three different scales: a Constant-Q scale, using 48 bins per octave, a Mel scale, and an ERB scale of 400 bins, defined over a range from 30 to 11000 Hz. As training is performed on spectral magnitude, we used Griffin-Lim procedures to reconstruct the corresponding phase.

**Models.** In this experiment we used a simple VAE with 64 latent dimensions, using 3-layer MLP with 2000 hidden units. The $\beta$ weight of the $D_{KL}$ regularization was set to 2 to encourage the rate of the representation, scheduled with a *warm-up* procedure on the first 100 epochs. The system is jointly trained with an ADAM gradient-descent algorithm, setting the learning rate to 1e-5. We adopted a two-step training procedure, such that the system is trained without perceptive regularization during the first 5000 epochs. This schedule enforces
the perceptual regularization to adapt to the unsupervised features extracted from the data, ensuring that the additional criterion does not interfere too much with the learning process. We further observed that this two-step procedure was mandatory to the success of the procedure.

### 3.4.3.2 Evaluations.

**Generative evaluations.** We evaluate the generation abilities of the trained systems with two different criteria: we first evaluate the likelihood of the encoded example over the generative model $p(x|z)$, and then the mean-squared error of the encoded signal with the mean of the generated distributions. While the first allows us to validate the probability of the generated signal given a certain uncertainty, the second allows us to evaluate the precision of the maximum a-posteriori estimation of the targeted signal. We also compared the obtained results with two baseline PCA and unregularized auto-encoder, whose architecture can be compared with the proposed model. For clarity, we divided the results table in two: in table 3.10 are written the results obtained across the different regularization strategies, while the results obtained across different transforms are shown fig. .

The results table 3.10 show that methods based on Bayesian formulations, that are the *prior* and *student-t* regularizations, achieve the better likelihood results. Conversely, the $\ell^2$ regularization method is providing the worst likelihood performance, but provide the best reconstruction over the mean-squared error criterion. As will be confirmed in the following sections, these results are coherent with the uncertainty trade-off of Bayesian methods: as $\ell^2$ does not model the uncertainty of the perceptual distances it performs better on individual examples, but has the worse performances when the variance of the generated distribution $p(x|z)$ is taken into account. However, we can see that the reconstructions obtained across the different regularization strategies are all satisfactory.

<table><caption>**Table 3.10:** Generative capabilities evaluated on the log likelihood and reconstruction error over the test set.</caption><thead><tr><th>Method</th><th>log p(x)</th><th>|x - ŷ|<sup>2</sup></th></tr></thead><tbody><tr><td>PCA</td><td>-</td><td>2.2570</td></tr><tr><td>AE</td><td>-1.2008</td><td>1.6223</td></tr><tr><td>VAE</td><td>-2.3443</td><td>0.1593</td></tr><tr><td>Prior</td><td>-2.7143</td><td>0.1883</td></tr><tr><td>ℓ<sup>2</sup></td><td>17.8960</td><td><strong>0.1223</strong></td></tr><tr><td>Gaussian</td><td>0.2894</td><td>0.1749</td></tr><tr><td>Student-T</td><td><strong>-2.8723</strong></td><td>0.1610</td></tr></tbody></table>

**Latent space evaluation** The observation of the respective latent for each representation confirm the assumptions developed in the prece-

<figure><img src="image_71.png" /><figcaption>**Figure 3.28:** PCAs of perceptually regularized latent spaces for the different regularization techniques (a) un-regularized (b) *prior* (b) $\ell^2$ (d) *Student-t*</figcaption></figure>
<table><caption>**Table 3.11:** Discriminative capabilities in classifying *family*, *instrument*, *pitch* and *dynamics* of the test set.</caption><thead><tr><th>Method</th><th>Family</th><th>Instrument</th><th>Pitch</th><th>Dynam.</th></tr></thead><tbody><tr><td>PCA</td><td>0.790</td><td>0.697</td><td>0.167</td><td>0.527</td></tr><tr><td>AE</td><td>0.973</td><td>0.957</td><td>0.936</td><td>0.597</td></tr><tr><td>VAE</td><td>0.978</td><td><strong>0.993</strong></td><td>0.963</td><td>0.941</td></tr><tr><td>Prior</td><td>0.975</td><td>0.991</td><td><strong>0.993</strong></td><td>0.936</td></tr><tr><td>Euclidean</td><td>0.972</td><td>0.990</td><td>0.990</td><td>0.943</td></tr><tr><td>Gaussian</td><td><strong>0.982</strong></td><td>0.991</td><td>0.989</td><td><strong>0.948</strong></td></tr></tbody></table>

dent paragraph. Indeed, the over-fitting behaviour of the $ℓ^2$ is clearly visible, the latent projections being shrunk towards the instruments centroid. Conversely, regularization methods modeling the uncertainty of the perceptual distances and are then providing smoother distributions, a the cost of slightly reducing the MSE performance. Moreover, the prior regularization indicates an *over-regularization* effect, that can be taken similar to the results obtained with vanilla VAEs taking non-informative prior. The results obtained with Gaussian regularization naturally lie between both methods, as it assumes a distribution over the distance spaces but without modeling the uncertainty of its parameters. However, regardless of the chosen regularization, we can see that the structure of the perceptual MDS is successfully enforced onto the latent spaces, proving the validity of our approach.

We further evaluate the representational properties of the extracted latent space by training low-capacity discriminators over a bunch of tasks corresponding to inner criteria of variation of the data : *pitch*, *dynamics*, *instrument*, and *family*. The discriminators all have the same architecture, composed from a single-layer network of 512 ReLU units with batch normalization and softmax non-linearity. Corresponding results are provided in table 3.11. We can see that the performances of every proposed methods are highly satisfactory, and significantly gain from the proposed baselines, especially on dynamics. This indicates that the VAE successfully use its latent representation to express the variety of data, while providing consistency across these different tasks to make the problem linearly separable.

**Perceptual inference.** Unlike pure MDS approaches to timbre space extraction, the proposed models allow to project signals unknown to the model onto the extracted perceptual space, allowing us to achieve *perceptual inference*. We then selected samples from the *acidsInstruments* database not present neither in the training nor test sets (Contrabass, Guitar, Harp, Piccolo, Viola), and projected them onto the perceptual space. We can see that the inferred centroids of the obtained projections are consistent, each instrument being close to ones with similar timbre properties : the piccolo lie in the wind part of the space, close to flute and oboe, viola and contra-bass are projected close, and the guitar is projected close to piano. The location of harp however seem a little hazardous, lying between trombone and english horn. Further investigations are then required to describe furthers the capacities of the model.

<figure><img src="image_72.png" /><figcaption>**Figure 3.29:** Perceptive features inferred by the model from unknown instruments.</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>**Figure 3.30:** Perceptual maps obtained from PCA planes of the model.</figcaption></figure>
<figure><img src="image_74.png" /><figcaption>Figure 3.31: Example of descriptor synthesis</figcaption></figure>

**Topology of audio descriptors.** As summarized section 3.4, some perceptual characteristics of sound have been successfully correlated with *audio descriptors*, that are DSP measures applied on a given slice of spectrum. We can then evaluate the global perceptual properties of the extracted representation by studying the repartition of these descriptors across the latent space. Similarly to section 3.2.2.2, we generate the spectral distributions from the space following a grid sampling procedure of the PCA showed in 3.29, and then compute the corresponding descriptors to obtain a *topology map* of the provided representation. Here we provide topology maps of three different descriptors, whose correlation with definite perceptual attributes is confirmed by human tests : *spectral centroid, spectral bandwidth, spectral flatness, and roll-off*.

The obtained maps are shown fig 3.30, where we split the space in four planes among the x axis on positions –.75, –.25, .25, .75, that we sample on a 50 × 50 grid between [-1; 1] over the y and z dimensions. We can see that the repartition of audio descriptors in the latent space is quite non-linear but smoothly evolving, such as the spectrum are coherently spread across the latent space. Interestingly, sound with high flatness seems to be concentrated around the origin of the representation, that correspond to the latent space with higher probability p(z). Furthermore, this non-linear repartition of the audio descriptors confirms the idea of a non-linear relationship between timbre spaces and the corresponding descriptors. We can see that, contrary to the maps fig. 3.19, the descriptors are more parsimoniously spread, and seem to correspond to the features of the target perceptive space.

**Descriptor-based synthesis.** The smoothness of topology maps previously observed points out that the descriptors are behaving locally in a linear way. Writing a descriptor function as a function of z such that $Desc[z] = Desc[\mu_p(z)]$, where $\mu_p$ is the mean function of the decoder, this observation comes back to assuming the local derivability of $Desc[z]$. Given a starting point $z_0$, we can obtain a rough estimation of this derivative by taking a neighbourhood N around $z_0$, evaluating the descriptor $Desc[z]$, and computing the local evolution
tions Desc[$z_0$] – Desc[$N$]. We can thus derive a path search algorithm to retrieve the optimal path that matches a target descriptor profile $t[n]$, from a point $z_0$. While quite rudimentary, this method allows us to perform descriptor-based synthesis, and provides a very intuitive way of sampling new sounds from the latent space. This algorithm is presented in Alg. 1, and an example is shown fig. 3.31.

<caption>**Algorithm 1:** Descriptor-based path synthesis</caption>

**Data:** space $z$, encoder $q_\phi(z|x)$, decoder $p_\theta(x|z)$
**Data:** origin spectrum $x_0$, target series $t_{1..N}$, descriptor $d$
**Result:** spectral distrib. $S \in \mathbb{R}^{N \times F}$
1 // Find origin position in latent space
2 $z_0 = q_\phi(x_0)$
3 // Evaluate origin descriptor
4 $d_0 = \text{evaluate}(x_0, d)$
5 **for** $i \in [1, N]$ **do**
6 // Latent 3-d neighbourhood of current point
7 $N_i = \text{neighborhood}(z_{i-1})$
8 // Sample and evaluate descriptors
9 $X_i = q_\phi(N_i)$
10 $D_i = \text{evaluate}(X_i, d)$
11 // Compute difference to target
12 $\Delta_i = |(D_i - d_{i-1}) - (t[i] - t[i-1])|^2$
13 // Find next latent point
14 $z_i = \text{argmin}(\Delta_i)$
15 // Decode distribution
16 $S[i] = p_\theta(z_i)$
17 **end**

### 3.4.4 Conclusion

In this section, we then investigated the use of variational auto-encoders and variants (Wasserstein Auto-Encoders, Adversarial Auto-Encoders, Renyi Auto-Encoders) for magnitude spectrum learning. Despite efficient reconstruction abilities, we saw that the unsupervised latent space extracted from spectral features were not organized to correspond main audio factors of data, while seeming to get a smooth (but not sufficient) perceptual organization. Unfortunately, main approaches proposed in the literature to increase the disentangling abilities of the latent representation, such as increasing the $\beta$ parameter, did not seem to be efficient with audio factors of variation.

Therefore, we proposed two different regularization strategies to make variational spaces more compliant to audio features: a method based on external symbol information, that used the latent representation as a translation space between audio data and musical labels, and a method based on perceptual regularization using human similarity experiments. Both approaches significantly increased the representative properties of the latent space for audio analysis, and can also provide smoother and controllable latent topologies for creative uses. Furthermore, we saw that the proposed processes allowed numerous applications, as perceptual inference, constrained generation, descriptor-based syntheses, direct interaction, instrument morphing,
providing a complete analysis-synthesis framework.

However, there is still an important aspect missing in this system : *time*.
Indeed, time is fundamental in sound and music, and so far we were
not able to model the temporal shape of the addressed signals. How-
ever, this spectral approach provided interesting results, such that we
would like to keep this process while modeling time at a higher level.
This is also coherent with the temporal *multi-scale* nature of musical
signals, that can be modelled as stacked levels of temporal hierarchies,
from the spectral bin to the full signal. This approach would also allow
us to directly perform *raw audio signals*, that is still a challenge for most
machine-learning methods. We address time and sequentiality in the
next chapter.
Time and prediction in
generative models

"I am sure that i lost my keys, but when...?"

(30 september 2012)

In the previous section, we investigated the capacity of Auto-Encoding
Variational Bayes to learn single spectral frames obtained from time-
frequency representations, that can be seen as a local snapshot of the
signal frequency content. This allowed us to model incoming signals
on a short time range (about 100ms, taking 2048 samples with a sample
rate of 22050Hz). However, these frames are learned individually, such
that the dynamical structure of each individual is not modeled by the
system. This limitation prevents not only to extract dynamical features
from the representation, but also to represent the true underlying man-
ifold involving the time dimension, that is yet mandatory if we target
to model audio information.

Hence, we have to develop AEVB-inspired models that are able to
model the dynamical aspect of the incoming audio files, such that the
obtained representation also reflect the temporal evolution of a given
signal. However, modeling time in AEVB models requires an addi-
tional step, as these systems are developed on data sets composed
from singletons x rather than sequences {x₁...xₜ}. In this work, we
model this input sequence in the latent space by a trajectory {z₁...zₜ},
hence using the latent representation as a state space. This state space,
that is independent to time, is then used to model the sequence dyn-
namics, corresponding to a higher-level features of much lower dimen-
sionality. Recalling the manifold hypothesis underlying variational
auto-encoders (see sec. 2.2.2), this means that a full sound x₁:T corre-
sponds to a trajectory curve in the manifold underlying all the possible
realizations of the modeled dataset.

However, we still need a way to encode the dynamics of the original
sequence in the system. We address this shortage with two original
contributions, targeting to model the sequence {x_t}t \in 0...T by only mod-
eling the corresponding {z_t}t \in 0...T. The first contribution is to model
the distribution p(z_t<\tau|z_{t-1},\tau), that can be seen as an extended version
of the transition model p(z_t|z_{t-1}) used in state-space models. The mod-
eling of this distribution, that we call the prediction task, allows to learn
the latent trajectories of the corresponding time series, and then to
perform the dynamical modeling of the data. As this prediction task is
performed jointly with inference and generation, the chosen prediction
method will shape the topology of the latent space, performing implic-
itly a temporal regularization of the representation. We then proposed
three different prediction methods, based on separate mathematical
backgrounds : a first method based on contrastive predictive coding, a
representation framework used for slow-feature analysis, a second
based on normalizing flows, modeling deformations of the latent space
across time index t, and a third based on Gaussian processes, a powerful
continuous regression method.

However, while modeling the predictive distribution $p(\mathbf{z}_{t<\tau}|\mathbf{z}_{t\ge\tau})$ allows to model the dynamics of the input sequence directly in the latent space, the extraction of dynamical features from this trajectory is still impossible. Hence, we propose a new model, called ShrubVAE, based on the progressive encoding/decoding of increasing time scales. The interest of the proposed model is twofold : first, it provides a new method to use AEVB for time series modeling, enabling to generate entire sequences using a single latent vector, while preventing the undesired degeneracies that occur with auto-regressive decoders (see ref. 2.3.2.3). Secondly, ShrubVAE leverages multi-layered latent architectures to model latent dynamics of increasing temporal scopes, enabling analysis and generation of multiple time granularities from the lower (single spectral frames) to the highest (in this work, up to 150 frames) level of temporal modeling. Moreover, this consecutive downsampling tasks during inference allow to greatly extend the scope of time steps model by a single latent state, hence also significantly extending the scope of the prediction methods proposed previously.

# 4.1 Time series analysis and prediction

Time series analysis is a seminal field of statistics, that have been extensively studied for decades and gave rise to numerous diverse approaches. As we will not be able to give an exhaustive state of the art of such a large research field, we will focus on models that are related to our purpose. First, we will shortly summarize stochastic processes, that are widely used in the audio signal processing (especially auto-regressive processes, that are related to linear systems modeling). Then, we will review dynamical Bayesian systems, leveraging statistical Bayesian learning (see sec. 2.1.1) for sequence analysis, that provided popular models such as Hidden Markov Models and State-Space Models. Finally, we will review Gaussian Processes, a powerful Bayesian modeling method that can be used for regression or classification, that we use in one of our prediction method.

## 4.1.1 Stochastic processes and model identification

Formally, a time series is a sequence of dimensional vectors $\mathbf{x} = [\mathbf{x}_1, \mathbf{x}_2, \dots, \mathbf{x}_N]$, here defined on $\mathbb{R}^D$, generally assumed to be taken over regular time steps. Time series analysis consists in extracting statistics or modeling the underlying process of one or multiple time series, generally assuming an *underlying structure* between time steps, that can be used for feature extraction or prediction. The complexity of this task gave rise to many different methods, targeting distinct applications and series types. By example, some approaches propose to model time series with a parametric function of $t$ such that $\mathbf{x}(t) = f_{\theta}(t)$, coming back to *regression* techniques (see 4.1.3). Alternatively, a flourishing class
<figure><img src="image_75.png" /><figcaption>**Figure 4.1:** Example of a continuous Wiener process, that can be obtained from a cumulative sum of centered normal distributions.</figcaption></figure>

of models rather model the transition function $x_t = f_\theta(x_{t-1})$, obtained from initial conditions. Such approaches are much more powerful, as they can be understood as the discrete identification of a differential equation of arbitrary order (more details are given next paragraph).

**Stochastic processes.** Time series can also be defined using probabilities with the notion of stochastic process, formally defined as a collection of random variables {$x(t)}_{t \in T}$ defined on an index set T, that can also be taken stochastic. In that case, a given signal is considered as a realization or sample path of the corresponding process. The notion of stochastic process is central in time series analysis, grounding the probabilistic approach to time series analysis. Particularly, the Wiener process has a special role in stochastic analysis, that is normally distributed noise continuous everywhere but differentiable nowhere (see fig. 4.1). A Wiener process, sampled with discrete index set T = $\mathbb{N}$, is then called gaussian noise, and is often added to deterministic system to model them using probabilities. (see fig. 4.1)

An interesting property of Wiener processes is that their mean $\mu = E[x(t)] = 0$ (or first-order moment) is constant, and that their correlation function $E[(x_t - \mu_t)(x_u - \mu_u)]$ (or second-order moment) only dependant from the distance $|u - t|$. Furthermore, a Wiener process verifies $p(x_{t1}, x_{t2}, x_{t_n}) = p(x_{t1+\tau}, x_{t2+\tau}, x_{t_n+\tau})$ for any $\tau \in \mathbb{R}$, verifying the stationarity of the process. Stationary is a seminal property of stochastic processes, meaning that its properties do not change over time, and therefore can be learned by a given model. As stationarity is a strong assertion, weaker definitions can be given, as the wide-sense stationarity that only requires the first- and second-order moments of the process to be constant. The concept of stationary is of central importance in time series analysis, as is it conditions the learnability of the process. However, it may be too restrictive to model complex signals, such that some stochastic models propose to bypass this assumption by proposing alternative formulations.

#### 4.1.1.1 Moving-Average and Auto-Regressive processes.

The two notions of moving-average auto-regressive process are seminal in time series analysis, and furthermore find a natural application in digital audio processing for system identification. A moving-average process can be described

$$x_t = \mu + c + \epsilon_t + \sum_{i=0}^{p} b_i \epsilon_{t-i} \quad (4.1)$$

<figure><img src="image_76.png" /><figcaption>**Figure 4.2:** Graph of a probabilistic auto-regressive process</figcaption></figure>

The model parameters are thus the p coefficients $b_i$, then defining a MA-process $MA(p)$ of order p. Analysing a time series with MA then implies to identify the parameters b, that can understood as a local linear regression of the new step under a fixed temporal scope p. A realization of such process can be then generated by iteratively applying
the equation 4.1 with a given set of initial conditions $[x_0...x_p]$. However, MA-processes have a limited expressiveness, and the identifying of the coefficients is non-trivial, as we cannot access the stochastic components $\epsilon_i$. Another family of processes, called *auto-regressive processes*, are defined

$$x_t = c + \epsilon_t + \sum_{i=0}^{p} a_i x_{t-i} \quad (4.2)$$

where the parameters are the coefficients $a_i$, defining an AR-process $AR(p)$ of order $p$. AR-processes are very expressive processes that are widely used in practical applications for time series modeling, the mean-squared error between the modelled and the target signals decreasing with order $p$. The coefficients **a** can be easily retrieved using Yule-Walker equations, or using the estimated auto-correlation function of the signal

$$R_{ff}(\tau) = \sum_{n \in \mathbb{Z}} x[n]x[n-l]$$

the describes the signal auto-correlation depending on lag $\tau$. Furthermore, the auto-correlation function is an efficient method to detect periodicities in the signal, and then to estimate the optimal order $p$ necessary to accurately model the signal $x[n]$ without over-fitting (see fig. 4.3). Alternative schemes, such as *Akaike-Information Criterion* or *Minimum Description Length*, provide information-theory criteria that can be used to optimally retrieve the optimal order $p$ of the model.

The two processes 4.1 and 4.2 can be combined, giving the following *Auto-Regressive Moving-Average* (ARMA) process

$$x_t = c + \epsilon_t + \sum_{i=0}^{p} a_i x_{t-i} + \sum_{i=0}^{m} b_i \epsilon_{t-i} \quad (4.3)$$

whose parameters **a** and **b** can be fully extracted with the Box-Jenkins method. ARMA processes still provide a strong baseline for time series analysis, and are still widely used in numerous domains. Further evolutions of this models can be found in literature, such that *Autoregressive Conditional Heteroskedasticity* models that updates the variance of stochastic components based on previous errors, or *auto-regressive modules with exogeneous inputs* that allows to model ARMA processes depending on external inputs.

An interesting property of ARMA models is their natural incidence in digital processing. Indeed, the equation 4.3 can be formulated as the *difference equation* of an infinite impulse response (IIR) digital filter, plus an additional source of noise. Hence, identifying the coefficients $\{a, b\}$ of an ARMA process can be thought as finding the coefficients of

<figure><img src="image_77.png" /><figcaption>**Figure 4.3:** Auto-correlation allow to give estimates of the intern regularities of the signal (here, original frequencies in red)</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>**Figure 4.4:** Digital filter equivalent to ARMA model 4.3, where $z^{-1}$ is a unit delay. Poles and zeros can be extracted from coefficients $a_0...a_N$ and $b_1...b_M$, representing the frequency response on the unit circle.</figcaption></figure>
the digital filter (see fig. 4.4)

$$H(z) = \frac{b_0 + \sum_{i=1}^{p} b_i z^{-i}}{a_0 + \sum_{i=1}^{q} a_i z^{-i}} \quad (4.4)$$

$$= \frac{\prod_{i=0}^{q} (1 - q_i z^{-1})}{\prod_{i=0}^{p} (1 - p_i z^{-1})} \quad (4.5)$$

where $q_i$ and $p_i$ are respectively the zeros and poles of the filter, the zeros being given by the MA coefficients and the poles by the AR coefficients. ARMA thus provides a stochastic framework for linear model identification, such as the signal was obtained by forwarding a white noise in the linear system identified by coefficients **a**, **b**. Model identification can be whether performed by averaging these coefficients on windows of data, assuming the ergodicity of the signal, or successively updated to provide a representation of the signal across time. The second method, called *Linear Predictive Coding* (LPC) [323, 324], is a well-known representation for audio source-filter separation, and can be used to speech processing [325] and harmonic-noise separation [326]. Furthermore, LPC can be understood as an estimation of the Discrete Fourier Transform of size $T$ whose approximation is optimal when $p \to T$, under some assumed noise and up to an amplitude coefficient proportional to the variance of the noise. Auto-regressive approaches are thus fundamental in digital audio processing, and still provide the basis of more recent algorithms as WaveNet (see sec. 4.2.1).

### 4.1.2 State-space models.

Auto-regressive approaches are powerful methods for time series analysis, performing model identification at the data level. However, the structure of the underlying process strongly conditions the model identification, and does not allow to model the uncertainty of the extracted parameters. Alternatively, some methods rather model the observed data $x_t$ as being generated from an unknown system with an hidden state $z_t$, taking inspiration from continuous state-spaces representations in automation theory that model

$$\begin{cases} \frac{d\mathbf{z}}{dt} = \mathbf{F}(t)\mathbf{z}(t) + \mathbf{B}(t)\mathbf{u}(t) \\ \mathbf{x}(t) = \mathbf{H}(t)\mathbf{z}(t) \end{cases}$$

where $\mathbf{u}(t))$ is an external action modifying the system state at time $t$. In discrete time, the continuous derivative can be replaced with a first-order difference, such we aim to model the transition from $\mathbf{z}_t$ to $\mathbf{z}_{t-1}$. If some normal noise $\epsilon_t$ is added to both transition and observation models, and all functions as linear transforms, we obtain the *Kalman filter*
$$ egin{cases} \mathbf{z}_t = \mathbf{F}_k \mathbf{z}_{t-1} + \mathbf{B}_t \mathbf{u}_t + \Delta_t^z \sigma_z^2 \epsilon_t \text{ transition model} \\ \mathbf{x}_t = \mathbf{H}_t \mathbf{z}_t + \Delta_t^z \sigma_z^2 \epsilon_t \text{ observation model} \end{cases} \quad (4.7) $$

that is a classical probabilistic formulation of state-space models (SSM, see fig. 4.5). Provided the linearity of transition and observation functions, the system is fully tractable with Bayesian inference, assuming the following generative model

$$ egin{align*} p(\mathbf{z}_0) &= \mathcal{N}(\mu_0, \sigma_0^2) \\ p(\mathbf{z}_t | \mathbf{z}_{t-1}) &= \mathcal{N}(\mu_{z,t}(\mathbf{z}_{t-1}, \mathbf{u}_t, \Delta_t), \sigma_{z,t}(\mathbf{z}_{t-1}, \mathbf{u}_t, \Delta_t^z)) \\ p(\mathbf{x}_t | \mathbf{z}_t) &= \mathcal{N}(\mu_{x,t}(\mathbf{z}_{t-1}), \Delta_t^x) \end{align*} $$

hence inferring the posterior $p(\mathbf{z}_t|\mathbf{z}_{t-1}, \mathbf{x}_t)$. Kalman filters are then an example of *dynamic Bayesian learning*, that aims to perform Bayesian inference to model the transition function. The availability of such models to accurately describe the target time series is based on Markov hypothesis, that assumes

$$ p(\mathbf{x}_t|\mathbf{x}_{t-1}, \mathbf{x}_{t-2}, \dots, \mathbf{x}_0) = p(\mathbf{x}_t|\mathbf{x}_{t-1}) \quad (4.8) $$

such that the corresponding process is memoryless. Such processes, then called *Markov chains*, ensure the Bayesian tractability of dynamical models. Despite this strong assumption, Bayesian dynamical models can be expressive enough to model complex relationships on probability graphs of arbitrary shape. *Hidden Markov Models* (HMM) [327] is an example of dynamical Bayesian learning, modeling the dependencies between the series time steps with higher-level latent variables (see fig. 4.6. Unlike state-space models, the underlying Bayesian graph of HMMs is not necessarily causal, allowing to infer more expressive models. However, the tractability of HMMs requires discrete latent variables, then estimating the probability $p(\mathbf{x}_t|\mathbf{z}_1 = z_1, \dots, \mathbf{z}_N = z_N) = p(\mathbf{x}_t|\mathbf{z}_N = z_N)$ and a set of transition probabilities. HMMs are then usually trained using Expectation-Maximization or MCMC techniques on the log-likelihood of the data. Similarly to non-dynamical Bayesian systems, variational approaches to HMMs have been proposed [328, 329], allowing to infer continuous latent spaces. However, inference in Hidden Markov Models is quite computationally heavy, and scale to long time series with difficulty.

<figure><img src="image_79.png" /><figcaption>**Figure 4.5:** General graph of a state-space model (SSM), with observed variables x, hidden variables z and actions u</figcaption></figure>

<figure><img src="image_80.png" /><figcaption>**Figure 4.6:** Example of a *Hidden Markov Model*, a generalized acyclic graph where temporality is driven in hidden discrete latent variables.</figcaption></figure>

### 4.1.3 Gaussian processes.

**Formulation.** Finally, *Gaussian Processes* is a commonly used method for continuous time regression over finite observations. While classical regression models use a finite number of K basis functions $\{\phi_1(x), \phi_2(x), \dots, \phi_K(x)\}$, Gaussian processes can be proven to be equivalent to a regression with an infinite basis of functions, and is thus appealing because of its expressive power. Given a collection of random variables $\mathbf{X} = \{(t_i, \mathbf{x}_i)\}_{i=1}^N$,
<figure><img src="image_81.png" /><figcaption>**Figure 4.7:** Example of Gaussian Process regression performed on a noisy curve. The input signal is the noisy black line, the bold red line is the MAP posterior of the GP, and all the other little lines are functions drawn from the GP.</figcaption></figure>

where $t$ is the continuous time index and $\mathbf{x}$ the corresponding output, Gaussian processes define a probability density over the function space $\mathbf{f}$ defined by $\mathbf{x} = f(t) + \sigma_n^2 + \epsilon_t$ [330]

$$p(\mathbf{f}|\mathbf{X}) = \mathcal{N}(\mathbf{f}|\mathbf{m}(\mathbf{X}), k(\mathbf{X}, \mathbf{X})) \quad (4.9)$$

where the mean vector can be taken as the empirical mean $\mathbf{m}(\mathbf{X}) = \frac{1}{N} \sum_i x_i$ and $k(\cdot, \cdot)$ is the kernel function defining the process. Hence, provided new observations $\{(t_i^*, \mathbf{x}_i^*)\}_{i=1}^N$, we can tractably obtain the Bayesian posterior distribution over $\mathbf{f}$, that is

$$\begin{cases} \mu^* = k(t^*, t)(K + \sigma_n^2 \mathbb{I})^{-1}\mathbf{x} \\ \sigma^* = k(t^*, t^*) - k(t^*, t)^\top(K + \sigma_n^2 \mathbb{I})^{-1}k(t^*, t) \end{cases} \quad (4.10)$$

such that we can update the prior $\mathcal{GP}(\mathbf{f})$ after each new arriving values (see fig. 4.7). The ability of Gaussian processes to perform continuous-time regression on $N$ basis on function is based on the kernel trick, that models the $N \times N$ correlation in input space using a kernel function, rather than in the $K \times K$ feature space. Several kernels $k(\cdot, \cdot)$ are provided in the literature, that can be also parametric and then learned during training on some criteria. In this work we use a simple radial basis kernel $k(\mathbf{x}, \mathbf{x}') = \frac{|\mathbf{x} - \mathbf{x}'|^2}{\ell^2}$ where the parameter $\ell$ is the scale of the kernel.

**Sparse methods.** While Gaussian processes is an efficient framework for continuous-time regression, it unfortunately suffers from scalability problems, as the complexity of the kernel matrix grows with the number $N$ of examples. As the posterior computation requires the inversion of a $N \times N$ matrix, the complexity of this regression is $\mathcal{O}(n^3)$ and can thus become prohibitive for big datasets. A method to alleviate this heavy computation is then to perform the regression on a reduced number of points, allowing to enlighten the heavy inversion of $K_{xx}$. Several methods have then been proposed to reduce this computational cost, such as *Projected Process Approximation* that approximates the covariance matrix with a matrix $Q_{xx} = \text{diag}[K_{xx} - K_{xu}K_{uu}^{-1}K_{ux}]$, where $u$ are $M$ inducing points [331]. Alternatively, *Sparse Posterior Gaussian Process* use a matrix $Q_{xx} = \text{diag}[K_{xx} - K_{xu}K_{uu}^{-1}K_{ux}] + K_{xu}K_{uu}^{-1}K_{ux}$ enforcing the Nystöm approximation to be exact on the diagonal. These methods allow to reduce the complexity up to $\mathcal{O}[nm^2]$, but do not ensure finding the optimal inducing points $u$ that minimizes the approximation error.

Alternatively, Titsias propose to use variational methods to approximate a light Gaussian process $\mathbf{f}_q$ depending on independent auxiliary inputs $\mathbf{u}$, whose variational distribution is defined $q(\mathbf{f}_q, \mathbf{f}) = p(\mathbf{f}|\mathbf{f}_q)q(\mathbf{f}_m)$ [332]. The full model is then trained to reduce the distance between approximated and true posteriors, that is $D_{KL}[q(\mathbf{f}, \mathbf{f}_q)\|q(\mathbf{f}, \mathbf{f}_q, \mathbf{x})]$. This variational approximation performs significantly better with a lower number of points than previous sparse methods, and is then a common scheme for scaling GP to heavy datasets. Sparse GP methods can also advantageously be used in Bayesian Neural Networks [333],
allowing to model a GP prior over the network parameters, or used to perform probabilistic non-linear PCA [334].

**GP for SSM modeling.** As the linearity of the functions involved in usual state-space formulations prevents them to model expressive relations between the observations **x** and the hidden states **z**, Gaussian processes can provide a solution to model richer relationships without dropping the tractability of the posterior. Wang & al. then propose a GP formulation for both observation and dynamical models, using a radial basis kernel for modeling $p(y|\mathbf{x}) \sim \mathcal{GP}(0, \mathbf{K}_{zz})$ and a transition model [335]

$$p(\mathbf{x}|\mathbf{f}) = p(\mathbf{x}_1) \int \prod_{t=2}^{N} p(\mathbf{x}_t|\mathbf{x}_{t-1}, \mathbf{f}_t)p(\mathbf{f}_t)d\mathbf{f}$$

using a linear + RBF kernel $k_z = k_{RBF}(\mathbf{z}, \mathbf{z}') + \alpha \mathbf{x}^T \mathbf{x}$ to add an auto-regressive component to the prediction. However, the method proposed by Wang & al. have the same drawbacks than standard GP when applied to big datasets, thus limiting their applications. A more flexible method is provided by Frigola & al., leveraging variational learning to infer the optimal GP inputs **u** [56]

$$p(\mathbf{x}, \mathbf{z}, \mathbf{f}, \mathbf{u}) = p(\mathbf{x}_0) \prod_{t=1}^{T} p(y_t|\mathbf{x}_t)p(\mathbf{x}_t|\mathbf{f}_t)p(\mathbf{f}_t|\mathbf{f}_{1:t-1}, \mathbf{x}_{0:t-1}, \mathbf{u})$$

$$q(\mathbf{z}, \mathbf{f}, \mathbf{u}) = q(\mathbf{u})q(\mathbf{z}) \prod_{t=1}^{T} p(\mathbf{f}_t|\mathbf{f}_{1:t-1}, \mathbf{x}_{0:t-1}, \mathbf{u})$$

that is fully tractable and trainable with stochastic variational inference, and can be sampled using MCMC particles. However, in this model the approximated transition function $q(\mathbf{z}|\mathbf{z}_{-1})$ is not explicitly modelled, and the MCMC-based particle sampling procedure can be heavy to run. Subsequent similar models alleviate this constraint : Eleftheriadis & al. explicitly model $q(\mathbf{z}|\mathbf{z}_{-1})$ using LSTM networks (see sec. 4.2.1) allowing to model expressive transition functions and to train the model using black-box VI methods [336]. Doerr & al., rather suggest to use the real transition model $p(\mathbf{z}|\mathbf{z}_{-1})$, and also use stochastic gradient descent to train the overall model [337]. Mattos & al. also propose a hierarchical state-space framework using Gaussian processes, also resorting to variational GP sparse approximations [338].

All of these recent advances in GP processing provided powerful models for Bayesian time series analysis, and have been proven to model complex behaviours. However, their scalability to large and multi-dimensional datasets is still very demanding, and GPs are still designed to model only one sequence, being deprived of a conditioning latent model. However, its appealing properties are still encouraging researchers to work develop similar approaches, and GPs will ground one of the proposed prediction methods.
## 4.2 Neural methods for time series modeling

The approaches developed in the previous sections are based on the modeling of time series according to explicit structures, whether based on model identification or likelihood methods using tractable Bayesian posteriors. However, the development of connectionist approaches also allowed to address time series analysis with specific neural networks using gradient descent methods, then providing an alternative approach to this domain. These methods have been now extensively investigated in the literature, and motivated the development of numerous frameworks that are whether entirely based on neural processing, or based on hybrid models that use both the flexibility of these methods and the constrained architectures of previous modeling approaches.

### 4.2.1 Recurrent Neural Networks

Neural networks, that have been briefly introduced section 2.2.2, are powerful function approximators that can be trained with given loss criterion, using a stochastic gradient-descent optimization algorithm. However, standard neural networks are not able to naturally handle time series, and have to be given a *memory* to keep information about previous states. Such networks, called *Recurrent Neural Networks*, are embedded with an hidden memory $\mathit{mathbf{fh}}_{t-1}$ updated at each time-step [339]. The classical formulation for RNNs is the one formulated by Elman [340]

$$h_t = \tanh(\mathbf{W}_{ih}x_t + \mathbf{W}_{hh}z_t) \qquad (4.11)$$

where the hidden representation given at time $t$ also depends on the last output by adding a transition matrix $\mathbf{W}_{hh}$ to the standard neuron (see fig. 4.8). RNNs can thus be thought as modeling a non-linear Markov Chain that learning its transition matrix. Despite the simplicity of this method, training RNNs over long time scales is very difficult, as the iterative update for each sequence step can whether make the gradient vanish, preventing the system to learn, whether to explode, providing degenerated functions. Indeed, performing an iterative map on high-capacity parameteric function such as RNNs can show a *chaotic* behaviour, and can then present hysteresis cycles or instability [341]. Several methods have then been developed to regularize RNNs during training such as $\ell_2$ weight regularization, teacher forcing [342], or leaky integration [343]. The professor forcing learning schedule also help to increase the consistency of RNN units, using an adversarial method to match the generations obtained by feedback with those directly generated from the sequence [344]. However, while these methods can help to alleviate their instability, they do not address the problem of learning long-term dependencies. Subsequent models, such as hierarchical RNNs, can be used to model longer timescales, limiting the scope of individual units to prevent the degeneracy of gradients [345]. Alternatively, a decisive milestone was

<figure><img src="image_82.png" /><figcaption>**Figure 4.8:** Example of an Recurrent Neural Network, that uses an hidden representation updated at each step with a specific operation.</figcaption></figure>
the *Long Short-term Memory* networks by Hochreiter & Schmidhuber [346], allowing the system to *forget* information by gating before and after each subsequent step (see fig. 4.9)

$$ i_t = \sigma_g (W_i x_t + U_i h_{t-1} + b_i) \qquad \text{input gate} \quad (4.12) $$

$$ o_t = \sigma_g (W_o x_t + U_o h_{t-1} + b_o) \qquad \text{output gate} \quad (4.13) $$

$$ c_t = i_t \odot \tanh(W_c x_t + U_c h_{t-1} + b_c) \qquad \text{gated recurrent input} \quad (4.14) $$

$$ h_t = o_t \odot c_t \qquad \text{gated recurrent output} \quad (4.15) $$

An additional *forget gate* can be added [347], such that $c_t^f = f_t \odot c_{t-1}^f + c_t$, where $f_t$ is an additional gating mechanism. These recurrent units are then able to forget data, allowing to free some memory in the retained representation to learn new features during iteration. The LSTM can also be added a reverse process, learning features $h_t$ from $h_{t+1}$, then also learning information from the backward sequence [348].

Gating mechanisms allow to accurately train recurrent neural units over large time spans, and are then commonly used for neural time series modeling. Their efficiency have then been extended to non-temporal purposes, providing efficient gradient selection methods as *attention mechanisms* [349, 350] and *transformers* [351, 352]. Alternative gating mechanisms have been proposed since such as *Gated Recurrent Units* (GRU) [353], that are used in this thesis. Some approaches also perform recurrent modeling in the spectral domain, that have been shown to accurately describe chaotic processes such as Lorentz attractors [354]. Finally, the recently proposed Neural Processes propose to combine neural networks and stochastic processes to infer probabilities directly in the function space, providing a very efficient framework for Bayesian learning of stochastic processes [355–357]. Time series modeling using neural networks is still an active field of research, that is likely to be further improved in the near future.

<figure><img src="image_83.png" /><figcaption>Figure 4.9: Graph of LSTM unit. Layer here is a sum of each input multiplied by a specific parameter matrix, and the gate is just a sigmoid unit allowing to mask given dimensions with Hadamard product ⊙.</figcaption></figure>

### 4.2.2 Dynamical Bayes and AR-models

Recurrent neural networks are then powerful units to model functions $y_t = f_t(x_t)$ with an independent recurrent mechanism. Hence, neural networks can be used in Dynamical Bayes to model observation or transition functions, empowering the expressiveness of these systems. In that case, neural networks can be used whether to model the transition functions of the latent variables, or directly generating new data samples, conditioned on a definite range of observed steps. Provided the *causality* of the system, meaning that the probability density of a step $x_t$ is conditioned over the *previous* ones $x_{t-1}$, these systems can be used to model the underlying auto-regressive process $x_t = f(x_{t-1}...x_{t-T})$ of the series.
<figure><img src="image_84.png" /><figcaption>Figure 4.10: The Neural Autoregressive Distribution Estimator [358]</figcaption></figure>

<figure><img src="image_85.png" /><figcaption>Figure 4.11: The WaveNet auto-regression module, that transform auto-regressive neural modules (top) using dilated convolutions (down) to extend the temporal scope using less connections.</figcaption></figure>

<figure><img src="image_86.png" /><figcaption>Figure 4.12: SampleRNN hierarchical modeling used to learn the underlying auto-regressive distribution p(x)</figcaption></figure>

The Neural Auto-regressive Distribution Estimator [358, 359] network, proposed by Larochelle & al., aims to model the auto-regressive distribution of a binary data vector $x \in [0; 1]^N$ (see fig. 4.10)

$$p(x_i = 1 | \mathbf{x}_{<i}) = \sigma(b_i + \mathbf{W}_i^\top \mathbf{h}_i) \\ \mathbf{h}_i = \sigma(\mathbf{c} + \mathbf{W}_{<i} \mathbf{x}_{<i})$$

that can be understood as an hybrid between recurrent neural networks, also using a memory $\mathbf{h}$, and an auto-encoder with hidden representation $\mathbf{h}_i$ (see fig. 4.10). As we summarized sec. 4.1.1.1, auto-regressive distributions are very expressive models, and are widely used outside the scope of temporal models. In that case, auto-regressive distributions are used to successively model the dimensions of a multi-dimensional vector $\mathbf{x}_n$, such that $\mathbf{x}_n = \prod_{i}^{D} p(x_i|\mathbf{x}_{<i})$. In that case, the ordering of dimensions $d_1, d_2...d_N$ does not have to be causal, and an arbitrary ordering can be chosen. Based on the NADE model, the Masked Auto-Regressive Distribution Estimator (MADE) rather uses the full matrix $\mathbf{W}$ rather than $\mathbf{W}_{<i}$ during generation, but uses mask to perform automatic dimension selection for the auto-regressive models.

However, auto-regressive models such as NADE and MADE can be slow to converge, and assume weight matrices whose size matches the sequence length, and can then become prohibitive with long-term sequences. Another method for auto-regressive modeling is the WaveNet model, proposed by van den Oord & al., that provides so far the baseline for raw waveform synthesis. WaveNet is based on stacked 1d dilated convolution networks with exponentially growing dilation factors, that significantly increase the receptive field of the context $\mathbf{x}_{i-1}...\mathbf{x}_{-1t}$ used to generate the next step $\mathbf{x}_t$ (see fig. 4.11). The operation processed by a dilated convolution unit is

$$\mathbf{h}_c = \mathbf{b}_c + \sum_{k=0}^{c_{in}-1} \mathbf{w}_{c_{out},k} \star_d \mathbf{x}_k \quad (4.16)$$

$$\mathbf{w}_n \star_d \mathbf{x}_n = \sum_{n} \mathbf{w}_n \mathbf{x}_{d \cdot n+m} \quad (4.17)$$

that can be described as a spaced convolution operation, such that a kernel of $m$ samples will be applied over a time range of $m \times d$ samples of the sequence. The use of dilated convolutions allowed WaveNet to extend its receptive field up to 300 milliseconds, that is sufficient to ensure the perceptual consistency of the generated signal. A similar method is the SampleRNN model proposed by Mehri & al., that uses structured hierarchical RNNs to predict windows of samples, thus fastening the generation process [360] (see fig. 4.12). Both methods are compatible with teacher forcing and professor forcing, that can help the robustness of prediction.
While WaveNet and SampleRNN are the first machine-learning models able to generate convincing audio waveforms, an impressive success regarding the complexity of the challenge, these models rely on heavy conditional architectures that makes their learning tedious. Indeed, these models requires weeks of training on single GPUs, and the generation is still far from being real-time. Haque & al. recently proposed reduce WaveNet architectures using attention LSTMs, fastening a little computation during generation time [361].

In addition to this performance issues, another weakness of such auto-regressive networks is the lack of latent representation inferred from the data, then deprived of higher-level space that could use the extracted features for audio analysis. While systems such as WaveNet can be conditioned during training to allowing class-dependent generation, the number of possible classes is limited, and the scale of the datasets needed to train such algorithms can make their annotation discouraging. A step towards controlling WaveNet decoders have been proposed by Engel & al., based on the extraction of latent features extracted from variational auto-encoders to condition a WaveNet decoder [362]. However, these features are obtained from spectral information, and the do not extract a representation based on the raw information itself. Hence, real-time raw signal synthesis is still a challenging task, and developing models that are able long-term dependencies with an acceptable computational cost is an on-going active research problem.

### 4.2.3 Variational auto-encoding for time-series

AEVB methods for inference/generation of time series were investigated shortly after the publication of the original article. However, addressing time series analysis with AEVB-inspired methods can be divided in two different approach (see fig. 4.13). The first consists in modeling the full sequence {$x_1...x_T$} with a single latent vector **z**, that we call *dynamical compression*. A first approach using VAEs with auto-regressive decoders was investigated by Bowman for language processing, targeting to generate full phrases from a single latent vector[183]. However, this approach unveiled some undesirable properties of AEVB, such as *latent over-pruning*, that showed the inherent trade-off between the capacity of the decoding module and the information represented in the latent space (see sec. 2.3.2.3). Hence, Chen & al. then proposed to rather encode *local* features of the input, that were decoded and then assembled by an additional auto-regressive module. The second approach is rather based on representing the input {$x_1...x_T$} with a latent sequence $z = z_1, ..., z_T$, modeling one different latent position for each step. Hence, this approach rather considers the variational auto-encoder as a state-space model, using the generative distribution $p(x|z)$ as the observation model (see sec. 4.5). Hence, the different solutions proposed in the literature for sequential AEVB will design different structures for the transition model $p(z_t|z_{t-1})$, whether by explicitly modeling this distribution in the variational process, or by modeling recurrence directly in the encoding and decoding modules