科学と家事とプログラミング (python を中心に)

python 温度計測 湿度計測 DS18B20 USB9097

多変数最適化の Nelder-Mead は、なかなか使い勝手が良いのです

多変数最適化のアルゴリズムの一種である Nelder-Mead は、なかなか使い勝手の良いアルゴリズムですね。 汎用的に使える(対象が最小二乗に限定されない)ことに加えて、 微分情報を使わないところも魅力です。 評価関数の性格が多少悪くても、上手に最適をみつけてくれる、便利屋さんといったところでしょうか。

Nelder-Mead のお勧めバージョン

Nelder-Mead には、いろんな改良/拡張が存在しますが、 場合分けの見た目の安心感? から、以下の 5 種類の動作を含んだ設計がお勧めです。

  1. reflection
  2. expansion
  3. inside-contraction
  4. outside-contraction
  5. shrinking

細かく言うと、最悪, 2番悪, 最良 の 3点の情報を使って 7 通りに場合分けしたあと、 動作が同じになるところをまとめると、最終的に 5 ケースとなります。 例えば http://adl.stanford.edu/aa222/lecture_notes_files/chapter6_gradfree.pdf

アルゴリズムの妥当性と試験について

この手のアルゴリズムは、経験則の寄せ集め的なところがあって、 絶対的な正解という考え方には馴染まない気がします。 適用限界を明示するのが難しい。 あと設計やコーディングに間違いを含んでいても、 なんとなく動いてしまって、しかも、本来と異なる経路で最適値に収束したりして、 外部仕様を定めにくいですね。つまり、ブラックボックス的な試験設計は難しい。 開発手法(管理)の観点からは、興味深い対象かもしれません(扱いが難しいということ)。 というわけで、ホワイトボックス的な試験に重点を置くのが良さそうです。とにかく全てのパスを通す。 このとき、見やすい図を作って直感的に動作を確認するのは、 そう悪くない方法かと思ったりもします。例を下図に示します。

f:id:sken20k:20180513141435j:plain

適当な評価関数(2次元)を用意して、乱数で初期値を選ぶようにして、 見やすい図が得られるまで乱数のシードをいろいろ試しました。。 スライダーなどの GUI で初期値を選べるようにすると、おもしろいデモが 作れそうですね。

パラメータの一部を固定して評価したいとき

最適化の問題をあつかっていると、一部のパラメータを固定して推定したい。 という場面に出くわすことがあります。 パラメータの組み合わせ毎に最適化対象の関数を定義しはじめると、 評価関数だらけになってしまって、何をやってるのか分からなくなりがちです。 最適化のルーチンの I/F に手を加えて、要素毎に推定する/しないを指定するフラグを追加 しても良いのですが、 python のようなインタープリターを使うのであれば、 固定パラメータを補う関数を動的に生成するのも悪くないですね。 こうすれば、ベーシックな最適化ルーチンに手を加えずシンプルに保ったまま、 目的を達することができます。

落穂ひろい

Nelder-Mead 以外で、実用の見地から押さえておきたい最適化アルゴリズムは、 やっぱり Levenberg-Marquardt でしょう(LM法)。 適用できる問題は最小二乗形式に限定されますが、 変数の数が多めのデリケートな問題でも、良好な結果が得られることが多いですね。 推定誤差や correlation を評価がしやすいところもポイントが高いです。 注意が必要なのは、ダンピングファクターの初期値の調整くらいでしょうか。 実装にもよりますが、ダンピングファクターと収束判定の関係を理解していないと、 使いこなせない場面が、時々あります。 最適化ルーチンをブラックボックス的に利用するのは良いのですが、こういう 特性を理解した上でないと、使いこなせないこともある。 ブラックボックスの表面に、そういう細かな特性を認識できているかどうかで話が大分違ってきそうです。 個々のアルゴリズムや理論の前提条件や適用限界に注意を促したり、それらを身に付けるための方法論が 求められる。そう思う次第です。

剛体回転に関するオイラーの定理は、絵的に理解したいのです

「剛体回転におけるオイラーの定理」について、 天球図を使った直感的な説明を紹介します。 代数的な説明だけだと、なんとなく分かった気がしない...という人におすすめです。

オイラーの定理

オイラーの名の付く定理や公式は、いくつもあって紛らわしいのですが、 その中の一つに「剛体回転におけるオイラーの定理」というのがあります。

【定理】 「点のまわりの回転変位は, その点を通る1つの軸のまわりの回転によって達せられる」

3次元座標の座標回転の言葉で言い直すと、

「原点まわりの任意の座標回転は、ある方向ベクトルを回転の軸とする1回の回転で達成できる」

1回の回転で変換を達成する回転の軸を「オイラー軸」あるいは「オイラー回転軸」と呼び、 その時の回転の角度を「オイラー回転角」と呼んだりします。

オイラー角」と名前が似てますが別物です。「オイラー角」が座標軸回りの3回の回転の組み合わせであるのに対して、「オイラー回転角」は任意方向回りの1回の回転です。

幾何学的な説明

f:id:sken20k:20180427211221j:plain

上の図は、座標系 X-Y-Z から X'-Y'-Z' への回転を表しています。 ひとまず Z を除いて、X と Y に注目した上で、回転軸の成分が保存される (回転軸とXの成す角は、回転軸とX'の成す角と等しい)ことに注目すれば、 回転軸の方向は、 X, X' を底辺とする二等辺(球面)三角形の頂点に存在するはずです。つまり、 X を X' に移す回転軸は X と X' の垂直二等分線のどこかに存在する(球面なので大円のどこか)。 同様に、Y を Y' に移す回転軸は Y と Y' の垂直二等分線のどこかに存在する(球面なので大円のどこか)。 これらを同時に満たす点は2つの大円の交点として求まるので、これを P とする。 (逆方向の解も存在するが、どちらを選んでも議論の一般性は失われないので見やすい方を選ぶ)。

これで、オイラー回転軸の候補 P が手に入りましたが、 X', Y' を同時に満たす回転(角度)が存在するかどうか分かりません。 以下、それを確認します。

XをX'に移す回転角をφx、YをY'に移す回転角をφyとする。 オイラーの定理が成り立つためには、φx とφy が一致すればよい。

Pの成分が座標回転に対して不変なので、 球面三角 P-X-Y と P-X'-Y' は合同であることが分かる。 Pを頂点とする角度から重なり部分を除いた角度がφxとφyなので、 φxとφyは確かに等しい。

Z, Z' についても同様の議論が成り立つ。

これで、オイラーの定例が成り立つことが確認できました。

代数的な説明

先に幾何学的な説明をしましたが、 オイラーの定理は代数的に説明されることの方が多いですね。 座標回転を表す行列の固有値固有ベクトルを使った議論は、 あいまいさなく簡潔に証明できるところが利点と思います。 こちらは、どの教科書にも載っているので説明は省略します。

おまけ

オイラーの定理の絵的な理解は、実践的(応用範囲が広い)知見と思います。 図解した教科書が少ない(気がするので)ここに紹介しました。

補間と言えば akima-spline は知らなきゃもぐりだろう、と言ってみたいのです

スプライン補間でもっともメジャーなのは cubic-spline でしょうか。 ある種の理解しやすさや微分も一緒に求まるという使い良さもあるのですが、 データにギャップがある場合に、補間結果が振動するという やっかいな特性を持っています。 下の図の赤線は、サンプル点(黒点)を cubic-spline で補間した結果です。 ギャップの前後で補間値があばれているのが分かります。

同じ図に、akima-spline で補間した結果を青線で示します。 ギャップの前後も補間値があばれることなく、スムースな補間結果が得られことが分かります。

f:id:sken20k:20180420234011j:plain

参考までに、別の例も載せておきます。 f:id:sken20k:20180420235709j:plain

補間方法には、それぞれ利点/欠点はあるので、 akima-spline が万能というわけではありませんが、 区間ごとに推定した軌道や姿勢をつなぐ場面では、つなぎ目に微小なギャップが 現れることが多いので、そういうとき akima-spline のなめらかさ(オーバーシュートの無さ)は とてもうれしい特性なのでした。

akima は、秋間浩さんという日本人の名前に由来します。

高級な補間公式が、万能とは限らないのです。

補間公式には、ラグランジュ、エルミート、スプライン...と、各種ありますが、 私がよく使うのは、線形補間(直線補間)とAkima-Spline の 2 種類くらいです(JPL の惑星歴のチェビシェフは除く)。 線形補間(直線補間)とは、随分と原始的なものを使うのですね...と思う人がいるかもしれませんが、 高級な補間公式で安全性を保障するのは、結構むずかしかったりします。

高級な補間公式の弱点

高級な補間公式は、(使い方を間違わなければ) 少ないサンプル点で高い精度を実現できますが、 それらは以下のような欠点ももっています。

  • 単調性が保証されない
  • 高階差分を適用するとき桁落ちが発生しやすい

これらの欠点について注意が払われることは、少ないように思います。 数値計算の教科書も、精度の観点から議論が多くを占めて、 単調性が保証されない点に関しては、たとえ記述されていたとしても さらりと一言述べられるだけの事が多いようです。 記述量が少ないため、読み飛ばされることもしばしば。 実運用や解析ツールの設計においては、とても大切な特性ですが、 痛い目に会わないと、なかなか身に付かない知見の一つかと思い、 ここに紹介します。

単調性の保証が大事な時

例えば、軌道情報が時系列データ(点列)として与えられたときに、 距離 R がある閾値を下回るタイミングを求めたいとする。 このとき、単調性が保証されない補間公式を使うと、 偽のイベントがジッターのように出現する可能性がある。 次の図では、本当のタイミングは time が 1.1 を超えあたりにありそうだけど、 高次の補間を使って求めようとしたために 0.9 ~ 1.0 の間に偽の事象が出現しています。

f:id:sken20k:20180417185817j:plain

補間の時に最後の Time = 1.2 の点を使わないのが原因ではないか? と思う人がいるかもしれませんが、Time = 1.0 までの点が揃ったところで、 直近のN 点から事象を判定する場面は、そう無理な仮定ではなかろうと思います。 高次の補間公式で単調性を保証するのは、いろいろとむずかしいので、 とにかく安全を保証したい場面でであれば、単調性が保証できる線形補間を使うことを お勧めする次第です。

上の図は、単調性が保たれない箇所を拡大しています。元のスケールで補間の様子全体を 表示すると以下のようになります(6点を使ったラグランジュ補間)。 f:id:sken20k:20180420183939j:plain

高階差分の桁落ち

高次の補間は、単調性の話以外にも精度の面で問題が生じる場合があります。 具体的には、補間公式には高階の差分演算が現れるため、似たような数値の減算の繰り返し が多くなり、次数が高くなるほど、つまり階数が深くなるほど、桁落ちが生じやすくなるからです。 これに対抗するには、最初のデータの有効桁を増やすしかありません。 時刻を例にとれば、準ユリウス日ユリウス日を時刻に使っていたりすると、 高階差分のために思わぬ誤差が生じることがあるので注意が必要です。 最近では、2000年からの時刻が使われたりすることもありますが、高い精 度が求められる場合には、解析期間の開始時刻を時刻原点に採用するのが無難です。 普段、あまり次数の高い補間公式を使わないようにしていても、DOP853 という 高精度の数値積分のルーチンに組み込まれた densout (稠密出力)などを利用する 場面で、あまり意識することなく、高階の差分を使ってしまうこともあります。

落穂ひろい

高次の補間公式が振動する話は有名ですが、単調性に焦点を当てた議論というのは 少ないように思います。教科書を書いている数値計算の専門家にとっては あたりまえのことであるので、ページ数が限られている教科書においては、 話が精度に偏ってしまうのも仕方ないことかもしれません。あるいは、 数値計算の専門家(研究者)には、実運用のシステムを設計/運用する経験はほとんどないと 思われるので、適用場面で重用になる特性に興味が薄いことも一因かと疑ったりしています。 というわけで、たとえあたりまえの話であっても、ある程度のボリュームを割いて、議論する意味 というのはあるのかなと思う次第です。もちろん、知識の効率的な吸収という観点とは相反するわけだけれども、 ある種の無駄は必要だろうという立場をとります。

少し一般的な話をすると、およそ公式や手法には、それぞれ利点だけでなく欠点もあるものです。 使える場面や条件を明晰に意識できるかどうかが、応用が利くかどうかの分かれ目ですね。 利点だけでなく欠点や制約を2つ3つ列挙できるかどうかが、その手法の理解度を測るバロメータ だと言う人もいます。

Akima Spline については、別記事にて...

緯度経度で分割された星カタログの選択は、なかなかやっかいなのです

カメラ視野を模擬したい

カメラの視野を模擬するとき、 手元の星カタログを全部スキャンすると星の数が多すぎて性能が出ない時があります。 こういう時は、星カタログをいくつかの部分に分けて対応します。 分けたものをサブカタログと呼んだりします。

経度緯度で分割したサブカタログ

サブカタログの分割方法で一番簡単なのは、緯度と経度で分ける方法でしょうか。 さて、分割したからには必要なサブカタログがどれかを調べなければいけません。 簡単なのは、各サブカタログの4隅(極だと3隅)と中心方向が、カメラの視野に入るかどうかを判定することでしょうか。 こういう方法だと、辺の一部に視野が重なる場合などはうまく判定できません。 特に困るのが極付近の判定です。その難しさは、下の図を見てもらえば分かってもらえるでしょうか。

力技による判定は、細かく面でサンプルする方法ですが、判定に時間がかかって残念な感じになってしまいます。 なんとかスマートに判定したいものです。

ビットマップを使ったサブカタログの選択

ここで、ビットマップを使った視野判定の方式を紹介します。

f:id:sken20k:20180330172056j:plain

手順はこうです

  1. 視野より大きめの描画エリアを用意して、カタログの境界線を描画する
  2. 背景色で視野を抜く
  3. 中心付近を視点として、塗りつぶす。
  4. 各サブカタログの中心(代表点)の色を調べる

細かな説明は不要と思いますが、これで視野の端がサブカタログにちょっとだけかかるときも、 ビットマップの分解能で正しく判定できます。実用上は、視野を少し大きく採ればよいでしょう。 けっきょく、面でサンプルしているといえなくもないのですが、サブカタログの境界線と境界線の間は、 ビットマップの塗りつぶし処理だけで代用できているので、ちょっとお得かなと思います。

他の方式

球座標のような特異点を避けるという観点からは、N 面体で分割するのも一案です。 用途によっては、排他的に分割するのではなく、重複を許すと便利そうです。 他にも分割の工夫はありそうですが、最近は使えるメモリーも桁違いに増えたので、 一括で扱った方が良い場面も多いかもしれません。 ビットマップを使った解法というのは、他にも応用が利きそうに思うので、一応紹介しておきます。

tanθに従う現象をべき級数で近似するのは、なんだか残念なのです

tan のべき級数展開には気をつけよう

べき級数で近似するのは手軽な方法ですが、対象とする現象が tan に従う場合には、 次数を上げたところで精度の向上は期待できません。

かりに次数を上げることで必要精度を達成できた場合でも、 それは単に tan を再現するために高次の項を必用とするだけかもしれません(tan の収束半径は有限)。 たとえば、θ が 60 deg (1.25 rad) くらいになると、次数を 9 まで増やしても 目に見えるスケールの誤差が残ります。 なにかしらの理想曲線からのズレや歪の表現としてべき級数を選んだのだとしたら、 ちょっと話が違うのではと思う次第です。

f:id:sken20k:20180324175558j:plain

実際の場面で思いつくのは、光学系の歪補正の補正式などでしょうか。 θの範囲はそれほど広くなくても、 要求精度が高くなると上のような状況が発生します。 理想的なピンホールカメラ(平面への投影)は、まさに tan(θ) ですから、 そこからの歪を表現するなら、θの級数ではなく、w = tan(θ) として w のべき級数を採用すると、少ない次数でよく近似できる場合があります。

θのべき級数が使われてしまう理由も、なんとなくは想像はできて、 教科書で扱われる演習問題などでは、θの一次式や二次式で定式化したものが多く、 精度や範囲を限定したものになっている。現場で、要求精度が厳しくなったりθの範囲の広く なっても、根本的に見直すことなく(改良と称して)、θの次数を上げて精度を確保しようとする。 そんなことが起こっている気がします。 ただし、組み込み分野などで tan などの超越関数を気軽に使えるようになったのは、 わりと最近のことなので tan が使えない環境であればべき級数による近似は、合理的 なのかもしれません。

sin, cos の場合

sin や cos の場合は、収束半径が無限大なので、次数を上げると精度も上がるし、 近似できる範囲もどんどん広がります。 次数を上げることで再現できる波の数が増えていくというのが、直感的には 理解できなかったりします ^^;

f:id:sken20k:20180324175629j:plain

座標変換ルーチンの爆発的な増殖は、なんとしても避けたいのです

問題設定

天文宇宙に関連したアプリケーションでは、多くの座標系が登場します。 黄道座標系、地球赤道座標系、地球固定座標系、月固定座標系、などなど。 さらに、赤道座標系や黄道座標系には、分点の違いによって、 基準エポック、平均分点(MOD)、真分点(TOD) といった区別もあります。

いくつも座標系がある状況で、相互の座標変換ルーチンを場当たり的に用意しはじめると、 座標変換ルーチンばかりが増えていきます。いわゆる、組み合わせ論的に爆発ですね。 こういう事態は、なんとしても避けねばなりません。

解決案

こういうときは、以下の指針が有効です。

  • 基準座標を定める
  • 各座標系と基準座標系との変換だけ用意する
  • 相互の変換は必ず基準座標経由で実施する

以下、オブジェクト指向的な実装の例を示します。 (座標系 R を基準座標として、座標系 A, B, を実装する)

class Frame:
    def __init__(self, t):
        """ 生成時に時刻を指定して座標変換行列を固定する """
        self.sys_fm_ref = None # sys_fm_ref は、基準座標からの座標変換行列
        
    def _from(self, other):                  # self が A, other が B なら
         """ 座標系 other から self への座標変換行列 """
        sys_fm_ref = self.sys_fm_ref         # A <== R
        ref_fm_oth = other.sys_fm_ref.tran() # R <== B  (B<==R を転置)
        sys_fm_oth = sys_fm_ref*ref_fm_oth
        return sys_fm_oth
        
class R:
    def __init__(self):
        # 基準座標から基準座標への変換は、単位行列
        self.sys_fm_ref = mat([[1, 0, 0], [0, 1, 0], [0, 0, 1]])

class A:
    def __init__(self, t):
        """ 基準座標からの座標変換行列 """
        self.sys_fm_ref = rotx(1.0*t) # 適当

class B:
    def __init__(self, t):
        """ 基準座標からの座標変換行列 """
        self.sys_fm_ref = roty(1.0*t) # 適当

def test():
    t = 0.0
    a_fm_b = A(t)._from(B(t)) # A <== B  座標回転行列

新たに座標系 C を追加したいときは、class C を準備するだけで、 すでに存在する全ての座標系との変換行列が手に入ります。

c_fm_a = C(t)._from(A(t))
c_fm_b = C(t)._from(B(t))
c_fm_r = C(t)._from(R())
a_fm_c = A(t)._from(C(t)) 
b_fm_c = B(t)._from(C(t)) 
r_fm_c = R()._from(C(t)) 

利点

これで、座標変換ルーチンの組み合わせ論的な爆発を防ぐことができます。 検証の面でも有利でしょう。

欠点

倍精度の限界のような高精度が求められる場合には、こういう方式では対応できない 場合があります。性能(実行速度)面は、まったく考慮していません。

捕捉

基準座標の R は、例えば J2000.0 EQ や ICRF などを割り当てます。 座標系の数が増えれば増えるほど、この設計の有難味が増します。 時刻系の変換ルーチンについても、同じような方針を適用できます。

f:id:sken20k:20180324202958j:plain

コード中、fm は from の略です。 座標回転の結合は、a_fm_c = a_fm_b * b_fm_c という感じ記述できます。何がよいかというと、 * の前後に現れる座標系が同じである確認した上で(斜線で)消去する操作を、プログラムコードで直接表現できます。 座標系の変換でよくある変換の向きの間違いを、ソースコードの形式的な工夫によって減らすことを狙ってます。