zoukankan      html  css  js  c++  java
  • 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;
    }
    
    
    
  • 相关阅读:
    C#深入浅出 修饰符(二)
    HDU 5785 Interesting
    HDU 5783 Divide the Sequence
    HDU 5781 ATM Mechine
    UVA 714 Copying Books
    uva 1471 Defense Lines
    UVA 11134 Fabled Rooks
    UVA 11572 Unique Snowflakes
    UVA 11093 Just Finish it up
    UVA 10954 Add All
  • 原文地址:https://www.cnblogs.com/theodoric008/p/9506015.html
Copyright © 2011-2022 走看看