Drineas P, Kannan R, Mahoney M W, et al. Fast Monte Carlo Algorithms for Matrices II: Computing a Low-Rank Approximation to a Matrix[J]. SIAM Journal on Computing, 2006, 36(1): 158-183.

问题

我们有一个矩阵，我们需要对其进行矩阵的分解，很完美很经典的一种方法就是SVD,但是这种方法 的缺憾在于，需要的计算量比较大。不妨设的奇异值分解为：

其中：，, 。
假设，那么，矩阵的零空间，矩阵的值域为
那么可以有下面的方法表示：

到这里，我们简单介绍了SVD。回到正题，为了避免计算量大的问题，这篇文章提出了一种基于蒙特卡洛采样的矩阵分解的算法。

算法

为什么可以这么采样，以及概率的选择，在FAST MONTE CARLO ALGORITHMS FOR MATRICES I中有介绍。算法的思想很朴素，但是通篇的证明让人抓耳挠腮。

LINEARTIMESVD 算法

在这里插入图片描述

CONSTANTTIMESVD 算法

在这里插入图片描述

理论

俩个算法，作者都给除了形如下的界（大概率）：

在这里插入图片描述

算法1的理论

作者先给出的是下面的证明，

在这里插入图片描述

我们先来看第一部分的证明，这部分只是简单地利用了的性质。

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

到此关于范数的一个理论就得到了，接下来作者给出了关于范数的性质。

在这里插入图片描述

第一部分的证明，不等式部分利用了三角不等式，及的性质。最后一个等式成立的原因是

在这里插入图片描述

在这里插入图片描述

在这里插入图片描述

接下来，还有定理4：

在这里插入图片描述

在这里插入图片描述

这部分的证明，需要利用FAST MONTE CARLO ALGORITHMS FOR MATRICES I 中的性质，这里便不讲了。

算法2 的理论

我们只给出了结果，证明实在有些长。

在这里插入图片描述

代码

import numpy as np


class FastSVD:

    def __init__(self, A):
        self.m, self.n = A.shape
        self.A = np.array(A, dtype=float)
        self.norm_F = FastSVD.forbenius(self.A)

    @classmethod
    def forbenius(cls, A):
        """矩阵A的F范数"""
        return np.sum(A ** 2)

    @classmethod
    def approx_h(cls, A):
        """A=UDV^T, 我们要U"""
        value, vector = np.linalg.eig(A.T @ A)
        U = []
        for i in range(len(value)):
            if value[i] < 1e-15:
                break
            else:
                U.append(A @ vector[:, i] / np.sqrt(value[i]))
        return np.array(U).T

    def fastSVD(self, c):
        """返回的H的每一列是我们所需要的"""
        assert isinstance(c, int), "{0} is not an integer"
        p = np.array([self.A[:, i] @ self.A[:, i] / self.norm_F for i in range(self.n)])
        lucky_dog = np.random.choice(np.arange(self.n), size=c, replace=True, p=p)
        C = np.zeros((self.m, c))
        for t, dog in enumerate(lucky_dog):
            C[:, t] = self.A[:, dog] / np.sqrt(c * p[dog])
        H = FastSVD.approx_h(C)
        return H