Algoritma yang kuat untuk SVD

26

Apa itu algoritma sederhana untuk menghitung SVD matriks?2×2

Idealnya, saya ingin algoritma yang kuat secara numerik, tetapi saya ingin melihat implementasi yang sederhana dan yang tidak terlalu sederhana. Kode C diterima.

Adakah referensi untuk makalah atau kode?

lhf
sumber
5
Wikipedia mencantumkan solusi formulir tertutup 2x2, tetapi saya tidak tahu sifat numeriknya.
Damien
Sebagai referensi, "Resep Numerik", Press et al., Cambridge Press. Buku yang cukup mahal tetapi bernilai setiap sen. Selain solusi SVD, Anda akan menemukan banyak algoritma bermanfaat lainnya.
Jan Hackenberg

Jawaban:

19

Lihat /math/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (maaf, saya akan menuliskannya di komentar tapi saya sudah mendaftar hanya untuk memposting ini jadi saya belum bisa memposting komentar).

Tetapi karena saya menulisnya sebagai jawaban, saya juga akan menulis metode:

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

Itu menguraikan matriks sebagai berikut:

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

Satu-satunya hal yang harus dijaga dengan metode ini adalah bahwa atau untuk atan2.G=F=0H=E=0Saya ragu itu bisa lebih kuat dari itu( Perbarui: lihat jawaban Alex Eftimiades!).

Referensi tersebut adalah: http://dx.doi.org/10.1109/38.486688 (diberikan oleh Rahul di sana) yang berasal dari bagian bawah posting blog ini: http://metamerist.blogspot.com/2006/10/linear-algebra -for-graphics-geeks-svd.html

Pembaruan: Sebagaimana dicatat oleh @VictorLiu dalam komentar, mungkin negatif. Itu terjadi jika dan hanya jika penentu matriks input negatif juga. Jika itu masalahnya dan Anda ingin nilai singular positif, ambil saja nilai absolut dari .sysy

Pedro Gimeno
sumber
1
Tampaknya dapat menjadi negatif jika . Ini seharusnya tidak mungkin. Q < RsyQ<R
Victor Liu
@ VictorLiu Jika matriks input membalik, satu-satunya tempat yang dapat tercermin dalam matriks penskalaan, karena matriks rotasi tidak mungkin flip. Hanya saja, jangan mengumpankannya masukan matriks yang flip. Saya belum melakukan matematika tetapi saya bertaruh bahwa tanda penentu matriks input akan menentukan apakah atau lebih besar. RQR
Pedro Gimeno
@ ViktorLiu Saya sudah melakukan matematika sekarang dan mengkonfirmasi bahwa memang, menyederhanakan menjadi yaitu penentu matriks input. m 00 m 11 - m 01 m 10Q2R2m00m11m01m10
Pedro Gimeno
9

@Pedro Gimeno

"Aku ragu itu bisa lebih kuat dari itu."

Tantangan diterima.

Saya perhatikan pendekatan yang biasa digunakan adalah menggunakan fungsi trigonometri seperti atan2. Secara intuitif, seharusnya tidak perlu menggunakan fungsi trigonometri. Memang, semua hasil berakhir sebagai sinus dan cosinus arctans - yang dapat disederhanakan menjadi fungsi aljabar. Butuh waktu cukup lama, tetapi saya berhasil menyederhanakan algoritma Pedro untuk hanya menggunakan fungsi aljabar.

Kode python berikut melakukan triknya.

dari asarray impor numpy, 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
sumber
1
Kode sepertinya salah. Pertimbangkan matriks identitas 2x2. Kemudian y1= 0, x1= 0, h1= 0, dan t1= 0/0 = NaN.
Pelukan
8

The GSL memiliki SVD 2-by-2 solver yang mendasari bagian QR dekomposisi algoritma SVD utama untuk gsl_linalg_SV_decomp. Lihat svdstep.cfile dan cari svd2fungsinya. Fungsi ini memiliki beberapa kasus khusus, tidak sepele, dan terlihat melakukan beberapa hal untuk berhati-hati secara numerik (misalnya, menggunakan hypotuntuk menghindari luapan).

horchler
sumber
1
Apakah fungsi ini memiliki dokumentasi? Saya ingin tahu apa parameter inputnya.
Victor Liu
@ ViktorLiu: Sayangnya saya belum melihat apa pun selain komentar yang sedikit dalam file itu sendiri. Ada sedikit ChangeLogfile jika Anda mengunduh GSL. Dan Anda dapat melihat svd.cdetail dari keseluruhan algoritma. Satu-satunya dokumentasi yang benar tampaknya untuk fungsi-fungsi yang dapat dipanggil pengguna tingkat tinggi, misalnya gsl_linalg_SV_decomp,.
horchler
7

Ketika kita mengatakan "kuat secara numerik" biasanya kita maksudkan algoritma di mana kita melakukan hal-hal seperti berputar untuk menghindari penyebaran kesalahan. Namun, untuk matriks 2x2, Anda dapat menuliskan hasilnya dalam bentuk rumus eksplisit - yaitu, menuliskan rumus untuk elemen SVD yang menyatakan hasil hanya dalam hal input , bukan dalam hal nilai menengah yang sebelumnya dihitung . Itu berarti bahwa Anda mungkin memiliki pembatalan tetapi tidak ada propagasi kesalahan.

Intinya adalah bahwa untuk sistem 2x2, tidak perlu khawatir tentang ketahanan.

Wolfgang Bangerth
sumber
Itu bisa bergantung pada matriks. Saya telah melihat metode yang menemukan sudut kiri dan kanan secara terpisah (masing-masing melalui arctan2 (y, x)) yang umumnya berfungsi dengan baik. Tetapi ketika nilai singular berdekatan, masing-masing arctans ini cenderung 0/0, sehingga hasilnya bisa tidak akurat. Dalam metode yang diberikan oleh Pedro Gimeno, perhitungan a2 akan didefinisikan dengan baik dalam kasus ini, sementara a1 menjadi tidak jelas; Anda masih memiliki hasil yang baik karena validitas dekomposisi hanya sensitif terhadap theta + phi ketika s.vals berdekatan, bukan theta-phi.
greggo
5

Kode ini didasarkan pada kertas Blinn ini , kertas Ellis , SVD kuliah , dan tambahan perhitungan. Algoritma cocok untuk matriks nyata reguler dan tunggal. Semua versi sebelumnya berfungsi 100% dan juga yang ini.

#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
sumber
5

Saya membutuhkan algoritma yang dimiliki

  • sedikit percabangan (semoga CMOVs)
  • tidak ada panggilan fungsi trigonometri
  • akurasi numerik tinggi bahkan dengan floats 32 bit

Kami ingin menghitung dan sebagai berikut:c1,s1,c2,s2,σ1σ2

A=USV , yang dapat diperluas seperti:

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

Gagasan utamanya adalah menemukan matriks rotasi yang mendiagonalisasi , yaitu adalah diagonal.VATAVATAVT=D

Ingat itu

USV=A

US=AV1=AVT (karena adalah orthogonal)V

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

Mengalikan kedua sisi dengan kita dapatkanS1

(STST)UTU(SS1)=UTU=STDS1

Karena adalah diagonal, pengaturan ke akan memberi kita , artinya adalah matriks rotasi, adalah matriks diagonal, adalah matriks rotasi dan , seperti yang kita cari untuk.DSDUTU=IdentityUSVUSV=A

Menghitung rotasi diagonalisasi dapat dilakukan dengan menyelesaikan persamaan berikut:

t22βαγt21=0

dimana

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

dan adalah tangen dari sudut . Ini dapat diturunkan dengan memperluas dan membuat elemen off-diagonal sama dengan nol (mereka sama satu sama lain).t2VVATAVT

Masalah dengan metode ini adalah bahwa ia kehilangan presisi floating point yang signifikan ketika menghitung dan untuk matriks tertentu, karena pengurangan dalam perhitungan. Solusi untuk ini adalah untuk melakukan dekomposisi RQ ( , segitiga atas dan orthogonal) pertama, kemudian menggunakan algoritma untuk pd . Hal ini memberikan . Perhatikan bagaimana pengaturan to 0 (seperti pada ) menghilangkan beberapa penambahan / pengurangan. (Dekomposisi RQ cukup sepele dari perluasan produk matriks).βαγA=RQRQUSV=RUSV=USVQ=RQ=AdR

Algoritme yang diimplementasikan secara naif dengan cara ini memiliki beberapa anomali numerik dan logis (misalnya atau ), yang saya perbaiki dalam kode di bawah ini.S +DD

Saya melemparkan sekitar 2000 juta matriks acak pada kode, dan kesalahan numerik terbesar yang dihasilkan adalah sekitar (dengan floats 32 bit, ). Algoritma ini berjalan di sekitar 340 siklus clock (MSVC 19, Ivy Bridge).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);
}

Gagasan dari:
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
sumber
3

Saya telah menggunakan deskripsi di http://www.lucidarme.me/?p=4624 untuk membuat kode C ++ ini. Matriksnya adalah dari perpustakaan Eigen, tetapi Anda dapat dengan mudah membuat struktur data Anda sendiri dari contoh ini:

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

Dengan fungsi tanda standar

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

Ini menghasilkan nilai yang persis sama dengan Eigen::JacobiSVD(lihat https://eigen.tuxfamily.org/dox-devel/classEigen_1_1JacobiSVD.html ).

Corbie
sumber
1
S2 = hypot( a*a + b*b - c*c - d*d, 2*(a*c + b*d))
greggo
2

I have pure C code for the 2x2 real SVD here. See line 559. It essentially computes the eigenvalues of ATA 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.

Victor Liu
sumber
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...
Carlos Scheidegger
2

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

Fifi
sumber
3
Since links are rarely permanent, it is important to summarize the approach rather than simply providing a link as an answer.
Paul
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.
Geoff Oxberry
1

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