Solidny algorytm dla

Co to jest prosty algorytm obliczania SVD macierzy ?$2 \times 2$

Idealnie chciałbym mieć solidny algorytm liczbowy, ale chciałbym zobaczyć zarówno proste, jak i nie tak proste implementacje. Kod C został zaakceptowany.

Wszelkie odniesienia do artykułów lub kodu?

Zobacz /math/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (przepraszam, umieściłbym to w komentarzu, ale zarejestrowałem się po prostu opublikować, żeby nie dodawać jeszcze komentarzy).

Ale ponieważ piszę to jako odpowiedź, napiszę również metodę:

To rozkłada macierz w następujący sposób:

Jedyną rzeczą, której należy się wystrzegać za pomocą tej metody jest to, że $G=F=0$ lub $H=E=0$ dla atan2.Wątpię, czy może być bardziej solidny( Aktualizacja: patrz odpowiedź Alexa Eftimiadesa!).

Odnośnik to: http://dx.doi.org/10.1109/38.486688 (podany tam przez Rahula), który pochodzi z dolnej części tego postu na blogu: http://metamerist.blogspot.com/2006/10/linear-algebra -for-graphics-geeks-svd.html

Aktualizacja: Jak zauważył @VictorLiu w komentarzu, $s_y$ może być ujemny. Dzieje się tak wtedy i tylko wtedy, gdy wyznacznik macierzy wejściowej jest również ujemny. Jeśli tak jest w przypadku i chcesz dodatnich wartości pojedynczych, po prostu weź bezwzględną wartość $s_y$ .

@Pedro Gimeno

„Wątpię, czy może być bardziej solidny.”

Wyzwanie przyjęte.

Zauważyłem, że typowym podejściem jest używanie funkcji trig, takich jak atan2. Intuicyjnie nie powinno być potrzeby używania funkcji trig. Rzeczywiście, wszystkie wyniki kończą się jako sinus i cosinus arctans - które można uprościć do funkcji algebraicznych. Trwało to dość długo, ale udało mi się uprościć algorytm Pedro, aby używał tylko funkcji algebraicznych.

Poniższy kod python załatwia sprawę.

z 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

GSL ma SVD 2-by-2 solver którym opiera się część rozkładu QR głównego algorytmu dla SVD gsl_linalg_SV_decomp. Zobacz svdstep.cplik i wyszukaj svd2funkcję. Funkcja ma kilka specjalnych przypadków, nie jest trywialna i wygląda na to, że robi kilka rzeczy, aby zachować ostrożność liczbową (np. Używać, hypotaby uniknąć przepełnienia).

Kiedy mówimy „odporny numerycznie”, zwykle mamy na myśli algorytm, w którym robimy takie rzeczy, jak obracanie, aby uniknąć propagacji błędów. Jednak w przypadku macierzy 2x2 można zapisać wynik w postaci jawnych wzorów - tj. Zapisać wzory dla elementów SVD, które podają wynik tylko w kategoriach danych wejściowych , a nie w odniesieniu do wcześniej obliczonych wartości pośrednich . Oznacza to, że możesz mieć anulowanie, ale brak propagacji błędów.

Chodzi o to, że w przypadku systemów 2x2 martwienie się o solidność nie jest konieczne.

Kod ten jest oparty na papierze Blinn jest , Ellis papieru , wykład SVD i dodatkowych obliczeń. Algorytm nadaje się do regularnych i pojedynczych rzeczywistych macierzy. Wszystkie poprzednie wersje działają w 100% tak dobrze, jak ta.

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

Potrzebowałem algorytmu, który ma

• małe rozgałęzienia (mam nadzieję, że CMOV)
• brak wywołań funkcji trygonometrycznych
• wysoka dokładność numeryczna nawet przy 32-bitowych liczbach zmiennoprzecinkowych

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

Odwołaj to

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

Obliczanie obrotu diagonalizującego można wykonać, rozwiązując następujące równanie:

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

gdzie

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

$t_2$$V$$VA^TAV^T$

$\beta-\alpha$$\gamma$$A=RQ$$R$$Q$$USV' = R$$USV=USV'Q=RQ=A$$d$$R$

$S$ $+\sqrt{D}$$-\sqrt{D}$

$6\cdot10^{-7}$$error = ||USV-M||/||M||$

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

Pomysły:
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/

Użyłem opisu pod adresem http://www.lucidarme.me/?p=4624 do stworzenia tego kodu C ++. Matryce pochodzą z biblioteki Eigen, ale z tego przykładu można łatwo utworzyć własną strukturę danych:

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

Ze standardową funkcją znaku

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

Daje to dokładnie takie same wartości jak Eigen::JacobiSVD(patrz https://eigen.tuxfamily.org/dox-devel/classEigen_1_1JacobiSVD.html ).

$A^TA$

Na własne potrzeby starałem się wyodrębnić minimalne obliczenia dla dysku DVD 2x2. Myślę, że jest to prawdopodobnie jedno z najprostszych i najszybszych rozwiązań. Szczegóły znajdziesz na moim osobistym blogu: http://lucidarme.me/?p=4624 .

Zalety: proste, szybkie i można obliczyć jedną lub dwie z trzech macierzy (S, U lub D) tylko wtedy, gdy nie potrzebujesz trzech macierzy.

Wadą jest atan2, który może być niedokładny i może wymagać biblioteki zewnętrznej (typ. Math.h).

Oto implementacja rozwiązania SVD 2x2. Oparłem go na kodzie Victora Liu. Jego kod nie działał dla niektórych matryc. Użyłem tych dwóch dokumentów jako matematycznego odniesienia do rozwiązania: pdf1 i pdf2 .

Metoda macierzowa setDatajest uporządkowana według rzędów. Wewnętrznie reprezentuję dane macierzy jako tablicę 2D podaną przez 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
            );
        }
    }
}