Alberto Bietti

Arshia Cont

Francis Bach

Audio segmentation is an essential problem in many audio signal processing tasks which tries to segment an audio signal into homogeneous chunks, or segments. Most current approaches rely on a change-point detection phase for finding segment boundaries, followed by a similarity matching phase which identifies similar segments. In this thesis, we focus instead on joint segmentation and clustering algorithms which solve both tasks simultaneously, through the use of unsupervised learning techniques in sequential models. Hidden Markov and semi-Markov models are a natural choice for this modeling task, and we present their use in the context of audio segmentation. We then explore the use of online learning techniques in sequential models and their application to real-time audio segmentation tasks. We present an existing online EM algorithm for hidden Markov models and extend it to hidden semi-Markov models by introducing a different parameterization of semi-Markov chains. Finally, we develop new online learning algorithms for sequential models based on incremental optimization of surrogate functions.

Résumé

Le problème de la segmentation audio, essentiel dans de nombreuses tâches de traitement du signal audio, cherche à décomposer un signal audio en courts segments de contenu homogène. La plupart des approches courantes en segmentation sont basées sur une phase de détection de rupture qui trouve les limites entre segments, suivie d'une phase de calcul de similarité qui identifie les segments similaires. Dans ce rapport, nous nous intéressons à une approche différente, qui cherche à effectuer les deux tâches -segmentation et clustering -simultanément, avec des méthodes d'apprentissage non supervisé dans des modèles séquentiels. Les modèles de Markov et de semi-Markov cachés sont des choix naturels dans ce contexte de modélisation, et nous présentons leur utilisation en segmentation audio. Nous nous intéressons ensuite à l'utilisation de méthodes d'apprentissage en ligne dans des modèles séquentiels, et leur application à la segmentation audio en temps réel. Nous présentons un modèle existant de online EM pour les modèles de Markov cachés, et l'étendons aux modèles de semi-Markov cachés grâce à une nouvelle paramétrisation des chaines de semi-Markov. Enfin, nous introduisons de nouveaux algorithmes en ligne pour les modèles séquentiels qui s'appuient sur une optimisation incrémentale de fonctions surrogées.

First, I would like to thank my advisors, Arshia Cont and Francis Bach, whose support and expertise have been of invaluable help throughout my internship. I also would like to thank them for giving me the opportunity to conciliate my two main passions -machine learning and music -with this project. Being a musician myself, working at IRCAM has been a wonderful experience. I enjoyed being exposed to the world of contemporary classical music as well as learning about state of the art scientific research on music and audio. I also learned a lot from the various seminars and group meetings I attended in the SIERRA group at INRIA, which has shown to be an incredible environment for machine learning research.

I would like to thank Philippe Cuvillier, with whom I shared many fruitful and insightful discussions. I also want to thank Olivier Cappé and Mathieu Lagrange for helpful discussions.

Chapter 1

Introduction

The task of audio segmentation has been an active area of research in the audio signal processing community since the 1980s. It aims to decompose an audio signal into homogeneous chunks, or segments, such that neighboring segments contain what a listener would consider different audio content. Perhaps the most common use case of audio segmentation is musical note segmentation, where an algorithm attempts to detect the start and end of each note in a musical piece. This can be extended to discovering musical structures, where segments can be more than just single notes, but rather short structures, such as chords or vibratos, which are homogeneous or close to being so. The problem of audio segmentation is however more general than musical segmentation, and could also be applied to segmenting auditory scenes in an audio file. Example applications of audio segmentation include music indexing for information retrieval, audio summarization (where the summary can then be synthesized), fingerprinting, or analyzing higher-level musical structures from the symbolic representation given by the segmentation.

Many approaches to audio segmentation rely on detecting abrupt changes in the audio content (such as musical onsets) which determine the boundaries of each segment, a task known as change detection. This is often accomplished by looking at specific audio features such as spectral centroid, spectral flux or zero crossings [START_REF] Tzanetakis | Multifeature audio segmentation for browsing and annotation[END_REF]. Another important technique for change-point detection is that of sequential hypothesis testing, using e.g. the CUSUM algorithm [START_REF] Tartakovsky | Sequential Analysis: Hypothesis Testing and Changepoint Detection[END_REF][START_REF] Basseville | Detection of abrupt changes: theory and application[END_REF]. Recently, [START_REF] Dessein | An information-geometric approach to real-time audio segmentation[END_REF]; [START_REF] Dessein | Computational Methods of Information Geometry with Real-Time Applications in Audio Signal Processing[END_REF] applied these techniques to real-time audio segmentation using the information geometry of exponential families. These approaches don't generally assign a label to each segment directly, and one usually needs to compute similarity matrices between segments in order to find similar segments and discover the underlying structure, which is not adapted to real-time settings. [START_REF] Lostanlen | Découverte automatique de structures musicales en temps réel par la géométrie de l'information[END_REF]; [START_REF] Cont | On the information geometry of audio streams with applications to similarity computing[END_REF] show how segmentations and similarities can be computed in an information-geometric context by finding the right centroids for each segment.

In many cases, the audio segmentation task is cast to a supervised learning task, in which detecting a change point or a specific segment is done using a classifier which needs to be trained from labeled data. This makes the algorithm depend heavily on the training examples used, and might not easily adapt to new auditory environments.

We will focus on unsupervised learning approaches to audio segmentation, where no training data is required and we rely on the expressiveness of the model to discover the correct segmentation. Rather than separately detecting changes and finding similarities, we wish to do both at the same time, by having different models for each class of segments and inferring the segment transitions from observations. This naturally leads to the framework of hidden Markov models (HMMs) [START_REF] Cappé | Inference in Hidden Markov Models[END_REF][START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF] and hidden semi-Markov models (HSMMs) [START_REF] Yu | Hidden semi-markov models[END_REF], which are powerful modeling tools and have had great success in audio signal processing, particularly speech recognition. In Chapter 2, we introduce the representation we use for audio signals, the use of Bregman divergences for computing similarities, and present the main algorithms for (offline) learning in various models of interest, including HMMs and HSMMs.

Online learning has proven to be an effective way to improve learning, especially in largescale settings [START_REF] Bottou | The tradeoffs of large scale learning[END_REF][START_REF] Bottou | Online learning and stochastic approximations[END_REF]. Most of its successes have been in the context of independent observations, e.g., [START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF][START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF], and little work has been done to apply online learning to sequential models, such as hidden Markov or semi-Markov models. O. Cappé proposed an online EM algorithm for HMMs based on a forward smoothing recursion [START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF]. [START_REF] Yildirim | An online expectation-maximization algorithm for changepoint models[END_REF] proposed an online EM algorithm for changepoint models (including HSMMs), but they rely on sequential Monte Carlo sampling techniques as in [START_REF] Cappé | Online sequential monte carlo EM algorithm[END_REF], which aren't always desirable.

In Chapter 3, we present various online learning algorithms for HMMs and HSMMs. We propose an extension of the online EM algorithm of [START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF] to HSMMs in Section 3.1.2, and various other online algorithms for HMMs (which can also be extended to HSMMs) based on incremental optimization of surrogate functions. Chapter 4 presents some experiments for the audio segmentation task on musical examples.

The online learning algorithms present in Chapter 3 could also be beneficial to the Antescofo score-following system built at IRCAM [START_REF] Cont | A coupled duration-focused architecture for real-time music-to-score alignment[END_REF], which is a real-time audio-to-score alignment system based on inference in a hidden semi-Markov model. At present time, Antescofo involves no learning and its observation model relies on fixed templates artificially created for each note or chord that is observed. Because this observation model is quite weak, the model relies heavily on the duration criteria of the notes on the score, and isn't very robust to changes in sound, due to different instruments or different environments. Our online learning algorithms could be used to learn these observations models adaptively over time, starting from the handcrafted templates as a prior model.

Our main contributions are the following:

• an extension of the online EM algorithm of [START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF] to HSMMs, thanks to a twovariable parameterization of the HSMM

• various algorithms for online learning in HMMs and HSMMs based on incremental optimization schemes, and a comparison of these models

• applications of these algorithms to audio segmentation, and a potential application to the Antescofo score-following system for adaptively learning the observation templates.

Chapter 2

Representations, models and offline algorithms

In this chapter, we will provide some background on the building blocks we use for audio segmentation. Our focus will be on unsupervised learning techniques, where no labeled data is provided and we rely on the richness of the model to correctly identity the segments and musical structures. We will start by discussing the choice of representation for the audio signals. We then present the modeling tools based on Bregman divergences that we use on our representations, and finally we describe Hidden Markov and semi-Markov models, which are well-suited for the sequential nature of our task, and derive appropriate learning algorithms.

Audio signal representation

An audio signal x(t) is commonly described by frequency representations, given by the Fourier transform.

Definition 2.1 (Fourier transform).

x(ω) = +∞ -∞

x(t)e -iωt dt (2.1)

The Fourier transform captures the frequency information of the entire signal, and therefore fails to capture local changes in frequency components. This problem is overcome by the windowed Fourier transform, or short-time Fourier transform (STFT), which uses a real and symmetric window g(t) (typically bell-shaped) to localize the Fourier transform around each time t. Definition 2.2 (Continuous-time STFT).

x(t, ω) = +∞ -∞ x(u)g(u -t)e -iωu du (2.2)
Common choices for the window function g are the Gaussian window, the Hamming window or the Hanning window. If we now consider a discrete signal x [t], where t ∈ Z, we get the following discrete-time transform:

Definition 2.3 (Discrete-time STFT). x(t, e iω ) = +∞ u=-∞ x[u]g[u -t]e -iωu
(2.3)

When g has compact support, each x(t, e iw ) for fixed t can be obtained by taking the discrete Fourier transform of the signal u → x[u]g[ut] on a time window centered around t (typically using an FFT algorithm), and is thus given by a fixed number p of coefficients xt,1 , . . . , xt,p ∈ C.

In the following, we consider the audio signal to be represented by the modulus of its STFT coefficients: x t ∈ R p = (|x t,1 |, . . . , |x t,p |) ⊤ , where the subscript t = 1, . . . , T denotes the index of the time windows -which have a constant offset and can have some overlap -and can potentially be unbounded in the online (streaming audio) case. In order to provide some invariance to the volume of the sound, we will consider the vectors x t to be normalized, i.e. i x t,i = c > 0 for all t.

Note that this representation is quite primitive and more informative representations could be used instead (e.g. MFCCs or scattering coefficients [START_REF] Andén | Deep scattering spectrum[END_REF]) to improve performance, however we will restrict ourselves to this simple representation, and will rely on the richness of our models as a means for improving results.

Bregman divergences

In order to cluster data points together, we need to be able to compare them using a similarity measure D(x, y). In many cases the chosen measure is the Euclidian distance, or its square: D(x, y) = xy 2 . However, the choice of similarity measure makes some implicit assumptions about the geometry of the data, and in particular using the Euclidian distance assumes Euclidian geometry. Empirical evidence shows that alternative measures, such as the Kullback-Leibler (KL) or Itakura-Saito (IS) divergences perform better in practice for audio signals [START_REF] Cont | On the information geometry of audio streams with applications to similarity computing[END_REF][START_REF] Stylianou | Perceptual and objective detection of discontinuities in concatenative speech synthesis[END_REF]. In particular, when dealing with non-negative normalized vectors, the KL divergence is a natural choice since it corresponds to a multinomial distribution, as we will see in Section 2.3.3 (see, e.g., [START_REF] Banerjee | Clustering with bregman divergences[END_REF]).

Definition and examples

Definition 2.4 (Bregman divergence). Let ψ : Ω → R be a strictly convex function defined on a convex set Ω ⊂ R p . The Bregman divergence associated with ψ is defined by

D ψ (x, y) = ψ(x) -ψ(y) -x -y, ∇ψ(y) .
(2.4)

It follows from the strict convexity of ψ that D ψ (x, y) ≥ 0, with equality if and only if x = y. Note that in general Bregman divergences are not symmetric and do not satisfy the triangle inequality, and hence are not distance metrics. We now give some classical examples of Bregman divergences relevant for audio.

Example 2.5. The squared Euclidian distance xy 2 is a Bregman divergence, with ψ(x) = x 2 . Indeed,

D ψ (x, y) = x 2 -y 2 -x -y, 2y = x 2 + y 2 -2 x, y = x -y 2 . Example 2.6. The KL divergence D KL (x y) = i x i log xi yi , with x i , y i ∈ R + and i x i = i y i = c (where c = 1 for probability distributions) is a Bregman divergence, with ψ(x) = i x i log x i . Indeed, we have D ψ (x, y) = i x i log x i - i y i log y i - i (x i -y i )(1 + log y i ) = i x i log x i y i - i x i + i y i = i x i log x i y i = D KL (x y).
Note that if we define ψ as a function of the p -

1-dimensional vector x = [x i ] i∈{1,...,p-1} by ψ(x) = p-1 i=1 x i log x i +(c- p-1 i=1 x i ) log(c- p-1 i=1
x i ), it can be easily shown that the associated Bregman divergence still takes the same form. This formulation has the benefit of having ψ defined on a set with non-empty interior, which will be useful for Legendre duality properties as we will see in §2.2.2 and §2.3.3.

Example 2.7. The Itakura-Saito divergence 1 D IS (x, y) = i ( xi yi -log xi yi -1) with x i , y i > 0 is a Bregman divergence with ψ(x) = -i log x i . We have indeed D ψ (x, y) = - i log x i + i log y i - i (x i -y i ) = i x i y i -log x i y i -1 = D IS (x, y).

Duality

When ψ is differentiable in the interior of Ω and diverges near its boundary, the function ψ and its Fenchel conjugate2 ψ * are so-called Legendre duals of each other [START_REF] Rockafellar | Convex analysis[END_REF][START_REF] Banerjee | Clustering with bregman divergences[END_REF][START_REF] Amari | Methods of Information Geometry[END_REF] and verify ∇ψ * = (∇ψ) -1 . If µ = ∇ψ(x), or equivalently x = ∇ψ * (µ), we have

ψ(x) + ψ * (µ) = x, µ . (2.5)
Then, if we let µ = ∇ψ(x) and µ = ∇ψ(y), we have the following duality relation between D ψ and D ψ * :

D ψ (x, y) = D ψ * (ν, µ).
(2.6)

Clustering with Bregman divergences

Now that we have defined a similarity measure between data points, our goal is to use it to cluster data points together. We will focus on centroid-based clustering methods, where each cluster of points is described by its centroid, the K-means algorithm being the simplest example. We will start by defining and characterizing centroids for Bregman divergences (see, e.g. [START_REF] Nielsen | Sided and symmetrized bregman centroids[END_REF]), then, based on [START_REF] Banerjee | Clustering with bregman divergences[END_REF] we will describe a generalization of K-means to Bregman divergences, as well as extensions to probabilistic models thanks to a bijection between Bregman divergences and exponential families. In this section, we will ignore the temporal dependency of data points, and will consider any set of points {x i } i=1..n .

Bregman centroids

A centroid is defined as the point µ minimizing the average distance to a set of other points {x i } i=1..n . If the distance we consider is the squared Euclidian distance, it is well known that the centroid is given by the average x = 1 n i x i . Since Bregman divergences are not symmetric, one can define three different types of centroids as follows.

Definition 2.8. The right-type, left-type and symmetrized Bregman centroids for the divergence D ψ are defined respectively by

µ R ψ = arg min c 1 n n i=1 D ψ (x i , c) µ L ψ = arg min c 1 n n i=1 D ψ (c, x i ) µ S ψ = arg min c 1 n n i=1 D ψ (x i , c) + D ψ (c, x i ) 2 .
The following proposition characterizes right-type Bregman centroids to be point averages and is a key result for extending K-means to Bregman divergences.

Proposition 2.1. Given a set of points {x i } i=1..n and a Bregman divergence D ψ , the right-type Bregman centroid is given by

µ R ψ = x := 1 n n i=1 x i .
(2.7)

Proof. We have

1 n n i=1 D ψ (x i , c) = 1 n n i=1 (ψ(x i ) -ψ(c) -x i -c, ∇ψ(c) ) = 1 n n i=1 ψ(x i ) -ψ(c) -x -c, ∇ψ(c) = 1 n n i=1 ψ(x i ) -ψ(x) + D ψ (x, c),
which is minimized for D ψ (x, c) = 0, i.e. for c = x.

The proof can easily be extended to the weighted case, and we obtain the following result.

Proposition 2.2. The unique minimizer of the weighted sum i w i D ψ (x i , µ), with w i ≥ 0 and

i w i = 1, is given by µ = i w i x i .
Interestingly, Banerjee et al. (2005a) give a reverse characterization of Bregman divergences as the only "loss" functions, under some mild assumptions, for which this proposition is true.

We will limit ourselves to right-type centroids in the following because of their simplicity. That said, it is possible to easily compute left-type centroids thanks to (2.6) by noting that the left-type centroid under D ψ corresponds to a right-type centroid under D ψ * in the dual space given by the one-to-one mapping x → ∇ψ(x):

µ L ψ = (∇ψ) -1 1 n n i=1 ∇ψ(x i ) .
(2.8)

By switching to the dual representation and changing the divergence to D ψ * , one can thus reformulate a problem with left-type centroids as one with right-type centroids. As for symmetrized Bregman centroids, [START_REF] Nielsen | Sided and symmetrized bregman centroids[END_REF] proposed a dichotomic geodesic-walk algorithm to approximate them, using the differential geometric properties of Bregman divergences.

K-means

The K-means algorithm clusters data points around K centroids by alternating between two steps: (i) assign each point to its closest centroid, (ii) update the centroids from the newly assigned points. Although the classical K-means algorithm uses the squared Euclidian distance xµ 2 as a way to measure distance between a point x and its centroid µ, [START_REF] Banerjee | Clustering with bregman divergences[END_REF] show that the algorithm can be extended to use any Bregman divergence, with the centroid on the right, D ψ (x, µ), to measure the closeness of x to µ. The algorithm is given in Algorithm 2.1, where the names E-step and M-step come from the close resemblance to the EM algorithm, which we will discuss below.

Algorithm 2.1: K-means algorithm.

Data:

x 1 , . . . , x n , initial centroids µ 1 , . . . , µ K , Bregman divergence D ψ . Result: final centroids µ k , assignments z i repeat z i ← arg min k D ψ (x i , µ k ) i = 1, . . . , n /* E-step */ µ k ← 1 |{i:zi=k}| i:zi=k x i k = 1, . . . , K /* M-step */ until convergence;
The following proposition justifies the algorithm and explains its convergence.

Proposition 2.3. The K-means algorithm decreases the following objective (also known as distortion) at each step by performing coordinate descent:

ℓ(µ, z) = n i=1 D ψ (x i , µ zi ).
(2.9)

Furthermore, convergence is reached in a finite number of iterations.

Proof. The E-step minimizes the objective w.r.t. z by minimizing each individual term D ψ (x i , µ zi ). That the M-step minimizes the objective w.r.t. µ follows from Proposition 2.1. Let ℓ 1 (z) = min µ ℓ(µ, z), that is the value of the objective after the M-step. ℓ 1 is non-negative and decreases with each iteration, and therefore converges. Since there is only a finite number of possible assignments, convergence is reached in a finite number of iterations.

Initialization One issue with the K-means objective -and with most of the problems we will deal with -is that it is non-convex, and thus relies heavily on a good initialization. A simple approach is to initialize the centroids to random points in the dataset, and keep the result with lowest distortion after multiple such initializations. K-means++ [START_REF] Arthur | K-means++: The advantages of careful seeding[END_REF] is an example of another initialization scheme.

Choice of K The number of clusters K needs to be chosen in order to run K-means, and the question comes up of how to pick a good value for it. Of course, if one knows in advance how many clusters are present in the data, then the choice is immediate, however such information isn't always available and it is common to run the algorithm for different values of K and evaluate which is best using some measure. In practice, this measure is often taken to be a penalized distortion measure, where a penalty term λK is added to the objective in (2.9). [START_REF] Kulis | Revisiting k-means: New algorithms via bayesian nonparametrics[END_REF] derive an algorithm, called DP-means, which incrementally adds new clusters when all existing clusters are further away than λ from the current point. Their algorithm is justified as an asymptotic Gibbs sampling algorithm in a Dirichlet process mixture of Gaussians (with fixed isotropic covariances), where the variance goes to zero, and is shown to decrease the penalized objective until convergence. Although the link with Bayesian nonparametrics made in [START_REF] Kulis | Revisiting k-means: New algorithms via bayesian nonparametrics[END_REF] relies on the Gaussian distribution and thus the Euclidian distance, the same algorithm with a Bregman divergence still satisfies the property of decreasing the penalized objective.

Bregman divergences and exponential families

Although K-means provides a simple and effective clustering algorithm, it lacks the flexibility of a probabilistic model and the ability to have soft-assignments of a data point to multiple clusters. Mixture models such as Gaussian mixture models are widely used to address these issues, and [START_REF] Banerjee | Clustering with bregman divergences[END_REF] show how one can use such models with Bregman divergences by establishing a bijection between Bregman divergences and exponential families (with some restrictions) and using the corresponding exponential family model as the emission distribution.

An exponential family distribution with (natural) parameter θ ∈ Θ ⊂ R p and sufficient statistic φ(x) ∈ R p is given by the density

p θ (x) = exp( φ(x), θ -a(θ)),
(2.10) with respect to some measure ν on the input space X (see, e.g., [START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF]; [START_REF] Barndorff-Nielsen | Information and exponential families in statistical theory[END_REF]). The log-partition (or cumulant) function a ensures normalization and is given by a(θ) = log X exp( φ(x), θ )ν(dx). The base measure ν is usually of the form ν(dx) = h(x)λ(dx), where λ is either the Lebesgue measure (for continuous distributions) or the counting measure (for discrete distributions) and h(x) is sometimes called ancillary statistic.

Most commonly used distributions, including normal, exponential, gamma, multinomial (with fixed number of trials), and many others, belong to this family. The sufficient statistic φ(x) is said to be minimal if there exists no constant c such that φ(x), θ = c ν-almost everywhere. An important property that we recall is that the derivatives of the log-partition function and the moments of the distribution are linked, and we have in particular

∇a(θ) = E X∼p θ [φ(X)].
(2.11)

In order to establish a bijection with Bregman divergence, we follow [START_REF] Banerjee | Clustering with bregman divergences[END_REF] and consider a restricted set of such distributions. Definition 2.9 (regular exponential family). A regular exponential family distribution is an exponential family distribution for which the parameter space Θ is open and the sufficient statistic x ∈ X = R p gives a minimal representation. Its density w.r.t. the Lebesgue or counting measure takes the form

p ψ,θ (x) = h(x) exp( x, θ -ψ(θ)),
(2.12) for x ∈ R p , where we denote by ψ the log-partition function.

The log-partition function is always convex, and with the additional properties of this definition, one can show that ψ and its conjugate ψ * satisfy Legendre-duality properties as in §2.2.2 (see [START_REF] Banerjee | Clustering with bregman divergences[END_REF]; [START_REF] Barndorff-Nielsen | Information and exponential families in statistical theory[END_REF] for details). If we write µ = ∇ψ(θ) -which corresponds to the mean of the distribution from (2.11) -, then θ = ∇ψ * (µ), and using (2.5), we have

p ψ,θ (x) = h(x) exp( x, θ -ψ(θ)) = h(x) exp( x, θ -µ, θ + ψ * (µ)) = h(x) exp(ψ(x) -ψ(x) + ψ * (µ) + x -µ, ∇ψ * (µ) ) = h ′ (x) exp(-D ψ * (x, µ)),
with h ′ (x) = h(x) exp(ψ(x)). Therefore, a regular exponential familiy can be expressed in terms of the Bregman divergence associated to the conjugate of the log-partition function. [START_REF] Banerjee | Clustering with bregman divergences[END_REF] show that the divergences obtained this way belong to a special family which they call regular Bregman divergences, and establish a bijection between regular exponential families and regular Bregman divergences.

In practice, this means that when D ψ is a regular Bregman divergence and we are considering distortions D ψ (•, µ) to a centroid µ, we can rely on the corresponding regular exponential family distribution p ψ * ,θ with θ = ∇ψ(µ) and with log-partition function ψ * , given by the density

p ψ * ,θ (x) = h(x) exp(-D ψ (x, µ))
(2.13) in a probabilistic context. Note that the mean of this distribution is µ, and the maximum likelihood estimate of µ given observations x 1 , . . . , x n is given by μ = n -1 i x i , which is the same as computing the right-type Bregman centroid. It is hence a natural choice of distribution to be used as an emission/observation distribution in a mixture model.

We now give two examples of one-to-one correspondences given by this bijection, following [START_REF] Banerjee | Clustering with bregman divergences[END_REF].

Example 2.10. If we consider a Gaussian distribution with mean µ and fixed covariance σ 2 I, taking the conjugate of the log-partition function gives ψ(x) = x 2 /2σ 2 and thus the corresponding Bregman divergence is D ψ (x, µ) = xµ 2 /2σ 2 , as we might have expected.

Example 2.11. In the case of a multinomial distribution with fixed number of trials N , we need to make sure to use a minimal representation in order to establish the correspondence. We can use x = [x j ] j∈{1,...,p-1} as the minimal sufficient statistic, where p j=1 x j = N , and θ = [log q j /q p ] j∈{1,...,p-1} as the natural parameter (which lies in the open set R p-1 ), where q j ≥ 0 are the event probabilities and are such that p j=1 q j = 1. The log-partition function is then ψ * (θ) = N log(1 + p-1 j=1 e θj ), and its conjugate is

ψ(µ) = N p j=1 xj N log xj N , leading to a form of KL-divergence D ψ (x, µ) = N p j=1 x j N log x j /N µ j /N .
Remark 2.12. Note that in a multinomial distribution, x is almost surely discrete since the base measure is discrete. The bijection implies that if D ψ is the KL divergence, the corresponding regular exponential family is a multinomial, hence also discrete, even though D ψ can potentially take any normalized vector x ∈ R p + . In practice, and in our experiments with the audio representation described in § 2.1, one can fix a "number of trials" integer N , normalize x ∈ R p + such that p j=1 x j = N and discretize each coordinate x j to a close integer while keeping the sum equal to N . Larger values of N can give a better, more granular approximation, but it is best to choose N empirically according to experimental results on data. Since our algorithms don't require computing the ancillary statistic h(x) thanks to normalization, the non-integer representation can also be used directly without discretization as an approximation.

Mixture models and the EM algorithm

A mixture model with K mixture components is a probabilistic latent variable model which can be described by the following generative process

z i ∼ π, i = 1, . . . , n x i |z i ∼ p µz i , i = 1, . . . , n,
where z i ∈ {1, . . . , K} is a latent random variable indicating which cluster observation x i belongs to. z ∼ π indicates that z follows a categorical distribution (or multinomial with 1 trial) with parameter π ∈ R K + , k π k = 1. p µ k denotes the observation/emission distribution of mixture

z i x i i = 1..n
µ k = (m k , Σ k ).
If we take our emission distributions to be in the regular exponential family corresponding to a Bregman divergence D ψ given by the density p (ψ * ,θ) , we can take a mean parameterization p µ k = p (ψ * ,∇ψ(µ k )) which takes the form

p µ k (x) = h(x) exp(-D ψ (x, µ k )).
(2.14)

We will call this a Bregman distribution associated with ψ.

A graphical model representation of the mixture model is given in figure 2.1. The joint distribution of the observed variables x = (x i ) i=1..n and hidden variables z = (z i ) i=1..n factorizes as follows

p(x, z; π, µ) = n i=1 p(z i ; π)p(x i |z i ; µ) = n i=1 K k=1 π ✶(zi=k) k K k=1 p(x i |z i = k; µ k ) ✶(zi=k) ,
(2.15)

where p(

x i |z i = k; µ k ) = p µ k (x i ).
In order to estimate parameters π and µ = (µ k ) k=1..K from the observed data, we can use a maximum likelihood estimator and try to maximize the (log-)likelihood of the observed data, given by

ℓ(π, µ) = log p(x; π, µ) = n i=1 log p(x i ; π, µ) = n i=1 log K k=1 p(x i , z i = k; π, µ).
(2.16)

Note that like minimizing the K-means objective, maximizing this log-likelihood is a nonconvex problem and thus makes it hard to find a global maximum. It is possible nonetheless to find local maxima using ascent methods, which try increase the value of the log-likelihood after each iteration. A standard technique used for maximum likelihood estimation in latent variable models is the Expectation-Maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF], which iteratively maximizes lower bounds on the likelihood.

Derivation of EM If we denote by x the set of observed variables, by z the hidden variables and by θ the set of parameters, we can use Jensen's inequality to find a lower bound on the log-likelihood in terms of any probability distribution q on z ℓ

(θ) = log z p(x, z; θ) = log z q(z) p(x, z; θ) q(z) ≥ z q(z) log p(x, z; θ) q(z) .
(2.17)

The EM algorithm proceeds by maximizing this lower bound with respect to q and θ, in a coordinate-ascent fasion. Maximizing it with respect to the distribution q can be done by having an equality in Jensen's inequality, thus having a tight bound. The equality case corresponds to having a constant in the expectation, i.e. q(z) ∝ p(x, z; θ), that is

q(z) = p(z|x; θ), (2.18)
which is known as the E-step, since it usually corresponds to computing expected sufficient statistics. If we now fix q in (2.17), we can now maximize this new lower bound as a function of θ, which, up to a constant term (the entropy of q), is equal to the expected complete-data likelihood under q

E z∼q [log p(x, z; θ)]. (2.19)
This is known as the M-step. This optimization problem can often be solved in closed form, as it resembles a complete-data maximum-likelihood problem. Given an initial parameter θ (0) , the EM algorithm thus repeats the following two steps until convergence

• (E-step) Compute q using current parameter θ q(z) = p(z|x; θ (t) )

• (M-step) Update parameter θ θ (t+1) = arg max θ E z∼q [log p(x, z; θ)].
Since the lower bound 2.17 increases after each step and is tight after each E-step, the loglikelihood increases after each iteration, that is, ℓ(θ (t+1) ) ≥ ℓ(θ (t) ), and thus converges, assuming it is bounded above3 . A good stopping criterion is then to see if the increase in likelihood is small. The rate of convergence of EM is linear4 , but is often preffered in practice to second-order methods because of its simplicity and because it generally suffices to quickly get a good fit to the sample data in typical machine learning applications [START_REF] Xu | On convergence properties of the EM algorithm for gaussian mixtures[END_REF].

EM for mixture models

In order to derive an EM algorithm for our mixture model, we can start by computing the expected complete-data log-likelihood used in the M-step, following Eq. (2.15)

E z∼q [log p(x, z; π, µ)] = i k E q [✶{z i = k}] log π k + i k E q [✶{z i = k}] log p(x i |z i = k; µ k ), (2.20)
where the τ ik := E q [✶{z i = k}] = p(z i = k|x i ) are computed using the previous parameter estimates in the E-step as follows

τ ik = p(z i = k|x i ; π, µ) = 1 Z p(z i = k; π)p(x i |z i = k; µ) = 1 Z π k e -D ψ (xi,µ k ) ,
where Z is a normalizing constant which ensures k τ ik = 1 and the last line holds in the specific case of Bregman emission distributions (2.14). Maximizing (2.20) with respect to π is straightforward, while maximizing with respect to µ k in the case of a Bregman distribution is equivalent to minimizing

i τ ik D ψ (x i , µ k ) i τ ik ,
which, by Proposition 2.2, gives

µ k = i τ ik x i i τ ik .
Note that this can also be seen as computing expected sufficient statistics and doing momentmatching (see, e.g., [START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF]). The entire algorithm for mixtures of Bregman distributions, i.e. Bregman soft-clustering, is given in Algorithm 2.2.

Algorithm 2.2: EM algorithm for soft clustering with Bregman divergences.

Data: x 1 , . . . , x n , initial parameters π, µ, Bregman divergence D ψ . Result: final parameters π, µ, soft-assignments τ ik repeat // E-step τ ik ← p(z i = k|x i ; π, µ) = 1 Z π k e -D ψ (xi,µ k ) // M-step π k ← 1 n i τ ik µ k ← i τ ik xi i τ ik until convergence;

Hidden Markov Models (HMMs)

Mixture models assume that the observations are i.i.d., and therefore using such a model in the case of our audio representation presented in Section 2.1 would ignore the temporal dependencies in the data. Hidden Markov models (HMMs) [START_REF] Cappé | Inference in Hidden Markov Models[END_REF][START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF]) are latent variable models in which the latent state variables evolve with Markovian dynamics, and thus are well-suited for modeling data with sequential structure, as is the case for audio.

Model

Let (x t ) t=1..T be our sequence of observations, with x t ∈ R p , and (z t ) t=1..T be our hidden state sequence, where each state is one of K states, that is z t ∈ {1, . . . , K}. An HMM is a generative latent variable model which can be described by the following generative process, for which a

z 1 z 2 z 3 . . . z T x 1 x 2 x 3 . . . x T Figure 2.2: Graphical model representation of an HMM. graphical model representation is given in Figure 2.2 z 1 ∼ π z t |z t-1 = i ∼ A i , t = 2, . . . , T x t |z t = i ∼ p µi , t = 1, . . . , T where π gives distribution of z 1 , A ∈ R K×K + is a transition matrix such that A ij = p(z t = j|z t-1 = i), A1 = 1 and A i = (A ij ) j
, where 1 = (1, . . . , 1) ⊤ . µ k is the parameter of the k-th emission distribution, which we will consider to be the distribution associated to a Bregman divergence D ψ , although one can easily extend what follows to other emission distributions such as Gaussians. The joint probability of a sequence z 1:T = (z 1 , . . . , z T ) and observations

x 1:T = (x 1 , . . . , x T ) is p(x 1:T , z 1:T ; π, A, µ) = p(z 1 ; π) T t=2 p(z t |z t-1 ; A) T t=1 p(x t |z t ; µ).
(2.21)

Inference

The goal of probabilistic inference is to infer hidden variables from observed variables, when the parameters θ of the model are fixed. This generally means computing the posterior distribution p(z Q |x; θ) over a set of query variables z Q among hidden variables. A different form of inference is maximum a posteriori (MAP) inference, which aims at finding an assignment over hidden variables with maximum posterior probability, z M AP = arg max z p(z|x; θ) = arg max z p(x, z; θ). In HMMs, there are a few inference tasks of interest:

• Smoothing: compute p(z t |x 1:T ) for t < T , uses all available observations to compute the marginal of a particular hidden state.

• Filtering: compute p(z t |x 1:t ), which is particularly useful when inference is done online, as is the case for the score-following system Antescofo.

• Prediction: compute p(z t |x 1:T ) for t > T , useful for predicting future states.

• MAP inference: compute the most likely sequence z M AP

1:T = arg max z 1:T p(z 1:T |x 1:T ).
We now describe the forward-backward and Viterbi algorithms, which are the standard algorithms for posterior and MAP inference, respectively.

Forward-Backward algorithm

The main difficulty of posterior inference comes from the complicated summations over many latent variables that appear for marginalization and for computing the normalization constants in the target posterior distribution. However, the factorized form of the joint probability in Eq. (2.21) makes it possible to propagate these sums in smaller factors and then propagate back, as is common in standard belief propation or message passing algorithms (also known as sum-product algorithms) in tree-structured graphical models, making the algorithm much more efficient. In the specific case of HMMs, it is usual to define

α t (i) = p(z t = i, x 1 , . . . , x t ) β t (i) = p(x t+1 , . . . , x T |z t = i).
The forward-backward algorithm computes these quantities recursively as follows. If we set α 1 (i) = π i p(x 1 |z 1 = i; µ i ), the other α t can be computed using the forward recursion

α t+1 (j) = p(z t+1 = j, x 1 , . . . , x t+1 ) = j p(z t = i, z t+1 = j, x 1 , . . . , x t+1 ) = j p(z t = i, x 1 , . . . , x t )p(z t+1 = j|z t = i)p(x t+1 |z t+1 = j) = i α t (i)A ij p(x t+1 |z t+1 = j; µ j ).
Similarly, if we let β T (i) = 1, the β t are computed using the backward recursion

β t (i) = j p(x t+1 , . . . , x T |z t = i, z t+1 = j)p(z t+1 = j|z t = i) = j A ij p(x t+1 |z t+1 = j; µ j )β t+1 (j).
The time complexity of the algorithm is O(T K 2 ). Note that these quantities can be very small, leading to numerical problems in the implementation. The usual way to deal with this problem is to use the logarithms of all the quantities involved, and use the log-sum-exp operation when a summation is involved.

Once the α and β quantities are computed, one can easily obtain the following useful quantities

p(z t = i|x 1:T ) = p(z t = i, x 1:t )p(x t+1:T |z t = i) p(x 1:T ) = 1 Z α t (i)β t (i) p(z t = i, z t+1 = j|x 1:T ) = 1 Z α t (i)A ij p(x t+1 |z t+1 = j; µ j )β t+1 (j) p(z t = i|x 1:t ) = 1 Z α t (i) p(x 1:T ) = i α T (i),
where Z is the appropriate normalization constant in each equation.

Viterbi algorithm

In order to compute the MAP estimate z M AP 1:T = arg max z 1:T p(z 1:T |x 1:T ), we use a similar recursion to the forward recursion, but where we replace i by max i . This is known as the max-product or Viterbi algorithm. Let's define

γ t (i) = max z1,...,zt-1 p(z 1 , . . . , z t-1 , z t = i, x 1 , . . . , x t ). If we let γ 1 (i) = π i p(x 1 |z 1 = i; µ i ), we have the recursion γ t+1 (j) = max i γ t (i)A ij p(x t+1 |z t+1 = j; µ j ).
(2.22)

By storing backpointers b t+1 (j) = arg max i γ t (i)A ij p(x t+1 |z t+1 = j; µ j ), we can then recover the MAP sequence z 1:T by taking z T = arg max i γ T (i) and z t = b t+1 (z t+1 ) for t = T -1, . . . , 1.

EM algorithm

We now derive an EM algorithm for HMMs, allowing us to estimate good parameters θ = (π, A, µ) from data. We start by computing the complete log-likelihood from Eq. (2.21)

ℓ(θ) = log p(x 1:T , z 1:T ; θ) = i ✶{z 1 = i} log π i + t≥2 i j ✶{z t-1 = i, z t = j} log A ij + t≥1 i ✶{z t = i} log p(x t |z t = i; µ i ).
The E-step thus consists in computing the expected sufficient statistics p(z t = i|x 1:T ; θ (k) ) and p(z t-1 = i, z t = j|x 1:T ; θ (k) ) given current parameter θ (k) , which is an inference task and can be solved with the forward-backward algorithm. The M-step then uses these expected sufficient statistics to maximize the expected complete log-likelihood and obtain a new parameter

θ (k+1) = arg max θ ′ E z [ℓ(θ ′ )|x 1:T ; θ (k) ].
The resulting algorithm is outlined in Algorithm 2.3. Note that because the likelihood is non-convex, initialization plays a big role, and emissions are often initialized using the results of K-means or a mixture model.

Algorithm 2.3: EM algorithm for HMMs. Data:

x 1:T , initial parameters θ = (π, A, µ), Bregman divergence D ψ . Result: final parameters θ = (π, A, µ). repeat // E-step α, β ← ForwardBackward(x 1:T , θ) τ t (i) ← p(z t = i|x 1:T ; θ) = 1 Z α t (i)β t (i) τ t (i, j) ← p(z t-1 = i, z t = j|x 1:T ; θ) = 1 Z α t-1 (i)A ij p(x t |z t = j; µ j )β t (j) // M-step π i ← τ 1 (i) A ij ← t≥2 τt(i,j) j ′ t≥2 τt(i,j ′ ) µ i ← t≥1 τt(i)xi t≥1 τt(i)
until convergence;

Hidden Semi-Markov Models (HSMMs)

One drawback of HMMs is their limited capability for expressing segment durations, that is the number of time steps during which the hidden state stays the same. In fact, if we consider the prior probability -ignoring observations -of staying in a given state i for exactly d time steps, it is equal to

A d-1 ii (1 -A ii )
, that is, the duration distribution for each state is geometric. Although the parameters A ii can be estimated to give a good fit to the data in many applications, the choice of a geometric distribution can be quite restrictive in some cases.

Hidden Semi-Markov Models (HSMMs) allow us to overcome this problem by explicitely modeling the segment durations with any distribution [START_REF] Murphy | Hidden semi-markov models (hsmms)[END_REF][START_REF] Guédon | Estimating hidden semi-markov chains from discrete sequences[END_REF].

Model

An HSMM5 can be thought of as a generalization of an HMM where each hidden state corresponds to a segment, with both a state variable z and a segment length l sampled from the

Distribution pmf parameter mean Poisson

λ k e -λ k! λ λ Negative Binomial k + r -1 k p r (1 -p) k r, p (1-p)r p
Table 2.1: Duration distributions duration distribution p z (l), and where l observations are drawn independently from the emission distribution in state z. The transition matrix A now represents transitions between segments, and we will consider it to only depend on state variables. This HMM representation isn't tractable, since it would require having KD possible hidden states, where D is the maximum segment length, and using segments as states stops us from processing observations one by one. Instead, we follow [START_REF] Murphy | Hidden semi-markov models (hsmms)[END_REF] by introducing deterministic transitions and "finish nodes": z D t will denote a counter of the number of transitions left until the end of the segment (for a segment of length d, z D t will go from d to 1, deterministically), and the finish variable f t will stay equal to 0 inside a segment, and turn to 1 when the end of the segment is reached, i.e. f t = 1 if and only if z D t = 1. We then have

p(z t = j|z t-1 = i, f t-1 = f ) = δ(i, j), if f = 0 A ij , if f = 1 (transition) p(z D t = d|z t = i, f t-1 = 1) = p i (d) p(z D t = d|z t = i, z D t-1 = d ′ ≥ 2) = δ(d, d ′ -1),
where δ is the Dirac delta (δ(i, j) = 1 if i = j and 0 otherwise). Note that having d ′ = 1 in the last identity would give the previous identity, since the events {z D t-1 = 1} and {f t-1 = 1} are the same. The rest of the model is similar to the HMM:

z 1 ∼ π x t |z t = i ∼ p µi .

Duration distributions

The distributions p i (d) can either be nonparametric tabular distributions described by each entry, or parametric distributions such as Poisson or negative binomial. Table 2.1 shows the probability mass function, parameters and means for these parametric distributions.

Complete likelihood If we consider a hidden sequence z 1:T given by a segment description (t i , l i , q i ) i=1..τ where τ is the number of segments, t i is the start of each segment, l i its length and q i its state, with

t 1 = 1, t i+1 = t i + l i , τ i=1 l i = T , the joint probability of z 1:T and x 1:T is p(x 1:T , z 1:T ; θ) = τ i=1 p qi (l i )p(q i |q i-1 ; π, A) ti+li-1 t=ti p(x t |z t = q i ; µ) , (2.23)
where θ = (π, A, µ, (p j ) j ) is the set of parameters.

Inference

Forward-backward algorithm As in the HMM case, we can derive a forward-backward algorithm using the following quantities:

α t (j) = p(z t = j, f t = 1, x 1:t ) α * t (j) = p(z t+1 = j, f t = 1, x 1:t ) β t (i) = p(x t+1:T |z t = i, f t = 1) β * t (i) = p(x t+1:T |z t+1 = i, f t = 1
). If we define α * 0 (j) = π j , the following forward recursions are obtained by marginalizing respectively on the length of the current segment and the state of the previous segment:

α t (j) = d p(x t-d+1:t |j)p j (d)α * t-d (j) α * t (j) = i α t (i)A ij ,
where we have, in our case, p(x t-d+1:t |j) = t u=t-d+1 p(x u |z u = j; µ j ). The backward recursions are similarly obtained by marginalizing respectively on the state of the next segment and the length of the current segment:

β t (i) = j β * t (j)A ij β * t (i) = d β t+d (i)p i (d)p(x t+1:t+d |i),
with β T (i) = 1. The time complexity of the algorithm is O(T K 2 D) when we force segments to be of length at most D. Note that in these recursions we assumed that the first segment starts exactly at time t = 1 and the last segment ends at time t = T . This is known as an uncensored formulation. In contrast, if we allow the last segment to end after t = T , this is known as right-censoring [START_REF] Guédon | Estimating hidden semi-markov chains from discrete sequences[END_REF], and the backward recursion on β * becomes

β * t (i) = T -t d=1 β t+d (i)p i (d)p(x t+1:t+d |i) + D i (T -t + 1)p(x t+1:T |i),
where

D i (d) := d ′ ≥d p i (d ′ ) is the survivor function of the segment length, that is D i (d) = p(l ≥ d|z = i).
The last term is thus a right-censoring term which accounts for the case of an unfinished segment at the end of the sequence. A left-censored formulation can be obtained by adding a similar term to the forward recursion, allowing initial segments to start before time t = 1.

Expected sufficient statistics

We can now compute the expected sufficient statistics for transitions across segment boundaries and durations

p(z t = i, z t+1 = j, f t = 1|x 1:T ) = 1 Z α t (i)A ij β * t (j) p(l t = d, z t = i, f t = 1|x 1:T ) = 1 Z α * t-d (i)p i (d)p(x t-d+1:t |i)β t (i)
In order to compute p(z t = i|x 1:T ), we follow [START_REF] Murphy | Hidden semi-markov models (hsmms)[END_REF] and define

γ t (i) = p(z t = i, f t = 1|x 1:T ) = 1 Z α t (i)β t (i) γ * t (i) = p(z t+1 = i, f t = 1|x 1:T ) = 1 Z α * t (i)β * t (i).
We then have

p(z t = i|x 1:T ) = τ <t (γ * τ (i) -γ τ (i)), (2.24)
which follows from the identity p(z t+1 = i|x 1:T ) = p(z t = i|x 1:T )+p(z t+1 = i, f t = 1|x 1:T )-p(z t = i, f t = 1|x 1:T ). See [START_REF] Murphy | Hidden semi-markov models (hsmms)[END_REF][START_REF] Guédon | Estimating hidden semi-markov chains from discrete sequences[END_REF][START_REF] Rabiner | A tutorial on hidden markov models and selected applications in speech recognition[END_REF] for details.

If we are in an online filtering task where the quantities p(z t = i|x 1:t ) are needed, we can use the following:

p(z t = i, x 1:t ) = d α * t-d (i)D i (d)p(x t-d+1:t |i),
which is obtained by marginalizing over the start of the segment, and where D i (d) appears because the segment is possibly unfinished [START_REF] Guédon | Estimating hidden semi-markov chains from discrete sequences[END_REF]. The p(z t = i|x 1:t ) are then obtained by normalization.

Viterbi algorithm If we define

δ t (j) = max z1:t-1 p(z 1:t-1 , z t = j, f t = 1, x 1:t ) δ * t (j) = max z1:t p(z 1:t , z t+1 = j, f t = 1, x 1:t ),
then if δ * 0 (j) = π j , we have the Viterbi recursion

δ t (j) = max d p(x t-d+1:t |j)p j (d)δ * t-d (j) b t (j) = arg max d p(x t-d+1:t |j)p j (d)δ * t-d (j) δ * t (j) = max i δ t (i)A ij b * t (j) = arg max i δ t (i)A ij ,
where b t , b * t are backpointers which can be used to recover the MAP sequence z 1:T segment by segment as follows: initialize t ← T , i ← arg max j δ t (j), then repeat the steps d ← b t (i), z t-d+1:t ← i, t ← td, i ← b * t (i) until the beginning of the sequence is reached.

EM algorithm

The E-step of the EM algorithm computes the expected sufficient statistics using the forwardbackward algorithm as described above: p(z t = i|x 1:T ) for emissions, p(z t = i, z t+1 = j, f t = 1|x 1:T ) for transitions, and p(l t = d, z t = i, f t = 1|x 1:T ) for duration distributions if they need to be learned. The M-step updates are as follows:

π i = p(z 1 = i|x 1:T ) A ij = t p(z t = i, z t+1 = j, f t = 1|x 1:T ) j ′ t p(z t = i, z t+1 = j ′ , f t = 1|x 1:T ) µ i = t p(z t = i|x 1:T )x t t p(z t = i|x 1:T )
A tabular duration distribution can be estimated easily from the sufficient statistics:

pi (d) = t p(l t = d, z t = i, f t = 1|x 1:T ) d ′ t p(l t = d ′ , z t = i, f t = 1|x 1:T )
.

However having such a nonparametric distribution can be impractical because of the large number of parameters that need to be learned. The parameter λ of a Poisson duration distribution can be estimated from the mean of the estimated distribution λi =

d dp i (d).
In the case of a negative binomial distribution with fixed r, the parameter p can be estimated as follows: pi = r r + d dp i (d) .

Summary and discussion

In this chapter, we introduced a framework for the representation of audio signals through shorttime Fourier coefficients, and their modeling with Bregman divergences, or the corresponding regular exponential family distributions in a probabilistic setting. We then developed algorithms for clustering these signals into homogeneous groups, starting with time-independent algorithms such as K-means and EM for mixture models, followed by sequential models such as HMMs and HSMMs, which explicitly model temporal dependencies and are thus more adapted to segmentation tasks.

One issue with these models is that the number of clusters K needs to be fixed in advance, which isn't always an obvious task, e.g. in an musical structure segmentation task where the number of structures present (like notes and chords) isn't known in advance. We briefly discussed a way to add new states in the K-means algorithm. In the probabilistic setting, this can be achieved by placing Bayesian nonparametric priors on the emission distributions, such as the Dirichlet process or the hierarchical Dirichlet process. These methods have been applied to HMMs [START_REF] Teh | Hierarchical dirichlet processes[END_REF][START_REF] Fox | An HDP-HMM for systems with state persistence[END_REF] and HSMMs [START_REF] Johnson | Bayesian nonparametric hidden semi-markov models[END_REF], but usually require expensive sampling procedures for Bayesian inference, for which convergence is hard to assess.

In the next chapter, we will develop online algorithms for these models, both as a way to accelerate learning and as a way to run the clustering and segmentation algorithms in a real-time context.

Chapter 3

Online algorithms

In this chapter, we will cover various algorithms for online learning in sequential models such as HMMs and HSMMs. Online learning has had many recent successes, and has been shown to often outperform batch learning, especially in large-scale settings (see, e.g., [START_REF] Bottou | The tradeoffs of large scale learning[END_REF]). In the context of latent variable models, online EM algorithms have been shown to improve performance, both in terms of speed and accuracy, compared to batch EM [START_REF] Liang | Online EM for unsupervised models[END_REF]. These algorithms have mainly been used for i.i.d. data. An extension of the online EM of [START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF] to HMMs has been given in [START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF], and we give a further extension to HSMMs in Section 3.1. We derive alternative online algorithms based on incremental minorization-maximization schemes (similar to, e.g., [START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF][START_REF]Incremental majorization-minimization optimization with application to large-scale machine learning[END_REF]), both in a non-probabilistic setting in Section 3.2 and a probabilistic one in Section 3.3, giving a new incremental EM similar to that of [START_REF] Neal | A view of the em algorithm that justifies incremental, sparse, and other variants[END_REF] for sequential models.

Online EM

In this section, we present the online EM algorithm of [START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF] for HMMs, which is based on similar ideas to the approach of [START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF] for independent observations, which we briefly present. Let's consider the case of a batch EM algorithm on a latent variable model with independent observations (x i ) i=1..n where the complete-data model p(x, z; θ) lies in the exponential family and takes the form p(x, z; θ) = h(x, z) exp( s(x, z), η(θ)a(θ)).

(3.1)

The batch EM algorithm at iteration t then takes the form

S t = 1 n n i=1 E z [s(x i , z i )|x i ; θ t-1 ] θ t = θ(S t ),
where θ(s) := arg max θ s, η(θ)a(θ) gives the complete-data maximum likelihood estimate from sufficient statistics s, which is assumed to be unique. In the case of Bregman emissions, the complete-data sufficient statistics s(x, z) are given by the vectors s 0 k (x, z) = ✶{z = k} and s 1 k (x, z) = ✶{z = k}x for every mixture component k, so that the parameter update θ is given by π k = S 0 t,k and µ k = S 1 t,k /S 0 t,k , where S 0/1 t,k are the sums in the E-step corresponding to s 0/1 k . By taking the limit n → ∞, the empirical average in the E-step converges to the expectation with respect to the true distribution P of the observations (assuming x 1 , . . . , x n iid ∼ P ), and we obtain the following limiting EM recursion

S t = E x∼P [E z [s(x, z)|x; θ t-1 ]] θ t = θ(S t ).
The online EM algorithm of [START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF] then tries to approach this limiting recursion by using a stochastic approximation (a.k.a. Robbins-Monro) procedure on the space of sufficient statistics, with an appropriate sequence of step sizes (γ t ) t , giving rise to the recursion

ŝt = (1 -γ t )ŝ t-1 + γ t E z [s(x t , z)|x t ; θt-1 ] θt = θ(ŝ t ).
We now describe the online EM algorithm of [START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF] for HMMs, which relies on a special form of smoothing called forward smoothing in order to compute sums of expected sufficient statistics in HMMs, and we will then adapt this algorithm to HSMMs by formulating them as a dynamic bayesian network with two hidden variables, which allows us to maintain a markovian assumption.

Online EM for HMMs

We consider an HMM, defined as in Section 2.4, where the final time T may be unbounded in the online case. The initial distribution π of z 0 is consider fixed as it cannot be estimated accurately from a single observation sequence (as explained in [START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF][START_REF] Cappé | Inference in Hidden Markov Models[END_REF]). We can then fully describe the model from the complete-data distribution at time t conditioned on the hidden state at time t -1, p(x t , z t |z t-1 ; θ), which we assume to be in the exponential family, just as in Eq. (3.1):

p(x t , z t |z t-1 ; θ) = h(z t , x t ) exp( s(z t-1 , z t , x t ), η(θ) -a(θ)).
(3.2)

The batch EM algorithm at iteration k then takes the form

S k = 1 T E z T t=1 s(z t-1 , z t , x t ) x 0:T ; θ k-1 θ k = θ(S k ),
where θ(s) := arg max θ s, η(θ)a(θ) is assumed unique. In the HMM with Bregman emissions, we can consider the vector s(z t-1 , z t , x t ) to be given by s ij (z t-1 , z t , x t ) = ✶{z t-1 = i, z t = j}, s 0 i (z t-1 , z t , x t ) = ✶{z t = i} and s 1 i (z t-1 , z t , x t ) = ✶{z t = i}x t for i, j ∈ {1, . . . , K} and the parameter update then becomes:

A ij = S k (i, j) j ′ S k (i, j ′ ) (3.3) µ i = S 1 k (i) S 0 k (i) , (3.4)
where S k (i, j), S 0/1 k (i) correspond to the E-step expected sums using s ij and s 0/1 i , respectively. [START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF] shows that under certain stationarity assumptions on the observed sequence and (strong) forgetting assumptions on the model, the sum in the E-step converges to an expectation, and we obtain a limiting EM recursion

S k = E x∼P [E z [s(z -1 , z 0 , x 0 )|x -∞:∞ ; θ k-1 ]] θ k = θ(S k ),
which can be approximated by the online EM algorithm we describe below, although a full proof of convergence is still lacking.

Forward smoothing recursion

The sums of expected sufficient statistics from the E-step are usually computed using a forwardbackward algorithm, which isn't suitable in an online setting. In contrast, the online EM approach of [START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF] is inspired by a forward-only smoothing recursion, which allows to recursively update the sum of interest as new observations come in. Let's define

S t = 1 t E z t t ′ =1 s(z t ′ -1 , z t ′ , x t ′ ) x 0:t ; θ φ t (i) = p(z t = i|x 0:t ) ρ t (i) = 1 t E z t t ′ =1 s(z t ′ -1 , z t ′ , x t ′ ) x 0:t , z t = i; θ .
We then have S t = i ρ t (i)φ t (i). If we initialize φ 0 (i) = 1 Z π i p(x 0 |z 0 = i) and ρ 0 (i) = 0, the forward smoothing recursion is given by the following:

φ t+1 (j) = 1 Z i φ t (i)A ij p(x t+1 |z t+1 = j) ρ t+1 (j) = 1 t + 1 E t+1 t ′ =1 s(z t ′ -1 , z t ′ , x t ′ ) x 0:t+1 , z t+1 = j; θ = 1 t + 1 i E t+1 t ′ =1 s(z t ′ -1 , z t ′ , x t ′ ) x 0:t+1 , z t+1 = j, z t = i; θ p(z t = i|z t+1 = j, x 0:t ) = i 1 t + 1 s(i, j, x t+1 ) + 1 - 1 t + 1 ρ t (i) p(z t = i|z t+1 = j, x 0:t ), (3.5)
where we used the conditional independence properties of HMMs: p(z t ′ -1 , z t ′ |z t = i, z t+1 = j, x 0:t+1 ) = p(z t ′ -1 , z t ′ |z t = i, x 0:t ) for t ′ ≤ t and p(z t = i|z t+1 = j, x 0:t+1 ) = p(z t = i|z t+1 = j, x 0:t ). This backward retrospective probability can be expressed in terms of φ as

p(z t = i|z t+1 = j, x 0:t ) = 1 Z φ t (i)A ij .

Online EM algorithm

The idea of [START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF] is to perform the M-step for updating the parameters after each iteration in the forward recursion (starting after a minimal number t min of iterations), leading to an online algorithm which isn't restricted to finite-length sequences. The recursion update (3.5) becomes a stochastic approximation update of sufficient statistics similar to the update in [START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF], and the step size 1/t can be replaced by different decreasing sequences (γ t ) t satisfying the usual requirements t γ t = ∞ and t γ 2 t < ∞. The online EM algorithm is as follows, given an initial parameter θ0 .

Initialization

φ0 (i) = 1 Z π i p(x 0 |z 0 = i; θ0 ) ρ0 (i) = 0,
The initialization of ρ should be replaced by the emission sufficient statistics of x 0 for the corresponding elements in ρ0 .

Online EM recursion For t = 1, 2, . . .

• Stochastic Approximation E-step φt+1 (j) = 1 Z i φt (i)A (t) ij p(x t+1 |z t+1 = j; θt ) ρt+1 (j) = i (γ t+1 s(i, j, x t+1 ) + (1 -γ t+1 )ρ t (i)) rt+1 (i|j), (3.6)
where rt+1 (i|j) = 1 Z φt (i)A

(t) ij approximates the backward retrospective probability.

• M-step θt+1 = θ i ρt+1 (i) φt+1 (i) , if t ≥ t min θt ,
otherwise.

Example 3.1. We saw in Eq. (3.3, 3.4) how parameters θ = (A, µ) can be updated in a batch EM algorithm (at iteration k) for HMMs with Bregman emission distributions using the quantities

S A T (i, j) = 1 T E z T t=1 ✶{z t-1 = i, z t = j} x 0:T ; Â(k) , μ(k) S µ,ν T (i) = 1 T E z T t=1 ✶{z t = i}s ν (x t ) x 0:T ; Â(k) , μ(k) ,
with s 0 (x) = 1 and s 1 (x) = x, and we will write s(x) = (s 0 (x), s 1 (x) ⊤ ) ⊤ . S A T and S µ T := (S µ,0 T , S µ,1⊤ T ) ⊤ can be obtained using the forward smoothing recursion by introducing the ρ quantities

ρ A T (i, j, k) = 1 T E z T t=1 ✶{z t-1 = i, z t = j} x 0:T , z T = k; Â(k) , μ(k) ρ µ T (i, k) = 1 T E z T t=1 ✶{z t = i}s(x t ) x 0:T , z T = k; Â(k) , μ(k) .
The sums of interest are then given by

S A T (i, j) = k ρ A T (i, j, k)φ T (k) S µ T (i, k) = k ρ µ T (i, k)φ T (k),
where φ T (k) = p(z T = k|x 0:T ). In the online EM algorithm, the ρ quantities at time t + 1 are estimated recursively from those at time t:

ρA t+1 (i, j, k ′ ) = γ t+1 k ✶{k = i, k ′ = j}r t+1 (k|k ′ ) + (1 -γ t+1 ) k ρt (i, j, k)r t+1 (k|k ′ ) = γ t+1 δ(j, k ′ )r t+1 (i|j) + (1 -γ t+1 ) k ρA t (i, j, k)r t+1 (k|k ′ ) ρµ t+1 (i, k ′ ) = γ t+1 δ(i, k ′ )s(x t+1 ) + (1 -γ t+1 ) k ρµ t (i, k)r t+1 (k|k ′ ),
where δ denotes the Kronecker delta and rt+1 the approximate backward retrospective probability. These updates require a time complexity of O(K 4 + K 3 p) for each observation. Estimates of the desired sums ŜA and Ŝµ can then be computed as before from ρ and φ, and can be used for online parameter updates (M-step)

z D t-1 z D t z D t+1 z t-1 z t z t+1 x t-1 x t x t+1
Â(t+1) ij = ŜA t (i, j) j ′ ŜA t (i, j ′ ) μ(t+1) i = Ŝµ,1 t (i) Ŝµ,0 t (i)
.

Online EM for HSMMs

Model

In order to derive an online EM algorithm for the HSMM, we introduce a new parameterization, different from the one presented in Section 2.5. We consider the HSMM as a dynamic Bayesian network with two hidden variables, z t and z D t , which correspond, respectively, to the current state and the length of the current segment up until time t (as opposed to the time left to the end of the segment in §2.5.1). That is, the event {z t = i, z D t = d} means that the chain was in state i for the past d states, and we must have z t-d+1:t = i and z D t-u = du for 0 ≤ u ≤ d -1. The transitions between segments are thus given by A ij = p(z t = j|z t-1 = i, z D t = 1). We recall the definition of the survivor function, D i (d) := d ′ ≥d p i (d ′ ), which is the probability of having a segment duration of at least d in state i.

The transition between time t -1 and time t can be encoded in the following conditional probability distributions:

p(z t = j|z t-1 = i, z D t = d) = A ij , if d = 1 δ(i, j), otherwise p(z D t = d ′ |z t-1 = i, z D t-1 = d) =      Di(d+1) Di(d) , if d ′ = d + 1 1 -Di(d+1) Di(d) , if d ′ = 1 0, otherwise.
A graphical model representation is given in Figure 3.1. The probability of having a segment in state i of length d is then:

D i (2) D i (1) D i (3) D i (2) • • • D i (d) D i (d -1) 1 - D i (d + 1) D i (d) = D i (d) -D i (d + 1) = p i (d), (3.7)
since D i (1) = 1. Note that this parameterization of the HSMM inherently includes the rightcensoring formulation of [START_REF] Guédon | Estimating hidden semi-markov chains from discrete sequences[END_REF]. Indeed, the probability that the last segment ends with

z D t = d in state i is D i (2) D i (1) D i (3) D i (2) • • • D i (d) D i (d -1) = D i (d).
The emission densities are assumed to only depend on the current state z t and are given by p(x t |z t = k). The distribution of the initial state z 0 is π and z D 0 = 1 with probability one. The complete likelihood of the model with finite horizon T is thus

p(z 0:T , z D 0:T , x 0:T ) = p(z 0 )p(z D 0 ) t≥1 p(z D t |z t-1 , z D t-1 )p(z t |z t-1 , z D t )p(x t |z t ).
The exponential family assumption for online EM is now on the model p(z t , z D t , x t |z t-1 , z D t-1 ), and we denote its sufficient statistics vector by s(z t-1 , z D t-1 , z t , z D t , x t ).

Forward smoothing recursion

Since we consider the joint state (z t , z D t ), our model satisfies the same conditional independence properties required for the forward smoothing recursion in the HMM. We have the following forward recursion:

α t (i, d) := p(z t = i, z D t = d, x 0:t ) α t+1 (j, 1) = d i α t (i, d) 1 - D i (d + 1) D i (d) A ij p(x t+1 |z t+1 = j) α t+1 (j, d) = α t (j, d -1) D j (d) D j (d -1) p(x t+1 |z t+1 = j), if d ≥ 2.
A similar recursion can be obtained for the forward filter

φ t (i, d) := p(z t = i, z D t = d|x 0:t ) = α t (i, d) j d ′ α t (j, d ′ )
by normalizing after each update, so that i d φ t (i, d) = 1. The goal is now to compute the sum

S t = 1 t E z t t ′ =1 s(z t ′ -1 , z D t ′ -1 , z t ′ , z D t ′ , x t ′ ) x 0:t ; θ (3.8) recursively. Let's define ρ t (i, d) = 1 t E z t t ′ =1 s(z t ′ -1 , z D t ′ -1 , z t ′ , z D t ′ , x t ′ ) x 0:t , z t = i, z D t = d; θ .
We then have

S t = i d ρ t (i, d)φ t (i, d).
The forward smoothing recursion is as follows.

Initialization

φ 0 (i, 1) = π i p(x 0 |z 0 = i) i π i p(x 0 |z 0 = i) φ 0 (i, d) = 0, if d ≥ 2 ρ 0 (i, d) = 0. Recursion φt+1 (j, 1) = d i φ t (i, d) 1 - D i (d + 1) D i (d) A ij p(x t+1 |z t+1 = j) φt+1 (j, d) = φ t (j, d -1) D j (d) D j (d -1) p(x t+1 |z t+1 = j), if d ≥ 2 φ t+1 (i, d) = φ(i, d) j d ′ φ(j, d ′ ) ρ t+1 (j, 1) = i d 1 t + 1 s(i, d, j, 1, x t+1 ) + 1 - 1 t + 1 ρ t (i, d) r t+1 (i, d|j, 1) ρ t+1 (j, d) = 1 t + 1 s(j, d -1, j, d, x t+1 ) + 1 - 1 t + 1 ρ t (j, d -1), if d ≥ 2,
where the backward retrospective probabilities are defined by

r t+1 (i, d|j, d ′ ) := p(z t = i, z D t = d|z t+1 = j, z t+1 = d ′ , x 0:t ).
Note that the summation in the ρ update disappears for d ≥ 2 since r t+1 (i,

d ′ |j, d) = 1 if i = j and d ′ = d -1 and r t+1 (i, d ′ |j, d) = 0 otherwise. For d = 1, we have r t+1 (i, d|j, 1) = 1 Z p(z t+1 = j, z t+1 = 1|z t = i, z D t = d)p(z t = i, z D t = d|x 0:t ) = 1 Z 1 - D i (d + 1) D i (d) A ij φ t (i, d),
with Z such that i d r t+1 (i, d|j, 1) = 1.

Online EM algorithm

Given a decreasing sequence (γ t ) t satisfying t γ t = ∞ and t γ 2 t < ∞, the minimum time before the first parameter update t min , and an initial parameter vector θ0 , the online EM algorithm for the previously described HSMM model is as follows.

Initialization φ0 (i, 1) = π i p(x 0 |z 0 = i; θ0 ) i π i p(x 0 |z 0 = i; θ0 ) φ0 (i, d) = 0, if d ≥ 2 ρ0 (i, d) = 0.
ρ0 can be initialized differently for emission sufficient statistics.

Online EM recursion For t = 1, 2, . . .

• SA E-step φt+1 (j, 1) = d i φt (i, d) 1 - D i (d + 1) D i (d) A ij p(x t+1 |z t+1 = j) φt+1 (j, d) = φt (j, d -1) D j (d) D j (d -1) p(x t+1 |z t+1 = j), if d ≥ 2 φt+1 (i, d) = φ(i, d) j d ′ φ(j, d ′ ) ρt+1 (j, 1) = i d (γ t+1 s(i, d, j, 1, x t+1 ) + (1 -γ t+1 )ρ t (i, d)) rt+1 (i, d|j, 1) ρt+1 (j, d) = γ t+1 s(j, d -1, j, d, x t+1 ) + (1 -γ t+1 )ρ t (j, d -1), if d ≥ 2,
where A, p(•|z t+1 = i) and possibly D i implicitly depend on the current parameter estimate θt , and

rt+1 (i, d|j, 1) = 1 Z 1 - D i (d + 1) D i (d) A ij φt (i, d). • M-step θt+1 = θ i d ρt+1 (i, d) φt+1 (i, d) , if t ≥ t min θt ,
otherwise.

Parameter estimation

Transition matrix The parameter update for the transition matrix A is similar to the case of regular HMMs. The expected sum of interest is

S A T (i, j) = 1 T E T t=1 ✶{z t-1 = i, z t = j, z D t = 1} x 0:T ,
and we define the corresponding forward smoothing ρ quantity

ρ A T (i, j, k, u) = 1 T E T t=1 ✶{z t-1 = i, z t = j, z D t = 1} x 0:T , z T = k, z D T = u .
The online EM recursion for estimating ρ A is as follows:

ρA t+1 (i, j, k ′ , 1) = γ t+1 δ(j, k ′ )r t+1 (i|j, 1) + (1 -γ t+1 ) k u ρA t (i, j, k, u)r t+1 (k, u|k ′ , 1) ρA t+1 (i, j, k ′ , u) = (1 -γ t+1 )ρ A t (i, j, k ′ , u -1), if u ≥ 2, where rt+1 (•, •|•, •)
is the approximate backward retrospective probability, and rt+1 (i|j, 1) = d rt+1 (i, d|j, 1). We then have an online EM approximation of the statistics

S A t : ŜA t (i, j) = k u ρA t (i, j, k, u)φ t (k, u),
and the parameter update on the transition matrix becomes

Â(t+1) ij = ŜA t (i, j) j ′ ŜA t (i, j ′ )
.

The updates for emission parameters are the same as the HMM case (with the additional duration state variable in the online EM recursion). The time complexity of the A and µ updates is now O(K 4 D + K 3 Dp), which only grows linearly with the maximum duration D thanks to the deterministic transitions in our model.

Duration distributions

In order to estimate duration distributions, we can consider the quantities λ i (d) = D i (d + 1)/D i (d) as our duration parameters. These parameters can then be estimated using the following expected sums of sufficient statistics:

S T (i, d) = 1 T E T t=1 ✶{z t-1 = i, z D t-1 = d, z D t = d + 1} x 0:T S ′ T (i, d) = 1 T E T t=1 ✶{z t-1 = i, z D t-1 = d, z D t = 1} x 0:T ,
or their online EM estimates. The maximum likelihood parameter is then given by λi

(d) = S T (i, d) S T (i, d) + S ′ T (i, d) .
The survivor function D i (d) can then be obtained by setting D i (1) = 1 and

D i (d + 1) = λ i (d)D i (d) for d ≥ 1.
One can also estimate a parametric duration distribution from the obtained tabular distribution. Note that it is hard to estimate duration distributions in an online manner at the beginning of the sequence since only few segments can be observed, thus it can help to add priors, either on the λ i or on the parameters of the parametric duration distribution, and perform MAP estimation rather than maximum likelihood (see §3.4).

Non-probabilistic models

If we consider a mixture model with emission distributions associated to the Bregman divergence D ψ , and with π k = 1/K fixed, then we have:

-log p(x, z; µ) = - i log p(z i ) - i log p(x i |z i ; µ zi ) = C + i D ψ (x i , µ zi ).
Thus, the negative comoplete-data log-likelihood is equal to the K-means objective (2.9), up to a constant term. The K-means algorithm can thus be seen as alternatively maximizing the complete-data log-likelihood with respect to the hidden variables z and the parameters µ. Or, conversely, if we have a probabilistic latent variable model, we can consider the complete-data likelihood as an objective function, and optimize it directly with respect to the latent variables and the parameters. More generally, we are replacing the the sum-product operation used to obtain expected sufficient statistics in the E-step of EM with a max-product operation, which finds a set of variables z with maximum joint probability (which is equivalent to MAP estimation). The two steps in the main loop of the algorithm are therefore:

z (k) = arg max z p(x, z; θ (k-1) ) θ (k) = arg max θ p(x, z (k) ; θ).
In the case of HMMs, the first step is given by the Viterbi algorithm, and the second steps computes the centroids of the points in each cluster. This has been called the segmental K-means algorithm [START_REF] Juang | The segmental k-means algorithm for estimating parameters of hidden markov models[END_REF], and a similar algorithm can be derived for HSMMs.

Online algorithm

If we consider an HMM with parameters π and A fixed and look at the form of the joint likelihood p(x 1:T , z 1:T ; µ) given in Eq. ( 2.21), we can derive a corresponding objective function to be minimized:

ℓ(z 1:T , µ) = 1 T t≥1 D ψ (x t , µ zt ) + λ 1 T t≥2 d(z t-1 , z t ), (3.9)
where d(z t-1 , z t ) is the cost of a transition from z t-1 to z t and λ 1 a parameter that we fix. An offline algorithm could minimize this objective with respect to the sequence z 1:T using a Viterbi algorithm, then with respect to µ by computing centroids. This can be seen as a majorizationminimization algorithm [START_REF]Incremental majorization-minimization optimization with application to large-scale machine learning[END_REF] on the objective

f T (µ) := min z 1:T ℓ(z 1:T , µ), (3.10)
where the majorizing surrogates (upper bounds on f T ) take the form µ → ℓ(z 1:T , µ) for a fixed sequence z 1:T , and are tight at the current parameter µ after the "E-step" (Viterbi step). In an online setting, we can follow [START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF] and at each time step t, minimize a new surrogate

ft (µ) := 1 t t i=1 D ψ (x i , µ zi ) + λ 1 t t i=2 d(z i-1 , z i ), (3.11)
where z 1:t-1 have been previously computed and z t is chosen greedily to minimize this surrogate with the current parameter µ. µ is then updated by minimizing ft , which corresponds to computing new centroids and can be done as an online update rule similar to the online K-means algorithm described in [START_REF] Bottou | Convergence properties of the k-means algorithms[END_REF]; [START_REF] Bottou | Online learning and stochastic approximations[END_REF]. The function ft is still a majorizing surrogate, but given that the previous sequence elements are not updated, the bound on f t isn't tight, as in [START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF]. The algorithm is given in Algorithm 3.1.

Algorithm 3.1: Online optimization algorithm for the non-probabilistic HMM.

Data: Input stream x t , t = 1, . . ., initial centroids µ 1 , . . . , µ K , Bregman divergence D ψ . Result: final centroids µ 1 , . . . , µ K n 1 = . . . = n K = 0; // cluster counts z 1 ← arg max j D ψ (x 1 , µ j ); µ z1 ← x 1 ; for t = 2, . . . do j ← arg min k D ψ (x t , µ k ) + λ 1 d(z t-1 , k) z t ← j; n j ← n j + 1; µ j ← µ j + 1 nj (x t -µ j ) end
Adding new states The algorithm can easily be modified to add new states when needed, by adding a penalty term on the number of clusters K to the objective, as we did with K-means:

ℓ(z 1:T , µ, K) = 1 T t≥1 D ψ (x t , µ zt ) + λ 1 T t≥2 d(z t-1 , z t ) + λ 2 T K. (3.12)
This gives us Algorithm 3.2, where we consider the transition cost to be d(j, k) = 0 if j = k and 1 otherwise.

Semi-Markov extension

Note that the algorithm can be extended to a semi-Markov model, by considering the parameterization described in §3.1.2, and taking a transition cost of the form

d(z t-1 , z D t-1 , z t , z D t ) = -log p(z t , z D t |z t-1 , z D t-1 ).
all i since the maximizer takes this form. The lower bound then takes the form fq (θ) = z1,...,zn

( i q i (z i )) log i p(x i , z i ; θ) i q i (z i ) = n i=1 z1,...,zn ( i q i (z i )) log p(x i , z i ; θ) q i (z i ) = n i=1 zi q i (z i ) log p(x i , z i ; θ) q i (z i ) .
The incremental EM algorithm proceeds by selecting a data point i, maximizing the lower bound with respect to q i , i.e. taking q i (z i ) = p(z i |x i ; θ), then maximizing the updated minorizing surrogate fq . In the online streaming case (where new observations come in), one can consider surrogates of the form

fn (θ) = 1 n n i=1 zi q i (z i ) log p(x i , z i ; θ) q i (z i ) , ( 3.13)
where q 1 , . . . , q n-1 are fixed from the past. q n (z n ) = p(z n |x n ; θ) is computed in the E-step using the current parameter estimate θ, and the M-step then maximizes this new surrogate fn , which is equal to (1/n) E q [log p(x, z; θ)] up to a constant entropy term on q independent of θ. When the complete-data model lies in the exponential family as in Eq. (3.1), we have

fn (θ) = C + 1 n n i=1 zi q i (z i )( s(x i , z i ), η(θ) -a(θ)) - 1 n n i=1 zi q i (z i ) log q i (z i ) = C + 1 n n i=1 E zi∼qi [s(x i , z i )], η(θ) -a(θ) + 1 n H(q),
where H(q) is the entropy of q and C a constant which only depends on observations. This leads to an online algorithm which is similar to the online EM of [START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF] described in Section 3.1, where sufficient statistics are updated with step sizes of the form γ n = 1/n and the parameter updates are obtained in closed form from sufficient statistics.

Incremental EM for HMMs

In the case of HMMs, the factorization q(z 1:T ) = t q t (z t ) doesn't hold anymore since observations aren't i.i.d., but using the chain rule and conditional independences in HMMs, we can factor the posterior distribution as follows:

p(z 1:T |x 1:T ; θ) = p(z 1 |x 1:T ) t≥1 p(z t |z t-1 , x 1:T ).
Hence, we can limit ourselves to distributions q of the form q(z 1:T ) = q 1 (z 1 ) t≥2 q t (z t |z t-1 ) (with j q t (j|i) = 1 for all i, t ≥ 2 and q 1 fixed to π) in the lower bound on the likelihood, since the maximizer takes this form. We then have (with the abuse of notation q 1 (z

1 |z 0 ) = q 1 (z 1 )) f T (θ) = 1 T log p(x 1:T ; θ) ≥ 1 T z 1:T q(z 1:T ) log p(x 1:T , z 1:T ; θ) q(z 1:T ) = 1 T z 1:T q(z 1:T ) log t p(x t , z t |z t-1 ; θ) t q t (z t |z t-1 ) = 1 T T t=1 z 1:T q(z 1:T ) log p(x t , z t |z t-1 ; θ) q t (z t |z t-1 ) = 1 T T t=1   zt-1,zt φ t-1 (z t-1 )q t (z t |z t-1 ) log p(x t , z t |z t-1 ; θ) q t (z t |z t-1 )   =: fT (θ),
where we defined φ t (z t ) := z1:t-1 q 1 (z 1 )q 2 (z 2 |z 1 ) . . . q t (z t |z t-1 ) = zt-1 φ t-1 (z t-1 )q(z t |z t-1 ) for t ≥ 2 and φ 1 (z 1 ) = q 1 (z 1 ). The quantity φ t (z t ) corresponds to the marginal of q on z t , while φ t-1 (z t-1 )q t (z t |z t-1 ) is the pairwise marginal on z t-1 and z t , which appears naturally in our lower bound. Thus, in the online streaming setting, we can consider at time T the surrogate fT , where q t and φ t are fixed from the past for t < T , and q T is obtained in the E-step to maximize the surrogate at the current parameter estimate, giving q T (z

T |z T -1 ) = p(z T |z T -1 , x T ; θ).
The M-step then maximizes this new surrogate. If we assume p(x t , z t |z t-1 ; θ) is in the exponential family as in Eq. 3.2, the surrogate takes the form

fT (θ) = C + 1 T T t=1   zt-1,zt φ t-1 (z t-1 )q t (z t |z t-1 )( s(z t-1 , z t , x t ), η(θ) -a(θ))   + 1 T H(q) = C + S T , η(θ) -a(θ) + 1 T H(q),
where

S T = 1 T T t=1
zt-1,zt φ t-1 (z t-1 )q t (z t |z t-1 )s(z t-1 , z t , x t ) and C is independent of θ. If we define θ(s) := arg max θ s, η(θ)a(θ) the maximizer of the complete-data likelihood (assumed unique as in §3.1.1), the surrogate fT is maximized by taking θ = θ(S T ). We obtain Algorithm 3.3. In order to get our definition of S T , the step sizes are taken to be γ t = 1/t, although in practice we often obtained improvements for slower step sizes such as γ t = t -0.6 . Algorithm 3.3: Incremental EM algorithm for HMMs.

Data: Input stream x t , t = 1, . . ., initial parameter θ, fixed π, first M-step time t min .

Result:

final parameter θ φ 1 (i) ← p(z 1 = i|x 1 ; θ) = 1 Z π i p(x 1 |z 1 = i; θ); S ← i φ 1 (i)s(i, x 1 ); for t = 2, . . . do // E-step q t (j|i) ← p(xt,zt=j|zt-1=i;θ) j ′ p(xt,zt=j ′ |zt-1=i;θ) ; φ t (j) ← i φ t-1 (i)q t (j|i); S ← (1 -γ t )S + γ t i j φ t-1 (i)q t (j|i)s(i, j, x t ); // M-step if t ≥ t min then θ ← θ(S); end end
Example 3.2. If we now consider the HMM with transition matrix A and emissions associated to Bregman divergence D ψ with centroids µ k , we have p(x t , z t = j|z t-1 = i; θ) = A ij p(x t |z t = j; µ j ).

The sufficient statistics s(z

t-1 , z t , x t ) are given as in §3.1.1 by s ij (z t-1 , z t , x t ) = ✶{z t-1 = i, z t = j}, s 0 i (z t-1 , z t , x t ) = ✶{z t = i}, s 1 i (z t-1 , z t , x t ) = ✶{z t = i}x t .
It follows that the updates of sufficient statistics in the E-step are simply

S A t (i, j) = (1 -γ t )S A t-1 (i, j) + γ t φ t-1 (i)q t (j|i) S µ,0 t (i) = (1 -γ t )S µ,0 t-1 (i) + γ t φ t (i) S µ,1 t (i) = (1 -γ t )S µ,2 t-1 (i) + γ t φ t (i)x t .
Note that these updates do not require summing over the previous state as in the ρ updates of the forward smoothing recursion in §3.1.1, and are thus less costly, with a time complexity of O(K 2 + Kp). The M-step is similar to the online EM algorithm in §3.1.1:

A (t) ij = S A t (i, j) j ′ S A t (i, j ′ ) µ (t) i = S µ,1 t (i) S µ,0 t (i)
.

Semi-Markov extension

Our incremental EM algorithm can be extended to HSMMs using the two-variable hidden chain described in Section 3.1.2. The quantities of interest are now q t (z t , z D t |z t-1 , z D t-1 ) and φ t (z t , z D t ). Thanks to the deterministic transitions, the quantities q t (j, d ′ |i, d) will be zero both if d ′ / ∈ {1, d + 1} or if i = j and d = d + 1, thus the E-step reduces to the following:

qt (j, 1|i, d) = (1 -λ i (d))A ij p(x t |z t = j) qt (i, d + 1|i, d) = λ i (d)p(x t |z t = i) q t (j, d ′ |i, d) = 1 Z qt (j, d ′ |i, d), if d ′ = 1 or d ′ = d + 1 and i = j 0, otherwise. φ t (j, 1) = i,d φ t-1 (i, d)q t (j, 1|i, d) φ t (j, d) = φ t-1 (j, d -1)q t (j, d|j, d -1), if d ≥ 2,
where λ i (d) = D i (d + 1)/D i (d) and the normalizing constant Z is given by Z = qt (i, d + 1|i, d) + j qt (j, 1|i, d). The complexity of these updates is thus O(K 2 D). Example 3.3. If we consider the HSMM version of the model in Example 3.2, the relevant sufficient statistics for estimating the transition matrix and the emission parameters are ✶{z t-1 = i, z t = j, z D t = 1}, ✶{z t = i} and ✶{z t = i}x t . The updates are then the following:

S A t (i, j) = (1 -γ t )S A t-1 (i, j) + γ t d φ t-1 (i, d)q t (j, 1|i, d) S µ,0 t (i) = (1 -γ t )S µ,0 t-1 (i) + γ t d φ t (i, d) S µ,1 t (i) = (1 -γ t )S µ,2 t-1 (i) + γ t d φ t (i, d) x t .
The M-step update are similar to the HMM (see Example 3.2). This gives us an overall time complexity per observation of O(K 2 D + KDp), where D is the maximal duration of a segment.

Note that one can also estimate the duration parameters λ i using the statistics ✶{z t-

1 = i, z D t-1 = d, z D t = d + 1} and ✶{z t-1 = i, z D t-1 = d, z D t = 1}
, as in §3.1.2.

Including prior knowledge with Bayesian priors

The online algorithms we presented so far can be noisy in the beginning and take some time before giving good results. In addition, one often needs to wait a given number t min of iterations before carrying M-step updates, so that the sufficient statistics are somewhat stabilized and the parameter update is numerically stable.

Often, we might have some reliable prior information about what the different clusters should look like. For example, if we are segmenting musical structures, such as notes and chords, we might know in advance the harmonic structure of these elements. In this case, it is desirable to encode this information into the model in order to improve the accuracy. A simple way to take this into account is to have good initial parameters which represent each cluster, but this initialization will be ignored after the first parameter update. A more suitable and natural way to encode prior information is to add custom Bayesian priors on the model parameters, and replace maximum likelihood updates of the parameters with maximum a posteriori (MAP) updates.

We will limit ourselves to conjugate priors, which will lead to a very similar form of the updates, by simply having an additional number of "virtual" observations that trade off against actual observations, and are forgotten as the number of observations increases. If we consider an exponential family form for p(x t , z t |z t-1 , θ):

p(x t , z t |z t-1 , θ) = h(z t , x t ) exp( s(z t-1 , z t , x t ), η(θ) -a(θ)),
the corresponding conjugate prior on θ takes the form

p(θ; κ, τ ) = exp( τ κ, η(θ) -τ a(θ) -b(τ, κ)), (3.14)
where b is the log-partition function. The product τ κ is often considered itself as a parameter, but our parameterization leads to a natural interpretation, where τ corresponds to an equivalent sample size, and κ a vector of sufficient statistics that represents what we consider to be a good observed sufficient statistic s(z t-1 , z t , x t ). To see why this is true, let's consider a completely observed HMM with state sequence z 1:T and observations x 1:T . Maximizing the posterior p(θ|x 1:T , z 1:T ) is equivalent to maximizing the product of the likelihood and the prior since by Bayes rule we have p(θ|x 1:T , z 1:T ) ∝ p(x 1:T , z 1:T |θ)p(θ). Taking the logarithm of this quantity gives us log p(θ; κ, τ )p(x 1:T , z 1:

T |θ) = C + τ κ, η(θ) -τ a(θ) + T t=1 s(z t-1 , z t , x t ), η(θ) -T a(θ) = C + τ κ + T S T , η(θ) -(τ + T )a(θ),
where S T = 1 T T t=1 s(z t-1 , z t , x t ) and C is a constant which doesn't depend on θ. If θ(s) := arg max θ s, η(θ)a(θ), the maximum a posteriori estimate of θ is then given by θMAP = θ τ κ + T S T τ + T .

(3.15)

We can see the desired trade-off between the statistics of the prior κ and of the observations S T , which leans towards the prior when T is small and towards observed data when T becomes large.

When the sequence z 1:T is unobserved, the same parameter update is obtained in the M-step of a (batch or online) EM algorithm, by considering S T to be the expected average sum of sufficient statistics, as we defined for the online EM algorithm of Section 3.1.1 and the incremental EM algorithm of Section 3.3.1.

Example 3.4. In our HMM with Bregman emissions, using the notations of Example 3.2, we can define priors on A using an equivalent number of transitions κ A ij that trade-off against the observed transitions, and for µ, an equivalent number of observations κ µ,0 i and a vector κ µ,1 i ∈ R p , which might correspond to an ideal short-time Fourier spectrum vector for a given note or chord in the context of musical note segmentation. For simplicity, we can take τ = 1 in this case, since the vector κ already encodes equivalent sample sizes. The M-step updates are then

A (t) ij = κ A ij + T S A t (i, j) j ′ (κ A ij ′ + T S A t (i, j ′ )) µ (t) i = κ µ,0 i κ µ,1 i + T S µ,1 t (i) κ µ,0 i + T S µ,0 t (i)
.

Note that the transitions are given by categorical distributions (multinomial with a single trial), hence the corresponding conjugate priors are Dirichlet distributions. Similarly, if the emissions are multinomial, that is, D ψ is the KL divergence, the conjugate prior is also Dirichlet.

In the case of HSMMs, similar priors can be defined, and additional priors can be defined on duration distributions, the conjugate priors being Gamma distributions for Poisson durations, and Beta distributions for the parameter p of Negative Binomial durations.

Experiments on synthetic data

In this section, we compare the behavior of online EM ( §3.1.1) and incremental EM ( §3.3.1) for HMMs on synthetic sequences generated from some fixed parameters, with either Gaussian or KL (multinomial) emissions. See [START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF] for a description of the sufficient statistics updates used for Gaussians. We used a training sequence of length 4000 and a test sequence of length 1000 to monitor how the test likelihood evolves. The step sizes were fixed to γ t = t -0.6 , and the time of the first M-step update to t min = 80.

Gaussian toy example

We first compared the two algorithms on two-dimensional Gaussians with K = 2:

µ 1 = 3 0 Σ 1 = 2 1 1 4 µ 2 = -2 3 Σ 2 = 3 -1 -1 2 .
Initial parameters were set to µ 0 1 = (1, 0) ⊤ , µ 0 2 = (-1, 0) ⊤ and covariances equal to identity. Figure 3.2 shows the evolution of parameter estimates and log-likelihood during the algorithm. We can see that both algorithms give very similar parameter estimates in this low-dimensional case, and that the test log-likelihood comes close to the batch value after a just a few hundred observations.

We noticed that the online approaches are not effective when considering Gaussians with full covariance matrices in higher dimensions (see Figure 3.3). This could be due to the fact that many observations are needed to correctly estimate covariance matrices, and it might be a bad idea to update them in an online manner without a strong prior.

Squared euclidian distance and KL divergence

We then compared batch and online algorithms for Bregman divergence-based emissions, in particular the squared Euclidian distance (corresponding to Gaussians with covariance matrix equal to identity) and the KL divergence (multinomial, for which we fixed the number of trials to N = 100), for different values of the number of clusters K and the dimensionality of the data p. Our Gaussian centroids were sampled according to N (0, I), and the means of our multinomials were sampled from a Dir(1, . . . , 1). Initial parameters were sampled from N (0, I) for Gaussians, and for the multinomial, they were taken to be slightly perturbed uniform distributions. 3.4, 3.5, 3.6 show results of the different algorithms for different values of K and p. The plots on the top show the evolution of the log-likelihood, both for the batch algorithm across each EM iteration, and for the online algorithms over time (plotted against the final value of batch EM). The bottom plots show the training sequence of observations, projected onto their first two principal components, where the color of each point gives the true or predicted state (for the batch algorithm, we show smoothing estimates, i.e. maximum posterior marginals, while for online algorithms we show the predicted state given by filtering during the algorithm).

We can see that the online algorithms give good predictions, and give a good value for the test log-likelihood quite quickly, especially when the dimensionality is relatively low. The convergence is usually faster in the KL divergence case, and the incremental EM approach generally performs better than the online EM algorithm of [START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF] when looking at the test likelihood, despite its lower computational complexity. For instance, in Figure 3.5b, we can see that after a few hundred observations incremental EM performs just as well (in terms of test likelihood) as batch EM after convergence, which required 5 iterations on the entire training sequence to converge. Of course, the non-convexity of the problem makes it highly sensitive to the random initializations, but these remarks remain consistent with the multiple simulations we ran.

Summary and discussion

In this chapter, we covered various algorithms for online learning in HMMs and HSMMs, mainly based on the EM algorithm. We explored algorithms based on stochastic approximations, following [START_REF] Cappé | Online expectation-maximization algorithm for latent data models[END_REF][START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF], as well as incremental majorization-minimization (or minorization-maximization) algorithms, with both probabilistic and non-probabilistic versions. We discussed how to add new states on the fly in some models, and how to encode prior information about the data in the model through the use of Bayesian conjugate priors. Some questions that can be explored in future work include the use of mini-batches in online algorithms for HMMs (see the use of blocks of increasing size in the Block Online EM algorithm in [START_REF] Corff | Online expectation maximization based algorithms for inference in hidden markov models[END_REF], or the frequent use of mini-batches in i.i.d. data, e.g., [START_REF] Liang | Online EM for unsupervised models[END_REF][START_REF] Mairal | Online learning for matrix factorization and sparse coding[END_REF], and an analysis of convergence of our incremental algorithms. Chapter 4

Audio segmentation experiments

In this chapter, we describe and discuss some of our experiments on the audio segmentation task using the algorithms presented in the previous chapters, both in the context of musical note/structure segmentation, and acoustic scene segmentation. For musical structure segmentation, we will use two main examples, one by M. Ravel for piano and one by J.S. Bach for violin. Their music score is shown in Figure 4.1. We use the short-time spectral representation described in Section 2.1, computed using Hamming windows of size 4096 with an offset of size 512 between each window (the sample rate of our audio files is 44.1KHz). Our emission distributions are multinomial distributions, that is, the regular exponential family associated to the KL divergence. The normalization constant of our short-time vectors (which corresponds to a number of trials N in the multinomial) is chosen depending on the example: we used N = 5 for the Ravel example and N = 20 for the Bach sonata. For the acoustic scenes, we used sounds from the Office Live Dataset [START_REF] Giannoulis | Detection and classification of acoustic scenes and events: An ieee aasp challenge[END_REF], with a similar representation, and we fixed N = 150.

Offline segmentation results

When running the segmentation algorithms offline, we have the entire audio signal at hand, and it is common in this case to initialize the emission parameters of the HMM and HSMM with centroids obtained using the K-means algorithm (itself run multiple times with different random initializations). Figure 4.2 shows results on the Ravel excerpt, with two different Kmeans initializations. In the case of HSMMs, we fixed the duration distributions to a Negative Binomial distribution with parameters r = 5, p = 0.95, up to a maximal duration D = 200. We can see that the initialization has a strong effect on the resulting sequences: the results on the left are very accurate for the HSMM, while the ones on the right show some segments which include multiple different notes. We can also see that the explicit duration distributions given by the HSMM help the model avoid short segments, such as attacks, which are obtained for all other models. Note that the N B(5, 0.95) distribution is actually quite similar to a geometric distribution, but we are forcing a slower decay, while the HMM learns a very steep duration distribution for the attack states through the transition matrix. Figure 4.3 shows results on the Bach example. In this case, the HMM and HSMM emissions were initialized with slightly perturbed uniform multinomial distributions (µ j ∝ 1 + ǫ j for j = 1, . . . , p, where ǫ j is a small noise term). We can see that non-sequential methods (K-means is initialized multiple times and EM uses K-means as initialization) perform well here, since the notes are quite distinctive in this example, and the onsets are not as disruptive on the violin as they are on the piano. We can see however that the orange segment in K-means and EM is split into two by the HSMM, thanks to the duration distribution which encourages segments of length close to 20. Note that despite the random initialization, the HSMM still manages to obtain accurate segments.

q = 70 ! ! " " # # #! # # ! $ # # #! # # ! $ # # # # # $ # # # # # $ A B C D E F G ∅ ∅ A B C D E F G ∅ (

Online EM for HMMs and HSMMs

Figure 4.4 compares the results of the online EM algorithms for HMMs and HSMMs described in Section 3.1 on the example by Bach. The duration distributions were chosen to be Negative Binomial with mean 20 and parameter r = 30. The first M-step time was t min = 80. We show both the filtered state given in real-time during online estimation, and the full viterbi sequence obtained with the final parameters, and in Figure 4.4b we repeat the sequence of notes twice to see how the results improve. We can see that the resulting segmentation is indeed improved after seeing the sequence twice, and we can see that the HSMM manages to discover more notes (clusters) than the HMM. The execution time can be quite slow due to the expensive forward smoothing recursions, especially in the HSMM case. In the HMM case, the segmentation took about 3s on a laptop (which is close to real-time) for the single sequence, while in the HSMM case it took about 40s, with a truncation of the duration distribution to D = 70, which is definitely too long for real-time applications, although it could be improved by using a compiled language rather than Python. 3.1, respectively. We show results for a single sequence or the same sequence repeated twice, with K = 10. We can see that both methods give similar results. The online EM approach seems to model transitions better than incremental EM, which gives many short segments, especially in the shorter sequence, while incremental EM seems to discover more states compared to online EM. Another important difference is in the execution time: incremental EM is faster by almost an order of magnitude compared to online EM (0.4s vs 1.5s), and the difference gets bigger as K increases, since the complexities per observation differ by a factor O(K 2 ).

Online vs incremental EM for HMMs

Segmentation of acoustic scenes

The previous examples we've shown considered the task of musical note segmentation, in which we wish to obtain a single segment for each note. In this situation, a single note or chord is mostly homogeneous and usually corresponds to a single segment (although the onset can be quite different). If we instead consider the task of detecting acoustic scenes or events -such as a door slam, a cough or the sound of keys being dropped -, the audio content isn't always homogeneous anymore, but might be described by a sequence of segments rather than a single segment. If we manage to get a similar sequence of segments for different instances of the acoustic scene, it is then easy to detect the scene with a higher-level algorithm on top of this segmentation. In Figure 4.6, we consider two audio signals, each of which contains three different instances of two acoustic scenes, taken from the Office Live Dataset [START_REF] Giannoulis | Detection and classification of acoustic scenes and events: An ieee aasp challenge[END_REF]. Emissions are initialized to slightly perturbed uniform distributions on spectra, and HSMM durations are fixed to N B(5, 0.2) (mean 20). For the incremental EM, we also show the real-time estimate ("filter"), in addition to the viterbi sequence given by the final parameters.

The top example alternates between a telephone ringing sound and a coughing sound. As we can see in the spectrogram, the telephone ringing instances look quite similar, while the middle coughing is quite different from the other two. We can see that all three algorithms give consistent segment sequences across the different instances. By encouraging longer segments, the HSMM has the added benefit of having fewer small segments for each scene, for instance the coughing sound is made of a single cluster (light blue), while the HMM breaks it into several clusters, which are different in the middle coughing. Finally, the segmentation given by the incremental EM algorithm for HMMs is quite satisfactory, given that it only observed the sequence once. Note that the background sound (silence) was detected accurately here in a single cluster, but in general one could also use a separate algorithm -based, e.g., on spectral flatness -for distinguishing it from actual content before performing the segmentation.

The bottom example alternates between sounds of keys dropping on a table and door slams. The second key drop is preceded by a sound of shaking keys that is clearly visible in the spectrogram. The background here is quite noisy and thus gives us several different clusters, but these are different from the ones obtained during the acoustic scenes, and it can be helpful in this case to rely on a separate silence detection algorithm for ignoring the background sound. If we focus on the final viterbi sequences around the acoustic events, we can see that all three algorithms give consistent segmentations across the different instances, and even manage to capture the shaking key sound (between t = 210 and t = 280) as a separate cluster. Once again, we can see that the HSMM uses a single segment to describe some scenes while the HMM might split it in two: for instance, the door slam sound uses a single cluster in the HSMM (dark blue), versus two clusters in the HMM (brown and red in the offline HMM).

These segmentation results on acoustic scenes are quite promising, and could eventually be used in systems such as Antescofo for detecting other sounds than musical notes, including some percussion instruments, or other complex acoustic events that are increasingly common in contemporary classical music, and which cannot be described by hand-crafted emission templates, as is currently done for musical notes in Antescofo.

Summary and discussion

In this chapter, we applied the learning algorithms presented in previous chapters to the task of audio segmentation and clustering. We've seen that offline algorithms can give very accurate segmentations of musical notes, and that the explicit modeling of durations in HSMMs can help in avoiding short segments or encouraging segments of specific lengths. We then experimented with online algorithms and their use for real-time audio segmentation. Although results in general are not as good as in the offline case, we saw that these algorithms can still give good segmentation results despite the initial lack of data, and how they can improve when more data is observed over time. Note that if prior information about the different clusters is at hand, encoding it into the model -as explained in Section 3.4 -can help getting better segmentation results from the start, for instance by having hand-designed prior templates for each different musical note. We also explored the use of our segmentation algorithms for acoustic scenes, and saw how different examples of a similar acoustic event produce similar sequences of segments, which can then be used at a higher level for detecting and recognizing these events. This could be a promising start towards the goal of having a real-time artificial listening machine, which might hear a few examples of a sound sequence and then recognize it later when it hears it again.

Figure 2 . 1 :

 21 Figure 2.1: Graphical model representation of a mixture model.

Figure 3 . 1 :

 31 Figure 3.1: Graphical model representation of the HSMM model used for online EM.

Figure 3

 3 Figure 3.2: Gaussian HMM toy example, online EM vs incremental EM. (a) Evolution of parameters µ 1,1 and µ 1,2 for both algorithms. (b) Test log-likelihood across EM iterations for the batch EM algorithm and online observations for the online algorithms (the x axis denotes every 100th iteration after the first M-step update).

Figures

  Figures 3.4, 3.5, 3.6 show results of the different algorithms for different values of K and p. The plots on the top show the evolution of the log-likelihood, both for the batch algorithm across each EM iteration, and for the online algorithms over time (plotted against the final value of batch EM). The bottom plots show the training sequence of observations, projected onto their first two principal components, where the color of each point gives the true or predicted state (for the batch algorithm, we show smoothing estimates, i.e. maximum posterior marginals, while for online algorithms we show the predicted state given by filtering during the algorithm).We can see that the online algorithms give good predictions, and give a good value for the test log-likelihood quite quickly, especially when the dimensionality is relatively low. The convergence is usually faster in the KL divergence case, and the incremental EM approach generally performs better than the online EM algorithm of[START_REF] Cappé | Online EM algorithm for hidden markov models[END_REF] when looking at the test likelihood, despite its lower computational complexity. For instance, in Figure3.5b, we can see that after a few hundred observations incremental EM performs just as well (in terms of test likelihood) as batch EM after convergence, which required 5 iterations on the entire training sequence to converge. Of course, the non-convexity of the problem makes it highly sensitive to the random initializations, but these remarks remain consistent with the multiple simulations we ran.

Figure 3 Figure 3 Figure 3 Figure 3

 3333 Figure 3.3: Results for Gaussians with full covariance matrices, K = 4, p = 5. (a) Test loglikelihood. (b) Projection of the observation sequence on the first two principal components. Color indicates predicted state.

Figure 4

 4 Figure 4.1: (a) Music score of the excerpt from M. Ravel's Ma Mère l'Oye. Source: (Lostanlen, 2013). (b) Score of one bar of Bach's violin sonata n. 2 (Allegro).

Figure 4

 4 Figure 4.2: Results on Ravel for offline algorithms with different K-means initializations. K = 9. HSMM duration distributions are fixed to N egBin(5, 0.95).

Figure 4

 4 Figure 4.3: Results on the Bach example for offline algorithms. K = 10. HSMM durations: N B(5, 0.2) (mean 20). The last segmentation shown is the ground truth.

Figure 4 .

 4 Figure4.5 compares the online and incremental EM algorithms for HMMs, described in §3.1.1 and §3.3.1, respectively. We show results for a single sequence or the same sequence repeated twice, with K = 10. We can see that both methods give similar results. The online EM approach seems to model transitions better than incremental EM, which gives many short segments, especially in the shorter sequence, while incremental EM seems to discover more states compared to online EM. Another important difference is in the execution time: incremental EM is faster by almost an order of magnitude compared to online EM (0.4s vs 1.5s), and the difference gets bigger as K increases, since the complexities per observation differ by a factor O(K 2 ).

Figure 4

 4 Figure 4.4: Online EM for HMM/HSMM on Bach. K = 10, N B(30, 0.6) (mean 20). (a) Single sequence. (b) Sequence repeated twice.

Figure 4

 4 Figure 4.5: Online vs incremental EM for HMMs. The Bach sequence is repeated twice in (b)

Figure 4 . 6 :

 46 Figure 4.6: Scenes segmentation results. In both examples, two different events are alternated and repeated 3 times. (a) phone ringing and coughing. (b) keys dropping and door slam.

In signal processing, the IS divergence is often used with the power spectrum F (w) = |x(w)|

of a signal x(t), which is a function rather than a vector, but we limit ourselves to vectors in the example for simplicity.2 The Fenchel conjugate of ψ is defined by ψ * (µ) = max x∈Ω µ, xψ(x).

Note that this doesn't necessarily imply convergence to a stationary point of the likelihood, nor does it imply the convergence of the iterates. For more details the convergence properties of EM, see, e.g.,[START_REF] Dempster | Maximum likelihood from incomplete data via the EM algorithm[END_REF];[START_REF] Wu | On the convergence properties of the EM algorithm[END_REF];[START_REF] Xu | On convergence properties of the EM algorithm for gaussian mixtures[END_REF].

This follows from the fact that EM can be seen as a fixed point problem, under a mapping M such that θ (t+1) = M (θ(t) ).

Several different types of HSMMs exist, see[START_REF] Yu | Hidden semi-markov models[END_REF] for details. We will focus on the "explicit duration HMM", as described in[START_REF] Murphy | Hidden semi-markov models (hsmms)[END_REF][START_REF] Guédon | Estimating hidden semi-markov chains from discrete sequences[END_REF]

Acknowledgements

Algorithm 3.2: Online optimization algorithm for the non-probabilistic HMM where new states can be added.

Data: Input stream x t , t = 1, . . ., Bregman divergence D ψ .

Result:

Because of the deterministic transitions in the model, most transitions have infinite cost, and the only possible choices are staying in the same state with z D t-1 = d and z D t = d + 1, or changing state and having z D t = 1.

Incremental EM

An alternative algorithm to perform online updates in EM, in the case of independent observations, is the incremental approach of [START_REF] Neal | A view of the em algorithm that justifies incremental, sparse, and other variants[END_REF]. Rather than relying on a stochastic approximation procedure as in the online EM of Section 3.1, the algorithm performs an incremental minorization-maximization of the log-likelihood just like batch EM, but where the E-step is only taken on a single observation. As we saw in Section 2.3.4, the E-step maximizes a lower bound on the log-likelihood with respect to a distribution q obtained with Jensen's inequality:

We saw that this lower bound is maximized (and made tight) for q(z) = p(z|x; θ) = i p(z i |x i ; θ). We can thus limit ourselves to distributions of the form q(z) = i q i (z i ) where k q i (k) = 1 for