非负矩阵分解NMF(2): 拟牛顿法及其他方法

程序员文章站 2024-03-17 20:34:40

...

文章目录

写在前面

在之前的一篇文章中非负矩阵分解（1）已经介绍了NMF的基本概念，数学基础以及代码实现，这里不再对NMF基本概念做更多讲解，建议先看完（1）再看这篇。

实现非负矩阵分解主要两种方法，一个是梯度下降法，一个是乘法迭代法，乘法迭代可以看作是梯度下降的一种变形。文章（1）中使用的就是乘法迭代，而这一篇的重点则是梯度下降

逆牛顿法(Quasi-Newton method)

公式推导

理论参考至《矩阵分析于应用（第二版》（张贤达著）第373页。

无约束优化问题： m i n D E ( R ∣ ∣ P Q ) = 1 2 ∣ ∣ R − P Q ∣ ∣ 2 minD_E(R||PQ)=\frac{1}{2}||R-PQ||^2 minDE(R∣∣PQ)=21∣∣R−PQ∣∣2

代价函数 D E ( R ∣ ∣ P Q ) D_E(R||PQ) DE(R∣∣PQ) 的黑塞(Hessian)矩阵 H A = ∇ P 2 ( D E ) = I I × I ⊗ Q Q T ∈ R I K × I K H_A=\nabla^2_P(D_E)=I_{I \times I}\otimes QQ^T\in R^{IK \times IK} HA=∇P2(DE)=II×I⊗QQT∈RIK×IK 是一个分块对角矩阵，对角线上的块矩阵为 Q Q T QQ^T QQT。

于是梯度下降更新法则为：

P = P − ∇ P ( D E ( R ∣ ∣ P Q ) ) H P − 1 P = P-\nabla_P(D_E(R||PQ))H_P^{-1} P=P−∇P(DE(R∣∣PQ))HP−1

以及

Q = Q − ∇ Q ( D E ( R ∣ ∣ P Q ) ) H Q − 1 Q = Q-\nabla_Q(D_E(R||PQ))H_Q^{-1} Q=Q−∇Q(DE(R∣∣PQ))HQ−1

其中 ∇ P ( D E ( R ∣ ∣ P Q ) ) = ( P Q − R ) Q T \nabla_P(D_E(R||PQ))=(PQ-R)Q^T ∇P(DE(R∣∣PQ))=(PQ−R)QT

可得：

P = P − ( P Q − R ) Q T ( Q Q T ) − 1 P = P-(PQ-R)Q^T(QQ^T)^{-1} P=P−(PQ−R)QT(QQT)−1

以及

Q = Q − ( P Q − R ) P T ( P P T ) − 1 Q = Q-(PQ-R)P^T(PP^T)^{-1} Q=Q−(PQ−R)PT(PPT)−1

为了防止矩阵 Q Q T QQ^T QQT奇异或者条件数很大，可以采用松弛法。于是最终更新公式为：

P = P − ( P Q − R ) Q T ( Q Q T + λ I K × K ) − 1 P=P-(PQ-R)Q^T(QQ^T+\lambda I_{K\times K})^{-1} P=P−(PQ−R)QT(QQT+λIK×K)−1

以及

Q = Q − ( P Q − R ) P T ( P P T + λ I K × K ) − 1 Q=Q-(PQ-R)P^T(PP^T+\lambda I_{K\times K})^{-1} Q=Q−(PQ−R)PT(PPT+λIK×K)−1

其中 I K × K I_{K \times K} IK×K 就是一个对角线全为1，大小为 K × K K \times K K×K的矩阵。

代码实现

导入package

!pip install update scikit-image

import numpy as np
from sklearn.datasets import fetch_olivetti_faces
import matplotlib.pyplot as plt
from time import time
from skimage.transform import resize

导入data

faces = fetch_olivetti_faces()
data_faces = faces.data
images_faces = faces['images']
number_of_train = 100
images_faces_train = images_faces[:number_of_train,:,:]

data_faces_train = data_faces[:number_of_train,:]

n_samples = len(images_faces_train)
image_shape = images_faces_train[0].shape

注意这里和文章（1）的不同是没有reshape图像成(16,16)的格式，而是直接使用原图像(64,64)

查看数据

data_faces_centered = data_faces_train - data_faces_train.mean(axis=0)

# local centering
data_faces_centered -= data_faces_centered.mean(axis=1).reshape(n_samples, -1)

# Let's show some centered faces
plt.figure(figsize=(20, 2))
plt.suptitle("Centered Olivetti Faces", size=16)
for i in range(10):
    plt.subplot(1, 10, i+1)
    plt.imshow(data_faces_centered[i].reshape(image_shape), cmap=plt.cm.gray)
    plt.xticks(())
    plt.yticks(())
    
R = data_faces_centered

显示效果如下
非负矩阵分解NMF(2): 拟牛顿法及其他方法

算法实现
注意代码中的Q=公式推导中的Q.T。

P的大小是(N,K), Q的大小是(M,K)

def newton(R, P, Q, lamb, steps=5000):
    
    
    for step in range(steps):
        Pu = P-(aaa@qq.com.T-R)@aaa@qq.com.linalg.inv(Q.aaa@qq.com+lamb*np.eye(K))
        Qu = Q-(aaa@qq.com.T-R).aaa@qq.com@np.linalg.inv(Pu.aaa@qq.com+lamb*np.eye(K))
        e_P = np.sqrt(np.sum((Pu-P)**2, axis=(0,1)))/P.size
        e_Q = np.sqrt(np.sum((Qu-Q)**2, axis=(0,1)))/Q.size
        if e_P<0.001 and e_Q<0.001:
            print("step is:",step)
            break
        P = Pu
        Q = Qu
    return P, Q

测试，这里的 λ \lambda λ 设置的0.001

N = len(R)
M = len(R[0])
K = 2
rng = np.random.RandomState(1)
P = rng.rand(N,K)
Q = rng.rand(M,K)

P_estimate, Q_estimate = newton(R, P, Q, 0.001) 

plt.figure(figsize=(6, 3))
plt.suptitle("learned basises from NMF", size=16)
for i in range(K):
    plt.subplot(1, K, i+1)
    plt.imshow(Q_estimate[:,i].reshape(64,64), cmap=plt.cm.gray)
    plt.xticks(())
    plt.yticks(())

结果
非负矩阵分解NMF(2): 拟牛顿法及其他方法
可以发现，这收敛速度非常快，只需6步即可完成，而文章（1）中则慢很多。

查看重组后的矩阵 R ^ \hat{R} R^中的图像

plt.imshow((aaa@qq.com_estimate.T)[0].reshape((64,64)), cmap=plt.cm.gray)

结果
非负矩阵分解NMF(2): 拟牛顿法及其他方法

其他方法

在《矩阵分析于应用（第二版》（张贤达著）第372页还介绍了交替非负最小二乘算法（alternating nonnegative least squares, ANLS）

中间的公式推导不再赘述，其最终得到的是加上了约束项的乘法迭代，跟文章（1）的方式较为相似。

公式推导

更新公式：

P i k = P i k [ R Q T ] i k [ P Q Q T + α P ] i k P_{ik}=P_{ik}\frac{[RQ^T]_{ik}}{[PQQ^T+\alpha P]_{ik}} Pik=Pik[PQQT+αP]ik[RQT]ik

以及

Q k j = Q k j [ P T R ] k j [ P T P Q + β Q ] k j Q_{kj}=Q_{kj}\frac{[P^TR]_{kj}}{[P^TPQ+\beta Q]_{kj}} Qkj=Qkj[PTPQ+βQ]kj[PTR]kj

代码实现

def ALS(R, P, Q, alpha=0.001, beta=0.001, steps=5000):

    for step in range(steps):
        

        Pu = P*(R.dot(Q))/(P.dot(Q.T).dot(Q)+alpha*P)+1e-7
        Qu = (Q.T*(Pu.T.dot(R))/(Pu.T.dot(Pu).dot(Q.T))+beta*Q.T).T+1e-7
        
        e_P = np.sqrt(np.sum((Pu-P)**2, axis=(0,1)))/P.size
        e_Q = np.sqrt(np.sum((Qu-Q)**2, axis=(0,1)))/Q.size
        if e_P<0.001 and e_Q<0.001:
            print("step is:",step)
            break
        
        P = Pu
        Q = Qu
    return P, Q