python中的复杂对称矩阵

Question

我正在尝试在 python 中对 complex 对称矩阵进行对角化。

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

我正在寻找的是某种方法来获得我的起始复杂对称矩阵的高木因式分解。这基本上是标准的特征分解 S = V D V^-1 但是，由于起始矩阵 S 是对称的，因此生成的 V 矩阵应该自动正交，即 V.T = 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)

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