100 numpy exercisesの解説 1~50 - minus9d’s diary の続きです。コードの難易度が高く、一問一問を理解するための調査量が多すぎて、この記事で完了させることができませんでした。 numpy-100/100 Numpy exercises.ipynb at master · rougier/numpy-100 · GitHub を片手にご覧ください。
76問目～100問目は 100 numpy exercisesの解説 76~100 - minus9d’s diary に書きました。
52. Consider a random vector with shape (100,2) representing coordinates, find point by point distances
100点間の距離をすべて求めるという問題です。ここでは簡単化のために3点で考えます。
>>> Z = np.array([[1,1],[3,1],[2,5]])
np.at_least_2d()は、入力されたndarrayの次元数が最低でも2になるように次元を拡張する関数のようです。
3点のXとYを分離するのに、普通だと以下のようになるところが、
>>> X, Y = Z[:,0], Z[:,1] >>> X.shape, Y.shape ((3,), (3,))
np.at_least_2d()を使うと、以下のようになります。
>>> X, Y = np.atleast_2d(Z[:,0], Z[:,1]) >>> X.shape, Y.shape ((1, 3), (1, 3))
次に、(X-X.T)**2の部分です。Xはarray([[1, 3, 2]])で1x3の形、X.Tはその転置なので
array([[1], [3], [2]])
で3x1の形です。
X-X.Tでは1x3の行列から3x1の行列を引くことになります。shapeが一致していないのに引き算なんてできるのかと思いますが、問題なくできます。これはNumPyのbroadcastingと呼ばれる性質に基づくものです。1.3.2. Numerical operations on arrays — Scipy lecture notes の “1.3.2.3. Broadcasting” 冒頭にあるイラストを見ると演算の中身が理解しやすいと思います。
結果、X-X.Tでは以下の結果が得られます。
array([[ 0, 2, 1], [-2, 0, -1], [-1, 1, 0]])
ここまでくればD = np.sqrt((X-X.T)**2 + (Y-Y.T)**2)の演算はできたも同然で、結果は
array([[ 0. , 2. , 4.12310563], [ 2. , 0. , 4.12310563], [ 4.12310563, 4.12310563, 0. ]])
となります。
この問題、自分だったら以下のように書きたくなりそうです。
>>> Z = np.array([[1,1],[3,1],[2,5]]) >>> X, Y = Z[:,0], Z[:,1] >>> X.shape, Y.shape ((3,), (3,)) >>> X_dist = np.subtract.outer(X, X) >>> Y_dist = np.subtract.outer(Y, Y) >>> X_dist.shape, Y_dist.shape ((3, 3), (3, 3)) >>> np.hypot(X_dist, Y_dist) array([[ 0. , 2. , 4.12310563], [ 2. , 0. , 4.12310563], [ 4.12310563, 4.12310563, 0. ]])
61. Find the nearest value from a given value in an array
以下の回答ではarrayとして長さ10の1次元arrayを例にとっています。flat()は任意次元のarrayを1次元とみなしてアクセスできるようにする関数ですが、この例の場合もとから1次元なので、Z.flatを使う必要性がなくなっています。
Z = np.random.uniform(0,1,10) z = 0.5 m = Z.flat[np.abs(Z - z).argmin()] print(m)
そこで、Zを
Z = np.random.uniform(0,1,(10,10))
とすると、Z.flatの有り難みがわかるようになります。
62. Considering two arrays with shape (1,3) and (3,1), how to compute their sum using an iterator?
np.nditer()は多次元の配列に対して要素をなめるときに使う関数で、numpy.nditer — NumPy v1.12 Manual を見てわかる通り、非常に多機能です。
まず基本形として、2x2の行列2つをイテレートする例を見てみます。
A = np.array([[1,2],[3,4]]) B = np.array([[5,6],[7,8]]) it = np.nditer([A,B]) for x, y in it: print("x y:", x, y)
を実行すると、
x y: 1 5 x y: 2 6 x y: 3 7 x y: 4 8
となります。
次に、上のコードで、Aが3x1の行列, Bが1x3の行列と変えることを考えます。このときも、Q52と同様、Broadcastingのルールが働きます（参考：NumPy Iterating Over Arrayの"Broadcasting Iteration"）。
A = np.array([[1],[2],[3]]) B = np.array([[5,6,7]]) print(A.shape) print(B.shape) it = np.nditer([A,B]) for x, y in it: print("x y:", x, y)
実行結果は以下です。
x y: 1 5 x y: 1 6 x y: 1 7 x y: 2 5 x y: 2 6 x y: 2 7 x y: 3 5 x y: 3 6 x y: 3 7
ここまで理解した上で、元のコードを見てみます。
A = np.arange(3).reshape(3,1) B = np.arange(3).reshape(1,3) it = np.nditer([A,B,None]) for x,y,z in it: z[...] = x + y print(it.operands[2])
A, Bに加えてNoneをも同時にイテレートして、値の代入までしてしまっているのがなんとも不思議なコードです。
63. Create an array class that has a name attribute
ndarrayをサブクラス化して自作クラスを作成する例です。Subclassing ndarray — NumPy v1.12 Manualに公式の説明がありますが、読んでいません。既存クラスのサブクラス化はハマりどころが多そうで、自分でやる気にはなれません。
64. Consider a given vector, how to add 1 to each element indexed by a second vector (be careful with repeated indices)?
np.bincount()は、値の生起回数をカウントする関数です。
まず基本形。np.bincount([3,5,3,8])を実行すると、array([0, 0, 0, 2, 0, 1, 0, 0, 1])が返ります。3が2回登場するので、3番目のindexが2になっていることがわかります。また、返る配列の長さは、登場する値の最大値で決まることもわかります。
これをふまえて問題文のコードを見てみます。I = np.random.randint(0,len(Z),20)でarray([8, 3, 8, 2, 0, 2, 4, 4, 8, 3, 8, 6, 7, 3, 5, 4, 0, 7, 8, 6])が生成されたものとします。0から9までの値をランダムに20個並べたものですが、今回はたまたま9が出なかったことに注意してください。
引数無しで単にnp.bincount(I)とすると、さきほど見たとおり、各数値をカウントした結果であるarray([2, 0, 2, 3, 3, 1, 2, 2, 5])が返ります。
ここでもしZ += np.bincount(I)とすると、Zの長さは10なのに対してnp.bincount(I)の長さは9のため、ValueErrorが起こります。
解答例のように、np.bincount(I, minlength=len(Z))とminlength=10と指定することで、np.bincount()が返す配列の長さが最低でも10確保されるようになり、ValueErrorを防げます。
次に、別解であるnp.add.at(Z, I, 1)を読み解きます。numpy.ufunc.at — NumPy v1.12 Manualによると、この関数呼び出しはZ[I] += 1と等価です。ここで、Iはインデックス列です。
65. How to accumulate elements of a vector (X) to an array (F) based on an index list (I)?
- に引き続き
bincount()に関する出題です。
X = [1,2,3,4,5,6] I = [1,3,9,3,4,1] F = np.bincount(I,X) print(F)
のnp.bincount(I,X)は、np.bincount(I, weights=X)と等価です。つまり、値の生起回数に重みを与えます。実行結果は[ 0. 7. 0. 6. 5. 0. 0. 0. 0. 3.]です。
66. Considering a (w,h,3) image of (dtype=ubyte), compute the number of unique colors
F = I[...,0]*(256*256) + I[...,1]*256 +I[...,2]により(R, G, B)の値をR * 256 * 256 + G * 256 + Bという単一の値に変換したあと、np.unique()で同一の色を持つ要素を削除して、残った要素数を数えています。
C++のunique()は隣り合う同一の数値しか削除しませんが、np.unique()の場合は隣り合わない同一の数値も削除するようです。
69. How to get the diagonal of a dot product?
第1の方法 np.diag(np.dot(A, B)) は愚直にAとBの積をとり、その対角成分を取得しています。対角成分以外の要素も計算しているので無駄が多いと言えます。
第2の方法 np.sum(A * B.T, axis=1) では、要素ごとの積をとる*を使用しています。紙と鉛筆で確かめるとわかりますが、AとBの積をとったときの対角成分に相当する部分のみを演算しているので高速化が期待できます。
第3の方法 np.einsum('ij,ji->i', A, B) は初めて見ました。アインシュタインの縮約記号にのっとって演算を行うための関数だそうです。詳細は【einsum】アインシュタインの縮約記法のように使えるnumpyの関数。性能と使い方を解説。 - プロクラシスト などをご覧ください。こちらも
実際に時間を計測した結果は以下の通りです。
>>> import timeit >>> test_string = "import numpy as np; size = 5; A = np.random.uniform(0,1,(size,size)); B = np.random.uniform(0,1,(size,size));" >>> timeit.timeit("np.diag(np.dot(A, B))", setup=test_string) 2.1652195840001696 >>> timeit.timeit("np.sum(A * B.T, axis=1)", setup=test_string) 3.46571562500003 >>> timeit.timeit("np.einsum('ij,ji->i', A, B)", setup=test_string) 0.8595781789999819
第3の方法がもっとも高速というのはよいとして、第1の方法より第2の方法が遅いのは予想外でした。そこで、行列のサイズを100x100に変えて再測定してみます。
>>> import timeit >>> test_string = "import numpy as np; size = 100; A = np.random.uniform(0,1,(size,size)); B = np.random.uniform(0,1,(size,size));" >>> timeit.timeit("np.diag(np.dot(A, B))", setup=test_string) 18.509561333999955 >>> timeit.timeit("np.sum(A * B.T, axis=1)", setup=test_string) 15.581413918999715 >>> timeit.timeit("np.einsum('ij,ji->i', A, B)", setup=test_string) 9.336061797000184
第2の方法が第1の方法より早くなりましたが、計算量から予想されるほどの差は出ませんでした。私の環境だと、第1の方法の実行中は4つの物理コアがすべてぶん回るのに対し、第2の方法と第3の方法では1つの物理コアしか回らないことがtopコマンドで確かめられたので、第1の方法には何らかの最適化が効いているのかもしれません。あくまで推測です。
71. Consider an array of dimension (5,5,3), how to mulitply it by an array with dimensions (5,5)?
(5,5)のshapeを持つBに対してB[:,:,None]とすると、(5,5,1)のshapeになるようです。
72. How to swap two rows of an array?
>>> A = np.arange(25).reshape(5,5) >>> A array([[ 0, 1, 2, 3, 4], [ 5, 6, 7, 8, 9], [10, 11, 12, 13, 14], [15, 16, 17, 18, 19], [20, 21, 22, 23, 24]])
上記のような(5,5)のshapeを持つarrayAに対して、A0,1というアクセスを行うと
>>> A[[0,1]] array([[0, 1, 2, 3, 4], [5, 6, 7, 8, 9]])
が返ります。このアクセス方法は"Advanced Indexing"と呼ばれる方法です（参考：Indexing — NumPy v1.12 Manual）。ふつうのindexingであるA[(0,1)]や、それの簡略系であるA[0,1]との違いに注意してください。対比のために書いておくと
>>> A[0,1] 1
です。Advanced Indexingについてまだ勉強不足ですが、公式リファレンスよりはとっつきやすい解説が NumPy Advanced Indexing にあります。slicingだと規則正しい
A0,1と一見ほぼ同等の結果を、以下のようにslicingを使って得ることもできます。
>>> A[0:2,:] array([[0, 1, 2, 3, 4], [5, 6, 7, 8, 9]])
ただし、advanced Indexingでは常にコピーが返るのに対し、slicingではビューが返るという違いがあるので注意です。
0行目と1行目を入れ替える、slicingを使った別解を以下に示します。
>>> A = np.arange(25).reshape(5,5) >>> A[0:2,:] = A[1::-1,:] >>> A array([[ 5, 6, 7, 8, 9], [ 0, 1, 2, 3, 4], [10, 11, 12, 13, 14], [15, 16, 17, 18, 19], [20, 21, 22, 23, 24]])
時間を計測すると、advanced indexingよりslicingのほうが高速でした。これがコピーとビューの差でしょうか？
>>> timeit.timeit("A[0:2,:] = A[1::-1,:]", setup="import numpy as np; A = np.arange(25).reshape(5,5)") 0.7406019989994093 >>> timeit.timeit("A[[0,1]] = A[[1,0]]", setup="import numpy as np; A = np.arange(25).reshape(5,5)") 3.9081818149988976
73. Consider a set of 10 triplets describing 10 triangles (with shared vertices), find the set of unique line segments composing all the triangles
告白すると問題の意味がよく分かりませんでした。問題文には「10個の三角形を構成する10組のtriplet」を考えよ、と述べていると思いますが、想定解
faces = np.random.randint(0,100,(10,3))
で生成している10組の数値はランダムで、当然ながら三角形を構成しないtripletも生成し得ます。
その後の"the set of unique line segments composing all the triangle"の意味もよく分かりませんでした。想定解から察するに、三角形の各角を構成する2辺を、重複を排除してprintしているようです。
75. How to compute averages using a sliding window over an array?
あるarrayが与えられたとき、arrayの上で固定ウィンドウを動かし、各位置でのウィンドウ内の数値の平均値を求める問題です。
公式回答は以下になりますが、順を追わないとなかなか理解が難しいです。
def moving_average(a, n=3) :
ret = np.cumsum(a, dtype=float)
ret[n:] = ret[n:] - ret[:-n]
return ret[n - 1:] / n
仮に、a = [A_0, A_1, A_2, ..., A_(K-1)]とします。
np.cumsum(a)は累積和をとる関数で、ret = [S_0, S_1, ..., S_(k-1)]です。ここで、S_i = A_0 + A_1 + ... + A_iです。
ret[n:] - ret[:-n]と、ちょうど累積和の数列をn個ずらして差をとることで、数列Aの隣り合うn項の和を順に並べたものを生成できています。
より具体的に書くと、
ret[n:]は[S_n, S_(n+1), ..., S_(K-1)]、ret[:-n]は[S_0, S_1, ..., S_(K-n)]
なので、
ret[n:] - ret[:-n]は[S_n - S_0, S_(n+1) - S_1, ..., S_(K-1) - S_(K-n)]
です。
ただ、最後の計算で得られた配列には、その最初に来るべきA_0 + A_1 + ... + A_(n-1)が抜けています。A_0 + A_1 + ... + A_(n-1)はret[n-1]に等しいので、ret[n-1]に、最後の計算で得られた配列ret[n:] - ret[:-n]を連結すると、無事数列aの上で固定ウィンドウを動かしてウィンドウ内の和をとった数列が求まります。あとはnで割れば平均になります。
