おっぱいそん！

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

Octave

とある論文のsupplemental materialにあったMATLAB（線形代数のパッケージを詰め込んだ有料ソフト）のコードを動かしたかったので、Octave（MATLABと互換性のある無料ソフト）をインストールして使ってみた。
MAC OSX（Yosemite）へのインストールはhttps://kiskeyix.org/articles/605を参考にした。
以下で、pythonとの比較と自分用につまづきやすいポイントのメモを。


pythonと同様にインターラクティブなモードで実行する方法と、ファイルに保存して実行する方法（拡張子は.mとする）がある。
インターラクティブな方の使い方は、pythonとほとんど同じ（terminalからpython or octaveと叩いて入る）。
ファイルに保存して実行する場合はpythonと同様に、terminalでファイルを置いているフォルダまで移動して"octave file_name.m"とする方法の他に、インターラクティブなモードでfile_nameとタイプして実行する方法がある。その際、インターラクティブなモードでcdでファイルを置いているフォルダまで移動する。


また、pythonと同様にあるファイル（pythonでいうモジュール）に保存した関数を別のファイル（スクリプトファイル）から読み込んで使うことができる。ただし、その仕様が結構違う。
pythonの場合は、あるファイルmy_module.pyに保存した関数my_func()を使うには、実行するファイルで

import my_module
my_module.my_func()

とします。
もちろん、一つのモジュールファイルの中に幾つかの関数を入れることも出来た。

一方、Octaveでは、スクリプトファイルはpythonとほとんど同じだが、呼び出す関数を保存しておくファイルは関数mファイルと呼ばれ、以下のような制約がある（どちらも拡張子は.m）。
関数名とファイル名は同じでないとダメで、その中身は関数の定義から始める。
関数名とファイル名の名前が一致していないと、

warning: function name 'func_name' does not agree with function file name '/Users/uesr_name/Desktop/test/file_name.m'

とエラーがでる。
また、関数mファイルを実行することはできない（pythonの場合は.pyファイルは実行されれば、スクリプトファイルとして扱われ、importされればモジュールとして扱われていた）。
関数mファイルなのか、スクリプトmファイルなのかは、ファイルの中身が関数の定義で始まるのか、そうでないのかによって決まる。
すなわち、

function y = func_name(x)
    y = x + 1;
end

みたいな感じで、functionで始まるかどうかで決まる。
従って、スクリプトファイルの中で関数を定義するのに、ファイルの一行目から関数を定義しようとすると、関数mファイルとみなされてしまって、実行する際にエラーがでる（僕はここでつまづいた）。
例えば、上のコードをfunc_name.mというファイルの保存して実行しようとすると、

error: 'x' undefined near line 2 column 13

とエラーがでる。

スクリプトmファイルから関数mファイルを呼び出すには、例えば、上の関数だと、

a = 1;
b = func_name(a);
disp(b);

とする。
ここで、どのファイルから関数を呼ぶか指定していないが、上に書いたように、関数名と関数mファイルの名前は一致させているので、スクリプトファイルと同じフォルダに関数mファイルをおいておけば自動的に探してくれる。

Ref:
http://www.arc.hokkai-s-u.ac.jp/~kusiyama/Inf1_Matlab/Grammar_4.htm
データ解析ツールoctaveを語ろう Part 2


Octaveのスクリプトmファイルと関数mファイルの仕様は、pythonの使用を知っていると、不便に感じてしまう（一つの関数mファイルに複数の関数を定義できないのと、同じ.mファイルでスクリプトmファイルと関数mファイルの両方の機能をもたせられない）。

行列のランクと固有値の関係

線形代数

行列のrankとnonzeroの固有値の関係について、すぐに忘れるのでメモしておく。
結論から先に書くと、一般の正方行列Aについて、

rank A ≧ Aの（重複度も含めた）nonzeroの固有値の数

が成り立つ。
特に、対角化可能な行列の場合には

rank = nonzeroの固有値の数

が成り立つ。
以下で、証明を書いておく。

続きを読む

pythonでフィッティングをする

numpy scipy

pythonでfittingをする方法。

例えば、
{ \displaystyle f(x) = a + b x + c x^2 }
という{\displaystyle a,b,c }をパラメータとする関数でデータ点{\displaystyle \{ x_i, y_i \}_{i=1,2,...N} }を
{ \displaystyle \sum_{i=1}^{N} ( y_i - f(x_i) )^2 }
が最小になるようにfittingしたいとする（最小二乗法）。
scipy.optimizeのcurve_fitを使うのが楽（scipy.optimizeにはleastsqという関数もあり、こちらでも同じことができるが、curve_fitの方が分かりやすい）。

import numpy as np
import scipy.optimize
import matplotlib.pylab as plt

# data which you want to fit
xdata = np.array([0.0,1.0,2.0,3.0,4.0,5.0,6.0])
ydata = np.array([0.1,0.9,2.2,2.8,4.2,5.9,7.4])
# initial guess for the parameters
parameter_initial = np.array([0.0, 0.0, 0.0]) #a, b, c
# function to fit
def func(x, a, b, c):
    return a + b*x + c*x*x

paramater_optimal, covariance = scipy.optimize.curve_fit(func, xdata, ydata, p0=parameter_initial)
print "paramater =", paramater_optimal

y = func(xdata,paramater_optimal[0],paramater_optimal[1],paramater_optimal[2])
plt.plot(xdata, ydata, 'o')
plt.plot(xdata, y, '-')
plt.show()

xdata, ydataはfittingしたいデータ点（上の例では手で作ってるが、普通は実験などで得られる）。
parameter_initialはパラメータの初期値（なくてもいい）。
funcがfittingする関数で、第一引数がxで、第二引数以降がパラメータ。返り値はy。
paramater_optimalがfittingで最適化されたパラメータで、covarianceは共分散。

結果は

paramater = [ 0.14761905  0.70357143  0.08452381]

f:id:nohzen:20150928224255p:plain

続きを読む

関数で複数の返り値を返す方法

pythonにおいて関数で複数の値を返すと、デフォルトではタプルになる。

def f_tuple(x):
  plus = x + 1
  cubic = x ** 3
  return (plus, cubic)
result = f_tuple(2)
print result[1]

returnの後の(plus, cubic)はカッコをのけても同じ動きをする。
また、受け取る時に、plus, cubic = f_tuple(2)のように受けても良い。


変数が少ない場合には上の書き方でよいが、返す値の数が多くなると、どれがどの返り値か分からなくなる。
上の例ではresult[1]で３乗した結果の値を取り出しているが、インデックスで指定するやり方だと、その要素が何の値なのか分かりにくい。
そこで、タプル以外で返り値を返すことを考えてみる。

続きを読む

Named tuple

pythonの組み込みのコンテナ型（=コレクション型）にはリスト、タプル、ディクショナリがある。
この記事では、名前付きタプルというコンテナ型を紹介したい。（名前つきタプルはクラスの書き方をふんだんに使っているので、クラスを知らない人はクラスを先に勉強した方がいいと思います。この記事でもクラスの基本的なことは前提知識としています。）
名前付きタプルはタプルの子クラスで、インデックスだけでなく、クラスのように変数名（ドット表記）でも要素にアクセスできるコンテナ型です。

使用例

例えば、以下のように使用できる。

from collections import namedtuple
from math import sqrt

Vector = namedtuple('Vector', 'x, y')
vec1 = Vector(x=1.0, y=5.0)
vec2 = Vector(2.5, 1.5)
distance = sqrt((vec1.x-vec2.x)**2 + (vec1.y-vec2.y)**2)
print "distance =", distance

使用するには、collectionsからimportする必要がある。
namedtuple()はクラスを作る関数で、引数はクラス名(慣習として、クラス名は大文字で始める)、要素（インスタンス変数）の名前。
要素の名前の部分は、'x y'や['x', 'y']のように書いてもよい。
また、要素の呼び出しはsqrt((vec1[0]-vec2[0])**2 + (vec1[1]-vec2[1])**2)と、普通のタプルのようにインデックスでも出来る。

長所

いいところ：

  • タプルより要素の取り出しが分かりやすい（ドットを使ったオブジェクト的な書き方が出来る）。
  • 辞書より高速でメモリが節約できる（メモリの使用量は普通のタプルと同じ）。

わるいところ：

  • import correctionがいる。

普通のタプルのような使い方もできるので、元々、普通タプルで書いていたコードを書き換えずにnamed tupleに、一部だけをnamed tupleにすることも出来ます。

続きを読む

numpyの型

numpy

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

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

ただし、（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.float64やnumpy.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次独立でなければ、対応するジョルダンブロックの非対角成分に１が現れるが、数値計算においてとても小さい誤差が入り、固有値の縮退がとけると、その瞬間に固有ベクトルは１次独立となり、非対角成分がゼロになってしまうからである。つまり、数値計算における小さい誤差で、ジョルダン標準形の値が急激に変わるので、ジョルダン標準形は数値計算的に安定でないのだ。
参考：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）ので十分条件にしかなっていないと思う。つまり、２つの行列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つの行列が相似であるかどうかは、小さい数値誤差で変わりうるので、浮動小数点数計算で確認するという考え自体がナンセンスな気がしてきた。
確認したい行列が誤差なく与えられているなら、代数計算でジョルダン標準形を使って確認できるし、確認したい行列が既に誤差を含んでいるなら、相似であるかどうかを確認するということ自体が意味がない（誤差で相似かどうかが変わりうるので）。