Mathematica 10への期待(とその後)


昨年ついに無料になったMathematicaですが、思ったほど盛り上がってはいないようです。盛る上がらない原因の一つに、性能の悪さがあるのは間違いないでしょう。特にGUIはかなり重く、Manipulateなどでインタラクティブに遊ぶのは無理です。(無料のMathematicaが動くのは、700MHzのCPU、最大でも512MBしか主記憶のないRaspberry Pi上。)

実際にベンチマークを取ってみると、Windows上のMathematica 9.01(CPUはCore i7 4700MQ、主記憶は16GB)で10秒で終わるMathematicaMark9が、Raspberry Pi上では約3300秒かかりました(800MHzにオーバークロックして約3000秒)。Core i7で動かすエミュレータ上でもおそらく同じ程度でしょう。

しかし不思議なことに、Raspberry Pi上のMathematicaのほうが早く解ける問題があります。

たとえば、A^2 + B^2 == 1かつP^2 + Q^2 == 1という条件で、以下のsの最小値を単純にMinValueで「解析的に」求めようとすると、上述のWindowsマシンではなんと約20000秒かかるのに対して、Raspberry Piでは約2200秒で終わります。追記:Raspberry Pi 2では約800秒です。

s = 1/24 (3 Sqrt[3] Abs[A P] + Sqrt[3] Abs[A P + 2 Sqrt[2] Q] + 
     Abs[Sqrt[3] A P - 3 Sqrt[2] B P - Sqrt[6] Q] + 
     Abs[Sqrt[3] A P + 3 Sqrt[2] B P - Sqrt[6] Q]);

MinValue[{s, And[A^2 + B^2 == 1, P^2 + Q^2 == 1]}, {A, B, P, Q}]

(これは比較のために作った問題ではなく、実際に私が遭遇した問題です。実は、単純にMinValueを使うというのが間違っているのですが、本来Mathematicaはそういうことを気にせず使うソフトウェアのはずです。)

Raspberry Piの方が速い原因は、Windows版Mathematicaのバージョンが9.01なのに対して、Raspberry Pi版Mathematicaのバージョンが10であり、バージョン9.01と10の間で、このあたりのアルゴリズムが改良されたためだと思われます。

非力なRaspberry Piでの話なので、比較的速いCPUがあれば、エミュレートされたRaspberry Piでも同程度の時間でできるでしょうし(Windowsの場合Ubuntuの場合)、ユーザーモードエミュレーションならこの半分くらいの時間でできるでしょう(ライセンス違反なので想像しているだけです)。

というわけで、Windows版Mathematica 10のリリースを待っています。

追記:Windows版Mathematica 10がリリースされたので試したところ54秒でした(CPUはCore i7 4700MQ、主記憶は16GB)。Raspberry Piは約2200秒なので、遅すぎです

フェルマーの最終定理の「反例」(MySQL編)


MySQLでは、a=5999856E0, b=99992800E0, c=100000000E0としたときに、a^3 + b^3 = c^3になります。

SELECT
  5999856E0*5999856E0*5999856E0
 +99992800E0*99992800E0*99992800E0
 -100000000E0*100000000E0*100000000E0\G
*************************** 1. row ***************************
5999856E0*5999856E0*5999856E0
 +99992800E0*99992800E0*99992800E0
 -100000000E0*100000000E0*100000000E0: 0
1 row in set (0.00 sec)

多くのプログラミング言語と違うのは、5999856.0のように書くと結果が変わることです。

SELECT
  5999856.0*5999856.0*5999856.0
 +99992800.0*99992800.0*99992800.0
 -100000000.0*100000000.0*100000000.0\G
*************************** 1. row ***************************
5999856.0*5999856.0*5999856.0
 +99992800.0*99992800.0*99992800.0
 -100000000.0*100000000.0*100000000.0: -2985984.000
1 row in set (0.00 sec)

5999856と5999856E0、5999856.0は、

  • 算数や数学では、5999856=5999856E0=5999856.0(すべて同じ)、
  • 自然科学では、5999856=5999856E0≠5999856.0(有効数字の意味で)、
  • 多くのプログラミング言語では、5999856≠5999856E0=5999856.0(整数と浮動小数点数)。

多くのプログラミング言語と異なり、MySQLではすべて別物です。次のように確認できます。

CREATE TEMPORARY TABLE t
  SELECT 5999856 AS a,5999856E0 AS b,5999856.0 AS c;

DESCRIBE t;

+-------+--------------+------+-----+---------+-------+
| Field | Type         | Null | Key | Default | Extra |
+-------+--------------+------+-----+---------+-------+
| a     | int(7)       | NO   |     | 0       |       |
| b     | double       | NO   |     | 0       |       |
| c     | decimal(8,1) | NO   |     | 0.0     |       |
+-------+--------------+------+-----+---------+-------+
3 rows in set (0.01 sec)

参考:精密計算の例

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

Toy problemsは役立たずか(Floydの問題)


U = {√1, √2, …, √50}を2つのグループに分け、グループ毎の和を求める。和の差がなるべく小さくなるようなグループ分けを、制限時間10秒で見つけよ。

goal = (√1 + √2 + … + √50) / 2 ≈ 119.517900301760392247020223113… としたときに、{√1, √2, …, √50}のサブセットで、和がgoalに最も近いものを求めればよい問題です。

1881526917Bob Floydによって提供されたというこの問題は、Knuthによって1977年に発表されました(Selected Papers on Computer Scienceに再掲されています。参考文献リストあり。索引あり。)。

学生だったときに研究室(の一部)でこの問題を話題にしたことがあります。遺伝的アルゴリズム(Genetic Algorithm, GA)と呼ばれる最適化手法を研究する人が多くいた研究室でしたが、GAではあまりよい解は見つかりませんでした。(エージェントアプローチ人工知能 第1版のp.624によれば、GAは何かをする上で3番目によい方法だそうです。)

4274200183伊庭先生の『進化論的計算手法』でこの問題が取り上げられているのは、そのような経緯のためと思われます(ちなみに、この本の第1版のまえがきによると、私は第5章の共著者ということになっていますが、第6章の間違いですね。)

すべてのグループ分けは2の49乗通りありますが、当時のコンピュータで10秒でこれを全探索するのはおそらく不可能だったはずです。当時のコンピュータでよい解を見つけるには、大量のデータを処理する方法と計算時間を見積もる方法を知る必要がありました。ですから、このようなtoy problemsからも、いろいろ学ぶことがあるというのがKnuthの主張でした。

学ぶことがあるという点は今日でも変わりませんが、コンピュータは圧倒的に速くなっているので、この問題から学べることは少なくなっているかもしれません。それでも、やったことがない人は、ちょっとプログラムを書いて試してみてください。楽しめることは間違いありません。

doubleで計算するとこんな感じになります(最良解は1つではありません)。

119.517900301760391812422312796: goal
119.517900301760320758148736786: {1,4,5,7,8,9,12,13,15,16,20,25,27,30,31,33,37,38,41,43,44,46,47,48,49}
119.517900301760462866695888806: {2,3,6,10,11,14,17,18,19,21,22,23,24,26,28,29,32,34,35,36,39,40,42,45,50}

goalの正確な値は119.517900301760392247020223113...で、doubleの計算結果は小数点以下14桁までしか一致しません。2つのグループの合計の一致は、(厳密に計算しても)小数点以下12桁までなので、けっこうぎりぎりのところで計算していることがわかります。

GCCでlong doubleで計算するとこんな感じです(最良解は1つではありません)。

119.517900301760392242633734838: goal
119.517900301760320737332055074: {1,4,5,7,8,9,12,13,15,16,20,25,27,30,31,33,37,38,41,43,44,46,47,48,49}
119.517900301760463740996520698: {2,3,6,10,11,14,17,18,19,21,22,23,24,26,28,29,32,34,35,36,39,40,42,45,50}

コードは後ほど。