線形回帰とロジスティック回帰をゼロから実装する

レッスン01では、モデルを空欄のあるプログラム、学習をその空欄を埋める数値の探索として説明しました。本レッスンでは、他のあらゆるモデルの比較対象となる2つのモデルについて、その探索を手作業で行います。線形回帰は数値を予測し、厳密な公式を持ちます。ロジスティック回帰は確率を予測し、公式を持たないので、勾配降下法を導入します。これはレッスン04が一般化し、あらゆるニューラルネットワークが使う探索です。各実装は数行のnumpyで、scikit-learnと照合するので、ライブラリが同じ計算をしていることが確認できます。

線形回帰: モデルと損失

モデルは重み付き和です。特徴量がxの1つなら、予測はw0 + w1 * x、つまり切片と傾きです。特徴量がd個なら、切片に特徴量ごとの重みを1つずつ加えたものです。パラメータは重みであり、学習がそれを選びます。

損失は平均二乗誤差、すなわち予測とラベルの差の2乗を事例について平均したものです。2乗することで、すべての誤差が正になり、大きな誤差が小さな誤差よりはるかに強く罰せられ、しかも最小値を探索ではなく代数で求められる損失になります。

ごく小さな例を頭の中で解いてみます。xが1、2、3で、ラベルが2、3、5の3事例です。xの平均は2、ラベルの平均は10÷3で約3.333です。特徴量が1つのとき、最良の傾きは、（xから平均を引いたもの）×（ラベルから平均を引いたもの）の和を、（xから平均を引いたもの）の2乗の和で割ったものです。分子は（−1×−1.333）＋0＋（1×1.667）で3.0。分母は1＋0＋1で2。傾きは1.5です。切片はラベルの平均から傾き×xの平均を引いたもの、3.333−3.0で0.333です。この直線は1.833、3.333、4.833を予測し、ラベルとのずれは0.167、0.333、0.167です。二乗損失でこれより良い直線はありません。この公式は、下の一般解の特徴量1つの場合です。

閉形式解

事例を、1事例1行、1特徴量1列の行列Xに積み上げ、切片を他の重みと同じように扱えるよう1の列を加えます。ラベルはベクトルyに積み上げます。すると二乗損失は重みベクトルwのお椀型の関数になり、お椀の底は、あらゆる方向で傾きがゼロになる場所です。勾配をゼロと置くと正規方程式が得られます。Xの転置×X×w＝Xの転置×y。これは小さな正方の線形連立方程式で、np.linalg.solveがwを直接求めます。

import numpy as np

rng = np.random.default_rng(0)
n = 50
x = rng.uniform(0, 10, size=n)              # hours of study, say
noise = rng.normal(0, 1.5, size=n)
y = 3.0 + 2.0 * x + noise                   # true line: intercept 3, slope 2

# Design matrix: a column of ones for the intercept, then x.
X = np.column_stack([np.ones(n), x])

# Normal equations: w = (X^T X)^-1 X^T y
w = np.linalg.solve(X.T @ X, X.T @ y)
print("closed-form intercept, slope:", np.round(w, 3))

residuals = y - X @ w
print("mean squared error:", round(float(np.mean(residuals ** 2)), 3))

closed-form intercept, slope: [2.636 2.079]
mean squared error: 2.235

データは切片3、傾き2に標準偏差1.5のノイズを加えて生成したので、当てはめた2.636と2.079は近いものの厳密ではありません。ノイズ入りの50点では真実を特定しきれず、平均二乗誤差2.235はノイズの分散（1.5の2乗は2.25）とほぼ同じです。これが、ノイズのあるデータにうまく当てはまったモデルの姿です。残差はノイズそのもので、それ以上取り除けるものはありません。

明示的に逆行列を計算するのではなくnp.linalg.solveを使うのは、身につける価値のある習慣です。連立方程式を解くほうが、行列を反転して掛けるより速く、数値的にも安全です。特徴量が非常に多い場合や、互いにほぼ複製の特徴量がある場合、行列は特異に近くなり、solveでも苦しくなります。ライブラリは代わりに分解を使い、重みに小さなペナルティを加えます。これが「リッジ回帰」の意味です。古典的機械学習のハブでその系統を扱います。

scikit-learnとの照合

import numpy as np
from sklearn.linear_model import LinearRegression

rng = np.random.default_rng(0)
n = 50
x = rng.uniform(0, 10, size=n)
y = 3.0 + 2.0 * x + rng.normal(0, 1.5, size=n)
X = np.column_stack([np.ones(n), x])
w = np.linalg.solve(X.T @ X, X.T @ y)

sk = LinearRegression().fit(x.reshape(-1, 1), y)
print("numpy  :", np.round(w, 6))
print("sklearn:", np.round([sk.intercept_, sk.coef_[0]], 6))
print("max difference:", float(np.max(np.abs(w - [sk.intercept_, sk.coef_[0]]))))

numpy  : [2.636207 2.078504]
sklearn: [2.636207 2.078504]
max difference: 1.9539925233402755e-14

両者は小数点以下14桁まで一致し、差は浮動小数点の丸め誤差です。LinearRegressionは切片の列を自動で加えるので、Xではなくxを列として渡しています。

分類: 二乗誤差が破綻する理由

次はラベルが「はい・いいえ」で、レッスン02の離脱の表のように1か0で符号化されています。0と1に直線を当てはめ、0.5より上を「はい」と呼ぶこともできます。これは2つの理由でうまくいきません。直線の出力には上下限がないので、極端な入力に対しては3.7や−2を予測し、答えとして意味をなさない数値になります。そして二乗誤差は、自信のある正解を罰します。ラベルが1で直線が4を出せば損失は9になり、学習はすでに正しく当てている事例の方へ直線を引き下げてしまいます。

修正は2つの部分からなります。重み付き和は、特徴量を組み合わせる最も安価な方法なので残します。ただしその出力を、任意の実数を0と1の間の区間に写す押しつぶし関数に通します。そして結果をラベルが1である確率として読み、目標の数値に近いことではなく、確率として正しいことに報いる損失を選びます。

シグモイド

押しつぶし関数はシグモイド、すなわち1÷（1＋eの−z乗）です。zが0のとき0.5を与えます。zが2なら約0.88、4なら約0.98、−2なら約0.12です。大きな正のスコアは1に近い確率に、大きな負のスコアは0に近い確率になり、遷移は滑らかです。探索がたどるべき傾きを必要とするので、この滑らかさは重要です。レッスン02で離脱ラベルを生成した隠れたルールは、まさにこれでした。重み付き和をシグモイドに通したものです。

対数損失

ラベルが1で予測確率がpの事例1つについて、損失はpの自然対数にマイナスを付けたものです。pが0.9なら損失は約0.105。0.5なら0.693。0.1なら2.303。損失は小さな間違いにはゆっくりと、自信のある誤答には急激に増えます。ラベルが0なら同じことが（1−p）に当てはまります。全事例の平均が対数損失で、二値交差エントロピーとも呼ばれます。これは予測されたpを較正された確率として振る舞わせる損失であり、あらゆる言語モデルの次トークン予測の背後にある損失を、数千の選択肢からなる語彙に拡大したものでもあります。

勾配

シグモイドと対数損失を組み合わせると、きれいなことが起こります。平均損失の重みベクトルに関する微分は、Xの転置×（p−y）をnで割ったものです。こう読んでください。各事例について、誤差は予測確率から真のラベルを引いた−1と1の間の数値であり、各重みの勾配は、その誤差と対応する特徴量の値との積の平均です。モデルが間違えるときに大きな値を取る特徴量に付いた重みが、最も大きく押されます。pはシグモイドを通してwに依存するので閉形式はなく、「勾配＝ゼロ」の方程式は代数的に解けません。代わりに、勾配と逆の方向へ繰り返し進みます。

import numpy as np

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

def log_loss(y, p):
    eps = 1e-12
    p = np.clip(p, eps, 1 - eps)
    return float(-np.mean(y * np.log(p) + (1 - y) * np.log(1 - p)))

rng = np.random.default_rng(0)
n = 200
# Two features. Class 1 tends to have larger x1 + x2.
X_raw = rng.normal(0, 1, size=(n, 2))
true_w = np.array([-0.5, 2.0, -1.5])          # bias, w1, w2
logits = true_w[0] + X_raw @ true_w[1:]
y = (rng.random(n) < sigmoid(logits)).astype(float)

X = np.column_stack([np.ones(n), X_raw])      # add the bias column
w = np.zeros(3)
lr = 0.5
for step in range(1, 501):
    p = sigmoid(X @ w)
    grad = X.T @ (p - y) / n                  # dL/dw for mean log loss
    w -= lr * grad
    if step in (1, 10, 50, 100, 500):
        print(f"step {step:4d}  loss {log_loss(y, p):.4f}  w {np.round(w, 3)}")

acc = float(np.mean((sigmoid(X @ w) >= 0.5) == y))
print("training accuracy:", round(acc, 3))

step    1  loss 0.6931  w [-0.04   0.134 -0.086]
step   10  loss 0.4789  w [-0.21   0.851 -0.555]
step   50  loss 0.3940  w [-0.369  1.765 -1.227]
step  100  loss 0.3861  w [-0.435  2.089 -1.485]
step  500  loss 0.3846  w [-0.488  2.319 -1.67 ]
training accuracy: 0.84

最初の損失は0.6931で、これは2の自然対数です。すべての重みがゼロのとき、モデルはすべてに0.5と答え、意見を持たないときの対数損失のコストがこれです。ステップ10までに重みは正しい符号になり、損失は3分の1下がっています。ステップ100までに、真の値−0.5、2.0、−1.5に近づいています。ステップ100から500の間、損失はほとんど動かない一方で重みは少しずつ大きくなり続けます。有限のノイズ入りデータでは、当てはめた重みが生成に使った重みをわずかに行き過ぎるからです。精度0.84はこのデータの上限近くです。ラベルは確率から引いたので、自分の確率と反対側に落ちた事例があり、どんな重みでも取り戻せません。

log_loss内のクリッピングは、ちょうど0の対数を取るのを防ぎます。numpyはそれを負の無限大として報告します。ライブラリは同じことを黙って行っています。

scikit-learnとの照合

import numpy as np
from sklearn.linear_model import LogisticRegression

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

rng = np.random.default_rng(0)
n = 200
X_raw = rng.normal(0, 1, size=(n, 2))
true_w = np.array([-0.5, 2.0, -1.5])
y = (rng.random(n) < sigmoid(true_w[0] + X_raw @ true_w[1:])).astype(float)

X = np.column_stack([np.ones(n), X_raw])
w = np.zeros(3)
for _ in range(5000):
    w -= 0.5 * X.T @ (sigmoid(X @ w) - y) / n

# C is the inverse regularisation strength. A huge C means "almost none",
# which is what our plain gradient descent is doing.
sk = LogisticRegression(C=1e6, max_iter=5000).fit(X_raw, y)
print("numpy  :", np.round(w, 3))
print("sklearn:", np.round(np.r_[sk.intercept_, sk.coef_[0]], 3))
print("true   :", true_w)
print("prediction agreement:",
      round(float(np.mean((sigmoid(X @ w) >= 0.5) == (sk.predict(X_raw) == 1))), 3))

numpy  : [-0.488  2.32  -1.671]
sklearn: [-0.488  2.32  -1.67 ]
true   : [-0.5  2.  -1.5]
prediction agreement: 1.0

重みは小数点以下3桁まで一致し、すべての予測が一致します。そこに至るために2つの設定が必要でした。第1に、numpyのループを5,000ステップ回して完全に収束させました。500ステップではまだ少しずつ動いていたからです。第2に、scikit-learnのLogisticRegressionはデフォルトで重みにL2ペナルティを課し、Cで制御され、重みをゼロの方へ引き寄せます。デフォルトのC=1.0では、素朴なループより小さな重みになっていたでしょう。C=1e6に設定すると、事実上ペナルティが切れます。このデフォルトは実務では妥当で、分離可能なデータで重みが暴走するのを止めてくれますが、ライブラリがこのページの損失とはわずかに異なる損失を最小化していることを意味し、数値を比べる前に知っておくべきことです。

線形モデルの位置づけ

線形回帰とロジスティック回帰は、どんな表形式の問題でも最初に試すモデルであり続けています。ミリ秒で学習でき、重みが読め（support_callsの重みが2.3なら、問い合わせ1回ごとに離脱スコアが2.3上がる）、より凝ったモデルが十分な差で上回らなければならないベースラインを設定します。限界は直線であることです。「支出の多さは新規顧客にだけ効く」といった相互作用は、その特徴量を手で作らない限り表現できません。レッスン05は、相互作用を自力で発見する木を導入し、レッスン07はシグモイドに似たユニットを積み重ねて、自ら特徴量を学習するネットワークを作ります。どちらも最後は、いま見た損失と勾配に行き着きます。

演習

1. 線形回帰のデータに2つ目の特徴量x2 = rng.uniform(0, 5, size=n)を真の重み−1.0で加え、Xを3列に拡張してください。np.linalg.solveが3つの重みすべてを復元し、LinearRegressionが一致することを確認してください。
2. ロジスティック回帰のループで、学習率0.05と5.0を試し、同じステップで損失を出力してください。それぞれ何が起きるかを記述し、その説明はレッスン04で読んでください。
3. LogisticRegressionをデフォルトのC=1.0で当てはめ、その重みをC=1e6のものと比べてください。ペナルティは重みをどちらの方向へ動かしますか。予測の一致率は変わりますか。

まとめ

• 線形回帰は、二乗誤差損失を持つ重み付き和です。損失はお椀型で、正規方程式がnp.linalg.solve1回でその底を見つけます。
• はい・いいえのラベルでは、シグモイドが重み付き和を確率に変え、対数損失がその確率を採点して、自信のある間違いを急激に罰します。
• 対数損失の勾配は、Xの転置×（予測−ラベル）の平均です。閉形式はないので、勾配降下法が損失が下がらなくなるまで勾配と逆向きに進みます。
• ゼロから書いた版は、正則化を切ってループを収束させれば、scikit-learnと一致します。ライブラリが加えるのはペナルティ、より速いソルバー、多クラス対応であって、別のモデルではありません。

さらに読む

• The Elements of Statistical Learning, Hastie, Tibshirani and Friedman: 3.2節が最小二乗解を導出し、4.4節がロジスティック回帰とその当てはめを扱います。
• Pattern Recognition and Machine Learning, Bishop: 4.3.2節が、シグモイドと交差エントロピー損失がこのページで使ったきれいな勾配を生む理由を示します。
• scikit-learn user guide, “Linear models”: LinearRegressionとLogisticRegressionの背後にある正則化オプションとソルバー。
• numpy reference, numpy.linalg.solve: 解くことが反転に勝る理由。