Algorithm for finding the symmetries of a linear operator $\DeclareMathOperator\Hom{Hom}\DeclareMathOperator\GL{GL}$Let $V, W$ be finite dimensional complex vector spaces and $M\in \Hom(V, W)$ a full rank linear map. I want to see if there exists a Lie group $G$ and representations $\pi: G \to \GL(V)$ and $\rho: G \to \GL(W)$ s.t. $\rho(g) M = M \pi(g)$. The problem is that I don't know what the group $G$ is. I'm looking for an algorithm that will find the group $G$ and the representations $\rho, \pi$ for any $M$. In case it makes things easier, I mostly care about finding unitary symmetries of $M$, meaning that $V,W$ have inner products and the representations are of the form $\pi: G \to U(V)$ and $\rho: G \to U(W)$.
Two questions:

*

*Is there a well defined notion of the maximal symmetry group of $M$? Obviously, the trivial group will always be a symmetry of $M$, I care about finding all the symmetries.

*Does there exist an algorithm to solve this problem? I'm imagining that the group $G$ could be constructed by finding generators. Starting with Lie algebra generators $X\in \mathfrak{gl}(V), Y\in \mathfrak{gl}(W)$ s.t. $YM = MX$. The process would start with the trivial Lie algebra and would extend it step by step until no more extensions can be found. Then, an analogous process should find all the discrete symmetries. The result should be the maximal symmetry group of $M$ together with the representations $\pi$ and $\rho$.

 A: Note that the identity $\rho M = M \pi$ (I'm omitting the $g$ argument here) implies that $\rho$ has to preserve the image of $M$ and can do anything it likes the cokernel of $M$, while $\pi$ has to preserve the kernel of $M$ and it can do anything it likes on the complement. Writing $V = \ker M \oplus Z$ and $W = Z \oplus \operatorname{coker} M$, where the subspaces $\operatorname{im} M = Z = \operatorname{coim} M$ are identified by putting $M$ in reduced row echelon form, the obvious most general parametrization of $\rho$ and $\pi$ in an adapted basis are
$$
  \rho = \begin{bmatrix} \zeta & * \\ 0 & GL(\operatorname{coker}M) \end{bmatrix}, \quad
  \pi = \begin{bmatrix} GL(\ker M) & * \\ 0 & \zeta \end{bmatrix}, \quad \text{with} \quad
  M = \begin{bmatrix} 0 & I \\ 0 & 0 \end{bmatrix} ,
$$
where $\zeta \in GL(Z)$, $GL(-)$ schematically denotes any elelement of the corresponding group, and $*$ denotes an arbitrary matrix block of the right dimension.
If you want to restrict $\rho$ and $\pi$ to be unitary, then it's a little bit more complicated. Instead of using reduced row echelon form, you have to use the singular value decomposition (SVD) to bring $M$ into positive diagonal form $\Sigma$ on $Z$. In an adapted basis, $\rho$ and $\pi$ are parametrized as
$$
  \rho = \begin{bmatrix} \kappa & 0 \\ 0 & U(\operatorname{coker}M) \end{bmatrix}, \quad
  \pi = \begin{bmatrix} U(\ker M) & 0 \\ 0 & \lambda \end{bmatrix}, \quad \text{with} \quad
  M = \begin{bmatrix} 0 & \Sigma \\ 0 & 0 \end{bmatrix} ,
$$
where $\kappa,\lambda \in U(Z)$ and $\kappa \Sigma = \Sigma \lambda$. The standard solution of the last equality is $\kappa = \lambda$ with the matrices in block-diagonal form, each diagonal block being unitary, the blocks corresponding to the groups of equal singular values in $\Sigma$.
In each case, you can reconstruct the group $G$ and the multiplication rule in it from the above parametrizations.
