*01 Apr 2016, 11:53 UTC*

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!

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;
}
```

**This is like listening to a bunch of southern whites complain about the "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 "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:

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

*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:

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

*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:

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

*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:

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

*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:

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

*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:

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

*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:

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

*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:

I'm on business

*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:

Punk not dead

*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:

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

*on 05 Sep 2016, 18:06 UTC*