26
$\begingroup$

What is a simple algorithm for computing the SVD of $2 \times 2$ matrices?

Ideally, I'd like a numerically robust algorithm, but I'll like to see both simple and not-so-simple implementations. C code accepted.

Any references to papers or code?

$\endgroup$
  • 5
    $\begingroup$ Wikipedia lists a 2x2 closed form solution, but I have no idea of its numerical properties. $\endgroup$ – Damien Oct 28 '13 at 22:57
  • $\begingroup$ As reference, "Numerical Recipes", Press et al., Cambridge Press. Quite expensive book but worth every cent. Besides SVD solutions you will find a lot of other useful algorithms. $\endgroup$ – Jan Hackenberg Aug 16 '16 at 15:46

10 Answers 10

19
$\begingroup$

See https://math.stackexchange.com/questions/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (sorry, I would have put that in a comment but I've registered just to post this so I can't post comments yet).

But since I'm writing it as an answer, I'll also write the method:

$$E=\frac{m_{00}+m_{11}}{2}; F=\frac{m_{00}-m_{11}}{2}; G=\frac{m_{10}+m_{01}}{2}; H=\frac{m_{10}-m_{01}}{2}\\ Q=\sqrt{E^2+H^2}; R=\sqrt{F^2+G^2}\\ s_x=Q+R; s_y=Q-R\\ a_1=\mathrm{atan2}(G,F); a_2=\mathrm{atan2}(H,E)\\ \theta=\frac{a_2-a_1}{2}; \phi=\frac{a_2+a_1}{2}$$

That decomposes the matrix as follows:

$$M=\pmatrix{m_{00}&m_{01}\\m_{10}&m_{11}}=\pmatrix{\cos\phi&-\sin\phi\\\sin\phi&\cos\phi}\pmatrix{s_x&0\\0&s_y}\pmatrix{\cos\theta&-\sin\theta\\\sin\theta&\cos\theta}$$

The only thing to guard against with this method is that $G=F=0$ or $H=E=0$ for atan2. I doubt it can be any more robust than that (Update: see Alex Eftimiades' answer!).

The reference is: http://dx.doi.org/10.1109/38.486688 (given by Rahul there) which comes from the bottom of this blog post: http://metamerist.blogspot.com/2006/10/linear-algebra-for-graphics-geeks-svd.html

Update: As noted by @VictorLiu in a comment, $s_y$ may be negative. That happens if and only if the determinant of the input matrix is negative as well. If that's the case and you want the positive singular values, just take the absolute value of $s_y$.

$\endgroup$
  • 1
    $\begingroup$ It seems that $s_y$ can be negative if $Q<R$. This should not be possible. $\endgroup$ – Victor Liu Jul 11 '14 at 17:27
  • $\begingroup$ @VictorLiu If the input matrix flips, the only place that can be reflected is in the scaling matrix, as the rotation matrices can't possibly flip. Just don't feed it input matrices that flip. I haven't done the math yet but I bet that the sign of the determinant of the input matrix will determine whether $Q$ or $R$ is greater. $\endgroup$ – Pedro Gimeno Jul 11 '14 at 22:47
  • $\begingroup$ @VictorLiu I've done the math now and confirmed that indeed, $Q^2-R^2$ simplifies to $m_{00}m_{11}-m_{01}m_{10}$ i.e. the determinant of the input matrix. $\endgroup$ – Pedro Gimeno Jul 12 '14 at 1:09
  • $\begingroup$ Oh hi. (+1) $ $ $\endgroup$ – Rahul Dec 15 '14 at 6:45
9
$\begingroup$

@Pedro Gimeno

"I doubt it can be any more robust than that."

Challenge accepted.

I noticed the usual approach is to use trig functions like atan2. Intuitively, there shouldn't be a need to use trig functions. Indeed, all the results end up as sines and cosines of arctans--which can be simplified to algebraic functions. It took quite a while, but I managed to simplified Pedro's algorithm to use only algebraic functions.

The following python code does the trick.

from numpy import asarray, diag

def svd2(m):

y1, x1 = (m[1, 0] + m[0, 1]), (m[0, 0] - m[1, 1]) y2, x2 = (m[1, 0] - m[0, 1]), (m[0, 0] + m[1, 1]) h1 = hypot(y1, x1) h2 = hypot(y2, x2) t1 = x1 / h1 t2 = x2 / h2 cc = sqrt((1 + t1) * (1 + t2)) ss = sqrt((1 - t1) * (1 - t2)) cs = sqrt((1 + t1) * (1 - t2)) sc = sqrt((1 - t1) * (1 + t2)) c1, s1 = (cc - ss) / 2, (sc + cs) / 2, u1 = asarray([[c1, -s1], [s1, c1]]) d = asarray([(h1 + h2) / 2, (h1 - h2) / 2]) sigma = diag(d) if h1 != h2: u2 = diag(1 / d).dot(u1.T).dot(m) else: u2 = diag([1 / d[0], 0]).dot(u1.T).dot(m) return u1, sigma, u2

$\endgroup$
  • 1
    $\begingroup$ The code seems incorrect. Consider the 2x2 identity matrix. Then y1=0, x1=0, h1=0, and t1=0/0=NaN. $\endgroup$ – Hugues Nov 16 '17 at 21:24
8
$\begingroup$

The GSL has a 2-by-2 SVD solver underlying the QR decomposition part of the main SVD algorithm for gsl_linalg_SV_decomp. See the svdstep.c file and look for the svd2 function. The function has a few special cases, isn't exactly trivial, and looks to be doing several things to be numerically careful (e.g., using hypot to avoid overflows).

$\endgroup$
  • 1
    $\begingroup$ Does this function have any documentation? I would like to know what its input parameters are. $\endgroup$ – Victor Liu Oct 29 '13 at 23:05
  • $\begingroup$ @VictorLiu: Sadly I haven't seen anything other than the meager comments in the file itself. There's a bit in the ChangeLog file if you download the GSL. And you can look at svd.c for details of the overall algorithm. The only true documentation seems to be for the high level user-callable functions, e.g., gsl_linalg_SV_decomp. $\endgroup$ – horchler Oct 30 '13 at 14:06
7
$\begingroup$

When we say "numerically robust" we usually mean an algorithm in which we do things like pivoting to avoid error propagation. However, for a 2x2 matrix, you can write the result down in terms of explicit formulas -- i.e., write down formulas for the SVD elements that state the result only in terms of the inputs, rather than in terms of intermediate values previously computed. That means that you may have cancellation but no error propagation.

The point simply is that for 2x2 systems, worrying about robustness is not necessary.

$\endgroup$
  • $\begingroup$ It can depend on the matrix. I've seen a method which finds the left and right angles separately (each via arctan2(y,x)) which generally works fine. But when the singular values are close together, each of these arctans tends to 0/0, so the result can be inaccurate. In the method given by Pedro Gimeno, the calculation of a2 will be well defined in this case, while a1 becomes ill-defined; you still have a good result because the validity of the decomposition is only sensitive to theta+phi when the s.vals are close together, not to theta-phi. $\endgroup$ – greggo Mar 3 '17 at 20:19
5
$\begingroup$

This code is based on Blinn's paper, Ellis paper, SVD lecture, and additional calculations. An algorithm is suitable for regular and singular real matrices. All previous versions works 100% as well as this one.

#include <stdio.h>
#include <math.h>

void svd22(const double a[4], double u[4], double s[2], double v[4]) {
    s[0] = (sqrt(pow(a[0] - a[3], 2) + pow(a[1] + a[2], 2)) + sqrt(pow(a[0] + a[3], 2) + pow(a[1] - a[2], 2))) / 2;
    s[1] = fabs(s[0] - sqrt(pow(a[0] - a[3], 2) + pow(a[1] + a[2], 2)));
    v[2] = (s[0] > s[1]) ? sin((atan2(2 * (a[0] * a[1] + a[2] * a[3]), a[0] * a[0] - a[1] * a[1] + a[2] * a[2] - a[3] * a[3])) / 2) : 0;
    v[0] = sqrt(1 - v[2] * v[2]);
    v[1] = -v[2];
    v[3] = v[0];
    u[0] = (s[0] != 0) ? (a[0] * v[0] + a[1] * v[2]) / s[0] : 1;
    u[2] = (s[0] != 0) ? (a[2] * v[0] + a[3] * v[2]) / s[0] : 0;
    u[1] = (s[1] != 0) ? (a[0] * v[1] + a[1] * v[3]) / s[1] : -u[2];
    u[3] = (s[1] != 0) ? (a[2] * v[1] + a[3] * v[3]) / s[1] : u[0];
}

int main() {
    double a[4] = {1, 2, 3, 6}, u[4], s[2], v[4];
    svd22(a, u, s, v);
    printf("Matrix A:\n%f %f\n%f %f\n\n", a[0], a[1], a[2], a[3]);
    printf("Matrix U:\n%f %f\n%f %f\n\n", u[0], u[1], u[2], u[3]);
    printf("Matrix S:\n%f %f\n%f %f\n\n", s[0], 0, 0, s[1]);
    printf("Matrix V:\n%f %f\n%f %f\n\n", v[0], v[1], v[2], v[3]);
}
$\endgroup$
5
$\begingroup$

I needed an algorithm that has

  • little branching (hopefully CMOVs)
  • no trigonometric function calls
  • high numerical accuracy even with 32 bit floats

We want to calculate $c_1, s_1, c_2, s_2, \sigma_1$ and $\sigma_2$ as follows:

$A = USV$, which can be expanded like:

$ \begin{bmatrix} a & b \\ c & d \end{bmatrix} = \begin{bmatrix} c_1 & s_1 \\ -s_1 & c_1 \end{bmatrix} \begin{bmatrix} \sigma_1 & 0 \\ 0 & \sigma_2 \end{bmatrix} \begin{bmatrix} c_2 & -s_2 \\ s_2 & c_2 \end{bmatrix} $

The main idea is to find a rotation matrix $V$ that diagonalizes $A^TA$, that is $VA^TAV^T=D$ is diagonal.

Recall that

$USV = A$

$US = AV^{-1} = AV^T$ (since $V$ is orthogonal)

$VA^TAV^T = (AV^T)^TAV^T = (US)^TUS = S^TU^TUS = D$

Multiplying both sides by $S^{-1}$ we get

$(S^{-T}S^T)U^TU(SS^{-1}) = U^TU = S^{-T}DS^{-1}$

Since $D$ is diagonal, setting $S$ to $\sqrt{D}$ will give us $U^TU=Identity$, meaning $U$ is a rotation matrix, $S$ is a diagonal matrix, $V$ is a rotation matrix and $USV = A$, just what we are looking for.

Calculating the diagonalizing rotation can be done by solving the following equation:

$t_2^2 - \frac{\beta-\alpha}{\gamma}t_2-1 = 0$

where

$ A^TA = \begin{bmatrix} a & c \\ b & d \end{bmatrix} \begin{bmatrix} a & b \\ c & d \end{bmatrix} = \begin{bmatrix} a^2+c^2 & ab+cd \\ ab+cd & b^2+d^2 \end{bmatrix} = \begin{bmatrix} \alpha & \gamma \\ \gamma & \beta \end{bmatrix} $

and $t_2$ is the tangent of angle of $V$. This can be derived by expanding $VA^TAV^T$ and making its off-diagonal elements equal to zero (they are equal to each other).

The problem with this method is that it loses significant floating point precision when calculating $\beta-\alpha$ and $\gamma$ for certain matrices, because of the subtractions in the calculations. The solution for this is to do an RQ decomposition ($A=RQ$, $R$ upper triangular and $Q$ orthogonal) first, then use the algorithm to factorize $USV' = R$. This gives $USV=USV'Q=RQ=A$. Notice how setting $d$ to 0 (as in $R$) eliminates some of the additions/subtractions. (The RQ decomposition is fairly trivial from the expansion of the matrix product).

The algorithm naively implemented this way has some numerical and logical anomalies (e.g. is $S$ $+\sqrt{D}$ or $-\sqrt{D}$), which I fixed in the code below.

I threw about 2000 million randomized matrices at the code, and the largest numerical error produced was around $6\cdot10^{-7}$ (with 32 bit floats, $error = ||USV-M||/||M||$). The algorithm runs in about 340 clock cycles (MSVC 19, Ivy Bridge).

template <class T>
void Rq2x2Helper(const Matrix<T, 2, 2>& A, T& x, T& y, T& z, T& c2, T& s2) {
    T a = A(0, 0);
    T b = A(0, 1);
    T c = A(1, 0);
    T d = A(1, 1);

    if (c == 0) {
        x = a;
        y = b;
        z = d;
        c2 = 1;
        s2 = 0;
        return;
    }
    T maxden = std::max(abs(c), abs(d));

    T rcmaxden = 1/maxden;
    c *= rcmaxden;
    d *= rcmaxden;

    T den = 1/sqrt(c*c + d*d);

    T numx = (-b*c + a*d);
    T numy = (a*c + b*d);
    x = numx * den;
    y = numy * den;
    z = maxden/den;

    s2 = -c * den;
    c2 = d * den;
}


template <class T>
void Svd2x2Helper(const Matrix<T, 2, 2>& A, T& c1, T& s1, T& c2, T& s2, T& d1, T& d2) {
    // Calculate RQ decomposition of A
    T x, y, z;
    Rq2x2Helper(A, x, y, z, c2, s2);

    // Calculate tangent of rotation on R[x,y;0,z] to diagonalize R^T*R
    T scaler = T(1)/std::max(abs(x), abs(y));
    T x_ = x*scaler, y_ = y*scaler, z_ = z*scaler;
    T numer = ((z_-x_)*(z_+x_)) + y_*y_;
    T gamma = x_*y_;
    gamma = numer == 0 ? 1 : gamma;
    T zeta = numer/gamma;

    T t = 2*impl::sign_nonzero(zeta)/(abs(zeta) + sqrt(zeta*zeta+4));

    // Calculate sines and cosines
    c1 = T(1) / sqrt(T(1) + t*t);
    s1 = c1*t;

    // Calculate U*S = R*R(c1,s1)
    T usa = c1*x - s1*y; 
    T usb = s1*x + c1*y;
    T usc = -s1*z;
    T usd = c1*z;

    // Update V = R(c1,s1)^T*Q
    t = c1*c2 + s1*s2;
    s2 = c2*s1 - c1*s2;
    c2 = t;

    // Separate U and S
    d1 = std::hypot(usa, usc);
    d2 = std::hypot(usb, usd);
    T dmax = std::max(d1, d2);
    T usmax1 = d2 > d1 ? usd : usa;
    T usmax2 = d2 > d1 ? usb : -usc;

    T signd1 = impl::sign_nonzero(x*z);
    dmax *= d2 > d1 ? signd1 : 1;
    d2 *= signd1;
    T rcpdmax = 1/dmax;

    c1 = dmax != T(0) ? usmax1 * rcpdmax : T(1);
    s1 = dmax != T(0) ? usmax2 * rcpdmax : T(0);
}

Ideas from:
http://www.cs.utexas.edu/users/inderjit/public_papers/HLA_SVD.pdf
http://www.math.pitt.edu/~sussmanm/2071Spring08/lab09/index.html
http://www.lucidarme.me/singular-value-decomposition-of-a-2x2-matrix/

$\endgroup$
3
$\begingroup$

I have used the description at http://www.lucidarme.me/?p=4624 to create this C++ code. The Matrices are those of the Eigen library, but you can easily create your own data structure from this example:

$A=U\Sigma V^T$

#include <cmath>
#include <Eigen/Core>
using namespace Eigen;

Matrix2d A;
// ... fill A

double a = A(0,0);
double b = A(0,1);
double c = A(1,0);
double d = A(1,1);

double Theta = 0.5 * atan2(2*a*c + 2*b*d,
                           a*a + b*b - c*c - d*d);
// calculate U
Matrix2d U;
U << cos(Theta), -sin(Theta), sin(Theta), cos(Theta);

double Phi = 0.5 * atan2(2*a*b + 2*c*d,
                         a*a - b*b + c*c - d*d);
double s11 = ( a*cos(Theta) + c*sin(Theta))*cos(Phi) +
             ( b*cos(Theta) + d*sin(Theta))*sin(Phi);
double s22 = ( a*sin(Theta) - c*cos(Theta))*sin(Phi) +
             (-b*sin(Theta) + d*cos(Theta))*cos(Phi);

// calculate S
S1 = a*a + b*b + c*c + d*d;
S2 = sqrt(pow(a*a + b*b - c*c - d*d, 2) + 4*pow(a*c + b*d, 2));

Matrix2d Sigma;
Sigma << sqrt((S1+S2) / 2), 0, 0, sqrt((S1-S2) / 2);

// calculate V
Matrix2d V;
V << signum(s11)*cos(Phi), -signum(s22)*sin(Phi),
     signum(s11)*sin(Phi),  signum(s22)*cos(Phi);

With the standard sign function

double signum(double value)
{
    if(value > 0)
        return 1;
    else if(value < 0)
        return -1;
    else
        return 0;
}

This results in exactly the same values as the Eigen::JacobiSVD (see https://eigen.tuxfamily.org/dox-devel/classEigen_1_1JacobiSVD.html).

$\endgroup$
  • 1
    $\begingroup$ S2 = hypot( a*a + b*b - c*c - d*d, 2*(a*c + b*d)) $\endgroup$ – greggo Mar 3 '17 at 20:27
2
$\begingroup$

I have pure C code for the 2x2 real SVD here. See line 559. It essentially computes the eigenvalues of $A^TA$ by solving a quadratic, so it's not necessarily the most robust, but it seems to work well in practice for not-too-pathological cases. It's relatively simple.

$\endgroup$
  • $\begingroup$ I don't think your code works when the eigenvalues of the matrix are negative. Try [[1 1] [1 0]], and u * s * vt is not equal to m... $\endgroup$ – Carlos Scheidegger Feb 20 '14 at 3:13
2
$\begingroup$

For my personal need, I tried to isolate the minimum computation for a 2x2 svd. I guess it is probably one of the simplest and fastest solution. You can find details on my personal blog : http://lucidarme.me/?p=4624.

Advantages : simple, fast and you can only calculate one or two of the three matrices (S, U or D) if you don't need the three matrices.

Drawback it uses atan2, which may be inacurate and may require an external library (typ. math.h).

$\endgroup$
  • 3
    $\begingroup$ Since links are rarely permanent, it is important to summarize the approach rather than simply providing a link as an answer. $\endgroup$ – Paul Dec 14 '14 at 22:59
  • $\begingroup$ Also, if you're going to post a link to your own blog, please (a) disclose that it's your blog, (b) even better would be to actually summarize or cut-and-paste your approach (the images of formulas can be translated into raw LaTeX and rendered using MathJax). The best answers for this sort of question state formulas, provide citations for said formulas, and then list things like drawbacks, edge cases, and potential alternatives. $\endgroup$ – Geoff Oxberry Dec 15 '14 at 5:51
1
$\begingroup$

Here is an implementation of a 2x2 SVD solve. I based it off of Victor Liu's code. His code was not working for some matrices. I used these two documents as mathematical reference for the solve: pdf1 and pdf2.

The matrix setData method is in row-major order. Internally, I represent the matrix data as a 2D array given by data[col][row].

void Matrix2f::svd(Matrix2f* w, Vector2f* e, Matrix2f* v) const{
    //If it is diagonal, SVD is trivial
    if (fabs(data[0][1] - data[1][0]) < EPSILON && fabs(data[0][1]) < EPSILON){
        w->setData(data[0][0] < 0 ? -1 : 1, 0, 0, data[1][1] < 0 ? -1 : 1);
        e->setData(fabs(data[0][0]), fabs(data[1][1]));
        v->loadIdentity();
    }
    //Otherwise, we need to compute A^T*A
    else{
        float j = data[0][0]*data[0][0] + data[0][1]*data[0][1],
            k = data[1][0]*data[1][0] + data[1][1]*data[1][1],
            v_c = data[0][0]*data[1][0] + data[0][1]*data[1][1];
        //Check to see if A^T*A is diagonal
        if (fabs(v_c) < EPSILON){
            float s1 = sqrt(j),
                s2 = fabs(j-k) < EPSILON ? s1 : sqrt(k);
            e->setData(s1, s2);
            v->loadIdentity();
            w->setData(
                data[0][0]/s1, data[1][0]/s2,
                data[0][1]/s1, data[1][1]/s2
            );
        }
        //Otherwise, solve quadratic for eigenvalues
        else{
            float jmk = j-k,
                jpk = j+k,
                root = sqrt(jmk*jmk + 4*v_c*v_c),
                eig = (jpk+root)/2,
                s1 = sqrt(eig),
                s2 = fabs(root) < EPSILON ? s1 : sqrt((jpk-root)/2);
            e->setData(s1, s2);
            //Use eigenvectors of A^T*A as V
            float v_s = eig-j,
                len = sqrt(v_s*v_s + v_c*v_c);
            v_c /= len;
            v_s /= len;
            v->setData(v_c, -v_s, v_s, v_c);
            //Compute w matrix as Av/s
            w->setData(
                (data[0][0]*v_c + data[1][0]*v_s)/s1,
                (data[1][0]*v_c - data[0][0]*v_s)/s2,
                (data[0][1]*v_c + data[1][1]*v_s)/s1,
                (data[1][1]*v_c - data[0][1]*v_s)/s2
            );
        }
    }
}
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.