Woodbury行列反転公式を使って、逆行列演算を高速化。

回帰の問題などでは(X'X)^-1をよく使う。
Xが横長であるときX'Xが非常に大きな行列になるので、
逆行列を求めるのが大変で高速化が重要になる。

そこで、Woodbury行列反転公式(Woodbury matrix inversion formula)
(A+BCD)^-1 = A^-1 - A^-1 B(C^-1+DA^-1 B)^-1DA^-1
という超有名公式を使って、どれくらい高速化できるか実験した。

コード

%% setting
nLoop = 10;
h = 100;    %1000;
w = 2000; %1000;
X = rand(h,w);

%% normal
tic;
for i = 1:nLoop
    Y = X'*X + 0.05*eye(w,w);
    Yi = inv(Y);
end
t1 = (toc/nLoop);

%% woodbury matrix inversion formula
%(A+BCD)^-1 = A^-1 - A^-1 B(C^-1+DA^-1 B)^-1DA^-1
tic;
for i = 1:nLoop
    Ai = 1/0.05*eye(w,w);
    XAi = X*Ai;
    Yi = Ai - Ai*X'/(XAi*X')*XAi;
end
t2 = (toc/nLoop);

%%
disp([t1 t2]);

結果

h=100、w= 100のとき、通常 0.0009s 高速化 0.0012s
h=100、w= 500のとき、通常 0.0320s 高速化 0.0205s
h=100、w=1000のとき、通常 0.1977s 高速化 0.0625s
h=100、w=2000のとき、通常 1.3256s 高速化 0.1901s

通常の場合に逆行列を求める際にinv()を使用しているが、実際の回帰の問題などではバックスラッシュ(/)などの演算子を使うので、もう少し早くなると思われる。

結論

行列の縦横比が5以上であれば、高速化の価値があると思われる。

追記(2013-02-25)

コメントで指摘をもらったので，訂正．
実験に使ったコードの21行目は

    Yi = Ai - Ai*X'/(XAi*X')*XAi;

ではなく，正しくは

    Yi = Ai - Ai*X'/(eye(h) + XAi*X')*XAi;

結果自体にはそれほど影響はないかと．