python中的复对称矩阵

Question

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

我看过numpy和scipy-linalg例程，但它们似乎都处理hermitian或实对称矩阵。在

我要找的是一些方法来获得我的起始复对称矩阵的Takagi因式分解。这基本上就是标准特征分解 S=vdv^-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)

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