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