Rémi Flamary
email: remi.flamary@unice.fr

Marco Cuturi
email: mcuturi@i.kyoto-u.ac.jp

Nicolas Courty
email: courty@univ-ubs.fr

Alain Rakotomamonjy

Wasserstein Discriminant Analysis

Wasserstein Discriminant Analysis (WDA) is a new supervised method that can improve classification of high-dimensional data by computing a suitable linear map onto a lower dimensional subspace. Following the blueprint of classical Linear Discriminant Analysis (LDA), WDA selects the projection matrix that maximizes the ratio of two quantities: the dispersion of projected points coming from different classes, divided by the dispersion of projected points coming from the same class. To quantify dispersion, WDA uses regularized Wasserstein distances, rather than cross-variance measures which have been usually considered, notably in LDA. Thanks to the the underlying principles of optimal transport, WDA is able to capture both global (at distribution scale) and local (at samples scale) interactions between classes. Regularized Wasserstein distances can be computed using the Sinkhorn matrix scaling algorithm; We show that the optimization of WDA can be tackled using automatic differentiation of Sinkhorn iterations. Numerical experiments show promising results both in terms of prediction and visualization on toy examples and real life datasets such as MNIST and on deep features obtained from a subset of the Caltech dataset.

Introduction

Dimensionality reduction techniques can convert high-dimensional data with potentially redundant features into low-dimensional vectors [START_REF] Van Der Maaten | Dimensionality reduction: a comparative review[END_REF][START_REF] Christopher | Dimension reduction: A guided tour[END_REF]. In doing so, dimensionality reduction can alleviate ensuing computations, and, ideally, extract relevant features out of noisy ones. Dimensionality reduction techniques come in all flavors: Some techniques only consider linear transformations of the data, whereas some extract features using nonlinear functions; Some of these techniques, such as PCA, are designed for unsupervised settings, whereas some can consider labeled data and fall in the supervised category. We consider in this paper linear and supervised techniques. Within that category, two families of methods stand out: Given a dataset of pairs of vectors and labels {(x i , y i )} i , the goal of Fisher Discriminant Analysis (FDA) and variants is to learn a linear map P : R d → R p , p d, such that the embeddings of these points Px i can be easily discriminated using linear classifiers. Mahalanobis metric learning (MML) methods proceed in the same way, except that the quality of the embedding P is judged by the ability of the k-nearest neighbor algorithm to obtain good classification accuracy. Although both FDA and MML learn a map P, their algorithms are shaped by the classifiers they use in their final phase.

FDA [3, §4.3] attempts to maximize w.r.t. P the sum of all distances ||Px i -Px j || between pairs of samples from different classes c, c while minimizing the sum of all distances ||Px i -Px j || between pairs of samples within the same class c. Because of this, it is well documented that the performance of FDA degrades when class distributions are multimodal. Several variants have been proposed to tackle this problem [3, §12.4]. For instance, a localized version of FDA was proposed in [START_REF] Sugiyama | Dimensionality reduction of multimodal labeled data by local fisher discriminant analysis[END_REF], which boils down to discarding in the computation all pairs of points that are not neighbors.

On the other hand, despite the fact that Mahalanobis metric learning (MML) was originally designed to operate with a k-nearest neighbor classifier, the first techniques that were proposed in that field [START_REF] Xing | Distance metric learning with application to clustering with side-information[END_REF] used a global criterion, based on all pairs of points. Later on, variations that focused instead exclusively on close-range interactions, such as LMNN [START_REF] Kilian | Distance metric learning for large margin nearest neighbor classification[END_REF], were shown to be far more efficient.

Supervised dimensionality approaches stemming either from FDA or MML consider thus either global or local interactions between points, but not both at the same time. We introduce in this work a novel approach that incorporates both global and local constraints. WDA can achieve this blend through the mathematics of optimal transport, which enforce that global and local aspects are taken into account when computing the weight t ij that regulates the influence of the distance ||Px i -Px j ||. Indeed, such weights are decided by (i) making sure that all points in one set are matched to all points in the other set (global constraint); (ii) making sure that points are only matched on average to nearby points, thanks to the optimality of the transportation matrix (local cost). Our method has the added flexibility that it can interpolate between an exclusively global viewpoint (identical, in that case, to FDA), to a more local viewpoint with a global matching constraint (different, in that sense, to that of purely local tools such as LMNN or Local-FDA). In mathematical terms, we adopt the ratio formulation of FDA to maximize the ratio of the regularized Wasserstein distances between inter class populations and between the intra-class population with itself, when these points are considered in their projected space:

max P∈∆ c,c >c W λ (PX c , PX c ) c W λ (PX c , PX c ) (1)
where

∆ = {P = [p 1 , . . . , p p ] | p i ∈ R d , p i 2 =
1 and p i p j = 0 for i = j} is the Stiefel manifold [START_REF] Absil | Optimization algorithms on matrix manifolds[END_REF], the set of orthogonal d × p matrices; PX c is the matrix of projected samples from class c and W λ is the regularized Wasserstein distance proposed in [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]. We recover FDA when the entropic regularization strength is very large (small λ). When, on the contrary, that regularization is weak, our approach tries to split two class distributions by, using loose language, making sure that their optimal matching distance is as large as possible. In that process, each point in one class is paired with another point in another class, and only that distance is maximized. Figure 1 illustrates how inter and intra-class distances are computed, and the difference in how data points interact locally when comparing Wasserstein discriminant analysis (WDA), a global approach such as FDA, and a purely local one such as Local-FDA [START_REF] Sugiyama | Dimensionality reduction of multimodal labeled data by local fisher discriminant analysis[END_REF].

Another strong feature brought by regularized Wasserstein distances is that relations between samples (as given by the optimal transport matrix T) are estimated in the projected space. This is an important difference compared to all previous local approaches which estimate local relations in the original space and make the hypothesis that these relations are unchanged after projection. Estimating local relations in the projected space leads to more robust results in the presence of noisy dimensions, as illustrated in Figure 2 only Gaussian noise. In this case, the 2D space estimated by our method is more discriminant than the local methods that rely on a noisy neighborhood estimation.

In what follows, we provide first background on regularized Wasserstein distances, present our Wasserstein-based discriminant analysis framework in §3 and discuss the optimization problem we tackle, which we solve using automatic differentiation. Numerical experiments that show the relevance of our approach are provided in § 4.

Background on Wasserstein distances

Wasserstein distances, also known as earth mover distances, define a geometry over the space of probability measures using principles from optimal transport theory [START_REF] Villani | Optimal transport: old and new[END_REF]. Recent computational advances [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF][START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF] have made it scale to dimensions relevant to machine learning applications.

Notations and Definitions

Let µ = 1 n i δ xi , ν = 1 m i δ zi be two empirical measures with locations in R d stored in the matrices X = [x 1 , • • • , x n ] and Z = [z 1 , • • • , z m ].
The squared Euclidean distance between samples from µ and ν is defined as M X,Z := [||x i -z j || 2 2 ] ij ∈ R n×m . Let U nm be the polytope of n × m nonnegative matrices such that their row and column marginals are equal to 1 n /n and 1 m /m respectively. Writing 1 n for the n-dimensional vector of ones,

U nm := {T ∈ R n×m + : T1 m = 1 n /n, T T 1 n = 1 m /m}.
Regularized Wassersein distance Let A, B := tr(A T B) be the Frobenius dot-product of matrices. For λ ≥ 0, the regularized Wasserstein distance we adopt in this paper between µ and ν is:

W λ (µ, ν) := W λ (X, Z) := T λ , M X,Z , (2)
where T λ is the solution of an entropy-smoothed optimal transport problem,

T λ := argmin T∈Unm λ T, M XZ -Ω(T), (3)
where Ω(T) is the entropy of T seen as a discrete joint probability distribution, namely Ω(T) := ij t ij log(t ij ). Note that problem (3) can be solved very efficiently using the Sinkhorn fixed point iterations [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]. The solution of the optimization problem can be expressed as:

T = diag(u)e -λM diag(v) = u1 T m • e -λM • 1 n v T , (4)
where • stands for elementwise multiplication of matrices and the exponential is also applied elementwise. The Sinkhorn iterations consist in updating left/right scaling vectors u k and v k of the matrix K = e -λM . These updates take the following form for iteration k:

v k = 1 m /m K T u k-1 , u k = 1 n /n Kv k (5)
with an initialization which will be fixed to u 0 = 1 n . Because it only involves matrix products, the Sinkhorn algorithm can be streamed efficiently on parallel architectures such as GPGPUs.

Optimal transport for machine learning The geometry of the space of probability measures endowed with the Wasserstein metric has been recently considered in applied settings, notably to compute barycenters [START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF][START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF], or generalizations of PCA [START_REF] Seguy | Principal geodesic analysis for probability measures under the optimal transport metric[END_REF]. Wasserstein distances have also been considered to carry out semi-supervised learning [START_REF] Solomon | Wasserstein propagation for semi-supervised learning[END_REF], domain adaptation [START_REF] Courty | Domain adaptation with regularized optimal transport[END_REF], or to define new loss functions for multiclass classifiers [START_REF] Frogner | Learning with a wasserstein loss[END_REF]. A recent application related to ours is that of identifying univariate projections of high-dimensional probability distributions [START_REF] Mueller | Principal differences analysis: Interpretable characterization of differences between distributions[END_REF]. Contrary to our approach, this reference builds upon the closed form formulation of 1D-Wasserstein distances and was not designed to find subspaces of dimension p > 1.

Wasserstein Discriminant Analysis

In this section we discuss the optimization problem (1) and its relation to FDA. Then we propose an efficient approach for computing the gradient of the objective function.

Optimization problem In the remaining, to simplify notations, we define a separate empirical measure for each of the C classes: The samples of each class c are stored in matrices X c ; the number of samples from class c is n c .

Using the definition (2) of the regularized Wasserstein distance, we can write the Wasserstein Discriminant Analysis optimization problem as

max P∈∆ J(P, T(P)) = c,c >c P T P, C c,c c P T P, C c,c (6)
s.t. C c,c = i,j T c,c i,j (x c i -x c j )(x c i -x c j ) T , ∀c, c and T c,c = argmin T∈Un c n c λ T, M PX c ,PX c -Ω(T),
which can be reformulated as max P∈∆ J(P, T(P))

s.t. T(P) = argmin T∈Un c n c E(T, P)

where T = {T c,c } c,c contains all the transport matrices between classes. The objective function J and the inner problem function E defined as

J(P, T(P)) = P T P, C b P T P, C w , E(T) = c,c>=c λ T c,c , M PX c ,PX c -Ω(T c,c ) (9)
where

C b = c,c >c C c,c and C w = c C c,c
are the between and within cross-covariance matrices that depend on T(P). Optimization problem ( 7)-( 8) is a classic bilevel optimization problem and it has a particular form that allows its resolution using gradient descent [START_REF] Colson | An overview of bilevel optimization[END_REF]. Indeed since T (P) is smooth and optimization problem ( 8) is strictly convex one can compute the gradient directly w.r.t. P using the chain rule as follows

∇ P J(P, T(P)) = ∂J(P, T) ∂P + c,c ≥c ∂J(P, T) ∂T c,c ∂T c,c ∂P (10)
The first term in gradient [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF] suppose that T is constant and can be computed (Eq. 94-95 [START_REF] Brandt Petersen | The matrix cookbook[END_REF]) as

∂J(P, T) ∂P = P 2 σ 2 w C b - 2σ 2 b σ 4 w C w (11)
with σ 2 w = P T P, C w and σ 2 b = P T P, C b . In order to compute the second term in (10), we will separate the cases when c = c and c = c since it corresponds to their position in the fraction. Their partial derivative is obtained directly from the scalar product and is a weighted vectorization of the transport cost matrix

∂J(P, T) ∂T c,c =c = vec 1 σ 2 b M PX c ,PX c , ∂J(P, T) ∂T c,c = -vec σ 2 w σ 4 b M PX c ,PX c (12)
We will see in the remaining that the main difficulty stands in computing the Jacobian ∂T c,c ∂P since the optimal transport matrix is not available as a closed form. In the sequel, we solve this problem by means of an automatic differentiation approach wrapped around the Sinkhorn fixed point iteration algorithm.

Once the gradient is computed, we can use classical manifold optimization tools such as projected gradient [START_REF] Schmidt | Minconf-projection methods for optimization with simple constraints in matlab[END_REF] or a trust region algorithm as implemented in Manopt [START_REF] Boumal | Manopt, a matlab toolbox for optimization on manifolds[END_REF]. The latter toolbox includes tools to optimize over the Stiefel manifold, notably automatic conversions from Euclidean to Riemannian gradients.

Relation to Fisher Discriminant Analysis We consider the limit behavior of WDA as λ approaches 0. In this case, we can see from Eq. ( 3) that the matrix T does not depends on the data. The solution T for each Wasserstein distance is the matrix that maximizes entropy, namely the uniform probability distribution T = 1 nm 1 n,m . The cross-covariance matrices become thus

C c,c = 1 n c n c i,j (x c i -x c j )(x c i -x c j ) T
and the matrices C w and C b correspond then to intra-and inter-class covariances as used in FDA. Since these matrices do not depend on P, the optimization problem (1) boils down to the usual Rayleigh quotient which can be solved using a generalized eigendecomposition of C -1 w C b . When λ > 0, the smoothed optimal transport promotes instead cross-covariance matrices that tend to favor local relations as illustrated in Figure 1.

Also note that the regularization parameter (λ < ∞) is required in the current formulation [START_REF] Van Der Maaten | Dimensionality reduction: a comparative review[END_REF]. Indeed, the classical Wasserstein distance of a discrete distribution with itself is 0. The entropic regularization ensures that the mass of each sample is split among neighbors hence promoting a local variance estimation and avoiding a division by 0 in [START_REF] Van Der Maaten | Dimensionality reduction: a comparative review[END_REF]. If one wants to use the non-regularized Wasserstein distance, then the intra-class distance should be computed between two non-overlapping splits of the examples of each class which would also leads to non-zero distance.

Automatic Differentiation A possible way to compute the Jacobian ∂T c,c /∂P is to use the implicit function theorem as proposed in [START_REF] Bengio | Gradient-based optimization of hyperparameters[END_REF] for hyperparameter estimation. Although we have tried it, as detailed in supplementary material, this computational approach fails to scale because it requires inverting a very large matrix, and also assumes that the exact optimal transport T can be obtained, despite the fact that one only gets it after a finite number of Sinkhorn iterations. Following the gist of [START_REF] Bonneel | Wasserstein barycentric coordinates: Histogram regression using optimal transport[END_REF] (which do not differentiate Sinkhorn iterations but a more complex fixed point iteration designed to compute Wasserstein barycenters), we propose in this section to differentiate the transportation matrices obtained after running exactly L Sinkhorn iterations, with a predefined L. Writing T k (P), for the solution obtained after k iterations as a function of P for a given c, c pair,

T k (P) = diag(u k )e -λM diag(v k )
where M is the distance matrix induced by P. T L (P) can then be directly differentiated:

∂T k ∂P = ∂[u k 1 T m ] ∂P • e -λM • 1 n v k T + u k 1 T m • ∂e -λM ∂P • 1 n v k T + u k 1 T m • e -λM • ∂[1 n v k T ] ∂P (13)
Note that the recursion occurs as u k depends on v k whose is also related to u k-1 . The Jacobians that we need can then be obtained from Equation [START_REF] Xing | Distance metric learning with application to clustering with side-information[END_REF]. For instance, the gradient of one component of u k 1 T m at the j-th line is

∂u k j ∂P = - 1/n [Kv k ] 2 j i ∂K j,i ∂P v k i + i K j,i ∂v k i ∂P ,
while for v k , we have

∂v k j ∂P = - 1/m [K T u k-1 ] 2 j i ∂K i,j ∂P u k-1 j + i K i,j ∂u k-1 j ∂P ,
∂K i,j ∂P = -2K i,j P(x i -x j )(x i -x j ) T .
The Jacobian ∂T k ∂P can be thus obtained by keeping track of all the Jacobians at each iteration and then by successive applying those equations. This approach is far cheaper than the implicit function theorem approach. Indeed, in this case, the computation of ∂T ∂P is dominated by the complexity of computing ∂K ∂P whose costs for one iteration is O(pn 2 d 2 ) for n = m. The complexity is then linear in L and quadratic in n.

Numerical experiments

In this section we illustrate WDA on three datasets. First, we evaluate our approach on a simple simulated dataset with a 2-dimensional discriminant subspace that is known a priori. Then we apply WDA on MNIST and Caltech datasets. Note that in the spirit of reproducible research all the Matlab code will be made available to the community on Github and mloss.org upon publication.

Practical implementation

In order to make the method less sensitive to the dimension and scaling of the data, we propose to use a pre-computed adaptive regularization parameter λ c,c for each Wasserstein distances in [START_REF] Van Der Maaten | Dimensionality reduction: a comparative review[END_REF]. In practice we initialize P with the PCA projection, and we use it to find the following λ c,c = λ( 1

ncn c i,j Px c i -Px c j 2
) -1 between class c and c . These values are computed a priori and fixed in the remaining iterations, buth they will promote a similar regularization between inter and intra-class distances.

Simulated dataset This dataset has been designed to evaluate the ability of a subspace method to estimate a linear subspace when the classes are non-linearly separable. It is a 3-class problem in dimension d = 10 where only the first two dimensions are discriminant and the remaining 8 contain Gaussian noise with the same variance as the discriminant dimensions. In the 2 discriminant features each class contains two modes that have centers positioned so that all the classes have the same mean at [0, 0]. An example of samples in the discriminant dimensions is shown in the left part of Figure 2.

First, we depict in Figure 2, the samples projected in several subspaces with different approaches. We can see that for this dataset WDA seems to estimate the subspace with the best generalization property. To verify this trend, we compute the prediction errors of all the linear subspace methods with a K-Nearest-Neighbors classifier (KNN) for n = 100 training examples and n t = 5000 test examples in Figure 3. In this simulation, all prediction errors are averaged over 20 data generations and the neighbors parameters of LMNN and LFDA have been selected empirically to maximize performances (K = 5 for LMNN and K = 1 for LFDA). We can see in the left part of the figure that WDA and to a lesser extent LMNN can estimate a relevant subspace that is robust to the choice of K. LFDA clearly outperforms PCA, Orig. and LDA since it can handle non-linearity but fails to attain the performances of LMNN and WDA. In the center plot, we illustrate the sensitivity of WDA w.r.t. the regularization parameter λ. WDA returns the best performance on almost a full order of magnitude which suggests that a coarse validation can be performed in practice. The right part of Figure 3 shows the performance of the WDA with fixed point gradient (WDA FP) for different values of the number of inner iterations L. We can see that even if this parameter leads to different performances on large values of λ it still leads to the best overall performances even for small L. The variability between the curves comes from the non-convexity of the objective function, since each descent direction might lead to a different stationary point.

MNIST dataset. Our aim on this dataset is to measure how robust our approach is when only few training samples are available in high-dimension. To this end, we draw n = 1000 samples for training and report the KNN prediction error as a function of k for the different subspace methods when projecting onto p = 2 and p = 10 dimensions (resp. left and middle plots of Figure 4). In order to leverage the impact of a random choice of the initial distribution, the reported scores are averages of 20 realizations of the same experiment. We also limit the analysis to L = 10 as the number of fixed point iterations. As expected, when p = 2 there is an important loss of information which leads to large classification errors compared to a KNN in the original space (Orig.). Nevertheless, WDA seems to be the method that finds the best 2D subspace. When p = 10, WDA finds a better subspace than the original space which suggests that all the discriminant information available in the training dataset has been extracted efficiently. Conversely, LFDA and LMNN struggle in finding a relevant subspace in this configuration. We believe that this is due to the small number of training samples and the fact that those approaches look at local relationship whereas the Wasserstein distance is a global measure that can be regularized to better avoid overfitting and to better capture local informations. In addition to better prediction performance, we want to emphasize that in this configuration WDA leads to a dramatic compression of the data from 784 to 10 features while keeping the discriminative information.

To gain a better understanding of the corresponding embedding, we project the data from the 10dimensional space to a 2-dimensional one using tSNE, that is a classic nonlinear projection method used for visualizing high dimensional samples for instance in deep architectures [START_REF] Van Der Maaten | Visualizing data using t-sne[END_REF]. In order to We finally compare in Figure 6 the different reconstructions induced by the adjoint projection operator. One can see that contrary to the other discriminative methods, our method preserves a coherent spatial information related to the class, which illustrates again the generalization capability of WDA.

The means that the estimated subspace is both generative as for PCA and discriminative thanks to the use of a proper divergence between discrete distributions.

Caltech dataset

In this experiment, we use a subset described in [START_REF] Donahue | DeCAF: a deep convolutional activation feature for generic visual recognition[END_REF] from the Caltech-256 image collection [START_REF] Griffin | Caltech-256 Object Category Dataset[END_REF]. The dataset uses features that are the output of the DeCAF deep learning architecture [START_REF] Donahue | DeCAF: a deep convolutional activation feature for generic visual recognition[END_REF]. More precisely, they are extracted as the sparse activation of the neurons from the 6th fully connected layer of a convolutional network trained on imageNet and then fine-tuned for the considered visual recognition task. As such, they form vectors of 4096 dimensions. In this setting, 500 images are considered for training, and the remaining portion of the dataset for testing (623 images). There are 9 different classes in this dataset. We examine in this experiment how the proposed dimensionality reduction performs when changing the subspace dimensionality (right plot in Figure 4). For this problem, the regularization parameter λ of WDA was empirically set to 10 -2 and we use early stopping at 10 iterations for WDA. The K from KNN was to set to 3 which is a common standard setting for this classifier. The reported results are averaged over 10 realizations of the same experiment. When p ≥ 5, WDA already finds a subspace which gathers all the discriminative information from the original space. In this experiment LMNN reaches a better subspace for small p values but WDA is the best performing method for p ≥ 6. Those results highlight the potential interest for using WDA as linear dimensionality reduction layers in neural-nets architecture.

Conclusion

This work presents the Wasserstein Discriminant Analysis, a new and original linear discriminant subspace estimation method. Based on the framework of Wasserstein distances, which measure a global similarity between empirical distributions, WDA operates by separating distributions of different classes in the subspace, while maintaining a coherent structure at a class level. To this extent, the use of regularization in the Wasserstein formulation allows to effectively bridge a gap between a global coherency and the local structure of the class manifold. This comes at a cost of a difficult optimization of a bilevel program, for which we showed that a naive computation of the gradient is achievable but computationally too costly in real settings. We proposed instead an efficient method that relies on differentiating the map obtained after a finite number of iterations of the Sinkhorn algorithm. Numerical experiments show that the method performs well on a variety of features, including those obtained with a deep neural architecture. Future work will consider stochastic versions of the same approach in order to enhance further the ability of the method to handle large volume of high-dimensional data.

Figure 1 :

 1 Figure 1: Weights used for inter/intra class variances for FDA, Local FDA and WDA for different regularizations λ. Only weights for two samples from class 1 are shown. The color of the link darkens as the weight grows. FDA computes a global variance with uniform weight on all pairwise distances, whereas LFDA focuses only on samples that lie close to each other. WDA relies on an optimal transport matrix T that matches all points in one class to all other points in another class. WDA has both a global (due to matching constraints) and local (due to transportation cost minimization) outlook on the problem, with a tradeoff controlled by the regularization strength λ.

2 WDAFigure 2 :

 22 Figure 2: Illustration of subspace learning methods on a nonlinearly separable 3-class toy example of dimension d = 10 with 2 discriminant features (shown on the left) and 8 Gaussian noise features. The projection onto p = 2 for training (up) and test data (down) are reported for several subspace estimation methods.

Figure 3 :

 3 Figure 3: Prediction error on the simulated dataset with projection dimension p = 2. (left) Error for a varying K in the KNN classifier (middle) Evolution of WDA performance for different regularization parameter values λ (right) comparison of WDA performances as a function of the number of fixed point iterations.

Figure 4 :Figure 5 :

 45 Figure 4: (left) Averaged prediction error on MNIST with projection dimension p = 2 and (middle) p = 10. (right) Averaged prediction error on the Caltech-256 dataset as a function of p .

Figure 6 :

 6 Figure 6: Image reconstruction after projection in the subspaces for test samples. The title on top of each column corresponds to the class of the sample.


Rémi Flamary
email: remi.flamary@unice.fr

Cédric Févotte
email: cedric.fevotte@irit.fr

Nicolas Courty
email: courty@univ-ubs.fr

Valentin Emiya
email: valentin.emiya@lif.univ-mrs.fr

Optimal spectral transportation with application to music transcription

Many spectral unmixing methods rely on the non-negative decomposition of spectral data onto a dictionary of spectral templates. In particular, state-of-the-art music transcription systems decompose the spectrogram of the input signal onto a dictionary of representative note spectra. The typical measures of fit used to quantify the adequacy of the decomposition compare the data and template entries frequency-wise. As such, small displacements of energy from a frequency bin to another as well as variations of timbre can disproportionally harm the fit. We address these issues by means of optimal transportation and propose a new measure of fit that treats the frequency distributions of energy holistically as opposed to frequency-wise. Building on the harmonic nature of sound, the new measure is invariant to shifts of energy to harmonically-related frequencies, as well as to small and local displacements of energy. Equipped with this new measure of fit, the dictionary of note templates can be considerably simplified to a set of Dirac vectors located at the target fundamental frequencies (musical pitch values). This in turns gives ground to a very fast and simple decomposition algorithm that achieves state-of-the-art performance on real musical data.

Context

Many of nowadays spectral unmixing techniques rely on non-negative matrix decompositions. This concerns for example hyperspectral remote sensing (with applications in Earth observation, astronomy, chemistry, etc.) or audio signal processing. The spectral sample v n (the spectrum of light observed at a given pixel n, or the audio spectrum in a given time frame n) is decomposed onto a dictionary W of elementary spectral templates, characteristic of pure materials or sound objects, such that v n ≈ Wh n . The composition of sample n can be inferred from the non-negative expansion coefficients h n . This paradigm has led to state-of-the-art results for various tasks (recognition, classification, denoising, separation) in the aforementioned areas, and in particular in music transcription, the central application of this paper.

In state-of-the-art music transcription systems, the spectrogram V (with columns v n ) of a musical signal is decomposed onto a dictionary of pure notes (in so-called multi-pitch estimation) or chords. V typically consists of (power-)magnitude values of a regular short-time Fourier transform [START_REF] Smaragdis | Non-negative matrix factorization for polyphonic music transcription[END_REF]. It may also consists of an audio-specific spectral transform such as the Melfrequency transform, like in [START_REF] Vincent | Adaptive harmonic spectral decomposition for multiple pitch estimation[END_REF], or the Q-constant based transform, like in [START_REF] Oudre | Chord recognition by fitting rescaled chroma vectors to chord templates[END_REF]. The success of the transcription system depends of course on the adequacy of the time-frequency transform & the dictionary to represent the data V. In particular, the matrix W must be able to accurately represent a diversity of real notes. It may be trained with individual notes using annotated data [START_REF] Boulanger-Lewandowski | Discriminative non-negative matrix factorization for multiple pitch estimation[END_REF], have a parametric form [START_REF] Rigaud | A parametric model and estimation techniques for the inharmonicity and tuning of the piano[END_REF] or be learnt from the data itself using a harmonic subspace constraint [START_REF] Vincent | Adaptive harmonic spectral decomposition for multiple pitch estimation[END_REF].

One important challenge of such methods lies in their ability to cope with the variability of real notes. A simplistic dictionary model will assume that one note characterised by fundamental frequency ν 0 (e.g., ν 0 = 440 Hz for note A 4 ) will be represented by a spectral template with non-zero coefficients placed at ν 0 and at its multiples (the harmonic frequencies). In reality, many instruments, such as the piano, produce musical notes with either slight frequency misalignments (so-called inharmonicities) with respect to the theoretical values of the fundamental and harmonic frequencies, or amplitude variations at the harmonic frequencies with respect to recording conditions or played instrument (variations of timbre). Handling these variabilities by increasing the dictionary with more templates is typically unrealistic and adaptive dictionaries have been considered in [START_REF] Vincent | Adaptive harmonic spectral decomposition for multiple pitch estimation[END_REF][START_REF] Rigaud | A parametric model and estimation techniques for the inharmonicity and tuning of the piano[END_REF]. In these papers, the spectral shape of the columns of W is adjusted to the data at hand, using specific time-invariant semi-parametric models. However, the note realisations may vary in time, something which is not handled by these approaches. This work presents a new spectral unmixing method based on optimal transportation (OT) that is fully flexible and remedies the latter difficulties. Note that [START_REF] Typke | Searching notated polyphonic music using transportation distances[END_REF] have previously applied OT to notated music (e.g., score sheets) for search-by-query in databases while we address here music transcription from audio spectral data.

A relevant baseline: PLCA

Before presenting our contributions, we start by introducing the PLCA method of [START_REF] Smaragdis | A probabilistic latent variable model for acoustic modeling[END_REF] which is heavily used in audio signal processing. It is based on the Probabilistic Latent Semantic Analysis (PLSA) of Hofmann (2001) (used in text retrieval) and is a particular form of nonnegative matrix factorisation (NMF). Simplifying a bit, in PLCA the columns of V are normalised to sum to one. Each vector v n is then treated as a discrete probability distribution of "frequency quanta" and is approximated as V ≈ WH. The matrices W and H are of size M × K and K × N , respectively, and their columns are constrained to sum to one. As a result, the columns of the approximate V = WH sum to one as well and each distribution vector v n is as such approximated by the counterpart distribution vn in V. Under the assumption that W is known, the approximation is found by solving the optimisation problem defined by

min H≥0 D KL (V|WH) s.t ∀n, h n 1 = 1, (1)
where

D KL (v|v) = i v i log(v i /v i )
is the KL divergence between discrete distributions, and by extension

D KL (V| V) = n D KL (v n |v n ).
An important characteristic of the KL divergence is its separability with respect to the entries of its arguments. It operates a frequency-wise comparison in the sense that, at every frame n, the spectral coefficient v in at frequency i is compared to its counterpart vin , and the results of the comparisons are summed over i. In particular, a small displacement in the frequency support of one observation may disproportionally harm the divergence value. For example, if v n is a pure note with fundamental frequency ν 0 , a small inharmonicity that shifts energy from ν 0 to an adjacent frequency bin will unreasonably increase the divergence value, when v n is compared with a purely harmonic spectral template with fundamental frequency ν 0 . As explained in Section 1 such local displacements of frequency energy are very common when dealing with real data. A measure of fit invariant to small perturbations of the frequency support would be desirable in such a setting, and this is precisely what OT can bring.

Elements of optimal transportation

Given a discrete probability distribution v (a non-negative real-valued column vector of dimension M and summing to one) and a target distribution v (with same properties), OT computes a transportation matrix T belonging to the set

Θ def = {T ∈ R M ×M + |∀i, j = 1, . . . , N, M j=1 t ij = v i , M i=1 t ij = vj }.
T establishes a bi-partite graph connecting the two distributions. In simple words, an amount (or, in typical OT parlance, a "mass") of every coefficient of vector v is transported to an entry of v. The sum of transported amounts to the j th entry of v must equal vj . The value of t ij is the amount transported from the i th entry of v to the j th entry of v. In our particular setting, the vector v is a distribution of spectral energies v 1 , . . . , v M at sampling frequencies f 1 , . . . , f M .

Without additional constraints, the problem of finding a non-negative matrix T ∈ Θ has an infinite number of solutions. As such, OT takes into account the cost of transporting an amount from the i th entry of v to the j th entry of v, denoted c ij (a non-negative real-valued number). Endorsed with this cost function, OT involves solving the optimisation problem defined by

min T J(T|v, v, C) = ij c ij t ij s.t T ∈ Θ, ( 2
)
where C is the non-negative square matrix of size M with elements c ij . Eq. ( 2) defines a convex linear program. The value of the function J(T|v, v, C) at its minimum is denoted D C (v|v). When C is a symmetric matrix such that c ij = f i -f j p p , where we recall that f i and f j are the frequencies in Hertz indexed by i and j, D C (v|v) defines a metric (i.e., a symmetric divergence that satisfies the triangle inequality) coined Wasserstein distance or earth mover's distance [START_REF] Rubner | A metric for distributions with applications to image databases[END_REF][START_REF] Villani | Optimal transport: old and new[END_REF]. In other cases, in particular when the matrix C is not even symmetric like in the next section, D C (v|v) is not a metric in general, but is still a valid measure of fit. For generality, we will refer to it as the "OT divergence".

By construction, the OT divergence can explicitly embed a form of invariance to displacements of support, as defined by the transportation cost matrix C. For example, in the spectral decomposition setting, the matrix with entries of the form c ij = (f i -f j ) 2 will increasingly penalise frequency displacements as the distance between frequency bins increases. This precisely remedies the limitation of the separable KL divergence presented in Section 2. As such, the next section addresses variants of spectral unmixing based on the Wasserstein distance.

Optimal spectral transportation (OST)

Unmixing with OT. In light of the above discussion, a direct solution to the sensibility of PLCA to small frequency displacements consists in replacing the KL divergence with the OT divergence. This amounts to solving the optimisation problem given by min

H≥0 D C (V|WH) s.t ∀n, h n 1 = 1, (3)
where

D C (V| V) = n D C (v n |v n ),
W is fixed and populated with pure note spectra and C penalises large displacements of frequency support. This approach is a particular case of NMF with the Wasserstein distance, which has been considered in a face recognition setting by [START_REF] Sandler | Nonnegative matrix factorization with earth mover's distance metric for image analysis[END_REF], with subsequent developments by [START_REF] Zen | Simultaneous ground metric learning and matrix factorization with earth mover's distance[END_REF] and [START_REF] Benamou | Iterative Bregman projections for regularized transportation problems[END_REF]. This approach is relevant to our spectral unmixing scenario but as will be discussed in Section 5 is on the downside computationally intensive. It also requires the columns of W to be set to realistic note templates, which is still constraining. The next two sections describes a computationally more friendly approach which additionally removes the difficulty of choosing W appropriately.

Harmonic-invariant transportation cost. In the approach above, the harmonic modelling is conveyed by the dictionary W (consisting of comb-like pure note spectra) and the invariance to small frequency displacements is introduced via the matrix C. In this section we propose to model both harmonicity and local invariance through the transportation cost matrix C. Loosely speaking, we want to define a class of equivalence between musical spectra, that takes into account their inherent harmonic nature. As such, we essentially impose that a harmonic frequency (i.e., a close multiple of its fundamental) can be considered equivalent to its fundamental, the only target of multi-pitch estimation. As such, we assume that a mass at one frequency can be transported to a divisor frequency with no cost. In other words, a mass at frequency f i can be transported with no cost to f i /2, f i /3, f i /4, and so on until sampling resolution. One possible cost matrix that embeds this property is

c ij = min q=1,...,qmax (f i -qf j ) 2 + δ q =1 , (4)
where q max is the ceiling of f i /f j and is a small value. The term δ q =1 favours the discrimination of octaves. Indeed, it penalises the transportation of a note of fundamental frequency 2ν 0 or ν 0 /2 to the spectral template with fundamental frequency ν 0 , which would be costless without this additive term. Let us denote by C h the transportation cost matrix defined by Eq. ( 4). to the more standard quadratic cost C 2 defined by c ij = (f i -f j ) 2 . With the quadratic cost, only local displacements are permissible. In contrast, the harmonic-invariant cost additionally permits larger displacements to divisor frequencies, improving robustness to variations of timbre besides to inharmonicities.

v v 1 v 2 v 3 Measure of fit D 2 DKL D C 2 D C h D(
Dictionary of Dirac vectors. Having designed an OT divergence that encodes inherent properties of musical signals, we still need to choose a dictionary W that will encode the fundamental frequencies of the notes to identify. Typically, these will consist of the physical frequencies of the 12 notes of the chromatic scale (from note A to note G, including half-tones), over several octaves. As mentioned in Section 1, one possible strategy is to populate W with spectral note templates. However, as also discussed, the performance of the resulting unmixing method will be capped by the representativeness of the chosen set of templates.

A most welcome consequence of using the OT divergence built on the harmonic-insensitive cost matrix C h is that we may use for W a mere set of Dirac vectors placed at the fundamental frequencies ν 1 , . . . , ν K of the notes to identify and separate. Indeed, under the proposed setting, a real note spectra (composed of one fundamental and multiple harmonic frequencies) can be transported with no cost to its fundamental. Similarly, a spectral sample composed of several notes can be transported to mixture of Dirac vectors placed at their fundamental frequencies. This simply eliminates the problem of choosing a representative dictionary! This very appealing property is illustrated in Fig. 2. Furthermore, the particularly simple structure of the dictionary leads to a very efficient unmixing algorithm, as explained in the next section. In the following, the unmixing method consisting of the combined use of the harmonic-invariant cost matrix C h and of the dictionary of Dirac vectors will be coined "optimal spectral transportation" (OST).

At this level, we assume for simplicity that the set of K fundamental frequencies {ν 1 , . . . , ν K } is contained in the set of sampled frequencies {f 1 , . . . , f M }. This means that w k (the k th column of W) is zero everywhere except at some entry i such that f i = ν k where w ik = 1. This is typically not the case in practice, where the sampled frequencies are fixed by the sampling rate, of the form f i = 0.5(i/T )f s , and where the fundamental frequencies ν k are fixed by music theory. Our approach can actually deal with such a discrepancy and this will be explained later in Section 5.

Optimisation

OT unmixing with linear programming. We start by describing optimisation for the state-of-theart OT unmixing problem described by Eq. ( 3) and proposed by [START_REF] Sandler | Nonnegative matrix factorization with earth mover's distance metric for image analysis[END_REF]. First, since the objective function is separable with respect to samples, the optimisation problem decouples with respect to the activation columns h n . Dropping the sample index n and combining Eqs. ( 2) and (3), optimisation thus reduces to solving for every sample a problem of the form

min h≥0,T≥0 T, C = ij t ij c ij s.t. T1 M = v, T 1 M = Wh, (5)
where 1 M is a vector of dimension M containing only ones and •, • is the Frobenius inner product.

Vectorising the variables T and h into a single vector of dimension M 2 + K, problem (5) can be turned into a canonical linear program. Because of the large dimension of the variable (typically in the order of 10 5 ), resolution can however be very demanding, as will be shown in experiments.

Optimisation for OST. We now assume that W is a set of Dirac vectors as explained at the end of Section 5. We also assume that K < M , which is the usual scenario. Indeed, K is typically in the order of a few tens, while M is in the order of a few hundreds. In such a setting v = Wh contains by design at most K non-zero coefficients, located at the entries such that f i = ν k . We denote this set of frequency indices by S. Hence, for j / ∈ S, we have vj = 0 and thus i t ij = 0, by the second constraint of Eq. ( 5). Additionally, by the non-negativity of T this also implies that T has only K non-zero columns, indexed by j ∈ S. Denoting by T this subset of columns, and by C the corresponding subset of columns of C, problem (5) reduces to

min h≥0, T≥0 T, C s.t. T1 K = v, T 1 M = h. (6)
This is an optimisation problem of significantly reduced dimension (M + 1)K. Even more appealing, the problem has a simple closed-form solution. Indeed, the variable h has a virtual role in problem (6). It only appears in the second constraint, which de facto becomes a free constraint. Thus problem ( 6) can be solved with respect to T regardless of h, and h is then simply obtained by summing the columns of T at the solution. Now, the problem

min T≥0 T, C s.t. T1 K = v (7)
decouples with respect to the rows ti of T, and becomes, ∀i = 1, . . . , M ,

min ti ≥0 k tik cik s.t. k tik = v i . ( 8
)
The solution is simply given by tik i = v i for k i = arg min k {c ik }, and tik = 0 for k = k i . Introducing the labelling matrix L which is everywhere zero except for indices (i, k i ) where it is equal to 1, the solution to OST is trivially given by ĥ = L v. Thus, under the specific assumption that W is a set of Dirac vectors, the challenging problem (5) has been reduced to an effortless assignment problem to solve for T and a simple sum to solve for h. Note that the algorithm is independent of the particular structure of C. In the end, the complexity per frame of OST reduces to O(M ), which starkly contrasts with the complexity of PLCA, in the order O(KM ) per iteration.

In Section 5, we assumed for simplicity that the set of fundamental frequencies {ν k } k was contained in the set of sampled frequencies {f i } i . As a matter of fact, this assumption can be trivially lifted in the proposed setting of OST. Indeed, we may construct the cost matrix C (of dimensions M × K) by replacing the target frequencies f j in Eq. ( 4) by the theoretical fundamental frequencies ν k . Namely, we may simply set the coefficients of C to be c ik = min q (f i -qν k ) 2 + δ q =1 , in the implementation. Then, the matrix T indicates how each sample v is transported to the Dirac vectors placed at fundamental frequencies {ν k } k , without the need for the actual Dirac vectors themselves, which elegantly solves the frequency sampling problem.

OST with entropic regularisation (OST e ). The procedure described above leads to a winnertakes-all transportation of all of v i to its cost-minimum target entry k i . We found it useful in practice to relax this hard assignment and distribute energies more evenly by using the entropic regularisation of [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transportation[END_REF]. It consists of penalising the fit T, C in Eq. ( 6) with an additional term Ω e ( T) = ik tik log( tik ), weighted by the hyper-parameter λ e . The negentropic term Ω e ( T) promotes the transportation of v i to several entries, leading to a smoother estimate of T. As explained in the supplementary material, one can show that the negentropy-regularised problem is a Bregman projection [START_REF] Benamou | Iterative Bregman projections for regularized transportation problems[END_REF] and has again a closed-form solution ĥ = L e v where L e is the M × K matrix with coefficients l ik = exp(-c ik /λ e )/ p exp(-c ip /λ e ). Limiting cases λ e = 0 and λ e = ∞ return the unregularised OST estimate and the maximum-entropy estimate h k = 1/K, respectively. Because L e becomes a full matrix, the complexity per frame of OST e becomes O(KM ).

OST with group regularisation (OST g ). We have explained above that the transportation matrix T has a strong group structure in the sense that it contains by construction M -K null columns, and that only the subset T needs to be considered. Because a small number of the K possible notes will be played at every time frame, the matrix T will additionally have a significant number of null columns. This heavily suggests using group-sparse regularisation in the estimation of T.

As such, we also consider problem ( 6) penalised by the additional term Ω g ( T) = k t k 1 which promotes group-sparsity at column level [START_REF] Huang | A group bridge approach for variable selection[END_REF]. Unlike OST or OST e , OST g does not offer a closed-form solution. Following [START_REF] Courty | Domain adaptation with regularized optimal transport[END_REF], a majorisation-minimisation procedure based on the local linearisation of Ω g ( T) can be employed and the details are given in the supplementary material. The resulting algorithm consists in iteratively applying unregularised OST, as of Eq. ( 6), with the iteration-dependent transportation cost matrix C (iter) = C + R (iter) , where R (iter) is the M × K matrix with coefficients r

(iter) ik = 1 2 t (iter) k -1 2 1 .
Note that the proposed group-regularisation of T corresponds to a sparse regularisation of h. This is because

h k = t k 1 and thus, Ω g ( T) = k √ h k .
Finally, note that OST e and OST g can be implemented simultaneously, leading to OST e+g , by considering the optimisation of the doubly-penalised objective function T, C + λ e Ω e ( T) + λ g Ω g ( T), addressed in the supplementary material.

Experiments

Toy experiments with simulated data. In this section we illustrate the robustness, the flexibility and the efficiency of OST on two simulated examples. The top plots of Fig. 3 display a synthetic dictionary of 8 harmonic spectral templates, referred to as the "harmonic dictionary". They have been generated as Gaussian kernels placed at a fundamental frequency and its multiples, and using exponential dampening of the amplitudes. As everywhere in the paper, the spectral templates are normalised to sum to one. Note that the 8th template is the upper octave of the first one. We compare the unmixing performance of five methods in two different scenarios. The five methods are as follows. PLCA is the method described in Section 2, where the dictionary W is the harmonic dictionary. Convergence is stopped when the relative difference of the objective function between two iterations falls below 10 -5 or the number of iterations (per frame) exceeds 1000. OT h is the unmixing method with the OT divergence, as in the first paragraph of Section 5, using the harmonic transportation cost matrix C h and the harmonic dictionary. OST is like OT h , but using a dictionary of Dirac vectors (placed at the 8 fundamental frequencies characterising the harmonic dictionary). OST e , OST g and OST e+g are the regularised variants of OST, described at the end of Section 5. The iterative procedure in the group-regularised variants is run for 10 iterations (per frame).

In the first experimental scenario, reported in Fig. 3 (a), the data sample is generated by mixing the 1st and 4th elements of the harmonic dictionary, but introducing a small shift of the true fundamental frequencies (with the shift being propagated to the harmonic frequencies). This mimics the effect of possible inharmonicities or of an ill-tuned instrument. The middle plot of Fig. 3 (a), displays the generated sample, together with the "theoretical sample", i.e., without the frequencies shift. This shows how a slight shift of the fundamental frequencies can greatly impact the overall spectral distribution. The bottom plot displays the true activation vector and the estimates returned by the five methods. The table reports the value of the (arbitrary) error measure ĥh true 1 together with the run time (on an average desktop PC using a MATLAB implementation) for every method. The results show that group-regularised variants of OST lead to best performance with very light computational burden, and without using the true harmonic dictionary. In the second experimental scenario, reported in Fig. 3 (b), the data sample is generated by mixing the 1st and 6th elements of the harmonic dictionary, with the right fundamental and harmonic frequencies, but where the spectral amplitudes at the latters do not follow the exponential dampening of the template dictionary (variation of timbre).

Here again the group-regularised variants of OST outperforms the state-of-the-art approaches, both in accuracy and run time.

Transcription of real musical data. We consider in this section the transcription of a selection of real piano recordings, obtained from the MAPS dataset [START_REF] Emiya | Multipitch estimation of piano sounds using a new probabilistic spectral smoothness principle[END_REF]. The data comes with a ground-truth binary "piano-roll" which indicates the active notes at every time. The note fundamental frequencies are given in MIDI, a standard musical integer-valued frequency scale that matches the keys of a piano, with 12 half-tones (i.e., piano keys) per octave. The spectrogram of each recording is computed with a Hann window of size 93-ms and 50% overlap (f s = 44.1Hz). The columns (time frames) are then normalised to produce V. Each recording is decomposed with PLCA, OST and OST e , with K = 60 notes (5 octaves). Half of the recording is used for validation of the hyper-parameters and the other half is used as test data. For PLCA, we validated 4 and 3 values of the width and amplitude dampening of the Gaussian kernels used to synthesise the dictionary. For OST, we set = q 0 in Eq. ( 4), which was found to satisfactorily improve the discrimination of octaves increasingly with frequency, and validated 5 orders of magnitude of 0 . For OST e , we additionally validated 4 orders of magnitude of λ e . Each of the three methods returns an estimate of H. The estimate is turned into a 0/1 piano-roll by only retaining the support of its P n maximum entries at every frame n, where P n is the ground-truth number of notes played in frame n. The estimated piano-roll is then numerically compared to its ground truth using the F-measure, a global recognition measure which accounts both for precision and recall and which is bounded between 0 (critically wrong) and 1 (perfect recognition). Our evaluation framework follows standard practice in music transcription evaluation, see for example [START_REF] Daniel | Perceptually-based evaluation of the errors usually made when automatically transcribing music[END_REF]. As detailed in the supplementary material, it can be shown that OST g and OST e+g do not change the location of the maximum entries in the estimates of H returned by OST and OST e , respectively, but only their amplitude. As such, they lead to the same F-measures than OST and OST e , and we did not include them in the experiments of this section.

We first illustrate the complexity of real-data spectra in Fig. 4, where the amplitudes of the first six partials (the components corresponding to the harmonic frequencies) of a single piano note are represented along time. Depending on the partial order q, the amplitude evolves with asynchronous beats and with various slopes. This behaviour is characteristic of piano sounds in which each note comes from the vibration of up to three coupled strings. As a consequence, the spectral envelope of such notes cannot be well modelled by a fixed amplitude pattern. Fig. 4 shows that, thanks to its flexibility, OST e can perfectly recover the true fundamental frequency (MIDI 50) while PLCA (a) Thresholded OSTe transcription is prone to octave errors (confusions between MIDI 50 and MIDI 62). Then, Table 1 reports the F-measures returned by the three competing approaches on seven 15-s extracts of pieces from Chopin, Beethoven, Mussorgski and Mozart. For each of the three methods, we have also included a variant that incorporates a flat component in the dictionary that can account for noise or non-harmonic components. In PLCA, this merely consists in adding a constant vector w f (K+1) = 1/M to W. In OST or OST e this consists in adding a constant column to C, whose amplitude has also been validated over 3 orders of magnitude. OST performs comparably or slightly inferiorly to PLCA but with an impressive gain in computational time (∼3000× speedup). Best overall performance is obtained with OST e +noise with an average ∼10% performance gain over PLCA and ∼750× speedup.

A Python implementation of OST and real-time demonstrator are available at https://github. com/rflamary/OST

Conclusions

In this paper we have introduced a new paradigm for spectral dictionary-based music transcription.

As compared to state-of-the-art approaches, we have proposed a holistic measure of fit which is robust to local and harmonically-related displacements of frequency energies. It is based on a new form of transportation cost matrix that takes into account the inherent harmonic structure of musical signals. The proposed transportation cost matrix allows in turn to use a simplistic dictionary composed of Dirac vectors placed at the target fundamental frequencies, eliminating the problem of choosing a meaningful dictionary. Experimental results have shown the robustness and accuracy of the proposed approach, which strikingly does not come at the price of computational efficiency. Instead, the particular structure of the dictionary allows for a simple algorithm that is way faster than state-of-the-art NMF-like approaches. The proposed approach offers new foundations, with promising results and room for improvement. In particular, we believe exciting avenues of research concern the learning of C h from examples and extensions to other areas such as in remote sensing, using application-specific forms of C.

Figure 1 :

 1 Figure 1: Comparison of transportation cost matrices C 2 and C h (full matrices and selected columns).

Figure 3 :

 3 Figure 3: Unmixing under model misspecification. See text for details.

Figure 4 :

 4 Figure 4: First 6 partials and transcription of a single piano note (note D3, ν 0 = 147 Hz, MIDI 50).

  Three example spectra v n compared to a given template v (left) and computed divergences (right). The template is a mere Dirac vector placed at a particular frequency ν 0 . D 2 denotes the standard quadratic error xy 2 2 . By construction of D C h , sample v 3 which is harmonically related to the template returns a very good fit with the latter OT divergence. Note that it does not make sense to compare output values of different divergences; only the relative comparison of output values of the same divergence for different input samples is meaningful.

	v1|v)	1.13 72.92	145.00	134.32
	D(v2|v)	1.13	5.42	10.00	10.00
	D(v3|v)	0.91	2.02	1042.67	1.00
	Figure 2:

Table 1 :

 1 Recognition performance (F-measure values) and average computational unmixing times.

	MAPS dataset file IDs	PLCA PLCA+noise	OST	OST+noise OSTe OSTe+noise
	chpn_op25_e4_ENSTDkAm	0.679	0.671	0.566	0.564	0.695	0.695
	mond_2_SptkBGAm	0.616	0.713	0.470	0.534	0.610	0.607
	mond_2_SptkBGCl	0.645	0.687	0.583	0.676	0.695	0.730
	muss_1_ENSTDkAm 4	0.613	0.478	0.513	0.550	0.671	0.667
	muss_2_AkPnCGdD	0.587	0.574	0.531	0.611	0.667	0.675
	mz_311_1_ENSTDkCl	0.561	0.593	0.580	0.628	0.625	0.665
	mz_311_1_StbgTGd2	0.663	0.617	0.701	0.718	0.747	0.747
	Average	0.624	0.619	0.563	0.612	0.673	0.684
	Time (s)	14.861	15.420	0.004	0.005	0.210	0.202

29th Conference on Neural Information Processing Systems (NIPS 2016), Barcelona, Spain.

Acknowledgments. This work is supported in part by the European Research Council (ERC) under the European Union's Horizon 2020 research & innovation programme (project FACTORY) and by the French ANR under the JCJC programme (project MAD). Many thanks to Antony Schutz for generating & providing some of the musical data.


Célia Da Costa Pereira

Mauro Dragoni
email: dragoni@fbk.eu

Andrea G B Tettamanzi
email: andrea.tettamanzi@unice.fr

Serena Villata
email: serena.villata@unice.fr

Fuzzy Labeling for Abstract Argumentation: An Empirical Evaluation

Argumentation frameworks have to be evaluated with respect to argumentation semantics to compute the set(s) of accepted arguments. In a previous approach, we proposed a fuzzy labeling algorithm for computing the (fuzzy) set of acceptable arguments, when the sources of the arguments in the argumentation framework are only partially trusted. The convergence of the algorithm was proved, and the convergence speed was estimated to be linear, as it is generally the case with iterative methods. In this paper, we provide an experimental validation of this algorithm with the aim of carrying out an empirical evaluation of its performance on a benchmark of argumentation graphs. Results show the satisfactory performance of our algorithm, even on complex graph structures as those present in our benchmark.

Introduction

In crisp argumentation, arguments are evaluated, following a specific semantics, as being acceptable or not acceptable, as shown by Dung [9]. Roughly, accepted arguments are those arguments which are not attacked by other accepted arguments, and rejected arguments are those attacked by accepted arguments. The set of accepted arguments, called extensions, represent consistent set(s) of arguments that can be accepted together. However, in many applications, such as decision making and agent-based recommendation systems, such a crisp evaluation of the arguments is not suitable to represent the complexity of the considered scenario. To address this issue, we [8] proposed to perform a fuzzy evaluation of the arguments of an argumentation framework. Such a fuzzy evaluation of the arguments is originated by the observation that some arguments may come from only partially trusted sources. To represent the degrees of trust, we rethink the usual crisp argument evaluation [9,6] by evaluating arguments in terms of fuzzy degrees of acceptability. In our previous contribution [8], we proved that the fuzzy labeling algorithm used to assign to the arguments fuzzy degrees of acceptability converges, and we discussed its convergence speed. However, no empirical evaluation was addressed to support this discussion.

In this paper, we face this issue by providing an extensive evaluation of the performance and scalability of the fuzzy labeling algorithm with respect to a benchmark of abstract argumentation frameworks. We first select three existing datasets for abstract argumentation tasks used in the literature, namely the datasets created by Bistarelli et al. [START_REF] Bistarelli | Benchmarking hard problems in random abstract afs: The stable semantics[END_REF], by Cerutti et al. [7], and by Vallati et al. [12]. Moreover, we generate our own dataset of abstract argumentation frameworks by randomly combining some well known graph patterns in argumentation theory into 20,000 bigger argumentation frameworks. Second, we study the behaviour of the algorithm with respect to the frameworks in the benchmark, to check whether its performance are satisfiable even considering huge and complex networks as those represented in the datasets, e.g., presenting an increasing number of strongly connected components.

The reminder of the paper is as follows. In Section 2, we provide the basics of abstract argumentation theory and fuzzy set theory. Section 3 firstly recalls the main concepts behind the definition of the fuzzy labeling algorithm presented in [8], and secondly, describes its current implementation. In Section 4, we describe the four datasets which compose the bechmark used to evaluate our algorithm, and we report about the obtained results. Conclusions end the paper.

Preliminaries

In this section, we provide some insights about abstract argumentation theory and fuzzy sets.

Abstract argumentation theory

We provide the basics of Dung's abstract argumentation [9]. Definition 1. (Abstract argumentation framework) An abstract argumentation framework is a pair A, → where A is a set of elements called arguments and →⊆ A×A is a binary relation called attack. We say that an argument A i attacks an argument A j if and only if (A i , A j ) ∈→.

Dung [9] presents several acceptability semantics which produce zero, one, or several sets of accepted arguments. These semantics are grounded on two main concepts, called conflict-freeness and defence. The concepts of Dung's semantics are originally stated in terms of sets of arguments. It is equal to express these concepts using argument labeling [10, 13, 6] In a reinstatement labeling [6], an argument is labeled "in" if all its attackers are labeled "out" and it is labeled "out" if it has at least an attacker which is labeled "in".

Definition 2. (Conflict-free, Defence) Let C ⊆ A. A set C is conflict-free if and only if there exist no A i , A j ∈ C such that A i → A j . A set C defends an argument A i if and only if for each argument A j ∈ A if A j attacks A i then there exists A k ∈ C such that A k attacks A j .
Definition 4. (AF-labeling [6]) Let A, → be an abstract argumentation framework. An AF-labeling is a total function lab : A → {in, out, undec}. We define in(lab)

= {A i ∈ A|lab(A i ) = in}, out(lab) = {A i ∈ A|lab(A i ) = out}, undec(lab) = {A i ∈ A|lab(A i ) = undec}.
Definition 5. (Reinstatement labeling [6]) Let lab be an AF-labeling.

a in-labeled argument is said to be legally in iff all its attackers are labeled out. a out-labeled argument is said to be legally out iff it has at least one attacker that is labeled in. an undec-labelled argument is said to be legally undec iff not all its attackers are labelled out and it does not have an attacker that is labelled in. A fuzzy set is completely defined by its membership function. Therefore, it is useful to define a few terms describing various features of this function. Given a fuzzy set A, its core is the (conventional) set of all elements x such that A(x) = 1; its support, supp(A), is the set of all x such that A(x) > 0. A fuzzy set is normal if its core is nonempty. The set of all elements x of A such that A(x) ≥ α, for a given α ∈ (0, 1], is called the α-cut of A, denoted A α .

The usual set-theoretic operations of union, intersection, and complement can be defined as a generalization of their counterparts on classical sets by introducing two families of operators, called triangular norms and triangular co-norms. In practice, it is usual to employ the min norm for intersection and the max co-norm for union. Given two fuzzy sets A and B, and an element x,

(A ∪ B)(x) = max{A(x), B(x)}; (1)
(A ∩ B)(x) = min{A(x), B(x)}; (2) Ā(x) = 1 -A(x). (3)
Finally, given two fuzzy sets A and B, A ⊆ B if and only if, for every element x, A(x) ≤ B(x).

Fuzzy Labeling for Abstract Argumentation

In this section, we recall the fuzzy labeling algorithm for abstract argumentation [8] we want to empirically evaluate over the available datasets for abstract argumentation to study its performance. For a complete description of the algorithm and its convergence theorem as well as the comparison with the related approaches we remind the reader to [8]. Moreover, we report about the implementation we develop to test the perfomances of the algorithm.

Algorithm

In order to account for the fact that arguments may originate from sources that are trusted only to a certain degree, the (crisp) abstract argumentation structure described in Section 2 may be extended by allowing gradual membership of arguments in the set of arguments A. We have that A is a fuzzy set of trustworthy arguments, and A(A), the membership degree of argument A in A, is given by the trust degree of the most reliable (i.e., trusted) source that offers argument

A, 3 A(A) = max s∈src(A) τ s , (4)
where src(A) is the set of sources proposing argument A and τ s is the degree to which source s ∈ src(A) is trusted. We do not make any further assumptions on the trust model, as it is out of the scope of this paper. However, we refer the interested reader to [11], where a more detailed description of how the source trustworthiness degree can be computed starting from elements, like the source sincerity and expertise, is provided.

Definition 7. (Fuzzy AF-labeling) Let A, → be an abstract argumentation framework. A fuzzy AF-labeling is a total function

α : A → [0, 1].
Such an α may also be regarded as (the membership function of) the fuzzy set of acceptable arguments: α(A) = 0 means the argument is outright unacceptable, α(A) = 1 means the argument is fully acceptable, and all cases inbetween are provided for.

Definition 8. (Fuzzy Reinstatement Labeling) Let α be a fuzzy AF-labeling. We say that α is a fuzzy reinstatement labeling iff, for all arguments A,

α(A) = min{A(A), 1 -max B:B→A α(B)}.
(5)

The above definition combines two intuitive postulates of fuzzy labeling: (1) the acceptability of an argument should not be greater than the degree to which the arguments attacking it are unacceptable and (2) an argument cannot be more acceptable than the degree to which its sources are trusted: α(A) ≤ A(A).

We can verify that the fuzzy reinstatement labeling is a generalization of the crisp reinstatement labeling of Definition 5, whose in and out labels are particular cases corresponding, respectively, to α(A) = 1 and α(A) = 0. The intermediate cases, 0 < α(A) < 1 correspond to a continuum of degrees of "undecidedness", of which 0.5 is but the most undecided representative.

We denote by α 0 = A the initial labeling, and by α t the labeling obtained after the t th iteration of the labeling algorithm. Definition 9. Let α t be a fuzzy labeling. An iteration in α t is carried out by computing a new labeling α t+1 for all arguments A as follows:

α t+1 (A) = 1 2 α t (A) + 1 2 min{A(A), 1 -max B:B→A α t (B)}. ( 6
)
This defines a sequence {α t } t=0,1,... of labelings which always converges to a limit fuzzy labeling, as proven in [8]. Moreover, the convergence speed is linear: in practice, a small number of iterations is enough to compute the limit up to the desired precision. The fuzzy labeling of a fuzzy argumentation framework is thus the limit of {α t } t=0,1,... . Definition 10. Let A, → be a fuzzy argumentation framework. A fuzzy reinstatement labeling for such argumentation framework is, for all arguments A, α(A) = lim t→∞ α t (A).

(7)

Implementation

The fuzzy labeling algorithm has been implemented by using the Java language without using any specific library. The first version of the algorithm was developed with the support of multi-threading where the update of each node was parallelized. However, preliminary tests run in multi-threading mode reported some computational overhead that make the observation of algorithm performance not consistent. For this reason and for measuring performance values easing the comparison with the obtained results, we opted for running all tests in a single-thread mode.

Each run was performed on a server equipped with a Xeon E5-2609 v2 @ 2.50 Ghz. According to official user-based benchmarks,4 the single thread mark of the CPU is 1229 points. This value can be used as reference for normalizing results concerning the timing of each run obtained on other machines.

Evaluation

In this section, we study the behaviour and the performances of the fuzzylabeling algorithm over a benchmark for abstract argumentation, and then we report about the obtained results.

The aim of our experimental analysis is to assess the scalability of the fuzzylabeling algorithm concerning two perspectives:

the number of iterations needed for convergence with respect to the number of the nodes in the graph, and the time needed for convergence with respect to the number of the nodes in the graph.

It must be stressed that the time needed for convergence depends on (i) the time needed for computing each iteration, (ii) the time needed to update the α of each single argument, and (iii) the number of iterations required for the labeling to converge.

Benchmark

The benchmark we used to evaluate the performances of the fuzzy labeling algorithm is composed of different datasets for abstract argumentation tasks used in the literature. More precisely, we have considered the following datasets:

-The Perugia dataset [2, 4, 3]:5 the dataset is composed of randomly generated directed-graphs. To generate random graphs, they adopted two different libraries. The first one is the Java Universal Network/Graph Framework (JUNG), a Java software library for the modeling, generation, analysis and visualization of graphs. The second library they used is NetworkX, a Python software package for the creation, manipulation, and study of the structure, dynamics, and functions of complex networks. Three kinds of networks are generated:

• In the Erdős-Rényi graph model, the graph is constructed by randomly connecting n nodes. Each edge is included in the graph with probability p independent from every other edge. For the generation of these argumentation graphs, they adopted p = c log n/n (with c empirically set to 2.5), which ensures the connectedness of such graphs.

• The Kleinberg graph model adds a number of directed long-range random links to an n × n lattice network, where vertices are the nodes of a grid with undirected edges between any two adjacent nodes. Links have a non-uniform distribution that favors edges to close nodes over more distant ones. • In the Barabási-Albert graph model, at each time step, a new vertex is created and connected to existing vertices according to the principle of "preferential attachment", such that vertices with higher degree have a higher probability of being selected for attachment. For more details about the generation of these networks as well as the graph models, we refer the reader to [START_REF] Bistarelli | A first comparison of abstract argumentation systems: A computational perspective[END_REF]4,[START_REF] Bistarelli | Benchmarking hard problems in random abstract afs: The stable semantics[END_REF].

-The dataset used by Cerutti et al. in their KR 2014 paper [7] (which we will call the KR dataset): the dataset has been generated to evaluate a metaalgorithm for the computation of preferred labelings, based on the general recursive schema for argumentation semantics called SCC-Recursiveness. The dataset is composed of three sets of argumentation frameworks, namely:

• 790 randomly generated argumentation frameworks where the number of strongly connected components (SCC) is 1, varying the number of arguments between 25 and 250 with a step of 25. • 720 randomly generated argumentation frameworks where the number of strongly connected components varies between 5 and 45 with a step of 5. The size of the SCCs is determined by normal distributions with means between 20 and 40 with a step of 5, and with a fixed standard deviation of 5. They similarly varied the probability of having attacks between arguments among SCCs. • 2800 randomly generated argumentation frameworks where the number of strongly connected components is between 50 and 80 with a step of 5. -The dataset presented by Vallati et al. at ECAI 2014 [12] (which we will call the ECAI dataset): the dataset was produced to study the features of argumentation frameworks. More precisely, it is composed of 10, 000 argumentation frameworks generated using a parametric random approach allowing to select (probabilistically average, standard deviation) the density of attacks for each strongly connected component, and how many arguments (probabilistically) in each SCC attack how many arguments (probabilistically) in how many (probabilistically) other SCCs. The number of arguments ranges between 10 and 40, 000, and they exploited a 10-fold cross-validation approach on a uniform random permutation of the instances.

The availability of real-world benchmarks for argumentation problems is quite limited, with some few exceptions like [5] or AIFdb.6 However, these bench- marks are tailored towards problems of argument mining and their representation as abstract argumentation frameworks usually leads to topologically simple graphs, such as cycle-free graphs. These kinds of graphs are not suitable for evaluating the computational performance of solvers for abstract argumentation problems. For this reason, we decided to use artificially generated graphs as benchmarks, in line with the preliminary performance evaluation of Bistarelli et al. [START_REF] Bistarelli | A first comparison of abstract argumentation systems: A computational perspective[END_REF].

a → b, b → a a → b, b → a, b → c a → b, b → a, b → c, c → d, d → c a → b, b → c, c → d, d → e, e → f , f → e a → b, b → c, c → d, d → c a → b, b → c, c → a a → b, b → a, b → c, c → c a → b, b → a, b → c, c → c, d → d a → b, b → a, b → c, a → c, c → d, d → c a → b, b → a, b → c, a → c, c → d a → b, b → c, c → a, b → d, a → d, c → d, d → e, e → d a → b, b → c, e → c, c → d a → b, b → a, b → c, c → d, d → e, e → c a → b, b → c, c → c a → b, b → c, c → a, a → d, b → d, c → d a → b, a → c, c → a, c → d, d → c, d → a, a → d, c → e, d → f
In order to ensure the consideration of all kinds of interesting "patterns" that could appear in argumentation frameworks (e.g., the abstract argumentation frameworks used to exemplify the behaviour of the semantics in [START_REF] Baroni | An introduction to argumentation semantics[END_REF]), we have generated further graphs by composing these basic well-known examples of interesting argumentation patterns (shown in Figure 1) into bigger frameworks.

Our generated dataset (which we will call the Sophia Antipolis dataset)7 consists of 20,000 argumentation graphs created through a random aggregation of the patterns shown above. This process has been executed with different settings in order to obtain complex graphs of specific sizes. In particular, a set of 1,000 argumentation graphs is generated for graph sizes from 5,000 to 100,000 nodes, with incremental steps of 5,000 nodes each. The aggregation of patterns has been done incrementally, and the connections (edges) between single patterns were generated randomly. The number of created graphs and their different sizes should support the evaluation of argumentation reasoning algorithms under a broad number of scenarios.

Results

Figures 2-11 summarize the behavior of the fuzzy-labeling algorithm on the four datasets we considered. For each dataset, we applied the algorithm to the argumentation graphs with all argument weights set to 1 (i.e., arguments coming from fully trusted sources) and with random weights (i.e., arguments coming from a variety of more or less trusted sources as it may be the case in application scenarios like multiagent systems). From a first inspection of the figures, it is clear that certain graph types are harder than others: the Sophia Antipolis appears to be the hardest, followed by the Erdős-Rényi, the Barabási-Albert, and the KR + EKAI datasets. The Kleinberg dataset appears to be the easiest. Furthermore, for all datasets, the graphs with random weights never require a smaller number of iterations for convergence than their counterparts with all weights fixed to 1.

Figures 2 and3 show the behaviour of the fuzzy labeling algorithm when applied to the Barabási-Albert dataset. In particular, Figure 2 (left-hand side) illustrates the evolution of the number of iterations needed to reach convergence when all the weights are equal to 1. We can notice that the curve follows a logarithmic rise with the increasing of the number of nodes. The figure illustrated through the (right-hand side) curve represents the evolution of the time needed to reach the convergence. It shows a behaviour rather linear. However, we can notice that the slope of the curve decreases with the increasing of the number of nodes. A similar behaviour is depicted in Figure 3 which illustrates the evolution of the quantity of time (in ms) needed to reach the convergence when the weights are assigned randomly. These two illustrations clearly show the capability of the fuzzy labeling algorithm to handle a growing amount of data.

In Figures 4 and5, we present the behaviour of the fuzzy labeling algorithm when applied to the Erdős-Rényi dataset. We can notice that when all the weights are equal to 1, the convergence is reached very quickly both when considering the number of iterations, and the quantity of time needed for convergence. However, while such a quantity is quite similar with respect to the case in which the weights are randomly assigned, we can notice that the number of iterations needed for convergence is higher with respect to the behaviour illustrated in Figure 4. This can be due to the fact that the Erdős-Rényi dataset is constructed by randomly connecting the nodes. As we can see in Figures 6 and7, the convergence with the Kleinberg dataset is even globally faster, either when all weights are equal to 1 or when the weights are randomly assigned. Instead, the behaviour on the Sophia Antipolis dataset, shown in Figures 10 and11, is quite similar to the behaviors obtained with the Barabási-Albert dataset.

It is less evident, but the fuzzy-labeling algorithm behaves on the KR + ECAI dataset (illustrated in Figures 8 and9) much like it does on the Barabási-Albert and Sophia Antipolis datasets, with the exception of a few small graphs which are outliers and which demand a relatively large number of iterations to converge. Nevertheless, the time behavior of Barabási-Albert, Sophia Antipolis and KR + ECAI is qualitatively identical. Despite the differences among the various graph types, we have a rate of increase in time which is at most log-linear for all graph types and for all weight assignments.

Conclusions

We have evaluated the performance of the fuzzy labeling algorithm proposed in [8] on a benchmark consisting of four datasets of argumentation graphs having widely different characteristics. The experimental results clearly indicate that the fuzzy labeling algorithm scales up nicely in all circumstances, and is thus a viable argumentation reasoning tool.

Definition 3 .

 3 (Acceptability semantics) Let C be a conflict-free set of arguments, and let D : 2 A → 2 A be a function such that D(C) = {A|C def ends A}. -C is admissible if and only if C ⊆ D(C). -C is a complete extension if and only if C = D(C).

Fig. 1 .

 1 Fig. 1. The "patterns" used for constructing the Sophia Antipolis dataset.

Fig. 2 .Fig. 3 .

 23 Fig. 2. Barabási-Albert dataset of the Perugia benchmark with all weights equal to 1.0.

Fig. 4 .Fig. 5 .Fig. 6 .Fig. 7 .Fig. 8 .Fig. 9 .Fig. 10 .

 45678910 Fig. 4. The Erdős-Rényi dataset of the Perugia benchmark with all weights equal to 1.0.

Fig. 11 .

 11 Fig. 11. The Sophia Antipolis benchmark with random weights.

  -C is a grounded extension if and only if it is the smallest (w.r.t. set inclusion) complete extension. -C is a preferred extension if and only if it is a maximal (w.r.t. set inclusion) complete extension. -C is a stable extension if and only if it is a preferred extension that attacks all arguments in A \ C.

Here, we suppose that the agent is optimistic. To represent a pessimistic behaviour, we should use the min operator, for example.

https://www.cpubenchmark.net/CPU_mega_page.html

The dataset is available at http://www.dmi.unipg.it/conarg/dwl/networks.tgz.

http://corpora.aifdb.org

The Sophia Antipolis dataset is available at https://goo.gl/pN1M9r.


Kartheek Medathati
email: kartheek.medathati@inria.fr

Andrew I Meso

Guillaume S Masson

Pierre Kornprobst

James Rankin

Understanding the impact of recurrent interactions on population tuning Application to MT cells characterization

A lot of theoretical work has been done in understanding neural computation and dynamics regulated by excitatory and inhibitory interactions [5,[START_REF] Raijmakers | [END_REF] but:

• Most studies are conducted without structured input.

• Their impact in terms of sensory information processing is often missing. The goal of this work is to bridge this gap by focusing on the impact of lateral interactions on population tuning. I Methodology

• A ring network model under neural fields formalism with a structured input.

• Bifurcation analysis to understand model's behavior depending on connectivity and input parameters. I Application to MT cells characterization Neurons in the middle temporal (MT) visual area have been linked to 2D motion perception [START_REF] Movshon | Pattern recognition mechanisms[END_REF]1]. Direction-selective cells encode either (local) component or (global) pattern motion. However, these properties depend on the spatiotemporal inputs (e.g. plaids vs. Random Dot Patterns (RDPs)) as well as on the recurrent local interactions. We aim to better understand how the local interactions in conjunction with driving stimuli lead to di erent tuning behaviors. • On the left, connectivity kernels with narrow to broad excitatory extent.

Model description

I
• On the right, connectivity kernels with low to high lateral inhibition strength.

I Solution space

• Stable and unstable solutions the network could converge to under the presence of driving input. • There are parameter regimes where multiple solutions co-exist.

• The stability of the solution also changes quite sharply w.r.t to parameters of the driving input. • All theoretically stable solutions are not reached unless there is a bias in the initial conditions.

Conclusion

Our model demonstrates that recurrent interactions can shape direction tuning of MT neurons. The properties of both excitation width and inhibition strength explain the di erent tuning classes and their dynamics with respect to spatiotemporal properties of the input. Moreover, prototypical tuning shapes (e.g. VA, SB and TP) correspond to di erent regimes of a single dynamical system, where sharp transitions can be identified.

I Experimental Observation MT cells show di erent tuning behaviors [6] when stimulated with two overlapped motion components such as RDPs or Plaids. u(0, t) u(0, t) u(0, t) u(0, t)

• Observation: Recurrently interacting homogeneous population of Direction Selective (DS) cells can reproduce di erent kinds of tuning observed experimentally.

• Fluctuations can occur at single cell level due to the multistability. I Input driven changes in tuning

• Model behavior changes from VA to SB or consistent SB to fluctuating component selection. These behavior changes are driven by changes in component separation and relative strength.

• Assymetries in the input could significantly stabilize MT single unit tuning from a fluctuating WTA to SB. Temporal dynamics I Experimental Observation Temporal dynamics of MT cells show tuning transitions [6] and develop their selectivity gradually. • Slow onset of inhibition could lead to a lack of tuning behavior in the initial time period, lateral competition then leads to observed tuning behaviors.

  Ring model with Excitatory/Inhibitory interactions u(✓, t) (✓ 2 [ ⇡, ⇡)) denotes a population of directionally tuned MT neurons with dynamics governed by: µu( , t))d + Iext(✓), where, J is the connectivity kernel, S is a sigmoid function (µ regulates the sigmoidal gain), Iext is the driving input. I Definition of Iext Iext is defined as a linear combination of gaussian bumps. It allows us to represent di erent kinds of stimuli used to probe motion estimation. • Peak width (PW) : Representative of uncertainity in local velocity estimates (standard dev. of the bump) • Peak separation (PS) : Directional difference between component motions. I Definition of Jg e, e,gi, i (✓) Let Jg e, e,gi, i (✓) denote weighted di erence of gaussians: Jg e, e,gi, i (✓) = geG(✓, e) g i G(✓, i ), where, G(✓, ): Gaussian function, ge: excitatory strength, e: extent of excitatory surround, g i : inhibitory strength, i : extent of inhibitory surround. Numerical study of the model I Exploration of the connectivity space with J↵, J↵, (✓) = J ge ↵ , e↵,gi↵+ , i (✓) • Assumption: We consider the case of a uniform lateral inhibition, so i is fixed to a large value ( i = 10⇡). • ↵ is a parameter to smoothly vary the extent of excitatory surround from narrow to broad while mainting constraints imposed on specific Fourier components. • is a free parameter to regulate the strength of lateral inhibition. Numerical study of the model (Cont.) I Connectivity kernels: J ↵,

I

  Likelihood of convergenceWe estimate the strength of the stable solutions using random noise as initial condition w.r.t.input and kernel parameters (100 trials)

1 •

 1 Vector average (VA) Side biased (SB) Side biased(II) Two peaked (TP) I Question Could recurrent local interactions in the feature space (direction of motion) lead to these observed tuning behaviors? In terms of modelling, can a single model reproduce these behaviours by considering local recurrent interactions in the feature space? I Model behavior: Population versus Single unit tuning

  component tuning classification has been traditionally obtained with plaid patterns [3] but several empirical studies have shown that such classfication cannot translate from one stimulus class to another [2, 6], for example from plaid stimuli to RDPs. Pattern Cells Component Cells I Question 2 Can our model capture the changes in the tuning behavior with respect to changes in the input and therefore explain the di culty in predicting the tuning behavior across stimuli type? I Model behavior: Connectivities w.r.t to Input stimuli • Extent of overlap between likelihood of di erent tunings vary across stimuli and certain transitions have a very weak likelihood.

  explained by temporal dynamics of inhibition? I Model behavior: slow onset of inhibition Time (ms)

This work was partially supported by the EC IP project FP7-ICT-2011-8 no. 318723 (MatheMACS)