Help us understand the problem. What is going on with this article?

僕の逆引きnumpy・scipyメモ

逆引きnumpyメモ

numpy? ムズいよね。使い慣れるまで好きになれないと思うよ。だけど、俺は負けないよ。
え〜、numだっ、numpyが躍動する俺のpythonを、皆さんに見せたいね。
※ この記事は少しずつ加筆修正すると思います。

俺的マイナーイディオム

下記は解説記事とかであまり見ない書き方だけど、
僕の中でほとんどイディオムなやつなのでこの記事では使うね。

# グラム行列ベースに何かする
# *は関数にぶち込む時に展開して貰うための修飾
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とかは今更なのでいいでしょう。

ユニバーサル関数を使う

ユニバーサル関数という、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

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

普通の連立方程式。例えばこう。
$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.lonalg.solve(X.T @ X, X.T @ Y)

なんだけど、既にちゃんと関数があるから上記ではなく下記が正解。

A = np.linalg.lstsq(X, Y)

最小ノルム解を解く(解が無限にある場合)

このような時は解が無限にある。
$8x_1 + x_2 + 6x_3 = 4$
$5x_1 + x_2 + 2X_3 = 3$
$$A = X^T{(XX^T)}^{-1}Y$$

A = np.lonalg.solve(X.T @ X, X.T @ Y)

ってなりそうじゃん?でも、これも下記が正解。

A = np.linalg.lstsq(X, Y)

この関数はすごいねー。

その他、いろんな最適化問題

https://docs.scipy.org/doc/scipy/reference/optimize.html#module-scipy.optimize

速くする

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

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
Comments
No comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  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
ユーザーは見つかりませんでした