Algoritmo robusto para SVD

26

¿Qué es un algoritmo simple para calcular la SVD de matrices ?2×2

Idealmente, me gustaría un algoritmo numéricamente robusto, pero me gustaría ver implementaciones simples y no tan simples. Código C aceptado.

¿Alguna referencia a documentos o código?

lhf
fuente
55
Wikipedia enumera una solución de forma cerrada de 2x2, pero no tengo idea de sus propiedades numéricas.
Damien
Como referencia, "Numerical Recipes", Press et al., Cambridge Press. Libro bastante caro pero vale cada centavo. Además de las soluciones SVD, encontrará muchos otros algoritmos útiles.
Jan Hackenberg

Respuestas:

19

Ver /math/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (lo siento, lo habría puesto en un comentario pero me he registrado solo para publicar esto para que no pueda publicar comentarios todavía).

Pero como lo escribo como respuesta, también escribiré el método:

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

Eso descompone la matriz de la siguiente manera:

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

Lo único contra lo que debe protegerse con este método es que o para atan2.G=F=0H=E=0Dudo que pueda ser más robusto que eso( Actualización: ¡ vea la respuesta de Alex Eftimiades!).

La referencia es: http://dx.doi.org/10.1109/38.486688 (dada por Rahul allí) que viene del final de esta publicación de blog: http://metamerist.blogspot.com/2006/10/linear-algebra -for-graphics-geeks-svd.html

Actualización: como señaló @VictorLiu en un comentario, puede ser negativo. Eso sucede si y solo si el determinante de la matriz de entrada también es negativo. Si ese es el caso y desea los valores singulares positivos, simplemente tome el valor absoluto de .sysy

Pedro Gimeno
fuente
1
Parece que puede ser negativo si . Esto no debería ser posible. Q < RsyQ<R
Victor Liu
@VictorLiu Si la matriz de entrada se voltea, el único lugar que puede reflejarse es en la matriz de escala, ya que las matrices de rotación no pueden voltearse. Simplemente no lo alimente con matrices de entrada que se invierten. Todavía no he hecho los cálculos, pero apuesto a que el signo del determinante de la matriz de entrada determinará si o es mayor. RQR
Pedro Gimeno
@VictorLiu Ya he hecho los cálculos y he confirmado que, de hecho, simplifica a es decir, el determinante de la matriz de entrada. m 00 m 11 - m 01 m 10Q2R2m00m11m01m10
Pedro Gimeno
9

@Pedro Gimeno

"Dudo que pueda ser más robusto que eso".

Desafío aceptado.

Noté que el enfoque habitual es utilizar funciones trigonométricas como atan2. Intuitivamente, no debería ser necesario usar funciones trigonométricas. De hecho, todos los resultados terminan como senos y cosenos de arctanos, que pueden simplificarse a funciones algebraicas. Me llevó bastante tiempo, pero logré simplificar el algoritmo de Pedro para usar solo funciones algebraicas.

El siguiente código de Python hace el truco.

de 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
fuente
1
El código parece incorrecto. Considere la matriz de identidad 2x2. Entonces y1= 0, x1= 0, h1= 0, y t1= 0/0 = NaN.
Hugues
8

El GSL tiene un solucionador SVD de 2 por 2 subyacente a la parte de descomposición QR del algoritmo SVD principal para gsl_linalg_SV_decomp. Vea el svdstep.carchivo y busque la svd2función. La función tiene algunos casos especiales, no es exactamente trivial y parece estar haciendo varias cosas para ser numéricamente cuidadoso (por ejemplo, usar hypotpara evitar desbordamientos).

horchler
fuente
1
¿Esta función tiene alguna documentación? Me gustaría saber cuáles son sus parámetros de entrada.
Victor Liu
@VictorLiu: Lamentablemente, no he visto nada más que los escasos comentarios en el archivo en sí. Hay un poco en el ChangeLogarchivo si descarga el GSL. Y puede consultar los svd.cdetalles del algoritmo general. La documentación sólo es cierto parece ser para las funciones se puede llamar de usuario de alto nivel, por ejemplo, gsl_linalg_SV_decomp.
horchler
7

Cuando decimos "numéricamente robusto" usualmente nos referimos a un algoritmo en el que hacemos cosas como pivotar para evitar la propagación de errores. Sin embargo, para una matriz de 2x2, puede escribir el resultado en términos de fórmulas explícitas, es decir, escribir fórmulas para los elementos SVD que expresan el resultado solo en términos de las entradas , en lugar de en términos de valores intermedios previamente calculados . Eso significa que puede tener cancelación pero no propagación de errores.

El punto simplemente es que para sistemas 2x2, no es necesario preocuparse por la robustez.

Wolfgang Bangerth
fuente
Puede depender de la matriz. He visto un método que encuentra los ángulos izquierdo y derecho por separado (cada uno a través de arctan2 (y, x)) que generalmente funciona bien. Pero cuando los valores singulares están muy juntos, cada uno de estos arctanos tiende a 0/0, por lo que el resultado puede ser inexacto. En el método proporcionado por Pedro Gimeno, el cálculo de a2 estará bien definido en este caso, mientras que a1 estará mal definido; todavía tiene un buen resultado porque la validez de la descomposición solo es sensible a theta + phi cuando los valores s están muy juntos, no a theta-phi.
Greg
5

Este código se basa en el artículo de Blinn , papel Ellis , conferencia SVD , y otros cálculos. Un algoritmo es adecuado para matrices reales regulares y singulares. Todas las versiones anteriores funcionan al 100% tan bien como esta.

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

Necesitaba un algoritmo que tenga

  • poca ramificación (con suerte CMOV)
  • sin llamadas a funciones trigonométricas
  • alta precisión numérica incluso con flotadores de 32 bits

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

A=USV

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

VATAVATAVT=D

Recordar que

USV=A

US=AV1=AVTV

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

S1

(STST)UTU(SS1)=UTU=STDS1

DSDUTU=IdentityUSVUSV=A

Se puede calcular la rotación de diagonalización resolviendo la siguiente ecuación:

t22βαγt21=0

dónde

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

Ideas 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/

petiaccja
fuente
3

He usado la descripción en http://www.lucidarme.me/?p=4624 para crear este código C ++. Las matrices son las de la biblioteca Eigen, pero puede crear fácilmente su propia estructura de datos a partir de este ejemplo:

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

Con la función de signo estándar

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

Esto da como resultado exactamente los mismos valores que el Eigen::JacobiSVD(consulte https://eigen.tuxfamily.org/dox-devel/classEigen_1_1JacobiSVD.html ).

Grajo negro
fuente
1
S2 = hypot( a*a + b*b - c*c - d*d, 2*(a*c + b*d))
Greggo
2

ATA

Victor Liu
fuente
No creo que su código funcione cuando los valores propios de la matriz son negativos. Pruebe [[1 1] [1 0]], y u * s * vt no es igual a m ...
Carlos Scheidegger
2

Para mi necesidad personal, traté de aislar el cálculo mínimo para un svd de 2x2. Supongo que es probablemente una de las soluciones más simples y rápidas. Puede encontrar detalles en mi blog personal: http://lucidarme.me/?p=4624 .

Ventajas: simple, rápido y solo puede calcular una o dos de las tres matrices (S, U o D) si no necesita las tres matrices.

La desventaja es que utiliza atan2, que puede ser impreciso y puede requerir una biblioteca externa (típ. Math.h).

Fifi
fuente
3
Como los enlaces rara vez son permanentes, es importante resumir el enfoque en lugar de simplemente proporcionar un enlace como respuesta.
Paul
Además, si va a publicar un enlace a su propio blog, (a) revele que es su blog, (b) aún mejor sería resumir o cortar y pegar su enfoque (las imágenes de las fórmulas pueden se traducirá a LaTeX sin formato y se representará con MathJax). Las mejores respuestas para este tipo de preguntas formulan fórmulas, proporcionan citas para dichas fórmulas y luego enumeran cosas como inconvenientes, casos extremos y posibles alternativas.
Geoff Oxberry
1

Aquí hay una implementación de una resolución SVD 2x2. Lo basé en el código de Victor Liu. Su código no funcionaba para algunas matrices. Usé estos dos documentos como referencia matemática para la resolución: pdf1 y pdf2 .

El setDatamétodo de la matriz está en orden de fila mayor. Internamente, represento los datos de la matriz como una matriz 2D dada por 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
fuente