Decomposition of a nonsquare affine matrix

You've written this somewhat unorthodoxly. To use that matrix for that transformation, one would more usually write

$$\pmatrix{x'\\y'\\1}=\pmatrix{a&b&c\\d&e&f\\0&0&1}\pmatrix{x\\y\\1}\;.$$

So the difference between a $2\times3$ matrix and a $4\times4$ matrix was only from your way of writing it; this works the same way as an affine transform in three dimensions, just with one fewer dimension. You can immediately factor out the translation,

$$\pmatrix{x'\\y'\\1}=\pmatrix{1&0&c\\0&1&f\\0&0&1}\pmatrix{a&b&0\\d&e&0\\0&0&1}\pmatrix{x\\y\\1}\;.$$

Then you just have to decompose $\pmatrix{a&b\\d&e}$ into shear, rotation and scaling in two dimensions.

[Edit in response to the comment:]

This isn't a unique decomposition, since you can do the shear, rotation and scaling in any order. Here's the decomposition I use:

$$A=\pmatrix{a&b\\d&e}=\pmatrix{p\\&r}\pmatrix{1\\q&1}\pmatrix{\cos\phi&\sin\phi\\-\sin\phi&\cos\phi}$$

with

$$ \begin{eqnarray} p&=&\sqrt{a^2+b^2}\;,\\ r&=&\frac{\det A}p=\frac{ae-bd}{\sqrt{a^2+b^2}}\;,\\ q&=&\frac{ad+be}{\det A}=\frac{ad+be}{ae-bd}\;,\\ \phi&=&\operatorname{atan}(b,a)\;, \end{eqnarray} $$

where $\operatorname{atan}$ is the two-argument arctangent function with operand order as in Java. This of course assumes $p\ne0$.

If $(x, y, 1)$ is a vector in homogeneous coordinates, we have, by decomposing $M$ into blocks, that

$$M \left[\begin{array}{c}x\\y\\1\end{array}\right] = \left[\begin{array}{cc} a& b\\ d&e\end{array}\right]\left[\begin{array}{c}x\\y\end{array}\right] + \left[\begin{array}{c}c\\f\end{array}\right].$$

Here $(c,f)$ is the translation component. We can decompose the 2x2 matrix into a composition of a rotation, shear, and scale by using the QR decomposition:

$$\begin{align*}\left[\begin{array}{cc}a & b\\d & e\end{array}\right] &= \left[\begin{array}{cc} \cos \theta &-\sin \theta \\ \sin\theta &\cos \theta\end{array}\right]\left[\begin{array}{cc} \sqrt{a^2+d^2} & b\cos \theta + e\sin \theta\\0 & e\cos \theta - b\sin \theta\end{array}\right]\\ &=\left[\begin{array}{cc} \cos \theta &-\sin \theta \\ \sin\theta &\cos \theta\end{array}\right]\left[\begin{array}{cc}1 & \frac{b\cos \theta + e\sin\theta}{e\cos \theta-b\sin\theta}\\0 & 1\end{array}\right]\left[\begin{array}{cc}\sqrt{a^2+d^2} & 0\\0 & e\cos\theta - b\sin\theta\end{array}\right],\end{align*}$$ where $\theta = \arctan\left(\frac{d}{a}\right).$