Restarting accelerated gradient methods with a rough strong convexity estimate

We propose new restarting strategies for accelerated gradient and accelerated coordinate descent methods. Our main contribution is to show that the restarted method has a geometric rate of convergence for any restarting frequency, and so it allows us to take profit of restarting even when we do not know the strong convexity coefficient. The scheme can be combined with adaptive restarting, leading to the first provable convergence for adaptive restarting schemes with accelerated gradient methods. Finally, we illustrate the properties of the algorithm on a regularized logistic regression problem and on a Lasso problem.

Motivation
The proximal gradient method aims at minimizing composite convex functions of the form where f is differentiable with Lipschitz gradient and ψ may be nonsmooth but has an easily computable proximal operator. For a mild additional computational cost, accelerated gradient methods transform the proximal gradient method, for which the optimality gap F (x k )−F (x * ) decreases as O(1/k), into an algorithm with "optimal" O(1/k 2 ) complexity [9]. Accelerated variants include the dual accelerated proximal gradient [10,12], the accelerated proximal gradient method (APG) [18] and FISTA [1]. Gradient-type methods, also called first-order methods, are often used to solve large-scale problems because of their good scalability and easiness of implementation that facilitates parallel and distributed computations.
In the case when the nonsmooth function ψ is separable, which means that it writes as coordinate descent methods are often considered thanks to the separability of the proximal operator of ψ.
These are optimization algorithms that update only one coordinate of the vector of variables at each iteration, hence using partial derivatives rather than the whole gradient. In [11], Nesterov introduced the randomized coordinate descent method with an improved guarantee on the iteration complexity. He also gave an accelerated coordinate descent method for smooth functions. Lee and Sidford [5] introduced an efficient implementation of the method and Fercoq and Richtárik [4] developed the accelerated proximal coordinate descent method (APPROX) for the minimization of composite functions.
When solving a strongly convex problem, classical (non-accelerated) gradient and coordinate descent methods automatically have a linear rate of convergence, i.e. F (x k ) − F (x * ) ∈ O((1 − µ) k ) for a problem dependent 0 < µ < 1, whereas one needs to know explicitly the strong convexity parameter in order to set accelerated gradient and accelerated coordinate descent methods to have a linear rate of convergence, see for instance [5,6,7,11,12]. Setting the algorithm with an incorrect parameter may result in a slower algorithm, sometimes even slower than if we had not tried to set an acceleration scheme [13]. This is a major drawback of the method because in general, the strong convexity parameter is difficult to estimate.
In the context of accelerated gradient method with unknown strong convexity parameter, Nesterov [12] proposed a restarting scheme which adaptively approximate the strong convexity parameter. The similar idea was exploited by Lin and Xiao [8] for sparse optimization. Nesterov [12] also showed that, instead of deriving a new method designed to work better for strongly convex functions, one can restart the accelerated gradient method and get a linear convergence rate. However, the restarting frequency he proposed still depends explicitly on the strong convexity of the function and so O'Donoghue and Candes [13] introduced some heuristics to adaptively restart the algorithm and obtain good results in practice.

Contributions
In this paper, we show that we can restart accelerated gradient and coordinate descent methods, including APG, FISTA and APPROX, at any frequency and get a linearly convergent algorithm. The rate depends on an estimate of the strong convexity and we show that for a wide range of this parameter, one obtains a faster rate than without acceleration. In particular, we do not require the estimate of the strong convexity coefficient to be smaller than the actual value. In this way, our result supports and explains the practical success of arbitrary periodic restart for accelerated gradient methods.
In order to obtain the improved theoretical rate, we need to define a novel point where the restart takes place, which is a convex combination of previous iterates. Our approach is radically different from the previous restarting schemes [8,12,13], for which the evaluation of the gradient or the objective value is needed in order to verify the restarting condition. In particular, our approach can be extended to a restarted APPROX, which admits the same theoretical complexity bound as the accelerated coordinate descent methods for strongly convex functions [7] and exhibits better performance in numerical experiments.
In Sections 2 and 3 we recall the main convergence results for accelerated gradient methods. In Section 4, we present our restarting rules: one for accelerated gradient and one for accelerated coordinate descent. Finally, we present numerical experiments on the lasso and logistic regression problem in Section 5.

Problem and assumptions
For simplicity we present the algorithm in coordinatewise form. The extension to blockwise setting follows naturally (see for instance [4]). We consider the following optimization problem: where f : R n → R is a differentiable convex function and ψ : R n → R ∪ {+∞} is a closed convex and separable function: Note that this implies that each function ψ i : R → R is closed and convex. Let x * denote a solution of (1). We further assume that for each positive vector v = (v 1 , . . . , v n ), there is a constant µ F (v) > 0 such that where · v denotes the weighted Euclidean norm in R n defined by: Remark 1. Note that (2) is weaker than the usual strong convexity assumption on F . In particular, (2) does not imply that the objective function F is strongly convex. However, (2) entails the uniqueness of solution to (1) and by abuse of language, we refer to (2) as the strong convexity assumption.

Accelerated gradient schemes
In this paper, we are going to restart accelerated gradient schemes. We will concentrated on three versions that we will call FISTA (Fast Iterative Soft Thresholding Algorithm) [1], APG (Accelerated Proximal Gradient) [18] and APPROX (Accelerated Parallel and PROXimal coordinate descent) [4]. In the following, ∇f (y k ) denotes the gradient of f at point y k and ∇ i f (y k ) denotes the partial derivative of f at point y k with respect to the ith coordinate.Ŝ is a random subset of We have written the algorithms in a unified framework to emphasize their similarities. Practical implementations usually consider only two variables: (x k , y k ) for FISTA, (y k , z k ) for APG and (z k , w k ) where w k = θ −2 k−1 (x k − z k ) for APPROX. One may also consider t k = θ −1 k instead of θ k . The update in FISTA, APG or APPROX employs a positive vector v ∈ R n . To guarantee the convergence of the algorithm, the positive vector v should satisfy the so-called expected separable overapproximation (ESO) assumption, developed in [3,15] for the study of parallel coordinate descent methods.
where for h = (h 1 , . . . , h n ) ∈ R n and S ⊂ [n], h [S] is defined as: with e i being the ith standard basis vectors in R n .
We require that the positive vector v used in APPROX satisfy (3) with respect to the samplingŜ used. Note that FISTA shares the same constant v with APG and APG can be seen as a special case of APPROX whenŜ = [n]. Therefore the positive vector v used in FISTA and APG should then satisfy: which is nothing but a Lipschitz condition on the gradient of f . In other words, the vector v used in FISTA and APG is just the Lipschitz constant of ∇f , given a diagonal scaling that may be chosen to improve the conditioning of the problem. When in each step we update only one coordinate, we have τ = 1 and (3) reduces to: It is easy to see that in this case the vector v corresponds to the coordinate-wise Lipschitz constants of ∇f , see e.g. [11]. Explicit formulas for computing admissible v with respect to more general samplingŜ can be found in [15,3,14].

Convergence results for accelerated gradients methods
In this section we review two basic convergence results of FISTA and APPROX, which will be used later to build restarted methods. We first recall the following properties on the sequence {θ k }.
Remark 2. The strength of these propositions is that they are independent of the strong convexity parameter. Indeed, APG, FISTA and APPROX work for non-strongly convex minimization.

Remark 3.
As APPROX generalizes APG, we have covered all three algorithms in the two propositions. A remarkable feature is that the result for FISTA and APG are exactly the same even though the algorithms are different.

Restarted gradient methods
The basic tool upon which we build our restarting rule is a contraction property. We first present two restarting rules that require a special condition in order to guarantee the linear convergence. Then we present new rules that are more complex but are always certified to give a linearly convergent algorithm.

Conditional restarting
The first rule is an extension of the "optimal fixed restart" of [12,13] to FISTA and APPROX.
Proposition 3 (Conditional restarting at x k ). Let (x k , z k ) be the iterates of FISTA or APPROX applied to (1). We have Proof. By (8) and (10), the following holds for the iterates of FISTA (θ 0 = 1) and APPROX: Condition (15) is equivalent to: and we have the contraction using (5). Proposition 4 (Conditional restarting at z k ). Let (x k , z k ) be the iterates of FISTA or APG applied to (1).
Proof. By (8) and (10), the following holds for the iterates of FISTA and APG (θ 0 = 1): Here, we do not need to know µ F (v) but we need to wait for F (z k ) to be smaller than F (x k ). This event does happen sometimes but there is no guarantee for it to happen when minimizing a given function. Hence we may wait for ever and never restart.

Unconditional restarting
In this section, we prove the main results of this paper: setting a convex combination of the past iterates as the restarting point leads to a linearly convergent restarted method.
We first show that for full-gradient accelerated methods, an arbitrary strict convex combination of the last iterates x k and z k works.
Theorem 1 (Restarting for FISTA and APG). Let (x k , z k ) be the iterates of FISTA or APG applied to (1).
Proof. By the definition ofx k : Next we apply (8) and (9) (the same holds for APG by taking θ 0 = 1 in (10) and (11)): For APPROX, we need a more complex restarting point.
Theorem 2 (Restarting for APPROX). Let γ i k be the coefficients defined in (12) and We have Proof. Note that by (7), Hencex k is a convex combination of {x 0 , . . . , x k }. By (11) and the definition ofx k , In view of the strong convexity assumption (2), Moreover, using (11) again, we can easily see that

Restarted APPROX
We describe in Algorithm 4 the restarted APPROX method and give the convergence result in Theorem 3.

Remark 5.
For this restarting rule to be useful, we need to be able to computex k efficiently, in particular without computing x i for i < k. A way to do this is to use the variable w i = θ −2 i−1 (x i −z i ), which is maintained up-to-date in the algorithm: Then, we can compute the sum using cumulative updates like in [2]. We develop this idea in Appendix A.
Theorem 3. Let us choose K ∈ N and σ ∈ (0, 1) as we wish. Using the notation defined in Theorem 2, the iterates of Algorithm 4 satisfy for any k ≥ K
Proof. Let us write the Euclidean division k = mK + r with r ∈ [0, K − 1]. Using (11) and Theorem 2, Theorem 3 shows that Algorithm 4 is linearly convergent with respect to arbitrary convex combination coefficient σ ∈ (0, 1) and arbitrary restarting period K ∈ N. This implies that we can always get linear convergence without any information on the strong convexity parameter of the objective function. The next proposition provides an estimation on the rate of convergence given a guess µ on the parameter µ F (v) and a particular choice of σ and K. and The iterates of Algorithm 4 satisfy for any k ≥ K Proof. This is a direct corollary of Proposition 6 by taking λ = 1 + µ, presented in Appendix B.
Let us have additional insight of Proposition 5.
Let µ ∈ (0, 1]. Choose K and σ as in (20) and (21). Then for The proof of Corollary 1 is deferred to Appendix B. We therefore showed that for any strong convexity estimator µ ∈ (0, 1], the iteration complexity of Algorithm 4 is on the order of where the O notation hides logarithms of problem dependent constants and universal constants. Recall that for coordinate descent methods [11,16], the iteration complexity bound is Therefore, if µ is an upper bound on µ F (v), our iteration complexity bound improves over that of the randomized coordinate descent method [7,17] by a factor of we also obtain an improved complexity bound; only if µ < µ 2 F is our bound worse. This observation will be illustrated in Section 5 Remark 6. The theorems presented in this section can easily be combined with an adaptive restart strategy. We just need to define an interval [K,K] and allow the adaptive restart only if k ∈ [K,K]. Then if k =K we force the restart. We obtain a linear convergence rate where the rate is given by the worst case in the interval.

Illustration of the theoretical bounds
We first illustrate the theoretical rate we have found. On Figure 1, we can see that restarted APPROX has a better rate of convergence than vanilla proximal coordinate descent for a wide range of estimates of the strong convexity. Indeed, in this example with n = 10 and µ F (v) = 10 −5 , one can take 1.6 10 −9 ≤ µ ≤ 0.04. Note that this shows that even if the estimate µ is much larger than the true strong convexity coefficient µ F (v), we already see an improved rate. Yet, of course the closer µ is to µ F (v), the faster the algorithm will be.
On Figure 2, we fix the estimate of the strong convexity as µ = 10 −3 and we plot the rate of convergence of the method for µ F (v) ∈ [10 −9 , 1].

Gradient methods
We then present experiments on accelerated gradient methods (the case n = τ ). We solve the L 1 -regularised least squares problem (Lasso) min x∈R N 1 2 Ax − b    on the Iris dataset where A ∈ R m×n is the design matrix and b ∈ R m is such that b j = 1 if the label is "Iris-setosa", b j = −1 otherwise. We chose λ = A T b ∞ /10. This dataset is rather small (n = 4 and m = 8124). As we can see on Table 1, non-accelerated proximal gradient (underlined numbers) is faster than accelerated variants (italics). Yet, accelerated gradient methods designed for strongly convex objectives may be faster than both vanilla proximal gradient and basic accelerated proximal gradient. The "true" strong convexity coefficient is around 5.3 10 −4 : we can see that taking µ est close to this value leads indeed to a faster algorithm but that the algorithms are rather stable to approximations. Dual APG is quite efficient on this dataset but as its restarting rule is based on a divergence detection scheme, some attention should be paid before using intermediate solutions. Also remark that when it is not restarted, APG exhibits the sublinear O(1/k 2 ) rate while FISTA seems to be take some profit of strong convexity even without restarting.
All restarting strategies seem to perform well. Note however that APG-µ and FISTA-µ are only proved to converge linearly when µ ≤ µ F (v) and that the adaptive restart of [13] is a heuristic restart. Indeed, with APG, the restart condition did not happen within the first 10,000 iterations, which shows that this adaptive restart is not always efficient.
The rule of Theorem 3 is more complex and seems to be slightly less efficient than the rule of Theorem 1. We still present this more complex rule because it is the only one that is proved to have a linear rate of convergence in the accelerated coordinate descent case n > τ .  Table 1: Number of iterations to reach F (x k ) − F (x * ) ≤ 10 −10 with various accelerated algorithms for the Lasso problem on the Iris dataset. In italics, no restart has taken place; underlined numbers means that the algorithm is equivalent to non-accelerated ISTA.

Coordinate descent
We solve the following logistic regression problem: We consider In particular, for serial sampling (τ = 1), (3) is satisfied for Then for the latter v, (2) is satisfied for Even if the logistic objective is not strongly convex, we expect that the local curvature around the optimum is nonzero and so, that taking µ > µ F (v) = µ ψ will be useful. We solve (22) for different values of λ 2 , using v and µ F (v) = µ ψ defined in (23) and (24).
We compare randomized coordinate descent (CD), APCG [7] and APPROX-restart (Algorithm 4) with K and σ given by Proposition 5, on the dataset rcv1. We run both APCG and APPROX-restart using four different values of µ: µ F (v), 10 µ F (v), 100 µ F (v) and 1,000 µ F (v). We stop the program when the duality gap is lower than 10 −10 or the running time is larger than 3,000s. The results are reported in Figure 3, where by µ F we refer to µ F (v) defined in (24).
Note that the convergence of APCG is only proved for µ ≤ µ F (v) in [7]. In our experiments, we observed numerical issues when running APCG for several cases when taking larger µ (we were not able to compute the ith partial derivative at y k = ρ k w k + z k because ρ k had reached the double precision float limit). Such cases can be identified in the plots if the line corresponding to APCG stops abruptly before the time limit (3000s) with a precision worse than 10 −10 . On all the experiments, restarted APPROX is faster or much faster than APCG. Moreover, it is stable for any restarting frequency while APCG may fail if one is too optimistic when setting the strong convexity estimate.

A Efficient implementation of the restart point
Computingx k using (18) will be inefficient for APPROX because it involves full-dimensional operations. We show in the following that a reformulation similar to the one in [5] can be done to avoid such expensive computations. Define: We first recall from [4] that the update for {w k } can be as efficient as for {z k }: Define: It is easy to deduce from the recursive equation (6) that  Figure 3: Plots of results for dataset RCV1; APCG and APPROX-restart run with: µ = µ F (v) (first row), µ = 10µ F (v) (second row), µ = 100µ F (v) (third row) and µ = 1000µ F (v) (fourth row). Note that each column corresponds to the same problem, only the parametrization of the algorithms differ. APCG failed for (µ ψ = 1/n, µ = 1000µ ψ ) and (µ ψ = 0.1/n, µ = 1000µ ψ ).
Next we show that the sum vector (28) can be obtained efficiently using auxiliary coefficients and vectors. Define: Then clearly the update for {g k } and {h k } are also as efficient as for {z k }: Moreover, it is easy to see that Hence,x Finally we plug in (25) and (31) to obtain: .
The above reasoning showed that Algorithm 5 is equivalent to Algorithm 4. More importantly, note that full dimensional operations are avoided in Algorithm 5.
Lemma 2. For any k ≥ 1 we have: Proof. Let any k ≥ 1. We first decompose the sum and use (12) to obtain: Then we get the recursive equation: which together with (6) and (12) yields (33).
Lemma 3. For any k ≥ 1 we have Proof. We proceed by induction on k. First for k = 1 we have: Now suppose that we have Thus we proved by recurrence the following upper bound: , ∀k = 1, 2, . . . .
In particular, Now suppose that we have Next we apply (35) to obtain: We thereby proved the lower bound for any k.
Thus we just need to prove: We first note that if µ > µ F (v), then m K (µ) > m K (µ F (v)).


extfieldsofstudy:
Medicine
Computer Science

text:
What to Do When K-Means Clustering Fails: A Simple yet Principled Alternative Algorithm

The K-means algorithm is one of the most popular clustering algorithms in current use as it is relatively fast yet simple to understand and deploy in practice. Nevertheless, its use entails certain restrictive assumptions about the data, the negative consequences of which are not always immediately apparent, as we demonstrate. While more flexible algorithms have been developed, their widespread use has been hindered by their computational and technical complexity. Motivated by these considerations, we present a flexible alternative to K-means that relaxes most of the assumptions, whilst remaining almost as fast and simple. This novel algorithm which we call MAP-DP (maximum a-posteriori Dirichlet process mixtures), is statistically rigorous as it is based on nonparametric Bayesian Dirichlet process mixture modeling. This approach allows us to overcome most of the limitations imposed by K-means. The number of clusters K is estimated from the data instead of being fixed a-priori as in K-means. In addition, while K-means is restricted to continuous data, the MAP-DP framework can be applied to many kinds of data, for example, binary, count or ordinal data. Also, it can efficiently separate outliers from the data. This additional flexibility does not incur a significant computational overhead compared to K-means with MAP-DP convergence typically achieved in the order of seconds for many practical problems. Finally, in contrast to K-means, since the algorithm is based on an underlying statistical model, the MAP-DP framework can deal with missing data and enables model testing such as cross validation in a principled way. We demonstrate the simplicity and effectiveness of this algorithm on the health informatics problem of clinical sub-typing in a cluster of diseases known as parkinsonism.

In the generalized MAP-DP algorithm (Algorithm 3), the computation of the variables d i,k and d i, K+1 (algorithm lines 8,9) requires the collapsed prior predictive distribution f (x|θ 0 ), and also the collapsed posterior predictive distribution f x|θ −i k . This predictive distribution requires the updated cluster posterior hyper parameters θ −i k (algorithm line 7). These updates depend upon the distribution, and the data type, of each data point x i . When the distribution is from the exponential family, the prior distribution over the parameters can be chosen to be conjugate: the prior over the parameters of the data distribution and the posterior have the same form of distribution. This simplifies the hyper parameter updates, and, furthermore, the form of the prior and posterior predictive distributions is the same and is available in closed form. The table below lists some possible data types and distributions, their conjugate prior/posterior distribution, the names given to the hyper parameters and the corresponding name of the predictive distributions. We discuss each case in more detail in the subsequent sections. Multivariate normal

Spherical normal data with known variance
This is the variant of MAP-DP described in Algorithm 2. When each data point x ∈ R D is assumed to be spherical Gaussian with known varianceσ 2 shared across dimensions, the conjugate prior distribution of the Gaussian mean vector parameter µ ∈ R D is also spherical normal with hyper parameters θ 0 = µ 0 , σ 2 0 . Then the posterior distribution for each cluster is also spherical normal with hyper parameters The hyper parameter updates (Algorithm 3, line 7) for each cluster are: The predictive distributions f (x|θ 0 ) and f x θ −i k are D-dimensional spherical normal distributions, whose negative logs are: Note that since the normalization term D 2 ln (2π) is common to both predictive distributions, it can be omitted when computing d i,k and d i,K+1 in the algorithm.

Multivariate normal data with known covariance
For data points x ∈ R D assumed to be multivariate Gaussian with known covariance matrixΣ, the conjugate prior distribution of the Gaussian mean vector parameter is also multivariate normal with hyper parameters θ 0 = (µ 0 , Σ 0 ). The posterior distribution for each cluster is also multivariate normal with hyper parameters The hyper parameter updates are: The predictive distributions f (x|θ 0 ) and f x θ −i k are D-dimensional normal distributions, whose negative logs are: Since the normalization term D 2 ln (2π) is common to both predictive distributions, it can be omitted when computing d i,k and d i,K+1 in the algorithm.

Multivariate Gaussian data
When each data point x ∈ R D is assumed to be multivariate Gaussian with unknown mean vector and covariance matrix, the conjugate prior distribution of the Gaussian parameters is Normal-Wishart, with hyper parameters θ 0 = (m 0 , c 0 , B 0 , a 0 ). Then, the posterior distribution for each cluster is also Normal-Wishart, with hyper parameters These are updated for each cluster according to: The predictive distributions f (x|θ 0 ) and f x θ −i k are D-dimensional multivariate Student-t distributions, whose negative log, written in terms of the parameters (µ, Λ, ν) is: where the Student-t parameters (µ, Λ, ν) are given in terms of the Normal-Wishart parameters µ = m, ν = a − D + 1 and Λ = cν c+1 B. We note that fast incremental updates of all these parameters are possible when including and then removing a single data point from a cluster, see Raykov et al. [1] for further details.

Exponential data
Given data points x ∈ R, x ≥ 0 assumed to be exponentially-distributed, the conjugate prior over the exponential rate parameter is the gamma distribution. This gamma distribution has hyper parameters θ 0 = (α, β) (shape, rate). So, the posterior probability of the rate parameter is also gamma, and the cluster hyper parameter are updated using: The predictive distributions f (x|θ 0 ) and f x θ −i k are the so-called Lomax distribution, with negative log:

Categorical data
For categorical data which can take on one of D > 1 possible values, x ∈ {1, 2, . . . D}, the conjugate prior over the D outcome probability parameters of this distribution are Dirichlet distributed. This Dirichlet distribution has hyper parameters θ 0 = (α 0,1 , . . . , α 0,D ). So, the posterior outcome probability parameters for each cluster are also Dirichlet, and for each cluster the D entries in the cluster hyper parameter θ −i k = α −i k are updated using: where δ (x, y) = 1 if x = y and 0 otherwise. The predictive distributions f (x|θ 0 ) and f x θ −i k are special cases of the Dirichlet-multinomial distribution, with negative log:

Binomial data
In the case of binomial data where the data can take on x ∈ {0, 1, . . . n} for n > 0, the conjugate prior over the binomial success probability parameter is beta distributed, with hyper parameters θ 0 = (α 0 , β 0 ). By conjugacy, the posterior cluster parameters are also beta distributed with hyper parameters θ −i k = α −i k , β −i k , and are updated according to: For such binomial data, the predictive distributions f (x|θ 0 ) and f x θ −i k are beta-binomial, with negative log: where B (·, ·) is the beta function:

Poisson data
For positive integer Poisson count data x ∈ Z, x ≥ 0, the conjugate prior over the single rate parameter is the gamma distribution with hyper parameters θ 0 = (α 0 , β 0 ) (shape and rate, respectively). The posterior cluster parameters are similarly gamma distributed with hyper parameters The updates for these hyper parameters are: For Poisson count data, the predictive distributions f (x|θ 0 ) and f x θ −i k are negative binomial distributed with negative log:

Geometric data
In the case of positive integer data x ∈ Z, x ≥ 0 which is assumed to be geometrically-distributed, the conjugate prior over the single success probability parameter is the beta distribution with hyper parameters θ 0 = (α 0 , β 0 ). The posterior cluster parameters are similarly beta distributed with hyper parameters θ −i k = α −i k , β −i k . The updates for these hyper parameters are: For geometric data, the predictive distributions f (x|θ 0 ) and f x θ −i k have negative log: − ln f (x|θ) = − ln B (α + 1, β + x) + ln B (α, β) where B (·, ·) is the beta function described above.


s2fieldsofstudy:
Computer Science
Medicine

text:
Stacked Autoencoders for Medical Image Search

Medical images can be a valuable resource for reliable information to support medical diagnosis. However, the large volume of medical images makes it challenging to retrieve relevant information given a particular scenario. To solve this challenge, content-based image retrieval (CBIR) attempts to characterize images (or image regions) with invariant content information in order to facilitate image search. This work presents a feature extraction technique for medical images using stacked autoencoders, which encode images to binary vectors. The technique is applied to the IRMA dataset, a collection of 14,410 x-ray images in order to demonstrate the ability of autoencoders to retrieve similar x-rays given test queries. Using IRMA dataset as a benchmark, it was found that stacked autoencoders gave excellent results with a retrieval error of 376 for 1,733 test images with a compression of 74.61%.

Introduction
Many physicians have experienced lawsuits due to perceived or actual medical malpractice and negligence in recent decades. Radiology, due to its diagnostic nature, has been one of most liable branches in medicine. Most claims and complaints deal with a correct diagnosis [11]. For instance, wrong interpretation of a malignant mass as a benign lesion accounts for approximately 45% of radiologists' errors [8]. In addition, medical claims of misdiagnosis from medical images in Canada cost taxpayers over 8.3 million per year [2].
Radiologists do not have a reliable tool to cross-reference their initial diagnosis instantaneously. Receiving a second opinion on the diagnosis requires consulting a peer who must be pyhsically or virtually available, which is not always possible and constitutes an insurmountable financial and personal challenge for hospitals and clinics around the globe. Building a search engine for fast and accurate content-based image retrieval (CBIR) can potentially revolutionize diagnostic radiology and ultimately decrease both the number of patients affected by misdiagnosis and the number of costly claims against clinicians and hospitals. There are untapped petabytes of data already available in digital archives of hospitals -constantly growing -, that can be leveraged to reduce misdiagnosis through retrieving similar cases to exploit available knowledge from the past.
The challenge CBIR is that the image must be examined at the pixel level to quantify the similarity between images which is not an inherently quantitative attribute, especially in medical imaging. Additionally, when dealing with "big image data", the computational costs of retrieval can become infeasible, if the image representation, i.e., image feature, is not chosen appropriately. Accuracy and speed of retrieval are understandably crucial in medical image applications. This paper looks at medical image search based on global image similarity. Single layered and stacked autoencoders are examined to search for similar images within the IRMA dataset that contains 14,410 x-ray images. Similarity is measured using the IRMA error score. The idea proposed is that given an input x-ray image, the system would search for the top n = 3, 5, . . . most similar images to display to the clinicians (ideally accompanied with other corresponding information such as biopsy reports, treatment plans, monitoring and follow-ups). This would allow practitioners to leverage the diagnosis, monitoring and treatment results of other patients with similar cases.

Background
The main idea proposed in this paper is feature extraction through the use of autoencoders. The main purpose of extracting features from x-ray images is to create a high-level description from low-level pixel values data. The accuracy of image retrieval is highly dependent on the quality of the features extracted. If the features do not represent the image content adequately, similar images cannot be retrieved.
Image search and retrieval has been studied over the past 20 years [7]. CBIR focuses on visual information as opposed to textual metadata, which is textbased search as we all know from daily Internet search when we type keywords to search for desired webpages [10]. However, implementing autoencoders as a feature extraction technique for x-rays has only been studied recently. CBIR is a valuable option for medical images. Text-based search may be limited for medical images due to insufficient text-based data or features, as well as lack of standardized software and procedure in clinical enviroments. Akgül et al. examine the various applications of CBIR towards medical imaging [1].
Feature extraction is the process of consolidating high signal information from the noisy input, in this case an image. These features are valuable in a model for training and retrieval. Some examples of visual features are color, shape, and texture, which can be extracted from low-level pixel information [6]. These visual features can range from low-level descriptors that are focused on a pixel or small group of pixels, to mid-level features that describe shapes, textures and colour, and finally to high-level features that describe the image and its components in rather general or abstract ways.
Tradition feature extractions have been extensively examined in computer vision. Using artificial neural networks to "represent" an image is rather a new development in computer vision. An autoencoder is a special type of neural network that automatically finds compressed encodings of images, while minimizing error. The compressed image can then be used to represent the image to search algorithms. A stacked autoencoder inputs the results from the previous autoencoder. Hinton et al. introduced the idea of using backpropagation networks to develop a complex unsupervised learning scheme [15]. The autoencoders were trained through Deep Belief Networks [15]. A denoising autoencoder can be trained on corrupted input versions in order to reduce the impact of future input noise [19]. Autoencoders have also been applied to mammograms for compression [17]. Ideas have been proposed to detect (and to ignore) irrelevant image blocks in each medical image class, by analyzing the error histogram of a autoencoder [5]. The purpose in such approaches is to reduce the dimensionality of features for image retrieval when dealing with a large number of images. The relevance of image blocks is directly proportional to the error of an autoencoder [5]. However, this approach only uses shallow autoencoders that mostly depend on features used to create the prediction model. On the other hand, the use of deep autoencoders allow to extract a better representation of the raw features of the image, therefore allowing to create a more accurate retrieval scheme. Some considerations in evaluating the suitability of autoencoders for this domain were the accuracy of retrieving similar images. Providing fast and accurate results to practitioners is necessary. Autoencoders may require longer training time dependent on the dimensionality of the problem, however this does not impact the retrieval time, once all images have been processed.

Proposed Approach
Fast and accurate CBIR is possible by extracting features from images that have the following two intuitive properties: two very similar images should produce two very similar feature vectors, and two very dissimilar images should produce two very different feature vectors. These two properties are especially critical for medical imaging where a region of interest needs to be analyzed (e.g., a tumour, an organ, a tissue type). A system is only valuable if it can return cases which contain regions similar to the selected region of interest.
Some CBIR solutions use local feature from feature descriptors such as SURF [18]. This technique creates many local feature descriptors for a given image and then combines all of them (up to a threshold) to create a pseudo-global image descriptor. This approach works well for general images but, as we verified in many preliminary experiments, fails in medical imaging, specifically for x-ray images. This could be due to x-rays not having many sharp or "good" features, e.g., corners, for the algorithms to use. To circumvent this problem, a global first approach is to be used to create one global feature descriptor (vector) for a given image. One such approach is to use autoencoders to generate a compressed feature vector.
Autoencoders, with n/p/n architecture, encodes n inputs into p positions, and then decode p positions back into n outputs. By setting p < n, the autoencoder functions as a compressor to reduce the dimensionality. Such an autoencoder is basically a shallow neural network with some level of error to reconstruct the input signal from its minimalistic representation in the deepest layer.
While it is possible to use many visual features such as color, shape and texture as part of a feature extraction pipeline, we propose a very simple scheme.
We use the unprocessed pixel values as features into the autoencoder. The autoencoder, by virtue of its dimensionality reduction capability, will extract the complex anatomical features using the reconstruction error as its guide. We hypothesize that these features will satisfy the two properties mentioned above since the features are constructed relative to the training set, a claim that we seek to experimentally establish in following sections. The architecture can be deepened by creating an architecture of form n/m/k/p/k/m/n (5 hidden layers) with n < m < k < p by stacking multiple single layer autoencoders serially. Stacked autoencoders (SAE) can capture highly nonlinear mapping between input and output from the interactions between the many hidden layers and thousands of trainable weights.

Experiments
In this section, we describe the IRMA dataset, report the details of training and testing the stacked autoencoder, and finally analyze the results of image retrieval using autoencoders applied on IRMA dataset.

Image Data
Our goal is to provide the n most similar image to an trained medical professional, the clinician, when she/he provides a query image. Images can be verified to be similar using many techniques such as RMS, and SSIM, however these metrics do not reflect medical domain information. To verify if two x-rays are indeed similar, the images would need to be shown to a medical professional. The IRMA data set contains images already annotated by medical professionals, which makes it a very useful benchmark dataset for retrieval purposes.
The IRMA dataset, supplied for the imageCLEF organization, had annotated codes were similar to domain expert knowledge [9,12,13,14]. This database has been used by many researchers and is comprised of preset test and training portions, which enable direct comparisons of metrics. The IRMA codes (manually created by several clinicians) contains information on technical, biological and diagnostic traits of the image in a structured manner: TTTT-DDD-AAA-BBB. Each section is hierarchical meaning there is a least significant bit and most significant digit (Table 1). Sample IRMA images are depicted in Fig. 1.  As mentioned before, the objective of any CBIR system is to retrieve the most similar image from a dataset given an input (query) image. To build the database, each image in the training dataset is encoded using a trained Autoencoder. Processing the image through the autoencoder results in a compressed (encoded) feature vector that is assumed to "describe" the image content. This feature vector is created for all 12,998 x-ray images of IRMA dataset and stored. To retrieve an image, the query image is put through the trained autoencoder and a query feature vector is generated. The feature vector is then compared with all other indexed feature vectors to locate the most similar image. Various autoencoder specifications were tried, their results are stated in following sections.
IRMA Error Score -Let I =l 1 ,l 2 , . . . ,l i , . . . ,l I be the classified IRMA code (for one axis: T,D,A or B) of an image; where l i is set precisely for every position, and inli is used when we 'don't know', marked by '*'. I may be different for different images. If the classification at positionl i is wrong, then all succeeding decisions are wrong and, given a not-specified position, all succeeding decisions are considered to be not specified. if the correct code is unspecified and the predicted code is a wildcard, then no error is counted. In such cases, all remaining positions are regarded as not specified. Wrong (easy) decisions (i.e., fewer possible choices) are penalized for wrong difficult decisions (many possible choices at that node). A decision at position l i is correct by chance with a probability of 1 bi if b i is the number of possible labels for position i. This assumes equal priors for each class at each position. Furthermore, wrong decisions at an early stage in the code (higher up in the hierarchy) are more penalized than wrong decisions at a later stage in the code (lower down on the hierarchy): i.e., l i is more important than l i+1 . Considering all these thoughts: with δ(I i ,Î i ) being 0, 0.5 or 1 for agreement, 'don't know' and disagreement, respectively. In order to normalize the error, the maximal possible error is calculated for every axis such that a completely wrong decision (i. The result is normalized on a scale 0 to 1. The IRMA error score is accumulated over the entire test set (1,734 images). The resulting number is the IRMa score reported in the results section. Since each IRMA error can range from 0 to 1, the sum of all error divided by the total possible error gives the error percentage for the test set.

Architecture of Autoencoder
We investigate both single layer autoencoders n/p/n where p n and stacked autoencoders where the output of the previous autoencoder is fed to the next and so on. For the actual search we use the k-NN algorithm.
Training single layer autoencoder -We take 12,998 training x-ray images from the IRMA dataset. Each image is gray-scaled, and downscaled to 32 × 32 and normalized in [0, 1]. It is then fed to the input layer of the autoencoder as a column vector of size 1024. The encoding and decoding layers are trained to reconstruct the input vector (down sampled and vectorized image) by minimizing the reconstruction error measured by cross entropy using stochastic gradient descent. Note that we are using tied weights for encoding and decoding weights which has a regularization effect. We use sigmoid function for nonlinearity. Finally, we obtain the feature vector for each image by obtaining the latent representation from the hidden layer for each image.
Training stacked autoencoder -We obtain stacked autoencoders by serially arranging autoencoders such that the output of the first layer is fed as input to the second and so on (Fig. 2). Each autoencoder layer is trained greedily. Latent features from the first autoencoder are used as input to the second autoencoder layer. The process is repeated for all subsequent hidden layers. Finally, the latent feature vector from the last encoding layer acts as the final feature vector for the input image. We use Theano [3,4] and various python libraries to implement both the autoencoder and stacked autoencoder [16].

Results
We report compression reduction, IRMA score and reconstruction error (root mean squared) for various AEs (autoencoders) and SAEs (stacked autoencoders).
We tried various configurations for a single layer autoencoder where p < n as reported in Table 2. We used Euclidean distance to compare feature vectors between two images. The first hit (1-NN) IRMA score as well as the reconstruction error is reported on the test set. Interestingly, better reconstruction error did not necessarily lead to lower IRMA score. For example, while the 1024/512/1024 configuration clearly had the lowest reconstruction error amongst the reported single layer architectures, it did not yield the lowest IRMA score. This implies that the features trained by minimizing the reconstruction error heuristic are clearly useful, we cannot easily predict the performance of the features on the test dataset from the reconstruction error alone. We note that the latent features produced significant compression while still maintaining low error. Each trial used the following settings: 30 epochs, a learning rate of 0.1, and a batch size of 20 for stochastic gradient descent.
We also include results from stacking two and three autoencoders sequentially (see Tables 3 and 4). The reconstruction errors are reported on the deepest layer. The 1024/600/1024, 600/500/600, 500/260/500 architecture yielded the lowest IRMA score of 376 while still maintaining relatively high compression (74.61% reduction). We believe stacked autoencoders, as opposed to a single layer autoencoder, produce a more useful higher-level representation from the lower-level representation output by the previous layer. To test for overfitting on the training data, the reconstruction error was calculated on the training set and compared to the testing images. For example, a triple stacked autoencoder 1024/600/1024, 600/500/600, 500/260/500, the error was 0.096200 and 0.101679, respectively. One can interpret the negligible difference in reconstruction errors as acceptable generalization. It should be noted that reconstruction error for stacked autoencoders is reported only for the first layer, and therefore the true reconstruction error is not simple to obtain. However, the differences in reconstruction for a simple non-stacked autoencoder (see Table 4) are also negligible.
Observing the distribution of data in the IRMA training and testing images, we noticed that there exists both inter-and intra-statistical imbalances between the data sets. That is, the training set is severely imbalanced such that various classes are not represented normally. A histogram of image classes is shown in Fig. 3. Despite this problem, the trained autoencoders have a low IRMA score and exhibit negligible differences in reconstruction error (between training and test data). This can further compound the notion that there is no overfitting, hence the autoencoders are generalizing well.

Conclusions
Content-based image retrieval is a challenging problem for medical images. It is possible to analyze one image in a short amount of time and classify it in a medically meaningful way. However, given a new image, comparing it to the entire collection of medical images, to find the similar images, can be computationally expensive. Such a task requires a technique that searches through the big image data efficiently while still finding highly similar images. One method to reduce the dimensionality of the problem is to encode the images in a smaller, compact representation to enable faster image comparisons during the retrieval task. Such technology has clearly enormous potential in aiding domains such as diagnostic radiology and pathology. Due to the scarcity of text-based features associated with medical images directly embedded in medical image formats, i.e., in DICOM files, as well as the large quantity of images available, a search focused on visual features was explored in this paper. A global, content-based image retrieval technique was successfully applied to the application of medical images. The best performing feature extraction technique examined was an stacked autoencoder. Of the examined architectures, the most accurate setup consisted of three hidden layers with configurations 1024/600/1024, 600/500/600, and 500/260/500. On the IRMA dataset, with 14,410 x-ray images, the proposed approach achieved an IRMA score of 376 with a compression of 74.61%.


Cross-domain Text Classification with Multiple Domains and Disparate Label Sets

Advances in transfer learning have let go the limitations of traditional supervised machine learning algorithms for being dependent on annotated training data for training new models for every new domain. However, several applications encounter scenarios where models need to transfer/adapt across domains when the label sets vary both in terms of count of labels as well as their connotations. This paper presents first-of-its-kind transfer learning algorithm for cross-domain classification with multiple source domains and disparate label sets. It starts with identifying transferable knowledge from across multiple domains that can be useful for learning the target domain task. This knowledge in the form of selective labeled instances from different domains is congregated to form an auxiliary training set which is used for learning the target domain task. Experimental results validate the efficacy of the proposed algorithm against strong baselines on a real world social media and the 20 Newsgroups datasets.

Introduction
A fundamental assumption in supervised statistical learning is that training and test data are independently and identically distributed (i.i.d.) samples drawn from a distribution. Otherwise, good performance on test data cannot be guaranteed even if the training error is low. On the other hand, transfer learning techniques allow domains, tasks, and distributions used in training and testing to be different, but related. It works in contrast to traditional supervised techniques on the principle of transferring learned knowledge across domains. Pan and Yang, in their survey paper (2010), de- scribed different transfer learning settings depending on if domains and tasks vary as well as labeled data is available in one/more/none of the domains. In this paper, we propose a generic solution for multi-source transfer learning where domains and tasks are different and no labeled data is available in the target domain. This is a relatively less chartered territory and arguably a more generic setting of transfer learning.
Motivating example: Consider a social media consulting company helping brands to monitor their social media channels. Two problems typically of interest are: (i) sentiment classification (is a post positive/negative/neutral?) and (ii) subject classification (what was the subject of a post?). While sentiment classification attempts to classify a post based on its polarity, subject classification is towards identifying the subject (or topic) of the post, as illustrated in Figure 1. The company has been using standard classification techniques from an off-the-shelf machine learning toolbox. While machine learning toolkit helps them to create and apply statistical models efficiently, the same model can not be applied on a new collection due to variations in data distributions across collections 1 . It requires a few hundreds of manually labeled posts for every task on every collec-tion. As social media are extremely high velocity and low retention channels, human labeling efforts act like that proverbial narrow bottleneck. Need of the hour was to reduce, if not eliminate, the human-intensive labeling stage while continue to use machine learning models for new collections.
Several transfer learning techniques exist in the literature which can reduce labeling efforts required for performing tasks in new collections. Tasks such as sentiment classification, named entity recognition (NER), part of speech (POS) tagging that have invariant label sets across domains, have shown to be greatly benefited from these works. On the other hand, tasks like subject classification that have disparate label sets across domains have not been able to gain at pace with the advances in transfer learning. Towards that we formulate the problem of Cross-domain classification with disparate label sets as learning an accurate model for the new unlabeled target domain given labeled data from multiple source domains where all domains have (possibly) different label sets.
Our contributions: To the best of our knowledge, this is the first work to explore the problem of cross-domain text classification with multiple source domains and disparate label sets. The other contributions of this work includes a simple yet efficient algorithm which starts with identifying transferable knowledge from across multiple source domains useful for learning the target domain task. Specifically, it identifies relevant class-labels from the source domains such that the instances in those classes can induce classseparability in the target domain. This transferable knowledge is accumulated as an auxiliary training set for an algorithm to learn the target domain classification task followed by suitable transformation of the auxiliary training instances.
Organization of the paper is as follows: Section 2 presents the preliminaries and notation, Section 3 summarizes the related work. Section 4 and 5 present the proposed algorithm and experimental results respectively. Section 6 concludes the paper.

Preliminaries and Notations
A domain D = {X , P (X)} is characterized by two components: a feature space X and a marginal probability distribution P (X), where X = {x 1 , x 2 , ...x n } ∈ X . A task T = {Y, f (·)} also consists of two components: a label space Y and an objective predictive function f (·).
In our settings for cross-domain classification with disparate label sets, we assume M source domains, denoted as D S i , where i = {1, 2, ..M }. Each source domain has different marginal distribution i.e. P (X S i ) = P (X S j ) and different label space i.e. Y S i = Y S j , ∀i, j ∈ M . The label space across domains vary both in terms of count of class-labels as well as their connotations; however, a finite set of labeled instances are available from each source domain. The target domain (D T ) consists of a finite set of unlabeled instances, denoted as t i where i = {1, .., N }. Let Y T be the target domain label space with K class-labels. We assume that the number of classes in the target domain i.e. K is known (analogous to clustering where the number of clusters is given). Table 1 summarizes different settings of transfer learning (Pan and  and how this work differentiates from the existing literature 2 . The first scenario represents the ideal settings of traditional machine learning (Mitchell, 1997) where a model is trained on a fraction of labeled data and performs well for the same task on the future unseen instances from the same domain.

Related Work
The second scenario where the domains vary while the tasks remain the same is referred to as transductive transfer learning. This is the most extensively studied settings in the transfer learning literature and can be broadly categorized as single and multi-source adaptation. Single source adaptation (Chen et al., 2009;Ando and Zhang, 2005;Daumé III, 2009) primarily aims at minimizing the divergence between the source and target domains either at instance or feature levels. The general idea being identifying a suitable low dimensional space where transformed source and target domains data follow similar distributions and hence, a standard supervised learning algorithm can be trained (Daumé III, 2009;Jiang and Zhai, 2007;Pan et al., 2011;Blitzer et al., 2007;Dai et al., 2007;Bhatt et al., 2015).
While several existing single source adaptation techniques can be extended to multi-source adaptation, the literature in multi-source adaptation can be broadly categorized as: 1) feature representation approaches (Chattopadhyay et al., 2012;Sun et al., 2011;Duan et al., 2009;Duan et al., 2012;

Transfer
Learning unlabeled data from the target domain P(X S ) = P(X T )

Multi-source adaptation
Classifier combination; efficient combination of information from multiple sources; Feature representation Intelligent selection of transferable knowledge from multiple sources for adaptation.

Inductive
Unlabeled data in source domain(s) and labeled data in target domain

Self-taught learning
Extracts higher level representations from unlabeled auxiliary data to learn instance-to-label mapping with labeled target instances Learns instance-to-label mapping in the unlabeled target domain using multiple labeled source domains having different data distributions and label spaces.
Labeled data is available in all domains

Multi-task learning
Simultaneously learns multiple tasks within (or across) domain(s) by exploiting the common feature subspace shared across the tasks Learns the optimal class distribution in an unlabeled target domain by minimizing the differences with multiple labeled source domains.
Labeled data in source and target domains Transfer learning with disparate label set Disparate fine grained label sets across domains, however, same coarse grained labels set can be invoked across domains No coarse-to-fine label mapping due to heterogeneity of label sets, Assumes no labelled data in target domain. Bollegala et al., 2013;Crammer et al., 2008;Mansour et al., 2009;Ben-David et al., 2010;Bhatt et al., 2016) and 2) combining pre-trained classifiers (Schweikert and Widmer, 2008;Sun and Shi, 2013;Xu and Sun, 2012;. Our work differentiates in intelligently exploiting selective transferable knowledge from multiple sources unlike existing approaches where multiple sources contribute in a brute-force manner. The third scenario where the tasks differ irrespective of the relationship among domains is referred to as inductive transfer learning. Self-taught learning (Raina et al., 2007) and multi-task (Jiang, 2009;Maurer et al., 2012;Xu et al., 2015;Kumar and Daume III, 2012) learning are the two main learning paradigms in this scenario and Table 1 differentiates our work from these.
This work closely relates to the fourth scenario where we allow domains to vary in the marginal probability distributions and the tasks to vary due to different label spaces 3 . The closest prior work by Kim et al. (2015) address a sequential labeling problem in NLU where the fine grained label sets across domains differ. However, they assume that there exists a bijective mapping between the coarse and fine-grained label sets across domains. They learn this mapping using labeled instances from the target domain to reduce the problem to a standard domain adaptation problem (Scenario 2). elaborated in the next sections.

Exploiting Multiple Domains
If we had the mappings between the source and target domain label sets, we could have leveraged existing transfer learning approaches. However, heterogeneity of label sets across domains and the unlabeled data from the target domain exacerbate the problem. Our objective is to leverage the knowledge from multiple source domains to induce class-separability in the target domain. Inducing class-separability refers to segregating the target domain into K classes using labeled instances from selective K source domain classes.
Towards this, the proposed algorithm divides each source domain into clusters/groups based on the class-labels such that instances with the same label are grouped in one cluster. All source domains are divided into Q clusters where Q = M m=1 ||Y m || represents the count of class-labels across all sources. ||Y m || being the count of classlabels in the m th source domain. C q denotes the q th cluster and µ q denotes its centroid computed as the average of all the members in the cluster. We assert that the target domain instances that have high similarity to a particular source domain cluster can be grouped together. Given N target domain instances and Q source domain clusters, a matrix R (dimension N × Q) is computed based on the similarity of the target domain instances with the source clusters. The i th row of the matrix captures the similarity of the i th target domain instance (t i ) with all the source domain clusters. It captures how different source domain class-labels are associated with the target domain instances and hence, can induce class-separability in the target domain.

Extracting Transferable Knowledge
The similarity matrix R associates target domain instances to the source domain clusters in proportion to their similarity. However, the objective is to select the optimal K source domain clusters that fit the maximum number of target domain instances. This problem is similar to the well-known combinatorial optimization problem of Maximum Coverage (Vazirani, 2003) where given a collection of P sets, we need to select A sets (A < P ) such that the size of the union of the selected sets is maximized. In this paper, we are given Q source domain clusters and need to select K clusters such that the corresponding number of associated tar-get domain instances is maximized. As the Maximum Coverage problem is NP-hard, we implement a greedy algorithm for selecting the k source domain clusters, as illustrated in Algorithm 1.

Algorithm 1 Selecting K Source Clusters
Input: A matrix R, K = number target domain classes, l= number of selected cluster. Initialize: l = 0, Normalize R such that each row sums up to 1. repeat: 1: Pick the column in R which has maximum sum of similarity scores for uncovered target domain instances. 2: Mark elements in the chosen column as covered. 3: l = l + 1 until: l = K Output: K source domain clusters.
A source domain contributes partially in terms of zero or more class-labels (clusters) identified using the Algorithm 1. Therefore, we refer to the labeled instances from the selected clusters of a source domain as the partial transferable knowledge from that domain. This partial transferable knowledge from across multiple source domains is congregated to form an auxiliary training set, referred to as (AU X).

Adapting to the Target Domain
The auxiliary training set comprises labeled instances from selected K source domain clusters 4 . Since, the auxiliary set is pulled out from multiple source domains, it follows different data distribution as compared to the target domain. For a classifier, trained on the K-class auxiliary training set, the distributional variations have to be normalized so that it can generalize well on the target domain.
In this research, we proposed to use an instance weighting technique (Jiang and Zhai, 2007) to minimize the distributional variations by deferentially weighting instances in the auxiliary set. Intuitively, the auxiliary training instances similar to the target domain are assigned higher weights while training the classifier and vice versa. The weight for the i th instance in the auxiliary set should be proportional to the ratio (Pt(x i )) (Pa(x i )) . However, since the actual probability distributions C: Adapting to target domain: 1: Minimize distributional variations using instance weighing technique. 2: Train a K-class classifier using AU X.
Output: K-class target domain classifier.
(P a (x) and P t (x) for the auxiliary set and target domain respectively) are unknown, the instance difference is approximated as (Pt(x i |d=target)) (Pa(x i |d=auxilliary)) , where d is a random variable used to represent whether x i came from the auxiliary set or the target domain. To calculate this ratio, a binary classifier is trained using the auxiliary set and target domain data with labels {-1} and {+1} respectively. The predicted probabilities from the classifier are used to estimate the ratio as the weight for the i th auxiliary instance x i . Finally, a K-class classifier is trained on the weighted auxiliary training set to perform classification on the target domain data.

Algorithm
As shown in Figure 2, the step-by-step flow of the proposed algorithm is summarized below: 1. Divide M source domains into Q clusters, each represented as C q , q = {1, 2, .., Q}.
2. Compute centroid of each cluster as the average of the cluster members, as shown in Eq. 1.
where µ q is the centroid, ||C q || is the membership count and x i is the i th member of C q .
3. For target instances t i ∀i ∈ N , compute cosine similarity with all the source domain cluster centroids to form the matrix R (dimensions: N × Q), as shown in Eq. 2 4. Run Algorithm 1 on R to select K optimal source clusters (i.e. columns of R).
5. Congregate labeled instances from the selected source domain clusters to form the Kclass auxiliary training set.
6. Minimize the divergence between the auxiliary set and target domain using the instance weighing technique, described in Section 4.3.
7. Finally, train a K-class classifier on deferentially weighted auxiliary training instances to perform classification in the target domain.
The K-class classifier trained on the auxiliary training set is an SVM classifier (Chih-Wei Hsu and Lin, 2003) with L2 − loss from the LIB-LINEAR library (Fan et al., 2008). The classifier used in the instance weighing technique is again an SVM classifier with RBF kernel. The proposed algorithm uses distributional embedding i.e. Doc2Vec (Le and Mikolov, 2014) to represent instances from the multiple source and target domains. We used an open-source implementation of Doc2Vec (Le and Mikolov, 2014) for learning 400 dimensional vector representation using DBoW.

Experimental Evaluation
Comprehensive experiments are performed to evaluate the efficacy of the proposed algorithm for cross-domain classification with disparate label sets across domains on two datasets.

Datasets
The first dataset is a real-world Online Social Media (OSM) dataset which consists of 74 collections. Each collection comprises comments/tweets that are collected based on user-defined keywords. These keywords are fed to a listening engine which crawls the social media (i.e. Twitter.com) and fetches comments matching the keywords. The task is to classify the comments in a collection  into user-defined categories. These user-defined categories may vary across collections in terms of count as well as their connotations. Table 2 shows an example of the user-defined categories for a few collections related to "Apple" products. In the experiments, one collection is used as unlabeled target collection and the remaining collections are used as the labeled source collections. We randomly selected 5 target collections to report the performance, as described in Table 3. The second dataset is the 20 Newsgroups (NG) (Lang, 1995) dataset which comprises 20, 000 news articles organized into 6 groups with different sub-groups both in terms of count as well as connotations, as shown in Figure 3(a). Two different experiments are performed on this dataset. In the first experiment ("Exp-1"), one group is considered as the target domain and the remaining 5 groups as the source domains. In the second experiment ("Exp-2"), one sub-group from each of the first five groups 5 is randomly selected to synthesize a target domain while all the groups (with the remaining sub-groups) are used as source domains. Figure 3(b) shows an example on how to synthesize target domains in "Exp-2". There are 720 possible target domains in this experiment and we report the average performance across all possible target domains, referred to as "Grp 7". The task in both the experiments is to categorize the target domain into its K categories (sub-groups) using labeled data from multiple source domains.

Evaluation Metric
The performance is reported in terms of classification accuracy on the target domain. There is no definite mapping between the actual class-labels in the target domain and the K categories (i.e. induced categories) in the auxiliary training set. Therefore, we sequentially evaluate all possible one-to-one mappings between the K categories in the auxiliary training set and target domain to report results for the best performing mapping.

Experimental Protocol
The performance of the proposed algorithm is skylined by the in-domain performance (Gold), i.e. a classifier trained and tested on the labeled target domain data. We also compared the performance with spherical K-means clustering (Dhillon and Modha, 2001) used to group the target domain data into K categories against the ground truth, referred to CL. Spherical K-means clustering is based on cosine similarity and performs better for high-dimensional sparse data such as text.
To compare with a baseline and an existing adaptation algorithm, we selected the most similar source domain 6 with exactly K number of classlabels and report the performance on the best possible mapping, as described in Section 5.2. To compute the baseline (BL), a classifier trained on the source domain is used to categorize the target domain. A widely used domain adaptation algorithm, namely structural correspondence learning (SCL) (Blitzer et al., 2007) is also applied using the selected source domain.

Results and Analysis
Key observations and analysis from the experimental evaluations are summarized below:

Results on the OSM Dataset
Results in Figure 4 and Table 4 show the efficacy of the proposed algorithm for cross-domain classification with disparate label sets as it outperforms other approaches by at least 15%. Coll ID(#) refers to the target collection and the corresponding count of class-labels. Results in Table  4 also compare the performance of the proposed technique without the distributional normalization of the auxiliary training set, referred to as "W/O". Results suggest that suitably weighing instances from the auxiliary training set mitigates the distributional variations and enhances the cross-domain performance by at least 3.3%.

Results on the 20Newsgroups Dataset
Results in Table 5 show that the proposed algorithm outperforms other techniques for both the experiments by at least 15 % and 18% respectively on the 20 Newsgroups dataset. In Table 5, "-" refers to the cases where a single source domain with the same number of class-labels as in the target domain is not available. In "Exp-1" where the source and target categories vary in terms of counts as well as their connotations, the proposed algorithm efficiently induces the classes in the unlabeled target domain using the partial transferable knowledge from multiple sources. For "Exp-2", it is observed that the performance of the proposed algorithm is better than the performance in "Exp-1" as the target categories have closely related categories (from the same group) in the source do- Figure 5: Effects of selected source collections on the OSM dataset. mains. Table 5 reports the average performance across all the 720 possible combinations of target domains with a standard deviation of 2.6. Table 6 validates our assertion that multiple sources are necessary to induce class-separability in the target domain as a single source is not sufficient to cater to the heterogeneity of class-labels across domains. It also suggests that the proposed algorithm can learn class-separability in the target domain by using arbitrary diverse class-labels from different sources and does not necessarily require class-labels to follow any sort of coarse-tofine mapping across domains.

Effect of Multiple Source Domains
To evaluate the effects of using multiple sources, further experiments were performed by varying the number of available source domains. For the OSM dataset, we varied the number of available source collections from 1 to 73 starting with the most similar source collection and repeatedly adding the next most similar collection in the pool of available collections. We observe that even the most similar collection was not independently sufficient to induce classes in the target collection and it was favorable to exploit multiple collections. Moreover, adding collections based on similarity to the target collection had a better likelihood of achieving higher performance as compared to adding random collections.
In another experiment, we first identified the source collections which contributed to learning the target task. We removed these collections and applied the proposed algorithm on the remaining source collections. Figure 5 shows the perfor- mance of the proposed algorithm on 5 such iterations of removing the contributing source collections from the previous iteration. We observed a significant drop in the performance with each iteration which signifies the effectiveness of the proposed algorithm in extracting highly discriminating transferable knowledge from multiple sources.

Comparing with Domain Adaptation
We applied domain adaptation techniques considering the auxiliary training set to be a single source domain with the same number of classes as that in the target domain. We applied two of the widely used domain adaptation techniques, namely SCL (Blitzer et al., 2007) and SFA  referred to as "AuxSCL" and "AuxSFA" respectively. Results in Table 7 suggest that the proposed algorithm significantly outperforms "AuxSCL" and "AuxSFA" on the two datasets. Generally, existing domain adaptation techniques are built on the co-occurrences of the common features with the domain specific features and hence, capture how domain specific features in one domain behaves w.r.t to the domain specific features in the other domain. They assume homogeneous labels and expect the aligned features across domains to behave similarly for the prediction task. However, these features are misaligned when the label set across domains vary in terms of their connotations.

Effect of Different Representations
The proposed algorithm uses Doc2Vec (Le and Mikolov, 2014) for representing instances from multiple domains. However, the proposed algorithm can build on different representations and hence, we compare its performance with traditional TF-IDF representation (including unigrams  and bigrams) and a dense representation using TF-IDF+PCA ( reduced to a dimension such that it covers 90% of the variance). We observe that Doc2Vec representation clearly outperforms the other two representations as it addresses the drawbacks of bag-of n-gram models in terms of implicitly inheriting the semantics of the words in a document and offering a more generalizable concise vector representation.

Conclusions
This paper presented the first study on crossdomain text classification in presence of multiple domains with disparate label sets and proposed a novel algorithm for the same. It proposed to extract partial transferable knowledge from across multiple source domains which was beneficial for inducing class-separability in the target domain. The transferable knowledge was assimilated in terms of selective labeled instances from different source domain to form a K-class auxiliary training set. Finally, a classifier was trained using this auxiliary training set, following a distribution normalizing instance weighing technique, to perform the classification task in the target domain. The efficacy of the proposed algorithm for cross-domain classification across disparate label sets will expand the horizon for ML-based algorithms to be more widely applicable in more general and practically observed scenarios.


extfieldsofstudy: Mathematics Physics

text:
Entropy-reducing dynamics of a double demon

We study the reduction in total entropy, and associated conversion of environmental heat into work, arising from the coupling and decoupling of two systems followed by processing determined by suitable mutual feedback. The scheme is based on the actions of Maxwell's demon, namely the performance of a measurement on a system followed by an exploitation of the outcome to extract work. When this is carried out in a symmetric fashion, with each system informing the exploitation of the other (and both therefore acting as a demon), it may be shown that the second law can be broken, a consequence of the self-sorting character of the system dynamics.

We study the reduction in total entropy, and associated conversion of environmental heat into work, arising from the coupling and decoupling of two systems followed by processing determined by suitable mutual feedback. The scheme is based on the actions of Maxwell's demon, namely the performance of a measurement on a system followed by an exploitation of the outcome to extract work. When this is carried out in a symmetric fashion, with each system informing the exploitation of the other (and both therefore acting as a demon), it may be shown that the second law can be broken, a consequence of the self-sorting character of the system dynamics.
Entropy production may be viewed as the irreversible growth in uncertainty of the microscopic state of a system, arising from the complexity of the underlying dynamics and the incompleteness of the model we employ to represent it. It can be quantified using a framework of stochastic rules describing the forward evolution in time of a system influenced by its coarse grained environment [1]. Nevertheless, it is possible to conceive of procedures, often involving a time-asymmetric feedback mechanism, that operate against the usual tendency for processes to operate irreversibly, a classic example being Maxwell's demon [2]. The argument is that a measurement of a system can reveal a route for exploitation that leads to a reduction in entropy; typically the conversion of environmental heat to potential energy represented by the raising of a weight.
The original demon exploited his observations to sort molecules of a gas into fast and slow groups, creating a resource for a heat engine without the expenditure of work. Szilard's later conception of a demon-operated heat engine made this explicit [3]. Maxwell saw no edict against breaking the second law by time-asymmetric dynamical processing [4][5][6][7], but much attention has been given to finding a way to 'exorcise the demon' and protect the second law in these circumstances. The majority view is that dissipative processes, operating either in the act of measurement [3,[8][9][10][11][12][13][14] or the act of restoring the initial condition of the demon prior to measurement [15][16][17][18], generate enough entropy to cancel out any possible gains [19,20].
But we can imagine dynamical schemes that emulate a successful demon and we are obliged either to accept that they are possible, or to find reasons to exclude them. Consider, for example, a particle tethered by a harmonic spring to a point and coupled to a heat bath. The tether point might be moved instantaneously towards the current position of the particle, relaxing the spring and harvesting potential energy to lift a weight. Subsequently, the system will evolve back towards equilibrium under the influence of the heat bath, with its mean potential energy replenished through the absorption of heat. There is no act of measurement and no demon here, just an autonomous dynamical system that employs feedback, breaking time reversal symmetry. We might call such dynamics self-sorting and this system a self-adjusting oscillator.
On the other hand, we might feel uneasy about a system that alters the dynamical rules that control its future behaviour, depending on its current state, and question whether such a system is admissible for thermodynamic consideration. We might prefer to channel the feedback instead through a second entity, a demon, such that the system dynamics cannot be described as self-sorting. The system would be coupled to the demon, or more prosaically a measuring device, in such a way that establishes a correlation between their coordinates, and then decoupled. The subsequent exploitation of the system would then be determined by the device coordinate. Since the effects of the feedback are felt by the system after the device is decoupled, there is no element of self-sorting in the dynamics.
Analysis of such a procedure indeed shows that the second law is preserved. The process of coupling and decoupling must be entropy generating if it is to yield an exploitable correlation between system and device [10,12,14]. This entropy production exceeds the reduction in entropy made possible by the measurement, and so overall the entropy is never decreased. This can be regarded as a satisfactory outcome.
However, if we accept that the state of a device can inform the subsequent exploitation of the system, it is natural to consider a symmetric arrangement where the state of the system is allowed to inform an exploitation of the device. Both parts act as a demon and we might therefore call this a realisation of a 'double demon'. Intuitively, we suspect that the combination might display self-sorting behaviour. The feedback mechanism preserves the second law when applied in a demon-system context, but when operated mutually the effect might be different.
Our purpose is to analyse the dynamics and thermodynamics of a double demon constructed using harmonic springs, and to investigate its irreversibility. We find that a suitable exploitation protocol makes possible an overall entropy reduction.

arXiv:1605.02202v1 [cond-mat.stat-mech] 7 May 2016
We consider two 1-d harmonic oscillators that can be coupled and decoupled through a further harmonic spring. Both oscillators are influenced by noise from the environment. We consider four intervals of time. In the period −∞ ≤ t ≤ −τ m the coupling spring strength K is zero, the oscillator spring strengths are unity. We take kT to equal unity so at t = −τ m an equilibrium state is established described by the probability density function (pdf) p(x, y, −τ m ) = (2π) −1 exp[−(x 2 + y 2 )/2] where x and y are the oscillator displacements.
In the next time interval −τ m ≤ t ≤ 0 the oscillators interact through a nonzero K(t). The situation is illustrated in Fig. 1, where the shading of the springs indicates their strength. We consider the coupling spring constant to evolve to a value K(0) = K 0 just prior to t = 0 (labelled t = 0 − ) and then to go to zero abruptly.
We model x and y in this period using overdamped stochastic differential equations representing the effects of the various spring forces, with white environmental noise described by increments dW x and dW y in separate Wiener processes. The work performed during this 'measurement' period is The pdf satisfying the Fokker-Planck equation in the measurement interval takes the form with the time-dependent parameterK(t) determined by withK(−τ m ) = 0. In the quasistatic limit,K(t) follows the evolution of the coupling strength K(t) exactly. The stochastic entropy production during measurement is given by [1,14,21,22] which implies a mean rate of production At t = 0, the joint pdf may be represented as p(x 0 , y 0 , 0) = p(x 0 , 0|y 0 )p(y 0 , 0) where x 0 = x(0) and y 0 = y(0), with the state of the first oscillator described by a conditional pdf Figure 1. Two oscillators in a noisy environment are coupled and decoupled in the period −τm ≤ t ≤ 0 using a harmonic spring of strength K(t). The spring strengths Kx,y are adjusted and tether points λx,y shifted, using feedback from the displacements x0 and y0 at t = 0, in preparation for an exploitation procedure in the period t > 0.
whereK 0 =K(0), and the second oscillator described by with κ = (1 + 2K 0 )/(1 +K 0 ). If the measurement process took the form of quasistatic coupling and instantaneous decoupling (denoted 'qi'), for whichK 0 would equal K 0 , the mean work of measurement would be a free energy of coupling minus the mean energy of the coupling spring at t = 0, which is In the third period 0 ≤ t ≤ τ e the outcome of the measurement is exploited through changes in the strengths K x (t) and K y (t) of the oscillator springs and the positions λ x and λ y of their tethering points. Optimal exploitation sequences for a single oscillator after measurement have been studied previously [10,12,19]. We similarly employ λ x = K 0 y 0 /(1 + K 0 ), λ y = K 0 x 0 /(1 + K 0 ), and K x (0) = K y (0) = 1+K 0 , the rationale for which will become clear, and K x (τ e ) = K y (τ e ) = 1. The changes at the start of the 'exploitation' interval are indicated in Fig. 1 (labelled t = 0 + ). The subsequent dynamics are modelled using and the work of exploitation is given by (11) If the pdf p(x 0 , y 0 , 0) is characterised byK 0 = K 0 and if instantaneous changes to the spring strengths and tether points are followed by a quasistatic evolution of K x and K y (a process denoted 'iq'), then the mean work of exploitation is given by for each oscillator. The mean work of measurement and exploitation for a protocol consisting of quasistatic coupling, instantaneous decoupling and changes in spring parameters, followed by further quasistatic processing (denoted 'qiiq') is therefore for exploitation of just the first oscillator (i.e. where K y and λ y are not modified for the exploitation interval). Since this mean work cannot be negative, the second law is preserved. However, if both oscillators are subjected to the exploitation procedure we get which offers the prospect of a violation of the law.
To check such a claim, we compute the stochastic entropy production during exploitation. We write the pdf as p(x, y, t) = p(x, t|y 0 )p(y, t|x 0 ), such that at t = 0 the conditional pdf p(x, t|y 0 ) takes the form in Eq. (7) and p(y, t|x 0 ) is given by Eq. (8). This asymmetric specification of the pdf provides a rationale for the chosen exploitation protocol. If a quasistatic-instantaneous (qi) measurement procedure has been performed, such that K 0 = K 0 , then the changes made to K x and λ x at t = 0 effectively place the first oscillator in a state of canonical equilibrium. A subsequent quasistatic evolution of the spring strength of the first oscillator could be carried out without the generation of entropy.
In general, however, the evolving pdf of the first oscillator would be written (15) with time-dependent parametersK x (t) andx(t) given by subject toK x (0) = 1 +K 0 andx(0) =K 0 y 0 /(1 +K 0 ). The stochastic entropy production associated with the first oscillator during exploitation is [21] For a qi measurement processK 0 = K 0 andx(0) = λ x , in which case d∆s x tot reduces to leading to an average rate of production Similarly, the pdf of the second oscillator may be written (21) for the exploitation interval, with but this time subject toK y (0) = (1 + 2K 0 )/(1 +K 0 ) and y(0) = 0. The stochastic entropy production associated with the second oscillator during exploitation is [21] d∆s y It is instructive to consider a special case where we as-sumeK 0 = K 0 , set K y equal to 1+K 0 for the initial part of the exploitation interval, and compute a stochastic entropy production d∆s rel,y tot associated with the relaxation of the parametersK y andȳ to K y and λ y (x 0 ), respectively. This is the key to understanding the breakage of the second law. The mean rate of production, averaged over x 0 and y 0 , is given by . Note that this mean rate of production is negative. After this relaxation, d ∆s y tot /dt takes the form of an analogue of Eq. (20) associated with the deviation betweenK y and K y over the remainder the exploitation interval.
In the final period τ e ≤ t ≤ ∞ the oscillator spring constants have returned to unity and the system relaxes until there is no further stochastic entropy production. The shifts in the tether points are irrelevant to the irreversibility, since a quasistatic process can take the λ x,y back to zero at no cost in mean work or mean entropy production. This completes the cycle.
The only contribution to the mean total entropy production for a qiiq process involving exploitation of both  Figure 2. Evolution of the mean work performed on the two oscillator system (solid lines), and the associated mean stochastic entropy production (dashed lines), obtained from numerical simulations of the process specified in the main text, conducted over a range of timescales τ = τm = τe. The quantities become equal when equilibrium is restored at the end of the process, and they approach a negative limit when the process timescale goes to infinity.
oscillators is the relaxation represented by Eq. (25). When integrated over time this gives ∆s rel,y tot = ∆s qiiq,2 which is consistent with the negative mean work W qiiq,2 m+e in Eq. (14) obtained from a cycle of measurement and exploitation under such conditions. If exploitation is not invoked for the second oscillator, its relaxation would take place while K y remained equal to unity and with λ y = 0. The mean total entropy production for this case can be shown to be consistent with the positive mean work W qiiq,1 m+e in Eq. (13). Thus the mechanical and entropic assessments of irreversibility are consistent with one another for the special case of the qiiq process.
If the double demon oscillator system is processed slowly enough, we therefore expect environmental heat to be converted into stored potential energy, when averaged over many realisations, with a maximum change in the total thermodynamic entropy per cycle given by Eq. (26). This arises from the exploitation of each oscillator in a manner informed by the state of the other at t = 0, a symmetrising of the feedback traditionally employed in a system-demon scenario. We next study how this might emerge in numerical simulations.
We choose K 0 = 4 and generate sets of realisations for various values of τ m and τ e . The averages of work done and stochastic entropy produced are shown in Fig.  2 using a timestep of 5 × 10 −4 and 10 6 trajectories for each case. The two quantities coincide at the end of the cycle, as the laws of thermodynamics suggest they should, and for slower processes they approach the negative limit given by Eqs. (14) and (26). The reduction in entropy is achieved through the feedback invoked at t = 0 and the partial harvesting of the potential energy of the springs. However, positive mean entropy production arises from the deviations ofK,K x andK y from K, K x and K y , respectively, and this occurs more strongly for faster processes.
We have shown that the stochastic evolution of two oscillators evolving according to a particular form of autonomous time-asymmetric dynamics can break the second law. Their interactions are conceived as an extension of a scenario where a measurement made by a device or demon is used to inform the exploitation of a system in order to convert environmental heat into work. In that case the second law is preserved because of a prior work of measurement, but passing feedback in both directions, thus allowing the device to be manipulated in a fashion informed by the system, extracts additional work. The total thermodynamic entropy can be reduced after completion of such a sequence of measurement and exploitation, and we attribute this to the self-sorting dynamics of the double demon.
We thank Stefan Grosskinsky and Rosemary J. Harris for helpful discussions, and acknowledge support from Engineering and Physical Sciences Research Council (EPSRC) Grant No. EP/I01358X/1, and the COST1209 network.