奇異值分解(singular value decomposition, SVD)：將矩陣分解為奇異向量(singular vector)和奇異值(singular value)。通過奇異值分解，我們會得到一些與特征分解相同類型的信息。然而，奇異值分解有更廣泛的應用。每個實數矩陣都有一個奇異值分解，但不一定都有特征分解。例如，非方陣的矩陣沒有特征分解，這是我們只能使用奇異值分解。

將矩陣A分解成三個矩陣的乘積：A=UDV^T。

假設A是一個m*n的矩陣，那么U是一個m*m的矩陣，D是一個m*n的矩陣，V是一個n*n矩陣。

這些矩陣中的每一個經定義后都擁有特殊的結構。矩陣U和V都被定義為正交矩陣，而矩陣D被定義為對角矩陣。注意，矩陣D不一定是方陣。對角矩陣D對角線上的元素被稱為矩陣A的奇異值(singular value)。矩陣U的列向量被稱為左奇異向量(left singular vector)，矩陣V的列向量被稱為右奇異向量(right singular vector)。A的左奇異向量是AA^T的特征向量。A的右奇異向量是A^TA的特征向量。A的非零奇異值是A^TA特征值的平方根，同時也是AA^T特征值的平方根。

奇異值分解(singular value decomposition)是線性代數中一種重要的矩陣分解，在信號處理、統計學等領域有重要應用。奇異值分解能夠用于任意m*n矩陣，而特征分解只能適用于特定類型的方陣，故奇異值分解的適用范圍更廣。

以上內容摘自：?《深度學習中文版》?和 ?維基百科

基于雅克比(Jacobi)方法對矩陣進行奇異值分解操作步驟(參考OpenCV中實現)：

(1)、對原始矩陣A(m*n)進行判斷，若m<n(即行數小于列數)，則交換m、n值；若m≥n，則對A進行轉置變換；

(2)、初始化臨時變量D、U、Vt：D為奇異值，為n行1列；U為左奇異向量，為m行m列；Vt為轉置后的右奇異向量，為n行n列；并初始化D、U、Vt值均為0；

(3)、初始化臨時變量A′：A′為m行m列，并將A的值賦值給A′，A′中多余元素賦初值為0；

(4)、由A′、D、Vt開始進行基于Jacobi方法的奇異值分解；

(5)、設置臨時變量W，長度為n，將A′中前n行中，每行元素的平方和賦值給W；

(6)、設置Vt為單位矩陣；

(7)、循環(huán)計算旋轉矩陣，并更新A′、W、Vt對應位置的值；最大循環(huán)次數為std::max(m, 30)；

(8)、重置W值為A′中前n行，每行元素平方和的開方；

(9)、將W中元素按照從大到小排序，排序后的W即為D中主對角線元素值；

(10)、按照(9)中對W的排序規(guī)則對A′和Vt進行排序；

(11)、計算A′中值；

(12)、最終的A′和Vt即為所求的左、右奇異向量。

以下是分別采用C++(參考opencv sources/modules/core/src/lapack.cpp)和OpenCV實現的矩陣奇異值分解：

#include "funset.hpp"
#include <math.h>
#include <iostream>
#include <string>
#include <vector>
#include <opencv2/opencv.hpp>
#include "common.hpp"// ================================= 矩陣奇異值分解 =================================
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;_Tp eps = (_Tp)(FLT_EPSILON * 2);const int m = At[0].size();const int n = _W.size();const int n1 = m; // urowsstd::vector<double> W(n, 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 _Tp val0 = (_Tp)(1. / m);for (int k = 0; k < m; k++) {unsigned 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++) {_Tp 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\n");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);elsetmp_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, matU);matVt = tmp_a_;}return 0;
}int test_SVD()
{//std::vector<std::vector<float>> vec{ { 1.2f, 2.5f, 5.6f, -2.5f },//				{ -3.6f, 9.2f, 0.5f, 7.2f },//				{ 4.3f, 1.3f, 9.4f, -3.4f },//				{ 6.4f, 0.1f, -3.7f, 0.9f } };//const int rows{ 4 }, cols{ 4 };//std::vector<std::vector<float>> vec{ { 1.2f, 2.5f, 5.6f, -2.5f },//				{ -3.6f, 9.2f, 0.5f, 7.2f },//				{ 4.3f, 1.3f, 9.4f, -3.4f } };//const int rows{ 3 }, cols{ 4 };std::vector<std::vector<float>> vec{ { 0.68f, 0.597f },{ -0.211f, 0.823f },{ 0.566f, -0.605f } };const int rows{ 3 }, cols{ 2 };fprintf(stderr, "source matrix:\n");print_matrix(vec);fprintf(stderr, "\nc++ implement singular value decomposition:\n");std::vector<std::vector<float>> matD, matU, matVt;if (svd(vec, matD, matU, matVt) != 0) {fprintf(stderr, "C++ implement singular value decomposition fail\n");return -1;}fprintf(stderr, "singular values:\n");print_matrix(matD);fprintf(stderr, "left singular vectors:\n");print_matrix(matU);fprintf(stderr, "transposed matrix of right singular values:\n");print_matrix(matVt);fprintf(stderr, "\nopencv singular value decomposition:\n");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);}}/*w calculated singular valuesu calculated left singular vectorsvt transposed matrix of right singular vectors*/cv::Mat w, u, vt, v;cv::SVD::compute(mat, w, u, vt, 4);//cv::transpose(vt, v);fprintf(stderr, "singular values:\n");print_matrix(w);fprintf(stderr, "left singular vectors:\n");print_matrix(u);fprintf(stderr, "transposed matrix of right singular values:\n");print_matrix(vt);return 0;
}

執(zhí)行結果如下：

以下是采用Eigen實現的矩陣奇異值分解code：

#include "funset.hpp"
#include <math.h>
#include <iostream>
#include <vector>
#include <string>
#include <opencv2/opencv.hpp>
#include <Eigen/Dense>
#include "common.hpp"int test_SVD()
{//std::vector<std::vector<float>> vec{ { 1.2f, 2.5f, 5.6f, -2.5f },//				{ -3.6f, 9.2f, 0.5f, 7.2f },//				{ 4.3f, 1.3f, 9.4f, -3.4f },//				{ 6.4f, 0.1f, -3.7f, 0.9f } };//const int rows{ 4 }, cols{ 4 };//std::vector<std::vector<float>> vec{ { 1.2f, 2.5f, 5.6f, -2.5f },//				{ -3.6f, 9.2f, 0.5f, 7.2f },//				{ 4.3f, 1.3f, 9.4f, -3.4f } };//const int rows{ 3 }, cols{ 4 };std::vector<std::vector<float>> vec{ { 0.68f, 0.597f },{ -0.211f, 0.823f },{ 0.566f, -0.605f } };const int rows{ 3 }, cols{ 2 };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);fprintf(stderr, "source matrix:\n");std::cout << m << std::endl;Eigen::JacobiSVD<Eigen::MatrixXf> svd(m, Eigen::ComputeFullV | Eigen::ComputeFullU); // ComputeThinU | ComputeThinVEigen::MatrixXf singular_values = svd.singularValues();Eigen::MatrixXf left_singular_vectors = svd.matrixU();Eigen::MatrixXf right_singular_vectors = svd.matrixV();fprintf(stderr, "singular values:\n");print_matrix(singular_values.data(), singular_values.rows(), singular_values.cols());fprintf(stderr, "left singular vectors:\n");print_matrix(left_singular_vectors.data(), left_singular_vectors.rows(), left_singular_vectors.cols());fprintf(stderr, "right singular vecotrs:\n");print_matrix(right_singular_vectors.data(), right_singular_vectors.rows(), right_singular_vectors.cols());return 0;
}

執(zhí)行結果如下：

由以上結果可見：C++、OpenCV、Eigen結果是相同的。

GitHub：

https://github.com/fengbingchun/NN_Test

https://github.com/fengbingchun/Eigen_Test

總結

以上是生活随笔為你收集整理的矩阵奇异值分解简介及C++/OpenCV/Eigen的三种实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。