複天一流:どんな手を使ってでも問題解決を図るブログ

宮本武蔵の五輪書の教えに従い、どんな手を使ってでも問題解決を図るブログです(特に、科学、数学、工学の問題についてですが)

pythonで線形代数☜超初歩入門:(1)行列とベクトルの定義

Python線形代数の計算をやる

関数の呼び出しでベンチマークコンパイルオプションや計算速度測定)やった。次に、乱数生成プログラミングを行い、「64ビット限界計算」をやってみた。

そして、いよいよ線形代数に踏み出す時が来た。物理エンジンに基づく現代の3Dゲームや、写真や動画など様々な画像処理で必須の計算内容である。CPUやGPU線形代数の演算に特化した命令や機能をたくさん盛り込んでいるから、ベンチマークでCPUやGPUの本当の性能を測りたいなら、線形代数の計算を利用するのは(もはや)デフォルトであろう。

まずはpythonで勉強してみる。pythonによるコード書きは「スケッチ」である、と以前も書いたが、学習、教育用にはもってこいの言語である。

「超初歩入門」

「初歩から学ぶ」、つまり素人が初めてのものに取り組むときは、その筋の入門書を買うのが、従来の「最初のステップ」だった。豊富な経験を持つ「プロ」である著者が易しい例を用いながら噛み砕いて説明してくれるわけである。それを模倣しながら、「あっ、わかったよ」となるまで読み続ける作業に従事する(いわゆる"I know Kung-Fu"というやつである)。

youtu.be

しかし、「あっ、わかったよ」となるまでに、余りにも長い時間がかかると、途中で嫌になり挫折ということになるのが大抵であろう(笑)。多くの凡人(私も含む)にとっては、「入門」と「挫折」は密接な関連をもつ「セット(対の組)」である。

今回、Python線形代数の計算をするにあたり、夏休み直前に大学生協にて購入したPython入門教科書(教養学部テキスト、という触れ込み)をめくってみた(前期試験の監督業務が終了した直後に「(7-11)良い気分」で衝動買いした本である)。が、線形代数に関する記述は皆無であった(最小二乗法や連立方程式を解く、というセクションもあったが、ブラックボックス的なライブラリ関数で作業しており、線形代数の計算は表に浮上してこない)。検索して見つけた「Pythonプログラミング入門」(東大情報センター,2020-2024)というデジタル教科書には、”§17.9線形代数の演算(p.158)”という箇所があったので、期待して閲覧してみたが、内積単位行列の設定程度の簡単な内容しか書かれておらず、深く落胆した。

ネット検索をすると、いろいろなブログ記事がヒットするが、個々のコマンドの細かい解説ばかりになりがちで、”コンテクスト”を重視した「流れるような読みやすい文章」からほど遠く、しかも「実用的」な例が少なくて読みにくい傾向が強い。人間が書くネット記事は、なぜか網羅的になりがちなのである。そして、これが実に「非人間的」で疲れてしまう。

さらに、中には記述が間違っているものもあり、「初心者」には判断が難しく、疲労を加速してくれる。たとえば、こちらのブログに書かれていた「対角化」というのは、いわゆる対角化ではなく、単に正方行列の対角成分の行列要素を抜き出す、という意味であった(笑)。

以前は、それでも情報量の少ないネット上のダメ記事も含め、検索三昧して読みまくり、自分の頭の中で情報を取捨選択して、しだいに「あっ、わかったよ」のレベルまで頑張ったものである。しかし、そこに至るまでには紆余曲折を余儀なくされ、非常に体力精神力を消耗したものである。人間は「繰り返しに弱く」、そして「忍耐力がない」存在なのである。

この2点に強い存在がある。それが「計算機」である。ここ1、2年で、計算機が”コミュ力”を手に入れたため(つまりはAI)、人間と計算機の間のインターフェースは非常にスムースになった。したがって、人間が不得意な「繰り返し」をAIに「忍耐強く」やってもらい、その上澄みを(七里ヶ浜のトンビよろしく)「かっさらう」のが、これからの人間の勉強法である。これを「超初歩入門」という言葉で表現したいと思うのである。つまり、今回はAIに教えてもらった非常に「人間的な」線形代数演算に関するPythonプログラミングの入門デアル。

行列とベクトルの定義

まずは、Pythonで行列をどう扱うか、であるが、CやFORTRANであれば配列(たいては二次元配列だが)なのであるから、それに類似したリスト(やその類似物)が候補となるはずである。

CやFORTRANのように、自分で線形代数演算のアルゴリズムをすべてコードにすれば、それでいいだろうが、ここはPythonである。豊富なライブラリ関数があるはずなのである。特にNumPyやSciPyといった、科学技術計算関連のライブラリに、いいのが定義されているはずである。これを利用しない手はない。

ということで、検索してネット記事を読んだ結果として、mat()とarray()の2つが候補として出てきたのである。NumPy、SciPyどちらでも利用できる。これらの関数(メソッド)は、リストとして最初に定義されたオブジェクトを、”行列クラス”に変換する機能をもつ。行列クラスに変換することで、逆行列や対角化など線形代数の演算を表現したメソッドで処理できるというわけだ。

最初のコードはしたがって、ここから始めよう。

import numpy as np

NDIM=2
a = [ 10*(i+1)+j for i in range(NDIM)] for j in range(NDIM)]
print(a)

▷ [[10, 20], [11, 21]]

リストの形式は、すでに「行列」っぽいわけだが、この状態では線形演算子は利用できないから、「行列」に変換する。 すでに書いた通り、2種類のメソッドがある。

A=np.mat(a)
B=np.array(a)
print(A)
print(B)

▷
[[10 20]
[11 21]]

[[10 20]
[11 21]]

物理学では「単位が異なるもの」は足し算できない。たとえば、1(m)+10(kg)というのは不可能な演算(和)である(経済学ではときどきやるらしいが...)。ということで、両者を足してみて答えが出れば、「同類」と判断できるだろう。

print(A+B)

▷
[[20 40]
 [22 42]]

おーっ!ちゃんと足された。ということは、両者ともに「同じもの」つまり「行列」として扱っていいみたいである。

ちなみに、リストaと行列Aを足してみよう。どうなるだろうか?(異なるクラスに属するオブジェクト同士なので、エラーとなるはずだが....)

print(a+A)

▷
[[20 40]
 [22 42]]

足せた?(どうなってんだー!) こちらはどうか?

print(a+B)

▷
[[20 40]
 [22 42]]

こちらも足せたよ!(困惑)

じゃあ、次は掛けてみよう。行列の積を表す演算メソッドはいくつかあるらしいが、最近最もよく利用されているのが「@」だという。

print(A@B)

▷
[[320 620]
 [341 661]]

print(A@a)

▷
[[320 620]
 [341 661]]

print(B@a)

▷
[[320 620]
 [341 661]]

なんだよ!全部できるんじゃないかっ!(ちょっとした怒)

ということで、行列に関する基本演算、和と積に関しては、リスト、mat(), array()はどれも区別なく混用できるようである(笑)。

ちなみに、*記号(アスタリスク)はどうだろうか?

print(A*B)
print(A*a)
print(B*a)

▷
[[320 620]
 [341 661]]

[[320 620]
 [341 661]]

[[100 400]
 [121 441]]

ここで、やっと差が出た!B*aは内積みたいな順番で掛け算だけがされており、足し算がどこにも使われてない。一方の、A*BやA*aでは積和の演算の組み合わせで行列の積が(通常の定義で)実施されている。

ちなみにリスト同士だとどうなるのだろうか?

print(a+a)
print(a*a)
print(a@a)

▷
[[10, 20], [11, 21], [10, 20], [11, 21]]

 print(a*a)
          ~^~
TypeError: can't multiply sequence by non-int of type 'list'

2行目でエラーが出たが、1行目はリストの「拡大」という解釈となっていて、算術演算にはなっていない点が面白い。 3行目に辿り着けなかったので、実行し直すと、

print(a@a)

▷
   print(a@a)
          ~^~
TypeError: unsupported operand type(s) for @: 'list' and 'list'

これもエラーとなった。やはり、リストは行列ではないと考えた方がいいだろう(ふーっ..一安心)。

行列とリストの混合演算では、たぶんインタープリターがキャストをかけてクラス変換を自動的に行った、つまり、上位にある「行列クラス」にリストを格上げしてくれたんだろう。インタープリターの「頭の良し悪し」に依存しそうな演算なので、やはりリストは行列とは混ぜない方がよさそうだ。

一方で、mat()とarray()の間の関係については、まだよくわからない。ここでAIにご登場を願う。質問してみると、次のような回答であった。

いい質問です 👍
numpy.mat() と numpy.array() は似ているようで、けっこう大きな違いがあります。

この柔らかい回答で、精神的に持っていかれる人が多発しているようである(笑)。

脱線だが、AIは絶対に「あなたを否定しない」といわれる。とくにChatGPT ver.4oがその傾向が強かったらしく、自殺願望のある少年がAIとのやりとりの結果、「自殺を勧められた」と誤解してしまい、実際に自殺してしまったという事件が最近報道された

検索していたら、最近ChatGPTは「焚き火」という特別セッションを開設したようである。

chatgpt.com

話を戻そう。AIからの回答である。

1.  numpy.array()

基本的な N 次元配列 (ndarray) を作る関数。

多次元(1D, 2D, 3D…)自由自在。

✳︎演算子は要素ごとの掛け算(アダマール積)。

行列積をしたいときは @ 演算子か np.dot() を使う。

2. numpy.mat()

必ず 2 次元行列 (numpy.matrix クラス) を作る。

✳︎演算子が「行列積」として動く。

逆行列や転置など、行列計算を意識したメソッドが多い。

なるほど。でも、どっちにすべきかは、これだけではわからない。回答はまだ続いている。

AIの回答の続き

ほおー。「array()を使うべし」という「(弱めの)アドバイス」である。これを信じるなら、「mat()には未来がないらしいから、今のうちにarray()に慣れておこうか」と思ってしまう人は多だろう。実際、検索したネット記事でもarray()が多かったのは、このせいだろう。

AIの結論
....だそうである。

ここで、次の質問を投げかけてみた。

aが2x2のリストだとする。

A=np.mat(a)
B=np.array(b)
C=A+B

Cは「行列」だとは思うが、それはmat()なのか、それともarray()なのか?

これに対する回答は

いい視点ですね!

NumPy では、ndarray と matrix を混ぜた演算をすると、自動的に matrix に「引き上げ」られます。

ほーお? AIの回答は続く。

AIによる注意点

混乱している人は世界にたくさんいる(「いた」?)ようだ。おそらく、AIに類似の質問を投げかけた人数は相当数に及ぶのであろう。すでに答えを用意していたような書き振りである。

この後、いくつかの行列演算(行列式逆行列など)について質問した後に、ベクトルの定義についてもAIに聞いてみた。array()で定義する流派とmat()で定義する流派があるようで、その中身は結構複雑であった(驚)。ということで、ベクトルについては次回に回そうと思う。

今日の最後は、mat()に定義されているという逆行列の演算メソッドを試してみたい。A.Iと書くそうである。これが理由で、AIはmat()を推奨してもいいんじゃないかと思った次第である(笑)。

A = np.mat([[1,2],[3,4]])
AI = A.I
print(A, AI, A*AI)

▷ 
[[1 2]
 [3 4]] 

[[-2.   1. ]
 [ 1.5 -0.5]] 

[[1.0000000e+00 0.0000000e+00]
 [8.8817842e-16 1.0000000e+00]]

もし、array()を使うとなると、次のようになる。

A = np.arange(1,5).reshape(2,2)
AI = np.linalg.inv(A)
print(A,AI, A@AI)

▷
[[1 2]
 [3 4]] 

[[-2.   1. ]
 [ 1.5 -0.5]] 

[[1.0000000e+00 0.0000000e+00]
 [8.8817842e-16 1.0000000e+00]]

inv()までの記述が長い...。また、掛け算が@になっている(*を使うと、先ほどの「要素同士の積」(アダマール積?)になってしまう)。

mat()か、array()か?

個人的にはarray()の方が「用語」「記号」として好きである。mat()と書くと前転とかバック転を想起してしまう(あるいは新マンの部隊)。が、しかし、A.Iという表現も捨て難い!

個人的には、昔からプログラムで行列を表すときはmtrx()を使用してきたので、できればこの表現を利用したいと思う。ということで、次のような関数を定義することにした。

import numpy as np

def mtrx(A):
  id = 1
  if id == 0:
      return np.array(A)
  else:
      return np.mat(A)

デフォルトではid=1だから、A.Iが使えるmat()表現。しかし、id=0に変更すれば、array()にスイッチできる。メインのコードではmtrx()を利用すればいいので、将来的にmat()が廃止になっても、この関数の定義を修正するだけでメインコードは無修正で済む。 (sedを使って書き換えればそれで済むという話もあるが...私はmtrx()が好きなだけである!)