Levenberg-Marquardt法あたりの実装練習

 非線形最小二乗問題を解く上で、

と整理できて、これらを実装してみる。

 サンプルの関数としてお馴染みのRosenbrock functionを使う。

 (1, 1)のところに極小値が存在する。

 これに対してPyhtonで各種最適化手法を実装した。Line Searchなどは使わず、Levenberg-Marquardt法でもλの混合具合の調整などはせず、かなりシンプルに関係を把握する目的のコードになる。

実装

import numpy as np
from copy import deepcopy
from collections import defaultdict
import matplotlib.pyplot as plt


def error1(x1: float, x2: float) -> float:
    return x1 - 1


def error2(x1: float, x2: float) -> float:
    return 10 * (x2 - x1**2)


def rosenbrock_function(x1: float, x2: float) -> float:
    return 1 / 2 * (error1(x1, x2) ** 2 + error2(x1, x2) ** 2)


def rosenbrock_gradient(x1: float, x2: float) -> np.ndarray:
    g1 = error1(x1, x2) * 1.0 + error2(x1, x2) * (-20 * x1)
    g2 = error1(x1, x2) * 0.0 + error2(x1, x2) * 10
    return np.array([g1, g2])


def hessian_matrix(x1: float, x2: float) -> np.ndarray:
    return np.array([[1 - 200 * x2 + 600 * x1**2, -200 * x1], [-200 * x1, 100]])


def jacobian_matrix(x1: float, x2: float) -> np.ndarray:
    row1 = np.array([1, 0])
    row2 = np.array([-20 * x1, 10])
    return np.vstack([row1, row2])


if __name__ == "__main__":
    np.set_printoptions(formatter={"float": "{:.3f}".format})

    initial_guess = np.array([-10.985, 50.940])

    # Newton method
    x = deepcopy(initial_guess)
    dict_newton = defaultdict(list)
    dict_newton["x1"].append(x[0])
    dict_newton["x2"].append(x[1])
    dict_newton["f(x)"].append(rosenbrock_function(x[0], x[1]))
    for i in range(100):
        g = rosenbrock_gradient(x[0], x[1])
        h = hessian_matrix(x[0], x[1])
        x -= np.linalg.inv(h) @ g
        dict_newton["x1"].append(x[0])
        dict_newton["x2"].append(x[1])
        dict_newton["f(x)"].append(rosenbrock_function(x[0], x[1]))
        if np.linalg.norm(g) < 1e-6:
            break
    plt.plot(dict_newton["x1"], dict_newton["x2"], marker="o", label="Newton")

    # Gauss-Newton method
    # print("Gauss-Newton method")
    x = deepcopy(initial_guess)
    dict_gauss_newton = defaultdict(list)
    dict_gauss_newton["x1"].append(x[0])
    dict_gauss_newton["x2"].append(x[1])
    dict_gauss_newton["f(x)"].append(rosenbrock_function(x[0], x[1]))
    for i in range(100):
        g = rosenbrock_gradient(x[0], x[1])
        j = jacobian_matrix(x[0], x[1])
        h_hat = j.T @ j
        x -= np.linalg.inv(h_hat) @ g
        dict_gauss_newton["x1"].append(x[0])
        dict_gauss_newton["x2"].append(x[1])
        dict_gauss_newton["f(x)"].append(rosenbrock_function(x[0], x[1]))
        if np.linalg.norm(g) < 1e-6:
            break
    plt.plot(dict_gauss_newton["x1"], dict_gauss_newton["x2"], marker="x", label="Gauss-Newton")

    # Gradient descent method
    x = deepcopy(initial_guess)
    dict_gradient_descent = defaultdict(list)
    dict_gradient_descent["x1"].append(x[0])
    dict_gradient_descent["x2"].append(x[1])
    dict_gradient_descent["f(x)"].append(rosenbrock_function(x[0], x[1]))
    for i in range(100):
        g = rosenbrock_gradient(x[0], x[1])
        x -= 0.000025 * g
        dict_gradient_descent["x1"].append(x[0])
        dict_gradient_descent["x2"].append(x[1])
        dict_gradient_descent["f(x)"].append(rosenbrock_function(x[0], x[1]))
        if np.linalg.norm(g) < 1e-6:
            break
    plt.plot(dict_gradient_descent["x1"], dict_gradient_descent["x2"], marker="s", label="Gradient descent")

    # Levenberg-Marquardt method
    x = deepcopy(initial_guess)
    dict_levenberg_marquardt = defaultdict(list)
    dict_levenberg_marquardt["x1"].append(x[0])
    dict_levenberg_marquardt["x2"].append(x[1])
    dict_levenberg_marquardt["f(x)"].append(rosenbrock_function(x[0], x[1]))
    for i in range(100):
        g = rosenbrock_gradient(x[0], x[1])
        j = jacobian_matrix(x[0], x[1])
        h_hat = j.T @ j
        x -= np.linalg.inv(h_hat + 0.1 * np.eye(2)) @ g
        dict_levenberg_marquardt["x1"].append(x[0])
        dict_levenberg_marquardt["x2"].append(x[1])
        dict_levenberg_marquardt["f(x)"].append(rosenbrock_function(x[0], x[1]))
        if np.linalg.norm(g) < 1e-6:
            break
    plt.plot(dict_levenberg_marquardt["x1"], dict_levenberg_marquardt["x2"], marker="^", label="Levenberg-Marquardt")

    plt.scatter(initial_guess[0], initial_guess[1], color="blue", label="Initial guess")
    plt.scatter(1, 1, color="red", label="Minimum")
    plt.legend()
    plt.savefig("optimization.png", bbox_inches="tight", pad_inches=0.05)

 (-10.985, 50.940)というかなり離れたところから始めて、ある程度は(1, 1)へ収束していくことがわかった。

 だいぶ大きく跳ね回っていたり最急降下法は途中で止まっていたりなど、かなり怪しいところは多いが、Levenberg-Marquardt法がとりあえず挙動としてそうなりそうというところが見えてきたのでこれでいったん打ち止め。

 本来、実装としては明に逆行列を求めるのではなく、np.linalg.solveを使うと良いとのこと。

おまけ

 ここに興味を持ったきっかけは、fast_gicpではNDTやGICPがlsq_registrationという形で抽象化されており、最適化手法としてGaussNewtonとLevenbergMarquardtが選択できるようになっていたのを発見したことにある。

 もとのPCLではNDTの最適化はかなりいろいろと工夫が凝らされており、なにをやっているのかまだあまり理解できていないが、こうやって整理できるならかなり良さそう。

ガウス過程の実装練習

 ベイズ最適化に興味がありガウス過程を学びたいと思ったので以下の本を読んだ。

 一周目なので細かい式変形は追っておらず、まず大枠の導出の流れを理解することで妥協したが、結局実装部分がわからないとピンとこないので、そこだけPythonで実装した。といっても本に疑似コードが載っており、それをほぼそのまま写すだけ。

 まず、1次元の関数 f(x) = x ^ 3 - 2 xが未知のものとしてあり、 x = -5, -4, \dots, 4, 5の11点についてノイズありで値が与えられるので関数を回帰する問題設定でやってみた。

main_1d.py

"""Gaussian Processの実装練習."""

import matplotlib.pyplot as plt
import numpy as np


def f(x: float) -> float:
    return x**3 - 2 * x


def radial_basis_function(x1: np.ndarray, x2: np.ndarray) -> float:
    diff = x1 - x2
    norm = np.linalg.norm(diff)
    theta1 = 100
    theta2 = 1
    return theta1 * np.exp(-norm / theta2)


def predict(
    x_test: np.ndarray, y_train: np.ndarray, x_train: np.ndarray
) -> tuple[np.ndarray, np.ndarray]:
    assert len(x_train) == len(y_train)
    n = len(x_train)
    m = len(x_test)
    K = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            K[i][j] = radial_basis_function(x_train[i], x_train[j])

    K += 25 * np.eye(n)

    K_inv = np.linalg.inv(K)
    yy = K_inv @ y_train
    mu = []
    var = []
    for i in range(m):
        k = np.zeros((n, 1))
        for j in range(n):
            k[j, 0] = radial_basis_function(x_train[j], x_test[i])
        curr_mu = k.T @ yy
        curr_var = radial_basis_function(x_test[i], x_test[i]) - k.T @ K_inv @ k
        mu.append(curr_mu[0])
        var.append(curr_var[0, 0])
    return np.array(mu), np.array(var)


if __name__ == "__main__":
    MIN_X = -5
    MAX_X = 5
    x = np.arange(MIN_X, MAX_X, 1, dtype=np.float32)
    n = len(x)
    y = f(x)

    # ノイズ付与
    y += np.random.normal(0, 20, n)

    WIDTH = (MAX_X - MIN_X) * 0.25
    x_list = np.arange(MIN_X - WIDTH, MAX_X + WIDTH, 0.1, dtype=np.float32)
    mean_list, variance_list = predict(x_list, y, x)
    gt_list = f(x_list)

    plt.plot(x_list, mean_list, label="predict_mean")
    plt.fill_between(
        x_list, mean_list - variance_list, mean_list + variance_list, alpha=0.3
    )
    plt.plot(x_list, gt_list, linestyle=":", label="GT(y=x**3 - 2 * x)")
    plt.scatter(x, y, marker="x", label="data", color="black")
    plt.legend()
    save_path = "./plot_1d.png"
    plt.savefig(save_path, bbox_inches="tight", pad_inches=0.05)
    print(f"Saved to {save_path}")
    plt.close()

 やってみるとカーネル関数(今回だとRadial Basis Function)のパラメータを上手く設定しないと想定していたようなグラフにならず、そこに気づくのが難しかった。上の図だと手動でいくつか試して \theta _ 1 = 100, \theta _ 2 = 1という値になっている。これがどの程度それっぽいのかはよくわからない。もちろん本でもハイパーパラメータの決め方についてはきちんとフォローされており、パラメータの推定方法まできちんと載っている。微分は大変そうに見えたので、今回は適当な山登り法である程度調整した。

 xがサンプリングされる区間が一定だと見た目が面白くないのでそれもちょっと変えて、パラメータは \theta _ 1 = 1079.53, \theta _ 2 = 1.25となった。 \theta _ 1はまだまだ大きいほうが良さそうなのでちょっと怖い。

 GTのオレンジ点線がある程度は 1\sigma区間に入っているのが正当と考えると、これくらい大きくするのが当たり前か。多分これで間違ってない、と思う。まぁ真剣に使うときになったらもっとよく調べるということで。

 自分の考えている問題だと2Dまでは最低限欲しいと思われるので、xyの2Dでも試した。

main_2d.py

"""Gaussian Processの実装練習."""

import matplotlib.pyplot as plt
import numpy as np


def f(x: np.ndarray) -> float:
    mu = np.array([[0, 0]])
    sigma = [
        [0.2, 0.0],
        [0.0, 10.0],
    ]
    det = np.linalg.det(sigma)
    inv = np.linalg.inv(sigma)
    const = 1.0 / (2.0 * np.pi * np.sqrt(det))
    diff = x - mu
    exp_term = np.exp(-0.5 * diff @ inv @ diff.T)
    exp_term = np.diag(exp_term)
    return const * exp_term


def radial_basis_function(x1: np.ndarray, x2: np.ndarray) -> float:
    diff = x1 - x2
    norm = np.linalg.norm(diff)
    theta1 = 100
    theta2 = 10
    return theta1 * np.exp(-norm / theta2)


def predict(
    x_test: np.ndarray, y_train: np.ndarray, x_train: np.ndarray
) -> tuple[np.ndarray, np.ndarray]:
    assert len(x_train) == len(y_train)
    n = len(x_train)
    m = len(x_test)
    K = np.zeros((n, n))
    for i in range(n):
        for j in range(n):
            K[i][j] = radial_basis_function(x_train[i], x_train[j])

    K += 25 * np.eye(n)

    K_inv = np.linalg.inv(K)
    yy = K_inv @ y_train
    mu = []
    var = []
    for i in range(m):
        k = np.zeros((n, 1))
        for j in range(n):
            k[j, 0] = radial_basis_function(x_train[j], x_test[i])
        curr_mu = k.T @ yy
        curr_var = radial_basis_function(x_test[i], x_test[i]) - k.T @ K_inv @ k
        mu.append(curr_mu[0])
        var.append(curr_var[0, 0])
    return np.array(mu), np.array(var)


if __name__ == "__main__":
    MIN_X = -5
    MAX_X = 5
    # [-5, 5]の範囲で2次元のランダムな点
    x = np.random.uniform(MIN_X, MAX_X, (20, 2))
    n = len(x)
    y = f(x)

    print(f"{x.shape=}")
    print(f"{y.shape=}")

    # ノイズ付与
    # y += np.random.normal(0, 20, n)

    fig = plt.figure(figsize=(20, 10))

    # GTをプロット
    x_list = np.arange(MIN_X, MAX_X, 0.05, dtype=np.float32)
    x_list = np.array(np.meshgrid(x_list, x_list)).T.reshape(-1, 2)
    gt_list = f(x_list)
    ax = fig.add_subplot(311, projection="3d")
    ax.scatter(x_list[:, 0], x_list[:, 1], gt_list, c=gt_list)
    ax.set_xlim(MIN_X, MAX_X)
    ax.set_ylim(MIN_X, MAX_X)
    ax.set_title("GT")

    # サンプル点をプロット
    ax = fig.add_subplot(312, projection="3d")
    ax.scatter(x[:, 0], x[:, 1], y, c=y)
    ax.set_xlim(MIN_X, MAX_X)
    ax.set_ylim(MIN_X, MAX_X)
    ax.set_title("Sample points")

    WIDTH = (MAX_X - MIN_X) * 0.25
    # [-5, 5]の範囲で0.1刻みのxの2dリスト
    x_list = np.arange(MIN_X, MAX_X, 0.1, dtype=np.float32)
    x_list = np.array(np.meshgrid(x_list, x_list)).T.reshape(-1, 2)
    mean_list, variance_list = predict(x_list, y, x)
    gt_list = f(x_list)

    print(f"{x_list.shape=}")
    print(f"{mean_list.shape=}")
    print(f"{variance_list.shape=}")
    print(f"{gt_list.shape=}")

    # 3dプロット
    ax = fig.add_subplot(322, projection="3d")
    ax.scatter(x_list[:, 0], x_list[:, 1], mean_list, c=mean_list)
    ax.set_xlim(MIN_X, MAX_X)
    ax.set_ylim(MIN_X, MAX_X)
    ax.set_title("Predict mean")

    ax = fig.add_subplot(324, projection="3d")
    ax.scatter(
        x_list[:, 0],
        x_list[:, 1],
        mean_list + variance_list,
        c=mean_list + variance_list,
    )
    ax.set_xlim(MIN_X, MAX_X)
    ax.set_ylim(MIN_X, MAX_X)
    ax.set_title("Predict mean + var")

    plt.tight_layout()

    save_path = "./plot_2d.png"
    plt.savefig(save_path, bbox_inches="tight", pad_inches=0.05)
    print(f"Saved to {save_path}")
    plt.close()

 GTとして、y方向にだけ分散値が大きい正規分布のようなものを考えて、それに対して20点分評価できるサンプル点がもらえるので、そこから関数を回帰する。

 また、ベイズ最適化をちょっと意識して、予測した平均と分散を足したもの(簡易UCBみたいな気持ち)もプロットした。

 3Dでのプロットにしてしまったのでちょっと見づらいところもあるが、まぁだいたい近似はできていそうで、mean + varだと手薄になっている奥のところがちゃんと優先されそうというのも見えたのでだいたい満足。

 しかし、ベイズ最適化の文脈になると、ここから次に調べるべき点をUCBだったりEIだったりを基準にして決めていくわけだが、結局そういう代理関数の最適化パートが挟まることになり、これの極大点自体は解析的には求まらない(と理解している)ので、勾配法やサンプリングで見つけていかないといけない。そこでまた計算コストがかかってしまうので、やっぱり大元の対象がとても高コストのときにこそ有効なのかとは感じた。

 簡単な適用はできないかもしれないということがわかっただけでも収穫なので、ガウス過程は手段の一つとして取り出せる程度に頭の一部に置いておくこととして、いったん切り上げか。

Gaussian Splattingを試す

概要理解

 流石に中身を全く理解しないままに動かして結果だけ得るのもつまらないと思って、簡単なガウス球3つをレンダリングする部分だけ簡単に実装して原理を確認した。

実装と結果

サンプル実装

from dataclasses import dataclass
import torch
import roma


@dataclass
class GaussianPoint:
    position = torch.zeros((3))
    rotation = torch.tensor([0.0, 0.0, 0.0, 1.0])
    scale = torch.ones((3))
    color = torch.zeros((3))
    opacity = torch.tensor(0.5)

    # when projected to the screen
    i = torch.tensor(0.0)
    j = torch.tensor(0.0)
    z = torch.tensor(0.0)
    sigma_2d = torch.zeros((2, 2))


def render_image(T: torch.Tensor, # extrinsic translation (3,)
                 W: torch.Tensor, # extrinsic orientation (3, 3)
                 K: torch.Tensor, # intrinsic (3, 3)
                 height: int,
                 width: int,
                 gaussian_points: list[GaussianPoint]
                 ) -> torch.Tensor:
    for gaussian in gaussian_points:
        # convert to camera coordinate
        p_in_world = gaussian.position
        p_from_camera = W.T @ (p_in_world - T)
        p_in_screen = K @ p_from_camera
        i = p_in_screen[1] / p_in_screen[2]
        j = p_in_screen[0] / p_in_screen[2]
        if p_from_camera[2] <= 0:
            gaussian.z = float("inf")
            continue

        R = roma.unitquat_to_rotmat(gaussian.rotation)
        S = torch.diag(gaussian.scale)
        sigma_3d = R @ S @ S.T @ R.T

        # ref. EWA Volume Splatting p.6
        # https://www.cs.umd.edu/~zwicker/publications/EWAVolumeSplatting-VIS01.pdf
        # (with scale of K)
        l = torch.norm(p_from_camera)
        u0 = p_from_camera[0]
        u1 = p_from_camera[1]
        u2 = p_from_camera[2]
        J = torch.zeros((3, 3))
        J[0, 0] = 1 / u2
        J[1, 1] = 1 / u2
        J[0, 2] = -u0 / (u2 ** 2)
        J[1, 2] = -u1 / (u2 ** 2)
        J[2, 0] = u0 / l
        J[2, 1] = u1 / l
        J[2, 2] = u2 / l

        sigma_2d = J @ W @ sigma_3d @ W.T @ J.T
        sigma_2d = sigma_2d[0:2, 0:2]
        sigma_2d = K[0:2, 0:2] @ sigma_2d

        # update
        gaussian.i = i
        gaussian.j = j
        gaussian.z = u2
        gaussian.sigma_2d = sigma_2d

    # sort by z
    gaussian_points = sorted(gaussian_points, key=lambda x: x.z)

    rem_alpha = torch.ones((height, width), dtype=torch.float32)
    pred_image = torch.zeros((height, width, 3), dtype=torch.float32)
    I, J = torch.meshgrid(torch.arange(
        height), torch.arange(width), indexing="ij")
    for gaussian in gaussian_points:
        if gaussian.z == float("inf"):
            continue
        i = gaussian.i
        j = gaussian.j
        u2 = gaussian.z
        sigma_2d = gaussian.sigma_2d

        # calc inv and det
        sigma_2d_inv = torch.inverse(sigma_2d)
        sigma_2d_det = torch.det(sigma_2d)

        # calc diff vec
        diff_i = I - i
        diff_j = J - j
        diff_vec = torch.stack([diff_j, diff_i], dim=2)  # [h, w, 2]
        diff_vec = diff_vec.unsqueeze(-2)  # [h, w, 2] -> [h, w, 1, 2]

        # calc 2d gauss
        # [h, w, 1, 2] @ [2, 2] @ [h, w, 2, 1] -> [h, w, 1, 1]
        exp = -0.5 * diff_vec @ sigma_2d_inv @ diff_vec.transpose(-1, -2)
        exp = exp.squeeze()  # [h, w, 1, 1] -> [h, w]
        factor = 1 / (2 * torch.pi * sigma_2d_det)
        gauss = factor * torch.exp(exp)

        # calc color
        curr_weight = rem_alpha * gaussian.opacity * gauss
        add_color = curr_weight.unsqueeze(-1) * gaussian.color * 50000
        pred_image += add_color
        rem_alpha *= 1.0 - curr_weight

    pred_image = torch.clamp(pred_image, 0, 1)
    pred_image *= 255
    pred_image = pred_image.to(torch.uint8)
    return pred_image

 簡単に3つのGaussianを浮かべて実行した。

  • 青:X方向(画像左右方向)だけ分散大
  • 緑:Y方向(画像上下方向)だけ分散大
  • 赤:Z方向(画像奥行き方向)だけ分散大

 赤いのが横に伸びて見えるのは、左前方にあるものを斜めから眺めるような形になっているからだと思うが、もしかしたらなにか間違えているかもしれない。自信はない。

知見

 射影変換が非線形なのでテイラー展開で近似するというのはなるほど納得だった。この実装で、大本のEWA Volume Splattingの論文だと、カメラ座標系から光線座標系への変換において、Z座標として座標までの距離としている(なので射影変換というか純粋に座標変換?)が、Gaussian Splattingのcuda_rasterizerの実装だとよくある射影変換を偏微分した行列になっており、それで良いらしい。自分の実装だとカメラの内部パラメータをいつかければ良いのかがわからなくて2x2行列にした後でかけてみていた(そうでないと2回かかってしまう気がした)。

 色の計算がよくわからず、そのままだとGaussianの値が小さすぎるので色が付かなかった。仕方なく50000などとんでもない係数をかけているが、これは妙だ。まだ間違っているところはあるのだろうが、あまりデバッグに時間かける気力が出なかったのでここで切り上げた。

疑問点

  • ガウス球なので、不透明度を1にしても、すでに中央は確率密度としては1にならない。なので、完全に不透明な赤い球を表現するのに、ガウス球が一個では足りないし、薄くなっていくことも考えるとなんだかとても不揃いなガウス球をたくさん浮かべないといけなさそう。多分実際は完全な球とかは少ないので問題ない? というか確率密度関数を使うので良いのだっけ? ピクセルが1x1の面積を持つので、確率密度が確率の積分値と一致するということ?

実装しなかった部分

 ガウス球のレンダリングというところまでは20年前にもあるような話ではあるので、Gaussian Splatting論文の貢献部分とは言い難いかもしれない。Gaussian Splatting論文だと、その他ガウス球の生成削除アルゴリズム部分と、レンダリングの高速化パートという概ね2つの貢献があると見なせると思う。

 特に高速化パートの実装が大変そうだ。必要な処理として

  1. 画像を小さなタイルに分割し、そのタイルごとにガウス球がその範囲に十分被っているかチェックする
  2. ガウス球がどのタイルに属するかという情報と奥行き情報を1つの変数にパッキングする
  3. パッキングした変数についてGPU基数ソートを行う
  4. ソートされた結果から、タイルごとに分割して前からαブレンディングをする

などがあり、さらに勾配を計算することなどを考えると実装がとんでもないことになりそう。論文だけ見せられて「再現実装をしなさい」と言われたら、まず半年専念することは覚悟して、それでも全然できないかもしれない。ところで公式実装はライセンスが商用利用不可なので、頑張る人は再現実装をするのだろうか……。

公式を動かす

 AWSIMで取得したデータで試した。PCDとして西新宿の点群マップデータがあるので、それを読み込んで初期値として利用するようにした。scene/dataset_readers.pyに追加で読み込むようなコードを書けばとりあえず動いた。

使用したコード片

def readAutowareSceneInfo(path):
    import pandas as pd
    from scipy.spatial.transform import Rotation
    import yaml
    from copy import deepcopy
    camera_info = yaml.load(open(os.path.join(path, "camera_info.yaml"), 'r'), Loader=yaml.FullLoader)
    width = camera_info["width"]
    height = camera_info["height"]
    K = np.array(camera_info["K"]).reshape(3, 3)
    fx = K[0, 0]
    fy = K[1, 1]
    FovX = 2 * np.arctan(width / (2 * fx))
    FovY = 2 * np.arctan(height / (2 * fy))
    df_cam_pose = pd.read_csv(os.path.join(path, "pose.tsv"), sep="\t")
    train_cam_infos = []
    T0 = None
    min_pos = np.min(df_cam_pose[["x", "y", "z"]].values, axis=0)
    max_pos = np.max(df_cam_pose[["x", "y", "z"]].values, axis=0)
    for i in range(0, len(df_cam_pose), 3):
        row = df_cam_pose.iloc[i]
        R = Rotation.from_quat([row.qx, row.qy, row.qz, row.qw]).as_matrix()
        T = np.array([row.x, row.y, row.z])
        if i == 0:
            T0 = deepcopy(T)
        T -= T0
        c2w = np.eye(4)
        c2w[0:3, 0:3] = R
        c2w[0:3, 3] = T
        w2c = np.linalg.inv(c2w)
        R = w2c[0:3, 0:3].T
        T = w2c[0:3, 3]
        image_path = os.path.join(path, "images", f"{i:08d}.png")
        image = Image.open(image_path)
        train_cam_infos.append(
            CameraInfo(uid=i,
                       R=R,
                       T=T,
                       FovY=FovY,
                       FovX=FovX,
                       image=image,
                       image_path=image_path,
                       image_name=f"{i:08d}",
                       width=width,
                       height=height,
                       )
        )
    test_cam_infos = train_cam_infos

    nerf_normalization = getNerfppNorm(train_cam_infos)

    ply_path = os.path.join(path, "points3D.ply")
    if not os.path.exists(ply_path):
        print("Converting point3d.bin to .ply, will happen only the first time you open the scene.")
        import open3d as o3d
        loaded_pcd = o3d.io.read_point_cloud(f"{path}/pointcloud_map.pcd")

        print(f"before downsample: {len(loaded_pcd.points)}")
        margin = 100  # [m]
        min_bound = min_pos - margin
        max_bound = max_pos + margin
        print(f"{min_bound=}, {max_bound=}")
        bbox = o3d.geometry.AxisAlignedBoundingBox(min_bound, max_bound)
        loaded_pcd = loaded_pcd.crop(bbox)
        print(f"after  downsample: {len(loaded_pcd.points)}")

        xyz = np.asarray(loaded_pcd.points)
        xyz -= T0
        num_pts = len(xyz)

        shs = np.random.random((num_pts, 3)) / 255.0
        storePly(ply_path, xyz, SH2RGB(shs) * 255)
    try:
        pcd = fetchPly(ply_path)
    except:
        print("Error reading ply file")
        exit(1)

    scene_info = SceneInfo(point_cloud=pcd,
                           train_cameras=train_cam_infos,
                           test_cameras=test_cam_infos,
                           nerf_normalization=nerf_normalization,
                           ply_path=ply_path)
    return scene_info

結果

 あまり良い結果にならなかった。特に直進している場合に、なにも存在していないはずの道路上に浮いているものが多い。学習時の表示ではPSNR25あたり出ているということだったが、途中の崩壊しているところではそこまであるとは思えない。序盤の静止中などはきれいなので、平均値としてはそこで底上げされてしまっているのかもしれない。

 流石に単眼カメラの、視野に重複が多くはないシーンでは難しかったか。6方向カメラとかがあればまだ違うかもしれないが、これ以上試すのは大変なので、とりあえずここで切り上げる。

Emergent Communication through Metropolis-Hastings Naming Game with Deep Generative Modelsを読んだメモ

 メトロポリスヘイスティングス名付けゲームが気になったので読んだ。

概要

 まずは図1を見て概要を把握する。

 二者間でだいたい同じものを見ている(共同注意)という前提のもと、明示的なフィードバックなしでサインの授受だけをして、メトロポリスヘイスティングス法に従って受け入れ/棄却を判断していくと、二者間で同じものに対しては同じサインが割り当てられていくというものだと理解した。人間の言語が社会的にちゃんと同じサインが同じものを指示することのモデルになっている。

実験と評価指標

 概要はわかったが、結局それはなにがどうできたら成功と見なせるのか。

 今回だとMNISTをラベルなしで教師なしクラスタリングするタスクを考える。10クラスという事前知識はあり、10個のサインが使える。共同で最適化した際に同じサインには同じクラスが対応するようなクラスタリングができていれば成功となる。(たとえば、個別でクラスタリングした場合、クラスタリング自体は上手くいくかもしれないが、サインとクラスの対応関係が揃うことは望めないのは直観的に理解できる)

 Adjusted Rand Index(ARI)という指標があるのでそれを用いる。

グラフィカルモデル

 手法に戻ってグラフィカルモデルを理解する。こういうグラフィカルモデルに慣れていないのでいろいろ混乱した。

  •  w _ dはlatent variableであるということがまず若干よくわからなかった。他者から見ると観測になりそうだが? でも内部的には確かに潜在変数なのかも
    • さらに実際の分布ではなくサンプリングされたものをやり取りすることになるわけなのでとても概念がややこしい
  •  \thetaから o _ dに矢印が伸びている意味がわからなかったが、VAEのデコーダのパラメータということらしい
  • 普通のVAEだと zとしては標準正規分布を考えるが、別にそうしないこともできるので、ここではGMMと w _ dから z _ dが導かれるというものを考えているらしい

アルゴリズム

 話し手側がサインを振り出して、受け手側がそれを採用するか、自分が思ったサインをまだ信じるかが確率的に決まる。この確率がメトロポリスヘイスティングス法の計算になっている。

 ( w _ d ^ {Sp}のサンプリングまでしてそれを引数として渡す方が、なんとなく依存関係がきれいなような? Communication過程と考えるとこのシグニチャが適切なのだろうか)


 各データについてMH-Communicationをして今回での暫定名付けを決めて、それで受け手側のパラメータを更新する。立場を逆転してMH-Commnicationをやって、元の話し手側のパラメータも更新する。なるほど明快。

コード

 さらに具体的に理解するにはコードを見た方が良い。main.py

 各iterationでのループが、

  1. それぞれのVAE学習
    • VAEのパラメータは毎回初期から学習開始されている
    • GMMのmuとvarがあるのでそれを事前分布として使うことになるのだと思う
  2. GMM更新パート
    1. VAEから得られる事後確率から各データのサインを振り出す
    2. それをMH-Communicationで受け入れ/棄却の判定を行う
    3. 得られたラベルでGMMのパラメータを推論する

となっていることがなんとなく理解できてくる。

疑問

 GMMは必須か?→ 「6. Conclusion and discussion」のところに書いてある通り、まぁ他の手法で置き換えることもたぶん可能ではありそう。ニューラルネットにしか馴染みのない人にとってはそこの部分が簡潔になるとわかりやすいなぁと感じる。

 最近の生成系の話題だとトークナイズとかいってVQ-VAE系を使うことが多いが、この手法のVAEをVQ-VAEにすることはできるか?→GMMと繋げるという事自体は難しくなってしまうかもしれないが、理屈として確率モデルとして成立している以上は同じことができそうな気がする。

AHC030復習

 やったこと

  • 解説ページにあるものを読んだ
  • 特にwriter解を参考にして(細かい高速化などは省いて)自分なりにRustからC++に書き直し、自分の本番時より良いスコアを得る提出を実装した

理解内容

 解説ページが充実しているのであえて繰り返すまでもないが、自分の理解をもう一度書いておくと

(1) 複数の占いが所与のとき、事前分布が一様であるという情報から、事後分布は尤度に比例する。尤度は問題で定義されている通りのものを、正規分布の累積分布関数を使って計算できる。ここまではよくあるベイズの定理使用法。

(1') ただ、ありえる配置を全て列挙するのは Mが大きい場合には難しいので、焼き鈍し法で有力な配置をサンプリングする。これはMCMC(の特にメトロポリス法)と近いものと見なせる。

 (AHCラジオでの言及としては、温度1にすれば遷移ルールが同じになり、MCMCだと出現頻度を分布として見たサンプリングとして使うが、今回尤度が求めるので尤度を正規化すれば良い)

(2) 次の占い範囲の決定は、「占い結果がyになる」という事象と「真の配置が x である」という事象の相互情報量を求めて、相互情報量/コストが高いものとしたい。厳密に求めるのは計算量的に難しいので各マスの相互情報量を計算して山登り法で範囲を作る。

 一般に、「複数の選択肢から一つ選んで行動が可能なときに、得られる情報量が一番大きい行動を選択したい」ということはそれなりにありそうで、条件によってどこまで使えるのかはわからないが、考える一つの方針にはなりそう。


 (1), (2)両方とも、今までは「上位陣はなにかよくわからないことをやっているなぁ」という感じだったのが一歩理屈が見えてきて取っ掛かりが掴めてきた、かもしれない。それを実際に使えるレベルとなるとまた遠いのだろうが。

実装面

 手法自体を大まかに理解しても、実装してみると躓くところがいくつかあった。しょうもないパートとしては

  • 積分布関数と確率密度関数を勘違いして意味不明な値が出る
  • 対数尤度から確率に直すときに、最大の対数尤度をベースラインとして引くことを忘れてnanが多発する
  • 何と何の相互情報量を求めるのか混乱する
  • 各クラスがどの情報を管理しているのかという関係がハチャメチャになる

などのことが起きる。今回はwriter解という指針があったので細かい実装のテクニックやクラスの切り方なども困ったらそこに立ち返ればよかったが、自分で考えて実装しているとまとめきるのは相当大変なことになりそう。

焼き鈍しパートがとても難しい

 最終的にはwriter解の遷移をそのまま真似することになったが、それでも精度が出きらなかった。seed 0019などでちゃんと収束していかない。また、処理の高速化もとても難しいのでかなり妥協して、制限時間を超えそうだったら諦めて[tex: N2]を全部掘るという処理にした。それでも間に合うケースでは自分の本番時よりとてもスコアが良いので総合的にも本番時超えには十分だった。

 操作回数が少なくなってきた場合に、できれば滑らかに1マス掘りへ移行させたかったのだが、焼き鈍しの最中で1マス掘りの確定情報と整合しない候補でトラップされて抜けられない、有効な遷移を作れないということになってそれも挫折した。

 何も考えずにやると遷移 O(1)で、評価関数 O(QK) Kはクエリの面積)になってしまうが、遷移 O(Q)、評価関数 O(Q)とした方が良いというのは当たり前なんだけどそう。

 遷移の種類として二つの油田をスワップするのはギリギリ思いつきたいとして、スワップするときにできるだけ被るマスが多くなるようにズラす量を事前計算しておいて適用するなどは厳しそう。

所感

 学びの多い回だった。余裕があればこの知識をもってAHC022とかに立ち戻ると良さそうだが、時間。

AtCoder Heuristic Contest 030

 暫定64位。

0086

考えたこと

ソフトマックス分布を使いたい

  M個ある各油田について、それぞれ左上のマス位置が i, jである確率を保持してみたい。これは i, jについてlogitを持つSoftmax分布として表現したくなった。

 この分布を持って、マップ全域について期待値を取っていけばどのマスにどの程度の量の油田がありそうかを予測することができる。

微分で分布を改善できそう

 この分布をクエリの結果から改善したい。各マスについて、予測した値と1マスクエリ結果の値の自乗誤差を損失関数として誤差を微分すれば良さそう。損失の式を考えて連鎖率を追っていくとそれっぽい値は出てくる。

 (しかし、結局Softmax関数の正しそうな勾配の式は遅いし精度も悪化するしで、謎の嘘勾配を使った。この式じゃ間違っているっけ? うーん)

0マスを確定すると大きい

 1マスクエリで0が確定すると、そこに油田の量が発生するような置き方をするのは不可能ということがわかり、かなり分布を削ることができる。それを狙っていく。

 どのマスで1マスクエリをかけるか。確定することで得られる情報が多いのは、0か1以上か不確かなところとなる。つまり、期待値が0.5くらいのマスを優先的に調べていく。この調べ方でかなりスコアが改善された。

おまけ程度に複数マスクエリの使いみちを探す

 各油田から配置パターンをサンプリングしてきて、その配置パターンでクエリをかける。やはり自乗誤差を考えるが、各マスではなくクエリマス全体の和を取ってから損失計算ということになるので精度はそこまで出ないようだった。一応最初の方にフィールドサイズの N回だけ入れた方がちょっとだけマシになったので入れた。


 ということをやると冒頭に貼ったGifのような挙動になる。

所感

 頭ニューラルネットという感じの解法しか思いつかない。

 確率的な要素が多少ある方が単に焼きなましゲームよりは得意な気がする。ここ最近のAHCで成績良かったのもAHC022で、やや近い匂いがするな?

STL-10データでFSQを試す

 最近の世界モデル系手法ではVQ-VAEが当たり前のように使われているので試してみたいが、生のVQ-VAEだとcommitment lossとかentropy lossとか、様々な工夫を入れなければいけないことが気になり、そういう工夫が要らないと主張されていて内容もシンプルなFSQを試してみた。

 公式実装がjaxで存在しているのでこれを利用した。

 タスクとしては単純な画像再構成で、CIFAR-10あたりだとちょっとトイ・タスクすぎるので、STL-10データを使った。

 STL-10データの特徴は冒頭に書いてある通り

  • 10 クラス: airplane, bird, car, cat, deer, dog, horse, monkey, ship, truck
  • 画像サイズ 96x96 pixels
  • 各クラス train500枚、test800枚の画像がある

となる。unlabeledなデータもあるが、今回はとりあえずtrain 500 x 10クラスの5000枚で学習を試した。

 VQ-VAEを使う上で、Encoder-Decoder部分のCNNも、あまり適当なものではなく実績のあるものを使いたかったため、MaskGITの実装を利用した。(両方Google Researchの研究なので)こちらもjaxでありFSQと繋げるのは簡単。

実装

 FSQをlevel3で要素数10としたので、コードブックサイズとしては 3 ^ {10} = 59049と結構大きくなっており、VQ-VAEの方は単純にコードブックサイズ1024としたので、フェアではないかもしれない。FSQのlevel2とするとnanが発生して上手く学習が進まなかった。

結果

比較実験

 まずFSQと補助損失がなにも入っていないvanilla VQ-VAEと比較した。

 再構成損失が明らかにFSQの方が良く、損失の推移も安定的だった。

FSQの長時間実験

 FSQを上記比較実験から10倍の8時間回してどの程度再構成ができるかを確認した。10倍回すとかなり損失も下がる。

学習データにあるもの
データID 結果
0
1
2
3
4

 学習データ内にあるものはまぁまぁ再現ができているといえばできているかもしれないが、No.2のチーターの画像はちょっと細かくてよくわからない感じになっている。

テストデータ
データID 結果
0
1
2
3
4

 テストデータにあるものだとまだ崩れがちである。学習データも学習量も足りていなさそうな中ではよくやっているか。

今後

  • BBFの次状態予測部分にFSQを導入してみる

余談

 FSQのLevel2だと学習が上手くいかないのは、元のコードがバグっているからのような気がした。自前で実装してみると学習は進んだ。しかし、損失の減りはLevel3(10次元)より悪かった。

 前述の通り、コードブックサイズが 3 ^ {10} = 59049 2 ^ {10} = 1024ではかなり違うので仕方がないことではありそう。

BBFで連続的に2つのゲームを実行する

 強化学習エージェントが環境の変化に対応できるかどうか、という点に興味があるため、atariゲームを連続的に切り替える設定での学習を試せるようにした。

 デフォルトのatari環境では各ゲームによって有効な行動数が異なる。たとえばBreakoutだと4種類(NOOP、FIRE、RIGHT、LEFT)だけが有効な行動である一方、Alienだと18行動が有効となる(8方向 + 停止の9通り動きがあり、FIREするかしないかで2倍の18通りになる)。ここが変わってしまうとネットワーク構造なり出力方策の処理を変えなければならなくて不便なため、必要ない行動も含めて全てを入力できるようにした。これはatariの環境を作るときにfull_action_space=Trueというオプションを入れれば良い。ドキュメント

 環境設定周りのBBFの実装がわりと複雑で、環境を作成する関数をRunnerに渡して、Runner内部でそれをAgentを作る関数に渡してAgentにActionNumを教えるという形になっているため、いろいろこねくり回す必要がある。

 簡単に実験したところ、Breakoutで行動が増えることによる性能の低下は、そこまで大きくないようだった。

 今回はまず実装を試すということで、atari 100kベンチマークで使われる26ゲームのうち、BBFが高いスコアを出せていたBreakoutとBoxingの2ゲームを選び、それをこの順番で10kずつ学習させた。

 学習中のHumanNormalizedScoreは以下の通り。

学習中のHumanNormalizedScore

 10kでゲームが切り替わったことで報酬が全然得られなくなっており、学習して持ち直していくことがわかる。

 それぞれのゲームの学習が終わったタイミングでBreakoutとBoxingの両方を100回テストした結果が以下の通り。

評価タイミング BreakoutのHumanNormalizedScore BoxingのHumanNormalizedScore
Breakout学習終了時(10k終了時点) 1.01 -0.72
Boxing学習終了時(20k終了時点) 0.18 1.12

 ある程度想定通りの結果が得られた。

  1. まず、Breakout学習終了時ではBoxingは全く学習されてないのでとても低いスコアになる。
  2. Boxingまで学習するとBoxingのスコアは良くなるが、Breakoutの学習結果が失われてBreakoutのスコアが低くなる。

 やや想定外のこととしては

  1. Breakoutのスコアの下がり方が思ったよりも大きくない
    • よく考えると20kではまだ1回目のパラメータリセットが走らないタイミングかもしれないので、その点で忘却が起きない
    • リプレイバッファは十分大きなサイズで作っており、ゲーム切り替えタイミングで特にリセットはしていないつもりなので、Boxing学習中もバッファ内にはBreakoutのデータが残っていて適宜リプレイされる

 やっていることが複雑化してくるので、各種リセットの挙動が想定通りに実装できているかどうかが不安になる。

 今後はAugmenting Replay in World Models for Continual Reinforcement Learningで示されている指標を計算できるようにして、忘却性と可塑性の評価をできるようにしてみたい。

読書メモ 西田洋平『人間非機械論』

 西田洋平『人間非機械論 サイバネティクスが開く未来』

 著者多数『未来社会と「意味」の境界: 記号創発システム論/ネオ・サイバネティクスプラグマティズム

を読み、面白かったので記録を残しておく。主に前者の本についてがメインで、部分的に後者を参照する。それぞれ有料の本なので、内容をまとめるということはしない。読んで浮かんだ疑問点・所感などを残す。興味がある人は買って読んでください。

面白かった点

 テクニカルな面として、相互再帰によって無限退行を解消できるという話は面白い。単なる再帰よりも性質的に良い点がどれくらい出てくるのだろうか。これは元論文とか、もっとの他の参考資料を辿らないとわからなさそう。場合によってはいろいろ応用が利く概念なのかもしれない。

 観測者の再帰性、あるいは説明の円環あたりも、完全に納得できているわけではないがそれらしい気もする。モデル化の中に自分を含めるかどうかというのは大きな差異になりそう。主観的な知覚刺激の中からどうやってそれを学習するのか。不思議だ。

 あとはHACSの階層性の見出し方は感触の良さを感じる。

(1) オートポイエーシスがどこまで成立するのかよくわからない

 細胞が、タンパク質等を構成素(正確には構造らしいが)オートポイエティック・システムであるというのは納得感がある(生物学について詳しくないので機序を説明せよと言われたら全く自信はないが)。また、生物個体が細胞を構成素とするオートポイエティック・システムであるというのも納得できる。

 一方で神経システムがオートポイエティック・システムであるのかどうかは(書内でも注意書きがあったように)、すぐにはよくわからない。構成素としてはやはり細胞レベルのものを考えるのか、それとも抽象的な神経ネットワーク的なものを考えるのだろうか。また心的システムがオートポイエティック・システムであるかどうかもさらによくわからない。構成素が何に相当するのか、あまり見当がつかない。

 HACSを前提とすると全て構成素はコミュニケーションになるはずなのだろうか。抽象度がかなり高まり、解釈にブレが生じてしまいそうな気がする。会議の例はわかりやすいが、その他に自分で応用して考えてみようとすると恣意的になんでも当てはめてしまうか、全然当てはめられないか、ということになる。

 工学的に、細胞レベルでオートポイエーシスを再現するのは難しそうなので、一歩抽象的な段階で実現できるかどうかというのが争点になりそうだというのはわかる(未来社会と「意味」の境界 第3章)。が、上の通り、それがオートポイエティック・システムなのかどうかそもそもよくわからないので、目指すべきものなのかも不明。正直、工学レベルとしては、とても複雑なノン・トリビアルマシンであるか、オートポイエティックなシステムであるか、という区別にはあまり意味がないのではないか、という方に同意したい。

(2) 性質は明確な違いがあるのか

 オートポイエーシスについても、閉鎖性・自律性についても、非連続な変化というよりも中間を考えたくなった。オートポイエティックである/ないといったものが、完全に非連続的な2つの極であると、応用を考えにくい気がする。

 たとえば、生物も無限に生きられるわけじゃないのでオートポイエーシスも完全な形式では成立していない。なので、オートポイエーシスが成立する時間的長さがあるわけで、それが十分に長いと生物っぽく見えるが、成立時間を徐々に短くしていくとどこかで非生物へ切り替わる瞬間が生じるだろうか。そのオートポイエーシスが成り立つ時間の長さで「オートポイエーシスらしさ」を連続的に考えたくなる。

 閉鎖性・自律性についても、結局構造は外部とやり取りするというのを考えて、どの程度意味が構造的になった情報を受け取るかというのは段階があるのではないかという気がした。「自律性」という言葉を身近なレベルで捉えて卑近な例を考えると

  • 手取り足取り全て教えなければいけない
  • 言葉でやってほしいことを伝えるだけで良い
  • なにも言わずにやる

などの段階が考えられる。相互作用(入力)される情報がリッチであるか、そういう連続的な段階で自律性を考えられないか、とは思った。

 入力から情報を希薄化していく。有益な情報自体も自分で獲得してくださいとなる。それは結局表現学習の一つでしかないということになってしまうか。

その他

 結局こういう考えを知って、では工学的にどうしますかというと、別にそこまでやることは変わらなくて、ある程度妥当そうなモデルを考えてあとは性能次第、というようにも感じる。その「妥当そうなモデル」が厳密な確率的生成モデルになるのか、ちょっと妥協していろいろ近似の入ったニューラルネットワークいろいろになるのかはよくわからないが。

 全然別の方向性としては、ライフゲーム的なものが超高度に進歩して、外界とも接点を持ち、かつ、意味があると錯覚するような振る舞いをできるようになると、それはオートポイエティックなシステムだしハッピー、ということになる道もあるのだろうか。やっぱりそのレベルのオートポイエーシス性は別にあまり興味ないなぁとは思ってしまう。

Bigger, Better, Faster: Human-level Atari with human-level efficiencyを読んだメモ

 実験の結果が

  • Atari 100k ベンチマークで、Human Normalized ScoreのIQM(26ゲーム中の上位25%と下位25%を除いた中間50%ゲームについての平均スコア)が1.045
  • 学習時間は6 hours on single GPU

とのことであり、魅力的。

 Max Schwarzer氏(Google DeepMind)が筆頭著者として関わっている一連の研究があるようだ。

  1. SPR(Self-Predictive Representations)
  2. SR-SPR(Scaled-by-Resetting SPR)
  3. BBF(Bigger, Better, Faster):この論文
    • SR-SPRに細かい工夫を詰め込んでネットワークサイズを大きくすることに成功し、性能向上

手法

 まずこれらの手法はQ値ベースの方策オフ型であり、Atari 100kベンチマークをやるに際しても、リプレイバッファに経験を溜めて定期的にサンプリングして学習をする。

 環境との相互作用1回(行動選択1回)に対して、「リプレイバッファからランダムサンプリングして損失計算して勾配法でネットワークのパラメータ更新」を何回やるか、というReplay Ratioを考える。これは一般的な手法だとRR = 0.25とか、せいぜい1とかになる。しかし、ネットワークの定期的なパラメータリセットを入れるとRR = 8とか16まで伸ばせるようになる。まずここまでがSR-SPRの主張。

 これに加えてBBFでは

  • Harder reset : SR-SPRでは最終層だけ完全リセット、それ以外はランダムパラメータと20%で線形補間、という形でリセットしていた。BBFではもっとリセットを強めて50%で線形補間する
  • Larger network : Impala-CNNを持ち出してきて、CNNのchを4倍などにする
  • Receding update horizon : 価値関数の学習でnステップTD誤差を計算するときのnを10から3まで小さくしていく。これはネットワークがリセットされてから最初の10kステップで行う。リセットは40kステップごとに行うので、要するに最初の学習の25%分で下げていく
    • これはhttps://www.learningtheory.org/colt2000/papers/KearnsSingh.pdfで主張されているらしい
    • nが大きいほど収束は速くなるが、最適価値関数に対して漸近誤差が大きくなるので、最初だけ大きくしておいてから下げて学習を速めることが狙い
    • nは大きければ大きいほど良いと思っていたのでこれは意外な話
  • Increasing discount factor:割引率を0.97から0.997まで上げていく。これもnと同じく学習の最初25%でやる [1512.02011] How to Discount Deep Reinforcement Learning: Towards New Dynamic Strategies
  • Weight decay:0.1という大きな値でAdamWを使って重み減衰を入れる
  • Removing noisy nets:SPRやSR-SPRではNoisy netを使っていたがこれを削除する

ということを行う。

 各工夫がどれだけ効いているか。

Figure 5

 Replay Ratio2と8で実験されている(横軸の範囲が揃っていないの見にくい……)

自分がグラフを見た印象では

  • とにかくNo Resetは論外。ネットワークのリセットは大事
    • しかもSR-SPRのリセットでも結構性能が落ちているので、強めのリセットが重要らしい
  • nを10に固定する場合でもかなり性能が落ちている。びっくり
  • その他はまぁほどほどか

 なので、「ネットワークのリセット」と「nステップTDのnの影響」について調べる。

ネットワークのリセット

 これについては一つ前の論文SR-SPRにいろいろ書かれている。

 どうもニューラルネットワークは学習を進めてしまうと一般化する能力を失うらしい。滑らかに2つのタスクを切り替えていくマルチタスク学習の設定でも、最初のタスクで長く学習しすぎると次のタスクでの性能が上がらないらしい。 ([1910.08475] On Warm-Starting Neural Network Training; [2106.00042] A study on the plasticity of neural networks)

 (確かに、学習によって全く可塑性を失わないのならばもっとFine-Tuningが大幅にできるようなことを意味するはずなので、そういう意味では当たり前なのかもしれない。しかし100kの定数倍程度の学習ステップ数でもそれが大きく効いてくることは意外)

 オンライン強化学習というのは、学習が進めば振る舞いが変わっていくので、ある意味では滑らかに連続したマルチタスク学習と言える(極端な例としてステージ制のゲームだったりすると明確)。Replay Rateを上げてしまうと、学習が進んでいない段階のもので長時間学習してしまうのと同じことになり、その結果次のタスク(それ以降のデータ)での性能が上がらなくなってしまうと捉える。そして一般化能力を回復させる手法がネットワークのリセットである。

 リセットタイミングでの精度低下が嫌な場合は、リセット直後にある程度(ほぼオフライン強化学習として)学習してしまえば良いらしい。なるほどではある。

 その他参考文献

nステップTDのnの影響

https://www.learningtheory.org/colt2000/papers/KearnsSingh.pdf

 だいぶ話が難しくてあまり理解はできていない。phased TDという解析用の条件設定(ただしそんなに非現実な設定というわけでもない)の下で、誤差の上限 \Delta

 
\Delta _ t \le \frac{1 - \gamma ^ k}{1 - \gamma} \sqrt{\frac{3 \log (k/\delta)}{n}} + \gamma ^ k \Delta _ {t - 1}

となるらしく、この右辺第一項はk(ややこしいが、この論文ではnステップTDのnをkとしている)が大きくなると大きくなる。右辺第二項はkが大きくなると小さくなる。なのでkがいい感じのところでトレードオフのバランスが取られるはずということになるらしい。

 若干自信はないが、直感的な解釈としては、nステップを大きくとれば当然nステップ分は実報酬を使うことになり、たいていその分散は大きいのでそれに影響されるということなのだと思う。分散が大きいだけでバイアスではないので長い目で見ると漸近誤差には影響ないのかと思っていたが、確率的な収束を考えると理論的にも考えられるのか。

 理解は浅いが、まぁそうかなという気分にはなってきたので、nをいじるのは必要なことだと感じられてきた。