如何避免处理大行数自变量矩阵时帽矩阵迹计算中的内存错误？

Question

我有一个独立变量矩阵，叫做X。现在我想要计算这个矩阵的“帽子矩阵”的迹，或者找到一种更简单的方法来获取这个迹，而不需要实际计算帽子矩阵。问题是，X有14826行。

res = glm_binom.fit()
YHatTemp = res.mu
HatMatTemp = X*res.pinv_wexog

（或者，可以把第三行替换成）

HatMatTemp = X*np.linalg.inv(np.transpose(X)*X)*np.transpose(X)

上面的内容给了我以下结果：

File "C:\Python27\lib\site-packages\numpy\matrixlib\defmatrix.py", line 341, in __mul__
return N.dot(self, asmatrix(other))
MemoryError

我之所以想要帽子矩阵的迹，是为了计算GCV标准，以便进行模型选择。所以一旦这个方法有效，我会多次重复这个过程。

我非常希望能找到解决这个问题的方法，或者干脆绕过它。谢谢！

Answer 1

对于你提到的帽子矩阵，hat = X.dot(np.linalg.inv(X.T.dot(X)).dot(X.T))，它的迹（trace）其实就是rank(X)。所以，如果X的列满秩，并且列数少于行数，那么这个值就是X.shape[1]。你在这里找到的是无偏线性估计器的自由度，也就是普通最小二乘法的自由度。

不过，你想要进行广义交叉验证（GCV），这就需要完全了解帽子矩阵的对角线。你可以不计算整个帽子矩阵，而是通过以下方式来计算对角线：

import numpy as np
rng = np.random.RandomState(42)
n_samples, n_features = 20, 5
X = rng.randn(n_samples, n_features)

hat = X.dot(np.linalg.inv(X.T.dot(X)).dot(X.T))

hat_diag = np.diagonal(hat)
trace = hat_diag.sum()  # this is equal to n_features == 5
print trace

only_diag = np.einsum('ij, ij -> j', X.T, np.linalg.inv(X.T.dot(X)).dot(X.T))

print (only_diag == hat_diag).all()  # this evaluates to True

从上面可以看出，only_diag只包含了对角线的元素。你也可以用一种更简单的方法来计算，不过可能会占用更多的内存，方法如下：

only_diag = (X.T * np.linalg.inv(X.T.dot(X)).dot(X.T)).sum(0)