SVD分解 opencv实现

头文件

#ifndef DEBUG_LRN_SVD_H
#define DEBUG_LRN_SVD_H
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
#include <opencv2/opencv.hpp>
#include </usr/local/include/eigen3/Eigen/Dense>

using namespace std;
using namespace cv;

int test_SVD();

int opencv_svd(vector<vector<float>>& vec);

int eigen_svd(vector<vector<float>>& vec);

int cpp_svd(vector<vector<float>>& vec);

void print_matrix(const vector<vector<float>>& vec);

#endif //DEBUG_LRN_SVD_H

源文件

#include "svd.h"

void print_matrix(const vector<vector<float>>& vec){
    if(vec.empty()){
        return;
    }
    for ( auto row: vec){
        if(row.empty()){
            return;
        }
        for ( auto elem: row){
            cout << elem << ", ";
        }
        cout << endl;
    }
    cout << endl;
}

int opencv_svd(vector<vector<float>>& vec){
    cout << "source matrix：" << endl;
    print_matrix(vec);

    cout << "opencv singular value decomposition:" << endl;
    const auto rows = (int)vec.size();
    const auto cols = (int)vec[0].size();
    cv::Mat mat(rows, cols, CV_32FC1);
    for (int y = 0; y < rows; ++y) {
        for (int x = 0; x < cols; ++x) {
            mat.at<float>(y, x) = vec.at(y).at(x);
        }
    }

    cv::Mat matD, matU, matVt;
    cv::SVD::compute(mat, matD, matU, matVt, 4);
    cout << "opencv singular values:" << endl;
    cout << matD << endl;
    cout << "left singular vectors:" << endl;
    cout << matU << endl;
    cout << "transposed matrix of right singular values:" << endl;
    cout << matVt << endl;
//    cv::Mat matUt, matV;
//    cv::transpose(matU, matUt);
//    cv::transpose(matVt, matV);
//
//    cout << "verify if the result is correct:" << endl;
//    cv::Mat reconMatD = matUt * mat * matV;
//    cout << reconMatD << endl;
    return 0;
}


int eigen_svd(vector<vector<float>>& vec){
//    cout << "source matrix：" << endl;
//    print_matrix(vec);

    cout << "opencv singular value decomposition:" << endl;
    const auto rows = (int)vec.size();
    const auto cols = (int)vec[0].size();

    std::vector<float> vec_;
    for (int i = 0; i < rows; ++i) {
        vec_.insert(vec_.begin() + i * cols, vec[i].begin(), vec[i].end());
    }

    Eigen::Map<Eigen::Matrix<float, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>> m(vec_.data(), rows, cols);

    cout << "source matrix:" << endl;
    std::cout << m << std::endl;

    Eigen::JacobiSVD<Eigen::MatrixXf> svd(m, Eigen::ComputeFullV | Eigen::ComputeFullU); // ComputeThinU | ComputeThinV
    Eigen::MatrixXf singular_values = svd.singularValues();
    Eigen::MatrixXf left_singular_vectors = svd.matrixU();
    Eigen::MatrixXf right_singular_vectors = svd.matrixV();

    cout << "eigen singular values:" << singular_values << endl;
    cout << "left singular vectors:" << left_singular_vectors << endl;
    cout << "right singular vectors:" << right_singular_vectors << endl;
    return 0;
}



// ================================= 矩阵奇异值分解 =================================
template<typename _Tp>
static void JacobiSVD(std::vector<std::vector<_Tp>>& At,
                      std::vector<std::vector<_Tp>>& _W, std::vector<std::vector<_Tp>>& Vt)
{
    double minval = FLT_MIN;
    auto eps = (_Tp)(FLT_EPSILON * 2);
    const int m = At[0].size();
    const auto n = (int)_W.size();
    const int n1 = m; // urows
    std::vector<double> W(_W.size(), 0.);

    for (int i = 0; i < n; i++) {
        double sd{0.};
        for (int k = 0; k < m; k++) {
            _Tp t = At[i][k];
            sd += (double)t*t;
        }
        W[i] = sd;

        for (int k = 0; k < n; k++)
            Vt[i][k] = 0;
        Vt[i][i] = 1;
    }

    int max_iter = std::max(m, 30);
    for (int iter = 0; iter < max_iter; iter++) {
        bool changed = false;
        _Tp c, s;

        for (int i = 0; i < n - 1; i++) {
            for (int j = i + 1; j < n; j++) {
                _Tp *Ai = At[i].data(), *Aj = At[j].data();
                double a = W[i], p = 0, b = W[j];

                for (int k = 0; k < m; k++)
                    p += (double)Ai[k] * Aj[k];

                if (std::abs(p) <= eps * std::sqrt((double)a*b))
                    continue;

                p *= 2;
                double beta = a - b, gamma = hypot((double)p, beta);
                if (beta < 0) {
                    double delta = (gamma - beta)*0.5;
                    s = (_Tp)std::sqrt(delta / gamma);
                    c = (_Tp)(p / (gamma*s * 2));
                } else {
                    c = (_Tp)std::sqrt((gamma + beta) / (gamma * 2));
                    s = (_Tp)(p / (gamma*c * 2));
                }

                a = b = 0;
                for (int k = 0; k < m; k++) {
                    _Tp t0 = c*Ai[k] + s*Aj[k];
                    _Tp t1 = -s*Ai[k] + c*Aj[k];
                    Ai[k] = t0; Aj[k] = t1;

                    a += (double)t0*t0; b += (double)t1*t1;
                }
                W[i] = a; W[j] = b;

                changed = true;

                _Tp *Vi = Vt[i].data(), *Vj = Vt[j].data();

                for (int k = 0; k < n; k++) {
                    _Tp t0 = c*Vi[k] + s*Vj[k];
                    _Tp t1 = -s*Vi[k] + c*Vj[k];
                    Vi[k] = t0; Vj[k] = t1;
                }
            }
        }

        if (!changed)
            break;
    }

    for (int i = 0; i < n; i++) {
        double sd{ 0. };
        for (int k = 0; k < m; k++) {
            _Tp t = At[i][k];
            sd += (double)t*t;
        }
        W[i] = std::sqrt(sd);
    }

    for (int i = 0; i < n - 1; i++) {
        int j = i;
        for (int k = i + 1; k < n; k++) {
            if (W[j] < W[k])
                j = k;
        }
        if (i != j) {
            std::swap(W[i], W[j]);

            for (int k = 0; k < m; k++)
                std::swap(At[i][k], At[j][k]);

            for (int k = 0; k < n; k++)
                std::swap(Vt[i][k], Vt[j][k]);
        }
    }

    for (int i = 0; i < n; i++)
        _W[i][0] = (_Tp)W[i];

    srand(time(nullptr));

    for (int i = 0; i < n1; i++) {
        double sd = i < n ? W[i] : 0;

        for (int ii = 0; ii < 100 && sd <= minval; ii++) {
            // if we got a zero singular value, then in order to get the corresponding left singular vector
            // we generate a random vector, project it to the previously computed left singular vectors,
            // subtract the projection and normalize the difference.
            const auto val0 = (_Tp)(1. / m);
            for (int k = 0; k < m; k++) {
                unsigned long int rng = rand() % 4294967295; // 2^32 - 1
                _Tp val = (rng & 256) != 0 ? val0 : -val0;
                At[i][k] = val;
            }
            for (int iter = 0; iter < 2; iter++) {
                for (int j = 0; j < i; j++) {
                    sd = 0;
                    for (int k = 0; k < m; k++)
                        sd += At[i][k] * At[j][k];
                    _Tp asum = 0;
                    for (int k = 0; k < m; k++) {
                        auto t = (_Tp)(At[i][k] - sd*At[j][k]);
                        At[i][k] = t;
                        asum += std::abs(t);
                    }
                    asum = asum > eps * 100 ? 1 / asum : 0;
                    for (int k = 0; k < m; k++)
                        At[i][k] *= asum;
                }
            }

            sd = 0;
            for (int k = 0; k < m; k++) {
                _Tp t = At[i][k];
                sd += (double)t*t;
            }
            sd = std::sqrt(sd);
        }

        _Tp s = (_Tp)(sd > minval ? 1 / sd : 0.);
        for (int k = 0; k < m; k++)
            At[i][k] *= s;
    }
}

// matSrc为原始矩阵，支持非方阵，matD存放奇异值，matU存放左奇异向量，matVt存放转置的右奇异向量
template<typename _Tp>
int svd(const std::vector<std::vector<_Tp>>& matSrc,
        std::vector<std::vector<_Tp>>& matD, std::vector<std::vector<_Tp>>& matU, std::vector<std::vector<_Tp>>& matVt)
{
    int m = matSrc.size();
    int n = matSrc[0].size();
    for (const auto& sz : matSrc) {
        if (n != sz.size()) {
            fprintf(stderr, "matrix dimension dismatch
");
            return -1;
        }
    }

    bool at = false;
    if (m < n) {
        std::swap(m, n);
        at = true;
    }

    matD.resize(n);
    for (int i = 0; i < n; ++i) {
        matD[i].resize(1, (_Tp)0);
    }
    matU.resize(m);
    for (int i = 0; i < m; ++i) {
        matU[i].resize(m, (_Tp)0);
    }
    matVt.resize(n);
    for (int i = 0; i < n; ++i) {
        matVt[i].resize(n, (_Tp)0);
    }
    std::vector<std::vector<_Tp>> tmp_u = matU, tmp_v = matVt;

    std::vector<std::vector<_Tp>> tmp_a, tmp_a_;
    if (!at)
        transpose(matSrc, tmp_a);
    else
        tmp_a = matSrc;

    if (m == n) {
        tmp_a_ = tmp_a;
    } else {
        tmp_a_.resize(m);
        for (int i = 0; i < m; ++i) {
            tmp_a_[i].resize(m, (_Tp)0);
        }
        for (int i = 0; i < n; ++i) {
            tmp_a_[i].assign(tmp_a[i].begin(), tmp_a[i].end());
        }
    }
    JacobiSVD(tmp_a_, matD, tmp_v);

    if (!at) {
        transpose(tmp_a_, matU);
        matVt = tmp_v;
    } else {
//        transpose(tmp_v, matVt);
//        matU = tmp_a_;
    }
    return 0;
}

int cpp_svd(vector<vector<float>>& vec){
    std::vector<std::vector<float>> matD, matU, matVt;
    if (svd(vec, matD, matU, matVt) != 0) {
        fprintf(stderr, "C++ implement singular value decomposition fail
");
        return -1;
    }
//    cout << "singular values:" << endl;
//    print_matrix(matD);
//    cout << "left singular vectors" << endl;
//    print_matrix(matU);
//    cout << "transposed matrix of right singular values:" << endl;
//    print_matrix(matVt);
    return 0;
}


int test_SVD()
{
    std::vector<std::vector<float>> vec{ { 0.68f, 0.597f, 0.2752f, 0.68f, 0.597f, 0.2752f },
                                         { -0.211f, 0.823f, -0.9273f, 0.68f, 0.597f, 0.2752f },
                                         { 0.566f, -0.605f, 0.1260f, 0.68f, 0.597f, 0.2752f } };
    auto time = (double)getTickCount();
    const int repeat = 100;
    for(int i = 0;i<repeat;i++){
        opencv_svd(vec); // 对比下来这个最快
    }
    auto time2 = ((double)getTickCount() - time) / (repeat * getTickFrequency());
    for(int i = 0;i<repeat;i++){
        eigen_svd(vec);
    }
    auto time3 = ((double)getTickCount() - time2) / (repeat * getTickFrequency());
    for(int i = 0;i<repeat;i++){
        cpp_svd(vec);
    }
    auto time4 = ((double)getTickCount() - time3) / (repeat * getTickFrequency());
    time = 0;
    return 0;
}