AdvML2024_ex02notebookA.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/takatakamanbou/AdvML/blob/2024/AdvML2024_ex02notebookA.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "2UzDNLmVS3po"
      },
      "source": [
        "# AdvML ex02notebookA\n",
        "\n",
        "<img width=72 src=\"https://www-tlab.math.ryukoku.ac.jp/~takataka/course/AdvML/AdvML-logo.png\"> [この授業のウェブページ](https://www-tlab.math.ryukoku.ac.jp/wiki/?AdvML)\n",
        "\n",
        "※ この notebook は，授業時間中の解説や板書と併用することを想定して作っていますので，説明が不十分なところが多々あります．\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "gQ0oDxG3iygx"
      },
      "source": [
        "----\n",
        "## 汎化と過適合\n",
        "----"
      ]
    },
    {
      "cell_type": "code",
      "source": [
        "# 準備あれこれ\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "import pandas as pd\n",
        "import seaborn\n",
        "seaborn.set()"
      ],
      "metadata": {
        "id": "RjLTO96x2EvK"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "### 多項式回帰\n",
        "\n",
        "「汎化と過適合」の説明で使う道具として，1入力1出力のデータに多項式を当てはめる **多項式回帰** を導入する．"
      ],
      "metadata": {
        "id": "2wnmuzmQ16t6"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "**［多項式回帰（最小二乗法による多項式当てはめ）の問題設定］**\n",
        "\n",
        "変数 $x$ と $y$ の値のペア $N$ 組から成るデータ\n",
        "$$\n",
        "(x_1, y_1), (x_2, y_2), \\ldots, (x_N, y_N)\n",
        "$$\n",
        "がある．変数 $x$ の値から $y$ の値が決まるものとし，そのモデルとして $(D+1)$ 個のパラメータをもつ $D$ 次多項式\n",
        "\n",
        "$$\n",
        "y = f(x) = w_0 + w_1 x + w_2x^2 + \\cdots + w_Dx^{D} \\qquad (1)\n",
        "$$\n",
        "\n",
        "を考える．\n",
        "\n",
        "\n",
        "このとき，モデルの出力 $f(x_n)$ とその正解の値 $y_n$ との間の二乗誤差の和\n",
        "$$\n",
        "\\sum_{n=1}^{N}(y_n - f(x_n))^2\n",
        "$$\n",
        "を最小にするパラメータ $w_0, w_1, \\ldots, w_D$ を求めたい．\n",
        "\n"
      ],
      "metadata": {
        "id": "BSD8USkbK3jD"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "この多項式回帰の問題は，線形回帰（平面当てはめ）の問題をちょこっと変形したものとみなせる．\n",
        "\n",
        "平面当てはめでは，$1, x_1, x_2, \\ldots, x_D$ という $(D+1)$ 個の値をならべた $(D+1)$ 次元ベクトル（$1\\times (D+1)$ 行列）\n",
        "$$\n",
        "\\mathbf{x} = (1, x_1, x_2, x_3, \\ldots, x_D)\n",
        "$$\n",
        "を考えたが，そのかわりに\n",
        "$$\n",
        "\\mathbf{x} = (1, x, x^2, x^3, \\ldots, x^D)\n",
        "$$\n",
        "としてやれば ok．こうすれば，平面当てはめと同様に$(D+1)$個のパラメータをならべたベクトル（$1\\times (D+1)$ 行列）を\n",
        "$$\n",
        "\\mathbf{w} = (w_0, w_1, w_2, \\ldots, w_D)\n",
        "$$\n",
        "として，式$(1)$を\n",
        "$$\n",
        "y = f(\\mathbf{x}) = \\mathbf{w}\\cdot\\mathbf{x} \\qquad (2)\n",
        "$$\n",
        "と表せる．以下，平面当てはめの定式化と同じなので，正規方程式も同じ形になる．\n",
        "\n",
        "つまり，$N\\times(D+1)$ 行列 $X$ と $N\\times 1$ 行列 $Y$ を\n",
        "\n",
        "$$\n",
        "X = \\begin{pmatrix}\n",
        "\\mathbf{x}_1\\\\\n",
        "\\mathbf{x}_2\\\\\n",
        "\\vdots\\\\\n",
        "\\mathbf{x}_N\\\\\n",
        "\\end{pmatrix}\n",
        "=\n",
        "\\begin{pmatrix}\n",
        "1 & x_{1} & x_{1}^{2} & \\cdots & x_{1}^{D}\\\\\n",
        "1 & x_{2} & x_{2}^{2} & \\cdots & x_{2}^{D}\\\\\n",
        "& & \\vdots\\\\\n",
        "1 & x_{N} & x_{N}^{2} & \\cdots & x_{N}^{D}\\\\\n",
        "\\end{pmatrix}\n",
        "\\qquad\n",
        "Y = \\begin{pmatrix}\n",
        "y_1\\\\ y_2 \\\\ \\vdots \\\\ y_N\n",
        "\\end{pmatrix}\n",
        "$$\n",
        "\n",
        "とおくと，二乗誤差を最小にするパラメータ $\\mathbf{w}$ は次の正規方程式の解である．\n",
        "\n",
        "$$\n",
        "X^{\\top}X\\mathbf{w}^{\\top}\n",
        "= X^{\\top}Y\n",
        "$$\n"
      ],
      "metadata": {
        "id": "xrIKO0QWGZtX"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "### 正弦波の多項式回帰の実験"
      ],
      "metadata": {
        "id": "xMP4BEJO4RjV"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### 問題"
      ],
      "metadata": {
        "id": "kcvzFrEp7IJJ"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "$x$ と $y$ の間の真の関係が $y = h(x) = \\sin{(2\\pi x)} + 1$ という式で表されるときに，区間 $[0, 1]$ で一様分布する $x$ の値が $N$ 個得られ，それぞれの値に対応する $y$ の値も同数得られるとする．それらを $(x_n, y_n)$ ($n = 1, 2, \\ldots, N$) とする．\n",
        "ただし，得られる $y_n$ にはノイズが乗っており，$y_n = h(x_n) + \\epsilon_n$ と表されるものとする．$\\epsilon_n$ は標準偏差 $0.1$ の正規乱数である．\n",
        "\n",
        "$(x_n, y_n)$  ($n = 1, 2, \\ldots, N$) のみが与えられる（$h(x)$ の式やどんなノイズが乗っているかは未知）ときに，$x$ から $y$ の値を予測するモデルを作りたい．\n"
      ],
      "metadata": {
        "id": "QwSFt2MR413-"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "### データを生成する関数\n",
        "#\n",
        "def generateData(N, trueFunc=False, seed=None, sigma=0.0):\n",
        "\n",
        "    if trueFunc:\n",
        "        x = np.linspace(-0.1, 1.1, num=N)\n",
        "        y = np.sin(2*np.pi*x) + 1\n",
        "    else:\n",
        "        # 乱数生成器\n",
        "        rng = np.random.default_rng(seed)\n",
        "        x = rng.uniform(0.0, 1.0, N)\n",
        "        y = np.sin(2*np.pi*x) + 1 + sigma*rng.standard_normal(N)\n",
        "\n",
        "    return x, y"
      ],
      "metadata": {
        "id": "MtBcrQJmV8_i"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "x_true, y_true  = generateData(100, trueFunc=True)\n",
        "x_noisy, y_noisy = generateData(25, sigma=0.1)\n",
        "fig, ax = plt.subplots()\n",
        "ax.set_xlim(-0.1, 1.1)\n",
        "ax.set_ylim(-0.2, 2.2)\n",
        "ax.plot(x_true, y_true, color='blue', label='true function')\n",
        "ax.scatter(x_noisy, y_noisy, color='red', label='noisy data')\n",
        "ax.legend()\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "3rd04d84mtHu"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### 多項式回帰をやってみよう"
      ],
      "metadata": {
        "id": "MpBSBfug7PR6"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "以下の「多項式次数」をいろいろ変えて実行してみよう．"
      ],
      "metadata": {
        "id": "OA5F9pNW7YNA"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "##### 1, x, x^2, x^3, ..., x^D をならべたデータ行列（N x (D+1)）をつくる\n",
        "#\n",
        "def makeDataMatrix(x, D):\n",
        "\n",
        "    N = x.shape[0]\n",
        "    X = np.zeros((N, D+1))\n",
        "    X[:, 0] = 1\n",
        "    for d in range(1, D+1):\n",
        "        X[:, d] = x**d\n",
        "\n",
        "    return X"
      ],
      "metadata": {
        "id": "UG1n6Z7V7evF"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "### 線形回帰\n",
        "#\n",
        "def LinearRegression(X, y):\n",
        "    # 正規方程式の左辺の行列を A とする\n",
        "    A = X.T @ X\n",
        "    # 正規方程式の右辺のベクトルを b とする\n",
        "    b = X.T @ y\n",
        "    # np.linalg.solve を使って正規方程式を解き，解を w とする\n",
        "    w = np.linalg.solve(A, b)\n",
        "    return w"
      ],
      "metadata": {
        "id": "fc_sd3cs7iBX"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# 実験条件\n",
        "NL, NT = 10, 1000 # 学習データ数，テストデータ数\n",
        "sig = 0.1 # ノイズの標準偏差\n",
        "D = 1 # 多項式次数\n",
        "\n",
        "# 学習データを作って多項式回帰\n",
        "xL, yL = generateData(NL, sigma=sig)\n",
        "XL = makeDataMatrix(xL, D)\n",
        "w = LinearRegression(XL, yL)\n",
        "yL_est = XL @ w\n",
        "msqeL = np.mean((yL - yL_est)**2)\n",
        "\n",
        "# テストデータを作って予測\n",
        "xT, yT = generateData(NT, sigma=sig)\n",
        "XT = makeDataMatrix(xT, D)\n",
        "yT_est = XT @ w\n",
        "msqeT = np.mean((yT - yT_est)**2)\n",
        "\n",
        "## グラフを描く\n",
        "#\n",
        "x_true, y_true = generateData(1000, trueFunc=True)\n",
        "X_true = makeDataMatrix(x_true, D)\n",
        "y_est = X_true @ w\n",
        "fig, ax = plt.subplots()\n",
        "ax.set_xlim(-0.1, 1.1)\n",
        "ax.set_ylim(-0.2, 2.2)\n",
        "ax.plot(x_true, y_true, color='blue', label='true function')\n",
        "ax.scatter(xL, yL, color='red', label='training data')\n",
        "ax.plot(x_true, y_est, color='green', label='fitted curve')\n",
        "ax.legend()\n",
        "plt.show()\n",
        "\n",
        "# 試行ごとの平均二乗誤差の平均を表示\n",
        "print(f'D = {D}, msqeL = {msqeL:.5f}, msqeT = {msqeT:.5f}')"
      ],
      "metadata": {
        "id": "76Eo4yDL7WM7"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "`msqeL` は，学習データに対する平均二乗誤差(mean squared error)である．モデルを $f(x)$ と表すとき，これは次式で与えられる．\n",
        "\n",
        "$$\n",
        "\\frac{1}{N}\\sum_{n=1}^{N}(y_n - f(x_n))^2\n",
        "$$\n",
        "\n",
        "`msqeT` は，学習データとは別に生成されたテストデータに対する平均二乗誤差である．ここでは，テストデータの生成条件は学習データと同じで，データ数 $N$ のみ異なっている．"
      ],
      "metadata": {
        "id": "cCjgwxk97toa"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "観察しよう:\n",
        "- 多項式の次数を変えると `msqeL` はどのように変化するか\n",
        "- 同じく `msqeT` はどのように変化するか"
      ],
      "metadata": {
        "id": "t9YZEySb8cy7"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "次のセルを実行すると，上記の実験をもう少し大規模に行うことができる．ただし，学習データの数を $10$ から $25$ に増やしており，多項式次数 $D$ の値を一つ決めるごとに学習データ（とテストデータ）を1000通り作って1000通りのモデルを学習させている．"
      ],
      "metadata": {
        "id": "aJl3ugF680uY"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# 実験条件\n",
        "Dlist = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10] # 多項式次数\n",
        "Ntrial = 1000 # 試行回数\n",
        "NL, NT = 25, 1000 # 学習データ数，テストデータ数\n",
        "sig = 0.1 # ノイズの標準偏差\n",
        "\n",
        "msqeL = np.zeros((len(Dlist), Ntrial))\n",
        "msqeT = np.zeros((len(Dlist), Ntrial))\n",
        "\n",
        "for i, D in enumerate(Dlist):\n",
        "    for n in range(Ntrial):\n",
        "        xL, yL = generateData(NL, sigma=0.1, seed=2*n)\n",
        "        XL = makeDataMatrix(xL, D)\n",
        "        w = LinearRegression(XL, yL)\n",
        "        yL_est = XL @ w\n",
        "        msqeL[i, n] = np.mean((yL - yL_est)**2)\n",
        "        xT, yT = generateData(NT, sigma=0.1, seed=2*n+1)\n",
        "        XT = makeDataMatrix(xT, D)\n",
        "        yT_est = XT @ w\n",
        "        msqeT[i, n] = np.mean((yT - yT_est)**2)\n",
        "    eL, eT = np.mean(msqeL[i, :]), np.mean(msqeT[i, :])\n",
        "    print(f'# D = {D}, msqeL = {eL:.4g}, msqeT = {eT:.4g}')\n",
        "\n",
        "fig, ax = plt.subplots(1, 2, figsize=(8, 4))\n",
        "ax[0].boxplot(msqeL.T)\n",
        "ax[0].set_xticklabels(Dlist)\n",
        "ax[0].set_ylim(0, 0.1)\n",
        "ax[0].set_xlabel('D')\n",
        "ax[0].set_title('msqeL')\n",
        "ax[1].boxplot(msqeT.T)\n",
        "ax[1].set_xticklabels(Dlist)\n",
        "ax[1].set_ylim(0, 0.1)\n",
        "ax[1].set_xlabel('D')\n",
        "ax[1].set_title('msqeT')\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "mTU0TGWK8z3P"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "グラフの横軸は多項式次数 $D$ を表す．左図の縦軸は，1000回の試行ごとの学習データに対する平均二乗誤差を表し，右図の縦軸は，テストデータに対する同様の値を表す．箱ひげ図なので，の箱の中の赤線が1000試行の平均）．"
      ],
      "metadata": {
        "id": "DVNNukv99pT0"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "観察しよう:\n",
        "- 多項式の次数を変えると `msqeL` はどのように変化するか\n",
        "- 同じく `msqeT` はどのように変化するか"
      ],
      "metadata": {
        "id": "b-X-1g3w-IIf"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "### 汎化と過適合"
      ],
      "metadata": {
        "id": "Nc5xTB7p-KYB"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "データが少ないときに次数の大きい多項式を当てはめようとすると，学習データにはよく当てはまっても，その背後にある真の関数関係をうまくとらえられるとは限らず，未知のデータ（テストデータ）に対する当てはまりは逆に大きくなってしまうことがある．\n",
        "\n",
        "**汎化** できない，**汎化能力** が低い，**過適合** (over-fitting, 過学習)している\n",
        "\n",
        "一般に，パラメータ数の多い複雑な学習機械ほど学習データによく当てはまるが，過適合によって未知データにうまく汎化できなくなりがちである．\n",
        "逆に言えば，複雑なモデルでもたくさんのデータを用いて学習させることができるならば，過適合は起こりにくくなる．\n",
        "\n"
      ],
      "metadata": {
        "id": "HR5kUZiK-VjG"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "以下は，上記と同じ問題で，学習データ数を増やしてみた例である．"
      ],
      "metadata": {
        "id": "vGbj8rr6xp4p"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# 実験条件\n",
        "NL, NT = 100, 1000 # 学習データ数，テストデータ数\n",
        "sig = 0.1 # ノイズの標準偏差\n",
        "D = 9 # 多項式次数\n",
        "\n",
        "# 学習データを作って多項式回帰\n",
        "xL, yL = generateData(NL, sigma=sig)\n",
        "XL = makeDataMatrix(xL, D)\n",
        "w = LinearRegression(XL, yL)\n",
        "yL_est = XL @ w\n",
        "msqeL = np.mean((yL - yL_est)**2)\n",
        "\n",
        "# テストデータを作って予測\n",
        "xT, yT = generateData(NT, sigma=sig)\n",
        "XT = makeDataMatrix(xT, D)\n",
        "yT_est = XT @ w\n",
        "msqeT = np.mean((yT - yT_est)**2)\n",
        "\n",
        "## グラフを描く\n",
        "#\n",
        "x_true, y_true = generateData(1000, trueFunc=True)\n",
        "X_true = makeDataMatrix(x_true, D)\n",
        "y_est = X_true @ w\n",
        "fig, ax = plt.subplots()\n",
        "ax.set_xlim(-0.1, 1.1)\n",
        "ax.set_ylim(-0.2, 2.2)\n",
        "ax.plot(x_true, y_true, color='blue', label='true function')\n",
        "ax.scatter(xL, yL, color='red', label='training data')\n",
        "ax.plot(x_true, y_est, color='green', label='fitted curve')\n",
        "ax.legend()\n",
        "plt.show()\n",
        "\n",
        "# 試行ごとの平均二乗誤差の平均を表示\n",
        "print(f'D = {D}, msqeL = {msqeL:.5f}, msqeT = {msqeT:.5f}')"
      ],
      "metadata": {
        "id": "vV5PhPEsw-Mg"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "## 正則化とリッジ回帰\n",
        "---"
      ],
      "metadata": {
        "id": "-PLcr3DDtLES"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "真の関数関係に合うようなモデルパラメータ数やモデルの複雑度があらかじめわかるならば，適切に設定して汎化性能の高いモデルが得られるかもしれない．しかし，そんなものは通常は分からない．どうしよう？\n"
      ],
      "metadata": {
        "id": "NZ5mF-FR_1Mz"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "上の多項式回帰の実験で，得られたパラメータの値を表示させてみると，次のようになっている．"
      ],
      "metadata": {
        "id": "V2plrAPMA8FV"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# 実験条件\n",
        "NL, NT = 10, 1000 # 学習データ数，テストデータ数\n",
        "sig = 0.1 # ノイズの標準偏差\n",
        "Dlist = [1, 3, 5, 7, 9] # 多項式次数\n",
        "\n",
        "for D in Dlist:\n",
        "    xL, yL = generateData(NL, sigma=sig)\n",
        "    XL = makeDataMatrix(xL, D)\n",
        "    w = LinearRegression(XL, yL)\n",
        "    print(f'D = {D}, w = {w}, ||w||^2 = {w @ w}')"
      ],
      "metadata": {
        "id": "rTwFez--A7bt"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "各行の末尾には，パラメータの $L_2$ ノルム，すなわちパラメータの値の二乗和\n",
        "\n",
        "$$\n",
        "||\\mathbf{w}||^2 = w_0^2 + w_1^2 + \\ldots + w_D^2\n",
        "$$\n",
        "\n",
        "の値も表示してある．"
      ],
      "metadata": {
        "id": "E7nD4ft4BAjR"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "次数の大きい（パラメータ数の多い）モデルでは，パラメータの値が大きくなりがち．そうならないように制約したらどうだろう？"
      ],
      "metadata": {
        "id": "WLNizdegB1ZU"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "### リッジ回帰\n",
        "\n"
      ],
      "metadata": {
        "id": "yrm8NbyBDr8f"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### リッジ回帰とは"
      ],
      "metadata": {
        "id": "kXsaCGW6GTzP"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "線形回帰（最小二乗法による平面当てはめ）の問題において，\n",
        "\n",
        "$$\n",
        "E(\\mathbf{w}) = \\frac{1}{2}\\sum_{n=1}^N (y_n - \\mathbf{w}\\cdot\\mathbf{x}_n)^2$$\n",
        "\n",
        "のかわりに\n",
        "\n",
        "$$\n",
        "F(\\mathbf{w}) = E(\\mathbf{w}) + \\frac{\\alpha}{2}||\\mathbf{w}||^2\n",
        "$$\n",
        "\n",
        "を最小化する．$\\alpha$ は $0$ 以上の定数とする．右辺第2項は，正則化項（$L_2$正則化項）と呼ばれる．二乗誤差に正則化項を加えた誤差関数を最小化する線形回帰を，\n",
        "**リッジ回帰** (Ridge Regression) または $L_2$正則化付き線形回帰という．\n",
        "\n",
        "パラメータの $L_2$ ノルムが過大にならないように学習させることで，モデルの実質的な自由度が抑えられ，過適合が抑制されることが期待できる．\n"
      ],
      "metadata": {
        "id": "iwPBKi9_D7rX"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "#### 正規方程式の導出"
      ],
      "metadata": {
        "id": "x3NpWQshGWbw"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "$F(\\mathbf{w})$ を最小にするパラメータ $\\mathbf{w}$ は次の正規方程式の解である．\n",
        "\n",
        "$$\n",
        "\\mathbf{w}\\left( X^{\\top}X + \\alpha I \\right)\n",
        "= Y^{\\top}X\n",
        "$$\n",
        "\n",
        "両辺を転置して書くと次のようになる．\n",
        "\n",
        "$$\n",
        "\\left( X^{\\top}X + \\alpha I \\right)\\mathbf{w}^{\\top}\n",
        "= X^{\\top}Y\n",
        "$$\n"
      ],
      "metadata": {
        "id": "YSeD61pYGigW"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "### （よだんだよん）正則化という言葉の意味"
      ],
      "metadata": {
        "id": "CtxkcJzpH99g"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "行列 $A$ を $A = X^{\\top}X$ とおく．$A^{\\top} = \\left(X^{\\top}X\\right)^{\\top} = X^{\\top}\\left(X^{\\top}\\right)^{\\top}  = X^{\\top}X = A$ より，$A$ が対称行列であることは明らか．対称行列の固有値は実数であり，固有ベクトルは互いに直交するようにとることができる．証明は省略するが，行列 $A$ は「準正定値」対称行列（固有値が $0$ 以上）となる．\n",
        "\n",
        "一方，任意の正方行列 $A$ の固有値を $\\lambda$，対応する固有ベクトルを（列ベクトルとして） $\\mathbf{u}$ とおくと，$A\\mathbf{u} = \\lambda\\mathbf{u}$ である．このとき，$B =A + \\alpha I$ に対して $B\\mathbf{u} = A\\mathbf{u} + \\alpha\\mathbf{u} = (\\lambda + \\alpha)\\mathbf{u}$ が成り立つ．したがって，$\\lambda + \\alpha$ が $B$ の固有値である．\n",
        "\n",
        "以上をあわせて考えると，パラメータ数に対してデータ数が少なくて $X^{\\top}X$ のランクが落ちている = 正則でない = いくつかの固有値が $0$ である，という場合でも，$X^{\\top}X + \\alpha I$ の固有値は $\\alpha$ 以上になる．したがって，$\\alpha > 0$ とすればこの行列は正則となる．"
      ],
      "metadata": {
        "id": "p7_QFuF2IHSG"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "---\n",
        "### リッジ回帰の実験\n",
        "\n"
      ],
      "metadata": {
        "id": "SjTLItejIp91"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "正弦波の多項式回帰の例でリッジ回帰をやってみよう．\n",
        "\n",
        "※注: ここでは，モデルのパラメータのうち，定数項を表す $w_0$ は正則化の対象外とするようにしています．線形回帰のモデルでは定数項には制約を加えない方が適切なことが多いので．"
      ],
      "metadata": {
        "id": "F4aBDY3BJEfF"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "### リッジ回帰\n",
        "#\n",
        "def RidgeRegression(X, y, lam, excludeBias=True):\n",
        "    # 正規方程式の左辺の行列を A とする\n",
        "    eyet = np.ones(X.shape[1])\n",
        "    if excludeBias:\n",
        "        eyet[0] = 0.0\n",
        "    A = X.T @ X + lam * np.diag(eyet)\n",
        "    # 正規方程式の右辺のベクトルを XTY とする\n",
        "    b = X.T @ y\n",
        "    #print(X.shape, A.shape, b.shape)\n",
        "    # np.linalg.solve を使って正規方程式を解き，解を w とする\n",
        "    w = np.linalg.solve(A, b)\n",
        "    return w"
      ],
      "metadata": {
        "id": "zmAeXyzkBrc_"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "source": [
        "次のセルの変数 `alpha` の値が，正則化項の係数 $\\alpha$ を表している．適当に変えて結果を観察しよう．"
      ],
      "metadata": {
        "id": "aSGcVfakJ4gJ"
      }
    },
    {
      "cell_type": "code",
      "source": [
        "# 実験条件\n",
        "NL, NT = 25, 1000 # 学習データ数，テストデータ数\n",
        "sig = 0.1 # ノイズの標準偏差\n",
        "D = 9 # 多項式次数\n",
        "alpha = 0.0 # 0.0, 0.0001, 0.001, 0.01, 0.1, 1.0\n",
        "\n",
        "# 学習データを作って多項式回帰\n",
        "xL, yL = generateData(NL, sigma=sig)\n",
        "XL = makeDataMatrix(xL, D)\n",
        "w = RidgeRegression(XL, yL, alpha)\n",
        "yL_est = XL @ w\n",
        "msqeL = np.mean((yL - yL_est)**2)\n",
        "\n",
        "# テストデータを作って予測\n",
        "xT, yT = generateData(NT, sigma=sig)\n",
        "XT = makeDataMatrix(xT, D)\n",
        "yT_est = XT @ w\n",
        "msqeT = np.mean((yT - yT_est)**2)\n",
        "\n",
        "## グラフを描く\n",
        "#\n",
        "x_true, y_true = generateData(1000, trueFunc=True)\n",
        "X_true = makeDataMatrix(x_true, D)\n",
        "y_est = X_true @ w\n",
        "fig, ax = plt.subplots()\n",
        "ax.set_xlim(-0.1, 1.1)\n",
        "ax.set_ylim(-0.2, 2.2)\n",
        "ax.plot(x_true, y_true, color='blue', label='true function')\n",
        "ax.scatter(xL, yL, color='red', label='training data')\n",
        "ax.plot(x_true, y_est, color='green', label='fitted curve')\n",
        "ax.legend()\n",
        "plt.show()\n",
        "\n",
        "# 試行ごとの平均二乗誤差の平均を表示\n",
        "print(f'D = {D}, alpha = {alpha}, msqeL = {msqeL:.5f}, msqeT = {msqeT:.5f}')"
      ],
      "metadata": {
        "id": "Ycm3q3yRI_A7"
      },
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "source": [
        "# 実験条件\n",
        "D = 9\n",
        "lamlist = [0.0, 0.00001, 0.0001, 0.001, 0.01] # 正則化パラメータ\n",
        "Ntrial = 1000 # 試行回数\n",
        "NL, NT = 25, 1000 # 学習データ数，テストデータ数\n",
        "sig = 0.1 # ノイズの標準偏差\n",
        "\n",
        "msqeL = np.zeros((len(lamlist), Ntrial))\n",
        "msqeT = np.zeros((len(lamlist), Ntrial))\n",
        "\n",
        "for i, lam in enumerate(lamlist):\n",
        "    for n in range(Ntrial):\n",
        "        xL, yL = generateData(NL, sigma=0.1, seed=2*n)\n",
        "        XL = makeDataMatrix(xL, D)\n",
        "        w = RidgeRegression(XL, yL, lam)\n",
        "        yL_est = XL @ w\n",
        "        msqeL[i, n] = np.mean((yL - yL_est)**2)\n",
        "        xT, yT = generateData(NT, sigma=0.1, seed=2*n+1)\n",
        "        XT = makeDataMatrix(xT, D)\n",
        "        yT_est = XT @ w\n",
        "        msqeT[i, n] = np.mean((yT - yT_est)**2)\n",
        "    eL, eT = np.mean(msqeL[i, :]), np.mean(msqeT[i, :])\n",
        "    print(f'# D = {D}, lam = {lam}, msqeL = {eL:.4g}, msqeT = {eT:.4g}')\n",
        "\n",
        "fig, ax = plt.subplots(1, 2, figsize=(8, 4))\n",
        "ax[0].boxplot(msqeL.T)\n",
        "ax[0].set_xticklabels(lamlist)\n",
        "ax[0].set_ylim(0, 0.1)\n",
        "ax[0].set_xlabel(r'$\\alpha$')\n",
        "ax[0].set_title('msqeL')\n",
        "ax[1].boxplot(msqeT.T)\n",
        "ax[1].set_xticklabels(lamlist)\n",
        "ax[1].set_ylim(0, 0.1)\n",
        "ax[1].set_xlabel(r'$\\alpha$')\n",
        "ax[1].set_title('msqeT')\n",
        "plt.show()"
      ],
      "metadata": {
        "id": "YLxBYP1bJBNK"
      },
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.9.7"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}