1.演算法描述

     最近在做AutoEncoder的一些探索，看到2016年的一篇論文，雖然不是最新的，但是思路和方法值得學習。論文原文連結 http://proceedings.mlr.press/v48/xieb16.pdf，論文有感於t-SNE演算法的t-分佈，先假設初始化K個聚類中心，然後資料距離中心的距離滿足t-分佈，可以用下面的公式表示:

其中 i表示第i樣本,j表示第j個聚類中心, z表示原始特徵分佈經過Encoder之後的表徵空間。

$q_{ij}$可以解釋為樣本i屬於聚類j的概率，屬於論文上說的"軟分配"的概念。那麼“硬分配”呢?那就是樣本一旦屬於一個聚類，其餘的聚類都不屬於了，也就是其餘聚類的概率為0。由於$\alpha$在有label的訓練計劃中，是在驗證集上進行確定的，在該論文中，全部設定成了常數1。

     然後神奇的事情發生了，作者發明了一個輔助分佈也用來衡量樣本屬於某個聚類的分佈，就是下面的公式了:

，其中$f_{j}=\sum_{i}q_{ij}$

也許你會疑問，上面這個玩意怎麼來的?作者的論文中說主要考慮一下三點:

1. 強化預測。q分佈為軟分配的概率，那麼p如果使用delta分佈來表示，顯得比較原始。
2. 置信度越高，屬於某個聚類概率越大。
3. 規範每個質心的損失貢獻，以防止大類扭曲隱藏的特徵空間。分子中那個$f_{j}$就是做這個的。

假設分佈有了，原始的資料分佈也有了，剩下衡量兩個分佈近似的方法，作者使用了KL散度，公式如下：

這個也是DEC聚類的損失函式。有了具體的公式，明確一下每次迭代更新需要Update的引數：

第一個公式是優化AE中的Encoder引數，第二個公式是優化聚類中心。也就是說作者同時優化了聚類和DNN的相關引數。

作者設計的網路概念圖如下:

   DEC演算法由兩部分組成，第一部分會預訓練一個AE模型；第二部分選取AE模型中的Encoder部分，加入聚類層，使用KL散度進行訓練聚類。 

2.實驗分析

    實驗部分比較了幾種演算法，比較的指標是ACC，對比表格如下:

DEC的效果還是比較不錯的，另外值得一提的是DEC w/o backprop演算法，是將第一部分Encoder的引數固定之後，不再參加訓練，只更新聚類中心的演算法，從結果上看，並沒有兩者同時訓練效果來的好。

3.原始碼分析

論文使用的Caffe寫的，對於我這種半路出家的和尚有點吃力，網上找了一個keras的實現程式碼，https://github.com/XifengGuo/DEC-keras/blob/master/DEC.py 

首先是DEC的預訓練的部分，預訓練的模型先儲存了起來方便訓練聚類使用:

def pretrain(self, x, y=None, optimizer='adam', epochs=200, batch_size=256, save_dir='results/temp'):
        print('...Pretraining...')
        self.autoencoder.compile(optimizer=optimizer, loss='mse')
        self.autoencoder.fit(x, x, batch_size=batch_size, epochs=epochs, callbacks=cb)
        self.autoencoder.save_weights(save_dir + '/ae_weights.h5')
        self.pretrained = True

在進行訓練之前，我們看一下作者構造的一個新的網路層

class ClusteringLayer(Layer):
    .....

    def build(self, input_shape):
        assert len(input_shape) == 2
        input_dim = input_shape[1]
        self.input_spec = InputSpec(dtype=K.floatx(), shape=(None, input_dim))
        //在這裡定義了需要訓練更新的權重，就是說把K個聚類當作了“權重”來進行更新了
        self.clusters = self.add_weight((self.n_clusters, input_dim), initializer='glorot_uniform', name='clusters')
        if self.initial_weights is not None:
            self.set_weights(self.initial_weights)
            del self.initial_weights
        self.built = True

    def call(self, inputs, **kwargs):
        """ student t-distribution, as same as used in t-SNE algorithm.
                 q_ij = 1/(1+dist(x_i, u_j)^2), then normalize it.
        Arguments:
            inputs: the variable containing data, shape=(n_samples, n_features)
        Return:
            q: student's t-distribution, or soft labels for each sample. shape=(n_samples, n_clusters)
        """
        q = 1.0 / (1.0 + (K.sum(K.square(K.expand_dims(inputs, axis=1) - self.clusters), axis=2) / self.alpha))
        q **= (self.alpha + 1.0) / 2.0
        q = K.transpose(K.transpose(q) / K.sum(q, axis=1))
        return q

    def compute_output_shape(self, input_shape):
        assert input_shape and len(input_shape) == 2
        return input_shape[0], self.n_clusters

    def get_config(self):
        config = {'n_clusters': self.n_clusters}
        base_config = super(ClusteringLayer, self).get_config()
        return dict(list(base_config.items()) + list(config.items()))

接下來就是第二部分fit的過程了

def fit(self, x, y=None, maxiter=2e4, batch_size=256, tol=1e-3,
            update_interval=140, save_dir='./results/temp'):

        # Step 1: initialize cluster centers using k-means
        kmeans = KMeans(n_clusters=self.n_clusters, n_init=20)
        y_pred = kmeans.fit_predict(self.encoder.predict(x))
        y_pred_last = np.copy(y_pred)
        self.model.get_layer(name='clustering').set_weights([kmeans.cluster_centers_])

        # Step 2: deep clustering
    
        loss = 0
        index = 0
        index_array = np.arange(x.shape[0])
        for ite in range(int(maxiter)):
            if ite % update_interval == 0:
                q = self.model.predict(x, verbose=0)
                p = self.target_distribution(q)  # update the auxiliary target distribution p

                # evaluate the clustering performance
                y_pred = q.argmax(1)
                if y is not None:
                    acc = np.round(metrics.acc(y, y_pred), 5)
                    nmi = np.round(metrics.nmi(y, y_pred), 5)
                    ari = np.round(metrics.ari(y, y_pred), 5)
                    loss = np.round(loss, 5)
                    logdict = dict(iter=ite, acc=acc, nmi=nmi, ari=ari, loss=loss)
                    logwriter.writerow(logdict)
                    print('Iter %d: acc = %.5f, nmi = %.5f, ari = %.5f' % (ite, acc, nmi, ari), ' ; loss=', loss)

                # check stop criterion
                delta_label = np.sum(y_pred != y_pred_last).astype(np.float32) / y_pred.shape[0]
                y_pred_last = np.copy(y_pred)
                if ite > 0 and delta_label < tol:
                    print('delta_label ', delta_label, '< tol ', tol)
                    print('Reached tolerance threshold. Stopping training.')
                    logfile.close()
                    break

            # train on batch
            # if index == 0:
            #     np.random.shuffle(index_array)
            idx = index_array[index * batch_size: min((index+1) * batch_size, x.shape[0])]
            loss = self.model.train_on_batch(x=x[idx], y=p[idx])
            index = index + 1 if (index + 1) * batch_size <= x.shape[0] else 0

            # save intermediate model
            if ite % save_interval == 0:
                print('saving model to:', save_dir + '/DEC_model_' + str(ite) + '.h5')
                self.model.save_weights(save_dir + '/DEC_model_' + str(ite) + '.h5')

            ite += 1

        return y_pred

文章只是簡單分析了一下，具體細節還是看原始碼來得實在。

想到的一些問題如下:

1.DEC的假設分佈從實驗效果上看起來不錯，是否存在其他的比較牛逼的分佈呢?

2.DEC聚類不能產生新的樣本，這也是VADE類似的聚類演算法的優勢，抽空再看看。

3.DEC的使用除了聚類，還有什麼呢？個人能想到的一點就是做離散化，相比於AE的那種Encoder的抽象降唯來說，DEC可以產生離散的變數，而不是多維的連續變數。後續可以在工程中嘗試一下。