行列演算だから速いだろうと油断してはいけない

回帰分析の最小二乗法による係数の推定は
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^Ty をかけて(X^Ty) としてから(X^TX)^{-1} にかけるのが高速である。
 
20180607追記
スーパーポスドクが計算量について教えてくれた。

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

p\times p行列Ap\times n行列Bn\times 1行列Cの掛け算を行う場合、

(AB)Cと計算すると、ABp\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)と計算すると、BCp\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