# Transductive Few-Shot Learning:
Clustering is All You Need?

Imtiaz Masud Ziko, Malik Boudiaf, Jose Dolz, Eric Granger and Ismail Ben Ayed
Imtiaz Masud Ziko is with Thales - Cortaix, Canada and ÉTS Montreal, Canada.
E-mail<EMAIL_ADDRESS>Malik Boudiaf, Jose Dolz, Eric Granger
and Ismail Ben Ayed are with ÉTS Montreal, Canada.

###### Abstract

We investigate a general formulation for clustering and transductive few-shot
learning, which integrates prototype-based objectives, Laplacian
regularization and supervision constraints from a few labeled data points. We
propose a concave-convex relaxation of the problem, and derive a
computationally efficient block-coordinate bound optimizer, with convergence
guarantee. At each iteration, our optimizer computes independent (parallel)
updates for each point-to-cluster assignment. Therefore, it could be trivially
distributed for large-scale clustering and few-shot tasks. Furthermore, we
provides a thorough convergence analysis based on point-to-set maps. We report
comprehensive clustering and few-shot learning experiments over various data
sets, showing that our method yields competitive performances, in term of
accuracy and optimization quality, while scaling up to large problems. Using
standard training on the base classes, without resorting to complex meta-
learning and episodic-training strategies, our approach outperforms state-of-
the-art few-shot methods by significant margins, across various models,
settings and data sets. Surprisingly, we found that even standard clustering
procedures (e.g., K-means), which correspond to particular, non-regularized
cases of our general model, already achieve competitive performances in
comparison to the state-of-the-art in few-shot learning. These surprising
results point to the limitations of the current few-shot benchmarks, and
question the viability of a large body of convoluted few-shot learning
techniques in the recent literature. Our code is publicly available at
https://github.com/imtiazziko/SLK-few-shot.

###### Index Terms:

Clustering, few-shot learning, transductive inference, bound optimization,
regularization, point-to-set maps.

##  1 Introduction

Deep learning models achieve outstanding performances when trained on large-
scale labeled data sets. However, they are seriously challenged when dealing
with new classes that were not observed during training, even if a handful of
labeled samples becomes available for these new classes. To address this
challenge, few-shot learning methods have recently attracted wide interest
within the machine learning and computer vision communities. In standard few-
shot learning settings, we assume that we have access to training data
containing labeled samples from an initial set of classes, commonly referred
to as the base classes. A fully supervised model is first trained on this base
data. Then, supervision for a new set of classes that are completely different
from those seen in the base training is restricted to just one or a few
labeled samples per class. Model generalization is assessed over few-shot
tasks: Each task involves a query set, which contains unlabeled samples for
evaluation, and is supervised by a support set, which contains a few labeled
samples per new class.

The recent and abundant few-shot literature is widely dominated by convoluted
meta-learning and episodic-training strategies, which create a set of few-shot
tasks (or episodes) during base training, each containing support and query
samples, so as to mimic the generalization challenges encountered at testing
times. Very popular methods within this framework include [57, 67, 15, 50].
Prototypical networks [57] use an embedding prototype for each class, and
optimize the log-posteriors of the query points, within each episode during
base training. Matching networks [67] model the predictions of query samples
as linear functions of the support labels, and use memory architectures along
with episodic-training strategies. MAML (Model-Agnostic Meta-Learning) [15]
attempts to enable a model to become easy to fine-tune, and the meta-learner
in [50] views optimization as a model for few-shot learning. These widely used
methods were followed by a large body of few-shot learning literature based on
meta-learning, e.g., [61, 47, 40, 53, 77, 25, 79], to list a few.

Transductive few-shot learning: Recently, a large body of works investigated
transductive inference for few-shot tasks, e.g., [77, 49, 25, 13, 26, 73, 23,
78, 90, 38, 39, 5]. The transductive setting leverages the statistics of the
unlabeled query set and, in general, yields substantial improvements in
performances over inductive inference. Such boosts in performances are in line
with established facts in classical transductive inference for semi-supervised
learning [66, 29, 87], well-known to outperform its inductive counterpart on
small labeled data sets. While transductive inference uses exactly the same
training and testing data as its inductive counterpart111In the transductive
setting, each few-shot task is treated independently of the other tasks.
Therefore, the setting does not assume to have access to additional unlabeled
data. This is different from semi-supervised few-shot learning [51], which
uses additional unlabeled data., it classifies all the unlabeled query points
of a few-shot task jointly, rather than one data point at a time as in
inductive inference.

As pointed out in [6], most of meta-learning methods rely critically on
transductive batch normalization (TBN) to achieve competitive performances,
for instance, the works in [15, 20, 86], among others. Adopted initially in
the widely used MAML [15], TBN performs normalization using the statistics of
the query set of a given few-shot task, and yields significant increases in
performances [6]. Therefore, due to the popularity of MAML, several meta-
learning techniques have used TBN. There are also several recent methods that
proposed specialized architectures and loss functions for transductive few-
shot inference. For instance, the meta-learning, graph-based method in [77]
propagates labels from labeled to unlabeled data points, following label-
propagation concepts, which are popular in the classical semi-supervised
learning literature [87]. The transductive inference in [13] encourages
confident class predictions at the unlabeled query points of a few-shot task
by minimizing the Shannon entropy, in conjunction with the standard cross-
entropy loss over the labeled support samples. It is worth noting that the
general principle of minimizing entropy to leverage unlabeled data is widely
used in semi-supervised learning [21, 41, 4]. The method in [5] maximized the
mutual information between the query samples of a task and their latent
labels, and could be viewed as an extension of the entropy loss in [13]. The
method in [90] regularized the inductive predictions of query samples with a
pairwise Laplacian term.

It is important to note that the computational times of transductive-inference
techniques are, in general, much heavier than inductive inference. For
example, the entropy minimization in [13] updates all the trainable parameters
of a deep network at inference, requiring several orders of magnitude more
computations than inductive inference. Another example is the transductive
label-propagation inference is [77], whose complexity is cubic in the number
of unlabeled query points, as it requires matrix inversion.

###  1.1 Contributions

We view transductive few-shot inference as a constrained clustering, and
investigate a general multi-term formulation of the problem, which integrates:
a prototype-based objective (e.g, K-means or K-modes); a pairwise Laplacian
regularizer; and supervision constraints from the support set. We propose a
concave-convex relaxation of the objective, for which we derive a
computationally efficient block-coordinate bound optimizer, with convergence
guarantee. At each iteration, our optimizer computes independent (parallel)
updates for each point-to-cluster assignment. Therefore, it could be trivially
distributed for large-scale clustering and few-shot tasks. Furthermore, we
provides a thorough convergence analysis based on point-to-set maps. We report
comprehensive clustering and few-shot experiments over various data sets,
showing that our method yields competitive performances, both in term of
accuracy and optimization quality (i.e., the objectives obtained at
convergence), while scaling up to large problems. Without resorting to complex
meta-learning and episodic-training schemes, and using a standard training on
the base classes, our approach outperforms state-of-the-art few-shot methods
by significant margins across different models, settings, and data sets.
Furthermore, and surprisingly, we found that even standard clustering
procedures (e.g., K-means), which correspond to particular, non-regularized
cases of our general model, already achieve competitive performances in
comparison to the the state-of-the-art in few-shot learning (even though these
are not the best-performing variants of our general model). These surprising
results point to the limitations of current few-shot benchmarks, and question
the viability of a large body of convoluted few-shot learning techniques in
the recent literature.

###  1.2 Related works

Our work integrates several powerful ideas in clustering: pairwise graph-
Laplacian regularization, cluster-representative objectives and density-mode
finding. Cluster-representative objectives, such as the standard K-means, are
very popular because they provide prototype representations of the clusters,
explicit cluster-assignment variables and straightforward out-of-sample
extensions. Furthermore, typically, they are amenable to optimizers that scale
up to large and high-dimensional data sets. Pairwise Laplacian regularization
encourages data points that are neighbors in the feature space to have
consistent latent assignments, and is widely used in the contexts of spectral
graph clustering [56, 69, 55, 62] and semi-supervised learning [3, 74, 27, 63,
10, 87]. Unlike standard prototype-based methods, Laplacian regularization can
separate non-convex (or manifold-structured) clusters. However, despite their
complementary benefits, graph-Laplacian and prototype-based terms are, in
general, used separately in the clustering literature, except a very few
exceptions [89, 71, 62]. This might be due to the significant differences in
the optimization algorithms used to tackle these two different types of
objectives.

On the one hand, prototype-based objectives are commonly tackled with block-
coordinate descent procedures, which alternate optimization with respect to
the cluster prototypes and discrete assignment variables. It is worth noting
that optimizing prototype-based objectives over discrete variables is NP-
hard222For instance, a proof of the NP-hardness of K-means could be found in
[1].. On the other hand, optimizing the Laplacian term over discrete
variables, which is also NP-hard [64], is commonly tackled by relaxing the
integer constraints. For instance, in the context of graph clustering,
spectral relaxation [69, 56] is widely used for optimizing the Laplacian term
subject to balancing constraints. The problem is expressed in the form of a
generalized Rayleigh quotient, which yields a closed-form solution in terms of
the $K$ largest eigenvectors of the affinity matrix ($K$ is the number of
clusters). However, spectral relaxation is not applicable to joint objective
integrating K-prototypes and Laplacian regularization, as in our formulation.
Furthermore, spectral relaxation has high computational and memory load for
large data sets: For $N$ data points, one has to store an $N\times N$ affinity
matrix and compute its eigenvalue decomposition. The complexity is cubic with
respect to $N$ for a straightforward implementation and, to our knowledge,
super-quadratic for fast implementations [64]. Another approach is to replaces
the integer constraint with a probability-simplex constraint, which results in
a convex relaxation of the Laplacian term [71]. Such a direct convex
relaxation could be integrated with prototype-based objectives [71]. However,
unfortunately, it requires solving for $N\times K$ variables jointly.
Furthermore, it requires additional projections onto the $K$-dimensional
simplex, with a quadratic complexity with respect to $K$. Therefore, as we
will see in our experiments, the relaxation in [71] does not scale up for
large-scale problems (i.e., when $N$ and $K$ are large).

Laplacian regularization is also widely used in semi-supervised learning [3,
74, 27, 63, 88, 36]. This includes, most notably, label-propagation methods
[88, 10, 87], which spread soft labels from a few labeled samples to unlabeled
ones. Popular label-propagation algorithms [88, 87] could be unified under a
general convex-relaxation framework [10], which minimizes a quadratic
objective function, integrating a Laplacian regularization term and a
supervised-learning term defined over the labeled data points, along with
additional terms to avoid degenerate solutions. While the relaxation is
convex, solving the ensuing necessary condition for an optimum has high
computational and memory load as it requires an $N\times N$ matrix inversion.
In general, such matrix inversions in the standard label-propagation
algorithms take O($N^{3}$) time and O($N^{2}$) memory [36]. This impedes the
use of label propagation for large-scale data sets.

Our concave-convex relaxation computes independent (parallel) pointwise
updates during each iteration, and has important computational and memory
advantages over the above-discussed methods. It does not require eigenvalue
decomposition/matrix inversion (unlike spectral relaxation and label
propagation), neither does it perform expensive projection steps and
Lagrangian-dual inner iterates for the simplex constraints of each point
(unlike convex relaxation). In fact, investigating the scalability of spectral
relaxation for large-scale problems is an active research subject, for both
laplacian regularization and clustering in general [55, 36, 64, 68]. For
instance, the studies in [55, 64] investigated deep-learning approaches to
spectral clustering, so as to ease the scalability issues for large data sets,
and the authors of [68] examined the variational Nyström method for large-
scale spectral problems, among many other efforts on the subject. Also, there
is an active line of research works to mitigate the scalability issue of label
propagation [36], e.g., via graph approximations or stochastic optimization
methods. In fact, such scalability issues are being actively investigated even
for basic clustering procedures such as K-means [18, 45].

Finally, we point to the preliminary clustering results of this work, which we
published at the NeurIPS conference [89]. This journal extension expands
significantly on [89], and generalizes the work in many different ways,
including: (i) a generalized and constrained version of our formulation, along
with its adaptation to transductive few-shot learning; (ii) a new convergence
analysis based on point-to-set maps and Cauchy sequences; (iii) a completely
new set of comprehensive experiments and comparisons in the context of few-
shot classification, along with surprising findings; and (iv) new discussions
of the recent few-shot learning literature.

##  2 Proposed Formulation

We advocate Laplacian K-prototypes for clustering and transductive few-shot
inference, and propose a general bound- and relaxation-based solution of the
problem, which can be used to regularize several prototype-based objective
(e.g, K-means, K-median or K-modes). Importantly, our concave-convex
relaxation yields a parallel-structure algorithm, with convergence guarantee,
and can handle large-scale data sets. First, we discuss a general formulation
for clustering. Then, we describe an adaptation to transductive few-shot
learning, which has recently attracted significant interest. In this case, our
model takes the form of a constrained clustering for a given few-shot task.

Let $\mathcal{X}=\\{{\mathbf{x}}_{p}\in\mathbb{R}^{d},p=1,\dots,N\\}$ denotes
a set of feature inputs, which need to be assigned to $K$ different clusters
$\mathcal{S}_{k},\,k=1,\dots,K$. For each point ${\mathbf{x}}_{p}$, we define
a latent one-hot assignment vector
${\mathbf{s}}_{p}=[s_{p,1},\dots,s_{p,K}]^{t}\in\\{0,1\\}^{K}$, which is
constrained to be within the $K$-dimensional simplex: $s_{p,k}=1$ if
${\mathbf{x}}_{p}$ belongs to cluster $\mathcal{S}_{k}$ and $s_{p,k}=0$
otherwise. Superscript $t$ denotes the transpose operator. Let ${\mathbf{S}}$
denote the $N\times K$ matrix whose rows are formed by assignment variables
${\mathbf{s}}_{p}$, $p=1,\dots,N$. Finally, each cluster $\mathcal{S}_{k}$ is
represented by a prototype variable ${\mathbf{m}}_{k}\in\mathbb{R}^{d}$ in the
feature space. Let ${\mathbf{M}}$ be the $K\times d$ prototype matrix. We
consider the following general Laplacian K-prototypes model for clustering:

$\displaystyle\min_{{\mathbf{S}},{\mathbf{M}}}$ $\displaystyle\quad{{\cal
E}}({\mathbf{S}},{\mathbf{M}})={{\cal
F}({\mathbf{S}},{\mathbf{M}})}+\frac{\lambda}{2}{{\cal L}({\mathbf{S}})}$ s.t.
$\displaystyle\quad{\mathbf{1}}^{t}{\mathbf{s}}_{p}=1;\,{\mathbf{s}}_{p}\in\\{0,1\\}^{K}\,\forall
p$ (1)

General model (2) integrates, within a single objective, the advantages of
several powerful and well-known ideas in clustering and SSL, which makes it
very convenient:

  * •

Prototype-based term: ${\cal F}({\mathbf{S}},{\mathbf{M}})$ is a standard
prototype-based objective such as K-means or K-modes; see Table I for detailed
expressions. It identifies cluster prototypes (or representatives)
${\mathbf{m}}_{k}$, for instance, modes or means. It is worth noting that the
K-modes term in Table I is closely related to density mode estimation methods,
for instance, the very popular mean-shift algorithm [12]. The value of
${\mathbf{m}}_{k}$ globally optimizing this term for a given fixed cluster
$\mathcal{S}_{k}$ is the mode of the kernel density of feature points within
the cluster [62]. Therefore, the K-modes term, as in [9], can be viewed as an
energy-based formulation of mean-shift algorithms with a fixed number of
clusters [62]. Cluster modes are valid data points in the input set and,
therefore, convenient representations for manifold-structured, high-
dimensional inputs such as images, where simple parametric prototypes such as
the means, as in K-means, may not be good representatives.

  * •

Laplacian term: ${\cal L}({\mathbf{S}})$ is the well-known graph Laplacian
regularizer:

$\displaystyle{{\cal L}({\mathbf{S}})}$ $\displaystyle=$
$\displaystyle\sum_{p,q}w_{\mathcal{L}}({\mathbf{x}}_{p},{\mathbf{x}}_{q})\|{\mathbf{s}}_{p}-{\mathbf{s}}_{q}\|^{2}$
(2)

This term can be equivalently written as
$\mathrm{tr}({\mathbf{S}}^{t}{\mathbf{L}}{\mathbf{S}})$, with ${\mathbf{L}}$
the Laplacian matrix corresponding to affinity matrix
${\mathbf{W}}_{\mathcal{L}}=[w_{\mathcal{L}}({\mathbf{x}}_{p},{\mathbf{x}}_{q})]$.
Laplacian regularization encourages nearby data points to have similar latent
assignments, and is widely used in spectral clustering [56, 69, 55, 62], as
well as in SSL [3, 74, 27, 63].

Model (2) is flexible: On the one hand, Laplacian term ${{\cal
L}({\mathbf{S}})}$ enables to handle non-convex (or manifold-structured)
clusters, unlike standard prototype-based clustering techniques such as
K-means. On the other hand, and unlike spectral clustering methods, term
${{\cal F}({\mathbf{S}},{\mathbf{M}})}$ yields prototype representations of
the clusters. Furthermore, as we will see later in the few-shot setting and in
our experiments, model (2) enables to easily and effectively embed supervision
information from a few labeled examples in the prototype terms. It is not
clear how to embed effectively such supervision constraints in spectral
clustering, as the latter takes the form of a generalized Rayleigh quotient
optimization.

###  2.1 Adaptation to transductive few-shot inference

In standard few-shot learning settings, we assume that we have access to a
training set $\mathcal{X}_{\text{base}}$ containing labeled data points from
an initial set of classes, referred to as the base classes. A fully supervised
deep-network model is first trained on this base data, yielding embedding
function ${\mathbf{x}}_{p}=f_{\theta}({\mathbf{i}}_{p})$ for a given data
input ${\mathbf{i}}_{p}$, with $\theta$ denoting the base-training parameters.
Then, supervision for a new set of classes, which are completely different
from those seen in the base training, is restricted to just one or a few
labeled samples per class. Model generalization is assessed over few-shot
tasks: Each task involves a query set, which contains unlabeled samples for
evaluation, and is supervised by a support set, which contains a few labeled
samples per new class. For a given few-shot task, let
$\mathcal{X}_{\text{s}}=\bigcup_{k=1}^{K}\mathcal{X}_{\text{s}}^{k}$ denote
the labeled support set, with $K$ the number of classes of the task. Each new
class $k$ has $|\mathcal{X}_{\text{s}}^{k}|$ labeled examples, for instance,
$|\mathcal{X}_{\text{s}}^{k}|=1$ for 1-shot and
$|\mathcal{X}_{\text{s}}^{k}|=5$ for 5-shot. The objective is, therefore, to
accurately classify the unlabeled samples of query set
$\mathcal{X}_{\text{q}}=\bigcup_{k=1}^{K}\mathcal{X}_{\text{q}}^{k}$ into one
of these $K$ test classes. This setting is referred to as the
$|\mathcal{X}_{\text{s}}^{k}|$-shot $K$-way few-shot learning. For each few-
shot task with concatenated features
${\mathcal{X}=\\{{\mathbf{x}}_{p}\in\mathbb{R}^{D},p=1,\dots,N\\}=\mathcal{X}_{\text{s}}\cup\mathcal{X}_{\text{q}}}$,
we pose the problem as optimizing a Laplacian K-prototypes objective, subject
to supervision constraints from the support set:

$\displaystyle\min_{{\mathbf{S}},{\mathbf{M}}}\ {{\cal
E}}({\mathbf{S}},{\mathbf{M}})$ $\displaystyle=$ $\displaystyle{{\cal
F}({\mathbf{S}},{\mathbf{M}})}+\frac{\lambda}{2}{{\cal
L}({\mathbf{S}})}\quad\mbox{s.t.}$ $\displaystyle s_{p,k}$ $\displaystyle=$
$\displaystyle 1~{}\mbox{if}~{}{\mathbf{x}}_{p}\in\mathcal{X}_{\text{s}}^{k}$
$\displaystyle{\mathbf{1}}^{t}{\mathbf{s}}_{p}$ $\displaystyle=$
$\displaystyle 1,\;{\mathbf{s}}_{p}\in\\{0,1\\}^{K}\;\forall p$ (3)

TABLE I: Auxiliary functions of several prototype-based clustering objectives. ${\mathbf{S}}_{k}=[s_{1,k},\dots,s_{N,k}]^{t}\in\\{0,1\\}^{N}$ is the binary indicator vector of cluster $\mathcal{S}_{k}$. When used for variables ${\mathbf{S}}_{k}$ and ${\mathbf{m}}_{k}$, superscript $i$ means the current solution at outer iteration indexed by $i$. ${\mathbf{X}}$ is the $N\times d$ matrix whose rows are formed by feature inputs ${\mathbf{x}}_{p}$. Clustering |  ${\cal F}({\mathbf{S}},{\mathbf{M}})$ |  ${\mathbf{a}}_{p}^{i}=[a_{p,1}^{i},\dots,a_{p,K}^{i}]^{t}$ |  Where
---|---|---|---
K-means |  $\sum_{p}\sum_{k}s_{p,k}({\mathbf{x}}_{p}-{\mathbf{m}}_{k})^{2}$ |  $({\mathbf{x}}_{p}-{\mathbf{m}}_{k}^{i})^{2}$ |  ${\mathbf{m}}_{k}^{i}=\frac{{\mathbf{X}}^{t}{\mathbf{S}}_{k}^{i}}{{\mathbf{1}}^{t}{\mathbf{S}}_{k}^{i}}$
K-modes |  $-\sum_{p}\sum_{k}s_{p,k}w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}}_{k})$ |  $-w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}}_{k}^{i})$ |  $w_{\mathcal{F}}$ is kernel-based affinity
${\mathbf{m}}_{k}^{i}=\operatorname*{arg\,max}_{{\mathbf{m}}}\sum_{p}s_{p,k}^{i}w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}})$

##  3 Concave-convex relaxation

It is easy to verify that, for binary simplex variables, the Laplacian term
${\cal L}({\mathbf{S}})$ can be written as follows:

$\displaystyle{\cal L}({\mathbf{S}})$
$\displaystyle=\mathrm{tr}({\mathbf{S}}^{t}{\mathbf{L}}{\mathbf{S}})$ (4)
$\displaystyle=2\left(\sum_{p}{\mathbf{s}}_{p}^{t}{\mathbf{s}}_{p}d_{p}-\sum_{p,q}w_{\mathcal{L}}({\mathbf{x}}_{p},{\mathbf{x}}_{q}){\mathbf{s}}_{p}^{t}{\mathbf{s}}_{q}\right)$
$\displaystyle=2\left(\sum_{p}d_{p}-\sum_{p,q}w_{\mathcal{L}}({\mathbf{x}}_{p},{\mathbf{x}}_{q}){\mathbf{s}}_{p}^{t}{\mathbf{s}}_{q}\right)$

where the last equality is valid only for binary (integer) variables, and
$d_{p}$ denotes the degree of data point $p$:

$d_{p}=\sum_{q}w_{\mathcal{L}}({\mathbf{x}}_{p},{\mathbf{x}}_{q})$

When we replace the integer constraints ${\mathbf{s}}_{p}\in\\{0,1\\}$ by
${\mathbf{s}}_{p}\in[0,1]$, our relaxation becomes different from direct
convex relaxations of the Laplacian [71], which optimizes
$\mathrm{tr}({\mathbf{S}}^{t}{\mathbf{L}}{\mathbf{S}})$ subject to
probabilistic simplex constraints. For relaxed variables,
$\mathrm{tr}({\mathbf{S}}^{t}{\mathbf{L}}{\mathbf{S}})$ is a convex function
because the Laplacian is always positive semi-definite (PSD). In contrast, our
relaxation of the Laplacian term is concave for PSD affinity matrix
${\mathbf{W}}_{\mathcal{L}}=[w_{\mathcal{L}}({\mathbf{x}}_{p},{\mathbf{x}}_{q})]$.
As we will see later, concavity yields a scalable (parallel) algorithm for
large $N$, which computes independent updates for assignment variables
${\mathbf{s}}_{p}$. Our updates can be trivially distributed, and do not
require storing a full $N\times N$ affinity matrix. These are important
computational and memory advantages over direct convex relaxations of the
Laplacian [71], which require solving for $N\times K$ variables all together
as well as expensive simplex projections, and over common spectral relaxations
[69], which require storing a full affinity matrix and computing its
eigenvalue decomposition.

We propose the following concave-convex relaxation of the objective in (2) by
relaxing the integer constraints and adding a convex negative-entropy barrier
for constraints ${\mathbf{s}}_{p}\geq 0$:

${{\cal R}}({\mathbf{S}},{\mathbf{M}})={\cal
E}({\mathbf{S}},{\mathbf{M}})+\mathrm{tr}({\mathbf{S}}\log{\mathbf{S}}^{t})$
(5)

${\cal R}({\mathbf{S}},{\mathbf{M}})$ is minimized over each
${\mathbf{s}}_{p}\in\nabla_{K}$, with $\nabla_{K}$ denoting the
$K$-dimensional probability simplex:

$\nabla_{K}=\\{{\mathbf{s}}\in[0,1]^{K}\;|\;{\mathbf{1}}^{t}{\mathbf{s}}=1\\}$

It is easy to check that, at the vertices of the simplex, our relaxation in
(5) is equivalent to the initial discrete objective in Eq. (2). The last term
we introduced in Eq. (5) is a convex negative-entropy barrier function, which
completely avoids expensive projection steps and Lagrangian-dual inner
iterations for the simplex constraints of each point. First, this entropy
barrier restricts the domain of each ${\mathbf{s}}_{p}$ to non-negative
values, which avoids extra dual variables for constraints
${\mathbf{s}}_{p}\geq 0$. Second, the presence of such a barrier function
yields closed-form updates for the dual variables of constraints
$\mathbf{1}^{t}{\mathbf{s}}_{p}=1$. In fact, entropy-like barriers are
commonly used in Bregman-proximal optimization [80], and have well-known
computational and memory advantages when dealing with the challenging simplex
constraints [80]. Surprisingly, to our knowledge, they are not common in the
context of clustering. In machine learning, such entropic barriers appear
frequently in the context of conditional random fields (CRFs) [30, 31], but
are not motivated from an optimization perspective; they result from standard
probabilistic and mean-field approximations of CRFs [30].

##  4 Optimization

In this section, we derive a block-coordinate bound optimization algorithm for
tackling relaxation (5). Our algorithm alternatively computes parallel updates
of assignment variables ${\mathbf{s}}_{p}$ (${\mathbf{s}}$-updates), for each
point $p$ and independently of the other points, along with prototype updates
(${\mathbf{M}}$-updates) at each outer iteration, with guaranteed convergence.
As we will see in our experiments, our bound optimizer yields consistently
lower values of objective ${\cal E}$ at convergence than the proximal
algorithm in [71], while being highly scalable to large-scale and high-
dimensional problems. Our SLK (Scalable Laplacian K-prototypes) algorithm
alternates the following two steps until convergence:

  * •

Optimizing $\cal{R}({\mathbf{S}},{\mathbf{M}})$ w.r.t ${\mathbf{S}}$, with
${\mathbf{M}}$ fixed (subsection 4.1).

  * •

Optimizing $\cal{R}({\mathbf{S}},{\mathbf{M}})$ w.r.t ${\mathbf{M}}$, with
${\mathbf{S}}$ fixed (subsection 4.2).

###  4.1 Bound optimization for the ${\mathbf{S}}$-updates

At outer iteration $i$, we keep ${\mathbf{M}}={\mathbf{M}}^{i}$ fixed and
minimize
$\mathcal{R}^{i}({\mathbf{S}})=\mathcal{R}({\mathbf{S}},{\mathbf{M}}^{i})$
with respect to ${\mathbf{S}}$. To tackle efficiently the minimization of our
relaxation $\mathcal{R}^{i}$, we minimize a series of auxiliary functions
${\cal A}^{i,n}$, indexed by inner counter $n\in\mathbb{N}$. As an auxiliary
function, ${\cal A}^{i,n}$ is an upper bound on ${\cal R}^{i}$, is tight at
the current inner solution ${{\mathbf{S}}}^{i,n}$ and is chosen to make
optimization easier. Formally:

######  Definition 1.

At outer iteration $i$, and inner iteration $n$, ${\cal
A}^{i,n}({\mathbf{S}})$ is an auxiliary function of
$\mathcal{R}^{i}({\mathbf{S}})$ at current solution ${\mathbf{S}}^{i,n}$ if it
satisfies:

$\displaystyle\mathcal{R}^{i}({\mathbf{S}})\,$ $\displaystyle\leq\,{\cal
A}^{i,n}({\mathbf{S}}),\,\forall~{}{\mathbf{S}}$ (6a)
$\displaystyle\mathcal{R}^{i}({\mathbf{S}}^{i,n})$ $\displaystyle={\cal
A}^{i,n}({\mathbf{S}}^{i,n})$ (6b)

In general, bound optimizers update current solution ${\mathbf{S}}^{i,n}$ to
the optimum of the auxiliary function:
${\mathbf{S}}^{i,n+1}=\arg\min_{{\mathbf{S}}}~{}{\cal A}^{i,n}({\mathbf{S}})$.
This guarantees that, at each inner iteration, the original objective function
does not increase:

$\displaystyle{\cal R}^{i}({\mathbf{S}}^{i,n+1})\leq{\cal
A}^{i,n}({\mathbf{S}}^{i,n+1})\leq{\cal A}^{i,n}({\mathbf{S}}^{i,n})={\cal
R}^{i}({\mathbf{S}}^{i,n})$ (7)

Bound optimizers can be very effective as they transform difficult problems
into easier ones [85]. Examples of well-known bound optimizers include the
concave-convex procedure (CCCP) [81], expectation maximization (EM) algorithms
and submodular-supermodular procedures (SSP) [44], among others. Furthermore,
bound optimizers are not restricted to differentiable functions333Our
objective is not differentiable with respect to the modes as each of these is
defined as the maximum of a function of the assignment variables., neither do
they depend on optimization parameters such as step sizes.

######  Proposition 1.

At outer iteration $i$, inner iteration $n$, and given the current solution
${\mathbf{S}}^{i,n}=[s_{p,k}^{i,n}]$, we have the following auxiliary function
(up to an additive constant) for concave-convex relaxation
$\mathcal{R}^{i}({\mathbf{S}})$ and psd444We can consider
${\mathbf{W}}_{\mathcal{L}}$ to be psd without loss of generality. When
${\mathbf{W}}_{\mathcal{L}}$ is not psd, we can use a diagonal shift for the
affinity matrix, i.e., we replace ${\mathbf{W}}_{\mathcal{L}}$ by
${\mathbf{W}}_{\mathcal{L}}+\delta{\mathbf{I}}_{N}$. Clearly,
${\mathbf{W}}_{\mathcal{L}}+\delta{\mathbf{I}}_{N}$ is psd for sufficiently
large $\delta$. For integer variables, this change does not alter the
structure of the minimum of discrete function ${\cal E}$. affinity matrix
${\mathbf{W}}_{\mathcal{L}}=[w_{\mathcal{L}}({\mathbf{x}}_{p},{\mathbf{x}}_{q})]_{p,q}$:

${\cal
A}^{i,n}({\mathbf{S}})=\sum_{p=1}^{N}{\mathbf{s}}_{p}^{t}\left(\log({\mathbf{s}}_{p})-{{\mathbf{a}}}_{p}^{i}-\lambda~{}{{\mathbf{b}}}_{p}^{i,n}\right)$
(8)

where ${{\mathbf{a}}}_{p}^{i}$ and ${{\mathbf{b}}}_{p}^{i,n}$ are the
following $K$-dimensional vectors:

$\displaystyle{{\mathbf{a}}}_{p}^{i}$
$\displaystyle=[a_{p,1}^{i},\dots,a_{p,K}^{i}]^{t},\,\mbox{\em
with}~{}a_{p,k}^{i}~{}\mbox{given in Table \ref{tab:a_p}}$ (9a)
$\displaystyle{{\mathbf{b}}}_{p}^{i,n}$
$\displaystyle=[b_{p,1}^{i,n},\dots,b_{p,K}^{i,n}]^{t},\,\mbox{\em
with}\,\,b_{p,k}^{i,n}=\sum_{q}w_{\mathcal{L}}({\mathbf{x}}_{p},{\mathbf{x}}_{q})s_{q,k}^{i,n}$
(9b)

###### Proof.

See Appendix A. ∎

Notice that the bound in Eq. (8) is the sum of independent functions, each
corresponding to a point $p$. As a result, both the bound and simplex
constraints ${\mathbf{s}}_{p}\in\nabla_{K}$ are separable over assignment
variables ${\mathbf{s}}_{p}$. Therefore, inner iteration $n$ corresponds to
minimizing the auxiliary function $\mathcal{A}^{i,n}$ by minimizing
independently each term in the sum over ${\mathbf{s}}_{p}$, subject to the
simplex constraint, while guaranteeing convergence:

$\min_{{\mathbf{s}}_{p}\in\nabla_{K}}{\mathbf{s}}_{p}^{t}(\log({\mathbf{s}}_{p})-{{\mathbf{a}}}_{p}^{i}-\lambda{{\mathbf{b}}}_{p}^{i,n}),\,\forall
p$ (10)

Note that, for each $p$, negative entropy
${\mathbf{s}}_{p}^{t}\log{\mathbf{s}}_{p}$ restricts ${\mathbf{s}}_{p}$ to be
non-negative, which removes the need for handling explicitly constraints
${\mathbf{s}}_{p}\geq 0$. This term is convex and, therefore, the problem in
(10) is convex: The objective is convex (sum of linear and convex functions)
and constraint ${\mathbf{s}}_{p}\in\nabla_{K}$ is affine. Hence, one can
minimize this constrained convex problem for each $p$ by solving the Karush-
Kuhn-Tucker (KKT) conditions555Note that strong duality holds since the
objectives are convex and the simplex constraints are affine. This means that
the solutions of the (KKT) conditions minimize the auxiliary function.. The
KKT conditions yield a closed-form solution for both primal variables
${\mathbf{s}}_{p}$ and the dual variables (Lagrange multipliers) corresponding
to simplex constraints ${\mathbf{1}}^{t}{\mathbf{s}}_{p}=1$. Each closed-form
update, which globally optimizes (10) and is within the simplex, is given by:

${\mathbf{s}}_{p}^{i,n+1}=\frac{\exp({{\mathbf{a}}}_{p}^{i}+\lambda{{\mathbf{b}}}_{p}^{i,n})}{{\mathbf{1}}^{t}\exp({{\mathbf{a}}}_{p}^{i}+\lambda{{\mathbf{b}}}_{p}^{i,n})}\,\,\forall\,p$
(11)

Inner updates (11) are repeated for each point $p$ until convergence. Letting
${\mathbf{s}}_{p}^{i,*}$ denote the stationary solution obtained for point $p$
at convergence, the outer update then reads:

$\displaystyle{\mathbf{s}}^{i}_{p}\leftarrow{\mathbf{s}}^{i,*}_{p}$ (12)

The complexity of each ${\mathbf{S}}$ inner iteration is $\mathcal{O}(N\rho
K)$, with $\rho$ the neighborhood size for the affinity matrix. Typically, we
use sparse matrices ($\rho<<N$). Note that the complexity becomes
$\mathcal{O}(N^{2}K)$ for dense matrices, i.e., the case where all the
affinities are non-zero. However, the update of each ${\mathbf{s}}_{p}$ can be
done for point $p$ independently of the other points, which enables parallel
implementations.

Input : ${\mathbf{X}}$, Initial prototypes ${\mathbf{M}}^{0}$

Output : ${\mathbf{S}}$ and ${\mathbf{M}}$

$i\leftarrow 1$

// outer iteration index

1 ${\mathbf{M}}^{i}\leftarrow{\mathbf{M}}^{0}$

2 repeat

// ${\mathbf{S}}$-updates

3 foreach  _$\mathbf{x}_{p}$_ do

4 Compute $\mathbf{a}^{i}_{p}$ from (9a)

5 Initialize
${\mathbf{s}}^{i,0}_{p}\leftarrow\frac{\exp\\{\mathbf{a}^{i}_{p}\\}}{\mathbf{1}^{t}\exp\\{\mathbf{a}^{i}_{p}\\}}$

6 Obtain ${\mathbf{s}}^{i,*}_{p}$ by repeating (11) until convergence

7 Update ${\mathbf{s}}^{i}_{p}$ using (18)

8

9 end foreach

// ${\mathbf{M}}$-updates

10 foreach  _$k$_ do

11 if  _SLK-MS_ then

12 Obtain ${\mathbf{m}}^{i,*}_{k}$ by repeating (17) until convergence

13 Update ${\mathbf{m}}^{i}_{k}$ using (18);

14 if  _SLK-Means_ then

15 Update ${\mathbf{m}}^{i}_{k}$ using (20);

16

17 end foreach

18
${\mathbf{M}}^{i+1}\leftarrow[{\mathbf{m}}^{i}_{1};\dots;{\mathbf{m}}^{i}_{K}]$

19 $i\leftarrow i+1$

20until  _convergence_

21${\mathbf{S}}\leftarrow[{\mathbf{s}}^{i}_{1};\dots;{\mathbf{s}}^{i}_{N}]$

22 ${\mathbf{M}}\leftarrow{\mathbf{M}}^{i}$

return ${\mathbf{S}}$, ${\mathbf{M}}$

Algorithm 1 SLK algorithm

###  4.2 Prototype updates

In this section, we give the update rules for prototype matrix ${\mathbf{M}}$.
Therefore, we consider ${\mathbf{S}}={\mathbf{S}}^{i}$ fixed to its value at
outer iteration $i$. These updates come in two different forms, one is closed-
form and the other iterative (based on convergent fixed-point iterations),
depending on the chosen prototype-based term $\cal
F({\mathbf{S}},{\mathbf{M}})$. In what follows, we detail the prototype
updates for the K-modes and K-means objectives.

####  4.2.1 SLK-MS

In this case, the goal is to find ${\mathbf{m}}^{i}_{k}$ that minimizes the
K-modes objective, with the assignment variables fixed, i.e.,

$\displaystyle{\mathbf{m}}^{i}_{k}=\operatorname*{arg\,min}_{{\mathbf{m}}}-\sum_{p}\sum_{k}s^{i}_{p,k}w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}})$
(13)

For the K-modes term, we use the well-known RBF kernel, defined as follows:

$\displaystyle
w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}})=\exp\left(-\frac{\lVert{\mathbf{x}}_{p}-{\mathbf{m}}\rVert^{2}}{2\sigma^{2}}\right)$
(14)

with $\sigma>0$ denoting a kernel width parameter. Setting the gradient of Eq.
(13) w.r.t ${\mathbf{m}}$ to ${\mathbf{0}}$ yields the following condition for
a stationary point of the K-modes objective in (13):

$\displaystyle{\mathbf{m}}-g^{i}_{k}({\mathbf{m}})={\mathbf{0}}$ (15) with
$\displaystyle
g^{i}_{k}({\mathbf{m}})=\frac{\sum_{p}s^{i}_{p,k}w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}}){\mathbf{x}}_{p}}{\sum_{p}s^{i}_{p,k}w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}})}$
(16)

The following proposition describes how to find the unique solution of Eq.
(15), via convergent fixed-point iterates:

######  Proposition 2.

For every outer iteration $i$, and every class $k$, the following fixed-point
updates:

$\displaystyle{\mathbf{m}}_{k}^{i,n+1}=g^{i}_{k}({\mathbf{m}}_{k}^{i,n})$ (17)

yield a Cauchy sequence $\\{{\mathbf{m}}_{k}^{i,n}\\}_{n\in\mathbb{N}}$. As
such, this sequence converges to a unique value
${\mathbf{m}}_{k}^{i,*}=\lim_{n\rightarrow\infty}{\mathbf{m}}_{k}^{i,n}$.
Furthermore, ${\mathbf{m}}_{k}^{i,*}$ is the unique solution of the
stationary-point condition in Eq. (15).

The proof of Proposition 2 is provided in Appendix B. K-modes prototypes’
updates therefore require inner loop iterations given by Proposition 2. Once
convergence is reached, the prototypes can be updated as:

$\displaystyle{\mathbf{m}}^{i}_{k}\leftarrow{\mathbf{m}}_{k}^{i,*}$ (18)

####  4.2.2 SLK-Means

In this case, for each class $k$, we compute ${\mathbf{m}}_{k}$ such that:

$\displaystyle{\mathbf{m}}^{i}_{k}=\operatorname*{arg\,min}_{{\mathbf{m}}}\sum_{p}\sum_{k}s^{i}_{p,k}\lVert{\mathbf{x}}_{p}-{\mathbf{m}}\rVert^{2}$
(19)

This problem is convex and can be solved in closed-form, yielding the class-
wise means as prototypes, as in the standard K-means procedure:

$\displaystyle{\mathbf{m}}^{i}_{k}\leftarrow\frac{{\mathbf{X}}^{t}{\mathbf{S}}^{i}_{k}}{{\mathbf{1}}^{t}{\mathbf{S}}^{i}_{k}}$
(20)

where ${\mathbf{S}}_{k}=[s^{i}_{1,k},\dots,s^{i}_{N,k}]^{t}\in\\{0,1\\}^{N}$
and ${\mathbf{X}}$ is the $N\times d$ matrix whose rows are formed by feature
inputs ${\mathbf{x}}_{p}$.

The pseudo-code for our full Scalable Laplacian K-modes (SLK) method is
provided in Algorithm 1.

##  5 Convergence study

Putting together the ${\mathbf{S}}$\- and ${\mathbf{M}}$-update steps, one
could notice that sequence
$\\{\mathcal{R}({\mathbf{S}}^{i},{\mathbf{M}}^{i})\\}_{i\in\mathbb{N}}$ is
non-increasing:

$\displaystyle\mathcal{R}({\mathbf{S}}^{i+1},{\mathbf{M}}^{i+1})\underbrace{\leq}_{{\mathbf{M}}\text{-update}}\mathcal{R}({\mathbf{S}}^{i+1},{\mathbf{M}}^{i})\underbrace{\leq}_{{\mathbf{S}}\text{-update}}\mathcal{R}({\mathbf{S}}^{i},{\mathbf{M}}^{i})$
(21)

Furthermore, this sequence is bounded from below:

$\displaystyle\mathcal{R}({\mathbf{S}},{\mathbf{M}})=\underbrace{\mathcal{F}({\mathbf{S}},{\mathbf{M}})}_{\leavevmode\resizebox{140.92792pt}{}{$\begin{cases}\geq
0&\text{SLK-Means}\\\ \geq-N&\text{SLK-
MS}\end{cases}$}}+\frac{\lambda}{2}\underbrace{\mathcal{L}({\mathbf{S}})}_{\geq
0}+\underbrace{\mathrm{tr}({\mathbf{S}}\log{\mathbf{S}}^{t})}_{-N\log K}$ (22)

This makes of
$\\{\mathcal{R}({\mathbf{S}}^{i},{\mathbf{M}}^{i})\\}_{i\in\mathbb{N}}$ a
converging sequence, regardless of initialization
$\\{{\mathbf{S}}^{0},{\mathbf{M}}^{0}\\}$. However, without additional
analysis, this does neither imply that sequence
$\\{{\mathbf{S}}^{i},{\mathbf{M}}^{i}\\}_{i\in\mathbb{N}}$ itself converges,
nor does it characterize its limit points.

In what follows, we analyze our proposed SLK algorithm through the lens of
Zangwill’s global convergence theory, which provides a simple but general
framework to study the convergence of iterative algorithms. Note that this
theory was already used to prove the convergence of the concave-convex [58]
and EM/GEM procedures [76], both of which could be viewed as special cases of
bound optimization. In particular, we show that all limits points of any
sequence $\\{{\mathbf{S}}^{i},{\mathbf{M}}^{i}\\}_{i\in\mathbb{N}}$ produced
by our algorithm are stationary points.

###  5.1 Background

In this section, we introduce the minimal set of required elements to
understand Zangwill’s theory, upon which our own convergence result is based.
We first introduce the central concept of a point-to-set map $\psi$, which
maps a point ${\boldsymbol{\theta}}\in\Theta$ to a set of points
$\psi({\boldsymbol{\theta}})\subset\Theta$. Map $\psi$ could be understood as
the iteration of an optimization algorithm that, from a point
${\boldsymbol{\theta}}^{i}$ in some space $\Theta$ of the optimization
variables, outputs a new point
${\boldsymbol{\theta}}^{i+1}\in\psi({\boldsymbol{\theta}}^{i})\subset\Theta$
from a set of local candidate solutions. We further recall the important
notion of closedness, which generalizes the concept of continuity in standard
point-to-point maps to point-to-set maps:

######  Definition 2.

(Closedness) Let us consider two converging sequences:

$\displaystyle{\boldsymbol{\theta}}^{n}\xrightarrow[n\to\infty]{}{\boldsymbol{\theta}}^{*}$
$\displaystyle\widetilde{{\boldsymbol{\theta}}}^{n}\xrightarrow[n\to\infty]{}\widetilde{{\boldsymbol{\theta}}}^{*}$

Let us further assume that:

$\displaystyle\forall
n\in\mathbb{N},~{}\widetilde{{\boldsymbol{\theta}}}^{n}\in\psi({\boldsymbol{\theta}}^{n})$
(23)

Then, the point-to-set map $\psi$ is said to be closed at point
${\boldsymbol{\theta}}^{*}$ if relation (23) extends to the limit when
$n\rightarrow\infty$, i.e:

$\displaystyle\tilde{{\boldsymbol{\theta}}}^{*}\in\psi({\boldsymbol{\theta}}^{*})$

The point-to-set map $\psi$ is said to be closed on the set $\Theta$ if it is
closed at every point of $\Theta$.

We are now ready to enunciate Zangwill’s theorem:

######  Theorem 1.

(p.29 in [83]) Consider a compact set $\Theta$, and let
$\psi:\Theta\rightarrow\mathcal{P}(\Theta)$ a point-to-set map, i.e., an
algorithm generating a sequence
$\\{{\boldsymbol{\theta}}^{i}\\}_{i\in\mathbb{N}}$ from an initial point
$\\{{\boldsymbol{\theta}}^{0}\\}$, via iteration
${\boldsymbol{\theta}}^{i+1}\in\psi({\boldsymbol{\theta}}^{i})$. Let
$\Gamma\subset\Theta$ denotes a given solution set, e.g., the set of fixed
points of the algorithm:
$\Gamma=\\{{\boldsymbol{\theta}}\in\Theta\,|\,\psi({\boldsymbol{\theta}})={\boldsymbol{\theta}}\\}$.
Assume that:
(1) $\psi({\boldsymbol{\theta}})$ is nonempty and $\psi$ is closed at
${\boldsymbol{\theta}}$, $\forall{\boldsymbol{\theta}}\in\Theta$
(2) There is a continuous function
$\mathcal{L}:{\boldsymbol{\theta}}\rightarrow\mathbb{R}$ such that:
$\forall{\boldsymbol{\theta}}\in\Theta$,
$\quad\forall~{}{\boldsymbol{\theta}}^{\prime}\in\psi({\boldsymbol{\theta}})$:

  * •

${\boldsymbol{\theta}}^{\prime}\notin\Gamma\implies\mathcal{L}({\boldsymbol{\theta}}^{\prime})<\mathcal{L}({\boldsymbol{\theta}})$

  * •

${\boldsymbol{\theta}}^{\prime}\in\Gamma\implies\mathcal{L}({\boldsymbol{\theta}}^{\prime})\leq\mathcal{L}({\boldsymbol{\theta}})$

Then any sequence $\\{{\boldsymbol{\theta}}^{i}\\}_{i\in\mathbb{N}}$ defined
by ${\boldsymbol{\theta}}^{i+1}\in\psi({\boldsymbol{\theta}}^{i})$ has all of
its limit points in $\Gamma$.

As noted in [58], the general idea to prove the convergence of an iterative
algorithm is to properly set $\Gamma$ and $\mathcal{L}$. The natural choice
for $\Gamma$ is to set it as the set of fixed points of the algorithm
$\Gamma=\\{{\boldsymbol{\theta}}\in\Theta\,|\,\psi({\boldsymbol{\theta}})={\boldsymbol{\theta}}\\}$.
$\mathcal{L}$ is the objective function that the algorithm minimizes. With
that in mind, assumption (2) in Theorem 1 simply ensures that, even though the
algorithm has not yet reached stationary points $\Gamma$, objective function
$\mathcal{L}$ is strictly decreasing. Theorem 1 is, unfortunately, not
straightforwardly applicable to our SLK algorithm, as it considers a single
algorithm update $\psi$, while SLK alternates between both ${\mathbf{S}}$\-
and ${\mathbf{M}}$-updates, which handle two completely different types of
optimization variables (assignments and prototypes). A corollary proposed in
[16] enables us to verify the assumptions on each of type of updates
independently while preserving the final result:

######  Corollary 1.

(Corollary 2. in [16]) Consider $\Theta$ a compact set, $\Gamma\subset\Theta$,
$J$ a finite set of indices, a set of point-to-set maps
$\\{\psi_{j}:\Theta\rightarrow\mathcal{P}(\Theta),~{}j\in J\\}$, and
$\mathcal{L}:{\boldsymbol{\theta}}\rightarrow\mathbb{R}$ a continuous
function. Using the sequence of indices $I:\mathbb{N}\rightarrow J$, assume
that for any $n\in\mathbb{N}$, for any
${\boldsymbol{\theta}}\in\Theta\setminus\Gamma$:
(1) $\psi_{I(n)}({\boldsymbol{\theta}})$ is nonempty and $\psi_{I(n)}$ is
closed at ${\boldsymbol{\theta}}$
(2)
${\boldsymbol{\theta}}^{\prime}\in\psi_{I(n)}({\boldsymbol{\theta}})\implies\mathcal{L}({\boldsymbol{\theta}}^{\prime})<\mathcal{L}({\boldsymbol{\theta}})$
Then any sequence $\\{{\boldsymbol{\theta}}^{n}\\}_{n\in\mathbb{N}}$ defined
by ${\boldsymbol{\theta}}^{n+1}\in\psi_{I(n)}({\boldsymbol{\theta}}^{n})$ has
all of its limit points in $\Gamma$.

###  5.2 Main result

In this section, we derive our convergence result for SLK-Means, but the
result and proof are very similar for SLK-modes. We frame an inner iteration
of the ${\mathbf{S}}$-update defined by (11) and the whole
${\mathbf{M}}$-update in (20) as two distinct point-to-set maps:

$\displaystyle\psi_{{\mathbf{S}}}({\mathbf{S}},{\mathbf{M}})$
$\displaystyle=\operatorname*{arg\,min}_{{\mathbf{S}}}\mathcal{A}({\mathbf{S}},{\mathbf{M}})\times\\{{\mathbf{M}}\\}$
(24) $\displaystyle\psi_{{\mathbf{M}}}({\mathbf{S}},{\mathbf{M}})$
$\displaystyle=\\{{\mathbf{S}}\\}\times\operatorname*{arg\,min}_{{\mathbf{M}}}\mathcal{R}({\mathbf{S}},{\mathbf{M}})$
(25)

$\psi_{{\mathbf{M}}}$ and $\psi_{{\mathbf{S}}}$ are the building blocks of our
SLK method: a full ${\mathbf{S}}$-update simply consists in sequentially
applying $\psi_{{\mathbf{S}}}$ until inner convergence. On the other hand, the
${\mathbf{M}}$-update is captured by a single application of
$\psi_{{\mathbf{M}}}$. Finally, we define $\Gamma$ as the intersection of the
sets of fixed points of $\psi_{\mathbf{S}}$ and $\psi_{\mathbf{M}}$:

$\displaystyle\Gamma=\\{({\mathbf{S}},{\mathbf{M}})~{}:~{}({\mathbf{S}},{\mathbf{M}})=\psi_{\mathbf{S}}({\mathbf{S}},{\mathbf{M}})=\psi_{\mathbf{M}}({\mathbf{S}},{\mathbf{M}})\\}$
(26)

######  Proposition 3.

Starting from some initial point $\\{{\mathbf{S}}^{0},{\mathbf{M}}^{0}\\}$
such that $\mathcal{R}({\mathbf{S}}^{0},{\mathbf{M}}^{0})<\infty$, every limit
point of the sequence $\\{{\mathbf{S}}^{i},{\mathbf{M}}^{i}\\}$ built using
the SLK-Means algorithm 1 is in $\Gamma$.

###### Proof.

The idea of the proof is to apply Corollary 1 with the right ingredients.
First, we define the set of point-to-set maps as
$\\{\psi_{{\mathbf{M}}},\psi_{{\mathbf{S}}}\\}$, and the index function $I$ as
the one given by the SLK method. We now define our optimization-variable
space, and show its compactness. Given that the objective keeps decreasing,
$({\mathbf{S}}^{i},{\mathbf{M}}^{i}),~{}i\in\mathbb{N}$ has to live in the
following space:

$\displaystyle\Theta=L_{0}\times\Delta_{K}$ (27)

where $L_{0}=\\{{\mathbf{M}}\in\mathbb{R}^{K\times
d}\,|\,\mathcal{E}({\mathbf{S}}^{0},{\mathbf{M}})\leq\mathcal{E}({\mathbf{S}}^{0},{\mathbf{M}}^{0})\\}$
represents a sublevel set of $\mathcal{E}({\mathbf{S}}^{0},.)$.

Compactness of $\Theta$: Provided the continuity of $\mathcal{E}$, $L_{0}$ is
closed as the inverse image of a continuous function. Furthermore, it is
bounded as
$\mathcal{E}({\mathbf{S}}^{0},{\mathbf{M}})\xrightarrow[\lVert{\mathbf{M}}\rVert\to\infty]{}\infty$.
Therefore, $L_{0}$ is compact. Furthermore, simplex $\Delta_{K}$ is a compact
subset of $\mathbb{R}^{K}$. Hence, as the product of two compact sets,
$\Theta$ is compact.

Assumption (1):  To prove the closedness assumption (1), we will use the
following lemma from [22]:

######  Lemma 1.

(Proposition 7 of Appendix A in [22]) For a continuous function $\phi:A\times
B\rightarrow B$ defined as

$\displaystyle\psi({\mathbf{a}})=\operatorname*{arg\,min}_{{\mathbf{b}}\in
B}\phi({\mathbf{a}},{\mathbf{b}})$ (28)

Then $\psi$ is closed at ${\mathbf{a}}$ if $\psi({\mathbf{a}})$ is nonempty.

Using the continuity of both $\mathcal{A}$ and $\mathcal{R}$ in both
${\mathbf{S}}$ and ${\mathbf{M}}$, along with Lemma 1, the closedness of
$\psi_{\mathbf{S}}$ and $\psi_{\mathbf{M}}$ follows.

Assumption (2): We now prove that the objective function strictly decreases
for non-stationnary points. The argument relies on the strict convexity of
$\mathcal{A}({\mathbf{S}},{\mathbf{M}})$ in ${\mathbf{S}}$ and of
$\mathcal{R}({\mathbf{S}},{\mathbf{M}})$ in ${\mathbf{S}}$ and ${\mathbf{M}}$.
Specifically, consider $n\in\mathbb{N}$ and
$\\{{\mathbf{S}},{\mathbf{M}}\\}\notin\Gamma$;
$\exists~{}\psi\in\\{\psi_{{\mathbf{S}}},\psi_{{\mathbf{M}}}\\}\,|\,\\{{\mathbf{S}},{\mathbf{M}}\\}\neq\psi({\mathbf{S}},{\mathbf{M}})$.
SLK pre-defined sequence ensures that such $\psi$ is chosen at iteration $n$.
Without loss of generality, assume $\psi=\psi_{I(n)}=\psi_{{\mathbf{S}}}$.
Then, recall that $\mathcal{A}({\mathbf{S}},{\mathbf{M}})$ is strictly convex
in ${\mathbf{S}}$, which implies that it has a unique minimum. Therefore,
either the current solution ${\mathbf{S}}^{i}$ is already a stationnary point,
with
${\mathbf{S}}^{n}=\operatorname*{arg\,min}_{\mathbf{S}}\mathcal{A}({\mathbf{S}}^{n},{\mathbf{M}}^{n})$,
which implies $\\{{\mathbf{S}}^{n},{\mathbf{M}}^{n}\\}\in\Gamma$ (impossible),
or the current solution is not a stationary point. In this case, for
${\mathbf{S}}^{n+1}\in\psi_{{\mathbf{S}}}({\mathbf{S}}^{n},{\mathbf{M}}^{n})$,
we inevitably have
$\mathcal{A}({\mathbf{S}}^{n+1},{\mathbf{M}}^{n})<\mathcal{A}({\mathbf{S}}^{n},{\mathbf{M}}^{n})$.

∎

TABLE II: Average accuracy (in %) for the few-shot learning experiments on the miniImageNet, tieredImageNet and CUB dataset. The best results are highlighted in bold font. The shaded areas are the results based on clustering approaches for transductive few-shot inference. |  |  |  | miniImageNet | tieredImageNet | CUB
---|---|---|---|---|---|---
Methods | Reference | Transductive | Network | 1-shot | 5-shot | 1-shot | 5-shot | 1-shot | 5-shot
Baseline | ICLR 2019[11] | ✗ | ResNet-18 | 51.87 | 75.68 | - | - | 67.02 | 83.58
RelationNet | CVPR 2018 [61] | ResNet-18 | 52.48 | 69.83 | - | - | 68.58 | 84.05
MatchingNet | NeurIPS 2016 [67] | ResNet-18 | 52.91 | 68.88 | - | - | 73.49 | 84.45
ProtoNet | NeurIPS 2017 [57] | ResNet-18 | 54.16 | 73.68 | - | - | 72.99 | 86.64
MTL | CVPR 2019 [60] | ResNet-12 | 61.20 | 75.50 | - | - | - | -
vFSL | ICCV 2019 [84] | ResNet-12 | 61.23 | 77.69 | - | - | - | -
MetaoptNet | CVPR 2019 [35] | ResNet-12 | 62.64 | 78.63 | 65.99 | 81.56 |  |
Neg-cosine | ECCV 2020 [37] | ResNet-18 | 62.33 | 80.94 | 69.68 | 84.56 | - | -
SimpleShot | arXiv 2019 [72] | ResNet-18 | 63.10 | 79.92 | 69.68 | 84.56 | 70.28 | 86.37
Distil | ECCV 2020 [65] | ResNet-18 | 64.82 | 82.14 | 71.52 | 86.03 | - | -
MAML | ICML 2017 [15] |  | ResNet-18 | 49.61 | 65.72 | - | - | 68.42 | 83.47
RelationNet + T | NeurIPS 2019 [25] |  | ResNet-12 | 52.40 | 65.36 | - | - | - | -
ProtoNet + T | NeurIPS 2019 [25] |  | ResNet-12 | 55.15 | 71.12 | - | - | - | -
MatchingNet + T | NeurIPS 2019 [25] |  | ResNet-12 | 56.31 | 69.80 | - | - | - | -
Team | ICCV 2019 [49] |  | ResNet-18 | 60.07 | 75.90 | - | - | - | -
TPN | ICLR 2019 [77] |  | ResNet-12 | 59.46 | 75.64 | - | - | - | -
Entropy-min | ICLR 2020 [13] |  | ResNet-12 | 62.35 | 74.53 | - | - | - | -
DPGN | CVPR 2020 [78] |  | ResNet-18 | 66.63 | 84.07 | 70.46 | 86.44 | - | -
CAN+T | NeurIPS 2019 [25] |  | ResNet-12 | 67.19 | 80.64 | 73.21 | 84.93 | - | -
LaplacianShot | ICML 2020 [90] |  | ResNet-18 | 70.74 | 82.33 | 77.60 | 86.25 | 79.93 | 88.59
K-means |  |  | ResNet-18 | 71.11 | 81.87 | 77.94 | 85.99 | 80.30 | 88.51
K-modes |  |  | ResNet-18 | 72.31 | 82.21 | 78.85 | 86.12 | 81.73 | 88.58
SLK-MEANS |  |  | ResNet-18 | 72.58 | 82.92 | 79.22 | 86.50 | 81.40 | 88.61
SLK-MS |  | ✓ | ResNet-18 | 73.10 | 82.82 | 79.99 | 86.55 | 81.88 | 88.55
LEO | ICLR 2019 [53] | ✗ | WRN | 61.76 | 77.59 | 66.33 | 81.44 | - | -
ProtoNet | NeurIPS 2017 [57] | WRN | 62.60 | 79.97 | - | - | - | -
MatchingNet | NeurIPS 2016 [67] | WRN | 64.03 | 76.32 | - | - | - | -
CC+rot | ICCV 2019 [17] | WRN | 62.93 | 79.87 | 70.53 | 84.98 | - | -
FEAT | CVPR 2020 [79] | WRN | 65.10 | 81.11 | 70.41 | 84.38 | - | -
SimpleShot | arXiv 2019 [72] | WRN | 65.87 | 82.09 | 70.90 | 85.76 | - | -
AWGIM | CVPR 2020 [23] |  | WRN | 63.12 | 78.40 | 67.69 | 82.82 | - | -
Entropy-min | ICLR 2020 [13] |  | WRN | 65.73 | 78.40 | 73.34 | 85.50 | - | -
SIB | ICLR 2020 [26] |  | WRN | 70.0 | 79.2 | 70.90 | 85.76 | - | -
BD-CSPN | ECCV 2020 [38] |  | WRN | 70.31 | 81.89 | 78.74 | 86.92 | - | -
SIB+E3BM | ECCV 2020 [39] |  | WRN | 71.4 | 81.2 | 75.6 | 84.3 | - | -
LaplacianShot | ICML 2020 [90] |  | WRN | 73.44 | 83.93 | 78.80 | 87.72 | - | -
K-means |  |  | WRN | 73.80 | 84.62 | 79.78 | 87.23 | - | -
K-modes |  |  | WRN | 74.78 | 84.45 | 80.67 | 87.23 | - | -
SLK-MEANS |  |  | WRN | 74.75 | 84.61 | 80.55 | 87.57 | - | -
SLK-MS |  | ✓ | WRN | 75.17 | 84.28 | 81.13 | 87.69 | - | -

##  6 Few-shot learning experiments

In the following, we describe our experimental setup and report the results of
various clustering approaches in the few-shot learning scenario.

###  6.1 Datasets

For this task, we resort to three widely used benchmarks for few-shot
classification: miniImageNet, tieredImageNet and CUB-200-2011. The details of
these standard few-shot benchmarks are given below.

miniImageNet is a subset of the larger ILSVRC-12 dataset [52]. It has a total
of 60,000 color images with 100 classes, where each class has 600 images of
size $84\times 84$, following [67]. We use the standard split of 64 classes
for base training, 16 for validation and 20 for testing [50, 72].

tieredImageNet [51] is also a subset of the ILSVRC-12 dataset, but with 608
classes instead. Following [72], we split the dataset into 351 classes for
base training, 97 for validation and 160 for testing. As in the miniImageNet
dataset, the images are resized to $84\times 84$ pixels.

CUB-200-2011 [70] is a fine-grained image classification dataset with 200
categories. Following the setting in [11] for few-shot classification, we
split the CUB dataset into 100 classes for base training, 50 for validation
and 50 for testing. Similarly to the previous datasets, the images are resized
to $84\times 84$ pixels.

###  6.2 Evaluation Protocol

We follow standard evaluation practices in few-shot learning, using 600 tasks
(5-way classification) randomly sampled from the test classes. The query set
of each task contains 15 images per class, and the average accuracy over 600
few-shot tasks is reported, for both the 1-shot and 5-shot supervision
scenarios.

###  6.3 Architectures

We evaluate all clustering methods on two different network models as feature
extractor function $f_{\theta}$, following the same settings used recently in
[90]. Particularly, we consider the following standard models for the few-shot
learning experiments:

ResNet-18 is based on the standard deep residual network, ResNet [24], whose
first two down-sampling layers are removed, by setting the stride to 1 in the
first convolutional layer and removing the first max-pooling layer.
Furthermore, the kernel size of the first convolutional layer is set to
$3\times 3$, instead of $7\times 7$ in the original model. The used ResNet-18
model has 8 basic residual blocks, and the dimension of the extracted features
is 512.

WRN [82] widens the residual blocks by adding more convolutional layers and
feature planes. In our particular case, we used an architecture with 28
convolutional layers, along with a widening factor of 10. This network yields
an extracted feature vector of dimension 640.

###  6.4 Implementation Details

Base training:  We utilize the publicly available trained models from
[90]666https://github.com/imtiazziko/LaplacianShot. The network models are
trained using the standard cross-entropy loss on the training set of $A$ base
classes, without resorting to meta-learning and episodic-training strategies.
Training starts from scratch and runs during 90 epochs. The models are
optimized using SGD. Similarly to [90, 72], the initial learning rate is set
to 0.1 and scaled by 0.1 at the 45th and 66th epoch for miniImageNet and CUB,
and after every 30 epochs for tieredImageNet. The mini-batch size is set to
256 for ResNet-18 and 128 for WRN.

Inference procedure: Before running any clustering approach in the few-shot
learning setting, we extract the features from both the support and query
samples with the initial $f_{\theta}$ trained on the base classes. Similarly
to [72, 90], we perform CL2 normalization on the features. This amounts to
computing the mean of the base-class features, followed by centering the
extracted support and query features via subtracting the mean. Then, we
perform an L2 normalization on the centered features. Similarly to [38, 90],
we also perform a bias correction by computing the cross-class bias from the
difference between the mean features of the support set
$\mathcal{X}_{\text{s}}$ and the query set $\mathcal{X}_{\text{q}}$. Then, we
add this bias to each query features vector before performing the clustering
for a given few-shot task. For each task, the initial prototypes
$\\{{\mathbf{m}}_{k}^{0}\\}_{k=1}^{K}$ for the clustering algorithms (see Alg.
1) are either the support examples
${\mathbf{x}}_{p}\in\mathcal{X}_{\text{s}}^{k}$ of each class in the 1-shot
setting, or the mean of the support examples within each class in the 5-shot
setting. Once the clustering of a task is performed, we label each query
sample according to the class label of the support examples belonging to the
same cluster.

SLK hyperparameters:  We define $\mathcal{N}_{p}^{\rho}$ as the set of the
$\rho$ nearest neighbors of data point ${\mathbf{x}}_{p}$. Mode estimation is
based on the Gaussian kernel as per (14), with $\sigma^{2}$ estimated as:

$\displaystyle\sigma^{2}=\frac{1}{N\rho}\sum_{\mathbf{x}_{p}\in\mathbf{X}}\sum_{\mathbf{x}_{q}\in\mathcal{N}_{p}^{\rho}}\|\mathbf{x}_{p}-\mathbf{x}_{q}\|^{2}$
(29)

The Laplacian affinities $w_{\mathcal{L}}$ are built using the $\rho$-nearest
neighbors:

$\displaystyle
w_{\mathcal{L}}({\mathbf{x}}_{p},{\mathbf{x}}_{q})=\begin{cases}1&\text{if
}\mathbf{x}_{q}\in\mathcal{N}_{p}^{\rho}\\\ 0&\text{otherwise}\end{cases}$
(30)

This yields a sparse affinity matrix, which is efficient in terms of memory
and computational complexity. In our experiments, $\rho$ is simply chosen from
three typical values (3, 5 or 10) tuned over training-data by sampling 500
few-shot tasks from the training set. Such procedure yielded a fixed $\rho=3$
during inference. Since no test data is involved in choosing the value of
$\rho$, this procedure is not biased towards the testing classes. Then, we
tune parameter $\lambda$ on the validation classes by sampling $500$ few-shot
tasks and evaluating a range of $\lambda$ values from
$[0.1,0.3,0.5,0.7,0.8,1.0]$. The best $\lambda$ values corresponding to the
best average 1-shot and 5-shot accuracy over the 500 tasks on the validation
classes are selected for inference on the test classes.

TABLE III:  Average inference time (in seconds) for the 5-shot tasks. The results are for miniImageNet dataset with WRN network. Methods | inference time
---|---
SimpleShot [72] | 0.009
Transductive tuning [13] | 20.7
K-means | 1.17
K-modes | 1.10
SLK-MEANS | 1.71
SLK-MS | 1.15

###  6.5 Results

The results on miniImagenet, tieredImagenet and CUB datasets are reported in
Table II. It is very interesting to observe that the clustering-based
approaches (shaded gray in the tables) consistently outperform transductive
few-shot learning methods, yielding state-of-the-art on these standard
benchmarks, without resorting to complex meta-learning and episodic-training
strategies. In the 1-shot settings, our SLK-MS clustering outperforms the
best-performing state-of-the-art few-shot method by important margins (2 to
$3\%$), across data sets and models. It also outperforms several convoluted
meta-learning methods by substantial margins, up to $10\%$ or more, including
many very recent tranductive methods. Surprisingly, our results show that even
the standard K-means clustering procedure already achieves competitive
performances in comparison to the state-of-the-art. The standard K-means
procedure is on par with the best-performing method in the tranductive
setting, and outperforms substantially several very recent methods based on
complex meta-learning and/or episodic training schemes. These surprising
results point to the limitations of the current few-shot benchmarks, and
question the viability of a large body of convoluted few-shot learning
techniques in the recent literature. It also suggests rethinking transductive
few-shot learning and using simple constrained clustering mechanisms as
baselines.

###  6.6 Inference time

Table III reports an evaluation in terms of efficiency. This is done by
measuring the average time required for inference on the 5-shot tasks in
miniImageNet. Note that all the methods listed in Table III are transductive,
except SimpleShot. The transductive network fine-tuning approach in [13], for
instance, is significantly slower than the other methods. This is due to the
gradient updates performed by [13] at inference over all the network
parameters. While the clustering methods are more efficient than the approach
in [13], they require longer inference times than inductive predictions [72].

##  7 Experiments on clustering

We conducted a comprehensive experimental analysis of our SLK approaches on 5
different data sets, each having different features. To evaluate the
clustering performance, we use two metrics widely adopted in the clustering
community: the Normalized Mutual Information (NMI) [59] and the Clustering
Accuracy (ACC). The optimal mapping of clustering assignments to the true
labels is determined by using the Kuhn-Munkres algorithm [43]. We further
report evaluations and comparisons in terms of optimization quality, which we
measure by the values obtained at convergence for the discrete-variable
objective ${\cal E}$ in (2).

TABLE IV: Clustering results as NMI/ACC. Algorithm | MNIST | MNIST (GAN) | LabelMe (Alexnet) | LabelMe (GIST) | YTF | Shuttle | Reuters
---|---|---|---|---|---|---|---
K-means | 0.53/0.55 | 0.68/0.75 | 0.81/0.90 | 0.57/0.69 | 0.77/0.58 | 0.22/0.41 | 0.48/0.73
K-modes | 0.56/0.60 | 0.69/0.80 | 0.81/0.91 | 0.58/0.68 | 0.79/0.62 | 0.33/0.47 | 0.48/0.72
NCUT | 0.74/0.61 | 0.77/0.67 | 0.81/0.91 | 0.58/0.61 | 0.74/0.54 |  0.47/0.46 | -
KK-means | 0.53/0.55 | 0.69/0.68 | 0.81/0.90 | 0.57/0.63 | 0.71/0.50 | 0.26/0.40 | -
LK | - | - | 0.81/0.91 | 0.59/0.61 | 0.77/0.59 | - | -
Spectralnet | - | - | - | - | - | - | 0.46/0.65
SLK-MS |  0.80/0.79 | 0.86/0.94 |  0.83/0.91 |  0.61/0.72 |  0.82/0.65 | 0.45/0.70 | 0.43/0.74
SLK-Means | 0.78/0.75 | 0.84/0.91 |  0.83/0.91 | 0.60/0.68 |  0.82/0.64 | 0.31/0.71 |  0.56/0.77

(a) MNIST (small)

(b) LabelMe (Alexnet)

Figure 1: Discrete-variable objective (2): Comparison of the objectives
obtained at convergence for SLK-MS (ours) and LK [71]. The objectives at
convergence are plotted versus different values of parameter $\lambda$.

We resort to well-known benchmarks to report comprehensive evaluations of the
proposed algorithm777Code is available at: https://github.com/imtiazziko/SLK,
as well as comparisons to the following related baseline methods: Laplacian
K-modes (LK) [71], K-means, NCUT [56], K-modes [54, 8], Kernel K-means (KK-
means) [14, 62] and Spectralnet [55].

###  7.1 Clustering datasets

TABLE V: Datasets used in the clustering experiments. Datasets | $\text{Samples}~{}(N)$ | $\text{Dimensions}~{}(D)$ | $\text{Clusters}~{}(L)$ | Imbalance
---|---|---|---|---
MNIST (small) | $2,000$ | $784$ | $10$ | $1$
MNIST | $70,000$ | $784$ | $10$ | $\sim 1$
MNIST (GAN) | $70,000$ | $256$ | $10$ | $\sim 1$
Shuttle | $58,000$ | $9$ | $7$ | $4,558$
LabelMe (Alexnet) | $2,688$ | $4,096$ | $8$ | $1$
LabelMe (GIST) | $2,688$ | $44,604$ | $8$ | $1$
YTF | $10,036$ | $9,075$ | $40$ | $13$
Reuters (code) | $685,071$ | $10$ | $4$ | $\sim 5$

We used a mix of image datasets (MNIST, LabelMe, YTF) and non-image datasets
(Reuters, Shuttle) for the clustering task. The overall summary of the
datasets is given in Table V. For each dataset, imbalance is defined as the
ratio between the size of the largest and smallest clusters. We use three
versions of the MNIST dataset [34], which are detailed below. MNIST is the
original dataset and contains the whole set of $70,000$ images. In order to
compare to LK [71], which does not scale up for large datasets, we use MNIST
(small). This reduced version includes only $2,000$ images by randomly
sampling $200$ images per class. To obtain the MNIST (GAN) dataset, we train a
GAN following [19] on $60,000$ training images and extract the
$256$-dimensional features from the discriminator network for the whole set of
$70,000$ images. We use the publicly available auto-encoder in [28] to extract
$10$-dimensional features to generate the Reuters (code) dataset. LabelMe [46]
consists of $2,688$ images divided into $8$ categories. Furthermore, we used
the pre-trained AlexNet [32] and extracted the 4096-dimensional features from
the fully-connected layer, leading to the LabelMe (AlexNet) dataset. To show
the performances on high-dimensional data, we extract $44604$-dimensional GIST
features [46] for the LabelMe dataset, which is referred to as LabelMe(GIST).
Finally, the Youtube Faces (YTF) [75] dataset consists of videos of faces with
$40$ different subjects.

###  7.2 Implementation details

Similarly to the few-shot experiments, we used $\rho$-nearest neighbor
affinities to build $w_{\mathcal{L}}$, with fixed $\rho=5$ for all the
datasets. For the large datasets, such as MNIST, Shuttle and Reuters, we used
the _Flann_ library [42] with the KD-tree algorithm, which finds approximate
nearest neighbors. For the other smaller datasets, we used an efficient
implementation of exact nearest-neighbor computations. We used the same sparse
${\mathbf{W}}_{\mathcal{L}}$ matrix for the pairwise-affinity algorithms we
compared with, i.e., NCUT, KK-means and Laplacian K-modes. Furthermore, for
each of these baseline methods, we evaluated the default setting of affinity
construction with tuned $\sigma$, and report the best result found. Initial
centers $\\{{\mathbf{m}}^{0}_{k}\\}_{k=1}^{K}$ are based on K-means++ seeds
[2]. We choose the best initial seed and regularization parameter $\lambda$
empirically based on the accuracy over a validation set (10% of the total
data). The $\lambda$ is determined by tuning over a small range from $1$ to
$4$. In Algorithm 1, all assignment variables ${\mathbf{s}}_{p}$ are updated
in parallel. We run the publicly released codes for K-means [48], NCUT [56],
Laplacian K-modes [7], KK-means and Spectralnet [55].

###  7.3 Clustering results

Table IV reports the clustering results, showing that, in most of the cases,
our algorithms SLK-MS and SLK-Means yielded the best NMI and ACC values. For
MNIST with the raw intensities as features, the proposed SLK-MS achieved
almost 80% NMI and ACC, outperforming all the standard clustering methods by
substantial margins. With better learned GAN features for MNIST (GAN), the
accuracy (ACC) increases up to 95%, also outperforming standard clustering
methods by substantial margins. For the Reuters (code) datasets, we used the
same features and Euclidean distance based affinity as Spectralnet [55], and
obtained better NMI/ACC performances. The Shuttle dataset is quite imbalanced
and, therefore, all the baseline clustering methods fail to achieve high
accuracy. Notice that, in regard to ACC for the Shuttle dataset, our SLK
methods outperformed all the other methods by a large margin (over $20\%$ of
difference).

###  7.4 Comparison in terms of optimization quality

To assess the optimization quality of our optimizer, we computed the values of
discrete-variable objective ${\cal E}$ in model (2) at convergence for our
concave-convex relaxation (SLK-MS) as well as for the convex relaxation in
[71] (LK). We compare the discrete-variable objectives for different values of
$\lambda$. For a fair comparison, we use the same initialization, $\sigma$,
$w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{x}}_{q})$, $\lambda$ and mean-shift
(MS) modes for both methods. As shown in the plots in Figure 1, and reported
in Table VI, our relaxation consistently obtained lower values of discrete-
variable objective ${\cal E}$ at convergence than the convex relaxation in
[71].

TABLE VI: Discrete-variable objectives at convergence for LK [71] and SLK-MS. Datasets |  LK [71] | SLK-MS (ours)
---|---|---
MNIST (small) | $273.25$ | $67.09$
LabelMe (Alexnet) | $-1.513\,84\text{\times}{10}^{03}$ | $-1.807\,77\text{\times}{10}^{03}$
LabelMe (GIST) | $-1.954\,90\text{\times}{10}^{03}$ | $-2.024\,10\text{\times}{10}^{03}$
YTF | $-1.000\,32\text{\times}{10}^{4}$ | $-1.000\,35\text{\times}{10}^{4}$

##  8 Conclusion

We presented a general formulation for clustering and transductive few-shot
learning, integrating Laplacian regularization, prototype-based terms and
supervision constraints. Based on a concave-convex relaxation, our bound
optimizer performs parallel point-wise updates, while guaranteeing
convergence. We reported comprehensive clustering and few-shot experiments
over various data sets, showing that our methods yield competitive
performances, in regard to both accuracy and optimization quality. Using
standard training on the base classes, and without resorting to complex meta-
learning and episodic-training strategies, our clustering methods yield state-
of-the-art results in the transductive few-shot setting, outperforming
significantly a large number of convoluted methods in the recent literature,
across various models, settings and data sets. In the 1-shot settings, our
SLK-MS clustering outperforms the best-performing state-of-the-art few-shot
method by up to $3\%$. It also outperforms several convoluted meta-learning
methods by substantial margins, up to $10\%$ or more. Surprisingly, we found
that even the basic K-means already achieves competitive performances in
comparison to the state-of-the-art in few-shot learning, outperforming
substantially several very recent methods. These surprising performances point
to the limitations of the current few-shot benchmarks, question the viability
of a large body of complex few-shot learning works in the literature, and
suggest simple clustering mechanisms as baselines.

## Acknowledgements

This research was supported by the National Science and Engineering Research
Council of Canada (NSERC), via its Discovery Grant program (grant CRSNG RGPIN
2019-05954). It was also supported by computational resources provided by
Compute Canada and Calcul Quebec.

## References

  * [1] D. Aloise, A. Deshpande, P. Hansen, and P. Popat.  Np-hardness of euclidean sum-of-squares clustering.  Machine Learning, 75(2):245–248, 2009.
  * [2] D. Arthur and S. Vassilvitskii.  k-means++: The advantages of careful seeding.  In ACM-SIAM symposium on Discrete algorithms, pages 1027–1035. Society for Industrial and Applied Mathematics, 2007.
  * [3] M. Belkin, P. Niyogi, and V. Sindhwani.  Manifold regularization: A geometric framework for learning from labeled and unlabeled examples.  Journal of Machine Learning Research, 7:2399–2434, 2006.
  * [4] D. Berthelot, N. Carlini, I. Goodfellow, N. Papernot, A. Oliver, and C. A. Raffel.  Mixmatch: A holistic approach to semi-supervised learning.  In Advances in Neural Information Processing Systems (NeurIPS), 2019\.
  * [5] M. Boudiaf, Z. I. Masud, J. Rony, J. Dolz, P. Piantanida, and I. B. Ayed.  Transductive information maximization for few-shot learning.  Neural Information Processing Systems (NeurIPS), 2020.
  * [6] J. Bronskill, J. Gordon, J. Requeima, S. Nowozin, and R. E. Turner.  Tasknorm: Rethinking batch normalization for meta-learning.  In International Conference on Machine Learning (ICML), 2020.
  * [7] M. Á. Carreira-Perpiñán.  Gaussian mean-shift is an em algorithm.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 29(5):767–776, 2007.
  * [8] M. A. Carreira-Perpinán.  A review of mean-shift algorithms for clustering.  arXiv preprint arXiv:1503.00687, 2015.
  * [9] M. Á. Carreira-Perpiñán and W. Wang.  The k-modes algorithm for clustering.  arXiv preprint arXiv:1304.6478, 2013.
  * [10] O. Chapelle, B. Schlkopf, and A. Zien.  Semi-Supervised Learning.  The MIT Press, 1st edition, 2010.
  * [11] W.-Y. Chen, Y.-C. Liu, Z. Kira, Y.-C. F. Wang, and J.-B. Huang.  A closer look at few-shot classification.  In International Conference on Learning Representations (ICLR), 2019\.
  * [12] D. Comaniciu and P. Meer.  Mean shift: A robust approach toward feature space analysis.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 24(5):603–619, 2002.
  * [13] G. S. Dhillon, P. Chaudhari, A. Ravichandran, and S. Soatto.  A baseline for few-shot image classification.  In International Conference on Learning Representations (ICLR), 2020\.
  * [14] I. S. Dhillon, Y. Guan, and B. Kulis.  Kernel k-means: spectral clustering and normalized cuts.  In International Conference on Knowledge Discovery and Data Mining (SIGKDD), 2004.
  * [15] C. Finn, P. Abbeel, and S. Levine.  Model-agnostic meta-learning for fast adaptation of deep networks.  In International Conference on Machine Learning (ICML), 2017.
  * [16] J. C. Fiorot and P. Huard.  Composition and union of general algorithms of optimization.  In Point-to-Set Maps and Mathematical Programming, pages 69–85. Springer, 1979.
  * [17] S. Gidaris, A. Bursuc, N. Komodakis, P. Pérez, and M. Cord.  Boosting few-shot visual learning with self-supervision.  In IEEE/CVF International Conference on Computer Vision (ICCV), 2019\.
  * [18] Y. Gong, M. Pawlowski, F. Yang, L. Brandy, L. Bourdev, and R. Fergus.  Web scale photo hash clustering on a single machine.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2015.
  * [19] I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio.  Generative adversarial nets.  In Neural Information Processing Systems (NeurIPS), 2014.
  * [20] J. Gordon, J. Bronskill, M. Bauer, S. Nowozin, and R. Turner.  Meta-learning probabilistic inference for prediction.  In International Conference on Learning Representations (ICLR), 2019\.
  * [21] Y. Grandvalet and Y. Bengio.  Semi-supervised learning by entropy minimization.  In Neural Information Processing Systems (NeurIPS), 2005.
  * [22] A. Gunawardana, W. Byrne, and M. I. Jordan.  Convergence theorems for generalized alternating minimization procedures.  Journal of Machine Learning Research, 6(12), 2005.
  * [23] Y. Guo and N.-M. Cheung.  Attentive weights generation for few shot learning via information maximization.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
  * [24] K. He, X. Zhang, S. Ren, and J. Sun.  Deep residual learning for image recognition.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition, 2016.
  * [25] R. Hou, H. Chang, M. Bingpeng, S. Shan, and X. Chen.  Cross attention network for few-shot classification.  In Advances in Neural Information Processing Systems (NeurIPS), 2019\.
  * [26] S. X. Hu, P. G. Moreno, Y. Xiao, X. Shen, G. Obozinski, N. D. Lawrence, and A. Damianou.  Empirical bayes transductive meta-learning with synthetic gradients.  In International Conference on Learning Representations (ICLR), 2020\.
  * [27] A. Iscen, G. Tolias, Y. Avrithis, and O. Chum.  Label propagation for deep semi-supervised learning.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition(CVPR), 2019.
  * [28] Z. Jiang, Y. Zheng, H. Tan, B. Tang, and H. Zhou.  Variational deep embedding: An unsupervised and generative approach to clustering.  In International Joint Conference on Artificial Intelligence (IJCAI), 2017.
  * [29] T. Joachims.  Transductive inference for text classification using support vector machines.  In International Conference on Machine Learning (ICML), 1999.
  * [30] P. Krähenbühl and V. Koltun.  Efficient inference in fully connected crfs with gaussian edge potentials.  In Neural Information Processing Systems (NeurIPS), 2011.
  * [31] P. Krähenbühl and V. Koltun.  Parameter learning and convergent inference for dense random fields.  In International Conference on Machine Learning (ICML), 2013.
  * [32] A. Krizhevsky, I. Sutskever, and G. E. Hinton.  Imagenet classification with deep convolutional neural networks.  In Neural Information Processing Systems (NeurIPS), 2012.
  * [33] K. Lange, D. R. Hunter, and I. Yang.  Optimization transfer using surrogate objective functions.  Journal of Computational and Graphical Statistics, 9(1):1–20, 2000\.
  * [34] Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [35] K. Lee, S. Maji, A. Ravichandran, and S. Soatto.  Meta-learning with differentiable convex optimization.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2019.
  * [36] D.-M. Liang and Y.-F. Li.  Lightweight label propagation for large-scale network data.  In International Joint Conference on Artificial Intelligence (IJCAI), 2018.
  * [37] B. Liu, Y. Cao, Y. Lin, Q. Li, Z. Zhang, M. Long, and H. Hu.  Negative margin matters: Understanding margin in few-shot classification.  In European Conference on Computer Vision (ECCV), 2020.
  * [38] J. Liu, L. Song, and Y. Qin.  Prototype rectification for few-shot learning.  In European Conference on Computer Vision (ECCV), 2020.
  * [39] Y. Liu, B. Schiele, and Q. Sun.  An ensemble of epoch-wise empirical bayes for few-shot learning.  In European Conference on Computer Vision (ECCV), 2020.
  * [40] N. Mishra, M. Rohaninejad, X. Chen, and P. Abbeel.  A simple neural attentive meta-learner.  In International Conference on Learning Representations (ICLR), 2018\.
  * [41] T. Miyato, S. Maeda, M. Koyama, and S. Ishii.  Virtual adversarial training: A regularization method for supervised and semi-supervised learning.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 41(8):1979–1993, 2019.
  * [42] M. Muja and D. G. Lowe.  Scalable nearest neighbor algorithms for high dimensional data.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 36(11):2227–2240, 2014.
  * [43] J. Munkres.  Algorithms for the assignment and transportation problems.  Journal of the Society for Industrial and Applied Mathematics, 5(1):32–38, 1957.
  * [44] M. Narasimhan and J. Bilmes.  A submodular-supermodular procedure with applications to discriminative structure learning.  In Conference on Uncertainty in Artificial Intelligence (UAI), pages 404–412, 2005.
  * [45] J. Newling and F. Fleuret.  Nested mini-batch k-means.  In Neural Information Processing Systems (NeurIPS), 2016.
  * [46] A. Oliva and A. Torralba.  Modeling the shape of the scene: A holistic representation of the spatial envelope.  International Journal of Computer Vision, 42(3):145–175, 2001.
  * [47] B. Oreshkin, P. R. López, and A. Lacoste.  Tadam: Task dependent adaptive metric for improved few-shot learning.  In Neural Information Processing Systems (NeurIPS), 2018.
  * [48] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  Journal of Machine Learning Research, 12:2825–2830, 2011.
  * [49] L. Qiao, Y. Shi, J. Li, Y. Wang, T. Huang, and Y. Tian.  Transductive episodic-wise adaptive metric for few-shot learning.  In IEEE/CVF International Conference on Computer Vision (ICCV), 2019\.
  * [50] S. Ravi and H. Larochelle.  Optimization as a model for few-shot learning.  In International Conference on Learning Representation (ICLR), 2017\.
  * [51] M. Ren, E. Triantafillou, S. Ravi, J. Snell, K. Swersky, J. B. Tenenbaum, H. Larochelle, and R. S. Zemel.  Meta-learning for semi-supervised few-shot classification.  In International Conference on Learning Representations ICLR, 2018\.
  * [52] O. Russakovsky, J. Deng, H. Su, J. Krause, S. Satheesh, S. Ma, Z. Huang, A. Karpathy, A. Khosla, M. Bernstein, A. C. Berg, and L. Fei-Fei.  ImageNet Large Scale Visual Recognition Challenge.  International Journal of Computer Vision (IJCV), 115(3):211–252, 2015.
  * [53] A. A. Rusu, D. Rao, J. Sygnowski, O. Vinyals, R. Pascanu, S. Osindero, and R. Hadsell.  Meta-learning with latent embedding optimization.  In International Conference on Learning Representations (ICLR), 2019\.
  * [54] M. B. Salah, I. B. Ayed, J. Yuan, and H. Zhang.  Convex-relaxed kernel mapping for image segmentation.  IEEE Transactions on Image Processing, 23(3):1143–1153, 2014.
  * [55] U. Shaham, K. Stanton, H. Li, R. Basri, B. Nadler, and Y. Kluger.  Spectralnet: Spectral clustering using deep neural networks.  In International Conference on Learning Representations (ICLR), 2018\.
  * [56] J. Shi and J. Malik.  Normalized cuts and image segmentation.  IEEE Transactions on Pattern Analysis and Machine Intelligence, 22(8):888–905, 2000.
  * [57] J. Snell, K. Swersky, and R. Zemel.  Prototypical networks for few-shot learning.  In Neural Information Processing Systems (NeurIPS), 2017.
  * [58] B. K. Sriperumbudur and G. R. Lanckriet.  On the convergence of the concave-convex procedure.  In Neural Information Processing Systems (NeurIPS), 2009.
  * [59] A. Strehl and J. Ghosh.  Cluster ensembles—a knowledge reuse framework for combining multiple partitions.  Journal of Machine Learning Research, 3(12):583–617, 2002.
  * [60] Q. Sun, Y. Liu, T. Chua, and B. Schiele.  Meta-transfer learning for few-shot learning.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2019.
  * [61] F. Sung, Y. Yang, L. Zhang, T. Xiang, P. H. Torr, and T. M. Hospedales.  Learning to compare: Relation network for few-shot learning.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2018.
  * [62] M. Tang, D. Marin, I. B. Ayed, and Y. Boykov.  Kernel cuts: Kernel and spectral clustering meet regularization.  International Journal of Computer Vision, 127:477–511, 2019.
  * [63] M. Tang, F. Perazzi, A. Djelouah, I. B. Ayed, C. Schroers, and Y. Boykov.  On regularized losses for weakly-supervised CNN segmentation.  In European Conference on Computer Vision (ECCV), 2018.
  * [64] F. Tian, B. Gao, Q. Cui, E. Chen, and T.-Y. Liu.  Learning deep representations for graph clustering.  In AAAI Conference on Artificial Intelligence, 2014.
  * [65] Y. Tian, Y. Wang, D. Krishnan, J. B. Tenenbaum, and P. Isola.  Rethinking few-shot image classification: a good embedding is all you need?  In European Conference on Computer Vision (ECCV), 2020.
  * [66] V. Vapnik.  An overview of statistical learning theory.  IEEE Transactions on Neural Networks, 10(5):988–999, 1999.
  * [67] O. Vinyals, C. Blundell, T. P. Lillicrap, K. Kavukcuoglu, and D. Wierstra.  Matching networks for one shot learning.  In Neural Information Processing Systems (NeurIPS), 2016.
  * [68] M. Vladymyrov and M. Carreira-Perpiñán.  The variational nystrom method for large-scale spectral problems.  In International Conference on Machine Learning (ICML), 2016.
  * [69] U. Von Luxburg.  A tutorial on spectral clustering.  Statistics and computing, 17(4):395–416, 2007.
  * [70] C. Wah, S. Branson, P. Welinder, P. Perona, and S. Belongie.  The caltech-ucsd birds-200-2011 dataset.  2011\.
  * [71] W. Wang and M. A. Carreira-Perpinán.  The laplacian k-modes algorithm for clustering.  arXiv preprint arXiv:1406.3895, 2014.
  * [72] Y. Wang, W.-L. Chao, K. Q. Weinberger, and L. van der Maaten.  Simpleshot: Revisiting nearest-neighbor classification for few-shot learning.  arXiv preprint arXiv:1911.04623, 2019.
  * [73] Y. Wang, C. Xu, C. Liu, L. Zhang, and Y. Fu.  Instance credibility inference for few-shot learning.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
  * [74] J. Weston, F. Ratle, H. Mobahi, and R. Collobert.  Deep learning via semi-supervised embedding.  In Neural networks: Tricks of the trade, pages 639–655. Springer, 2012.
  * [75] L. Wolf, T. Hassner, and I. Maoz.  Face recognition in unconstrained videos with matched background similarity.  In IEEE/CVF Computer Vision and Pattern Recognition (CVPR), 2011\.
  * [76] C. J. Wu.  On the convergence properties of the em algorithm.  The Annals of statistics, pages 95–103, 1983.
  * [77] L. Yanbin, J. Lee, M. Park, S. Kim, E. Yang, S. Hwang, and Y. Yang.  Learning to propagate labels: Transductive propagation network for few-shot learning.  In International Conference on Learning Representations (ICLR), 2019\.
  * [78] L. Yang, L. Li, Z. Zhang, X. Zhou, E. Zhou, and Y. Liu.  Dpgn: Distribution propagation graph network for few-shot learning.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
  * [79] H.-J. Ye, H. Hu, D.-C. Zhan, and F. Sha.  Few-shot learning via embedding adaptation with set-to-set functions.  In IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
  * [80] J. Yuan, K. Yin, Y. Bai, X. Feng, and X. Tai.  Bregman-proximal augmented lagrangian approach to multiphase image segmentation.  In Scale Space and Variational Methods in Computer Vision (SSVM), pages 524–534, 2017.
  * [81] A. L. Yuille and A. Rangarajan.  The concave-convex procedure (CCCP).  In Neural Information Processing Systems (NeurIPS), 2001.
  * [82] S. Zagoruyko and N. Komodakis.  Wide residual networks.  In British Machine Vision Conference (BMVC), 2016.
  * [83] W. I. Zangwill.  Nonlinear programming: a unified approach, volume 52.  Prentice-Hall Englewood Cliffs, NJ, 1969.
  * [84] J. Zhang, C. Zhao, B. Ni, M. Xu, and X. Yang.  Variational few-shot learning.  In IEEE/CVF International Conference on Computer Vision (ICCV), 2019\.
  * [85] Z. Zhang, J. T. Kwok, and D.-Y. Yeung.  Surrogate maximization/minimization algorithms and extensions.  Machine Learning, 69:1–33, 2007.
  * [86] X. Zhen, H. Sun, Y. Du, J. Xu, Y. Yin, L. Shao, and C. Snoek.  Learning to learn kernels with variational random features.  In International Conference on Machine Learning (ICML), 2020.
  * [87] D. Zhou, O. Bousquet, T. N. Lal, J. Weston, and B. Schölkopf.  Learning with local and global consistency.  In Advances in Neural Information Processing Systems (NeurIPS), 2004\.
  * [88] X. Zhu, Z. Ghahramani, and J. Lafferty.  Semi-supervised learning using gaussian fields and harmonic functions.  In International Conference on Machine Learning (ICML), 2003.
  * [89] I. Ziko, E. Granger, and I. Ben Ayed.  Scalable laplacian k-modes.  In Neural Information Processing Systems (NeurIPS), 2018.
  * [90] I. M. Ziko, J. Dolz, E. Granger, and I. B. Ayed.  Laplacian regularized few-shot learning.  In International Conference on Machine Learning (ICML), 2020.

##  Appendix A Proof of Proposition 1

Proposition 1 states that, given current solution
${\mathbf{S}}^{i,n}=[s_{p,k}^{i,n}]$ at outer iteration $i$, we have the
following auxiliary function (up to an additive constant) for concave-convex
relaxation (5) and psd affinity matrix ${\mathbf{W}}_{\mathcal{L}}$:

${\cal
A}^{i,n}({\mathbf{S}})=\sum_{p=1}^{N}{\mathbf{s}}_{p}^{t}(\log({\mathbf{s}}_{p})-{\mathbf{a}}_{p}^{i}-\lambda{{\mathbf{b}}}_{p}^{i,n})$
(31)

where ${{\mathbf{a}}}_{p}^{i}$ and ${{\mathbf{b}}}_{p}^{i}$ are the following
$K$ dimensional vectors:

$\displaystyle{{\mathbf{a}}}_{p}^{i}$
$\displaystyle=[a_{p,1}^{i},\dots,a_{p,K}^{i}]^{t},\,\mbox{\em
with}\,\,a_{p,k}^{i}=w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}}_{k}^{i})$
(32a) $\displaystyle{{\mathbf{b}}}_{p}^{i}$
$\displaystyle=[b_{p,1}^{i},\dots,b_{p,K}^{i}]^{t},\,\mbox{\em
with}\,\,b_{p,k}^{i}=\sum_{q}w_{\mathcal{L}}({\mathbf{x}}_{p},{\mathbf{x}}_{q})s_{q,k}^{i}$
(32b)

###### Proof.

We want to study the sign of the following difference. Discarding the constant
$2\sum_{p}d_{p}$ in (4) for the Laplacian term ${\cal L}({\mathbf{S}})$, one
can write:

$\displaystyle\mathcal{R}^{i}({\mathbf{S}})-{\cal A}^{i,n}({\mathbf{S}})$
$\displaystyle=\lambda~{}\sum_{p}{\mathbf{s}}_{p}^{t}({\mathbf{b}}_{p}^{i,n}-\sum_{q}w_{\mathcal{L}}({\mathbf{x}}_{p},{\mathbf{x}}_{q}){\mathbf{s}}_{q})$
(33)
$\displaystyle=\lambda~{}\sum_{p,q}{\mathbf{s}}_{p}^{t}w_{\mathcal{L}}({\mathbf{x}}_{p},{\mathbf{x}}_{q})({\mathbf{s}}_{q}^{i}-{\mathbf{s}}_{q})$
(34)

First, instead of considering the $N\times K$ matrix ${\mathbf{S}}$, let us
represent our assignment variables by a flattened vector
${\mathbf{s}}\in[0,1]^{KN}$, that takes the form
${\mathbf{s}}=[{\mathbf{s}}_{1},{\mathbf{s}}_{2},\dots,{\mathbf{s}}_{N}]$.
Recall that each ${\mathbf{s}}_{p}$ is a vector of dimension $K$ containing
the probability variables of all labels for point $p$:
${\mathbf{s}}_{p}=[s_{p,1},\dots,s_{p,K}]^{t}$. Let
$\Psi=-{\mathbf{W}}_{\mathcal{L}}\otimes\mathbf{I}_{N}$, where $\otimes$
denotes the Kronecker product and $\mathbf{I}_{N}$ the $N\times N$ identity
matrix. One can rewrite (5) in the following convenient form:

$-\sum_{p,q}w_{\mathcal{L}}({\mathbf{x}}_{p},{\mathbf{x}}_{q}){\mathbf{s}}_{p}^{t}{\mathbf{s}}_{q}={\mathbf{s}}^{t}\Psi{\mathbf{s}}$
(35)

Which allows us to write the difference (33) as:

$\displaystyle\mathcal{R}^{i}({\mathbf{S}})-{\cal
A}^{i,n}({\mathbf{S}})=\lambda({\mathbf{s}}^{t}\Psi{\mathbf{s}}-{\mathbf{s}}^{t}\Psi{\mathbf{s}}^{i,n})$
(36)

Now, notice that the Kronecker product $\Psi$ is negative semi-definite when
${\mathbf{W}}_{\mathcal{L}}$ is positive semi-definite. Hence, the function
${\mathbf{s}}\rightarrow{\mathbf{s}}^{t}\Psi{\mathbf{s}}$ is concave and,
therefore, is upper bounded by its first-order approximation at current
solution ${\mathbf{S}}^{i,n}$ (outer iteration $i$, inner iteration $n$). In
fact, concavity arguments are standard in deriving auxiliary functions for
bound-optimization algorithms [33]. With this condition, we have the following
auxiliary function for the Laplacian-term relaxation:

${\mathbf{s}}^{t}\Psi{\mathbf{s}}\leq({\mathbf{s}}^{i,n})^{t}\Psi{\mathbf{s}}^{i,n}+(\Psi{\mathbf{s}}^{i,n})^{t}({\mathbf{s}}-{\mathbf{s}}^{i,n})$
(37)

Equation (37) directly implies that:

$\displaystyle{\mathbf{s}}^{t}\Psi{\mathbf{s}}-{\mathbf{s}}^{t}\Psi{\mathbf{s}}^{i,n}\leq
0$ (38)

Which implies that:

$\displaystyle\mathcal{R}^{i}({\mathbf{S}})\leq{\cal A}^{i,n}({\mathbf{S}})$
(39)

Finally note that for ${\mathbf{s}}={\mathbf{s}}^{i,n}$, equation (37) becomes
an equality, and so do following equations (38) and (39), which finalizes the
proof. ∎

##  Appendix B Proof of Proposition 2

Let us consider outer iteration $i$, with the associated current assignements
${\mathbf{S}}_{i}$. For each class $k$, let us prove that
$\\{{\mathbf{m}}_{k}^{i,n}\\}_{n\in\mathbb{N}}$ is a Cauchy sequence. Recall
the recursive relation:

$\displaystyle{\mathbf{m}}_{k}^{i,n}$
$\displaystyle=\frac{\sum_{p}s_{p,k}w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}}_{k}^{i,n}){\mathbf{x}}_{p}}{\sum_{p}s_{p,k}w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}}_{k}^{i,n})}$
(40)

with
$w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}}_{k}^{i,n})=\exp(-\lVert\frac{{\mathbf{x}}_{p}-{\mathbf{m}}_{k}^{i,n}}{\sigma}\rVert^{2})$,
for some $\sigma>0$. To simplify notations in what follows, we define:

$\displaystyle h(x)$ $\displaystyle=\exp(-x)$ (42) $\displaystyle u^{n}$
$\displaystyle=\sum_{p}s_{p,k}w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}}_{k}^{i,n})$
(43) $\displaystyle v^{n}$
$\displaystyle=\sum_{p}s_{p,k}w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}}_{k}^{i,n}){\mathbf{x}}_{p}$
(44)

Step 1:  First, let us prove that $\\{u^{n}\\}_{n\in\mathbb{N}}$ is a Cauchy
sequence. Recall that in a metric space, a convergent sequence is necessarily
a Cauchy sequence. Therefore, we only need to show that $u^{n}$ is convergent
(i.e bounded and strictly monotonic).
Notice that for $x>0$, $0\leq h(x)\leq 1$. Therefore:

$\displaystyle u^{n}$
$\displaystyle=\sum_{p}s_{p,k}h(\lVert\frac{{\mathbf{x}}_{p}-{\mathbf{m}}_{k}^{i,n}}{\sigma}\rVert^{2})$
(45) $\displaystyle\leq\sum_{p}s_{p,k}\leq N$ (46)

Therefore, $u^{n}$ is bounded between 0 and N. Now, let us study the
consecutive differences $\Delta^{n}=u^{n+1}-u^{n}$:

$\displaystyle\Delta^{n}=\sum_{p}s_{p,k}\left[h(\frac{\lVert{\mathbf{x}}_{p}-{\mathbf{m}}_{k}^{i,n+1}\rVert^{2}}{\sigma^{2}})-h(\frac{\lVert{\mathbf{x}}_{p}-{\mathbf{m}}_{k}^{i,n}\rVert^{2}}{\sigma^{2}})\right]$
(47)

Because $h$ is convex, one can say that $\forall a,b\in\mathbb{R}$:

$\displaystyle h(a)-h(b)\geq h^{\prime}(b)(a-b)$ (48)

And because $h^{\prime}(b)=-h(b)$ in our case, one ends up with:

$\displaystyle h(a)-h(b)\geq h(b)(b-a)$ (49)

Applied with
$a=\frac{\lVert{\mathbf{x}}_{p}-{\mathbf{m}}_{k}^{i,n+1}\rVert^{2}}{\sigma^{2}}$
and
$b=\frac{\lVert{\mathbf{x}}_{p}-{\mathbf{m}}_{k}^{i,n}\rVert^{2}}{\sigma^{2}}$,
one can obtain:

$\begin{split}\Delta^{n}\geq&\sum_{p}s_{p,k}w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}}_{k}^{i,n})\left[\frac{\lVert{\mathbf{x}}_{p}-{\mathbf{m}}_{k}^{i,n}\rVert^{2}}{\sigma^{2}}-\frac{\lVert{\mathbf{x}}_{p}-{\mathbf{m}}_{k}^{i,n+1}\rVert^{2}}{\sigma^{2}}\right]\\\
=&\frac{1}{\sigma^{2}}\sum_{p}s_{p,k}w_{\mathcal{F}}({\mathbf{x}}_{p},{\mathbf{m}}_{k}^{i,n})\left[\lVert{\mathbf{m}}_{k}^{i,n}\rVert^{2}-\lVert{\mathbf{m}}_{k}^{i,n+1}\rVert^{2}\right.\\\
&\left.-2<{\mathbf{m}}_{k}^{i,n},{\mathbf{x}}_{p}>+2<{\mathbf{m}}_{k}^{i,n+1},{\mathbf{x}}_{p}>\right]\\\
\end{split}$

Now is time to recall recursive relation
${\mathbf{m}}_{k}^{i,n+1}=\displaystyle\frac{v^{n}}{u^{n}}$. By simply
expanding, one can end up with:

$\displaystyle\Delta^{n}\geq$
$\displaystyle\frac{1}{\sigma^{2}}\left[\lVert{\mathbf{m}}_{k}^{i,n}\rVert^{2}u^{n}-\frac{\lVert
v^{n}\rVert^{2}}{u^{n}}-2<{\mathbf{m}}_{k}^{i,n},v^{n}>+2\frac{\lVert
v^{n}\rVert^{2}}{u^{n}}\right]$ $\displaystyle=$
$\displaystyle\frac{1}{\sigma^{2}}u^{n}\left[\lVert{\mathbf{m}}_{k}^{i,n}\rVert^{2}-2<{\mathbf{m}}_{k}^{i,n},{\mathbf{m}}_{k}^{i,n+1}>+\lVert{\mathbf{m}}_{k}^{i,n+1}\rVert^{2}\right]$
$\displaystyle=$
$\displaystyle\frac{1}{\sigma^{2}}u^{n}\lVert{\mathbf{m}}_{k}^{i,n}-{\mathbf{m}}_{k}^{i,n+1}\rVert^{2}$
(50)

Therefore, $\Delta^{n}>0$, which shows that $\\{u^{n}\\}_{n\in\mathbb{N}}$ is
strictly increasing. This concludes the proof that $u^{n}$ is a convergent
sequence, and therefore a Cauchy one.

Step 2:  Now, on top of concluding the proof that
$\\{u^{n}\\}_{n\in\mathbb{N}}$ is a Cauchy sequence, Eq. (50) also offers an
interesting relation between $\\{\Delta^{n}\\}_{n\in\mathbb{N}}$ and the
sequence of interest $\\{{\mathbf{m}}_{k}^{i,n}\\}_{n\in\mathbb{N}}$, which we
can use. Indeed, for any $n_{0},m\in\mathbb{N}$, we can sum Eq. (50):

$\displaystyle\sum_{n=n_{0}}^{n_{0}+m}\Delta^{n}\geq$
$\displaystyle\frac{1}{\sigma^{2}}\sum_{n=n_{0}}^{n_{0}+m}u^{n}\lVert{\mathbf{m}}_{k}^{i,n}-{\mathbf{m}}_{k}^{i,n+1}\rVert^{2}$
(51) $\displaystyle\geq$
$\displaystyle\frac{u_{0}}{\sigma^{2}}\sum_{n=n_{0}}^{n_{0}+m}\lVert{\mathbf{m}}_{k}^{i,n}-{\mathbf{m}}_{k}^{i,n+1}\rVert^{2}$
(52) $\displaystyle\geq$
$\displaystyle\frac{u_{0}}{\sigma^{2}}\lVert{\mathbf{m}}_{k}^{i,n_{0}+m}-{\mathbf{m}}_{k}^{i,n_{0}}\rVert^{2}$
(53)

Where Eq. (52) follows because $\\{u^{n}\\}_{n\in\mathbb{N}}$ is strictly
increasing, and Eq. (53) follows from the triangle inequality. Now, the left-
hand side of Eq. () can be reduced to
$\sum_{n=n_{0}}^{n_{0}+m}\Delta^{n}=u^{n_{0}+M+1}-u^{n_{0}}$. But because we
proved in Step 1 that $\\{u^{n}\\}_{n\in\mathbb{N}}$ was a Cauchy sequence,
this difference is bounded by a constant. This concludes the proof that
$\\{{\mathbf{m}}_{k}^{i,n}\\}_{n\in\mathbb{N}}$ is itself a Cauchy sequence in
the Euclidean space.

Step 3:  We just proved that $\\{{\mathbf{m}}_{k}^{i,n}\\}_{n\in\mathbb{N}}$
was a Cauchy sequence. Therefore
$\\{{\mathbf{m}}_{k}^{i,n}\\}_{n\in\mathbb{N}}$ can only converge to a single
value ${\mathbf{m}}_{k}^{i,*}$. We now use the continuity of function
$g_{k}^{i}$ to conclude that ${\mathbf{m}}_{k}^{i,*}$ has to be a solution of
the initial equation (15):

$\displaystyle{\mathbf{m}}_{k}^{i,*}$
$\displaystyle=\lim_{n\rightarrow\infty}{\mathbf{m}}_{k}^{i,n+1}=\lim_{n\rightarrow\infty}g_{i}^{k}({\mathbf{m}}_{k}^{i,n})$
(54)
$\displaystyle=g_{i}^{k}(\lim_{n\rightarrow\infty}{\mathbf{m}}_{k}^{i,n})=g_{i}^{k}({\mathbf{m}}_{k}^{i,*})$
(55)