Solidny algorytm dla

26

Co to jest prosty algorytm obliczania SVD macierzy ?2×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?

lhf
źródło
5
Wikipedia wymienia rozwiązanie w formie zamkniętej 2x2, ale nie mam pojęcia o jego liczbowych właściwościach.
Damien
Jako odniesienie, „Numerical Recipes”, Press i in., Cambridge Press. Książka dość droga, ale warta każdego centa. Oprócz rozwiązań SVD znajdziesz wiele innych przydatnych algorytmów.
Jan Hackenberg,

Odpowiedzi:

19

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

E=m00+m112;F=m00m112;G=m10+m012;H=m10m012Q=E2+H2;R=F2+G2sx=Q+R;sy=QRa1=atan2(G,F);a2=atan2(H,E)θ=a2a12;ϕ=a2+a12

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

M=(m00m01m10m11)=(cosϕsinϕsinϕcosϕ)(sx00sy)(cosθsinθsinθcosθ)

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, sy 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ść sy .

Pedro Gimeno
źródło
1
Wydaje się, że może być ujemny, jeśli Q < R . To nie powinno być możliwe. syQ<R
Victor Liu
@VictorLiu Jeśli matryca wejściowa się odwraca, jedynym miejscem, które można odbić, jest matryca skalowania, ponieważ matryce obrotowe prawdopodobnie nie mogą się odwrócić. Tylko nie karm tego, że matryce wejściowe się odwracają. Jeszcze nie zrobiłem matematyki, ale założę się, że znak wyznacznika macierzy wejściowej określi, czy lub R jest większy. QR
Pedro Gimeno,
@VictorLiu Zrobiłem teraz matematykę i potwierdziłem, że rzeczywiście, upraszcza do m 00 m 11 - m 01 m 10, tj. Wyznacznik macierzy wejściowej. Q2R2m00m11m01m10
Pedro Gimeno,
9

@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

Alex Eftimiades
źródło
1
Kod wydaje się niepoprawny. Rozważ macierz tożsamości 2x2. Następnie y1= 0, x1= 0, h1= 0 i t1= 0/0 = NaN.
Hugues
8

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).

horchler
źródło
1
Czy ta funkcja ma jakąś dokumentację? Chciałbym wiedzieć, jakie są jego parametry wejściowe.
Victor Liu
@VictorLiu: Niestety nie widziałem nic poza skromnymi komentarzami w samym pliku. W ChangeLogpliku GSL jest trochę pliku. I możesz spojrzeć na svd.cszczegóły dotyczące całego algorytmu. Jedyną prawdziwą dokumentacją wydaje się być funkcja wysokiego poziomu, którą można wywoływać, np gsl_linalg_SV_decomp.
horchler
7

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.

Wolfgang Bangerth
źródło
Może zależeć od matrycy. Widziałem metodę, która odnajduje osobno lewy i prawy kąt (każdy poprzez arctan2 (y, x)), który ogólnie działa dobrze. Ale gdy pojedyncze wartości są blisko siebie, każdy z tych arktanów ma tendencję do 0/0, więc wynik może być niedokładny. W metodzie podanej przez Pedro Gimeno obliczenia a2 będą w tym przypadku dobrze określone, podczas gdy a1 będzie źle zdefiniowane; nadal masz dobry wynik, ponieważ ważność rozkładu jest wrażliwa na theta + phi tylko wtedy, gdy s.vals są blisko siebie, a nie na theta-phi.
greggo,
5

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]);
}
Martynas Sabaliauskas
źródło
5

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

c1,s1,c2,s2,σ1σ2

A=USV

[abcd]=[c1s1s1c1][σ100σ2][c2s2s2c2]

VATAVATAVT=D

Odwołaj to

USV=A

US=AV1=AVTV

VATAVT=(AVT)TAVT=(US)TUS=STUTUS=D

S1

(STST)UTU(SS1)=UTU=STDS1

DSDUTU=IdentityUSVUSV=A

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

t22βαγt21=0

gdzie

ATA=[acbd][abcd]=[a2+c2ab+cdab+cdb2+d2]=[αγγβ]

t2VVATAVT

βαγA=RQRQUSV=RUSV=USVQ=RQ=AdR

S +DD

6107error=||USVM||/||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/

petiaccja
źródło
3

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ΣVT

#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 ).

Corbie
źródło
1
S2 = hypot( a*a + b*b - c*c - d*d, 2*(a*c + b*d))
greggo,
2

ATA

Victor Liu
źródło
Nie sądzę, aby twój kod działał, gdy wartości własne macierzy są ujemne. Spróbuj [[1 1] [1 0]], a u * s * vt nie jest równe m ...
Carlos Scheidegger
2

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).

Fifi
źródło
3
Ponieważ linki rzadko są trwałe, ważne jest podsumowanie podejścia, a nie tylko podanie linku jako odpowiedzi.
Paweł
Ponadto, jeśli zamierzasz opublikować link do własnego bloga, (a) ujawnij, że to Twój blog, (b) jeszcze lepiej byłoby podsumować lub wyciąć i wkleić swoje podejście (obrazy formuł mogą zostać przetłumaczone na surowy LaTeX i renderowane przy użyciu MathJax). Najlepsze odpowiedzi dla tego rodzaju formuł stanów pytań, cytowanie tych formuł, a następnie lista rzeczy takich jak wady, przypadki krawędzi i potencjalne alternatywy.
Geoff Oxberry
1

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
            );
        }
    }
}
Azmisow
źródło