Algoritmo robusto para

O que é um algoritmo simples para calcular o SVD de matrizes?$2 \times 2$

Idealmente, eu gostaria de um algoritmo numericamente robusto, mas gostaria de ver implementações simples e não tão simples. Código C aceito.

Alguma referência a documentos ou código?

Consulte /math/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (desculpe, eu teria colocado isso em um comentário, mas me registrei apenas para postar isso, então eu não posso postar comentários ainda).

Mas, como estou escrevendo como resposta, também escreverei o método:

Isso decompõe a matriz da seguinte maneira:

A única coisa a ser evitada com esse método é que $G=F=0$ ou $H=E=0$ para atan2.Duvido que possa ser mais robusto que isso( Atualização: veja a resposta de Alex Eftimiades!).

A referência é: http://dx.doi.org/10.1109/38.486688 (fornecida por Rahul), que vem da parte inferior desta postagem do blog: http://metamerist.blogspot.com/2006/10/linear-algebra -para-gráficos-geeks-svd.html

Update: Como observado por @VictorLiu em um comentário, $s_y$ pode ser negativo. Isso acontece se, e somente se, o determinante da matriz de entrada também for negativo. Se for esse o caso e você desejar os valores singulares positivos, basta pegar o valor absoluto de $s_y$ .

@Pedro Gimeno

"Duvido que possa ser mais robusto que isso".

Desafio aceito.

Percebi que a abordagem usual é usar funções trigonométricas como atan2. Intuitivamente, não deve haver necessidade de usar funções trigonométricas. De fato, todos os resultados acabam em senos e cossenos de arctans - que podem ser simplificados para funções algébricas. Demorou um pouco, mas consegui simplificar o algoritmo de Pedro para usar apenas funções algébricas.

O seguinte código python faz o truque.

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

O GSL possui um solucionador SVD 2 por 2 subjacente à parte de decomposição QR do algoritmo SVD principal para gsl_linalg_SV_decomp. Veja o svdstep.carquivo e procure a svd2função. A função tem alguns casos especiais, não é exatamente trivial e parece estar fazendo várias coisas para ser numericamente cuidadosa (por exemplo, usar hypotpara evitar estouros).

Quando dizemos "numericamente robusto", geralmente queremos dizer um algoritmo no qual fazemos coisas como rotação para evitar a propagação de erros. No entanto, para uma matriz 2x2, é possível anotar o resultado em termos de fórmulas explícitas - ou seja, anotar fórmulas para os elementos SVD que indicam o resultado apenas em termos de entradas , e não em valores intermediários anteriormente calculados . Isso significa que você pode ter cancelamento, mas nenhuma propagação de erro.

O ponto simplesmente é que, para sistemas 2x2, não é necessário se preocupar com robustez.

Este código é baseado no papel de Blinn , papel Ellis , palestra SVD e adicionais cálculos. Um algoritmo é adequado para matrizes reais regulares e singulares. Todas as versões anteriores funcionam 100%, bem 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]);
}

Eu precisava de um algoritmo que tenha

• pouca ramificação (espero CMOVs)
• nenhuma chamada de função trigonométrica
• alta precisão numérica, mesmo com flutuadores de 32 bits

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

Lembre-se que

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

O cálculo da rotação da diagonalização pode ser feito resolvendo a seguinte equação:

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

Onde

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

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

Eu usei a descrição em http://www.lucidarme.me/?p=4624 para criar esse código C ++. As matrizes são as da biblioteca Eigen, mas você pode criar facilmente sua própria estrutura de dados a partir deste exemplo:

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

Com a função de sinal padrão

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

Isso resulta exatamente nos mesmos valores que o Eigen::JacobiSVD(consulte https://eigen.tuxfamily.org/dox-devel/classEigen_1_1JacobiSVD.html ).

$A^TA$

Para minha necessidade pessoal, tentei isolar o cálculo mínimo para um svd 2x2. Eu acho que é provavelmente uma das soluções mais simples e rápidas. Você pode encontrar detalhes no meu blog pessoal: http://lucidarme.me/?p=4624 .

Vantagens: simples, rápido e você só pode calcular uma ou duas das três matrizes (S, U ou D) se não precisar das três matrizes.

A desvantagem é o atan2, que pode ser impreciso e pode exigir uma biblioteca externa (tipo. Math.h).

Aqui está uma implementação de uma solução SVD 2x2. Baseei-o no código de Victor Liu. Seu código não estava funcionando para algumas matrizes. Usei esses dois documentos como referência matemática para a solução: pdf1 e pdf2 .

O setDatamétodo da matriz está na ordem das principais linhas. Internamente, represento os dados da matriz como uma 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
            );
        }
    }
}