Error-Controlled Simultaneous Block-Diagonalization Algorithm


Problem

The Simultaneous Block-Diagonalization Problem is the following problem:

Given: n x n real (or complex) matrices A_1, ..., A_N.

Find: An orthogonal (or a unitary) matrix P such that P' A_1 P, ..., P' A_N P become block-diagonal matrices with a common diagonal structure (simultaneous block-diagonal form).

The problem is widely studied in many areas such as physics, numerical computation, optimization, signal processing.

In many applications, the input matrices are contaminated with numerical errors (observation errors, approximation errors, etc.) So the error-controlled variant of the problem, Error-Controlled Simultaneous Block-Diagonalization Problem, is more useful in practice:

Given: n x n real (or complex) matrices A_1, ..., A_N and error-controlling parameter ε ≥ 0.

Find: An orthogonal (or a unitary) matrix P such that P' A_1 P, ..., P' A_N P are ε-approximated simultaneous block-diagonal form, i.e., P' A_j P = (Block-Diagonal Matrix) + ( O(ε) Matrix ).

See [1] for more introduction and precise definition.

Algorithm

In [1], we have proposed the following simple algorithm, "MM-algorithm", for Error-Controlled Simultaneous Block-Diagonalization Problem.

  1. Sample randomly from { X: n x n symmetric (or Hermitian) matrix   |  |X A_j - A_j X| ≤ ε (for all j) }.
  2. Return an orthogonal (or a unitary) matrix P which diagonalizes X.

Step 1 is reduced to eigenvalue computation of a large matrix S (n^2 x n^2), and the Lanczos method can be applied.

Remark: The algorithm is a dual (in the sense of the commutant of the matrix *-algebra) version of the algorithm in [2,3], the so-called "MKKKM algorithm".

Source Code (in Matlab)

We have implemented the algorithm in Matlab.

The program is provided without warranty of any kind. You can use/modify/redistribute the program for any purpose.

We welcome your feedback; Please send e-mail to: maehara@nii.ac.jp

Remark: Since this program is implemented for a demonstration of the algorithm, you may have to adjust the program for your use. I believe that the code is easy to read.

Example of usage

Here we show the usage of the program. Put "commde.m" to your working directory.

First, build block-diagonal matrices B{:} which consist of 4 blocks of size 1x1, 2x2, 3x3, and 4x4.

>>for j = 1 : 5
B{j} = blkdiag(randn(1),randn(2),randn(3),randn(4)) + 0.01*randn(10);
end

>>colormap gray;
>>imagesc( -abs(B{1}) );
density of B{1}

The image shows the density of matrix B{1}. (larger value corresponds to blacker cell.)

Next, build input matrices A{:}

>>Q = orth(rand(10)); 
>>for j = 1 : 5
A{j} = Q' * B{j} * Q;
end;

>>imagesc( -abs(A{1}) );
density of A{1}

The problem is to recover the block-diagonal structure of B{:} from A{:}.

Use the function "commdec" in "commdec.m". Then we obtain the matrix P which recovers the block-diagonal structure.

>>P = commdec(A);
err = 1.907e-02
>>imagesc( -abs(P' * A{1} * P) );
density plot of P' A{j} P

References


Version History

2012.11.13
Increase the number of iterations in the Lanczos method
Sparse matrix implementation
Automatic parameter selection
2012.10.20
First published

Takanori MAEHARA (maehara@nii.ac.jp).

Last Modified: 2012.11.13. / Published: 2012.10.20.