逆引きnumpyメモ
numpy? ムズいよね。使い慣れるまで好きになれないと思うよ。だけど、俺は負けないよ。
え〜、numだっ、numpyが躍動する俺のpythonを、皆さんに見せたいね。
※ この記事は少しずつ加筆修正すると思います。
俺的マイナーイディオム
下記は解説記事とかであまり見ない書き方だけど、
僕の中でほとんどイディオムなやつなのでこの記事では使うね。
あと、xとyは全部np.ndarrayの行列とします。
# グラム行列ベースに何かする
# この書き方なら、例えばカーネル関数を明確に分けて定義できる。
# *は関数にぶち込む時に展開して貰うための修飾
fuga = hoge(*np.meshgrid(x, y))
# Aという行列のx軸に対してhogeを行ったものを取り出す
# map的な何か。何故か皆for文を使うけど、numpyは関数型っぽく書いたほうがわかり易くない？
fuga = np.apply_along_axis(hoge, x, A)
# イディオムというか、普通にdot積なんだけど、
# 何故皆これ使わずにnp.dotばかりなんですか？サクッと書く時は楽でいいんだけど…
hoge @ fuga
初期化系
Iを作る
この辺は流石に要らないかと思ったけど、意外とググってる自分が居るんだな。
np.eye(5)
ただの等差数列をつくる
0から100までの大きさの、互いの差が1の等差数列。
np.arange(0, 100, 1)
といっても、最後の1とかは黙ってれば1にしてくれるので
np.arange(100)
とか。ま、基本やね。
等比数列を作る
これは、1から10の3乗までを1000分割するよという意味。
適当に読み替えてくださいな。
np.logspace(1, 3, 1000, base=10.0)
基本的なところ
一番大きい数を取ってくる
np.nanmaxが良い。
この関数はNanを無視してくれる。
まぁ、Nan引っ掛けて欲しい場合もあるやろうけどね。
逆順にする
np.flipでおｋ
連結
並列というか、縦につなげるやつはnp.vstack
直列というか、横につなげるやつnp.hstack
なんだけど、この関数はnp.vstack(hoge, fuga)ではなくて
np.vstack((hoge, fuga))みたいに使うのです。
あとは、np.concatenateとかもあります。
グラム行列を作る
普通にグラム行列を作る。こんな時はmeshgridがいい。
…なんだけど、meshgrid自体はnumpy芸とも言いたくなるような
技巧たっぷりな書き方になってしまうのは僕だけ？
x = np.arange(10)
y = np.arange(5)
xg, yg = np.meshgrid(x, y)
gram = xg * yg # ここがやりたい演算
だけど、こう書くことも出来る。この場合、下記のカーネル法によく似ているという意味で
わかりやすいので、この方が僕は好き。(カーネル法は一般化内積ですしね)
from operator import mul
x = np.arange(10)
y = np.arange(5)
gram = mul(*np.meshgrid(x, y)) # ここではまだ大差ない
回帰とか
カーネル法
取り敢えずガウシアンカーネルで書いてみた。
上記グラム行列のmulを適当にカーネルで置き換えてあげればいいだけ。
ちなみに、この書き方なら色々なカーネル関数…
例えば、ガウシアンや線形等を簡単に入れ替えられるので便利。
x = np.arange(10)
y = np.arange(5)
def kernel(a: np.ndarray, b: np.ndarray, v: float) -> np.ndarray:
return np.exp(- v * np.square(a - b))
gram = kernel(*np.meshgrid(x, y), 3)
二次正則化付きの場合は
lmd = 0.1 # 適当
gram = kernel(*np.meshgrid(x, x)) + lmd * np.eye(x.shape[0])
というのを応用して、以前、ふざけて書いたカーネル回帰ゴルフ。
def solve_kernel(x: np.ndarray, y: np.ndarray,
kernel: Callable, lmd: float) -> np.ndarray:
K = kernel(*np.meshgrid(x, x)) + lmd * np.eye(x.shape[0])
alpha = np.linalg.solve(K, y)
return np.vectorize(lambda xs: np.sum(alpha * kernel(x, xs)))
カーネル回帰はPEP8守っても実質たった4行で書ける。numpyはえらいなぁ…。
三次元で山っぽいのを作ってみる
100行100列で、ガウシアンな山を作ってみる。
これにもmeshgridが使える。若干黒魔術っぽいのが悲しいが、大したことはない。
def mount(x: int, y: int) -> np.ndarray:
return np.exp(-(np.square(x) + np.square(y)) / 1000)
x, y = 100, 100
g = mount(*np.meshgrid(*(np.arange(-n/2, n/2, 1) for n in [x, y])))
何故こんなの作ったかって？クラスタリングの実験用に書いたやつです。
meshgridと*の組み合わせは黄金かも知れないし、コーディング規約違反と言われて殴られるのかも知れない。
関数型的な何か
for文が嫌な時
複数やり方がある。まぁ、mapとかは今更なのでいいでしょう。
ちなみに、mapに限らず、内包表記もpythonなので糞遅いことを自覚すること。
そんなものに頼っているようではnumpy使いとして二流である。
(ごく稀に、そのほうが速いこともありはする)
ユニバーサル関数を使う
ユニバーサル関数という、numpyっぽい挙動の関数を作れるので、それを使う方法が一つ。
ただし、速度が速くなるわけではない。
A = [[2, 4, 8], [1, 1, 1]]
def test(x):
return x * 2
np.vectorize(test)(A)
> array([[4, 8, 16], [2, 2, 2]])
ちなみに、デコレータとしても使えるけど、コーディング規約的な意味で許されるかは分からん。
ちゃんと関数になるためには引数がただの数じゃないと駄目とか、色々あるっぽい。
np.apply_along_axisを使う
この方法は軸を指定できるのでかなり応用が効く。
しかも、公式曰く速いらしい。vectorizeとどっちのほうが速いかは知らない。
# 上記と同様にして
np.apply_along_axis(test, 0, A)
ちゃんと動くためには次元の制約があったりもするらしい。
他にも関数はいくつかある。https://numpy.org/devdocs/reference/routines.functional.html
np.tileの応用
例えば、あるarrayの列一つ一つに、違う事をしていきたいときがあるとする。
こんな感じで。
matrix = np.reshape(np.arange(70), (7, 10))
nums = np.arange(10)
result = []
for r, num in zip(matrix, nums):
result.append(r * num)
このコードは愚の骨頂である。非常に遅いからだ。だけど、apply_along hogehogeもどうも思い浮かばなかった。
そこで、こうする。
matrix = np.reshape(np.arange(70), (7, 10))
nums = np.arange(10)
num_arr = tile(i, (matrix.shape[0], 1))
result = matrix * num_arr
勿論、圧倒的に速い。しかし、どうもnumpy芸である…
そこで、関数にしてしまうのはどうか。
ここでは、先に書いたnp.meshgridグラム行列風のインターフェイスを目指してみる。
def for_each_row(x, y):
return np.tile(x, (y.shape[0], 1)), y
matrix = np.reshape(np.arange(70), (7, 10))
nums = np.arange(10)
def func(x, y): return x * y
result = func(*for_each_row(nums, matrix))
ちなみに、実はこのケースではこれで良い。
matrix = np.reshape(np.arange(70), (7, 10))
nums = np.arange(10)
matrix * nums
上記のやり方は、より複雑なものを実装しないといけなくなった時に役立つ…かもしれない。
逆問題
連立方程式を解く(解ける場合)
普通の連立方程式。例えばこう。
$5x + 3 = y$
$x + 1 = y$
行列表現ではこう。
$$AX=Y$$
これを解く場合は下記のようにする。
X = np.linalg.solve(A, Y)
この場合、未知の数は変数Xだった。
しかし、実用的には測定値Xから定数Aを知りたい事がほとんどと思う。
だから、最小二乗法なんかではここのAとXが逆の立場になる。
逆行列は明示的に出さずにこの関数を使うほうが良いらしい。
最小二乗法を解く(解がない場合)
解が存在しないパターンの場合は最小二乗法を使うことになる。
「解が存在しない」とは例えば下記の様なやつ。
$8x + 4 = y$
$5x + 3 = y$
$x + 1 = y$
この連立方程式は互いに矛盾する解となるため、解がない。
最小二乗法の解は以下の式なので
$$X^TXA = X^TY$$
A = np.linalg.solve(X.T @ X, X.T @ Y)
疑似逆行列を使うことも出来る。
A = pinv(x) @ y
なんだけど、既にちゃんと関数があるから上記ではなく下記が正解。
A = np.linalg.lstsq(X, Y)[0]
正則化付き最小二乗法
正則化付きの場合はどうするのかよく分からない。
ググったけど、上記のような正解を見つけられなかったから、適当に自分で書くことになるのかな？
$$(X^TX + \lambda I)A = X^TY$$
lmd = 1
seisoku = lmd * np.eye(min(X.shape))
A = np.linalg.solve(X.T @ X + seisoku , X.T @ Y)
最小ノルム解を解く(解が無限にある場合)
このような時は解が無限にある。
$8x_1 + x_2 + 6x_3 = 4$
$5x_1 + x_2 + 2X_3 = 3$
この方程式の最適解を求めてみる。
$$A = X^T{(XX^T)}^{-1}Y$$
A = X.T @ np.linalg.inv(X @ X.T) @ Y
ってなりそうじゃん？
だけど、ムーアペンローズの逆行列は凄いので、普通にいける。
A = pinv(X) @ Y
でも、これも正解は下記。
A = np.linalg.lstsq(X, Y)[0]
この関数はすごいねー。
正則化付き最小二乗法(解が無限にある場合)
これも正攻法が分からん。だけど、式は簡単に作れる。
$$A = X^T{(XX^T + \lambda I)}^{-1}Y$$
lmd = 1
seisoku = lmd * np.eye(min(X.shape))
A = X.T @ np.linalg.inv(X @ X.T + seisoku) @ Y
うーん…(･ัω･ั)
正則化付き探せなかったのは僕の勉強不足か？
大して難しくないから良いんだけどさ…
その他、いろんな最適化問題
速くする
numpyでは**とかの演算子は比較的遅い。e**hogeとかやるべきではなくて、np.expとかを使うほうが速い。
僕の書いたWavelet変換のパッケージではそうだった。(なので可読性が悪い…)
hoge = np.arange(1000)
e ** hoge # 遅い
np.exp(hoge) # 速い
hoge ** 2 # 遅い
np.square(hoge) # 速い
hoge ** (0.5) # 遅い
np.sqart(hoge) # 速い
hoge ** (3) # 遅い
np.power(hoge, 3) # 速い
hoge ** (3.4) # 遅い
np.float_power(hoge, 3.4) # 速い
この他、np.powerとか色々あるけど、基本、最適化された関数を使うほうが速い。
楽をするために
行列の積
何故皆@演算子使わないのかよくわからない。np.dotとかより読みやすいと思うんだよな…。
遅かったりするんですか？
np.arange(10).T@np.arange(10)
複素数関連
複素共役
np.conjでおｋ
実軸虚軸関連
hoge.real, hoge.imag
という感じで書けますから、こいつを代入する時は
fuga.real, fuga.imag = hoge.real, hoge.imag
何故これが欲しいかと言うと、虚軸だけ、実軸だけに何らかの処理をすることが意外とあるから。
フーリエ変換
フーリエ変換はどっちかと言うとscipyを使うことが多いと思う。numpyにもあるんだけどね。
scipy.fftpack.fft(wave)
逆変換
comp # フーリエ変換後の複素数の行列とする
scipy.fftpack.ifft(comp)
cupy関連
cupyへ変換
cp.asarray(numpy_array)
numpyへ変換
cp.asnumpy(cupy_array)
cupy自体、numpyと違う所もそれなりにあるので、色々気をつけようね。
場合分け
場合分け
np.whereを使う。
もしもxが条件式に合致すれば第二引数を、合致しなければ第三引数を出す。
x = np.arange(10)
np.where(x > 5, x * 10, x)
> array([0, 1, 2, 3, 4, 5, 60, 70, 80, 90])
条件文
条件式をぶっこむとRっぽい挙動になる。
x = np.arange(10)
x > 5
> array([False, False, False, False, False, False, True, True, True, True])
数学の基礎っぽいやつ
行列に関して
np.det(hoge) # 行列式
np.norm(hoge) # ノルム
積分
これもscipy
これ、マジカヨすげーってなった。
from scipy.integrate import quad
def test(x):
return x ** 2 - x
print(quad, test, -5, 5)
他に、2回積分するのとかもある。それはdblquad
特殊関数
さては君は数学好きだね？
このリンクを見ると良い。
https://docs.scipy.org/doc/scipy/reference/special.html#module-scipy.special
ゼロ埋め
hogeは配列
次のは3番目以前と5番目移行を埋めるということ。
constantは定数で埋めるという意味。
最後のは左右で埋める定数。
他にも色々モードがある。
np.pad(hoge, [3, 5], 'constant' ,[0, 0])
応用
画像のくり抜き
値の大きな所だけくり抜いてみる。
さっきのようにTrue, Falseの配列を作って、条件式でくり抜く。
例ではnumpyで山みたいな物を作った上で、山のてっぺんを切り取っている。
この時に、np.nanが必要になったりする。matplotlibのimshowではnanが無視されるからだ。
def mount(x: int, y: int) -> np.ndarray:
return np.exp(-(np.square(x) + np.square(y)) / 1000)
# くり抜く関数。
# valuesが普通の行列で、clusterがboolの行列
def make_single_cluster(values, cluster):
img = np.nan * np.zeros_like(values)
img[cluster] = tvalues[cluster]
return img
x_y = 100, 100
# 今回は一寸ノイズを入れてみた
X = np.random.rand(*x_y) / 100
# 山を作る
X += mount(*np.meshgrid(*(np.arange(-n/2, n/2, 1) for n in x_y)))
plt.imshow(make_single_cluster(X, X > 0.5))
plt.show()
互換性
matlab読み書き
scipyにmatlab読み書き機能がある。
import scipy.io as sio
sio.loadmat('hoge.mat')
data = {'fuga', np.arange(100)}
sio.savemat(fuga.mat, data)
残念ながら僕はmatlab使ったことないからよく分からん。