1

首先让我澄清一下，这里的“稀疏 PCA”是指具有 L1 惩罚和稀疏负载的 PCA，而不是稀疏矩阵上的 PCA。

我已经阅读了 Zou 和 Hastie 的关于稀疏 PCA 的论文，我已经阅读了关于 sklearn.decomposition.SparsePCA 的文档，并且我知道如何使用 PCA，但我似乎无法从 SparsePCA 获得正确的结果。

也就是说，当 L1 惩罚为 0 时，SparsePCA 的结果应该与 PCA 一致，但负载差异很大。为了确保我没有弄乱任何超参数，我在 R 中使用了相同的超参数（收敛容差、最大迭代次数、岭惩罚、套索惩罚……）和来自“elasticnet”的“spca”，R 给了我正确的结果。如果有人有使用此功能的经验并且可以让我知道我是否犯了任何错误，我宁愿不必通过 SparsePCA 的源代码。

以下是我生成数据集的方式。这有点令人费解，因为我想要一个特定的马尔可夫决策过程来测试一些强化学习算法。只需将其视为一些非稀疏数据集。

import numpy as np
from sklearn.decomposition import PCA, SparsePCA
import numpy.random as nr

def transform(data, TranType=None):
    if TranType == 'quad':
        data = np.minimum(np.square(data), 3)
    if TranType == 'cubic':
        data = np.maximum(np.minimum(np.power(data, 3), 3), -3)
    if TranType == 'exp':
        data = np.minimum(np.exp(data), 3)
    if TranType == 'abslog':
        data = np.minimum(np.log(abs(data)), 3)
    return data

def NewStateGen(OldS, A, TranType, m=0, sd=0.5, nsd=0.1, dim=64):
    # dim needs to be a multiple of 4, and preferably a multiple of 16.
    assert (dim == len(OldS) and dim % 4 == 0)
    TrueDim = dim / 4
    NewS = np.zeros(dim)
    # Generate new state according to action
    if A == 0:
        NewS[range(0, dim, 4)] = transform(OldS[0:TrueDim], TranType) + \
            nr.normal(scale=nsd, size=TrueDim)
        NewS[range(1, dim, 4)] = transform(OldS[0:TrueDim], TranType) + \
            nr.normal(scale=nsd, size=TrueDim)
        NewS[range(2, dim, 4)] = nr.normal(m, sd, size=TrueDim)
        NewS[range(3, dim, 4)] = nr.normal(m, sd, size=TrueDim)
        R = 2 * np.sum(transform(OldS[0:int(np.ceil(dim / 32.0))], TranType)) - \
            np.sum(transform(OldS[int(np.ceil(dim / 32.0)):(dim / 16)], TranType)) + \
            nr.normal(scale=nsd)
    if A == 1:
        NewS[range(0, dim, 4)] = nr.normal(m, sd, size=TrueDim)
        NewS[range(1, dim, 4)] = nr.normal(m, sd, size=TrueDim)
        NewS[range(2, dim, 4)] = transform(OldS[0:TrueDim], TranType) + \
            nr.normal(scale=nsd, size=TrueDim)
        NewS[range(3, dim, 4)] = transform(OldS[0:TrueDim], TranType) + \
            nr.normal(scale=nsd, size=TrueDim)
        R = 2 * np.sum(transform(OldS[int(np.floor(dim / 32.0)):(dim / 16)], TranType)) - \
            np.sum(transform(OldS[0:int(np.floor(dim / 32.0))], TranType)) + \
            nr.normal(scale=nsd)

    return NewS, R

def MDPGen(dim=64, rep=1, n=30, T=100, m=0, sd=0.5, nsd=0.1, TranType=None):
    X_all = np.zeros(shape=(rep*n*T, dim))
    Y_all = np.zeros(shape=(rep*n*T, dim+1))
    A_all = np.zeros(rep*n*T)
    R_all = np.zeros(rep*n*T)

    for j in xrange(rep*n):
        # Data for a single subject
        X = np.zeros(shape=(T+1, dim))
        A = np.zeros(T)
        R = np.zeros(T)
        NewS = np.zeros(dim)

        X[0] = nr.normal(m, sd, size=dim)
        for i in xrange(T):
            OldS = X[i]
            # Pick a random action
            A[i] = nr.randint(2)
            # Generate new state according to action
            X[i+1], R[i] = NewStateGen(OldS, A[i], TranType, m, sd, nsd, dim)

        Y = np.concatenate((X[1:(T+1)], R.reshape(T, 1)), axis=1)
        X = X[0:T]

        X_all[(j*T):((j+1)*T)] = X
        Y_all[(j*T):((j+1)*T)] = Y
        A_all[(j*T):((j+1)*T)] = A
        R_all[(j*T):((j+1)*T)] = R

    return {'X': X_all, 'Y': Y_all, 'A': A_all, 'R': R_all, 'rep': rep, 'n': n, 'T': T}    

nr.seed(1)
MDP = MDPGen(dim=64, rep=1, n=30, T=90, sd=0.5, nsd=0.1, TranType=None)
X = MDP.get('X').astype(np.float32)

现在我运行 PCA 和 SparsePCA。当 lasso 惩罚 'alpha' 为 0 时，SparsePCA 应该给出与 PCA 相同的结果，但事实并非如此。其他超参数使用 R 中 elasticnet 的默认值设置。如果我使用 SparsePCA 的默认值，结果仍然不正确。

PCA_model = PCA(n_components=64)
PCA_model.fit(X)
Z = PCA_model.transform(X)

SPCA_model = SparsePCA(n_components=64, alpha=0, ridge_alpha=1e-6, max_iter=200, tol=1e-3)
SPCA_model.fit(X)
SZ = SPCA_model.transform(X)

# Check the first 2 loadings from PCA and SPCA. They are supposed to agree.
print PCA_model.components_[0:2]
print SPCA_model.components_[0:2]

# Check the first 2 observations of transformed data. They are supposed to agree.
print Z[0:2]
print SZ[0:2]

当 lasso 惩罚大于 0 时，SparsePCA 的结果仍然与 R 给我的结果有很大不同，后者基于人工检查和我从原始论文中学到的东西是正确的。那么，是 SparsePCA 坏了，还是我错过了什么？

4

2 回答 2

3

通常：有许多不同的公式和实现。

sklearn 正在使用具有不同特征的不同实现。

让我们看看它们有什么不同：

所以看起来 sklearn 至少在基于 l2-norm 的组件方面做了一些不同的事情（它不见了）。

这是设计使然，因为这是字典学习领域内的基本形式：（由 sklearn 链接的算法论文用于实现）。

在此处输入图像描述

很有可能，当稀疏参数为零时，这种替代公式不能保证（或根本不关心）模拟经典 PCA（这并不奇怪，因为这些问题在优化理论和sparsePCA 必须驻留在一些基于启发式的算法中，因为问题本身是 NP-hard，ref )。通过这里对等价定理的描述，这个想法得到了加强：

在此处输入图像描述

于 2017-01-14T22:18:16.560 回答
0

答案不一样。首先，我认为它可能是求解器，但检查不同的求解器，我得到几乎相同的载荷。看到这个：

nr.seed(1)
MDP = MDPGen(dim=16, rep=1, n=30, T=90, sd=0.5, nsd=0.1, TranType=None)
X = MDP.get('X').astype(np.float32)

PCA_model = PCA(n_components=10,svd_solver='auto',tol=1e-6)
PCA_model.fit(X)

SPCA_model = SparsePCA(n_components=10, alpha=0, ridge_alpha=0)                       
SPCA_model.fit(X)

PC1 = PCA_model.components_[0]/np.linalg.norm(PCA_model.components_[0])
SPC1 = SPCA_model.components_[0].T/np.linalg.norm(SPCA_model.components_[0])
print(np.dot(PC1,SPC1))

import pylab
pylab.plot(PC1)
pylab.plot(SPC1)
pylab.show()
于 2018-04-04T23:02:22.277 回答