在说明前,我也是查了大量文档,弄清楚各个名词的意思,才写下这篇博客。。。
特征值和特征向量:
根据公式:
A是n*n的方阵(必须是方阵),x是特征向量,λ是特征值。
一般情况下,会有n个特征值,和n个特征向量。(这里的一般是指方阵是满秩,各行各列都是线性无关)。
引出问题:如果不是n*n维的矩阵,怎么求?
假设矩阵X是m*n,则可以求
或
矩阵的特征分解。
一般的非方阵的特征分解称为奇异值分解。
奇异值分解(Singular Value Decomposition)是线性代数中一种重要的矩阵分解,奇异值分解则是特征分解在任意矩阵上的推广。在信号处理、统计学等领域有重要应用。
奇异值分解SVD:
A是m*n大小的矩阵
根据公式:
U是左奇异矩阵
同理V是右奇异矩阵
剩下一个S:奇异值矩阵,理论上S的大小应该是m*n才能满足矩阵乘法运算。
=》
=》
理解:
, V是单位正交基,不知道的可百度。理解:
,S是奇异值,一般的是n行一列的矩阵,为了方便计算,我们设计S变成对角矩阵,即n*1 =》n*n。
右奇异矩阵举例:
其中V是矩阵
的特征向量, n*n大小
得到奇异值
,且S是n列矩阵,V是n*n维,反求左奇异值
=》
最后得到U是m*n维即
左奇异矩阵举例:
其中U是矩阵
的特征向量 m*m大小
得到奇异值
,且S是m列矩阵,U是m*m维,反求右奇异值
=》
得到V是n*m维(这里是U的转置)即
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:比率)
公式:
(公式1)=》
(公式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格式。