手計算で学ぶ誤差逆伝播法、そしてautogradへ

本レッスンは、loss.backward()でネットワークを学習したことがあり、その1行が実際に何をしているのか知りたい人のためのものです。鉛筆で扱える大きさの数値で1つの入力を2層ネットワークに通し、連鎖律だけを使って損失をすべての重みで微分し、同じ計算をnumpyで書き、有限差分でそれが正しいことを証明し、それから同一のネットワークをPyTorchのautogradに通します。勾配はdoubleの最後の数ビットまで一致し、読み終えるころには、なぜ一致しなければならないのかを理解しているはずです。

学習ループそのものは損失と勾配に関する入門レッスンで、小さなネットワークの形は最初のニューラルネットワークで扱いました。ここでは勾配だけに焦点を当てます。

ネットワークと数値

ネットワークは、入力2つ、tanh非線形性を持つ隠れユニット2つ、線形出力1つです。損失はスカラーの目標値に対する二乗誤差の半分です。半分にしているのは導関数をきれいにするためです。

順伝播の計算を4行で書くと次のようになります。

z1 = W1 x + b1          W1 is 2 by 2, b1 has 2 entries
h  = tanh(z1)
ŷ  = W2 h + b2          W2 has 2 entries, b2 is a scalar
L  = ½ (ŷ − y)²

具体的な値は、最初のステップが切りのよい数になるように選んであります。

記号	値
x	[0.5, −1.0]
y	1.0
W1	行 [0.2, −0.4] と [0.6, 0.1]
b1	[0.1, −0.3]
W2	[0.8, −0.5]
b2	0.2

学習可能な数値は全部で9つです。W1に4つ、b1に2つ、W2に2つ、b2に1つ。誤差逆伝播法の仕事は、パラメータごとに1つ、9つの偏微分を生み出すことです。

順伝播

上の値で4行を順に計算します。

1つ目の隠れユニットが見るのは0.2 × 0.5 + (−0.4) × (−1.0) + 0.1で、これは0.1 + 0.4 + 0.1なので、このユニットのz1は0.6です。2つ目が見るのは0.6 × 0.5 + 0.1 × (−1.0) − 0.3で、これは0.3 − 0.1 − 0.3なので−0.1です。活性化前のベクトルはz1 = [0.6, −0.1]です。

tanhを適用するとh = [0.5370, −0.0997]になります。出力は0.8 × 0.5370 + (−0.5) × (−0.0997) + 0.2 = 0.4296 + 0.0498 + 0.2 = 0.6795です。

誤差はŷ − y = −0.3205で、損失はその二乗の半分、0.0514です。この誤差をeとして書き留めてください。逆向きに進むときに最初に必要になるもので、すべての勾配に現れます。

連鎖律を逆向きに

連鎖律が言うのは、Lがŷに依存し、ŷがW2に依存するなら、LのW2に関する微分は、Lのŷに関する微分にŷのW2に関する微分を掛けたものだ、ということです。誤差逆伝播法とは、この規則を演算ごとに1回ずつ、損失から始めて入力に向かって適用し、途中の微分を再計算せずにすべて再利用すること以上のものではありません。

出力層

½ (ŷ − y)²のŷに関する微分はŷ − y、つまりe = −0.3205です。この1つの数値が、他のすべての種になります。

出力ŷ = W2 h + b2はW2について線形なので、ŷのW2の各要素に関する微分は、hの対応する要素です。種を掛けます。

dL/dW2 = e · h  = −0.3205 × [0.5370, −0.0997] = [−0.1721, 0.0319]
dL/db2 = e      = −0.3205

9つの勾配のうち2つが終わりました。出力はhについても線形で、ŷのhに関する微分はW2なので、隠れ層に手渡さなければならない勾配は次のとおりです。

dL/dh = e · W2 = −0.3205 × [0.8, −0.5] = [−0.2564, 0.1603]

隠れ層

次はh = tanh(z1)で、要素ごとに適用されます。tanh(z)の微分は1 − tanh(z)²で、tanh(z1)はすでに手元にあるhです。つまり各隠れユニットの局所微分は1 − h²です。2つのユニットについてそれは1 − 0.2884 = 0.7116と1 − 0.0099 = 0.9901です。要素ごとに掛けます。

dL/dz1 = dL/dh ⊙ (1 − h²) = [−0.2564 × 0.7116, 0.1603 × 0.9901] = [−0.1825, 0.1587]

最後に、z1 = W1 x + b1はW1とb1について線形で、出力層がW2とb2について線形だったのとまったく同じです。W1の各行は1つの隠れユニットに入力するので、行iの勾配はユニットiのdL/dz1に入力ベクトルxを掛けたものです。これは外積です。

dL/dW1 = dL/dz1 ⊗ x = rows [−0.1825 × 0.5, −0.1825 × −1.0] and [0.1587 × 0.5, 0.1587 × −1.0]
       = rows [−0.0912, 0.1825] and [0.0793, −0.1587]
dL/db1 = dL/dz1 = [−0.1825, 0.1587]

9つの数値がすべて揃いました。パターンに注目してください。すべての線形層は、入ってくる勾配にその入力を掛けたものを勾配として寄与し、すべての非線形性は要素ごとの係数を寄与し、ある層に入ってくる勾配は常に、前の層に入ってきた勾配にその層の重みを掛けたものです。何も2回微分していません。

同じことをnumpyで

下のコードは手計算を1行ずつそのまま書いたものです。逆伝播の部分は順伝播の逆順に実行され、各行はすでに手元にある値だけを使います。

import numpy as np

np.set_printoptions(precision=4, suppress=True)

x = np.array([0.5, -1.0])
y = 1.0
W1 = np.array([[0.2, -0.4],
               [0.6, 0.1]])
b1 = np.array([0.1, -0.3])
W2 = np.array([0.8, -0.5])
b2 = 0.2

# forward
z1 = W1 @ x + b1
h = np.tanh(z1)
y_hat = W2 @ h + b2
loss = 0.5 * (y_hat - y) ** 2

# backward, in reverse order of the forward pass
d_y_hat = y_hat - y
d_W2 = d_y_hat * h
d_b2 = d_y_hat
d_h = d_y_hat * W2
d_z1 = d_h * (1.0 - h ** 2)
d_W1 = np.outer(d_z1, x)
d_b1 = d_z1

print("z1    ", z1)
print("h     ", h)
print("y_hat ", round(float(y_hat), 4))
print("loss  ", round(float(loss), 4))
print("dW2   ", d_W2)
print("db2   ", round(float(d_b2), 4))
print("dz1   ", d_z1)
print("dW1   ", d_W1)
print("db1   ", d_b1)

表示される値は手計算と一致します（numpyが列を揃えるので、空白は異なることがあります）。

z1     [ 0.6 -0.1]
h      [ 0.537  -0.0997]
y_hat  0.6795
loss   0.0514
dW2    [-0.1721  0.0319]
db2    -0.3205
dz1    [-0.1825  0.1587]
dW1    [[-0.0912  0.1825]
 [ 0.0793 -0.1587]]
db1    [-0.1825  0.1587]

これをバッチに拡張するとき、バグの定番の原因は形状です。規則はテンソル形状とブロードキャストのリファレンスで扱っており、本レッスン後半のバッチ版では、外積がどのように行列積になるかを示します。

有限差分による勾配チェック

手計算の導出は符号や転置を間違えていても、コードは動いてしまいます。唯一の独立した検査は、各パラメータをごくわずかに動かし、損失がどれだけ変わるかを測り、その傾きを解析的な勾配と比較することです。中心差分(L(p + ε) − L(p − ε)) / 2εはε²のオーダーまで正確なので、倍精度でε = 1e-5なら、おおよそ小数点以下10桁までの一致が期待できます。

import numpy as np

x = np.array([0.5, -1.0])
y = 1.0
W1 = np.array([[0.2, -0.4], [0.6, 0.1]])
b1 = np.array([0.1, -0.3])
W2 = np.array([0.8, -0.5])
b2 = np.array([0.2])


def forward_backward(W1, b1, W2, b2):
    z1 = W1 @ x + b1
    h = np.tanh(z1)
    y_hat = W2 @ h + b2[0]
    loss = 0.5 * (y_hat - y) ** 2
    d_y_hat = y_hat - y
    grads = {"W2": d_y_hat * h, "b2": np.array([d_y_hat])}
    d_z1 = d_y_hat * W2 * (1.0 - h ** 2)
    grads["W1"] = np.outer(d_z1, x)
    grads["b1"] = d_z1
    return loss, grads


params = {"W1": W1, "b1": b1, "W2": W2, "b2": b2}
loss, grads = forward_backward(**params)

eps = 1e-5
worst = 0.0
for name, p in params.items():
    numeric = np.zeros_like(p)
    for idx in np.ndindex(p.shape):
        saved = p[idx]
        p[idx] = saved + eps
        loss_plus, _ = forward_backward(**params)
        p[idx] = saved - eps
        loss_minus, _ = forward_backward(**params)
        p[idx] = saved
        numeric[idx] = (loss_plus - loss_minus) / (2 * eps)
    rel = np.abs(numeric - grads[name]) / (np.abs(numeric) + np.abs(grads[name]) + 1e-12)
    worst = max(worst, float(rel.max()))
    print(f"{name}: analytic {grads[name].ravel().round(6)}  numeric {numeric.ravel().round(6)}")

print("worst relative error:", f"{worst:.1e}")

各行で解析的勾配と数値勾配が小数点以下6桁まで一致し、最悪の相対誤差は1e-10のオーダーの数値として表示されます。隠れ層の符号誤りならこれが1のオーダーに、(1 − h²)の係数の欠落なら0.1のオーダーに跳ね上がるので、この検査は微妙なものではありません。

同じネットワークをPyTorchのautogradで

次に、同一のネットワークをテンソルで構築し、PyTorchに微分させます。唯一の新しい要素は、4つのパラメータに付けるrequires_grad=Trueで、これがautogradにそれらに触れるすべての演算を記録するよう指示します。

import numpy as np
import torch

torch.set_printoptions(precision=4, sci_mode=False)

x = torch.tensor([0.5, -1.0], dtype=torch.float64)
y = torch.tensor(1.0, dtype=torch.float64)
W1 = torch.tensor([[0.2, -0.4], [0.6, 0.1]], dtype=torch.float64, requires_grad=True)
b1 = torch.tensor([0.1, -0.3], dtype=torch.float64, requires_grad=True)
W2 = torch.tensor([0.8, -0.5], dtype=torch.float64, requires_grad=True)
b2 = torch.tensor(0.2, dtype=torch.float64, requires_grad=True)

z1 = W1 @ x + b1
h = torch.tanh(z1)
y_hat = W2 @ h + b2
loss = 0.5 * (y_hat - y) ** 2

print("loss           ", round(loss.item(), 4))
print("grad_fn of loss", type(loss.grad_fn).__name__)
print("grad_fn of h   ", type(h.grad_fn).__name__)

loss.backward()

print("dW1", W1.grad)
print("db1", b1.grad)
print("dW2", W2.grad)
print("db2", b2.grad)

# the numpy backward pass from earlier, for a direct comparison
d_y_hat = y_hat.item() - 1.0
h_np = h.detach().numpy()
d_z1 = d_y_hat * np.array([0.8, -0.5]) * (1.0 - h_np ** 2)
print("W1 gradients match numpy:", np.allclose(W1.grad.numpy(), np.outer(d_z1, [0.5, -1.0]), atol=1e-12))

損失は0.0514と表示され、2つのgrad_fnの名前はMulBackward0とTanhBackward0と表示され、4つの.gradテンソルはnumpyのブロックと同じ9つの数値を保持し、最後の行はTrueと表示されます。

backward()を呼ぶ前から、すべての中間テンソルはすでにgrad_fnを持っています。この属性が、順伝播の実行中にautogradが残した記録です。損失は乗算で作られたのでMulBackward0、hはtanhで作られたのでTanhBackward0です。それぞれのオブジェクトは自身の入力を知っており、自身の局所微分の計算方法を知っています。

なぜ数値が一致するのか

autogradは手計算より賢いことをしているわけではありません。まさに手計算そのものを、帳簿付けを自動化して行っているのです。

順伝播の間、勾配を必要とするテンソルに対するすべての演算がグラフにノードを追加します。ノードは逆伝播の式が必要とするものを保存します。tanhなら微分が1 − h²なので出力hを保存します。行列積ならもう一方のオペランドを保存します。これは、順方向に進みながら紙に書き留めたのと同じ情報です。

backward()が実行されると、損失から種となる勾配1で始め、グラフを逆トポロジカル順にたどります。各ノードでそのノードの逆伝播関数を呼び、それは上から届いた勾配を受け取って各入力に関する勾配を返します。tanhノードはdL/dhを受け取り、dL/dh ⊙ (1 − h²)を返します。行列積ノードはdL/dz1を受け取り、重みの入力に対してはxとの外積を、ベクトルの入力に対してはW1との積を返します。勾配がrequires_grad=Trueの葉テンソルに到達すると、そのテンソルの.gradフィールドに加算されます。

この説明から3つの帰結が導かれ、3つとも実務で重要です。

• リバースモードである。 パラメータがいくつあろうと、作業量は1回の逆向き掃引に比例します。掃引は勾配の単一のベクトルを運び、それに局所ヤコビアンを掛けていくからです。完全なヤコビアンを計算すると出力ごとに1回の掃引が必要になり、だからこそフォワードモードの微分は数百万の重みに対するスカラー損失には不向きな道具です。
• 勾配は累積する。 .gradフィールドは上書きではなく加算されます。1つのパラメータが複数のノードに入力することがあり、それらの寄与を合計しなければならないからです。学習ループが各ステップの前にoptimizer.zero_grad()を呼ぶかgradをNoneにするのもそのためです。忘れると今回のバッチの勾配が前回のバッチの勾配に加算され、損失曲線が不可解にでこぼこになります。標準的な順序はPyTorch学習ループのチートシートにあります。
• グラフは使用後に解放される。 デフォルトではbackward()は保存されたテンソルを解放するので、同じ損失に対して2回呼ぶと、retain_graph=Trueを渡さない限りエラーになります。順伝播のたびにグラフを作り直すからこそ、ifやforのような制御フローが自然に機能します。グラフは、コードが実際に実行したものそのものです。

上のテンソルはfloat64で作ったので、2つの実装は約15桁の有効数字まで一致します。デフォルトのfloat32でも一致しますが約6桁までで、それはどちらかの方法の誤りではなく、想定どおりの丸め誤差です。

2つを並べて学習する

手計算の例では入力は1つでした。実際の学習ではバッチを使い、逆伝播の唯一の変更点は、外積が行列積になり、バイアスの勾配がバッチにわたって合計されることです。下のブロックは、同じ64個の合成データ点に対して、同じ初期重みから、素朴な勾配降下法でnumpyのネットワークとautogradのコピーを学習します。numpyの逆伝播が正しければ、2つの損失の列は同一になります。

import numpy as np
import torch

rng = np.random.default_rng(0)
X = rng.uniform(-1.0, 1.0, size=(64, 2))
Y = X[:, 0] * X[:, 1] + 0.05 * rng.normal(size=64)

hidden = 8
W1 = rng.normal(scale=0.5, size=(hidden, 2))
b1 = np.zeros(hidden)
W2 = rng.normal(scale=0.5, size=hidden)
b2 = 0.0
lr = 0.1

# torch copies of the same starting point, in float64 so the arithmetic agrees
tW1 = torch.tensor(W1, requires_grad=True)
tb1 = torch.tensor(b1, requires_grad=True)
tW2 = torch.tensor(W2, requires_grad=True)
tb2 = torch.tensor(b2, dtype=torch.float64, requires_grad=True)
tX = torch.tensor(X)
tY = torch.tensor(Y)

for step in range(1001):
    # numpy: forward, backward, update
    Z1 = X @ W1.T + b1                      # (64, hidden)
    H = np.tanh(Z1)
    y_hat = H @ W2 + b2                     # (64,)
    loss = np.mean((y_hat - Y) ** 2)
    d_y_hat = 2.0 * (y_hat - Y) / len(Y)    # derivative of the mean squared error
    dW2 = H.T @ d_y_hat
    db2 = d_y_hat.sum()
    dZ1 = np.outer(d_y_hat, W2) * (1.0 - H ** 2)
    dW1 = dZ1.T @ X
    db1 = dZ1.sum(axis=0)
    W1 -= lr * dW1
    b1 -= lr * db1
    W2 -= lr * dW2
    b2 -= lr * db2

    # torch: same forward, gradients from autograd, same update
    t_loss = torch.mean((torch.tanh(tX @ tW1.T + tb1) @ tW2 + tb2 - tY) ** 2)
    t_loss.backward()
    with torch.no_grad():
        for p in (tW1, tb1, tW2, tb2):
            p -= lr * p.grad
            p.grad = None

    if step % 200 == 0:
        print(f"step {step:4d}  numpy {loss:.6f}  torch {t_loss.item():.6f}")

print("max |W1 difference|:", f"{np.abs(W1 - tW1.detach().numpy()).max():.1e}")

2つの列はどの行でも同じ小数点以下6桁の損失を表示し、損失は実行を通じて着実に下がり、最終的な重みの差はdoubleの最後のビットのノイズで、1e-10をはるかに下回ります。numpyのコードにはgrad_fnもグラフもbackward()もありませんが、同じ連鎖律を手で符号化しているので、同じ重みに到達します。

d_y_hatを計算する行を見てください。単一の例のバージョンと異なるのは、損失が平均になったので係数2/Nが現れる点だけです。それ以外のすべての行は、先ほどの行にバッチ次元を加えたものです。深層学習ハブが「層の逆伝播は2つの行列積で定義される」と言うのは、まさにこのことです。PyTorchでは、この積こそがnn.Linearが逆伝播関数として登録するものであり、50層のモデルは本レッスンの2層の計算を50回適用したものです。

演習

numpy版とPyTorch版の両方で、隠れ層の非線形性をシグモイド関数σ(z) = 1 / (1 + exp(−z))に置き換えてください。出力を使った局所微分を導出し（σ(z) (1 − σ(z))になります）、逆伝播を更新し、有限差分チェックで最悪の相対誤差が依然として1e-10のオーダーであることを確認してください。次に、損失を目標値1.0に対する二値交差エントロピーに切り替えて繰り返してください。最後に、すべてのテンソルをfloat32、ε = 1e-5で勾配チェックを試し、なぜもう通らないのかを表示された数値から説明してください。

まとめ

• 誤差逆伝播法は、損失から入力に向かって演算ごとに1回ずつ連鎖律を適用し、各中間勾配を再計算せずに再利用するものです。
• 線形層の重みの勾配は、入ってくる勾配にその層の入力を掛けたもので、後ろに渡す勾配は、入ってくる勾配にその重みを掛けたものです。非線形性は要素ごとの係数を寄与します。
• float64での中心有限差分は、解析的勾配を小数点以下約10桁まで検査し、符号や転置の誤りを即座に捕まえます。
• autogradは順伝播の各演算を、その逆伝播の式が必要とするものを保存するノードとして記録し、グラフを逆向きにたどって.gradに累積します。
• numpyとPyTorchの勾配が一致するのは、同じ算術を同じ順序で実行しているからで、差があるとすれば浮動小数点の丸め誤差です。

さらに読む

• D. E. Rumelhart, G. E. Hinton and R. J. Williams, “Learning representations by back-propagating errors”, Nature, 1986: https://www.nature.com/articles/323533a0
• A. G. Baydin, B. A. Pearlmutter, A. A. Radul and J. M. Siskind, “Automatic differentiation in machine learning: a survey”, JMLR, 2018: https://jmlr.org/papers/v18/17-468.html
• PyTorch documentation, “Autograd mechanics”: https://pytorch.org/docs/stable/notes/autograd.html
• I. Goodfellow, Y. Bengio and A. Courville, Deep Learning, chapter 6 on back-propagation: https://www.deeplearningbook.org/contents/mlp.html