Numpy¶
Numpyとは¶
Pythonは、科学技術計算及び機械学習の分野で広く使われていますが、その理由のひとつは使いやすいパッケージが整備されているからです。 ここでは、その中でも最も重要なNumpyの使用法を簡単に説明します。
物理では、ベクトルや行列が用いられますが、その計算は、基本的には繰り返しです。 例えば、C言語では、ベクトル $(1.0,2.0,3.0)$の長さの2乗は
double vec[3];
vec[0]=1.0; vec[1]=2.0; vec[2]=3.0;
double sum = 0;
for(int i=0; i<3; i++){
sum += vec[i]*vec[i];
}
printf("%f\n",sum);
と書くことで計算ができます。
これをそのままPythonで書くと、
vec = [1.0, 2.0, 3.0]
sum=0.0
for i in range(3):
sum += vec[i]**2
print(sum)
となります。
この書き方でも問題はありませんが、Pythonの短所のひとつは繰り返し文が非常に遅いことです。 そのため、for文などはできるだけ避けることが、高速計算へと繋がります。
上記でやっていることは、単純に配列の各要素を2乗し、得られた配列の要素の和をとれば良いだけです。 Numpyを使えば、その処理は1行で以下のように書けます。
import numpy as np
vec = np.array([1.0, 2.0, 3.0])
sum = np.sum(vec**2)
print(sum)
配列はnp.array()関数で用意し、vec**2で各要素を2乗し、np.sum()という関数で和をとることになります。
それぞれの要素を2乗したり、和をとったりするところのnp.sum()関数の中ではfor文を回しているから同じでは?と思われますが、この部分は、別個にC言語などで書かれた高速なアプリケーションを用いて実行されていますので、C言語と遜色ない速度が出ます。
高速で見やすいプログラムを書くためには、ブロードキャスト及びスライスは特に重要です。
詳細は、公式ページ を参考にしてください。
Numpyを使うには¶
まず、Numpyを使うには、そのパッケージがインストールされている必要があります。 Anacondaを使っているのであれば、すでに導入されているはずです。
プログラムで使用するには、コードの最初に
import numpy as np
と書いてください。以下では、特に断りがない場合は、上記が書かれているとして進めます。
配列の作成¶
基本的な生成方法¶
まずは、配列をどのように作るかを説明します。 最も簡単なのは、リストから作成する方法で、
vec = np.array([1.0, 2.0, 3.0])
と書きます。
len(vec)とすれば、配列の長さがわかります。
リストの使用法も参考にしてください。
Numpyの配列(np.ndarray)では、特にsizeを用いて要素数を取り出せます。
N = vec.size
とすると、配列のサイズがNに入ります。
2次元配列(行列)の場合は、
mat = np.array([[1.0, 2.0, 3.0], [4.0, 5.0, 6.0]])
です。この場合、2x3の行列になります。
要素へのアクセスは、例えば、1行3列目はmat[0,2]となります(0始まりなので)。
行列のサイズの取得方法ですが、mat.sizeとすると6が返ってきます。
行の数と列の数を別個にほしければ、
M, N = mat.shape
とすれば、行数がMに列数がNに入ります。
注釈
C言語での最も致命的なエラーは配列の数とfor文などで配列の要素にアクセスする数が合っていない場合です。
定義した配列の長さを超えてアクセスしようとすると、コンパイルが通ったとしても変な挙動を示します。
それは、自分で定義した要素数と、実際に定義されている配列の数が合っていないからです。
Pythonでは、配列を定義した後で、 vec.size や mat.shape などを用いるか、for文をできるだけ使わないようにして、これを回避するようにしてください。
ある範囲を等間隔に分割¶
np.linspace関数を使います。
例えば、0から10までを11個の要素に等分割(10分割)する場合には、
np.linspace(0, 10, 11)
と書きます(答えは[0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])。
np.linspace(0, 10, 10)
とすると、[ 0., 1.11111111, 2.22222222, 3.33333333, 4.44444444, 5.55555556, 6.66666667, 7.77777778, 8.88888889, 10.]のように配列が構築されますので注意してください。
注釈
最後の要素を入れたくない場合もあると思いますが、そのときは、 np.linspace(0, 10, 10, endpoint=False) とすればよいです。
詳しい使い方は公式ページを参考にしてください。
等差数列¶
配列を等差数列的に生成したい場合もあります。
このときは、np.arange関数を使います。
例えば、0以上1未満の範囲で0.1ずつ増加するような要素を持つ配列は
np.arange(0, 1, 0.1)
として作成できます。
詳しい使い方は公式ページを参考にしてください。
空の配列¶
空の配列は、np.empty関数で作成できます。C言語で配列を定義することに相当します。
C言語で、
double vec[10];
に相当するのは、Pythonでは
vec = np.empty(10)
となります。 この場合、中身に何が入っているかはわかりませんので、必ずプログラムのどこかで値を入れる必要があります。
Numpyは暗黙的に倍精度としています。もし、例えば、複素数の配列にしたければ、
vec = np.empty(10, dtype='complex')
としてください。
行列の場合は、例えば、2x3のときは、
mat = np.empty((2, 3))
とします。
詳しい使い方は公式ページを参考にしてください。
すべて特定の数字を入れた配列¶
すべての要素に特定の数字を入れた配列を作る場合は、np.full()関数を使います。
値を1.0を持つ要素を10個持つ配列は
np.full(10, 1.0)
で作成できます。
詳しい使い方は公式ページを参考にしてください。
注釈
すべての要素が0の配列は np.zeros() 関数、1の配列は np.ones() 関数で作成できます。デフォルトは float 型(Cでいうdouble型)で作成します。
np.zeros(10) は np.full(10, 0.0) と、np.ones(10) は np.full(10, 1.0) と同じです。
ブロードキャスト¶
Numpyでは配列に対する要素の演算を容易に記述できます。
配列と数の演算¶
配列と普通の数(スカラー)同士の演算は、配列のすべての要素に作用します。 例えば、
vec = np.array([1.0, 2.0, 3.0])
vec + 2.0
とすると、[3., 4., 5.]となります。
配列同士の演算¶
配列同士の演算は、その長さが合致している場合のみ計算ができます。
vec1 = np.array([1.0, 2.0, 3.0])
vec2 = np.array([-1.0, 4.0, 0.0])
vec1 * vec2
とすると、配列の要素同士のかけ算になるので、結果は[-1., 8., 0.]となります。
注釈
ベクトルの内積は、 np.dot を使います(複素数の場合は np.vdot )。例えば、 np.dot(vec1, vec2) のように使います。これは np.sum(vec1 * vec2) と同じです。
配列の合成¶
2つの1次元配列をくっつけることができます。
1次元配列を縦に繋げるには
vec1 = np.array([1.0, 2.0, 3.0])
vec2 = np.array([-1.0, 4.0, 0.0])
np.r_[vec1,vec2]
とすると、[ 1., 2., 3., -1., 4., 0.]が得られます。
また、
vec1 = np.array([1.0, 2.0, 3.0])
vec2 = np.array([-1.0, 4.0, 0.0])
np.c_[vec1,vec2]
とすると、3x2の2次元配列になります。
[[ 1., -1.], [ 2., 4.], [ 3., 0.]]
スライス¶
スライスの基礎¶
スライスを用いると、配列の一部を取り出したり、書き換えたりすることが容易になります。
ある配列に対して、vec[2:]とすると、2番目(0番を最初として)からのものになります。
例えば、
vec = np.array([1.0, 5.0, 3.0, -1.0, 4.1, 9.6])
vec[2:]
とすると、[ 3. , -1. , 4.1, 9.6]になります。
vec[2:5]とすると、2番目(0番を最初として)から5番目未満のものを取り出します。
vec = np.array([1.0, 5.0, 3.0, -1.0, 4.1, 9.6])
vec[2:5]
とすると、[ 3. , -1. , 4.1]になります。
最初を省略しvec[:5]とすると最初(0番目)から5番目未満のものを取り出します(vec[0:5]と同じ)。
さらにvec[:5:2]とすれば、最初(0番目)から5番目未満のもののうち、2個おき(1個飛ばし)で取り出します。
vec = np.array([1.0, 5.0, 3.0, -1.0, 4.1, 9.6])
vec[:5:2]
とすると、[1. , 3. , 4.1]となります。
スライスを用いた書き換え¶
スライスを用いて配列の中身を書き換えることもできます。 例えば最初の3つだけ書き換えたい場合、
vec = np.array([1.0, 5.0, 3.0, -1.0, 4.1, 9.6])
vec[:3] = [0,1,2]
とすれば良いです。
vec = np.array([1.0, 5.0, 3.0, -1.0, 4.1, 9.6])
vec[:3] = -vec[:3]
とすれば、最初の3つの要素の符号が反転します。
多次元配列におけるスライス¶
行列の一部だけ取り出したり、書き換えたいときがよくあります。 例えば行列の0行目(一番最初の行)だけ取り出したいときは、
A=np.array([[1., 2., 3.], [4., 5., 6.]])
A[0,:]
とします。一次元配列[1., 2., 3.]が得られます。列に関しても同様です。
A=np.array([[1., 2., 3.], [4., 5., 6.]])
A[:2,:2]
とすれば、左側の2x2の部分を取り出せます。
Numpyの便利な関数¶
よく使う簡単な関数¶
np.sum(): 配列の要素の和をとります。np.average(): 配列の要素の平均をとります。np.diff(): 近接する配列の要素の間の差の配列を求めます (オプションを指定しなれば、(もとの要素数-1)の配列になります)。np.max()、np.min(): 配列の要素の最大値、および最小値を求めます。np.argmax()、np.argmin(): 配列の要素の最大値、および最小値となるインデックス(何番目の要素か)を求めます。
数学定数¶
np.pi: 円周率np.e: ネイピア数np.euler_gamma: オイラーの定数np.inf: 無限大np.nan: 非数値
数学関数¶
Numpyの数学関数はブロードキャストが有効です。もちろん配列でなく普通の数に対しても使えます。
例えば、np.sin(vec)とすれば、各要素に対してsin関数を適用した配列が得られます。
vec = np.pi * np.array([0.0, 0.25, 0.5, 0.75, 1.0])
np.sin(vec)
は、[0.0 0.7071, 1.0, 0.7071, 0.0]となります。
np.sqrt(): 平方根np.exp(): 指数関数np.log(): 自然対数np.sin()、np.cos()、np.tan(): 三角関数np.arcsin()、np.arccos()、np.arctan(): 三角関数の逆関数np.sinh()、np.cosh()、np.tanh(): 双曲線関数np.real()、np.imag()、np.conj()、np.abs()、np.arg(): 複素数の実部、虚部、複素共役、絶対値、偏角
など。
複素数に対しても有効ですが、多価関数となる場合は値域に注意してください。
ファイルの読み書き¶
ファイルの読み書きもNumpyを使うことで容易にできるようになります。
ファイルの読み込み¶
例えば、以下のようなdata.txtを読み込む場合を考えます。
1.0 2.0
2.0 4.0
3.0 6.0
4.0 8.0
以下からダウンロードできます。
1列目をxs[n]、2行目をys[n]という配列に入れる場合、C言語では例えば以下のように書きます。
読み込むまで配列の大きさがわからないので、とりあえず100個からなる配列を用意しています。
配列の大きさはnに入るようになっています。
#include <stdio.h>
#include <stdlib.h>
int main(){
double x, y , xs[100], ys[100];
int n=0;
FILE *fp;
fp = fopen("data.txt", "r"); // ファイルをオープン
while (fscanf(fp, "%lf %lf", &x, &y) != EOF){ // ファイルの最後まで1行ずつの読み込みを繰り返す
xs[n]=x;
ys[n]=y;
n++;
}
fclose(fp); //ファイルを閉じる
return 0;
}
これはpythonでは、
xs, ys = np.loadtxt('data.txt', unpack = True)
の1行でできます。
注釈
unpack = True オプションをつけない場合、 np.loadtxt は対応する2次元配列(上の data.txt の場合は、4x2の行列)として返します。
そのため、 xs, ys = np.loadtxt('data.txt', unpack = True) は、
data = np.loadtxt('data.txt')
xs = data[:, 0]
ys = data[:, 1]
と同じです。
ファイルの書き込み¶
配列xs[n]とys[n]が与えられたときに、それをファイルに書き出す方法を説明します。
C言語では、
#include <stdio.h>
#include <stdlib.h>
int main(){
int n=4;
double xs[4] = {1.0, 2.0, 3.0, 4.0};
double ys[4] = {2.0, 4.0, 6.0, 8.0};
FILE *fp;
fp = fopen("out.txt", "w"); // ファイルをオープン
for(i = 0; i < n; i++){
fprintf(fp, "%lf %lf\n", xs[i], ys[i]);
}
fclose(fp); //ファイルを閉じる
return 0;
}
となりますが、Pythonだと、np.savetxt()関数を使います。
xs = np.array([1.0, 2.0, 3.0, 4.0]);
ys = np.array([2.0, 4.0, 6.0, 8.0]);
np.savetxt('out.txt', np.c_[xs, ys])
とシンプルに書けます。
np.c_[xs, ys]は、配列の合成にあるとおり、xsとysを結合し、4x2の2次元配列にします。
この2次元配列をファイルに出力しています。
行列計算¶
Numpyを使うと物理でよく出てくる行列計算が簡単に記述できます。
行列積¶
例えば、 $ A= \begin{pmatrix} 1&2&3\cr 4&5&6 \end{pmatrix} $ と $ B= \begin{pmatrix} 1&-2\cr 5&-3\cr -2&0 \end{pmatrix} $ の間の行列積を計算します。当然、Aの列の数とBの行の数が合致していなければいけません。
python3.5(2015年リリース)以降では、行列の積を計算する@が導入されています。
A=np.array([[1., 2., 3.], [4., 5., 6.]])
B=np.array([[1., -2.], [5., -3.], [-2., 0.]])
A @ B
答えは、[[ 5., -8.], [ 17., -23.]]となります。
ちなみにベクトル(一次元配列)に使うと内積になります。
転置、エルミート共役¶
転置はT、エルミート共役は転置してからnp.conj()で複素共役を取れば良いです。
A=np.array([[1., 2., 3.], [4., 5., 6.]])
A.T
とすると、[[1., 4.], [2., 5.], [3., 6.]]となり、
A=np.array([[1.-0.5j, 2.+1.j], [4.+3.j, 5.+7.j]])
np.conj(A.T)
とすると[[1.+0.5j, 4.-3.j ], [2.-1.j , 5.-7.j ]]が得られます。
逆行列¶
逆行列はnp.linalg.inv()を使います。
A = np.array([[1.-0.5j, 2.+1.j], [4.+3.j, 5.+7.j]])
np.linalg.inv(A)
とすると、
[[-0.49411765+1.22352941j, -0.03529412-0.34117647j], [ 0.05882353-0.76470588j, 0.14705882+0.08823529j]]
が得られます。
実際、
A = np.array([[1.-0.5j, 2.+1.j], [4.+3.j, 5.+7.j]])
B = np.linalg.inv(A)
A @ B
とすれば、単位行列が得られることが確認できます。
対角化¶
固有値と固有ベクトルを求める対角化には
np.linalg.eig()を用います。すなわち、
$$
P^{-1} A P=D=
\begin{pmatrix}
\lambda_1 & & \cr
& \lambda_2 & & \cr
& &\lambda_3 & \cr
& & & \ddots \cr
\end{pmatrix}
$$
の固有値$\{\lambda_1,\lambda_2,\cdots \}$と変換行列$P=(\bm{v}_1,\bm{v}_2, \cdots)$を求めます。
A = np.array([[1.,-0.2, 2.], [4., 2., -1.], [-3., -1., 2.]])
w, P = np.linalg.eig(A)
P @ np.diag(w) @ np.linalg.inv(P)
この一番下の式は、$PDP^{-1}$を計算しています。np.diag(w)は配列wを対角に持ち他の要素がゼロの行列で、$D$に対応します。
最後の式が、Aと合致することを確認してみてください。
特に対象とする行列が対称行列もしくはエルミート行列の場合は、np.linalg.eig()を用います。
詳しい使い方は公式ページ(eig)と公式ページ(eigh)を参考にしてください。
また、scipyという科学技術計算用のパッケージを使っても同様のことができます。
from scipy import linalg
A = np.array([[1.,2, 1.], [2., 1., 2.], [1., 2., 1.]])
w, U = linalg.eigh(A)
U @ np.diag(w) @ U.T
最初のfrom scipy import linalgで線形代数に関係する関数を使えるようにしています。ここでは、
対称行列
$
A=\begin{pmatrix}
1&2&1 \cr
2&1&2 \cr
1&2&1
\end{pmatrix}
$を直交行列$U$を用いて$U^T A U=D$のように対角化しています。
最後の行のU @ np.diag(w) @ U.Tで$U D U^T$を計算しているので、これがAと一致することを確認してください。
詳しい使い方は公式ページ(eig)と公式ページ(eigh)を参考にしてください。
3重対角行列の対角化¶
3重対角行列の対角化はnumpyではできないのでscipyを使います。
ここでは特に対称な3重対角行列を考えます。
このとき行列は、
$$
A=
\begin{pmatrix}
a_1 & b_1 & & \cr
b_1& a_2 & b_2 & & \cr
& b_2& \ddots &\ddots & \cr
& & \ddots& \ddots & b_{N-1} \cr
& & & b_{N-1} &a_{N} \cr
\end{pmatrix}
$$
となります(例によって空白の行列要素はゼロです)。
対角部分$\{a_i\}$とそこからひとつずれた部分$\{b_i\}$のみ有限で、それ以外の要素はゼロです。
$N\times N$の行列だと、$\{a_i\}$は$N$個、$\{b_i\}$は$N-1$個です。
以下は、 $ A= \begin{pmatrix} 1 & -0.5 & 0 \cr -0.5& 2 & -0.8 \cr 0& -0.8& 3 \end{pmatrix} $ の行列を対角化します。
from scipy import linalg
a = [1,2,3]
b = [-0.5,-0.8]
w, U = linalg.eigh_tridiagonal(a, b)
U @ np.diag(w) @ U.T
最後の行の答えがもとの行列と合致していることを確認してください。
慣れてきたら¶
Pythonでやりたいことがあれば、まずはGoogleなどで検索しながら調べるのをおすすめします。
たとえば、numpy 配列 和をとるなどで検索すれば、多くのwebページが見つかると思います。
但し、内容は玉石混交です。
Pythonに慣れてきたら、Numpyの公式のページを参照した方が、最新の情報と正しい関数の使い方がわかります。