Convex Covariate Clustering for Classification

TL;DR

This paper introduces a convex covariate clustering method that incorporates class label information to improve interpretability and compression in classification tasks, with an efficient optimization algorithm and a new model selection criterion.

Contribution

It presents a novel convex clustering approach that integrates class labels, along with a scalable ADMM solver and a marginal likelihood-based model selection criterion.

Findings

Effective clustering on synthetic data

Improved interpretability on real datasets

Scalable to thousands of variables

Abstract

Clustering, like covariate selection for classification, is an important step to compress and interpret the data. However, clustering of covariates is often performed independently of the classification step, which can lead to undesirable clustering results that harm interpretability and compression rate. Therefore, we propose a method that can cluster covariates while taking into account class label information of samples. We formulate the problem as a convex optimization problem which uses both, a-priori similarity information between covariates, and information from class-labeled samples. Like ordinary convex clustering [Chi and Lange, 2015], the proposed method offers a unique global minima making it insensitive to initialization. In order to solve the convex problem, we propose a specialized alternating direction method of multipliers (ADMM), which scales up to several thousands of variables. Furthermore, in order to circumvent computationally expensive cross-validation, we propose a model selection criterion based on approximating the marginal likelihood. Experiments on synthetic and real data confirm the usefulness of the proposed clustering method and the selection criterion.

Tables

Table 1. Table 1: Shows the ANMI score of the proposed method, convex clustering, and k 𝑘 k -means clustering on synthetic data. Standard deviation in brackets. For hyper-parameter selection the proposed marginal likelihood criteria is used. Small number shows the ANMI score when cross-validation is used instead of the marginal likelihood.

Proposed Model
		n
		40	400	4000
d	40	0.84 (0.11) 0.89 (0.03)	0.71 (0.2) 0.9 (0.0)	1.0 (0.0)
	200	0.87 (0.16) 0.87 (0.13)	0.95 (0.03) 0.93 (0.04)	1.0 (0.0)
	1000	0.94 (0.0)	0.94 (0.0)	1.0 (0.0)
Convex Clustering
		n
		40	400	4000
d	40	0.59 (0.04) 0.58 (0.05)	0.54 (0.18) 0.48 (0.23)	0.38 (0.27)
	200	0.69 (0.0) 0.69 (0.0)	0.69 (0.0) 0.69 (0.0)	0.48 (0.31)
	1000	0.01 (0.01)	0.01 (0.0)	0.0 (0.0)
$k$-means Clustering
		n
		40	400	4000
d	40	0.55 (0.09) 0.52 (0.12)	0.43 (0.14) 0.43 (0.14)	0.35 (0.21)
	200	0.56 (0.2) 0.58 (0.09)	0.57 (0.15) 0.53 (0.14)	0.54 (0.13)
	1000	0.65 (0.06)	0.7 (0.0)	0.41 (0.08)

Table 2. Table 2: Average runtime (in minutes with standard deviation in brackets) of the proposed ADMM solution, and two standard solvers ECOS (Domahidi et al., 2013 ) and SCS (O’Donoghue et al., 2016 ) . When a method was too slow to perform 10 repeated experiments, we either ran only 3 experiments, or estimated the expected runtime based on a few iterations, marked as ( ∗ ) (*) .

Solving with proposed ADMM
		n
		40	400	4000
d	40	5.314 (0.503)	12.271 (0.488)	63.318 (1.339)
	200	19.331 (4.031)	86.815 (4.246)	187.734 (4.101)
	1000	102.721 (7.896)	552.763 (51.104)	982.877 (12.316)
Solving with ECOS
		n
		40	400	4000
d	40	6.75 (0.841)	18.936 (0.186)	360.677 (11.68)
	200	147.6 (35.325)	196.136 (8.52)	3800(*)
	1000	15400(*)121212For this setting ECOS did not converge properly, which led to exceptional long runtime.	6990(*)	65100(*)
Solving with SCS
		n
		40	400	4000
d	40	11.333 (1.159)	155.13 (5.767)	1590 (*)
	200	163.145(33.505)	360.96(20.762)	1680(*)
	1000	15300(*)	5730(*)	23100(*)

Table 3. Table 3: Results for Newsgroup20 and IMDB dataset. Shows the result of each clustering method, when hyper-parameters are selected with the approximate marginal likelihood. For reference, also the results when performing no clustering, i.e. the full model, are shown.

Newsgroup20
Method	Marginal Likelihood	Nr. Clusters	Accuracy
Proposed	-9487.1	1173	0.86
Convex	-10695.8	1538	0.88
$k$-means	-10152.4	1036	0.85
No Clustering	-10695.8	1538	0.88
IMDB
Proposed	-578.4	140	0.85
Convex	-3637.6	2604	0.85
$k$-means	-3589.7	2528	0.85
No Clustering	-3637.6	2604	0.85

Table 4. Table 4: Results for Newsgroup20 and IMDB dataset. Comparison of held-out accuracy of the proposed method, convex clustering and k 𝑘 k -means for around the same number of clusters: 100, 500, 1000.

Newsgroup20
Nr. Clusters	Proposed	Convex	$k$-means
100	0.67	0.44	0.67
500	0.80	0.70	0.79
1000	0.84	0.82	0.83
IMDB
Nr. Clusters	Proposed	Convex	$k$-means
100	0.85	0.59	0.8
500	0.84	0.76	0.82
1000	0.84	0.8	0.82

Table 5. Table 5: Shows the ANMI score of the proposed method, convex clustering, and k 𝑘 k -means clustering on synthetic data: prior covariate similarity is fully agreeing as shown in Figure 2 . Standard deviation in brackets. For hyper-parameter selection the proposed marginal likelihood criteria is used. Small number shows the ANMI score when cross-validation is used instead of the marginal likelihood.

Proposed Model
		n
		40	400	4000
d	40	0.7 (0.12) 0.0 (0.0)	0.83 (0.16) 1.0 (0.0)	1.0 (0.0)
	200	0.88 (0.07) 1.0 (0.0)	0.92 (0.1) 1.0 (0.0)	1.0 (0.0)
	1000	0.95 (0.07)	0.96 (0.06)	1.0 (0.0)
Convex Clustering
		n
		40	400	4000
d	40	0.93 (0.13) 0.93 (0.13)	1.0 (0.0) 1.0 (0.0)	1.0 (0.0)
	200	0.99 (0.04) 0.98 (0.04)	1.0 (0.0) 1.0 (0.0)	1.0 (0.0)
	1000	0.04 (0.02)	0.0 (0.0)	0.0 (0.0)
$k$-means Clustering
		n
		40	400	4000
d	40	0.72 (0.21) 0.66 (0.22)	0.84 (0.2) 0.77 (0.25)	0.99 (0.03)
	200	0.92 (0.11) 0.88 (0.16)	0.93 (0.07) 0.88 (0.09)	0.99 (0.02)
	1000	0.98 (0.03)	1.0 (0.0)	0.95 (0.07)

Equations

$$f(y|\mathbf{x},B,{\mathbf{\beta }}_0)=\frac{\text{exp}(B_{y,\cdot }\mathbf{x}+{\mathbf{\beta }}_0(y))}{{\sum }_{y^{\prime }}\text{exp}(B_{y^{\prime },\cdot }\mathbf{x}+{\mathbf{\beta }}_0(y^{\prime }))}.$$

$$\underset{B,{\mathbf{\beta }}_0}{\mathrm{minimize}}$$

$$+\nu \sum _{i_1<i_2}{S}_{i_1,i_2}||B_{\cdot ,i_1}-B_{\cdot ,i_2}|{|}_2,$$

$$\underset{B,{\mathbf{\beta }}_0}{\mathrm{minimize}}$$

$$+\nu \frac{1}{2}\sum _{i=1}^d\sum _{j=1}^{d_i}{S}_{i,e_i(j)}||B_{\cdot ,i}-B_{\cdot ,e_i(j)}|{|}_2,$$

$$\tilde{f}(y_s|{\mathbf{x}}_s,B,{\mathbf{\beta }}_0):=f(y_s|{\mathbf{x}}_s,B,{\mathbf{\beta }}_0)\cdot \text{exp}(-\frac{\lambda }{n}||B|{|}_F^2).$$

$$\underset{B,Z,{\mathbf{\beta }}_0}{\mathrm{minimize}}$$

$$+\nu \frac{1}{2}\sum _{i=1}^d\sum _{j=1}^{d_i}{S}_{i,e_i(j)}||{\mathbf{z}}_{i\to e_i(j)}-{\mathbf{z}}_{e_i(j)\to i}|{|}_2$$

$$\forall i\in \{1,\dots d\},j\in \{1,\dots ,d_i\}:{\mathbf{z}}_{i\to e_i(j)}-{\mathbf{b}}_i=0,$$

$$B^{k+1}:=\underset{B,{\mathbf{\beta }}_0}{\mathrm{argmin}}g(B,{\mathbf{\beta }}_0),$$

$$Z_{\cdot \to \cdot }^{k+1}:=\underset{Z_{\cdot \to \cdot }}{\mathrm{argmin}}\nu \frac{1}{2}\sum _{i=1}^d\sum _{j=1}^{d_i}{S}_{i,e_i(j)}||{\mathbf{z}}_{i\to e_i(j)}-{\mathbf{z}}_{e_i(j)\to i}|{|}_2$$

$$+\frac{\rho }{2}\sum _{i=1}^d\sum _{j=1}^{d_i}||{\mathbf{z}}_{i\to e_i(j)}-{\mathbf{b}}_i^{k+1}+{\mathbf{u}}_{i\to e_i(j)}^k|{|}_2^2,$$

$$\forall \{i,j\}:i\in \{1,\dots d\},j\in \{1,\dots ,d_i\}:$$

$${\mathbf{u}}_{i\to j}^{k+1}:={\mathbf{u}}_{i\to j}^k+{\mathbf{z}}_{i\to j}^{k+1}-{\mathbf{b}}_i^{k+1},$$

$$g(B,{\mathbf{\beta }}_0):=$$

$$+\frac{\rho }{2}\sum _{i=1}^d\sum _{j=1}^{d_i}||{\mathbf{z}}_{i\to e_i(j)}^k-{\mathbf{b}}_i+{\mathbf{u}}_{i\to e_i(j)}^k|{|}_2^2.$$

$$\sum _{i=1}^d\sum _{j=1}^{d_i}||{\mathbf{z}}_{i\to e_i(j)}^k-{\mathbf{b}}_i+{\mathbf{u}}_{i\to e_i(j)}^k|{|}_2^2$$

$$\displaystyle=q_{all}-2\sum_{i=1}^{d}\mathbf{b}_{i}^{T}\mathbf{q}\uparrow_{i}+\sum_{i=1}^{d}d_{i}\Big{(}\mathbf{b}_{i}^{T}\mathbf{b}_{i}\Big{)}\,,$$

$${\nabla }_{{\mathbf{b}}_i}g(B,{\mathbf{\beta }}_0)$$

$$\forall \{i,j\}:i\in \{1,\dots d\},j\in \{1,\dots ,d_i\}:$$

$${\mathbf{z}}_{i\to j}^{k+1},{\mathbf{z}}_{j\to i}^{k+1}:=\underset{{\mathbf{z}}_{i\to j},{\mathbf{z}}_{j\to i}}{\mathrm{argmin}}\nu S_{i,j}||{\mathbf{z}}_{i\to j}-{\mathbf{z}}_{j\to i}|{|}_2+\frac{\rho }{2}||{\mathbf{z}}_{i\to j}-{\mathbf{b}}_i^{k+1}+{\mathbf{u}}_{i\to j}^k|{|}_2^2$$

$$+\frac{\rho }{2}||{\mathbf{z}}_{j\to i}-{\mathbf{b}}_j^{k+1}+{\mathbf{u}}_{j\to i}^k|{|}_2^2.$$

$${\mathbf{z}}_{i\to j}^{k+1}:=\theta ({\mathbf{b}}_i^{k+1}-{\mathbf{u}}_{i\to j}^k)+(1-\theta )({\mathbf{b}}_j^{k+1}-{\mathbf{u}}_{j\to i}^k),$$

$${\mathbf{z}}_{j\to i}^{k+1}:=(1-\theta )({\mathbf{b}}_i^{k+1}-{\mathbf{u}}_{i\to j}^k)+\theta ({\mathbf{b}}_j^{k+1}-{\mathbf{u}}_{j\to i}^k),$$

$$\displaystyle\theta:=\begin{cases}0.5&\text{if }m=0,\\ \max\Big{(}1-\frac{\nu S_{i,j}}{\rho h}\,,\,0.5\Big{)}&\text{otherwise. }\end{cases}$$

$$p(\mathbf{y},B|X,\mathcal{C},{\mathbf{\beta }}_0)=\pi (B)\prod _{s}f(y_s|{\mathbf{x}}_s^{\mathcal{C}},B,{\mathbf{\beta }}_0)$$

$${\mathbf{x}}^{\mathcal{C}}=T^{\mathcal{C}}\mathbf{x},$$

$$\displaystyle T^{\mathscr{C}}_{ij}=\left\{\begin{array}[]{rl}1&\text{if }j\in\mathcal{C}_{i}\,,\\ 0&\text{else }\,.\end{array}\right.$$

$$\log f(y|\mathbf{x},B,{\mathbf{\beta }}_0)+\frac{1}{n}\log \pi (B)=$$

$$\displaystyle\quad B_{y,\cdot}\mathbf{x}+\boldsymbol{\beta}_{0}(y)-\log\Big{(}\sum_{y^{\prime}}\text{exp}(B_{y^{\prime},\cdot}\mathbf{x}+\boldsymbol{\beta}_{0}(y^{\prime}))\Big{)}$$

$$-\frac{1}{n}(\frac{1}{2{\sigma }^2}||B|{|}_F^2+\frac{cm}{2}\log (2\pi {\sigma }^2)).$$

$$S(B{)}_{ij}:=\sum _{s=1}^n{\nabla }_{{\mathbf{r}}_i,{\mathbf{r}}_j}^2\{\log f(y_s|{\mathbf{x}}_s,B,{\mathbf{\beta }}_0)+\frac{1}{n}\log \pi (B)\},$$

$$p(X|{\mathbf{\beta }}_0)$$

$$\sum _{i=1}^d\sum _{j=1}^{d_i}||{\mathbf{z}}_{i\to e_i(j)}^k-{\mathbf{b}}_i+{\mathbf{u}}_{i\to e_i(j)}^k|{|}_2^2$$

$$=\sum _{i=1}^d\sum _{j=1}^{d_i}||{\mathbf{q}}_{ij}-{\mathbf{b}}_i|{|}_2^2$$

$$\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{d_{i}}\Big{(}\mathbf{q}_{ij}^{T}\mathbf{q}_{ij}-2\cdot\mathbf{b}_{i}^{T}\mathbf{q}_{ij}+\mathbf{b}_{i}^{T}\mathbf{b}_{i}\Big{)}$$

$$\displaystyle=\sum_{i=1}^{d}\sum_{j=1}^{d_{i}}\Big{(}\mathbf{q}_{ij}^{T}\mathbf{q}_{ij}\Big{)}-2\sum_{i=1}^{d}\mathbf{b}_{i}^{T}\Big{(}\sum_{j=1}^{d_{i}}\mathbf{q}_{ij}\Big{)}+\sum_{i=1}^{d}\sum_{j=1}^{d_{i}}\Big{(}\mathbf{b}_{i}^{T}\mathbf{b}_{i}\Big{)}$$

$$\displaystyle=q_{all}-2\sum_{i=1}^{d}\mathbf{b}_{i}^{T}\mathbf{q}\uparrow_{i}+\sum_{i=1}^{d}d_{i}\Big{(}\mathbf{b}_{i}^{T}\mathbf{b}_{i}\Big{)}\,,$$

$${\nabla }_{{\mathbf{b}}_i}g(B,{\mathbf{\beta }}_0)$$

$$\displaystyle\quad+\rho\sum_{j=1}^{d_{i}}\Big{(}-\mathbf{z}_{i\rightarrow e_{i}(j)}^{k}+\mathbf{b}_{i}-\mathbf{u}_{i\rightarrow e_{i}(j)}^{k}\Big{)}$$

Full text

1 Introduction

Interpretability is paramount to communicate classification and regression results to the domain expert. Especially in high-dimensional problems, sparsity of the solution is often assumed to increase interpretability. As a consequence, various previous work in statistics focuses on covariate selection, with the $\ell_{1}$-penalty being a particular popular choice (Hastie et al., 2015).

However, even after successful covariate selection, the remaining solution might still contain several hundreds or thousands of covariates. We therefore suggest to further cluster the covariates. In particular, we propose a clustering method that uses a-priori similarity information between covariates while respecting the response variable of the classification problem. The resulting covariate clusters can help to identify meaningful groups of covariates and this way, increase the interpretability of the solution.

As a motivating example, consider the situation of document classification, where we are interested in classifying into positive and negative movie reviews. Assuming, a logistic regression model, words with similar regression coefficients might convey similar sentiment, e.g. ”wonderful” and ”great”, and therefore can be considered as one cluster. On the other hand, we might also have prior word similarity information from word embeddings (Mikolov et al., 2013) that can help us to identify meaningful clusters.

In order to incorporate, both types of information, class labels and a-priori covariate similarities, we propose to formulate the clustering problem as a joint classification and covariate clustering problem. In particular, we use a logistic regression loss with a pair-wise group lasso penalty for each pair of covariate regression coefficient vector. Our formulation leads to a convex optimization problem, like convex clustering (Chi and Lange, 2015; Hocking et al., 2011). As a consequence, we have the desirable properties of a unique global minima, and the ease to plot a clustering hierarchy, instead of just a single clustering.

Our proposed method is conceptually related to joint convex covariate clustering and linear regression as proposed in (She, 2010). However, the change from linear regression to logistic regression is computationally challenging. The main reason is that the objective function is not decomposable for each covariate anymore. Therefore, we propose a specialized alternating direction method of multipliers (ADMM) with an efficient gradient descent step for the non-decomposable primal update. Our solution allows us to scale the covariate clustering to problems with several 1000 covariates.

Since we often want to decide on one clustering result, we also propose a model selection criterion using an approximation to the marginal likelihood. The motivation for this criterion is similar to the Bayesian information criterion (BIC) (Schwarz, 1978), but prevents the issue of a singular likelihood function. By using the marginal likelihood, we can circumvent the need for hyper-parameter selection with cross-validation, which is computationally infeasible when the number of covariates is large.

The outline of this article is as follows. In the next section, we introduce our proposed objective function, and describe an efficient ADMM algorithm for solving the convex optimization problem. In Section 3, we describe our model selection criterion for selecting a plausible clustering result out of all clustering results that were found with the proposed method. In our experiments, in Sections 4 and 5, we compare the proposed method to a two-step approach, first step: covariate clustering with $k$-means clustering or convex clustering (Chi and Lange, 2015); second step: classification.

Finally, in Section 6, we summarize our conclusions.

A note on our notation: we denote a matrix by capital letter, e.g. $B\in\mathbb{R}^{c\times d}$, and a column vector by bold font e.g. $\mathbf{x}\in\mathbb{R}^{d}$. Furthermore, the i-th row of $B$ is denoted by $B_{i,\cdot}$ and is a row vector. The j-th column of $B$ is denoted by $B_{\cdot,j}$ or simply $\mathbf{b}_{j}$, and is a column vector. For a vector $\mathbf{x}$, we denote by $||\mathbf{x}||_{2}$ the $\ell_{2}$ norm, and for a matrix $B$, $||B||_{F}$ denotes the Frobenius norm.

2 Proposed Method

Let $B\in\mathbb{R}^{c\times d}$, where $c$ is the number of classes, and $d$ is the number of covariates. $B_{l,\cdot}$ is the weight vector for class $l$. Furthermore, $\boldsymbol{\beta}_{0}\in\mathbb{R}^{c}$ contains the intercepts. We assume the multi-class logistic regression classifier defined by

[TABLE]

We propose the following formulation for jointly classifying samples $\boldsymbol{x}_{s}$ and clustering the covariates:

[TABLE]

where $S_{i_{1},i_{2}}\geq 0$ defines a similarity measure between covariate $i_{1}$ and $i_{2}$ and is assumed to be given a-priori. The last term is a group lasso penalty on the class weights for any pair of two covariates $i_{1}$ and $i_{2}$. The penalty is large for similar covariates, and therefore encourages that $B_{\cdot,i_{1}}-B_{\cdot,i_{2}}$ is $\mathbf{0}$, that means that $B_{\cdot,i_{1}}$ and $B_{\cdot,i_{2}}$ are equal. The clustering of the covariates can be found by grouping two covariates $i_{1}$ and $i_{2}$ together if $B_{\cdot,i_{1}}$ and $B_{\cdot,i_{2}}$ are equal.

The advantage of this formulation is that, as long as $\lambda>0$, the model is identifiable, and the optimization problem (2) is strongly convex, and we are therefore guaranteed to find the unique global minima.111On the other hand, if $\lambda=0$, and, for example, the optimal solution for $B$ has all entries equal to some constant, then the problem in Equation (2) is not strongly convex, but only convex. In practice, a value $\lambda>0$ can potentially help to speed up convergence, and improve numerical stability. For all our experiments we set $\lambda$ to $0.1$.

Note that this penalty shares some similarity to convex clustering as in (Chi and Lange, 2015; Hocking et al., 2011). However, one major difference is that we do not introduce latent vectors for each data point, and our method can jointly learn the classifier and the clustering.

We remark that it is straight forward to additionally add a sparsity penalty on the columns of $B$ to jointly select and cluster all covariates. We omit in the following such extensions and focus on the computationally difficult part, the clustering. Our implementation222Released here https://github.com/andrade-stats/convexCovariateClusteringClassification. allows to perform joint or sequential covariate selection and clustering.

Finally, we remark that the similarity matrix can be represented as an undirected weighted graph with an edge between covariates $i_{1}$ and $i_{2}$ iff $S_{i_{1},i_{2}}>0$. In practice, $S$ might be manually crafted from a domain expert (e.g. given by some ontology), or learned from data (as we do in Sections 4 and 5).

2.1 Optimization using ADMM

Here, we focus on the optimization problem from Equation (2), which we can rewrite as

[TABLE]

where $d_{i}$ denotes the number of adjacent covariates of $i$, and

[TABLE]

Assuming that the adjacent covariates of $i$ are ordered from 1 to $d_{i}$, the function $e_{i}(j)$ returns the global covariate id of the $j$-th adjacent covariate of $i$. We can then formulate the problem as

[TABLE]

where we denote $\mathbf{b}_{j}:=B_{\cdot,j}\in\mathbb{R}^{c}$. Therefore, $\mathbf{z}_{i\rightarrow a}$ can be read as “a copy of $\mathbf{b}_{i}$ for the comparison with $\mathbf{b}_{a}$”, where $a$ is an adjacent node of $i$.

Using ADMM this can be optimized with the following sequence:

[TABLE]

where $k$ denotes the current iteration; $\mathbf{u}_{i\rightarrow j}$ denotes the scaled dual variables for $\mathbf{z}_{i\rightarrow j}$; $Z^{k+1}_{\cdot\rightarrow\cdot}$ denotes the set of variables $\{\mathbf{z}_{i\rightarrow j}\,|\,i\in\{1,\ldots d\},j\in\{1,\ldots,d_{i}\}\}$.333Note that our formulation is a scaled ADMM, see e.g. (Boyd et al., 2011) in Section 3.1.1. Therefore we have $\mathbf{u}_{i\rightarrow j}^{k+1}:=\mathbf{u}_{i\rightarrow j}^{k}+\mathbf{z}_{i\rightarrow j}^{k+1}-\mathbf{b}_{i}^{k+1}$ and not $\mathbf{u}_{i\rightarrow j}^{k+1}:=\mathbf{u}_{i\rightarrow j}^{k}+\rho(\mathbf{z}_{i\rightarrow j}^{k+1}-\mathbf{b}_{i}^{k+1})$. Furthermore, we defined

[TABLE]

Update of primal variables

The update of $B$ can be solved with an approximate gradient Newton method (Byrd et al., 1995), where fast calculation of the gradient and function evaluation is key to an efficient implementation.

Due to the second term, the calculation of $g(B,\boldsymbol{\beta}_{0})$ is in $O(d\cdot\max_{i}d_{i})$, and therefore, $O(d^{2})$ for dense graphs. However, the double sum in the second term can be expressed as follows:444Details in Appendix A.

[TABLE]

where we defined $\mathbf{q}_{ij}:=\mathbf{z}_{i\rightarrow e_{i}(j)}^{k}+\mathbf{u}_{i\rightarrow e_{i}(j)}^{k}$, and $q_{all}:=\sum_{i=1}^{d}\sum_{j=1}^{d_{i}}\mathbf{q}_{ij}^{T}\mathbf{q}_{ij}$, and $\mathbf{q}\uparrow_{i}:=\sum_{j=1}^{d_{i}}\mathbf{q}_{ij}$. Since $q_{all}$ and $\mathbf{q}\uparrow_{i}$ can be precalculated, each repeated calculation of $g(B,\boldsymbol{\beta}_{0})$ is in $O(d)$.

Furthermore, we note that

[TABLE]

Update of auxiliary variables

The update of $\mathbf{z}_{i\rightarrow j}$ and $\mathbf{z}_{j\rightarrow i}$, for each unordered pair $\{i,j\}$ can be performed independently, i.e.:

[TABLE]

This optimization problem has a closed form solution, which was proven in a different context in (Hallac et al., 2015), with

[TABLE]

where

[TABLE]

with $h:=||\mathbf{b}_{i}^{k+1}-\mathbf{u}_{i\rightarrow j}^{k}-\mathbf{b}_{j}^{k+1}+\mathbf{u}_{j\rightarrow i}^{k}||_{2}$.

2.2 Identifying the covariate clusters

Although in theory the optimization problem ensures that certain columns of $B$ are exactly $\mathbf{0}$, due to the use of ADMM this is not true in practice.555Details of the stopping criteria for ADMM are provided in Appendix A. We therefore follow the strategy proposed in (Chi and Lange, 2015): after convergence of the ADMM, we investigate $Z$, and place an edge between node $i$ and $j$, iff $\mathbf{z}_{i\rightarrow j}$ equals $\mathbf{z}_{j\rightarrow i}$. Note that we can test this equality exactly (without numerical difficulties) due to the thresholding of $\theta$ to 0.5 in Equation (2).

In the resulting graph, there are two possible ways to identify clusters:

• •

identify the connected components as clusters.666A component is in general not fully connected.

• •

consider only fully connected components as clusters.

Of course, in theory, since the optimization problem is convex, after complete convergence, we must have that the two ways result into the same clustering. However, in practice, we find that the latter leads to too many covariates not being clustered (i.e. each covariate is in a single cluster). The latter is also computationally difficult, since identifying the largest fully connected component is NP-hard. Therefore, we proceed here with the former strategy.

We denote the identified clusters as $\mathcal{C}_{1},\ldots,\mathcal{C}_{m}$, where $m$ is the number of clusters and $\mathscr{C}:=\{\mathcal{C}_{1},\ldots,\mathcal{C}_{m}\}$ is a partition of the set of covariates.

3 Approximate Bayesian Model Selection

Note that different hyper-parameter settings for $\nu$ will result in different clusterings. For our experiments, we consider $\nu\in\{n\cdot 2^{-0.1a}\,|\,a\in 0,1,2,\ldots,299\}$. This way, we get a range of clusterings. Like convex clustering, this allows us to plot a clustering hierarchy. However, in most situations, we are interested in finding the most plausible clustering, or ranking the clusterings according to some criterion.

One obvious way is to use cross-validation on the training data to estimate held-out classification accuracy for each setting of $\nu$. However, the computational costs are enormous. Another, more subtle issue is that the group lasso terms in Equation (2) jointly perform clustering and shrinkage, but controlled by only one parameter $\nu$. However, the optimal clustering and the optimal shrinkage might not be achieved by the same value of $\nu$, an issue well known for the lasso penalty (Meinshausen, 2007). Therefore, we consider here a different model selection strategy: an approximate Bayesian model selection with low computational costs.

After we have found a clustering, we train a new logistic regression classifier with the parameter space for $B$ limited to the clustering. Assuming some prior over $B$, the marginal likelihood $p(\mathbf{y})$ provides a trade-off between model fit and number of parameters (= number of clusters). A popular choice is the Bayesian information criterion (BIC) (Schwarz, 1978), which assumes that the prior can be asymptotically ignored777That is, the prior $\log p(B)$ does not increase with $n$.. However, BIC requires that the unpenalized maximum-likelihood estimate is defined. Unfortunately, this is not the case for logistic regression where linearly separable data will lead to infinite weights. For this reason, we suggest here to use a Gaussian prior on $B$ and then estimate the marginal likelihood with a Laplace approximation.

In detail, in order to evaluate the quality of a clustering, we suggest to use the marginal likelihood $p(\mathbf{y}|X,\mathscr{C},\boldsymbol{\beta}_{0})$, where we treat the intercept vector $\boldsymbol{\beta}_{0}$ as hyper-parameter. $X\in\mathbb{R}^{n\times d}$ denotes the design matrix, and $\mathbf{y}\in\{1,\ldots,c\}^{n}$ the responses.

We define the following Bayesian model

[TABLE]

where $\pi$ denotes our prior on $B$, and $\mathbf{x}^{\mathscr{C}}$ denotes the projection of covariates $\mathbf{x}$ on the clustering defined by $\mathscr{C}$. This means

[TABLE]

where we define the matrix $T^{\mathscr{C}}\in\mathbb{R}^{m\times d}$ as follows

[TABLE]

Assuming a Gaussian prior $N(0,\sigma^{2})$ on each entry of $B\in\mathbb{R}^{m\times m}$, the log joint distribution for one data point $(y,\mathbf{x})$ is given by

[TABLE]

For simplicity, let us denote $\mathbf{r}_{i}:=(B_{i,\cdot})^{T}$. Note that $\mathbf{r}_{i}\in\mathbb{R}^{m}$.

Let us denote by $S(B)$ the Hessian of the log joint probability of the observed data $X=\{\mathbf{x}_{1},\ldots,\mathbf{x}_{n}\}$ and prior $B$, i.e.

[TABLE]

where $S(B)_{ij}\in\mathbb{R}^{m\times m}$ denotes the $i,j$-th block of $S(B)\in\mathbb{R}^{cm\times cm}$. By the Bayesian central limit theorem, we then have that the posterior is approximately normal with covariance matrix $-S(\hat{B})^{-1}$. Then applying the Laplace approximation (see e.g. Ando (2010)), we get the following approximation for the marginal likelihood $p(X|\boldsymbol{\beta}_{0})$:

[TABLE]

where $\hat{B}$ and $\hat{\boldsymbol{\beta}}_{0}$ are the maximum a-posteriori (MAP) estimates from model (3). For large number of clusters $m$ the calculation of $|-S(\hat{B})^{-1}|$ is computationally expensive. Instead, we suggest to use only the diagonal of the Hessian.888The diagonal of the Hessian can be calculated for one sample $(y,\mathbf{x})$ as follows $\text{Diag}(\nabla_{\mathbf{r}_{i},\mathbf{r}_{j}}^{2}\{\log f(y|\mathbf{x},B,\boldsymbol{\beta}_{0})+\frac{1}{n}\log\pi(B)\})_{u}=(p_{i}-1)p_{i}\cdot x_{z}^{2}-\frac{1}{\sigma^{2}n}\,,$ with $i:=\lceil\frac{u}{m}\rceil$, and $z:=u-(i-1)m$.

We consider the intercept terms $\boldsymbol{\beta}_{0}$ as hyper-parameters and use empirical Bayes to estimate them, i.e. we set them to $\hat{\boldsymbol{\beta}}_{0}$. The hyper-parameter $\sigma$ cannot be estimated with empirical Bayes, since the unpenalized maximum-likelihood might not be defined, and this would lead to $\sigma\rightarrow\infty$. Therefore, we estimate $\sigma$ using cross-validation on the full model, i.e. no clustering. This is computationally feasible since the MAP is an ordinary logistic regression with $\ell_{2}$ penalty and cross-validation needs to be done only once and not for every clustering.

4 Synthetic Data Experiments

Here, in this section, we investigate the performance of the proposed covariate clustering method, as well as the proposed model selection criterion on synthetic data.

For all synthetic datasets, we set the number of classes $c$ to 4. The intercept vector $\boldsymbol{\beta}_{0}$ is set to the zero vector. We group the covariates evenly into 10 clusters. The weight vector $\mathbf{b}_{i}$ for covariate $i$ is set to the all zero vector except one position which is set to $0.5\cdot\text{ground truth cluster id}$. That means each covariate is associated with exactly one class. If two covariates $i$ and $j$ belong to the same cluster (ground truth), then $\mathbf{b}_{i}$ and $\mathbf{b}_{j}$ are equal. Finally, we generate samples from a multivariate normal distribution as follows: given a positive definite covariate similarity matrix $S^{*}$, we generate a sample $\mathbf{x}$ from class $y$ by $\mathbf{x}\sim N(B_{y,\cdot},S^{*})$. For each class we generate the same number of samples. We consider $d\in\{40,200,1000\}$ and $n\in\{40,400,4000\}$.

In practice, the a-priori similarity information between covariates will not fully agree with the associated classes. For that reason, we define $S^{*}$ such that the cluster structure implied by $S^{*}$ covers covariates that are associated to two different class labels. An example is shown in Figure 1.

Estimation of $S$

Since in practice, we only have access to a noisy estimate of $S^{*}$, rather than setting $S:=S^{*}$, we estimate $S$ from the data, by first normalizing the data to have zero mean for each class, and then pooling all samples from all classes for acquiring an estimate of the covariance matrix $\hat{\Sigma}$ using the Ledoit-Wolf shrinkage estimator (Ledoit and Wolf, 2004). Finally, we hard-threshold all negative entries in $\hat{\Sigma}$ to acquire $S$.

Baseline Methods

We compare the proposed method to $k$-means with the initialization proposed in (Arthur and Vassilvitskii, 2007) and convex clustering (Chi and Lange, 2015), which both use the same similarity matrix $S$ for clustering. For $k$-means, we consider $k\in\{1,2,\ldots,d\}$. For convex clustering the hyper-parameter $\gamma$ is tested in the range 0 to 5 evenly spaced with step size $1/300$, and weights $w_{ij}$ are set to the exponential kernel with $\phi=0.5$.999For $d>1000$ the convex clustering method led to memory overflow. Therefore, we needed to limit the weights to 5-nearest neighbors. After clustering with then train an ordinary multinomial logistic regression model. 101010On the other hand, first classification and then covariate clustering makes it difficult to incorporate the similarity information in $S$.

Clustering Evaluation

We repeat each experiment 10 times and report average and standard deviation of the adjusted normalized mutual information (ANMI) (Vinh et al., 2010) when compared to the true clustering. The ANMI score ranges from 0.0 (agreement with true clustering at pure chance level) to 1.0 (complete agreement with true clustering). All results are summarized in Table 5. For selecting a clustering we use the approximate marginal likelihood selection criterion described in Section 3, we also show the results when using cross-validation111111We use 5-fold cross-validation selecting the smallest model that is within one standard deviation from the best model, as suggested in (Hastie et al., 2009). (small font size), where it is computationally feasible.

Quality of marginal likelihood model selection criterion

Comparing the results when using cross-validation with the proposed approximate marginal likelihood criteria, we see that in most cases, there is no significant difference. However, compared to $k$-fold cross-validation, the approximate marginal likelihood criteria is $k$ times more computationally efficient.

Quality of clustering result

Even for small number of samples, $n\leq 400$, the proposed method’s clustering result is close to the ground truth, and fully recovers the correct clustering given sufficient number of samples, $n=4000$. This is neither the case for convex clustering, nor $k$-means clustering. Moreover, we find that for large number of variables, $d=1000$, convex clustering performed poorly. Finally, we note that in the more unrealistic scenario, where the a-priori similarity information fully agrees with the class association, convex clustering and $k$-means perform similar to the proposed method (see results in Appendix A).

4.1 Runtime Experiments

In order to check the efficiency of our proposed ADMM solution, we compare it to two standard solvers for convex optimization, namely ECOS (Domahidi et al., 2013) and SCS (O’Donoghue et al., 2016) using the CVXPY interface (Diamond and Boyd, 2016; Akshay Agrawal and Boyd, 2018). Note that SCS uses a generic ADMM method for which we use the same stopping criteria as our method. We run all experiments on a cluster with 88 Intel(R) Xeon(R) CPUs with 2.20GHz. We repeated each experiment 10 times and report the average runtime (wall-clock time) in in Table 2.

Our evaluation shows that the proposed ADMM scales well, both in terms of number of samples $n$, and in terms of number of variables $d$. Moreover, the proposed ADMM is considerably faster than the standard solvers when $n$ or $d$ are large, enabling us to apply the proposed method to real data.

5 Experiments on Real Data

For our experiments on real data, we used the movie review corpus IMDB (Maas et al., 2011), which consists of in total 100k (labeled and unlabeled) documents. IMDB is a balanced corpus with 50% of the movies being reviewed as “good movie”, and the remaining as “bad movie”. As the second dataset, we used the 20 Newsgroups corpus (Newsgroup20)131313http://people.csail.mit.edu/jrennie/20newsgroups/ with around 19k documents categorized into 20 classes (topics like “hockey”, “computer”,…). In order to check whether our model selection criterion correlates well with accuracy on held-out data, we use 10000 documents for training and clustering selection, and the remaining documents as held-out data.14141439k and 9k held-out documents for IMDB and Newsgroup20, respectively. We removed all duplicate documents, and performed tokenization and stemming with Senna (Collobert et al., 2011). Furthermore, we removed irrelevant covariates (details in Appendix A), leading to 2604 and 1538 covariates for IMDB and Newsgroup20, respectively.

5.1 Covariate Similarity Measure

From additional unlabeled documents, we determined the similarity between two covariates $i$ and $j$ as follows.

First, using the unlabeled datasets, we created for each covariate $i$ a word embedding $\mathbf{w}_{i}$. For IMDB, we created 50-dimensional word embeddings with word2vec (Mikolov et al., 2013) using 75k documents from the IMDB corpus.151515https://code.google.com/p/word2vec/ word2vec was used with the default settings. For Newsgroup20, since, the number of samples is rather small, we used the 300 dimensional word embeddings from GloVe (Pennington et al., 2014) that were trained on Wikipedia + Gigaword 5. Finally, the similarity between covariate $i$ and $j$ was calculated using $S_{ij}=e^{-\frac{1}{2}||\mathbf{w}_{i}-\mathbf{w}_{j}||_{2}^{2}}$. Using the similarity matrix $S$, the baseline methods are trained as in Section 4.

5.2 Quantitative Results

For real data, no ground truth for clustering is available. Therefore, we suggest to use the number of clusters vs the classification accuracy on held-out data. An ideal covariate clustering method should achieve high classification accuracy, while having only few clusters. Note that fewer number of clusters leads to more compact models, and potentially eases interpretability.

For selecting a clustering we use the proposed marginal likelihood criterion (Section 3), and the held-out data is only used for final evaluation. The results for IMDB and Newsgroup20 are shown in Table 3. We find that the marginal likelihood criterion tends to select models which accuracy on held-out data is similar to the full model, but with fewer number of parameters. Furthermore, for IMDB, we find that the proposed method leads to considerably more compact models than convex clustering and $k$-means, while having similar held-out accuracy. On the other hand, for Newsgroup20, the accuracies of the proposed method and $k$-means clustering are similar, indicating that there is good agreement between the similarity measure $S$ and the classification task.

In order to confirm that these conclusions are true independently of the model selection criterion, we also show in Table 4 the held-out accuracy of clusterings with number of clusters being around 100, 500 and 1000.161616Note that, in contrast to $k$-means, the proposed method and convex clustering can only control the number of clusters indirectly through their regularization hyper-parameter. The results confirm that the proposed method can lead to better covariate clusterings than $k$-means and convex clustering.

5.3 Qualitative Results

Like convex clustering, our proposed method can be used to create a hierarchical clustering structure by varying the hyper-parameter $\nu$. In particular, for large enough $\nu$, all covariates collapse to one cluster. On the other side, for $\nu=0$, each covariate is assigned to a separate cluster. Visualizations of the results from the proposed method, convex clustering and $k$-means are given in Appendix A. As expected, we observe that $k$-means and convex clustering tend to produce clusters that include covariates that are associated to different classes. For example, in IMDB, with $k$-means and convex clustering, all number ratings are clustered together. However, these have very different meaning, since for example 7/10, 8/10, 9/10, 10/10 rate good movies and 0, 1/10, *1/, 2/10,… rate bad movies. In contrast, our proposed method distinguishes them correctly by assigning them to different clusters associated with different classes.

6 Conclusions

We presented a new method for covariate clustering that uses an a-priori similarity measure between two covariates, and additionally, class label information. In contrast to $k$-means and convex clustering, the proposed method creates clusters that are a-posteriori plausible in the sense that they help to explain the observed data (samples with class label information). Like ordinary convex clustering (Chi and Lange, 2015), the proposed objective function is convex, and therefore insensitive to heuristics for initializations (as needed by $k$-means clustering).

Solving the convex objective function is computationally challenging. Therefore, we proposed an efficient ADMM algorithm which allows us to scale to several 1000 variables. Furthermore, in order to prevent computationally expensive cross-validation, we proposed a marginal likelihood criterion similar to BIC. For synthetic and real data, we confirmed the usefulness of the proposed clustering method and the marginal likelihood criterion.

Appendix A Appendices

A.1 Detailed Derivations of ADMM Update Equations

Update of primal variables

The double sum in the second term of $g(B,\boldsymbol{\beta}_{0})$ can be expressed as follows:

[TABLE]

where we defined $\mathbf{q}_{ij}:=\mathbf{z}_{i\rightarrow e_{i}(j)}^{k}+\mathbf{u}_{i\rightarrow e_{i}(j)}^{k}$, and $q_{all}:=\sum_{i=1}^{d}\sum_{j=1}^{d_{i}}\mathbf{q}_{ij}^{T}\mathbf{q}_{ij}$, and $\mathbf{q}\uparrow_{i}:=\sum_{j=1}^{d_{i}}\mathbf{q}_{ij}$. Since $q_{all}$ and $\mathbf{q}\uparrow_{i}$ can be precalculated, each repeated calculation of $g(B,\boldsymbol{\beta}_{0})$ is in $O(d)$.

Furthermore, we note that

[TABLE]

Update of auxiliary variables

The update of $\mathbf{z}_{i\rightarrow j}$ and $\mathbf{z}_{j\rightarrow i}$, for each unordered pair $\{i,j\}$ can be performed independently, i.e.:

[TABLE]

A.2 Stopping criteria for ADMM

As stopping criteria, we use, as suggested in (Boyd et al., 2011), the $\ell_{2}$ norm of the primal $\mathbf{r}^{k}$ and dual residual $\mathbf{s}^{k}$ at iteration $k$, defined here as

[TABLE]

We stop after iteration $k$, if

[TABLE]

where $l$ is the total number of edges in our graph, and $\epsilon$ is set to 0.00001. We also stop, in the cases where $k\geq 1000$.

A.3 Details of Approximate Bayesian Model Selection

The gradient with respect to $\mathbf{r}_{i}$ is

[TABLE]

with $p_{i}:=\frac{exp(\mathbf{r}_{i}^{T}\mathbf{x}+\beta_{i})}{\sum_{y^{\prime}}exp(\mathbf{r}_{y^{\prime}}^{T}\mathbf{x}+\beta_{y^{\prime}})}$.

The Hessian with respect to $\mathbf{r}_{i}$, $\mathbf{r}_{j}$ is

[TABLE]

A.4 Additional Experiments on Synthetic Data

$S^{*}$ agrees with class label information

If covariates $i$ and $j$ are in different clusters, then $S^{*}_{ij}=0$, otherwise $S^{*}_{ij}=0.9$. $S^{*}_{ii}$ is set to $1.0$. An example, where each cluster has four covariates is show in Figure 2.

A.5 Additional Details and Visualization of Real Data Results

A.5.1 Covariate Selection

First, for each dataset we extract the 10000 most frequent words as covariates, and represent each document $s$ by a vector $\mathbf{x}_{s}$ where each dimension contains the tf-idf score of each word in the document.171717We also perform $\ell_{2}$ scaling of each $\mathbf{x}_{s}$ which is known to improve classification performance. This is performed before the normalization of the covariates. Finally, we normalize the covariates to have mean 0 and variance 1.

For such high dimensional problems, many covariates are only negligibly relevant. Therefore, in order to further reduce the dimension, we apply a group lasso penalty on the columns of $B$ like in (Vincent and Hansen, 2014)181818We use also an ADMM algorithm to perform this optimization. We omit the details since the optimization function is the same as in Definition 1 in (Vincent and Hansen, 2014) with $\alpha=1$ and $\gamma\in\{2^{-a}\,|\,a\in 0,1,2,\ldots,29\}$. and select the model with the highest approximate marginal likelihood. This leaves us with 2604 and 1538 covariates for IMDB and Newsgroup20, respectively.

A.5.2 Qualitative Results

Like convex clustering, our proposed method can be used to create a hierarchical clustering structure by varying the hyper-parameter $\nu$. In particular, for large enough $\nu$, all covariates collapse to one cluster. On the other side, for $\nu=0$, each covariate is assigned to a separate cluster. In order to ease interpretability, we limit the visualization of the hierarchical clustering structure to all clusterings with $\nu\leq\nu_{0}$, where $\nu_{0}$ corresponds to the clustering that was selected with the marginal likelihood criterion. Furthermore, for the visualization of IMDB, we only show the covariates with odds ratios larger than 1.1. Part of the clustering results of the proposed method, convex clustering and $k$-means are shown in Figures 5, 6, 7, 8, and Figures 9, 10, for IMDB and Newsgroup20, respectively.191919For comparison, for $k$-means and convex clustering we choose (around) the same number of clusters as the proposed method.

Each node represents one cluster. If the cluster is a leaf node, we show the covariate, otherwise we show the size of the cluster (in the following called cluster description). From the logistic regression model of each clustering, we associate each cluster to the class with highest entry in weight matrix $B$, and then calculate the odds ratio to the second largest weight. The color of a node represents its associated class. All classes with coloring are explained in Figure 3 and 4 for IMDB and Newsgroup20, respectively. The odds ratio of each cluster is shown below the cluster description.

For IMDB, the qualitative differences between the proposed method and $k$-means are quite obvious. $k$-means clustering tends to produce clusters that include covariates that are associated to different classes. An example is shown in Figure 5, where all number ratings are clustered together. However, these have very different meaning, since for example 7/10, 8/10, 9/10, 10/10 rate good movies and 0, 1/10, *1/, 2/10,… rate bad movies. We observe a similar result for convex clustering as shown in Figure 6. In contrast, our proposed method distinguishes them correctly by assigning them to different clusters associated with different classes (see Figure 7 and 8).

For Newsgroup20, the qualitative differences of the proposed method with $k$-means are more subtle, but still we can identify some interesting differences. For example, $k$-means assigns “apple” and “cherry” to the same cluster which is associated with the class “Macintosh” (see Figure 9). In contrast, our proposed method correctly distinguishes these two concepts in the context of “Macintosh” (see Figure 10).