Sinkhorn Regression

1 Introduction

Regression analysis is an important statistical technique frequently applied in machine learning for a large variety of tasks such as image classification [Naseem et al., 2010] and subspace segmentation [Lu et al., 2012]. However, traditional linear regression models (e.g., least square regression) are mainly based on the normal population, which cannot provide the robustness against outlier. To address this problem, robust regression has emerged as a powerful tool to handle the data with non-Gaussian noises. The representative models include Robust Sparse Coding (RSC) [Yang et al., 2011]

∗Corresponding Author. H. Huang was partially supported by U.S. NSF IIS 1836945, IIS 1836938, IIS 1845666, IIS 1852606, IIS 1838627, IIS 1837956.

and Correntropy based Sparse Representation (CESR) [He et al., 2011]. They can be viewed as M-estimation problems induced by different influence functions which can make outlying data have a relatively smaller influence on the regression estimator.

However, those robust regression methods ignore the spatial structure information of real noises since they usually assume that noise pixels are independently generated. To alleviate this problem, [Li et al., 2013] proposed the Structured Sparse Error Coding (SSEC) model by exploring the intrinsic structure of continuous occlusion. [Jia et al., 2012] used a class of structured sparsity inducing norms to fit structural noises. [Yang et al., 2017] presented the two-dimensional image-matrix-based error models by employing nuclear norm to characterize the practical noise matrix. Although these methods have shown great potential in handling structural noises such as occlusion and illumination, they cannot characterize the geometric knowledge that often exists in real data and is helpful for improving the model performance.

More recently, Wasserstein distance, derived from the Optimal Transport (OT) theory, has drawn ample attention in many machine learning tasks. Differing from Lp distances (p ≥ 1) or Kullback-Leibler and other f-divergences, which require distributions to be absolutely continuous with respect to each other or to a base measure, Wasserstein distance is well-defined between any pair of probability distributions over a sample space equipped with a metric. Thus, it provides a meaningful notion of closeness (i.e., distance) for distributions supported on non-overlapping low dimensional manifolds. Due to this advantage, Wasserstein distance has been successfully applied to cancer detection [Ozolek et al., 2014] and super-resolution [Kolouri and Rohde, 2015] problems. In addition, Kusner et al. [Kusner et al., 2015] proposed the Word Mover's Distance (WMD), an implementation of the Wasserstein distance for textual data. To capture both global (at distribution scale) and local (at samples' scale) interactions between classes, Flamary et al. [Flamary et al., 2016] described a new Wasserstein Discriminant Analysis by using a mechanism that can induce neighborhood preservation.

Motivated by the success of Wasserstein distances, in this paper, we propose a novel efficient and robust Matrix Regression method by employing joint Sinkhorn distances (i.e., Wasserstein distances with Entropic constraints) minimization on both loss function and regularization. This is the first time for exploiting Sinkhorn distances to characterize both error and coefficient matrix in the regression model. To provide the high flexibility for non-distribution features (e.g., original data features), we relax the proposed new model into its unbalanced formulation by virtue of Kullback-Leibler divergence. According to the OT geometry [Janati et al., 2018], we know that using Sinkhorn distance as regularizer can promote parameters that are close since it takes into account a prior geometric knowledge on the regressor variables. As a result, our method is more reliable and applicable than traditional regression method using Lp-norm regularizer. To solve SMR, we derive an efficient algorithm based on Sinkhorn iteration, which iterates through applications of alternating optimization. Moreover, as the theoretical contribution of this paper, we provide a statistical bound for the proposed new model by leveraging the Rademacher complexity. We apply the proposed algorithm to feature selection problem, and make comparisons with some existing approaches. The experimental results show that the proposed method is more effective than the state-of-the-art methods.

2 Preliminaries

We summarize the notations and definitions used in this paper. Matrices are written as boldface uppercase letters. Vectors are written as boldface lowercase letters. ‖ · ‖F and ‖ · ‖∗ denote Frobenius norm and nuclear norm [Liu et al., 2010], respectively. 〈·, ·〉 is the inner product operation. eNS ∈ RNS is a column vector of ones. R+ denotes the set of all non-negative real number. RNS×NG + denotes the set of all positive semi-define matrices on RNS×NG . For matrix M ∈ Rn×m, its i-th row, j-th column are denoted by mi , mj, respectively. The L2,1-norm of M is defined as ‖ M ‖2,1= ∑ni=1 √∑mj=1 m2ij = ∑ni=1 ‖ mi ‖2, where ‖ · ‖2 denotes the L2-norm. We define the Kullback-Leibler (KL) divergence [Janati et al., 2018] between two positive vectors x and y by

The KL divergence between two matrices A ∈ RNS×NS + (i.e., aij ≥ 0 for all (i, j)) and B ∈ RNS×NS ++ (i.e., bij > 0 for all (i, j)) is defined as

Optimal transport: OT theory, originally used to study the problem of resource allocation, provides a powerful geometrical tool for comparing probability distributions.

In a more formal way, given access to two sets of points XS = {xSi ∈ Rd }NS i=1 and XG = {xGi ∈ Rd }NG i=1 , we construct two discrete empirical probability distributions as:

where pSi and pGi are probabilities associated to xSi and xGi respectively and δx is a Dirac measure that can be interpreted as an indicator function taking value 1 at the position of x and 0 elsewhere. For these two distributions, the optimal transport (or Monge-Kantorovich) problem consists in finding a probabilistic coupling defined as a joint probability measure over XS × XG with marginals ˆµS and ˆµG that minimizes the cost of transport with respect to some ground metric D : Xs ×Xg →R+ :

where UˆµS,ˆµG = {P ∈ RNS×NG + : PeNG = pS , PT eNS = pG} is a set of doubly stochastic matrices, pS = [pS1 , pS2 , · · · , pSNS ]T and pG = [pG1 , pG2 , · · · , pGNG ]T . C ∈ RNS×NG + is a given cost matrix, i.e., ci,j = D(xSi , xGj ), defining the energy needed to move a probability mass from xSi to xGj . This problem admits a unique solution P∗ and defines a metric on the space of probability measures (called the Wasserstein distance) as

Entropic regularization: Solving problems with Wasserstein distance fitting errors can require solving several costly optimal transport problems. As a minimum of affine functions, the Wasserstein distance itself is not a smooth function. To avoid both of these issues, [Cuturi, 2013] proposed to smooth the optimal transport problem with an entropic term:

where γ > 0 and e is the (strictly concave) entropy function:

(6) allows to solve the optimal transportation problem efficiently using Sinkhorn-Knopp matrix scaling algorithm [Cuturi, 2013]. Thus, (6) is also called Sinkhorn distance, which preserve all advantages of Wasserstein distance. Compared to the original case (5), it can obtain smoother and more numerically stable solutions.

Matrix regression: In data mining and machine learning, a common paradigm for matrix regression is to minimize the penalized empirical loss:

where λ > 0 is the balance parameter, Z ∈ Rm×n is the parameter to be estimated from the training sample matrix A ∈ Rm×l and the response matrix Y ∈ Rl×n , L(Y-AZ) is the empirical loss on the training set, and Ω(Z) is the regularization term that encodes feature relatedness. Different assumptions on the loss Y - AZ and variate Z can lead to different models. For example, if both L(·) and Ω(·) are L2,1-norm, (6) becomes the feature selection model based on L2,1-norms [Nie et al., 2010]; if L(·) is L1-norm [Yang et al., 2012; Yang et al., 2013] and Ω(·) is nuclear norm, (6) leads to the Low Rank Representation (LRR) model [Liu et al., 2010]; if both L(·) and Ω(·) are the square of Frobenius norm, (6) is the well-known Least Squares Regression (LSR) model. These models have been widely applied to many tasks such as feature selection, subspace segmentation and image classification.

3 Matrix Regression Based on Joint Sinkhorn Distances

In the above formulations, the loss term and estimated variate are characterized via the simple matrix norm. Thus, these models can be easily solved by conventional convex optimization methods (e.g., ADMM [Liu et al., 2010] and reweighted iterative methods [Nie et al., 2010]). However, they suffer from two limitations. First, these matrix norm cannot offer the exibility in adapting them to various data types since they are nonparametric. Secondly, they do not taking into account the geometry of the data through the pairwise distances between the distributions' points. Accordingly, these models often achieve the suboptimal results in practical applications, especially when real data is corrupted by noise.

Sinkhorn Regression. Comparing with matrix norm, Sinkhorn distance can circumvent the above limitations (as stated in the introduction). Therefore, we propose to use Sinkhorn distance to jointly characterize loss term and estimated variate Z, which is formulated as

where h(·) denotes the histogram operator and the latent variable d consists in estimating the barycenter of {h(z1T ), h(z2T ), · · · , h(zmT )}.

In model (9), the histogram operator is used to constrain ZT ai, YT and ZT , respectively. Thus, we do not need to extract the histogram features of data in the preprocessing stage. Meanwhile, since histogram features are invariant in movement, rotation and zoom, our method is more robust to the real noisy data than the other methods using matrix norm over original data.

Remark 1. Sinkhorn distance is used to compare row vectors of AT Z and Y in model (9). In fact, Sinkhorn distance can also be imposed on the corresponding column vectors, which induces the following model:

The algorithm for solving the above model is similar to Algorithm 1. Due to the space limitations, we only focus on model (9) in this paper.

The proposed algorithm. Solving problem (9) is extremely challenging since it not only contains the composition of h(·) and W(·, ·), but also the computations of Sinkhorn distances with regard to different terms. Some existing algorithms [Rolet et al., 2016; Sommerfeld et al., 2018] are only suitable for solving Wasserstein loss minimization with matrix norm regularizer (e.g., L1- and L2,1-norm). To cope with this challenge, we relax the marginal constraints UˆµS,ˆµG in (6) using a Kullback-Leibler divergence from the matrix to target marginals ˆµS and ˆµG is [Frogner et al., 2015], Algorithm 1 Solving (13) via Alternating Optimization

Input: data matrix A and Y, parameters λ, γ, µ and ρ

Initialization: P0 and ˆP0

Output: model parameter Z

Repeat

for i = 1 to m do

Update each zi with proximal coordinate descent. end for

for i = 1 to l do

Update each ui with proximal coordinate descent.

end for

Update d with proximal coordinate descent.

Update P1, · · ·, Pl, ˆP1, · · ·, ˆPm via Sinkhorn iteration (Algorithms 2)

until convergence

i.e., (5) is converted as

where K = exp(-C/γ) and C is defined in Eq.(4). Let

and assume A, Y > 0. Then model (9) ultimately becomes:

where each zi > 0, i = 1, 2, · · · ,m.

To facilitate the design of the algorithm, we bring into an auxiliary variable U = AT Z and rewrite model (12) as:

where ρ > 0.

The block coordinate descent method can be used to solve (13). This method alternates the minimization with respect to the variables {P1, · · ·, Pl, ˆP1, · · ·, ˆPm}, Z, U and d. It is easy see that these variables are independent. Thus, they can be updated in parallel. The main iterative procedure for solving model (13) is summarized in Algorithm 1. In the following, we give the details for optimizing each variable.

Updating coefficient matrix Z. For each zi , fixing other variables, problem (13) is simplified as:

Considering the definition of KL divergence, (14) becomes

Algorithm 2 Sinkhorn iteration for optimizing P1, · · · , Pl

Input: data matrix A, coefficient matrix Z and auxiliary variable U

Output: P1, · · · ,Pl

for i = 1 to l do

Ki = exp(-Ci/γ), where Ci is the ground metric matrix between uiT and yiT .

repeat

until convergence

Pi ← (p(i)jt)n×n, where the (j, t)-th element of Pi is p(i)jt = wi(j)k(i)jtvi(t) (wi(j) is the j-th element of wi). end for

To solve problem (15), we introduce the following Lemma: Lemma 1. [Janati et al., 2018] Let a, b, α ∈ R+. The function η : x → (x - alog(x)) + bx is convex on R+ and its proximal operator can be obtained in closed form, i.e.,

Let

By Lemma 1, we can easily obtain the proximal operator proxς(·) . Therefore, using proximal coordinate descent, each ziT can be updated by

where δ > 0 is the step size and (AAT )(i) is the i-column of AAT .

Updating auxiliary variable U. For each iteration, we obtain the optimal ui by solving

which can be further written as

Using the similar approach as in (15-18), uTi is updated by

Updating auxiliary variable d. We can use the similar method as in (18) or (19) to update d. Here we omit it.

Updating parameter set {P1, ···, Pl, ˆP1, ···, ˆPm}. For fixed Z and U, the update of each Pi or ˆPi boils down to an OT problem, which can be solved via Sinkhorn iteration [Cuturi, 2013]. These steps are summarized in Algorithm 2, where we list the detailed iteration process for each Pi. For each ˆPi, we can use the similar method to optimize it. The detailed iterative process is listed in Algorithm 3.

Convergence analysis. As pointed out by [Sandler and Lindenbaum, 2011], the alternate optimization process (Algorithm 1) generates a sequence of lower bounded nonincreasing values for the objective of Problem (12), so the

Algorithm 3 Sinkhorn iteration for optimizing ˆP1, · · · , ˆPm

Input: data matrix A and coefficient matrix Z

Output: ˆP1, · · · , ˆPm

for i = 1 to m do

Ki = exp(-Ci/γ), where Ci is the ground metric matrix between ziT and εen T , ε > 0.

repeat

until convergence

ˆPi ← (ˆ p(i)jt)n×n, where the (j, t)-th element of ˆPi is ˆ p(i)jt = wi(j)k(i)jtvi(t) (wi(j) is the j-th element of wi). end for

sequence of objectives converges. Meanwhile, according to [Janati et al., 2018], we know that every accumulation point of the sequences of iterates of Z, U and Pi (or ˆPi) is a generalized fixed point. For the convergence of Algorithm 2, it can be guaranteed by [Cuturi, 2013]. Here we omit it.

4 Theoretical Guarantee for Sinkhorn Regression

Tremendous studies have already been done on the statistical properties of Wasserstein Distances. For example, [Arras et al., 2017] provided the bound to estimate the 2Wasserstein distance between random variables which can be represented as linear combinations of independent random variables. [Frogner et al., 2015] showed that Wasserstein loss can encourage smoothness of the predictions with respect to a chosen metric on the output space. However, theoretical Guarantee of the Sinkhorn distance loss minimization is still deficient. To bridge this gap, we will establish risk bounds of JWMR in this section. For the convenience, as in (8), let

For a sequence S = {(Z1, y1 ), · · · , (ZN, yN )} of i.i.d. training samples, we denote the empirical error and its expectation of L(Z) as

where E(·) is the expectation operator.

To derive the bound of model (9), we need to introduce the following definition and Lemmas.

Definition 1. [Bartlett and Mendelson, 2002] Let G be a family of mapping from Z to R, and S = (z1, · · · , zN) a fixed sample from Z. The empirical Rademacher complexity of G with regard to S is defined as

where σ = (σ1, σ2, · · · , σN), with σ′i s independent uniform random variables taking values in {+1,-1}. σ′i s are called the Rademacher random variables. The Rademacher complexity is defines by taking expectation with respect to the samples S,

We define a space of loss function induced by the hypothesis space Z as

Lemma 2. [Frogner et al., 2015] For any δ > 0, with probability at least 1 - δ, the following holds for all ι ∈ L,

with the constant MC = maxi,jci,j .

Assume that Z∗ is the optimal solution of model (9), we have the following theorems:

Theorem 1. For any δ > 0, with probability at least 1 - δ, it holds that

where the constant MC = maxi,jci,j and RN(H0 ) is the Rademacher complexity measure the complexity of the hypothesis space Ho .

proof. Assume that the loss function Wγ(h(ZT a), h(y)) is preceded with a softmax layer

We know that the function space can be expressed as:

with I a singleton function space with only the identity map. Since ˆRS(I) = 0 (see (25)), it holds

Therefore, connecting Lemma 1, we have

This theorem shows that the Sinkhorn Regression estimator, as well as the empirical risk minimizer asymptotically reaches a √1/2N speed of convergence under very weak hypotheses.

5 Experiments

In Section 3, we provided a general framework for Joint Sinkhorn matrix regression, which can be applied to many tasks such as feature selection, multi-task learning and image classification, etc. In this paper, we evaluate the performance of Sinkhorn regression using the feature selection task, which is important for many real-world applications, such as bioinformatics, text mining, etc. Here, A = [ai, a2, · · ·, an] ∈ Rm×l is a data matrix and Y = [y1 , y2 , · · ·, yn ]T ∈ Rn×c is a label matrix. Our goal is to choose a subset of features for improving prediction accuracy or decreasing the size of the structure without significantly decreasing prediction accuracy of the classifier built using only the selected features.

5.1 Data Descriptions

We use five publicly available microarray data sets and one Mass Spectrometry (MS) datasets: ALLAML data set [Fodor, 1997], the malignant glioma (GLIOMA) data set [Nutt et al., 2003], the human lung carcinomas (LUNG) data set [Bhattacharjee et al., 2001], Human Carcinomas (Carcinomas) data set [Yang et al., 2006], Prostate Cancer gene expression (Prostate-GE) data set [Singh et al., 2002] for microarray data; and Prostate Cancer (Prostate-MS) [Petricoin III et al., 2002] for MS data. To be fair, the Support Vector Machine (SVM) classifier is employed to these data sets, using 5-fold cross-validation for all compared methods.

5.2 ComparedMethods

To validate the effectiveness of Sinkhorn regression for feature selection, we compare it with several classical feature selection methods. The compared algorithms are enumerated as follows.

1. F-statistic(F-s) [Ding and Peng, 2005]: it describes the statistically expected level of (usually) heterozygosity in a population; more specifically the expected degree of a reduction in heterozygosity when compared to Hardy-Weinberg expectation.
2. ReliefF(RF) [Kira and Rendell, 1992]: It relies entirely on statistical analysis and employs few heuristics and is less often foled.
3. mRMR [Sulaiman and Labadin, 2015]: A two-stage feature selection algorithm by combining mRMR and other more sophisticated feature selectors. 4) T-test: it determines if two sets of data are significantly different from each other
4. Information Gain (IG) [Raileanu and Stoffel, 2004]: A formal comparison of the behavior of two of the most popular split functions
5. Robust Feature Selection Based on L2,1-Norms (RFS) [Nie et al., 2010]: The L2,1-norm is emphasizing on both loss function and regularization.
6. Clustering-Guided Sparse Structural Learning (CGSSL) [Li et al., 2014]: An unsupervised feature selection framework by exploiting the cluster analysis and structural analysis with sparsity simultaneously.
7. Regularized Self-Representation (RSR) [Zhu et al., 2015]: An unsupervised features election method by exploiting the self-representation ability of features.

5.3 Classification Accuracy Comparisons

In the experiments, all data sets are standardized to be zeromean and normalized by standard deviation. SVM classifier has been individually performed on all data sets using 5-fold cross-validation. Our feature selection method is compared to several popularly used feature selection methods in bioinformatics, such as F-statistic, reliefF, mRMR, t-test, and information gain. Since the above data sets are for multiclass classification problem, L1-SVM, HHSVM and other methods that were designed for binary classification are not compared Table 1 and 2 show the detailed experimental results using SVM classifier. The average accuracy for each feature selection approach is calculated using the top 20 and top 80 features. It can be seen that our approaches obviously outperform most of methods significantly. Specifically, with top 20 features, our method is around 0%-2% better than other methods all six data sets. With top 80 features, our method achieve an improvement of 0.45% than the second best method (RFS) on average accuracy.

Like many other feature selection algorithms, our proposed Sinkhorn regression also requires several parameters λ, γ and µ to be set in advance. For the results reported in the above subsection, we do not tune the parameter γ and µ and only set them as: γ = 0.01 and µ = 0.1. Better results may be achieved with tuning it. In this subsection, we will discuss sensitivity of parameter λ. Here, we take the top 20 features as an example on all data sets. The detailed results are shown in Table 3. It can be found the best results of our method mainly lie in the interval [0.01 0.5]. But in Table 1 and 2, we choose λ = 0.1.

6 Conclusions

In this paper, we proposed a new robust matrix regression method with emphasizing joint Sinkhorn distances minimization on both loss function and regularization. The Sinkhorn distances based loss function is robust to noise in data points. Meanwhile, Sinkhorn based regularizer can promote parameters that are close. In addition, the proposed Sinkhorn regression is extended to the unbalanced formulation which does not rely on distribution features of data. We provided an efficient algorithm to solve the proposed new model and described its generalization bound from the statistical viewpoint. Our method has been applied into feature selection task. Extensive empirical studies on six standard data sets demonstrated that the proposed algorithm works much more robustly than some existing methods.