4

I am trying to diagonalise a complex symmetric matrix in python.

I had a look at numpy and scipy linalg routines but they all seem to deal with either hermitian or real symmetric matrices.

What I am looking for is some way of obtaining the Takagi factorisation of my starting complex and symmetric matrix. This basically is the standard eigendecomposition S = V D V^-1 but, as the starting matrix S is symmetric, the resulting V matrix should automatically be orthogonal, i.e. V.T = V^-1.

any help?

Thanks

purpleshift
  • 170
  • 1
  • 7
  • 4
    A slight clarification: the Takagi factorization is a special case of the singular value decomposition. The values in D are singular values of S, not eigenvalues of S. – Warren Weckesser Oct 21 '13 at 13:05
  • numpy provides two sets of eigenvalue related routines: one for general matrices, another for hermitian matrices (which include real, symmetric matrices). Your matrices are not hermitian, so you cannot take advantage of those special routines, but the general case ones will work with no problems. On the other hand, as @WarrenWeckesser pointed out, the Takagi factorization is not simply the diagonalization of a complex symmetric matrix, and cannot in general be constructed from the eigendecomposition or SVD of the matrix. – Jaime Oct 21 '13 at 13:35
  • Yes, sorry, that depends on the starting matrix being normal. – purpleshift Oct 22 '13 at 08:25

1 Answers1

2

Here is some code for calculating a Takagi factorization. It uses the eigenvalue decomposition of a Hermitian matrix. It is not intended to be efficient, fault tolerant, numerically stable, nor guaranteed correct for all possible matrices. An algorithm designed for this factorization is preferable, particularly if large matrices need to be factored. Even so, if you just need to factor some matrices and get on with your life, then using mathematical tricks such as this can be useful.

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)

To understand the algorithm it is convenient to write out each step and see how it leads to the desired factorization. The code can then be made more efficient if need be.

Craig J Copi
  • 1,754
  • 1
  • 15
  • 12
  • Hi Craig, thanks. I thought of using S^\dg S too, however I need something reliable as I am doing precise simulations. My questions was actually about the existence of dedicated routines for this particular factorisation. For instance, I believe matLab offers such possibility. – purpleshift Oct 22 '13 at 08:30