こんにちは｜こんばんは。カエルのアイコンで活動しております @kyamaz 
です。
はじめに
本稿では、SymPyを使って例題を解きます。最初にsympyを読み込みます。
>>> import sympy as sym
本稿の環境
本稿のために使用した環境は以下となります。macOS: Sonoma 14.2.1 (chip: Apple M1)
Homebrew: 4.2.9
anyenv: 1.1.5 (homebrewにてインストール)
Python: 3.12.1
SymPy: 1.12
ピタゴラス数を生成する行列
$a^2+b^2=c^2$を満たす正の整数の組$(a,b,c)$のことをピタゴラス数と言います。その中で，$a,b,c$の最大公約数が$1$であるものを原始ピタゴラス数と言います。$(3,4,5)$に対して次の３つの行列$U,A,D$をかける操作を何度か繰り返すことで、すべての原始ピタゴラス数を生成できます。1
$
A=\displaystyle \left(\begin{matrix}1 & 2 & 2\\2 & 1 & 2\\2 & 2 & 3\end{matrix}\right)
,U=\displaystyle \left(\begin{matrix}1 & -2 & 2\\2 & -1 & 2\\2 & -2 & 3\end{matrix}\right)
,D=\displaystyle \left(\begin{matrix}-1 & 2 & 2\\-2 & 1 & 2\\-2 & 2 & 3\end{matrix}\right)
$
次のサイトに詳しく証明等が紹介されています。
本稿では、行列$A$の対角化と行列$U,D$のジョルダン標準形を求めてみます。数学的な意義を述べるものではなく、行列の対角化やジョルダン標準形の例題としてSymPyを使った演習問題として解いてみます。
行列 A の対角化
行列$A$の固有値は以下のように求められます。
>>> A=sym.Matrix([[1,2,2],[2,1,2],[2,2,3]])
>>> A.eigenvals()
{-1: 1, 3 - 2*sqrt(2): 1, 2*sqrt(2) + 3: 1}
固有値は$-1,3\pm2\sqrt{2}$の３つです。
そして、対角化はdiagonalize()を使って次のように求めます。ここで天下り的に、対角化行列PAを$\frac{1}{\sqrt{2}}$倍したPPAを代わりに対角化行列$P_{A}$として用います。
>>> PA, DA = A.diagonalize()
>>> DA
Matrix([
[-1, 0, 0],
[ 0, 3 - 2*sqrt(2), 0],
[ 0, 0, 2*sqrt(2) + 3]])
>>> PPA=sym.simplify(PA/sym.sqrt(2))
>>> PPA
Matrix([
[-sqrt(2)/2, -1/2, 1/2],
[ sqrt(2)/2, -1/2, 1/2],
[ 0, sqrt(2)/2, sqrt(2)/2]])
>>> PPA.inv()
Matrix([
[-sqrt(2)/2, sqrt(2)/2, 0],
[ -1/2, -1/2, sqrt(2)/2],
[ 1/2, 1/2, sqrt(2)/2]])
$
P_{A}=\left(\begin{matrix}- \frac{\sqrt{2}}{2} & - \frac{1}{2} & \frac{1}{2}\\\frac{\sqrt{2}}{2} & - \frac{1}{2} & \frac{1}{2}\\0 & \frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2}\end{matrix}\right)
\ ,P_{A}^{-1}=\left(\begin{matrix}- \frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} & 0\\- \frac{1}{2} & - \frac{1}{2} & \frac{\sqrt{2}}{2}\\\frac{1}{2} & \frac{1}{2} & \frac{\sqrt{2}}{2}\end{matrix}\right)
$
対角化の公式より$P_{A}^{-1}AP_{A}=D_{A}$ですので$A=P_{A}D_{A}P_{A}^{-1}$となります。検算してみましょう。
>>> sym.simplify(PPA*DA*PPA.inv())
Matrix([
[1, 2, 2],
[2, 1, 2],
[2, 2, 3]])
行列 U,D のジョルダン標準形
行列$U,D$の固有値は以下のようにどちらとも重解になります。
>>> U=sym.Matrix([[1,-2,2],[2,-1,2],[2,-2,3]])
>>> U.eigenvals()
{1: 3}
>>> D=sym.Matrix([[-1,2,2],[-2,1,2],[-2,2,3]])
>>> D.eigenvals()
{1: 3}
重解の場合は、対角化ではなくジョルダン標準形を計算します。ジョルダン標準形はjordan_form()を使うと計算できます。なお、こちらも天下り的に、それぞれの変換行列$P_{U},P_{D}$はPU,PDを$\frac{1}{2}$倍したPPU,PPDを用います。
>>> PU, JU = U.jordan_form()
>>> JU
Matrix([
[1, 1, 0],
[0, 1, 1],
[0, 0, 1]])
>>> PPU=sym.simplify(PU/2)
>>> PPU
Matrix([
[ 0, -1, 0],
[-2, -1, 1/2],
[-2, -1, 0]])
>>> PD, JD = D.jordan_form()
>>> JD
Matrix([
[1, 1, 0],
[0, 1, 1],
[0, 0, 1]])
>>> PPD=sym.simplify(PD/2)
>>> PPD
Matrix([
[-2, -1, 1/2],
[ 0, -1, 0],
[-2, -1, 0]])
ジョルダン標準形は、行列$U,D$とも同じで次のようになります。
$
J_{U}=J_{D}=\left(\begin{matrix}1 & 1 & 0\\0 & 1 & 1\\0 & 0 & 1\end{matrix}\right)
$
$U$の変換行列とその逆行列は次のようになります。
$
P_{U}=\left(\begin{matrix}0 & -1 & 0\\-2 & -1 & \frac{1}{2}\\-2 & -1 & 0\end{matrix}\right)
\ ,P_{U}^{-1}=\left(\begin{matrix}\frac{1}{2} & 0 & - \frac{1}{2}\\-1 & 0 & 0\\0 & 2 & -2\end{matrix}\right)
$
$D$の変換行列とその逆行列は次のようになります。
$
P_{D}=\left(\begin{matrix}-2 & -1 & \frac{1}{2}\\0 & -1 & 0\\-2 & -1 & 0\end{matrix}\right)
\ ,P_{D}^{-1}=\left(\begin{matrix}0 & \frac{1}{2} & - \frac{1}{2}\\0 & -1 & 0\\2 & 0 & -2\end{matrix}\right)
$
$U=P_{U}J_{U}P_{U}^{-1}$, $D=P_{D}J_{D}P_{D}^{-1}$を計算して検算してみましょう。
>>> PPU*JU*PPU.inv()
Matrix([
[1, -2, 2],
[2, -1, 2],
[2, -2, 3]])
>>> PPU*JU*PPU.inv()-U
Matrix([
[0, 0, 0],
[0, 0, 0],
[0, 0, 0]])
>>> PPD*JD*PPD.inv()
Matrix([
[-1, 2, 2],
[-2, 1, 2],
[-2, 2, 3]])
>>> PPD*JD*PPD.inv()-D
Matrix([
[0, 0, 0],
[0, 0, 0],
[0, 0, 0]])
おわりに
本稿では、pythonのコマンドラインでの例で進めましたが、Jupyter Notebook で実行すると行列の表示を見やすくするなります。例のプロンプト>>>のあとの操作を入力すると実行できます。各自でお試し頂けると幸いです。
ピタゴラス数を生成する$3\times3$行列で見てきましたが、参考文献1にもあるように次の$2\times2$行列も知られています。
$
A=\left(\begin{matrix}2 & 1\\1 & 0\end{matrix}\right)
\ ,U=\left(\begin{matrix}2 & -1\\1 & 0\end{matrix}\right)
\ ,D=\left(\begin{matrix}1 & 2\\0 & 1\end{matrix}\right)
$
本稿と同じように、行列$A$は対角化できますが、行列$U,D$は固有値が重解となり、ジョルダン標準形で表される行列です。
ピタゴラス数を生成する３つの行列は興味深い行列ではありますが、本稿は数学的な話題を抜きに単なる計算問題とさせていただきました。