建立问题

我们知道Pytorch工具可以很方便的对深度神经网络进行训练和建模，但是你是否想过，把Pytorch当成一种解决凸优化问题的工具？今天就带大家简单尝试用Pytorch来解决一个聚类的凸优化问题，希望给大家带来启发。

开始之前，首先给出一个聚类问题的凸优化表达式。对于给定的数据集$\mathbf{X}\in {\mathbb{R}}^{N\times D}$，其中$N$表示数据规模，$D$表示数据维度。我们给定一些聚类中心为$\mathbf{C}\in {\mathbb{R}}^{Z\times D}$，其中$Z$表示聚类中心的个数。最后给定相似度矩阵$\mathbf{S}\in {\mathbb{R}}^{N\times Z}$，举例而言，$s_{ij}$就表示对于第$i$个数据点到第$j$个聚类中心的连接概率，这个数值越接近于1，则表明这个数据点越大概率属于这个聚类中心，否则就说明这个数据点不属于这个聚类中心。

非常显然的，我们给出这个聚类问题的凸优化表达形式：
$$\begin{gathered} \underset{\mathbf{C},\mathbf{S}}{\min }\sum _{i=1}^N\sum _{j=1}^Z{s}_{ij}{\left|\left|{\mathbf{x}}_i-{\mathbf{c}}_j\right|\right|}_2^2 \\ s.t.{\mathbf{s}}_i^{T}1=1,s_{ij}\ge 0 \end{gathered}$$

但是这个问题的解是非常“硬”（Hard）的，也就是每个样本点都只会去选择离它最近的那个样本中心点，换句话说就是得到的结果${\mathbf{s}}_i$和One-Hot的差不多 （sij∈{0,1}s_{ij} \in \{0,1\}sij​∈{0,1}）。简单论证一下为什么会出现这样的情况，不妨对第${\mathbf{x}}_i$个样本而言其最近的聚类中心是${\mathbf{c}}_m$，其次近的聚类中心是${\mathbf{c}}_n$。可以证明：
$${\left|\left|{\mathbf{x}}_i-{\mathbf{c}}_m\right|\right|}_2^2\le \alpha {\left|\left|{\mathbf{x}}_i-{\mathbf{c}}_m\right|\right|}_2^2+(1-\alpha ){\left|\left|{\mathbf{x}}_i-{\mathbf{c}}_n\right|\right|}_2^2$$
其中$0\le \alpha \le 1$。为了避免这种情况导致的解的塌陷，一般在其后面加一个正则约束项即可：
$$\begin{gathered} \underset{\mathbf{C},\mathbf{S}}{\min }\sum _{i=1}^N\sum _{j=1}^Z{s}_{ij}{\left|\left|{\mathbf{x}}_i-{\mathbf{c}}_j\right|\right|}_2^2+\gamma {\left|\left|\mathbf{S}\right|\right|}_2^2 \\ s.t.{\mathbf{s}}_i^{T}1=1,s_{ij}\ge 0 \end{gathered}$$

该优化问题是有闭式解的，但是，今天我们用Pytorch来强行解决它。

代码

下面是喜闻乐见的代码时间

'''
batch_size := 数据样本个数
cluster := 聚类中心个数
feat_dim := 特征的维度
'''
s = torch.rand([batch_size,cluster],requires_grad=True)
c = torch.rand([cluster,feat_dim],requires_grad=True)

上面给出了两个变量的初始化定义。这两个变量在后续过程中需要求取梯度并更新。

for epoch in range(1,50000):'''data: [Batch * Dim]dis_mat: [Batch * Cluster]out: [Batch * Cluster]'''data = Variable(torch_pack(moon_data))dis_mat = get_euclid(data,c)out = torch.nn.functional.softmax(s,dim=1)loss_means = torch.sum(out * dis_mat)loss_reg = torch.norm(out,p=2)loss = loss_means + gamma * loss_reg if epoch % 1000 == 0:print(epoch, loss)loss.backward()s.data.sub_(lr * s.grad.data)c.data.sub_(lr * c.grad.data)s.grad.data.zero_()c.grad.data.zero_()

这里我采用了对${\mathbf{s}}_i$的softmax归一化，其实用其他的归一化应该也是可以的，操作类似。在优化结束之后，得到的$out$就是我们想要的结果了。