SVD in 2D

01 Apr 2016, 11:53 UTC

The SVD

The SVD is a matrix decomposition of an arbitrary matrix M into two orthonormal matrices U,V and a diagonal matrix D:

M = UDV

In addition, U and V have determinant +1 (they are rotations, not reflections). Definitions vary but if M has negative determinant we want that represented by negative elements of D.

In 2-dimensions the SVD can be computed very fast in just a few lines of code! This article shows you how!

Initial Observations

U,V are orthonormal, so their transposes are their inverses. U' and V' are these.

So we have

U'MV' = D

Let's break this up into two parts

E = U'M
EV' = D

Now, V' is a simple rotation matrix which operates on the rows of E, rotating each row into rows of D. Since the rows of D are orthogonal (D is diagonal) therefore the rows of E must also be orthogonal. V' then just transforms those rows into the axes. (Note that the *columns* of E are generally not orthogonal.)

Let's express E using variables c3=cos(x), s3=sin(x), a and b as follows (this is justified by the row-orthogonality):

E = (a.c3 -a.s3)
    (b.s3  b.c3)

Then we can just write down V':

(a.c3 -a.s3)( c3  s3) = (a  0)
(b.s3  b.c3)(-s3  c3) = (0  b)

So if we have a matrix E written as follows, and we know that |EF| >= |GH| (i.e. E*E+F*F >= G*G+H*H) then we can easily extract c3/s3:

E = (E  F)
    (G  H)

a.c3 = E
-a.s3 = F

a = sqrt(E*E + F*F)
c3 = E/a
s3 = -F/a

We can now compute an improved estimate of a and b:

(E F)( c3  s3) = (a  0)
(G H)(-s3  c3) = (0  b)

So overall,

V = (c3 -s3)  D = (a  0)
    (s3  c3)      (0  b)

So given E, we can solve for D and V'. We just need to find E, the matrix whose rows are orthogonal, whose top row has largest magnitude, and which can be formed from M as:

E = U'M

So let's express U' using variables c=cos(t), s=sin(t):

E = U'M

M = (A  B)
    (C  D)

(E  F) = ( c s)(A  B)
(G  H) = (-s c)(C  D)

E =  Ac + Cs
F =  Bc + Ds
G = -As + Cc
H = -Bs + Dc

And our orthogonality condition is just:

EG + FH = 0

So we'll just expand this all out as follows:

(Ac + Cs)(-As + Cc) + (Bc + Ds)(-Bs + Dc) = 0

-AAcs + ACcc - ACss + CCsc - BBcs + BDcc - BDss + DDsc = 0

(AC + BD)cc - (AC + BD)ss + ((DD + CC) - (AA + BB))sc = 0

(AC + BD)(cc - ss) + 1/2.((DD + CC) - (AA + BB)).2sc = 0

Now, this turns out to be quite easy to solve, using the double-angle formulae. Recall that

cos(2t) = cos(t)*cos(t) - sin(t)*sin(t)
sin(2t) = 2*cos(t)*sin(t)

Which means we simply have:

(AC + BD)c2 + 1/2.((DD + CC) - (AA + BB)).s2 = 0

This is easily solved by setting

r.c2 = ((AA + BB) - (CC + DD))/2
r.s2 = (AC + BD)

Where r is chosen so c2*c2 + s2*s2 = 1.

Problems occur if r is zero, but it's amazing how that happens.

The bottom value (r.s2) is zero iff the rows are already orthogonal - that is, this step of the computation isn't needed. In this case s2 = 0, c2 = +/-1, so the initial rotation can either do nothing, or it can swap the rows over, which is what we want to give us |EF| >= |GH|.

The top value (r.c2) is zero iff the row vectors have equal length. But if the rows vectors have equal length and the rows are orthogonal, this is a necessary condition for the entire matrix to be orthogonal.

So if we find r = 0, we just set c2 = 1, s2 = 0. We don't need to swap |EF| and |GH| since |EF| = |GH|.

But all this means that r is a measure of the non-orthogonality of the matrix, which is kind of interesting.

Note that we can choose r to be positive or negative. It turns out we want r to be non-negative; this means that c2 is > 0 if |AB| > |CD| and < 0 if |CD| > |AB|

Once we find c2,s2 we can find c and s. There are two ways to do this.

We can find the results from c2:

c2 = c*c - s*s
   = c*c - s*s + (c*c + s*s - 1)
   = 2*c*c - 1
c = sqrt((1 + c2) / 2)

c2 = c*c - s*s
   = c*c - s*s - (c*c + s*s - 1)
   = -2*s*s + 1
s = sqrt((1 - c2) / 2)

note: we choose the signs so that c is >= 0, and then since s2 = 2.c.s, we choose
the sign of s such that 2.c.s has the correct sign

We can also find the results from s2:

s2 = 2cs
1 + s2 = cc + 2cs + ss = (c+s)^2
1 - s2 = cc - 2cs + ss = (c-s)^2

c+s = sqrt(1 + s2)
c-s = sqrt(1 - s2)

note: c2 = cc - ss, so we must choose one sign here such that cc > ss (if c2 > 0) or such that cc < ss (if c2 < 0) - this sign choice amounts to swapping c and s.
then, to match the above, we choose the other sign so that c is positive; this sign choice amounts to negating c and s.

Which formula to use? Well, if c or s is nearly 1, and the other one is not quite zero, there's actually more information in the smaller value. For instance you will find that cosf(epsilon) = 1.0f exactly, while sinf(epsilon) ~= epsilon. The cos is indistinguishable from the cos of other small values, but the sin retains the information about epsilon. This matters, so we compute c and s from c2 or s2 depending on which has the smallest magnitude (fixing this issue took the implementation below from a maximum error of about 0.0002 down to a maximum error of about 1E-6).

So we end up with a value for c and s, with c >= 0, and (here's the tricky bit) if c2 > 0 then |s| < |c|, otherwise |s| > |c| (this follows from c2 = cc - ss).

But we already know that c2 is > 0 if |AB| > |CD| (i.e. A*A + B*B > C*C + D*D), and < 0 if |AB| < |CD|.

Therefore |AB| > |CD| <=> |c| > |s|.

Now we are ready to construct matrices E and U:

( c  s)(A  B) = (E  F)
(-s  c)(C  D)   (G  H)

E =  Ac + Cs
F =  Bc + Ds
G = -As + Cc
H = -Bs + Dc

U = (c -s)
    (s  c)

Now, suppose |AB| > |CD|, then we have c >= 0, and |s| < |c|, so EF is some amount of AB plus some smaller amount of CD, while GH is some amount of CD and some smaller amount of AB, which gives us that |EF| > |GH|. The converse is also true. So we have automatically constructed the larger axis into the top row, which means that in our diagonal matrix the larger value will appear first! (This explanation is a bit hand-wavy, but in the code below in fact |EF| >= |GH| for all millions of tested examples.)

And this constructs the 2D SVD. Here is the actual code. This has been tested on millions of random matrices, and a combinatorial set of matrices from a fixed set of numbers, which includes many nearly-singular and truly singular matrices.

#include <math.h>

// SVD 2D
//
// given an *arbitrary* 2D matrix M, find orthonormal matrices U,V and diagonal matrix D s.t.
//   M = UDV
//   D[0] >= 0
//   D[0] >= fabsf(D[3])

// matrices are stored as follows:
//   0 1
//   2 3

void SVD_2D(const float M[4], float outU[4], float outD[4], float outV[4])
{
    float   A = M[0], B = M[1], C = M[2], D = M[3];
    float   E,F,G,H;

    // precondition A,B,C,D
    // - this helps if all values are very small, because otherwise
    //   we can suffer exponent underflow when squaring
    float   largest = fabsf(M[0]);
    for (unsigned int ii = 0; ii < 4; ii++) {
        if (fabsf(M[ii]) > largest) {
            largest = fabsf(M[ii]);
        }
    }
    if (largest && (largest < 1)) {
        A /= largest;
        B /= largest;
        C /= largest;
        D /= largest;
    }

    // determine l*cos(2t) and l*sin(2t)
    float   lc2 = ((A*A + B*B) - (C*C + D*D)) * 0.5f;
    float   ls2 = A*C + B*D;

    float   magsq2 = lc2*lc2 + ls2*ls2;
    if (magsq2 == 0) {
        // matrix is precisely orthogonal already - U = I, EFGH = M
        outU[0] = 1; outU[1] = 0; outU[2] = 0; outU[3] = 1;
        E = A; F = B; G = C; H = D;
    } else {
        // determine cos(2t) and sin(2t)
        float   mag = sqrtf(magsq2);
        float   c2 = lc2 / mag;     // c2 is +ve when |AB| > |CD|, -ve when |AB| < |CD|
        float   s2 = ls2 / mag;

        float   c,s;

        if (fabsf(s2) < fabsf(c2)) {
            // a smaller magnitude has more precision
            // - compute c,s from s2 rather than c2 in this case

            // s2 = 2cs
            // 1 + s2 = cc + 2cs + ss = (c+s)^2
            // 1 - s2 = cc - 2cs + ss = (c-s)^2
            float   cps = sqrtf(1 + s2);    // +/-
            float   cms = sqrtf(1 - s2);    // +/-
          
            c = (cps + cms) * 0.5f;
            s = (cps - cms) * 0.5f;
            if (((c*c - s*s > 0) && (c2 < 0)) ||
                ((c*c - s*s < 0) && (c2 > 0))) {
                // swap c & s - corresponds to choosing -sign for cms
                float tmp = c; c = s; s = tmp;
            }
            // make c positive
            if (c < 0) {
                // negate c & s - corresponds to choosing -sign for cps
                c = -c;
                s = -s;
            }
        } else {
            // determine cos(t) and sin(t)
            c = sqrtf((1 + c2) * 0.5f); // c > s <=> |AB| > |CD|
            s = sqrtf((1 - c2) * 0.5f);

            // swap sign of s to match sin(2t) = 2*sin(t)*cos(t) (cos(t) >= 0)
            if (s2 < 0) {
                s = -s;
            }
        }

        outU[0] = c; outU[1] = -s; outU[2] = s; outU[3] = c;
        E = c*A + s*C;      // c > s <=> EF contains more AB and less CD
        F = c*B + s*D;      // c > s <=> GH contains more CD and less AB
        G = -s*A + c*C;     // so |EF| > |GH| here
        H = -s*B + c*D;
    }

    // since EF is the largest vector, just scale by this length
    float   mag = sqrtf(E*E + F*F);
    float   c3,s3;
    if (mag > 0) {
        c3 = E / mag;
        s3 = -F / mag;
    } else {
        // the matrix is zero
        c3 = 1; s3 = 0;
    }

    float   D0 = E*c3 - F*s3;
    float   D3 = G*s3 + H*c3;

    // undo preconditioning
    if (largest && (largest < 1)) {
        D0 *= largest;
        D3 *= largest;
    }

    outV[0] = c3; outV[1] = -s3; outV[2] = s3; outV[3] = c3;
    outD[0] = D0; outD[1] = 0; outD[2] = 0; outD[3] = D3;
}