各種ライブラリを利用した特異値分解(SVD)のC/C++&Python&MATLABサンプルソースコード

Tweet

特異値分解のサンプルソースコードをいくつか紹介します．といっても，特異値分解そのものはライブラリの関数を呼び出すだけです．MATLABによるSVDの実装，PythonとnumpyによるSVDの実装などを載せています．また，C/C++のプログラムも載せています．C/C++は，Eigenライブラリによる特異値分解のソースコード，OpenCVによる特異値分解のソースコードなどを載せています．

MATLABによる特異値分解(SVD)のサンプル ソースコード

ROW = randi([2 10]);
COL = randi([2 10]);
A = 99.9*rand(ROW, COL);
[U,W,V] = svd(A,'econ');
VT = V.';
disp("SVD result")
disp("A =")
disp(A)
disp("U =")
disp(U)
disp("W =")
disp(W)
disp("V^T =")
disp(VT)
check = U * W * VT;
disp("Check")
disp("U W V^T =")
disp(check)
disp("Coincide with A")

Pythonとnumpyによる特異値分解(SVD)のサンプルソースコード

import numpy as np

def main():
    ROW = np.random.randint(2, 10)
    COL = np.random.randint(2, 10)
    A = np.random.uniform(0.0, 99.9, (ROW, COL))
    U,W,VT = np.linalg.svd(A, full_matrices=False)
    print("SVD result")
    print("A=", A)
    print("U=", U)
    print("W=", W)
    print("V^T=", VT)
    matw = np.diag(W)
    check = np.dot(np.dot(U, matw), VT)
    print("Check")
    print("U W V^T=", check)
    print("Coincide with A")

if __name__ == '__main__':
    main()

PythonとOpenCVによる特異値分解(SVD)のサンプルソースコード

import numpy as np
import cv2

def main():
    ROW = np.random.randint(2, 10)
    COL = np.random.randint(2, 10)
    A = np.random.uniform(0.0, 99.9, (ROW, COL))
    W,U,VT = cv2.SVDecomp(A)
    print("SVD result")
    print("A=", A)
    print("U=", U)
    print("W=", W)
    print("V^T=", VT)
    matw = np.diagflat(W)
    check = np.dot(np.dot(U, matw), VT)
    print("Check")
    print("U W V^T=", check)
    print("Coincide with A")

if __name__ == '__main__':
    main()

C/C++とEigenライブラリによる特異値分解(SVD)のサンプルソースコード

#include <stdlib.h>
#include <time.h>

#include <iostream>
using namespace std;

#include "Eigen/Core"
#include "Eigen/SVD"
using namespace Eigen;

int main()
{
    srand((unsigned)time(NULL));
    int ROW = rand() % 8 + 2;
    int COL = rand() % 8 + 2;
    MatrixXd A = (MatrixXd::Random(ROW, COL) + MatrixXd::Ones(ROW, COL)) / 2.0 * 99.9;
    JacobiSVD<MatrixXd> svd(A, ComputeThinU | ComputeThinV);
    MatrixXd VT = svd.matrixV().transpose();
    cout << "SVD result" << endl << endl;
    cout << "A=" << endl << A << endl << endl;
    cout << "U=" << endl << svd.matrixU() << endl << endl;
    cout << "W=" << endl << svd.singularValues() << endl << endl;
    cout << "V^T=" << endl << VT << endl << endl;
    cout << endl;
    MatrixXd W = svd.singularValues().asDiagonal();
    MatrixXd check = svd.matrixU() * W * VT;
    cout << "Check" << endl << endl;
    cout << "U W V^T=" << endl << kenzan << endl << endl;
    cout << "Coincide with A" << endl << endl;
    cout << "Push enter" << endl;
    cin.get();
    return 0;
}

C/C++とEigenライブラリによる特異値分解(SVD)のサンプルソースコード

#include <stdlib.h>
#include <time.h>

#include <iostream>
using namespace std;

#include "Eigen/Core"
#include "Eigen/SVD"
using namespace Eigen;

int main()
{
    srand((unsigned)time(NULL));
    int ROW = rand() % 8 + 2;
    int COL = rand() % 8 + 2;
    MatrixXd A = (MatrixXd::Random(ROW, COL) + MatrixXd::Ones(ROW, COL)) / 2.0 * 99.9;
    BDCSVD<MatrixXd> svd(A, ComputeThinU | ComputeThinV);
    MatrixXd VT = svd.matrixV().transpose();
    cout << "SVD result" << endl << endl;
    cout << "A=" << endl << A << endl << endl;
    cout << "U=" << endl << svd.matrixU() << endl << endl;
    cout << "W=" << endl << svd.singularValues() << endl << endl;
    cout << "V^T=" << endl << VT << endl << endl;
    cout << endl;
    MatrixXd W = svd.singularValues().asDiagonal();
    MatrixXd check = svd.matrixU() * W * VT;
    cout << "Check" << endl << endl;
    cout << "U W V^T=" << endl << kenzan << endl << endl;
    cout << "Coincide with A" << endl << endl;
    cout << "Push enter" << endl;
    cin.get();
    return 0;
}

C/C++とOpenCVによる特異値分解(SVD)のサンプルソースコード

// Include directory C:\OpenCV3.1\build\include
// Library directory (x64) C:\OpenCV3.1\build\x64\vc14\lib
// PATH (x64) C:\OpenCV3.1\build\x64\vc14\bin
#if _DEBUG
#pragma comment(lib, "opencv_world310d.lib")
#else
#pragma comment(lib, "opencv_world310.lib")
#endif

#include <time.h>
#include <iostream>
using namespace std;

#include <opencv2/opencv.hpp>
using namespace cv;

int main()
{
    RNG gen(getTickCount());
    int ROW = gen.uniform(2, 10);
    int COL = gen.uniform(2, 10);
    Mat A(ROW, COL, CV_64F);
    gen.fill(A, RNG::UNIFORM, 0.0, 99.9);
    SVD svd(A);
    cout << "SVD result" << endl << endl;
    cout << "A=" << A << endl << endl;
    cout << "U=" << svd.u << endl << endl;
    cout << "W=" << svd.w << endl << endl;
    cout << "V^T=" << svd.vt << endl << endl;
    cout << endl;
    Mat matw = Mat::diag(svd.w);
    Mat check = svd.u * matw * svd.vt;
    cout << "Check" << endl << endl;
    cout << "U W V^T=" << check << endl << endl;
    cout << "Coincide with A" << endl << endl;
    cout << "Push enter" << endl;
    cin.get();
    return 0;
}

C/C++とOpenCVによる特異値分解(SVD)のサンプルソースコード

// Include directory C:\OpenCV2.4.11\build\include
// Library directory (x86) C:\OpenCV2.4.11\build\x86\vc12\lib
// Library directory (x64) C:\OpenCV2.4.11\build\x64\vc12\lib
// PATH (x86) C:\OpenCV2.4.11\build\x86\vc12\bin
// PATH (x64) C:\OpenCV2.4.11\build\x64\vc12\bin
#if _DEBUG
#pragma comment(lib, "opencv_core2411d.lib")
#pragma comment(lib, "opencv_highgui2411d.lib")
#pragma comment(lib, "opencv_imgproc2411d.lib")
#else
#pragma comment(lib, "opencv_core2411.lib")
#pragma comment(lib, "opencv_highgui2411.lib")
#pragma comment(lib, "opencv_imgproc2411.lib")
#endif

#include <opencv2/opencv.hpp>
using namespace cv;

#include <stdio.h>
#include <conio.h>
#include <time.h>

int main()
{
    srand((unsigned)time(NULL));
    int ROW = rand() % 8 + 2;
    int COL = rand() % 8 + 2;
    Mat A(ROW, COL, CV_64F);
    for (int i = 0; i < A.rows; i++) {
        for (int j = 0; j < A.cols; j++) {
            A.at<double>(i, j) = 99.9 * (double)rand() / (double)RAND_MAX;
        }
    }
    SVD svd(A);
    printf("SVD result\n");
    printf("\n");
    printf("A=\n");
    for (int i = 0; i < A.rows; i++) {
        for (int j = 0; j < A.cols; j++) {
            printf("%1.2f ", A.at<double>(i, j));
        }
        printf("\n");
    }
    printf("U=\n");
    for (int i = 0; i < svd.u.rows; i++) {
        for (int j = 0; j < svd.u.cols; j++) {
            printf("%1.2f ", svd.u.at<double>(i, j));
        }
        printf("\n");
    }
    printf("W=\n");
    for (int i = 0; i < svd.w.rows; i++) {
        for (int j = 0; j < svd.w.cols; j++) {
            printf("%1.2f ", svd.w.at<double>(i, j));
        }
        printf("\n");
    }
    printf("V^T=\n");
    for (int i = 0; i < svd.vt.rows; i++) {
        for (int j = 0; j < svd.vt.cols; j++) {
            printf("%1.2f ", svd.vt.at<double>(i, j));
        }
        printf("\n");
    }
    printf("\n");
    Mat matw = Mat::zeros(svd.w.rows, svd.w.rows, CV_64F);
    for (int i = 0; i < svd.w.rows; i++) {
        matw.at<double>(i, i) = svd.w.at<double>(i, 0);
    }
    Mat check = svd.u * matw * svd.vt;
    printf("Check\n");
    printf("U W V^T=\n");
    for (int i = 0; i < check.rows; i++) {
        for (int j = 0; j < check.cols; j++) {
            printf("%1.2f ", check.at<double>(i, j));
        }
        printf("\n");
    }
    printf("Coincide with A\n");
    _getch();
    return 0;
}

C/C++とNumerical Recipes in C 3rd editionによる特異値分解(SVD)のサンプルソースコード

#include <stdio.h>
#include <conio.h>
#include <time.h>

#include "nr3/nr3.h"
#include "nr3/svd.h"

int main()
{
    srand((unsigned)time(NULL));
    int ROW = rand() % 8 + 2;
    int COL = rand() % 8 + 2;
    MatDoub A(ROW, COL);
    for (int i = 0; i < ROW; i++) {
        for (int j = 0;j < COL;j++) {
            A[i][j] = 99.9 * (double)rand() / (double)RAND_MAX;
        }
    }
    SVD svd(A);
    printf("SVD result\n");
    printf("\n");
    printf("A=\n");
    for (int i = 0; i < ROW; i++) {
        for (int j = 0; j < COL; j++) {
            printf("%1.2f ", A[i][j]);
        }
        printf("\n");
    }
    printf("U=\n");
    for (int i = 0; i < svd.u.nrows(); i++) {
        for (int j = 0; j < svd.u.ncols(); j++) {
            printf("%1.2f ", svd.u[i][j]);
        }
        printf("\n");
    }
    printf("W=\n");
    for (int i = 0; i < svd.w.size(); i++) {
        printf("%1.2f\n", svd.w[i]);
    }
    printf("V^T=\n");
    for (int i = 0; i < svd.v.ncols(); i++) {
        for (int j = 0; j < svd.v.nrows(); j++) {
            printf("%1.2f ", svd.v[j][i]);
        }
        printf("\n");
    }
    printf("\n");
    MatDoub check(ROW, COL);
    for (int i = 0; i < ROW; i++) {
        for (int j = 0; j < COL; j++) {
            double val = 0.0;
            for (int k = 0; k < svd.u.ncols(); k++) {
                val += svd.u[i][k] * svd.w[k] * svd.v[j][k];
            }
            check[i][j] = val;
        }
    }
    printf("Check\n");
    printf("U W V^T=\n");
    for (int i = 0; i < check.nrows(); i++) {
        for (int j = 0; j < check.ncols(); j++) {
            printf("%1.2f ", check[i][j]);
        }
        printf("\n");
    }
    printf("Coincide with A\n");
    printf("Push any key\n");
    _getch();
    return 0;
}

もどる