行列演算だから速いだろうと油断してはいけない - 驚異のアニヲタ社会復帰の予備

回帰分析の最小二乗法による係数の推定は
$b=(X^TX)^{-1}X^Ty$
で求められる。ここで、標本数 $n$ 、パラメータ数 $p$ とすると
$\underbrace{b}_{p\times 1}=(\underbrace{X^TX}_{p\times p})^{-1}\underbrace{X^T}_{p \times n}\underbrace{y}_{n\times 1}$
である。
普通に演算すると、 $(X^TX)^{-1}$ に $X^T$ をかけて $y$ をかける、という順番だが、 $(X^TX)^{-1}(X^Ty)$ として先に $X^T$ と $y$ をかけて $(X^Ty)$ としてから $(X^TX)^{-1}$ にかけるのが高速である。
　
20180607追記
スーパーポスドクが計算量について教えてくれた。

$l\times m$ 行列と $m\times n$ 行列の掛け算では、 $l\times n$ 個の成分を計算することになり、各成分について $m$ 回の掛け算と $m-1$ 回の足し算が必要になる。

$p\times p$ 行列 $A$ 、 $p\times n$ 行列 $B$ 、 $n\times 1$ 行列 $C$ の掛け算を行う場合、

$(AB)C$ と計算すると、 $AB$ は $p\times n$ 行列だから、
掛け算は $p\times p \times n + n\times p\times 1 = p(np+n)$ 回
足し算は $(p-1)\times p\times n + (n-1)\times p\times 1 = p(np-1)$ 回
　
$A(BC)$ と計算すると、 $BC$ は $p\times 1$ 行列だから、
掛け算は $n\times p\times 1 + p\times p\times 1 = p(n+p)$ 回
足し算は $(n-1)\times p\times 1 + (p-1)\times p\times 1 = p(n+p-2)$ 回

前者だと3次の項が出てくるが、後者だと2次まで。

n <- 1000
p <- 100
X <- matrix(runif(n*p), n, p)
y <- matrix(runif(n), n, 1)
iter <- 3000

前からかけていく場合は

system.time(replicate(iter, solve(t(X) %*% X) %*% t(X) %*% y))

   ユーザ   システム       経過  
    64.515     54.919     10.156

先に $(X^Ty)$ を行う場合は

system.time(replicate(iter, solve(t(X) %*% X) %*% (t(X) %*% y)))

   ユーザ   システム       経過  
    41.333     37.563      6.672