python で行列ベクトル演算が可能です。でも、実際に行列ベクトル計算をしようとしたとき戸惑わされました。python での行列ベクトル演算について手頃な解説がありませんでした。コード例も殆どなく、試行錯誤で使う必要がありました。回り道をしました。特に Matrix と array の使い分けに戸惑いました。結論は「慣れるまでは Matrix を使わずに array の範囲だけで使っとけ。」です。慣れた後でも Matrix を使うメリットは限られます。array だけで済ましたほうが余分なことを考えずに済みます。
このような遠回りをすることなく python での数値計算を手っ取り早く始められるようにように、この Web page を書きました。C 言語や数値計算についての素養はあるが python は使い始めの方、早急に行列 ベクトル演算を行う必要がある方を対象に、python で行列・ベクトル演算をてっとり早く使いこなすための解説を行います。
python の線形演算モジュールには Numeric と numarray の二つが並存しています。両者の機能は殆ど同じです。numarray のほうが Numeric よりも新しいモジュールです。でも、入門的な使い方しかしないユーザーにとっては、大差ありません。どちらでも構いません。例えば Numeric では浮動小数点演算がオーバー・フローしたときなど NaN(Not a Number を意味します。) などを出力させていました。これを numarray では、例外を投げるかなど幾つかの選択をさせられるように進化させています。
新しい numarray を使えば良いとも限りません。vtk や mayavi などのモジュールが古い Numeric を使っていて、まだ新しい numarray に切り替わっていないためです。現在でも両方が並行して使われています。こちらから、どちらもダウンロードできます。
しかし直ぐにウンロードしないでください。、まだ python をインストールしていない方、まだ python を使いこんでない方は numarray や Numeric を直接インストールするのではなく、次に述べる Enthought パッケージを使い、他の応用モジュールも含め一度に纏めてインストールすることを勧めます。
python はパッケージとしても配布されています。python は多数の開発者によって多様な分野のライブラリが開発されています。多数のモジュールと共にソフト資産が蓄積され続けています。でも、この活発さは悪いことも引き起こします。多数の開発者が互いに関連し合うモジュールを独立して開発しているため、相互のモジュールがパージョンによって正常に動作する保証が取れなくなっています。そこで、それら複数のアプリケーション・モジュールの整合性を保ったパッケージとして python パッケージが作られ配布されています。
私自身は Enthought 社のパッケージを使っています。python が 2.3 の古いバージョンになってしまうのですが、次のアプリケーション群を一括してインストールできる手軽さを優先しています。
とくに科学技術分野の数値計算で ScyPy の部分が役立ちます。下の数値計算パッケージ纏めて入手できるのはありがたいはずです。
fft 特殊関数 Matlab 互換プロットパッケージ 統計処理パッケージ 積分 常微分方程式ソルバー Lapack 線形演算
Ehthought のパッケージが大きすぎると思う方は、こちらを使って SciPy パッケージだけをインストールすることもありだと思います。こちらから入手できる SciPy Tutrorial だけでも価値があります。python での科学技術計算が 42 page に要領よくまとめられています。
scipy の tutorial は良くできていますが、行列演算については、mat(Matrix に同じ) を中心に書いてあり惑わされました。mat や Matrix の array ラッパー・インターフェースを介さずに array のまま行列ベクトルを扱うべきです。
Enthought 社のパッケージで python をインストールすると、古い Nemeric がインストールされてしまいます。新しい numarray はインストールされません。
他の numarray を使っている python コードを使うためにも、numarray もインストールしておきましょう。こちらからダウンロードして追加インストールしましょう。Numeric と numarray の両方をインストールしても、混在して使わない限りは問題も起きないようです。
numarray/Numeric は array 型を使えるようにするモジュールです。python での線形演算のために tuple, list とは別に array データー型が numarray/Numeric で追加されました。numarray/Numeric では Fortran で開発された線形演算パッケージ Lapack を python に移植していますが、それは array 型のデータを対象にしています。
python での行列ベクトル計算のためには、tuple / array / Matrix / list と良く似たデータ型の違いを理解し使い分けることが必要となります。この四つのデータ型の関係を あえて図示すると次のようになります
tuple ⊂ ⊂ list
array
△
Matrix
list は tuple を包含します。array は tuple と list の中間にあります。中間と言っても「少しずれた中間」です。array は行列やベクトルという巨大な数値の集まりを効率的に処理するために設けられたデータ型です。C などの配列と同様に、数値データを連続したメモリブロックに配置します。tuple や list が参照を要素とするのに対し array は数値自体を要素とするので「少しずれた中間」と表現しました。
Matrx は array を継承したラッパー・クラスです。 *, ** 演算子がオーバーロードされています。行列どうしや行列とベクタの積を * 演算子で書けるようになります。行列のマイナスも含んだ べき乗を ** 演算子で記述できるようになります。
Enthought/SciPy パッケージでは、さらに Matrix を mat でラップします。mat と Matrix との違いは良く分かりません。mat コンストラクタは Matrix コンストラクタを呼び出しています。私は Matrix と殆ど同じだと推測しています。
でも、Matrix は array に追加する形で設けられたクラスです。Lapack の移植は Matrix ではなく array に対して行われています。この結果 Matrix は lapack の線形代数や array 関数との整合性が不完全です。Matrix インスタンスに lapack 関数を働かせると array に変わったりします。Matrix は + - * ** 演算子の範囲でのみの計算に限って使うもののようです。
最初は Matrix 型の使用を避けるべきです。Python での線形代数処理の理解を混乱させるだけです。Matrix 型の使用は python での線形代数処理の全体が見渡せるようになるまで控えるべきです。Matrix 型を使うにしても、その用途を限定して使うべきです。
array は多量の数値を扱うための型です。lapack はlist に移植することもできるとは思いますが、list は汎用的過ぎ、高速な数値計算に適しません。数値専用のコンテナ型 array を導入したほうが、数値計算プログラムを効率的に実装できます。
tuple,list,array の違いを、メモリに配置されているデータの視点から見てみましょう。
tuple は参照が並んだものです。tuple はメモリ内で参照ポインタを連続して配置することで実装されています。
┌──┬──┬──┬──┬ ┬──┬──┬──┬──┐ │参照│参照│参照│参照│ ・・・・・・│参照│参照│参照│参照│ └┼─┴┼─┴┼─┴┼─┴ ┴┼─┴┼─┴┼─┴┼─┘ ↓ ↓ ↓ ↓ ↓ ↓ ↓ ↓
list は、参照と list pointer のペアがリンクされたものです。tuple と同様に list 要素には異なった型のインスタンスを混在させられます。メモリ上にポインタを介して飛び飛びにリンク配置させるので、途中への挿入や抹消 sort などが可能です。でも tuple よりは複雑な処理が必要となります。遅くなります。
● が list pointer
┌──┬─┐ ┌──┬─┐ ┌──┬─┐ ┌──┬─┐
│参照│●┼→│参照│●┼→ ........... │参照│●┼→│参照│●┼→
└┼─┴─┘ └┼─┴─┘ └┼─┴─┘ └┼─┴─┘
↓ ↓ ↓ ↓
array はメモリ上に数値を直接並べたものです。処理は高速になりますが、整数/実数/複素数などサイズの異なる要素を混在させられません。
┌──┬──┬──┬──┬ ┬──┬──┬──┬──┐ │数値│数値│数値│数値│ ・・・・・・│数値│数値│数値│数値│ └──┴──┴──┴──┴ ┴──┴──┴──┴──┘
//@@
from numarray import *
#from Numeric import *
mtAt = array([ [1,2],
[3,5] ])
mtAt2 = array([[1,2],
[3,4, 5] ])
//@@@
# numarray のエラー
Traceback (most recent call last):
File "temp.py", line 5, in ?
mtAt2 = array([ [1,2],
File "C:\lng\Python\Lib\site-packages\numarray\numarraycore.py", line 395, in array
return fromlist(sequence,type,shape)
File "C:\lng\Python\Lib\site-packages\numarray\numarraycore.py", line 258, in fromlist
arr.fromlist(seq)
ValueError: Nested sequences with different lengths.
# Numeric のエラー
Traceback (most recent call last):
File "temp.py", line 5, in ?
mtAt2 = array([ [1,2],
TypeError: an integer is required
python 実行マクロお使いのエディタ・マクロで python コードを下のようにブロック実行させるマクロを作る事を奨めます。私自身は WZ エディタを使っており WZ の shell マクロと、それに行、ブロック実行マクロを組み合わせています。これについてもソースと一緒に公開してあります。宜しければ御利用ください。
この python 実行マクロは、 //@@ から //@@@ までの間にはさまれた文字列をカレント・ディレクトリの temp.py に書き出し「python temp.py」を実行させるものです。
//@@ ↑ │<== ctrl O + P 操作されたとき、この範囲をカレント ディレクトリの temp.py に ↓ 書き出し python temp.py を実行する //@@@以後 python code //@@ ... //@@@ で囲まれたコードは、このように実行される python コードと解釈してください。
numarray/Numeric には、下のように数多くの array を操作する関数が備わっています。でも行列・ベクトル計算をさせる分には ○ の付いた関数たちを使うだけで済みます。その他の関数は、python での線形代数計算に慣れてから再度調べれば充分でしょう。
Numeric.py にある最初のドキュメンテーション文字列より
Functions
○ array - NumPy Array construction
○ zeros - Return an array of all zeros
- shape - Return shape of sequence or array
- rank - Return number of dimensions
- size - Return number of elements in entire array or a
certain dimension
- fromstring - Construct array from (byte) string
- take - Select sub-arrays using sequence of indices
- put - Set sub-arrays using sequence of 1-D indices
- putmask - Set portion of arrays using a mask
○ reshape - Return array with new shape
- repeat - Repeat elements of array
- choose - Construct new array from indexed array tuple
- cross_correlate - Correlate two 1-d arrays
- searchsorted - Search for element in 1-d array
- sum - Total sum over a specified dimension
- average - Average, possibly weighted, over axis or array.
- cumsum - Cumulative sum over a specified dimension
- product - Total product over a specified dimension
- cumproduct - Cumulative product over a specified dimension
- alltrue - Logical and over an entire axis
- sometrue - Logical or over an entire axis
- allclose - Tests if sequences are essentially equal
More Functions:
- arrayrange (arange) - Return regularly spaced array
- asarray - Guarantee NumPy array
- sarray - Guarantee a NumPy array that keeps precision
- convolve - Convolve two 1-d arrays
- swapaxes - Exchange axes
- concatenate - Join arrays together
○ transpose - Permute axes
- sort - Sort elements of array
- argsort - Indices of sorted array
- argmax - Index of largest value
- argmin - Index of smallest value
- innerproduct - Innerproduct of two arrays
○ dot - Dot product (matrix multiplication)
- outerproduct - Outerproduct of two arrays
○ resize - Return array with arbitrary new shape
- indices - Tuple of indices
- fromfunction - Construct array from universal function
- diagonal - Return diagonal array
○ trace - Trace of array
- dump - Dump array to file object (pickle)
- dumps - Return pickled string representing data
- load - Return array stored in file object
- loads - Return array from pickled string
- ravel - Return array as 1-D
- nonzero - Indices of nonzero elements for 1-D array
- shape - Shape of array
- where - Construct array from binary result
- compress - Elements of array where condition is true
- clip - Clip array between two values
○ zeros - Array of all zeros
○ ones - Array of all ones
○ identity - 2-D identity array (matrix)
(Universal) Math Functions
add logical_or exp
subtract logical_xor log
multiply logical_not log10
divide maximum sin
divide_safe minimum sinh
conjugate bitwise_and sqrt
power bitwise_or tan
absolute bitwise_xor tanh
negative invert ceil
greater left_shift fabs
greater_equal right_shift floor
less arccos arctan2
less_equal arcsin fmod
equal arctan hypot
not_equal cos around
logical_and cosh sign
arccosh arcsinh arctanh
各関数の概略機能を知りたいとき、次のコードでドキュメンテーション文字列を表示させられます。
//@@
import Numeric as Nm
print dir(Nm)
help(Nm.zeros)
# other helps of commands
//@@@
['ArrayType', 'Character', 'Complex', 'Complex0', 'Complex16', 'Complex32', 'Complex64',
'Complex8', 'DumpArray', 'Float', 'Float0', 'Float16', 'Float32', 'Float64', 'Float8', 'Int',
'Int0', 'Int16', 'Int32', 'Int8', 'LittleEndian', 'LoadArray', 'NewAxis', 'Pickler', 'PrecisionError',
'PyObject', 'StringIO', 'UInt', 'UInt16', 'UInt32', 'UInt8', 'Unpickler', 'UnsignedInt16',
'UnsignedInt32', 'UnsignedInt8', 'UnsignedInteger', '__builtins__', '__doc__', '__file__',
'__name__', '__version__', '_numpy', 'absolute', 'add', 'allclose', 'alltrue', 'arange',
'arccos', 'arccosh', 'arcsin', 'arcsinh', 'arctan', 'arctan2', 'arctanh', 'argmax', 'argmin',
'argsort', 'around', 'array', 'array2string', 'array_constructor', 'array_repr', 'array_str',
'arrayrange', 'arraytype', 'asarray', 'average', 'bitwise_and', 'bitwise_or', 'bitwise_xor',
'ceil', 'choose', 'clip', 'compress', 'concatenate', 'conjugate', 'convolve', 'copy', 'copy_reg',
'cos', 'cosh', 'cross_correlate', 'cumproduct', 'cumsum', 'diagonal', 'divide', 'divide_safe',
'dot', 'dump', 'dumps', 'e', 'equal', 'exp', 'fabs', 'floor', 'floor_divide', 'fmod', 'fromfunction',
'fromstring', 'greater', 'greater_equal', 'hypot', 'identity', 'indices', 'innerproduct', 'invert',
'left_shift', 'less', 'less_equal', 'load', 'loads', 'log', 'log10', 'logical_and', 'logical_not',
'logical_or', 'logical_xor', 'math', 'matrixmultiply', 'maximum', 'minimum', 'multiarray', 'multiply',
'negative', 'nonzero', 'not_equal', 'ones', 'outerproduct', 'pi', 'pickle', 'pickle_array', 'power',
'product', 'put', 'putmask', 'rank', 'ravel', 'remainder', 'repeat', 'reshape', 'resize', 'right_shift',
'sarray', 'searchsorted', 'shape', 'sign', 'sin', 'sinh', 'size', 'sometrue', 'sort', 'sqrt', 'string',
'subtract', 'sum', 'swapaxes', 'take', 'tan', 'tanh', 'trace', 'transpose', 'true_divide', 'typecodes',
'types', 'vdot', 'where', 'zeros']
Help on built-in function zeros:
zeros(...)
zeros((d1,...,dn),typecode='l',savespace=0) will return a new array of shape (d1,...,dn) and type
typecode with all it's entries initialized to zero. If savespace is nonzero the array will be a
spacesaver array.
Fortran で実績のある線形演算パッケージ lapack が python に移植されています。次の行列関係の関数が用意されています
# numarray\linear_algebra\LinearAlgebra2.py または # Numeric\LinearAlgebra.py に実装されています。 def solve_linear_equations(a, b): def inverse(a): def cholesky_decomposition(a): def qr_decomposition(a, mode='full'): def eigenvalues(a): def eigenvectors(a): def singular_value_decomposition(a, full_matrices = 0): def generalized_inverse(a, rcond = 1.e-10): def determinant(a): def linear_least_squares(a, b, rcond=1.e-10):
これらの関数は array 型の行列に対して実装されています。array を継承した Matrix ではありません。Lapack の移植を行った方たちは Matrix モジュールの開発者とは別の方のようです。array や lapack の実装が終わった後に Matlab との互換を求めて Matrix wrapper を実装したように見えます。
現在の段階では、なまじ Matrix 型を使わずに array の範囲だけで計算プログラムを記述した方が簡単になります。行列の積を dot(,) で記述しようが * operator で記述しようが大差ありません。*, ** 演算子による見た目は数式に近くなりますが、そうすると Matrix と array 型が混在することになり色々と惑わされます。日常の数式に近い記述を望むなら、後で説明する sf を使うべきです。
tuple, list の数値一次元シーケンス・データから、一次元ベクトル：一次元 array を次のように生成できます。またrange(..) でリストを作ると同様に、start/end[/step] 引数を使って arange(..) でベクトル：一次元 array を作ることもできます。arange(..) は浮動小数点も引数にできます。(range(..) では int 引数に限定されます。)。でも arange でも複素数までは扱えません。複素数では大小関係が定義でないので、上限/下限に到達したか判定できないためでしょう。(sf ならば複素数まで扱えます。start, size, stride による slice array 形式のパラメータ設定の方が優れています。)
//@@ from Numeric import * # generage array from tuple, list print array([1,2,3]) print array( (4,5,6,7) ) print array(range(1,10)) # arange(..) function print arange(3,12) print arange(-0.4, 0.5, 0.1) #TypeError: can't convert complex to float; use abs(z) #print arange(-0.4+0j, 0.5+0j, 0.1j) //@@@ [1 2 3] [4 5 6 7] [1 2 3 4 5 6 7 8 9] [ 3 4 5 6 7 8 9 10 11] [-0.4 -0.3 -0.2 -0.1 0. 0.1 0.2 0.3 0.4]
array 要素が保持できる数値型が次のように定められています。
bool Int8, UInt8, # Uint8 は unsigned integer 8bit の意味です Int16, UInt16, Int32, UInt32, Int64, UInt64, Float32, Float64 Complex32, Complex64
これらの数値型は array 自体の型です。個別 array 要素の型ではありません。array の型は .type() で確認できます。tuple または list のシーケンス型から作られる array 型は、デフォルトでは Int32, Int32, Float64, Complex64 の何れかです。array(.) にシーケンス型データを与えて一次元ベクタを生成するとき,二番目の引数に array 型を指定できます。array 型を指定しないときは、デフォルト型と指定されている Int32 または Foat64, Complex64 型になります。
//@@ from numarray import * vctAt=array([1,2,3]) print type(vctAt[0]) print vctAt.type() vctAt=array([1.01,2,3]) print type(vctAt[0]) print vctAt.type() vctAt=array([1.01,2,3], Float32) print vctAt.type() vctAt=array([1+0.1j,2,3]) print vctAt.type() print arange(3,12).type() //@@@ < type 'int'> Int32 < type 'float'> Float64 Float32 Complex64 Int32
一次元 arrayについては list や tuple と同様に、[..] や ':' を使い、要素へのインデックスを指定することでアクセスします。[start: end: step] の三つのパラメータによる部分列の指定は重宝します。特に行列の斜め要素群の指定に便利です。慣れておくべきです。
//@@ # -*- encoding: cp932 -*- from Numeric import * vctAt = array(range(1,10)) print vctAt print vctAt[3] print vctAt[3:7] print vctAt[3:] print vctAt[:3] print vctAt[-2] # step:2 を指定した要素の取り出し print vctAt[3:7:2] vctAt[0]=100 print vctAt # ステップ step:2 を使った要素への代入 # Int32 型の array に浮動小数点を指定しても Int32 整数になってしまう vctAt[3:7:2] = [4.1, 8.2] print vctAt //@@@ [1 2 3 4 5 6 7 8 9] 4 [4 5 6 7] [4 5 6 7 8 9] [1 2 3] 8 [4 6] [100 2 3 4 5 6 7 8 9] [100 2 3 4 5 8 7 8 9]
サイズが同じベクタ:一次元 array どおしでは + - * / 演算子を使った四則演算が可能です。
//@@ from scipy import * vctLeftAt = array(range(0,5),Float64) vctRightAt = array(range(5,10)) print vctLeftAt print vctRightAt print vctLeftAt + vctRightAt print vctLeftAt - vctRightAt print vctLeftAt * vctRightAt print vctLeftAt / vctRightAt //@@@ [ 0. 1. 2. 3. 4.] [5 6 7 8 9] [ 5. 7. 9. 11. 13.] [-5. -5. -5. -5. -5.] [ 0. 6. 14. 24. 36.] [ 0. 0.16666667 0.28571429 0.375 0.44444444]
サイズの揃わないベクタ：一次元 array どおしの加減乗除算はエラーになります。
//@@ from numarray import * # ValueError: Arrays have incompatible shapes print array([1,2,3]) + array([1,2]) //@@@
scipy モジュールEnthought または scipy パッケージをインストールしているときは、Numeric の変わりに scipy モジュールを import するだけで、Numeric 機能を使えるようになります。scipy モジュールは array 以外に fft や常微分方程式ソルバーなど、パッケージが用意している数値演算一般を纏めて import／使えるようにします。これも重宝します。
また ベクタ:一次元 array と scalar の + - * / 四則演算が可能です。
//@@ from Numeric import * vctLeftAt = array(range(0,5)) print vctLeftAt * 3.3 print vctLeftAt / 3.3 print vctLeftAt + 3.3 print vctLeftAt - 3.3 //@@@ [ 0. 3.3 6.6 9.9 13.2] [ 0. 0.3030303 0.60606061 0.90909091 1.21212121] [ 3.3 4.3 5.3 6.3 7.3] [-3.3 -2.3 -1.3 -0.3 0.7]
最後のベクタ：一次元 array と scalar との加減算は数学では定義されていませんが、array と scalar の + - 演算として実装されています。でも、このような数学では使わない一次元 array と scalar の +,- 演算は避けるべきです。(array での行列とベクタとの積では、もっと多くの数学的には定義されていない演算が実装されています。)
//@@ from numarray import * print dot(array([1,2,3]),array([4,5,6]) ) //@@@ 32
複素数ベクトルのときは dot 関数は内積ではありません。 要素毎の掛け算と足し合せにすぎないからです。conjugate を取ってからの掛け算と足し合わせではありません。複素数ベクトルの内積を計算させるには自分で conjugate(.) 関数を働かせます。
//@@ from Numeric import * print dot(array([1+2j,2,3]),array([1+2j,5,6]) ) print dot(array([1+2j,2,3]),conjugate(array([1+2j,5,6])) ) //@@@ (25+4j) (33+0j)
Numeric モジュールには innerproduct(.) 関数も有ります。でも残念ながら dot(.) と同じでした。複素ベクトルの内積を計算させるには conjugate() が必要でした。
array は二次元以上の行列も扱えます。matrixInstance[i,j] のような行列形式でのインデックス指定が可能です。二次元 array インスタンスに対して lapack 関数を働かせられます。
array 配列の : によるアクセスも多すぎるほどに様々の使い方が追加されています。でもとりあえずは「 行列[index, :] 」で index 行のベクトル array を、「 行列[:, index] 」で index 列のベクトル array を表すことを理解しておけば、大体の行列処理は可能です。
tuple, list の二次元シーケンスから行列：二次元 array を次のように生成します。
//@@
from numarray import *
mtAt = array([[1,2,3],
[4,5,6],
[7,8,10]])
print mtAt
print mtAt.type()
//@@@
[[ 1 2 3]
[ 4 5 6]
[ 7 8 10]]
Int32
Numeric の二次元 array では .type() が使えません。
//@@
from Numeric import *
mtAt = array([[1,2,3],
[4,5,6],
[7,8,10]])
print mtAt
# AttributeError: type
print mtAt.type()
//@@@
行列：二次元 array では全ての行、列のサイズが同じでなければなりません。そうでないとエラーになります。
//@@
from numarray import *
# ValueError: Nested sequences with different lengths.
mtAt = array([[1,2,3],
[4,5],
[7,8,10]])
//@@@
reshape(.) を使うことで、一次元ベクタ array から二次元行列 array に変形することも可能です。
//@@ # -*- encoding: cp932 -*- from numarray import * # [1,2,..,10] のベクトルを作り 2x5 行列に reshape する print reshape( arange(10), (2,5) ) //@@@ [[0 1 2 3 4] [5 6 7 8 9]]
zeros(,) によって要素全てが 0 の行列を、ones(,) によって、要素全てが 1 の行列を生成できます。 identity(.) によって、指定サイズの単位行列を作れます。
//@@ # -*- encoding: cp932 -*- from scipy import * # 2x5 の全要素が 0 の整数行列を作ります print zeros( (2,5) ) # 3x5 の全要素が 1 の整数行列を作ります print ones( (3,5) ) # 5x5 の浮動小数点の単位行列を作ります print identity(5, Float64 ) //@@@ [[0 0 0 0 0] [0 0 0 0 0]] [[1 1 1 1 1] [1 1 1 1 1] [1 1 1 1 1]] [[ 1. 0. 0. 0. 0.] [ 0. 1. 0. 0. 0.] [ 0. 0. 1. 0. 0.] [ 0. 0. 0. 1. 0.] [ 0. 0. 0. 0. 1.]]
二次元 array については matrixInstance[i,j] のように、二つの index で、行列要素を指定できます。matrixInstance[i,:] により i 行要素を、 matrixInstance[:,j] により j 列要素を指定できます。matrixInstance[start:end:step] で部分要素を指定できます。
//@@ # -*- encoding: cp932 -*- from Numeric import * mtAt = reshape( arange(4*4), (4,4) ) print mtAt mtAt[0,0] = 100 print mtAt print mtAt[0,:] print mtAt[:,1] #行列の対角要素部分を指定します print mtAt[0:5:6] //@@@ [[ 0 1 2 3] [ 4 5 6 7] [ 8 9 10 11] [12 13 14 15]] [[100 1 2 3] [ 4 5 6 7] [ 8 9 10 11] [ 12 13 14 15]] [100 1 2 3] [ 1 5 9 13] [ [100 1 2 3]]
サイズの同じ行列どおしの四則演算は要素毎の演算を作用させるだけです。
//@@ # -*- encoding: cp932 -*- from Numeric import * mtLeftAt = reshape( arange(3*3), (3,3) ) mtRightAt = reshape( arange(1.0, 3*3+1), (3,3), ) print mtLeftAt print mtRightAt print mtLeftAt + mtRightAt print mtLeftAt - mtRightAt print mtLeftAt * mtRightAt print mtLeftAt / mtRightAt //@@@ [[0 1 2] [3 4 5] [6 7 8]] [[ 1. 2. 3.] [ 4. 5. 6.] [ 7. 8. 9.]] [[ 1. 3. 5.] [ 7. 9. 11.] [ 13. 15. 17.]] [[-1. -1. -1.] [-1. -1. -1.] [-1. -1. -1.]] [[ 0. 2. 6.] [ 12. 20. 30.] [ 42. 56. 72.]] [[ 0. 0.5 0.66666667] [ 0.75 0.8 0.83333333] [ 0.85714286 0.875 0.88888889]]
行列：二次元 array と scalar との乗除算は数学と同じ結果になります。でも scalar との加減算は行列要素全部に対して同じ scalar 量を足したり引いたりします。エラーにはならないのでご注意ください。
//@@ # -*- encoding: cp932 -*- from numarray import * mtAt = reshape( arange(3*3), (3,3) ) print mtAt*3.4 print mtAt/3.4 print mtAt+3.4 print mtAt-3.4 //@@@ [[ 0. 3.4 6.8] [ 10.2 13.6 17. ] [ 20.4 23.8 27.2]] [[ 0. 0.29411765 0.58823529] [ 0.88235294 1.17647059 1.47058824] [ 1.76470588 2.05882353 2.35294118]] [[ 3.4 4.4 5.4] [ 6.4 7.4 8.4] [ 9.4 10.4 11.4]] [[-3.4 -2.4 -1.4] [-0.4 0.6 1.6] [ 2.6 3.6 4.6]]
ベクトル一次元 array との積は判りにくい結果となります。数学における行列とベクタの積にはなりません。使わないほうが無難です。数学の意味での行列とベクタの積を計算させるには次に述べる dot(,) 関数を使います。
//@@ # -*- encoding: cp932 -*- from scipy import * mtLeftAt = reshape( arange(3*3), (3,3) ) mtRightAt = reshape( arange(1.0, 3*3+1), (3,3), ) print dot( mtLeftAt, mtRightAt) print dot( mtLeftAt, [1,2,3]) print dot( (4,5,6), [1,2,3]) //@@@ [[ 18. 21. 24.] [ 54. 66. 78.] [ 90. 111. 132.]] [ 8 26 44] 32
行列とベクタの dot 演算ではベクトルを縦行列に直さなくても、正しく計算してくれます。後で述べる Matrix を使うと縦行列に修正してやらないと計算してくれません。array での行列とベクトルの dot 演算を使っておいたほうが記述が簡単です。
Fortran で開発された lapack が phtyhon の linear algebra モジュールに移植されています。ここにある関数を array 行列インスタンスに適用します。まず一番使う頻度の高い、逆行列を求める inverse(.) 関数から始めます。
//@@ # -*- encoding: cp932 -*- from numarray import * mtAt = reshape( arange(3*3), (3,3) ) # 逆行列が存在する値にする mtAt[2,2] = 10 from numarray.linear_algebra import * print inverse(mtAt) # 確認計算 print dot(inverse(mtAt), mtAt) //@@@ [[-0.83333333 -0.66666667 0.5 ] [ 0. 2. -1. ] [ 0.5 -1. 0.5 ]] [[ 1. 0. 0.] [ 0. 1. 0.] [ 0. 0. 1.]]
linear algebra モジュールdot などの lapack による線形代数演算を行うときiには、
import します。
- numarray では numarray.linear_algebra モジュールを
- Numeric や scipy では LinearAlgebra モジュールを
逆行列が存在しない行列に inverse(.) を働かせるとエラーになってくれます。
//@@ # -*- encoding: cp932 -*- from Numeric import * # mtAt の逆行列が存在しません。 mtAt = reshape( arange(3*3), (3,3) ) from LinearAlgebra import * # LinearAlgebra.LinAlgError: Singular matrix print inverse(mtAt) //@@@
その他の lapack 関数を実行させます。
//@@
# -*- encoding: cp932 -*-
from numarray import *
mtAt = reshape( arange(3*3), (3,3) )
# 逆行列が存在する値にする
mtAt[2,2] = 10
from numarray.linear_algebra import *
#from LinearAlgebra import *
#print cholesky_decomposition(mtAt)
print qr_decomposition(mtAt)
print eigenvalues(mtAt)
print eigenvectors(mtAt)
print singular_value_decomposition(mtAt)
print generalized_inverse(mtAt)
print determinant(mtAt)
//@@@
(array([[ 0. , 0.91287093, 0.40824829],
[-0.4472136 , 0.36514837, -0.81649658],
[-0.89442719, -0.18257419, 0.40824829]]), array([[ -6.70820393, -8.04984472, -11.18033989],
[ 0. , 1.09544512, 1.82574186],
[ 0. , 0. , 0.81649658]]))
[ 14.6544476 -1.04590831 0.39146071]
(array([ 14.6544476 , -1.04590831, 0.39146071]), array([[ 0.15082592, 0.45445878, 0.87790589],
[ 0.90425716, -0.13682593, -0.40447207],
[ 0.23673828, -0.84967597, 0.47117485]]))
(array([[-0.12899165, 0.9914205 , -0.02113178],
[-0.45702173, -0.07834688, -0.88599825],
[-0.88005244, -0.1046287 , 0.46320681]]), array([ 15.45352414, 1.02187585, 0.37994924]), array([[-0.43041184, -0.5252812 , -0.73404721],
[-0.84434213, -0.05320404, 0.53315635],
[ 0.31911128, -0.8492638 , 0.4206174 ]]))
[[ -8.33333333e-01 -6.66666667e-01 5.00000000e-01]
[ 1.20736754e-15 2.00000000e+00 -1.00000000e+00]
[ 5.00000000e-01 -1.00000000e+00 5.00000000e-01]]
-6.0
//@@
from numarray import *
from numarray.matrix import *
mtAt = Matrix([[1,2,3],
[4,5,6],
[7,8,10]])
print mtAt*mtAt
print mtAt**2
print mtAt**-1
print mtAt**-1 * mtAt
//@@@
[[ 30 36 45]
[ 66 81 102]
[109 134 169]]
[[ 30 36 45]
[ 66 81 102]
[109 134 169]]
[[-0.66666667 -1.33333333 1. ]
[-0.66666667 3.66666667 -2. ]
[ 1. -2. 1. ]]
[[ 1.00000000e+00 0.00000000e+00 0.00000000e+00]
[ 0.00000000e+00 1.00000000e+00 -3.55271368e-15]
[ 0.00000000e+00 0.00000000e+00 1.00000000e+00]]
でも逆行列が存在しないはずなのに計算してくれます
//@@
from numarray import *
from numarray.matrix import *
mtAt = Matrix([[1,2,3],
[4,5,6],
[7,8,9]])
print mtAt*mtAt
print mtAt**2
print mtAt**-1
print mtAt**-1 * mtAt
//@@@
[[ 30 36 42]
[ 66 81 96]
[102 126 150]]
[[ 30 36 42]
[ 66 81 96]
[102 126 150]]
[[ -4.50359963e+15 9.00719925e+15 -4.50359963e+15]
[ 9.00719925e+15 -1.80143985e+16 9.00719925e+15]
[ -4.50359963e+15 9.00719925e+15 -4.50359963e+15]]
[[ 0. -4. -8.]
[-8. 0. 0.]
[ 4. 0. 0.]]
Matrix * array を計算させると Matrix になってしまいます。 dot(array, array) は array の範囲で閉じており無用な型の変換が発生しません。dot(.) の方が使いやすく感じます。
Matrix モジュールは、array を継承して Matrix クラスを作って無理に * ** 演算子を実装しただけに見えます。Matrix の完成度は低く感じます。
実際、全世界の web page でも Matrix については殆ど解説されていません。Matrix を使うと計算速度も何倍も遅くなります。Numerx\Matrix.py を見ても、実装してあるのは演算子のオーバーロードだけです。
さらには、上にあるような逆行列が存在しないはずなのに計算してしまう不都合が入り込んでいます。
私自身は Matrix を使わないようにしています。どうしても最初は Matrix を使ったほうが便利にみえます。でも、これで計算しようして、いろいろとトラブルに遭遇しました。始めは Matrix を使わないことを推奨します。Python での線形代数処理全体が見渡せるようになってから用途を限定して Matrix 型を使うべきだと主張します。そのほうが効率的に python での線形代数計算を使いこなせるようになります。
sf は行列/ベクトル計算に特化した CUI アプリケーションです。 Windows コンソールからも使えますが、エディタに組み込むことで、より効果的に使えます。私は Wz マクロ に組み込んでいます。ctrl O + N または ctrl O + S 操作により、エディタ画面のカーソル位置にある sf式 または sf ブロック式を実行させています。
その他に次のような特徴を備えています
# ε0 は真空の誘電率です。`oku は 1 億の意味です。 # 1`coulomb の電荷から 1`meter の距離の電位は 89.9 億ボルトです。 q@=1`coulomb, (q/(4`π ε0 1`meter))/`oku < 89.8755 >
mtRes = `σx `σy^-1 + `σz^3 (2+i) < 2+2i, 0 > < 0, -2-2i >
type mtRes.val % BaseComplexDoubleMatrix 2, 2 < 2 +2i, 0 +0i > < 0 +0i, -2 -2i >
mtRes^-3 < -0.03125-0.03125i, 0 > < 0, 0.03125+0.03125i >
# パイソンでのリストの生成 //@@ start=0; end=10; step=2; print range(start,end,step) //@@@ [0, 2, 4, 6, 8]
# 上と同じ操作を sf のベクトル生成で行います。 start=0, end=10, step=2, << start,(end-start)/step,step >> < 0, 2, 4, 6, 8 > # 同じことを、より sf らしく記述します start=0, size=5, step=2, << start,size,step >> < 0, 2, 4, 6, 8 >
# パイソンでの内包リスト操作 //@@ print [x**2 + 2*x + 3 for x in range(0,10,2)] //@@@ [3, 11, 27, 51, 83]
# 上と同じ内包リスト操作を、次の sf ベクトル操作で行います。 <<0,5,2 @x| x^2 + 2x + 3>> < 3, 11, 27, 51, 83 >
//@@
<<0,3,1 @j|\
<<0,4,3 @k| j+k^2>>
>>
//@@@
< 0, 9, 36, 81 >
< 1, 10, 37, 82 >
< 2, 11, 38, 83 >
sf /?
・
・
basic functions
!sin !cos !tan !sinh !cosh !tanh !exp !log !asin !acos !atan !log10
other functions
!abs: absolute for each element
!norm: square root norm
!fft: Fast Fourier Transformation
!rft: Reverse Fast Fourier Transformation
!dggr: dagger for matrix, conjugate for vector and scalar
!image, !real: get image/real from complex number
!size: get vector size
!print: print argment string into left side value file and on console
!sqrt : square root
!sum : sum up all vector element
!prdct: product of all vector element
!print(....): print #g .... to file variable
jacobi.exe 対称行列/エルミート行列の対角化と対角化させる直交化行列を求めます。 expM.exe 行列の exponetial 指数関数を求めます random.exe ランダム データを保持するベクタを生成します。 shift.exe ベクタを shift します。次のような使い方ができます。
# 5x5 のゼロ正方行列 mtSample を作る mtSample =[[5]] < 0, 0, 0, 0, 0 > < 0, 0, 0, 0, 0 > < 0, 0, 0, 0, 0 > < 0, 0, 0, 0, 0 > < 0, 0, 0, 0, 0 > # mtSample の対角要素にランダム値を seed 0,最大値 100 で設定する mtSample[*,*] = ~random(5,0,100) < 3, 0, 0, 0, 0 > < 0, 42, 0, 0, 0 > < 0, 0, 40, 0, 0 > < 0, 0, 0, 69, 0 > < 0, 0, 0, 0, 29 > # 対角要素の一つ上の斜めにランダム値を seed 1,最大値 100 で設定する mtSample[<1,4,6>] = (~shift(mtSample[*,*],-2))[<0,4,1>] < 3, 40, 0, 0, 0 > < 0, 42, 69, 0, 0 > < 0, 0, 40, 29, 0 > < 0, 0, 0, 69, 0 > < 0, 0, 0, 0, 29 > # 対角要素の一つ下のの斜めに上、対角要素を 2 つぶんシフトさせたものを設定する mtSample[<5,4,6>] = mtSample[<1,4,6>] < 3, 40, 0, 0, 0 > < 40, 42, 69, 0, 0 > < 0, 69, 40, 29, 0 > < 0, 0, 29, 69, 0 > < 0, 0, 0, 0, 29 > # mtSample の t = 0.01 秒あとの指数行列値を計算し mtExp に設定する t@=0.01, mtExp = ~expM(t mtSample) < 1.13029, 0.563328, 0.195688, 0.0204324, 0 > < 0.563328, 2.0171, 1.16757, 0.175587, 0 > < 0.195688, 1.16757, 1.94361, 0.54758, 0 > < 0.0204324, 0.175587, 0.54758, 2.07341, 0 > < 0, 0, 0, 0, 1.33643 > # mtExp の対角成分：固有値を求める ~Jacobi(mtExp) # mtExp を対角化させる直交行列 < 0.506821, -0.802171, 0, -0.232383, 0.213662 > < -0.647721, -0.126369, 0, -0.381128, 0.647479 > < 0.552174, 0.524912, 0, -0.0122326, 0.647627 > < -0.136723, -0.254988, 0, 0.894757, 0.340144 > < 0, 0, 1, 0, 0 > # mtExp の固有値 < 0.618039, 1.09748, 1.33643, 1.98583, 3.46307 >
sf を使えば python の 1/2 から 1/3 程度で行列 ベクタ計算処理を記述できます。sf でプロトタイプを作って動作を確認してから、python に落として詳細・高速計算をさせる方法が使えます。実際に独楽の運動のシミュレーションで、 sf によりプロトタイプを作り、python で詳細・高速計算させています。
gdsp.bat gspl.bat により sf ベクタ変数／行列変数を gnuplot にグラフ表示させられます。このバッチ プログラムは sf ファイル変数を gnuplot で表示できるフォーマットに変換し gnuplot.ini に書き出してやり、gnuplot.exe を起動します。gdsp.bat, gdsp.bat の中身は次のようなものです。
type gdsp.bat gnpltMDt.exe %1 > gnuplot.ini start wgnuplot.exe type gspl.bat gnSplt.exe %1 > gnuplot.ini start wgnuplot.exe
gnpltMDt.exe, gnSplt.exe が sf ファイル変数を gnuplot 向けのデータに変換するフィルター プログラムです。これらはこちら にある sf12?zip:sf の配布 zip ファイルに含まれています。
例えば [0,2π] での sin(1+x^2) の関数グラフを次のように表示させられます
N@=128, <<0,N,2`π/N @x| !sin(1+x^2)>>
gdsp
.gif)
N@=128, temp=<<0,N,2`π/N @x| !sin(1+x^2)>>
gdsp temp
gnuplot などを起動する手間を含めて考えると、sf と gdsp.bat によることが、関数グラフを表示させる一番単純な方法だと主張します。
極端なまでに単純化しているため、上の gnuplot のグラフ表示で x 軸のレンジは [0,140] になっています。[0, 2π] になっていません。でもユーザーは横軸の範囲が [0,2π] であることは解っています。[0,140] で良いと主張します。!print(表示制御コマンド文字列) sf コマンドにより、gnuplot に表示制御コマンド文字列を引き渡せます。gnuplot グラフの横軸を [0,2π] にできます。でも通常は、それだけの手間をかける意味がありません。
scipy と pylab モジュールを使えば python で次のように、matlab に近い形でグラフを表示させられます。( Enthought のパッケージでインストールしていれば pylab もインストールされています。)
//@@ from scipy import * from pylab import * _N = 128 _arSin = [sin(1+x**2) for x in arange(0,2*pi, 2*pi/_N)] plot(arange(_N)*2*pi/_N, _arSin) # pair connect with line show() //@@@
.gif)
gspl.bat は三次元メッシュ表示をさせるときに使います。下に例を示します。
//@@
/s
N@=30
_temp@=[[N]]
<<0,N,2`π/N @x,j|\
<<0,N,2`π/N @y,k|\
_temp[j,k]=!sin(x^2+y^2)
>>
>>
temp = _temp
//@@@
gspl temp
.gif)
gdsp, gspl のより詳細な説明については、こちらを参照ください。
sf と python をミックスして使えます。sf のユーザー関数として python コードで書かれた関数を呼び出せます。逆に python コードで計算した行列／ベクトル データを sf ファイル変数に変換し、sf で ちまちまと操作・処理する方法も使えます。このような sf と python の共用を可能にするのが sfVal.py モジュールです。
sfVal.py モジュールは python の array クラスで記述した行列／ベクトル変数と sf ファイル変数との間で交互に変換することを可能にします。
sfVal の getSf() を使えば sf ファイル変数を次のように python の array 変数に変換できます。sf での計算結果を python の array 変数に移せます。
N@=128, << 0,N,2`π/N @x|!exp(!sin(x^2))>>
//@@
import scipy
import pylab
import sfVal
_N = 128
_ar = sfVal.getSf('_dt')
# pair connect with line
pylab.plot(scipy.arange(_N)*2*scipy.pi/_N, _ar)
pylab.show()
//@@@
逆に sfVal の putSf() を使えば python の array 変数を次のように sf ファイル変数に変換できます。python での計算結果を sf でも使えます。
# 複素平面での Γ 関数の複素数値分布を行列データにする
//@@
import scipy
import sfVal
_N = 128
_ar = scipy.array(
[[ scipy.special.gamma(x+y*1j) for x in scipy.arange(-1+0.01,1,0.02)]
for y in scipy.arange(1-0.01,-1,-0.02)
]
)
sfVal.putSf(_ar,'rs')
//@@@
# 実数軸に近い位置で [-1,1]の範囲のデータを積分してみる !sum(rs[49,*])0.01 < -2.27873-0.713216i > # trace をとってみる !sum(rs[*,*]) < -29.9937+24.9733i >
sf には Γ 関数などの特殊関数がありません。でも python の Γ 関数を、下の python ルーチンを介して Γ 関数を呼び出せます。_arg{1}.val ファイル変数に呼び出し時の引数が入ること、_rs.val に結果を入れて返す規則を守って python コードを書くだけです。
# _arg[1}.val 変数で与えられたスカラー値に対応する Γ 関数値を返す。
//@@
import scipy
import sfVal
_ar = sfVal.getSf('_arg{1}')
sfVal.putSf(scipy.array([scipy.special.gamma(_ar[0])]),'_rs')
//@@@
//copy \#####.### gamma.py /y
~gamma.py(2) < 1 > ~gamma.py(4) < 6 > !exp(~gamma.py(2+3i)^3) < 1.00152+0.00109786i >
gamma.py() 引数には複素数も可能です。分野によっては、この gamma.py() sf ユーザー関数は重宝するでしょう。
sf からの python コードの呼び出しは特殊関数に限りません。常微分方程式の求解など使い切れないほどの数値演算処理パッケージが python モジュールとして公開されています。ユーザーの必要に応じて、getSf(),putSf() 介して sf からも呼び出してください。。
sfVal.putSf(), sfVal.getsf() の機能はカレント ディレクトリへの sf テキスト ファイル形式でのデータの書き出しと読み出しです。シリアライズを行っています。 array 専用の pickles であるとも言えます。
sfVal.putSf(), sfVal.getsf() ならば、ファイルを open(.) する必要も有りません。記述が簡単になります。シリアライズした結果できるカレント ディレクトリの sf ファイル変数データは sf により、グラフ表示させるなど、そのままでも数学的な処理ができます。より便利に使いまわせます。
array の pickles としては sfVal.py のほうが重宝するはずです。この意味でも putSf()/getSf() を御利用ください。
小規模な計算は sf のほうが手っ取り早く行えます。~gamma.py(3.4i) 4/3 などと数式に近い記述ができます。でも規模が大きくなり、沢山のループ処理を行うようになると、CPU の計算時間が何十分とかかるようになります。そのときは python の利用を考えるべきです。python の方が二桁近く高速です。
python ならば、数値演算のためのモジュールが既に存在しています。線形演算モジュールや特殊関数モジュール、積分／微分方程式ソルバーモジュール、プロット モジュールなどが無償で配布されています。これらを使えば C プログラムのときより簡単に、教科書に出てくる数式の具体例を作れます。でもプログラム作業を必要とします。
逆に多く使われるが単純な加減乗除などの処理は、プログラムの手間がない sf で行わせた方が簡単です。何度も使われる python 計算は sf から利用できる python コードで書かれたユーザー関数にしてやれば、sf 計算式に python 計算を組み込んでしまえます。
適宜、sf, python の処理を組み合わせて両方を活用ください。
負整数の割り算による丸め方が c 言語とは 0.5 ぶんだけずれています。下のコードを参照ください。
//@@ print -1/2 0xff print -2/2 0xfe print -3/2 0xfd print -4/2 0xfc //@@@ -1 -1 -2 -2
//05.07.22 test -1/2 #includeusing namespace std; void main (void) { cout << -1/2 << endl; cout << -2/2 << endl; cout << -3/2 << endl; cout << -4/2 << endl; } //@@@ >copy c:\#####.### a.cpp /y >cl a.cpp /MLd /W3 /Od /I.\vrf /D"_CONSOLE" /YX /GX /GR /Fp#.pch netapi32.lib 0 -1 -1 -2
バイナリ ビットのデータ配列整数の割り算回路で、MSB 符号を溜まったままに右シフトによって 1/2 演算を行わせることを考えると pyhon での結果も納得できます。
c でも C99 が決まるまでは負整数の割り算が決まっていなかったそうです。python が負の整数の割り算で上のようになるのは、歴史的経緯のためです。