Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
5
Help us understand the problem. What is going on with this article?
@uesseu

僕の逆引きnumpy・scipyメモ

逆引き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)

とか。ま、基本やね。

逆順にする

np.flipでおk

連結

並列というか、縦につなげるやつはnp.vstack
直列というか、横につなげるやつnp.hstack
なんだけど、この関数はnp.vstack(hoge, fuga)ではなくて
np.vstack((hoge, fuga))みたいに使うのです。
あとは、np.concatenateとかもあります。

グラム行列を作る

普通にグラム行列を作る。こんな時はmeshgridがいい。
…なんだけど、meshgrid自体はnumpy芸とも言いたくなるような
技巧たっぷりな書き方になってしまうのは僕だけ?

xg, yg = np.meshgrid(x, y)
gram = xg * yg

だけど、こう書くことも出来る。この場合、下記のカーネル法によく似ているという意味で
わかりやすいので、この方が僕は好き。(カーネル法は一般化内積ですしね)

from operator import mul
gram = mul(*np.meshgrid(x, y))

カーネル法

取り敢えずガウシアンカーネルで書いてみた。
上記グラム行列のmulを適当にカーネルで置き換えてあげればいいだけ。

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, x))

二次正則化付きの場合は

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))

今の僕はこれを使っている。

連立方程式を解く(解ける場合)

普通の連立方程式。例えばこう。
$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でおk

実軸虚軸関連

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)  # ノルム

ゼロ埋め

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()

Figure_1.png

matlab読み書き

scipyにmatlab読み書き機能がある。

import scipy.io as sio
sio.loadmat('hoge.mat')

data = {'fuga', np.arange(100)}
sio.savemat(fuga.mat, data)

残念ながら僕はmatlab使ったことないからよく分からん。

積分

これも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

5
Help us understand the problem. What is going on with this article?
Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
uesseu
「君は絶対にエンジニアどころか、工学部の学生にすら勝てないのだ。何故ならば、プログラミングを習っていないから」私はそう仰せつかりました。

Comments

No comments
Sign up for free and join this conversation.
Sign Up
If you already have a Qiita account Login
5
Help us understand the problem. What is going on with this article?