MathematicaのNMinimize, NMinValueのバグ

Ver.11で修正されました。

Mathematica 9.0, 10.0, 10.1, 10.2, 10.3, 10.4.1 for Microsoft Windows (64-bit)と10.0.0 for Linux ARM (32-bit)でのことです。

a2 + c2 == 1, b2 + d2 == 1という条件の下で、次のような関数の最小値を求めたいとします。

Mathematicaには最大・最小を求めるためのさまざまな関数が用意されていますが、NMinimizeMinimizeNMinValueMinValueは、頭にNをつけるかどうかで、数値的な方法と解析的な方法を切り替えられて便利です。最小値とその時の各変数の値を知りたいときは[N]Minimize、最小値だけを知りたいときは[N]MinValueを使います。

しかし、ちょっとうまくいかない例に遭遇しました。バグが修正されたバージョンを受け取るために、プレミアユーザになっています。

問題その1

f1 = Sqrt[
 1 + 4 Sqrt[3] b c d + d^2 - 5 c^2 (-1 + d^2) + 
  2 a (-Sqrt[2] b d + Sqrt[6] c (-1 + d^2))];

NMinimize[{f1, a^2 + c^2 == 1 && b^2 + d^2 == 1}, {a, b, c, d}]

このように実行すると、

関数の値4.35271 +1.45672\ Iは{a,b,c,d} = {0.79784,1.52261,-0.444508,0.634251}において実数ではありません.

という警告とともに結果{1.22618, {a -> -0.82781, b -> 0.492467, c -> 0.224354, d -> -0.333602}}が返ります。このとき、a2 + c2は約0.74、b2 + d2は約0.35なので、与えた制約条件から大きくずれています。

問題その2

f2の最小値を求めようとしたときにも問題が起こります。

f2 = Sqrt[
  1 - 4 Sqrt[3] b c d + d^2 - 5 c^2 (-1 + d^2) + 
   2 a (-Sqrt[2] b d - Sqrt[6] c (-1 + d^2))];

NMinimize[{f2, a^2 + c^2 == 1 && b^2 + d^2 == 1}, {a, b, c, d}]

このように実行しても、入力がそのまま返ってくるだけで、まとも(?)な結果が得られません(間違った答えが返ってくるよりはましかもしれません)。

この問題は、WolframAlphaでも起こります(2015年5月17日確認)。

f1に関しては間違った結果が返ってきます(下はNMinValueの場合。NMinimizeでも同様)。

f2に関しては、珍しいことに、WolframAlphaは完全に沈黙します(下はNMinValueの場合。NMinimizeは違う結果になります)。

MathematicaからWolframAlphaに問い合わせたときはまた違う結果で、f2は解釈できないとのことでした。ほとんど同じf1は解釈できたので、そんなことはないと思うのですが。

そもそも、こういう関数の最小値を何の工夫もせずに得ようとするのが間違いなのかもしれませんが、あえて工夫をしないのにも理由があるのです。それについては別の機会に書きます。

病的な計算式?(『ソフト・エッジ』から)

4621053833中島 震・みわ よしこ『ソフト・エッジ: ソフトウェア開発の科学を求めて』(丸善出版, 2013)で、数値計算の面白い例が紹介されていました。(参考文献リストあり。索引なし。)

精度を向上させれば近似はよくなるというのが、自然な考え方でしょう。でも、この直観に当てはまらない「病的な計算式」の例があります。(p.132)

ということで、f(a, b)=(333.75-a**2)*b**6+a**2*(11*a**2*b**2-121*b**4-2)+5.5*b**8+a/(2*b)としたときの、f(77617, 33096)が紹介されています。単精度・倍精度・4倍精度の浮動小数点数での計算結果は約1.17で、正しい値(約-0.827)とは符号も違ってしまうとのことです。

これだけ読むと、精度を向上させても近似はよくならないと思うかもしれませんが、そんなことはないはずです。実際、精度をある程度高くすれば、正しく計算できます。

Pythonで試します(ideone.comでも動きます)。

from decimal import getcontext, Decimal
for n in range(16, 40):
  getcontext().prec = n
  a = Decimal("77617")
  b = Decimal("33096")
  print(n, (Decimal("333.75")-a**2)*b**6+a**2*(11*a**2*b**2-121*b**4-2)+Decimal("5.5")*b**8+a/(2*b))
16 1.000000000000000E+21
17 -4.0000000000000000E+20
18 -2.00000000000000000E+19
19 2000000000000000001
20 -99999999999999998.827
21 1.17260394005317863186
22 -999999999999998.8273961
23 100000000000001.17260394
24 10000000000001.1726039401
25 -3999999999998.827396059947
26 -99999999998.827396059946821
27 -19999999998.8273960599468214
28 -999999998.8273960599468213681
29 100000001.17260394005317863186
30 20000001.1726039400531786318588
31 -999998.8273960599468213681411651
32 300001.17260394005317863185883490
33 -9998.82739605994682136814116509548
34 -1998.827396059946821368141165095480
35 -198.82739605994682136814116509547982
36 21.1726039400531786318588349045201837
37 -0.827396059946821368141165095479816292
38 -0.8273960599468213681411650954798162920
39 -0.82739605994682136814116509547981629200

このようにPythonでは、getcontext().precの値が37以上なら適切な結果が得られます。

Rubyでも、ちょっと工夫すれば適切な結果が得られます(ideone.comで確認)

Windows付属の電卓では、何の工夫も要りません(WIndows 7で確認)。xyのショートカットキーが「Y」なので、「(333.75-77617Y2)*33096Y6+77617Y2*(11*77617Y2*33096Y2-121*33096Y4-2)+5.5*33096Y8+77617/(2*33096)=」と打てば、適切な結果が得られます。

Mac OS Xでも計算できます(10.8.4で確認。できるようになったのは比較的最近のことだと思いますが)。「計算機」で計算する方法はよくわかりませんが、Spotlightに「(333.75-77617^2)*33096^6+77617^2*(11*77617^2*33096^2-121*33096^4-2)+5.5*33096^8+77617/(2*33096)」と入力すれば、適切な結果が得られます。

Windowsの電卓は32桁までの入力にしか対応していないので、230928922463004537535516678003369の平方根を求めたいというようなときに不便ですが、Spotlightなら、sqrt(230928922463004537535516678003369)とすればとりあえず計算できます。

参考:電卓に求められるコト

本稿執筆時点のWolframAlphaでは、まったく違う結果が2つ返ってきて、片方はあっていました

Googleはだめでした

Googleの例は別として、精度を上げてもうまくいかないようにするには、もっと「病的」でなければなりません。

米国製の計算機と日本製の計算機の違い

AppleとSonyの一番の違い part IIという記事を見て、「なるほど」と思ったのですが、こういう「たとえ」を使わなくても、そのままずばりのわかりやすい例がありました。

「高さが3の正三角形の面積は?」と訊かれたら、たいていの人は答えられると思いますが、ウェブで検索するより早く答えられる人は少ないでしょう。実際に、ウェブで解決しようとするとき、米国製の計算機と日本製の計算機、どっちを使いますか?

米国製の計算機

日本製の計算機

米国製の計算機、WolframAlphaの場合

「equilateral triangle height=3」と入力すれば答えが得られます。簡単に答えがわかるのはもちろん、結果のURLが得られるのもいいですね。

日本製の計算機、高精度計算サイトkeisanの場合

  1. 「数学・物理の計算」をクリックする
  2. 「三角形」をクリックする
  3. 「正三角形」をクリックする
  4. 入力指定を「高さ」にする
  5. 高さを「3」にする
  6. 「計算」ボタンをクリックする

これでやっと答えが得られます。初めて行く人は、どこに何があるかがわからないので、「手で計算した方が早い」ということになるでしょう。

keisanがWolframAlphaになる必要はありませんが、「英語がわからなくても大丈夫」以外のメリットを打ち出せるようになるといいですね。

スリーマイル島の裏側はチャイナではないしフクシマの裏側もブラジルではない (WolframAlpha)

まあ、常識なんでしょうが、みんな大好きWolframAlphaで確認しておきましょう。

「location “three mile island”」で検索
スリーマイル島は40.15°N, 76.72°W
「40.15°N, 76.72°W antipode」で問い合わせ

チャイナからはかなり遠いですね。

スリーマイル島よりレベルが高いせいか、フクシマは「antipode fukushima」という1回の問い合わせでできました。ブラジルというよりは、ブラジルの南の海ですね。



Google Mapsで調べる方法:ブラジル・シンドローム

長さを変えた15個の振り子を一斉に揺らす実験

長さを変えた15個の振り子を一斉に揺らす実験(動画)をMathematicaでシミュレートする。

2次元

3次元

手順は以下のとおり。

  1. 運動方程式の確認(g=1、m=1、糸はたるまないことを仮定)
  2. 1分間の振動数が整数になるような振り子の長さを求める
  3. 運動方程式を解いて結果をアニメーションにする

ラグランジの運動方程式は次のように確認できる。

lagrangian = (r \[Theta]'[t] )^2/2 - r Cos[\[Theta][t]];
D[D[lagrangian, \[Theta]'[t]], t] + D[lagrangian, \[Theta][t]]

r Sin[theta[t]] + r^2 theta''[t]

これを使って、シミュレーションを作る。

(* 1分間の振動数が整数になる長さを求める *)
tMax = 60; (* maximum time *)
n = 15; (* number of objects *)
theta = Table[Unique[], {n}];
initialAngle = Pi/6;
r = Table[
   First[x /. Solve[
      tMax/(4 Sqrt[x] EllipticK[Sin[initialAngle/2]^2]) == (20 + i), x]
    ], {i, 1, n}];

(* 運動方程式を解く *)
equations = Flatten[Table[{
     (* ラグランジの運動方程式 *)
     r[[i]] Sin[theta[[i]][t]] + r[[i]]^2 theta[[i]]''[t] == 0,
     theta[[i]][0] == Pi/6,
     theta[[i]]'[0] == 0},
    {i, 1, n}]];
solution = NDSolve[equations, theta, {t, 0, 1.1 tMax}];

(* 2次元アニメーション *)
Animate[
 points = Table[First[
    {r[[i]]  Sin[theta[[i]][t]], -r[[i]]  Cos[theta[[i]][t]]}
      /. solution /. t -> T], {i, 1, n}];
 Graphics[{
     Black, Disk[#, 0.005],
     Gray, Line[{{0, 0}, #}]} &
   /@ points,
  PlotRange -> {{-0.8 Max[r], 0.8 Max[r]}, {0, -1.1 Max[r]}}],
 {T, 0, 1.1 tMax}, DefaultDuration -> 1.1 tMax, 
 SaveDefinitions -> True]

実行結果は冒頭の動画の通り。

最後の部分を次のようにすれば3次元になる。

(* 3次元アニメーション *)
Animate[
 points = Table[First[
    {i Max@r/n, r[[i]] Sin[theta[[i]][t]], -r[[i]] Cos[theta[[i]][t]]}
      /. solution /. t -> T], {i, 1, n}];
 Graphics3D[{
     Black, Sphere[#, 0.005],
     Gray, Line[{{#[[1]], 0, 0}, #}]} &
   /@ points,
  PlotRange -> {{0, 1.1 Max@r}, {-0.7 Max@r, 
     0.7 Max@r}, {0, -1.1 Max@r}}],
 {T, 0, 1.1 tMax}, DefaultDuration -> 1.1 tMax, 
 SaveDefinitions -> True]

参考