python实现矩阵补全和奇异值软阈值算法改进模型

这段代码实现了一个低秩矩阵补全 (Low-Rank Matrix Completion, LRMC) 算法，专门用于填补稀疏矩阵（缺失值矩阵）中的空缺条目。例如，协同过滤推荐系统、图像修复或任何需要从观测数据中恢复缺失数据模式的场景。使用RSE指标可以评估算法预测缺失值的性能。我就不过多描述，有问题可以私信我。

import numpy as np
import pandas as pd

def LRMC(sparse_mat, dense_mat, alpha, rho, maxiter):
    pos_train = np.where(sparse_mat != 0)
    pos_test = np.where((sparse_mat == 0) & (dense_mat != 0))
    binary_mat = sparse_mat.copy()
    binary_mat[pos_train] = 1
    print("pos_test", pos_test)

    X = sparse_mat.copy()
    Z = np.zeros(sparse_mat.shape)
    T = np.zeros(sparse_mat.shape)
    rse = np.zeros(maxiter)

    for it in range(maxiter):
        u, s, v = np.linalg.svd(X + T / rho, full_matrices=False)
        vec = s - alpha / rho
        vec[np.where(vec < 0)] = 0
        Z = np.matmul(np.matmul(u, np.diag(vec)), v)
        X = Z - T / rho
        X[pos_train] = sparse_mat[pos_train]
        T = T - rho * (Z - X)
        rse[it] = (np.linalg.norm(X[pos_test] - dense_mat[pos_test], 2)
                   / np.linalg.norm(dense_mat[pos_test], 2))
        print(f"Iteration {it + 1}: rse = {rse[it]}")
    return X, rse


def RM(missing_rate, dense_mat):
    dim = dense_mat.shape
    sparse_tensor = dense_mat * np.round(np.random.rand(*dim) + 0.5 - missing_rate)
    return sparse_tensor


# 示例用法
dense_mat = pd.read_csv('sz_speed.csv', header=None).values  # 假设原始矩阵是一个随机矩阵
missing_rate = 0.2
sparse_mat = RM(missing_rate, dense_mat)
alpha = 0.1
rho = 0.005
maxiter = 200
mat_hat, rse_lr = LRMC(sparse_mat, dense_mat, alpha, rho, maxiter)

# 打印MAE结果
print("X", mat_hat)
print("Final MAE for completed values:", rse_lr[-1])
result_df = pd.DataFrame(mat_hat)

result_df.to_csv('result_matrix.csv', index=False)

print("Result matrix saved to result_matrix.csv")

改进后的模型如下(主要是加入L正则项)：

新代码实现了一种使用交替方向乘子法 (ADMM) 的矩阵补全方法，与之前的低秩矩阵补全 (LRMC) 方法有一些区别。在LRMC中，只用SVD进行奇异值的软阈值化。在ADMM中，同样有奇异值处理，但是它将矩阵补全问题分为多个子问题。此外ADMM算法引入了图的邻接矩阵 A 和由此生成的拉普拉斯矩阵 L（利用图结构信息进行约束），这在处理如图结构数据时，例如社交网络数据，可以通过反映节点之间的关系来辅助补全任务。

# 随机缺失的处理
 def RM(self, missing_rate, dense_mat):
     dim = dense_mat.shape
     sparse_tensor = dense_mat * np.round(np.random.rand(*dim) + 0.5 - missing_rate)
     return sparse_tensor

def admm_matrix_completion(self, dense_mat, A, gamma, lambd, maxiter, missing_rate, epsilon=1e-6):
        """
        ADMM矩阵补全算法
        """
        sparse_mat = self.RM(missing_rate, dense_mat)
        pos_train = np.where(sparse_mat != 0)
        pos_test = np.where((sparse_mat == 0) & (dense_mat != 0))
        X = sparse_mat.copy()
        Z = np.zeros(sparse_mat.shape)
        W = np.zeros(sparse_mat.shape)
        rse = np.zeros(maxiter)
        # 拉普拉斯矩阵（L=D-A)
        L = np.diag(A.sum(axis=0)) - A
        

首先调用随机缺失处理方法 RM 生成一个稀疏矩阵 sparse_mat，其中一些元素根据 missing_rate 被随机置为零。

其中pos_train: 稀疏矩阵中非零元素的索引，用于训练。
      pos_test: 稀疏矩阵中为零且稠密矩阵中非零元素的索引，用于测试，计算重构误差。

然后初始化变量，X为当前估计的矩阵，最初设置为稀疏矩阵。Z为中间变量，W为拉格朗日乘子。rse用于记录每次迭代的相对重构误差。L = np.diag(A.sum(axis=0)) - A是为了创建拉普拉斯矩阵 L，其中 D 是度矩阵，A 是邻接矩阵，L = D - A。

        for it in range(maxiter):
            # X更新公式
            X = self.singular_value_thresholding(Z, W, 1 / lambd)
            X[pos_train] = sparse_mat[pos_train]

            # Z更新公式
            # Z =（γL.T*L+λI）-1 * (λX+W) + Y   (对应文中的公式)
            inv_term = np.linalg.inv(gamma * np.dot(L.T, L) + lambd * 
            np.identity(L.shape[0]))
            Z_Omega_c = np.dot(inv_term, lambd * X + W)
            Z_Omega = sparse_mat
            Z = self.l1_projection(Z_Omega + Z_Omega_c, epsilon)

            # W更新公式
            W = W + lambd * (X - Z)
            rse[it] = (np.linalg.norm(X[pos_test] - dense_mat[pos_test], 2)
                       / np.linalg.norm(dense_mat[pos_test], 2))
            print(f"Iteration {it + 1}: rse = {rse[it]}")
        return X, rse

主迭代循环使用奇异值阈值法对 Z 和 W 进行处理以更新 X，并且保留稀疏矩阵中训练位置的原始值。然后计算更新公式中涉及的逆矩阵更新 Z，通过将当前估计的矩阵 (X 和 W) 输入到 L1 投影中来进行稀疏性约束。接着更新拉格朗日乘子 W 用于调整 Z 的变化。并且加入计算记录当前迭代的相对重构误差（rse），用于监控算法收敛情况。最后返回补全后的矩阵 X 和每次迭代的重构误差 rse。

    def l1_projection(self, Z, epsilon):
        Z[Z < -epsilon] += epsilon
        Z[Z > epsilon] -= epsilon
        Z[(Z >= -epsilon) & (Z <= epsilon)] = 0
        return Z

l1_projection函数实现了L1正则化

    def singular_value_thresholding(self, Z, W, lambd):
        """
        奇异值阈值处理函数
        """
        # 对矩阵 Z - W / lambd 进行奇异值分解（SVD）
        U, S, VT = svd(Z - W / lambd, full_matrices=False)
        # 将奇异值进行软阈值处理，将小于 lambd / 2 的奇异值置零。
        S_thresh = np.maximum(S - lambd / 2, 0)
        # 用处理后的奇异值重新构建矩阵
        return U @ np.diag(S_thresh) @ VT

singular_value_thresholding函数实现了奇异值阈值法。