Robuster Algorithmus für

Was ist ein einfacher Algorithmus zur Berechnung der SVD von $2 \times 2$ Matrizen?

Idealerweise hätte ich gerne einen numerisch robusten Algorithmus, aber ich würde gerne sowohl einfache als auch nicht so einfache Implementierungen sehen. C-Code akzeptiert.

Verweise auf Papiere oder Code?

Siehe https://math.stackexchange.com/questions/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (Entschuldigung, das hätte ich in einen Kommentar eingefügt, aber ich habe mich registriert nur um dies zu posten, damit ich noch keine Kommentare posten kann).

Da ich es aber als Antwort schreibe, schreibe ich auch die Methode:

Das zerlegt die Matrix wie folgt:

Das einzige, wovor man sich bei dieser Methode hüten muss, ist, dass $G=F=0$ oder $H=E=0$ für atan2 ist.Ich bezweifle, dass es robuster sein kann( Update: siehe Antwort von Alex Eftimiades!)

Die Referenz lautet: http://dx.doi.org/10.1109/38.486688 (von Rahul dort angegeben) und befindet sich am Ende dieses Blogposts: http://metamerist.blogspot.com/2006/10/linear-algebra -für-Grafik-Geeks-svd.html

Update: Wie von @VictorLiu in einem Kommentar vermerkt, kann $s_y$ negativ sein. Dies geschieht genau dann, wenn die Determinante der Eingangsmatrix ebenfalls negativ ist. Wenn dies der Fall ist und Sie die positiven Singularwerte wollen, nehmen Sie einfach den absoluten Wert von $s_y$ .

@ Pedro Gimeno

"Ich bezweifle, dass es robuster sein kann."

Herausforderung angenommen.

Mir ist aufgefallen, dass der übliche Ansatz darin besteht, Triggerfunktionen wie atan2 zu verwenden. Intuitiv sollte es nicht erforderlich sein, Triggerfunktionen zu verwenden. In der Tat enden alle Ergebnisse als Sinus und Cosinus von Arctanen - was zu algebraischen Funktionen vereinfacht werden kann. Es hat eine Weile gedauert, aber ich habe es geschafft, den Algorithmus von Pedro so zu vereinfachen, dass er nur algebraische Funktionen verwendet.

Der folgende Python-Code erledigt den Trick.

von 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

Die GSL verfügt über einen 2-mal-2-SVD-Löser, der dem QR-Zerlegungsteil des SVD-Hauptalgorithmus für zugrunde liegt gsl_linalg_SV_decomp. Sehen Sie sich die svdstep.cDatei an und suchen Sie nach der svd2Funktion. Die Funktion hat einige Sonderfälle, ist nicht gerade trivial und scheint verschiedene Dinge zu tun, um numerisch vorsichtig zu sein (z. B. hypotum Überläufe zu vermeiden).

Wenn wir "numerisch robust" sagen, meinen wir normalerweise einen Algorithmus, bei dem wir Dinge wie das Schwenken tun, um eine Fehlerausbreitung zu vermeiden. Bei einer 2x2-Matrix können Sie das Ergebnis jedoch in Form expliziter Formeln aufschreiben, dh Formeln für die SVD-Elemente, die das Ergebnis nur anhand der Eingaben und nicht anhand der zuvor berechneten Zwischenwerte angeben . Das bedeutet, dass Sie möglicherweise eine Stornierung, aber keine Fehlerweitergabe haben.

Der Punkt ist einfach, dass bei 2x2-Systemen die Sorge um die Robustheit nicht erforderlich ist.

Dieser Code basiert auf Blinns Artikel , Ellis-Artikel , SVD-Vorlesung und zusätzlichen Berechnungen. Ein Algorithmus eignet sich für reguläre und singuläre reelle Matrizen. Alle vorherigen Versionen funktionieren zu 100% so gut wie diese.

#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]);
}

Ich brauchte einen Algorithmus, der hat

• kleine Verzweigung (hoffentlich CMOVs)
• keine trigonometrischen Funktionsaufrufe
• Hohe numerische Genauigkeit auch bei 32-Bit-Floats

Wir wollen berechnen$c_1, s_1, c_2, s_2, \sigma_1$$\sigma_2$

$A = USV$

$\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}$

$V$$A^TA$$VA^TAV^T=D$

Erinnere dich daran

$USV = A$

$US = AV^{-1} = AV^T$$V$

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

$S^{-1}$

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

$D$$S$$\sqrt{D}$$U^TU=Identity$$U$$S$$V$$USV = A$

Die Berechnung der Diagonalisierungsdrehung kann durch Lösen der folgenden Gleichung erfolgen:

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

woher

$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);
}

Ideen von:
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/

Ich habe die Beschreibung unter http://www.lucidarme.me/?p=4624 verwendet , um diesen C ++ - Code zu erstellen. Die Matrizen sind die der Eigen-Bibliothek, aber Sie können aus diesem Beispiel leicht Ihre eigene Datenstruktur erstellen:

$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);

Mit der Standardzeichenfunktion

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

Daraus ergeben sich genau dieselben Werte wie für Eigen::JacobiSVD(siehe https://eigen.tuxfamily.org/dox-devel/classEigen_1_1JacobiSVD.html ).

Ich habe reinen C - Code für die 2x2 echte SVD hier . Siehe Zeile 559. Sie berechnet im Wesentlichen die Eigenwerte von$A^TA$durch das Lösen eines Quadrats ist es nicht unbedingt das robusteste, aber es scheint in der Praxis für nicht zu pathologische Fälle gut zu funktionieren. Es ist relativ einfach.

Aus persönlichen Gründen habe ich versucht, die Mindestberechnung für eine 2x2-DVD zu isolieren. Ich denke, es ist wahrscheinlich eine der einfachsten und schnellsten Lösungen. Details finden Sie in meinem persönlichen Blog: http://lucidarme.me/?p=4624 .

Vorteile: einfach, schnell und Sie können nur eine oder zwei der drei Matrizen (S, U oder D) berechnen, wenn Sie die drei Matrizen nicht benötigen.

Nachteil: Es wird atan2 verwendet, was möglicherweise ungenau ist und eine externe Bibliothek erfordert (typ. Math.h).

Hier ist eine Implementierung eines 2x2 SVD zu lösen. Ich habe es auf Victor Lius Code aufgebaut. Sein Code funktionierte für einige Matrizen nicht. Ich habe diese beiden Dokumente als mathematische Referenz für die Lösung verwendet: pdf1 und pdf2 .

Die Matrixmethode setDataist in der Hauptreihenfolge. Intern repräsentiere ich die Matrixdaten als 2D-Array von 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
            );
        }
    }
}