Python（SymPy）で学ぶVandermonde行列式

数学の具体的な計算にPython（SymPy）を使って、数学もPython（SymPy）も同時に学んでしまいましょう。今回はPython（SymPy）を使ってVandermonde（ヴァンデルモンド）行列式の計算をしてみたいと思います。

本記事はPython（SymPy）で行列の扱い方を解説した過去記事

pianofisica.hatenablog.com

の続編という位置付けになりますので、あわせて読んでみてください。

Vandermonde（ヴァンデルモンド）行列式
公式
具体例：４次正方行列の場合
具体例：少しスマートなコーディング

Vandermonde（ヴァンデルモンド）行列式

Vandermonde（ヴァンデルモンド）行列式は、次のヴァンデルモンド行列

$\qquad\displaystyle{V=\left[\begin{array}{ccccc} 1 & x_{1} & x_{1}^{2} & \cdots & x_{1}^{n-1} \\ 1 & x_{2} & x_{2}^{2} & \cdots & x_{2}^{n-1} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_{n} & x_{n}^{2} & \cdots & x_{n}^{n-1} \end{array}\right]}$

の行列式のことをいいます。行列サイズが小さい場合（ $n=2,3$ ）に具体的に求めてみましょう。

$\qquad\displaystyle{\begin{aligned} \det V_2&=\det\left[\begin{array}{cc} 1 & x_{1} \\ 1 & x_{2} \end{array}\right]=x_2-x_1 \\ \det V_3&=\det\left[\begin{array}{ccc} 1 & x_{1} & x_{1}^{2} \\ 1 & x_{2} & x_{2}^{2} \\ 1 & x_{3} & x_{3}^{2} \end{array}\right]=x_1x_2^2+x_2x_3^2+x_3x_1^2-x_1^2x_2-x_2^2x_3-x_3^2x_1\\ &=x_1x_2(x_2-x_1)+x_3^2(x_2-x_1)-x_3(x_2^2-x_1^2)\\ &= (x_2-x_1)\{x_3^2-(x_1+x_2)x_3+x_1x_2\}\\ &=(x_2-x_1)(x_3-x_1)(x_3-x_2) \end{aligned}}$

以上のように因数分解されることがわかります。これは以下の事実を反映したものです。つまり、行列式の変数 $x_i$ と $x_j$ （ $i\neq j$ ）が一致すると、Vandermonde行列式がゼロになることが行列式の一般論からわかるため、 $\det V$ は $(x_i-x_j)$ の因数を持ちます。より一般に次の公式が成り立ちます。

公式

$\qquad\displaystyle{ \det V=\prod_{1 \leq i {<} j \leq n} \left(x_{j}-x_{i}\right) }$

その数学的な証明は線形代数学の教科書に譲るとして、この公式が成り立っていることをSymPyで計算して確かめてみましょう。というのも、３次正方行列くらいまでであれば、手計算で行列式を展開・整理して因数分解することもできるでしょうが、４次、５次、６次となってくると大変になってきます。そういった計算にこそ数式処理ライブラリのSymPyが役に立ちます。

具体例：４次正方行列の場合

SymPyで４次のVandermonde行列式

$\qquad \displaystyle{ \det V_4=\det\left[\begin{array}{cccc} 1 & x_{1} & x_{1}^{2} & x_{1}^{3} \\ 1 & x_{2} & x_{2}^{2} & x_{2}^{3} \\ 1 & x_{3} & x_{3}^{2} & x_{3}^{3} \\ 1 & x_{4} & x_{4}^{2} & x_{4}^{3} \end{array}\right] }$

を考えてみます。行列をPythonに入力するには

from sympy import *
x1, x2, x3, x4 = var('x1, x2, x3, x4')

V4 = Matrix([
[1,x1,x1**2,x1**3],
[1,x2,x2**2,x2**3],
[1,x3,x3**2,x3**3],
[1,x4,x4**2,x4**3]])

とします。その行列式は

VD4 = V4.det()

によって計算できますが、このままでは展開式が表示されるだけですので、factor関数によって因数分解しましょう。

factor(VD4)

すると、行列式が

$\qquad\displaystyle{\det V_4=(x_1-x_2)(x_1-x_3)(x_1-x_4)(x_2-x_3)(x_2-x_4)(x_3-x_4)}$

と因数分解され、公式と一致することが確かめられます。

具体例：少しスマートなコーディング

さて、以上のようにして行列の成分を具体的に一つひとつ入力することで求めたい量の計算はできるわけですが、あまり汎用性の高い書き方であるとはいえません。というのは、５次、６次と拡張して同様の計算をしようとしたときに、その都度全く同じものをコピーし、必要な成分を書き加えないといけないためです。この不満足な点を改良するべく、自然数を与えればそれに対応する次数のVandermonde行列式を与えるような関数をSymPyで定義してみましょう。まず、Vandermonde行列の $(i,j)$ -成分が

$\qquad\displaystyle{V_{ij}=x_i^{j-1} }$

で与えられることに注意しておきます。また、SymPyのMatrix関数がリストから行列に変換する関数であることに留意しておきましょう。そこで以下のようにして、列と行それぞれについてfor文を回してMatrix関数に入力するリスト L(n) を作成します。

from sympy import * 

x1, x2, x3, x4, x5, x6 = var('x1, x2, x3, x4, x5, x6')
X = [1, x1, x2, x3, x4, x5, x6]

def L(n):
    L = []
    R = []
    for i in range(1,n+1):
        for j in range(1, n+1):
            R.append(X[i]**(j-1))
        L.append(R)
        R = []
    return L

このリストをMatrix関数に入力して、その行列式を（因数分解したうえで）出力する関数を定義してみます。

def VD(n):
    return factor(Matrix(L(n)).det())

すると

VD(2)

VD(3)

VD(4)

などとして、上で次数ごとに個別に求めていた計算を一つのプログラムで実現することができました。ただし、６次以上になると計算が重くなってきますので注意してください。

キーワード：Python、行列、Vandermonde行列式