はじめに

3次方程式の解の公式は調べれば出てくるし、33な行列の固有値を手計算で求める方法も出てきます。
しかし、33な行列の固有値を解の公式を用いて数値計算しているものが見つけられなかったので、ここに残しておきます。

なお、前半では数式を、後半ではJavaScriptによるプログラムを、ぞれぞれ説明します。

数式で求める

3次方程式の解の公式

3次方程式

a'_3 x^3 + a'_2 x^2 + a'_1 x + a'_0 = 0

は、全体を$a'_3$で除算することにより

x^3 + a_2 x^2 + a_1 x + a_0 = 0

とすることができます。この代数的な解はWikipediaによると、カルダノの公式より

x = \omega^k \sqrt[3]{-\frac{q}{2} + \sqrt{\left(\frac{q}{2}\right)^2 + \left(\frac{p}{2}\right)^3}} + \omega^{3 - k} \sqrt[3]{-\frac{q}{2} - \sqrt{\left(\frac{q}{2}\right)^2 + \left(\frac{p}{2}\right)^3}} - \frac{a_2}{3}
\\
p = a_1 - \frac{a_2^2}{3}
\\
q = a_0 - \frac{1}{3}a_1 a_2 + \frac{2}{27}a_2^3

となります。

$\omega$は$1$の立方根の一つで、

\begin{eqnarray}
\omega^0 &=& \omega^3 = 1 \\
\omega^1 &=& \frac{-1 + \sqrt 3 i}{2} \\
\omega^2 &=& \frac{-1 - \sqrt 3 i}{2} \\
\end{eqnarray}

といった感じです。なおこれは明示的には使いません。後の複素数の立方根を求める段階で、それぞれに対応する立方根が計算されます。

固有値

固有値は等式

\boldsymbol A \boldsymbol x = \lambda \boldsymbol x

が成り立つ$\lambda$のことです。なお、$\boldsymbol x$は固有ベクトルとなりますが、本記事では扱いません。

この時、固有値を求める方程式は

| \lambda \boldsymbol I - \boldsymbol A | = 0

となります。$| \cdot |$は行列式を表します。

特に3行3列の場合、上の式は

\begin{eqnarray}
\left | \lambda \boldsymbol I -
\begin{pmatrix}
a_{00} & a_{01} & a_{02} \\
a_{10} & a_{11} & a_{12} \\
a_{20} & a_{21} & a_{22} \\
\end{pmatrix}
\right | &=& 0
\\
\begin{vmatrix}
\lambda - a_{00} & -a_{01} & -a_{02} \\
-a_{10} & \lambda - a_{11} & -a_{12} \\
-a_{20} & -a_{21} & \lambda - a_{22} \\
\end{vmatrix} &=& 0
\\
(\lambda - a_{00})(\lambda - a_{11})(\lambda - a_{22}) - a_{01}a_{12}a_{20} - a_{02}a_{10}a_{21} - (\lambda - a_{00})a_{12}a_{21} - (\lambda - a_{11})a_{02}a_{20} - (\lambda - a_{22})a_{01}a_{10} &=& 0
\end{eqnarray}

となり、$\lambda$について整理すると

\lambda^3 + a_2\lambda^2 + a_1\lambda + a_0 = 0

が得られます。ただし、$a_2, a_1, a_0$はそれぞれ以下のとおりです。

\begin{eqnarray}
a_2 &=& -(a_{00} + a_{11} + a_{22}) \\
a_1 &=& a_{00}a_{11} + a_{11}a_{22} + a_{22}a_{00} - a_{12}a_{21} - a_{02}a_{20} - a_{01}a_{10} \\
a_0 &=& -| \boldsymbol A | = - 
\begin{vmatrix}
a_{00} & a_{01} & a_{02} \\
a_{10} & a_{11} & a_{12} \\
a_{20} & a_{21} & a_{22} \\
\end{vmatrix}
\end{eqnarray}

この式を先述した3次方程式の解に適用することで、3行3列の行列の固有値を代数的に求めることができます。

複素数上の計算

全体の大まかな処理が導けたのですが、計算途中で複素数を扱う必要があります。
残念ながらJavaScriptには複素数を扱う型はありません。（いくつかのライブラリには存在します）

行列が実数の場合、複素数の範囲で計算が必要となるのは、加算、平方根、立方根の三つです。
この三つのうち加算は難しくないと思いますので、複素数の平方根および立方根の計算方法を示します。
※減算もありますが、引く数が実数なので$-1$を乗算することで加算扱いします。

このサイトを参考に、複素数$z = a + bi$の偏角を$\theta$、絶対値を$|z| = \sqrt{a^2 + b^2}$とすると、平方根および立方根はそれぞれ

\begin{array}{l}
\sqrt {|z|} \left(\cos \left(\frac{\theta}{2} + \pi k \right) + i\sin \left(\frac{\theta}{2} + \pi k \right) \right) & (k = 0, 1)
\\
\sqrt[3] {|z|} \left(\cos \left(\frac{\theta}{3} + \frac{2\pi}{3}k \right) + i\sin \left(\frac{\theta}{3} + \frac{2\pi}{3}k \right) \right) & (k = 0, 1, 2)
\end{array}

となります。

なおここで、三角関数$\sin$および$\cos$、複素数の偏角、実数の平方根および立方根を計算する必要がありますが、これはJavaScriptではMath.sin、Math.cos、Math.atan2、Math.sqrtおよびMath.cbrtで計算できます。他の多くの言語でも同様の関数が用意されていますので、これらの計算の方法については割愛します。

プログラムで求める

とりあえず、実数の行列から求める方法を示します。
複素数の行列の場合は、四則演算を複素数のものに置き換えてもらえれば問題ないかと思います。

複素数のクラス

以下の通り複素数を扱うクラスを定義します。

class Complex {
    constructor(real = 0, imag = 0) {
        this.real = real
        this.imag = imag
    }

    abs() {
        return Math.sqrt(this.real ** 2 + this.imag ** 2)
    }

    add(other) {
        if (typeof other === 'number') {
            return new Complex(this.real + other, this.imag)
        }
        return new Complex(this.real + other.real, this.imag + other.imag)
    }

    sqrt() {
        const th = Math.atan2(this.imag, this.real) / 2
        const s = Math.sqrt(this.abs())
        return [
            new Complex(Math.cos(th) * s, Math.sin(th) * s),
            new Complex(Math.cos(th + Math.PI) * s, Math.sin(th + Math.PI) * s),
        ]
    }

    cbrt() {
        const th = Math.atan2(this.imag, this.real) / 3
        const s = Math.cbrt(this.abs())
        return [
            new Complex(Math.cos(th) * s, Math.sin(th) * s),
            new Complex(Math.cos(th + (2 * Math.PI) / 3) * s, Math.sin(th + (2 * Math.PI) / 3) * s),
            new Complex(Math.cos(th + (4 * Math.PI) / 3) * s, Math.sin(th + (4 * Math.PI) / 3) * s),
        ]
    }
}

このような、イミュータブルのような雰囲気のクラスです。
realが実部、imagが虚部の値となり、絶対値、（数値および複素数との）加法、平方根、立方根をそれぞれ定義しています。

※本当はrealとimagをgetterにしてあげた方がイミュータブルな意識が出るとは思いますが、簡略化のためこのようにしています。

係数の計算

まずは行列となる適当な変数を用意し、そこから三次方程式の係数$a_2, a_1, a_0$を計算します。

const mat = []
for (let i = 0; i < 3; i++) {
    mat[i] = []
    for (let j = 0; j < 3; j++) {
        mat[i][j] = Math.random()
    }
}

const a2 = -mat[0][0] - mat[1][1] - mat[2][2]
const a1 = mat[0][0] * mat[1][1]
         + mat[0][0] * mat[2][2]
         + mat[1][1] * mat[2][2]
         - mat[1][2] * mat[2][1]
         - mat[0][2] * mat[2][0]
         - mat[0][1] * mat[1][0]
const a0 = -mat[0][0] * mat[1][1] * mat[2][2]
         - mat[0][1] * mat[1][2] * mat[2][0]
         - mat[0][2] * mat[1][0] * mat[2][1]
         + mat[0][0] * mat[1][2] * mat[2][1]
         + mat[1][1] * mat[0][2] * mat[2][0]
         + mat[2][2] * mat[0][1] * mat[1][0]

媒介変数の計算

媒介変数$p, q$および、$t = \left(\frac{q}{2}\right)^2 + \left(\frac{p}{2}\right)^3$をそれぞれ計算します。

const p = a1 - a2 ** 2 / 3
const q = a0 - (a1 * a2) / 3 + (a2 ** 3 * 2) / 27
const t = (q / 2) ** 2 + (p / 3) ** 3

各項の計算

解の三つある項のうち、前の二つについて計算します。
ここからは複素数で計算する必要があります。

let [u3, v3] = new Complex(t).sqrt()
u3 = u3.add(-q / 2)
v3 = v3.add(-q / 2)

const [u0, u1, u2] = u3.cbrt()
const [v0, v1, v2] = v3.cbrt()

解の計算

求められた解を適切に組み合わせ、$-\frac{a_2}{3}$を加算して解を計算します。

const e0 = u0.add(v0).add(-a2 / 3)
const e1 = u1.add(v2).add(-a2 / 3)
const e2 = u2.add(v1).add(-a2 / 3)

ここで得られたe0、e1、e2が解になります。

実数への置き換え

実数のみが必要な場合は、適当な許容範囲を設定し、以下のようにするとよいと思います。

const tol = 1.0e-12

const re0 = Math.abs(e0.imag) < tol ? e0.real : Number.NaN
const re1 = Math.abs(e1.imag) < tol ? e1.real : Number.NaN
const re2 = Math.abs(e2.imag) < tol ? e2.real : Number.NaN

おわりに

これは、個人的に作成している機械学習ライブラリにおいて実装したものを、抜粋・簡略化したものです。

もともとのコードは、固有値計算はここ、複素数計算はここにあります。