矩阵的SVD分解

矩阵的SVD分解

#include<opencv2/opencv.hpp>
/*
用对角阵左乘一个矩阵，就是用对角阵的对角元分别乘这个矩阵的对应各行；
用对角阵右乘一个矩阵，就是用对角阵的对角元分别乘这个矩阵的对应各列。

只有对角线上有非0元素的矩阵称为对角矩阵，或说若一个方阵除了主对角线上的元素外，
其余元素都等于零。

对角线上的元素相等的对角矩阵称为数量矩阵，对角线上的元素都为1的n阶对角（矩）阵
称为单位（矩）阵。
*/

cv::Mat mysvd(cv::Mat mat) {

	cv::Mat U, V;
	cv::Mat tar;
	//转置矩阵T
	cv::Mat t = mat.t();

	cv::Mat mValuesMat;		//特征值
	cv::Mat mVectorsMat;	//特征向量
	cv::eigen(mat * t, mValuesMat, mVectorsMat);
	U = mVectorsMat;

	cv::Mat tValuesMat;
	cv::Mat tVectorsMat;
	cv::eigen(t * mat, tValuesMat, tVectorsMat);
	V = tValuesMat;

	cv::Mat S = cv::Mat(cv::Size(cv::countNonZero(mValuesMat) + cv::countNonZero(tValuesMat), 
		cv::countNonZero(mValuesMat) + cv::countNonZero(tValuesMat)),
		mValuesMat.type());

	for (int i = 0; i < cv::countNonZero(mValuesMat) - 1; i++) {
		if (sqrt(mValuesMat.at<double>(i) >= 0)) {
			S.at<double>(i, i) = sqrt(mValuesMat.at<double>(i));
		}
	}
	for (int j = 0; j < cv::countNonZero(tValuesMat); j++) {
		if (sqrt(tValuesMat.at<double>(j)) >= 0) {
			S.at<double>(cv::countNonZero(mValuesMat) + j, cv::countNonZero(mValuesMat) + j) =
				sqrt(tValuesMat.at<double>(j));
		}
	}
	for (int i = 0; i < S.rows; i++) {
		double* data = S.ptr<double>(i);
		for (int j = 0; j < S.cols; j++) {
			if (data[j] < 0) {
				data[j] = 0;
			}
		}

	}
	//std::cout << U.type() << '\t' << S.type() << '\t' << V.type() << std::endl;
	//std::cout << U.cols << '\t' << S.rows << '\t';
	return S;
}

//求矩阵的SVD分解
int main() {
	cv::Mat C = (cv::Mat_<double>(3, 3) <<
		1, 2, 3, 4, 5, 6, 7, 9, 10);		//直接赋初始值的方法
	mysvd(C);
	std::cout << mysvd(C) << std::endl;
	cv::Mat w, wt, u, vt;
	cv::SVDecomp(C, w, u, vt);
	cv::SVDecomp(C.t(), wt, u, vt);
	std::cout << w << std::endl;
	std::cout << wt << std::endl;
	return 0;
}