# Joint Approximate Diagonalization under Orthogonality Constraints

Ronald de Vlaming111School of Business and Economics, Vrije Universiteit
Amsterdam, Amsterdam, The Netherlands<EMAIL_ADDRESS>https://research.vu.nl/en/persons/ronald-de-vlaming). and Eric A.W. Slob222MRC
Biostatistics Unit, School of Clinical Medicine, University of Cambridge,
Cambridge, United Kingdom (mailto:eric.slob@mrc-bsu.cam.ac.uk).

## Abstract

Joint diagonalization of a set of positive (semi)-definite matrices has a wide
range of analytical applications, such as estimation of common principal
components, estimation of multiple variance components, and blind signal
separation. However, when the eigenvectors of involved matrices are not the
same, joint diagonalization is a computationally challenging problem. To the
best of our knowledge, currently existing methods require at least $O(KN^{3})$
time per iteration, when $K$ different $N\times N$ matrices are considered. We
reformulate this optimization problem by applying orthogonality constraints
and dimensionality reduction techniques. In doing so, we reduce the
computational complexity for joint diagonalization to $O(N^{3})$ per quasi-
Newton iteration. This approach we refer to as JADOC: Joint Approximate
Diagonalization under Orthogonality Constraints. We compare our algorithm to
two important existing methods and show JADOC has superior runtime while
yielding a highly similar degree of diagonalization. The JADOC algorithm is
implemented as open-source Python code, available at
https://github.com/devlaming/jadoc.

Key words. diagonalization, blind signal separation, variance component
estimation

AMS subject classifications. 65F25, 90C53

##  1 Introduction

Given a set of $N\times N$ symmetric matrices, $\textbf{C}_{k}$ for
$k=1,\ldots,K$, that are (at least) positive semidefinite (PSD) the objective
of finding an $N\times N$ matrix B such that
$\textbf{B}\textbf{C}_{k}{\textbf{B}}^{\top}$ is ‘as diagonal as possible’ for
$k=1,\ldots,K$ arises naturally in estimation of common principal components
[Pham, 2001].

Assuming matrices $n_{k}\textbf{C}_{k}$, for $k=1,\ldots,K$, follow Wishart
distributions with $n_{k}$ degrees of freedom and covariance matrix
$\boldsymbol{\Gamma}_{k}$ respectively, such that (_i_)
$\boldsymbol{\Gamma}_{k}$ have the same eigenvectors for $k=1,\ldots,K$ and
(_ii_) $n_{k}\textbf{C}_{k}$ and $n_{l}\textbf{C}_{l}$ are independent for
$k\neq l$, one can apply maximum likelihood estimation (MLE) to obtain
estimates of these common eigenvectors. The resulting optimization problem can
be rewritten as finding a matrix B that minimizes the following criterion
[Pham, 2001]:

$\mathcal{L}\left(\textbf{B}\right)=\frac{1}{2}\sum_{k=1}^{K}n_{k}\left(\log\det\operatorname{diag}\left(\textbf{B}\textbf{C}_{k}{\textbf{B}}^{\top}\right)-\log\det\left(\textbf{B}\textbf{C}_{k}{\textbf{B}}^{\top}\right)\right).$
(1)

This criterion provides a reasonable overall measure of deviation of
$\textbf{B}\textbf{C}_{k}{\textbf{B}}^{\top}$ from diagonality for
$k=1,\ldots,K$. Thus, this MLE-based criterion, as proposed by Pham [Pham,
2001], is a natural starting point for finding a matrix B that aims to
diagonalize $\textbf{B}\textbf{C}_{k}{\textbf{B}}^{\top}$ for $k=1,\ldots,K$.

Subsequent research on this general approach to find B has been fairly
limited, however, possibly owing to the fact that this optimization problem is
computationally challenging and has poor scalability [Mesloub et al., 2014].
Existing approach typically use so-called sweeps [Pham, 2001]. Within a sweep,
an optimization is performed for all pairwise combinations of rows of B.

Importantly, relatively recently, a quasi-Newton method has been proposed to
find a matrix B that minimizes Equation 1 more efficiently [Ablin et al.,
2018]. The great advantage of that particular approach is that all elements of
B are updated jointly in each iteration. However, that approach requires at
least $K$ matrix multiplications of $N\times N$ matrices in order to calculate
the gradient of the loss function with respect to B, thus requiring
$O(KN^{3})$ time per iteration.

Finally, that approach does not utilize further computational speed-ups that
are possible by enforcing a simple orthogonality constraint on B, viz., where
we require B to be orthonormal (i.e.,
$\textbf{B}{\textbf{B}}^{\top}=\textbf{I}_{N}$). Here we note that although
this so-called orthogonal joint diagonalization (OJD) property [Mesloub et
al., 2014] is built into several existing methods, most notably JADE [Cardoso
and Souloumiac, 1996], those methods still use sweeps, where each sweep also
requires $O(KN^{3})$ time.

Instead, we (_i_) follow a quasi-Newton approach to jointly update all
elements B in each iteration, (_ii_) exploit orthonormality to further
simplify the criterion, (_iii_) parametrize the optimization problem such that
we can use unconstrained optimization while guaranteeing orthonormality of B,
(_iv_) use a low-dimensional approximation of $\textbf{C}_{k}$ to avoid
calculation of matrix products of $N\times N$ matrices for $k=1,\ldots,K$ in
each iteration, (_v_) apply a mild regularization to the approximations of
matrices $\textbf{C}_{k}$ such that it does not perturb our solution for B,
and (_vi_) linearize updates to efficiently implement a golden-section search
[Kiefer, 1953] within each iteration.

This overall approach, incorporating these six points, we refer to as JADOC:
joint approximate diagonalization under orthogonality constraints. JADOC
requires $O(N^{3})$ time per iteration, constituting an $O(K)$ reduction in
runtime compared to existing methods. Moreover, the marginal cost of
performing the golden-section search within each iteration is negligible.

Our main derivations are provided in Section 2, our JADOC algorithm is
presented in Section 3 together with the algorithm that we use to simulate
data, simulation results are discussed in Section 4, and the conclusions
follow in Section 5.

##  2 Main derivations

We follow Pham’s notation [Pham, 2001], while using the criterion as proposed
by Ablin et al. [Ablin et al., 2018], which is defined as

$\mathcal{L}\left(\textbf{B}\right)=\frac{1}{2K}\sum_{k=1}^{K}\left(\log\det\operatorname{diag}\left(\textbf{B}\textbf{C}_{k}{\textbf{B}}^{\top}\right)-\log\det\left(\textbf{B}\textbf{C}_{k}{\textbf{B}}^{\top}\right)\right).$
(2)

This criterion is in fact equivalent to that in Equation 1, except for the
scaling by $n_{k}$, which effectively weights $\textbf{C}_{k}$ by its degrees
of freedom. We ignore this scaling, as we do not necessarily assume
$\textbf{C}_{k}$ to follow a Wishart distribution. We simply want to find B
that diagonalizes symmetric, PSD matrices $\textbf{C}_{k}$ for $k=1,\ldots,K$
as much as possible, using the criterion in Equation 2 (i.e., without
assigning particular weights to the individual matrices $\textbf{C}_{k}$).

Given that (_i_)
$\det(\textbf{A}_{1}\textbf{A}_{2})=\det(\textbf{A}_{1})\det(\textbf{A}_{2})$
for square matrices $\textbf{A}_{1}$ and $\textbf{A}_{2}$, and (_ii_)
$\textbf{B}{\textbf{B}}^{\top}=\textbf{I}_{N}$, we can rewrite this criterion
as follows

$\displaystyle\mathcal{L}\left(\textbf{B}\right)$
$\displaystyle=\frac{1}{2K}\sum_{k=1}^{K}\left(\log\det\operatorname{diag}\left(\textbf{B}\textbf{C}_{k}{\textbf{B}}^{\top}\right)-\log\left(\det\left(\textbf{B}\right)\det\left(\textbf{C}_{k}\right)\det\left({\textbf{B}}^{\top}\right)\right)\right)$
$\displaystyle=\frac{1}{2K}\sum_{k=1}^{K}\left(\log\det\operatorname{diag}\left(\textbf{B}\textbf{C}_{k}{\textbf{B}}^{\top}\right)-\log\left(\det\left(\textbf{C}_{k}\right)\det\left(\textbf{B}{\textbf{B}}^{\top}\right)\right)\right)$
$\displaystyle=\frac{1}{2K}\sum_{k=1}^{K}\left(\log\det\operatorname{diag}\left(\textbf{B}\textbf{C}_{k}{\textbf{B}}^{\top}\right)-\log\left(\det\left(\textbf{C}_{k}\right)\right)\right).$

Clearly, the second log-determinant in the summand is independent of B.
Consequently, this term merely contributes to the constant. Thus, we can
safely drop the second term from the loss function. Moreover, let the
eigenvalue decomposition of $\textbf{C}_{k}$ be given by

$\textbf{C}_{k}=\textbf{P}_{k}\textbf{D}_{k}^{2}{\textbf{P}}^{\top}_{k},$ (3)

where $\textbf{P}_{k}$ denotes the orthonormal matrix of eigenvectors and
$\textbf{D}_{k}^{2}$ the diagonal matrix of eigenvalues. Now, further,
consider a rank-$S$ approximation to $\textbf{C}_{k}$, defined as

$\widetilde{\textbf{C}}_{k}=\textbf{L}_{k}{\textbf{L}}^{\top}_{k},\textrm{
where }\textbf{L}_{k}=\widetilde{\textbf{P}}_{k}\widetilde{\textbf{D}}_{k},$
(4)

where $\widetilde{\textbf{D}}_{k}$ is an $S\times S$ diagonal matrix, with the
square root of the $S$ leading eigenvalues of $\textbf{C}_{k}$ as its diagonal
elements, and $\widetilde{\textbf{P}}_{k}$ denotes the $N\times S$ matrix of
corresponding eigenvectors.

JADOC aims to find a matrix B that diagonalizes this rank-$S$ approximation,
$\widetilde{\textbf{C}}_{k}$, for $k=1,\ldots,K$. Importantly, when
$\widetilde{\textbf{C}}_{k}$ is rank-deficient, for at least one
$k\in\\{1,\ldots,K\\}$, the loss function is ill-defined. However, observe
that a matrix B that diagonalizes $\widetilde{\textbf{C}}_{k}$ also
diagonalizes $\widetilde{\textbf{C}}_{k}+\lambda\textbf{I}$ for some fixed
value of $\lambda$.

Moreover, by setting $\lambda>0$, matrix
$\widetilde{\textbf{C}}_{k}+\lambda\textbf{I}$ is positive definite for
$k=1,\ldots,K$. Diagonalization of these matrices is, thus, a valid approach
to find a matrix B that also (approximately) diagonalizes
$\widetilde{\textbf{C}}_{k}$, with the additional advantage that the loss
function is properly defined for all orthonormal real matrices B. Thus, more
concretely, the aim of JADOC is to find an orthonormal, real matrix B that
diagonalizes the rank-$S$ approximations of $\textbf{C}_{k}$ for
$k=1,\ldots,K$, by minimizing:

$\displaystyle\mathcal{L}\left(\textbf{B}\right)$
$\displaystyle=\frac{1}{2K}\sum_{k=1}^{K}\log\det\operatorname{diag}\left(\textbf{B}\left(\widetilde{\textbf{C}}_{k}+\lambda\textbf{I}\right){\textbf{B}}^{\top}\right),\textrm{
where}$ (5) $\displaystyle\lambda$
$\displaystyle=\lambda_{0}+\frac{1}{NK}\sum_{k=1}^{K}\operatorname{tr}\left(\textbf{C}_{k}-\widetilde{\textbf{C}}_{k}\right),\textrm{
where }\lambda_{0}=1.$ (6)

The first term in our definition of $\lambda$ (i.e., $\lambda_{0}$) ensures
the loss function is well-defined even though $\widetilde{\textbf{C}}_{k}$ is
rank deficient for sure when $S<N$ and irrespective of whether the original
input matrices $\textbf{C}_{k}$ are rank deficient or not. The second term
ensures the scale of the diagonal elements in
$\widetilde{\textbf{C}}_{k}+\lambda\textbf{I}$ is on average roughly the same
as the diagonal elements of $\textbf{C}_{k}$. Further, we note that the trace
in Equation 6 simply equals the trace of $\textbf{C}_{k}$ minus the sum of its
$S$ leading eigenvalues.

We can now simplify the preceding criterion to the main JADOC criterion:

$\mathcal{L}\left(\textbf{B}\right)=\frac{1}{2K}\sum_{k=1}^{K}\sum_{i=1}^{N}\log\left(\lambda+\sum_{j=1}^{S}\left(\textbf{A}_{k}\right)_{ij}^{2}\right),\textrm{
where }\textbf{A}_{k}=\textbf{B}\textbf{L}_{k},$ (7)

where $\lambda$ is as defined in Equation 6. Now, by setting

$S=\left\lceil\frac{N}{K}\right\rceil,$ (8)

and given (_i_) we pre-calculate $\textbf{L}_{k}$, requiring $O(KN^{3})$ time
in total, and (_ii_) we want to evaluate the loss function for some
orthonormal B, we can calculate $\textbf{A}_{k}$ for $k=1,\ldots,K$ jointly in
$O(N^{3})$ time. This time complexity is the bottleneck for calculations in
each iteration.

We consider the choice of $S$ in Equation 8 to be fair: (_i_) it increases
with $N$, (_ii_) it extracts the most salient dimensions from each matrix
$\textbf{C}_{k}$, for $k=1,\ldots,K$, and (_iii_) it is the lowest $S$ such
that there may exist convex combinations of $\widetilde{\textbf{C}}_{k}$, for
$k=1,\ldots,K$, that are positive definite. Notice that users of the JADOC
tool may specify $S$ different from the recommended value seen in Equation 8.

###  2.1 Parametrization

Given our loss function in Equation 7 and the fact that we will use an
iterative procedure to obtain B, let us assume that in iteration $t$ we have B
from the directly preceding iteration, denoted by $\textbf{B}_{(t-1)}$. Here,
we further assume that $\textbf{B}_{(t-1)}$ is orthonormal. That is,
$\textbf{B}_{(t-1)}{\textbf{B}}^{\top}_{(t-1)}=\textbf{I}_{N}$.

Bearing these considerations in mind, the goal of iteration $t$ is to ‘update’
$\textbf{B}_{(t-1)}$ such that (_i_) the loss function is further decreased
and (_ii_) orthonormality of B is preserved. The easiest way to ensure
orthogonality is indeed maintained, is by using a rotation matrix R (i.e., a
matrix such that $\textbf{R}{\textbf{R}}^{\top}=\textbf{I}_{N}$) to set
$\textbf{B}_{(t)}=\textbf{R}\textbf{B}_{(t-1)}$, as clearly under this
definition of $\textbf{B}_{(t)}$, we have that
$\textbf{B}_{(t)}{\textbf{B}}^{\top}_{(t)}=\textbf{I}_{N}$.

However, defining such an $N\times N$ matrix R ‘straight up’ is far from
trivial. Fortunately, for any rotation matrix R there exists a skew-symmetric
matrix S such that $e^{\textbf{S}}=\textbf{R}$ [Gallier and Xu, 2002, Bröcker
and Tom Dieck, 1985], where $e^{c\textbf{M}}$ denotes the matrix exponential
of square matrix M, which is defined in terms of the following power series
[Marsden and Ratiu, 1999]:

$e^{c\textbf{M}}=\textbf{I}+\sum_{k=1}^{\infty}\frac{c^{k}}{k!}\textbf{M}^{k}.$
(9)

Thus, we can achieve our second aim in any iteration by defining a strictly
lower-triangular matrix E (i.e., with diagonal elements equal to zero), and
considering each of its free elements
$\varepsilon_{lm}=\left(\textbf{E}\right)_{lm}\in\mathbb{R}$ for
$l,m=1,\ldots,N$ where $l>m$ as a parameter of interest. Under this
parametrization, (_i_) $\textbf{S}=\textbf{E}-{\textbf{E}}^{\top}$ defines a
skew-symmetric matrix (i.e., such that ${\textbf{S}}^{\top}=-\textbf{S}$), for
which it holds that $e^{\textbf{S}}$ is a rotation matrix R (i.e., such that
$\textbf{R}{\textbf{R}}^{\top}=\textbf{I}_{N}$), and (_ii_) any B that meets
the orthonormality requirement can be written as follows

$\textbf{B}_{(t)}=e^{\textbf{E}-{\textbf{E}}^{\top}}\textbf{B}_{(t-1)}.$ (10)

Under these definitions, letting $\textbf{A}_{(t)k}$ denote $\textbf{A}_{k}$
in iteration $(t)$, notice that

$\textbf{A}_{(t)k}=\textbf{B}_{(t)}\textbf{L}_{k}=e^{\textbf{E}-{\textbf{E}}^{\top}}\textbf{B}_{(t-1)}\textbf{L}_{k}=e^{\textbf{E}-{\textbf{E}}^{\top}}\textbf{A}_{(t-1)k},$
(11)

implying we can define both B and $\textbf{A}_{k}$ in a recursive manner
across iterations.

Bearing all these considerations in mind, we can write our loss function in
terms of the free parameters in E and $\textbf{A}_{k}$ from the preceding
iteration as follows:

$\displaystyle\mathcal{L}\left(\textbf{E}\right)$
$\displaystyle=\frac{1}{2K}\sum_{k=1}^{K}\sum_{i=1}^{N}\log\left(d_{ik}\left(\textbf{E}\right)\right),\textrm{
where }$ (12) $\displaystyle d_{ik}\left(\textbf{E}\right)$
$\displaystyle=\lambda+\sum_{j=1}^{S}\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\textbf{A}_{(t-1)k}\right)_{ij}^{2}.$
(13)

###  2.2 First derivatives

Here, let us consider the partial derivative of this loss function with
respect to $\varepsilon_{lm}$ for $l>m$. This derivative can be written as

$\frac{\partial\mathcal{L}\left(\textbf{E}\right)}{\partial\varepsilon_{lm}}=\frac{1}{K}\sum_{k=1}^{K}\sum_{i=1}^{N}\frac{1}{d_{ik}\left(\textbf{E}\right)}\sum_{j=1}^{N}\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\textbf{A}_{(t-1)k}\right)_{ij}\frac{\partial\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\textbf{A}_{(t-1)k}\right)_{ij}}{\partial\varepsilon_{lm}},$
(14)

where

$\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\textbf{A}_{(t-1)k}\right)_{ij}=\sum_{p=1}^{N}\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\right)_{ip}\left(\textbf{A}_{(t-1)k}\right)_{pj}.$

Furthermore, using the power series in Equation 9, it is easy to show that

$\frac{\partial
e^{\textbf{E}-{\textbf{E}}^{\top}}}{\partial\varepsilon_{lm}}=e^{\textbf{E}-{\textbf{E}}^{\top}}\left(\textbf{J}_{lm}-\textbf{J}_{ml}\right),$

where $\textbf{J}_{lm}$ is the single-entry matrix, which equals one for
element $l,m$ and zero elsewhere. Thus, we have that

$\displaystyle\frac{\partial\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\textbf{A}_{(t-1)k}\right)_{ij}}{\partial\varepsilon_{lm}}$
$\displaystyle=\sum_{p=1}^{N}\frac{\partial\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\right)_{ip}}{\partial\varepsilon_{lm}}\left(\textbf{A}_{(t-1)k}\right)_{pj}$
$\displaystyle=\sum_{p=1}^{N}\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\left(\textbf{J}_{lm}-\textbf{J}_{ml}\right)\right)_{ip}\left(\textbf{A}_{(t-1)k}\right)_{pj}$

Substituting this expression in Equation 14 yields the following expression
for the partial derivative of the loss function with respect to
$\varepsilon_{lm}$:

$\displaystyle\frac{\partial\mathcal{L}\left(\textbf{E}\right)}{\partial\varepsilon_{lm}}$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\sum_{i=1}^{N}\frac{1}{d_{ik}\left(\textbf{E}\right)}\sum_{j=1}^{N}\left[\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\textbf{A}_{(t-1)k}\right)_{ij}\right.$
$\displaystyle\hphantom{=\frac{1}{K}}\left.\sum_{p=1}^{N}\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\left(\textbf{J}_{lm}-\textbf{J}_{ml}\right)\right)_{ip}\left(\textbf{A}_{(t-1)k}\right)_{pj}\right]$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\sum_{i=1}^{N}\frac{1}{d_{ik}\left(\textbf{E}\right)}\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\textbf{A}_{(t-1)k}{\textbf{A}}^{\top}_{(t-1)k}\left(\textbf{J}_{ml}-\textbf{J}_{lm}\right)e^{{\textbf{E}}^{\top}-\textbf{E}}\right)_{ii}$

Evaluation of the loss function and gradient at current value of B, which
corresponds to $\textbf{E}=\textbf{0}$ and
$e^{\textbf{E}-{\textbf{E}}^{\top}}=\textbf{I}_{N}$, yields

$\displaystyle{\left.\frac{\partial\mathcal{L}\left(\textbf{E}\right)}{\partial\varepsilon_{lm}}\right|}_{\textbf{E}=\textbf{0}}$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\sum_{i=1}^{N}\frac{1}{d_{ik}}\left(\textbf{A}_{(t-1)k}{\textbf{A}}^{\top}_{(t-1)k}\left(\textbf{J}_{ml}-\textbf{J}_{lm}\right)\right)_{ii},\textrm{
where }$ (15) $\displaystyle d_{ik}$
$\displaystyle=d_{ik}\left(\textbf{0}\right)\equiv\lambda+\sum_{j=1}^{S}\left(\textbf{A}_{(t-1)k}\right)^{2}_{ij}$
(16)

Using properties of single-entry matrices, this expression can be simplified
as follows:

${\left.\frac{\partial\mathcal{L}\left(\textbf{E}\right)}{\partial\varepsilon_{lm}}\right|}_{\textbf{E}=\textbf{0}}=\frac{1}{K}\sum_{k=1}^{K}\left[\left(\frac{1}{d_{lk}}-\frac{1}{d_{mk}}\right)\left(\textbf{A}_{(t-1)k}{\textbf{A}}^{\top}_{(t-1)k}\right)_{lm}\right].$
(17)

Now, the gradient can be expressed in matrix notation as follows:

G
$\displaystyle=\frac{\partial\mathcal{L}\left(\textbf{E}\right)}{\partial\textbf{E}}=\operatorname{stril}\left(\textbf{F}-{\textbf{F}}^{\top}\right),\textrm{
where}$ (18) F
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\operatorname{diag}\left(\frac{1}{d_{1k}},\ldots,\frac{1}{d_{Nk}}\right)\left(\textbf{A}_{(t-1)k}{\textbf{A}}^{\top}_{(t-1)k}\right),$
(19)

where $\operatorname{stril}(\textbf{F}-{\textbf{F}}^{\top})$ is the strict
lower-triangular submatrix of square matrix $\textbf{F}-{\textbf{F}}^{\top}$.

###  2.3 Second derivatives

Let us now consider second partial derivatives, which are, in general, defined
by the following expression for $l>m$ and $p>q$:

$\displaystyle\frac{\partial^{2}\mathcal{L}\left(\textbf{E}\right)}{\partial\varepsilon_{lm}\partial\varepsilon_{pq}}$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\sum_{i=1}^{N}\frac{\partial\left(\frac{1}{d_{ik}\left(\textbf{E}\right)}\left(\textbf{T}_{k}\left(\textbf{E}\right)\right)_{ii}\right)}{\partial\varepsilon_{pq}},\textrm{
where}$ (20) $\displaystyle\textbf{T}_{k}\left(\textbf{E}\right)$
$\displaystyle=e^{\textbf{E}-{\textbf{E}}^{\top}}\textbf{A}_{(t-1)k}{\textbf{A}}^{\top}_{(t-1)k}\left(\textbf{J}_{ml}-\textbf{J}_{lm}\right)e^{{\textbf{E}}^{\top}-\textbf{E}}.$
(21)

Using the chain rule, the partial derivative on the right-hand side of
Equation 20 can be written as follows:

$\frac{\partial\left(\frac{1}{d_{ik}\left(\textbf{E}\right)}\left(\textbf{T}_{k}\left(\textbf{E}\right)\right)_{ii}\right)}{\partial\varepsilon_{pq}}=\frac{1}{d_{ik}\left(\textbf{E}\right)}\frac{\partial\left(\textbf{T}_{k}\left(\textbf{E}\right)\right)_{ii}}{\partial\varepsilon_{pq}}-\frac{1}{\left(d_{ik}\left(\textbf{E}\right)\right)^{2}}\frac{\partial
d_{ik}\left(\textbf{E}\right)}{\partial\varepsilon_{pq}}\left(\textbf{T}_{k}\left(\textbf{E}\right)\right)_{ii}.$
(22)

Immediately evaluating
$\left(\textbf{T}_{k}\left(\textbf{E}\right)\right)_{ii}$ in the second term
at $\textbf{E}=\textbf{0}$, yields

$\left(\textbf{T}_{k}\left(\textbf{0}\right)\right)_{ii}=\left\\{\begin{array}[]{ll}\left(\textbf{A}_{(t-1)k}{\textbf{A}}^{\top}_{(t-1)k}\right)_{lm}&\textrm{
if }i=l,\\\
-\left(\textbf{A}_{(t-1)k}{\textbf{A}}^{\top}_{(t-1)k}\right)_{lm}&\textrm{ if
}i=m,\\\ 0&\textrm{ otherwise. }\end{array}\right.$ (23)

Effectively, this implies the second term on the right-hand side of Equation
22 is zero when $\textbf{A}_{(t-1)k}{\textbf{A}}^{\top}_{(t-1)k}$ is diagonal.
Thus, assuming (near) diagonalization is possible and has already been
(almost) achieved in iteration $t-1$, this last term can be ignored.

Consequently, the most important part of the second partial derivatives is
finding an efficient expression for

$\frac{\partial\left(\textbf{T}_{k}\left(\textbf{E}\right)\right)_{ii}}{\partial\varepsilon_{pq}}=\frac{\partial\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\textbf{A}_{(t-1)k}{\textbf{A}}^{\top}_{(t-1)k}\left(\textbf{J}_{ml}-\textbf{J}_{lm}\right)e^{{\textbf{E}}^{\top}-\textbf{E}}\right)_{ii}}{\partial\varepsilon_{pq}}.$
(24)

We here, thus, assume that $\textbf{A}_{(t-1)k}{\textbf{A}}^{\top}_{(t-1)k}$
is diagonal. That is, we can replace
$\textbf{A}_{(t-1)k}{\textbf{A}}^{\top}_{(t-1)k}$ by diagonal matrix
$\textbf{D}_{k}$ of which element $i,i$ is given by $d_{ik}-\lambda$.
Therefore, we have that

$\frac{\partial\left(\textbf{T}_{k}\left(\textbf{E}\right)\right)_{ii}}{\partial\varepsilon_{pq}}=\sum_{s=1}^{N}\sum_{t=1}^{N}\frac{\partial\left(\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\right)_{it}\left(\textbf{D}_{k}\left(\textbf{J}_{ml}-\textbf{J}_{lm}\right)\right)_{ts}\left(e^{{\textbf{E}}^{\top}-\textbf{E}}\right)_{si}\right)}{\partial\varepsilon_{pq}}.$
(25)

Using the product rule and previously discussed properties of the partial
derivative of the matrix exponential with respect to a skew symmetric matrix,
it then follows that

$\displaystyle\frac{\partial\left(\textbf{T}_{k}\left(\textbf{E}\right)\right)_{ii}}{\partial\varepsilon_{pq}}$
$\displaystyle=\sum_{s=1}^{N}\sum_{t=1}^{N}\left[\frac{\partial\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\right)_{it}}{\partial\varepsilon_{pq}}\left(\textbf{D}_{k}\left(\textbf{J}_{ml}-\textbf{J}_{lm}\right)\right)_{ts}\left(e^{{\textbf{E}}^{\top}-\textbf{E}}\right)_{si}\right.$
$\displaystyle\hphantom{=}\left.+\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\right)_{it}\left(\textbf{D}_{k}\left(\textbf{J}_{ml}-\textbf{J}_{lm}\right)\right)_{ts}\frac{\partial\left(e^{{\textbf{E}}^{\top}-\textbf{E}}\right)_{si}}{\partial\varepsilon_{pq}}\right]$
$\displaystyle=\sum_{s=1}^{N}\sum_{t=1}^{N}\left[\left(\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\right)\left(\textbf{J}_{pq}-\textbf{J}_{qp}\right)\right)_{it}\left(\textbf{D}_{k}\left(\textbf{J}_{ml}-\textbf{J}_{lm}\right)\right)_{ts}\left(e^{{\textbf{E}}^{\top}-\textbf{E}}\right)_{si}\right.$
$\displaystyle\hphantom{=}\left.+\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\right)_{it}\left(\textbf{D}_{k}\left(\textbf{J}_{ml}-\textbf{J}_{lm}\right)\right)_{ts}\left(\left(\textbf{J}_{qp}-\textbf{J}_{pq}\right)\left(e^{{\textbf{E}}^{\top}-\textbf{E}}\right)\right)_{si}\right]$
$\displaystyle=\left(\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\right)\left(\textbf{J}_{pq}-\textbf{J}_{qp}\right)\textbf{D}_{k}\left(\textbf{J}_{ml}-\textbf{J}_{lm}\right)\left(e^{{\textbf{E}}^{\top}-\textbf{E}}\right)\right)_{ii}$
$\displaystyle\hphantom{=}+\left(\left(e^{\textbf{E}-{\textbf{E}}^{\top}}\right)\textbf{D}_{k}\left(\textbf{J}_{ml}-\textbf{J}_{lm}\right)\left(\textbf{J}_{qp}-\textbf{J}_{pq}\right)\left(e^{{\textbf{E}}^{\top}-\textbf{E}}\right)\right)_{ii}.$

Evaluating this expression at $\textbf{E}=\textbf{0}$ and substituting in our
expression for the second-order derivative under complete diagonality, we have
that

$\displaystyle\left.\frac{\partial^{2}\mathcal{L}\left(\textbf{E}\right)}{\partial\varepsilon_{lm}\partial\varepsilon_{pq}}\right|_{\textbf{E}=\textbf{0}}$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\sum_{i=1}^{N}\frac{1}{d_{ik}}\left[\left(\left(\textbf{J}_{pq}-\textbf{J}_{qp}\right)\textbf{D}_{k}\left(\textbf{J}_{ml}-\textbf{J}_{lm}\right)\right)_{ii}\right.$
$\displaystyle\hphantom{=\frac{1}{K}}\left.+\left(\textbf{D}_{k}\left(\textbf{J}_{ml}-\textbf{J}_{lm}\right)\left(\textbf{J}_{qp}-\textbf{J}_{pq}\right)\right)_{ii}\right].$

The second term in the inner summand is easiest to simplify further. Notice
that by properties of the single-entry matrix, this term is zero, except when
$l=p$ and $m=q$. That is, when we consider the second-order derivative of the
loss function with respect to $\varepsilon_{lm}$. In that particular case, the
second term is given by
$-\left(\textbf{D}_{k}\textbf{J}_{mm}+\textbf{D}_{k}\textbf{J}_{ll}\right)_{ii}$.
For the first term, things are slightly more involved. Yet, here, we can show
that it equals zero, again except when $l=p$ and $m=q$, in which case it
equals the following expression:
$\left(\left(d_{mk}-\lambda\right)\textbf{J}_{ll}+\left(d_{lk}-\lambda\right)\textbf{J}_{mm}\right)_{ii}$.

Now, considering index $i$, for all values of $i$ other than $m$ and $l$, both
terms are also zero, even if $l=p$ and $m=q$. Consequently, we can write the
second-order derivative with respect to $\varepsilon_{lm}$ as follows:

$\displaystyle\left.\frac{\partial^{2}\mathcal{L}\left(\textbf{E}\right)}{\partial^{2}\varepsilon_{lm}}\right|_{\textbf{E}=\textbf{0}}$
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\sum_{i=1}^{N}\frac{1}{d_{ik}}\left[\left(\left(d_{mk}-\lambda\right)\textbf{J}_{ll}+\left(d_{lk}-\lambda\right)\textbf{J}_{mm}\right)_{ii}\right.$
(26)
$\displaystyle\hphantom{=\frac{1}{K}}\left.-\left(\textbf{D}_{k}\textbf{J}_{mm}+\textbf{D}_{k}\textbf{J}_{ll}\right)_{ii}\right]$
(27)
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\frac{\left(d_{mk}-\lambda\right)-\left(d_{lk}-\lambda\right)}{d_{lk}}+\frac{\left(d_{lk}-\lambda\right)-\left(d_{mk}-\lambda\right)}{d_{mk}}$
(28)
$\displaystyle=\frac{1}{K}\sum_{k=1}^{K}\frac{d_{mk}}{d_{lk}}+\frac{d_{lk}}{d_{mk}}-2.$
(29)

In short, when the solution is very close to being diagonal, the Hessian
becomes a diagonal matrix, with the diagonal element that corresponds to the
second derivative with respect to parameter $\varepsilon_{lm}$ as described in
the last equation. Thus, by dividing the corresponding elements in the
gradient by these elements, we obtain a quasi-Newton update.

Defining matrix H with element $l,m$ as seen in Equation 29, the search
direction in this minimization problem can now be set as follows:

$\textbf{E}=-\textbf{G}\oslash\textbf{H},$ (30)

where ‘$\oslash$’ denotes the elementwise division, of $N\times N$ matrices G
and H, comprising corresponding elements of the gradient and Hessian
respectively.

Each element of H can be described as a function of some scalar
$x\in[0,+\infty)$, viz., $x+x^{-1}-2$. It is easy to show this function
attains its minimum at $x=1$, where then the value is exactly zero. Thus, the
elements of H are non-negative by definition. To ensure numerical stability,
we set elements of H below threshold $\tau_{H}=0.01$ to be equal to
$\tau_{H}$. This modification to H is equivalent to requiring that some
element in E can have a magnitude that is at most 100 times as large as the
corresponding element in G.

This overall approach for calculating G, H, and E is similar to what is used
in an undocumented part of the Python package qndiag [Ablin et al., 2018].
However, our derivations (_i_) formalize and prove the correctness of
expressions that are similar to those utilized in that code, (_ii_) provide
expressions that are formulated such that greater computational speedups are
possible, and (_iii_) take our low-rank approximation with regularization
fully into account.

###  2.4 Line search

The preceding quasi-Newton approach formulates the update in terms of matrix
E. In a standard line search, one considers $\alpha\in[0,1]$, where $\alpha=0$
would correspond to no update and $\alpha=1$ would correspond to the full
update. In our formulation, the update has a multiplicative nature.

That is, define

$\textbf{R}^{*}=e^{\textbf{E}-{\textbf{E}}^{\top}}.$ (31)

Now, considering $\textbf{A}_{(t-1)k}$ for $k=1,\ldots,K$ (readily available
from the previous iteration), if $\alpha=1$, then $\textbf{A}_{(t)k}$ would be
given by $\textbf{R}^{*}\textbf{A}_{(t-1)k}$. Thus, given E and
$\textbf{A}_{(t-1)k}$, we have that $\textbf{R}^{*}$ and
$\textbf{R}^{*}\textbf{A}_{(t-1)1},\ldots,\textbf{R}^{*}\textbf{A}_{(t-1)K}$
altogether can be calculated in $O(N^{3})$ time. Moreover, if $\alpha=0$, we
have that $\textbf{A}_{(t)k}$ is simply given by its value in the previous
iteration.

By taking a convex combination of these two end-points, we ‘linearize’ the
update during the line search. That is, we set

$\textbf{A}^{*}_{k}=\alpha\textbf{R}^{*}\textbf{A}_{(t-1)k}+(1-\alpha)\textbf{A}_{(t-1)k}.$
(32)

Now, given $\textbf{R}^{*}\textbf{A}_{(t-1)k}$ and $\textbf{A}_{(t-1)k}$ for
$k=1,\ldots,K$, we can calculate the loss function as defined in Equation 7 in
$O(KN^{2})$ time for each value of $\alpha$ considered. This approach allows
us to efficiently apply a golden-section search to find the best value of
$\alpha$.

Finally, we observe that as a result of the linearization, for
$\alpha\in(0,1)$ its value is slightly underestimated. Therefore, given
$\alpha$ after the linearization, we cast it back to the exponential scale by
setting

$\alpha^{*}=\log\left(1+\left(\alpha(e-1)\right)\right).$ (33)

Finally, we set

R $\displaystyle=e^{\alpha^{*}\left(\textbf{E}-{\textbf{E}}^{\top}\right)}$
(34) $\displaystyle\textbf{B}_{(t)}$
$\displaystyle=\textbf{R}\textbf{B}_{(t-1)}$ (35)
$\displaystyle\textbf{A}_{(t)k}$
$\displaystyle=\textbf{R}\textbf{A}_{(t-1)k},$ (36)

after which we move to the next iteration, provided convergence has not yet
occurred.

###  2.5 Initialization

JADOC initializes by setting $\textbf{B}_{(0)}=\textbf{I}_{N}$, which clearly
is a valid orthonormal transformation matrix.

###  2.6 Convergence

Convergence is defined in terms of the root-means-square deviation (RMSD) from
zero of the gradient. By default, a value below $10^{-4}$, in conjunction with
JADOC having made at least $T_{min}=10$ iterations, is considered to
constitute convergence. The latter requirement is to deal with the fact that
the gradient can be quite small near the starting point (i.e., in the first
few iterations) due to poor starting values. After $T=100$ iterations JADOC
will terminate, irrespective of whether convergence has occurred.

##  3 Algorithms

Algorithm 1 shows an overview of the JADOC algorithm as implemented in our
Python code, available at https://github.com/devlaming/jadoc. Within the
iterations (i.e., the last for loop), updating B and $\textbf{A}_{k}$ for
$k=1,\ldots,K$ (i.e., the last step) is the most demanding step,
computationally speaking, causing the algorithm to require $O(N^{3})$ per
iteration. This, however, still constitutes an $O(K)$ reduction in time
complexity, compared to existing algorithms for joint diagonalization.

Algorithm 2 shows an overview of the data-generating process we use in our
simulations, to compare performance of JADOC to existing algorithms. This
algorithm is implemented in the SimulateData(iK,iN,iR,dAlpha) function, as
implemented in the Python package on GitHub.

Algorithm 1 JADOC algorithm.

Input: $N\times N$ positive (semi)-definite matrices, $\textbf{C}_{k}$, for
$k=1,\ldots,K$
Output: $N\times N$ orthonormal diagonalizer matrix

Set $\tau_{conv}=10^{-4}$, $\tau_{H}=10^{-2}$,
$S=\left\lceil\frac{N}{K}\right\rceil$, $\lambda_{0}=1$,
$\textbf{B}_{(0)}=\textbf{I}_{N}$, $T=100$, $T_{min}=10$

for $k=1,\ldots,K$ do

Compute $\textbf{L}_{k}$ in line with Equation 4 and store as
$\textbf{A}_{(0)k}$

Update $\lambda$ in line with Equation 6

end for

for $t=1,\ldots,T$ do

Compute gradient in accordance with Equation 18 $\rightarrow\textbf{G}$

if root-mean-square deviation of gradient $<\tau_{conv}$ and $t>T_{min}$ then

break for loop

end if

Compute diagonal of approximated Hessian using Equation 29
$\rightarrow\textbf{H}$

Set elements of H below $\tau_{H}$ equal to $\tau_{H}$

Compute update using Equation 30 $\rightarrow\textbf{E}$

Apply line search described in Subsection 2.4 $\rightarrow\alpha^{*}$

Computation rotation matrix using Equation 34 $\rightarrow\textbf{R}$

Update B and $\textbf{A}_{k}$ using Equations 35 & 36
$\rightarrow\textbf{B}_{(t)},\textbf{A}_{(t)k}$

end for

return $\textbf{B}_{(t)}$

Algorithm 2 Algorithm for the data-generating process of our simulations.

Input: $K$, $N$, replicate number $r$, degree of similarity eigenvectors
$\alpha\in[0,1]$
Output: $K$ different $N\times N$ PSD matrices

Set random seed based on $r$

Get $N^{2}$ iid draws from $\mathcal{N}\left(0,1\right)$ and store in $N\times
N$ matrix $\rightarrow\textbf{X}$

for $k=1,\ldots,K$ do

Get $N^{2}$ iid draws from $\mathcal{N}\left(0,1\right)$ and store in $N\times
N$ matrix $\rightarrow\textbf{X}_{k}$

Compute
$\alpha\textbf{X}+\left(1-\alpha\right)\textbf{X}_{k}\rightarrow\textbf{X}_{k}$

Compute $e^{\textbf{X}_{k}-{\textbf{X}}^{\top}_{k}}\rightarrow\textbf{R}$

Get matrix of zeros $\rightarrow\textbf{D}$

for $n=1,\ldots,N$ do

Draw from $\chi^{2}\left(1\right)\rightarrow d$

Set $n$-th diagonal element of D equal to $d$

end for

Compute $\textbf{R}\textbf{D}{\textbf{R}}^{\top}\rightarrow\textbf{C}_{k}$

end for

return $\textbf{C}_{1},\ldots,\textbf{C}_{K}$

##  4 Simulation results

We consider two simulation designs. In Design 1, we set set $K=10$ and
consider $N=100,200,300,400,500$. In Design 2, we set $N=256$ and consider
$K=2,4,8,16,32$. In each design, for a given value of $N$ and $K$, we consider
$R=10$ replicates.

For each of these designs, we consider four values of the parameter $\alpha$
defined in Algorithm 2, viz., $\alpha=0,0.25,0.5,0.75$. Here, $\alpha=0$
corresponds to $\textbf{C}_{1},\ldots,\textbf{C}_{K}$ having completely
independent eigenvectors. Values of $\alpha$ between zero and one correspond
to intermediate scenarios, where eigenvectors across matrices are correlated;
the higher $\alpha$, the stronger this correlation.

In total, we have eight scenarios: two designs and four different values of
$\alpha$. In each scenario, we apply

  1. 1.

JADE [Cardoso and Souloumiac, 1996] as implemented on
https://github.com/gabrieldernbach/approximate_joint_diagonalization,

  2. 2.

qndiag [Ablin et al., 2018] as implemented on
https://github.com/pierreablin/qndiag, and

  3. 3.

our JADOC algorithm as implemented on
https://github.com/devlaming/jadoc.

For each of the ten replicates and each of the three methods, we run the
relevant Python code on separate machines with identical configurations (i.e.,
each machine has 24 cores with a clock speed of 2.6GHz and 64GB of memory).
That is, for a given replicate and a given method, all ten scenarios are
carried out on a single, separate machine. For JADOC and JADE, we make use of
the built-in parallelized functionality. Finally, for fair comparison, for
qndiag, we use the built-in option to enforce orthonormality of B.

For each scenario and for each method, we report the median (across the ten
replicates) of the runtimes (solid lines) and the median (across the ten
replicates) of the root-mean-square deviation from zero of the off-diagonal
elements after transformation using the matrix B (dotted lines).

$\alpha=0.00$
---
|
$\alpha=0.25$
|
$\alpha=0.50$
|
$\alpha=0.75$
|
Figure 1: Runtime in seconds (solid lines with circles) and root-mean-square
deviation (RMSD) from zero of off-diagonal elements (dotted lines with
diamonds). Shown is the median across 10 replicates, for various sample sizes
$N$ with $K=10$ (left column) and for various number of matrices to
diagonalize $K$ with $N=256$ (right column), for various degrees of
similarities in eigenvectors across matrices $\alpha$ (rows). Each plot shows
lines for JADOC (black), JADE (dark gray), and qndiag (light gray).

Figure 1 comprises eight line plots (one for each scenario), where each plot
shows the median runtime of each algorithm on the left $y$-axis and the median
RMSD of each algorithm on the right $y$-axis. Each plot has three lines,
corresponding to JADOC, JADE, and qndiag. In the plots on the left (Design 1)
$N$ is shown on the $x$-axis and in the plots on the right (Design 2) $K$ is
shown on the $x$-axis. In the consecutive rows, results are shown for
$\alpha=0,0.25,0.5,0.75$.

We see that JADOC outperforms qndiag by several orders of magnitude in terms
of CPU time both as $N$ and as $K$ increase. In turn, qndiag outperforms JADE
by another order of magnitude timewise. Overall, JADOC is much fast than
qndiag which in turn is much faster than JADE. Also, as expected, for fixed
$N$ we see the runtime of JADOC does not increase $K$. In fact, the runtime of
JADOC even seems to decrease slightly as $K$ increases.

In terms of the degree of diagonalization, as measured by the RMSD, we see
that JADE performs best, very closely followed by JADOC, and finally followed
at some distance by qndiag. Thus, in terms of the degree of diagonalization,
JADOC closely trails the golden standard provided by JADE, whereas qndiag
provides a considerably poorer result.

Overall, the improved runtime afforded by JADOC simply dwarfs the minor
increase in the degree of diagonalization that is achieved by JADE.

##  5 Conclusions

We have derived a new framework for finding a square matrix B such that
$\textbf{B}\textbf{C}_{k}{\textbf{B}}^{\top}$ is as diagonal as possible for
$k=1,\ldots,K$, where $\textbf{C}_{1},\ldots,\textbf{C}_{K}$ are symmetric,
positive (semi)-definite matrices, where
$\textbf{B}{\textbf{B}}^{\top}=\textbf{I}_{N}$ (i.e., B is orthonormal). Our
approach for finding B we refer to as JADOC: Joint Approximate Diagonalization
under Orthogonality Constraints.

In a nutshell, JADOC (_i_) applies a simple dimensionality-reduction technique
combined with mild regularization, (_ii_) uses a quasi-Newton approach which
is combined with a golden section, to jointly update B in each iteration,
(_iii_) guarantees an orthonormal B, and (_iv_) requires only $O(N^{3})$ time
per iteration. JADOC is implemented as an open-access tool for Python 3.$x$
and is available on GitHub (see: https://github.com/devlaming/jadoc).

In most real-life scenarios, positive (semi)-definite matrices
$\textbf{C}_{1},\ldots,\textbf{C}_{K}$ will not be jointly diagonalizable. In
other words, their eigenvectors will not be completely identical. Under this
scenario, we find that JADOC improves upon two important approaches, viz., the
classical JADE method and a recent quasi-Newton approach with a Python package
called qndiag. JADOC has far superior runtime compared to both. Moreover, in
terms of the degree of diagonalization, JADOC closely follows the golden
standard set by JADE.

Finally, we observe, in line with our theory, that for fixed $N$, increasing
$K$ does not increase the computational complexity of JADOC. Its time
performance remains steady for all values of $K$ considered here. In fact, we
even see a slight decrease in CPU time as $K$ increases, as the
parallelization implemented in JADOC is tailored towards higher values of $K$.

We conclude that JADOC is an efficient and accurate tool to diagonalize $K$
different $N\times N$, symmetric, positive (semi)-definite matrices, requiring
only $O(N^{3})$ time per iteration, which constitutes an $O(K)$ reduction in
complexity, compared to the $O(KN^{3})$ time required in each iteration by
competing methods. This reduction in computational complexity makes JADOC far
more scalable and, therefore, much more useful than its competitors in the age
of big data.

## Acknowledgments

We would like to acknowledge the fruitful discussions we have had about this
topic with Niels Rietveld, Robert Kirkpatrick, and Patrick Groenen. This work
was carried out on the Dutch national e-infrastructure with the support of
SURF Cooperative (NWO Call for Compute Time EINF-403 to E.A.W.S.).

## References

  * [Ablin et al., 2018] Ablin, P., Cardoso, J.-F., and Gramfort, A. (2018).  Beyond pham’s algorithm for joint diagonalization.  arXiv preprint, arXiv:1811.11433.
  * [Bröcker and Tom Dieck, 1985] Bröcker, T. and Tom Dieck, T. (1985).  Representations of compact Lie groups, volume 98.  Springer-Verlag Berlin Heidelberg.
  * [Cardoso and Souloumiac, 1996] Cardoso, J.-F. and Souloumiac, A. (1996).  Jacobi angles for simultaneous diagonalization.  SIAM journal on matrix analysis and applications, 17(1):161–164.
  * [Gallier and Xu, 2002] Gallier, J. and Xu, D. (2002).  Computing exponentials of skew-symmetric matrices and logarithms of orthogonal matrices.  International Journal of Robotics and Automation, 17(4):1–11.
  * [Kiefer, 1953] Kiefer, J. (1953).  Sequential minimax search for a maximum.  Proceedings of the American mathematical society, 4(3):502–506.
  * [Marsden and Ratiu, 1999] Marsden, J. E. and Ratiu, T. S. (1999).  Introduction to mechanics and symmetry: a basic exposition of classical mechanical systems, volume 17.  Springer-Verlag New York.
  * [Mesloub et al., 2014] Mesloub, A., Abed-Meraim, K., and Belouchrani, A. (2014).  A new algorithm for complex non-orthogonal joint diagonalization based on shear and givens rotations.  IEEE Transactions on Signal Processing, 62(8):1913–1925.
  * [Pham, 2001] Pham, D. T. (2001).  Joint approximate diagonalization of positive definite hermitian matrices.  SIAM Journal on Matrix Analysis and Applications, 22(4):1136–1152.