EigenGame：将主成份分析（PCA）作为一个博弈游戏

EigenGame由两部分组成，“Eigen”意为特征，也是主成份分析（PCA）方法的核心。而"Game"则意为博弈论，是一种研究理性决策者之间的冲突与合作的数学模型。与本篇博客要介绍的EigenGame相关的论文主要有发表于ICLR2021的《EIGENGAME: PCA AS A NASH EQUILIBRIUM》与arXiv.org上公开预发表的《EigenGame Unloaded When playing games is better than optimizing》。这两篇论文的作者相同。

一、PCA简介

本人与PCA相关的博客见：

• 主成分分析|PCA算法大全
• 从主成分分析到发育网络的核心算法（PCA-CCIPCA-CCILCA）

此处还想强调一点：PCA并没有把握一定能提高后续机器学习任务的效果，也没有把握能够防止过拟合问题。

PCA可以用来解决的问题【Andrew Ng曾在讲PCA时提到过】：

1）减少数据因为存储而造成的内存和硬盘的占用；
2）加速训练过程；
3）高维数据可视化。

假设原始数据集为$\mathbf{X}\in {\mathbb{R}}^{m\times n}$，我们的目标是找到最佳的投影空间$\mathbf{W}=[{\mathbf{w}}_1,{\mathbf{w}}_2,...,{\mathbf{w}}_k]\in {\mathbb{R}}^{n\times k}$，其中${\mathbf{w}}_i$是单位向量且${\mathbf{w}}_i$是单位向量且${\mathbf{w}}_i$与${\mathbf{w}}_j(i\ne j)$正交，何为最佳的$\mathbf{W}$？就是原始样本点投影到$W$上之后，使得投影后的样本点方差最大。

我们先将原始数据集$\mathbf{X}$零均值化为 U = { u 1 , u 2 , . . . , u n } \mathbf{U}=\{\mathbf{u_1}, \mathbf{u_2}, ..., \mathbf{u_n}\} U={u1​,u2​,...,un​},其中，${\mathbf{u}}_{\mathbf{i}}={\mathbf{x}}_{\mathbf{i}}-{\sum }_{j=1}^n{\mathbf{x}}_{\mathbf{j}}$。

由于投影后均值为$0$，因此数据集$\mathbf{U}$沿某一单位向量$\mathbf{w}$投影后的总方差为：
$$\frac{1}{n}\sum _{i=1}^n({\mathbf{u}}_i^T\mathbf{w}{)}^2=\frac{1}{n}\sum _{i=1}^n{\mathbf{w}}^T{\mathbf{u}}_i{\mathbf{u}}_i^T\mathbf{w}={\mathbf{w}}^T\left[\frac{1}{n}\sum _{i=1}^n{\mathbf{u}}_i{\mathbf{u}}_i^T\right]\mathbf{w}$$

其中$\frac{1}{n}{\sum }_{i=1}^n{\mathbf{u}}_i{\mathbf{u}}_i^T=\frac{1}{n}{\sum }_{i=1}^n({\mathbf{x}}_i-\frac{1}{n}{\sum }_{j=1}^n{\mathbf{x}}_j)({\mathbf{x}}_i-\frac{1}{n}{\sum }_{j=1}^n{\mathbf{x}}_j{)}^T$就是原始数据集$\mathbf{X}$的协方差矩阵（因为无偏估计的原因，一般协方差矩阵除以$n-1$，这是用$n$）。

$$\lambda ={\mathbf{w}}^T\mathbf{\Sigma }\mathbf{w}$$
其中，$\lambda =\frac{1}{n}{\sum }_{i=1}^n({\mathbf{x}}_i^T\mathbf{w}{)}^2,\mathbf{\Sigma }=\frac{1}{n}{\sum }_{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^T$。

上式两边同时左乘$\mathbf{w}$，注意到$\mathbf{w}{\mathbf{w}}^T=1$（单位向量），则有

$$\lambda \mathbf{w}=\mathbf{\Sigma }\mathbf{w}$$
所以$\mathbf{w}$是矩阵$\mathbf{\Sigma }$的特征值所对应的特征向量。

欲使投影后的总方差最大，即$\lambda$最大， 可知最佳的投影向量$w$是特征值$\lambda$最大时对应的特征向量，因此，当我们将$w$设置为与具有最大的特征值向量相等时，方差会达到最大值。这个特征向量被称为第一主成分。

二、EigenGame

论文《EigenGame Unloaded When playing games is better than optimizing》中第一个公式给出EigenGame的核心思想：通过设计一个多智能体的博弈学习系统，达到最终的纳什均衡（Nash equilibrium）；此时，每个智能体则是按特征值大小排序的特征向量。

上面公式右边第一项${\hat{v}}_i^T\Sigma {\hat{v}}_i={\lambda }_i$，也即第$i$大的特征值。第二项度量了第$i$个特征向量与前面第$j$个特征向量的方向一致程度，第$i$个向量与前面第$j$个向量越一致，第二项的值越大，垂直时为最小值0。整个公式右边，期望第$i$个特征向量在尽量与它前面特征向量垂直的同时，使得它对应的特征值最大。每个智能体被赋予这样一条规则，将公式(1)作为学习目标函数，最后就能学到数据集的top-k特征向量。

两篇论文总共给出三个相关的EigenGame算法。三个算法实现功能一样，主要是计算过程串行（Algorithm1 EIgenGame R-sequential）并行(Algorithm 2 EigenGame R)，以及Alogrithm 1 $\mu -$EigenGame R实现了将整个数据切分成小batch来训练。

注意：经过验证，按照上面算法伪代码实现与sklean库里的PCA结果不一致。以Algorithm 1 $\mu -$EigenGame R算法为例，$\text{reward}\leftarrow X_{tm}^T{X}_{tm}{\hat{v}}_i$中的$X_{tm}^T{X}_{tm}$表示一个batch大小为$n^{\prime \prime }$数据集的方差，准确公式应该为${\Sigma }_m=\frac{1}{n^{\prime \prime }}(X_{tm}-X_{mean}{)}^T(X_{tm}-X_{mean})$。然后，$\text{reward}\leftarrow \Sigma {\hat{v}}_i$,$\text{penalties}\leftarrow {\sum }_{j<i}<{\hat{v}}_i,\Sigma {\hat{v}}_j>{\hat{v}}_j$。总之一句，PCA标准过程：1）数据去中心化；2）利用中心化的数据计算协方差（要除以数据的个数）。

三、python实现注意

3.1 sklearn PCA结果

• 加载手写数字数据集

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasets

digits = datasets.load_digits()
X = digits.data
y = digits.target

• 利用sklearn.decomposition.PCA计算基准结果

from sklearn.model_selection import train_test_split
X_train,X_test,y_train,y_test = train_test_split(X,y,random_state=666)
from sklearn.decomposition import PCA

pca = PCA(n_components=10)
pca.fit(X_train)
X_train_reduction = pca.transform(X_train)
X_test_reduction = pca.transform(X_test)
pca.explained_variance_

前10个特征向量对应的特征值如下：

3.2 EigenGame结果

此处给出三个关于EigenGame的结果

3.2.1 EigenGame基本代码

本博客以最简的串行计算方式，来实现EigenGame

import numpy as np

def sample_spherical(npoints, ndim=3):
    vec = np.random.randn(ndim, npoints)
    vec /= np.linalg.norm(vec, axis=0)
    return vec

def update_eigenvector(M, parents, d):
    """
    Follow Algorithm 1
    """
    v_i = sample_spherical(1, d)  # want d x 1
    lr = 0.001
    t_i = 10000 # TEMP
    # precompute rewards upfront for each vector of parents
    rewards_j = [np.matmul(M, v_j) for v_j in parents]
    for i in range(t_i):
        reward_i = np.matmul(M, v_i)  # n x 1
        penalty = np.zeros((reward_i.shape))
        for r_j in rewards_j:
            p = float(np.dot(reward_i.T, r_j) / np.dot(r_j.T, r_j))
            penalty += p * r_j

        delta_vi = 2.0 * (reward_i - penalty)

        reimann_projection = delta_vi - float(np.dot(delta_vi.T, v_i)) * v_i

        v_prime = v_i + lr * reimann_projection
        
        v_i = v_prime / np.linalg.norm(v_prime, axis=0)

    return v_i, parents

3.2.2 论文算法中原始作法:$X^{T}X$

M = np.dot(X_train.T, X_train)
parents = []
for i in range(10):
    v_new, parents = update_eigenvector(M, parents, d=64)
    print(np.dot(np.dot(v_new.T, M), v_new))
    parents.append(v_new)

发现结果与sklearn基准结果相差很远。

3.2.3 第二种：$\frac{1}{n}{X}^{T}X$

M = np.dot(X_train.T, X_train)/len(X_train)
parents = []
for i in range(10):
    v_new, parents = update_eigenvector(M, parents, d=64)
    print(np.dot(np.dot(v_new.T, M), v_new))
    parents.append(v_new)

我们发现，第一位貌似多余了，从第二个开始，与sklearn基准结果对上了。

3.2.3 第三种：$\frac{1}{n}(X-X_{mean}{)}^T(X-X_{mean})$

X_mean = np.mean(X_train, axis=0)
X_zero = X_train - X_mean
M = np.dot(X_zero.T, X_tzero)/len(X_zero)
parents = []
for i in range(10):
    v_new, parents = update_eigenvector(M, parents, d=64)
    print(np.dot(np.dot(v_new.T, M), v_new))
    parents.append(v_new)

我们发现，结果与sklearn基准结果一致。

最后把特征向量Plot出来与sklearn的基准结果进行比较，发现除了有些向量正负号相反外，基本重合。

by windSeS 2021-5-27