Efficient algorithm to find a symmetric matrix that commute with a given set of matrices

Question

Given a set of $K$ real $n \times n$ matrices $\{B_k\}$, what's the most efficient way to find a symmetric matrix $X$ that commute with all of $\{B_k\}$, i.e.,

$[X,B_k] = XB_k - B_kX = 0 \quad \forall \; k=1,\dots, K.$

The matrix $X$ needs to be a generic solution in the sense that it has no unnecessary eigenvalue degeneracy (e.g., identity matrix $I_n$ is always a solution but not a generic one).

I can map this problem to an eigenvector problem of $n^2 \times n^2$ matrix and solve it with method of $\mathcal{O}(n^4)$ complexity. But I wonder whether more efficient algorithms exist. Approximate solutions with $[X,B_k] < \epsilon$ are acceptable and any practical methods with high success rate are welcome.

Answer 1

I'm not sure how well this will work, but have you tried the following? (I'll be talking about the simultaneous block-diagonalization problem immediately). Note that if the matrices $B_j$ are simultaneously block-diagonalizable, then, for any linear combination $B=\sum_j \xi_j B_j$, the top eigenvector of $B$ is lying in one of the blocks. So, start with taking a random linear combination $B$ and finding the top eigenvector $v$ by some fast method (see something like this though I surmise you know such stuff better than I do). That is (almost and on average) $O(Kn^2)$ (at least it is as much $O(Kn^2)$ and your $O(n^4)$ is really $O(n^4)$. Now notice that that block should be (nearly, if you are fine with an approximate diagonalization) closed under the action of all $B^j$ and $B_j^T$. So, start running Gram-Schmidt on $v_1=v$, $B_j v, B_j^Tv$ ignoring the vectors that test almost linearly dependent with the previous ones. The orthogonalization itself is fast until the block grows really large and even then it is just $n^2$ per vector, so I'll not even count it. The real quadratic cost comes from finding the images. Now assume that you obtained $s$ vectors $v_1,v_2,\dots,v_s$ in those $2Kn^2$ steps. Mix them into a random combination (Gaussian randomness is the best, as usual) and try to see if every image of that combination is (nearly) in the span of $v_1,\dots,v_s$ under $B_j$ and $B_j^T$ (sam G-S). If it is, declare a complete block (you may want to test 2-5 combinations to be more sure). If not, expand the system accordingly and repeat. Every time your cost is $2Kn^2$ and you increase the dimension of the block by at least $1$.

Suppose you have a chunk of your basis now and you know that it is (nearly) invariant. Replace your operators by $(I-P)B_j(I-P)$ where $P$ is the projection to the chunk (only do not multiply matrices but operate on vectors instead when iterating for eigenvalues) and repeat. You'll get a second block and so on. The total cost of constructing the basis is $O(Kn^3)$. Since the cost of the single final rotation of the operators to the basis is $Kn^3$ (unless you use really fancy matrix multiplication schemes, but they become effective only in the range of thousands), so you can hardly do better. The final check that the rotation worked is about $Kn^2$, which is nothing again.

Does that make sense to you? If it does, try it and let me know what happens.

Answer 2

Another perspective:

What you're really doing is solving a system of Sylvester equations. In particular, you're looking for the nullspace of the $(n^2 K) \times n^2$ matrix given by $$ \pmatrix{I \otimes B_1 - B_1^T \otimes I \\ I \otimes B_2 - B_2^T \otimes I\\ \vdots\\ I \otimes B_K - B_K^T \otimes I} $$ With the further constraint that $X^T = X$. You might be able to reduce that down a bit using half-vectorization in some way. There are plenty of linear systems algorithms around. That being said, I'm not convinced any of them beat $O(n^4)$ in this case.

Answer 3

$X$ must map the generalized eigenspaces of each $B_k$ into themselves. That is, if $(B_k - \lambda)^m v = 0$ then $(B_k - \lambda)^m Xv = 0$. Depending on the structure of these eigenspaces, that can be a very helpful constraint.