numpy - 高效计算帽子矩阵: inv(X'WX )'X' 的对角线

Question

作为疾病风险模型的一部分，我试图从一篇论文(在 Python/numpy 中)实现计算，其中一部分是以下矩阵计算:



在哪里:

X an [n.m]。 m 的大小合理 (1..200)，但 n 相当大 (>500k)

W 是 [n.n]，并且过大，但在对角线上只有非零值

另外，我只需要获取 Q 的对角线元素作为输出。

是否有一些 numpy 矩阵魔法可以让我有效地计算它？

笔记:

paper有一个 implementation in R我(相信)这样做如下:

Qdiag <- lm.influence(lm(y ~ X-1, weights=W))$hat/W

R的 docs for lm.influence$hat说这给出了“一个包含‘帽子’矩阵对角线的向量。”虽然 Wikipedia's definition 听起来有点像我想要的的帽子矩阵(==影响或投影矩阵)，看起来略有不同。

——

我认为以下是一个有效的(幼稚的)实现。对于大 n，它的内存不足

m = 3
n = 20 # 500000 -- out of memory for large n

np.random.seed(123)
X = np.random.random((n,m))
W = np.random.random(n)
W = np.diag(W)
xtwx = X.T.dot(W.dot(X))
xtwxi = np.linalg.pinv(xtwx)
xtwxixt = xtwxi.dot(X.T)
Q = X.dot(xtwxixt)
Qdiag = np.diag(Q)
print Qdiag.shape, Qdiag.sum() # Checksum of output 
print Qdiag

Answer 1

(在现已删除的评论中，我说这是不可能的，因为某些密度和硬件假设将其视为黑盒。但似乎可以做到。这并不意味着这是正确的方法!)

因此，在不分析这个公式的背景的情况下，我们可以在给定最小假设和经典规则的情况下做一些基本方法，例如:

A:矩阵乘法的结合性

B:求解线性方程组而不是求逆

我们假设 XtWX 是非奇异的

C:识别 A*W(仅 W 对角线)是具有对角线向量

的逐行元素乘积

D:仅计算 Q 的对角线条目(否则我们将为您的数字生成 N*N = 2.5e8 个条目的结果矩阵)

我用过 this

代码:

import numpy as np
from time import perf_counter as pc     # python 3 only

m = 200
n = 500000

np.random.seed(123)
X = np.random.random((n,m))
W_diag = np.random.random(n)            # C -> dense vector

start_time = pc()

lhs = np.multiply(X.T, W_diag).dot(X)   # C (+A)
x = np.linalg.solve(lhs, X.T)           # B

# EDIT: Paul Panzer recommends the inverse in his comment based on the rhs-dims!

# if you know something about lhs (looks symmetric; maybe even PSD)
# use one of the following for speedups and more robustness
# i highly recommend this research: feels PSD
# import scipy.linalg as slin
# x = slin.solve(lhs, X.T, assume_a='sym')
# x = slin.solve(lhs, X.T, assume_a='pos')

Q_ = np.einsum('ij,ji->i', X,x)         # D most important -> N*N elsewise
print(Q_)

end_time = pc()
print(end_time - start_time)

出去:

[ 0.00068103  0.00083676  0.00080945 ...,  0.00077864  0.00078945
  0.0007804 ]
3.1077745566331165  # seconds

与您为单个测试用例给出的代码相比，结果是相同的!

一般来说，我会建议通过基础数学而不是提取公式本身。由于论文基本上说更完整的问题是加权最小二乘问题，因此对于某些研究来说，这是一个良好的开端。