As a refresher, the transformation matrix works like so:
|a d 0|
|x y 1| * |b e 0| = |x' y' 1|
|c f 1|
When two matrices are interpolated, each one is broken down like so:
M = S * R * T
The S is for scale and shear combined, R is for rotation, and T is transformation. There's a routine in the core code called MatrixDeconstruct() that handles this, and it outputs S*T (because the T is downright trivial to pull out) and an angle.
Step 1: Pull out T. This is super easy, because all you have to do is keep track of what the c and f values are. A translation matrix doesn't set anything else.
Step 2: Analyze the x and y axes. The x axis is (a,d), because those are the two values x gets multiplied by; the y axis is (b,e). In a non-rotated matrix, the x axis should be a multiple of (1,0) and the y axis is a multiple of (0,1). Find the cosine and sine of the angle by which the axis has been rotated from normal. The angle with the greater cosine (smallest angle) is the rotation angle.
For the x axis, the cosine is a/sqrt(a*a+d*d), and the sine is -d/sqrt(a*a+d*d). For y, it's e/sqrt(b*b+e*e) and b/sqrt(b*b+e*e). Note: If the square root term for either axis is near 0, the matrix is treated as degenerate; the rotation angle is then 0, and S*T is the same as M, so all the following steps are skipped.
Step 3: Multiply by the opposite rotation, R', to get S.
M = S * R * T
M/T = S * R
S = S * (R * R') = M/T * R'
Step 4: Discard any translation components (c and f) from S, and replace them with the original values (T), to get S * T. These are kept together in the calculation just because you can store them in a single matrix, and they're easy to separate again.
Once you have the rotation angle, S, and T for each matrix, you can determine where to go from there.
1) If either matrix is degnerate, linear interpolation is the only way forward.
2) If the difference in rotation angles is 0, linear interpolation is easier.
3) If the difference in rotation angles is 180°, linear interpolation is preferred. (A flip was probably the intended result, and anyway it's not clear which direction was intended for rotation.)
4) All other cases get done the long way.
For the long form, the final interpolated matrix Mi should be:
Mi = Si * Ri * Ti
Si and Ti can be interpolated linearly. Ri is simply the rotation matrix for angle1 + (angle_diff * t).