RIMS 共同研究 (公開型)
数値解析が拓く次世代情報社会〜エッジから富岳まで〜

連続最適化に対する数値解析学的アプローチ

佐藤 峻，牛山寛生，松尾宇泰
東京大学2022 年 10 月 12 日

1 / 27

はじめの例
x(k+1) = x(k) − hk ∇f x(k)連続最適化の視点

数値解析の視点

最小化問題 min f (x) に対する
最急降下法

勾配流 ẋ = −∇f (x) に対する
陽的 Euler 法• 目的：
最適解の発見• hk ：
関数値を基準に選ぶ• キーワード：
収束条件，収束レート，...

• 目的：
解軌道の忠実な再現• hk ：
誤差を基準に選ぶ• キーワード：
収束性，安定性，精度，...

目的は違うが，全く同じ更新式が使われている→ 数値解析の視点は連続最適化に有用である (と期待される)
2 / 27

最適化手法と微分方程式の関係の例
陽的 Euler 法

最急降下法ẋ = −∇f (x)陰的 Euler 法

近接点法[Scieur et al. ’17]
線形多段階法

加速勾配法 (強凸)[Polyak ’64]

Polyak’s heavy ballẍ + γ ẋ + ∇f (x) = 0[野沢 et al. ’22]
可変刻み線形多段階法

加速勾配法[Su et al. ’14]ẍ + 3t ẋ + ∇f (x) = 0

右向きの矢印：連続極限
左向きの矢印：既存の数値解法としての解釈

3 / 27

連続最適化と微分方程式の数値解析

最適化問題

微分方程式

例：無制約最適化問題

例：勾配流

min f (x)ẋ = −∇f (x)x∈Rn

反復による最適解の探索

連続極限最適化手法
例：最急降下法

x(k+1) = x(k) − hk ∇f x(k)

離散化
数値解法

 対応

Nesterov の加速勾配法，Newton 法，...

例：陽的 Euler 法

x(k+1) = x(k) − hk ∇f x(k)Runge–Kutta 法，線形多段階法，...

4 / 27

連続最適化と微分方程式の数値解析
最適化に適したダイナミクス
最適化問題
微分方程式
例：無制約最適化問題

例：勾配流

min f (x)ẋ = −∇f (x)x∈Rn

反復による最適解の探索

連続極限

最適化手法
例：最急降下法

x(k+1) = x(k) − hk ∇f x(k)離散化

最適化に適した離散化手法

数値解法

 対応

Nesterov の加速勾配法，Newton 法，...

例：陽的 Euler 法

x(k+1) = x(k) − hk ∇f x(k)Runge–Kutta 法，線形多段階法，...

4 / 27

我々の周辺で得た成果(A) 最適化に適したダイナミクス
最適化問題

微分方程式

例：無制約最適化問題

例：勾配流

min f (x)ẋ = −∇f (x)x∈Rn

反復による最適解の探索

連続極限最適化手法
例：最急降下法

x(k+1) = x(k) − hk ∇f x(k)離散化(B) 最適化に適した離散化手法

数値解法

 対応

Nesterov の加速勾配法，Newton 法，...

例：陽的 Euler 法

x(k+1) = x(k) − hk ∇f x(k)Runge–Kutta 法，線形多段階法，...

(B) 構造保存数値解法の適用 [Onuma–S. ’22]
Lagrange 未定乗数法 アプローチ [Cheng–Liu–Shen ’20] の勾配流への適用(A) 本質的な収束レート [Ushiyama–S.–Matsuo ’22]
連続時間ダイナミクスにおける収束レートの不定性を除去した本質的収束レート(B) 安定な陽解法の適用 [Ushiyama–S.–Matsuo ’22]：今日の話題(B) 弱い離散勾配を用いた手法の構成と統一的収束証明 [牛山–佐藤–松尾 ’22]
離散勾配の定義を最適化用に拡張，収束レートの証明(B) 加速勾配法の可変刻み線形多段階法解釈 [野沢–松尾–佐藤 ’22]
加速勾配法が可変刻み線形多段階法とみなせることの指摘，手法の拡張

5 / 27

Table of contents

1

最適化における一次法と最適レート

2【本研究】 安定な陽的解法

3

既存研究との関係

4

まとめ

6 / 27

準備扱う問題：無制約連続最適化問題

min f (x)

x∈Rn

仮定：

• f が L-平滑 (L-smooth) ：勾配 ∇f が L-Lipschitz 連続．
• argmin f ̸= ∅：最適解 x⋆ が存在 (f ⋆ := f (x⋆ ))よく使う用語や記法：• ⟨·, ·⟩, ∥·∥：Rn の通常の内積と 2-ノルム．
• (関数 f の) 凸性：
f ((1 − λ)x + λy) ≤ (1 − λ)f (x) + λf (y)(∀ x, y ∈ Rn , λ ∈ [0, 1])• µ-強凸性：f (x) − µ2 ∥x∥2 が凸関数．
(µ-強凸 ⇒ 凸)
7 / 27

最適化手法における「一次法」の枠組
直感的な説明一次法：現時点までの勾配のみを用いて反復を行う手法の総称
一次法の例：
• 最急降下法：x(k+1) = x(k) − hk ∇f x(k)• Nesterov の加速勾配法 (強凸) [Nesterov ’83]：
x(k+1) = y (k) − 1 ∇f y (k)L √√y (k+1) = x(k+1) + √L−√µ x(k+1) − x(k)

L+ µ• 非線形共役勾配法：

d(k) = −∇f x(k) + βk d(k−1)
 (k+1)x= x(k) + hk d(k)2∥∇f (x(k+1) )∥βk+1 =2∥∇f (x(k) )∥8 / 27

最適化手法における「一次法」の枠組
直感的な説明一次法：現時点までの勾配のみを用いて反復を行う手法の総称
一次法の例：• 最急降下法：x(k+1) = x(k) − hk ∇f x(k)• Nesterov の加速勾配法 (強凸) [Nesterov ’83]：x(k+1) = y (k) − 1 ∇f y (k)L√√y (k+1) = x(k+1) + √L−√µ x(k+1) − x(k)L+ µ一次法ではない例：Newton 法−1 
x(k+1) = x(k) − ∇2 f x(k)∇f x(k)8 / 27

最適化手法における「一次法」の枠組
直感的な説明一次法：現時点までの勾配のみを用いて反復を行う手法の総称
一次法の例：
• 最急降下法：x(k+1) = x(k) − hk ∇f x(k)• Nesterov の加速勾配法 (強凸) [Nesterov ’83]：
x(k+1) = y (k) − 1 ∇f y (k)L √√y (k+1) = x(k+1) + √L−√µ x(k+1) − x(k)

L+ µ定義 (一次法：ﬁrst order methods) cf. [Nesterov ’04]

一次法 ⇐⇒ 点列 x(k)が以下の性質を満たす：

n 


ox(k+1) −x(0) ∈ span ∇f x(0) , ∇f x(1) , . . . , ∇f x(k)(∀k ≥ 0)8 / 27µ-強凸関数に対する最適レート
最急降下法は，hk = 2/(µ + L) の下で以下を満たす (κ := L/µ：条件数)：

f x(k)L− f⋆ ≤2

κ−1
κ+1!2kx(0) − x⋆2.

この収束レートは満足できるものか？

9 / 27µ-強凸関数に対する最適レート
最急降下法は，hk = 2/(µ + L) の下で以下を満たす (κ := L/µ：条件数)：

f x(k)L−f ≤2⋆

κ−1
κ+1!2kx(0) − x⋆2.

この収束レートは満足できるものか？
→ No! 下界が知られている (達成する手法も存在)：µ-強凸関数に対する一次法の最適レート
適切な仮定の下，(2 次) cf. [Nesterov ’04]

f x(k)µ−f ≥2⋆µ(一般) [Drori–Taylor ’22] f x(k) − f ⋆ ≥2!2k√2κ−1√x(0) − x⋆ ,κ+1!2k√2κ−1√x(0) − x⋆ .
κ9 / 27

Table of contents

1

最適化における一次法と最適レート

2【本研究】 安定な陽的解法

3

既存研究との関係

4まとめ

10 / 27

最急降下法の収束条件の数値解析的解釈
凸関数に対する最急降下法の収束条件 cf. [Nesterov ’04, Thm. 2.1.14]



最急降下法は，hk < 2/L のとき f x(k) → f ⋆ を満たす．

↓ 数値解析の言葉で言い換え
勾配流に対する陽的 Euler 法は hk < 2/L のもとで平衡点に収束する．
数値解析の手法で同様の条件は出るか？ → 線形安定性解析
観察• 勾配流 ẋ = −∇f (x) の線形化：δ̇ = −∇2 f (x)δ• f が L-平滑な凸関数 ⇒ ∇2 f (x) の固有値は [0, L] に含まれる．
11 / 27

Dahlquist のテスト方程式と安定領域 cf. [Hairer–Wanner ’96]
安定領域

Dahlquist のテスト方程式 (λ ∈ C は定数)ẏ = λy,y(0) = 1

に対して 1 ステップ数値解法を適用した場合の出力を R(hλ) とする．
このとき，R を安定性関数，{z ∈ C | |R(z)| ≤ 1} を安定領域という．
陽的 Euler 法の安定性関数は R(z) = 1 + z1 Im

安定領域：(|λh + 1| ≤ 1)

Re−2−11

前頁の設定では λ ∈ [−L, 0]→ h ∈ [0, 2/L] で安定(最急降下法の収束条件と一致)−112 / 27

より広い安定領域をもつ数値解法へ

安定領域：(|λh + 1| ≤ 1)
−21 Im
今の設定では λ ∈ [−L, 0][−hL, 0] Re → h ∈ [0, 2/L] で安定(最急降下法の収束条件と一致)−1−1→ より広い安定領域をもつ数値解法は，最適化に適していそう．

13 / 27

より広い安定領域をもつ数値解法へ

安定領域：(|λh + 1| ≤ 1)
−21 Im
今の設定では λ ∈ [−L, 0][−hL, 0] Re → h ∈ [0, 2/L] で安定(最急降下法の収束条件と一致)−1−1→ より広い安定領域をもつ数値解法は，最適化に適していそう．

[Hairer–Wanner ’96] より引用．
13 / 27

広い安定領域をもつ陽的解法の構成
通常のアプローチ：
数値解法の構成 (精度などを基準とする) → 安定領域を調べる
広い安定領域をもつ解法を構成するアプローチ：
安定領域を作る → それに合うような数値解法を作る．
具体例 (前ページの図)s 次 Chebyshev 多項式 Ts を「うまく」シフト，スケーリング：

Rs (z) :=Ts (w0 + w1 z)
Ts (w0 )w0 := 1 +η
s2

w1 :=

Ts (w0 )Ts′ (w0 )• 安定領域の左端 ≈ −2Ts′ (1) = −2s2
• 陽的 Euler 法 s 個の合成で実現可能
陽的 Eulers 次の場合

最大ステップ幅

勾配の評価回数比22s21s22s
14 / 27

s 段 RKC (Runge–Kutta–Chebyshev) 法のイメージ(s)Φh (x)：初期値 x に対して s 段 RKC 法を 1 回適用した際の出力．
s = 1, 2, 3, 5 のそれぞれで，h を動かした曲線：
f (x) = 21 (x1 /5)2 + (x2 )2厳密解∇f (x)(2)
Φh (x)xx⋆(5)

Φh (x)(3)

Φh (x)(1)Φh (x)• s を大きくしても (精度は上がっていないので) 厳密解には近づかない．
• s が大きい RKC 法は非常に大きな h を使える (2s2 /L 程度)．

15 / 27

RKC 法の実装について同じ安定性関数をもつ複数の Runge–Kutta 法が存在

• 微分方程式が線形であれば，出力は同じ (安定性関数のみで定まる)• 所望の安定性関数が s 次多項式 ⇒ s 段陽的 Runge–Kutta 法
実装の方法：• 陽的 Euler 法の合成 [Saul’ev ’60],[Guillou–Lago ’61]：• Chebyshev 多項式の根を用いる．
• 根を用いる順番に任意性があり，注意が必要．

• 2 段 RK 法の合成 [Lebedev ’89]：• Chebyshev 多項式の根を 2 つづつペアにして用いる．
• やはり順番には任意性があり，注意が必要．

• Chebyshev 多項式の漸化式を用いる方法 [van der Houwer–Sommeijer ’80]：• Chebyshev 多項式の 3 項間漸化式を用いて定める．
• 内部段も安定に計算でき，計算手順に任意性がない.

16 / 27

数値実験
目的関数 (強凸) (d = 128, L ≈ 512 + α, µ ≈ 0.0759, κ ≈ 6760)：

d 2 X
f (x) =x +(xi − xi−1 )2 + x2d2 1d!

+ α log

i=2dX!

exp(xi ) .

i=1

勾配計算の回数

勾配ノルムが 10−10 以下になるまでの勾配計算の回数：105α=0α = 0.5
α = 1.0
参考 (1/s)

104103 010101s102

強凸ではない (凸でもない) 場合の数値実験でもよく動く：
K. Ushiyama, S., T. Matsuo: Deriving eﬃcient optimization methods based on stable
explicit numerical methods, JSIAM Letters, 14:29–32, 2022.

17 / 27

Table of contents

1

最適化における一次法と最適レート

2【本研究】 安定な陽的解法

3

既存研究との関係

4まとめ

18 / 27

関連する既存研究
目的関数が強凸 2 次関数の場合 (=線形方程式の場合) の研究• RKC 法を勾配流に適用 [Eftekhari–Vandereycken–Vilmart–Zygalakis ’21]:• [van der Houwer–Sommeijer ’80] の実装を利用• スキームパラメータを上手く選ぶと，s → ∞ で (ほぼ) 最適レート• 「強凸 2 次関数+小さい凸関数」の場合への拡張• Chebyshev iterative method：• 強凸 2 次関数専用の最適化手法• 最適なレートを達成• 線形計算における Chebyshev 加速法 (cf. [杉原–室田 ’09])(再掲) µ-強凸関数に対する一次法の最適レート
適切な仮定の下(k)µ−f ≥2⋆(2 次) cf. [Nesterov ’04]

f x(一般) [Drori–Taylor ’22]
µf x(k) − f ⋆ ≥2!2k
√
2
κ−1
√
x(0) − x⋆ ,
κ+1!2k√2κ−1√x(0) − x⋆ .
κ19 / 27

強凸 2 次関数に対する RKC 法 (1/2) [Eftekhari et al. ’21]
RKC 法 (漸化式を利用した実装)
事前に定めるパラメータ：s, η, h.
上のパラメータから事前に計算しておく量：

w0 := 1 +2w0 Ti−1 (w0 )Ts (w0 )η, νi :=, w1 := ′2s
Ts (w0 )Ti (w0 )(i = 1, . . . , s)

反復：


w1
(k)(k)(k)
X1 = X0 − h ∇f X0 ,w0
w1(k)(k)(k)(k)Xi = −hνi ∇f Xi−1 + νi Xi−1 − (νi − 1)Xi−2
w0x(k+1) = Xs(k) .
(k)

X0= x(k) ,(i = 2, 3, ..., s),

20 / 27強凸 2 次関数に対する RKC 法 (2/2) [Eftekhari et al. ’21]
方針 η を固定して，s と h を最適に選ぶ：&rs='(κ − 1)η
,2h=w0 − 1.
w1 µ1
• C(η) := T (w) とすると (1 ステップで勾配を s 回計算していることに注意)，
s0



 

f x(k+1) − f ⋆ ≤ (C(η))2s f x(k) − f ⋆ .
• 収束レートを表す C(η) について，
lim C(η) =η→∞

√
κ−1
√
κ+1

+ O κ−3/2 .

最適レート

s を固定してやりなおせば少し改善できる (次ページ)．
21 / 27

簡単な改善方針 s を固定して，η と h を最適に選ぶ：

η=2s2,
κ−1

h=w0 − 1.
w1 µ1• C(s) := T (w) とすると (1 ステップで勾配を s 回計算していることに注意)，
s0 
f x(k+1) − f ⋆ ≤ (C(s))2s f x(k) − f ⋆ .
• 収束レートを表す C(s) について，
lim C(s) =s→∞

√
κ−1
√κ+1.

最適レート
前述の数値実験はこの実装に基づく．
22 / 27

Chebyshev iterative method (1/2)1
ある程度一般的な 1 次法x(k) = x(k−1) +k−2
X(k−1)(k−1)
cix(i+1) − x(i) + ck−1 ∇f x(k−1)
i=0の中で，2 次関数に対する収束レートが最適になるようなパラメータの決定⇐⇒ 上の手法に関係する多項式 P について，min max |P (λ)|.
P

λ∈[µ,L]→ P を Chebyshev 多項式を用いて定めるのが最適．
ここから遡ってパラメータを定めたものが Chebyshev iterative method．
(RKC 法のときと似た導出になっている．)
1cf. ”F. Pedregosa, Residual Polynomials and the Chebyshev method,
URL: http://fa.bianp.net/blog/2020/polyopt/”

23 / 27

Chebyshev iterative method (2/2)2
Chebyshev iterative method
事前に計算しておく量：

ρ=

κ−1,
κ+1

ν1 = 2,x(1) = −2∇f x(0) .
L+µ

反復 (k = 2, . . . )：

νk =

ρ2
1 − νk−1
4

−1,


x(k) = x(k−1) + (νk − 1) x(k−1) − x(k−2) − νk2∇f x(k−1) .
L+µcf. RKC 法の内部段の計算 (を整理したもの)：
w1(k)(k)(k)(k)(k)Xi = Xi−1 + (νi − 1) Xi−1 − Xi−2 − νi h ∇f Xi−1 .
w02cf. ”F. Pedregosa, Residual Polynomials and the Chebyshev method,
URL: http://fa.bianp.net/blog/2020/polyopt/”

24 / 27

RKC 法と Chebyshev iterative method
RKC 法
陽的 Euler 法の合成

2 段法の合成
漸化式を用いた実装

η を固定した最適スキーム
s を固定した最適スキーム

s → ∞ かつ h → ∞ として「1 回」で終わらせる
Chebyshev iterative methods 回ごとにリスタート• 両者「最善」を目指した結果なので，同じ漸化式が生まれるのはある意味で必然．

• 2 次関数であり，µ と L が精密に分かっている場合には，リスタートは無益．

• RKC 法をうまく使うことで，2 次関数以外の関数に拡張できることを期待．
25 / 27

Table of contents

1

最適化における一次法と最適レート

2【本研究】 安定な陽的解法

3

既存研究との関係

4まとめ

26 / 27

まとめ(B) 構造保存数値解法の適用 [Onuma–S. ’22]
(A) 本質的な収束レート [Ushiyama–S.–Matsuo ’22](B) 安定な陽解法 (RKC 法) の適用 [Ushiyama–S.–Matsuo ’22]：今日の話題• 線形安定性解析の観点から，勾配流に適した数値解法 (RKC 法) を選択• 既存手法との関係の整理(B) 弱い離散勾配を用いた手法の構成と統一的収束証明 [牛山–佐藤–松尾 ’22](B) 加速勾配法の可変刻み線形多段階法解釈 [野沢–松尾–佐藤 ’22]
今後の課題：(A) 最適化に適したダイナミクス
最適化問題

微分方程式

例：無制約最適化問題

反復による最適解の探索(B) 一般の非線形の場合のRKC 法の収束レートの証明

例：勾配流

min f (x)ẋ = −∇f (x)x∈Rn連続極限

最適化手法
例：最急降下法

x(k+1) = x(k) − hk ∇f x(k)離散化(B) 最適化に適した離散化手法

数値解法

 対応

Nesterov の加速勾配法，Newton 法，...

例：陽的 Euler 法

x(k+1) = x(k) − hk ∇f x(k)Runge–Kutta 法，線形多段階法，...

(A) より良いダイナミクス(ODE) の提案(B) より複雑な最適化 ODE に対する離散化の検討
27 / 27