1. ニュートン法の基本原理
非線形方程式 $F(x) = 0$ を反復で解く。反復回数 $k$ での解の推定値を $x_k \in \mathbb{R}^n$、ヤコビ行列を $F'(x_k) \in \mathbb{R}^{n \times n}$ とすると、ニュートン法は以下のように定式化される：
x_{k+1} = x_k - [F'(x_k)]^{-1}F(x_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系）に対してニュートン法を適用する。$\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$ として：
\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}
2. 制約付き最適化問題の定式化
十分に滑らかな目的関数 $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$）とする。最小化問題は：
\min_{x \in \mathbb{R}^n} f(x)
subject to:
g(x) = 0, \quad g: \mathbb{R}^n \to \mathbb{R}^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} \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}\frac{\partial x}{\partial p} + \frac{\partial g}{\partial p} = 0
$\frac{\partial g}{\partial x}$ は行フルランク（$m \leq n$）なので、その擬似逆行列 $\left(\frac{\partial g}{\partial x}\right)^{+}$ を用いて感度行列は：
S = \frac{\partial x}{\partial p} = -\left(\frac{\partial g}{\partial x}\right)^{+} \frac{\partial g}{\partial p}
以下では $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}$ は、$p$ が $f$ に直接与える影響 $\frac{\partial f}{\partial p} \in \mathbb{R}^{1 \times k}$ と、$p$ が $x$ を通して間接的に $f$ に与える影響の和になる：
\frac{df}{dp} = \frac{\partial f}{\partial p} + \left(\frac{\partial x}{\partial p}\right)^T \frac{\partial f}{\partial x}
以下では、$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} =: A^{-T}$ より：
\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}
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$）、感度行列を経由すると $k$ 本の線形方程式が要るが、Adjoint法は $\lambda$ の1本で勾配が求まる。
10. トポロジー最適化への応用
コンプライアンス最小化問題の例を考察する。変位場を $u \in \mathbb{R}^n$、密度場を $\rho \in \mathbb{R}^m$、外力ベクトルを $F \in \mathbb{R}^n$、密度に直接依存する正則化項を $R(\rho)$ とする：
f(u, \rho) = \int_\Omega F^T u \, d\Omega + \alpha R(\rho)
対称正定値の剛性マトリクスを $K(\rho) \in \mathbb{R}^{n \times n}$ とすると、線形弾性体の平衡方程式は：
g(u,\rho) = K(\rho)u - F = 0
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)u = F$ と一致し、$\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{\partial f}{\partial \rho} = \alpha \frac{\partial R(\rho)}{\partial \rho}$、感度行列 $\frac{\partial u}{\partial \rho}$ は陰関数定理から得られる：
\frac{df}{d\rho} = \frac{\partial f}{\partial \rho} + \left(\frac{\partial u}{\partial \rho}\right)^T \frac{\partial f}{\partial u}
Adjoint法により、感度行列を明示的に計算することなく $\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{\partial T}{\partial x} \in \mathbb{R}^{1 \times n}$, $\frac{\partial T}{\partial p} \in \mathbb{R}^{1 \times k}$ を用いて一階全導関数を表現すると：
\frac{dT}{dp} = \frac{\partial T}{\partial p} + \frac{\partial T}{\partial x}S
二階全導関数（ヘッセ行列）は連鎖律により次のように導出される：
\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}$ の更新が必要