スカイフィールドで移動せずにscipy.optimize.minimize（デフォルト）が成功を報告するのはなぜですか？

Question

scipy.optimize.minimizeデフォルトメソッドを使用すると、エラーまたは警告メッセージなしで結果として初期値が返されます。 this answerにより示唆されるようネルダー・ミードの方法を使用して問題を解決するが、私は理解したいと思います：ん

は、デフォルトの方法は、答えとして出発点を警告なし間違った答えを返すのはなぜ - とあります私はを "警告なしで間違った答え"から守ることができますこの場合、この現象は避けてください。

注意、機能separationは、シンプレックスがここに優れている理由とすることができる、滑らかで保証されていない最小化する値を生成するためのpythonパッケージSkyfieldを使用しています。

結果：

試験結果：[2.14159739] '正しい'：初期2.14159265359：0.0

既定の結果：[] '正しい'：13054初期：

ネルダ - ミードの結果：[13053.81011963] '正しい'：初期：ここ

FULL OUTPUT using DEFAULT METHOD: 
    status: 0 
    success: True 
    njev: 1 
    nfev: 3 
hess_inv: array([[1]]) 
     fun: 1694.98753895812 
     x: array([ 10000.]) 
    message: 'Optimization terminated successfully.' 
     jac: array([ 0.]) 
     nit: 0 

FULL OUTPUT using Nelder-Mead METHOD: 
    status: 0 
    nfev: 63 
success: True 
    fun: 3.2179306044608054 
     x: array([ 13053.81011963]) 
message: 'Optimization terminated successfully.' 
    nit: 28 

10000がいっぱいスクリプトは次のとおりです。

def g(x, a, b): 
    return np.cos(a*x + b) 

def separation(seconds, lat, lon): 
    lat, lon, seconds = float(lat), float(lon), float(seconds) # necessary it seems 
    place = earth.topos(lat, lon) 
    jd = JulianDate(utc=(2016, 3, 9, 0, 0, seconds)) 
    mpos = place.at(jd).observe(moon).apparent().position.km 
    spos = place.at(jd).observe(sun).apparent().position.km 
    mlen = np.sqrt((mpos**2).sum()) 
    slen = np.sqrt((spos**2).sum()) 
    sepa = ((3600.*180./np.pi) * 
      np.arccos(np.dot(mpos, spos)/(mlen*slen))) 
    return sepa 

from skyfield.api import load, now, JulianDate 
import numpy as np 
from scipy.optimize import minimize 

data = load('de421.bsp') 

sun = data['sun'] 
earth = data['earth'] 
moon = data['moon'] 

x_init = 0.0 
out_g = minimize(g, x_init, args=(1, 1)) 
print "test result: ", out_g.x, "'correct': ", np.pi-1, "initial: ", x_init # gives right answer 

sec_init = 10000 
out_s_def = minimize(separation, sec_init, args=(32.5, 215.1)) 
print "default result: ", out_s_def.x, "'correct': ", 13054, "initial: ", sec_init 

sec_init = 10000 
out_s_NM = minimize(separation, sec_init, args=(32.5, 215.1), 
       method = "Nelder-Mead") 
print "Nelder-Mead result: ", out_s_NM.x, "'correct': ", 13054, "initial: ", sec_init 

print "" 
print "FULL OUTPUT using DEFAULT METHOD:" 
print out_s_def 
print "" 
print "FULL OUTPUT using Nelder-Mead METHOD:" 
print out_s_NM 

Answer 1

1）

あなたの関数は、区分的定数（ある小規模な "階段" パターンを有しています）。 どこでも微分可能ではありません。

最初の推測での関数の勾配はゼロです。

デフォルトのBFGSオプティマイザはゼロ勾配を見て、それが基準（微分可能性などのこの場合には当てはまらない入力関数に関する仮定に基づく）によって極小であると判断します。

基本的に、正確に平坦な領域はオプティマイザを爆撃します。オプティマイザは、最初の点のまわりの小さな正確に平坦な領域の関数を調べます。関数の定数はすべて定数であるため、関数に定数を与えたと考えます。あなたの関数はどこでも微分可能ではないので、ほとんどすべての初期点がそのような平坦な領域の内側にある可能性があります。そのため、初期点の選択に不運が生じることはありません。

Nelder-Meadはではなく、であることに注意してください。最初のシンプレックスは階段のサイズよりも大きいので、より大きなスポットの周りの関数を調べます。初期シンプレックスが階段サイズより小さい場合、オプティマイザはBFGSと同様に動作します。

2）

一般的な答えは：地元のオプティマイザは、局所最適解を返します。これらが真の最適値と一致するかどうかは、関数の性質に依存する。

一般的に、あなたがローカルの最適状態に陥っていないかどうかを確認するには、初期の推測値を変えてみてください。

また、微分不可能な関数で微分ベースのオプティマイザを使用することはお勧めできません。関数が「大」スケールで微分可能な場合は、数値微分のステップサイズを調整できます。

関数がどこでも微分可能であるかどうか数値的にチェックする安価で一般的な方法はないので、そのようなチェックは行われません。代わりに、目的関数を入力した者が保証する必要がある最適化方法最適化方法を選択します。

Answer 2

@pvによる受け入れ可能な回答。スカイフィールドには「階段」の応答があると説明されています。つまり、スカラーフィールドから返される値の一部は離散ジャンプを除いてローカルにフラットです。

私は最初のステップで時間をJulianDateオブジェクトに変換して少し実験をしましたが、実際には1増分あたり約40マイクロ秒、つまり約5E-10日間に見えます。 JPLデータベースが何千年もの間に広がっていることを考えると、合理的です。これはほぼすべての一般的な天文学的規模のアプリケーションではおそらく問題ありませんが、実際には滑らかではありません。答えが指摘しているように、局所平坦性は、（おそらく多くの）ミニマイザで「成功」を引き起こします。これは予期され、合理的であり、決してその方法の失敗ではない。

from skyfield.api import load, now, JulianDate 
import numpy as np 
import matplotlib.pyplot as plt 

t = 10000 + np.logspace(-10, 2, 25)  # logarithmic spacing 
jd = JulianDate(utc=(2016, 3, 9, 0, 0, t)) 

dt = t[1:] - t[:-1] 
djd = jd.tt[1:] - jd.tt[:-1] 

t = 10000 + np.linspace(0, 0.001, 1001)  # linear spacing 
jd = JulianDate(utc=(2016, 3, 9, 0, 0, t)) 

plt.figure() 

plt.subplot(1,2,1) 

plt.plot(dt, djd) 
plt.xscale('log') 
plt.yscale('log') 

plt.subplot(1,2,2) 

plt.plot(t, jd.tt-jd.tt[0]) 

plt.show() 

Answer 3

私はアルゴリズムが時間を通じてどのように動作しているか見るためにあまりにも非常にprint文の値をほめることはできません。

def separation(seconds, lat, lon): 
    print seconds 
    ... 

は、この行を追加すると、あなたはネルダー・ミードの方法がないことを確認できるようになります：あなたはあなたのseparation()関数の先頭に1を追加してみてください場合は、最小化ルーチンが答えに向けて自分の道を働く見るようになるだろう500秒単位で前方跨ぐ秒の範囲の徹底的な検索は、それが再生を開始する前に：それはこれらが原因でこのようなソルバに、500秒ごとです知らない、もちろん

[ 10000.] 
[ 10500.] 
[ 11000.] 
[ 11500.] 
[ 12500.] 
... 

、問題には単位がありません。これらの調整は、500メートル、または500オングストローム、または500年とすることができる。しかし、それは盲目的に前方につまずき、Nelder-Meadの場合は、あなたが好きな答えを入力するために出力がどのように変化するかが十分に分かります。

[ 10000.] 
[ 10000.00000001] 
[ 10000.] 

それだ：ここで

は、コントラストのために、既定のアルゴリズムによって行われた全体の検索です。わずか1〜8秒ずつステップアップを試み、それが得られる答えを見分けることができず、他の答えのいくつかが正しくアサーションされているので、あきらめます。

時には、（a）より大きいステップを開始し、（b）ステップサイズが小さくなるとテストを中止するようにアルゴリズムに指示することで解決できます。ミリ秒。まだデフォルトの最小化技術は建設この助けを与えても、最小を見つけることが、私たちは何かを行うことができますあまりにも脆弱であると思われる。この場合

out_s_def = minimize(separation, sec_init, args=(32.5, 215.1), 
        tol=1e-3, options={'eps': 500}) 

：：私たちは最小化を伝えることができるようにあなたが何かを試してみてくださいそれが実際にどれくらい多くのビットを再生するかを機能させます。

これらの最小化ルーチンは、64ビット浮動小数点の精度が利用可能になるまでにどのくらい正確にプッシュでき、そのポイントの前で停止するように設計されているかについて、かなり明示的に書かれていることがよくあります。しかし、あなたは精度を隠しています。あなたはルーチンに「秒数をつけてください」と言っています。実際には秒が組み合わさっている秒値の非常に小さいローエンドの数字でさえも、時間と日だけではなく数年で、このプロセスの中で、秒の最下位にある小さな精度は失われます。最小化はわかりません！

実際の浮動小数点時間をアルゴリズムに公開しましょう。

1. は、あなたがやっているfloat()操縦の必要性を回避してみましょう：プロセスでは、私は3つのことをやります。私たちのprintステートメントは、問題を示しています。あなたは、スカラ浮動小数点数を提供していても、最小化はとにかくnumpyの配列に変換しますこと：

   (array([ 10000.]), 32.5, 215.1) 
   

   しかし、それは修正するのは簡単です：今Skyfieldはそれができるに建てseparation_from()を持っていることうまく配列を扱う、我々はそれを使用します。

   sepa = mpos.separation_from(spos) 
   return sepa.degrees 
   

2. 私はそれが1.0に向けてヘッドとしてSkyfieldが採用していること、日付を作成するための新しい構文に切り替わります。

私たちのようなものを与える

（しかし、あなたは一度だけtoposを構築し、それが渡された場合、これは代わりにそれを再構築し、それがすべての時間をかけてその数学をやらせで、速くなることに注意してください）：

ts = load.timescale() 

... 

def separation(tt_jd, lat, lon): 
    place = earth.topos(lat, lon) 
    t = ts.tt(jd=tt_jd) 
    mpos = place.at(t).observe(moon).apparent() 
    spos = place.at(t).observe(sun).apparent() 
    return mpos.separation_from(spos).degrees 

... 

sec_init = 10000.0 
jd_init = ts.utc(2016, 3, 9, 0, 0, sec_init).tt 
out_s_def = minimize(separation, jd_init, args=(32.5, 215.1)) 

結果は、成功した分娩であり、私がここで私を二重にチェックできるかどうかは分かりますか？私はSkyfieldで事前に構築された縮小ルーチンの数をバンドルする間もなく期待し

print ts.tt(jd=out_s_def.x).utc_jpl() 

['A.D. 2016-Mar-09 03:37:33.8224 UT'] 

- 実際には、PyEphemを交換するためにそれを書くための大きな理由のことができるようにしたいと思った： - あなたが探している答え強力なSciPyオプティマイザを発揮し、PyEhemがCで実装しているような貧弱なものを放棄できるようにする必要があります。ここでは、オプティマイザに浮動小数点数を与える必要があります。途中。

多分私は、Timeオブジェクトが2つの浮動小数点オブジェクトから時間を構成できるようにする必要があります。そうすれば、より多くの秒数を表現できます。私はAstroPyがこれをやったと思うし、天文学のプログラミングで伝統的です。