Робастный алгоритм для

Что такое простой алгоритм для вычисления SVD матрицы?$2 \times 2$

В идеале мне нужен численно устойчивый алгоритм, но я бы хотел увидеть как простые, так и не очень простые реализации. Код C принят.

Любые ссылки на документы или код?

См. Https://math.stackexchange.com/questions/861674/decompose-a-2d-arbitrary-transform-into-only-scaling-and-rotation (извините, я бы добавил это в комментарий, но я зарегистрировался просто чтобы опубликовать это, я пока не могу оставлять комментарии).

Но так как я пишу это как ответ, я также напишу метод:

Это разбивает матрицу следующим образом:

Единственное, на что следует обратить внимание при использовании этого метода, это то, что или для atan2.$G=F=0$$H=E=0$Я сомневаюсь, что это может быть более надежным, чем это( Обновление: см. Ответ Алекса Эфтимиадеса!).

Ссылка: http://dx.doi.org/10.1109/38.486688 (предоставлена ​​там Рахулом), которая находится в нижней части этого сообщения в блоге: http://metamerist.blogspot.com/2006/10/linear-algebra -дль-график-вундеркинды-svd.html

Обновление: Как отметил @VictorLiu в комментарии, может быть отрицательным. Это происходит тогда и только тогда, когда определитель входной матрицы также отрицателен. Если это так, и вы хотите положительные значения единственного числа, просто возьмите абсолютное значение .$s_y$$s_y$

@Pedro Gimeno

«Я сомневаюсь, что это может быть более надежным, чем это».

Вызов принят.

Я заметил, что обычный подход заключается в использовании тригонометрических функций, таких как atan2. Интуитивно понятно, что не должно быть необходимости использовать функции триггера. Действительно, все результаты заканчиваются синусами и косинусами арктанов, которые можно упростить до алгебраических функций. Это заняло много времени, но мне удалось упростить алгоритм Педро, чтобы использовать только алгебраические функции.

Следующий код Python делает свое дело.

от numy import asarray, diag


def svd2 (м):

    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 имеет СВД 2-на-2 решателя , лежащая в основе часть QR - разложения основного алгоритма SVD для gsl_linalg_SV_decomp. Смотрите svdstep.cфайл и ищите svd2функцию. У функции есть несколько особых случаев, она не совсем тривиальна, и, похоже, она делает несколько вещей, чтобы быть численно осторожной (например, использовать, hypotчтобы избежать переполнений).

Когда мы говорим «численно устойчивый», мы обычно имеем в виду алгоритм, в котором мы делаем такие вещи, как поворот, чтобы избежать распространения ошибок. Однако для матрицы 2x2 вы можете записать результат в терминах явных формул - то есть записать формулы для элементов SVD, которые указывают результат только в терминах входных данных , а не в терминах ранее вычисленных промежуточных значений . Это означает, что у вас может быть отмена, но нет распространения ошибки.

Дело просто в том, что для систем 2x2 беспокоиться о надежности не нужно.

Этот код основан на бумаге Blinn в , Эллис бумаги , СВД лекции и дополнительных расчетов. Алгоритм подходит для регулярных и сингулярных вещественных матриц. Все предыдущие версии работают на 100% так же, как и эта.

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

Мне нужен алгоритм, который имеет

• немного ветвления (надеюсь, CMOV)
• нет вызовов тригонометрических функций
• высокая точность вычислений даже с 32-разрядными числами

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

Напомним, что

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

Вычисление диагонализационного вращения может быть выполнено путем решения следующего уравнения:

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

где

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

Идеи из:
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/

Я использовал описание на http://www.lucidarme.me/?p=4624 для создания этого кода C ++. Матрицы принадлежат библиотеке Eigen, но вы можете легко создать свою собственную структуру данных из этого примера:

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

Со стандартной функцией знака

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

Это приводит к точно таким же значениям, что и Eigen::JacobiSVD(см. Https://eigen.tuxfamily.org/dox-devel/classEigen_1_1JacobiSVD.html ).

$A^TA$

Для моих личных нужд я попытался выделить минимальные вычисления для 2x2 SVD. Я думаю, что это, вероятно, одно из самых простых и быстрых решений. Вы можете найти подробности в моем личном блоге: http://lucidarme.me/?p=4624 .

Преимущества: просто, быстро, и вы можете рассчитать одну или две из трех матриц (S, U или D), если вам не нужны эти три матрицы.

Недостатком является использование atan2, которое может быть неточным и может потребовать внешней библиотеки (тип. Math.h).

Вот реализация решения 2x2 SVD. Я основал это на коде Виктора Лю. Его код не работал для некоторых матриц. Я использовал эти два документа в качестве математического справочника для решения: pdf1 и pdf2 .

Матричный setDataметод находится в главном порядке строк. Внутренне я представляю данные матрицы в виде двумерного массива, заданного выражением 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
            );
        }
    }
}