制約付き最適化・LP・凸最適化

概要

制約付き最適化は実数ベクトル $\mathbf{x} \in \mathbb{R}^n$ に対し

$$\min_{\mathbf{x}} f(\mathbf{x}) \quad \text{subject to} \quad \mathbf{g}(\mathbf{x}) \leq 0, \;\; \mathbf{h}(\mathbf{x}) = 0$$

を解く問題で、目的関数・制約の線形性・凸性・微分情報の有無に応じて多様なアルゴリズムが使い分けられる。 sangi は線形計画 (LP)、非線形最小二乗、境界制約 L-BFGS-B、ペナルティ法、拡張ラグランジュ法、凸最適化 (FISTA, ADMM) を optimization_nd.hpp / linear_programming.hpp / ConvexOptimization.hpp で提供する。

無制約最適化は Optimization-unconstrained、微分不要法は Optimization-derivative-free を参照。

線形計画 (LP)

LP は目的関数・制約がすべて線形である問題:

$$\min \mathbf{c}^T \mathbf{x} \quad \text{subject to} \quad A \mathbf{x} \;\{\leq, =, \geq\}\; \mathbf{b}, \quad \mathbf{x} \geq 0$$

輸送計画・生産計画・ポートフォリオ最適化・ネットワークフローなどの古典的応用がある。

シンプレックス法

Dantzig (1947) のシンプレックス法は実行可能領域 (凸多面体) の頂点を辿って最適頂点を探す。 各反復で:

1. 現在の基底解に対し縮約コストを計算 (双対変数経由)。
2. 負の縮約コストを持つ非基底変数を選び (入る変数)、上限制約から基底から出る変数を決める (pivot 規則: Bland 規則、最大改善規則等)。
3. 基底を更新して次の頂点へ移る。

最悪計算量は指数 (Klee-Minty cube)、平均的には多項式時間で動く実用的アルゴリズム。 改訂シンプレックス法は基底逆行列を陽に保持せず LU 分解を更新することで疎な大規模問題でも効率的に動作する。

2 フェーズ法

実行可能基底解が直接見つからない場合、フェーズ I で人工変数を導入した補助 LP を解いて初期実行可能解を取得し、フェーズ II で本来の目的を最小化する。 sangi の simplex は 2 フェーズ改訂シンプレックス法を採用し、$\leq, =, \geq$ すべての制約型に対応する。

双対性

主問題 $\min \mathbf{c}^T \mathbf{x}$ s.t. $A \mathbf{x} \leq \mathbf{b}, \mathbf{x} \geq 0$ に対し 双対問題 $\max \mathbf{b}^T \mathbf{y}$ s.t. $A^T \mathbf{y} \leq \mathbf{c}, \mathbf{y} \geq 0$ が定義され、強双対定理から両者の最適値は一致する。 感度分析や経済学的解釈 (シャドープライス) に重要。

内点法と凸計画

シンプレックス法が境界を辿るのに対し、内点法は実行可能領域の内側を通って最適点に近づく。

主双対パスフォロワー

不等式制約 $\mathbf{x} \geq 0$ を対数バリア $-\mu \sum_i \log x_i$ で柔軟化:

$$\min \mathbf{c}^T \mathbf{x} - \mu \sum_i \log x_i \quad \text{subject to} \quad A \mathbf{x} = \mathbf{b}$$

$\mu \to 0$ の極限で元問題の解に収束する。中心パス上を Newton 法で近似的に追跡する。 Karmarkar (1984) で多項式時間が確立、Mehrotra (1992) の predictor-corrector 法で実用性能が大きく向上した。 密で大規模な LP では内点法が圧倒的に速い。

凸計画と SOCP / SDP

内点法は LP の枠を超えて、より一般の凸計画にも拡張される:

• QP (二次計画): $\min \tfrac{1}{2} \mathbf{x}^T Q \mathbf{x} + \mathbf{c}^T \mathbf{x}$, $A \mathbf{x} \leq \mathbf{b}$。SVM 等に多用。
• SOCP (二次錐計画): $\|A_i \mathbf{x} + \mathbf{b}_i\| \leq \mathbf{c}_i^T \mathbf{x} + d_i$ の形の二次錐制約。ロバスト LP・ポートフォリオ最適化。
• SDP (半正定値計画): 行列変数の半正定値制約。LMI 制御・組合せ緩和。

これらは self-concordant barrier の理論 (Nesterov-Nemirovski 1994) で統一的に扱える。 sangi 現バージョンは LP のシンプレックス法のみを提供し、QP / SOCP / SDP は今後の拡張対象である。

KKT 条件とラグランジュ乗数

等式制約 $\mathbf{h}(\mathbf{x}) = 0$ と不等式制約 $\mathbf{g}(\mathbf{x}) \leq 0$ のもとでの最適化問題に対し、 ラグランジュ関数

$$\mathcal{L}(\mathbf{x}, \boldsymbol{\lambda}, \boldsymbol{\mu}) = f(\mathbf{x}) + \boldsymbol{\lambda}^T \mathbf{h}(\mathbf{x}) + \boldsymbol{\mu}^T \mathbf{g}(\mathbf{x})$$

を導入する。$\mathbf{x}^*$ が局所最小点で制約が適切な regularity (LICQ 等) を満たすとき、 ある $\boldsymbol{\lambda}^*, \boldsymbol{\mu}^*$ が存在して以下の KKT 条件を満たす:

凸性 ($f, g_i$ 凸, $h_j$ アフィン) のもとでは KKT は十分条件にもなる。 多くの制約付き最適化アルゴリズム (SQP、内点法、active set) は KKT 系を反復的に解く。

SQP と active set 法

SQP (Sequential Quadratic Programming)

非線形制約付き最適化を反復的に QP 部分問題に置き換えて解く手法。各反復で:

$$\min_{\mathbf{p}} \; \nabla f(\mathbf{x}_k)^T \mathbf{p} + \tfrac{1}{2} \mathbf{p}^T H_k \mathbf{p}$$ $$\text{s.t.} \;\; \nabla h_j(\mathbf{x}_k)^T \mathbf{p} + h_j(\mathbf{x}_k) = 0, \;\; \nabla g_i(\mathbf{x}_k)^T \mathbf{p} + g_i(\mathbf{x}_k) \leq 0$$

$H_k$ はラグランジュ関数の Hessian (またはその BFGS 近似)。 二次収束する強力な手法で、SNOPT, NLPQL 等の商用ソルバの基盤。sangi 現バージョンでは未提供 (拡張候補)。

active set 法

不等式制約のうち現在の解で等号成立しているもの (active set) を仮定し、その上で等式制約問題を解く。 Lagrange 乗数の符号や KKT 違反で active set を更新する。 QP では古典的、$n$ が小〜中規模で疎な制約構造を持つ問題に有効。

ペナルティ法と拡張ラグランジュ

ペナルティ法

制約違反を目的関数のペナルティ項に変換し、無制約問題として解く:

$$\Phi_\mu(\mathbf{x}) = f(\mathbf{x}) + \mu \sum_i \max(0, g_i(\mathbf{x}))^2 + \mu \sum_j h_j(\mathbf{x})^2$$

$\mu \to \infty$ で制約満足に近づく。実装が単純で勾配ベース無制約ソルバを再利用できるが、 $\mu$ が大きいと条件数が悪化して内部最適化が困難になる。

sangi の penaltyMinimize は不等式制約版で、$\mu$ を段階的に増加させる。

拡張ラグランジュ法

ペナルティ項にラグランジュ乗数推定を加えた目的関数:

$$\mathcal{L}_A(\mathbf{x}, \boldsymbol{\lambda}, \boldsymbol{\mu}; \rho) = f(\mathbf{x}) + \boldsymbol{\lambda}^T \mathbf{h}(\mathbf{x}) + \tfrac{\rho}{2} \|\mathbf{h}(\mathbf{x})\|^2 + \cdots$$

外側ループで $\boldsymbol{\lambda} \leftarrow \boldsymbol{\lambda} + \rho \, \mathbf{h}(\mathbf{x}_k)$ と更新することで、$\rho$ を中程度に保ったまま制約満足を達成できる。 Powell-Hestenes-Rockafellar 法とも呼ばれる。 ペナルティ法の ill-conditioning を回避し、SQP と並ぶ標準的手法である。

sangi の augmentedLagrangianMinimize は等式・不等式両対応の実装。

境界制約 (L-BFGS-B)

各変数に独立な上下限 $l_i \leq x_i \leq u_i$ を持つ問題は、勾配射影法と L-BFGS の組み合わせで効率的に解ける。 Byrd, Lu, Nocedal, Zhu (1995) の L-BFGS-B が標準実装。

各反復で:

1. Cauchy point を勾配射影で計算し、active set を確定。
2. 自由変数 ($l_i < x_i < u_i$) のみで L-BFGS ステップを計算。
3. 射影ライン探索で境界を尊重しつつ更新。

パラメータが物理的に有界 (確率は $[0, 1]$ 等) な問題で広く使われる。 無制約変数には $l_i = -\infty, u_i = +\infty$ を渡せばよい。sangi の lbfgsbMinimize がこの実装。

非線形最小二乗 (Gauss-Newton, LM)

残差関数 $\mathbf{r}: \mathbb{R}^n \to \mathbb{R}^m$ ($m \geq n$) に対し $\min \tfrac{1}{2} \|\mathbf{r}(\mathbf{x})\|^2$ を解く問題。 曲線フィッティング・ロボット工学・コンピュータビジョン (bundle adjustment) で頻出する。

Gauss-Newton

残差の線形化 $\mathbf{r}(\mathbf{x} + \mathbf{p}) \approx \mathbf{r}(\mathbf{x}) + J \mathbf{p}$ から正規方程式:

$$J^T J \, \mathbf{p} = -J^T \mathbf{r}$$

$J$ がフルランクなら 1 反復ごとに大幅な減少を期待できる。 Hessian の二次項 $\sum r_i \nabla^2 r_i$ を無視するため、$\mathbf{r}$ が小さい近傍では Newton 法と等価になる。 $J^T J$ が悪条件・特異な場合に発散しやすい。

Levenberg-Marquardt

正規方程式に damping を加える:

$$(J^T J + \lambda I) \mathbf{p} = -J^T \mathbf{r}$$

$\lambda$ が大きいと最急降下方向に近づき、小さいと Gauss-Newton 方向に近づく。 gain ratio $\rho = \frac{\Delta f_{\text{actual}}}{\Delta f_{\text{predicted}}}$ で $\lambda$ を信頼領域的に更新する (Marquardt 1963)。

sangi の gaussNewtonMinimize は LM 正則化付き Gauss-Newton (固定 $\lambda$ + Armijo ライン探索)、 leastSquaresFit は gain ratio で適応的に $\lambda$ を制御する標準 LM 実装。 curveFit はモデル $y = f(x, \mathbf{p})$ のフィッティングを高レベル API で提供する。

共分散と統計

収束後、パラメータ共分散行列は $\mathrm{Cov}(\hat{\mathbf{p}}) = s^2 (J^T J)^{-1}$ で推定される ($s^2 = \|\mathbf{r}\|^2 / (m - n)$)。 各パラメータの標準誤差・$\chi^2$・reduced $\chi^2$ も LeastSquaresFitResult から取得できる。

FISTA と ADMM

機械学習・信号処理で頻出する非滑らかな凸最適化 (L1 正則化、TV 正則化等) には専用の高速アルゴリズムがある。

FISTA

Beck and Teboulle (2009) の Fast ISTA。次の合成問題を扱う:

$$\min_{\mathbf{x}} \; f(\mathbf{x}) + g(\mathbf{x})$$

$f$ は微分可能凸 (Lipschitz 定数 $L$ の勾配)、$g$ は近接作用素 $\mathrm{prox}_{\alpha g}(v) = \arg\min_u \{g(u) + \frac{1}{2\alpha} \|u - v\|^2\}$ が安価に計算できる凸関数。

反復式:

$$\mathbf{y}_k = \mathbf{x}_k + \frac{t_{k-1} - 1}{t_k} (\mathbf{x}_k - \mathbf{x}_{k-1})$$ $$\mathbf{x}_{k+1} = \mathrm{prox}_{(1/L) g}\bigl(\mathbf{y}_k - \tfrac{1}{L} \nabla f(\mathbf{y}_k)\bigr)$$ $$t_{k+1} = \tfrac{1 + \sqrt{1 + 4 t_k^2}}{2}$$

Nesterov の運動量で $O(1/k^2)$ 収束を達成。Lasso ($g = \lambda \|\cdot\|_1$, prox = soft thresholding)、TV 復元、低ランク行列復元等で標準。

sangi の fista が一般 FISTA、lassoFISTA が Lasso 特化版 (Lipschitz $L$ をパワーメソッドで自動推定)。

ADMM

Glowinski-Marrocco (1975) と Gabay-Mercier (1976) に遡る分割技法。 次の形を扱う:

$$\min \; f(\mathbf{x}) + g(\mathbf{z}) \quad \text{subject to} \quad A \mathbf{x} + B \mathbf{z} = \mathbf{c}$$

拡張ラグランジュを $\mathbf{x}, \mathbf{z}$ で交互最小化、双対変数 $\mathbf{y}$ を更新:

$$\mathbf{x}_{k+1} = \arg\min_{\mathbf{x}} \; f(\mathbf{x}) + \tfrac{\rho}{2} \|A \mathbf{x} + B \mathbf{z}_k - \mathbf{c} + \mathbf{u}_k\|^2$$ $$\mathbf{z}_{k+1} = \arg\min_{\mathbf{z}} \; g(\mathbf{z}) + \tfrac{\rho}{2} \|A \mathbf{x}_{k+1} + B \mathbf{z} - \mathbf{c} + \mathbf{u}_k\|^2$$ $$\mathbf{u}_{k+1} = \mathbf{u}_k + A \mathbf{x}_{k+1} + B \mathbf{z}_{k+1} - \mathbf{c}$$

変数分離で複雑な問題が 2 つの近接問題に分解されるため、$f, g$ の構造を活用できる。 分散最適化・consensus 最適化・複合正則化に強い。sangi の admmLasso は Lasso 特化版。

FISTA / ADMM の包括的解説は Boyd et al. (2011) の Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers を参照。

sangi の dispatch

スケーリングの重要性

制約付き最適化では変数・制約のスケールが結果に大きく影響する。 変数を $[-1, 1]$ 程度に正規化、制約も同オーダーに揃えると ill-conditioning が大きく改善する。 LP では並列化された行・列スケーリング、非線形最小二乗では Jacobian 列正規化が標準的対策。

参考文献

• Nocedal, J. and Wright, S. J. (2006). Numerical Optimization, 2nd ed. Springer.
• Boyd, S. and Vandenberghe, L. (2004). Convex Optimization. Cambridge University Press.
• Bertsekas, D. P. (1999). Nonlinear Programming, 2nd ed. Athena Scientific.
• Marquardt, D. W. (1963). "An Algorithm for Least-Squares Estimation of Nonlinear Parameters". SIAM J. Appl. Math., 11, 431–441.
• Beck, A. and Teboulle, M. (2009). "A Fast Iterative Shrinkage-Thresholding Algorithm for Linear Inverse Problems". SIAM J. Imaging Sci., 2, 183–202.
• Boyd, S., Parikh, N., Chu, E., Peleato, B., and Eckstein, J. (2011). "Distributed Optimization and Statistical Learning via the Alternating Direction Method of Multipliers". Foundations and Trends in Machine Learning, 3, 1–122.