python - python中的复杂对称矩阵

Question

我试图在 python中对角化一个复杂的对称矩阵。

我查看了 numpy 和 scipy linalg 例程，但它们似乎都处理厄米特矩阵或实对称矩阵。

我正在寻找的是某种方法来获得我的起始复杂和对称矩阵的高木分解。这基本上是标准的特征分解 S = VDV^-1 但是，由于起始矩阵 S 是对称的，因此生成的 V 矩阵应该自动正交，即 VT = V^-1。

有什么帮助吗？

谢谢

Answer 1

这是一些用于计算高木分解的代码。它使用 Hermitian 矩阵的特征值分解。它的目的不是高效、容错、数值稳定，也不能保证所有可能的矩阵都是正确的。为这种因式分解而设计的算法更可取，特别是在需要对大型矩阵进行因式分解时。即便如此，如果您只需要分解一些矩阵并继续您的生活，那么使用诸如此类的数学技巧可能会很有用。

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)

为了理解算法，写出每个步骤并查看它如何导致所需的分解是很方便的。如果需要，可以使代码更有效率。