概要

SciPy

Wikipedia（2021/10/06 現在）より抜粋

SciPyはPythonのための科学的ツールのオープンソース・ライブラリとして開発されている。SciPyは配列の高速な操作のためのすべてのライブラリを含んでおり、人気のNumericモジュールを置き換え、ひとつのパッケージとして高レベルな科学と工学のモジュールを集めたもの。

SciPyは、配列オブジェクトとその他の基本的な機能を備えたNumPyを基礎にしている。SciPy は統計、最適化、積分、線形代数、フーリエ変換、信号・イメージ処理、遺伝的アルゴリズム、ODE (常微分方程式) ソルバ、特殊関数、その他のモジュールを提供する。

scipy.optimize.fsolve

optimize は最適化を意味し，ここにはその名の通り最適化問題を解くためのツールがまとまっている．

その中の一つ，fsolve は与えられた（数学的）関数の根を求めるための（プログラムにおける）関数である．

公式ドキュメントに簡単な利用例として次のようなものが掲載されている．

from scipy.optimize import fsolve
def func(x):
    return [x[0] * np.cos(x[1]) - 4,
            x[1] * x[0] - x[1] - 5]
root = fsolve(func, [1, 1])
root
array([6.50409711, 0.90841421])
np.isclose(func(root), [0.0, 0.0])  # func(root) should be almost 0.0.
array([ True,  True])

特に説明の必要はなさそうだが，np は numpy のことである．

この例では，次の連立方程式を満たす $x, y$ を求めている．

$$
\left\{
\begin{array}{rcl}
x \cos y - y &=& 0 \newline
xy - y - 5 &=& 0
\end{array}\right.
$$

この記事の目指すところ

fsolve がどのような数学的根拠の元で設計されているのかを探る．
とくに，xtol によってコントロールされる収束判定条件について，条件は本当に有効であるのか，関数値を用いた誤差判定はできないのかに興味がある．

掘り下げ

fsolve の公式ドキュメントには以下の記述がある．

fsolve is a wrapper around MINPACK’s hybrd and hybrj algorithms.

MINPACK は1980年頃，プログラミング言語 FORTRAN の上で開発されたライブラリの一つで，非線形方程式の解の計算を可能とするものである．（Wikipedia 参照）
MINPACK の公式ドキュメント P. 57 には，hybrd の実装について以下の通り記述されている．

DeepL によって和訳した結果は以下の通り．

HYBRDの目的は、Powellハイブリッド法の修正により、N変数のN個の非線形関数の系のゼロを求めることである。ユーザーは関数を計算するサブルーチンを提供しなければならない。ヤコビアンは前方差分近似により計算される。

同ドキュメント P.71 には，hybrdj の実装について以下の通り記述されている．

すなわち，こちらの関数では Jacobi 行列も自前で用意して渡す必要がある．

いずれの関数も A modification of Powell hybrid method なる手法を実装していることがわかった．（MINPACK に記載の元論文：Powell, M. J. (1970). A hybrid method for nonlinear equations. Numerical methods for nonlinear algebraic equations.）
- Powell hybrid method は通称，Powell's dog leg method と呼ばれるもので，steepest descent method という大域的な収束性に特化した手法と，Broyden 法という局所的な収束性に特化した手法を並列的に組み合わせた（いいとこ取りをした）ものである．(H. S. Chen の論文 (1980) に説明あり: Chen, H. S., & Stadtherr, M. A. (1981). A modification of Powell's dogleg method for solving systems of nonlinear equations. Computers & Chemical Engineering, 5(3), 143-150.)
- Modification という点が公式ドキュメントでは不明瞭であるが，年代を考えても H. S. Chen の論文で提案されている手法とみなすのが自然か？とはいえ，氏は開発者でもない上に当該論文が Reference にすら含まれていないことを鑑みると，MINPACK のところでアルゴリズム的な細かい修正を入れたというだけの記述にも思える．
ひとまずは上記の論文あたりが，fsolve の妥当性の数学的根拠といえる．
ここまで調べてわかったことは，ある程度の精度で初期値が求まっているのなら，Broyden 法を用いるのみで事足りるということだ．無理にハイブリッドなやり方を使う必要は無いのである．

ここまで書いて思い出したこと

fsolve の正体に Broyden 法が含まれることを見て，「そういえば，SciPy には Broyden 法の実装もあったな」と思い出し，公式ドキュメントを覗いてみる．すると中のリンクから scipy.optimize.root へ飛べるようになっており（実は fsolve からも飛べたことはあとになって気づいた...），ふと，「そういえば fsolve と root は何が違うのだろう」と思い至る．調べてみると，Quora という掲示板？で以下のコメントを見つけた．

要は，「古い関数を残してあるだけなので，これからの人は新しいの（scipy.optimize.root）を使いましょう」ということだ．

もうここまで書いてしまったので，タイトルはそのままで公開するが，本記事の興味の対象は fsolve から root(method=="hybr") へと移った．（中身は同じであるし，これまでの調査で得られた事実は全く無駄にならない．）

root(method=="broyden1") の恩恵

root を使用するなら，用いるメソッドによっては xtol 以外にも様々な収束判定条件を付与できる．筆者は fsolve における xtol の条件よりも，root(method=="broyden1") で設定できる fatol を好ましく考えているため，この発見は大きな成果である．

そして何より，xtol の有効性如何はまだ検討できていないが，fatol が設定できることがわかった以上，当面をこちらを用いれば問題ないと考えるに至った．

xtol の謎は MINPACK の公式ドキュメントにも細かな説明が見えたので，後日にでもあらためることにする．

今日のところはこのあたりで筆を置く．