NIPS-2019《Learning Feature-Sparse Principal Subspace》

2. 核心思想

这篇论文的核心思想是解决特征稀疏主子空间估计（Feature Sparse PCA, FSPCA）这一具有挑战性的问题。其核心创新点可以概括为以下几点：

1. 从优化角度提供确定性保证：与许多现有方法不同，这些方法通常需要对数据生成过程（如尖峰模型、次高斯分布）做出假设，并且只能提供高概率下的统计保证，本文完全从优化角度出发，提供确定性（deterministic）的分析。这意味着其算法的性能保证不依赖于特定的数据分布，使其更具普适性。
2. 利用协方差矩阵的低秩结构：作者敏锐地观察到，当总体协方差矩阵$A$是低秩（$rank(A)\le m$）时，一个看似NP-hard的FSPCA问题竟然可以被全局求解。这是论文最深刻的洞察之一。
3. 设计代理矩阵的迭代策略：对于一般的高秩协方差矩阵，作者提出了一个巧妙的迭代算法。其核心思想是，在每一步迭代中，构建一个精心设计的、低秩的代理协方差矩阵$P_t$。然后，利用在低秩情况下可以全局求解的算法（Algorithm 1）来求解这个代理问题，从而逐步逼近原问题的解。
4. 提供理论保证：论文不仅提出了新算法，还为这些算法提供了坚实的理论基础，包括近似比界（approximation bound）和收敛性保证（convergence guarantee），这在稀疏PCA领域是难得的。

3. 目标函数

论文的目标是求解特征稀疏主子空间，其目标函数是一个标准的主成分分析（PCA）方差最大化问题，但增加了一个关键的行稀疏性约束：

$$\underset{W\in {\mathbb{R}}^{d\times m}}{\max }Tr(W^{T}AW)\text{s.t.}{W}^{T}W=I_{m\times m},\Vert W{\Vert }_{2,0}\le k$$

其中：

• $A\in {\mathbb{R}}^{d\times d}$ 是给定的正半定协方差矩阵。
• $W\in {\mathbb{R}}^{d\times m}$ 是待求解的投影矩阵，其列向量是前$m$个主成分。
• $W^{T}W=I_{m\times m}$ 是正交性约束，确保主成分是正交的。
• $\Vert W{\Vert }_{2,0}$ 是行稀疏性范数（row-sparsity norm），定义为$W$中非零行的数量。这个约束强制$W$只有$k$个非零行，意味着只有$k$个原始特征被用于构建主子空间，实现了特征选择和降维的同步进行。

4. 目标函数的优化过程

论文提出了两种互补的优化策略，分别针对不同情况：

1. 全局最优策略 (GO - Algorithm 1)：适用于$rank(A)\le m$的情况。

   • 核心洞察：当协方差矩阵$A$的秩不超过主子空间维度$m$时，问题变得可解。
   • 优化步骤：
     • 选择$A$的对角线元素中最大的$k$个，记其索引为$I$。这一步基于一个关键观察：对于低秩$A$，最重要的特征（行）往往对应着较大的对角线元素。
     • 构建一个行选择矩阵$S$，根据索引$I$从$d$维中选出$k$维。
     • 在选出的$k$维子空间上，计算子矩阵$A_{I,I}$的前$m$个特征向量$V$。
     • 最终的解为$W=SV$，即把$m$个特征向量$V$“嵌入”回原始的$d$维空间中，其余未被选中的行置零。
   • 关键点：此过程是非迭代的，论文证明了在这种特定情况下，该算法能得到全局最优解。

2. 迭代代理更新策略 (IPU - Algorithm 2)：适用于一般的$rank(A)>m$情况。

   • 核心思想：将一个困难的高秩问题，通过迭代的方式，转化为一系列容易求解的低秩代理问题。
   • 优化步骤（迭代过程）：
     • 代理构建 (Proxy Construction)：在第$t$步，利用当前的解$W_t$，构建一个低秩代理矩阵：
       $$P_t=A{W}_t(W_t^{T}A{W}_t{)}^{†}{W}_t^{T}A$$
       这个$P_t$的设计非常巧妙，它保留了$A$在当前解$W_t$方向上的信息，且其秩不超过$m$。
     • 求解代理问题：将$P_t$作为新的协方差矩阵，调用Algorithm 1（GO）来求解代理问题，得到一个新的候选解${\tilde{W}}_{t+1}$。
     • 特征向量精炼 (Eigenvectors Refinement)：为了加速收敛，不直接使用${\tilde{W}}_{t+1}$，而是固定其非零行的索引（即行选择矩阵$S_{t+1}$），然后在原始协方差矩阵$A$上，重新计算这些$k$个特征所对应的$m$个最优特征向量：
       $$V_{t+1}\leftarrow \arg \underset{V^{T}V=I_{m\times m}}{\max }Tr(V^T{S}_{t+1}^{T}A{S}_{t+1}V)$$
     • 更新：将精炼后的特征向量$V_{t+1}$与行选择矩阵$S_{t+1}$结合，得到下一步的解$W_{t+1}=S_{t+1}{V}_{t+1}$。
   • 关键点：该过程是迭代的，并且论文证明了目标函数值$Tr(W^{T}AW)$在每一步都是单调递增的，这保证了算法的收敛性。

5. 主要贡献点

1. 理论突破：首次证明了在$rank(A)\le m$的条件下，FSPCA问题可以被全局求解（Algorithm 1），这是对NP-hard问题的一个重要特例分析。
2. 创新算法：提出了一个全新的迭代代理更新（IPU）框架（Algorithm 2）来解决一般的FSPCA问题，其核心是通过构建低秩代理矩阵将复杂问题简化。
3. 严格的理论保证：
   • 为Algorithm 1提供了近似比界（Theorem 5.1），当$A$的特征值衰减很快时（如指数分布、Zipf分布），该近似界非常紧，意味着解的质量很高。
   • 为Algorithm 2证明了单调递增性（Theorem 5.5）和收敛性（Theorem 5.8），这是现有许多算法所缺乏的。
4. 高效的计算复杂度：两种算法的计算复杂度都较低。Algorithm 1为$O(d+k^3)$，Algorithm 2为O(max⁡{dkm,k3}T)O(\max\{dkm, k^3\}T)O(max{dkm,k3}T)，其中$T$是迭代次数，且实验表明$T$通常很小（$\le 10$），优于或媲美现有方法。
5. 优越的实验表现：在合成和真实数据集上的实验结果表明，新算法在解释方差（Explained Variance）和命中频率（Hit Frequency）等指标上均优于或匹配现有的最先进方法（如SOAP, SRT）。

6. 实验结果

论文的实验设计严谨，结果有力地支持了其理论贡献：

1. 合成数据（Table 2）：

   • 在低秩协方差（Scheme C）上，Algorithm 1（Go）和Algorithm 2（IPU）均达到了100%的命中频率（HF=1.00），验证了其全局最优性。
   • 在多种设置下，新算法（尤其是IPU）在交集比（Intersection Ratio, IR）和相对误差（Relative Error, RE）上普遍优于或匹配SOAP和SRT。
   • 基于Algorithm 1（Low Rank Approx.）的初始化效果优于凸松弛（Convex Relaxation），且计算成本更低。

2. 真实数据（Figure 2）：

   • 在Lymphoma、NUS-WIDE和Numerical Numbers三个真实数据集上，新算法（Go和IPU）的归一化解释方差（Normalized Explained Variance）持续优于SOAP和SRT。
   • 虽然TPower(D)的解释方差最高，但其采用逐个计算特征向量的“消元法”（deflation），导致不同特征向量的稀疏模式不一致，因此不适用于特征稀疏主子空间估计。

3. 收敛性与效率（Figure 3 & 4）：

   • 收敛性（Figure 4）：实验验证了IPU算法的目标函数值单调递增，而SOAP和SRT则不具备此性质。
   • 计算时间（Figure 3）：新算法在大规模数据上计算效率高，扩展性良好。

7. 算法实现过程详解

下面结合伪代码和关键步骤，详细解释两个算法的实现：

Algorithm 1: Go for rank(A) ≤ m

procedure GO(A, m, k, d)
    % 步骤1: 选择对角线元素最大的k个特征
    % 计算协方差矩阵A的对角线元素
    diag_A = diag(A); 
    % 找到最大的k个元素的索引I (如果并列，优先选择索引小的)
    [~, I] = sort(diag_A, 'descend'); 
    I = I(1:k); % 取前k个
    
    % 步骤2: 构建行选择矩阵S
    % S是一个 d x k 的矩阵，根据索引I从d维中选择k维
    S = zeros(d, k);
    for j = 1:k
        S(I(j), j) = 1; % 第j列只有在第I(j)行为1
    end
    
    % 步骤3: 在k维子空间上计算特征向量
    % 提取A的k x k子矩阵
    A_II = A(I, I); 
    % 计算A_II的前m个特征向量V
    [~, V] = eigs(A_II, m, 'largestabs'); 
    
    % 步骤4: 将特征向量映射回d维空间
    W = S * V; % W的大小为 d x m，只有I中的行非零
    
    return W;
end

Algorithm 2: IPU for general A

procedure IPU(A, m, k, d, W0)
    t = 0;
    W_t = W0; % 初始化
    
    repeat
        t = t + 1;
        
        % 步骤1: 构建代理矩阵P_t
        % 计算中间矩阵: W_t^T * A * W_t (k x k)
        WtAWt = W_t' * A * W_t; 
        % 求其Moore-Penrose伪逆
        WtAWt_inv = pinv(WtAWt); 
        % 构建代理矩阵P_t = A * W_t * (W_t^T * A * W_t)^{-1} * W_t^T * A
        P_t = A * W_t * WtAWt_inv * W_t' * A; 
        
        % 步骤2: 调用Algorithm 1求解代理问题
        % 这里[GO(P_t, m, k, d)]会返回W_{t+1}，但论文伪代码第5行返回了[S,I]
        % 实际上，我们更关心的是通过GO得到的行索引I和特征向量V
        % 因此，我们分解GO的执行过程：
        [~, I] = sort(diag(P_t), 'descend'); 
        I = I(1:k);
        % 构建S
        S = zeros(d, k);
        for j = 1:k
            S(I(j), j) = 1;
        end
        % 计算P_t在I上的子矩阵的特征向量V_tilde
        P_II = P_t(I, I);
        [~, V_tilde] = eigs(P_II, m, 'largestabs');
        W_tilde = S * V_tilde; % 代理问题的解
        
        % 步骤3: 特征向量精炼
        % 固定行选择矩阵S（即索引I不变），在原始矩阵A上重新计算最优特征向量
        A_II = A(I, I); % 提取原始A在索引I上的子矩阵
        [~, V] = eigs(A_II, m, 'largestabs'); % 计算A_II的前m个特征向量
        
        % 步骤4: 更新W
        W_{t+1} = S * V; % 更新解
        
    until W_{t+1} == W_t (up to rotation) % 收敛条件：解不再变化
    return W_t;
end

总结：这篇论文通过深刻的洞察（低秩可解）和巧妙的设计（代理迭代），为特征稀疏PCA这一难题提供了兼具理论保证和实践效率的解决方案，是优化与统计结合的一个优秀范例。

代理矩阵 $ P_t $ 的构造原理分析

在论文中，作者提出了一个巧妙的代理矩阵 $ P_t $ 的构建方法，其公式如下：

$$P_t=A{W}_t(W_t^{T}A{W}_t{)}^{†}{W}_t^{T}A$$

这个代理矩阵的设计非常关键，因为它不仅保留了原始协方差矩阵 $ A $ 在当前解 $ W_t $ 方向上的信息，而且其秩不超过 $ m $。下面我们将详细分析这一构造的数学原理。

1. 代理矩阵 $ P_t $ 的设计目标

• 目标 1：低秩性
  $ P_t $ 的秩不超过 $ m $，即 $ \text{rank}(P_t) \leq m $。这是为了将原问题（可能高秩）转化为一个低秩子问题，从而利用低秩情况下的全局最优算法（Algorithm 1）求解。

• 目标 2：信息保留
  $ P_t $ 应该尽可能保留 $ A $ 在当前解 $ W_t $ 方向上的信息，以便迭代过程中逐步逼近原问题的解。

• 目标 3：单调性
  通过迭代更新 $ P_t $，确保目标函数值 $ \text{Tr}(W_t^T A W_t) $ 单调递增，从而保证算法的收敛性。

2. 代理矩阵 $ P_t $ 的构造细节

代理矩阵 $ P_t $ 的公式为：

$$P_t=A{W}_t(W_t^{T}A{W}_t{)}^{†}{W}_t^{T}A$$

其中：

• $ A $ 是原始的协方差矩阵。
• $ W_t $ 是第 $ t $ 步的当前解，满足正交性约束 $ W_t^T W_t = I_{m \times m} $。
• $ (W_t^T A W_t)^\dagger $ 是矩阵 $ W_t^T A W_t $ 的 Moore-Penrose 广义逆。

3. 代理矩阵 $ P_t $ 的秩分析

(1) 矩阵 $ W_t^T A W_t $ 的性质

• $ W_t $ 是一个 $ d \times m $ 的矩阵，且满足 $ W_t^T W_t = I_{m \times m} $，因此 $ W_t $ 的列是正交的。
• $ W_t^T A W_t $ 是一个 $ m \times m $ 的矩阵，其秩最多为 $ m $，因为 $ W_t $ 的列数为 $ m $。

(2) 广义逆 $ (W_t^T A W_t)^\dagger $ 的性质

• $ (W_t^T A W_t)^\dagger $ 是 $ W_t^T A W_t $ 的 Moore-Penrose 广义逆，它是一个 $ m \times m $ 的矩阵。
• 如果 $ W_t^T A W_t $ 是满秩的，则 $ (W_t^T A W_t)^\dagger = (W_t^T A W_t)^{-1} $。
• 如果 $ W_t^T A W_t $ 不满秩，则 $ (W_t^T A W_t)^\dagger $ 仍然保持 $ m \times m $ 的大小，但秩不会超过 $ m $。

(3) $ P_t $ 的秩

• $ P_t = A W_t (W_t^T A W_t)^\dagger W_t^T A $ 可以分解为：
  $$P_t=A{W}_{t}Q{W}_t^{T}A$$
  其中 $ Q = (W_t^T A W_t)^\dagger $。
• $ W_t $ 是 $ d \times m $ 的矩阵，因此 $ W_t Q W_t^T $ 的秩最多为 $ m $。
• $ A $ 是一个 $ d \times d $ 的矩阵，但它与 $ W_t Q W_t^T $ 相乘后，$ P_t $ 的秩仍然不会超过 $ m $，因为 $ W_t Q W_t^T $ 的秩限制了整个表达式的秩。

因此，$ P_t $ 的秩不超过 $ m $，即 $ \text{rank}(P_t) \leq m $。

4. 代理矩阵 $ P_t $ 的信息保留特性

(1) $P_t$ 的定义

$$P_t=A{W}_t(W_t^{T}A{W}_t{)}^{†}{W}_t^{T}A$$
可以理解为：

• $W_t^{T}A{W}_t$ 是 $A$ 在 $W_t$ 列空间上的投影。
• $ (W_t^T A W_t)^\dagger $ 是对 $W_t^{T}A{W}_t$ 的逆操作，用于恢复部分信息。
• 最终，$P_t$ 是通过 $A$ 和 $W_t$ 的相互作用构建的，保留了 $A$ 在 $W_t$ 方向上的信息。

(2) $ P_t $ 的投影性质

• $W_t$ 的列是正交基，因此 $W_t{W}_t^T$ 是一个投影矩阵，将任意向量投影到 $W_t$ 的列空间。
• $ P_t $ 的形式表明，它主要关注 $A$ 在 $W_t$ 列空间上的行为，同时通过 $(W_t^{T}A{W}_t{)}^{†}$ 对 $W_t$ 方向上的信息进行调整和优化。

(3) 与原矩阵 $ A $ 的关系

• $P_t$ 是 $A$ 的一个低秩近似，但它并不是简单的截断奇异值分解（SVD）结果，而是基于当前解 $W_t$ 的动态近似。
• $P_t$ 的设计确保了它在 $W_t$ 方向上的表现尽可能接近 $A$，同时保持低秩性质。

5. 代理矩阵 $ P_t $ 的单调性

(1) 目标函数

原问题的目标函数为：
$$\text{Tr}(W^{T}AW)$$
迭代过程中，我们希望目标函数值单调递增。

(2) 迭代过程中的单调性

• 在第 $t$ 步，使用 $P_t$ 作为代理矩阵，调用 Algorithm 1 求解代理问题，得到新的解 ${\tilde{W}}_{t+1}$。
• 为了加速收敛，进一步对 ${\tilde{W}}_{t+1}$ 的特征向量进行精炼，得到最终的 $W_{t+1}$。
• 论文证明了每一步迭代中，目标函数值 $\text{Tr}(W_t^{T}A{W}_t)$ 都会单调递增，即：
  $$\text{Tr}(W_t^{T}A{W}_t)<\text{Tr}(W_{t+1}^{T}A{W}_{t+1})$$
  这一性质保证了算法的收敛性。

6. 总结

代理矩阵 $P_t$ 的构造公式 $P_t=A{W}_t(W_t^{T}A{W}_t{)}^{†}{W}_t^{T}A$ 的数学原理可以总结为以下几点：

1. 低秩性：通过 $W_t$ 和 $(W_t^{T}A{W}_t{)}^{†}$ 的设计，确保 $P_t$ 的秩不超过 $m$。
2. 信息保留：$P_t$ 保留了 $A$ 在 $W_t$ 方向上的信息，同时通过广义逆 $(W_t^{T}A{W}_t{)}^{†}$ 对信息进行优化。
3. 单调性：迭代过程中，$P_t$ 的设计保证了目标函数值 $\text{Tr}(W_t^{T}A{W}_t)$ 单调递增，从而确保算法收敛。

这种设计使得代理矩阵 $P_t$ 成为连接低秩情况和高秩情况的关键桥梁，使得全局最优算法（Algorithm 1）能够应用于一般高秩协方差矩阵的情况，同时保证了算法的效率和理论保证。

代理矩阵 $P_t$ 的物理意义

在论文中，代理矩阵 $ P_t $ 的公式为：

$$P_t=A{W}_t(W_t^{T}A{W}_t{)}^{†}{W}_t^{T}A$$

这一公式是算法设计中的核心部分，其物理意义可以从多个角度进行解释。以下是对其物理意义的详细分析：

1. $P_t$ 的作用：构建低秩代理矩阵

• 目标：
  $P_t$ 是一个低秩代理矩阵，其秩不超过 $m$，用于近似原始协方差矩阵 $A$。通过将高秩的 $A$ 转化为低秩的 $P_t$，可以利用低秩情况下的全局最优算法（Algorithm 1）求解代理问题。

• 物理意义：
  $P_t$ 的设计确保了它保留了 $A$ 在当前解 $W_t$ 方向上的关键信息，同时保持低秩性质。这种低秩近似使得代理问题变得可解，从而逐步逼近原问题的解。

2. $A{W}_t$ 的作用：关注 $W_t$ 方向

• $A{W}_t$ 的含义：
  $W_t$ 是当前迭代的解，代表了当前对主子空间的估计。$A{W}_t$ 表示 $A$ 在 $W_t$ 列空间上的投影，反映了 $A$ 在当前解方向上的行为。

• 物理意义：
  $A{W}_t$ 将 $A$ 的注意力集中在 $W_t$ 方向上，确保 $P_t$ 主要关注 $A$ 在当前解方向上的表现，而不是全局的 $A$ 行为。这使得 $P_t$ 能够逐步逼近原问题的解。

3. $(W_t^{T}A{W}_t{)}^{†}$ 的作用：恢复部分信息

• $ W_t^T A W_t $ 的含义：
  $W_t^{T}A{W}_t$ 是 $A$ 在 $W_t$ 列空间上的投影，是一个 $m\times m$ 的矩阵，反映了 $A$ 在 $W_t$ 方向上的行为。

• 广义逆 $ (W_t^T A W_t)^\dagger $ 的作用：
  $(W_t^{T}A{W}_t{)}^{†}$ 是 $W_t^{T}A{W}_t$ 的 Moore-Penrose 广义逆，用于“恢复”部分信息。具体来说：

  • 如果 $W_t^{T}A{W}_t$ 是满秩的，则 $(W_t^{T}A{W}_t{)}^{†}=(W_t^{T}A{W}_t{)}^{-1}$，直接计算其逆。
  • 如果 $W_t^{T}A{W}_t$ 不满秩，则 $(W_t^{T}A{W}_t{)}^{†}$ 是对 $W_t^{T}A{W}_t$ 的一种广义逆，能够处理秩不足的情况。

• 物理意义：
  $(W_t^{T}A{W}_t{)}^{†}$ 的设计确保了 $P_t$ 在 $W_t$ 方向上尽可能接近 $A$。即使 $W_t^{T}A{W}_t$ 不满秩，通过引入 $(W_t^{T}A{W}_t{)}^{†}$，可以部分恢复丢失的信息，从而保证 $P_t$ 更加接近 $A$ 在 $W_t$ 方向上的真实表现。

4. $W_t^{T}A$ 的作用：返回到全局空间

• $W_t^{T}A$ 的含义：
  $W_t^{T}A$ 是 $A$ 在 $W_t$ 列空间上的投影矩阵的转置，用于将代理矩阵 $P_t$ 的结果映射回全局空间。

• 物理意义：
  $W_t^{T}A$ 确保了 $P_t$ 的构造不仅关注 $W_t$ 方向，还能够与全局协方差矩阵 $A$ 的其他部分相兼容。这种设计使得 $P_t$ 在保留 $A$ 关键信息的同时，能够逐步逼近原问题的解。

5. 整体物理意义：动态逼近原问题

• 动态更新：
  每次迭代后，$W_t$ 更新为 $W_{t+1}$，从而动态调整 $P_t$ 的构造。这种动态更新确保了 $ P_t$ 随着迭代逐步逼近原问题的解，同时保留了 $A$ 在当前解方向上的关键信息。

• 单调性：
  $P_t$ 的设计确保了每一步迭代的目标函数值 $\text{Tr}(W_t^{T}A{W}_t)$ 单调递增，从而保证了算法的收敛性。这种单调性表明 $P_t$ 在逐步逼近原问题的解时，不会退化或偏离正确的方向。

总结

代理矩阵 $P_t$ 的物理意义可以总结为以下几点：

1. 关注 $ W_t $ 方向：
   $P_t$ 的设计确保了它主要关注 $A$ 在当前解 $W_t$ 方向上的行为，而不是全局的 $A$ 行为。这使得 $P_t$ 能够逐步逼近原问题的解。

2. 信息恢复：
   $(W_t^{T}A{W}_t{)}^{†}$ 是对 $W_t^{T}A{W}_t$ 的逆操作，用于“恢复”部分信息，确保 $P_t$ 在 $W_t$ 方向上尽可能接近 $A$，即使 $W_t^{T}A{W}_t$ 不满秩。

3. 低秩性：
   $P_t$ 的秩不超过 $m$，使其适合使用低秩情况下的全局最优算法（Algorithm 1）求解代理问题。

4. 动态更新：
   $P_t$ 的设计随着迭代逐步逼近原问题的解，同时保留了 $A$ 在当前解方向上的关键信息。

通过这些设计，$P_t$ 成为连接低秩情况和高秩情况的关键桥梁，既保证了低秩性，又保留了 $A$ 的关键信息，从而确保了代理问题的有效性和收敛性。

$$\boxed{\text{Pt 的物理意义是通过动态构建低秩代理矩阵，逐步逼近原问题的解，同时保留 A 在当前解方向上的关键信息。}}$$