Python♪FEM:行列の積や転置行列のコード記述のコツ

行列の積や転置行列についてNumpyの関数を使えば簡単ですが、関数を使わない方法も紹介します。PythonのFEMの書籍はまだ少ないため、CやFORTRANのコードを読むときに役に立ちます。熟語を覚えるような感覚でコードを覚えれば簡単です。

1.行列の積

X=ABという行列の積を考えてみます。例としてAは3行3列、Bは3行2列です。行数、列数は1から始まりますが、Pythonのfor文が0から始まることを考慮して、添え字を00から始めています。

慣れない方は以下のような図で計算をイメージするとわかりやすいです。図のポイントは以下の通りです。

  • 行列の計算では下図のkaとkbの大きさは同じです。赤や青の点線で囲った部分同士の内積を求めるのですから、同じ長さでなければ計算できません。
  • 図下側の計算式をを参考に、AやBのどの位置のベクトルを使うのか確認してください。
  • 計算結果である行列Xの行数、列数はそれぞれi行、j列となります。

コード01は、X= ABの行列の積の計算例です。11行目の.dot()を使えば行列計算も1行で終了です。Pythonはfor文の処理が遅いので、行列処理はできるだけnumpyの関数を使うべきです。

#コード01
import numpy as np

a = np.array([[0., 1., 2.],
              [3., 4., 5.],
              [6., 7., 8.]])
b =  np.array([[0., 1.],
               [2., 3.],
               [4., 5.]])

x = np.dot(a, b)
print(x)
#出力01
[[10. 13.]
 [28. 40.]
 [46. 67.]]

コード02はnumpyの関数を使わない方法です。15~18行目がx = np.dot(a, b)に相当する部分です。この部分のコードは熟語を覚えるようなつもりでコードを覚えておいても損はありません。Pythonでは使うことはないかもしれませんが、他のプログラム言語ででこの形を見たときに、行列の積であるとわかれば読みやすくなります。

#コード02
import numpy as np

a = np.array([[0., 1., 2.],
              [3., 4., 5.],
              [6., 7., 8.]])
b =  np.array([[0., 1.],
               [2., 3.],
               [4., 5.]])
x = np.zeros((3, 2))
ii = 3
kk = 3
jj = 2

for i in range(ii):
    for j in range(jj):
        for k in range(kk):
            x[i][j] = x[i][j] + a[i][k] * b[k][j]
print(x)
#出力02
[[10. 13.]
 [28. 40.]
 [46. 67.]]

2.行列、ベクトルの積と内積

コード03は1行3列の行列と、3行1列の行列の積です。出力03のように両側を[ ]で囲まれた[14.]が出力されます。

#コード03
import numpy as np

a = np.array([0., 1., 2.])
b =  np.array([[3.],
               [4.],
               [5.]])

x = np.dot(a, b)
print(x)
#出力03
[14.]

だたし、numpyのdot関数では、コード04のようにベクトルの内積を求めることもできます。行列計算であれば、aの行数とbの列数が異なるので計算できないことになりますが、エラーにはならず内積が求められます。dot関数では要素数が同じ1次元配列同士の計算では行列の積ではなく内積が計算されます。出力04のように両側に[ ]はなく、14.0という値が出力されます。

#コード04
import numpy as np

a = np.array([0., 1., 2.])
b = np.array([3., 4., 5.])

x = np.dot(a, b)
print(x)
#出力04
14.0

なお、書籍やサイトで「行列の積」を「行列の内積」と表現しているケースが散見されますので注意してください。内積は対応する成分の積の全ての和のことであり1つの値となります。また、式も積はX=AB、内積はX=A・Bです。

例えば「ゼロから作るでDeep Learning」でも初版第8刷で「3.3.2 行列の内積」となっていますが、これは「内積」ではなく「積」です。

3.転置行列

転置行列は行と列を入れ替えた行列であり、下図のようになります。

転置行列は行と列を入れ替えた行列であり、下図のようになります。行列A, Bの転置行列がそれぞれAT, BTです。赤色の点線でひっくりかえしたような行列になります。

では、実際に下図の計算を行い、行列XをPythonで求めてみましょう。

コード05はnumpyの.Tメソッドを使った計算方法です。10行目でcの転置行列としてc.Tを用いています。

#コード05
import numpy as np

a = np.array([[0., 1., 2.],
              [3., 4., 5.],
              [6., 7., 8.]])
c = np.array([[0., 2., 4.],
              [1., 3., 5.]])

x = np.dot(a, c.T)
print(x)
#出力05
[[10. 13.]
 [28. 40.]
 [46. 67.]]

コード06はnumpyのメソッドを使わない方法です。17行目をよく見ると、コード02とは異なりc[k][j]ではなく、 c[j][k]となっていることに注目してください。転置行列は行と列の入れ替えなので、jとkを入れ替えているのです。

#コード06
import numpy as np

a = np.array([[0., 1., 2.],
              [3., 4., 5.],
              [6., 7., 8.]])
c =  np.array([[0., 2., 4.],
               [1., 3., 5.]])
x = np.zeros((3, 2))
ii = 3
kk = 3
jj = 2

for i in range(ii):
    for j in range(jj):
        for k in range(kk):
            x[i][j] = x[i][j] + a[i][k] * c[j][k]
print(x)
#出力06
[[10. 13.]
 [28. 40.]
 [46. 67.]]

線形代数の知識としては、「積」「転置行列」「内積」の他にも「逆行列」「単位行列」ぐらいは、その概要を理解しておきましょう。

私が実際に購入した教材のご紹介

以下、私が実際に購入したPythonの教材をまとめてみました。 Pythonを学習する上で、少しでもお役に立つことができればうれしいです。

Python♪私が購入したPythonの書籍のレビュー

UdemyのPythonの動画講座を書籍を買う感覚で購入してみた