第 13 章: 主成分分析による次元削減

13.1 はじめに

ここまでの章では、正解ラベルのあるデータで予測をしてきました。この章から扱う 教師なし学習 では、正解ラベルを使わずに、データそのものの構造を調べます。

最初に扱うのは 主成分分析（PCA） です。住宅価格のデータには 15 個の列がありますが、列どうしは互いに関係しあっていて、実際にはもっと少ない数の「軸」で大部分を説明できることがあります。主成分分析は、データのばらつきを最もよく説明する新しい軸を探し、たくさんの列を少数の軸に要約します。

Haskell 版の第 13 章 は、固有値分解を自作の Jacobi 法（ヤコビ法）と hmatrix の eigSH（LAPACK の DSYEV）の 2 通りで用意し、15 列の実データで突き合わせました。なでしこ3 版には突き合わせる相手がいません。 行列の命令は転置の 表行列交換 だけで、固有値分解はもちろん、行列の積もありません（ADR 015）。この章で書くヤコビ法が、そのまま最終実装になります。

突き合わせる相手がいないとき、「答えが正しい」ことは何で確かめればよいでしょうか。固有値分解には、答えが満たすべき性質がはっきり決まっています。

• 固有値 λ と固有ベクトル v の組は、A v = λ v を満たす
• 対称行列の固有ベクトルは、互いに直交する単位ベクトルになる
• 固有値の和は、対角成分の和（トレース）に等しい

この 3 つを満たす λ と v の組が行列の大きさだけそろえば、それは固有値分解そのものです。比べる相手の代わりに、性質をテストにします。 加えて、これまでの章と同じく、ほかの言語版の数値と、自作のデータでの Python の参照実装とも突き合わせます。

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

• ヤコビ法を配列の書き換えで書く。 Haskell 版は不変のリストを Map に写して「変更の集合」を重ねました。なでしこ3 の配列はその場で書き換えられるので、教科書どおりの形になります。その代わり、配列が参照で渡ることに気をつけます
• 並べ替えの安定性を自分で作る。 固有値の降順と係数の絶対値の降順を、同点なら元の順を残す挿入ソートで書きます
• 助詞の違う呼び出しが黙って通る。 関数を定義したときと違う助詞で呼ぶと、エラーにならず別の値が引数に入ることがありました

13.2 主成分分析とは

ばらつきが最も大きい方向を探す

データの点がばらついている空間に、新しい軸を 1 本引くことを考えます。第 1 主成分 は、その軸にデータを射影したときに、ばらつき（分散）が最も大きくなる向きです。第 2 主成分は、第 1 主成分と直交する向きのうち、ばらつきが最も大きい向きです。以下、列の数だけ同じように決めていきます。

分散共分散行列の固有ベクトル

この「ばらつきが最大になる向き」は、分散共分散行列の固有ベクトル として求まります。固有値が、その向きでのばらつきの大きさ（分散）です。

分散共分散行列は、対角に各列の分散、非対角に 2 列ずつの共分散を並べた正方行列で、必ず対称になります。対称行列の固有値分解では、実数の固有値と、互いに直交する固有ベクトルが得られます。13.1 節で挙げた性質のうち 2 つ目は、この事実から来ています。

寄与率

固有値をすべて足すと、データ全体のばらつきの合計になります。ある主成分の固有値を合計で割った値が 寄与率 で、「この軸だけで全体の何割を説明できるか」を表します。寄与率を大きいほうから足していったものが 累積寄与率 です。累積寄与率が 0.8 に届くまでの主成分を採る、といった目安の決め方をします。

13.3 題材とデータ

ボストンの住宅価格（Boston.csv、100 件）を使います。列は CRIME・ZN・INDUS・CHAS・NOX・RM・AGE・DIS・RAD・TAX・PTRATIO・B・LSTAT・PRICE の 14 個で、CRIME はカテゴリ値（high・low・very_low）、NOX と RAD に欠損値があります。

この章では正解ラベルを使わないので、PRICE も含めたすべての列を主成分分析にかけます。前処理は Haskell 版と同じです。

1. CRIME をダミー変数にする。カテゴリを文字の順に並べて最初の high を基準として落とし、CRIME_low・CRIME_very_low の 2 列を末尾に加える
2. 欠損値を列の平均値で補完する（第 2 章の 列平均算出・欠損補完 をそのまま呼ぶ）
3. すべての列を平均 0・標準偏差 1 にそろえる（標準偏差は件数で割る。0 なら 1 とみなす）

これで 13 列と 2 列を合わせた 15 列 になります。主成分分析は「ばらつきの大きさ」を見るので、単位の大きい列（TAX は数百、NOX は 0.5 前後）がそのままでは支配的になります。標準化はそれを防ぐためです。

13.4 TODO リストの作成

TODO リスト:

• 列ごとの平均を求める
• 各列から平均を引く（中心化）
• 分散共分散行列を求める
• ヤコビ法で固有値分解を自作する
• 対称でない行列・正方でない行列を断る
• A v = λ v・固有ベクトルの正規直交・固有値の和とトレースの一致を確かめる
• 固有ベクトルの符号をそろえる
• 寄与率と累積寄与率を求める
• データを主成分の向きに射影する
• 累積寄与率がしきい値に届く主成分の数を求める
• 主成分への影響が大きい列を求める
• Boston を前処理して実データで要約する

Haskell 版の TODO リストにあった「hmatrix の eigSH の振る舞いを学習用テストで確かめる」「自作と hmatrix を 15 列の実データで突き合わせる」が、なでしこ3 版にはありません。代わりに、固有値分解の性質を小さな行列と 15 列の行列で確かめます。

13.5 分散共分散行列を求める

Red: 分散と共分散

行列は「長さのそろった数値の配列」を並べた配列で表します。2 列目が 1 列目のちょうど 2 倍になる、完全に相関する 3 件を架空のデータにします。分散は 1 と 4、共分散は 2 になるはずです。

相関二列 = [[1, 2], [2, 4], [3, 6]]

「列ごとの平均を求める」で(相関二列の主成分列平均)と[2, 4]が一致確認
「各列から平均を引く」で(相関二列を[2, 4]で主成分中心化)と[[-1, -2], [0, 0], [1, 2]]が一致確認
「中心化しても元の行列は変わらない」で相関二列と[[1, 2], [2, 4], [3, 6]]が一致確認
「分散共分散行列は対称で、対角に件数から 1 を引いた数で割った分散が並ぶ」で((相関二列の分散共分散行列算出)の検査用平坦化)と[1, 2, 2, 4]を0.000000001で近似確認

テスト補助の 近似確認 は数値か数値の配列を受け取るので、行列はテストの中で 1 本の配列に並べてから比べます（検査用平坦化）。

3 つ目のテストは、なでしこ3 に固有のものです。配列は参照で渡るので、中心化で元の行列の値を引いてしまうと、呼び出し元の行列が書き換わります。 新しい配列を作って返すことを、元の行列が変わっていないことで確かめます。

Green: 中心化してから内積を取る

●(データ行列の)分散共分散行列算出とは
    件数 = データ行列の要素数
    もし、件数 < 2ならば
        「主成分分析には 2 件以上のデータが必要です（{件数} 件）」のエラー発生
    ここまで
    中心化行列 = データ行列を(データ行列の主成分列平均)で主成分中心化
    # 列を 1 本ずつ取り出して、列どうしの内積を取る。
    列一覧 = []
    Jを 0から (中心化行列[0]の要素数) - 1まで繰り返す
        列 = []
        Iを 0から 件数 - 1まで繰り返す
            列に中心化行列[I][J]を配列追加
        ここまで
        列一覧に列を配列追加
    ここまで
    結果 = []
    Aを 0から (列一覧の要素数) - 1まで繰り返す
        行 = []
        Bを 0から (列一覧の要素数) - 1まで繰り返す
            行に((列一覧[A]と列一覧[B]の主成分内積) ÷ (件数 - 1))を配列追加
        ここまで
        結果に行を配列追加
    ここまで
    結果で戻る
ここまで

Haskell 版が transpose と内包表記の 2 行で書いたところが、三重の繰り返しになります。列を取り出す部分は転置にあたります。何をしているかが読めるように、転置も書き下しました。

割る数を n−1 にするか n にするかは、この章の結論に影響しません。 寄与率は「固有値 ÷ 固有値の合計」なので、行列全体を定数倍しても変わらないからです。それでも n−1 を選んだのは、ほかの言語版と分散（固有値）そのものの値をそろえるためです。

内積は左から順に足します。足す順番が変わると下の桁が動くので、Haskell 版の foldl' と同じ順にしています。

●(左ベクトルと右ベクトルの)主成分内積とは
    和 = 0
    Iを 0から (左ベクトルの要素数) - 1まで繰り返す
        和 = 和 + 左ベクトル[I] × 右ベクトル[I]
    ここまで
    和で戻る
ここまで

関数の名前には、どれも章の題材の名詞（主成分・固有・ボストン）を含めています。なでしこ3 では取り込んだファイルをまたいで関数名が 1 つの空間に並ぶので、内積 や 列平均 のような汎用の名前は、線形代数を自作するほかの章とぶつかるからです。

13.6 固有値分解を自作する

対称かどうかは自分で確かめる

ヤコビ法は対称行列にしか使えません。Haskell 版は hmatrix の trustSym が「対称であることを確かめずに信じ、上三角だけを見る」ことを学習用テストで確かめ、手前に自分の検査を置きました。なでしこ3 版には信じてしまうライブラリがありませんが、自作のヤコビ法も対称であることを前提にしています。回転を決めるのは上三角の (P, Q) 成分で、回転したあとの値を (K, P) と (P, K) の両方に書き込むからです。前提を破った入力は、手前で断ります。

非対称分解 = 関数()
    [[1, 2], [3, 4]]のヤコビ固有値分解
ここまで
「対称でない行列は固有値分解できない」で非対称分解を「固有値分解できません（対称行列ではありません）」でエラー確認
非正方分解 = 関数()
    [[1, 2, 3], [2, 4, 5]]のヤコビ固有値分解
ここまで
「正方行列でなければ固有値分解できない」で非正方分解を「正方行列ではありません: 2 行」でエラー確認

エラーの文言は Haskell 版と同じにしています。

●(対象行列の)固有値分解前検査とは
    次数 = 対象行列の要素数
    Iを 0から 次数 - 1まで繰り返す
        もし、(対象行列[I]の要素数) ≠ 次数ならば
            「正方行列ではありません: {次数} 行」のエラー発生
        ここまで
    ここまで
    Iを 0から 次数 - 1まで繰り返す
        Jを 0から 次数 - 1まで繰り返す
            もし、対象行列[I][J] ≠ 対象行列[J][I]ならば
                「固有値分解できません（対称行列ではありません）」のエラー発生
            ここまで
        ここまで
    ここまで
ここまで

Haskell 版は m /= transpose m の 1 行でした。なでしこ3 でも 一致確認 のように JSON にして比べれば 1 行になりますが、数値の比較を文字列の比較に置き換えることになるので、成分ごとに比べています。

ヤコビ法

ヤコビ法は、対称行列の非対角成分を 2 次元の回転で 1 つずつ 0 にしていく手続きです。1 つ 0 にすると別のところがわずかに復活するので、全体が対角行列に十分近づくまで何巡もします。回転行列を掛け合わせたものが、そのまま固有ベクトルになります。固有値と固有ベクトルが同時に出るので、「固有値を求めてから、それを使って固有ベクトルを探す」という 2 段構えになりません。 PHP 版 でライブラリが 15 列の実データで落ちたのは、その 2 段目でした。

1 回の回転は次のとおりです。

●(作業行列と固有列行列をPとQで)ヤコビ回転とは
    APQ = 作業行列[P][Q]
    もし、APQ = 0ならば
        戻る
    ここまで
    APP = 作業行列[P][P]
    AQQ = 作業行列[Q][Q]
    角度比 = (AQQ - APP) ÷ (2 × APQ)
    回転符号 = 1
    もし、角度比 < 0ならば
        回転符号 = -1
    ここまで
    正接 = 回転符号 ÷ ((角度比の絶対値) + ((角度比 × 角度比 + 1)の平方根))
    余弦 = 1 ÷ ((正接 × 正接 + 1)の平方根)
    正弦 = 正接 × 余弦
    次数 = 作業行列の要素数
    Kを 0から 次数 - 1まで繰り返す
        もし、(K ≠ P) かつ (K ≠ Q)ならば
            KP値 = 余弦 × 作業行列[K][P] - 正弦 × 作業行列[K][Q]
            KQ値 = 正弦 × 作業行列[K][P] + 余弦 × 作業行列[K][Q]
            作業行列[K][P] = KP値
            作業行列[P][K] = KP値
            作業行列[K][Q] = KQ値
            作業行列[Q][K] = KQ値
        ここまで
    ここまで
    作業行列[P][P] = APP - 正接 × APQ
    作業行列[Q][Q] = AQQ + 正接 × APQ
    作業行列[P][Q] = 0
    作業行列[Q][P] = 0
    Kを 0から 次数 - 1まで繰り返す
        VKP値 = 余弦 × 固有列行列[K][P] - 正弦 × 固有列行列[K][Q]
        VKQ値 = 正弦 × 固有列行列[K][P] + 余弦 × 固有列行列[K][Q]
        固有列行列[K][P] = VKP値
        固有列行列[K][Q] = VKQ値
    ここまで
ここまで

回転符号 ÷ ((角度比の絶対値) + …の平方根) は、2 次方程式の根のうち絶対値の小さいほうを、引き算を通さずに求める書き方です。素直に解の公式を書くと、角度比が大きいときに近い数どうしの引き算になって桁落ちします。式は Haskell 版と演算の順までそろえました。

Haskell 版は、更新する成分の Map を作って元の行列に重ねました。なでしこ3 の配列はその場で書き換えられるので、作業行列[K][P] = KP値 と教科書どおりに書けます。気をつけるのは 2 点です。

• 新しい値を求めてから書き換える。 KQ値 の計算には書き換える前の 作業行列[K][P] が要るので、2 つの値を先に求めてから代入します。固有ベクトルの行列でも同じです
• 呼ぶ側で複製したものを渡す。 配列は参照で渡るので、ヤコビ回転 は渡された行列そのものを書き換えます。ヤコビ固有値分解 は最初に行列を 1 成分ずつ写した作業用の行列を作り、テストの「分解しても元の行列は変わらない」で守っています

巡回は、非対角成分の 2 乗和が打ち切りの値を下回るまで続けます。

    巡回数 = 0
    (巡回数 < ヤコビ最大巡回数)の間
        もし、(作業行列の非対角二乗和) < ヤコビ打切値ならば
            抜ける
        ここまで
        Pを 0から 次数 - 2まで繰り返す
            Qを (P + 1)から 次数 - 1まで繰り返す
                作業行列と固有列行列をPとQでヤコビ回転
            ここまで
        ここまで
        巡回数 = 巡回数 + 1
    ここまで

打ち切りの条件は Haskell 版の自作ヤコビ法と同じにしました。

設定	値	意味
ヤコビ最大巡回数	100	掃き出し（上三角のすべての組を 1 回ずつ回転する 1 巡）の回数の上限
ヤコビ打切値	1e-30	非対角成分の 2 乗和がこれを下回ったら対角行列とみなす

繰り返しの条件はかっこで囲みます。第 3 章で踏んだとおり、囲まないと文法の検査は通るのに実行で「『while』はまだ実行に対応していません。」になります。Qを (P + 1)から のかっこも同じ理由です。

Haskell 版は「再帰の引数に残り回数を持たせる」形で上限を書きました。なでしこ3 版は数を数えて 抜ける で抜ける、命令型の形です。そのぶん掃き出しの回数を数えるのが簡単なので、回数と打ち切ったときの 2 乗和も結果の辞書に入れて返します（13.10 節で使います）。

大きい順に並べる

eigSH は固有値を大きい順に返しました。自作では回転を終えた対角成分の順なので、並べ替えが要ります。

なでしこ3 の 表数値ソート は安定ですが、昇順で、元の表を書き換えます。固有値と固有ベクトルの組を降順に、同じ固有値なら元の順を残して並べたいので、挿入ソートを書きました。

●(固有値一覧と固有ベクトル一覧の)固有値降順並替とは
    順番 = []
    Iを 0から (固有値一覧の要素数) - 1まで繰り返す
        位置 = 順番の要素数
        (位置 > 0)の間
            もし、固有値一覧[I] > 固有値一覧[順番[位置 - 1]]ならば
                位置 = 位置 - 1
            違えば
                抜ける
            ここまで
        ここまで
        順番の位置にIを配列挿入
    ここまで
    …（順番どおりに固有値と固有ベクトルを並べ直す）
ここまで

より大きいときだけ前へ動かす ので、同じ固有値は後から来たものが後ろに残ります。Haskell 版の sortBy (comparing (Down . fst)) が安定であることに寄りかかった部分を、自分で書いたことになります。

同値分解 = [[2, 0], [0, 2]]のヤコビ固有値分解
「同じ固有値なら元の順を残す」で(同値分解["固有ベクトル"])と[[1, 0], [0, 1]]が一致確認

> を ≧ にすると、このテストが落ちることを確かめました。

性質をテストにする

2×2 の小さな行列から始めます。[[2, 1], [1, 2]] の固有値は 3 と 1 で、1 回の回転で対角になります。

二次分解 = [[2, 1], [1, 2]]のヤコビ固有値分解
「対称行列の固有値を大きい順に返す」で(二次分解["固有値"])と[3, 1]を0.000000001で近似確認
「固有ベクトルは単位ベクトルで、互いに直交する」で((二次分解["固有ベクトル"])の検査用グラム行列の検査用平坦化)と(2の検査用単位行列)を0.000000001で近似確認
「一度の回転で対角になれば一巡で止まる」で(二次分解["掃き出し回数"])と1が一致確認
対角分解 = [[3, 0], [0, 5]]のヤコビ固有値分解
「対角行列はそのまま固有値になり、大きい順に並べ替える」で(対角分解["固有値"])と[5, 3]が一致確認
「対角行列なら一巡もしない」で(対角分解["掃き出し回数"])と0が一致確認

検査用グラム行列 は、固有ベクトルどうしの内積をすべて並べた行列です。固有ベクトルが互いに直交する単位ベクトルなら、単位行列になります。

答えを手で求めにくい 3×3 と 5×5 では、13.1 節の 3 つの性質を確かめます。

対称五次 = [[6, 2, -1, 0, 3], [2, 5, 1, -2, 0], [-1, 1, 4, 1, 1], [0, -2, 1, 3, -1], [3, 0, 1, -1, 7]]
五次分解 = 対称五次のヤコビ固有値分解
「五次で A v = λ v が成り立つ」で(対称五次と五次分解の固有方程式残差)と0を0.000000001で近似確認
「五次の固有ベクトルが正規直交する」で((五次分解["固有ベクトル"])の検査用グラム行列の検査用平坦化)と(5の検査用単位行列)を0.000000001で近似確認
「五次の固有値の和がトレースに等しい」で((五次分解["固有値"])の検査用総和)と(対称五次の検査用対角和)を0.000000001で近似確認
「非対角成分が打ち切りのしきい値を下回るまで回す」で(((五次分解["非対角二乗和"]) < 0.000000000000000000000000000001) かつ ((五次分解["掃き出し回数"]) ≦ 100))と真が一致確認

固有方程式残差 は、すべての組について A v − λ v の成分の絶対値の最大を返すテストの中の関数です。行列とベクトルの積もテストの中で書いています。

性質のテストは、書いたときから通ったので、実装をわざと壊して落ちることを確かめました。

壊し方	落ちたテスト
正弦の符号を逆にする	A v = λ v（3×3・5×5）、主成分の向き、射影、報告
固有ベクトルの行列の Q 列を回転しない	正規直交（2×2・3×3・5×5）、A v = λ v、主成分の向き、射影、報告
打ち切りの値を 1 にする	A v = λ v（3×3・5×5）、打ち切りのしきい値、報告
降順の比較を ≧ にする	同じ固有値なら元の順を残す

正弦の符号を逆にしても、固有値の和とトレースの一致は落ちませんでした。 回転の向きを間違えても、回転は対角成分の和を変えないからです。トレースの一致は必要条件でしかなく、A v = λ v と正規直交があってはじめて固有値分解だと言えます。3 つの性質を並べたのはこのためです。

符号の規則を決める

固有ベクトルは、すべての要素の符号を反転しても同じ向きを表します。どちらが出てくるかは実装しだいなので、この連載では「絶対値が最大の要素を正にする」で統一しています。

●(ベクトル一覧の)固有ベクトル符号整列とは
    結果 = []
    Iを 0から (ベクトル一覧の要素数) - 1まで繰り返す
        ベクトル = ベクトル一覧[I]
        最大要素 = 0
        Jを 0から (ベクトルの要素数) - 1まで繰り返す
            もし、(ベクトル[J]の絶対値) > (最大要素の絶対値)ならば
                最大要素 = ベクトル[J]
            ここまで
        ここまで
        新ベクトル = []
        Jを 0から (ベクトルの要素数) - 1まで繰り返す
            もし、最大要素 < 0ならば
                新ベクトルに(0 - ベクトル[J])を配列追加
            違えば
                新ベクトルにベクトル[J]を配列追加
            ここまで
        ここまで
        結果に新ベクトルを配列追加
    ここまで
    結果で戻る
ここまで

絶対値が同じ要素が並んだときにどちらを見るかで結果が変わるので、≧ ではなく厳密な > にします。≧ にすると「絶対値が同じなら前の要素を見る」のテストが落ちることを確かめました。

符号を反転するところは 0 - ベクトル[J] と書いています。なでしこ3 では -(式) のようにかっこの前に負号を付けると「(...)の解析エラー」になります。 -1 のような数の直前の負号は書けます。

13.7 主成分を求める

固有値分解を呼んで寄与率を出す

●(データ行列を主成分数で)主成分分析学習とは
    共分散結果 = データ行列の分散共分散行列算出
    分解結果 = 共分散結果のヤコビ固有値分解
    固有値一覧 = 分解結果["固有値"]
    列数 = 固有値一覧の要素数
    もし、(主成分数 < 1) または (主成分数 > 列数)ならば
        「主成分の数は 1 以上 {列数} 以下にしてください: {主成分数}」のエラー発生
    ここまで
    総和 = 0
    Iを 0から 列数 - 1まで繰り返す
        総和 = 総和 + 固有値一覧[I]
    ここまで
    …（先頭から主成分数だけ、固有ベクトル・固有値・固有値 ÷ 総和を集める）
    {"平均": データ行列の主成分列平均, "主成分": 選択ベクトルの固有ベクトル符号整列, "分散": 分散一覧, "寄与率": 寄与率一覧, "掃き出し回数": 分解結果["掃き出し回数"]}で戻る
ここまで

Haskell 版は fitWith に固有値分解を引数で渡し、自作と hmatrix を同じ手順で回しました。なでしこ3 でも関数を値として渡せますが（B76 で確かめました）、差し替える相手がいないので、引数にはしていません。 学習の結果は、第 3 章の決定木と同じく辞書で返します。

完全に相関する 2 列

完全に相関する 2 列なら、1 本の軸ですべてを説明できます。第 1 主成分の寄与率が 1、第 2 主成分が 0 になり、第 1 主成分は 2 列を 1 対 2 で混ぜる向き（1/√5, 2/√5）になるはずです。

相関モデル = 相関二列を2で主成分分析学習
「完全に相関する 2 列なら第 1 主成分の寄与率が 1 になる」で(相関モデル["寄与率"])と[1, 0]を0.000000001で近似確認
「主成分の分散は固有値そのもの」で(相関モデル["分散"])と[5, 0]を0.000000001で近似確認
「第 1 主成分は 2 列を 1 対 2 で混ぜる向き」で((相関モデル["主成分"])[0])と[1 ÷ (5の平方根), 2 ÷ (5の平方根)]を0.000000001で近似確認

Haskell 版は「hmatrix で学習しても同じ寄与率になる」を比べました。なでしこ3 版では、答えが手で求まる例で主成分の向きそのものを確かめています。

13.8 データを主成分の向きに射影する

●(モデルでデータ行列を)主成分射影とは
    中心化行列 = データ行列を(モデル["平均"])で主成分中心化
    主成分一覧 = モデル["主成分"]
    結果 = []
    Iを 0から (中心化行列の要素数) - 1まで繰り返す
        行 = []
        Kを 0から (主成分一覧の要素数) - 1まで繰り返す
            行に(中心化行列[I]と主成分一覧[K]の主成分内積)を配列追加
        ここまで
        結果に行を配列追加
    ここまで
    結果で戻る
ここまで

完全に相関する 3 件を 1 次元に落とすと、真ん中の点が 0、前後の点が同じ大きさで符号だけ逆になります。符号は 13.6 節の規則で決まるので、値まで書けます。

射影結果 = (相関二列を1で主成分分析学習)で相関二列を主成分射影
「主成分の向きに射影すると、完全に相関する 2 列が 1 列で表せる」で(射影結果の検査用平坦化)と[0 - (5の平方根), 0, 5の平方根]を0.000000001で近似確認

13.9 必要な主成分の数と、影響が大きい列

累積寄与率がしきい値に届くまでの数

●(寄与率一覧を閾値で)必要主成分数算出とは
    累積 = 0
    Iを 0から (寄与率一覧の要素数) - 1まで繰り返す
        累積 = 累積 + 寄与率一覧[I]
        もし、累積 ≧ 閾値ならば
            (I + 1)で戻る
        ここまで
    ここまで
    寄与率一覧の要素数で戻る
ここまで

見つけたところで 戻る で抜け、しきい値に届かないまま終わったら、すべての主成分の数を返します。Haskell 版の scanl1 と遅延評価が担った「見つけたら止める」を、繰り返しの途中の 戻る で書いた形です。

名前の「しきい値」は 閾値 と漢字で書きました。ひらがなを含む名前は助詞で割られるおそれがあるからです（第 1 章）。

主成分への影響が大きい列

主成分の「意味」は、係数（ローディング）の大きい列から読みます。係数の絶対値の降順に、同じなら元の列の順を残して並べるので、固有値と同じく挿入ソートで書きます。

「係数の絶対値が大きい順に列名と係数を返す」で([0.1, -0.9, 0.5]と[「a」, 「b」, 「c」]を2で主成分影響列抽出)と[[「b」, -0.9], [「c」, 0.5]]が一致確認
「係数の絶対値が同じなら元の列の順を残す」で([0.5, -0.5, 0.1]と[「a」, 「b」, 「c」]を2で主成分影響列抽出)と[[「a」, 0.5], [「b」, -0.5]]が一致確認

列名と係数の組は、Haskell 版のタプルの代わりに 2 要素の配列で返します。

13.10 Boston を前処理して突き合わせる

ダミー変数化と標準化

Haskell 版は第 2 章と第 9 章の関数をそのまま呼びました。なでしこ3 版の第 9 章（特徴量エンジニアリング）はこの章より後に書くので、この章の時点ではダミー変数化と標準化を章の中に書いています。規則は Haskell 版の第 9 章と同じです。

●(CSV文字列の)ボストン行列読込とは
    読込表 = CSV文字列を表読込
    変換表 = 読込表を主成分カテゴリ列でボストンダミー変数化
    列名一覧 = 変換表["列名一覧"]
    行一覧 = 変換表["行一覧"]
    平均一覧 = 行一覧と列名一覧の列平均算出
    補完済 = 行一覧と列名一覧を平均一覧で欠損補完
    標準化済 = 補完済と列名一覧のボストン標準化
    {"列名一覧": 列名一覧, "行列": 標準化済と列名一覧の主成分行列変換}で戻る
ここまで

欠損値の平均と補完は第 2 章の 列平均算出・欠損補完 をそのまま呼びます。ダミー変数の列には数値の 1 と 0 を入れるので、第 2 章の関数がそのまま使えます。

カテゴリは空欄を除いて 配列ソート で並べ、最初の 1 つを落とします。配列ソート は文字列として比べる安定な並べ替えで、処理系のソースで Go の文字列の比較（バイトの順）であることを確かめました。Haskell の Text の比較と同じ順になります。

架空の 4 件（CRIME が 3 種類、RM に欠損値が 1 件）で、列名と標準化の結果を確かめます。

ボストン風 = 「CRIME,ZN,RM,PRICE{改行}low,18,6.5,24{改行}high,0,,21{改行}very_low,12,7.1,35{改行}low,0,5.9,18{改行}」
風行列 = ボストン風のボストン行列読込
「CRIME を除き、ダミー変数を末尾に加える」で(風行列["列名一覧"])と[「ZN」, 「RM」, 「PRICE」, 「CRIME_low」, 「CRIME_very_low」]が一致確認
「欠損値を補完してから各列の平均を 0 にそろえる」で((風行列["行列"])の主成分列平均)と[0, 0, 0, 0, 0]を0.000000001で近似確認

列の順は配列で持ち回る

特徴量は列名から数値への辞書です。なでしこ3 の辞書はキーを入れた順を保ちますが（第 2 章）、行列にするときは列名の配列の順で並べ、辞書の順には頼りません。 Haskell 版の Map が鍵の順に並び替えてしまうのと事情は逆ですが、「列の順を 1 か所で決める」ことは同じです。

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

実データと同じ形の「ボストン風」のデータ（100 行・14 列、CRIME が 3 種類、NOX に 3 か所・RAD に 4 か所の欠損。値の範囲を実データに似せた架空のデータで、配布データではありません）を作り、Python で書いた参照実装と突き合わせました。参照実装は、Haskell 版と同じ規則（ダミー変数・平均値での補完・件数で割る標準化・n−1 で割る分散共分散行列・同じ順の回転・同じ打ち切り）を、演算の順までそろえて書いたものです。

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

データ件数: 100, 列数: 15
寄与率: PC1 0.5225, PC2 0.1217, PC3 0.0912, PC4 0.0661
累積寄与率が 0.8 に届く主成分の数: 4（累積寄与率 0.8014）
第 1 主成分で影響の大きい列: NOX 0.328, TAX 0.323, DIS -0.322
第 2 主成分で影響の大きい列: RM 0.567, PRICE 0.482, LSTAT -0.345
ヤコビ法の掃き出し: 8 巡

報告は参照実装と 1 文字も違いませんでした。 さらに 15 個の固有値を最短の 10 進表記で並べたものも、打ち切ったときの非対角成分の 2 乗和（3.2702036702512856e-63）も、参照実装と一致しました。float64 の四則演算と平方根はどちらも IEEE 754 のとおりに丸めるので、演算の順がそろえば最後のビットまで同じになります。これは自作のデータでの値で、実データの値ではありません。

同じデータで、15 列の分散共分散行列にも 13.6 節の 3 つの性質が成り立つこと（A v = λ v の残差・正規直交のずれ・固有値の和とトレースの差がいずれも 1e-9 未満）を、実データ用のテストの分岐を自作のデータで走らせて確かめました。

打ち切りまでの回数と時間

打ち切りの条件（上限 100 巡、非対角成分の 2 乗和 1e-30）で、掃き出しの回数と時間を測りました。

行列	掃き出しの回数	打ち切ったときの非対角成分の 2 乗和
2×2（[[2, 1], [1, 2]]）	1	0
3×3（テストの 対称三次）	4	1.46e-109
5×5（テストの 対称五次）	5	7.17e-63
7×7（13.12 節の自作の小さなデータ、10 件）	6	5.71e-63
15×15（ボストン風のデータ、100 件）	8	3.27e-63

どの値も Python の参照実装と同じでした。

15×15 でも 8 巡で打ち切りに届き、上限の 100 巡には遠く届きません。 ヤコビ法は収束が進むと非対角成分がほぼ 2 乗の速さで小さくなるので、1e-30 という厳しい値でも数巡で下回ります。

時間は外から date で測りました。なでしこ3 の システム時間ミリ秒 は実行中に進まないからです（13.13 節）。

計測	時間
章を取り込むだけ（処理系の起動と構文解析）	約 0.3 秒
章の実行（読み込みから報告まで、起動を含む）	約 0.35〜0.45 秒
15×15 の固有値分解 1 回（20 回の差分から）	約 0.04 秒

15×15 の掃き出し 8 巡は 105 組 × 8 = 840 回の回転です。インタプリタで配列の添字を読み書きしても、1 回の固有値分解は 0.04 秒ほどでした。時間のほとんどは処理系の起動が占めています。

データが無ければ保留にする

実データのテストは、Haskell 版と同じ期待値で書き、データが無ければ保留にします。

./bin/gonako main.nako3 chapter13

期待する出力の先頭 5 行は次のとおりです（6 行目は掃き出しの回数です）。

データ件数: 100, 列数: 15
寄与率: PC1 0.4110, PC2 0.1448, PC3 0.1019, PC4 0.0645, PC5 0.0623, PC6 0.0581
累積寄与率が 0.8 に届く主成分の数: 6（累積寄与率 0.8427）
第 1 主成分で影響の大きい列: INDUS 0.359, NOX 0.350, TAX 0.328
第 2 主成分で影響の大きい列: PRICE 0.444, CRIME_low 0.423, RM 0.405

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

この章は学習データの無い環境で書いたため、上の出力はまだなでしこ3 版の実装で実測していません。実データでの掃き出しの回数も未測定です。学習データを配置した環境で npx gulp apps:check:nadesiko3 を実行し、保留が消えて通ることを確かめてから、この注記を外します。自作のボストン風データでは、Python の参照実装と報告が 1 文字も違わず、固有値も最後の桁まで一致することを確かめています。自作のデータを ML_DATA_DIR で指して tools/check.sh を実行すると、実データのテストは値の不一致（上の期待値と自作のデータの値の違い）だけで落ち、書き間違いのエラーでは止まりませんでした。

Haskell 版の報告の 6 行目は「hmatrix の eigSH との差」でした。なでしこ3 版は差を測る相手がいないので、代わりに掃き出しの回数を出しています。実装によって動く行なので、テストでは Haskell 版と同じく先頭 5 行だけを比べます。

結果を読む

期待する出力（ほかの言語版の実測）から読めることは次のとおりです。

• 第 1 主成分だけで 4 割を説明する: 15 列のばらつきのうち 41.1% を、1 本の軸で説明できます
• 15 列を 6 本の軸に要約できる: 累積寄与率 0.8 を目安にすると、6 つの主成分で全体の 84.3% を説明できます
• 第 1 主成分は「産業化・都市化」の軸と読める: 非小売業の土地の割合（INDUS）、窒素酸化物の濃度（NOX）、固定資産税率（TAX）が同じ向きに強く効いています
• 第 2 主成分は「住環境のよさ」の軸と読める: 住宅価格（PRICE）、犯罪率が低い地区（CRIME_low）、部屋数（RM）が同じ向きに効いています

主成分の「意味」は計算で決まるものではなく、係数を見た人間の解釈です。符号の向きも 13.6 節の規則で決めたものなので、係数の符号とあわせて読む必要があります。

13.11 可視化について

Notebook による探索と可視化の節は設けません。主成分の散布図、累積寄与率の折れ線（スクリープロット）、ローディングのヒートマップは、Python 版の第 13 章 と Kotlin 版の第 13 章 の「Notebook による探索と可視化」の節を参照してください。

13.12 突き合わせる相手がいないということ

Haskell 版は「何が置き換えられて、何が置き換えられないのか」を表にし、置き換えられるのは固有値分解だけだと結論しました。なでしこ3 版では、その固有値分解も含めて、置き換える相手がいません。

代わりに、次の 4 つで支えています。

支え	確かめられること
固有値分解の性質（2×2・3×3・5×5、データがあれば 15×15）	A v = λ v、固有ベクトルの正規直交、固有値の和とトレースの一致。3 つそろえば固有値分解そのもの
手で答えが求まる例	完全に相関する 2 列の寄与率・主成分の向き・射影、対角行列、同じ固有値の並び
自作データと参照実装	10 件（下）と 100 件の自作データで、報告の文字列と固有値が最後の桁まで一致すること
ほかの言語版の数値	実データでの寄与率と主成分の係数（データのある環境で確かめる）

報告を組み立てる関数は、データの無い環境でも走るよう、自作の 10 件のデータで直接呼ぶテストも書きました。期待値は Python の参照実装で求めたものです。

データ件数: 10, 列数: 7
寄与率: PC1 0.6751, PC2 0.2110
累積寄与率が 0.8 に届く主成分の数: 2（累積寄与率 0.8861）
第 1 主成分で影響の大きい列: RM 0.450, PRICE 0.449, NOX -0.447
第 2 主成分で影響の大きい列: CRIME_low 0.808, CRIME_very_low -0.406, RAD -0.379
ヤコビ法の掃き出し: 6 巡

性質のテストは、ライブラリとの突き合わせより強い面もあります。 ライブラリと一致しても、両方が同じ前提で同じ間違いをしている可能性は残ります。A v = λ v は定義そのものなので、その可能性がありません。弱い面は、固有値分解が正しくても、並べ方と符号の規則がほかの言語版と同じかどうかは分からない ことです。そこはほかの言語版の数値と、同じ規則で書いた参照実装で確かめます。

13.13 品質チェック

./tools/check.sh

文法・整形・テスト・テストが最後の報告まで届いたか・関数の網羅の 5 つをまとめて確かめます（第 6 章）。この章のテストは 44 件（保留 3 件）で、関数の網羅は次のとおりでした（この章を書いた作業ツリーでの値です）。

関数の網羅: 70 件中 70 件（直接 53 件、間接 17 件）

この章で新しく出会ったなでしこ3 の振る舞いは 3 つです。

関数を定義したときと違う助詞で呼ぶと、エラーにならないことがある。 テストの補助の関数を ●(行列を)検査用平坦化とは と「を」で定義し、(…)の検査用平坦化 と「の」で呼んだところ、文法の検査は通り、行列ではなく直前に作った JSON の文字列が引数に入って、1 文字ずつに分解された配列が返りました。

失敗: 分散共分散行列は対称で、対角に件数から 1 を引いた数で割った分散が並ぶ
  期待: [1,2,2,4]（誤差 1e-9）
  実際: ["[","[","1",",","2","]",",","[","2",",","4","]",",","[","3",",","6","]","]"]

「の」は定義に無い助詞なので引数として読まれず、足りない引数が直前に計算した値で埋められたと考えられます。第 5 章で見た「命令の引数が足りないと、直前の値で埋められる」と同じ仕組みが、自分で定義した関数でも起きます。同じ間違いを関数の定義の中で書いたとき（●(対称行列を)ヤコビ固有値分解とは を別の関数の中から「の」で呼んだとき）は、文法の検査が「未解決の単語があります」で止めました。止まるかどうかは書いた場所しだいです。 定義と呼び出しの助詞をそろえ、テストで結果を確かめるしかありません。

組み込みの 符号 に代入できない。 回転の向きを表す変数を 符号 と名付けたら「関数『符号』に代入できません」で止まりました。回転符号 にしています。

システム時間ミリ秒 は実行中に進まない。 固有値分解の前後で読んだ差が 0 ミリ秒でした。処理系のソース（internal/event/loop.go）を読むと、時刻はイベントループが持つ値で、タイマーの予定時刻へ飛ぶときにしか進みません。計算の時間は外から測ります。

13.14 まとめ

この章では、主成分分析をなでしこ3 の TDD で自作し、固有値分解をヤコビ法で書きました。

1. 固有値分解も自作が最終実装になる — 行列の命令が無いので、突き合わせる相手がいない
2. 性質をテストにする — A v = λ v、固有ベクトルの正規直交、固有値の和とトレースの一致。トレースの一致だけでは回転の向きの誤りを見逃すので、3 つを並べる。実装をわざと壊して、それぞれが落ちることを確かめた
3. 打ち切りの条件は Haskell 版と同じ — 上限 100 巡、非対角成分の 2 乗和 1e-30。自作の 100 件・15 列のデータで 8 巡、1 回の固有値分解は約 0.04 秒
4. 参照実装と最後の桁まで一致した — 演算の順をそろえた Python の参照実装と、報告の文字列も 15 個の固有値も同じになった。実データの値（寄与率 PC1 0.4110、必要な主成分 6 個、累積寄与率 0.8427）はデータのある環境で確かめる

なでしこ3 版ならではの学びもありました。

• 配列をその場で書き換えられる — ヤコビ法が教科書どおりの形で書ける。その代わり、配列は参照で渡るので、作業用に複製してから回し、元の行列が変わらないことをテストで守る
• 安定な降順の並べ替えを自分で書く — 「より大きいときだけ前へ動かす」挿入ソートで、同じ値なら元の順を残す
• 助詞をそろえる — 定義と違う助詞で呼ぶと、エラーにならず別の値が渡ることがある
• -(式) は書けない — 符号の反転は 0 - 式 と書く

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

• 列ごとの平均を求める
• 各列から平均を引く（中心化）
• 分散共分散行列を求める
• ヤコビ法で固有値分解を自作する
• 固有ベクトルの符号をそろえる
• 寄与率と累積寄与率を求める
• データを主成分の向きに射影する
• 累積寄与率がしきい値に届く主成分の数を求める
• 主成分への影響が大きい列を求める
• Boston を前処理して実データで要約する（実データでの実測は未確認。13.10 節の注記を参照）

次の章では、もう 1 つの教師なし学習である K-means によるクラスタリングを扱います。

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