Polar Decomposition

To decompose \(X\) so that \(X = QS\) we need set \(Q\) equal to \(X\) and then average \(Q\) with its inverse transpose repeatedly. Doing so over and over will eventually converge on a polar rotation, so long as the determinant of \(Q\) was positive. To find \(S\) you just need to multiply \(X\) by the inverse of \(Q\).

To actually implement this, set Q equal to X

$$ Q = X $$

Then, repeatedly average \(Q\) with its inverse transpose

$$ Q_{i + 1} = \frac{1}{2}(Q_{i} + (Q_{i}^{-1})^{T}) $$

Keep doing this until \(Q_{i + 1} - Q_{i} \approx 0\). Alternatley you can do a fixed number of iterations, 20 steps should be more than enough. After you know what \(Q\) is, finding \(S\) is trivial, just multiply the input matrix \(X\) by the inverse of \(Q\):

$$ S = Q^{-1} X $$

Below is an interactive example. The Red, Green and Blue lines represent the right, up and forward basis vectors of a transform matrix. This transform matrix contains both rotation and scale, causing the object to look skewed. The value of the slider in the upper left corner represents the current iteration out of a maximum of 9. Noice how past iteration 5 you can not visually tell that the matrix is being adjusted.

Canvas support required WebGL support required

Understanding how the iterative algorithm works, we can actually implement polar decompositon in code.

struct PolarDecompResult {
    Matrix Q; // Q is a rotation, or the negative of a rotation
    Matrix S; // Scale and skew matrix
}

bool PolarDecompositionEarlyOut(Matrix Q, Matrix Qprev) {
  Matrix res = Q - Qprev;

  for (var i = 0; i < 16; ++i) {
    if (!EpsilonCompare(res[i], 0)) {
      return false;
    }
  }

  return true;
}

PolarDecompResult PolarDecomposition(Matrix X) {
    PolarDecompResult result;

    Matrix Q = X;

    for (int i = 0; i < 20; ++i) {
        Matrix Qprev = Q;

        Matrix Qit = Inverse(Transpose(Q));
        Q = MulMatFloat(AddMatMat(Q, Qit), 0.5)

        if (PolarDecompositionEarlyOut(Q, Qprev)) {
          break;
      }
    }


    result.Q = Q;
    result.S = Inverse(Q) * X;

    return result;
}

Remember, \(Q\) can be further factored into a flip and a rotation matrix. If the determinant of \(Q\) is negative, it contains a flip.

struct FactorRotationResult {
    Matrix F; // Positive or negative identity (flip if negative)
    Matrix R; // Rotation matrix
}

FactorRotationResult FactorRotation(Matrix Q) {
    FactorRotationResult result;
    float f = 1.0;

    if (Det(Q) < 0) {
        f = -1;
        Q = Q * -1; // Component wise multiplication
    }

    result.F[0] = f;
    result.F[5] = f;
    result.F[10] = f;

    result.R = Q;

    return result;
}