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 anderror-controlling parameter ε ≥ 0.

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

See [1] for more introduction and precise definition.

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

- Sample randomly from { X: n x n symmetric (or Hermitian) matrix | |X A_j - A_j X| ≤ ε (for all j) }.
- 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".

We have implemented the algorithm in Matlab.

- commdec.m (Last update: 2012.11.13)
- older versions

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.

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}) );

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}) );

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) );

- [1] T. Maehara and K. Murota (2011): Algorithm for error-controlled simultaneous block-diagonalization of matrices, SIAM Journal on Matrix Analysis and Applications, 32, 605-620.
- [2] K. Murota, Y. Kanno, M. Kojima and S. Kojima (2010): A numerical algorithm for block-diagonal decomposition of matrix *-algebras with application to semidefinite programming, Japan Journal of Industrial and Applied Mathematics, 27, 125-160.
- [3] T. Maehara and K. Murota (2010): A numerical algorithm for block-diagonal decomposition of matrix *-algebras with general irreducible components, Japan Journal of Industrial and Applied Mathematics, 27, 263-293.

- 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.