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

写在前面

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

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

拟牛顿法(Quasi-Newton method)

公式推导

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

无约束优化问题：$$min{D}_E(R||PQ)=\frac{1}{2}||R-PQ|{|}^2$$

代价函数$D_E(R||PQ)$ 的黑塞(Hessian)矩阵$H_A={\nabla }_P^2(D_E)=I_{I\times I}\otimes Q{Q}^T\in R^{IK\times IK}$ 是一个分块对角矩阵，对角线上的块矩阵为$Q{Q}^T$。

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

$$P=P-{\nabla }_P(D_E(R||PQ))H_P^{-1}$$

以及

$$Q=Q-{\nabla }_Q(D_E(R||PQ))H_Q^{-1}$$

其中${\nabla }_P(D_E(R||PQ))=(PQ-R)Q^T$

可得：

$$P=P-(PQ-R)Q^T(Q{Q}^T{)}^{-1}$$

以及

$$Q=Q-(PQ-R)P^T(P{P}^T{)}^{-1}$$

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

$$P=P-(PQ-R)Q^T(Q{Q}^T+\lambda I_{K\times K}{)}^{-1}$$

以及

$$Q=Q-(PQ-R)P^T(P{P}^T+\lambda I_{K\times K}{)}^{-1}$$

其中$I_{K\times K}$ 就是一个对角线全为1，大小为$K\times 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

显示效果如下

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

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

def newton(R, P, Q, lamb, steps=5000):
    
    
    for step in range(steps):
        Pu = P-(P@Q.T-R)@Q@np.linalg.inv(Q.T@Q+lamb*np.eye(K))
        Qu = Q-(Pu@Q.T-R).T@Pu@np.linalg.inv(Pu.T@Pu+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(())

结果

可以发现，这收敛速度非常快，只需6步即可完成，而文章（1）中则慢很多。

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

plt.imshow((P_estimate@Q_estimate.T)[0].reshape((64,64)), cmap=plt.cm.gray)

结果

其他方法

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

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

公式推导

更新公式：

$$P_{ik}=P_{ik}\frac{[R{Q}^T{]}_{ik}}{[PQ{Q}^T+\alpha P{]}_{ik}}$$

以及

$$Q_{kj}=Q_{kj}\frac{[P^{T}R{]}_{kj}}{[P^{T}PQ+\beta Q{]}_{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

书中还提到过一些不错的方法，比如KL散度，AB散度等，由于水平有限没有写出code。想要对算法进行修改和优化，还可以从修改目标函数（Object function）入手。