# 4X4 Matrix Decomposition

I was in need to decompose transformations given as 4X4 matrices into simpler ones. What I found on the web wasn’t suiting my need for understanding what I’m doing. So I had to figure it out myself.

## The Math of Matrix Decomposition

The lots of math involved and my self-set goal to keep this site as free from Javascript as possible do not allow to use common Javascript libraries to typeset math formulars here. Therefore I exhumed my LaTeX knowledge and wrote a paper what describes the math.

The paper describes how to decompose 4X4 matrices, and on the way the handled matrices
are becoming smaller, so when you are interested in **4X3** or **3X3 matrix decomposition**
you’ll find it in the paper, too. Just skip the first step(s).

A slightly abbreviated section for handling **2D** cases like **3X2** and, on the way,
**2X2 matrix decomposition** is also included.

## Source Code

Source code on the other hand can be shown here perfectly well.

Obviously the sources depend a lot on the internal representation of matrices, and there are many of them. As I’m lazy I just copied over code from my library. Although that is not open and you cannot use it directly if you are not a customer the naming and comments should definitely help to understand what is necessary to port it to your needs.

### Decomposition of a 4X4 Matrix

This is a default method in my `Transformation3D`

interface.

`Vector3D`

is a 3D column vector. `OrderedPair`

is a typed pair of values.

```
/**
* Decompose this transformation into simpler ones.
*
* Be prepared that there is not only one way to decompose a matrix.
*
* The algorithm will produced the following order (where some or all
* of the simple transformations might be missing if they are just
* the identity):
* <ol>
* <li>{@link TransformFlags#Projection}</li>
* <li>{@link TransformFlags#Translation}</li>
* <li>{@link TransformFlags#Rotation}</li>
* <li>{@link TransformFlags#Shearing}</li>
* <li>{@link TransformFlags#Scaling} (including possible mirroring)</li>
* </ol>
*
* @return a chain of tagged simple transformations like projection, translation, rotation, shearing
* and scaling that when multiplied in the given order recreate this
* transformation (within floating point accuracy).
* The order to apply this transformations to one object step by step is inverse!
* If this represents identity the
* returned list will be empty, indicating that no transformation needs to be applied.
* If this transformation has an 0 determinant it is not
* considered decomposable and {@code null} is returned.
*/
@Nullable
default List<OrderedPair<TransformFlags, Transformation3D>> decompose()
{
// Column vectors of intermediate 4X3 matrix
final Vector3D columnX = new Vector3D(getXX(), getYX(), getZX());
final Vector3D columnY = new Vector3D(getXY(), getYY(), getZY());
final Vector3D columnZ = new Vector3D(getXZ(), getYZ(), getZZ());
final Vector3D columnW = new Vector3D(getXW(), getYW(), getZW());
// 3X3 matrix used for various purposes
final Matrix3X3 m = Matrix3X3.fromColumns(columnX, columnY, columnZ);
final List<OrderedPair<TransformFlags, Transformation3D>> result = new LinkedList<>();
final double wx = getWX();
final double wy = getWY();
final double wz = getWZ();
final double ww = getWW();
if (wx != 0 || wy != 0 || wz != 0 || ww != 1) {
// The trick here is to assume that the transformation is a product
// of a projective matrix where only the last row differs from an identity
// matrix and a 4X3 matrix (in that order).
// Then the upper 3 rows come out unharmed, and the last row can be considered
// the transformation of a vector w. We only use a 3X3 sub matrix here for faster
// inversion, as the ww coordinate is easily calculated directly when we
// know w.
final Matrix3X3 subTransposedInverted = m.getTransposed().getInverse();
if (subTransposedInverted == null) {
Debug.warn("Decomposition of transformation with projective parts and uninvertable sub matrix: %0",
this);
return null;
}
final Vector3D w = subTransposedInverted.multVector(wx, wy, wz);
result.add(OrderedPair.create(TransformFlags.Projection,
new Projection3D(w, ww - w.mult(getXW(), getYW(), getZW()))));
}
// Translation is easy now when it is considered to be applied last (before a possible projection), which we do here.
if (!TrafoUtil.equals(columnW, IVector3D.ZERO)) {
result.add(OrderedPair.create(TransformFlags.Translation,
createTranslation(columnW)));
}
// Some early checks
if (m.isIdentity()) {
return result;
}
final double det = m.getDeterminant();
if (det == 0) {
Debug.warn("Decomposition of uninvertable subtransformation: %0", m);
return null;
}
if (det == 1 && isIdentity()) {
return result;
}
// Now what is left is a 3X3 matrix which is considered to be a combination
// of rotation, shearing and scaling.
// We make use of the fact that for a rotation R the following is always valid:
// R^T = R^{-1}
// We now consider the remaining 3X3 matrix M to be the combination of
// a rotation R and a shear/scale matrix S:
// M = R * S
// If we calculate
// M^T * M = (R * S)^T * (R * S)
// = (S^T * R^T) * R * S
// = S^T * R^{-1} * R * S // because R^T = R^{-1}
// it follows that
// M^T * M = S^T * S
// This means that we can derive S from the last equation. For this we
// assume that S is a triangle matrix in the form
// / a d e \
// S = | 0 b f |
// \ 0 0 c /
// Multiplying S^T * S results in
// / a² ad ae \
// S^T * S = | ad b²+d² bf+ed | = M^T * M
// \ ae bf+ed c²+e²+f² /
// Resolving this leads to the formulas for a,b,c,... below.
// The squares allow to take plus or minus values for the
// scaling coefficients on the diagonal. For symmetry we
// use either always + or always -, although switching just one
// of them would suffice. This selects 1 of the 8 possible solutions
// and takes care that any mirroring operation is put into the
// scaling, so all other transformations in the chain have det=1.
final Matrix3X3 mTm = m.getTransposed().mult(m);
final boolean minus = det < 0;
final double a = minus // a cannot be 0 because of det != 0 (and m^T already was invertable when we had projection)
? -Math.sqrt(mTm.getXX())
: Math.sqrt(mTm.getXX());
final double d = mTm.getXY() / a;
final double e = mTm.getXZ() / a;
final double b = minus // cannot be 0, see a
? -Math.sqrt(mTm.getYY() - d*d)
: Math.sqrt(mTm.getYY() - d*d);
final double f = (mTm.getYZ() - d * e) / b;
final double c = minus // cannot be 0, see a
? -Math.sqrt(mTm.getZZ() - e*e - f*f)
: Math.sqrt(mTm.getZZ() - e*e - f*f);
final Matrix3X3 shearAndScale = new Matrix3X3(a, d, e,
0, b, f,
0, 0, c);
// Having shear-and-scale we can calculate the rotation from
// R * S = M <=> R = M * S^{-1}
final Transformation3D sasInverse = shearAndScale.getInverse();
assert sasInverse != null; // as shearAndScale's det != 0
final Transformation3D rotation = m.mult(sasInverse);
if (!rotation.isIdentity()) {
result.add(OrderedPair.create(TransformFlags.Rotation,
new Matrix3X3(rotation)));
}
// Last step: decompose the scaling-and-shearing S.
// Here we assume scaling happening first, which makes S:
// S = shear * scale.
// The scale matrix is obvious and uses a,b,c. For the shearing
// we assume
// / 1 i j \
// shear = | 0 1 k |
// \ 0 0 1 /
// Multiplying gives
// / a bi cj \ / a d e \
// shear * scale = | 0 b ck | = S = | 0 b f |
// \ 0 0 c / \ 0 0 c /
// therefore:
final double i = d / b;
final double j = e / c;
final double k = f / c;
if (i != 0 || j != 0 || k != 0) {
result.add(OrderedPair.create(TransformFlags.Shearing,
new Matrix3X3(1, i, j,
0, 1, k,
0, 0, 1)));
}
if (a != 1 || b != 1 || c != 1) {
result.add(OrderedPair.create(TransformFlags.Scaling,
Transformation3D.createScaling(a, b, c)));
}
return result;
}
```

### Decomposition of a 3X2 Matrix

This is a default method in my `Transformation3D`

interface.

`Vector2D`

is a 2D column vector. `OrderedPair`

is a typed pair of values.

```
/**
* Decompose this transformation into simpler ones.
*
* Be prepared that there is not only one way to decompose a matrix.
*
* The algorithm will produced the following order (where some or all
* of the simple transformations might be missing if they are just
* the identity):
* <ol>
* <li>{@link TransformFlags#Translation}</li>
* <li>{@link TransformFlags#Rotation}</li>
* <li>{@link TransformFlags#Shearing}</li>
* <li>{@link TransformFlags#Scaling} (including possible mirroring)</li>
* </ol>
* @return a chain of simple transformations like translation, rotation, shearing and scaling
* (but no projection) that when multiplied in the given order recreate this
* transformation. If this represents identity the returned list will be empty.
* If the matrix has an 0 determinant it is not considered decomposable
* and {@code null} is returned.
*/
@Nullable
default List<OrderedPair<TransformFlags, Transformation2D>> decompose()
{
final double det = getDeterminant();
if (det == 0) {
Debug.warn("Decomposition of uninvertable transformation: %0", this);
return null;
}
if (det == 1 && isIdentity()) {
return Collections.emptyList();
}
final Vector2D columnX = getColumnX();
final Vector2D columnY = getColumnY();
final Vector2D columnW = getColumnW();
final List<OrderedPair<TransformFlags, Transformation2D>> result = new LinkedList<>();
// Translation is easy when it is considered to be applied last, which we do here.
if (!TrafoUtil.equals(columnW, IVector2D.ZERO)) {
result.add(OrderedPair.create(TransformFlags.Translation,
createTranslation(columnW)));
}
// Now what is left is a 2X2 matrix which is considered to be a combination
// of rotation, shearing and scaling.
// We make use of the fact that for a rotation R the following is always valid:
// R^T = R^{-1}
// We now consider the remaining 2X2 matrix M to be the combination of
// a rotation R and a shear/scale matrix S:
// M = R * S
// If we calculate
// M^T * M = (R * S)^T * (R * S)
// = (S^T * R^T) * R * S
// = S^T * R^{-1} * R * S // because R^T = R^{-1}
// M^T * M = S^T * S
// This means that we can derive S from the last equation. For this we
// assume that S is a triangle matrix in the form
// / a d \
// S = \ 0 b /
// Multiplying S^T * S results in
// / a² ad \
// S^T * S = \ ad b²+d² | = M^T * M
// Resolving this leads to the formulas for a,b, and d below.
// The squares allow to take plus or minus values for the
// scaling coefficients on the diagonal. Here we put a
// the sign (+ or -) into factor a, b is always positive.
// This selects 1 of the 4 possible solutions
// and takes care that any mirroring operation is put into the
// scaling, so all other transformations in the chain have det=1.
final Matrix2X2 m = Matrix2X2.fromColumns(columnX, columnY);
if (det == 1 && m.isIdentity()) {
return result;
}
final Matrix2X2 mTm = m.getTransposed().mult(m);
final double a = det < 0 // a cannot be 0 because of det != 0
? -Math.sqrt(mTm.getXX())
: Math.sqrt(mTm.getXX());
final double d = mTm.getXY() / a;
final double b = Math.sqrt(mTm.getYY() - d*d);
final Matrix2X2 shearAndScale = new Matrix2X2(a, d,
0, b);
// Having shear-and-scale we can calculate the rotation from
// R * S = M <=> R = M * S^{-1}
final Transformation2D sasInverse = shearAndScale.getInverse();
assert sasInverse != null; // as shearAndScale's det != 0
final Transformation2D rotation = m.mult(sasInverse);
if (!rotation.isIdentity()) {
result.add(OrderedPair.create(TransformFlags.Rotation,
new Matrix2X2(rotation)));
}
// Last step: decompose the scaling-and-shearing S.
// Here we assume scaling happening first, which makes S:
// S = shear * scale.
// The scale matrix is obvious and uses a and b. For the shearing
// we assume
// / 1 k \
// shear = | 0 1 |
// Multiplying gives
// / a bk \ / a d \
// shear * scale = | 0 b | = S = | 0 b |
// therefore:
final double k = d / b;
if (k != 0) {
result.add(OrderedPair.create(TransformFlags.Shearing,
new Matrix2X2(1, k,
0, 1)));
}
if (a != 1 || b != 1) {
result.add(OrderedPair.create(TransformFlags.Scaling,
Transformation2D.createScaling(a, b)));
}
return result;
}
```