線形回帰をPython標準ライブラリだけで：勾配降下法が発散する瞬間を実測
numpyなしのPython標準ライブラリだけで線形回帰を閉形式・バッチGD・ミニバッチ・SGDの4通りで実装。学習率を上げて発散する損失の実測値、スケーリングなしで140,538エポックかかった記録、SGDのseed依存まで載せた。
scikit-learn の fit() を呼べば線形回帰は1行で終わります。ただ、その中で何が起きているかを説明できないと、学習率を変えたら損失が nan になった、という場面で手が止まります。この記事は、そういう「なんとなく使っている」状態から抜けたいエンジニア向けに、numpy を使わず Python の標準ライブラリだけで線形回帰を4通りに実装し、実際に壊れる様子を数字で見せます。結論を先に書くと、閉形式は0.47ミリ秒で真の傾きを誤差0.01以内で当て、勾配降下法はスケーリングありなら2エポック、なしだと140,538エポックかかりました。学習率を限界より6%大きくしただけで、損失は100エポックで5,014から9.3×10¹²に跳ね上がりました。
この記事はシリーズ「機械学習をゼロから（Python標準ライブラリのみ）」の第1回です。続く回でも、損失・ループ・止める条件という同じ考え方を使います：k-meansクラスタリング、ナイーブベイズによるテキスト分類、パーセプトロンからニューラルネットワークへ、TF-IDF/BM25で作る検索エンジン。
実行環境： Python 3.9.6 標準ライブラリのみ（math、random、time、sys）、Apple Silicon Mac（M1 Pro）。非推奨のAPIは使っていないので、Python 3.12 以降でもそのまま動くはずです。乱数はすべて seed を固定しています。
データ：正解がわかっている生成器
学習アルゴリズムを検証するには、正解がわかっているデータが要ります。賃貸物件サイトで見かける3つの項目を真似て、専有面積（20〜150㎡）、築年数（0〜40年）、駅からの距離（0.1〜3.0km）を持つ600件の物件と、その価格を作ります。
価格 = 5.0 + 0.9 × 面積 − 0.4 × 築年数 − 6.0 × 駅距離 ＋ ガウスノイズ（σ = 4）
行は最初からランダムな順番なので、先頭480件を学習用、残り120件をテスト用にします。σ = 4 なのでノイズだけで分散が16あり、どんなモデルでもMSEは16前後より大きくは下がりません。これが判断の基準線になります。
"""Linear regression from scratch: closed form vs gradient descent (Python stdlib only)."""
import math, random, sys, time
TRUE_B, TRUE_W = 5.0, [0.9, -0.4, -6.0] # the answer we hope to recover
def make_data(n=600, seed=42):
rng = random.Random(seed)
X, y = [], []
for _ in range(n):
x = [rng.uniform(20, 150), rng.uniform(0, 40), rng.uniform(0.1, 3.0)] # area, age, km
X.append(x)
y.append(TRUE_B + sum(w * v for w, v in zip(TRUE_W, x)) + rng.gauss(0, 4))
return X, y
def add_bias(X):
return [[1.0] + row for row in X]
def predict(row, w):
return sum(a * b for a, b in zip(row, w))
def mse(A, y, w):
return sum((predict(r, w) - t) ** 2 for r, t in zip(A, y)) / len(y)
add_bias は各行の先頭に定数1.0を足します。切片を「0番目の重み」として扱えるので、以降の式で切片を特別扱いせずに済みます。
閉形式：連立方程式を1回解くだけ
損失を L(w) = (1/n)·Σ(xᵢ·w − yᵢ)² と書き、重み wⱼ で偏微分します。
∂L/∂wⱼ = (2/n)·Σ(xᵢ·w − yᵢ)·xᵢⱼ
全部の j をまとめると ∇L = (2/n)·Xᵀ(Xw − y) です。この1行が勾配降下法のすべてでもあるので覚えておいてください。最小点では勾配が0なので、正規方程式 XᵀX w = Xᵀy が得られます。重みは4つなので XᵀX は4×4行列で、ガウスの消去法で解きます。逆行列は作りません。解く方が計算量も少なく、精度もよいからです。
# ---------- 1. closed form: solve (X^T X) w = X^T y ----------
def solve(M, v):
"""Gaussian elimination with partial pivoting."""
d = len(v)
M = [row[:] + [v[i]] for i, row in enumerate(M)]
for c in range(d):
p = max(range(c, d), key=lambda r: abs(M[r][c]))
M[c], M[p] = M[p], M[c]
for r in range(c + 1, d):
f = M[r][c] / M[c][c]
for k in range(c, d + 1):
M[r][k] -= f * M[c][k]
w = [0.0] * d
for r in range(d - 1, -1, -1):
w[r] = (M[r][d] - sum(M[r][k] * w[k] for k in range(r + 1, d))) / M[r][r]
return w
def closed_form(A, y):
d = len(A[0])
XtX = [[sum(r[i] * r[j] for r in A) for j in range(d)] for i in range(d)]
Xty = [sum(r[i] * t for r, t in zip(A, y)) for i in range(d)]
return solve(XtX, Xty)
出力（重みの並びは [切片, 面積, 築年数, 駅距離]）：
Python 3.9.6 train 480 test 120
true [5.000, 0.900, -0.400, -6.000]
closed [4.063, 0.909, -0.394, -6.005] train 16.6759 test 14.6532
true w train 16.8045 test 14.6847
target: train MSE <= 16.6926
3つの傾きは真の値から0.01以内に入りました。切片は1近くずれていますが、バグではありません。切片は「面積0㎡・築0年・駅0km」の物件の予測価格で、データには20㎡未満の物件が1件もありません。そこまで外挿すると、傾きのわずかな誤差が拡大されます。回帰は予測で評価するもので、切片がそれらしいかで評価するものではありません。
学習MSE 16.68 はノイズの下限16のすぐ上です。テストMSEの方が低い（14.65）理由は true w の行でわかります。生成器の真の重みそのものでも、学習用では16.80、テスト用では14.68でした。テストの120件にたまたま小さいノイズが乗っただけです。また、学習データでは当てはめた重みが真の重みより良い値（16.68 < 16.80）を出しています。最小二乗法がノイズを少し拾っているからで、テストではほぼ同点（14.65対14.68）に戻ります。
target の行は、以降のすべての勾配降下法の目標です。学習MSEが閉形式の最適値から0.1%以内に入ったら「到達」とします。
バッチ勾配降下法：同じ微分から作る
勾配降下法は、さっきの勾配を0と置いて解く代わりに、坂を少しずつ下ります。
w ← w − η·(2/m)·Xᵀ(Xw − y)
η は学習率、m は1回の更新に使う行数です。全行を使えばバッチGD、ランダムに選んだ32行ならミニバッチ、1行ずつなら確率的勾配降下法（SGD）です。つまり batch 引数が違うだけの1つの関数で書けます。
# ---------- 2/3. batch, mini-batch and stochastic gradient descent ----------
def gd(A, y, lr, max_epochs, batch=None, target=None, seed=0, trace=None):
rng = random.Random(seed)
n, d = len(A), len(A[0])
w, idx = [0.0] * d, list(range(n))
bs = batch or n # batch=None -> full batch
for ep in range(1, max_epochs + 1):
if batch:
rng.shuffle(idx)
for s in range(0, n, bs):
rows = idx[s:s + bs]
g = [0.0] * d
for i in rows:
e = predict(A[i], w) - y[i] # residual
for j in range(d):
g[j] += e * A[i][j]
k = 2.0 / len(rows) # d/dw of mean squared error
w = [wj - lr * k * gj for wj, gj in zip(w, g)]
try:
loss = mse(A, y, w)
except OverflowError:
loss = math.inf
if trace is not None and ep in trace:
print(f" epoch {ep:>5} loss {loss:.4g}")
if not math.isfinite(loss):
return w, ep, loss, "diverged"
if target is not None and loss <= target:
return w, ep, loss, "reached"
return w, max_epochs, loss, "not reached"
try/except OverflowError は実際にクラッシュしたので入れました。純粋な Python では、発散した重みで1e308近いfloatを2乗すると OverflowError: (34, 'Result too large') が出て止まります。numpy のように黙って inf になってはくれません。エポックは学習データを1周することで、1エポックあたりの更新回数はバッチGDが1回、ミニバッチ32が15回、SGDが480回です。
失敗その1：学習率が大きすぎる
次の節で説明する標準化済みの特徴量に、フルバッチGDを6種類の学習率で回しました。
[learning-rate sweep, scaled features, full batch]
lr 0.01 reached epoch 316 loss 16.69
lr 0.1 reached epoch 29 loss 16.69
lr 0.5 reached epoch 2 loss 16.69
lr 0.9 reached epoch 57 loss 16.69
lr 1.0 diverged epoch 2917 loss inf
lr 1.1 diverged epoch 1193 loss inf
目立つ点が2つあります。大きければ速いわけではなく、0.5は2エポック、0.9は57エポックかかりました。そして0.9と1.0の間に崖があります。trace 引数で発散した回の損失を表示するとこうなります。
[lr 1.1]
epoch 1 loss 7219
epoch 2 loss 1.037e+04
epoch 3 loss 1.507e+04
epoch 5 loss 3.29e+04
epoch 10 loss 2.858e+05
epoch 20 loss 5.169e+07
epoch 30 loss 1.616e+10
[lr 1.0]
epoch 1 loss 5014
epoch 10 loss 8174
epoch 100 loss 9.325e+12
epoch 500 loss 3.354e+54
epoch 1000 loss 2.954e+106
epoch 2000 loss 2.292e+210
epoch 2500 loss 2.019e+262
上の表が隠していたことがあります。表では学習率1.0の発散は「2917エポック目」でしたが、損失は100エポック目ですでに9.3×10¹²でした。inf や NaN だけを見る検査は、何千エポックも遅れて反応します。 実際の学習コードでは「1エポック目の損失を超えたら止める」も入れておくべきです。この検査なら、どちらの回も10エポック以内に止まっていました。
崖はどこにあるのか。MSEの損失はお椀型で、その曲率はヘッセ行列 H = (2/n)XᵀX です。最大固有値 λmax の方向では、1ステップごとに誤差が (1 − η·λmax) 倍になります。この値が −1 を下回ると、誤差は符号を変えながら大きくなります。つまり η > 2/λmax で勾配降下法は発散します。べき乗法で λ を求めました。
# ---------- why a learning rate is "too high" ----------
def curvature(A):
"""largest / smallest eigenvalue of the MSE Hessian (2/n) X^T X (power iteration)"""
n, d = len(A), len(A[0])
H = [[2 / n * sum(r[i] * r[j] for r in A) for j in range(d)] for i in range(d)]
out = []
for step in (lambda v: [predict(row, v) for row in H], lambda v: solve(H, v)):
v = [1.0] * d
for _ in range(300):
u = step(v)
nrm = math.sqrt(sum(x * x for x in u))
v = [x / nrm for x in u]
out.append(predict(v, [predict(row, v) for row in H]))
return out # [lambda_max, lambda_min]
unscaled lambda_max 1.745e+04 lambda_min 0.1563 condition 111,625.8 max stable lr 2/lambda_max = 0.0001146
scaled lambda_max 2.127 lambda_min 1.878 condition 1.1 max stable lr 2/lambda_max = 0.9403
理論と実験がきれいに一致しました。0.9 < 0.9403 は収束し、1.0 > 0.9403 は発散です。0.9が遅かった理由もわかります。1 − 0.9 × 2.127 = −0.91 なので、重みは谷を左右に行き来しながら、誤差を1エポックあたり9%しか縮めません。最適な固定学習率は 2/(λmax + λmin) = 0.499 で、0.5が2エポックで終わったのはこのためです。
失敗その2：特徴量のスケールがばらばら
上の出力の unscaled の行が2つ目の失敗です。面積は最大150、切片の列は常に1なので、お椀はある方向で別の方向の111,626倍も急です。学習率は急な方向に合わせて0.000115未満にする必要があり、その学習率では緩やかな方向がほとんど動きません。限界を少し超えた1.2e-4では、生の特徴量でも発散します。
[unscaled, lr 1.2e-4]
epoch 1 loss 6047
epoch 2 loss 7201
epoch 3 loss 8584
epoch 5 loss 1.222e+04
epoch 10 loss 2.972e+04
epoch 50 loss 3.782e+07
epoch 100 loss 2.883e+11
epoch 200 loss 1.676e+19
epoch 400 loss 5.661e+34
対策は、各特徴量を平均0・標準偏差1にそろえる標準化です。平均と標準偏差は学習データだけで計算し、テストデータにも同じ値を使います。600件全体で計算すると、テストの情報が学習に漏れます（データリーク）。学習後は重みを元の単位に戻します。
# ---------- 4. feature scaling (standardize with TRAIN statistics only) ----------
def fit_scaler(X):
cols = list(zip(*X))
mean = [sum(c) / len(c) for c in cols]
std = [math.sqrt(sum((v - m) ** 2 for v in c) / len(c)) for c, m in zip(cols, mean)]
return mean, std
def scale(X, mean, std):
return [[(v - m) / s for v, m, s in zip(r, mean, std)] for r in X]
def unscale(w, mean, std):
"""weights learned on scaled features -> weights in original units"""
ws = [wj / s for wj, s in zip(w[1:], std)]
return [w[0] - sum(wj * m for wj, m in zip(ws, mean))] + ws
同じ目標までのエポック数です。スケーリングなしは、発散しない限界のすぐ下の1.1e-4で回しました。
[epochs to target: scaled vs unscaled, full batch]
scaled lr 0.5 reached at epoch 2 [4.312, 0.906, -0.394, -5.980]
scaled lr 0.9 reached at epoch 57
epoch 1000 loss 34.73
epoch 10000 loss 18.55
epoch 20000 loss 17.73
epoch 100000 loss 16.74
unscaled lr 1.1e-4 reached at epoch 140538 loss 16.6926 100.7 s [3.610, 0.911, -0.389, -5.913]
前処理1つで2エポック対140,538エポックです。スケーリングなしの回にはもう1つ落とし穴があります。--long を付けずに実行すると20,000エポックで打ち切られ、そのときの損失17.73は最適値より6%高いだけで、ほぼ収束したように見えます。ところが重みは別の話をしています。
unscaled lr 1.1e-4 not reached at epoch 20000 loss 17.7341 14.1 s [0.483, 0.925, -0.357, -5.196]
駅距離の重みは−6.0のはずが−5.20、切片は約4のはずが0.48です。0.1%の目標に届いた140,538エポック後でも、駅距離の重みは−5.913で、閉形式の−6.005とはまだ差がありました。損失がほとんど下がらなくなっていても、緩やかな方向の重みには大きな誤差が残り得ます。「損失が動かなくなったら止める」というやり方では、間違ったモデルのまま止める可能性があります。
ミニバッチとSGD
同じ関数で batch=32 と batch=1、標準化済み、シャッフルの seed は1、学習率をいくつか試しました。
[mini-batch and SGD, scaled, shuffle seed 1]
batch 32 lr 0.02 reached epoch 11 loss 16.6874 test 14.7088 [4.066, 0.908, -0.395, -5.974]
batch 32 lr 0.05 reached epoch 4 loss 16.6845 test 14.6414 [4.045, 0.908, -0.388, -6.005]
batch 32 lr 0.1 reached epoch 9 loss 16.6906 test 14.7560 [3.942, 0.910, -0.399, -5.992]
batch 32 lr 0.2 reached epoch 290 loss 16.6893 test 14.6619 [4.195, 0.907, -0.390, -6.058]
batch 1 lr 0.001 reached epoch 7 loss 16.6910 test 14.7085 [4.080, 0.908, -0.392, -6.007]
batch 1 lr 0.002 reached epoch 4 loss 16.6845 test 14.6064 [4.116, 0.910, -0.391, -6.045]
batch 1 lr 0.005 reached epoch 184 loss 16.6837 test 14.6878 [4.111, 0.909, -0.392, -6.053]
batch 1 lr 0.01 reached epoch 3 loss 16.6805 test 14.6572 [3.868, 0.910, -0.389, -5.972]
batch 1 lr 0.02 not reached epoch 2000 loss 17.7942 test 15.5188 [7.181, 0.882, -0.411, -6.256]
batch 1 lr 0.05 not reached epoch 2000 loss 17.8677 test 15.5023 [7.453, 0.885, -0.399, -6.600]
[SGD lr 0.005, shuffle seeds 1-5]
epochs: [184, 15, 17, 163, 68]
予想と違う結果でしたが、この節でいちばん役に立つのはここです。
- エポック数が学習率の順に並びません。 SGDでは0.002が4エポック、0.005が184エポック、0.01が3エポックでした。1行だけの勾配は全体の勾配のノイズ入りの推定なので、重みは最適点に落ち着かず、その周りで揺れ続けます。「目標に到達」とは、揺れのどれかがたまたま0.1%の幅に入ったという意味です。同じ0.005でも、シャッフルの seed を変えるだけで184、15、17、163、68エポックとばらつきました。SGDの設定を1つの seed で比べてはいけません。
- 学習率が高すぎるSGDは爆発せず、高止まりします。 0.02と0.05では損失は有限のまま17.8付近で止まり、切片は7.2〜7.5までずれました。揺れの幅が目標の幅より大きいので、いつまでも入りません。発散より静かで、見落としやすい失敗です。
- ミニバッチはノイズを平均で打ち消します。 1回32行（1エポック15回更新）なら、0.02〜0.1のどの学習率でも4〜11エポックで届きました。0.2では揺れが大きくなり、290エポックかかりました。
480行程度ではSGDに利点はありません。SGDが活きるのは、データ全体をメモリに載せられないときや、新しい行が届くたびに学習を続けたいときです。実務のSGDは学習率を時間とともに下げて揺れを小さくしますが、ここでは効果を見せるために固定しています。
実行時間とテストデータ
それぞれ5回実行した最速値で、開始から目標の損失に届くまでの時間です。SGDの行は学習率0.005・seed 1、つまり184エポックかかった運の悪いケースです。
[wall time, best of 5]
closed form 0.47 ms
batch GD lr 0.5 1.39 ms
batch GD lr 0.9 39.69 ms
mini-batch 32 3.28 ms
SGD (batch 1) 212.59 ms
test MSE closed 14.6532 batch GD 14.6333
GDの時間には、止めどきを判断するために毎エポック計算している学習データ全体の損失も含みます。この規模の問題では閉形式の圧勝です。4×4行列を1回作って小さな方程式を1回解くだけだからです。テストMSEの差は0.02でノイズの範囲内です。学習の最適値の0.1%手前で止めても悪影響はなく、今回のテストデータではむしろわずかに良い値でした。
| 閉形式 | GD（バッチ／ミニバッチ） | SGD | |
|---|---|---|---|
| 計算量 | XᵀX の作成 O(n·d²) ＋ 解く O(d³) | 1エポック O(n·d) × エポック数 | 1回の更新 O(d) |
| 向いている場面 | 特徴量 d が数千程度まで、データがメモリに載る | d が大きい（d³ が重い）、L1のように閉形式のない正則化 | データが流れてくる、メモリに載らない、オンラインで更新し続ける |
| スケーリング | 不要 | 必須（今回2対140,538エポック） | 必須 |
| ハイパーパラメータ | なし | 学習率、止める条件 | 学習率（と減衰）、バッチサイズ、seed |
| 失敗する場面 | XᵀX が特異（同じ特徴量の重複など）。リッジ回帰か最小二乗ソルバーへ | η > 2/λmax で発散。条件数が悪いと遅い | 学習率が高いと目標の上で揺れ続ける |
正直な限界
- numpy がないので遅いです。 掛け算1回ごとに Python のインタプリタが動きます。480行×重み4つなら問題ありませんが、100万行×100特徴量は現実的ではありません。スケーリングなしの140,538エポックは100.7秒かかりましたが、同じ問題を閉形式は1ミリ秒未満で解きます。ベクトル化されたライブラリなら1エポックがはるかに軽くなります。
- ライブラリは正規方程式では解きません。 XᵀX を作ると問題の条件数が2乗になります。float64 で条件数10⁵程度なら害はありませんが、ほぼ共線な特徴量では桁落ちします。scikit-learn なら代わりに
LinearRegression（XᵀX を作らず行列分解で解く最小二乗ソルバー）、スケーリングにStandardScaler、SGD版にSGDRegressorを使うところです。このシリーズの目的は仕組みを見ることなので、名前を挙げるだけで実行はしていません。 - 合成データ1種類だけです。 特徴量は独立で一様分布、ノイズはガウス分布なので、問題としてはとても素直です。相関のある特徴量は標準化しても条件数を押し上げますし、外れ値は二乗誤差を大きく狂わせます。そもそも線形でないものは、ここにあるどの方法でも直せません。
- 学習率の探索は小規模です。 試した学習率は数種類、seed は1〜5個で、学習率の減衰も使っていません。実際の学習コードではスケジュールを組み、損失が上がったら早めに止めます。
自分で動かす
ここまでのコードを上から順につなげると linreg.py になります。最後に次のドライバーを足してください。python3 linreg.py は手元で27秒、python3 linreg.py --long はスケーリングなしの回を目標まで走らせるので、さらに100秒ほどかかります。
def fmt(w):
return "[" + ", ".join(f"{v:.3f}" for v in w) + "]"
def best_of(fn, k=5):
ts = []
for _ in range(k):
t = time.perf_counter(); fn(); ts.append(time.perf_counter() - t)
return min(ts)
if __name__ == "__main__":
LONG = "--long" in sys.argv # unscaled run until it converges (minutes)
X, y = make_data()
cut = int(len(X) * 0.8) # rows are already random: 480 / 120
Xtr, ytr, Xte, yte = X[:cut], y[:cut], X[cut:], y[cut:]
Atr, Ate = add_bias(Xtr), add_bias(Xte)
print(f"Python {sys.version.split()[0]} train {len(Xtr)} test {len(Xte)}")
print("true ", fmt([TRUE_B] + TRUE_W))
w_cf = closed_form(Atr, ytr)
best = mse(Atr, ytr, w_cf)
print(f"closed {fmt(w_cf)} train {best:.4f} test {mse(Ate, yte, w_cf):.4f}")
w_true = [TRUE_B] + TRUE_W # noise only: the floor for any model
print(f"true w train {mse(Atr, ytr, w_true):.4f} test {mse(Ate, yte, w_true):.4f}")
target = best * 1.001 # within 0.1% of the optimum
print(f"target: train MSE <= {target:.4f}")
mean, std = fit_scaler(Xtr)
Str, Ste = add_bias(scale(Xtr, mean, std)), add_bias(scale(Xte, mean, std))
for name, A in (("unscaled", Atr), ("scaled", Str)):
hi, lo = curvature(A)
print(f"{name:<9} lambda_max {hi:.4g} lambda_min {lo:.4g} "
f"condition {hi/lo:,.1f} max stable lr 2/lambda_max = {2/hi:.4g}")
print("\n[learning-rate sweep, scaled features, full batch]")
for lr in (0.01, 0.1, 0.5, 0.9, 1.0, 1.1):
w, ep, loss, st = gd(Str, ytr, lr, 5000, target=target)
print(f"lr {lr:<5} {st:<9} epoch {ep:>5} loss {loss:.4g}")
print("[lr 1.1]")
gd(Str, ytr, 1.1, 30, trace={1, 2, 3, 5, 10, 20, 30})
print("[lr 1.0]")
gd(Str, ytr, 1.0, 5000, trace={1, 10, 100, 500, 1000, 2000, 2500})
print("[unscaled, lr 1.2e-4]")
gd(Atr, ytr, 1.2e-4, 400, trace={1, 2, 3, 5, 10, 50, 100, 200, 400})
print("\n[epochs to target: scaled vs unscaled, full batch]")
w_s, ep, loss, st = gd(Str, ytr, 0.5, 5000, target=target)
print(f"scaled lr 0.5 {st} at epoch {ep} {fmt(unscale(w_s, mean, std))}")
w_s9, ep, loss, st = gd(Str, ytr, 0.9, 5000, target=target)
print(f"scaled lr 0.9 {st} at epoch {ep}")
t = time.perf_counter()
w_u, ep, loss, st = gd(Atr, ytr, 1.1e-4, 600000 if LONG else 20000, target=target,
trace={1000, 10000, 20000, 100000, 200000, 300000})
print(f"unscaled lr 1.1e-4 {st} at epoch {ep} loss {loss:.4f} "
f"{time.perf_counter()-t:.1f} s {fmt(w_u)}")
print("\n[mini-batch and SGD, scaled, shuffle seed 1]")
for bs, lrs in ((32, (0.02, 0.05, 0.1, 0.2)), (1, (0.001, 0.002, 0.005, 0.01, 0.02, 0.05))):
for lr in lrs:
w, ep, loss, st = gd(Str, ytr, lr, 2000, batch=bs, target=target, seed=1)
print(f"batch {bs:<3} lr {lr:<6} {st:<11} epoch {ep:>4} loss {loss:.4f} "
f"test {mse(Ste, yte, w):.4f} {fmt(unscale(w, mean, std))}")
print("[SGD lr 0.005, shuffle seeds 1-5]")
print(" epochs:", [gd(Str, ytr, 0.005, 2000, batch=1, target=target, seed=sd)[1]
for sd in range(1, 6)])
print("\n[wall time, best of 5]")
print(f"closed form {best_of(lambda: closed_form(Atr, ytr))*1000:8.2f} ms")
for lbl, lr, kw in (("batch GD lr 0.5", 0.5, {}), ("batch GD lr 0.9", 0.9, {}),
("mini-batch 32", 0.05, {"batch": 32, "seed": 1}),
("SGD (batch 1)", 0.005, {"batch": 1, "seed": 1})):
ms = best_of(lambda: gd(Str, ytr, lr, 5000, target=target, **kw)) * 1000
print(f"{lbl:<21} {ms:8.2f} ms")
print(f"\ntest MSE closed {mse(Ate, yte, w_cf):.4f} batch GD {mse(Ste, yte, w_s):.4f}")
損失・エポック数・重みは seed が同じなら毎回同じになります。実行時間はマシンによって変わりますが、比率は大きくは変わらないはずです。
次回： k-meansクラスタリングをゼロから。「勾配を求めて進む」の代わりに「割り当ててから更新する」ループです。同じように実測で比べる記事としては、C言語のソートアルゴリズムの実測もどうぞ。
よくある質問
学習率を大きくすると、なぜ勾配降下法は発散するのですか？
平均二乗誤差の損失はお椀型で、その一番急な方向の曲率は (2/n)XᵀX の最大固有値 λmax です。学習率が 2/λmax を超えると、毎ステップ谷を飛び越して誤差が大きくなります。今回の標準化データでは 2/λmax = 0.94 で、0.9は収束し、1.0は損失が100エポックで5,014から9.3×10¹²まで増えました。
特徴量のスケーリングで、勾配降下法はどれくらい速くなりますか？
同じ480行で、標準化した特徴量は学習率0.5で2エポックで目標の損失に届きました。生の特徴量は発散しない最大級の学習率でも140,538エポック必要でした。方向ごとの曲率が111,626倍違っていたためです。
線形回帰は閉形式と勾配降下法のどちらを使うべきですか？
特徴量が少なくデータがメモリに載るなら閉形式（最小二乗ソルバー）です。今回は0.47ミリ秒で、一番速く調整した勾配降下法の1.39ミリ秒より速く終わりました。特徴量が非常に多い、データが流れてくる、閉形式のない正則化を加える、といった場合は勾配降下法やSGDを使います。