simplestarの技術ブログ

目的を書いて、思想と試行、結果と考察、そして具体的な手段を記録します。

Math:Unityで対称行列の対角化と同時に固有値・固有ベクトルを求める方法

ゲーム業界でも固有値固有ベクトルという単語はよく出てきますが、対称行列を与えられた状態から具体的な解き方までご存知の方はそう多くいません。
答えを先に言うと、対称行列を対角化すると、その時点で対角行列の各要素は固有値、左右の回転行列はその固有値に対する固有ベクトルで構成されることになります。
数学カテゴリ第二回目の今回は、Unity での対称行列の対角化と、固有値固有ベクトルの求め方を示したいと思います。

今回示すプログラムを走らせると、次の動作を確認できます。
f:id:simplestar_tech:20170412084653g:plain
元の行列の形をまったく変えずに、左下の行列が対角行列に相似変換されていく様子と、それを実現する左右の回転行列(相似変換行列)が更新されていく様子が確認できます。

さて、固有値固有ベクトルの求め方だけでネットで調べると、固有方程式を立てて、高次方程式の解として固有値を得るというものが多く見つかります。
私が固有方程式を解く方法を嫌う理由としては、与えられた行列が解に近い形をしているとまず解けませんし、立てた高次方程式が結局解けないという問題や、重解となったときに対応する複数の固有ベクトル固有値から求められないという問題があるからです。(固有方程式は筆算でよく使われますが、数値計算プログラムとしては欠陥だらけですね)

もう少し資料を探してみると、ギブンス回転による相似変換(ヤコビ法)、ハウスホルダー変換、三重対角行列変換とQR法(ウィルキンソンの移動法)などが見つかります。
速度の速い順に手法を紹介しましたが、対照的に精度が上がっていきます。
ゲーム業界で使われるだろう次元数は最大で4くらいですので、この中で最速のギブンス回転による相似変換(ヤコビ法)が、ゲーム業界向けに最適な手法と言えます。

そこで、今回はギブンス回転による相似変換(ヤコビ法)によって、対称行列の対角化と固有値固有ベクトルを求めるプログラムを示します。

以下の通り

        /// <summary>
        /// 対称行列の固有値・固有ベクトル分解
        /// </summary>
        /// <param name="A">[in]分解する対象のNxNの実対称行列</param>
        /// <param name="D">[out]実対称行列に対応する、対角成分を固有値とする行列</param>
        /// <param name="L">[out]その列の固有値に対する固有列ベクトルで構成される行列</param>
        /// <param name="e">収束判定用最小誤差</param>
        /// <param name="maxCount">最大ギブンス回転回数</param>
        /// <returns>最大回数以下で収束した時trueを返す</returns>
        public static bool EigenDecomposition(float[][] A, float[][] D, float[][] L, float e, int maxCount)
        {
            int sz = A.Length;
            for (int i = 0; i < sz; i++)
            {
                for (int j = 0; j < sz; j++)
                {
                    D[i][j] = A[i][j];
                    if (j != i)
                    {
                        L[i][j] = 0;
                    }
                }
                L[i][i] = 1;
            }

            int count = 0;
            while (count++ < maxCount)
            {
                int row = 0;
                int col = 0;
                AbsoluteMaxOffDiagonalElementPosition(D, out row, out col);
                if (e > Mathf.Abs(D[row][col]))
                    break;
                // pp = qq = (pp - qq) * sin2θ + 2 * pq * cos2θ = 0
                // αsin2θ + βcos2θ = 0
                float alpha = (D[row][row] - D[col][col]);
                float beta = - 2 * D[row][col];
                // (αcos2θ) * (αcos2θ) + (αsin2θ) * (αsin2θ) = α * α 
                // (α * α + β * β)cos2θ^2 = α * α
                // cos2θ = Abs(α) / Sqrt(α * α + β * β)
                float cos2 = Mathf.Abs(alpha) / Mathf.Sqrt(alpha * alpha + beta * beta);
                // cosθ
                float cos1 = Mathf.Sqrt((1 + cos2) / 2);
                // sinθ
                float sin1 = Mathf.Sqrt((1 - cos2) / 2);
                if (0 > alpha * beta)
                    sin1 *= -1.0f;

                GivensRotation(D, row, col, sin1, cos1);

                SimilarityTransformationLeft(L, row, col, sin1, cos1);
            }

            if (count >= maxCount)
                return false;

            return true;
        }

        /// <summary>
        /// ギブンス回転
        /// </summary>
        /// <param name="A">入出力の対角化する行列</param>
        /// <param name="row">非対角成分における絶対値最大要素の行</param>
        /// <param name="col">非対角成分における絶対値最大要素の列</param>
        /// <param name="sin1">その要素=0としたときの回転量sinθ</param>
        /// <param name="cos1">その要素=0としたときの回転量cosθ</param>
        private static void GivensRotation(float[][] A, int row, int col, float sin1, float cos1)
        {
            float sin1cos1 = sin1 * cos1; // sinθcosθ
            float sin1_2 = sin1 * sin1; // sinθ^2
            float cos1_2 = cos1 * cos1; // cosθ^2

            // pp * cosθ^2 + qq * sinθ^2 - 2 pq * sinθcosθ
            float tmp1 = A[row][row] * cos1_2 + A[col][col] * sin1_2 - 2 * A[row][col] * sin1cos1;
            // pp * sinθ^2 + qq * cosθ^2 + 2 pq * sinθcosθ
            float tmp2 = A[row][row] * sin1_2 + A[col][col] * cos1_2 + 2 * A[row][col] * sin1cos1;

            A[row][col] = A[col][row] = 0; // pq, qp = (pp - qq) / 2 * sin2θ + pq * cos2θ = 0
            A[row][row] = tmp1; // pp
            A[col][col] = tmp2; // qq

            int sz = A.Length;

            // A[*][p], A[*][q]
            for (int r = 0; r < sz; ++r)
            {
                if (r == row || r == col) continue;
                tmp1 = A[r][row] * cos1 - A[r][col] * sin1;
                tmp2 = A[r][row] * sin1 + A[r][col] * cos1;
                A[r][row] = tmp1;
                A[r][col] = tmp2;
            }
            // A[p][*], A[q][*]
            for (int c = 0; c < sz; ++c)
            {
                if (c == row || c == col) continue;
                tmp1 = A[row][c] * cos1 - A[col][c] * sin1;
                tmp2 = A[row][c] * sin1 + A[col][c] * cos1;
                A[row][c] = tmp1;
                A[col][c] = tmp2;
            }
        }

        /// <summary>
        /// 行列の相似変換(右側変換行列の積み上げ)
        /// </summary>
        /// <param name="L">入出力の相似変換行列</param>
        /// <param name="row">変換行</param>
        /// <param name="col">変換列</param>
        /// <param name="sin1">sinθ</param>
        /// <param name="cos1">cosθ</param>
        private static void SimilarityTransformationLeft(float[][] L, int row, int col, float sin1, float cos1)
        {
            int sz = L.Length;
            // for Right Matrix
            //for (int c = 0; c < sz; ++c)
            //{
            //    float d1 = cos1 * R[row][c] - sin1 * R[col][c];
            //    float d2 = sin1 * R[row][c] + cos1 * R[col][c];
            //    R[row][c] = d1;
            //    R[col][c] = d2;
            //}

            // for Left Matrix
            for (int r = 0; r < sz; ++r)
            {
                float d1 = L[r][row] * cos1 - L[r][col] * sin1;
                float d2 = L[r][row] * sin1 + L[r][col] * cos1;
                L[r][row] = d1;
                L[r][col] = d2;
            }
        }

        /// <summary>
        /// 対称行列の非対角成分における絶対値最大の要素の位置を返す
        /// </summary>
        /// <param name="A">対称行列</param>
        /// <param name="row">行の位置</param>
        /// <param name="col">列の位置</param>
        private static void AbsoluteMaxOffDiagonalElementPosition(float[][] A, out int row, out int col)
        {
            int sz = A.Length;
            row = 0;
            col = 1;
            float maxValue = Mathf.Abs(A[row][col]);

            for (int r = 0; r < sz; ++r)
            {
                for (int c = r + 1; c < sz; ++c)
                {
                    float val = Mathf.Abs(A[r][c]);
                    if (maxValue < val)
                    {
                        row = r;
                        col = c;
                        maxValue = val;
                    }
                }
            }
        }

この関数の使い方は次の通り

        int n = 6;
        SimpleMath.CreateMatrix(n, n, out _A);
        int row = 0;
        int col = 0;
        foreach (InputField child in panelA.GetComponentsInChildren<InputField>())
        {
            float random = Random.Range(-3, 3);
            child.text = random.ToString("0.00");
            _A[row][col++] = float.Parse(child.text);
            if (n == col)
            {
                col = 0;
                row++;
                if (n == row)
                {
                    break;
                }
            }
        }
        SimpleMath.CreateMatrix(n, n, out _AtA);
        SimpleMath.GetAutocorrelation(_A, _AtA);

        float[][] L;
        SimpleMath.CreateMatrix(n, n, out L);
        float[][] D;
        SimpleMath.CreateMatrix(n, n, out D);
        SimpleMath.EigenDecomposition(_AtA, D, L, 0.001f, n * n);

動画はギブンス回転の最大回数を n * n から 1 にして、一回転ずつの変換の様子を観察したものです。
だいたい 30~40 回転させると収束するようです。
これで n x n の対称行列の対角化、固有値固有ベクトルが得られるようになりました。

優秀な数値計算ライブラリに頼るのもアリですが、ちょこっとした固有値問題を解くケースで、上記コードを使ってみたいと思いまして作ってみました。
いつでも自分が参照できるように、ここに公開します。

以下、参考文献です。

ヤコビ法のロジックはここの資料を見て書きました。ギブンス回転により必ず対角化に向かうことの証明も書かれています。(大変助かりました、ありがとうございます!)
http://www.na.scitec.kobe-u.ac.jp/~yamamoto/lectures/appliedmathematics2004/chapter11.PDF

相似変換は、ランク、トレース、行列式固有値、固有方程式を保存するという説明向けに紹介したいページです。
そう、ギブンス回転を何回行っても固有値は変わりません。(だから、今回の方法で固有値を求めることができます)
線形代数II/相似変換に対するトレース、行列式、固有値の保存 - 武内@筑波大

ギブンス回転の実装の参考にしたコードが書かれているページです。
式の説明に落とし穴がなく、すばらしい資料になっています。
ただ、ギブンス回転のコードには誤りがあるようで、そのまま利用するとなかなか対角行列に収束しませんでした。
fussy.web.fc2.com

ヤコビ法やハウスホルダー変換が簡潔に説明されている資料です。
上の資料との重複もあり、理解の確認用に使わせてもらいました。
http://www.aics.riken.jp/aicssite/wp-content/uploads/2014/03/sps14_enshu2-1.pdf