Discrete-Continuous ADMM for Transductive Inference in Higher-Order MRFs

TL;DR

This paper presents a new ADMM-based algorithm for transductive inference in higher-order MRFs that jointly optimizes continuous classifier parameters and discrete labels, ensuring global convergence and integrality.

Contribution

It introduces a decoupled optimization method combining discrete and continuous subproblems with ADMM, improving inference in higher-order MRFs over prior relaxations.

Findings

Effective in video object segmentation on DAVIS dataset

Preserves label integrality and guarantees convergence

Outperforms existing relaxation-based methods

Abstract

This paper introduces a novel algorithm for transductive inference in higher-order MRFs, where the unary energies are parameterized by a variable classifier. The considered task is posed as a joint optimization problem in the continuous classifier parameters and the discrete label variables. In contrast to prior approaches such as convex relaxations, we propose an advantageous decoupling of the objective function into discrete and continuous subproblems and a novel, efficient optimization method related to ADMM. This approach preserves integrality of the discrete label variables and guarantees global convergence to a critical point. We demonstrate the advantages of our approach in several experiments including video object segmentation on the DAVIS data set and interactive image segmentation.

Tables

Table 1. Table 1 : Comparison with the method of [ 23 ] on the SSL benchmark [ 9 ] . Reported are the average label-accuracy (in %) and variance over the splits. Our evaluation suggests that our method performs better for a standard hyperparameter setting except for three out of 20 settings. Moreover, it produces more consistent results, i.e. lower variances over the splits.

	Linear Kernel		RBF Kernel
Dataset	SDP	Ours	SDP	Ours
Digit1,10l	69.27$\pm$27.56	82.20$\pm$4.54	53.93$\pm$9.43	78.18$\pm$8.43
USPS,10l	57.72$\pm$13.73	64.58$\pm$3.37	40.10$\pm$11.67	48.19$\pm$5.84
BCI,10l	50.44$\pm$3.16	50.62$\pm$2.08	50.00$\pm$3.06	51.67$\pm$0.44
g241c,10l	49.88$\pm$38.92	55.42$\pm$3.95	62.33$\pm$36.84	89.98$\pm$0.32
g241n,10l	52.77$\pm$34.37	57.61$\pm$4.44	50.13$\pm$0.53	51.13$\pm$0.13
Digit1,100l	75.74$\pm$29.73	85.60$\pm$2.91	88.65$\pm$0.49	87.61$\pm$3.44
USPS,100l	63.44$\pm$9.97	72.14$\pm$0.84	39.83$\pm$12.63	56.54$\pm$3.31
BCI,100l	60.58$\pm$6.87	65.23$\pm$1.25	64.19$\pm$1.23	62.62$\pm$1.00
g241c,100l	64.92$\pm$17.47	86.31$\pm$0.91	85.63$\pm$0.76	89.34$\pm$1.07
g241n,100l	54.14$\pm$17.13	54.11$\pm$0.64	52.23$\pm$1.61	53.98$\pm$0.38

Equations

$$\underset{y\in \mathcal{Y}}{\min }\sum _{i\in \mathcal{V}}{E}_i(y_i)+\sum _{\begin{array}{c}C\in \mathcal{C}\\ |C|>1\end{array}}{E}_C(y_C).$$

$$\underset{\begin{array}{c}y\in \mathcal{Y},\\ W\in {\mathbb{R}}^{|\mathcal{L}|\times d^{\prime }}\end{array}}{\min }$$

$$W=\underset{W\in {\mathbb{R}}^{|\mathcal{L}|\times d^{\prime }}}{\mathrm{argmin}}\sum _{i\in \mathcal{V}}\ell (y_i;W\phi (x_i))+g(W).$$

$$g(W)=h(\Vert W{\Vert }_F),$$

$$f(\alpha ):=h(\Vert \Phi (X)\alpha {\Vert }_F),$$

$$\underset{\begin{array}{c}y\in \mathcal{Y},\\ \alpha \in {\mathbb{R}}^{|\mathcal{V}|\times |\mathcal{L}|}\end{array}}{\min }\sum _{i\in \mathcal{V}}\ell (y_i;K_i\alpha )+f(\alpha )+\sum _{C\in \mathcal{C}}{E}_C(y).$$

$$\underset{\begin{array}{c}y\in \mathcal{Y},\\ \alpha ,\beta \in {\mathbb{R}}^{|\mathcal{V}|\times |\mathcal{L}|}\end{array}}{\min }$$

$$K\alpha =\beta .$$

$$\underset{\begin{array}{c}y\in \mathcal{Y},\\ \alpha ,\beta \in {\mathbb{R}}^{|\mathcal{V}|\times |\mathcal{L}|}\end{array}}{\min }\underset{\lambda \in {\mathbb{R}}^{|\mathcal{V}|\times |\mathcal{L}|}}{\max }$$

$${\mathfrak{L}}_{\rho }(\alpha ,\beta ,\lambda ,y):=\sum _{i\in \mathcal{V}}\ell (y_i;{\beta }_i)+f(\alpha )$$

$$+\sum _{C\in \mathcal{C}}{E}_C(y)+⟨\lambda ,K\alpha -\beta ⟩+\frac{\rho }{2}\Vert K\alpha -\beta {\Vert }_F^2.$$

$$({\beta }^{t+1},y^{t+1})$$

$${\alpha }^{t+1}$$

$${\lambda }^{t+1}$$

$$\underset{{\beta }_i,y_i}{\mathrm{argmin}}\underset{{\psi }_i({\beta }_i,y_i;{\alpha }^t,{\lambda }_i^t)}{\underbrace{\ell (y_i;{\beta }_i)+\frac{\rho }{2}\Vert {\beta }_i-K_i{\alpha }^t-\text{nicefrac}{\lambda }_i^t\rho {\Vert }_F^2}},$$

$$B_{i,y_i}^{t+1}:=\underset{{\beta }_i}{\mathrm{argmin}}{\psi }_i({\beta }_i;y_i,{\alpha }^t,{\lambda }_i^t),$$

$$u_{i,y_i}^{t+1}:=\underset{{\beta }_i}{\min }{\psi }_i({\beta }_i;y_i,{\alpha }^t,{\lambda }_i^t).$$

$$y^{t+1}=\underset{y\in \mathcal{Y}}{\mathrm{argmin}}\sum _{i\in \mathcal{V}}{u}_{i,y_i}^{t+1}+\sum _{C\in \mathcal{C}}{E}_C(y).$$

$${\beta }_i^{t+1}=B_{i,y_i^{t+1}}^{t+1}.$$

$${\mathfrak{L}}_{\rho }({\alpha }^t,{\beta }^{t+1},{\lambda }^t,y^{t+1})-{\mathfrak{L}}_{\rho }({\alpha }^t,B_{:,y^t}^{t+1},{\lambda }^t,y^t)\le -\delta .$$

$${\alpha }^{t+1}=\underset{\alpha }{\mathrm{argmin}}\nu ⟨\alpha ,K\alpha ⟩+\frac{\rho }{2}\Vert K\alpha -{\beta }^{t+1}+\text{nicefrac}{\lambda }^t\rho {\Vert }_F^2.$$

$$\frac{L^2}{\rho {\sigma }_{\min }(K^{\top }K)}+\frac{m-\rho {\sigma }_{\min }(K^{\top }K)}{2}<0.$$

$$\begin{gathered} {\mathfrak{L}}_{\rho }({\alpha }^{t+1},{\beta }^{t+1},{\lambda }^{t+1},y^{t+1})-{\mathfrak{L}}_{\rho }({\alpha }^t,{\beta }^t,{\lambda }^t,y^t) \\ \le \left(\frac{L^2}{\rho {\sigma }_{\min }(K^{\top }K)}+\frac{m-\rho {\sigma }_{\min }(K^{\top }K)}{2}\right)\Vert {\alpha }^{t+1}-{\alpha }^t{\Vert }_F^2 \\ -\delta [[y^{t+1}\ne y^t]], \end{gathered}$$

$$\Vert {\alpha }^{t+1}-{\alpha }^t{\Vert }_F\to 0,$$

$$\Vert {\beta }^{t+1}-{\beta }^t{\Vert }_F\to 0,$$

$$\Vert {\lambda }^{t+1}-{\lambda }^t{\Vert }_F\to 0,$$

$$\Vert K{\alpha }^{t+1}-{\beta }^{t+1}{\Vert }_F\to 0.$$

$$0\in \partial (\ell (y_i^{\ast };\cdot ))({\beta }_i^{\ast })-{\lambda }_i^{\ast },$$

$$0\in \partial g({\alpha }^{\ast })+K^{\top }{\lambda }^{\ast }$$

$$K{\alpha }^{\ast }={\beta }^{\ast },$$

$${\alpha }^{\ast }=\underset{\alpha }{\mathrm{argmin}}\sum _{i\in \mathcal{V}}\ell (y_i^{\ast };K_i\alpha )+f(\alpha ).$$

$$f(y)\le f(x)+⟨\nabla f(x),y-x⟩+\frac{L}{2}\Vert x-y{\Vert }_F^2.$$

$$f(y)\ge f(x)+⟨\nabla f(x),y-x⟩-\frac{m}{2}\Vert x-y{\Vert }_F^2.$$

$${\mathfrak{L}}_{\rho }({\alpha }^{t+1},{\beta }^{t+1},{\lambda }^{t+1},y^{t+1})-{\mathfrak{L}}_{\rho }({\alpha }^t,B_{:,y^t}^t,{\lambda }^t,y^t)$$

$$={\mathfrak{L}}_{\rho }({\alpha }^t,B_{:,y^t}^{t+1},{\lambda }^t,y^t)-{\mathfrak{L}}_{\rho }({\alpha }^t,B_{:,y^t}^t,{\lambda }^t,y^t)$$

$$+{\mathfrak{L}}_{\rho }({\alpha }^t,{\beta }^{t+1},{\lambda }^t,y^{t+1})-{\mathfrak{L}}_{\rho }({\alpha }^t,B_{:,y^t}^{t+1},{\lambda }^t,y^t)$$

$$+{\mathfrak{L}}_{\rho }({\alpha }^{t+1},{\beta }^{t+1},{\lambda }^t,y^{t+1})-{\mathfrak{L}}_{\rho }({\alpha }^t,{\beta }^{t+1},{\lambda }^t,y^{t+1})$$

Full text

1 Introduction

Various problems in computer vision, computer graphics and machine learning can be formulated as MAP inference in a (possibly higher order) Markov random field (MRF) [34, 27, 42, 29, 31, 13, 33, 45]. The resulting optimization problem is defined over a hypergraph $(\mathcal{V},\mathcal{C})$ and a finite label set $\mathcal{L}$ as:

[TABLE]

The optimization variable $y\in\mathcal{Y}:=\mathcal{L}^{|\mathcal{V}|}$ corresponds to a labeling of the vertices $\mathcal{V}$ and assigns a label $y_{i}\in\mathcal{L}$ to each vertex $i\in\mathcal{V}$. For convenience, we make a distinction between the singleton clique energies (unaries) $E_{i}(y_{i})$ and the higher order energies $E_{C}(y_{C})$, $|C|>1$.

In computer vision tasks, where the image is interpreted as a (higher order) pixelgrid $(\mathcal{V},\mathcal{C})$, the higher order potentials often correspond to priors favoring spatially smooth solutions. In semantic image segmentation, for instance, inference in MRFs is widely used as a post-processing step to introduce spatial smoothness on the labeling $y$ [11]. In this sense, the overall task of semantic segmentation is subdivided into two tasks: First, a classifier, parameterized by $W$, is trained in a supervised fashion on a sufficiently large labeled training set, which assigns to each pixel in a test image a class (probability) score. Second, to enforce spatial smoothness of the labeling of the test image, MAP-inference in an MRF is performed in a post-processing, in which the class (probability) scores are interpreted as the unary energies $E_{i}(y_{i};W)$ for each $y_{i}$. We argue that it is advantageous to merge such a two-step approach into a joint approach, as the training of the classifier profits from both, (i) the distribution of the unlabeled pixels in the color or feature space, and (ii) the available structural information about the unlabeled pixels, namely the spatial smoothness prior. Conversely, the segmentation will also profit from an improved classifier.

A joint formulation has two interpretations; on the one hand, it is a semi-supervised learning method that makes use of structural knowledge about the training data to learn a classifier. Such knowledge may take the form of higher-order clique energies $E_{C}$ [59, 3] on the labeling and acts as weak supervision in the training process. This approach helps to mitigate the need for large amounts of annotated training data in typical modern machine learning applications. We, on the other hand, focus on its interpretation as a transductive inference method [58], i.e. the approach to directly infer the labels of specific test data given specific training data, which we accomplish by incorporating a variable classifier in the inference process. Transductive inference stands in contrast to inductive inference, which refers to first learning a general model from training data and subsequently applying the model to predict labels of a-priori unknown test data. We show the benefits of using transductive inference for the tasks of video object segmentation (cf. Fig. 1) as well as scribble-based segmentation (cf. Fig. 4).

1.1 Contributions

We propose a general joint model that assumes the unaries not be fixed for inference in the MRF, but rather optimized jointly with the labeling. Let for each $i\in\mathcal{V}$, $x_{i}\in\mathbb{R}^{d}$ denote the $d$-dimensional feature vector associated to the $i$th vertex and let $\varphi:\mathbb{R}^{d}\to\mathbb{R}^{d^{\prime}}$ be a feature map. Then, mathematically, such a task can be naturally formulated in terms of a bilevel optimization problem:

[TABLE]

Here, the upper-level task is inference in a MRF with additional unaries $E_{i}(y_{i};W):=\ell(y_{i};W\varphi(x_{i}))$, parameterized by linear classifier weights $W\in\mathbb{R}^{|\mathcal{L}|\times d^{\prime}}$, and a loss function $\ell\colon\mathcal{L}\times\mathbb{R}^{|\mathcal{L}|}\to\mathbb{R}$. The lower-level task associates to each given set of labels $y$ the optimal parameters $W$. For instance, if $\ell$ is the hinge loss and $g(W)=\|W\|_{F}^{2}$, then the lower-level optimization problem amounts to the training of a classical SVM. Note that in a semi-supervised learning context it might be more convenient to swap upper and lower level tasks, since the primary interest is the estimation of the classifier that minimizes the generalization error and not the inferred labeling. However, mathematically, both viewpoints are equivalent.

The model (2) suggests a simple alternating optimization scheme as in Lloyd’s algorithm [38] for $k$-means to compute a local optimum. However, such an approach has two major drawbacks: (i) The lower-lever subproblems are expensive, which is prohibitive for large scale applications. (ii) The optimization is prone to poor local optima and therefore sensitive to initialization [63]. Motivated by the good practical performance of the alternating direction method of multipliers (ADMM) in nonconvex optimization, we propose to generalize vanilla ADMM (commonly applied in continuous optimization) to discrete-continuous problems of the form (2), while preserving integrality of the discrete variables. Since our method serves as a general algorithmic framework to tackle such problems, it is also relevant to semi-supervised and transductive learning in a broader sense.

The main contributions of this work can be summarized as follows:

• •

We devise a decomposition of the model into simple, purely discrete and purely continuous subproblems within the framework of proximal splitting. The subproblems can be solved in a distributed fashion.

• •

We devise a tailored ADMM-inspired algorithm, discrete-continuous ADMM, to compute a local optimum of (2). In contrast to vanilla ADMM, our algorithm allows us to obtain sub-optimal solutions of the MAP inference problem so that also computationally more challenging MRFs can be considered.

• •

We generalize the convergence of nonconvex vanilla ADMM to the presented inexact discrete-continuous ADMM.

• •

In diverse experiments we demonstrate the relevance and generality of our model and the efficiency of our method: In contrast to standard $k$-means, our model integrates well with deep features. In contrast to a tailored SDP relaxation approach for transductive logistic regression, our method produces more consistent results, while being more efficient in terms of both runtime and memory consumption.

1.2 Related work

To improve image segmentation results it is common practice to treat the unary terms $E_{i}$ as additional variables in the optimization [5, 64, 8, 49, 55, 56, 57]. More recently, [55, 56] revealed the equivalence of $k$-means clustering with pairwise constraints [59, 3] and the Chan-Vese [8] approach, where the average foreground and background intensities (corresponding to the centroids in $k$-means) are not assumed to be fixed, but are rather treated as additional variables. The goal of this approach is to jointly cluster the pixels in the color space and regularize the cluster-assignment (the segmentation) in the image space. The clustering viewpoint suggests the application of the “kernel trick”, which allows us to separate more complicated, possibly nonlinearly deformed color clusters [59, 3, 55, 56]. Experimentally, it has been shown that this approach integrates well with color or even depth pixel features [55, 56]. Due to the enormous success of deep convolutional neural networks on computer vision tasks, it is tempting to replace the color features by more sophisticated deep features that are capable of compactly representing complicated semantic information [32, 11]. However, high-dimensional deep features are in general not “$k$-means friendly” [61] and without further preprocessing of the features as in [61] the plain Chan-Vese $k$-means approach does not generalize very well to deep features, despite of (almost) linear separability of the data. Since deep neural network classifiers can be viewed as a linear model on top of a deep feature extractor, we propose to alter the approaches from related work by using a (multiclass) SVM or a (multinomial) logistic regression model along with deep features. Under the absence of general higher order terms (i.e. $E_{C}=0$, for all $C\in\mathcal{C}$) our model (2) is closely related to transductive SVMs [58, 22, 4] and transductive logistic regression [23]. In such a setting, optimization schemes that alternately optimize w.r.t. to labels and model parameters as in Lloyd’s algorithm [38] are ineffective [63]. Other approaches, such as SDP relaxations are computationally expensive [23].

Instead, we propose an algorithm related to ADMM, which has recently been successfully applied to many nonconvex continuous optimization problems [10, 44, 52, 40, 35]. ADMM appears similar in form to message passing and subgradient descent schemes applied to the Lagrangian dual problem (dual decomposition) [30, 6, 39, 62, 53]. The latter is a Lagrangian relaxation approach, so that in difficult nonconvex cases the linear equality constraints may remain violated in the limit [30]. In contrast, ADMM attempts to solve the problem exactly and enforce the linear equality constraints strictly via additional quadratic penalty terms. In order to make mixed discrete-continuous problems such as (2) amenable to ADMM, related approaches often relax the discrete variable and perform rounding operations [21, 54]. In contrast, we propose a generalization of vanilla ADMM that preserves the integrality of the label variable and admits a theoretical convergence guarantee under affordable conditions. In the traditional convex and continuous setting, ADMM [19, 18] converges under mild conditions [17, 14]. For more restrictive nonconvex problems, its convergence has only been established recently [20, 37]. In this case, however, the required assumptions are fairly strong.

2 Discrete-Continuous ADMM

The coupling of the discrete labeling variable $y$ and the continuous variable $W$ renders problem (2) hard to solve. This is not surprising since the related $k$-means clustering problem is known to be np-hard. A common approach is to compute a local optimum by a simple discrete-continuous coordinate descent approach as in Lloyd’s algorithm [38]. Instead, we propose an advantageous decoupling into purely discrete and purely continuous subproblems, which allows us to compute a local optimum by updating the continuous and discrete variables jointly and efficiently.

2.1 Variable decoupling via ADMM

To this end, we employ a change of representation to make the proposed problem amenable to the “kernel trick”. Note that, for any fixed labeling $y$, the lower-level task in (2) amounts to supervised SVM training (resp. supervised logistic regression). Thus, we can apply the representer theorem [50]: Let $\Phi(X)$ be the feature matrix for a (possibly infinte-dimensional) matrix feature map $\Phi:\mathbb{R}^{d\times|\mathcal{V}|}\to\mathbb{R}^{d^{\prime}\times|\mathcal{V}|}$ and let

[TABLE]

for $h\colon[0,\infty)\to\mathbb{R}$ strictly monotonically increasing. Then, the weights $W^{\top}=\Phi(X)\alpha$ can be substituted via their representation $\alpha\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{L}|}$ in terms of the features. More precisely, we replace the scalar products $W\varphi(x_{i})$, up to transposition, by $K_{i}\alpha=(W\varphi(x_{i}))^{\top}$ where $K:=\Phi(X)^{\top}\Phi(X)$ denotes the Gram or kernel matrix.

For $f:\mathbb{R}^{|\mathcal{V}|\times|\mathcal{L}|}\to\mathbb{R}$, being defined as

[TABLE]

this substitution leaves us with the following equivalent mixed integer nonlinear program formulation of (2):

[TABLE]

In order to decompose problem (5) into simple subproblems associated with each $i\in\mathcal{V}$, we introduce auxiliary variables $\beta_{i}=K_{i}\alpha$, which yields

[TABLE]

Note that the objective of (6) is a separable function over the $\beta_{i}$. This suggests to relax the linear constraint $K\alpha=\beta$ and consider the equivalent saddle point problem:

[TABLE]

where $\lambda\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{L}|}$ are the Lagrange multipliers corresponding to $K\alpha=\beta$ and $\mathfrak{L}_{\rho}$ denotes the “discrete-continuous” augmented Lagrangian, that for some penalty parameter $\rho>0$ is defined as

[TABLE]

We show in Sec. 2.2 that, for fixed $\lambda$ and $\alpha$, the function $\mathfrak{L}_{\rho}(\alpha,\cdot,\lambda,\cdot)$ can be minimized (not necessarily to global optimality) efficiently and jointly over $y$ and $\beta$. This central observation and the good practical performance of ADMM in nonconvex optimization motivates the following generalization to discrete-continuous problems of the form (7).

We propose an algorithm that, similar to continuous ADMM, updates the discrete-continuous variable-pair $(\beta^{t+1},y^{t+1})$ via joint (and possibly suboptimal) minimization of $\mathfrak{L}_{\rho}(\alpha^{t},\cdot,\lambda^{t},\cdot)$. Subsequently, it updates $\alpha^{t+1}$ via minimization of $\mathfrak{L}_{\rho}(\cdot,\beta^{t+1},\lambda^{t},y^{t+1})$ and the Lagrange multiplier $\lambda^{t}$ by performing one iteration of gradient ascent on $\mathfrak{L}_{\rho}(\alpha^{t+1},\beta^{t+1},\cdot,y^{t+1})$ with step size $\rho>0$. In summary, the update steps at iteration $t$ are given as

[TABLE]

In practice, we choose the step size adaptively, as this often leads to better solutions in terms of objective value: For finitely many iterations, the penalty parameter $\rho$ is increased according to the schedule $\rho_{t+1}=\min\;\{\rho_{\max},\tau\rho_{t}\}$ with $\tau>1$ and some $\rho_{\max}>0$ that guarantees theoretical convergence of the algorithm (cf. Sec. 3).

2.2 Distributed solution of the subproblems

In this section, we describe the implementation of update steps (9)–(11) in our algorithm. In principle, (9) could be solved by minimization over $\beta$ for every feasible labeling $y\in\mathcal{Y}$. Obviously, this is not a viable approach, as it implies performing exhaustive search over the set $\mathcal{Y}$, which has size $\lvert\mathcal{L}\rvert^{\lvert\mathcal{V}\rvert}$. Instead, we pursue the following more efficient strategy.

Solution via lookup-tables.

Assume first the absence of any higher order energies, i.e. $E_{C}=0$, for all $C\in\mathcal{C}$. Then, since $\mathfrak{L}_{\rho}(\alpha^{t},\beta,\lambda^{t},y)$ is separable w.r.t. $\beta_{i}$ and $y_{i}$, we can decompose problem (9) into $|\mathcal{V}|$ independent problems of the form

[TABLE]

which can thus be solved in parallel. In the presence of higher-order energies, however, the problems (12) are not completely independent, because the variables $y_{i}$ are coupled via the energies $E_{C}$ in which they appear. In this case, we first solve (12) w.r.t. only the continuous variables $\beta_{i}$ for every possible label $y_{i}\in\mathcal{L}$ and store the results in a lookup-table $(u^{t+1},B^{t+1})$.

Precisely, for each $1\leq i\leq|\mathcal{V}|$ and each $y_{i}\in\mathcal{L}$ we create an entry $(u_{i,y_{i}}^{t+1},B_{i,y_{i}}^{t+1})$ according to

[TABLE]

In a second step, we determine the discrete variable update $y^{t+1}$ as the (possibly suboptimal) solution of the MRF

[TABLE]

Afterwards, the continuous variable updates $\beta_{i}^{t+1}$ can be read off from the solution of (13) via

[TABLE]

Note that there is an abundance of algorithms available to tackle problems of the form (14) such as graph cuts for binary submodular MRFs [29], move making and message passing algorithms [13, 28], primal-dual algorithms [31, 15] and more. For an overview, see also [24].

The matrix $u$ specifies the unary energies in problem (14) that pushes the MRF to attain a labeling which corresponds to a more suitable classifier. The latter is determined by a tradeoff between minimizing the distance of $\beta_{i}$ to the current consensus parameters $K_{i}\alpha^{t}+\nicefrac{{\lambda_{i}^{t}}}{{\rho}}$ and minimizing the loss term corresponding to sample $i$.

In case of suboptimality of $y^{t+1}$ we require that $y^{t+1}$, for some $\delta\geq 0$, satisfies a (sufficient) descent condition

[TABLE]

If this condition is violated, then we keep the previous iterate $y^{t+1}=y^{t}$. Under condition (16), the overall convergence of our algorithm is guaranteed (cf. Prop. 1 and Prop. 2 in Sec. 3). We summarize our method in Alg. 1.

Note that if the discrete subproblem (14) is solved to global optimality, our method specializes to classical nonconvex ADMM applied to a purely continuous problem $\min_{\alpha}E(K\alpha)+f(\alpha)$. The function $E(\beta)=\min_{y\in\mathcal{Y}}\sum_{i\in\mathcal{V}}\ell(y_{i};\beta)+\sum_{C\in\mathcal{C}}E_{C}(y)$ encapsulates the minimization over the discrete labelings $y$. This results in a pointwise minimum over exponentially many functions.

Distributed optimization.

Distributed optimization is considered one of the main advantages of ADMM in supervised learning [16, 6]. In our method, the $(\beta,y)$ update requires the solution of only $\lvert\mathcal{L}\rvert\cdot|\mathcal{V}|$ many (instead of $\lvert\mathcal{L}\rvert^{\lvert\mathcal{V}\rvert}$ for the naive approach) independent and small-scale continuous minimization problems of the form (13) and one additional discrete problem (14). This suggests the distributed solution of the subproblems (13), for instance on a GPU. Subsequently, the optimization of the MRF (14) and the update of the consensus variable is carried out after gathering the solutions of the subproblems. Since (14) need not be solved to optimality, the MRF solver may be stopped early to speed up computation. This is particularly useful if a primal-dual algorithm for solving the LP-relaxation is used [31, 15].

Exploit duality.

If the loss terms $\ell(y_{i};\cdot)$ are convex and lower semicontinuous, then the independent subproblems (13) can be solved efficiently via duality as follows. For all loss functions we consider, it is convenient to solve the dual problem as it scales linearly with the number of training samples (which is equal to one in our case). For the Crammer and Singer multiclass SVM loss [12], for instance, there exists an efficient variable fixing algorithm [26] for solving the dual problem. For the softmax loss the dual problem reduces to a one-dimensional nonlinear equation via the Lambert-$W$ function [36] and may be solved by performing a few iterations of Newton’s or Halley’s method. For the special case of the one-vs.-all hinge loss, (13) can be solved in closed form. In any case, each subproblem involves only a small number of instructions, which is important for a GPU-based implementation.

2.3 Consensus update

For a quadratic regularizer $h(x)=\nu x^{2}$, where $\nu$ is the regularization parameter, the update step (10) is equivalent to

[TABLE]

This is a quadratic problem that can be solved via a normal equation using either a cached eigenvalue decomposition of the kernel matrix, or an iterative algorithm such as conjugate gradient (CG). The latter is preferred for large scale applications, as each CG iteration involves a kernel-matrix-vector multiplication $Kv$. For the linear kernel $K=X^{\top}X$, this guarantees efficiency of our method, since $K$ does not have to be stored explicitly. For general kernels such as the RBF kernel, a low rank approximation to the kernel matrix $K\approx GG^{\top}$, for some $G\in\mathbb{R}^{|\mathcal{V}|\times l}$ with $l\ll|\mathcal{V}|$ can be obtained for instance via the Nyström method [43, 60] or random features [47]. Furthermore, in practice, often only a small number of conjugate gradient iterations are necessary.

3 Convergence analysis

In this section, we provide a complete convergence analysis of the proposed algorithm. To this end, we make the following assumptions:

• •

The function $f$ is $L$-smooth, $m$-semiconvex and lower-bounded, i.e. $f$ is differentiable and $\nabla f$ is Lipschitz-continuous with modulus $L$ and there exists $m>0$ sufficiently large so that $f+\frac{m}{2}\|\cdot\|_{F}^{2}$ is convex.

• •

For all $y_{i}\in\mathcal{L}$, $\ell(y_{i};\cdot)$ is lower-bounded.

• •

The kernel matrix $K\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{V}|}$ is surjective, i.e. the smallest eigenvalue $\sigma_{\min}(K^{\top}K)>0$ is positive.

• •

After finitely many iterations $t$ the penalty parameter $\rho$ is sufficiently large and kept fixed such that

[TABLE]

When the MRF subproblem is solved to global optimality, convergence can be guaranteed by considering a pointwise minimum over exponentially many augmented Lagrangians and applying existing theory [37, 20]. For the general case, however, the theory needs to be extended. Our convergence proof borrows arguments from [37, 20], where the convergence of ADMM in the nonconvex setting is shown via a monotonic decrease of the augmented Lagrangian. In our case, for a sufficiently large penalty parameter $\rho$, we can achieve a monotonic decrease of the “discrete-continuous” augmented Lagrangian (8), even if the MRF subproblem (14) is not solved globally optimal. This allows us to stop exact MRF solvers early or to apply heuristic solvers if computing global optima is intractable.

For the complete proofs of all the theoretical results, presented in this section, cf. the Appendix A.

Lemma 1**.**

Let $K\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{V}|}$ be surjective and $\delta\geq 0$. For $\rho$ meeting condition (18) we have that

The “discrete-continuous” augmented Lagrangian (8) decreases monotonically with the iterates $(\alpha^{t},\beta^{t},\lambda^{t},y^{t})$:

[TABLE]

where $\llbracket\cdot\rrbracket$ denotes the Iverson bracket. 2. 2.

$\{\mathfrak{L}_{\rho}(\alpha^{t+1},\beta^{t+1},\lambda^{t+1},y^{t+1})\}_{t\in\mathbb{N}}$* is lower bounded.* 3. 3.

$\{\mathfrak{L}_{\rho}(\alpha^{t+1},\beta^{t+1},\lambda^{t+1},y^{t+1})\}_{t\in\mathbb{N}}$* converges.*

We are now able to guarantee that feasibility is achieved in the limit. This is in contrast to a dual-decomposition approach [30, 62, 53] or a Gauss-Seidel quadratic penalty method (with finite penalty parameter $\rho$), used for instance in [51, 25], where a violation of the consensus constraint remains in the limit. Moreover, if $\delta>0$ is chosen strictly positive, then the discrete variable is guaranteed to converge, i.e. for $T$ sufficiently large, we have $y^{t+1}=y^{t}$ for all $t>T$.

Lemma 2**.**

Let $\{(\alpha^{t},\beta^{t},\lambda^{t},y^{t})\}_{t\in\mathbb{N}}$ be the iterates produced by Alg. 1. Then $\{(\alpha^{t},\beta^{t},\lambda^{t},y^{t})\}_{t\in\mathbb{N}}$ is a bounded sequence. Furthermore, for $t\to\infty$ the distance of two consecutive continuous iterates vanishes, and feasibility is achieved in the limit:

[TABLE]

Finally, if $\delta>0$ is chosen strictly positive, then there exists some $T\in\mathbb{N}$ such $y^{t+1}=y^{t}$ for all $t>T$.

The limit points of our algorithm correspond to “discrete-continuous” critical points of the augmented Lagrangian.

Definition 1** (“Discrete-continuous” critical point).**

We call $(\alpha^{*},\beta^{*},\lambda^{*},y^{*})$ a “discrete-continuous” critical point of the “discrete-continuous” augmented Lagrangian (8) if it satisfies

[TABLE]

for $y^{*}$ with $E_{C}(y^{*})<\infty$ for all $C\in\mathcal{C}$. Here, $\partial f(x)$ denotes the “limiting” subdifferential [48, Definition 8.3] of the function $f$ at $x$ with $f(x)<\infty$.

Proposition 1**.**

Let $\delta\geq 0$. Then any limit point $(\alpha^{*},\beta^{*},\lambda^{*},y^{*})$ of the sequence $\{(\alpha^{t},\beta^{t},\lambda^{t},y^{t})\}_{t\in\mathbb{N}}$ is a “discrete-continuous” critical point.

Finally, under convexity of $f$ and $\ell(y_{i};\cdot)$, for all $y_{i}\in\mathcal{L}$ and strictly positive $\delta>0$, we can guarantee that the sequence of iterates produced by Alg. 1 globally converges to a point $(\alpha^{*},\beta^{*},\lambda^{*},y^{*})$ which has the following property: $\alpha^{*}$ is the global optimum of the supervised learning problem w.r.t. the estimated training labels $y^{*}$:

[TABLE]

Proposition 2**.**

Let $\ell(y_{i};\cdot)$ and $g$ be proper, convex and lower-semicontinuous and let $\delta>0$. Then the sequence $\{(\alpha^{t},\beta^{t},\lambda^{t},y^{t})\}_{t\in\mathbb{N}}$ produced by Alg. 1 converges to a “discrete-continuous” critical point $(\alpha^{*},\beta^{*},\lambda^{*},y^{*})$ of (8) and $\alpha^{*}$ solves the problem (27) to global optimality.

Discussion of the assumptions.

Note that in general the kernel matrix $K$ is not surjective. However, for the strictly positive definite RBF kernel, $K$ is strictly positive definite so that convergence can be achieved for finite $\rho$. In order to enforce theoretical convergence for general kernels, we may add a small constant to the diagonal of the kernel matrix $K:=K+\gamma I$ that alters the model only slightly. In fact, for the binary SVM, this change is equivalent to replacing the hinge loss with its square.

4 Experiments

In this section, we present the experimental results of our method on several transductive learning tasks. First, we compare our method to an SDP relaxation method for transductive multinomial logistic regression by [23]. Second, we use our model and solver for the tasks of object video segmentation as well as image segmentation with user interaction, showing improvements on the false positive rate of object pixels.

4.1 Comparison with SDP relaxation for transductive learning

In this experiment, we consider the standard SSL benchmark [9] for a comparison with the SDP relaxation method for transductive multinomial logistic regression by [23]. The benchmark is a collection of several datasets, with varying feature dimensions and number of classes. Each dataset is provided with 12 splits into $l=10$ or $l=100$ labeled and $N-l$ unlabeled samples. We introduce additional unary energies $E_{C}$ with $|C|=1$ for all the labeled examples, to constrain their label to be fixed during optimization. While [23] incorporates an entropy prior on the labeling which favors an equal balance distribution, we introduce a higher order potential $E_{C}$, with $C=\mathcal{V}$, that restricts the solution to deviate at most 10 percent from the equal balance distribution. We solve the LP-relaxation of the higher-order MRF subproblem (14) with the dual-simplex method and round the solution. The baseline results are computed with a MATLAB implementation that is provided by the authors. For these experiments, we use the softmax loss and set the regularization parameter $\nu=0.05$ for the linear kernel. For the RBF kernel we manually chose the variance parameter $\sigma=0.5477$ and the regularization parameter $\nu=0.0025$. We chose the initial penalty parameter $\rho_{0}=0.001$ and $\tau=1.003$. All values are averaged over 12 different splits. The evaluation in Tab. 1 suggests, that our method performs better for a standard hyperparameter setting except for three out of 20 settings. Moreover, it produces more consistent results, i.e. lower variances over the splits, which suggests that our method is more robust towards noise and poorly labeled data.

4.2 Video object segmentation

In this experiment, we evaluated our method on video object segmentation. Here, the task is to segment an object throughout a video, given its mask in the first frame. This problem has been successfully approached by [7], using end-to-end deep learning with fully convolutional neural networks. At test time, their classifier is fine-tuned on the appearance of the object and the background in the first frame and predicts the object pixels of individual later frames. However, this method struggles with drastic appearance changes of the object, which have not been learned in advance. These include pose changes, sharp lighting and background changes or severe occlusions as shown in Fig. 2.

We propose to use a transductive approach instead. More precisely, we use the pre-trained (not fine-tuned) OSVOS parent network [7] as a deep feature extractor and a MRF model with a variable classifier in the form of (2). We use a simple linear kernel SVM in our model, as the extracted deep features are almost linearly separable. Further, we introduce unary indicator energies $E_{i}$ to fix the labels of the user-annotated pixels in the first frame and pairwise energies $E_{ij}$ for adjacent pixels in any frame to favor spatially smooth solutions. Similar to [7], we do not use any temporal consistency terms. To reduce the number of examples, we apply our method on a superpixel level and extract 6000 super-pixels [1] for each frame and apply average pooling over the superpixels. We compare the proposed transductive approach to OSVOS [7] and the classical (inductive) MRF inference approach (where the classifier is learned with the first frame only) on the DAVIS benchmark [46]. For both the inductive and the transductive approach the used linear kernel SVM model, the higher order energies $E_{C}$ and the extracted superpixels are the same. The results are shown in Fig. 2. It can be seen that [7] works well as long as the appearance of the object and background are sufficiently similar to the first frame (first column). In frames 60 to 68, where the object is occluded, both the inductive MRF inference approach and OSVOS produce a large number of false positive object pixels. In this experiment the intersection-over-union scores (the higher the better) are 0.7087 for our method, 0.6452 for OSVOS and 0.5063 for the inductive approach. Similarly in the car-shadow sequence, OSVOS and the inductive approach mask additionally the other car and the motorbike in frame 39 (cf. Fig. 3). In contrast our method masks the correct car only. Here, the intersection-over-union scores are 0.9262 for OSVOS, 0.9196 for our method and 0.8844 for the inductive approach.

4.3 Image segmentation with user interaction

We evaluated our method on the task of interactive foreground-background segmentation with deep features. Like in the previous experiment we used OSVOS as a deep feature extractor. On this task we compare our method to the Chan-Vese kernel $k$-means approach proposed in [55, 56] as a baseline method. Since the features are almost linearly separable, we use a simple linear kernel for both our model and the baseline model. As it is shown in Fig. 4, the $k$-means approach often fails to find a good cluster-center assignment, despite of strong supervision (provided in the form of user-scribbles) and richness of the features. This is due to the fact that deep high dimensional features are in general not $k$-means friendly [61], which means further pre-processing or a $k$-means suited kernel would be required. In contrast, our method provides a reasonable result for all cases, without the need for feature-pre-processing or kernel-parameter tuning.

5 Conclusion

We considered the joint solution of MAP-inference in MRFs and parameter learning, which can be viewed as a transductive inference problem. To solve this task, we proposed a novel algorithm that jointly optimizes over the discrete label variables and the continuous model parameters. The proposed method is related to classical ADMM from continuous optimization and admits a convergence proof under suitable assumptions even though the objective function is discrete-continuous and nonconvex. Our algorithm makes use of a decoupling of the problem into purely discrete and purely continuous subproblems and can be implemented in a distributed fashion. We evaluated our approach in several experiments including video object segmentation and interactive image segmentation. Our results suggest that the proposed optimization method performs favorable compared to alternating optimization (as in $k$-means) and convex relaxations. In particular, this indicates that the presented method also serves as an alternative approach to optimization problems arising in semi-supervised or transductive learning, e.g., in the case of SVMs. Furthermore, the visual results show that the transductive inference model is able to reduce the hallucination of false object pixels in image and video segmentation tasks.

Appendix A Theoretical Results

In the remainder of this section we make use of the following properties of $L$-smooth functions (known as the descent-lemma) and $m$-semiconvexity [2, 41], which are standard results and therefore stated without proof:

Lemma 3**.**

Let $f:\mathbb{R}^{k\times n}\to\mathbb{R}$ be continuously differentiable and let $x,y\in\mathbb{R}^{k\times n}$.

• •

If $f$ is $L$-smooth (meaning that $\nabla f$ is Lipschitz continuous with modulus $L$), then

[TABLE]

• •

If $f$ is $m$-semiconvex (meaning that $f+\frac{m}{2}\|\cdot\|_{F}^{2}$ is convex), then

[TABLE]

For showing convergence we make the following assumptions on our problem:

• •

The function $f$ is $L$-smooth, $m$-semiconvex and lower-bounded.

• •

For all $y_{i}\in\mathcal{L}$, $\ell(y_{i};\cdot)$ is lower-bounded.

• •

The kernel matrix $K\in\mathbb{R}^{|\mathcal{V}|\times|\mathcal{V}|}$ is surjective, i.e. the smallest eigenvalue $\sigma_{\min}(K^{\top}K)>0$ is positive.

• •

After finitely many iterations the penalty parameter $\rho$ is sufficiently large and kept fixed such that condition (18) holds.

A.1 Proof of Lemma 1

In [37, 20], to show convergence of nonconvex ADMM, a monotonic decrease of the augmented Lagrangian is guaranteed. Following a similar line of argument, we show that the “discrete-continuous” augmented Lagrangian (8) monotonically decreases with the iterates. Whereas its value decreases with the primal and discrete variable updates, the dual update yields a positive contribution to the overall estimate. Yet, for $\rho>0$ chosen large enough, $K$ surjective and $f$ being $L$-smooth, this ascent can be dominated by a sufficiently large descent in the primal block $\alpha$, updated last.

We need the following notation. Let $B_{:,y^{t}}^{t+1}$ denote the matrix whose $i$-th row is given by $B^{t+1}_{i,y^{t}_{i}}$. In particular, by definition of the $\beta$ update, this means $\beta^{t+1}=B^{t+1}_{:,y^{t+1}}$.

Proof.

We rewrite the difference of two consecutive “discrete-continuous” augmented Lagrangians as

[TABLE]

We now bound each of the four differences separately:

Since the augmented Lagrangian is separable in $\beta$ and we solve for any $y_{i}$ a minimization problem in $\beta_{y_{i}}$ globally optimal we have that

[TABLE]

A similar estimate holds for the the discrete variable $y^{t+1}$ due to the update in the algorithm:

[TABLE]

Now we devise a bound for the third term given by

[TABLE]

We apply the identity $\|a+c\|_{F}^{2}-\|b+c\|_{F}^{2}=-\|b-a\|_{F}^{2}+2\langle a+c,a-b\rangle$ with $a:=K\alpha^{t+1}$, $b:=K\alpha^{t}$ and $c=-\beta^{t+1}$ and obtain

[TABLE]

The optimality condition for the update of the variable $\alpha$ is given as

[TABLE]

We replace the term $\langle K\alpha^{t+1}-K\alpha^{t},\lambda^{t}+\rho(K\alpha^{t+1}-\beta^{t+1})\rangle=\langle\alpha^{t+1}-\alpha^{t},K^{\top}(\lambda^{t}+\rho(K\alpha^{t+1}-\beta^{t+1}))\rangle$ and obtain from the optimality condition of the $\alpha$ update that

[TABLE]

Moreover, due to the $m$-semiconvexity of the $f$ we know that

[TABLE]

Overall we can bound

[TABLE]

Since by assumption $K$ is surjective, the smallest eigenvalue of $K^{\top}K$ is greater than zero: $\sigma_{\min}(K^{\top}K)>0$. This means there exists some $\rho>0$ large enough so that $\frac{m-\rho\sigma_{\min}(K^{\top}K)}{2}<0$.

Finally, we estimate the last term:

[TABLE]

From the update of the dual variable and the optimality condition for the $\alpha$ update (32) it follows that

[TABLE]

Further, since $f$ is $L$-smooth we know that

[TABLE]

Overall, we obtain

[TABLE]

This gives the bound for the last term:

[TABLE]

Then, by merging the four estimates we obtain the desired result.

We proceed showing the lower boundedness of $\{\mathfrak{L}_{\rho}(\alpha^{t+1},\beta^{t+1},\lambda^{t+1},y^{t+1})\}_{t\in\mathbb{N}}$. Since $K$ is surjective, there exists $\alpha^{\prime}$ such that $K\alpha^{\prime}=\beta^{t+1}$ and it holds that

[TABLE]

Let $\rho>\frac{L}{\sigma_{\min}(K^{\top}K)}$. Then, since $f$ is $L$-smooth, we have

[TABLE]

Overall, since by assumption $f$ and $\ell(y_{i};\cdot)$ are bounded from below (for all $y_{i}\in\mathcal{L}$), this means

[TABLE]

is bounded from below.

Since $\{\mathfrak{L}_{\rho}(\alpha^{t+1},\beta^{t+1},\lambda^{t+1},y^{t+1})\}_{t\in\mathbb{N}}$ is monotonically decreasing and bounded from below, $\{\mathfrak{L}_{\rho}(\alpha^{t+1},\beta^{t+1},\lambda^{t+1},y^{t+1})\}_{t\in\mathbb{N}}$ converges. This completes the proof. ∎

A.2 Proof of Lemma 2

Proof.

We sum over the estimate (19) which yields

[TABLE]

Due to the lowerboundedness, the infinite sums have to converge. This yields that $\|\alpha^{t+1}-\alpha^{t}\|_{F}\to 0$. Since $0\leq\sigma_{\min}(K^{\top}K)\|\lambda^{t+1}-\lambda^{t}\|_{F}^{2}\leq\frac{L^{2}}{\rho\sigma_{\min}(K^{\top}K)}\|\alpha^{t+1}-\alpha^{t}\|_{F}^{2}$ and $\sigma_{\min}(K^{\top}K)>0$ also $\|\lambda^{t+1}-\lambda^{t}\|_{F}\to 0$. Since due to the dual update $\lambda^{t+1}-\lambda^{t}=\rho(K\alpha^{t+1}-\beta^{t+1})$, also $\|K\alpha^{t+1}-\beta^{t+1}\|_{F}\to 0$. Moreover, it holds that

[TABLE]

for $t\to\infty$.

Finally, suppose that there exists an infinite subsequence $\{t_{j}\}_{j=1}^{\infty}\subset\{t\}_{t=1}^{\infty}$ so that $y^{t_{j}+1}\neq y^{t_{j}}$. The last sum rewrites as,

[TABLE]

which diverges for $\delta>0$ positive. This however contradicts the lower boundedness of $\mathfrak{L}_{\rho}(\alpha^{t},\beta^{t},\lambda^{t},y^{t})$. ∎

A.3 Proof of Proposition 1

Proof.

Let $(\alpha^{*},\beta^{*},\lambda^{*},y^{*})$ be a limit point of $\{(\alpha^{t},\beta^{t},\lambda^{t},y^{t})\}_{t\in\mathbb{N}}$, and let $\{t_{j}\}_{j=1}^{\infty}\subset\{t\}_{t=1}^{\infty}$ be the corresponding subsequence of indices. The optimality conditions for the update of the variables $\beta_{i}$ (for any $i$) and $\alpha$ are given as:

[TABLE]

Passing the limit $j\to\infty$ and applying Lemma 2 we arrive at condittions (24)–(26). This completes the proof. ∎

A.4 Proof of Proposition 2

Proof.

Let $\delta>0$. Then, due to Lemma 2 the discrete variable converges, i.e. there is $T>0$ so that for all $t>T$

[TABLE]

Then, since $f$ and $\ell(y_{i};\cdot)$ are convex proper and lsc., after finitely many iterations our scheme Alg. 1 reduces to convex ADMM and the global convergence is a direct consequence of [17, 14, 6]. This completes the proof. ∎

Appendix B Additional Experimental Results

As a proof of concept we conduct a synthetic experiment with data sampled from 2D moon-shape distributions (600 samples, 4 classes, 150 per class). We sample 25 (possibly overlapping) cliques $C\subset\mathcal{V}$ of cardinality 25 from the set of examples. The synthetic labeling prior in this experiment is given in terms of constraints, that balance the label assignment within each clique. More precisely, it restricts the maximal deviation of the determined labeling from the true labeling to a given bound within each clique $C\in\mathcal{C}$. Mathematically, the higher order energies $E_{C}$ in the MRF are defined so that $E_{C}(y_{C})=0$ if $L_{C}^{j}\leq|\{i\in C:y_{C}^{i}=j\}|\leq U_{C}^{j}$, and $\infty$ otherwise. The bounds $L_{C}^{j}$ and $U_{C}^{j}$ are fixed and chosen a-priori, such that the number of samples $i\in C$ assigned to class $j$ deviates by at most 3 from the true number within clique $C$. This means that we do not provide any exact labels to the algorithm.

The overall task is to infer the correct labels from both, the distribution of the examples in the feature space, and the combinatorial prior encoded within the higher order energies. Within the algorithm, we solve the LP-relaxation of the higher order MRF-subproblem (14) with the dual-simplex method and threshold the solution. On this task, we compare our method to constrained kernel $k$-means and plain discrete-continuous coordinate descent on (5) with an RBF kernel and an SVM-loss (see Figure 5). Like [55, 56, 59, 3], we apply $k$-means in the RBF kernel space and solve the E-step w.r.t. to (14). It can be seen that both, discrete-continuous coordinate descent on the SVM-based model (Figure 5(c)) and constrained kernel $k$-means (Figure 5(b)) get stuck in poor local minima. In contrast, our method is able to infer the correct labels of most examples and finds a reasonable classifier, even for a trivial initialization of the parameters, cf. Figure 5(d). The label errors are 66.6% for constrained RBF kernel $k$-means, 68.5% for coordinate descent and 2.5% for our method.