# commuting matrices

I have a set $\{A_1, A_2, .. A_k\}$ of $n$ by $n$ real matrices and I know that they are 'perturbated versions' of a set of commuting matrices : $\{P_1,..,P_k\}$, by perturbated versions I mean that I have a bound on $\sum_i\|A_i-P_i\|\leq \epsilon$. My problem is that I want to find an epsilon perturbation of $\{A_1,...,A_k\}$ so that they commute. It is clear that at least one such perturbation exists since they are perturbed versions of a commuting set. But my question is if there is any algorithm to find a permutation that makes the matrices in $A$ commute.

-
Well, $A_1,\ldots,A_n$ span a commutative sub-algebra. Just take $P_1,\ldots,P_n$ in the same subspace. –  Denis Serre Jan 6 '11 at 20:10
@Denis If $A_1$ and $A_2$ are 2 by 2 with constant diagonals, one upper triangular and the other lower, then they span all 2 by 2 matrices. They might be perturbed diagonal matrices. –  Aaron Meyerowitz Jan 6 '11 at 20:35
The problem isn't well stated. Are $A_1,\dots,A_k$ commuting? What happens if one simply takes $P_j=(1+\delta)A_j$ for a suitable infinitesimal $\delta$? –  Wadim Zudilin Jan 6 '11 at 23:01
"permutation"?? –  darij grinberg Jan 7 '11 at 12:30

I think it really comes down to looking at the lattice of invariant subspaces for each $A_i$. I believe this is what Igor basically is saying.

For a single real linear transformation $V \rightarrow V$ with distinct characteristic roots, $V$ decomposes as a sum of 1 and 2-dimensional invariant subspaces, and these generate the lattice of invariant subspaces. When you perturb the matrix by a known amount, the subspaces can wander around over a certain range --- how much they shift depends on how close together are the characteristic roots, and also on the angles between the subspaces. For simplicity, think of the case with all real distinct eigenvalues. The projectivized action, on $\mathbb{RP}^{n-1}$, then has $n$ fixed points. The first derivative of the projectivized map is simple to determine, given the eigenvectors and eigenvalues --- the eigenvalues of the derivative at a point for eigenvalue $\lambda$ are the ratios of the other eigenvalues to $\lambda$. Therefore, if you change the map by adding multiples of the other eigenvectors as a linear function of the projection to the eigenspace under consideration you shift your given eigenvector by a known linear function of the coefficients.

If you have several matrices with all real distinct eigenvalues that are near enough to a commuting set of matrices, then presumably the eigenvectors are also near to each other, so you can tell which eigenvectors need to line up with each other. You have a good measurement of how hard it is to move each eigenvector in a given direction, so you can find a point that's approximately equally hard for each (do this separately, for each matching collection of eigenvectors.) Then adjust each of the eigenvectors to go to that point, and repeat on the other eigenspaces.

When there are complex roots, there are 2-dimensional invariant subspaces, which appear as invariant circles in $\RP^{n-1}$. You can do much the same thing with these.

When there are multiple real or complex roots, then it's necessary to look at larger-dimensional invariant subspaces. One strategy is to first fix up invariant subspaces associated with the product of all factors of the characteristic polynomial whose roots are in a small neighborhood, or a pair of complex conjugate neighborhoods. The goal is to get such subspaces to agree with, to contain, or to be contained in similar suspaces for the other matrices. If there are repeated roots, then first see if there's a small perturbation that makes multiple eigenvectors (or multiple 2-dimensional invariant subspaces) --- when this is possible, it gives a bigger target for matching commuting matrices. If not, look at the largest invariant subspace it contains.

I doubt if there's an easy way other than looking at these invariant subspace structures. Perhaps in the actual applications you have in mind, they're not so complicated.

-
All that merely restates the question. It is totally non-obvious to me how to efficiently find the closest in the sum of operator norms metric commuting family even if we restrict ourselves to orthogonal projections to 1D subspaces for both the original operators and their perturbations. –  fedja Jan 7 '11 at 5:49
Dear fedja, I don't see why you say that this simply restates the question. It suggests an algorithm for moving efficiently as possible to a commuting family: i.e., move nearby eigenlines to a common eigenline in way "that's approximately equally hard for each" (i.e. tries to uniformly bound the perturbation of each $A_i$ as much as possible), where the calculation implicit in "approximately" is made using the discussion of derivatives in the previous paragraph. This seems to be more information about the problem than is in the original question. –  Emerton Jan 7 '11 at 6:50
Dear Emerton, my point merely was that we all know that the commuting families are (essentially) families with common eigenspaces. The devil is in the details. The phrase "move in a way that is approximately equally hard for each" is neither very instructive, nor even always a correct way to go. Algorithm stability is another important issue. You are right that something along these lines is worth trying and may work but even in the example I mentioned (where all eigenspaces are totally explicit), I cannot say that I'm ready to sit and write a code after reading this post. –  fedja Jan 7 '11 at 14:57
@fedja: In the initial part of my discussion, I was assuming that the perturbation is small enough compared to the differences of eigenvalues that the eigenspaces are clustered, so the pairing between them is obvious. The L^2 penalty for moving each a certain distance, in local projective coordinates, is given by a readily computed quadratic function (you have the derivative explicitly, you need to invert to calculate the penalty). If you want the smallest perturbation measured as the L^2 norm of the sequence of perturbations, this is minimum of the sum. Minimizing the max is the median. –  Bill Thurston Jan 7 '11 at 15:19
@fedja: the complexity comes from cases when there are multiple eigenvalues, or nearly equal eigenvalues. In this case, the pairing might not be obvious, but the complexity can be reduced by first finding subspaces of any dimension where the projective action is close to the identity --- by first clumping eigenvalues that are close together, then looking at the sizes of the shear effects (caused either by Jordan blocks, or eigenspaces that are projectively close to each other). It's a bit complicated to describe, but for moderate dimension, I think it should be computationally effeective. –  Bill Thurston Jan 7 '11 at 15:29

Well, the following is a not-too-fast algorithm. Reduce all matrices to Jordan canonical form, since the matrices almost commute, their right eigenvectors are almost the same. Adjust the matrices of eigenvectors to actually BE the same (and in case of nontrivial jordan blacks perturb them so the matrix is actually diagonal. Jordan blocks and sets of nearby eigenvalues makes implementing this a nuisance, since the eigenvectors of a slightly perturbed diagonal matrix are very unstable -- you will need to look at which subspaces those blocks span.

-

First, let us consider a simpler setting. Let $A_1$ and $A_2$ be given as inputs. The aim is to find a matrix $B$ such that $A_1+B$ commutes with $A_2+B$, and that $B$ has small norm.

The key-constraint is thus: $(A_1+B)(A_2+B) = (A_2+B)(A_1+B)$, which simplifies to the worst-case (since no unique solution is possible) Sylvester equation

$$(A_1-A_2)B + B(A_2-A_1) = A_2A_1-A_1A_2$$

Thus, a formulation as an optimization problem could be

$$\min\quad \|B\|^2\qquad \text{s.t.}\quad (A_1-A_2)B + B(A_2-A_1) = A_2A_1-A_1A_2.$$

In principle (if we use the Frobenius norm for $B$), this is just a least-norm problem (similar to least-squares).

In general, if we have $A_1,\ldots,A_k$, then we will have a huge-number of constraints (for all pairs of commutativity requirements), but the overall problem should still be essentially a least-norm problem, and can be thus solved easily.''

-

I think there is enough here to qualify as an answer. Plus it is too long for a comment.

1. Hey, it's just calculus! The input and output are lists of $kn^2$ values and we want to minimize the square of the distance subject to a list of constraints (that is my norm). Use Lagrange multipliers. Alternately one could try some method as above which depends on eyeballing the eigenvectors and then take the result as an initial point for some iterative solver. (OK that is probably totally impractical. I was just thinking about how close a given answer is to the absolute global minimum.)

2. If the perturbation is random (say noise of mean $\epsilon$ in each variable) then it might be "easy" to find a solution better than the starting data.

3. I used this combined method for a (super small, baby) problem to try the eyeball method and then see how it did relative to fedja's request for the absolute minimum (choosing a norm I could work with). I worked in MAPLE starting with the 2 by 2 symmetric matrix.

$$P_1=\left[ \begin {array}{cc} 20000&10000\\10000&10000 \end {array} \right]$$

and $P_2$ its square (actually, $1/10000$ times its square). Then I added a random integer in the range $-100..100$ to each entry preserving symmetry (I told you it was a baby problem). Then I got numerical eigenvalues and eigenvectors for each matrix, averaged the eigenvectors to get the compromise eigenbasis, and transformed back. The six random pertubation integers came out to be $-11,-83,70$ for the first matrix and $15,82,8$ for the second. The result of this process was a pair of commuting matrices which differed from $A_1$ and $A_2$ by approximately $5.35,-2.67,-5.35$ and $-16.23,8.1,16.23$. Then I asked MAPLE to do a MINIMZE problem in $\mathbb{R}^6$ finding the point on the surface $(r_1-r_3)r_5=r_2(r_4-r_6)$ closest to a given one (from the $A_i$). It failed to find an initial point but when I used my previous answer as an initial point it returned commuting matrices which differed from $A_1$ and $A_2$ by approximately $9.65,-4.81,-9.65$ and $-3.18,1.59,3.18$

Disclaimer: I realize that the constraints would be a nightmare for a pair (let alone $k$) large matrices. I don't have any idea of the state of the art for numerical solvers.

-