Block-diagonalization with unitary similarity transformations: ($A \rightarrow U B U^\dagger $, $B$ block-diagonal) The problem
Given a matrix $A$, find a unitary matrix $U$ such that $U^\dagger A U=B$, where $B$ is approximately block-diagonal (when possible).
Explanation
In other words, assume I'm given a matrix with some unknown approximate block-diagonal structure which has been scrambled by an unknown unitary similarity transformation. I want to be able to find and undo the similarity transformation to recover as much block-diagonal structure as I can.
Obviously I do not want the trivial solution of one block. I want to find two or more blocks, if they can be found. I also want to allow for the blocks to be imperfect: the found $B$ can have off-block elements, so long as there are relatively few with large magnitudes.
All matrices here can have complex entries, though we can consider more restricted cases if that helps.
My thinking
If $U$ was not constrained to be unitary, this would be easy: I could just use the Jordan normal form. Though in practice this would be pretty trivial because my $A$ matrix is noisy and I don't expect any of its eigenvalues to coincide, so the resulting $B$ matrix would just be diagonal.
If $U$ is further constrained to be a permutation matrix, then I've had good results treating (the elementwise magnitudes of) $A$ as a weighted graph, and using graph-community clustering algorithms (specifically Louvain partitioning) to extract the blocks. But this sort of method certainly doesn't seem suited  to cases where $U$ is allowed to be any unitary.
I'm happy to use heuristic / numerical / iterative methods, or (I guess) even to resort to some kind of gradient descent penalizing the L1 norm of $UAU^\dagger$, favoring sparsity.
 A: Let $u_1,\dots,u_n$ denote the columns of the unitary matrix $U$. $U^\dagger AU$ will have a block diagonal structure with $m$ block and each block "connecting" the $k_j$th and $(k_{j+1} - 1)$th diagonal entries if and only if the columns $u_{k_j},\dots,u_{k_{j+1}-1}$ span an invariant subspace of $A$ for each $j = 1,\dots,m$. Notably, these invariant subspaces would have to be mutually orthogonal because $U$ is unitary.
It turns out to be particularly easy to systematically consider the invariant subspaces of $A$ in the case that $A$ has distinct eigenvalues: each invariant subspace is simply a combination of one-dimensional eigenspaces.
With all that established: given a matrix $A$ with distinct eigenvalues that has an (exact) non-trivial block-diagonalization, the following procedure will necessarily block-diagonalize $A$ with the largest possible number of blocks.

*

*Find an eigenbasis $v_1,\dots,v_n$ of $A$. Let $V$ denote the matrix with $v_1,\dots, v_n$ as its columns.


*Compute the matrix $H = V^\dagger V$. For some (small) threshold $\epsilon > 0$, generate the adjacency matrix $J$ as follows. For each $1 \leq i,j \leq n$, take
$$
J_{ij} = \begin{cases}
1 & i \neq j \text{ and } |H_{i,j}| \geq \epsilon,\\
0 & \text{otherwise.}
\end{cases}
$$
The graph $G$ over $\{1,\dots,n\}$ associated with this adjacency matrix is such that for $i \neq j$, $i \sim j$ if and only if $v_i$ fails to be orthogonal to $v_j$ (up to some numerical threshold). A non-trivial block-diagonalization of $A$ exists if and only if the graph $G$ is disconnected.


*Decompose $G$ into connected components, let $m$ denote the number of connected components. If $m=1$, then no non-trivial block-diagonalization exists. Otherwise, let $S_j = \{k_{j,1}, \dots, k_{j,p_j}\}$ denote the connected components of $G$ for $j = 1,\dots,m$. Let $W_j$ denote the matrix whose columns are $v_k$ for all $k \in S_j$ (the order doesn't matter). Let $W$ denote the matrix $W = [W_1 \ \ \cdots \ \ W_m]$


*Compute a unitary $QR$ decomposition of $W$, which is to say we find a unitary $Q$ and upper-triangular $R$ such that $W = QR$.
Taking $U = Q$ gives us our block-diagonalization. That is, $Q^\dagger AQ$ will be block-diagonal with $m$ blocks.

Notes:

*

*The most computationally expensive of the above steps is the eigenvector computation, which has complexity $O(n^3)$. In other words, the above describes a polynomial-time algorithm.


*All of the blocks of $B$ attained in this fashion will be upper-triangular


*I am not sure how good this procedure will be at detecting all possible "imperfect" block diagonalizations. The $\epsilon$ threshold gives the algorithm the flexibility to find subspaces that are "exactly" invariant and "approximately" orthogonal, but it isn't able to use "approximately" invariant subspaces. This means that the entries within the same "block-row" as a given block in $B = U^\dagger AU$ and to the right or the same column and above are allowed to be non-zero up to a tunable threshold, but the elements within the same block-row and to the left or the same block-column and below are not. In other words, $B$ is "approximately" block-diagonal but "exactly" block-upper-triangular.


*I claim a lot without proof in the above. I'm happy to provide justification, if there's interest, when I have time.
