逆引き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使ったことないからよく分からん。