Simultaneous Clustering and Optimization for Evolving Datasets

TL;DR

This paper introduces a new approach for simultaneous clustering and optimization tailored for evolving datasets, utilizing a variant of ADMM to maintain model accuracy efficiently over time.

Contribution

It proposes a novel formulation of SCO for dynamic data and develops an ADMM-based method with theoretical guarantees for specific tasks.

Findings

Effective in handling evolving datasets

Theoretical guarantees for ridge regression and convex clustering

Empirical results confirm efficiency and accuracy

Abstract

Simultaneous clustering and optimization (SCO) has recently drawn much attention due to its wide range of practical applications. Many methods have been previously proposed to solve this problem and obtain the optimal model. However, when a dataset evolves over time, those existing methods have to update the model frequently to guarantee accuracy; such updating is computationally infeasible. In this paper, we propose a new formulation of SCO to handle evolving datasets. Specifically, we propose a new variant of the alternating direction method of multipliers (ADMM) to solve this problem efficiently. The guarantee of model accuracy is analyzed theoretically for two specific tasks: ridge regression and convex clustering. Extensive empirical studies confirm the effectiveness of our method.

Tables

Table 1. TABLE I: l p subscript 𝑙 𝑝 l_{p} norm and its dual norm l q subscript 𝑙 𝑞 l_{q} .

${\parallel \cdot \parallel }_p$	${\parallel \cdot \parallel }_1$	${\parallel \cdot \parallel }_2$	${\parallel \cdot \parallel }_{\infty }$
${\parallel \cdot \parallel }_q$	${\parallel \cdot \parallel }_{\infty }$	${\parallel \cdot \parallel }_2$	${\parallel \cdot \parallel }_1$
constraints for ${\parallel \cdot \parallel }_q$	box	ball	simplex

Table 2. TABLE II: Average CPU time, in seconds, of running the proposed REG-CC method by varying constraints. The numbers in parentheses represent the respective variances.

dataset/algos	REG-CC$\left(l_1\right)$	REG-CC$\left(l_2\right)$	REG-CC$\left(l_{\infty }\right)$	parallel REG-CC$\left(l_{\infty }\right)$
iris	$5.53\left(0.05\right)$	$24.83\left(0.81\right)$	$29.86\left(4.04\right)$	$1.20\left(0.0004\right)$
moon	$4.26\left(0.02\right)$	$33.25\left(0.85\right)$	$29.50\left(0.71\right)$	$2.36\left(0.0137\right)$
segment	$54.29\left(11.94\right)$	$106.05\left(4.60\right)$	$147.16\left(42.03\right)$	$6.44\left(0.06\right)$
svm-guide	$318.58\left(24.70\right)$	$731.83\left(38.28\right)$	$1209.35\left(5.44\right)$	$20.32\left(0.04\right)$

Table 3. TABLE III: Average CPU time, in seconds, needed to execute the proposed REG-RG method by varying the constraints. The numbers in the parentheses represent variance.

dataset/algos	REG-RG$\left(l_1\right)$	REG-RG$\left(l_2\right)$	REG-RG$\left(l_{\infty }\right)$	parallel REG-RG$\left(l_{\infty }\right)$
space-ga	$26.41\left(0.145\right)$	$25.25\left(1.586\right)$	$26.63\left(0.225\right)$	$9.43\left(0.34\right)$
airfoil	$29.43\left(0.52\right)$	$29.61\left(0.004\right)$	$29.40\left(0.13\right)$	$5.68\left(0.051\right)$

Equations

$$\underset{\mathbf{X}\in {\mathbb{R}}^{n\times d}}{\min }\underset{\text{the empirical loss}}{\underbrace{\sum _{v_i\in \mathcal{V}}{f}_i({\mathbf{X}}_i)}}+\alpha \underset{\text{the sum-of-norms regularizer}}{\underbrace{\sum _{e_{ij}\in \mathcal{E}}{w}_{ij}{\Vert {\mathbf{X}}_i-{\mathbf{X}}_j\Vert }_p}}.$$

$$\sum _{v_i\in \mathcal{V}}{f}_i({\mathbf{X}}_i)=\sum _{v_i\in \mathcal{V}}{\Vert {\mathbf{X}}_i-{\mathbf{A}}_i\Vert }_2^2={\Vert \mathbf{X}-\mathbf{A}\Vert }_F^2$$

$$\sum _{v_i\in \mathcal{V}}{f}_i({\mathbf{X}}_i)=\sum _{i=1}^n\left({\Vert {\mathbf{A}}_i{\mathbf{X}}_i^{\mathrm{T}}-{\mathbf{y}}_i\Vert }_2^2+\gamma {\Vert {\mathbf{X}}_i\Vert }_2^2\right)$$

$$\underset{\mathbf{X}\in {\mathbb{R}}^{n\times d}}{\min }f(\mathbf{X};\mathbf{A})+\alpha \sum _{e_{ij}\in \mathcal{E}}{w}_{ij}{\Vert {\mathbf{X}}_i-{\mathbf{X}}_j\Vert }_p.$$

$${\mathbf{Q}}_k=(0,\underset{{\mathbf{Q}}_{ki}:=\alpha w_{ij}}{\underbrace{{\mathbf{Q}}_{ki}}},0,\underset{{\mathbf{Q}}_{kj}:=-\alpha w_{ij}}{\underbrace{{\mathbf{Q}}_{kj}}},0).$$

$$\underset{\mathbf{X}\in {\mathbb{R}}^{n\times d}}{\min }f(\mathbf{X};\mathbf{A})+{\Vert \mathbf{Q}\mathbf{X}\Vert }_{1,p}.$$

$$\underset{\mathbf{\lambda }\in {\mathbb{R}}^{m\times d}}{\min }$$

$${\Vert {\mathbf{\lambda }}_i\Vert }_q\le 1,1\le i\le m.$$

$$\nabla f({\mathbf{X}}_{\ast };\mathbf{A})+({\mathbf{I}}_d\otimes \mathbf{Q}{)}^{\mathrm{T}}\mathrm{vec}({\mathbf{\lambda }}_{\ast })=0.$$

$$\underset{\mathbf{\lambda }\in {\mathbb{R}}^{m\times d}}{\min }$$

$${\Vert {\mathbf{\lambda }}_i\Vert }_q\le 1,1\le i\le m.$$

$$\Vert \nabla f({\mathbf{X}}_{\ast };\mathbf{A})\Vert =$$

$${\Delta }_{f^{\ast }}:=$$

$$-f^{\ast }(-{\mathrm{vec}}^{\mathrm{T}}({\mathbf{\lambda }}_{\ast })({\mathbf{I}}_d\otimes \mathbf{Q});\mathbf{A})|.$$

$$\underset{\mathbf{X}\in {\mathbb{R}}^{n\times d}}{\min }{\Vert \mathbf{X}-\mathbf{A}\Vert }_F^2+{\Vert \mathbf{Q}\mathbf{X}\Vert }_{1,p}.$$

$$\underset{\mathbf{\lambda }\in {\mathbb{R}}^{m\times d}}{\min }$$

$$+\frac{1}{4}{\mathrm{vec}}^{\mathrm{T}}(\mathbf{\lambda })({\mathbf{I}}_d\otimes \mathbf{Q})({\mathbf{I}}_d\otimes \mathbf{Q}{)}^{\mathrm{T}}\mathrm{vec}(\mathbf{\lambda })$$

$${\Vert {\mathbf{\lambda }}_i\Vert }_q\le 1,1\le i\le m.$$

$$\underset{\mathbf{\lambda }\in {\mathbb{R}}^{m\times d}}{\min }$$

$$+\frac{1}{4}{\mathrm{vec}}^{\mathrm{T}}(\mathbf{\lambda })({\mathbf{I}}_d\otimes \mathbf{Q})({\mathbf{I}}_d\otimes \mathbf{Q}{)}^{\mathrm{T}}\mathrm{vec}(\mathbf{\lambda })$$

$$+\beta {\Vert ({\mathbf{I}}_d\otimes \mathbf{Q}{)}^{\mathrm{T}}\mathrm{vec}(\mathbf{\lambda })\Vert }_s$$

$${\Vert {\mathbf{\lambda }}_i\Vert }_q\le 1,1\le i\le m.$$

$$\underset{\mathbf{X}\in {\mathbb{R}}^{n\times d}}{\min }$$

$$\mathbf{\Omega }:=\text{diag}(\mathrm{vec}(\mathbf{A}))\text{diag}(\mathrm{vec}(\mathbf{A}))+\gamma {\mathbf{I}}_{nd},$$

$$\mathbf{\Lambda }:=(1_{1\times d}\otimes {\mathbf{I}}_n)\text{diag}(\mathrm{vec}(\mathbf{A})).$$

$$\underset{\mathbf{\lambda }\in {\mathbb{R}}^{m\times d}}{\min }$$

$$+\frac{1}{2}{\mathrm{vec}}^{\mathrm{T}}(\mathbf{\lambda })({\mathbf{I}}_d\otimes \mathbf{Q}){\mathbf{\Omega }}^{-1}\mathbf{\Lambda }\mathbf{y}+\beta {\Vert ({\mathbf{I}}_d\otimes \mathbf{Q}{)}^{\mathrm{T}}\mathrm{vec}(\mathbf{\lambda })\Vert }_s$$

$${\Vert {\mathbf{\lambda }}_i\Vert }_q\le 1,1\le i\le m.$$

$$\mathcal{C}=\{\mathbf{\lambda }|{\Vert {\mathbf{\lambda }}_i\Vert }_q\le 1,1\le i\le m\}.$$

$$\mathbf{u}:=({\mathbf{I}}_d\otimes \mathbf{Q}{)}^{\mathrm{T}}\mathrm{vec}(\mathbf{\lambda }).$$

$$\underset{\mathbf{\lambda }\in \mathcal{C}}{\min }$$

$$({\mathbf{I}}_d\otimes \mathbf{Q}{)}^{\mathrm{T}}\mathrm{vec}(\mathbf{\lambda })-\mathbf{u}=0.$$

Full text

1 Introduction

Simultaneous clustering and optimization (SCO) has recently drawn much attention in the machine learning and data mining community [1, 2]. Let us consider an example to explain this. Suppose that we want to predict the price of houses in New York City. The prices of houses located in the same region should be predicted by using similar prediction models. The prices of houses located in different regions should be predicted by using different prediction models. Traditional methods usually involve two separate steps. Such methods first learn prediction models for every house and then use clustering methods such as k-means clustering [3, 4, 5, 6] to determine the similarity among the obtained prediction models. However, the purpose of the SCO task is to perform prediction and identify the similarity among prediction models simultaneously; this approach usually outperforms traditional solutions.

Previous methods such as network lasso [1, 7, 2] formulate the SCO problem as that of a convex objective function with a sum-of-norms regularizer, which usually leads to a high computational cost for to following reasons:

• •

The number of optimization variables increases linearly with the number of instances and features.

• •

The optimization variables are highly nonseparable due to the sum-of-norms regularization.

• •

The objective function is extremely nonsmooth.

Due to those challenges, many optimization methods such as the alternating minimization algorithm (AMA) [8] and the alternating direction method of multipliers (ADMM) [1, 2] have been developed to reduce the computational cost. Although these existing methods obtain great efficiency improvement for a static dataset, they are unable to handle an evolving dataset directly due to exceptionally high computational complexity. The reason is that when a dataset evolves over time, the optimal model of SCO has to be frequently updated over time. Otherwise, the model accuracy cannot be guaranteed. However, it is impractical to update the optimal model frequently due to the high computational cost. Therefore, finding an effective method of performing SCO on an evolving dataset with a guarantee of model accuracy is a meaningful problem.

Let us consider an example to explain our motivation. Determining a cluster path [8, 9, 10] is one of SCO tasks. This is done by performing convex clustering over multiple rounds. As illustrated in Figure 1, the cluster paths for images $1$ and $2$ are very similar due to few changes of pixels. However, the cluster path of image $3$, namely, Figure 1(f), is significantly different from them. It is impractical to update the optimal cluster path for every image timely because it takes at least $18$ seconds to obtain a cluster path. Additionally, compared with Figure 1(d) and Figure 1(e), we observe that a cluster path for image $1$, i.e., Figure 1(d), can be used for a similar image $2$ with some slight loss of accuracy. Compared with Figure 1(d) and Figure 1(f), we observe that a cluster path should be updated if the image is changed significantly. In a general scenario, we are thus motivated by the following two nontrivial and challenging problems:

• •

How to obtain a model with a guaranteed accuracy that can be used for a dataset undergoing slight changes, and

• •

When to update the model if the dataset has evolved to become sufficiently different.

In this paper, we aim to answer the above two challenging questions. We reformulate the problem of SCO as a convex problem with cone constraints in the dual space. Compared with the formulation in the previous studies, the new formulation does not contain the sum-of-norms regularizer and is thus much easier to solve. Then, a new regularizer is proposed to allow for the optimal model to be insensitive to slightly evolving data. A new metric is proposed for deciding when to update the model for a significant change of a dataset. Additionally, we propose a new variant of the alternating direction method of multipliers (ADMM) to solve the proposed problem efficiently. Furthermore, the guarantee of the model accuracy is analyzed theoretically for convex clustering and ridge regression tasks. Extensive empirical studies show the advantages of the proposed method. In brief, our contributions are summarized as follows:

• •

A new formulation of the SCO problem is proposed to handle an evolving dataset.

• •

An efficient variant of ADMM is proposed to solve the proposed problem efficiently.

• •

The accuracy of the model is analyzed theoretically for convex clustering and ridge regression tasks.

• •

Extensive empirical studies show the superiority of the proposed method.

The paper is organized as follows. Section 2 reviews the related studies. Section 3 introduces the preliminaries. Section 4 describes a new formulation of SCO for an evolving dataset. Section 5 discusses a novel ADMM method for solving the problem. Section 6 presents a theoretical analysis. Section 7 details the experiments. Section 8 concludes the paper.

2 Related studies

In this section, we briefly review the related literature.

2.1 Network lasso

Network lasso was the first method proposed to solve the SCO problem [1]. It was subsequently implemented and published as a general optimization tool for graph analysis in [11]. Recently, network lasso has drawn much attention in various application scenarios [7, 12, 2]. It was used in [7] to predict the location of a shared bike. The study [12] examined a sufficient condition for the network topology to yield an accurate solution. Network lasso was extended in [2] to handle noisy and missing data. The cited existing studies extended network lasso for static datasets; however, such approaches are unsuitable for handling evolving datasets due to a high computational cost.

2.2 Convex clustering

The special case of SCO, i.e., convex clustering [13], has been extensively studied for several years. Compared to the traditional clustering methods, it has a convex formulation, leading to a robust clustering result. For example, the clustering result of k-means is sensitive to the seeds, and picking good seeds is challenging [3, 14]. However, due to a convex objective function, the result of convex clustering can be determined. It is not impacted by any heuristic rules used in traditional clustering methods. A formulation of convex clustering was proposed in [13] by relaxing the formulation of k-means clustering. Subsequently, [15] and [16] provided several sufficient conditions for recovering the clustering membership theoretically. Other studies, e.g., [8, 17], focus on improving the efficiency of convex clustering. Although those previous studies attained great improvement of convex clustering for static datasets, they are unsuitable for handling evolving datasets due to a high computational cost. The method proposed in the paper reduces such computational cost and makes a good tradeoff between efficiency and accuracy.

3 Preliminaries

In this section, several important notations are introduced. Then, the problem of simultaneous clustering and optimization is presented.

3.1 SCO: simultaneous clustering and optimization

SCO is an optimization framework with a formulation that is usually presented as a convex objective function with a sum-of-norms regularizer. Before discussing its formulation, we introduce the data model – a graph abstraction of a dataset.

The basic data model consists of two graphs – a data graph ${\mathcal{G}}^{\text{data}}$ and a variable graph ${\mathcal{G}}^{\text{variable}}$ – that are generated as follows.

• •

Data graph ${\mathcal{G}}^{\text{data}}$. Every instance in a dataset is represented by a vertex. Generally, the $K$-nearest neighbors ($K$-NN) method is performed on the dataset. If an instance is one of the $K$ nearest neighbors of another instance, then an edge is generated to connect them. Thus, we obtain a graph ${\mathcal{G}}^{\text{data}}$ that measures the similarity among instances.

• •

**Variable graph ${\mathcal{G}}^{\text{variable}}$. ** Variable graph ${\mathcal{G}}^{\text{variable}}$ is generated based on graph ${\mathcal{G}}^{\text{data}}$. An optimization variable, e.g., ${\bf X}_{i}$, is represented by a vertex $v_{i}$. If two instances, e.g., ${\bf A}_{i}$ and ${\bf A}_{j}$ in graph ${\mathcal{G}}^{\text{data}}$, are connected by an edge, then vertices $v_{i}$ and $v_{j}$ are connected by edge $e_{ij}$.

Let us consider an illustrative example as a further explanation. Figure 2(a) shows a data graph for a dataset consisting of $6$ instances. Figure 2(b) shows the variable graph for the data graph ${\mathcal{G}}^{\text{data}}$ in the left panel. Since SCO is formulated based on the variable graph ${\mathcal{G}}^{\text{variable}}$, we denote the variable graph ${\mathcal{G}}^{\text{variable}}$ by ${\mathcal{G}}$ by default in the following content. The vertex set of ${\mathcal{G}}$ is denoted by $\mathcal{V}$, and the edge set of ${\mathcal{G}}$ is denoted by $\mathcal{E}$. Thus, SCO [1] is formulated as

[TABLE]

Here, $p$ is picked from $1,2,\infty$. Parameter $\alpha$ is a hyperparameter that should be specified. Parameter $w_{ij}$ is the weight of the edge $e_{ij}$. When the variable graph ${\mathcal{G}}$ is obtained, $w_{ij}$ is usually set to be inversely proportional to the distance between ${\bf A}_{i}$ and ${\bf A}_{j}$ [8]. The optimal value of ${\bf X}_{i}$, denoted by ${\bf X}_{i\ast}$, represents the prediction model for the $i$-th instance, i.e., ${\bf A}_{i}$.

The complete objective function of (1) consists of two parts explained as follows.

• •

Empirical loss. SCO can be used on various data analysis tasks by instantiating the empirical loss of $f_{i}({\bf X}_{i})$. Function $f_{i}(\cdot)$ usually represents the empirical loss due to the $i$-th instance, which is determined by a specific machine learning task, such as regression, classification, and clustering. For example, if

[TABLE]

holds, it represents the empirical loss of convex clustering. If the number of vertices in $\mathcal{V}$ is $n$, and

[TABLE]

holds, it represents the empirical loss of ridge regression.

• •

Sum-of-norms regularizer. The sum-of-norms regularizer $\left\lVert{\bf X}_{i}-{\bf X}_{j}\right\rVert_{p}$ is essential for performing simultaneous clustering and optimization. As illustrated in (1), SCO learns a model, e.g., ${\bf X}_{i}$, for every instance, e.g., ${\bf A}_{i}$. The most noticeable difference between SCO and the traditional machine learning task is the sum-of-norms regularizer. The regularizer controls the similarity among those models. If $\alpha=0$, every instance has a different model from those of others. With the increase in $\alpha$, some instances tend to have similar or even identical models. If $\alpha$ is sufficiently large, all instances may have the same model. We refer to an illustrative example [2] in Figure 3 to provide further explanations about the regularizer.

3.2 Notation

Suppose that graph ${\mathcal{G}}$ consists of $n$ vertices and $m$ edges. The $i$-th vertex is represented by $v_{i}$, and its corresponding optimization variable is ${\bf X}_{i}\in\mathbb{R}^{1\times d}$. The edge connecting $v_{i}$ and $v_{j}$ is represented by $e_{ij}$. The weight of $e_{ij}$ is denoted by $w_{ij}$. ${\bf A}\in\mathbb{R}^{n\times d}$ represents the data matrix consisting of $n$ instances, and every instance is characterized by $d$ features. ${\bf X}\in\mathbb{R}^{n\times d}$ represents the matrix of optimization variables. The other important notations are shown as follows.

• •

Ordinary lowercase letters, e.g., $\alpha$ and $\beta$, represent constant scalars. Bold lowercase letters, e.g., ${\bf y}$, represent vectors. Bold capital letters, e.g., ${\bf A}\in\mathbb{R}^{n\times d}$ and ${\bf Q}\in\mathbb{R}^{m\times n}$, represent matrices.

• •

A bold capital letter with a subscript, e.g., ${\bf A}_{i}\in\mathbb{R}^{1\times d}$, represents the $i$-th row of a matrix. A bold capital letter with two subscripts, e.g., ${\bf Q}_{ij}$, represents the element located at the $i$-th row and the $j$-th column.

• •

${\rm vec}(\cdot)$ represents the column stacking vectorization of a matrix.

• •

$\text{diag}({\bf v})$ represents the diagonal matrix consisting of the elements of vector ${\bf v}$.

• •

$\otimes$ represents the Kronecker product. $\odot$ represents the Hadamard product.

• •

$\left\lVert\cdot\right\rVert$ represents a norm of a vector, and $\left\lVert\cdot\right\rVert_{\ast}$ represents its dual norm.

• •

${\bf I}_{d}$ represents the identity matrix. $\mathbf{0}$ and $\mathbf{1}$ represent constant matrices with elements of [math] and $1$, respectively.

4 Simultaneous clustering and optimization for evolving datasets

In this section, we propose a new formulation for performing SCO. Then, several examples are provided to present more details about the formulation.

4.1 Formulation

Since the dataset is evolving, the change in the data matrix ${\bf A}$ will impact the solution. Without a loss of generality, we reformulate (1) as follows:

[TABLE]

Suppose that the number of vertices in $\mathcal{V}$ is n, i.e., $n=|\mathcal{V}|$, and the number of edges in $\mathcal{E}$ is $m$, i.e., $m=|\mathcal{E}|$. We introduce an auxiliary matrix ${\bf Q}\in\mathbb{R}^{m\times n}$ to help reformulate the optimization’s objective function. Any row of ${\bf Q}$ consists of one positive element, one negative element and the remaining zero elements. Any row of ${\bf Q}$ corresponds to an edge of graph ${\mathcal{G}}$. If the $k$-th row of ${\bf Q}$, i.e., ${\bf Q}_{k}$, has a positive value at the $i$-th element and a negative value at the $j$-th element, it corresponds to edge $e_{ij}$. In this case, the positive element of ${\bf Q}_{k}$ is ${\bf Q}_{ki}:=\alpha w_{ij}$, and the negative element of ${\bf Q}_{k}$ is ${\bf Q}_{kj}:=-\alpha w_{ij}$. In other words,

[TABLE]

Additionally, the sum-of-norms regularization is usually denoted by the $l_{1,p}$ norm, i.e., $\left\lVert\cdot\right\rVert_{1,p}$ [1, 8, 16, 17]. It is defined by $\left\lVert\cdot\right\rVert_{1,p}:=\sum_{e_{ij}\in\mathcal{E}}\left\lVert\cdot\right\rVert_{p}$. Using ${\bf Q}$ to reformulate the optimization objective function of (2), we thus obtain

[TABLE]

Theorem 1**.**

The dual formulation of (2) is equivalent to

[TABLE]

subject to

[TABLE]

Once the quantity $\boldsymbol{\lambda}_{\ast}$ that minimizes (4) has been obtained, the quantity ${\bf X}_{\ast}$ that minimizes the primal problem (2) is obtained by

[TABLE]

Here, $p\in\{1,2,\infty\}$, and $\frac{1}{p}+\frac{1}{q}=1$. Table I shows the specific values of $p$ and $q$.

Therefore, our new formulation of the dual problem is

[TABLE]

subject to

[TABLE]

Here, $s$ can be $1$, $2$ or $\infty$. Parameter $\beta$ controls the accuracy of the solution for the evolving dataset. A small $\beta$ means that the solution of (6) is sensitive to a perturbation of the data matrix ${\bf A}$. A large $\beta$ causes the solution of (6) to be robust to a perturbation of the data matrix.

Intuitive idea. When ${\bf A}$ evolves to be ${\bf A}+\boldsymbol{\Delta}$, the minimizer $\boldsymbol{\lambda}_{\ast}$ is not the true minimizer $\tilde{\boldsymbol{\lambda}}_{\ast}$ of (4) based on the evolving dataset ${\bf A}+\boldsymbol{\Delta}$. However, if the difference between ${\bf A}$ and ${\bf A}+\boldsymbol{\Delta}$ is not significant, it is reasonable to use $\boldsymbol{\lambda}_{\ast}$ as an approximation of the true minimizer of (4). Note that

[TABLE]

Consider a fixed ${\bf X}_{\ast}$; then, $\left\lVert\nabla f({\bf X}_{\ast};{\bf A})\right\rVert$ changes when ${\bf A}$ evolves into ${\bf A}+\boldsymbol{\Delta}$. The smaller the change is between $\left\lVert\nabla f({\bf X}_{\ast};{\bf A})\right\rVert$ and $\left\lVert\nabla f({\bf X}_{\ast};{\bf A}+\boldsymbol{\Delta})\right\rVert$, the closer ${\bf X}_{\ast}$ and the true minimizer $\tilde{{\bf X}}_{\ast}$ based on ${\bf A}+\boldsymbol{\Delta}$. Making $\left\lVert\nabla f({\bf X}_{\ast};{\bf A})\right\rVert$ insensitive to $\boldsymbol{\Delta}$ is equivalent to making $\left\lVert({\bf I}_{d}\otimes{\bf Q})^{\rm T}{\rm vec}(\boldsymbol{\lambda}_{\ast})\right\rVert$ insensitive to $\boldsymbol{\Delta}$. The reason is that $\boldsymbol{\lambda}_{\ast}$ determines ${\bf X}_{\ast}$ according to (5). We add a regularizer $\left\lVert({\bf I}_{d}\otimes{\bf Q})^{\rm T}{\rm vec}(\boldsymbol{\lambda})\right\rVert_{s}$ into the objective function of (4), which penalizes the objective function for a large $\left\lVert({\bf I}_{d}\otimes{\bf Q})^{\rm T}{\rm vec}(\boldsymbol{\lambda})\right\rVert_{s}$. When ${\bf A}$ evolves, the regularizer guarantees that the change in $\left\lVert({\bf I}_{d}\otimes{\bf Q})^{\rm T}{\rm vec}(\boldsymbol{\lambda})\right\rVert_{s}$ is bounded within the user’s control.

We define a new metric to decide when to update the model for the evolving dataset.

[TABLE]

If this metric exceeds a threshold, then the model needs to be updated. Finally, we propose Algorithm 1 to perform SCO for evolving datasets.

4.2 Examples

We provide two examples – convex clustering [13] and ridge regression [18] – to present additional explanations of the proposed formulation.

Convex clustering. The optimization objective function of convex clustering is

[TABLE]

In this case, $f({\bf X};{\bf A})=\left\lVert{\bf X}-{\bf A}\right\rVert^{2}_{F}$, and $f^{\ast}(-{\rm vec}^{\rm T}(\boldsymbol{\lambda};{\bf A})({\bf I}_{d}\otimes{\bf Q})^{\rm T})=-{\rm vec}^{\rm T}({\bf A})({\bf I}_{d}\otimes{\bf Q})^{\rm T}{\rm vec}(\boldsymbol{\lambda})+\frac{1}{4}{\rm vec}^{\rm T}(\boldsymbol{\lambda})({\bf I}_{d}\otimes{\bf Q})({\bf I}_{d}\otimes{\bf Q})^{\rm T}{\rm vec}(\boldsymbol{\lambda})$. We thus obtain the formulation of the dual problem of convex clustering as

[TABLE]

subject to

[TABLE]

The dual problem of convex clustering is

[TABLE]

subject to

[TABLE]

Once the minimizer $\boldsymbol{\lambda}_{\ast}$ of (6) has been obtained, the minimizer of the primal problem (2) is obtained according to (5), i.e., $2\left({\rm vec}({\bf X})-{\rm vec}({\bf A})\right)+({\bf I}_{d}\otimes{\bf Q})^{\rm T}{\rm vec}(\boldsymbol{\lambda}_{\ast})=0$ in this case.

Ridge regression. The objective function of ridge regression with the sum-of-norms regularizer is

[TABLE]

After the constant item is discarded, the above becomes equivalent to

[TABLE]

where $\boldsymbol{\Omega}\in\mathbb{R}^{(nd)\times(nd)}$ is defined by

[TABLE]

and $\mathbf{\Lambda}\in\mathbb{R}^{n\times(nd)}$ is defined by

[TABLE]

Therefore, the final dual problem of ridge regression is formulated as

[TABLE]

subject to

[TABLE]

Once the minimizer $\boldsymbol{\lambda}_{\ast}$ of (6) has been obtained, the minimizer of the primal problem (2) is obtained according to (5), i.e., $2\boldsymbol{\Omega}{\rm vec}({\bf X}_{\ast})-2\mathbf{\Lambda}^{\rm T}{\bf y}+({\bf I}_{d}\otimes{\bf Q})^{\rm T}{\rm vec}(\boldsymbol{\lambda}_{\ast})=0$ in this case.

5 New variant of ADMM for solving (6)

In this section, we describe a novel ADMM method for efficiently solving the proposed formulation, i.e., (6).

5.1 Algorithmic framework

The constraint set $\mathcal{C}$ of (6) is defined by

[TABLE]

Similarly, a new notation ${\bf u}\in\mathbb{R}^{md\times 1}$ is defined by

[TABLE]

Thus, (6) is reformulated as

[TABLE]

subject to

[TABLE]

Denote $h(\boldsymbol{\lambda}):=f^{\ast}(-({\bf I}_{d}\otimes{\bf Q})^{\rm T}{\rm vec}(\boldsymbol{\lambda});{\bf A})$. Recall that $g({\bf u})=\left\lVert{\bf u}\right\rVert_{1,p}$. Its augmented Lagrangian multiplier is

[TABLE]

Update $\boldsymbol{\lambda}$. The $(t+1)$-th update of $\boldsymbol{\lambda}$ is

[TABLE]

It is worth noting that most of the previous studies [8, 1, 13, 17, 16] investigate the case of $p=q=2$. However, we observe that the case of $p=1$ and $q=\infty$ is more efficient than that of $p=q=2$. If $p=1$ and $q=\infty$, the update of $\boldsymbol{\lambda}$ can be naturally parallelized. We present the details in the appendix.

Remark 1**.**

The parallelization speed-up of the update of $\boldsymbol{\lambda}$ is up to $O(d)$ if $p=1$ and $q=\infty$.

Update ${\bf u}$. The $(t+1)$-th update of ${\bf u}$ is formulated as follows:

Theorem 2**.**

[TABLE]

where

[TABLE]

According to Theorem 2, the update of ${\bf u}$ has a closed form that can be computed very efficiently.

Update $\boldsymbol{\mu}$. The $(t+1)$-th update of $\boldsymbol{\mu}$ is simply

[TABLE]

The new variant of ADMM is presented in Algorithm 2.

5.2 Complexity analysis

If $f(\cdot)$ is convex, ADMM has been proven to be convergent [19]. A new norm $\left\lVert{\bf w}\right\rVert_{\bf H}$ is defined by

[TABLE]

where ${\bf w}\in\mathbb{R}^{(nd+2md)}$ is defined by

[TABLE]

and ${\bf H}$ is defined by

[TABLE]

Therefore, we obtain the following theorem [20]:

Theorem 3** (Theorem $5.1$ in [20]).**

Assume that $f(\cdot)$ is convex. If the ADMM is run for $T$ iterations, we have

[TABLE]

where ${\bf w}_{\ast}=\left({\rm vec}^{\rm T}(\boldsymbol{\lambda}_{\ast});{\bf u}_{\ast};\boldsymbol{\mu}_{\ast}\right)$.

It means that the convergence rate of the ADMM method is $O(\frac{1}{T})$ [20]. Formally, the quality of the solution obtained by the ADMM method is presented as follows.

Remark 2**.**

If we want to obtain an $\epsilon$-accurate ($\epsilon>0$) solution, i.e., $\left\lVert{\bf w}^{T}-{\bf w}^{T+1}\right\rVert_{\bf H}^{2}\leq\epsilon$, the ADMM method needs to perform $O(\frac{1}{\epsilon})$ iterations.

6 Theoretical analysis

In this section, we analyze the accuracy of the approximate solution obtained by Algorithm 1 for two specific tasks: convex clustering and ridge regression.

6.1 Meta framework for analysis

We define the constraints of (6), i.e., $\mathcal{C}$ as

[TABLE]

Define the minimizer of (6) based on the evolving dataset ${\bf A}+\boldsymbol{\Delta}$ as

[TABLE]

Recall that $\boldsymbol{\lambda}_{\ast}$ is defined by

[TABLE]

According to Algorithm 1, we have

[TABLE]

This inequality provides a method for bounding the difference between $\boldsymbol{\lambda}_{\ast}$ and $\tilde{\boldsymbol{\lambda}}_{\ast}$. Additionally, according to (5), we have

[TABLE]

This equality transforms the bound on the difference between $\tilde{{\bf X}}_{\ast}$ and ${\bf X}_{\ast}$ into the bound on the difference between $\boldsymbol{\lambda}_{\ast}$ and $\tilde{\boldsymbol{\lambda}}_{\ast}$. Therefore, by combining (11) and (12), we are able to bound the difference between $\tilde{{\bf X}}_{\ast}$ and ${\bf X}_{\ast}$.

Remark 3**.**

The difference between $\tilde{{\bf X}}_{\ast}$ and ${\bf X}_{\ast}$ may be bounded by (11) and (12).

6.2 Guarantee of model accuracy for convex clustering

In convex clustering, $f({\bf X};{\bf A})=\left\lVert{\bf X}-{\bf A}\right\rVert_{F}^{2}$, $\nabla f({\bf X};{\bf A})=2({\rm vec}({\bf X})-{\rm vec}({\bf A}))$, and $f^{\ast}(-{\rm vec}^{\rm T}(\boldsymbol{\lambda})({\bf I}_{d}\otimes{\bf Q})^{\rm T};{\bf A})=-{\rm vec}^{\rm T}({\bf A})({\bf I}_{d}\otimes{\bf Q})^{\rm T}{\rm vec}(\boldsymbol{\lambda})+\frac{1}{4}{\rm vec}^{\rm T}(\boldsymbol{\lambda})({\bf I}_{d}\otimes{\bf Q})({\bf I}_{d}\otimes{\bf Q})^{\rm T}{\rm vec}(\boldsymbol{\lambda})$. The difference between $\tilde{{\bf X}}_{\ast}$ and ${\bf X}_{\ast}$ is bounded according to the following theorem:

Theorem 4**.**

When the data matrix ${\bf A}$ evolves to be ${\bf A}+\boldsymbol{\Delta}$, and Algorithm 1 is run to perform convex clustering, the difference between ${\bf X}_{\ast}$ obtained by Algorithm 1 based on ${\bf A}$ and the true solution $\tilde{{\bf X}}_{\ast}$ based on ${\bf A}+\boldsymbol{\Delta}$ is bounded as

[TABLE]

6.3 Guarantee of model accuracy for ridge regression

In Section 4.2, we showed that $f({\bf X};{\bf A})={\rm vec}^{\rm T}({\bf X})\boldsymbol{\Omega}{\rm vec}({\bf X})-2{\bf y}^{\rm T}\mathbf{\Lambda}{\bf X}$ and $f^{\ast}(-{\rm vec}^{\rm T}(\boldsymbol{\lambda})({\bf I}_{d}\otimes{\bf Q});{\bf A})=\frac{1}{4}{\rm vec}^{\rm T}(\boldsymbol{\lambda})({\bf I}_{d}\otimes{\bf Q})\boldsymbol{\Omega}^{-1}({\bf I}_{d}\otimes{\bf Q})^{\rm T}{\rm vec}(\boldsymbol{\lambda})+\frac{1}{2}{\rm vec}^{\rm T}(\boldsymbol{\lambda})({\bf I}_{d}\otimes{\bf Q})\boldsymbol{\Omega}^{-1}\mathbf{\Lambda}{\bf y}$ for the ridge regression task. Define two new notations, $\tilde{\boldsymbol{\Omega}}\in\mathbb{R}^{(nd)\times(nd)}$ and $\tilde{\mathbf{\Lambda}}\in\mathbb{R}^{n\times(nd)}$, as

[TABLE]

and

[TABLE]

According to (11), we have

[TABLE]

Define a new notation $\boldsymbol{\Phi}\in\mathbb{R}^{(nd)\times(nd)}$ as

[TABLE]

Theorem 5**.**

When the data matrix ${\bf A}$ evolves to be ${\bf A}+\boldsymbol{\Delta}$, and Algorithm 1 is run to perform ridge regression, the difference between ${\bf X}_{\ast}$ obtained by Algorithm 1 based on ${\bf A}$ and the true solution $\tilde{{\bf X}}_{\ast}$ based on ${\bf A}+\boldsymbol{\Delta}$ is bounded as

[TABLE]

Remark 4**.**

Given ${\bf A}$ and $\boldsymbol{\Delta}$, it is possible to obtain a relatively accurate solution for convex clustering and ridge regression by setting a small $c$ and a large $\beta$ in Algorithm 1.

7 Empirical studies

In the section, we first present the experimental settings. Then, we perform convex clustering, obtain the cluster path for video datasets, and perform the ridge regression tasks to evaluate our method.

7.1 Experimental settings

The experiments are performed on a server equipped with an Intel Xeon(R) 32-core E5-2610 CPU and 48 GB of RAM. We implement all the algorithms by using Matlab 2016b and the CVX solver [21].

We use six datasets: moon111https://cs.joensuu.fi/sipu/datasets/jain.txt, *iris222https://www.csie.ntu.edu.tw/$\sim$ cjlin/libsvmtools/datasets/multi

class.html#iris*, *segment333https://www.csie.ntu.edu.tw/$\sim$ cjlin/libsvmtools/datasets/multi

class.html#segment*, *svm-guide444https://www.csie.ntu.edu.tw/$\sim$ cjlin/libsvmtools/datasets/binar

y.html#svmguide1*, *space_ga555https://www.csie.ntu.edu.tw/$\sim$ cjlin/libsvmtools/datasets/regre

ssion.html#space_ga*, and airfoil666http://archive.ics.uci.edu/ml/datasets/Airfoil+Self-Noise. Given graph ${\mathcal{G}}$, the weight for every edge is set to be inversely proportional to the distance between the vertices; this approach is widely used in [8, 9, 22]. Additionally, our method is compared with many existing methods. These existing methods are presented briefly as follows:

• •

PRIMAL-CC [13] has been proposed to perform convex clustering. It solves the primal problem, i.e., (2), directly.

• •

SSNAL-CC [17] is a semi-smooth Newton-based augmented Lagrangian method. It has been proposed to perform convex clustering efficiently.

• •

NET-LASSO [1] has been proposed as a general framework for performing SCO for various tasks, e.g., convex clustering and ridge regression.

• •

AMA-CC [8] is an alternating minimization algorithm. It has been proposed to perform convex clustering efficiently.

We follow the above notation, whereby the minimizer of (6) for ${\bf A}$ is denoted by ${\bf X}_{\ast}$, and the true minimizer of (6) for ${\bf A}+\boldsymbol{\Delta}$ is denoted by $\tilde{{\bf X}}_{\ast}$. Their difference, i.e., $\left\lVert{\bf X}_{\ast}-\tilde{{\bf X}}_{\ast}\right\rVert_{F}^{2}$, is used to measure the approximation of ${\bf X}_{\ast}$ against $\tilde{{\bf X}}_{\ast}$. Here, $\left\lVert\cdot\right\rVert_{F}$ represents the Frobenius norm. The threshold $c$ in Algorithm 1 is set to be $c=10$ by default.

7.2 Convex clustering

As we have shown in the previous section, $f({\bf X};{\bf A})=\left\lVert{\bf X}-{\bf A}\right\rVert_{F}^{2}$ holds for convex clustering. Given a value $\alpha$, convex clustering recovers the clustering membership. If an instance contains more than $2$ features, we present the cluster path by using its first and second principal components. The edges in graph ${\mathcal{G}}$ are generated by running the $K$-NN method with $K=10$ for moon and iris and $K=5$ for segment and svm-guide. For datasets segment and svm-guide, we pick $200$ instances randomly and draw their cluster paths as illustrative examples. We increase $\alpha$ gradually and obtain a cluster path for every dataset. Additionally, we use the Gaussian distribution $N(0,0.1^{2})$ to generate the perturbed matrix $\boldsymbol{\Delta}$ that is used to simulate the evolving dataset. Parameter $\beta$ is set to $5$ by default in the evaluation. Our method is denoted by REG-CC for convex clustering. Note that the existing NET-LASSO method is unable to handle the datasets segment and svm-guide due to running out of memory.

As illustrated in Figure 4, every color represents a true cluster from the ground truth. If $\alpha$ is small, local communities of instances are detected. As $\alpha$ increases, different local communities are fused into a large community. The number of clusters declines for a large $\alpha$. Furthermore, as illustrated in Figure 5, our REG-CC method yields more accurate clustering results than do the existing methods for various values of $\alpha$. Due to the regularized item in the dual problem of (6), our REG-CC method penalizes a large fluctuation in ${\bf X}_{\ast}$. Thus, with the control of accuracy, the regularized item makes ${\bf X}_{\ast}$ robust to a perturbation of the data matrix. In other words, ${\bf X}_{\ast}$ is a satisfactory approximation of $\hat{{\bf X}}_{\ast}$ for the evolving dataset ${\bf A}+\boldsymbol{\Delta}$. Additionally, when we vary $\beta$, Figure 6 shows that a large $\beta$ yields more accurate clustering results than does a small $\beta$. The reason is that a large $\beta$ imposes a greater penalty on the fluctuation in ${\bf X}_{\ast}$, which makes ${\bf X}_{\ast}$ insensitive to the change in ${\bf A}$.

Additionally, we evaluate the performance of our method by varying $s$ in the regularizer. We test it in three kinds of evolving environments. In the first evolving environment, the perturbed matrix $\boldsymbol{\Delta}$ is sparse, and its nonzero values are generated from a Gaussian distribution, i.e., $N(0,0.1^{2})$. This leads to a sparse change in the data matrix. In this case, the number of nonzero elements in the perturbed matrix is set to be $0.2n$, where $n$ is the number of instances in the dataset. In the second evolving environment, the perturbed matrix $\boldsymbol{\Delta}$ is dense, and its values are generated from the Gaussian distribution $N(0,0.01^{2})$. This leads to a dense and relatively insignificant change in the data matrix. The third evolving environment generates the mixed perturbed matrix $\boldsymbol{\Delta}$ with $0.2n$ elements generated from the first case and the other elements generated from the second case. As illustrated in Figure 7, our method performs best in the setting of $s=1$ in all of these cases. In this setting, the regularizer favors obtaining a sparse $({\bf I}_{d}\otimes{\bf Q})^{\rm T}{\rm vec}(\boldsymbol{\lambda}_{\ast})$ in the evolving environment. Although the data matrix changes in the evolving environment, the sparsity makes the minimizer $\boldsymbol{\lambda}_{\ast}$ more robust to the change in the data matrix.

Finally, we test the efficiency of our method by varying $q$. Various values of $q$ imply different types of constraints in the dual problem (6). Variable $q$ is set to $1$, $2$ and $\infty$. In particular, if $q=\infty$, our method supports parallel computing naturally. It is implemented by using multiple cores. According to Figure 7, performance at $s=1$ is usually better than at $s=2$ and at $\infty$. Thus, we set $s=1$ in this experiment. For all settings, our algorithm is executed three times. We show the mean and the variance (shown as the number in parentheses) of the used CPU time in seconds in Table II. As shown in Table II, our method is most efficient at $q=\infty$, when it is executed in parallel.

7.3 Cluster paths for video datasets

We obtain cluster paths for the DAVIS 2017 video dataset777https://data.vision.ee.ethz.ch/csergi/share/davis/DAVIS-2017-test-challenge-480p.zip. The videos in the dataset are at 480 p resolution. In the experiment, each raw image is transformed into the corresponding grayscale image, and we determine the cluster path by using the grayscale image. Specifically, we partition every image into $256$ blocks according to the size of the image and assign every block a unique ID number. The number near a point in the cluster path represents the corresponding block in the grayscale image. Together with the position representing location in the image, every block is represented by a $3\times 1$ vector. Additionally, we set the number of neighbors to $K=4$ when $K$-NN is used to construct the sparse graph ${\mathcal{G}}$ for every image. Other parameters such as $q$ and $s$ are set to $q=1$ and $s=\infty$ by default.

Figure 8 shows the cluster paths. We obtain several interesting observations from cluster paths. At a low level, the cluster path shows the clustering membership among blocks of the grayscale image. For example, the blocks with id of $1$, $17$, and $33$ (blocks in the top-left region) belong to a cluster due to their similar gray values. Blocks $112$, $109$, and $127$ (those in the center-right region) belong to a cluster. This captures the similarity of the blocks at a low level. At a high level, we compare the four images and their corresponding cluster paths. The event of the plane crossing is detected by analyzing the change of the structures of those cluster paths. In particular, the lower-right part of cluster paths changes from sparse paths to dense paths. The change in the cluster path reflects the event of the plane crossing. The above low-level and high-level observations provide significant insights for understanding the video.

Additionally, as we have shown in the previous section, if $f({\bf X};{\bf A})=\left\lVert{\bf X}-{\bf A}\right\rVert_{F}^{2}$ holds, our method can be used to perform object segmentation for a video. Similarly, we partition every image into $256$ blocks. Every block is represented by its average gray value. The illustrative examples are presented in Figure 9. As we observe, it is significant to detect the moving object in the video. Although clustering is a traditional method of performing image segmentation, convex clustering is more robust than the classic clustering method such as k-means clustering. The clustering result of convex clustering is determined and is not impacted by the seeds or heuristic rules used in k-means clustering.

7.4 Ridge regression

We further evaluate the performance of the proposed method by conducting the ridge regression task on the space-ga and airfoil datasets. As we have shown, $f({\bf X};{\bf A})$ in (2) is instantiated to be

[TABLE]

for the ridge regression task. Here, $\gamma=5$ in the experiment, and we construct the network ${\mathcal{G}}$ by running the $K$-NN method on the chosen datasets with $K=5$. Additionally, in this analysis, the proposed method is denoted by REG-RG for the ridge regression task. The state-of-the-art method is network lasso [1], denoted by NET-LASSO. The proposed method is also compared with the method of solving the primal problem (2) directly, which is denoted by PRIMAL-RG.

As illustrated in Figure 10, the proposed REG-RG method yields more accurate prediction models than do the existing methods when datasets evolve. The superiority becomes significant with the increase in $\beta$, as shown in Figure 11. The reason is that the proposed REG-RG method tends to yield the solution that is relatively robust to the evolving data. The large $\beta$ means that REG-RG applies a larger penalty to the solution, which encourages it to be insensitive to the evolving data.

Next, we evaluate our method on three kinds of evolving data, including the sparse evolving data, the dense evolving data, and a mixture of both. As we have shown, the sparse evolving data consist of $0.2n$ nonzero values generated from the Gaussian distribution $N(0,0.1^{2})$. The dense evolving data are generated from the Gaussian distribution $N(0,0.01^{2})$. A mixture of them is obtained by generating $0.2n$ values from the sparse case, while the other values are generated from the dense case. As shown in Figure 12, it is more effective to obtain an approximate solution for the proposed method at $s=1$ than at other values of $s$. Finally, we evaluate the efficiency of the proposed method by varying the constraints. We execute every algorithm three times and record the average and the variance (the number in the parentheses) of CPU time, in seconds. Table III shows that the efficiency varies slightly when different types of constraints are used. However, due to the use of parallel computing, the proposed method is most efficient at $q=\infty$.

8 Conclusions

We investigate SCO in an evolving environment. We first reformulate the problem into a convex problem with cone constraints in a dual space. Then, a new regularizer is proposed to obtain an approximate solution for the evolving dataset. A novel ADMM method is proposed to solve the problem efficiently. Afterward, we analyze the quality of the solution theoretically for the cases of the convex clustering and ridge regression tasks. Finally, extensive empirical studies show the advantages of the proposed method.

Acknowledgments

This work was supported by the National Key R & D Program of China (Grant No. 2018YFB1003203), the National Natural Science Foundation of China (Grant Nos. 61672528, 61773392, 61701451, 61671463, and 61772544), and the National Basic Research Program (the 973 program) under Grant No. 2014CB347800.