CA Tech Lounge 1期生のちわはるです。普段はパターン認識や機械学習を専門としている修士一年生です。
CA Tech Loungeはエンジニアを志す人たちが集うCyberAgent主催の学習者コミュニティで、私は機械学習エンジニア・データサイエンティスト職種として参加しています。
そのため、この記事でも機械学習に関する話題で一本投稿させていただきます。
また、こちらはCA Tech Lounge Advent Calendar 2023の12日目の記事となります。
概要
本記事では掲題の通り、受験数学を機械学習likeに最適化して解いてみようという挑戦（と見せかけたニューラルネット崇拝主義へのちょっとしたアンチテーゼ）です。これは流行りのLLMに答案作ってもらおうなどの試みではなく、むしろ簡易的な最適化を実装することで機械学習が一種の関数機能を模倣する様子を見てみることが目的です。
ただ問題を解いて終わりというのも味気ないので、はじめに機械学習を学び始める学生に向けたポエムを書きつつ、それを踏まえて本題に移りたいと思います。
前説：飾らない機械学習の導入
Stable diffusionやChatGPTなどの生成系AIの認知度が飛躍的に高まったここ2年、機械学習に興味や関心を持ち、私の所属する研究室の門戸を叩く学生も増えてきました。しかし、多くの学生は（当時の私もそうだったように）、機械学習の華々しい側面しか知りません。近年の素晴らしいAI技術のイメージやバイアスが却ってじっくり基礎から学習する際の障害になっている感が否めません。
出来上がったサービスを使うだけであればそれで構わないのですが、しっかり学問知識としての機械学習に入門したい学生には「機械学習とは関数近似のための数学的な最適化手法」という一側面があることを常に心に置いてほしいものです。
実際、大抵の場合は機械学習を用いるまでの流れは以下のような流れになります。
- 解きたいタスクAがある
- タスクAを数理的に考えると、最適化問題Bに言い換えることができる
- 与えられたデータを元にBを機械学習手法で最適化しよう
当然、機械に問題を解かせる際にはプログラミング言語という角ばった表現にまで問題を落とし込む必要がありますから、数学的な問題に言い換えないと機械学習は始まりません。
すでに数学的に明確な目的がある入試問題を用いてこの流れをイメージし、道具の使い方を見てみましょう、というのが本記事のモチベーションなわけです。
機械学習で入試問題を解いてみよう
それでは本題に移って問題を解いてみましょう。今回の題材は2018年度の中央大学の問題（一部改題）です。
$q\ を\ 0<q<1\ を満たす実定数とする。関数\ f(x)\ に対し、F(x)\ を以下で定める$
$$ \displaystyle F(x)= \frac{f(x)-f(qx)}{(1-q)x} $$
$このとき、次の恒等式(*)を満たす実定数\ c\ と二次関数\ f(x)=x^2+ax+b\ $を全て求めよ。
$$x(x-1)F(x)-(1+q)xf(x)=cf(x) \ \cdots ( * )$$
実際には、二次式以外ではこの恒等式を満たす解がないことを証明させる面白い小問がありますが、今回は$q$の値に対して$c$と$f(x)$が求まれば良いとします。興味のある人はその証明も含めて普通に解いてみてください、良問です。
さて、通常なら解の次数を絞り込み、そこから係数比較によって解くのが定石でしょうが、その過程を全て機械でモデリングするのは複雑です。では機械に近似的に解いて妥協するとして、どのように近似解を得るのが良いでしょうか。
今回解きたいタスクは任意の実数$x$に対して$(*)$を成立させる$c$および多項式の係数$a,b$の決定ですが、これを別の最適化問題に言い換えることができれば、機械学習likeなアプローチが可能になるでしょう。
この問題で最適化するべきはa,b,cをパラメータとして恒等式を成立させることですから、左辺と右辺が常に同じになることが求めるべき条件です。
つまり、さまざまなxが与えられたときに、左辺と右辺の差を最小化（できれば0に）すれば、この問題は近似的に解けたことになります。誤差の測り方はL2ノルムでもL1ノルムでも構いませんが、今回は単純にL1（絶対値）を取ることにしましょう。最適化問題の言い換えとしては次のようになります。
$$\underset{a, b, c}{\mathrm{argmin}}\ L(x; a,b,c)= \underset{a, b, c}{\mathrm{argmin}}\ \Biggr |x(x-1)F(x;a,b)-(1+q)xf(x;a,b)-cf(x;a,b) \Biggr|$$
これが決まればあとは最適化をするだけです。この問題を一種の回帰問題として捉え、Pytorchを用いてモデリングし、解いていきましょう。
今回は2種類のアプローチを試してみます。
- ネットワークを定義せず、直接的に$a,b,c$を最適化するアプローチ
- ネットワークを定義し、その最適化によって間接的に$a,b,c$を最適化するアプローチ
さて、どちらがうまくいきそうでしょうか？
直接的に変数を最適化するアプローチ
手順は次のようになります。
- 初期化:
- $a, b, c$の初期値を選び、PyTorchのパラメータとして定義します
- これらの値は最適化過程で更新されます
- データ生成と損失関数の計算:
- $x$ の値をランダムに生成し、各$x=x'$に対しての損失$L(x';a,b,c)$を計算します
- 今回は[-10, 10]の範囲で100点をランダムに取って学習用の入力としました。
- 最適化:
- 損失関数を最小化するように$a,b,c$を更新します
import torch
import numpy as np
from torch.autograd import Variable
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 今回はqをハイパーパラメータとします
q = 0.5
# a,b,cの初期化
a = Variable(torch.randn(1)-1, requires_grad=True)
b = Variable(torch.randn(1)+1, requires_grad=True)
c = Variable(torch.randn(1), requires_grad=True)
# 関数fの定義
def f(x, a, b):
return x**2 + a * x + b
#関数Fの定義
def F(x, a, b, q):
return (f(x, a, b) - f(q * x, a, b)) / ((1 - q) * x)
# オプティマイザ
optimizer = torch.optim.Adam([a, b, c], lr=0.1)
# 損失関数 L
def loss_fn(x_values, a, b, c, q):
loss = 0
for x in x_values:
x = torch.tensor(x)
left_side = x * (x - 1) * F(x, a, b, q) - (1 + q) * x * f(x, a, b)
right_side = c * f(x, a, b)
loss += torch.abs(left_side - right_side)
return loss
# 学習エポック数
epochs = 10000
# エポックごとにa,b,cを保存
a_values, b_values, c_values, loss_values = [], [], [], []
# 学習
for i, epoch in enumerate(range(epochs)):
x_values = np.random.uniform(-10, 10, 100) # 入力x
optimizer.zero_grad()
loss = loss_fn(x_values, a, b, c, q)
loss.backward()
optimizer.step()
a_values.append(a.item())
b_values.append(b.item())
c_values.append(c.item())
loss_values.append(loss.item())
# 100epochごとに損失とa,b,cを表示
if epoch % 100 == 99:
print(f"Epoch {epoch+1}/{epochs}, Loss: {loss.item()}")
print(f" a:{a.item():.4f}, b:{b.item():.4f}, c:{c.item():.4f}")
# 最終的な回帰値
final_a = a.item()
final_b = b.item()
final_c = c.item()
final_loss = loss.item()
print(final_a, final_b, final_c, final_loss)
# epochごとの学習の可視化
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
scatter = ax.scatter(a_values, b_values, c_values, c=loss_values, cmap='plasma', marker='o')
ax.plot(a_values, b_values, c_values, color='black', linewidth=0.5)
ax.set_xlabel('a values')
ax.set_ylabel('b values')
ax.set_zlabel('c values')
ax.set_title('Changes in a, b, and c over Epochs')
plt.show()
実際に何度も学習を回すと、$q=0.5$では以下の3種類のパターンのどれかに近い値が主な収束値となりました。
- $(a,b,c)=(0.0, 0.0, -1.5)$
- $(a,b,c)=(-1.0, 0.0, -1.0)$
- $(a,b,c)=(-3.0, 2.0, 0.0)$
今回の問題の解としては、$\displaystyle (a, b, c)=(0,0,-1-q)(-1,0,-1),(-1-\frac{1}{q},\frac{1}{q},0)$
ですので、$q=0.5$の場合の3つの解にうまく収束したことがわかります。
どの解に収束するかは軒並み初期値に依存しており、勾配降下法などの最適化アルゴリズムの場合には初期値が重要な役割を果たすことが視覚的に理解できます。
$$学習結果が(a,b,c)=(0.0065, -0.0011, -1.4940)となった学習過程の可視化$$
$$学習結果が(a,b,c)=(-1.0233, 0.0040, -0.9948)となった学習過程の可視化$$
$$学習結果が(a,b,c)=(-3.0121, 1.9978, -0.0129)となった学習過程の可視化$$
ネットワークの最適化によって間接的に変数を最適化するアプローチ
手順は次のようになります。
- 初期化:
- 3次元出力の3層ニューラルネットを定義します
- $a, b, c$は明示的に定義せず、ネットワークの出力を$a,b,c$と考えます
- ネットワークが持つパラメータが最適化過程で更新されます
- データ生成と損失関数の計算:
- これは先ほどとほぼ同様です
- 係数の爆発を防ぐため、正則化項をつけています
- 最適化:
- これも同様です。損失関数を最小化するように$a,b,c$を更新します
import torch
import torch.nn as nn
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# 今回はqをハイパーパラメータとします
q = 0.5
# ネットワークの定義
class RegressionNet(nn.Module):
def __init__(self):
super(RegressionNet, self).__init__()
# Define layers
self.fc1 = nn.Linear(1, 5) # Input layer
self.fc2 = nn.Linear(5, 5) # Hidden layer
self.fc3 = nn.Linear(5, 3) # Output layer (a, b, c)
def forward(self, x):
x = torch.relu(self.fc1(x))
x = torch.relu(self.fc2(x))
x = self.fc3(x)
return x
# モデルのインスタンス化
net = RegressionNet()
# オプティマイザ
optimizer = torch.optim.Adam(net.parameters(), lr=0.001)
def loss_function(output, x, q):
a, b, c = output[:, 0], output[:, 1], output[:, 2]
left_side = x * (x - 1) * F(x, a, b, q) - (1 + q) * x * f(x, a, b)
right_side = c * f(x, a, b)
return torch.mean(torch.abs(left_side - right_side)) + torch.sum(a**2 + b**2 + c**2)*
# 学習epoch数
epochs = 200000
# エポックごとにa,b,cを保存
a_values, b_values, c_values, loss_values = [], [], [], []
for epoch in range(epochs):
x_values = torch.tensor(np.random.uniform(-10, 10, 100).reshape(-1, 1), dtype=torch.float32)
optimizer.zero_grad()
outputs = net(x_values)
loss = loss_function(outputs, x_values, q)
loss.backward()
optimizer.step()
a_values.append(torch.mean(outputs[:,0]).item())
b_values.append(torch.mean(outputs[:,1]).item())
c_values.append(torch.mean(outputs[:,2]).item())
loss_values.append(loss.item())
if epoch % 100 == 99:
print(f"Epoch {epoch+1}/{epochs}, Loss: {loss.item()}")
print(f" a:{torch.mean(outputs[:,0]):.4f}, b:{torch.mean(outputs[:,1]):.4f}, c:{torch.mean(outputs[:,2]):.4f}")
# epochごとの学習の可視化
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')
scatter = ax.scatter(a_values, b_values, c_values, c=loss_values, cmap='plasma', marker='o')
ax.plot(a_values, b_values, c_values, color='black', linewidth=0.5)
ax.set_xlabel('a values')
ax.set_ylabel('b values')
ax.set_zlabel('c values')
ax.set_title('Changes in a, b, and c over Epochs')
plt.show()
結果としては、学習がかなり不安定になり、先ほどのようにはうまく正解の解には収束しないことがほとんどでした。
以下の学習過程の画像の通り、局所解にハマって動けなくなっています。
ニューラルネットがうまくいかなかったのはなぜ？
この結果を意外に思った方もいるかも知れません。「ニューラルネットの方がすごそう」だと考える人も多かったでしょうが、今回はシンプルな変数の直接的な最適化の方がうまくいったのです。なぜでしょうか？
今回のニューラルネットベースでのアプローチがうまくいかなかった一番の要因は、ずばり探索するパラメータ空間が広いために局所解に陥りやすかったことでしょう。
直接パラメータを最適化する場合は、学習過程のプロットに表示した通りの3次元空間を動きながら損失を最小化する点を探します。
一方で、ニューラルネットのパラメータは$a,b,c$そのものではなく、ネットワークが持つパラメータであり、今回の非常に単純なモデルですら合計58個の学習パラメータを持っています。つまり、こちらのアプローチでは58次元空間を彷徨いながら損失を最小化する点を探していたのです。
今回扱った入試問題は、数学的には完全に恒等式を満たす解は3種類しかありませんが、局所的におよそ恒等式を満たすような「惜しい解」はたくさんあります。また、$a,b$などが少しでも変われば、$f(x)$の挙動が大きく変化し、損失に強い影響を与える問題でもあります。パラメータに対して損失が綺麗な凸性を持っておらず、でこぼこになっているような問題だったため、非凸性の損失空間に弱いニューラルネットにとっては不向きだったと言えます。
この結果はシンプルな手法が侮れない良い例を与えてくれていると思います。むしろ、ディープラーニングが流行しているからといって、何も疑わずニューラルネットワークを過信してはいけないとも言えるでしょう。目的に合った手法の選択はいつだって大事なのです。
終わりに
大学入試を乗り越えてきた学生ならば、華々しいAIの話題よりも、実はこのような数理パズルを題材にした方が、データサイエンスや機械学習の大事なポイントを効率よく摂取できるのではないか？と思い筆をとりました。少しでも楽しんでいただけたら幸いです。