Newton-like root-finding algorithm that does not use derivatives
数値解析 において 、 ステフェンセン法 は ヨハン・フレデリック・ステフェンセン にちなんで名付けられた 、数値 根を求める 反復法であり、 セカント法 や ニュートン法 に似ています 。 ステフェンセン法は 導関数 を使用せずに 二次 収束 を達成しますが、より一般的なニュートン法も二次収束しますが導関数を必要とします。セカント法は導関数を必要としませんが、収束速度は二次収束よりも遅くなります。
ステフェンセン法は、ステップごとに2回の関数評価を必要とするという欠点がある。一方、セカント法はステップごとに1回の評価で済む。そのため、 計算コスト の点では必ずしも最も効率的とは言えない。これは、それぞれの反復回数に依存する。ニュートン法もまた、ステップごとに関数とその導関数の2つの関数を評価する必要があり、その計算コストはせいぜいセカント法と同程度だが、最悪の場合、ステフェンセン法と同程度になる。ほとんどの関数において、導関数の計算コストは元の関数の計算コストと同程度であるため、通常はニュートン法とステフェンセン法の計算コストは同程度である。 [a]
ステフェンセン法は、エイトケンのデルタ二乗過程を 固定小数点反復 に適用し たものとして導出できる 。このように考えると、ステフェンセン法は、固定小数点が存在することが保証され、固定小数点反復がバナッハの固定小数点定理によって(遅くなる可能性はあるものの)収束することが保証される限り、一般 バナッハ空間における効率的な 固定小数点計算 に自然に一般化される 。
簡単な説明 ステフェンセン法の最も単純な公式は、実関数 の 零点を 求める際に用いられます 。つまり、次式 を満たす実数値を求めることです。 解の近傍では 、関数の導関数は [b] を正確に、あるいは非常に近い値で満たす必要があります。 一部の関数では、この条件が満たされなくてもステフェンセン法は機能しますが、そのような場合、開始値は 実際の解に 非常に 近い値でなければならないため 、解への収束が遅くなる可能性があります。後述するように、この方法の中間ステップのサイズを調整することで、これらのケースの一部で収束を改善できます。 f {\displaystyle f} x ⋆ {\displaystyle \ x_{\star }\ } f ( x ⋆ ) = 0 . {\displaystyle \ f(x_{\star })=0~.} x ⋆ , {\displaystyle \ x_{\star }\ ,} f ′ , {\displaystyle \ f'\ ,} − 1 < f ′ ( x ⋆ ) < 0 . {\displaystyle -1<f'(x_{\star })<0~.} x 0 {\displaystyle \ x_{0}\ } x ⋆ , {\displaystyle \ x_{\star }\ ,}
適切な初期値が与えられれば、以下の式を用いて 値の列 を生成できる。この式が成立する場合、列の各値は 前の値よりも解に非常に近くなる。現在のステップの値は 、式 [1] を用いて次のステップの値を生成する。 x 0 , {\displaystyle \ x_{0}\ ,} x 0 , x 1 , x 2 , … , x n , … {\displaystyle \ x_{0},\ \ x_{1},\ x_{2},\ \dots ,\ x_{n},\ \dots \ } x ⋆ {\displaystyle \ x_{\star }\ } x n {\displaystyle \ x_{n}\ } x n + 1 {\displaystyle \ x_{n+1}\ }
x n + 1 = x n − f ( x n ) g ( x n ) {\displaystyle x_{n+1}=x_{n}-{\frac {f(x_{n})}{g(x_{n})}}} ここで、 傾き関数は、 式で与えられる 元の関数の合成である。 n = 0 , 1 , 2 , 3 , . . . , {\displaystyle \ n=0,1,2,3,...\ ,} g ( x ) {\displaystyle \ g(x)\ } f {\displaystyle \ f\ }
g ( x ) = f ( x + f ( x ) ) f ( x ) − 1 {\displaystyle g(x)={\frac {f{\bigl (}x+f(x){\bigr )}}{f(x)}}-1} あるいはもっと明確に言えば、
g ( x ) = f ( x + h ) − f ( x ) h ≈ d f ( x ) d x ≡ f ′ ( x ) , {\displaystyle g(x)={\frac {f(x+h)-f(x)}{h}}\qquad \approx \quad {\frac {\operatorname {d} f(x)}{\operatorname {d} x}}\equiv f'(x),} ここで、最後の反復点 と補助点 との間のステップサイズは、 h ≡ f ( x ) {\displaystyle \ h\equiv f(x)\ } x , {\displaystyle \ x\ ,} x + h . {\displaystyle \ x+h~.}
技術的には、この関数は 2点間 の 1次 差と呼ばれます [c] 実際には、最後のシーケンスポイント と補助点の間 の関数の 傾きの平均値であり、 中間ステップのサイズ(およびその方向)は次のように与えられます。 g {\displaystyle \ g\ } f {\displaystyle \ f\ } f ′ {\displaystyle f'} f {\displaystyle \ f\ } ( x , y ) = ( x n , f ( x n ) ) {\displaystyle \left(x,y\right)={\bigl (}x_{n},f\left(x_{n}\right){\bigr )}} ( x , y ) = ( x n + h , f ( x n + h ) ) , {\displaystyle \ {\bigl (}x,y{\bigr )}={\bigl (}x_{n}+h,f\left(x_{n}+h\right){\bigr )}\ ,} h = f ( x n ) . {\displaystyle \ h=f(x_{n})~.}
の値は の 近似値であるため 、その値がステフェンセンのアルゴリズムの収束を保証するために必要な条件を満たすかどうかをオプションで確認することができます 。わずかな不適合は必ずしも重大な結果をもたらすとは限りませんが、条件から大きく逸脱した場合は、ステフェンセンの方法が失敗する可能性が高いことを警告しており、一時的に何らかのフォールバックアルゴリズム(例えば、より堅牢な イリノイアルゴリズム や、単純な regula falsi )を使用することが正当化されます。 g {\displaystyle \ g\ } f ′ , {\displaystyle \ f'\ ,} − 1 < g < 0 , {\displaystyle \ -1<g<0\ ,}
この補助点を求める目的のためだけに 、関数の値は [b] の要件を満たす必要がある。 計算の他の部分では、ステフェンセン法では関数が 連続であり、実際に近傍解を持つことのみが求められる。 [1] 傾きの式で使用される ステップに は、 h {\displaystyle \ h\ } f {\displaystyle \ f\ } − 1 < f ′ ( x ⋆ ) < 0 . {\displaystyle \ -1<f'(x_{\star })<0~.} f {\displaystyle \ f\ } h {\displaystyle \ h\ } g {\displaystyle \ g\ } 1 / 2 または 3 / 4 、 要件を完全に満たしていない機能に対応するためです。 f {\displaystyle \ f\ }
利点と欠点 ステフェンセン法の主な利点は、 ニュートン法と 同様に 二次収束性 [1] を示すことです。つまり、どちらの方法も方程式の根を 同じように「速く」求めます。この場合、「 速く 」とは、どちらの方法においても、解の正しい桁数がステップごとに倍増することを意味します。しかし、ニュートン法の式は関数だけでなくその導関数の評価も必要としますが、 ステフェンセン法では 関数 自身の評価のみで済みます。これは、導関数が容易に、あるいは効率的に得られない場合に重要です。 f {\displaystyle \ f\ } f ′ {\displaystyle \ f'\ } f , {\displaystyle \ f\ ,} f {\displaystyle \ f\ }
高速収束の代償として、関数評価が2回必要になります。つまり、と の両方 を計算する必要があり、が 複雑な場合は時間がかかる可能性があります。比較すると、 regula falsi 法と セカント法 はどちらも1ステップあたり1回の関数評価で済みます。セカント法は1ステップあたり正解桁数を「わずか」約1.6倍しか増やしませんが、一定時間内にセカント法の2倍のステップ数を実行できます。セカント法はSteffensen法と同じ時間で2倍のステップ数を実行できるため、 [d] 実用上、両方のアルゴリズムが成功した場合、セカント法の方がSteffensen法よりも速く収束します。セカント法は2ステップ(2回の関数評価)ごとに約 (1.6) 2 ≈ 2.6 倍 の桁数を達成しますが、Steffensen法は1ステップ(2回の関数評価)ごとに 2 倍の桁数を達成します。 f ( x n ) {\displaystyle \ f(x_{n})\ } f ( x n + h ) {\displaystyle \ f(x_{n}+h)\ } f {\displaystyle \ f\ }
他のほとんどの反復根探索アルゴリズム と同様に 、ステフェンセン法の決定的な弱点は「十分に近い」開始値を選択することです。 の値が 実際の解に「十分近く」ない場合 、この方法は失敗する可能性があり、値のシーケンスは 2 つ (またはそれ以上) の極端な値の間で不規則に反転したり、無限大に発散したり、またはその両方が発生する可能性があります。 x 0 . {\displaystyle \ x_{0}~.} x 0 {\displaystyle \ x_{0}\ } x ⋆ , {\displaystyle \ x_{\star }\ ,} x 0 , x 1 , x 2 , x 3 , … {\displaystyle \ x_{0},\,x_{1},\,x_{2},\,x_{3},\,\dots \ }
エイトケンのデルタ二乗過程を用いた導出 以下に示すMATLAB コードに実装されているSteffensen法のバージョンは、 収束加速 のための Aitkenのデルタ2乗法 を用いて見つけることができます 。以下の式を上のセクションの式と比較すると、 であることに留意してください 。この方法は、線形収束するシーケンスから開始することを前提とし、そのシーケンスの収束速度を高めます。 の符号が 一致し、 シーケンスの望ましい限界に「十分に近い」場合 、次のように仮定できます。 x n = p − p n {\displaystyle x_{n}=p-p_{n}} p n , p n + 1 , p n + 2 {\displaystyle p_{n},\,p_{n+1},\,p_{n+2}} p n {\displaystyle p_{n}} p {\displaystyle p}
p n + 1 − p p n − p ≈ p n + 2 − p p n + 1 − p , {\displaystyle {\frac {p_{n+1}-p}{p_{n}-p}}\approx {\frac {p_{n+2}-p}{p_{n+1}-p}},} となることによって
( p n + 2 − 2 p n + 1 + p n ) p ≈ p n + 2 p n − p n + 1 2 . {\displaystyle (p_{n+2}-2p_{n+1}+p_{n})p\approx p_{n+2}p_{n}-p_{n+1}^{2}.} シーケンスの望ましい限界を解くと 次のようになります。 p {\displaystyle p}
p ≈ p n + 2 p n − p n + 1 2 p n + 2 − 2 p n + 1 + p n {\displaystyle p\approx {\frac {p_{n+2}p_{n}-p_{n+1}^{2}}{p_{n+2}-2p_{n+1}+p_{n}}}} = ( p n 2 + p n p n + 2 − 2 p n p n + 1 ) − ( p n 2 − 2 p n p n + 1 + p n + 1 2 ) p n + 2 − 2 p n + 1 + p n {\displaystyle =~{\frac {\,(\,p_{n}^{2}+p_{n}\,p_{n+2}-2\,p_{n}\,p_{n+1}\,)-(\,p_{n}^{2}-2\,p_{n}\,p_{n+1}+p_{n+1}^{2}\,)\,}{\,p_{n+2}-2\,p_{n+1}+p_{n}\,}}} = p n − ( p n + 1 − p n ) 2 p n + 2 − 2 p n + 1 + p n , {\displaystyle =p_{n}-{\frac {(p_{n+1}-p_{n})^{2}}{p_{n+2}-2p_{n+1}+p_{n}}},} その結果、より急速に収束するシーケンスが得られます。
p ≈ p n + 3 = p n − ( p n + 1 − p n ) 2 p n + 2 − 2 p n + 1 + p n . {\displaystyle p\approx p_{n+3}=p_{n}-{\frac {(p_{n+1}-p_{n})^{2}}{p_{n+2}-2p_{n+1}+p_{n}}}.}
コード例
Matlabで 以下は、 MATLAB での Steffensen メソッドの実装のソースです 。
function Steffensen ( f, p0, tol ) % この関数は、固定小数点反復関数 f、 固定小数点への初期推定値 p0、および許容値 tol を入力として受け取ります。 % 固定小数点反復関数は、インライン関数として入力されると想定されています 。 % この関数は、式 f(x) = p が目的の 許容値 tol 内で真となる 固定小数点 p を計算して返します。 format compact % 出力を短くします。 format long % 小数点以下の桁数を多く出力します。 for i = 1 : 1000 % 大規模だが有限な回数の反復処理を実行する準備をします。 % これは、メソッドが収束に失敗した場合、 無限ループに陥らないようにするためです。 p1 = f ( p0 ) + p0 ; % 固定点の次の 2 つの推定値を計算します。 p2 = f ( p1 ) + p1 ; p = p0 - ( p1 - p0 ) ^ 2 / ( p2 - 2 * p1 + p0 ) % Aitken のデルタ 2 乗法を使用して、 % p0 のより適切な近似値を見つけます。 if abs ( p - p0 ) < tol % 許容範囲内かどうかをテストします。 break % 範囲内であれば、反復処理を停止します。答えが得られます。 end p0 = p ; % 次の反復処理のために p0 を更新します。 end if abs ( p - p0 ) > tol % 許容値を満たさない場合は、 失敗のメッセージを出力します。 % 「1000回の反復で収束できませんでした。」 終了
Pythonの場合 以下はPython での Steffensen メソッドの実装のソースです 。
import Callable , Iterator Func = Callable [[ float ], float , float ] と 入力して def g ( f : Func , x : float , fx : float ) -> Func : """1 階差分商関数。 引数: f: gへの関数入力 x: gを評価する点 fx: xで評価される関数f """ return f ( x + fx ) / fx - 1 def steff ( f : Func , x : float , tol : float ) -> Iterator [ float ]: """根を見つけるためのステッフェンセンアルゴリズム。 この再帰ジェネレータは最初に x_{n+1} 値を生成し、次にジェネレータが反復処理されるときに、 次の再帰レベルから x_{n+2} を生成します。 引数: f: ルートを検索する関数 x: 最初の呼び出し時の開始値、関数が再帰する各レベル n x は x_n """ n = 0 while True : if n > 1000 : print ( "1000回の反復で収束しませんでした" ) break else : n = n + 1 fx = f ( x ) if abs ( fx ) < tol : break else : gx = g ( f , x , fx ) x = x - fx / gx # x_{n+1} に更新 yield x # 値を返す
バナッハ空間への一般化 ステフェンセン法は、入力と同じ出力を生成する 異なる種類の関数の 入力を求めるのにも使用できます。 例えば、特殊な値 のような解は 固定点 と呼ばれます 。これらの関数の多くは、結果を繰り返し入力として再利用することで、自身の解を求めることができますが、収束速度が遅くなったり、関数によっては収束に全く至らない場合があります。ステフェンセン法は、この収束を加速し、 2次収束 へと導きます。 x = x ⋆ {\displaystyle \ x=x_{\star }\ } F {\displaystyle \ F\ } x ⋆ = F ( x ⋆ ) {\displaystyle \ x_{\star }=F(x_{\star })\ } x ⋆ . {\displaystyle \ x_{\star }~.} x ⋆ {\displaystyle \ x_{\star }\ }
例として、より一般的な バナッハ空間 と基本 実数 の問題を一時的に無視します。読者を前のセクションに再び向けると、 任意のルート関数を使用した 、固定小数点関数の単純な おもちゃのモデルは 、次のように作成できます。 ここでは、反復 処理で安定するのに十分小さい値でありながら、 関数の 非線形性 が顕著になる
のに十分な 大きさの適切な符号を持つ定数です。 F ~ , {\displaystyle \ {\tilde {F}}\ ,} f , {\displaystyle \ f\ ,} F ~ ( x ) = x + ε f ( x ) . {\displaystyle \ {\tilde {F}}(x)=x+\varepsilon \ f(x)~.} ε {\displaystyle \ \varepsilon \ } F ~ {\displaystyle \ {\tilde {F}}\ } f {\displaystyle \ f\ }
実数値関数 の不動点を求めるこの方法は、 バナッハ空間を それ自身に写す関数 、あるいはより一般的には、ある バナッハ空間 から別の バナッハ空間 に写す 関数に対して一般化されている。 この一般化された方法は、 および に関連付けられた 有界 線形作用素 の 族が、(局所的に)条件 [2] を満たすように考案できることを前提としている。 F : X → X {\displaystyle \ F:X\to X\ } X {\displaystyle \ X\ } F : X → Y {\displaystyle \ F:X\to Y\ } X {\displaystyle X} Y . {\displaystyle \ Y~.} { G ( u , v ) : u , v ∈ X } {\displaystyle \ {\bigl \{}\ G(u,v):u,v\in X\ {\bigr \}}\ } u {\displaystyle \ u\ } v {\displaystyle \ v\ }
F ( u ) − F ( v ) = G ( u , v ) ( u − v ) {\displaystyle F\left(u\right)-F\left(v\right)=G\left(u,v\right)\ {\bigl (}\ u-v\ {\bigr )}\quad } 1
演算子は 、すべての要素が ベクトル 引数と の関数である 行列 とほぼ等価です 。最初のセクションで示した 単純な関数 をもう一度参照してください 。この関数は実数を単に入力して出力するだけです。ここでは、関数は 商差 です 。ここでの一般化された形式では、演算子は バナッハ空間 で使用される商差の類似物です 。 G {\displaystyle \ G\ } u {\displaystyle \ u\ } v . {\displaystyle \ v~.} f , {\displaystyle \ f\ ,} g {\displaystyle \ g\ } G {\displaystyle \ G\ }
バナッハ空間 で除算が可能な場合 、線形演算子は 次のように得られる。 G {\displaystyle \ G\ }
G ( u , v ) = [ F ( u ) − F ( v ) ] ( u − v ) − 1 , {\displaystyle G\left(u,v\right)={\bigl [}\ F\left(u\right)-F\left(v\right)\ {\bigr ]}\ {\bigl (}\ u-v\ {\bigr )}^{-1}\ ,} これは、いくつかの洞察を与えるかもしれません。このように表現すると、線形演算子は、上記の最初のセクションで議論した 商差 の複雑なバージョンであることがより容易にわかります 。商形式は、ここでは説明のためにのみ示されており、それ 自体 は必須では ありません 。また、バナッハ空間内での除算は、この詳細なステフェンセン法の実行に必ずしも必要ではないことにも注意してください。唯一の要件は、演算子が ( 1 )を満たすことです。 G {\displaystyle \ G\ } g {\displaystyle \ g\ } G {\displaystyle \ G\ }
ステフェンセン法は、微分の代わりに 差分商を使用する点を除けば、ニュートン法と非常によく似ています。固定 点関数 とそれらの線形演算子が 条件( 1 )を満たす いくつかの固定点に近い引数については、恒等演算子 である ことに注意してください 。 G ( F ( x ) , x ) {\displaystyle \ G{\bigl (}F\left(x\right),x{\bigr )}\ } F ′ ( x ) . {\displaystyle \ F'(x)~.} x {\displaystyle \ x\ } x ⋆ , {\displaystyle \ x_{\star }\ ,} F {\displaystyle \ F\ } G {\displaystyle \ G\ } F ′ ( x ) ≈ G ( F ( x ) , x ) ≈ I , {\displaystyle \ F'(x)\ \approx \ G{\bigl (}F\left(x\right),x{\bigr )}\ \approx \ I\ ,} I {\displaystyle \ I\ }
バナッハ空間で分割が可能な場合、一般化された反復公式は次のように与えられる。
x n + 1 = x n + [ I − G ( F ( x n ) , x n ) ] − 1 [ F ( x n ) − x n ] , {\displaystyle x_{n+1}=x_{n}+{\Bigl [}\ I-G{\bigl (}F\left(x_{n}\right),x_{n}{\bigr )}\ {\Bigr ]}^{-1}{\Bigl [}\ F\left(x_{n}\right)-x_{n}\ {\Bigr ]}\ ,} より一般的な場合 、 割り算が不可能な場合、反復公式は 、 n = 1 , 2 , 3 , . . . . {\displaystyle \ n=1,\ 2,\ 3,\ ...~.} x n + 1 {\displaystyle \ x_{n+1}\ } x n {\displaystyle \ x_{n}\ }
[ I − G ( F ( x n ) , x n ) ] ( x n + 1 − x n ) = F ( x n ) − x n . {\displaystyle {\Bigl [}\ I-G{\bigl (}F\left(x_{n}\right),x_{n}{\bigr )}\ {\Bigr ]}{\bigl (}\ x_{n+1}-x_{n}\ {\bigr )}=F\left(x_{n}\right)-x_{n}~.} 同様に、ある程度簡約された形の 解を求めることもできる。 x n + 1 {\displaystyle \ x_{n+1}\ }
[ I − G ( F ( x n ) , x n ) ] x n + 1 = [ F ( x n ) − G ( F ( x n ) , x n ) x n ] , {\displaystyle {\Bigl [}\ I-G{\bigl (}F\left(x_{n}\right),x_{n}{\bigr )}\ {\Bigr ]}\ x_{n+1}={\Bigl [}\ F\left(x_{n}\right)-G{\bigl (}F\left(x_{n}\right),x_{n}{\bigr )}\ x_{n}\ {\Bigr ]}\ ,} 角括弧内の値はすべて に依存しない。 括弧内の項はすべて のみに依存する 。しかし、2番目の形式は最初の形式ほど 数値的に安定して いない可能性がある。最初の形式は(願わくば)小さな差の値を求めるため、反復値の過度に大きな変化や不規則な変化を回避する可能性が数値的に高い可能性がある。 x n + 1 : {\displaystyle \ x_{n+1}\ :} x n {\displaystyle \ x_{n}\ } x n . {\displaystyle \ x_{n}~.}
線形演算子 が G {\displaystyle \ G\ }
‖ G ( u , v ) − G ( x , y ) ‖ ≤ k ( ‖ u − x ‖ + ‖ v − y ‖ ) {\displaystyle {\Bigl \|}G\left(u,v\right)-G\left(x,y\right){\Bigr \|}\leq k{\biggl (}{\Bigl \|}u-x{\Bigr \|}+{\Bigr \|}v-y{\Bigr \|}{\biggr )}} ある正の実定数に対して、 初期近似値が次の式を 満たす 所望の解に「十分近い」場合 、この方法は2乗的に固定点に収束する。 k , {\displaystyle \ k\ ,} F {\displaystyle \ F\ } x 0 {\displaystyle \ x_{0}\ } x ⋆ {\displaystyle \ x_{\star }\ } x ⋆ = F ( x ⋆ ) . {\displaystyle \ x_{\star }=F(x_{\star })~.}
注記 ^ 稀な特殊なケースの関数については、主関数の評価から節約された部分を用いることで、ニュートン法の微分計算は無視できるコストで実行できます。このように最適化すると、ニュートン法はセカント法よりもステップあたりのコストがわずかに高くなるだけで、収束速度がわずかに速くなるという利点があります。 ^ ab 条件は、 を 自身の 解を求める ための補正関数として使用し た 場合 、 が解の方向( )に移動し、新しい値は解と以前の値( )の間に位置する傾向があることを保証する。ただし、 は 原理上は 自己補正関数に過ぎない点に注意する必要がある 。実際にはそのような目的で使用されることはなく、たとえ使用されるとしても効率的である必要はない。 − 1 < f ′ ( x ⋆ ) < 0 {\displaystyle -1<f'(x_{\star })<0\ } f {\displaystyle \ f\ } x , {\displaystyle \ x\ ,} f ′ < 0 {\displaystyle \ f'<0\ } − 1 < f ′ {\displaystyle -1\ <f'\ } f {\displaystyle \ f\ } ^ 商差は、 の 符号に応じて、 前方商 差または 後方 商 差の いずれかになります 。 g {\displaystyle \ g\ } h {\displaystyle \ h\ } ^ 2つの評価の事前計算を順番に実行する必要がある ため 、関数評価を並列に実行してもアルゴリズム 自体 を高速化することはできません。これはステフェンセン法のもう一つの欠点です。 f ( x n + h ) {\displaystyle \ f(x_{n}+h)\ } h ≡ f ( x n ) , {\displaystyle \ h\equiv f(x_{n})\ ,}
参考文献