ハット行列の対角を効率的に計算する：inv（X'WX） 'X'

Question

私は病気のリスクモデルの一部として、（Python/numpyの）論文から計算を実装しようとしています。次の行列演算：

：

• X [NM]。 mは（1..200）適度な大きさであるが、nは非常に大きく（> 500K）である
• Wは[NN]であり、過度に大部分であるが、唯一の対角線

に非ゼロ値を有します

また、Qの対角要素を出力するだけでよいです。

私はこれを効率的に計算することができるいくつかの魔法の行列の魔法がありますか？

注：

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：ある行方向要素単位（W対角のみを有する）* Wを認識する対角ベクトルとの積
• D：Qの対角成分を計算する（または、N * N = 2の結果行列を得る）。私はthis

を使用し、あなたの番号の5E8エントリ）

• はコード：

  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 
  

  結果が与えられたあなたのコードと比較して、同じですあなたの単一のテストケースのために！

  一般に、抽出された数式自体ではなく、基礎となる数学を推薦します。紙は基本的に、より完全な問題は加重最小二乗問題であると述べているため、これはいくつかの研究にとっては良いスタートです。