Skip to content

第 12 章: 正則化とモデル選択

12.1 はじめに

第 7 章の線形回帰は、訓練データの誤差だけを小さくしました。特徴量が増えて互いに似てくると、係数が極端に大きくなり、訓練データにはよく当たるのに未知のデータには当たらない 過学習 が起きます。

この章では、係数の大きさそのものに罰則を加える 正則化 を扱います。

  • リッジ回帰 — 係数の 2 乗の合計に罰則をかける。係数は 0 に近づくが 0 にはならない
  • ラッソ回帰 — 係数の絶対値の合計に罰則をかける。係数がちょうど 0 になるので、要らない特徴量が分かる

罰則の強さ alpha は自分で決める値（ハイパーパラメータ）なので、検証データ を用意して選びます。

Haskell 版の第 12 章 は、リッジ回帰の正規方程式を hmatrix にも解かせて自作と突き合わせ（差は 8.0e-15）、ラッソ回帰には突き合わせる相手がいないと書きました。なでしこ3 版には、リッジ回帰にも突き合わせる相手がいません。 行列の命令が無いので、正規方程式を解くところから自作です（ADR 015）。

そこで、第 13 章の固有値分解と同じく、答えが満たすべき性質をテストにします。

  • alpha が 0 のリッジ回帰は、第 7 章の線形回帰と同じ解 になる
  • alpha を大きくするほど、係数のノルムが小さくなる
  • リッジ回帰の解は Xᵀ (t − X w) = alpha w を満たす
  • ラッソ回帰の解は 軟しきい値の条件（KKT 条件） を満たす

なでしこ3 版では、次の 3 点に注目してください。

  • 性質のテストで両側から挟む。 2 つの正則化を、それぞれの解が満たす式で確かめます。わざと壊した実装が落ちることも確かめました
  • 新しく書いた前処理は外れ値の除去だけ。 標準化・2 次の項・3 分割は、第 2 章と第 9 章の関数を呼ぶだけで済みます
  • 座標降下法は素朴に書いても命令の上限に遠い。 報告 1 回が約 223 万命令でした。ただし上限の 1 万巡まで回ると、それだけで約 1.5 億命令になります

12.2 過学習と正則化

係数が大きくなりすぎる

特徴量どうしが強く相関していると、「片方を大きく足して、もう片方を大きく引く」ような解でも訓練データには当てはまってしまいます。係数の絶対値が跳ね上がり、入力がわずかに変わるだけで予測が大きく動くようになります。

係数の大きさに罰則を加える

最小二乗法が最小にするのは残差の 2 乗和だけです。

\[\lVert t - Xw \rVert^2\]

正則化はここに項を足します。

  • リッジ回帰: \(\lVert t - Xw \rVert^2 + \alpha \lVert w \rVert^2\)
  • ラッソ回帰: \(\frac{1}{2} \lVert t - Xw \rVert^2 + \alpha \lVert w \rVert_1\)

alpha が 0 なら最小二乗法に戻ります。大きくするほど係数は 0 のほうへ押されます。

リッジ回帰の解き方

リッジ回帰は微分して 0 と置くだけで解けます。

\[(X^\top X + \alpha I) w = X^\top t\]

第 7 章の正規方程式の左辺に、対角の alpha を足しただけです。解くところは第 7 章の 連立一次方程式求解（部分ピボット付きのガウスの消去法）がそのまま使えます。

ラッソ回帰の解き方

L1 の罰則は原点で折れているので微分できません。行列を解いて終わりにはできず、係数を 1 つずつ順に動かす 座標降下法を使います。1 つの係数を動かすときは、ほかを固定して最適な値を求め、そこに 軟しきい値作用素 をかけて 0 に寄せます。

12.3 題材とデータ

Boston.csv

ボストンの住宅価格のデータ（100 件）を使います。第 9 章・第 13 章でも使ったものです。この章で使う列は 4 つに絞ります。

列 意味
RM 1 戸あたりの平均部屋数
PTRATIO 生徒と教師の比率
LSTAT 低所得者層の割合
PRICE 住宅価格（正解） 
# 使う特徴量の列。2 次の項を作る順もこの並びで決まる（ほかの言語版と同じ）。
正則化特徴量列 = [「RM」, 「PTRATIO」, 「LSTAT」]

正解の列の名前は、第 9 章の 住宅目的列（「PRICE」）をそのまま読みます。取り込んだファイルのトップレベルの変数は、読めても書き換えられませんが、読むだけなら問題ありません。

過学習が起きやすい状況を作る

3 列だけでは過学習しません。そこで、標準化したあとに 2 次の項 を作ります。3 列から 2 乗と交互作用を足すと 9 列になります。2 次の項は元の列と強く相関するので、正則化の効き目が見える状況になります。

訓練・検証・テストの 3 つに分ける

alpha を選ぶのにテストデータを使ってしまうと、「テストデータに合わせて選んだ」ことになり、テストが未知のデータの代わりになりません。そこで 3 つに分けます。

  • 訓練データ — モデルを学習する
  • 検証データ — alpha を選ぶ
  • テストデータ — 選び終わったあと、一度だけ測る

12.4 TODO リストの作成

TODO リスト:

  • リッジ回帰を自作する
  • ラッソ回帰を座標降下法で自作する
  • 係数の絶対値の合計と、0 になった特徴量名を求める
  • alpha ごとの実験結果を記録して、検証データで選ぶ
  • z スコアで外れ値を除く
  • 標準化してから 2 次の項を作り、3 つに分ける
  • 性質のテストと参照実装で数値を確かめる
  • 実データで線形回帰とリッジ回帰を比べる

12.5 リッジ回帰を自作する

Red: 第 7 章の線形回帰と同じになるテスト

Haskell 版と同じ「答えが分かっているデータ」（切片 1、係数 2 と −3）から始めます。

正方特徴量 = [{"a": 1, "b": 0}, {"a": 0, "b": 1}, {"a": 1, "b": 1}]
正方正解 = [3, -2, 0]
二列 = [「a」, 「b」]

# リッジ回帰
最小二乗 = 正方特徴量と正方正解を二列で線形回帰学習
零リッジ = 正方特徴量と正方正解を二列と0でリッジ回帰学習
「alpha が 0 なら第 7 章の線形回帰と同じ切片になる」で(零リッジ["切片"])と(最小二乗["切片"])を0.000000001で近似確認
「alpha が 0 なら第 7 章の線形回帰と同じ係数になる」で(零リッジ["係数一覧"])と(最小二乗["係数一覧"])を0.000000001で近似確認

alpha = 0 なら罰則が消えるので、第 7 章の 線形回帰学習 と同じ答えになるはずです。すでに動いているものを答えに使える ので、期待値を手で計算する必要がありません。ここまでは Haskell 版と同じです。

ただし、なでしこ3 版ではこのテストの意味が少し重くなります。Haskell 版には hmatrix という別の解き方があり、「自作の正規方程式が正しい」ことはそちらでも確かめられました。なでしこ3 版では、第 7 章の線形回帰が、正規方程式の組み立て方の違う唯一の比較相手 です。第 7 章は切片の列（すべて 1 の列）を足した計画行列で解き、この章は平均を引いてから解きます。組み立て方が違う 2 つの式が同じ答えになるので、どちらかの組み立てを間違えれば落ちます。

Green: 中心化して解く

●(特徴量一覧と正解一覧を列名一覧と強度で)リッジ回帰学習とは
    問題 = 特徴量一覧と正解一覧を列名一覧で正則化問題中心化
    列一覧 = 問題["中心化列一覧"]
    正規行列 = []
    右辺一覧 = []
    Iを 0から (列一覧の要素数) - 1まで繰り返す
        行 = []
        Jを 0から (列一覧の要素数) - 1まで繰り返す
            値 = (列一覧[I])と(列一覧[J])の正則化内積算出
            # 対角に alpha を足す。切片は中心化で外してあるので、罰則は係数だけにかかる。
            もし、I = Jならば
                値 = 値 + 強度
            ここまで
            行に値を配列追加
        ここまで
        正規行列に行を配列追加
        右辺一覧に((列一覧[I])と(問題["中心化正解"])の正則化内積算出)を配列追加
    ここまで
    重み一覧 = 正規行列と右辺一覧を連立一次方程式求解
    問題と列名一覧を重み一覧で正則化モデル作成して戻る
ここまで

第 7 章との違いは 2 つだけです。

  1. 平均を引いてから解く — 切片に罰則をかけないため。切片まで 0 に押されると、予測全体が下に引っ張られます
  2. 対角に alpha を足す — もし、I = Jならば の 3 行

\(X^\top X\) の成分は「中心化した X の列」どうしの内積なので、列の配列を持っておけば、行列の積を書かずに済みます。 前準備は、リッジ回帰とラッソ回帰で同じなので 1 つの関数にまとめました。

●(特徴量一覧と正解一覧を列名一覧で)正則化問題中心化とは
    もし、(特徴量一覧の要素数) = 0ならば
        「訓練データが空です」のエラー発生
    ここまで
    もし、(特徴量一覧の要素数) ≠ (正解一覧の要素数)ならば
        「特徴量と正解の件数が違います: {特徴量一覧の要素数} と {正解一覧の要素数}」のエラー発生
    ここまで
    列一覧 = (特徴量一覧と列名一覧の特徴量行列変換)を行列転置
    平均一覧 = []
    中心化列一覧 = []
    Jを 0から (列一覧の要素数) - 1まで繰り返す
        平均 = (列一覧[J])の正則化平均算出
        新列 = []
        Iを 0から ((列一覧[J])の要素数) - 1まで繰り返す
            新列に(列一覧[J][I] - 平均)を配列追加
        ここまで
        平均一覧に平均を配列追加
        中心化列一覧に新列を配列追加
    ここまで
    正解平均 = 正解一覧の正則化平均算出
    中心化正解 = []
    Iを 0から (正解一覧の要素数) - 1まで繰り返す
        中心化正解に(正解一覧[I] - 正解平均)を配列追加
    ここまで
    {"平均一覧": 平均一覧, "正解平均": 正解平均, "中心化列一覧": 中心化列一覧, "中心化正解": 中心化正解}で戻る
ここまで

特徴量を数値の行に並べるのは第 9 章の 特徴量行列変換、転置は linear_algebra.nako3 の 行列転置 です。どちらも書き換えずに呼んでいます。

モデルは第 7 章のものを使い回す

戻り値は第 7 章と同じ辞書 {"切片", "列名一覧", "係数一覧"} です。切片は「正解の平均 − 列の平均と係数の内積」で、第 7 章の 線形回帰モデル作成 に渡して作ります。

●(問題と列名一覧を重み一覧で)正則化モデル作成とは
    切片 = 問題["正解平均"] - ((問題["平均一覧"])と重み一覧の正則化内積算出)
    切片と列名一覧を重み一覧で線形回帰モデル作成して戻る
ここまで

この章では新しいモデルの形を作っていません。 予測は第 7 章の 線形回帰予測、評価は 決定係数算出 をそのまま使います。

性質でも押さえる

Haskell 版と同じく、切片に罰則がかからないことを確かめます。

極リッジ = 正方特徴量と正方正解を二列と100000000でリッジ回帰学習
「切片には罰則がかからないので、強く正則化すると予測は正解の平均に寄る」で(極リッジで正方特徴量を線形回帰予測)と[1 ÷ 3, 1 ÷ 3, 1 ÷ 3]を0.000001で近似確認

alpha を極端に大きくすると係数は 0 に押し切られ、残るのは切片だけです。その切片が正解の平均（1/3）になる ことが、中心化が効いている証拠になります。

突き合わせる相手がいない分、解が満たす式も確かめます。リッジ回帰の目的関数を微分して 0 と置いた式を、並べ替えると次の形になります。

\[X^\top (t - Xw) = \alpha w\]

テストの側に、この式の両辺のずれの最大を求める関数を置きました。

●(問題と重み一覧を強度で)検査用リッジ条件違反とは
    残差一覧 = 問題と重み一覧の検査用残差算出
    列一覧 = 問題["中心化列一覧"]
    最大 = 0
    Jを 0から (列一覧の要素数) - 1まで繰り返す
        逸脱 = (((列一覧[J])と残差一覧の正則化内積算出) - 強度 × 重み一覧[J])の絶対値
        もし、逸脱 > 最大ならば
            最大 = 逸脱
        ここまで
    ここまで
    最大で戻る
ここまで

小リッジ = 小訓練と小正解を小列名と10でリッジ回帰学習
「リッジ回帰の解は Xᵀ (t − X w) = alpha w を満たす」で(小問題と(小リッジ["係数一覧"])を10で検査用リッジ条件違反)と0を0.000000001で近似確認
「わざとずらした係数はリッジ回帰の式を満たさない」で((小問題と(小リッジ["係数一覧"])を11で検査用リッジ条件違反) > 0.001)と真が一致確認

2 本目は、検査の関数そのものが空振りしていないことの確かめです。alpha を 11 と偽って渡すと、ずれは 0.001 を超えます（小問題 などは 12.11 節の自作のデータです）。

最後に、alpha を大きくするほど係数の 2 乗の合計の平方根（L2 ノルム）が小さくなることを、候補の 5 つの alpha すべてで確かめました。リッジ回帰が小さくするのは L2 ノルムなので、絶対値の合計ではなくこちらで単調性を見ます。

12.6 ラッソ回帰を座標降下法で自作する

軟しきい値作用素

●(値と閾値の)軟閾値算出とは
    もし、値 > 閾値ならば
        値 - 閾値で戻る
    ここまで
    もし、値 < (0 - 閾値)ならば
        値 + 閾値で戻る
    ここまで
    0で戻る
ここまで

しきい値の外なら 0 のほうへ縮め、中なら 0 にします。符号を反転した値は、第 13 章と同じく 0 - 式 の形で (0 - 閾値) と書きました（-(式) は解析エラーになるため）。Haskell 版のテスト 3 本に、しきい値ちょうど（2と2・(0 - 2)と2）が 0 になるテストを足しました。

係数と残差をその場で書き換える

Haskell 版は再代入が書けないので、係数を Map に持ち、「係数・残差・いちばん大きな動き」の 3 つを組にして畳み込みで持ち回りました。なでしこ3 版は配列をその場で書き換えられるので、PHP 版に近い 1 つの繰り返しになります。

●(特徴量一覧と正解一覧を列名一覧と強度で)ラッソ回帰学習とは
    問題 = 特徴量一覧と正解一覧を列名一覧で正則化問題中心化
    列一覧 = 問題["中心化列一覧"]
    二乗和一覧 = []
    重み一覧 = []
    Jを 0から (列一覧の要素数) - 1まで繰り返す
        二乗和一覧に((列一覧[J])と(列一覧[J])の正則化内積算出)を配列追加
        重み一覧に0を配列追加
    ここまで
    残差一覧 = (問題["中心化正解"])を配列複製
    件数 = 残差一覧の要素数
    周数 = 0
    (周数 < ラッソ反復上限)の間
        最大変化 = 0
        Jを 0から (列一覧の要素数) - 1まで繰り返す
            # 値がすべて同じ列は動かしようがないので 0 のままにする（0 で割らない）。
            もし、(二乗和一覧[J]) = 0ならば
                続ける
            ここまで
            現列 = 列一覧[J]
            旧値 = 重み一覧[J]
            # いったん列の寄与を残差に戻してから、残差との相関で係数を決め直す。
            復元一覧 = []
            Iを 0から 件数 - 1まで繰り返す
                復元一覧に(残差一覧[I] + 旧値 × 現列[I])を配列追加
            ここまで
            新値 = ((現列と復元一覧の正則化内積算出)と強度の軟閾値算出) ÷ 二乗和一覧[J]
            残差一覧 = []
            Iを 0から 件数 - 1まで繰り返す
                残差一覧に(復元一覧[I] - 新値 × 現列[I])を配列追加
            ここまで
            重み一覧[J] = 新値
            変化 = (新値 - 旧値)の絶対値
            もし、変化 > 最大変化ならば
                最大変化 = 変化
            ここまで
        ここまで
        周数 = 周数 + 1
        もし、最大変化 < ラッソ収束幅ならば
            抜ける
        ここまで
    ここまで
    問題と列名一覧を重み一覧で正則化モデル作成して戻る
ここまで

打ち切りの条件・更新の順・軟しきい値のかけ方は、Haskell 版の lassoFit にそろえました。

項目 決め方
打ち切り 1 巡でいちばん大きな係数の動きが 1e-10（ラッソ収束幅）より小さくなったら、その巡の結果で止める。上限は 1 万巡（ラッソ反復上限）
更新の順 列の順に 1 つずつ。動かした係数はすぐ残差に反映し、次の列はその残差で決める
残差の直し方 古い係数の寄与を足し戻した 復元一覧 を作り、新しい係数の寄与を引いて新しい残差にする（差分だけを引く形にはしない）
軟しきい値 列と 復元一覧 の内積に、しきい値 alpha のままかける（2 で割らない）

最後の 2 行は、計算の結果を Haskell 版とビット単位でそろえるための選択です。Go 版 は残差を「新旧の差 × 列」だけ直し、しきい値を alpha/2 にしています（目的関数の定数倍が違う）。どちらも正しい座標降下法ですが、浮動小数点の丸めが変わり、alpha の尺度も変わります。

「書き換わる変数」を組にして持ち回る Haskell 版に比べると、なでしこ3 版は 書き換える変数がそのまま名前で並ぶ ので、数式との対応は追いやすくなります。その代わり、配列は参照で渡るので、問題["中心化正解"] を書き換えないよう、最初に 配列複製 してから残差として使います。

名前に使えなかった 2 つの言葉

この関数を書いたとき、文法の検査が次のように止まりました。

[文法エラー]src/chapter12.nako3(135行目): 関数『chapter12__ラッソ回帰学習』の定義で以下のエラーがありました。
[文法エラー]src/chapter12.nako3(135行目): 『ここまで』がありません。関数定義の末尾に必要です。『=』の前に『ここまで』を記述してください。

原因は、巡回の回数を数える変数を 巡回 と名付けたことでした。回 は 3回 のように繰り返しを表す語なので、巡回 = 0 の 回 が繰り返しとして読まれ、そのあとの = の置き場所が無くなったと考えられます。周数 に変えると通りました。テストの側でも、ずれを表す変数を 違反 と名付けたら「『違』の前に『ここまで』を記述してください」で止まりました。違えば の 違 と読まれたようです。逸脱 に変えています。

どちらも、使い捨てのファイル（関数の中に 巡回 = 0 や 違反 = 0 とだけ書いたもの）で同じエラーになることを確かめました。エラーは定義の先頭の行を指し、原因の行を指しません。 これまでに見た「助詞やひらがなを名前に入れない」「組み込みの名前と重ねない」に加えて、回 や 違 のように構文の語と同じ漢字で始まる、または終わる名前も避けます。

Red: 0 になることを確かめる

強ラッソ = 正方特徴量と正方正解を二列と1000でラッソ回帰学習
「alpha が大きければ係数はちょうど 0 になる」で(強ラッソ["係数一覧"])と[0, 0]が一致確認
「係数が 0 になった特徴量の名前を返す」で(強ラッソの零係数列名抽出)と[「a」, 「b」]が一致確認

リッジ回帰の同じテストは 近似確認 でした。一致確認 と 近似確認 の使い分けが、2 つの正則化の違いそのもの になっています。Haskell 版の shouldBe と近似比較の使い分けと同じです。

もう片方の端も押さえます。

弱ラッソ = 正方特徴量と正方正解を二列と0.00000001でラッソ回帰学習
「alpha が十分小さければ最小二乗法に近づく」で(弱ラッソ["係数一覧"])と(最小二乗["係数一覧"])を0.000001で近似確認

値がすべて同じ列（二乗和が 0）の道もテストしました。その列の係数は 0 のまま残り、ほかの列は学習します（b の係数が 1 − 0.0005 になることを手で計算して確かめました）。

KKT 条件をテストにする

両端を挟むだけでは、途中の alpha で正しいかは分かりません。Haskell 版はそこを「ほかの言語版と 0 になった 3 列が一致すること」で支えました。なでしこ3 版は、それに加えて ラッソ回帰の解が満たすべき条件 をテストにします。

勾配を g = Xᵀ (t − X w) とすると、ラッソ回帰の最適な解では、列ごとに次が成り立ちます。

  • 係数が正なら g = alpha、負なら g = −alpha
  • 係数が 0 なら |g| ≦ alpha

2 つ目が「軟しきい値の中なら 0 にする」に対応します。この条件（KKT 条件）からのずれの最大を求める関数をテストに置きました。

●(問題と重み一覧を強度で)検査用ラッソ条件違反とは
    残差一覧 = 問題と重み一覧の検査用残差算出
    列一覧 = 問題["中心化列一覧"]
    最大 = 0
    Jを 0から (列一覧の要素数) - 1まで繰り返す
        勾配 = (列一覧[J])と残差一覧の正則化内積算出
        逸脱 = (勾配 + 強度)の絶対値
        もし、(重み一覧[J]) > 0ならば
            逸脱 = (勾配 - 強度)の絶対値
        ここまで
        もし、(重み一覧[J]) = 0ならば
            逸脱 = ((勾配)の絶対値) - 強度
        ここまで
        もし、逸脱 > 最大ならば
            最大 = 逸脱
        ここまで
    ここまで
    最大で戻る
ここまで

小ラッソ = 小訓練と小正解を小列名と2でラッソ回帰学習
零件数 = (小ラッソの零係数列名抽出)の要素数
「alpha 2 のラッソ回帰は 0 の係数と 0 でない係数を両方持つ」で[(零件数 > 0), (零件数 < 9)]と[真, 真]が一致確認
「ラッソ回帰の解は KKT 条件（軟しきい値の条件）を満たす」で(小問題と(小ラッソ["係数一覧"])を2で検査用ラッソ条件違反)と0を0.000001で近似確認
「わざとずらした強度では KKT 条件を満たさない」で((小問題と(小ラッソ["係数一覧"])を3で検査用ラッソ条件違反) > 0.1)と真が一致確認

1 本目は、条件の両方の場合（0 の係数と 0 でない係数）が実際に試されていることを確かめるためのものです。どちらかしか無ければ、KKT 条件のテストは片方の場合しか見ていないことになります。

わざと壊して落ちることも確かめました。 1 か所ずつ壊してテストを走らせた結果です。

壊し方 落ちたテスト
軟しきい値のしきい値を 強度 ÷ 2 にする（Go 版の尺度） KKT 条件・値がすべて同じ列・報告の 3 本
係数を決め直す前に、古い係数の寄与を残差に戻さない 最小二乗法に近づく・値がすべて同じ列・KKT 条件・報告の 4 本
リッジ回帰の対角に 強度 × 2 を足す 正規方程式・実験表とテストデータの決定係数の突き合わせ・報告の 4 本
最良の実験を、同じ値でも入れ替える（> を ≧ に） 同じ値なら先の実験を選ぶ、の 1 本
外れ値の標準偏差を n で割る 12.9 節で足した境界のテストの 1 本

12.7 実験結果を記録して選ぶ

alpha ごとの記録は辞書にしました。

●(分割データと強度一覧の)リッジ実験実行とは
    訓練特徴量 = 分割データ["訓練特徴量"]
    訓練正解 = 分割データ["訓練正解"]
    列名一覧 = 分割データ["列名一覧"]
    結果 = []
    Iを 0から (強度一覧の要素数) - 1まで繰り返す
        モデル = 訓練特徴量と訓練正解を列名一覧と(強度一覧[I])でリッジ回帰学習
        訓練決定係数 = 訓練正解と(モデルで訓練特徴量を線形回帰予測)の決定係数算出
        検証決定係数 = (分割データ["検証正解"])と(モデルで(分割データ["検証特徴量"])を線形回帰予測)の決定係数算出
        結果に{"強度": 強度一覧[I], "訓練決定係数": 訓練決定係数, "検証決定係数": 検証決定係数, "係数絶対値合計": モデルの係数絶対値合計算出}を配列追加
    ここまで
    結果で戻る
ここまで

選ぶところは「検証データの決定係数が最大」です。

●(実験一覧の)最良実験選択とは
    もし、(実験一覧の要素数) = 0ならば
        「実験結果が 1 件もありません」のエラー発生
    ここまで
    最良 = 実験一覧[0]
    Iを 1から (実験一覧の要素数) - 1まで繰り返す
        もし、(実験一覧[I]["検証決定係数"]) > (最良["検証決定係数"])ならば
            最良 = 実験一覧[I]
        ここまで
    ここまで
    最良で戻る
ここまで

Haskell 版は maximumBy が同値なら後ろを返すので畳み込みを自作しました。なでしこ3 には最大の要素を選ぶ命令が無いので、最初から自分で書きます。同じ値なら先の実験（小さい alpha） を選ぶよう、厳密に大きいときだけ入れ替えます。> を ≧ に変えると、「同じ値なら先の実験を選ぶ」テストだけが落ちました。

12.8 0 になった係数の特徴量名を返す

●(モデルの)零係数列名抽出とは
    結果 = []
    列名一覧 = モデル["列名一覧"]
    Iを 0から (列名一覧の要素数) - 1まで繰り返す
        もし、モデル["係数一覧"][I] === 0ならば
            結果に(列名一覧[I])を配列追加
        ここまで
    ここまで
    結果で戻る
ここまで

浮動小数点を等号で比べるのは普通なら赤信号ですが、ここは軟しきい値が 0 そのものを返すことが分かっています。「0 に近い」ではなく「0 である」を見たいので、これが正しい比較です。なでしこ3 の = は JavaScript の == と同じゆるい比較（0 = 「」 が真）なので、型まで比べる === を使いました（第 2 章）。

12.9 最小限の前処理

標準化してから 2 次の項を作り、3 つに分ける

●(CSV文字列を外側割合と内側割合とシードで)正則化住宅前処理とは
    除去表 = (CSV文字列を表読込)と[住宅目的列, 正則化特徴量列[0], 正則化特徴量列[1], 正則化特徴量列[2]]を正則化外れ値閾値で標準得点外れ値除去
    行一覧 = 除去表["行一覧"]
    特徴量一覧 = []
    ラベル一覧 = []
    Iを 0から (行一覧の要素数) - 1まで繰り返す
        特徴量 = {}
        Jを 0から (正則化特徴量列の要素数) - 1まで繰り返す
            特徴量[正則化特徴量列[J]] = 行一覧[I]の(正則化特徴量列[J])を欠損許容数値取得
        ここまで
        特徴量一覧に特徴量を配列追加
        ラベル一覧に(行一覧[I]の住宅目的列を列取得)を配列追加
    ここまで
    外側 = 特徴量一覧とラベル一覧を外側割合とシードで訓練テスト分割
    内側 = (外側["訓練特徴量"])と(外側["訓練正解"])を内側割合とシードで訓練テスト分割
    標準化器 = (内側["訓練特徴量"])と正則化特徴量列で標準化器学習
    {"訓練特徴量": (標準化器で(内側["訓練特徴量"])を標準化器適用)と正則化特徴量列で多項式特徴量展開, …（検証とテストも同じ）…, "列名一覧": 正則化特徴量列の多項式列名一覧作成, "残件数": 行一覧の要素数}で戻る
ここまで

（最後の行は、検証データとテストデータの 4 つの欄を省いて載せています。）

使っている関数の出どころは次のとおりです。この章で新しく書いた前処理は、外れ値の除去（標準得点外れ値除去）だけです。

関数 出どころ 役割
訓練テスト分割 第 2 章 java.util.Random と同じ線形合同法で並べ替えて分ける。2 回呼んで 3 つに分ける
標準化器学習・標準化器適用 第 9 章 件数 n で割る標準偏差で標準化する。平均と標準偏差は訓練データだけから求める
多項式特徴量展開・多項式列名一覧作成 第 9 章 2 乗と交互作用の項を、scikit-learn の PolynomialFeatures と同じ名前と順で作る
住宅目的値変換 第 9 章 正解の値を数値にする（空欄ならエラー）

外側で訓練用とテスト用に分け、外側の訓練用をもう一度分けます。内側の「テストデータ」が検証データになるので、返す辞書で名前を付け替えます。第 2 章の関数を作り直さずに済みました。

第 9 章の 住宅前処理 はそのまま使えませんでした。第 9 章は CRIME をダミー変数にし、2 つに分けてから欠損値を平均で補完する前処理で、この章は外れ値を除いて 3 つに分けるからです。前処理の手順は別に書き、部品の関数だけを共有しています。

z スコアで外れ値を除く

第 9 章は四分位範囲で外れ値を判定しました。この章は z スコア で、標準偏差は 件数 n − 1 で割る標本標準偏差 を使います（標準化の n とは違います）。

        統計一覧に[値一覧, 平均, (二乗和 ÷ ((値一覧の要素数) - 1))の平方根]を配列追加

Haskell 版が書いたとおり、n 件のデータの z スコアには (n − 1) ÷ √n という上限があり、10 件では 2.846 までしか届きません。Haskell 版の 2 本のテスト（20 件なら 1000 が除かれる、10 件なら除かれない）をそのまま書きました。

「z スコアがしきい値を超える行を除く」で((((二十件と[「v」]を3で標準得点外れ値除去)["行一覧"])の要素数))と19が一致確認
「件数が少ないと z スコアの上限 (n − 1) ÷ √n に届かず 1 件も除かれない」で((((十件と[「v」]を3で標準得点外れ値除去)["行一覧"])の要素数))と10が一致確認

Haskell 版の 2 本は、標準偏差を n で割っても通ってしまいました（12.6 節の表）。10 件の例では、n で割ったときの上限 √(n − 1) がちょうど 3 で、しきい値を超えないからです。そこで、n − 1 で割れば残り、n で割れば除かれる値を探して、境界のテストを 1 本足しました。

境界 = 「v{改行}0{改行}1{改行}2{改行}3{改行}4{改行}0{改行}1{改行}2{改行}3{改行}4{改行}0{改行}13{改行}」を表読込
「標準偏差は n − 1 で割る（13 の z スコアは 2.89 で残る。n で割ると 3.02 になり除かれてしまう）」で((((境界と[「v」]を3で標準得点外れ値除去)["行一覧"])の要素数))と12が一致確認

空欄のテストも Haskell 版と同じ形です。Haskell 版は「1 列の CSV では空欄の行が空の行として読み飛ばされる」ことに刺されて列を 2 つにしました。なでしこ3 版の 表読込 も、すべてのセルが空の行を読み飛ばします（第 1 章）。同じ理由で、同じ 2 列の CSV（v,w{改行}1,1{改行},2{改行}2,3{改行}）を使っています。

12.10 ほかの言語版の数値と突き合わせる

突き合わせる相手の代わり

Haskell 版はこの節で hmatrix と突き合わせました。なでしこ3 版の支えは次の 4 つです。

支え 確かめられること
第 7 章の線形回帰 alpha が 0 のリッジ回帰が、組み立て方の違う正規方程式と同じ解になる（12.5 節）
解が満たす式 リッジ回帰は Xᵀ (t − X w) = alpha w、ラッソ回帰は KKT 条件。係数のノルムは alpha とともに小さくなる（12.5・12.6 節）
自作データと参照実装 自作のボストン風データで、報告の文字列と完全な精度の値が Python の参照実装と一致する（下）
ほかの言語版の数値 実データの実験表・選ばれる alpha・テストデータの決定係数・ラッソ回帰が 0 にした 3 列（12.11 節。データのある環境で確かめる）

Haskell 版の報告の「hmatrix のリッジ回帰との係数の最大の差」の行は、なでしこ3 版では「第 7 章の線形回帰との係数の最大の差（alpha=0.0）」に置き換えました。差そのものは組み立て方の違いで 1e-13 前後になり、最後の桁が実装に左右されるので、1e-9 より小さければ 1.0e-9 未満 と表示します。

正則化の強さの尺度をそろえる

自作の目的関数は件数で割りません。Haskell 版と同じく、件数で割る実装（Java 版・Clojure 版が使う Tribuo の ElasticNetCDTrainer）の alpha = 0.5 を、訓練データ 47 件にあわせて 0.5 × 47 = 23.5 と読み替えます。

ラッソ強度 = 23.5

自作のデータと参照実装で突き合わせる

実データと同じ形の「ボストン風」のデータ（100 行・14 列、CRIME が 3 種類、NOX に 3 か所・RAD に 4 か所の欠損、RM・PTRATIO・LSTAT・PRICE は欠損なし。z スコアが 3 を超える行を 2 つ入れた架空のデータで、配布データではありません）を作り、Python で書いた参照実装と突き合わせました。参照実装は Haskell 版と同じ規則（n − 1 の z スコア・java.util.Random と同じ分割を 2 回・n の標準化・2 次の項の順・中心化したリッジ回帰・同じ打ち切りの座標降下法）を、演算の順までそろえて書いたものです。連立一次方程式は、なでしこ3 版と同じ部分ピボット付きのガウスの消去法で解きます。

なでしこ3 版の出力です。

データ件数: 98（外れ値 2 件を除外）
訓練データ: 47 件, 検証データ: 21 件, テストデータ: 30 件
特徴量: RM, PTRATIO, LSTAT, RM^2, RM PTRATIO, RM LSTAT, PTRATIO^2, PTRATIO LSTAT, LSTAT^2
alpha  訓練 R²  検証 R²  係数の絶対値の合計
  0.0  0.8109  0.4917  12.433
  0.1  0.8109  0.5100  12.181
  1.0  0.8088  0.6054  10.762
 10.0  0.7859  0.7214  7.360
100.0  0.5746  0.4819  3.617
検証データで選んだ alpha: 10.0
テストデータの決定係数: 線形回帰 0.7735, リッジ回帰 0.7983
第 7 章の線形回帰との係数の最大の差（alpha=0.0）: 1.0e-9 未満
ラッソ回帰（alpha=23.5）で係数が 0 になった特徴量: RM^2, RM PTRATIO, RM LSTAT, PTRATIO^2, PTRATIO LSTAT

報告は参照実装と 1 文字も違いませんでした。 完全な精度の値 37 個（実験表の 15 個、テストデータの決定係数 2 個、alpha 10 のリッジ回帰と alpha 23.5 のラッソ回帰の切片と係数 20 個）も、参照実装とビット単位で一致しました。乱数の種を変えた自作データをさらに 5 通り作り（除かれる外れ値が 3 件のもの、ラッソ回帰が 0 にする列が 3〜5 列と変わるものを含む）、6 通りすべてで報告も 37 個の値も一致しました。ラッソ回帰の巡回は 51〜118 巡で、上限の 1 万巡には遠く届きませんでした。

これは自作のデータでの値で、実データの値ではありません。

命令の数

gonako は 1 回の実行で 2 億命令を超えると止まります（ADR 015）。計測用に作り直した処理系で、自作のボストン風データでの命令の数を測りました。

計測 命令の数
第12章実行 1 回（読み込みから報告まで） 約 223 万
うち、ラッソ回帰 1 回（77 巡） 約 116 万（1 巡あたり約 1.5 万）
テストファイル 1 つ（データなし） 約 111 万
テストファイル 1 つ（自作データを実データの場所に置いたとき） 約 448 万（乱数の種を変えた別の自作データで約 508 万）

目安の 1.5 億命令には遠く、第 8・10・14 章のような書き直しは要りませんでした。ただし、ラッソ回帰が上限の 1 万巡まで回ると、1 巡あたり約 1.5 万命令で、ラッソ回帰だけで約 1.5 億命令になります。 実データで収束が遅ければ、報告全体が上限に近づきます。自作データで 51〜118 巡だったこと、Haskell 版が同じ打ち切りの条件で 0 になった 3 列をほかの言語版とそろえていることから、実データでも上限には届かないと見込んでいますが、実データでの巡回の数はまだ測っていません。時間は 第12章実行 1 回が起動を含めて約 0.8 秒でした。

12.11 実データで比べる

報告を組み立てる

第12章実行 は、ファイルを読んで文字列を 正則化報告作成 に渡すだけです。

●第12章実行とは
    (「Boston.csv」のデータパス取得)を開いて正則化報告作成で戻る
ここまで

報告を組み立てる関数は文字列を受け取るので、30 行の自作のデータ（NOX・RM・PTRATIO・LSTAT・PRICE の 5 列、NOX に空欄、z スコアが 3 を超える行を 2 つ）で直接呼ぶテストを置いています。実データの分岐は、データの無い環境では一度も実行されないためです。期待値は Python の参照実装で求めました。

データ件数: 28（外れ値 2 件を除外）
訓練データ: 13 件, 検証データ: 6 件, テストデータ: 9 件
特徴量: RM, PTRATIO, LSTAT, RM^2, RM PTRATIO, RM LSTAT, PTRATIO^2, PTRATIO LSTAT, LSTAT^2
alpha  訓練 R²  検証 R²  係数の絶対値の合計
  0.0  0.9545  0.2119  10.005
  0.1  0.9543  0.2644  9.740
  1.0  0.9441  0.5009  8.135
 10.0  0.7524  0.4595  4.184
100.0  0.2242  0.0823  0.956
検証データで選んだ alpha: 1.0
テストデータの決定係数: 線形回帰 0.2746, リッジ回帰 0.6674
第 7 章の線形回帰との係数の最大の差（alpha=0.0）: 1.0e-9 未満
ラッソ回帰（alpha=23.5）で係数が 0 になった特徴量: PTRATIO, LSTAT, RM^2, RM PTRATIO, RM LSTAT, PTRATIO^2, PTRATIO LSTAT, LSTAT^2

訓練データが 13 件しかないので、alpha 23.5 のラッソ回帰は RM 以外の 8 列を 0 にしました。この小さなデータでは alpha 1.0 が選ばれ、100 件のデータとは違う alpha を選ぶ道もテストで通っています。実験表の 15 個の値とテストデータの決定係数 2 個も、参照実装と 1e-12 以内で一致することを別のテストで確かめています。

alpha の列は 小数表示 で 1 桁にしてから、組み込みの 空白埋 で 5 文字の幅に右寄せしました（「0.0」を5で空白埋 が 「 0.0」 になることを確かめました）。Haskell 版の printf "%5s" にあたります。

実データの期待する出力

実データのテストは、Haskell 版と同じ期待値で書き、データが無ければ保留にします。期待する出力は次のとおりです。12 行目（第 7 章の線形回帰との差）だけがなでしこ3 版で置き換えた行で、ほかの 12 行は Haskell 版の実測と同じです。

データ件数: 98（外れ値 2 件を除外）
訓練データ: 47 件, 検証データ: 21 件, テストデータ: 30 件
特徴量: RM, PTRATIO, LSTAT, RM^2, RM PTRATIO, RM LSTAT, PTRATIO^2, PTRATIO LSTAT, LSTAT^2
alpha  訓練 R²  検証 R²  係数の絶対値の合計
  0.0  0.8827  0.7272  14.187
  0.1  0.8827  0.7274  14.104
  1.0  0.8823  0.7288  13.594
 10.0  0.8681  0.7349  11.573
100.0  0.6583  0.5985  5.684
検証データで選んだ alpha: 10.0
テストデータの決定係数: 線形回帰 0.5224, リッジ回帰 0.6243
第 7 章の線形回帰との係数の最大の差（alpha=0.0）: 1.0e-9 未満
ラッソ回帰（alpha=23.5）で係数が 0 になった特徴量: PTRATIO^2, PTRATIO LSTAT, LSTAT^2

実データでの実測について

この章は学習データの無い環境で書いたため、上の出力はまだなでしこ3 版の実装で実測していません。実データでのラッソ回帰の巡回の数も未測定です。学習データを配置した環境で npx gulp apps:check:nadesiko3 を実行し、保留が消えて通ることを確かめてから、この注記を外します。自作のボストン風データ 6 通りでは、報告の文字列と完全な精度の 37 個の値が Python の参照実装とビット単位で一致することを確かめています（12.10 節）。また、自作データを ML_DATA_DIR で指して検査を走らせ、実データの分岐が書き間違いのエラーで止まらず、値の不一致（上の期待値と自作データの値の違い）だけで落ちることを確かめました。

実データの分岐では、報告の文字列のほかに、Haskell 版の実データのテストと同じ性質も確かめます。

    「検証データで選ぶ alpha がほかの言語版と一致する」で(実最良["強度"])と10が一致確認
    「検証データで選んだリッジ回帰はテストデータで線形回帰を上回る」で((実データと(実最良["強度"])の正則化テスト決定係数算出) > (実データと0の正則化テスト決定係数算出))と真が一致確認
    実訓練零 = 実正解と(実零リッジで実訓練を線形回帰予測)の決定係数算出
    実訓練十 = 実正解と((実訓練と実正解を実列名と10でリッジ回帰学習)で実訓練を線形回帰予測)の決定係数算出
    「訓練データでは線形回帰がリッジ回帰を上回る」で(実訓練零 > 実訓練十)と真が一致確認

過学習とは「訓練でよくてテストで悪い」ことだ、という定義そのもの が、2 本のテストで言葉になります。Haskell 版と同じ組み立てです。

結果を読む

期待する出力（ほかの言語版の実測）から読めることは、Haskell 版と同じです。

  • alpha を大きくするほど、係数の絶対値の合計は 14.187 から 5.684 へ小さくなり、訓練データの決定係数は下がり続けます
  • 検証データの決定係数は alpha = 10.0 で最も高く（0.7349）、100.0 では訓練・検証とも大きく下がります。正則化が強すぎて学習不足になっています
  • 検証データで選んだ alpha = 10.0 のリッジ回帰は、テストデータの決定係数が 0.6243 で、線形回帰の 0.5224 を上回ります。訓練データでは線形回帰のほうが高い（0.8827）のに、未知のデータでは正則化したモデルのほうがよく当たっています
  • ラッソ回帰は、9 列のうち PTRATIO^2・PTRATIO LSTAT・LSTAT^2 の 3 列の係数をちょうど 0 にします

100 件ほどのデータでは、分け方やデータのわずかな違いで結論まで変わりえます。Haskell 版の 12.11 節は、Kotlin 版で分割の乱数が違ったためにリッジ回帰がテストで線形回帰を下回った例を挙げています。12.10 節の自作データ 6 通りでも、検証データで選ばれた alpha は 10.0 が 4 通り、1.0 と 0.0 が 1 通りずつで、選んだリッジ回帰がテストデータで線形回帰を上回ったのは 4 通り、同じ（alpha 0 が選ばれた）が 1 通り、下回ったのが 1 通りでした。1 回の分け方に頼らない方法として、第 11 章の交差検証を組み合わせられます。

12.12 品質チェック

./tools/check.sh

文法・整形・テスト・テストが最後の報告まで届いたか・関数の網羅の 5 つをまとめて確かめます（第 6 章）。この章のテストは 52 件（保留 3 件）で、学習データがあれば保留の 3 件に代わって 6 件の検査が走ります。関数の網羅は次のとおりでした（この章を書いた作業ツリーでの値です）。

関数の網羅: 239 件中 239 件（直接 198 件、間接 41 件）

この章で新しく出会ったなでしこ3 の振る舞いは、12.6 節の「構文の語と同じ漢字を含む名前」です。巡回（回）と 違反（違）は、どちらも変数名に使えず、定義の先頭の行を指す分かりにくいエラーになりました。

関数の名前は、取り込んだファイルをまたいで一意にする必要があります（ADR 015）。この章の関数には 正則化・リッジ・ラッソ などの題材の名詞を入れ、第 13 章の 主成分内積 のような既存の名前と重ならないようにしました。内積と平均は第 13 章にも似た関数がありますが、この章では足す順をこの章の中で決めておきたかったので、正則化内積算出・正則化平均算出 を別に書いています。

12.13 可視化について

alpha と決定係数のグラフ、係数の軌跡（正則化パス）の描き方は Python 版 と Kotlin 版 を参照してください。この章では 12.11 節の表が同じ役割を果たしています。

12.14 まとめ

この章では、リッジ回帰とラッソ回帰をなでしこ3 の TDD で自作し、検証データでモデルを選びました。

作ったもの 支え
リッジ回帰 第 7 章の線形回帰（alpha 0 で一致）、Xᵀ (t − X w) = alpha w、係数のノルムの単調な減少
ラッソ回帰（座標降下法） 両端のテスト（0 に押し切る・最小二乗法に戻る）、KKT 条件、ほかの言語版が 0 にした 3 列
z スコアによる外れ値の除去 性質と上限 (n − 1) ÷ √n のテスト
標準化・2 次の項・3 分割 第 2 章・第 9 章の関数をそのまま再利用

なでしこ3 版ならではの学びは次のとおりです。

  1. リッジ回帰にも突き合わせる相手がいない — 代わりに、組み立て方の違う第 7 章の正規方程式と、解が満たす式で確かめた。わざと壊した実装が落ちることも確かめた
  2. ラッソ回帰は KKT 条件で途中を確かめる — 両端のテストだけでは、途中の alpha での正しさは分からない。0 の係数と 0 でない係数の両方があることも確かめて、条件の両方の場合を試した
  3. 配列をその場で書き換えられる — 座標降下法が 1 つの繰り返しで書ける。その代わり、配列は参照で渡るので、残差にする前に複製する
  4. 命令の上限には遠いが、上限の巡回数では近づく — 報告 1 回が約 223 万命令。ラッソ回帰が 1 万巡回ると約 1.5 億命令になる
  5. 構文の語と同じ漢字を名前に入れない — 巡回 と 違反 は変数名に使えなかった

TODO リスト（この章の完了時点）:

  • リッジ回帰を自作する
  • ラッソ回帰を座標降下法で自作する
  • 係数の絶対値の合計と、0 になった特徴量名を求める
  • alpha ごとの実験結果を記録して、検証データで選ぶ
  • z スコアで外れ値を除く
  • 標準化してから 2 次の項を作り、3 つに分ける
  • 性質のテストと参照実装で数値を確かめる
  • 実データで線形回帰とリッジ回帰を比べる（実データでの実測は未確認。12.11 節の注記を参照）

次の章では、特徴量そのものを減らす主成分分析を扱います。行列の固有値分解も、この章の正規方程式と同じく自作します。

なでしこ3 版のほかの章は なでしこ3 版のトップ から辿れます。