Joint Approximate Diagonalization of Positive Definite Hermitian Matrices.
This paper provides an iterative algorithm to jointly approximately diagonalize K Hermitian positive definite matrices ${\bf\Gamma}_1$, \dots, ${\bf\Gamma}_K$. Specifically, it calculates the matrix B which minimizes the criterion $\sum_{k=1}^K n_k [\log \det\diag(\B\C_k\B^*) - \log\det(\B\C_k\B^*)]$, nk being positive numbers, which is a measure of the deviation from diagonality of the matrices BCkB*$. The convergence of the algorithm is discussed and some numerical experiments are performed showing the good performance of the algorithm.
Introduction
The problem of diagonalizing jointly approximately several positive definite
matrices has arisen in (at least) two different contexts. The first one is the statistical
problem of common principal components in k group introduced by Flury (1984).
He considers k populations of multivariate observations of size n 1 obeying
the Gaussian distribution with zero means and covariance matrices
and assume that \Gamma k can be written as B     k B for some orthogonal matrix B and
some diagonal matrices   k (the symbol   denotes the transpose). The problem is to
estimate the matrix B (the colum of which are the common principal components)
from the sample covariance matrices C k of the populations. As it is well
known that n k C k , are distributed independently according to the
Wishart distribution of n k degrees of freedom and covariance matrices \Gamma k (see for
example Seber, 1984), the log likelihood function of based on them is
where C is a constant term and tr denotes the trace. Therefore the log likelihood
method for estimating B and the   k amounts to minimizing
For fixed B, it is not hard to see that the above expression is minimized with respect
to the   k when   the notation diag(M) denoting the diagonal
matrix with the same diagonal as M. Thus, one is led to the minimization (with
respect to B) of
which is the same as that of
since B has unit determinant. But (1.1) is precisely a measure of the global deviation
of the matrices  from diagonality, since, from the Hadamard inequality (Noble
and Daniel, 1977, exercise 11.51), det M - det diag(M) with equality if and only if
M is diagonal. Thus minimizing (1.1) can be viewed as trying to find a matrix B
which diagonalizes jointly the matrices C much as it can.
More recently, several authors (Cardoso and Souloumiac, 1993, Belouchrami
et al., 1977, Pham and Garat, 1997) have introduced the joint approximate diagonalization
as a method for the separation of source problem. In this problem, there
are K sensors which record each a linear mixture of K sources, so that denoting
by X(t) and S(t) the vectors of measurements and of sources at time t, one has
A. The goal is to extract the sources from
the observations and in the so called blind separation one does not have any specific
knowledge about the sources other than that they are statistically independent. Thus
a sensible method is to try to find a matrix B such that the components of BX(t)
(which represent the reconstructed sources) are as independent as possible. As it
is easier to work with non correlation rather than independence, a simple method
would be to try to make the cross-correlation, eventually lagged, between the sources,
vanish. This would lead to the joint approximate diagonalization of a certain set of
covariance matrices, as proposed in Belouchrami et al. (1977). Note that Pham
and Garat (1997) also consider joint diagonalization but they use only two matrices
and then the diagonalization can be exact (see for ex. Golub and Van Loan, 1989).
Cardoso and Souloumiac (1993), on the other hand, do not consider lagged covariance
but use higher order cumulants between the sources instead. They construct
certain set of matrices in which such cumulants appear as off diagonal elements. The
separation of source is then solved through a joint approximate diagonalization of
these matrices.
It should be pointed out that the above authors use a different measure
of deviation to diagonality than that of Flury. Their measure is simply the sum
of squares of the off-diagonal elements of the considered matrices. But there is a
common feature in all above works in that the diagonalizing matrix B is taken to be
orthogonal. In this work we shall drop this restriction. The orthogonality condition
is part of the assumption of Flury (1984) but there is no clear reason that it should
be satisfied. This condition is justified in the works of Cardoso and Souloumiac
(1993) and Belouchrami et al. (1997), since these authors have pre-normalized their
observations to be uncorrelated and have unit variance. We want to avoid this pre-
normalizing stage in the separation of source procedure, which can adversely affect its
performance since the statistical error committed in this stage cannot be corrected in
the following "effective separation" stage. By dropping the orthogonality restriction,
we obtain a single-stage separation procedure which is simpler and can perform
better. Note that without the orthogonality restriction, exact joint diagonalization
is possible for two matrices (see for ex. Golub and Van Loan, 1989). But for
more than two one can only achieve approximate joint diagonalization, relative to
some measure of deviation to diagonality. We take this measure to be (1.1) for two
following reasons. Firstly, it can be traced back to the likelihood criterion, widely
used in statistics. Secondly, this criterion is invariant with respect to scale change: it
remains the same if the matrices to be diagonalized are pre- and post-multiplied by
a diagonal matrix. The other measure which consists in taking the sum of squares of
the off-diagonal elements of the matrices, does not have this nice invariant property.
Of course, one can introduce this property by first normalizing the matrices so that
they have unit diagonal element, but then the resulting criterion would be very hard
to manipulated.
The main result of this paper is the derivation of an algorithm to perform
the joint approximate diagonalization in the sense of the criterion (1.1) and without
the restriction that the diagonalizing matrix be orthogonal. Our algorithm has
some similarity with that of Cardoso and Souloumiac (1993) and even more with
that of Flury and Gautschi (1996): it also operates through successive transformations
on pairs of rows and columns of the matrices to be diagonalized. However the
convergence proof is completely different since we can no longer rely on the orthogonality
property. Incidentally our method of proof can be easily adapted to prove
the convergence result in Flury and Gautschi (1996), in a much simpler way.
2 The algorithm
As one frequently encounters complex data in signal processing applications,
we shall consider complex Hermitian (instead of real symmetric) positive definite
matrices . (Note that Cardoso and Soulomiac (1993) and Bellouchrani
et al. (1997) also work in a complex setting.) The goal is to find a complex matrix
B such that the matrices are as close to diagonal as possible,
the notation   now denoting the transpose complex conjugated. The measure of
deviation to diagonality is taken to be (1.1) where the n k are positive weighs (they
need not be integers). Note that since the C k do not depend of B, the minimization
criterion (1.1) can be reduced to
The algorithm consists in performing successive transformations, each time
on a pair of rows of B, the i-th row B i\Delta and the j-th row say, according to
a b
c d
and in such a way that the criterion is sufficiently decreased. Whether the decrease
is sufficient is a question which we shall returned in next section. Once this is done,
the procedure is repeated with another pair of indices until convergence is achieved.
Denote by u
ij the general element of the matrix the decrease of
the criterion associated with the transformation (2.1) is
du
ii )=(u (k)
ii u (k)
A natural idea is to chose a, b, c, d to maximize the above decrease. However, this
maximization cannot be done analytically. Our idea is to maximize a lower bound of
it instead. Since the logarithm function is convex, for any two sequences of positive
1, one has
Applying this inequality, the above decrease is bounded below by
ii
ii )=u (k)
. Introducing the matrices
ii
=n
this lower bound can be rewritten as
[a b]P
- a
log
[c d]Q
d
The maximization of (2.3) can be done analytically as it will be shown below. Since
(2.3) vanishes at a = clear that its maximum
is non negative and can be zero only if this maximum is also attained a
Thus the decrease of the criterion, associated with the transformation
using the values of a, b, c, d realizing the maximum of (2.3), is positive unless
attains it maximum at a =
Let us return to the maximization of (2.3). For any given a
can always parameterize a; b; c; d as
a b
c d
a
provided that the last matrix is non singular. Put
a
a
a
a
one can express (2.3) as
log ja 0 d
The first term of the above right hand side is no other than (2.3) evaluated at the
point (a Therefore, a necessary and sufficient condition that this point
realizes the maximum of (2.3) is that the last term in the above right side is non
negative, for all -, ffl, j. But for this term can be seen to
be equivalent to n[(-fflp 0
Hence a necessary condition for
this condition to hold is that This is the same that the matrices P and
Q be jointly diagonalized by
a
. Further, for the right hand side
of (2.6) reduces to
Again, for 0, the last term in the above expression can be
seen to be equivalent to jfflj 2 p 0
j. Therefore it is also necessary
that this quadratic form in the variables ffl; - j be non negative. This requirement
is satisfied if and only if p 0
2 . On the other hand, using the inequality
can be seen to be bounded above by
log ja 0 d
But for p 0
2 , the quadratic form j-fflj
(in
the variables -
-ffl; -j) is non negative, entailing that the above expression is bounded
above by n log[ja
choices of -; ffl; j. Thus we have proved
that a necessary and sufficient condition for a to realize the maximum
of (2.3) is that
a
jointly diagonalizes P and Q and that the diagonal terms
2 of the diagonalized matrix satisfy p 0
. (Note that if the last
condition is not satisfied, one simply needs to permute a
In the appendix, the problem of jointly diagonalizing two Hermitian matrices
P and Q of size two (not necessarily positive definite) is completely solved. It is
shown there that if P and Q are positive definite and not proportional, the solution
exists and is unique up to a permutation and a scaling, that is all solutions can be
obtained from a representative one by pre-multiplying it by a diagonal matrix and a
permutation matrix. Denoting by the diagonal and the upper
off diagonal elements of P and Q and putting
a representative solution is given by
2ff D
and is positive.
We shall prove below that (i) p 2 q 1 - 1 with equality if and only if the
matrices P and Q are proportional and (ii) the p 0
defined by (2.5) with
a are such that p 0
1 has the same sign as
As obviously the above results show that the solution
to the maximization of (2.3) is given by [a b] / [D 2fl], [c d] / [2ff D], the
meaning proportional to. Of course, this result doesn't apply in the case
where the two matrices P and Q are proportional. But in this case, diagonalize one
would diagonalize the other. As it has been proved before, this would be enough to
ensure that (2.3) be maximized. Thus in this case, there exists an infinite number
of solutions (even after eliminating the ambiguity associated with scaling).
Note
The Flury and Gautschi (1996) algorithm operates on a similar principle.
However, these authors iterate the transformation (2.1) with a fixed pair (i;
convergence and only then he changes to another pair. We feel that this is less
efficient, because by using the same pair, the decrease of the criterion tends to be
smaller each time while by changing it one can get big decrease in the first few
iterations. Our algorithm is also simpler to program.
We now proved the results (i) and (ii) announced above. By formula (2.2),
one has
l
ii
ii
l
ii
ii
ii
ii
l
hi u (k)
ii
ii
ii
ii
It follows that p 2 q 1 ? 1 with equality if and only if u (1)
ii =u (1)
ii =u (K)
jj . In
the last case P and Q are proportional. This proves the result (i).
On the other hand, for a one gets from (2.5)
Hence the product p 0
denoting the real part. The product p 0
2 can be obtained from the above formula
by interchanging Therefore, the difference p 0
The last three terms of the above expression may be regrouped as
D)]:
putting
\Delta, one has
and hence, noting that
fi)=2:
Further
Therefore, combining the above result, one gets
fi)]g:
As
fi is purely imaginary, the last term in the above expression equals
One thus obtains finally p 0
which has the
same sign as
Note
The computation of ff, fi and fl could be subjected to large relative error if
the matrices P and Q are nearly proportional. But this doesn't matter as long as
the matrices P and Q are diagonalized with sufficient accuracy so that the criterion
be adequately decreased. To this end, note that the solution to the problem remains
the same if one replaces Q by arbitrary ae. One can chose ae
such that when P and Q are nearly proportional, R is almost zero. The numerical
calculation of R will then be subjected to large relative error, but this is the only
main source of error since there will be no near cancellation is subsequent calculation.
More precisely, let ~
R be the calculated R, the algorithm would diagonalize P and ~
R.
Hence it would diagonalize ~
the absolute error ~
small (it is the relative error which is large), P and Q are still accurately diagonalized.
The above argument can be repeated with the role of P and Q interchanged, that
is one takes jointly diagonalizes Q and R. This alternative would
be preferable, from the numerical accuracy point of view, if it leads to a smaller (in
absolute value) of ae. A simple rule to chose ae is to require that it results in a zero
diagonal element of R and is as small as possible. Thus
is defined as is defined as
we want ae to be small, we chose the first possibility if q 1 and the second
otherwise. We shall assume that as it is the case here. (But the case
can be handled in a similar way). Thus the diagonal and upper
off diagonal elements and r of R are
and one jointly diagonalizes P and R in the first case and Q and R in the second
case. Hence in the first case, we compute ff, fi, fl as
and in the second case as
The second right hand sides in the above formula provide more efficient computations
without loss of accuracy. Further, they remain applicable even in if the matrices
P and Q are proportional, by replacing
by an arbitrary non zero number (ff in the first case or fl in the second case can
be taken arbitrary too). Indeed the above computation ensures that in the first
case
this is sufficient to entail that the matrix
P is diagonalized (see the Proof of Proposition A1), and hence so is Q as it is
proportional to P. Similarly, in the second case, the above computation ensures
that Q is diagonalized and hence so is P.
3 Convergence of the algorithm
We have shown in previous section that our algorithm decreases the criterion
at each step, unless (2.5) is maximized at a = But from the results
of this section, this implies that the matrix P and Q are diagonal (we already proved
that 1). If this occurs for a pair of indexes (i; j) then one would skip
this pair and continue the algorithm with another pair. But if this occurs for all
pairs, then the algorithm stops. Explicitly, it stops when
One may recognize that the above condition is no other than the condition that B be
a stationary point of the criterion (1.1). Indeed, consider a small change in B of the
form ffiB (hence the matrix ffi represents a relative change), then the corresponding
change of (1.1) is
log[(u
ri
r
rs -
denoting the general elements of ffi. Expanding this expression with respect to ffi,
up to the first order, one gets
Thus, the vector with components 2n- ij , i 6= j can be viewed as the relative gradient
vector of the criterion (1.1).
We now prove that the decrease in the criterion at each step of the algorithm
is sufficient to ensure the convergence to zero of the above gradient vector. As it
has been proved in previous section, by parameterizing a; b; c; d as in (2.4) with
a being the point realizing the maximum of (2.3), one can express (2.3)
as
are defined by (2.5). Take -; ffl; j such that the
left hand side or (2.4) is the identity matrix. Then (2.3) vanishes. Thus (2:6 0 ) must
also vanish, hence its upper bound derived in section 2 must be positive. Therefore,
noting that p 0
, one has
log ja 0 d
where -; ffl; j; -, are the elements of the inverse of the matrix
a
. The left
hand side of (3.4) is no other than the value of (2.3) at a . Further, the
transformation (2.1) for this step of the algorithm uses precisely these a
Thus the decrease of the criterion associated with this transformation, which is
bounded below by (2.3), must be at least as large as the left hand side of (3.4). But
from the definition of -; ffl; -; j and (2.5), one has
hence, noting that the right hand side of (3.4) can be seen to equal
The last expression can be rewritten as
Therefore, noting that the middle matrix in this quadratic form has eigenvalues 1 \Sigma ae
and 0 - ae - 1, it is bounded below by n[(q 0
This is also
the lower bound of the decrease of our criterion at this step.
Since the criterion is always decreased during our algorithm, it must converge
to a limit. Therefore the decrease of the criterion at each step of the algorithm must
converge to zero, implying that (q 0
tends to zero. Note that
p, q are no other than - ij , - ji defined in (3.1) and 2n- ij are the components of
the relative gradient vector at this step of the algorithm. Still, the above result
proved the convergence to zero of this vector. The difficulty is due
to the lack of normalization. Indeed, our algorithm constructs the transformation
only up to a scaling of its rows, hence a row of B can be arbitrary large or
arbitrary small and this has an effect on the gradient, even when relative gradient
is considered. To avoid this, we shall renormalize the transformation matrices B
after each step of the algorithm. Any reasonable normalization procedure will do
but for simplicity and definiteness, we will consider the normalization which makes
the rows of B having unit norm. Then u
ii will be bounded between the smallest
and the largest eigenvalue of C k . Therefore, letting m and M be the minimum of
the smallest eigenvalues and the maximum of the largest eigenvalues, of C
respectively, one has m - u
for all for all i all k. Note that m ? 0 since the
matrices are all positive definite. Therefore, from (2.2) and (2.5),
[a
mn
a 0
Mn
a 0
and the same inequalities hold for q 0
1 , and similar inequalities, with a 0 , b 0 replaced by
2 . Thus both p 0
2 can be bounded above by M=m and
below by m=M . It follows that the relative gradient vector of the criterion evaluated
at each step of the algorithm, which has components 2n- ij , converges to zero.
The above result shows that if the algorithm converges, then the limit must
be a stationary point of the criterion. Further, since the algorithm always decreases
the criterion, this point is actually a local minimum, unless the algorithm is started
at a stationary point, in which case it stops immediately. Note that, the sequence of
transformation matrices constructed by the algorithm, being normalized and hence
lying on a compact set, will admit a convergent subsequence and this in fact also
holds for any of its subsequence. Therefore, if the criterion admit an unique local
minimum, the algorithm will converge to it. However, Flury and Gautschi (1996)
has shown that in some extreme cases, the minimization of the criterion (1.1) under
the orthogonality constraint, admits more than one local minimum. Therefore, it
seems likely that in our problem, where the maximization is without constraint, the
uniqueness of the local minimum is also not satisfied in all cases. Nevertheless, if
there are only several local minima, one can still expect that the algorithm converge
to one of them. Indeed, if this is not so, then since we have proved that the gradient
vector converge to 0, the algorithm must jump continually from one local minimum
to another, a quite implausible thing. The existence of a finite number of local
clearly does not hold if the matrices C are all proportional to a
single matrix, but this is a very extreme case.
We conclude this section by showing that our algorithm behaves near the
solution very much like to the Newton-Rhapson iteration, provided that the matrices
can be nearly jointly diagonalized. To derive the Newton-Rhapson iter-
ation, one makes a second order Taylor expansion of the criterion around the current
point, then minimizes this expansion (instead of the true criterion) to obtain the new
point. As we have already computed the change of the criterion corresponding to
a change ffiB of B, resulting in the formula (3.2), we need only to expand it up to
second order in ffi. Note that the first order expansion has already been given by
(3.3), we need only to pursue the expansion up to second order, yielding
r
s
rs =u (k)
ii
r
s
ir =u (k)
ii )!(ffi is -
is =u (k)
Assume that the matrices C can be nearly jointly diagonalized, then near
the solution, the off diagonal term u (k)
rs , r 6= s and u (k)
ir , r 6= i, of the matrices
would be small relative to the diagonal term u (k)
ii . Hence we may neglect,
in the above expression the term of second order in ffi containg the factor u (k)
rs =u (k)
ii ,
r 6= s or u (k)
ir =u (k)
ii , r 6= i. With this approximations, the above expansion for (3.2)
(up to second order in reduces to
The (approximate) Newton-Rhapson algorithm consists in minimizing (3.2')
with respect to ffi, then change B into B being the solution to the above
minimization. Note that (3:2 0 ) can be written as the sum over all (unordered) pairs
The minimization of this expression with respect to [ffi ij
can be easily done,
yielding -
It is worthwhile to note that the diagonal elements of ffi do not appear
in (3:2 0 ) so they can be anything as long as they are small. For convenience, we
take them to be 0. This is further justified by the fact that by dividing the i-th row
of B+ ffiB by 1 one is led to the matrix
0 has zero diagonal
element and (i; diagonal element which is about the same as
Note also that because the ffi ij are small, the matrix B + ffiB can be obtained
through successive transformations of the form (2.1) with a
associated with all distinct pairs of indexes (i; j), i 6= j. Reverting to the notation
defined by (2.2), are no other than p 2 , q 1 , p and -
q. Thus
On the other hand, for small p and
q, fi and D, as defined in (2.7) and (2.8), can both be approximated by p 1
Hence, since one gets that b and c equal approximately 2fi=(fi +D) and
+D). The above results show that the new matrix B resulting from one step
Newton-Rhapson algorithm is about the same as that resulting from a "sweep" of the
algorithm of section 2. constituted by successive steps associated with
all distinct pairs of indexes.) Threfore, our algorithm has about the same quadratic
convergence speed of the Newton-Rhapson iteration, near the solution. But the
Newton-Rhapson method may converge badly, even not at all, if it is started at a
point far from the solution. Our algorithm would have better convergence behavior
in this case since it always decreases the criterion.
4 Some numerical examples
We consider the same example as in Flury and Gautschi (1996). The following
6 \Theta 6 matrices are to be diagonalized
\Gamma12:5 27:5 \Gamma4:5 \Gamma:5 2:04 \Gamma3:72
\Gamma:5 \Gamma4:5 24:5 \Gamma9:5 \Gamma3:72 \Gamma2:04
\Gamma4:5 \Gamma:5 \Gamma9:5 24:5 3:72 2:04
\Gamma2:04 2:04 \Gamma3:72 3:72 54:76 \Gamma4:68
3:72 \Gamma3:72 \Gamma2:04 2:04 \Gamma4:68 51:247 7 7 7 7 5
We take our algorithm with B being the identity matrix.
The following table reports the values of criterion after each sweep, that is after
steps associated with each of the 15 possible pairs of indexes.
Criterion 0:809676 0:189367 0:00562301
The last sweep produces a zero value of the criterion up to machine precision, the
slightly negative value we have got comes from the rounding errors. Note that since
there are only 2 matrices, exact joint diagonalization can be achieved. Actually, after
3 sweeps (sweep 0 corresponds to the initial matrices) the diagonalization is already
quite good. We have
50:0000 \Gamma0:0198 \Gamma0:0013 \Gamma0:0001 0:0000 0:0000
\Gamma0:0013 0:0001 29:8099 0:0000 \Gamma0:0000 \Gamma0:0000
\Gamma0:0001 0:0000 0:0000 39:0333 \Gamma0:0000 \Gamma0:0000
0:0000 \Gamma0:0000 \Gamma0:0000 \Gamma0:0000 20:1550 \Gamma0:0000
0:0000 \Gamma0:0000 \Gamma0:0000 \Gamma0:0000 \Gamma0:0000 10:04497 7 7 7 7 5
20:0000 0:0177 0:0052 0:0000 \Gamma0:0000 \Gamma0:0000
0:0177 10:0000 \Gamma0:0003 \Gamma0:0000 0:0000 0:0000
0:0052 \Gamma0:0003 40:0120 \Gamma0:0000 0:0000 0:0000
0:0000 \Gamma0:0000 \Gamma0:0000 30:7912 0:0000 0:0000
\Gamma0:0000 0:0000 0:0000 0:0000 59:1714 0:0000
\Gamma0:0000 0:0000 0:0000 0:0000 0:0000 48:47167 7 7 7 7 5
which corresponds to the transformation matrix
(For definiteness, the rows of B have been normalized to have unit norm.) The fourth
sweep zeros all off diagonal elements of C 1 and C 2 (at least up to 4 digits after the
decimal point) without changing their diagonal elements. The transformation matrix
B is also almost unchanged:
0:5000 0:5000 \Gamma0:5000 \Gamma0:5000 0:0000 \Gamma0:0000
0:5000 0:5000 0:5000 0:5000 \Gamma0:0000 0:0000
0:5527 \Gamma0:5527 0:4206 \Gamma0:4206 0:1688 \Gamma0:0817
0:3977 \Gamma0:3977 \Gamma0:5752 0:5752 \Gamma0:0664 \Gamma0:1324
0:0073 \Gamma0:0073 \Gamma0:0272 0:0272 0:6082 0:79287 7 7 7 7 5
One can see that our algorithm converges quite fast. The Flury and Gautschi
needs 4 to 5 sweeps to converge. Moreover it makes several iterations
for each pair of indexes while we make only one. However, our algorithm does
not solves the same problem, since we do not require the transformation matrix to
be orthogonal. A simple way to implement the orthogonality constraint, at least
approximately, is to add another matrix C 3 which is the identity matrix and give it
a large weigh n 3 . For 1), the values of the criterion after each
sweep are given below
Criterion 0:809676 0:226183 0:0291083 0:0290463 0:0290454 0:0290454
The criterion does not decrease further after 4 sweeps. The change in the transformation
produced by the fifth sweep is also very slight, affecting only the
last digit and never more than 2 units. This matrix, after sweep 5, is
and the corresponding matrices C 1 , C 2 are
50:0000 0:0000 0:0000 0:0000 \Gamma0:0000 \Gamma0:0000
0:0000 29:9224 0:0000 \Gamma1:8497 2:2318 0:1111
0:0000 0:0000 60:0000 0:0000 \Gamma0:0000 \Gamma0:0000
0:0000 \Gamma1:8497 0:0000 39:7221 \Gamma0:7727 1:0432
\Gamma0:0000 2:2318 \Gamma0:0000 \Gamma0:7727 20:2390 \Gamma0:0385
\Gamma0:0000 0:1111 \Gamma0:0000 1:0432 \Gamma0:0385 10:02407 7 7 7 7 5
20:0000 \Gamma0:0000 \Gamma0:0000 \Gamma0:0000 \Gamma0:0000 \Gamma0:0000
\Gamma0:0000 40:2097 \Gamma0:0000 \Gamma2:3088 4:4428 1:7605
\Gamma0:0000 \Gamma0:0000 10:0000 \Gamma0:0000 \Gamma0:0000 \Gamma0:0000
\Gamma0:0000 \Gamma2:3088 \Gamma0:0000 31:7746 \Gamma1:2795 7:2126
\Gamma0:0000 4:4428 \Gamma0:0000 \Gamma1:2795 59:3949 0:5032
\Gamma0:0000 1:7605 \Gamma0:0000 7:2126 0:5032 48:34577 7 7 7 7 5
These results are very similar to that of Flurry and Gautsch (1986). (Note that our
matrix B is the transposed of their.) Of course, the orthogonality constraint is not
exactly satisfied here. We have
1:0000 \Gamma0:0000 0:0000 \Gamma0:0000 0:0000 0:0000
\Gamma0:0000 1:0000 \Gamma0:0000 0:0119 \Gamma0:0185 \Gamma0:0047
0:0000 \Gamma0:0000 1:0000 \Gamma0:0000 0:0000 0:0000
\Gamma0:0000 0:0119 \Gamma0:0000 1:0000 0:0060 \Gamma0:0253
0:0000 \Gamma0:0185 0:0000 0:0060 1:0000 \Gamma0:0007
0:0000 \Gamma0:0047 0:0000 \Gamma0:0253 \Gamma0:0007 1:00007 7 7 7 7 5
but the difference of this matrix from the identity matrix is slight. We should mention
here that our algorithm is not designed to enforce orthogonality, the above numerical
results are given only as examples showing its good convergence property.


Appendix

Joint diagonalization of two Hermitian matrices of size two
The following results provide the explicit and complete solutions to the problem
of joint diagonalization of two non proportional Hermitian matrices of order two.
(Proportionality should here be understood in the large sense so that a null matrix
is proportional to any other one.) Note that if the matrices are proportional, then
the problem degenerates to the diagonalization of a single matrix.
Proposition
Let P and Q be two non proportional Hermitian matrices of order two,
with diagonal elements diagonal elements
respectively. Then
are not all zero and real and the matrix
a b
c d
diagonalizes
simultaneously P and Q if and only if (i) [a b] or [c d] vanishes or (ii)
[a and proportional to [2ff fi + ffi] or to
[c d] 6= [0 0] and proportional to [2ff
ffi] or to [fi
where ffi is any one of the two square roots of \Delta.
Remark:
proportional to [2ff fi implies that it is
proportional to conversely, since
then [a b] proportional to [2ff fi reduces to it being proportional to [0 1],
which is the same as it being proportional to [fi \Gamma ffi 2fl] only if ffi has been chosen
to be fi and not \Gammafi . Similar conclusion holds if Also in the condition (ii) of
the Proposition, one has to exclude the case [a since then [c d] needs not
be proportional to anything. Similarly, the case [c has to be excluded.
hence [c d] / [a b]. Thus P and Q are not
only be diagonalized but are actually transformed by a singular matrix into the null
matrix.
Proof
Since P and Q are non proportional, the vectors [p 1
are linearly independent. This entails that ff, fi, fl are not all zero.
To prove that real, we expand it as
Consider now the solutions to the joint diagonalization problem. The condition
that the transformed matrices be diagonal can be written as
[a b]P
[a b]Q
d
or
[c d]P
[c d]Q
- a
If we exclude the trivial solutions [a then the above equations
imply that the matrices in their left hand side have zero determinants. Thus
and the same equation holds with a, b replaced by c, d. After expansion, one gets
the equations ffb
The solution [a b] to the equation ffb determined
only up to a multiplicative factor. Thus, if ff 6= 0 then either [a
a 6= 0 and b=a is the root of the quadratic polynomial ffz
is a square root of \Delta. Therefore [a b] is proportional to
Similarly, if fl 6= 0, then [a b] must be proportional to We have
chosen here the minus sign before ffi so that this solution is the same as the previous
one (with the same ffi) in the case where ff 6= 0. Note that these solutions still apply
in the case since then they reduce to that [a b] being proportional to
[0 fi] or [fi 0] which is the same as ab = 0.
Similar calculations apply to the equation ffd solutions
[c d] must be proportional to [2ff is a square root of
(not necessarily the same as ffi).
We have shown that the rows [a b] and [c d] of a matrix jointly diagonalizing
P and Q must have the above form. But we haven't proved the converse. Further,
our results yield two choices for [a b] (modulo a constant multiple), depending on the
choice of the square root ffi of \Delta, and similarly there are two other choices for [c d].
Hence one must determine which choice of the later must be associated with one of
the former. For this purpose, we shall consider, for each choice [a
[a the possible corresponding choices for [c d] of the form [ff fi \Sigma - ffi]
or [fi \Upsilon - ffi 2fl] and see if there is one for which [a b]P[c d]
will be more convenient here to write the two square roots of \Delta in the form \Sigma - ffi.)
Suppose that ff 6= 0, then we need to consider the choice [a
since the other is proportional. For [c ffi], we has
[a b]P[c d]
d
The last expression can be expanded as
Note that p 1 -
hence the first term in the above expression
using the equality -
fi)=2 or 2-ff-p+ -
reduces to
This expression vanishes if the minus sign is used in \Sigma. Hence [a b]Q[c d]
for the choice [c
ffi]. A similar calculation, with q 1 , q 2 , q in place of
shows that the same choice leads to [a b]Q[c d]  One needs not
consider further choices. Indeed if there exists another choice [~c ~
d] not proportional
to [2ff fi such that [a b]P[~c ~
since one already
has [a b]P[2ff and the vectors [~c ~
d] and [2ff fi are linearly
independent, one must have [a By a similar argument, one also must
have [a these two equalities [a can be easily
seen to imply that P and Q are proportional, which contradicts our assumption.
The calculations are similar in the case fl 6= 0, interchanging a with b, c with
d, ff with fl, p 1 with
reversing the sign
of ffi. Finally, if the possible choice for [a b] is either proportional to
[0 1] or to [1 0] according to ffi equal plus or minus fi. But it can be easily seen that
this case can happen if and only if either
p-q is real, as P and Q are non proportional. Then, it can be checked that the only
possible corresponding choice for [c d] is, in the first case, proportional to [1 0] or
to [0 1] and in the second case, to [0 1] or to [1 0]. This completes the proof of the
Proposition.
Corollary
With the same notation and assumption of the Proposition and assume
further that at least one of the matrices P and Q has positive determinant. Then
and the matrix
a b
c d
jointly diagonalizes P and Q if only if it equals
denotes the sign function) pre-multiplied by a permutation and a diagonal
matrix.
Proof
We first show that \Delta ? 0 if det P ? 0 or det Q ? 0. We shall prove
the result only for the case det Q ? 0, the proof for the other case is similar. As
the last two terms in the above right hand side
can be written as
Thus
the first condition implies that -. Then the second
condition implies that Thus P is proportional to Q, contradicting our
assumption.
As its roots are real and we denote as usual
\Delta the positive one.
\Delta so that fi since
one obtains from Proposition A1 that
a b
c d
equals
pre-multiplied
by a diagonal matrix. On the other hand, for
\Delta so that
by a similar calculation, the Proposition shows that
a b
c d
equals a
diagonal matrix times
. This yields the result of the corollary.



--R

A blind source separation technique using second-order statistics

Matrix computations.
An algorithm for the simultaneous orthogonal transformation of several positive definite symmetric matrices to nearly orthogonal form.
Common principal components in k groups.
Applied linear Algebra.
Multivariate Observations.
--TR

--CTR
Ale Holobar , Milan Ojsterek , Damjan Zazula, Distributed Jacobi joint diagonalization on clusters of personal computers, International Journal of Parallel Programming, v.34 n.6, p.509-530, December 2006
Y. Moudden , J.-F. Cardoso , J.-L. Starck , J. Delabrouille, Blind component separation in wavelet space: application to CMB analysis, EURASIP Journal on Applied Signal Processing, v.2005 n.1, p.2437-2454, 1 January 2005
Ch. Servire , D. T. Pham, Permutation correction in the frequency domain in blind separation of speech mixtures, EURASIP Journal on Applied Signal Processing, v.2006 n.1, p.177-177, 01 January
Andreas Ziehe , Motoaki Kawanabe , Stefan Harmeling , Klaus-Robert Mller, Blind separation of post-nonlinear mixtures using linearizing transformations and temporal decorrelation, The Journal of Machine Learning Research, v.4 n.7-8, p.1319-1338, October 1 - November 15, 2004
Andreas Ziehe , Motoaki Kawanabe , Stefan Harmeling , Klaus-Robert Mller, Blind separation of post-nonlinear mixtures using linearizing transformations and temporal decorrelation, The Journal of Machine Learning Research, 4, 12/1/2003
Andreas Ziehe , Pavel Laskov , Guido Nolte , Klaus-Robert Mller, A Fast Algorithm for Joint Diagonalization with Non-orthogonal Transformations and its Application to Blind Source Separation, The Journal of Machine Learning Research, 5, p.777-800, 12/1/2004
principal components;separation of sources;diagonalization