How to find the Takagi decomposition of a symmetric (unitary) matrix?

Question

The Takagi decomposition is a special case of the singular value decomposition for symmetric matrices. More exactly:

Let $U$ be a symmetric matrix, then Takagi tells us there is a unitary $V$ such that $U = VDV^T$ (with $D>0$ diagonal).

My question is basically: how to construct this $V$? Preferably I am looking for the `easiest'/most straight-forward way (which probably won't be the most efficient way!)

Note: For the case I am interested in, $U$ is in fact unitary (in which case Takagi gives $U = VV^T$). I'm happy to specialize to that special case if that makes the algorithm easier.

Answer 1

There is this 2012 paper outlining a strategy to decompose a symmetric unitary matrix $A$ into $A = VV^T$. Interestingly, Ikramov notes that

[For matrices of dimension larger than $5$,] Takagi's decomposition of a general symmetric matrix cannot be obtained by performing a finite number of arithmetic operations and using (a finite number of) root extractions.

Luckily there is such a finite algorithm in case $A$ is unitary, as he outlines in his paper.

I transcribed it to Python and I'm including the code here in case it can be of help to anyone:

import numpy as np

def takagi_for_unitary(A):
    ###################################################
    ### takes a unitary matrix A such that A^T = A, ###
    ###    returns unitary U such that A = U U^T    ###
    ###################################################

    N = A.shape[0]
    U = np.eye(N)

    while A.shape[0] > 0:   
        n = A.shape[0]

        # constructing a vector v such that A*ybar = y
        x = np.random.rand(n)+1j*np.random.rand(n)
        x /= np.sqrt(np.dot(x.conj(),x))
        Axbar = np.tensordot(A,x.conj(),(1,0))
        mu = np.dot(x.conj(),Axbar)

        if np.abs(mu)>1-1e-8: y = mu**.5*x
        else: y = Axbar + x; y /= np.sqrt(np.dot(y.conj(),y))

        # constructing a unitary V with first column = y
        V = np.random.random((n,n))+1j*np.random.random((n,n))
        V[:,0] = y
        V,__ = np.linalg.qr(V)

        # update A
        A = np.tensordot(np.tensordot(V.conj().T,A,(1,0)),V.conj(),(1,0))[1:,1:]

        # update U
        Vi = np.eye(N,dtype=complex)
        Vi[N-n:,N-n:] = V
        U = np.tensordot(U,Vi,(1,0))

    return U

Answer 2

This is a general algorithm for Takagi that reduces to the Schur decomposition with only minimal preprocessing and postprocessing. It's numerically stable.

We seek to obtain the Takagi decomposition of a complex symmetric matrix $M$, i.e. one for which $M = M^T$. This decomposition takes the form $M = U\Sigma U^T$ where $U^{-1}=U^*$ and $\Sigma$ is a diagonal matrix with only non-negative real entries.

Algorithm:

Take the orthogonal eigendecomposition of $M'=\begin{pmatrix}-\Re(M) & \Im(M) \\ \Im(M) & \Re(M)\end{pmatrix}$, and write it as $PDP^T$. This can be obtained using the Schur decomposition of $M'$ since $M'$ is real symmetric. Now consider only the eigenvectors (columns of $P$) where the corresponding eigenvalues (that is, entries along the diagonal of $D$) are positive. Reduce $P$ to $P'$ by keeping only those columns, and obtain $U$ by $\begin{pmatrix}\Im(U) \\\Re(U) \end{pmatrix} = P'$. Define $\Sigma$ to be a diagonal matrix whose entries are the eigenvalues of those eigenvectors of $M'$ that are columns of $P'$, in the same order as they occur in $P'$. Observe that $M = U \Sigma U^T$ where $UU^* = U^* U = I$.

Scipy code:

from numpy import block, diag, real, imag
from scipy.linalg import schur

def takagi(M):
   n = M.shape[0]
   D, P = schur(block([[-real(M),imag(M)],[imag(M),real(M)]]))
   pos_eigenval_positions = diag(D) > 0
   Sigma = diag(D[pos_eigenval_positions,pos_eigenval_positions])
   # Note: The arithmetic below is technically not necessary
   U = P[n:,pos_eigenval_positions] + 1j*P[:n,pos_eigenval_positions]
   return U, Sigma