Part 4: より高度な数学

行列と線形代数

行列は、階数２の配列によって表現されます。一方、ベクトル（幾何ベクトルではなく、線形代数におけるより一般的な意味でのベクトル）は、階数１のベクトルで表現されます。 Numeric モジュールには、行列やベクトルを扱うのに実用的な演算がいくつかあります。

dot(a, b) は、a と b のドット積を返します。これは、二つのベクトルのスカラー積、行列とベクトルの積、二つの行列の積、に使えます。 a の最後の次元と b の最初の次元の長さは等しくなくてはなりません。
transpose(a) は、行列 a の転置行列を返します。実際には、この関数は任意の階数の配列に適用でき、一般的には次元の順序を反転させます。
identity(n) は、n 行 n 列の単位行列を返します。
ravel(a) は、配列 a の全ての要素からなるベクトルを返します。要素の順序は、最初の添字が最も緩やかに変化するように決められます。すなわち、２行２列の配列要素は、(0,0), (0,1), (1,0), (1,1) の順序で返されます。
reshape(a, shape) は、配列 a の要素からなり、（タプルで）指定された形状を持つような配列を返します。（配列要素の順序は ravel(a) と同じです。）これは、特殊な行列を作る際にしばしば便利です。この制限は resize(a, shape) 関数では任意の形状へと緩和されています。a が持つよりも多くの要素が必要な場合には、最初の要素に戻って繰り返されます。

LinearAlgebra モジュールには、様々な慣用の行列演算があります。行列は二次元の配列で表されます。正方行列でなければならないために入力のチェックがなされるような演算もあります。全ての演算は、実行列にも複素行列にも使えます。これらの関数は、十分にテストされ効率の良い LAPACK ルーチンを基礎に作られています。

solve_linear_equations(a, b) 関数は、正方正則行列 a と右側ベクトル b からなる線形方程式系の解きます。 b を二次元配列（すなわちベクトルの列）とすることによって、複数の右ベクトルを同時に扱うことができます。 inverse(a) 関数は、適切な b を用いて solve_linear_equations(a, b) を呼び出すことにより、正方正則行列 a の逆行列を計算します。

eigenvalues(a) 関数は、正方行列 a の固有値を返します。 eigenvectors(a) 関数は、固有値と固有ベクトルを返し、後者は二次元配列（すなわち、ベクトルの列）です。

determinant(a) 関数は、正方行列 a の行列式を返します。

singular_value_decomposition(a) 関数は、 V, S, W^T の三つの配列を返します。これらの行列積は元の行列 a です。 V と W^T はユニタリ行列（階数２の配列）であり、一方 S は特異値行列の対角要素からなるベクトル（階数１の配列）です。この関数は、行列が悪条件にないかどうか（あるとしたらどのようにか）をチェックするのに主に使われます。

generalized_inverse(a) 関数は、一般化逆行列（偽逆行列またはムーア・ペンローズ逆行列としても知られています）を返します。これには、線形方程式や最小二乗法の問題に関連した応用が沢山あります。

linear_least_squares(a, b) 関数は、データ過剰な線形方程式系の最小二乗解を返します。第三の引数はオプションで、特異値の範囲を決めるカットオフを与えます（デフォルトは 10^-10）。返り値は四つあります：降順で、最小二乗解そのもの、残渣の二乗和（すなわち、解が最小にする量）、行列 a の階数、そして a の特異値です。

例：

import Numeric; N = Numeric
import LinearAlgebra; LA = LinearAlgebra

# 行列の生成
a = N.reshape(N.arrayrange(25.), (5, 5)) + N.identity(5)

# 逆行列の計算
inv_a = LA.inverse(a)

# a * a^{-1} - identity(5) の最大絶対値の要素を印字し、逆行列を検証
print "Inversion error:", \
      N.maximum.reduce(N.fabs(N.ravel(N.dot(a, inv_a)-N.identity(5))))

# 行列の固有値を印字
print "Eigenvalues: ", LA.eigenvalues(a)

モジュール名の略記名を作る方法に注意して下さい。モジュールは特殊なデータオブジェクトに過ぎないので、変数に割り付けることも可能なのです。

フーリエ変換

FFT モジュールには、一および二次元の種々の高速フーリエ変換を行う関数があります。FFTPACK というパッケージを使っています。全ての関数は任意の階数の配列に適用できます；変換すべき次元のデフォルトは最後の一つあるいは二つですが、それ（ら）を指定することも可能です。

関数 fft(a) は、標準的な複素１次元 FFT を、inverse_fft(a) は逆変換を行います。両者とも、オプションでデータ列の長さ（デフォルトは配列の長さ）を指定する第二引数を受け取ります；それが配列長よりも大きい場合、適当数のゼロを追加します。オプションの第三引数は、変換を実行する際の変数となるべき次元を指定します。

関数 fft2d(a) は、二次元 FFT を行います。オプションの第二引数により、変換すべき部分配列の形状 (shape)（長さ２のタプル）を指定します；デフォルトでは、配列全体が変換されます。オプションの第三引数（これも長さ２のタプル）は、変換の二変数の次元を指定します。

パラメータ最適化

LeastSquares モジュールには、関数 leastSquaresFit(model, parameter, data) が提供されています。これは、Levenberg-Marqardt アルゴリズムを用いた一般的な非線形最小二乗法最適化を行います。

引数 model は、当てはめるべき関数を指定します。これは二つの引数を渡して呼び出します：一番目は最適化すべきパラメータからなるタプル、二番目はデータ点の第一要素（下記参照）です。返り値は一つの数です。

引数 parameter は、最適化するパラメータの初期値からなるタプルです。この方法の収束速度は、初期値の性質に依存するので、注意して選ばなくてはなりません。

引数 data は、モデル関数を適合すべきデータ点の列です。各々のデータ点は長さ２または３の列からなります。その第一要素は、モデルの独立変数を指定するもので、モデル関数に第一パラメータとして渡され、それ以外の用途には使われません。各データ点のタプルの第二要素は、モデル関数の返値が出来るだけ一致すべき対象となる数値です。第三要素（デフォルト値は１）はデータ点の統計的分散、すなわち、最適化の手続きにおけるその点の統計的重みの逆数です。

例：

from LeastSquares import leastSquaresFit
import Numeric

# 数学的モデル
def exponential(parameters, x):
    a = parameters[0]
    b = parameters[1]
    return a*Numeric.exp(-b/x)

# 適合化されるべきデータ
data = [(100, 4.999e-8),
	(200, 5.307e+2),
	(300, 1.289e+6),
	(400, 6.559e+7)]

fit = leastSquaresFit(exponential, (1e13, 4700), data)

print "Fitted parameters:", fit[0]
print "Fit error:", fit[1]

格子点で与えられた関数

数値計算では、関数の多くは格子点上で計算された一組の値で定義されます。格子点以外でこのような関数の値を求めるには、内挿が必要です。 Scientific.Functions.Interpolation というモジュールは、このような関数を定義し、よく使われる演算を提供します。格子は直交でなくてはなりませんが、等間隔である必要はありません。格子は各点の座標を指定する一組の軸で定義されます。一次元あたり軸が一つです。

関数は、InterpolatingFunction(grid, value) で定義します。ここで、grid は軸からなる列で、その長さによって格子の次元が決まります。引数 values は、格子の次元と同じまたはより高い階数を持つ配列です。階数がより高い場合には、関数は多価です（例えば、ベクトル場）。第三引数はオプションで、格子の外側の関数値を指定します。なにも与えられない場合、格子の外側で関数値を求めるとエラーになります。

内挿関数は、他の任意の関数と同様に呼び出すことができます。引数の個数は、格子の次元と等しくなくてはなりません。また、引数の型は整数または実数です。

内挿関数には、単純な求値だけでなく、幾つかの便利な演算操作が提供されています：

function.derivative(variable) は、variable 番目の変数（デフォルトはゼロ）に関する（数値的）微分を表す別の内挿関数を返します。
function.integral(variable) は、variable 番目の変数（デフォルトはゼロ）に関する（数値的）不定積分を表す別の内挿関数を返します。
function.definiteIntegral(variable) は、variable 番目の変数（デフォルトはゼロ）に関する（数値的）定積分を表す別の内挿関数を返します。この操作は、次元を一つ減らします；一次元関数の場合は、一個の数値が返されます。
function.selectInterval(first, last, variable) は、variable 番目の次元に沿って first から last までの間隔に縮小された、より小さな格子で定義された内挿関数を返します。これは、定積分の計算において最も有用です。

例：

from Scientific.Functions.Interpolation import InterpolatingFunction
import Numeric; N = Numeric

# 軸 (axis) を１０個の点で定義
axis = N.arrayrange(0, 1.1, 0.1)

# 値の計算（ここでは単純な平方根）
values = N.sqrt(axis)

# 内挿関数の生成
f = InterpolatingFunction((axis,), values)

# 比較のために値の表を印字
for x in N.arrayrange(0, 1., 0.13):
    print "%5.2f: %10.5f %10.5f" % (x, N.sqrt(x), f(x))

可視化 (VRML)

複雑なデータ（特に三次元）を描出するには、単純な図示では必ずしも十分ではありません。 Scientific.Visualization.VRML モジュールを用いれば、（立方体や球などの）標準的なオブジェクトからなる三次元的情景（シーン）を直ちに構築することができます。これらのシーンは、VRML (Virtual Reality Modelling Language) フォーマットでファイルに書き出し、VRML 表示ソフトで見ることができます。 VRML 表示ソフトは大抵のシステムで使用できます。 VRML と VRML 表示ソフトについての詳細は、サンディエゴ・スーパーコンピュータセンターの VRML Repository を参照してください。

一々ファイルを書き出して VRML 表示ソフトを起動する代りに、Python プログラムから直接 VRML シーンを表示したければ、インストールした VRML 表示ソフトの名前を Unix の環境変数 VRMLVIEWER に設定する必要があります。

注：このコースでは、VRML モジュールの基本的使用法のみ扱います。詳細については、ソースコードを見てください。

典型的な一連の操作は次のようになります：

VRML シーンを Scene(objects) によって生成。ここで、objects は、VRML オブジェクトのリスト。空のリストを渡して、空のシーンを生成することも可能。
VRML オブジェクトを生成（下参照）し、scene.addObject(object) でシーンに追加。
scene.writeToFile(filename) でシーンをファイルに書き出す。または、直接 scene.view() によって描画。

VRML オブジェクトは、幾何的データによって定義されます。例えば、球は中心を表す一つのベクトルと半径を表す一つの数で定義されます。さらに、オブジェクトはその外観に係わる一つ以上の性質 (properties) を持つことができます。最も重要な性質は、素材 (material) で、色、透明度、反射度などを定義します。このコースでは、予め決められた種々の色を持つ広汎的素材 (diffuse materials) 、すなわち、透明でもなく反射もない素材のみを使います。全ての可能な性質に関しては、VRML に関する書籍を参照して下さい。

属性は常にオプションで、キーワードパラメータによって指定します。例えば、素材は material=the_material のように指定します。指定の無い場合のデフォルトは、白の広汎的素材になります。

素材は、DiffuseMaterial(color) によって生成されます。ここで、color は色の名前です。使用できる色は、black, white, grey, red, green, blue, yellow, magenta, cyan です。例えば 'light blue' のように、これら全色の先頭に 'light ' または 'dark '（空白に注意）を付けることができます。

使用できる VRML オブジェクトは以下の通りです：

Sphere(center, radius) は、球を生成します。
Cube(center, edge) は、与えられた点を中心とする立方体を生成します。各辺は、x, y, z 座標に沿っています。
Cylinder(point1, point2, radius) は、軸が point1 から point2 までで、指定された半径を持つ円柱を生成します。
Cone(point1, point2, radius) は、point1 を中心とする半径 radius の円を底面とし、先端が point2 にある円錐を生成します。
Line(point1, point2) は、point1 から point2 までの線分を生成します。
Arrow(point1, point2, radius) は、軸半径が radius で、point1 から point2 を指す矢印を生成します。

次のプログラム例は、PDB ファイルを読み込み、各々が黒い球で表され、隣り合うものが赤い線で結ばれたα炭素からなるシーンを構築します。

from Scientific.IO.TextFile import TextFile
from Scientific.IO.FortranFormat import FortranFormat, FortranLine
from Scientific.Geometry import Vector
import string
from Scientific.Visualization import VRML

generic_format = FortranFormat('A6')
atom_format = FortranFormat('A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,' +
                            '3X,3F8.3,2F6.2,7X,A4,2A2')

# PDB ファイルを読み込み、全てのα炭素原子の位置からなるリストを作る。
c_alpha_positions = []
for line in TextFile('protein.pdb'):
    record_type = FortranLine(line, generic_format)[0]
    if record_type == 'ATOM  ' or record_type == 'HETATM':
        data = FortranLine(line, atom_format)
        atom_name = string.strip(data[2])
        position = Vector(data[8:11])
	if atom_name == 'CA':
	    c_alpha_positions.append(position)

# VRML シーン、球の素材と半径、線の素材を生成。
scene = VRML.Scene([])
c_alpha_material = VRML.DiffuseMaterial('black')
c_alpha_radius = 0.2
line_material = VRML.DiffuseMaterial('red')

# 球を生成し、シーンに配置。
for point in c_alpha_positions:
    scene.addObject(VRML.Sphere(point, c_alpha_radius,
				material=c_alpha_material))

# 線を生成し、シーンに配置。
for i in range(len(c_alpha_positions)-1):
    point1 = c_alpha_positions[i]
    point2 = c_alpha_positions[i+1]
    scene.addObject(VRML.Line(point1, point2, material=line_material))

# シーンを描画。
scene.view()

その他の科学計算モジュール

このコースで網羅されていない科学計算モジュールは、まだまだあります。マニュアルを参照して下さい。

練習問題

対角化法により、正方行列に算術関数 (exp, sqrt など) を適用する関数を書きなさい。手順は次の通りです： 1) 行列の固有値と固有ベクトルを見つける。 2) 固有値に算術関数を適用する。 3) 固有ベクトル行列を両側から掛け、対角化の逆操作を行う。
テーラー展開による別法で行列の指数関数を計算し、このプログラムをテストしなさい。
同一のタンパク質の異なる二配置を表す二つの PDB ファイルを読み込み、各α炭素原子について配置１から配置２へ向う矢印を描いた VRML シーンを制作するプログラムを書きなさい。
ヒント：ファイルからα炭素の位置を取り出すプログラムを関数化し、これを二回呼ぶようにしましょう。