当'n'很大时计算帽子矩阵（60000）答案

【问题标题】：Calculating the hat matrix when 'n' is large (60000)当'n'很大时计算帽子矩阵（60000）
【发布时间】：2022-01-26 19:16:43
【问题描述】：

我正在尝试在 python 中计算“帽子矩阵”。我正在使用以下公式。当 X 的长度很大（比如 60,000）时，我会遇到内存不足的问题。 H = X*inv(X'X)*X'

有没有一种计算上更有效的方法来做到这一点？代码片段如下

import numpy as np
X = np.random.rand(60000,1)
hat = X.dot(np.linalg.inv(X.T.dot(X)).dot(X.T))

【问题讨论】：

我没有任何解决方案，但对于初学者，您是否考虑过如果使用 float64 (60_000**2 * 8 / 1e9)，hat 矩阵的大小将是 28.8GB。也许你可以看看使用 SciPy 稀疏矩阵
你确定你需要明确的矩阵吗？人们通常可以在不明确形成对象的情况下求解线性系统，例如，使用np.linalg.solve。
我不喜欢测试内存错误，但如果我解决了这个问题，我会开始逐步评估它，以便清楚地了解错误出现的位置。换句话说，以X.T.dot(X) 开头。然后是inv，以此类推。但简单地估计每一步的数组大小可能就足够了。
可能存在特定类型的矩阵 X，例如疏。您需要分享您对 X 的所有了解，以便人们能够提供帮助。 Bnaecker 的问题非常中肯：你真的需要所有的帽子矩阵吗？做什么的？例如，如果您只需要不需要计算帽子矩阵的状态向量，并且这样做确实可能会失去准确性。

标签： python numpy linear-algebra

【解决方案1】：

我同意cmets中所说的，可能这是其他一些计算的中间结果。

您可以在自定义类中对其进行抽象。感谢操作符覆盖，您可以非常自然地实现使用它。

这里我给你一个类，它允许你对结果矩阵进行矩阵切片，或者乘以左边的向量。

class HatMatrix:
    def __init__(self, X):
        m,n = X.shape;
        self.X = X;
        assert m >= n;
    
    def __array__(self):
        X = self.X;
        return X @ np.linalg.inv(X.T.conj() @ X) @ X.T.conj()
    
    def __getitem__(self, indices):
        if not isinstance(indices, tuple):
            I = indices;
            J = slice(None,None);
        else:
            I,J = indices;
        X = self.X;
        XJ = np.atleast_2d(X[J,:]);
        XI = np.atleast_2d(X[I,:]);
        return XI @ np.linalg.inv(X.T.conj() @ X) @ XJ.T.conj()
    def __matmul__(self, v):
        X = self.X
        return (X @ np.linalg.inv(X.T.conj() @ X)) @ (X.T.conj() @ v)

这是一个如何使用它的示例

X = np.random.rand(60000,1)
hat = HatMatrix(X)
hat[0] # the first row
hat[30:40, 10:40] # a submatrix
hat[50000, 50000] # a single element
hat @ (X**4); # multiply hat by a 60000 x 1 vector

上面的例子在这里运行5ms。

【讨论】：