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）高维数据可视化。

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

我们先将原始数据集X\mathbf{X}X零均值化为U={u1,u2,...,un}\mathbf{U}=\{\mathbf{u_1}, \mathbf{u_2}, ..., \mathbf{u_n}\}U={u1,u2,...,un},其中，ui=xi−∑j=1nxj\mathbf{u_i}=\mathbf{x_i}-\sum_{j=1}^{n} \mathbf{x_j}ui=xi−∑j=1nxj。

由于投影后均值为0\mathbf{0}0，因此数据集U\mathbf{U}U沿某一单位向量w\mathbf{w}w投影后的总方差为：
1n∑i=1n(uiTw)2=1n∑i=1nwTuiuiTw=wT[1n∑i=1nuiuiT]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} n1i=1∑n(uiTw)2=n1i=1∑nwTuiuiTw=wT[n1i=1∑nuiuiT]w

其中1n∑i=1nuiuiT=1n∑i=1n(xi−1n∑j=1nxj)(xi−1n∑j=1nxj)T\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)^Tn1∑i=1nuiuiT=n1∑i=1n(xi−n1∑j=1nxj)(xi−n1∑j=1nxj)T就是原始数据集X\mathbf{X}X的协方差矩阵（因为无偏估计的原因，一般协方差矩阵除以n−1n-1n−1，这是用nnn）。

λ=wTΣw\lambda=\mathbf{w}^T\mathbf{\Sigma} \mathbf{w} λ=wTΣw
其中，λ=1n∑i=1n(xiTw)2,Σ=1n∑i=1nxixiT\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}λ=n1∑i=1n(xiTw)2,Σ=n1∑i=1nxixiT。

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

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

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

二、EigenGame

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

上面公式右边第一项v^iTΣv^i=λi\hat{v}_i^T\Sigma \hat{v}_i=\lambda_iv^iTΣv^i=λi，也即第iii大的特征值。第二项度量了第iii个特征向量与前面第jjj个特征向量的方向一致程度，第iii个向量与前面第jjj个向量越一致，第二项的值越大，垂直时为最小值0。整个公式右边，期望第iii个特征向量在尽量与它前面特征向量垂直的同时，使得它对应的特征值最大。每个智能体被赋予这样一条规则，将公式(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算法为例，reward←XtmTXtmv^i\text{reward} \leftarrow X_{tm}^{T} X_{tm}\hat{v}_ireward←XtmTXtmv^i中的XtmTXtmX_{tm}^{T}X_{tm}XtmTXtm表示一个batch大小为n′′n''n′′数据集的方差，准确公式应该为Σm=1n′′(Xtm−Xmean)T(Xtm−Xmean)\Sigma_{m}=\frac{1}{n''}(X_{tm}-X_{mean})^{T}(X_{tm}-X_{mean})Σm=n′′1(Xtm−Xmean)T(Xtm−Xmean)。然后，reward←Σv^i\text{reward} \leftarrow \Sigma\hat{v}_ireward←Σv^i,penalties←∑j<i<v^i,Σv^j>v^j\text{penalties} \leftarrow \sum_{j<i}<\hat{v}_i, \Sigma \hat{v}_j>\hat{v}_jpenalties←∑j<i<v^i,Σv^j>v^j。总之一句，PCA标准过程：1）数据去中心化；2）利用中心化的数据计算协方差（要除以数据的个数）。

三、python实现注意

3.1 sklearn PCA结果

加载手写数字数据集

import numpy as np
import matplotlib.pyplot as plt
from sklearn import datasetsdigits = 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 PCApca = 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 npdef sample_spherical(npoints, ndim=3):vec = np.random.randn(ndim, npoints)vec /= np.linalg.norm(vec, axis=0)return vecdef update_eigenvector(M, parents, d):"""Follow Algorithm 1"""v_i = sample_spherical(1, d)  # want d x 1lr = 0.001t_i = 10000 # TEMP# precompute rewards upfront for each vector of parentsrewards_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 1penalty = 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_jdelta_vi = 2.0 * (reward_i - penalty)reimann_projection = delta_vi - float(np.dot(delta_vi.T, v_i)) * v_iv_prime = v_i + lr * reimann_projectionv_i = v_prime / np.linalg.norm(v_prime, axis=0)return v_i, parents

3.2.2 论文算法中原始作法:XTXX^TXXTX

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 第二种：1nXTX\frac{1}{n}X^TXn1XTX

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 第三种：1n(X−Xmean)T(X−Xmean)\frac{1}{n}(X-X_{mean})^T(X-X_{mean})n1(X−Xmean)T(X−Xmean)

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