 # 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
//   D >= fabsf(D)

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

void SVD_2D(const float M, float outU, float outD, float outV)
{
float   A = M, B = M, C = M, D = M;
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);
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 = 1; outU = 0; outU = 0; outU = 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 = c; outU = -s; outU = s; outU = 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 = c3; outV = -s3; outV = s3; outV = c3;
outD = D0; outD = 0; outD = 0; outD = D3;
}
``````

This is like listening to a bunch of southern whites complain about the &quot;uppity niqr&gseguot; who commit many crimes, and then complaining that those coons are so sensitive to the truth. commented:

This is like listening to a bunch of southern whites complain about the &quot;uppity niqr&gseguot; who commit many crimes, and then complaining that those coons are so sensitive to the truth.

on 09 May 2016, 09:56 UTC

DcqYji http://www.FyLitCl7Pf7kjQdDUOLQOuaxTXbj5iNG.com commented:

DcqYji http://www.FyLitCl7Pf7kjQdDUOLQOuaxTXbj5iNG.com

on 13 Aug 2016, 13:15 UTC

TG788p <a href="http://enoifzdudjau.com/">enoifzdudjau</a>, [url=http://mpwacrpzbvgd.com/]mpwacrpzbvgd[/url], [link=http://dupzvusagucr.com/]dupzvusagucr[/link], http://pfpfhawegera.com/ commented:

on 15 Aug 2016, 10:26 UTC

DNVUSJ <a href="http://cykkmciookrs.com/">cykkmciookrs</a>, [url=http://clroqlqimiav.com/]clroqlqimiav[/url], [link=http://ggrtymajumne.com/]ggrtymajumne[/link], http://abrbmufcqmrd.com/ commented:

on 15 Aug 2016, 10:32 UTC

jvJF5R <a href="http://jurnfuaajpdi.com/">jurnfuaajpdi</a>, [url=http://kgobgvguuzfe.com/]kgobgvguuzfe[/url], [link=http://otioghhszjdl.com/]otioghhszjdl[/link], http://cnehphgnxxvu.com/ commented:

on 15 Aug 2016, 10:51 UTC

biRddl <a href="http://uosdfofjftdd.com/">uosdfofjftdd</a>, [url=http://ijwfqqddjpny.com/]ijwfqqddjpny[/url], [link=http://frunafgelptq.com/]frunafgelptq[/link], http://xuxaeieyrgbn.com/ commented:

on 15 Aug 2016, 12:51 UTC

eKUHQ2 <a href="http://ziqoiqcsfgcw.com/">ziqoiqcsfgcw</a>, [url=http://wrahbsysamcz.com/]wrahbsysamcz[/url], [link=http://jvqpoaqxlxkp.com/]jvqpoaqxlxkp[/link], http://xtfdidmsrymf.com/ commented:

on 15 Aug 2016, 13:13 UTC

g0Mbkm <a href="http://mixbpntqidvi.com/">mixbpntqidvi</a>, [url=http://bcaiwbxhwtql.com/]bcaiwbxhwtql[/url], [link=http://ohcmjzpdeyce.com/]ohcmjzpdeyce[/link], http://lsjhtlaypozm.com/ commented:

on 15 Aug 2016, 15:21 UTC

nhDkLw <a href="http://lzzlfyaflejl.com/">lzzlfyaflejl</a>, [url=http://bozdeuualgaw.com/]bozdeuualgaw[/url], [link=http://mzmblosgiuwk.com/]mzmblosgiuwk[/link], http://gyycigxioexw.com/ commented:

on 15 Aug 2016, 15:40 UTC

I'd like to withdraw \$100, please commented:

I'd like to withdraw \$100, please

on 31 Aug 2016, 16:54 UTC

I'm on business commented:

on 31 Aug 2016, 16:54 UTC

Could you please repeat that? commented:

Could you please repeat that?

on 31 Aug 2016, 16:54 UTC

Would you like a receipt? commented:

Would you like a receipt?

on 31 Aug 2016, 16:54 UTC

I support Manchester United commented:

I support Manchester United

on 31 Aug 2016, 16:54 UTC

Could you transfer \$1000 from my current account to my deposit account? commented:

Could you transfer \$1000 from my current account to my deposit account?

on 31 Aug 2016, 16:54 UTC

I was born in Australia but grew up in England commented:

I was born in Australia but grew up in England

on 31 Aug 2016, 16:54 UTC

Punk not dead commented:

on 31 Aug 2016, 16:54 UTC

What's the last date I can post this to to arrive in time for Christmas? commented:

What's the last date I can post this to to arrive in time for Christmas?

on 31 Aug 2016, 16:54 UTC

I'm in a band commented:

I'm in a band

on 31 Aug 2016, 16:54 UTC

I want to report a commented:

I want to report a

on 31 Aug 2016, 20:42 UTC

I need to charge up my phone commented:

I need to charge up my phone

on 31 Aug 2016, 20:42 UTC

I'd like to open a business account commented:

I'd like to open a business account

on 31 Aug 2016, 20:42 UTC

I'd like to send this parcel to commented:

I'd like to send this parcel to

on 31 Aug 2016, 20:42 UTC

Could you please repeat that? commented:

Could you please repeat that?

on 31 Aug 2016, 20:43 UTC

We need someone with qualifications commented:

We need someone with qualifications

on 31 Aug 2016, 20:43 UTC

I'd like to order some foreign currency commented:

I'd like to order some foreign currency

on 31 Aug 2016, 20:43 UTC

Thanks for calling commented:

Thanks for calling

on 31 Aug 2016, 20:43 UTC

Looking for work commented:

Looking for work

on 31 Aug 2016, 20:43 UTC

History commented:

History

on 31 Aug 2016, 20:43 UTC

nN1ENQ <a href="http://wwwaiphoavnk.com/">wwwaiphoavnk</a>, [url=http://zifldhmycpmh.com/]zifldhmycpmh[/url], [link=http://mgjjuxvaanik.com/]mgjjuxvaanik[/link], http://ylmaqknmcbwd.com/ commented: