矩阵奇异值分解简介及C++/OpenCV/Eigen的三种实现
Posted fengbingchun
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了矩阵奇异值分解简介及C++/OpenCV/Eigen的三种实现相关的知识，希望对你有一定的参考价值。
奇异值分解(singular value decomposition, SVD)：将矩阵分解为奇异向量(singular vector)和奇异值(singular value)。通过奇异值分解，我们会得到一些与特征分解相同类型的信息。然而，奇异值分解有更广泛的应用。每个实数矩阵都有一个奇异值分解，但不一定都有特征分解。例如，非方阵的矩阵没有特征分解，这是我们只能使用奇异值分解。
将矩阵A分解成三个矩阵的乘积：A=UDVT。
假设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的左奇异向量是AAT的特征向量。A的右奇异向量是ATA的特征向量。A的非零奇异值是ATA特征值的平方根，同时也是AAT特征值的平方根。
奇异值分解(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)、循环计算旋转矩阵，并更新A′、W、Vt对应位置的值；最大循环次数为std::max(m, 30)；
(8)、重置W值为A′中前n行，每行元素平方和的开方；
(9)、将W中元素按照从大到小排序，排序后的W即为D中主对角线元素值；
(10)、按照(9)中对W的排序规则对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; // urows
std::vector<double> W(n, 0.);
for (int i = 0; i < n; i++)
double sd0.;
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);
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, 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 values
u calculated left singular vectors
vt 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;
以下是采用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 | ComputeThinV
Eigen::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;
由以上结果可见：C++、OpenCV、Eigen结果是相同的。
GitHub：
https://github.com/fengbingchun/NN_Test
https://github.com/fengbingchun/Eigen_Test
以上是关于矩阵奇异值分解简介及C++/OpenCV/Eigen的三种实现的主要内容，如果未能解决你的问题，请参考以下文章