おっぱいそん!

pythonを数値計算に使うときのテクニックとかをまとめていきたい。

numpyの型

自分用のメモにまとめておく。

numpyの基本的な型は以下の5つ。

ただし、(bool型を除いた)データ型はそれぞれ異なるサイズがああります。
たとえば、int型なら、

  • numpy.int8
  • numpy.int16
  • numpy.int32
  • numpy.int64

4種類があります(intの後の数字はbit数を表しています)。単にnumpy.intとすると、環境に依って、numpy.int32もしくはnumpy.int64を表します(最近のPCだと、ほとんどnumpy.int64になっていると思います)。
同様に、uint,float,complexにも、異なるサイズのものがあり、単にnumpy.floatやnumpy.comolexとすると、それぞれ、numpy.float64numpy.complex128を表します。

Pythonのビルトインにもint, float, complexがありますが、numpy.int(=numpy.int32 or numpy.int64), numpy.float(=numpy.float64), numpy.complex(=numpy.complex128)はこれらと同じもののようです。
一方、numpy.int16やnumpy.float128などは、対応するpythonの型はありません。
numpyの中では、int, float, complexと書いても、numpy.int, numpy.float, numpy.complexと書いても問題ありませんが、たぶん、前者のように書いた方が良いと思われます。もちろん、ビルトインと違うサイズを使いたいときは、numpy.float128のように書く必要があります(pure pythonのfloat128は存在しません)。

import numpy as np

print np.array([1], dtype=int).dtype #int64 or int32
print np.iinfo(np.int64)

print np.array([1], dtype=float).dtype #float64
print np.finfo(np.float64)

print np.array([1], dtype=complex).dtype #complex128

dtype=floatの部分は上記の説明のように、dtype=np.floatと書いても、dtype=np.float64と書いても同じです。

続きを読む

2つの行列が相似であることの確認

2つの正方行列A,Bが与えられた時に、この2つの行列が相似、つまり、ある正則行列P(Pはユニタリ行列でなくてもよい)が存在して、
A = P B P^(-1)
と書けるかどうかを、数値的に確かめたい。
(結論から言うと、僕はどうすればいいか分かりません。いい解決法を知っている方は教えて下さい)
[2018.05.03 記事の末尾に、自分の結論を追加した]

解析的に

解析的には、2つの正方行列A, Bのジョルダン標準系を計算して、比べれば良い。
これは、「A, Bのジョルダン標準形が一致すること(ただし、ジョルダンブロックの並べ方の不定性があるので、適当に並べ替える)」と「A, Bが相似であること」は同値であるからである。
(証明は、ジョルダン標準形の一意性から簡単に示せるので略。)
(ちなみに、ユニタリ同値であるかどうかを確かめたい時にも、ジョルダン標準形が使え、ジョルダン標準形が一致することに加えて、A, Bをジョルダン標準形にする正則行列P_A*P_B^(-1)がユニタリーがどうかを確かめればいい。)

数値的に

しかしながら、ジョルダン標準形と数値計算は相性が悪い。
というのも、ある考えている行列に縮退した固有値があり、固有ベクトルが1次独立でなければ、対応するジョルダンブロックの非対角成分に1が現れるが、数値計算においてとても小さい誤差が入り、固有値の縮退がとけると、その瞬間に固有ベクトルは1次独立となり、非対角成分がゼロになってしまうからである。つまり、数値計算における小さい誤差で、ジョルダン標準形の値が急激に変わるので、ジョルダン標準形は数値計算的に安定でないのだ。
参考:http://en.wikipedia.org/wiki/Jordan_normal_form#Numerical_analysis
実際、ジョルダン標準形は数値的に安定でないため、numpyなどに関数は用意されていません。代数計算をしてくれるSympyにはジョルダン標準形を求めてくれる関数があるみたいです。
参考:http://stackoverflow.com/questions/20312216/compute-jordan-normal-form-of-matrix-in-python-numpy


そこで、数値計算においては他の方法を使いたい。

まず、特性多項式(charactaristic polynominal)であるが、これは相似であることの必要条件であるが、十分条件でない。つまり、相似であれば、特性多項式が一致する(証明はdetの性質より簡単)ことは言えるが、特性多項式が一致したからといって、相似であるとは限らないので、使えない。

http://math.stackexchange.com/questions/14075/how-do-i-tell-if-matrices-are-similar
上の知恵袋の記事や、Wikipediaの記事によると、Schur decompositionなるものが数値的には安定で使えると書かれているが、Schur decompositionは分解した際の、上三角行列が一意的でない(http://ja.wikipedia.org/wiki/%E3%82%B7%E3%83%A5%E3%83%BC%E3%83%AB%E5%88%86%E8%A7%A3)ので十分条件にしかなっていないと思う。つまり、2つの行列A,BをSchur decompositionをした時に、上三角行列が一致していれば相似であるが、逆は言えない。この上三角行列をさらに何らかの一意的な標準系みたいなものに持っていければ、使えるのだが...
参考:http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.linalg.schur.html

行列の標準形にはジョルダン標準形以外のものもあって、Frobenius normal formやWeyr canonical formというものもあるらしいのですが(Wikipediaによると、標準形は一意的なので使えそう)、数値的に安定かどうかは知りません。また、numpyなどに関数が用意されていないようです。
参考:http://en.wikipedia.org/wiki/Canonical_form#Linear_algebra

ということで、行列の相似を数値的に確かめる方法を調べた結果、諦めました…
いい方法を知っている方は教えて下さい。

[2018.05.03追記]
よく考えてみたら、2つの行列が相似であるかどうかは、小さい数値誤差で変わりうるので、浮動小数点数計算で確認するという考え自体がナンセンスな気がしてきた。
確認したい行列が誤差なく与えられているなら、代数計算でジョルダン標準形を使って確認できるし、確認したい行列が既に誤差を含んでいるなら、相似であるかどうかを確認するということ自体が意味がない(誤差で相似かどうかが変わりうるので)。

MacBook Pro (Mid 2012, 13インチ)のSATA cableの交換

2週間ほど前に、MACのOSが立ち上がらなくなりました。
結論から言うと、SATAケーブルを交換したら治ったのですが、その記録をここに書いておきます。
*特にMACBookPro Mid 2012モデルの方向けの追記をしました。

時系列

  • MACを久々に再起動させようとすると、OSが起動しなくなってしまいました。最初はソフト的な問題かと思い、PRAMのリセットやディスクの修復などWebにある対処法を一通りやってみましたが、ダメでした。
  • 次に、データやアプリケーションは消さずに(HDDをフォーマットせずに)OSXのみの再インストールを試みましたが、これも失敗しました。そこで、諦めて、HDDをフォーマットしてしまうことにしました。重要なデータはDropBox等に置いているのですが、バックアップを取っていないデータもあったので、救出しました。HDDから、データを救出する方法は、例えば、下のようなやり方があると思います。僕の場合は、1.の方法では外付けHDDとして認識されなかったので、方法2.でデータを取り出しました。
  1. MACの中のHDDを外付けのHDDとして、動かす(具体的には、他のMACとThunderboltケーブルでつないで、ターゲットディスクモードで他のMACからHDDのデータを取り出す)。詳しくは「ターゲットディスクモード Thunderbolt」などでググってください。
  2. リカバリーモードで使えるターミナルで、MACの中のHDDのデータをUSBメモリなどに移す。http://techracho.bpsinc.jp/morimorihoge/2014_07_09/18182を参照。上のリンクでは、ネットワーク経由でデータを吸いだしていますが、USBメモリや外付けHDDなどを使うほうが楽だと思います。また、最初、ターミナルからdfコマンドで、Macintosh HDが見えなくても、暗号化されている場合、ディスクユーティリティから、パスワードを入れることで、マウントできます(ディスクユーティリティからHDDが見れない場合はこの方法は無理だと思います)。この方法では、ターミナルからコマンドを使って、データを移動させないとダメです。シェルコマンドに慣れてない人はググる
  3. MacBookからHDDを物理的に取り出して、SATA-USBケーブルなど(内蔵のHDDを外付けにできるケースを買うのが楽かもしれません)を使って、他のMACに差して、データを取り出す。
続きを読む

2つのarrayが近似的に同じか比べる

2つのNdarrayが近い値かどうかを比べるには、numpy.allcloseを使う。

  • numpy.allclose(a, b, rtol=1e-05, atol=1e-08)

2つのarray a,bを受け取って(a,bは同じshapeのarrayでないといけない)、すべての成分が
absolute(a - b) ≦ (atol + rtol * absolute(b))
を満たせばTrueを、満たさない成分があればFalseを返す。

成分ごとに、比較したい場合はnumpy.iscloseを使う。
また、近い値かどうかではなく、exactに一致しているかどうかを調べたい時には、numpy.array_equalもしくはnumpy.array_equivを使う(numpy.array_equalはshapeが一致している時のみTrueで、numpy.array_equivはbroadcastedして一致するならTrue)。

Ref:
http://ibisforest.org/index.php?python%2Fnumpy#q0f5fdb7
http://docs.scipy.org/doc/numpy/reference/generated/numpy.allclose.html
http://docs.scipy.org/doc/numpy/reference/generated/numpy.isclose.html
http://docs.scipy.org/doc/numpy/reference/generated/numpy.array_equal.html
http://docs.scipy.org/doc/numpy/reference/generated/numpy.array_equiv.html

例)

import numpy as np
A = np.array([1.0, 1.0])
B = np.array([1.0000001, 1.0])
print np.allclose(A,B)
print np.isclose(A,B)
print np.array_equal(A,B)

True
[ True True]
False

複素数型が実数に近かったら、実数型にする

complex型のndarrayが与えられた時に、もし複素数の部分が小さいなら、floatなどの実数の型で近似する方法。
そういう関数、ありそうだと思ってググったらNumpyの関数にnumpy.real_if_closeというのがあった。

numpy.real_if_close(a, tol=100)

この関数はNdarray aの虚部a_imagを見て、それがa_imag≦マシンイプシロン*tol(Tolerance)なら、虚部を無視した実部だけのNdarrayを返す。
a = a_real + i*a_imag → a = a_real

マシンイプシロン
np.finfo(np.float).eps
で調べることができ、Pythonでは普通、2.22044604925e-16になっているみたいです。
tolはデフォルトでは100に設定されているので、a_imag≦2.22044604925e-14であれば虚部を無視することになります(tolはオプションで変更できます)。


このように、np.real_if_closeは虚部の絶対値をみて、小さいかどうかを判断する仕様になっている。
しかし、普通aの絶対値が大きくなると、誤差の絶対値も大きくなるので、虚部/絶対値が小さいかどうかを見るほうが、良い場合が多いと思う。:

print np.finfo(np.float).eps
a_real = np.array([1.0e4])
a_imag =  np.array([3e-12j])
a = a_real + a_imag
print a.dtype
a = np.real_if_close(a, tol=100*np.absolute(a))
print a.dtype

2.22044604925e-16
complex128
float64

Numpyは欲しい関数がたいてい用意されているのですごい!
http://docs.scipy.org/doc/numpy/reference/generated/numpy.real_if_close.html
http://docs.scipy.org/doc/numpy/reference/generated/numpy.absolute.html

2015.6.26 追記
オプションtolで1より小さい値を設定すると、real_if_closeは配列を常に実数型にしてしまうようです(本来は複素数型のままにして欲しい)。
そもそも、tolが1より小さいとa_imagがゼロでない値を持っていたら、最小値はマシンイプシロンなので、常に複素数のままになるべきですが、たぶんそういうtolを想定して関数が作られておらず、バグってしまうのだと思います。
上のtol=100*np.absolute(a)という設定はaの絶対値が大きい場合を想定していますが、逆に絶対値が小さい場合は使うときに注意してください。

Pythonの代入の挙動

Pythonで代入をした時に、各変数がどのオブジェクトを指すかをまとめておく。特に、list(リスト)とNdarray(配列)の違いをば。
前回の記事も参照 : http://oppython.hatenablog.com/entry/2015/01/15/014729

例1

前回の記事でも書いたようにPythonでは、代入(=)は変数に値(=オブジェクト)を格納しているのではなく(値渡し)、値(=オブジェクト)への参照を格納している(参照渡し)。

x = np.array([1])
y = x
print 'id(x)_before =', id(x)
print 'id(y)_berore =',  id(y)
x = np.array([2])
print 'id(x)_after =', id(x)
print 'id(y)_after =',  id(y)
print y

id(x)_before = 4300100832
id(y)_berore = 4300100832
id(x)_after = 4300101152
id(y)_after = 4300100832
[1]

y = xとすると、変数yもxと同じオブジェクトnp.array([1])を指すようになる(参照渡し)。
変数xに新しいオブジェクトnp.array([2])を代入すると、xは新しいオブジェクトを参照するようになる。
もちろん、yには何もしていないので、始めのオブジェクトを参照したままである。
上の挙動は、ndarrayだけでなく、タプル、文字列、リスト、数値などでも共通(型が変更可能かどうかに依らない)。

続きを読む

Pythonの型の分類と代入

変数が指すオブジェクトがどれかなどの振る舞いをすぐ忘れてしまうのでメモしておく。

組み込み型の分類

まず、Pythonの組み込み型についてまとめておく。

管理方式

管理方式の違いで分類すると、以下のようになる。

  • コレクション型(=コンテナ型):複数のデータをまとめたもの

+ シーケンス型:要素を順番に並べて管理。整数のindexでデータを管理できる。
例)文字列、リスト、タプル
+ マップ型:キーと値のペアで要素を管理。
例)ディクショナリ

  • 非コレクション

例)数値型(int, float, long, complex)、ブール型、None

文字型がシーケンス型なのは面白い。シーケンスなので、たとえば、

a = "oppai"
print a[4] #i

のようにスライスを使って、要素の指定とかができる。

変更可能性

変更可能性で分類すると、以下のようになる。

  • 変更可能 (mutable)

例)リスト、ディクショナリ

  • 変更不能 (immutable)

例)数値型、文字列、タプル

基本的に、代入なりなんなりをした時の挙動は、変更可能かどうかで決まる。
オブジェクトはメモリに保存されているが、変更可能ならば、そのメモリに書かれているオブジェクトは書き換え可能になる。

続きを読む