ニュートン法とAdjoint法
1. ニュートン法の基本原理
非線形方程式 $F(x) = 0$ を解く反復解法として、ニュートン法は以下のように定式化される：
x_{k+1} = x_k - [F'(x_k)]^{-1}F(x_k)
ここで：
- $x_k \in \mathbb{R}^n$ は現在の解の推定値
- $F'(x_k) \in \mathbb{R}^{n \times n}$ はヤコビ行列
- $k$ は反復回数
最適化問題への適用では、目的関数 $f(x)$ の最小化に対して：
x_{k+1} = x_k - [\nabla^2 f(x_k)]^{-1}\nabla f(x_k)
となり、ヘッセ行列 $\nabla^2 f(x_k) \in \mathbb{R}^{n \times n}$ を使用する。
1.1 制約付き最適化問題でのニュートン法
制約付き最適化問題では、ラグランジュ関数の2次導関数（KKT系）に対してニュートン法を適用する：
\begin{bmatrix}
\nabla^2_{xx} L & \nabla g^T \\
\nabla g & 0
\end{bmatrix}
\begin{bmatrix}
\Delta x \\
\Delta \lambda
\end{bmatrix} =
-\begin{bmatrix}
\nabla_x L \\
g(x)
\end{bmatrix}
ここで、$\nabla^2_{xx} L \in \mathbb{R}^{n \times n}$, $\nabla g \in \mathbb{R}^{m \times n}$, $\Delta x \in \mathbb{R}^n$, $\Delta \lambda \in \mathbb{R}^m$である。
2. 制約付き最適化問題の定式化
以下の最小化問題を考える：
\min_{x \in \mathbb{R}^n} f(x)
subject to:
g(x) = 0, \quad g: \mathbb{R}^n \to \mathbb{R}^m
ここで、$f: \mathbb{R}^n \to \mathbb{R}$ は十分に滑らかな目的関数、$g: \mathbb{R}^n \to \mathbb{R}^m$ は制約関数である。以下では、$g$ は微分可能で、制約点においてヤコビアン $\frac{\partial g}{\partial x} \in \mathbb{R}^{m \times n}$ は行フルランク（$\text{rank}\left(\frac{\partial g}{\partial x}\right) = m$）と仮定する。この条件は陰関数定理の適用のために必要である。
3. ラグランジュの未定乗数法
ラグランジュ関数を以下で定義する:
L(x, \lambda) = f(x) + \lambda^T g(x)
最適性の必要条件は:
\begin{align}
\frac{\partial L}{\partial x} &= \frac{\partial f}{\partial x} + \lambda^T \frac{\partial g}{\partial x} = 0 \\
\frac{\partial L}{\partial \lambda} &= g(x) = 0
\end{align}
4. パラメータ化による問題の簡略化
状態 $x$ をパラメータ $p \in \mathbb{R}^k$ の関数として表現する:
g(x(p)) = 0
5. 感度解析
陰関数定理を拘束条件 $g(x(p), p) = 0$ に適用すると、$\frac{\partial g}{\partial x}$ が行フルランクという条件下で、$p$ の近傍での状態変数 $x$ の変化を解析できる：
\frac{\partial g}{\partial x}\frac{\partial x}{\partial p} + \frac{\partial g}{\partial p} = 0
ここで、$\frac{\partial g}{\partial x} \in \mathbb{R}^{m \times n}$、$\frac{\partial x}{\partial p} \in \mathbb{R}^{n \times k}$、$\frac{\partial g}{\partial p} \in \mathbb{R}^{m \times k}$ である。
$\frac{\partial g}{\partial x}$ が行フルランク（$m \leq n$）と仮定すると、感度行列は：
S = \frac{\partial x}{\partial p} = -\left(\frac{\partial g}{\partial x}\right)^{+} \frac{\partial g}{\partial p}
ここで、$\left(\frac{\partial g}{\partial x}\right)^{+}$ は $\frac{\partial g}{\partial x}$ の擬似逆行列である。$m = n$ の場合は通常の逆行列となる。以下では簡単のため $m = n$ のケースを考える。
6. 目的関数の導関数
目的関数の導関数を連鎖律を用いて導出する。
初めに、全微分と部分微分の違いについて述べる。一般に、状態 $x \in \mathbb{R}^n$ がパラメータ $p \in \mathbb{R}^k$ の関数として陰的に定義されている場合（$g(x(p),p)=0$）、目的関数 $f(x(p),p)$ の全導関数 $\frac{df}{dp} \in \mathbb{R}^{1 \times k}$ は以下のように表現される：
\frac{df}{dp} = \frac{\partial f}{\partial p} + \left(\frac{\partial x}{\partial p}\right)^T \frac{\partial f}{\partial x}
ここで第一項 $\frac{\partial f}{\partial p} \in \mathbb{R}^{1 \times k}$ は $p$ が $f$ に直接与える影響を表し、第二項は $p$ が $x$ を通して間接的に $f$ に与える影響を表す。以下では一般的な場合として $f$ が $p$ に直接依存する状況を含めて考察する。
- まず、連鎖律より全導関数は：
\frac{df}{dp} = \frac{\partial f}{\partial p} + \left(\frac{\partial x}{\partial p}\right)^T \frac{\partial f}{\partial x}
- 感度行列 $\frac{\partial x}{\partial p}$ は陰関数定理により：
\frac{\partial x}{\partial p} = -\left(\frac{\partial g}{\partial x}\right)^{-1} \frac{\partial g}{\partial p}
- これを代入すると：
\frac{df}{dp} = \frac{\partial f}{\partial p} - \left(\frac{\partial g}{\partial p}\right)^T \left(\left(\frac{\partial g}{\partial x}\right)^{-1}\right)^T \frac{\partial f}{\partial x}
- 転置の性質 $(AB)^T = B^T A^T$ と逆行列の転置 $(A^{-1})^T = (A^T)^{-1}$ より：
\frac{df}{dp} = \frac{\partial f}{\partial p} - \left(\frac{\partial g}{\partial p}\right)^T \left(\frac{\partial g}{\partial x}\right)^{-T} \frac{\partial f}{\partial x}
ここで、$\left(\frac{\partial g}{\partial x}\right)^{-T}$ は $\left(\frac{\partial g}{\partial x}\right)^T$ の逆行列を表す。
7. Adjoint法の導入
目的関数の導関数の式を観察すると：
\frac{df}{dp} = \frac{\partial f}{\partial p} - \left(\frac{\partial g}{\partial p}\right)^T \left[\left(\frac{\partial g}{\partial x}\right)^{-T} \frac{\partial f}{\partial x}\right]
ここで、括弧内の項に着目し、随伴変数（adjoint variable）$\lambda \in \mathbb{R}^m$ を以下で定義する:
\lambda = \left(\frac{\partial g}{\partial x}\right)^{-T} \frac{\partial f}{\partial x}
これは以下の線形方程式（随伴方程式、adjoint equation）を解くことと等価である：
\left(\frac{\partial g}{\partial x}\right)^T \lambda = \frac{\partial f}{\partial x}
この線形方程式を解くことで $\lambda$ が得られる。
8. 最終的な導関数の表現
- 元の式に $\lambda$ を代入すると：
\frac{df}{dp} = \frac{\partial f}{\partial p} - \left(\frac{\partial g}{\partial p}\right)^T \lambda
- 転置の性質より最終形は：
\frac{df}{dp} = \frac{\partial f}{\partial p} - \lambda^T \frac{\partial g}{\partial p}
この形は、ラグランジュの未定乗数法の式と対応している。実際、ラグランジュ関数：
L(x, \lambda, p) = f(x, p) + \lambda^T g(x, p)
の最適性条件：
\frac{\partial L}{\partial x} = \frac{\partial f}{\partial x} + \lambda^T \frac{\partial g}{\partial x} = 0
を変形すると：
\frac{\partial f}{\partial x} = -\lambda^T \frac{\partial g}{\partial x}
となり、これは先ほどの随伴方程式と同じ構造を持つ。このことから、Adjoint法の $\lambda$ はラグランジュ乗数と本質的に同じ役割を果たしていることが理解できる。
この形により、導関数の計算が効率的になる利点がある：
- $\lambda$ を求める線形方程式を1回解くだけでよい（目的関数が複数のパラメータに依存する場合に特に効率的）
- 大規模な逆行列計算を避けることができる
- $\frac{\partial f}{\partial p}$ の項は直接計算可能である
9. 計算効率の比較
パラメータ $p \in \mathbb{R}^k$ に対する導関数 $\frac{df}{dp} \in \mathbb{R}^{1 \times k}$ の計算において、状態変数が $x \in \mathbb{R}^n$ であるとき：
-
感度行列による方法（Forward方向）:
- $\frac{\partial x}{\partial p} \in \mathbb{R}^{n \times k}$ の計算: $O(mn^2k)$ の計算量（$m = n$ を仮定）
- 全体の計算量: $O(n^2k + nk)$
-
Adjoint法（Backward方向）:
- 随伴変数 $\lambda \in \mathbb{R}^m$ の計算: $O(n^2m)$ の計算量（$m = n$ を仮定）
- 全体の計算量: $O(n^2 + nmk)$
パラメータの数 $k$ が状態変数の数 $n$ よりも大きい場合（$k \gg n$）、Adjoint法は大幅な計算効率の向上をもたらす。これは最適化問題や機械学習など、高次元パラメータ空間での勾配計算において特に重要な性質である。
10. トポロジー最適化への応用
コンプライアンス最小化問題の例を考察する：
f(u, \rho) = \int_\Omega F^T u \, d\Omega + \alpha R(\rho)
ここで、$u \in \mathbb{R}^n$ は変位場、$\rho \in \mathbb{R}^m$ は密度場、$R(\rho)$ は密度 $\rho$ に直接依存する正則化項を表す。
線形弾性体の平衡方程式は：
g(u,\rho) = K(\rho)u - F = 0
ここで、$K(\rho) \in \mathbb{R}^{n \times n}$ は対称正定値の剛性マトリクス、$F \in \mathbb{R}^n$ は外力ベクトルである。
Adjoint法を適用すると：
- 目的関数と拘束条件の対応：
\begin{align}
f(u, \rho) &= \int_\Omega F^T u \, d\Omega + \alpha R(\rho) & \text{(目的関数)} \\
g(u,\rho) &= K(\rho)u - F = 0 & \text{(拘束条件)}
\end{align}
- 随伴方程式の導出：
\begin{align}
\left(\frac{\partial g}{\partial u}\right)^T \lambda &= \frac{\partial f}{\partial u} \\
K(\rho)^T \lambda &= F
\end{align}
ここで、$\frac{\partial g}{\partial u} = K(\rho)$ であり、$\frac{\partial f}{\partial u} = F$ となる。$K(\rho)$ は対称行列なので $K(\rho)^T = K(\rho)$ であり、実際には随伴方程式は元の状態方程式と同じ線形システムとなる。したがって、$\lambda = u$ となる特殊なケースである点に注意が必要である。
- 最終的な感度：
\begin{align}
\frac{df}{d\rho} &= \frac{\partial f}{\partial \rho} - \lambda^T \frac{\partial g}{\partial \rho} \\
&= \alpha \frac{\partial R(\rho)}{\partial \rho} - \lambda^T \frac{\partial}{\partial \rho}(K(\rho)u - F) \\
&= \alpha \frac{\partial R(\rho)}{\partial \rho} - \lambda^T \frac{\partial K(\rho)}{\partial \rho}u
\end{align}
この結果は、$f(u(\rho),\rho)$ のような目的関数の全導関数の一般形式と整合している：
\frac{df}{d\rho} = \frac{\partial f}{\partial \rho} + \left(\frac{\partial u}{\partial \rho}\right)^T \frac{\partial f}{\partial u}
ここで、$\frac{\partial f}{\partial \rho} = \alpha \frac{\partial R(\rho)}{\partial \rho}$ は $\rho$ が $f$ に直接与える影響を表し、感度行列 $\frac{\partial u}{\partial \rho}$ は陰関数定理から導出される。この方法により、感度行列を明示的に計算することなく $\frac{df}{d\rho}$ を効率的に求めることが可能である。
この形式は一般的な Adjoint法の構造と完全に一致しており：
- $u$ が状態変数 $x$ に対応
- $\rho$ がパラメータ $p$ に対応
- $K(\rho)$ が拘束条件のヤコビアン $\frac{\partial g}{\partial x}$ に対応
- $F$ が目的関数の偏導関数 $\frac{\partial f}{\partial u}$ に対応
11. 二次の感度解析とニュートン法的アプローチ
目的関数 $T(x(p),p)$ のパラメータ $p \in \mathbb{R}^k$ に関する二次の感度解析を考察する。ここでの $\frac{dT}{dp} \in \mathbb{R}^{1 \times k}$ と $\frac{d^2 T}{dp^2} \in \mathbb{R}^{k \times k}$ は全導関数を表し、パラメータ $p$ の変化が状態変数 $x$ を通じて間接的に $T$ に与える影響と、$p$ が $T$ に直接与える影響の両方を含む。テイラー展開により：
\frac{dT}{dp}(p + \Delta p) \approx \frac{dT}{dp}(p) + \frac{d^2 T}{dp^2}\Delta p
ニュートン法的アプローチでは、$\frac{dT}{dp}(p + \Delta p) = 0$ を解くことで更新量 $\Delta p \in \mathbb{R}^k$ を得る：
\frac{d^2 T}{dp^2}\Delta p = -\frac{dT}{dp}^T
ここで、ヤコビアン $S = \frac{\partial x}{\partial p} \in \mathbb{R}^{n \times k}$ を用いて一階全導関数を表現すると：
\frac{dT}{dp} = \frac{\partial T}{\partial p} + \frac{\partial T}{\partial x}S
ここで、$\frac{\partial T}{\partial x} \in \mathbb{R}^{1 \times n}$, $\frac{\partial T}{\partial p} \in \mathbb{R}^{1 \times k}$ である。
二階全導関数（ヘッセ行列）は連鎖律により次のように導出される：
\frac{d^2 T}{dp^2} = \frac{\partial^2 T}{\partial p^2} + \frac{\partial^2 T}{\partial p \partial x}S + S^T\frac{\partial^2 T}{\partial x \partial p} + S^T \frac{\partial^2 T}{\partial x^2}S + \sum_{i=1}^n \frac{\partial T}{\partial x_i}\frac{\partial^2 x_i}{\partial p^2}
ここで、最後の項は各状態変数 $x_i$ に関する二階感度を含んでいる。実用上は、この項の計算コストが高いため、しばしば無視されるか近似される。
この式は、全導関数 $\frac{d^2 T}{dp^2}$ が直接二階項 $\frac{\partial^2 T}{\partial p^2}$、交差項 $\frac{\partial^2 T}{\partial p \partial x}S + S^T\frac{\partial^2 T}{\partial x \partial p}$、そして間接二階項 $S^T \frac{\partial^2 T}{\partial x^2}S + \sum_{i=1}^n \frac{\partial T}{\partial x_i}\frac{\partial^2 x_i}{\partial p^2}$ の和として表現されることを示している。ニュートン法では、この二階全導関数を用いて効率的に最適解へと収束させることが可能である。
12. ニュートン-Adjoint法の最終形式
前節までの議論を統合すると、ニュートン法とAdjoint法を組み合わせた効率的な最適化アルゴリズムが以下のように定式化される：
- 状態方程式を解く：
g(x, p) = 0
- Adjoint方程式を解いて $\lambda$ を得る：
\left(\frac{\partial g}{\partial x}\right)^T \lambda = \frac{\partial T}{\partial x}^T
- 感度行列（必要な場合）を計算する：
S = -\left(\frac{\partial g}{\partial x}\right)^{-1} \frac{\partial g}{\partial p}
- 目的関数の導関数を計算する：
\frac{dT}{dp} = \frac{\partial T}{\partial p} - \lambda^T \frac{\partial g}{\partial p}
- 近似ヘッセ行列を構築する（簡易版）：
H \approx \frac{\partial^2 T}{\partial p^2} + S^T \frac{\partial^2 T}{\partial x^2}S
- パラメータを更新する：
p_{k+1} = p_k - \alpha_k H^{-1}\frac{dT}{dp}^T
このアルゴリズムは以下の特徴を有する：
- Adjoint法による効率的な勾配計算
- ニュートン法的な二次収束性
- 大規模問題への適用可能性
- 実装の現実性
実装上の注意点として以下が挙げられる：
- 状態方程式とAdjoint方程式の効率的な解法の選択（線形・非線形システムによって異なる）
- ヘッセ行列の近似精度と計算コストのバランス（BFGS等の準ニュートン法の利用も検討）
- 適切なステップ幅 $\alpha_k$ の選択（ライン検索法やトラスト領域法）
- 収束判定基準の設定
- 拘束条件が非線形の場合、ヤコビアン $\frac{\partial g}{\partial x}$ の更新が必要