python - Pythonの複雑な対称行列

Question

Python で複雑な対称行列を対角化しようとしています。

numpy と scipy の linalg ルーチンを見てみましたが、それらはすべてエルミート行列または実対称行列のいずれかを扱っているようです。

私が探しているのは、最初の複雑な対称行列の高木因数分解を取得する方法です。これは基本的に標準の固有分解 S = VDV^-1 ですが、開始行列 S が対称であるため、結果の V 行列は自動​​的に直交する必要があります。つまり、VT = V^-1 です。

助けはありますか？

ありがとう

Answer 1

以下は、高木因数分解を計算するためのコードです。エルミート行列の固有値分解を使用します。効率的で、フォールトトレラントで、数値的に安定しており、考えられるすべての行列が正しいことを保証するものではありません。この因数分解用に設計されたアルゴリズムは、特に大きな行列を因数分解する必要がある場合に適しています。それでも、いくつかの行列を因数分解して生活を続ける必要がある場合は、このような数学的トリックを使用すると便利です.

import numpy as np
import scipy.linalg as la

def takagi(A) :
    """Extremely simple and inefficient Takagi factorization of a 
    symmetric, complex matrix A.  Here we take this to mean A = U D U^T
    where D is a real, diagonal matrix and U is a unitary matrix.  There
    is no guarantee that it will always work. """
    # Construct a Hermitian matrix.
    H = np.dot(A.T.conj(),A)
    # Calculate the eigenvalue decomposition of the Hermitian matrix.
    # The diagonal matrix in the Takagi factorization is the square
    # root of the eigenvalues of this matrix.
    (lam, u) = la.eigh(H)
    # The "almost" Takagi factorization.  There is a conjugate here
    # so the final form is as given in the doc string.
    T = np.dot(u.T, np.dot(A,u)).conj()
    # T is diagonal but not real.  That is easy to fix by a
    # simple transformation which removes the complex phases
    # from the resulting diagonal matrix.
    c = np.diag(np.exp(-1j*np.angle(np.diag(T))/2))
    U = np.dot(u,c)
    # Now A = np.dot(U, np.dot(np.diag(np.sqrt(lam)),U.T))
    return (np.sqrt(lam), U)

アルゴリズムを理解するには、各ステップを書き出して、目的の因数分解にどのようにつながるかを確認すると便利です。必要に応じて、コードをより効率的にすることができます。