SPCArt算法,利用旋转(正交变换更为恰当,因为没有体现出旋转这个过程),交替迭代求解sparse PCA。
对以往一些SPCA算法复杂度的总结
在这里插入图片描述
注:
Notation
在这里插入图片描述
论文概述
就是普通PCA的前
个载荷向量(loadings,按照特征值降序排列)
也是彼此正交的,张成同一子空间的向量组。
原始问题
在这里插入图片描述
如果能解出来,当然好,可是这是一个很难求解的问题,所以需要改进。
问题的变种
直接用
表示了,为了符号的简洁。
在这里插入图片描述
变成这个问题之后,我们所追求的便是
1.我们希望
2.
算法
在这里插入图片描述
这是一个交替迭代算法,我们来分别讨论。
固定
,计算
当固定,之后,问题就退化为:
在这里插入图片描述
这个问题在Sparse Principal Component Analysis(Zou 06)这篇论文里面也有提到。
上述最小化问题,可以变换为
若
就是要最大化:
当(注意
是正交矩阵)。
固定
,求解
(
)
1-范数
在这里插入图片描述
注意:
经过转换,上述问题还等价于:
通过分析(蛮简单的,但是不好表述),可以得到:
(新的初始问题)
在这里插入图片描述
等价于:
显然,若,
,此时函数值为
若,值为
,所以,为了最小化值,取:
,也就是说,
否则,
T-sp 考虑稀疏度的初始问题
在这里插入图片描述
等价于:
在条件不变的情况下。
证明挺简单的,但不好表述,就此别过吧。
最优解是:
T-en 考虑Energy的问题
文章到此并没有结束,还提及了一些衡量算法优劣的指标,但是这里就不提了。大体的思想就在上面,我认为这篇论文好在,能够把各种截断方法和实际优化问题结合在一起,很不错。
代码
def Compute_R(X, V):
W, D, Q_T = np.linalg.svd(X.T @ V)
return W @ Q_T
def T_S(V, R, k): #k in [0,1)
Z = V @ R.T
sign = np.where(Z < 0, -1, 1)
truncate = np.where(np.abs(Z) - k < 0, 0, np.abs(Z) - k)
X = sign * truncate
X = X / np.sqrt((np.sum(X ** 2, 0)))
return X
def T_H(V, R, k): #k in [0,1) 没有测试过这个函数
Z = V @ R.T
X = np.where(np.abs(Z) > k, Z, 0)
X = X / np.sqrt((np.sum(X ** 2, 0)))
return X
def T_P(V, R, k): #k belongs to {0, 1, 2, ..., (p-1)} 没有测试过这个函数
Z = V @ R.T
Z[np.argsort(np.abs(Z), 0)[:k], np.arange(Z.shape[1])] = 0
X = Z / np.sqrt((np.sum(Z ** 2, 0)))
return X
def Main(C, r, Max_iter, k): #用T_S截断 可以用F范数判断是否收敛,为了简单直接限定次数
value, V_T = np.linalg.eig(C)
V = V_T[:r].T
R = np.eye(r)
while Max_iter > 0:
Max_iter -= 1
X = T_S(V, R, k)
R = Compute_R(X, V)
return X.T
结果,稀疏的程度大点,反而效果还好点。