---
lang: ja
dir: ltr
breaks: true
slideOptions:
theme: white
transition: slide
---
<style>
.reveal, .reveal h1, .reveal h2, .reveal h3, .reveal h4, .reveal h5, .reveal h6 {
font-family: "Source Sans Pro", "Helvetica Neue", Helvetica, Arial, "Microsoft JhengHei", Meiryo, sans-serif;
}
h1, h2, h3, h4, h5, h6 {
text-transform: none !important;
}
.color-yellow{
color: yellow;
}
.alert {
padding: 15px;
margin-bottom: 20px;
border: 1px solid transparent;
border-radius: 4px;
text-align: left;
padding: 10px 0;
}
.alert-info {
color: #31708f;
background-color: #d9edf7;
border-color: #bce8f1;
}
.alert-success {
color: #3c763d;
background-color: #dff0d8;
border-color: #d6e9c6;
}
.alert-danger {
color: #a94442;
background-color: #f2dede;
border-color: #ebccd1;
}
.reveal .slides span {
text-align: left;
display: inline-block;
}
p, li {
font-size: 0.88em !important;
}
li>p {
font-size: 1em !important;
}
</style>

# 2次補間法を<br>Juliaで実装してみた

[堀川 由人, ほりたみゅ, @Hyrodium](https://hyrodium.github.io/ja)

---

### おしながき

- 自己紹介
- 数理最適化のざっっくり入門
- 2次補間法の導入
- 2次補間法のJulia実装 (1次元)
- 2次補間法のJulia実装 (2次元)
- まとめ

---

### 自己紹介

- 名前
- 堀川由人
- ほりたみゅ
- Hyrodium
- リンク
- 個人サイト: https://hyrodium.github.io/
- GitHub: https://github.com/hyrodium
- bluesky: https://bsky.app/profile/hyrodium.bsky.social
- JuliaTokaiやJuliaLangJaの運営
- 最近活動できてなくてすみません :pray:
- 近況報告
- CVPR2024に行ったついでにセミと戯れてきました
- [221年に一度のセミ祭り](https://hackmd.io/@hyrodium/S1efQbjUA#/)

---

### 数理最適化のざっっくり入門 (1)

(入門というよりは用語の整理)

関数 $f: D \to \mathbb{R} \ (D\subseteq \mathbb{R}^n)$ が与えられたときに
$f(x)$を最小にする$x$を求めたい場面は色々ある

- $f$: 目的関数
- $D$: 定義域
- $x$: 変数
- $x^*$: 最適値

----

### 数理最適化のざっっくり入門 (2)

- 任意の最適化問題に対して最適値$x^*$を求めたい！
- しかし問題の性質で使えるアルゴリズムは変わる
- $D$ は有限集合?
- $f$ は線形?
- $f$ は滑らか?
- $f$ の微分が計算可能?
- $f$ は凸?
- $D$ は凸?
- $f$ は単峰?
- 最適値の近くに初期値を配置できる?
- etc.

以降では非線形連続最適化を扱います

----

### 数理最適化のざっっくり入門 (3)

非線形連続最適化に使えるアルゴリズムの例

- Newton法
- 2階微分まで参照
- 2次近似して放物線の頂点を次の探索点とする
- 最急降下法
- 1階微分まで参照
- $\operatorname{grad}(f)$の向きに下って次の探索点とする
- 焼き鈍し法
- $\exp(-\beta f(x))$を正規化して確率密度関数を定義
- 逆温度$\beta$を上げつつMCMCで探索点を更新

----

### 数理最適化のざっっくり入門 (4)

解きたい問題

<iframe src="https://www.desmos.com/calculator/p5jhl7f2ss" width="600" height="400" frameborder="0" style="border:0" allowfullscreen></iframe>

- 非線形滑らか関数
- 多峰性がある

----

### 数理最適化のざっっくり入門 (4)

Newton法

<iframe src="https://www.desmos.com/calculator/ucvmgqzjub" width="600" height="400" frameborder="0" style="border:0" allowfullscreen></iframe>

収束は早いが初期値依存性が強い

----

### 数理最適化のざっっくり入門 (4)

最急降下法

<iframe src="https://www.desmos.com/calculator/qpd5pzylgs" width="600" height="400" frameborder="0" style="border:0" allowfullscreen></iframe>

初期値依存性が弱いが局所最適解に陥る

----

### 数理最適化のざっっくり入門 (5)

焼き鈍し法

<iframe src="https://www.desmos.com/calculator/rtgnmxu4fv" width="600" height="400" frameborder="0" style="border:0" allowfullscreen></iframe>

多峰性を考慮できるが、計算量が多い
(Desmos力が足りず、デモが不十分…:pray:)

----

### 数理最適化のざっっくり入門 (4)

- Newton法
- 😵‍💫 2階微分の計算が必要
- 最急降下法
- 😔 1階微分の計算が必要
- 🙄 Newton法に比べると遅い
- 🧐 学習率の決め方に工夫が必要
- 焼き鈍し法
- 🧐 逆温度の下げ方やサンプリングに工夫が必要
- 😵 サンプリングベースなので最適値への収束が遅い

微分の計算が不要で収束が早くハイパラの工夫が不要なアルゴリズムが欲しい！
(∵ $f$がブラックボックスなら微分が取れない)

---

### 2次補間法の導入

<iframe src="https://www.desmos.com/calculator/rjb7dlcymi" width="600" height="400" frameborder="0" style="border:0" allowfullscreen></iframe>

- 初期値を3点取る
- 3点を通る放物線の頂点を次の探索点とする

----

### 2次補間法による最適化 (詳細)

3点$(x_1, f(x_1)), (x_2, f(x_2)), (x_3, f(x_3))$が与えられれば
これらを通る放物線$\tilde{f}$が一意に定まる ([Lagrange補間](https://ja.wikipedia.org/wiki/%E3%83%A9%E3%82%B0%E3%83%A9%E3%83%B3%E3%82%B8%E3%83%A5%E8%A3%9C%E9%96%93))

\begin{align}
\scriptsize
\tilde{f}(x) = f(x_1)\frac{(x-x_2)(x-x_3)}{(x_1-x_2)(x_1-x_3)} + f(x_2)\frac{(x-x_1)(x-x_3)}{(x_2-x_1)(x_2-x_3)} + f(x_3)\frac{(x-x_1)(x-x_2)}{(x_3-x_1)(x_3-x_2)}
\end{align}

この係数を整理して次の探索点が計算できる

----

### ところで記法の提案

- Newton法 $\newcommand\iter[2]{\underset{#1}{\underline{#2}}}$
- $\iter{n}{x} = \iter{n-1}{x} - \iter{n-1}{\frac{f''(x)}{f'(x)}}$
- 最急降下法
- $\iter{n}{x} = \iter{n-1}{x} - \alpha f'(\iter{n-1}{x})$
- 記法のメリット
- 添字の上下が入り乱れたときに便利
- 手書き数式にも使いやすい

この記法を広めたい！！

----

### 2次補間法による最適化 (詳細2)

以下の式で探索を繰り返す！ $\newcommand\iter[2]{\underset{#1}{\underline{#2}}} \newcommand\Pare[1]{\left(#1\right)}$

\begin{align}
\small
\iter{n}{x} = \frac{1}{2} \cdot \frac{\iter{n-1}{f(x)}(\iter{n-2}{x^2}-\iter{n-3}{x^2}) + \iter{n-2}{f(x)}(\iter{n-3}{x^2}-\iter{n-1}{x^2}) + \iter{n-3}{f(x)}(\iter{n-1}{x^2}-\iter{n-2}{x^2})}{\iter{n-1}{f(x)}(\iter{n-2}{x}-\iter{n-3}{x}) + \iter{n-2}{f(x)}(\iter{n-3}{x}-\iter{n-1}{x}) + \iter{n-3}{f(x)}(\iter{n-1}{x}-\iter{n-2}{x})}
\end{align}

---

### 2次補間法のJulia実装 (1次元)

function rec1D(x₁, x₂, x₃)
    a₁ = f(x₁)*(x₂-x₃)
    a₂ = f(x₂)*(x₃-x₁)
    a₃ = f(x₃)*(x₁-x₂)
    return (a₁*(x₂+x₃)+a₂*(x₃+x₁)+a₃*(x₁+x₂))/2(a₁+a₂+a₃)
end

f(x) = sin(x) + x^2/10
xs_init = [1.2, 0.1, -2.2]
xs = copy(xs_init)
for _ in 1:100
    push!(xs, rec1D(xs[end], xs[end-1], xs[end-2]))
end

----

### 1次元問題, 収束の様子

using Plots
plot(f; xlims=(-5,5), color=:red3, label="objective")
plot!(xs, f.(xs); color=:blue3, label="iteration")
scatter!(xs_init, f.(xs_init); color=:blue3, label="initial points")


----

### 1次元問題, 目的関数の値の減少

目的関数に対して単調減少に探索できてるか確認

plot(f.(xs); label="objectve")


----

### 1次元問題, 目的関数の値の減少 (log)

logスケールで確認
→探索が進むにつれて放物線近似の数値誤差が拡大

plot(f.(xs).-minimum(f.(xs)).+eps(Float64); label="objectve", yscale=:log10)


---

### 2次補間法のJulia実装 (2次元)

やり方は同じ

- 多項式$ax^2+bxy+cy^2+dx+ey+f$で補間
- 頂点を次の探索点とする
- 初期値が6つ必要 ($n$次元なら$\frac{(n+1)(n+2)}{2}$つ)

using LinearAlgebra
function rec2D(xs, ys)
    zs = f.(xs,ys)
    A = hcat([[x^2, x*y, y^2, x, y, 1] for (x,y) in zip(xs, ys)]...)'
    a,b,c,d,e,_ = pinv(A)*zs
    p,q = -[2a b;b 2c]\[d;e]
    return p,q
end

----

### 最適化問題の設定

\begin{align}
f\left(x,y\right)=x^{2}+\sin\left(x\right)+1.5y^{2}+\sinh\left(y\right)-\frac{xy}{5}
\end{align}

<iframe src="https://www.desmos.com/calculator/ppkvx1ldkn" width="600" height="400" frameborder="0" style="border:0" allowfullscreen></iframe>

----

### 解いてみる

using Random
Random.seed!(42)
f(x,y) = x^2 + sin(x) + 1.5y^2 + sinh(y) - x*y/5
xs_init = rand(6)
ys_init = rand(6)
xs = copy(xs_init)
ys = copy(ys_init)
for _ in 1:100
    x, y = rec2D(xs[end-5:end], ys[end-5:end])
    push!(xs,x)
    push!(ys,y)
end

----

### 2次元問題, 収束の様子

xs_plot = -3:0.1:3
ys_plot = -5:0.1:3
zs_plot = f.(xs_plot', ys_plot)
plot(xs_plot, ys_plot, zs_plot; levels=-40:40, label="objective")
plot!(xs, ys; color=:blue2, label="iteration")
scatter!(xs_init, ys_init, label="initial points")


----

### 2次元問題, 目的関数の値の減少

plot(f.(xs, ys); label="objective")


----

### 2次元問題, 目的関数の値の減少 (log)

plot(f.(xs, ys) .- minimum(f.(xs, ys)) .+ eps(Float64); label="objective", yscale=:log10)


----

### FAQ

- なぜ2次近似なんですか？
- 1次近似
- 勾配は分かるが学習率を決める必要がある(すこし面倒)
- 3次近似
- 頂点が存在しなかったり
- 非一意的だったり
- 有名な手法なんですよね？
- わからない
- [しっかり学ぶ数理最適化](https://www.kspub.co.jp/book/detail/5212707.html)に記載なし
- ちょっと探して日本語文献が見つからなかった

---

### まとめ

- 微分しなくても早く最適化できるアルゴリズム！
- 最適化の後半では数値誤差が支配的
- 1次元ではLagrange補間が使えて便利
- 2次元以上はもっと効率的に実装できそう
- 今後の展望
- 文献を探す
- 収束性能について実験継続
- Juliaパッケージ化

----

### 参考文献

見つかった英語文献 (ちゃんと読んでない)

- [Quadratic Interpolation Optimization (QIO): A new optimization algorithm based on generalized quadratic interpolation and its applications to real-world engineering problems](https://www.sciencedirect.com/science/article/abs/pii/S0045782523005704)
- [Quadratic interpolation method of 1D minimization](https://www.youtube.com/watch?v=9Zejl2YzaYY)
- [An Algorithm using Quadratic Interpolation for Unconstrained Derivative Free Optimization](https://link.springer.com/chapter/10.1007/978-1-4899-0289-4_3)