在说明前,我也是查了大量文档,弄清楚各个名词的意思,才写下这篇博客。。。

特征值和特征向量:

根据公式:

opencv 矩阵分解 opencv svd分解_SVD奇异值分解

A是n*n的方阵(必须是方阵),x是特征向量,λ是特征值。

一般情况下,会有n个特征值,和n个特征向量。(这里的一般是指方阵是满秩,各行各列都是线性无关)。

引出问题:如果不是n*n维的矩阵,怎么求?

假设矩阵X是m*n,则可以求

opencv 矩阵分解 opencv svd分解_特征向量_02


opencv 矩阵分解 opencv svd分解_i++_03

矩阵的特征分解。

一般的非方阵的特征分解称为奇异值分解。

奇异值分解(Singular Value Decomposition)是线性代数中一种重要的矩阵分解,奇异值分解则是特征分解在任意矩阵上的推广。在信号处理、统计学等领域有重要应用。

奇异值分解SVD:

A是m*n大小的矩阵

根据公式:

opencv 矩阵分解 opencv svd分解_SVD奇异值分解_04

U是左奇异矩阵

同理V是右奇异矩阵

剩下一个S:奇异值矩阵,理论上S的大小应该是m*n才能满足矩阵乘法运算。

opencv 矩阵分解 opencv svd分解_特征向量_05

=》

opencv 矩阵分解 opencv svd分解_#include_06

=》

opencv 矩阵分解 opencv svd分解_SVD奇异值分解_07

理解:

opencv 矩阵分解 opencv svd分解_#include_08

, V是单位正交基,不知道的可百度。理解:

opencv 矩阵分解 opencv svd分解_SVD奇异值分解_09

,S是奇异值,一般的是n行一列的矩阵,为了方便计算,我们设计S变成对角矩阵,即n*1 =》n*n。

右奇异矩阵举例:

其中V是矩阵

opencv 矩阵分解 opencv svd分解_特征向量_10

的特征向量,   n*n大小 

opencv 矩阵分解 opencv svd分解_SVD奇异值分解_11

得到奇异值

opencv 矩阵分解 opencv svd分解_opencv 矩阵分解_12

,且S是n列矩阵,V是n*n维,反求左奇异值

opencv 矩阵分解 opencv svd分解_特征向量_13

  =》

opencv 矩阵分解 opencv svd分解_SVD奇异值分解_14

最后得到U是m*n维即

opencv 矩阵分解 opencv svd分解_#include_15

左奇异矩阵举例:

其中U是矩阵

opencv 矩阵分解 opencv svd分解_SVD奇异值分解_16

的特征向量 m*m大小

opencv 矩阵分解 opencv svd分解_i++_17

得到奇异值

opencv 矩阵分解 opencv svd分解_opencv 矩阵分解_18

,且S是m列矩阵,U是m*m维,反求右奇异值

opencv 矩阵分解 opencv svd分解_i++_19

  =》

opencv 矩阵分解 opencv svd分解_#include_20

得到V是n*m维(这里是U的转置)即

opencv 矩阵分解 opencv svd分解_SVD奇异值分解_21

opencv求SVD:

#include <opencv2/opencv.hpp>
#include <iostream>

using namespace cv;
using namespace std;
void main()
{
	Mat src = imread("1.jpg", IMREAD_GRAYSCALE);
	Mat src_ = src.clone();
	src.convertTo(src, CV_64FC1);
	Mat U, W, V;
	SVD::compute(src, W, U, V);
	int set_dim = min(src.rows, src.cols);
	//cout << "W >> <" << W.rows << "," << W.cols << ">" << endl;
	//cout << "U >> <" << U.rows << "," << U.cols << ">" << endl;
	//cout << "V >> <" << V.rows << "," << V.cols << ">" << endl;
	cout << "W value >> " << W << endl;

	Mat W_ = Mat(set_dim, set_dim, CV_64FC1, Scalar(0));
	double ratio = 0.01;
	int set_rows = set_dim * ratio;
	for (int i = 0; i < set_rows; i++)
	{
		W_.at<double>(i, i) = W.at<double>(i, 0);
	}
	Mat dst = U * W_ * V;
	system("pause");
}

opencv里的W,U,V分别对应前文的 S,U,V^T 且维数计算是左奇异矩阵举例的模式。

利用eigen自己实现svd:

详细看注释

void main()
{

	Mat m = imread("1.jpg", IMREAD_GRAYSCALE);
	Mat src = m.clone();
	Mat S, U;
	m.convertTo(m, CV_64FC1);
	double t1 = getTickCount();
	eigen(m * m.t(), S, U); //求出奇异值矩阵 和 左奇异矩阵
	U = U.t();  //根据公式 U变成转置后的
	int set_dim = min(src.rows, src.cols);
	
	for (int i = 0; i < S.rows; i++)
	{
		double val = S.at<double>(i, 0);
		double sq_val = sqrt(val);
		S.at<double>(i, 0) = sq_val; //λ = sqrt(S)
	}
	
	Mat V;
	V = m.t() * U; //根据公式
	for (int i = 0; i < V.cols; i++)
	{
		for (int j = 0; j < V.rows; j++)
		{
			V.at<double>(j, i) = V.at<double>(j, i) / S.at<double>(i, 0);
		}
	}
	

	Mat Vt = V.t();

	Mat S_ = Mat(S.rows, S.rows, CV_64FC1, Scalar(0));

	for (int i = 0; i < S.rows; i++)//列向量 变成 对角矩阵
	{
		S_.at<double>(i, i) = S.at<double>(i, 0);
	}

	dst = U * S_ * Vt;
	
	imshow("src", src);
	imshow("dst", dst);
	waitKey();

}

降维分析:

一直S是对角矩阵,且按降序排列,根据实际测试前10%甚至1%的奇异值就能涵盖99%以上的特征。

可以设置r = ratio * n (ratio:比率)

公式:

opencv 矩阵分解 opencv svd分解_SVD奇异值分解_22

  (公式1)=》

opencv 矩阵分解 opencv svd分解_特征向量_23

(公式2)

如果取r = 0.01,降维后公式2维度是公式1 的 亿分之一。

#include <opencv2/opencv.hpp>
#include <iostream>

using namespace cv;
using namespace std;


Mat getForeMatrix(int rows, int cols, Mat & img, int type)
{
	Mat tmp = Mat(rows, cols, type, Scalar(0));
	for (int i = 0; i < rows; i++)
	{
		for (int j = 0; j < cols; j++)
		{
			if (type == CV_8UC1)
				tmp.at<uchar>(i, j) = img.at<uchar>(i, j);
			else if (type == CV_32FC1)
				tmp.at<float>(i, j) = img.at<float>(i, j);
			else
				tmp.at<double>(i, j) = img.at<double>(i, j);
		}
	}
	return tmp;
}
void main()
{
    Mat src = imread(filename, IMREAD_GRAYSCALE);
    Mat src_ = src.clone();
    src.convertTo(src, CV_64FC1);
    Mat U, W, V;
    SVD::compute(src, W, U, V, SVD::Flags::MODIFY_A);

    double t2 = (getTickCount() - t1) / getTickFrequency();
    cout << "cost time >> " << t2 << endl;
    int set_dim = min(src.rows, src.cols);
    //cout << "W >> <" << W.rows << "," << W.cols << ">" << endl;
    //cout << "U >> <" << U.rows << "," << U.cols << ">" << endl;
    //cout << "V >> <" << V.rows << "," << V.cols << ">" << endl;
    //cout << "W value >> " << W << endl;
    
    double ratio = 0.01;
    int set_rows = set_dim * ratio;
    Mat W_ = Mat(set_rows , set_rows , CV_64FC1, Scalar(0));
    for (int i = 0; i < set_rows ; i++)
    {
	W_.at<double>(i, i) = W.at<double>(i, 0);
    }
    Mat U_ = getForeMatrix(set_dim, set_rows , U, CV_64FC1);
    Mat V_ = getForeMatrix(set_rows , V.cols, V, CV_64FC1);

    Mat dst = U_ * W_ * V_;
    system("pause")
}

不正确的希望指出。。。不胜感激

更新错误日志2019.12.23,看了下左奇异矩阵举例:不需要将U转置。且更改了一些定义问题

最后:若要显示最终图像,需要将CV_64FC1图像转成CV_8UC1格式。