scipy.optimize.fsolve の中身を考える
概要
SciPy
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.
fsolve
の公式ドキュメントには以下の記述がある.fsolve is a wrapper around MINPACK’s hybrd and hybrj algorithms.
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
の公式ドキュメントにも細かな説明が見えたので,後日にでもあらためることにする.
今日のところはこのあたりで筆を置く.
Author And Source
この問題について(scipy.optimize.fsolve の中身を考える), 我々は、より多くの情報をここで見つけました https://qiita.com/Yuu-Miino/items/c38449a1c6e03e8b5c89著者帰属:元の著者の情報は、元のURLに含まれています。著作権は原作者に属する。
Content is automatically searched and collected through network algorithms . If there is a violation . Please contact us . We will adjust (correct author information ,or delete content ) as soon as possible .