\documentclass[twoside,11pt]{article}
\usepackage{jmlr2e}

\usepackage {graphicx}

\usepackage{amssymb}
\usepackage{amsmath}

\usepackage{float}


\usepackage{algorithmic}
\usepackage{algorithm}


\newcommand{\knote}[1]{{\footnote{K: #1}}}


\newcommand{\MD}{\mathfrak{M}_D} 
\newcommand{\FFD}{\textsc{FFDiag}}
\newcommand{\FFS}{\textsc{FFSEP}}
\newcommand{\rel}{I\!\!R} 
\newcommand{\lbr}[2]{\underset{#1}{\underbrace{#2}}}
\renewcommand{\d}[2]{\ensuremath{\frac{\partial #1}{\partial #2}}}
\newcommand{\for}{\hspace{5mm}{\rm for}\hspace{2mm}}
\newcommand{\jac}{\mathsf{D}}
\newcommand{\ddef}{\overset{\text{def}}{=}}

\DeclareMathOperator{\opvec}{vec}
\DeclareMathOperator{\opvecs}{vecs}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\dg}{diag}
\newcommand{\RR}{\mathbb{R}}
\newcommand{\diff}{\text{d}}
\newcommand{\s}[1]{\ensuremath{\mathcal{#1}}}%
\newcommand{\x}{{\mathbf x}}

\newcommand{\id}{I}
\newcommand{\dett}{\mbox{det}}
\newcommand{\nin}{\hspace{2mm}/\!\!\!\!\!\in} 

\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\bea}{\begin{eqnarray}} 
\newcommand{\eea}{\end{eqnarray}}
\newcommand{\off}[1]{\text{off}(#1)}
\newcommand{\grad}{\nabla}

\newcommand{\pe}{\psi}
\def\d{\delta} 
\def\ds{\displaystyle} 
\def\enorm#1{\|#1\|_2} 
\def\Fp{F^\prime}  
\def\fishpack{{FISHPACK}} 
\def\fortran{{FORTRAN}} 
\def\MATALB{{MATLAB}} 
\def\norm#1{\|#1\|} 
\def\wb{{\bar w}} 
\def\zb{{\bar z}} 




\title{A Fast Algorithm for Joint Diagonalization with Non-orthogonal
  Transformations and its Application to\\ Blind Source Separation}

\author{%
  \name Andreas Ziehe
   \email ziehe@first.fhg.de \\
          \name Pavel Laskov      
  \email laskov@first.fhg.de \\
  \addr Fraunhofer FIRST.IDA\\ 
Kekul\'estrasse 7\\
  12489 Berlin, Germany
  \AND
  \name Guido Nolte 
  \email NolteG@ninds.nih.gov \\
  \addr
  National Institutes of Health \\
  10 Center Drive MSC 1428\\
  Bethesda, MD 20892, USA
  \AND
  \name Klaus-Robert M\"uller
  \email klaus@first.fhg.de  \\
  \addr Fraunhofer FIRST.IDA\\ 
Kekul\'estrasse 7\\
  12489 Berlin, Germany\\ and \\
  \addr  University of Potsdam, Department of Computer
  Science\\
  August-Bebel-Strasse 89\\
 14482 Potsdam, Germany
}

%submitted{12/2003}; revised{05/2004}; published{??/2004}


\jmlrheading{5}{2004}{777--800}{12/03; Revised 5/04}{7/04}{Andreas Ziehe, Pavel Laskov, Guido  Nolte, and  Klaus-Robert M\"uller} 
\ShortHeadings{Fast Non-orthogonal Joint Diagonalization}{Ziehe et al.}
\editor{Michael Jordan}
\firstpageno{777}
\begin{document}



\maketitle


\begin{abstract}
A new efficient algorithm is presented for joint diagonalization of
several matrices. The algorithm is based on the Frobenius-norm
  formulation of the joint diagonalization problem, and addresses
  diagonalization with a general, non-orthogonal transformation. The
  iterative scheme of the algorithm is based on a multiplicative
  update which ensures the invertibility of the diagonalizer.  The
  algorithm's efficiency stems from the special approximation of the
  cost function resulting in a sparse, block-diagonal Hessian to be
  used in the computation of the quasi-Newton update step.  Extensive
  numerical simulations illustrate the performance of the algorithm
  and provide a comparison to other leading diagonalization methods.
  The results of such comparison demonstrate that the proposed
  algorithm is a viable alternative to existing state-of-the-art joint
  diagonalization algorithms.  The practical use of our algorithm is
  shown for blind source separation problems.
\end{abstract}

\begin{keywords} 
  joint diagonalization, common principle component analysis,
  independent component analysis, blind source separation,
  nonlinear least squares, Newton method, Levenberg-Marquardt algorithm
\end{keywords}



\section{Introduction}

Joint diagonalization of square matrices is an important problem of
numeric computation. Many applications make use of joint
diagonalization techniques as their main algorithmic tool, for example
independent component analysis (ICA) and blind source separation (BSS)
\citep{Comon1994, MolgedeySchuster1994, Belo97, ica99:wu, JADE:NC,
  Ziehe98,ZieLasMueNol03, Cardoso00b, ZiTBME2000, Yer02, Hay00,
  HyKaOj01}, common principal component analysis (CPC)
\citep{Flury1988,AirFlu88,FenHaeVil01}, various signal processing
applications \citep{Veen92,Veen98} and, more recently, kernel-based
nonlinear BSS \citep{HarZieKawMue03}.

This paper pursues two goals. First, we propose a new efficient
algorithm for joint approximate matrix diagonalization.  Our algorithm
is based on the second-order approximation of a cost function for the
simultaneous diagonalization problem. 
Second, we demonstrate an application of our algorithm to BSS, which
allows to perform BSS without pre-whitening the data.

Let us begin by defining the notion of joint diagonalization.  It
is well known that exact joint diagonalization is in general only
possible for two matrices and amounts to the generalized eigenvalue
problem.  Extensive literature exists on this topic
\citep[e.g.][and references therein]{NobDan77,Golub:1989,BunseGerstner93, vandervorst97years}. When more than two matrices are to be
diagonalized, exact diagonalization may also be possible if the
matrices possess a certain common structure. Otherwise one can only speak
of approximate joint diagonalization. Our paper focuses on the
investigation of  algorithms for exact---whenever this is
possible---or otherwise approximate diagonalization of more than two
matrices. In the remainder of the paper we will refer to such problems
as  ``joint diagonalization'' problems.

A number of algorithms for joint diagonalization have been previously
proposed in the literature
\citep{Flury:1986:ASO,Car93,SC-siam,Hori99,Pham2001c,Vee01,Yer02,Joho2002}.
To understand the challenges of the joint diagonalization problem, as
well as the need for further improvement of currently known
algorithms and possible directions of such improvement some insight
into the main issues of joint diagonalization is now provided.
%

Let us consider a set $\{C^1, \ldots, C^K\}$ of real-valued symmetric
matrices of size $N
\times N$.\footnote{The formulations and the proposed algorithm will be
  presented for symmetric matrices. Extensions to the unsymmetric or
  complex-valued case can be obtained in a similar manner.} 
The goal of a joint diagonalization algorithm is to find a
transformation $V$ that in some sense ``diagonalizes'' all the given
matrices. The notion of diagonality and the corresponding formal
statement of the joint diagonalization problem can be defined in at
least three different ways:
%
\begin{enumerate}
\item \emph{Frobenius norm formulation}. This formulation is used in
  \citet{Car93,SC-siam,Joho2002} and, in a generalized form, in \citet{Hori99}. Let 
  \begin{equation}
    \label{eq:Fk}
    F^k = V C^k V^T
  \end{equation}
  denote the result of applying transformation $V$ to matrix
  $C^k$. Joint diagonalization is defined as the following
  optimization problem:
  \begin{equation}
    \min_{V \in \rel^{N \times N}} \sum\limits_{k=1}^{K} \MD(F^k),
    \label{gen_opt2}
  \end{equation}
  where the diagonality measure $\MD$ is the Frobenius norm of the
  off-diagonal elements in $F^k$:\ 
  \be
  \MD(F^k)= \off{F^k} = \sum_{i\neq j} (F^k_{ij})^2. 
  \label{sum_square_diag} 
  \ee
%
  A more careful look at the cost function in Equation (\ref{gen_opt2}) reveals a serious problem with the Frobenius
  norm formulation: this cost function is obviously minimized by the
  trivial solution $V = 0$.  The problem can be avoided by
  additionally requiring orthogonality of $V$. In fact, this
  assumption is very natural if the joint diagonalization problem is
  seen as an extension of the eigenvalue problem to several matrices.
  However, restricting $V$ to the group of orthogonal matrices
  may limit the applicability and unduly degrade the performance of the method.\footnote{
  In principle, a pre-sphering step could be applied to alleviate this problem,
nevertheless a performance degradation is to be expected in this case, especially in the context of blind source separation \citep{Car94,Yer02}.}
\smallskip
\item \emph{Positive definite formulation}. Another reasonable
  assumption on the initial problem is the positive-definiteness of
  the matrices $C^k$. This assumption is motivated by the fact that in
  many applications matrices $C^k$ are covariance matrices of some
  random variables. In this case, as proposed in
  \citet{Matsuoka95,Pham2001c}\, the criterion\footnote{Here the operator $\text{ddiag}(F)$ returns a diagonal matrix containing only the diagonal entries of $F$.}
  \begin{equation}
    \MD(F^k)=  \mbox{log det(ddiag}(F^k))- \mbox{log det}(F^k)
    \label{pham_diag_crit}
  \end{equation}
  can be used in the cost function (\ref{gen_opt2}) instead of the
  criterion (\ref{sum_square_diag}). The additional advantage of this
  criterion is that it allows for super-efficient estimation
  \citep{Cardoso01}. However in certain applications, such as blind
  source separation based on time-delayed decorrelation
  \citep{Belo97,Ziehe98}, correlation matrices are no longer guaranteed
  to be positive-definite, and diagonalization based on this
  criterion may fail.
\smallskip
\item \emph{Subspace fitting formulation}. The fact that exact joint
  diagonalization may not be possible can be explicitly accounted for in the
  problem formulation. This is to say that, instead of applying the
  transformation directly to the matrices $C^k$, another set of diagonal
  matrices $\Lambda^k$ is sought for, along with the transformation so as
  to best approximate the target matrices. The optimization problem
  resulting from this approach
  \begin{equation}
    \label{eq:subspace}
    \min_{A,\Lambda^1, \ldots, \Lambda^K} \sum_{k=1}^K||C^k - A \Lambda^k A^T||_F^2
  \end{equation}
  constitutes an instance of a subspace fitting problem  \citep{Vee01,Yer02}. 
  
  Compared to the previous approaches, the algorithms based on
  subspace fitting have two advantages: they do not require
  orthogonality, positive-definiteness or any other normalizing assumptions on the
  matrices $A$ and $C^k$, and they are able to handle non-square
  mixture matrices.  These advantages, however, come at a high
  computational cost: the algorithm of \cite{Vee01} has quadratic
  convergence in the vicinity of the minimum but its running time per
  iteration is $O(KN^6)$, whereas the AC-DC algorithm of \cite{Yer02}
  converges linearly with a running time per iteration of order
  $O(KN^3)$.
\end{enumerate}
%
As a short resume of the above mentioned approaches we notice the
following. The algorithms using the Frobenius norm formulation are
efficient but rely on the orthogonality assumption to prevent
convergence to the trivial solution. The algorithms using the
positive-definiteness assumption are also quite efficient but they may
fail if this assumption is not satisfied. 
Subspace fitting algorithms, which do not require such strong prior
assumptions, are computationally much more  demanding.
%
A natural question arises: could a single algorithm combine the
positive and avoid the negative features of the previous joint
diagonalization algorithms? In this contribution we present an
algorithm using the Frobenius norm formulation that strives towards
this goal. In particular, the algorithm, to be called {\FFD} (Fast
Frobenius Diagonalization), possesses the following features:
\begin{itemize}
\item computational efficiency: quadratic convergence (in the
  neighborhood of the solution) and $O(KN^2)$ running time per
  iteration,
\item guaranteed avoidance of the trivial solution,
\item no orthogonality and no positive-definiteness assumptions;
  nonetheless, orthogonality can be used to constrain the
  solution, which further reduces the computational complexity by a factor of two.
\end{itemize}
On top of that, the algorithm is  simple and easy to implement. 
\smallskip

The remainder of the paper is organized as follows. In Section
\ref{sec2} the main idea of the {\FFD} algorithm is
proposed. The computational details regarding the algorithm's update rule  are derived in
Section \ref{SparseLSProblem}. Section \ref{sec:lmalgo}, in a slight
digression from the main topic of the article, presents a connection
of our algorithm to the classical Levenberg-Marquardt
algorithm, and points out the main differences between the two. The
application of the {\FFD} algorithm to blind source separation is
developed in Section \ref{sec:bss}. Extensive numerical simulations
are presented in Section \ref{sec:sim}. Finally, Section
\ref{sec:conclu} is devoted to the discussion and conclusions.


\section{General Structure of the Algorithm} 
\label{sec2}

The  {\FFD} algorithm is an iterative scheme to approximate the solution of the following
optimization problem:

% insert cost function here
\begin{equation}
    \min_{V \in \rel^{N \times N}} \sum\limits_{k=1}^{K} \sum_{i\neq j} ((V C^k V^T)_{ij})^2.
    \label{gen_opt}
  \end{equation}

The basic  idea is to use the invertibility of
the matrix $V$ as a constraint preventing convergence of the
minimizer of the cost function in Equation (\ref{gen_opt}) to the zero
solution. Invertibility is tacitly assumed in many applications of
diagonalization algorithms, e.g.~in blind source separation, therefore
making use of such constraint is very natural and does not 
 limit the generality from the practical point of view.

Invertibility can be enforced by carrying out the 
update of $V$ in multiplicative form as
\begin{equation}
  \label{eq:update}
  V_{(n+1)} \leftarrow (I+W_{(n)})\, V_{(n)},
\end{equation}
where $I$ denotes the identity matrix, the update matrix $W_{(n)}$ is
constrained to have zeros on the main diagonal, and $n$ is the
iteration number.  Such update is rarely used in classical
unconstrained optimization algorithms; however, it is common  
for many successful BSS algorithms, such as relative-gradient \citep{LahCar96,AmariNC00a},
relative Newton \citep{AkuMura2001,Zib03}, as well as for joint diagonalization  
\citep{Pham2001c}. The off-diagonal component
$W_{(n)}$ of the update multiplier is to be determined so as to
minimize the cost function (\ref{gen_opt}). In order to maintain
invertibility of $V$ it clearly suffices to enforce invertibility of
$I+W_{(n)}$. The latter can be carried out using the following
well-known results of matrix analysis \citep{HJ85}.

\begin{definition}
  An $n \times n$ matrix $A$ is said to be strictly diagonally dominant
  if $$ |a_{ii}| > \sum_{j \neq i} |a_{ij}|, \quad\text{\emph{for all} $i = 1,\ldots,n.$}$$
\end{definition}

\begin{theorem}[Levi-Desplanques]
  If an $n \times n$ matrix $A$ is strictly diagonally-dominant, then it is invertible.
\end{theorem}

The Levi-Desplanques theorem can be used to control invertibility of
$I+W_{(n)}$ in a straightforward way. Observe that the diagonal
entries in  $I+W_{(n)}$  are all equal to 1; therefore, it suffices to
ensure that $$\max_{i}{\sum_{j \neq i} |W_{ij}|} = ||W_{(n)}||_{\infty} <
1.$$ This can be done by dividing $W_{(n)}$ by its infinity norm
whenever the latter exceeds some fixed $\theta < 1$. An even stricter
condition can be imposed by using a Frobenius norm in the same way:
\begin{equation}
  W_{(n)} \leftarrow \frac{\theta} {||W_{(n)}||_F} W_{(n)}.
  \label{eq:heuristic}
\end{equation}

To determine the optimal updates $W_{(n)}$ at each iteration,
first-order optimality constraints for the objective 
(\ref{gen_opt}) are used. A special approximation of the objective
function will enable us to efficiently compute $W_{(n)}$.
%
For this reason, we consider the expression for updating the matrices to be diagonalized
\begin{equation}
  \label{eq:target}
  C^k_{(n+1)} \leftarrow F^k = (I+W_{(n)}) \,\, C^k_{(n)} \, (I+W_{(n)})^T.
\end{equation}

Let $D^k_{(n)}$ and $E^k_{(n)}$ denote the diagonal and off-diagonal
parts of $C^k_{(n)}$, respectively. In order to simplify the
optimization problem we assume that the norms of $W_{(n)}$ and $E^k_{(n)}$ are
small, i.e.~quadratic terms in the expression for the new
set of matrices can be ignored
\begin{align}
C^k_{(n+1)} &= (I+W_{(n)}) (D^k_{(n)}+E^k_{(n)}) (I+W_{(n)})^T\nonumber\\
  &\approx D^k_{(n)}+W_{(n)} D^k_{(n)} + D^k_{(n)} W_{(n)}^T +
  E^k_{(n)}. \label{eq:simplify}
\end{align}
With these simplifications,
 and ignoring already diagonal terms $D^k$, the
diagonality measure (\ref{sum_square_diag}) can be
computed using expressions linear in $W$.\footnote{The iteration indices
  will be dropped in the following if all quantities refer to the
  same iteration.}

 \be F^k \approx
\tilde{F}^k = WD^k+D^kW^T+E^k.
\label{eq:approx}
\ee
The linearity of terms (\ref{eq:approx}) allows us to explicitly
compute the optimal update matrix $W_{(n)}$ minimizing the approximated
diagonality criterion
\begin{equation}
  \label{eq:min_update}
  \min_W \sum_{k=1}^K  \sum_{i\neq j} {((WD^k+D^kW^T+E^k)_{ij})^2}.
\end{equation}
Details of the efficient solution of problem (\ref{eq:min_update}) are
presented in Section~\ref{SparseLSProblem}.

The simplifying assumptions used in (\ref{eq:simplify}) require some
further discussion. The motivation behind them is that in the
neighborhood of the optimal solution, the optimal steps $W$ that take
us to the optimum are small and the matrices $C^k$ are almost
diagonal. Hence, in the neighborhood of the optimal solution the
algorithm is expected to behave similarly to Newton's algorithm and
converge quadratically. The assumption of small $E^k$ is potentially
problematic, especially in the case where exact diagonalization is
impossible. A similar derivation can be carried out with $E^k$ fully
accounted for, which leads to additional $WE^k$ and $E^k W^T$ terms in
the expression~(\ref{eq:min_update}). However, the resulting
algorithm, will not give rise to a computationally efficient solution
of problem (\ref{eq:min_update}).  As for the assumption of small $W$,
it is crucial for the convergence of the algorithm and needs to be
carefully controlled. The latter is done by the
normalization~(\ref{eq:heuristic}).

Pseudo-code describing the {\FFD} method is outlined in Algorithm~\ref{algo:GNdiag}.\footnote{MATLAB code for {\FFD} can be obtained at {\tt http://www.first.fhg.de/\~{}ziehe/research/FFDiag}.}
\begin{algorithm}
\begin{algorithmic} 
\STATE 
\STATE {INPUT: $C^k$} \COMMENT{ Matrices to be diagonalized}
\STATE  $W_{(1)} \leftarrow 0$, $V_{(1)} \leftarrow I$, $n \leftarrow 
1$ \COMMENT{ $V_{(1)}$ could also be initialized by a more clever guess.}
\STATE $C^k_{(1)} \leftarrow V_{(1)} \,\, C^k \, V_{(1)}^T$
\smallskip
\REPEAT
 \STATE compute $W_{(n)}$ from $C_{(n)}^k$ according to Equation (\ref{Wup}) or (\ref{W_ortho})
\smallskip
\IF{$||W_{(n)}||_F > \theta$}
\STATE $W_{(n)} \leftarrow \frac{\theta} {||W_{(n)}||_F} W_{(n)}$\ENDIF 
\smallskip
\STATE $V_{(n+1)} \leftarrow (I+W_{(n)})\, V_{(n)}$
\smallskip
\STATE $C^k_{(n+1)} \leftarrow (I+W_{(n)}) \,\, C^k_{(n)}\, (I+W_{(n)})^T$
\smallskip
\STATE $n \leftarrow n+1$
\UNTIL convergence
\STATE OUTPUT: $V , C^k$
\end{algorithmic}
\caption{{\FFD} \label{algo:GNdiag}}
\label{algo:main}
\end{algorithm}

Some remarks on convergence properties of the proposed algorithmic
scheme are due at this point.  In general, Newton-like algorithms are
known to converge only in the neighborhood of the optimal solution;
however, when they converge, the rate of convergence is quadratic
\citep[e.g.][]{Kan49}. Since the essential components of our algorithm---the second-order approximation of the objective function and the
computation of optimal steps by solving the linear system arising from
first-order optimality conditions---are inherited from Newton's
method, the same convergence behavior can be expected.  In practice,
however, the known theoretical estimates of convergence regions of
Newton's method, such as the ones provided, e.g.,~in Theorems 1 and 2
in \cite{Kan49}, are of little utility since they provide no guidance
how to reach the convergence region from an arbitrary starting point.


\section{Computation of the Update Matrix}
\label{SparseLSProblem}

The key to computational efficiency of the {\FFD} algorithm lies in 
exploiting the sparseness introduced by the approximation (\ref{eq:approx}).
The special structure of the problem can be best seen in the matrix-vector
notation presented next.

The  $N(N-1)$ off-diagonal entries of the update matrix $W$ are arranged as
\be
w=(W_{12},W_{21},\ldots,W_{ij},W_{ji},\ldots)^T.
\ee
Notice that this is \emph{not} the usual vectorization operation $\opvec W$,
as the order of elements in $w$ reflects the pairwise relationship of the
elements in $W$.
In a similar way the $KN(N-1)$  off-diagonal entries of the matrices $E^k$ are
arranged as
\be
e=(E^1_{12},E^1_{21},\ldots,E^1_{ij},E^1_{ji},\ldots,E^k_{ij},E^k_{ji},\ldots)^T.
\ee
Finally, a large but very sparse,  $KN(N-1)\times N(N-1)$ matrix $J$ is built, in
the form:
%
$$
J=\left(
\begin{array}{c}
J_1\\
\vdots\\
J_K
\end{array}\right)
\hspace{1.5mm}\mbox{with}\hspace{1.5mm} 
J_k=\left(
\begin{array}{cccc}
\s{D}^k_{12} &        &          &         \\
         & \ddots &          &         \\
         &        & \s{D}^k_{ij} &         \\
         &        &          &  \ddots 
\end{array}
\right),
$$
where each $J_k$ is block-diagonal, containing $N(N-1)/2$  matrices of dimension $2\times 2$  
$$
%
%
\s{D}^k_{ij}=\left(
\begin{array}{cc}
D^k_{j} &D^k_{i} \\
D^k_{j} &D^k_{i}
  \end{array}\right), \quad i, j = 1,\ldots,N,\; i \neq j,
%
$$
where $D^k_{i}$ is a short-hand notation for the $ii$-th entry of a diagonal matrix $D^k$.
Now the approximate cost function can be re-written as the linear least-squares problem 
\be
\s{L}(w)=\sum_k \sum_{i\neq j} (\tilde{F}_{ij}^k)^2 = (Jw+e)^T (Jw+e). \nonumber%=(Jw+e)^T (Jw+e). \nonumber
\ee
The well-known solution of this problem \citep{NR} reads 
\be
w=-(J^TJ)^{-1} J^T e.
\label{eq:sol}
\ee
%
We can now make use of the sparseness of $J$ and $e$ to enable the direct
computation of the elements of $w$ in (\ref{eq:sol}). Writing out the matrix
product $J^T J$ yields a block-diagonal matrix
$$
J^T J=
\left(\begin{array}{cccc}
\sum_k (\s{D}^k_{12})^T \s{D}^k_{12}&         &                     &         \\
                    &  \ddots &                     &         \\
                    &         & \sum_k (\s{D}^k_{ij})^T  \s{D}^k_{ij} &         \\
                    &         &                     & \ddots 
\end{array}
\right)
$$
whose blocks are $2 \times 2$ matrices. 
Thus the system (\ref{eq:sol}) actually
consists of decoupled equations
\begin{equation*}
\begin{pmatrix} W_{ij} \\
                W_{ji} 
\end{pmatrix}\,\, = \, 
-\begin{pmatrix} z_{jj} & z_{ij} \\
                z_{ij} & z_{ii}
\end{pmatrix}^{-1}
 \begin{pmatrix} y_{ij} \\
                 y_{ji}
 \end{pmatrix}, \quad i, j = 1,\ldots,N,\; i \neq j,
\end{equation*}
where
\begin{equation*}
  \begin{aligned}
    z_{ij} &= \sum_k D_i^k D_j^k\\
    y_{ij} &= \sum_k D_j^k \frac{E_{ij}^k+E_{ji}^k}{2} = \sum_k D_j^k E_{ij}^k.
  \end{aligned}
\end{equation*}
The matrix inverse can be computed  in closed form,
leading to  the following expressions for the
update of the entries of $W$:
\begin{equation}
\begin{aligned}
W_{ij}&=\frac{z_{ij}y_{ji} - z_{ii} y_{ij}}{z_{jj} z_{ii} - z_{ij}^2 } \\
W_{ji}&=\frac{z_{ij}y_{ij} - z_{jj} y_{ji}}{z_{jj} z_{ii} - z_{ij}^2 }.
\end{aligned}
\label{Wup}
\end{equation}
(Here, only the off-diagonal elements ($i\neq j$) need to be computed and the diagonal terms of $W$ are set to zero.)
Thus, instead of performing inversion and multiplication of large
matrices, which would have brought us to the same $O(KN^6)$ complexity
as in \citet{Vee01}, computation of the optimal $W_{(n)}$ leads to a
simple formula~(\ref{Wup}) which has to be evaluated for each of $N(N-1)$
components of $W_{(n)}$. Since the computation of $z_{ij}$ and
$y_{ij}$ also involves a loop over $K$, the overall complexity of the
update step is $O(KN^2)$. 

An even simpler solution can be obtained if the diagonalization matrix
$V$ is assumed to be orthogonal from the very beginning.
Orthogonality of $V$ can be preserved to the first order by requiring  $W$ to be skew-symmetric, i.e., $W =-W^T$. 
 Hence only one of each pair of its entries needs to be
computed.
 In this case the structure of the problem is already
apparent in the scalar notation, and one can easily obtain the partial
derivatives of the cost function. Equating the latter to zero yields
the following expression for the update of $W$:  
\begin{equation}
  W_{ij} = \frac{\sum_k E_{ij}^k (D_i^k - D_j^k)} {\sum_k (D_i^k - D_j^k)^2}, \quad i, j = 1,\ldots,N,\; i \neq j,
\label{W_ortho}
\end{equation}
which agrees with the result of \cite{PertDJ}.
%
To ensure orthogonality of $V$ beyond the first order the update (\ref{eq:update})
should be replaced by the matrix exponential update 
 \begin{equation}
  \label{eq:exp_update}
  V_{(n+1)} \leftarrow \exp(W_{(n)}) V_{(n)},
\nonumber
\end{equation}
where $W_{(n)}$ is skew-symmetric \cite[cf.][]{AkuMura2001}.

\section{Comparison with the Levenberg-Marquardt Algorithm}
\label{sec:lmalgo}

The Leven\-berg-Marquardt (LM) algorithm \citep{Lev44, Mar63} is one
of the powerful and popular algorithms for solving nonlinear
least-squares problems. Interestingly, the motivation in the original
article of \cite{Mar63} was somewhat similar to ours: he knew that
quadratic convergence of Newton's method was attainable only in the
neighborhood of the solution, and therefore he looked for an efficient
means of steering the algorithm to the area of quadratic convergence.
The particular mechanism used in the LM algorithm consists of a
controllable trade-off between Newton and gradient steps.

Although the problem of simultaneous diagonalization is essentially a
nonlinear (quadra\-tic) least-squares problem, the LM algorithm cannot
be directly applied to it. An implicit assumption of the simultaneous
diagonalization problem is the invertibility of the diagonalizing matrix,
and the classical LM algorithm does not provide for incorporation of
additional constraints. In what follows we present a modification
which allows one to incorporate the additional structure of our
problem into the LM algorithm.

A general problem of the nonlinear regression to be solved by the LM
algorithm is usually formulated as
\begin{equation*}
  \min_p \, \sum_k ||y_k - f_p(x_k)||^2.
\end{equation*}
The goal of the optimization is to find the parameters $p$ of the
regression function $f$ so as to minimize the squared deviations
between $f_p(x_k)$ and $y_k$ for all data points $1, \ldots, k$.
The signature of the function $f$ can be arbitrary, with an appropriate norm
to be chosen.
The cost function~(\ref{gen_opt}) can be seen as a
nonlinear regression problem over the vector-valued function $f_V(C)$,
 parameterized
by $V$, of matrix argument $C$, with zero target values:
\begin{equation*}
  \min_{V \in \rel^{n \times n}} \, \sum_k ||0 - f_V(C^k)||^2.
\end{equation*}
The construction and dimensionality of the function $f$ is explained
below, and the zero vector is set to have the appropriate dimension. 

To enforce invertibility, the diagonal entries of $V$ are set to 1 and
only off-diagonal entries are considered as free parameters. As in Section \ref{SparseLSProblem} such a
representation of $V$ is constructed  by means of symmetric
vectorization  $\opvecs V \ddef [V_{21}, V_{12}, V_{31}, V_{13},
\ldots]^T$.\footnote{The symmetric vectorization $\opvecs$ is related
  to column vectorization $\opvec$ by the special permutation matrix
  $S_{NN}$ such that $\opvecs X = S_{NN} \opvec X$.} 

The same vectorization is applied to construct the regression function
$$f_V(C) : \rel^{N \times N} \rightarrow \rel^{N(N-1) \times 1} \ddef
\opvecs V C V^T.$$
As a result of such vectorization, the diagonal entries of $V C V^T$
are discarded, and the Euclidean norm of this vector is equivalent to the
``off'' function.

The LM algorithm requires computation of the Jacobian
matrix of the regression function (w.r.t.~parameters $\opvecs\! V $) at all data points:
  \begin{equation*}
    J_{LM} \ddef [\jac_{\opvecs\! V} f_V(C^1),\ldots,\jac_{\opvecs\! V} f_V(C^K)]^T.
  \end{equation*}
The Jacobian matrices (of dimension $N(N-1)\times N(N-1)$) at
individual data points can be computed as:
\begin{equation*}
  \jac_{\opvecs\! V} f_V(C^k) = S_{NN} (I_{N^2} + K_{NN}) (V C^k \otimes
  I_N) S_{NN}^T,
\end{equation*}
where $K_{NN}$ is the commutation matrix \citep{Lue96}, $I$ is the
identity matrix of the appropriate size, and $\otimes$ is the
Kronecker matrix product. 

Denoting $f =
[f^T(C^1),\ldots,f^T(C^K)]^T$, the main step of the LM algorithm consists of
solving the following linear system:
\begin{equation}
  \label{eq:lmmain}
  ((J_{LM})^T J_{LM} + \lambda I) \opvecs V = -(J_{LM})^T f.
\end{equation}
The parameter $\lambda$ controls the trade-off between Newton-like
and gradient-based strategies: $\lambda = 0$ results in the pure
Newton-direction, whereas with large $\lambda$, the steps
 approach the gradient direction. We use the original heuristic of
Marquardt to choose $\lambda$: if the value of the cost function
provided by the current step $\opvecs V$ decreases, $\lambda$ can be
decreased by a factor of 10 while descent is maintained; if the value of the cost
function increases, increase $\lambda$ by a factor of 10 until descent is
achieved. This heuristic is very intuitive and easy to implement;
however, since it doesn't involve any line search, the algorithm may
fail to converge. More sophisticated strategies with convergence
guarantees of \cite{Osb76} and \cite{Mor78} can also be deployed.

From the theoretical point of view, one can draw the following
parallels between the {\FFD} and LM algorithms:
\begin{itemize}
\item Both algorithms pursue a Newton direction (cf.~equations
  (\ref{eq:sol}) and (\ref{eq:lmmain})) to compute the update matrix
  $V$. Whereas the LM algorithms computes the update step directly on
  $V$, the update of the {\FFD} is performed in a multiplicative way
  by computing $W$ to be used in the update rule (\ref{eq:update}).
\item Unlike the LM algorithm using the Hessian of the original cost
  function, the Newton direction in {\FFD} is computed based on the  Hessian of
  the second-order approximation of the cost function.\footnote{In fact, in both algorithms, the Hessians are approximated by the product of the Jacobian matrices.}  Taking advantage of
  the resulting special structure, this computation can be carried out
  very efficiently in {\FFD}.
\item Regularization in the LM algorithm results in a gradual shift from the
  gradient to the Newton directions (and back when
  necessary). Regularization in the {\FFD} algorithm is of quite
  different flavor. Since the computed direction is only approximately
  Newton, one cannot fully trust it, and therefore the update
  heuristic (\ref{eq:heuristic}) limits the impact of inaccurate
  computation of $W$. On the other hand, when {\FFD} converges
  to a close neighborhood of the optimal solution, the heuristic is
  turned off, and Newton-like convergence is no longer impeded.
\end{itemize}

It is interesting to compare performance of the LM and {\FFD}
algorithms experimentally. We use two criteria: the cost function and
the convergence ratio
\begin{equation*}
  \text{convergence ratio} = \frac{||f_{(n+1)} - f^{*}||}{||f_{(n)} - f^{*}||}.
\end{equation*}

The zero value of the convergence ratio indicates super-linear
convergence. The evolution of our criteria in two runs of the
algorithms are shown in Figure \ref{fig:lmfig}. In the cost function
plot one can see that convergence of both algorithms is linear for the
most part of their operation, with a gradual shift to quadratic
convergence as the optimal solution is approached. The same conclusion
can be drawn from the convergence ratio plot, in which one can see
that this criterion approaches zero in the neighborhood of the optimal
solution. Thus one can conclude that, similarly to the LM algorithm,
the heuristic~(\ref{eq:heuristic}) steers the algorithm to the area
where Newton-like convergence is achieved.  Furthermore, we note that due to the use
of the special structure, the per-iteration complexity of the {\FFD}
algorithm is significantly lower than that of the LM algorithm.

\begin{figure}[tbp]
  \centering
  \includegraphics[width=0.7\linewidth]{expt-lm}
  \caption{Comparison of the LM and the {\FFD} algorithms. The data
  matrices are generated as described in Section \ref{sec:sanity} with $K =
  10$, $N = 10$. Two illustrative runs are shown.}
  \label{fig:lmfig}
\end{figure}


\section{Application to Blind Source Separation}
\label{sec:bss}

First, we recall the definition of the blind source separation (BSS) problem \citep{Jutten+Herault1991}.
We are given  the instantaneous linear mixtures
 $x_i(t)$ of a number of source signals $s_j(t)$,
obeying the model
 \begin{equation}
    x_i(t)=\sum\limits_{j=1}^m a_{ij} s_j(t),\,\,\, (i = 1,\ldots,n,
\, j=1,\ldots,m),
      \label{icamodel}
 \end{equation}
 with $A$ being non-singular and $s_j(t)$  statistically
 independent. The goal is to estimate both $A$ and $s(t)$ from $x(t)$. 
 
 Linear BSS methods have been successfully applied to a variety
 of problems. An example of such an application is the reduction of
 artifacts in electroencephalographic (EEG) and magnetoencephalographic
 (MEG) measurements.  Due to the fact that the electromagnetic waves
 superimpose linearly and virtually instantaneously (because of the
 relatively small distance from sources to sensors) the model
 (\ref{icamodel}) is valid
 \citep{nips-8:Makeig+Bell+Jung:1996,nips-10:Vigario+Jousmaki+Hamalainen:1998,WueZieMacMueTraCur00,ZiTBME2000}.
 Note, however, that in other applications, such as the so called
 cocktail-party problem in auditory perception \citep{masc86}, the model
 from Equation (\ref{icamodel}) may be too simplistic, since time-delays in
 the signal propagation are no longer negligible. Extended models to
 deal with such convolutive mixtures have been considered
 \citep[e.g.][]{Parra98,icassp98:Lee+Ziehe+Orglmeister,MuIkeZie2001}.
 We will in the following only discuss how to solve the linear,
 instantaneous BSS problem stated in Equation (\ref{icamodel}).
% 
 The usual approach is to define an appropriate cost function that can
 subsequently be optimized.  Here our goal is to use the general joint
 diagonalization cost function (\ref{gen_opt}) and to construct
 certain matrices in such a way that their approximate joint diagonalizer
 is an estimate for the demixing matrix $V$ (up to an arbitrary
 permutation and scaling of its rows).
%
 
Let us consider for example the spatial covariance matrix of the mixed
signals $x(t)$,
\begin{equation}
\label{eq:cov}
C_{(x)} \ddef E\{x(t) x(t)^T\} = E\{(As(t))(As(t))^T\} = A \,E\{s(t) s(t)^T \}\, A^T,  \nonumber
\end{equation}
where the expectation is taken over $t$.
We see that theoretically  $C_{(x)}=A \, C_{(s)} \, A^T$ is similar to a diagonal matrix,
because the cross-correlation terms that form the off-diagonal part of
$ C_{(s)}$ are zero for independent signals.
%
%
There are many more possibilities to define matrices that
 have the same property as the covariance matrix, namely that
they are diagonal for the source signals and `similar to diagonal'
for the observed mixtures and, most important, that the inverse $V=A^{-1}$ of the
 mixing matrix $A$ diagonalizes them all simultaneously.
 Examples are time-lagged covariances
 \citep{MolgedeySchuster1994,Belo97,Ziehe98}, covariance matrices of
 different segments of the data \citep{Cardoso00b,Choi01}, 
 matrices obtained from spatial time-frequency distributions
 \citep{Pham2001c}, slices of the cumulant tensor \citep{JADE:NC} 
 or Hessians of the characteristic function \citep{Yeredor00}.
 Generally, for stationary and temporally correlated signals we can
 define a set
%
of matrices $C^k$ 
with entries 
\be
(C_{(x)})^k_{ij}=\frac{1}{2}\sum_{t=1}^{T}x_i(t)(\Phi^k\star x_j)(t) +
x_j(t)(\Phi^{k}\star x_i)(t),
\label{eq2}
\ee
where $\star$ denotes convolution and $\Phi^k(t)$, $k=1,\ldots,K$ are
linear filters \citep{Ziehe00b}.

We note that in the popular special case where the 
$\Phi^k$
are  simple time-shift  operators
 $\Phi^{\tau}(t)=\delta_{t\tau}$ 
\citep[cf.][]{Tong91,MolgedeySchuster1994,Belo97} 
the matrices  defined by Equation (\ref{eq2}) may become indefinite  for certain choices of $\tau$.
Furthermore, in practice, the above target matrices have always to be estimated from the available data
which inevitably  gives rise to estimation errors. Hence the best we can do is to
find the matrix which diagonalizes the estimated target set ``as good as possible''.
Since we are able to perform the approximate joint diagonalization with a
non-orthogonal transformation, we  avoid the problematic pre-whitening step and
 obtain an estimate of the mixing matrix $A=V^{-1}$ by applying
our {\FFD} algorithm directly to the empirical matrices (\ref{eq2}).
Algorithm \ref{algo:FFSEP} summarizes the typical steps in an application  to BSS.
\begin{algorithm}
\begin{algorithmic} 
\STATE
\STATE INPUT:  $x(t)$, $\Phi^k$
%\STATE
\STATE $ C^k=\ldots $    \COMMENT{Estimate a number of matrices $C^k$ according to
Equation (\ref{eq2})}
\STATE $V=\mbox{\FFD}(C^k)$ \COMMENT{Apply joint diagonalization method}
\STATE $u(t) =V x(t)$  \COMMENT{unmix signals}
%\STATE
\STATE OUTPUT: u(t),\, V  
\end{algorithmic}
\caption{The {\FFS} algorithm \label{algo:FFSEP}}
\end{algorithm}


\section{Numerical Simulations}
\label{sec:sim}

The experiments in this Section are intended to compare the {\FFD}
algorithm with state-of-the-art algorithms for simultaneous
diagonalization and to illustrate the performance of our algorithm in 
BSS applications. As we mentioned in the introduction, there exist
at least three alternative formulations of the simultaneous
diagonalization problem. The most successful
algorithms representing the respective approaches were chosen for comparison. 

We present the results of five progressively
more complex experiments. First, we perform a ``sanity check''
experiment on a relatively easy set of perfectly diagonalizable
matrices. This experiment is intended to emphasize that for 
small-size diagonalizable matrices the
algorithm's performance matches the expected quadratic convergence.
In the second experiment
 we compare the {\FFD} algorithm
with the extended Jacobi method as used in the JADE algorithm of \citet{Car93} (orthogonal Frobenius
norm formulation), Pham's algorithm for
positive-definite matrices \citep{Pham2001c} and Yeredor's AC-DC
algorithm \citep{Yer02} (non-orthogonal, subspace fitting
formulation). In the third experiment we investigate the scaling behavior
of our algorithm as compared to AC-DC. Furthermore, the
performance of the {\FFD} algorithm is tested and compared with the AC-DC
algorithm on noisy, non-diagonalizable matrices.
Finally, the application of our algorithm to BSS is illustrated.
%
\subsection{``Sanity Check'' Experiment}
\label{sec:sanity}
The test data in this experiment is generated as follows. We use $K=15$ diagonal
matrices $D^k$ of size $5 \times 5$ where the elements on the
diagonal are drawn from a uniform distribution in the range $[-1,
\ldots,1]$ \citep[cf.][]{Joho2002}. These matrices are `mixed' by an
orthogonal
 matrix $A$ 
according to $A\,D^kA^T$ to generate the set of target
matrices $\{C^k\}$ to be diagonalized.\footnote{ The orthogonal matrix  was obtained from  a singular value decomposition of a random  $5 \times 5$ matrix, where the entries are drawn from a standard normal distribution.}  
The {\FFD} algorithm is initialized with the identity matrix
$V_{(0)}=I$, and the skew-symmetric update rule (\ref{W_ortho}) is used. 

The convergence behavior of the algorithm in 10 runs is shown in
Figure \ref{Fig1}. The diagonalization error is measured by the $\off{\cdot}$
function. One can see that the algorithm has converged to the correct solution
after less than 10 iterations in all trials. The quadratic
convergence is observed from early iterations.
\begin{figure}[H]
\centerline{\includegraphics[width=0.54\linewidth]{fig1}}
\caption{Diagonalization errors of the {\FFD} algorithm for a diagonalizable problem.}
       \label{Fig1}
\end{figure}
%

 
\subsection{Comparison with the State-of-the-Art Algorithms}
\label{sec:comp}
Two scenarios are considered for a comparison of the four selected algorithms: {\FFD}, the extended Jacobi method, Pham's algorithm and AC-DC. First, we test these 
algorithms on diagonalizable matrices under the conditions
satisfying the assumptions of all of them. Such conditions are:
positive-definiteness of the target matrices $C^k$ and orthogonality
of the true transformation $A$ used to generate those matrices. 
These conditions are met by generating the target matrices $C^k = A\, D^k A^T$
where $D^k$ are diagonal matrices with positive entries on the main
diagonal. The data set consists of 100 random matrices of size $10
\times 10$ satisfying the conditions above.

A comparison of the four algorithms on orthogonal positive-definite
matrices is shown in Figure \ref{Fig2}. Two runs of the algorithms are
presented, for the AC-DC algorithm 5 AC steps were interlaced with 1
DC step at each iteration. Although the algorithms optimize different
objective functions, the $\off{\cdot}$ function is still an adequate
evaluation criterion provided that the arbitrary scale is properly
normalized.

To achieve this, we evaluate $\sum_k\off{\hat{A}^{-1} C^k\hat{A}^{-T}
} $ where $\hat{A}$ is the normalized estimated mixing matrix.  At the
true solution the criterion must attain zero.  One can see that the
convergence of Pham's algorithm, the extended Jacobi method and {\FFD}
is quadratic, whereas the AC-DC algorithm converges linearly.
\begin{figure}[ht]
 \centerline{\includegraphics[width=0.5\linewidth]{fig2.eps}}
 \caption{Comparison of the  {\FFD}, the extended Jacobi method, Pham's algorithm and AC-DC
   in the orthogonal, positive-definite case: diagonalization error per iteration  measured by the $\off{\cdot}$ criterion.}        \label{Fig2}  
\end{figure}
The average iteration complexity of the four algorithms is shown in
Table~\ref{foo}. It follows from this table that the {\FFD} algorithm
indeed lives up to its name: its running time per iteration is
superior to both Pham's algorithm and AC-DC, and is comparable to the
extended Jacobi method algorithm.\footnote{In all experiments, MATLAB
  implementations of the algorithms were run on a standard PC with a
  750MHz clock.}
\begin{table}[ht]
\begin{center}   
\begin{tabular}{|c|c|c|c|}
\hline 
 {\FFD} & ext.~Jacobi &  Pham's & AC-DC \\ \hline\hline
   0.025  &0.030& 0.168 & 2.430 \\\hline
\end{tabular}
\end{center}
\caption{Comparison of the  {\FFD}, ext.~Jacobi, Pham's and AC-DC
   algorithms in the orthogonal, positive-definite case: average running time
   per iteration in seconds.}   
\label{foo}
\end{table}
%

In the second scenario, the comparison of the {\FFD} and the AC-DC
algorithms is repeated for non-positive-definite matrices obtained
from a non-orthogonal mixing matrix. This case cannot be handled by
the other two algorithms, therefore they are omitted from the
comparison. The convergence plots are shown in Figure \ref{Fig2b};
average running time per iteration is reported in Table~\ref{bar}.
Convergence behavior of the two algorithms is the same as in the
orthogonal, positive-definite case; the running time per iteration of
{\FFD} increases due to the use of non-skew-symmetric updates.
%
\begin{figure}[ht]
 \centerline{\includegraphics[width=0.54\linewidth]{fig2b.eps}}
 \caption{Comparison of the {\FFD} and AC-DC algorithms in the
   non-orthogonal, non-positive-definite case: diagonalization error per iteration measured by the $\off{\cdot}$ criterion.}  
\label{Fig2b}  
\end{figure}
\begin{table}[ht]
\begin{center}   
\begin{tabular}{|c|c|}
\hline 
 {\FFD} & AC-DC \\ \hline\hline
   0.034  &2.64 \\\hline
\end{tabular}
\caption{Comparison of the {\FFD} and AC-DC algorithms in the
   non-orthogonal, non-positive-definite case: average running time
  per iteration in seconds.}   
\label{bar}
\end{center}
\end{table}
%

\subsection{Scaling Behavior of {\FFD}}

Scalability is essential for application of an algorithm to real-life
problems. The most important parameter of the simultaneous
diagonalization problem affecting the scalability of an algorithm is
the size of the matrices. Figure \ref{Fig3} shows the running time per
iteration of the {\FFD} and the AC-DC algorithms for 
problems with increasing matrix sizes, plotted at logarithmic scale. 
One can see that both algorithm exhibit running times of $O(N^2)$; however,
in  absolute terms the {\FFD } algorithm is almost two orders of
magnitude faster.\footnote{This seemingly controversial result---theoretically expected scaling factor of AC-DC is $O(N^3)$---is
  due to high constants hidden in the setup phase of AC-DC. The setup
  phase has $O(N^2)$ complexity, but because of the constants it
  outweighs the main part of the algorithm in our experiment.}
\begin{figure}[h]
 \centerline{\includegraphics[width=0.54\linewidth]{fig3.eps}}
 \caption{Scaling of the {\FFD} and AC-DC algorithms with respect to the
 matrix size. Two repetitions of the experiment have been performed.}  
\label{Fig3}  
\end{figure}

\subsection{Non-diagonalizable Matrices}
We now investigate the impact of non-diagonalizability of the set of
matrices on the performance of the {\FFD} algorithm. Again, two
scenarios are considered: the one of the ``sanity check'' experiment
and the comparative analysis against the established algorithms. 
Non-diagonalizability is modeled by adding a random non-diagonal
symmetric ``noise'' matrix to each of the input matrices:
 $$C^k = A\,D^kA^T + \sigma^2 (R^k) (R^k)^T, $$ where the elements of $R^k$ are drawn from a
standard normal  distribution. The parameter $\sigma$ allows one to
control the impact of the non-diagonalizable component.
Another example, with a more realistic noise model, will be presented in 
subsection \ref{ssec:bss}. 

Figure \ref{Fig1a} shows the convergence plots of {\FFD} for various
values of $\sigma$.  The experimental setup is the same as in
Section \ref{sec:sanity}, apart from the additive noise. The impact of the
latter can be quantified by computing the $\off{\cdot}$ function on
the noise terms only (averaged over all runs), which is shown by the
dotted line in Figure \ref{Fig1a}.
One can see that the algorithm converges quadratically to the level
determined by the noise factor.
%
\begin{figure}[h]
       \centerline{\includegraphics[width=1\linewidth]{fig1a.eps}}
       \caption{Diagonalization errors of the {\FFD} algorithm on
         non-diagonalizable matrices.}
       \label{Fig1a}
\end{figure}

Similar to the second scenario in Section \ref{sec:comp}, the previously mentioned
algorithms are tested on the problem of approximate joint
diagonalization with non-orthogonal transforms. 
(Only the extended Jacobi algorithm had to be
excluded from the comparison since it is not designed to work with
non-orthogonal diagonalizers.)
 However, in contrast to Section \ref{sec:comp}, 
positive-definite target matrices were generated in order to enable a comparison with
Pham's method. 

Furthermore, we introduce another measure to assess the algorithms' performance 
in the non-orthogonal, non-diagonalizable case.
In synthetic experiments with artifical data the distance from the
true solution is a good evaluation criterion.
To be meaningful, this
distance has to be invariant w.r.t.~the irrelevant scaling and permutation
ambiguities. For this reason, we choose a performance index that is
 commonly used in the context of ICA/BSS where the same invariances
exist \citep[see e.g.~in][]{Amari98,JADE:NC}. Following the formulation of \citet{Moreau01}
a suitable performance index is defined on the normalized ``global'' matrix $G \ddef VA$
according to 

\begin{equation}\label{eq:score}
\mbox{score}(G)= \frac{1}{2} \left[ \sum_i \left( \sum_j \frac{|G_{ij}|^2}{\max\limits_l |G_{il}|^2}-1 \right) 
+\sum_j \left( \sum_i \frac{|G_{ij}|^2}{\max\limits_l |G_{lj}|^2}-1 \right) 
  \right].
\end{equation}
Clearly, this non-negative index attains zero iff $G$ is a product
of an invertible diagonal matrix $D$ and of a permutation matrix $P$,
i.e.,~$G=DP$.

The results of the comparison of the {\FFD}, Pham's and AC-DC
algorithms on a non-orthogonal positive-definite problem (5 matrices of
dimension $5 \times 5$) at various noise
levels are shown in Figure \ref{Fig5} for three typical runs.
The graphs illustrate some interesting ascpects of the convergence behavior of the algorithms.
Both the {\FFD} and Pham's algorithm
converge within a small number of iterations to approximately the same
error level. The AC-DC algorithm converges linearly, and
occasionally convergence can be very slow, as can be seen in each of
the plots in Figure \ref{Fig5}. However, when AC-DC converges, it
exhibits better performance as measured by the score
function; the higher the noise level, the stronger the difference.

\begin{figure}[h]
       \centerline{\includegraphics[width=0.99\linewidth]{exp_noiselevel_pd_many}}
       \caption{Comparison of the {\FFD}, Pham's and AC-DC algorithms
         in the non-diagonalizable, non-orthogonal, positive-definite
         case at various noise levels: performance index as
         measured by the score function (\ref{eq:score}).}
       \label{Fig5}
\end{figure}

\subsection{Blind Source Separation}
\label{ssec:bss}
Finally, we apply our method to a blind source separation task. The
signal matrix $S$ contains seven
audio signals containing 10000 points recorded at 8kHz and one Gaussian
noise source of the same length. These signals are mixed by a $8
\times 8$ Hadamard matrix,
$$A= \tiny \left(
%\begin{smallmatrix}
\begin{array}{rrrrrrrr}
      1 &   1 &   1 &   1  &  1 &   1 &   1 &   1 \\
      1 &  -1 &   1 &  -1  &  1 &  -1 &   1 &  -1 \\
      1 &   1 &  -1 &  -1  &  1 &   1 &  -1 &  -1 \\
      1 &  -1 &  -1 &   1  &  1 &  -1 &  -1 &   1 \\
      1 &   1 &   1 &   1  & -1 &  -1 &  -1 &  -1 \\
      1 &  -1 &   1 &  -1  & -1  &  1 &  -1 &   1 \\
      1 &   1 &  -1 &  -1  & -1  & -1 &   1 &   1 \\
      1 & -1 & -1 & 1 & -1 & 1 & 1 & -1
\end{array} 
%\end{smallmatrix}
\right).$$
This scaled orthogonal matrix produces a complete mixture, in the sense that each
observation contains a  maximal contribution from each source.
%
We compute 50 symmetrized, time-lagged correlation matrices according
to Equation (\ref{eq2}) with $\Phi^{\tau}(t)=\delta_{t\tau}$ and apply the
{\FFD} algorithm with $V_{(0)}=I$.
Figure \ref{Fig2o} shows the evolution of performance measure defined
in (\ref{eq:score}) and of the entries of the (normalized) global
system  $V_{(n)} A $. One can see that the difference from the true
solution, in terms of the score function, approaches zero and that $V_{(n)} A$
converges to a permutation matrix (as shown in the middle and the right panels).  
 \begin{figure}[hbt]
    \centerline{\includegraphics[width=0.9\linewidth]{convergence_horz_exp2.eps}}
       \caption{Convergence progress of the {\FFS} algorithm on the
         BSS task. The middle panel shows the evolution of the entries of the normalized global matrix $G$. The right panel shows those entries for the final (8th)  iteration  in matrix form and indicates that the cross-talk is minimized since the matrix G resembles a scaled and permutated identity matrix. Here, black, white and gray squares correspond to values -1, 1 and 0, respectively. }
       \label{Fig2o}
     \end{figure}

In order to study the behavior of the \FFD~algorithm in a more realistic noisy
scenario the following experiment is conducted.
The data is generated by mixing three stationary, time-correlated sources
with the fixed matrix  $A= \left(
\begin{smallmatrix}
  8  &   1  &   6\\
     3 &    5 &    7\\
     4 &    9 &    2\\
\end{smallmatrix}  \right).$
%
The sources are generated by feeding an i.i.d.~random noise signal
into a randomly chosen, auto-regressive (AR) model of order 5 whose 
coefficients are drawn from a standard normal distribution
and are sorted in decreasing order (to ensure stability).
The generated signals have a total length of $50000$ samples.
To separate the sources we estimate 10 symmetrized, time-lagged
correlation matrices of the mixed signals according to Equation (\ref{eq2})
with $\Phi^{\tau}(t)=\delta_{t\tau}$ and perform simultaneous diagonalization
of these matrices.

Clearly, the quality of the estimation depends on the number $T$ of
samples used to estimate these correlation matrices.  By varying $T$
we can simulate different noise levels corresponding to the variance
of the estimates, which is more realistic than corrupting the target
matrices with small i.i.d.~additive noise.

The results of the experiment are shown in Figure \ref{Fig6}.
Performance of the {\FFD} and the AC-DC algorithms, as measured by the
score (\ref{eq:score}), is displayed for four different sample sizes,
the smaller samples corresponding to the higher noise level. 100
repetitions are performed for each sample size, and the 25\%, 50\% and
75\% quantiles of the log-score are shown in the plots. Two
observations can be made from Figure \ref{Fig6}: 
{\FFD} converges much faster than AC-DC, and when converged,
{\FFD} yields a better score (on average), with the difference  more pronounced
for samples sizes 10000 and 30000 in our experiment.

\begin{figure}[h]
\centerline{
\includegraphics[width=0.25\linewidth]{tube_plots_1}
\includegraphics[width=0.25\linewidth]{tube_plots_2}
\includegraphics[width=0.25\linewidth]{tube_plots_3}
\includegraphics[width=0.25\linewidth]{tube_plots_4}
}
       \caption{ Performance of {\FFD} and AC-DC measured by the log of the score (\ref{eq:score})
         for different sample sizes and 100 trials each. 25\% (lower edge of the shaded region), 50\% (thick line in the middle) and 75\% quantiles (upper edge of the shaded region) are shown. }
       \label{Fig6}
\end{figure}
%\newpage

\section{Discussion and Conclusions}
\label{sec:conclu}

We have presented a new algorithm {\FFD} for simultaneous
diagonalization of a set of matrices. The algorithm is based on the
Frobenius norm formulation of the simultaneous diagonalization problem
and provides an efficient means of diagonalization in the absence of
additional constraints, such as orthogonality or
positive-definiteness. The important feature of our algorithm is the
direct enforcement of invertibility of the diagonalizer; in previous
work this was usually achieved by an orthogonality constraint which
reduces the space of solutions.

The efficiency of the {\FFD} algorithm lies in the special second-order 
approximation of the cost function, which yields a block-diagonal Hessian
and thus allows for highly efficient computation of the Newton update step.
%
Although, theoretically, such approximation can be seen as a weakness of
the approach---and raise the question of whether the point of the
algorithm's convergence is indeed an optimizer of the full cost
function---we have empirically observed that the solution found by
the algorithm is of good quality for practical applications.

A series of comparisons of the {\FFD} algorithm with 
state-of-the-art diagonalization algorithms is presented under a
number of conditions that can or cannot be handled by other
algorithms. The main conclusions of this comparative evaluation is
that our algorithm is competitive with the best algorithms
(i.e.~Jacobi-based and Pham's algorithm) that impose additional
constraints either on the class of solutions or the type of input
data. {\FFD} is significantly more efficient---as far as both the
scaling factors and the absolute constants are concerned---than the
AC-DC algorithm, the only general algorithm applicable under the same
conditions as ours. The {\FFD} algorithm can be applied to matrices of
dimensions in the hundreds of rows/columns, under no additional
assumptions. It also performs reliably on non-diagonalizable data, for
which only an approximate solution is possible.


Several interesting research topics can be anticipated.  From a
theoretical point of view, convergence analysis could yield further
insights into the numerical behavior of {\FFD} as well as a better
understanding of the general techniques for optimization over
nonholonomic manifolds that the algorithm belongs to.
Further investigation of the robustness of joint diagonalization
algorithms in the presence of various forms of noise is a very
interesting practical issue.  Numerous applications of the algorithm
to real-life problems can be clearly foreseen. 

\acks{
The authors thank Benjamin Blankertz, Steven Lemm, Christin Sch\"afer,  Sebastian Mika,
Stefan {Har\-me\-ling}, Frank Meinecke, Guido Dornhege, Motoaki
Kawanabe, David Tax, Julian Laub, Matthias Krauledat, Marcel Joho, Michael
Zibulevsky and Arie Yeredor for sharing their insights and expertise in many 
fruitful discussions.
Furthermore, the in-depth comments and valuable
suggestions of the three anonymous reviewers are highly appreciated. 
This helped us to improve the initial version of the manuscript.

A.Z., P.L. and K.-R.M. acknowledge partial funding in the EU project
(IST-1999-14190 -- BLISS), by BMBF under contract FKZ~01IBB02A, the
SFB 618 and the PASCAL Network of
Excellence  (EU~\#506778). G.N. has been supported by a grant
from the National Foundation for Functional Brain Imaging.
}

\bibliography{lmdiag,BSS_all}


\end{document}

