SVT矩阵补全算法的实现及应用

1 前言

  最近在学习矩阵补全的方法主要是用来做药物重定位。入门矩阵补全方法后才发现这个坑有点大，需要太多的数学基础了。对于数学知识严重不足的我欲哭无泪，搞了两周之后对这个方法的现实意义跟数学背景有了一定的了解。在这里做个总结并用经典的SVT矩阵补全算法（《A singular value thresholding algorithm for matrix completion》）作为一个Demo（Python实现），解释一下矩阵补全的具体做法，感受一下矩阵补全到底是干啥的。

2 低秩矩阵补全的现实意义

  这一部分内容参考我之前的博客（矩阵分解方法概述）
  
时间仓促没大块时间，码字先把代码贴出来，以后慢慢补解释。

""
@Date   :2020/11/15 19:54
@Source 《A singular value thresholding algorithm for MC&RW》
"""
import numpy as np

#  creating data
svt_data = np.array([[0, 3, 0, 4],
                    [3, 0, 4, 0],
                    [0, 0, 2, 0],
                    [5, 0, 3, 4],
                    [0, 0, 4, 0],
                    [0, 3, 3, 0]])


#  changing data type into float
svt_data = svt_data.astype(float)
# print(svt_data)

#  generating Omega :0 denotes None 1 denotes true
shape = svt_data.shape
Omega = np.zeros(shape)
for i in range(0, shape[0]):
    for j in range(0, shape[1]):
        if svt_data[i, j] > 0:
            Omega[i, j] = 1
# print(Omega)

def svt_solve(A, Omega, tau=None, delta=None, epslion=1e-2, max_iterations=1000):
    #  矩阵初始化，生成一个和矩阵A形状一样的0矩阵
    Y = np.zeros_like(A)

    if not tau:
        tau = 5 * np.sum(A.shape) / 2
    if not delta:
        #  确定步长初始值
        delta = 1.2 * np.prod(A.shape) / np.sum(Omega)
    for _ in range(max_iterations):
        #  对矩阵Y进行奇异值分解
        U, S, V = np.linalg.svd(Y, full_matrices=False)
        #  soft-thresholding operator
        print(type(S))
        print(type(tau))
        print(tau)
        S = np.maximum(S - tau, 0)
        #  singular value shrinkage
        X = np.linalg.multi_dot([U, np.diag(S), V])
        #  Y的迭代
        Y += delta * Omega * (A-X)
        #  误差计算
        rel_recon_error = np.linalg.norm(Omega * (X-A)) / np.linalg.norm(Omega*A)
        if rel_recon_error < epslion:
            break
    return X


result = svt_data_hat = svt_solve(svt_data, Omega)
print(result)