共役勾配法をＲで

たまには R のコード書いとかないと忘れる。
ただでさえ R はいろいろ特殊だってのに。

というわけで、勉強中の共役勾配法(conjugate gradient method)を R で書いてみた。といっても、pseudo code をそのまま落とし込んだだけなのだが。しかも線形。
読んでいるのはこれ。

• Jonathan Richard Shewchuk, "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain"

「苦痛を伴わない共役勾配法入門」、略して「サルでもわかるＣＧ法」。
ほんとわかりやすくて、びっくり。

読者が知りたい内容にたどり着くにはこの順番で読んだらいいよ、というダイアグラムまで付いている。
例えば(線形)共役勾配法は、4. Steepest Descent(最急降下法) → 7. Conjugate Directions → 8. Conjugate Gradient と読めばいいよ、という感じ*1。
親切すぎる。

そして Appendix には各種 pseudo code が載っていて。
そして R なら、この pseudo code をそのまま逐語訳すれば、動く。

# "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain" Appendix B2.

# sample
A <- matrix(c(3, 2, 2, 6), 2);
b <- c(2, -8);
x <- runif(2, -100, 100);  # random start vector
print(x);

# CG method
r <- b - A %*% x;
d <- r;
delta_new <- sum(r ^ 2);
delta_0 <- delta_new;
for (i in 1:nrow(A)) {
	q <- A %*% d;
	alpha <- delta_new / sum(d * q);
	x <- x + alpha * d;
	r <- r - alpha * q;
	delta_old <- delta_new;
	delta_new <- sum(r ^ 2);
	beta <- delta_new / delta_old;
	d <- r + beta * d;
	print(i);
	print(as.vector(d));
	print(as.vector(x));
}

これでどうして最小値を与えるベクトルが正確に求まってしまうのか、しかもでたらめな開始点から始めても高々 (空間の次元)回の繰り返しで、というあたりは "An Introduction to the Conjugate Gradient Method Without the Agonizing Pain" よりわかりやすく書く自信はとてもないので、そちらを見ていただくとして。
まず Steepest Descent を説明して、その流れで「どうして『共役性』という概念を考え出したくなったか」という動機を導いてくれているので、とても納得させてくれるあたりがポイント。

*1:Notation とか Quadratic Form とか、当たり前のは除く