SVD of $2\times 2$ matrix where entries have different units I have the following matrix:
$$
M = \begin{pmatrix}
1 + xy & y \\
x & 1
\end{pmatrix}
$$
where $x$ has the unit $m^{-1}$ (per meter) and $y$ has the unit $m$ (meter). This matrix acts on real vectors with units $([m], [1])$.
Now I need to compute $M^n$ for $n \in \mathbb{N}$. My idea was to compute the SVD of the matrix $M$ in order to use the property $(U\Sigma U^\ast)^n = U\Sigma^nU^\ast$. As shown here the SVD of a 2x2 matrix can be obtained analytically by rewriting $M$ as a sum of the Pauli matrices $\sigma_i$ and the identity matrix $I$:
$$
M = z_0I + z_1\sigma_1 + z_2\sigma_2 + z_3\sigma_3
$$
This determines the coefficients $z_i$ as follows:
$$
\begin{aligned}
z_0 &= 1 + \frac{xy}{2} \\
z_1 &= \frac{x + y}{2} \\
z_2 &= i\frac{y - x}{2} \\
z_3 &= \frac{xy}{2}
\end{aligned}
$$
This however requires the two quantities $x$ and $y$ to be added together which doesn't work since they have incompatible units. So I'm wondering if that is the right approach, or if SVD is applicable at all in this case. Can SVD applied to matrices with (different) physical dimensions?
 A: This is a stark demonstration of the mess inflicted by failure to nondimensionalize.  If your quantities are dimensional, write their units next to them, in your case:
$$
M = \begin{pmatrix}
1 + xy & ym \\
x/m & 1
\end{pmatrix},
$$
which acts on vectors $( am, b)^T$, to produce like-dimensioned vectors, for numerical x,y,a,b.
Having ensured dimensional consistency of this map, you may now observe
$$
M = \begin{pmatrix}
1 + xy & ym \\
x/m & 1
\end{pmatrix}\\ =  \begin{pmatrix}
1  &0 \\
 0 & 1/m
\end{pmatrix} 
 \begin{pmatrix}
1 + xy & y \\
x & 1
\end{pmatrix}  \begin{pmatrix}
1  &0 \\
 0 &  m
\end{pmatrix} \equiv S^{-1} N S,
$$
where now, any and all dimensional information is carried by the similarity matrix S; while the equivalent matrix N is completely devoid of dimension, mapping numerical 2-vectors to such, and can be dealt with in the conventional manner. If you chose to use different units, like light-years, you only scale your Ss, and leave N and all functions thereof alone.
You may now take powers of N, strictly equivalent to the same powers of M, in the conventional manner -- but your $z_3$ needs a sign reversal! :
$$
N= z_0 1\!\! 1 + z ~~\hat {z} \cdot \vec {\sigma}, \\
z\equiv \sqrt{z_1^2+z_2^2+z_3^2}, ~~\hat{z}\cdot \hat{z}=1, ~~~\leadsto \\
N^2=(z_0^2+ z^2) + 2 z_0z~~\hat {z} \cdot \vec {\sigma},
$$
etc. I would assume you need no help with the recursion. But you could do so much better with Euler's formula for the exponential of the Pauli vector...

NB An easier way... basically SVD
You may further similarity transform N into a symmetric matrix, so orthogonally diagonalizable:
$$
N=  
 \begin{pmatrix}
1 + xy & y  \\
x  & 1
\end{pmatrix}\\ =  \begin{pmatrix}
1  &0 \\
 0 &  \sqrt{x/y}
\end{pmatrix} 
 \begin{pmatrix}
1 + xy & \sqrt{xy} \\
\sqrt{xy} & 1
\end{pmatrix}  \begin{pmatrix}
1  &0 \\
 0 &   \sqrt{y/x}
\end{pmatrix} \equiv T^{-1} K T,
$$
and K is symmetric and orthogonally diagonalizable via O, with eigenvalues λ, 1/λ, so that, at the end of the day,
$$
M^n= S^{-1}T^{-1} O^{-1}~~\operatorname{diag}(\lambda^n,1/\lambda^n)  ~~OTS. ~.
$$
