Algorithme robuste pour SVD

Qu'est-ce qu'un algorithme simple pour calculer la SVD de matrices?$2 \times 2$

Idéalement, j'aimerais un algorithme numériquement robuste, mais j'aimerais voir à la fois des implémentations simples et pas si simples. Code C accepté.

Des références à des articles ou à du code?

Voir /math/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (désolé, j'aurais mis cela dans un commentaire mais je me suis inscrit juste pour poster ceci donc je ne peux pas encore poster de commentaires).

Mais comme je l'écris comme réponse, j'écrirai également la méthode:

Cela décompose la matrice comme suit:

La seule chose à éviter avec cette méthode est que $G=F=0$ ou $H=E=0$ pour atan2.Je doute qu'il puisse être plus robuste que ça( Mise à jour: voir la réponse d'Alex Eftimiades!).

La référence est: http://dx.doi.org/10.1109/38.486688 (donnée par Rahul là-bas) qui vient du bas de cet article de blog: http://metamerist.blogspot.com/2006/10/linear-algebra -for-graphics-geeks-svd.html

Mise à jour: Comme indiqué par @VictorLiu dans un commentaire, $s_y$ peut être négatif. Cela se produit si et seulement si le déterminant de la matrice d'entrée est également négatif. Si c'est le cas et que vous voulez les valeurs singulières positives, prenez simplement la valeur absolue de $s_y$ .

@Pedro Gimeno

"Je doute qu'il puisse être plus robuste que cela."

Défi accepté.

J'ai remarqué que l'approche habituelle consiste à utiliser des fonctions trigonométriques comme atan2. Intuitivement, il ne devrait pas être nécessaire d'utiliser des fonctions trigonométriques. En effet, tous les résultats se retrouvent sous forme de sinus et cosinus d'arctans - ce qui peut être simplifié en fonctions algébriques. Cela a pris un certain temps, mais j'ai réussi à simplifier l'algorithme de Pedro pour utiliser uniquement les fonctions algébriques.

Le code python suivant fait l'affaire.

depuis 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

Le GSL a un solveur SVD 2 par 2 sous-jacent à la partie de décomposition QR de l'algorithme SVD principal pour gsl_linalg_SV_decomp. Consultez le svdstep.cfichier et recherchez la svd2fonction. La fonction a quelques cas particuliers, n'est pas vraiment triviale et semble faire plusieurs choses pour être prudent numériquement (par exemple, utiliser hypotpour éviter les débordements).

Lorsque nous disons "numériquement robuste", nous entendons généralement un algorithme dans lequel nous faisons des choses comme le pivotement pour éviter la propagation des erreurs. Cependant, pour une matrice 2x2, vous pouvez écrire le résultat en termes de formules explicites - c'est-à-dire, écrire des formules pour les éléments SVD qui énoncent le résultat uniquement en termes d'entrées , plutôt qu'en termes de valeurs intermédiaires précédemment calculées . Cela signifie que vous pouvez avoir une annulation mais pas de propagation d'erreur.

Le fait est simplement que pour les systèmes 2x2, il n'est pas nécessaire de se soucier de la robustesse.

Ce code est basé sur le papier de Blinn , papier Ellis , conférence SVD , et d' autres calculs. Un algorithme convient aux matrices réelles régulières et singulières. Toutes les versions précédentes fonctionnent à 100% ainsi que celle-ci.

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

J'avais besoin d'un algorithme qui a

• peu de ramification (espérons-le CMOV)
• aucun appel de fonction trigonométrique
• haute précision numérique même avec des flotteurs 32 bits

Nous voulons calculer $c_1, s_1, c_2, s_2, \sigma_1$ et $\sigma_2$ comme suit:

$A = USV$, qui peut être développé comme:

$\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}$ ou $-\sqrt{D}$), que j'ai corrigé dans le code ci-dessous.

J'ai jeté environ 2000 millions de matrices randomisées sur le code, et la plus grande erreur numérique produite était d'environ $6\cdot10^{-7}$ (avec flotteurs 32 bits, $error = ||USV-M||/||M||$). L'algorithme s'exécute sur environ 340 cycles d'horloge (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);
}

Idées de:
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/

J'ai utilisé la description à http://www.lucidarme.me/?p=4624 pour créer ce code C ++. Les matrices sont celles de la bibliothèque Eigen, mais vous pouvez facilement créer votre propre structure de données à partir de cet exemple:

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

Avec la fonction de signe standard

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

Il en résulte exactement les mêmes valeurs que le Eigen::JacobiSVD(voir https://eigen.tuxfamily.org/dox-devel/classEigen_1_1JacobiSVD.html ).

J'ai du code C pur pour le vrai SVD 2x2 ici . Voir ligne 559. Il calcule essentiellement les valeurs propres de$A^TA$en résolvant un quadratique, ce n'est donc pas forcément le plus robuste, mais il semble bien fonctionner en pratique pour les cas pas trop pathologiques. C'est relativement simple.

Pour mes besoins personnels, j'ai essayé d'isoler le calcul minimum pour un svd 2x2. Je suppose que c'est probablement l'une des solutions les plus simples et les plus rapides. Vous pouvez trouver des détails sur mon blog personnel: http://lucidarme.me/?p=4624 .

Avantages: simple, rapide et vous ne pouvez calculer qu'une ou deux des trois matrices (S, U ou D) si vous n'avez pas besoin des trois matrices.

Inconvénient, il utilise atan2, qui peut être inexact et peut nécessiter une bibliothèque externe (typ. Math.h).

Voici une implémentation d'une résolution SVD 2x2. Je l'ai basé sur le code de Victor Liu. Son code ne fonctionnait pas pour certaines matrices. J'ai utilisé ces deux documents comme référence mathématique pour la résolution: pdf1 et pdf2 .

La setDataméthode matricielle est dans l'ordre des lignes principales. En interne, je représente les données de la matrice comme un tableau 2D donné par 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
            );
        }
    }
}