C++矩阵求逆及求数据动态增长的最小二乘中的(X‘X)^-1(一种容易实现的方式,包含详细的证明)

这是一篇关于用C++计算矩阵的逆的方法,是自己几天奋斗的结果,从有想法开始到详细的数学证明再到计算机程序的实现以及验证。废话就不多说了,按照我的习惯,这里还是先列出这篇文中的结果,方便对号入座,如不相关那就赶紧换一个把:

          1.一种容易实现的矩阵求逆算法,且可以判断矩阵是否可逆;

          2.这种算法支持数据的动态增长,例如,最小二乘中计算(X’X)^-1时当来一个新数据时,可以通过简单地计算得到最后的逆矩阵,而不必从头开始。

          3.数学上论证了正确性,并给处计算和测试程序

          4.关于该问题在矩阵乘法上的扩展

这里如果读者知道一点OenpCV或者matlab会很有裨益的,当然不知道也没有关系,只是如果你要实现它的话你还要实现矩阵的乘法和数乘两种运算就好,有人可能会问数乘是什么,数乘就是矩阵和一个实数相乘,也就是矩阵的每一个元素都乘上该实数,当然一种结构来存储该矩阵也是必要的,最好写个类。

        那么这里可能会问,我都会使用OpenCV(matlab)求个逆还不就是inv(),也就是一行代码的事而已,何必还要在这里看你在这里罗嗦!注意这里如下几个原因让你看下去:

          首先这里只使用OpenCV中的Mat 的较好的封装结构,而这不只必要的,正如前面介绍的那么样逆大可自己实现拥有那两个操作的数据结构,而不必借助于OpenCV或matlab.这里我使用它只要是因为我是做计算机视觉的电脑上自然安装了OpenCV库,而且也是为了方便。

          再来就是,这里本着的是,知其然,还要知其所以然的原则(当然我并不知道OpenCV和matlab到底是这样实现求逆的,但是对于本文中的方法我是清楚的),还因该方法实现简单,编码量少。

          这里,也是为了交流的目的,后面会提到本文的中的一些扩展可自己的一些想法,希望与各位交流。

          好了!不多说了,直接先上程序:

bool invMat(Mat origMat, Mat& invert)//对origMat求逆,并将结果返回到invert中
{
	Size sz=origMat.size();
	if(sz.height!=sz.width||origMat.depth()!=CV_64F) //这里只处理double类型的方阵
	  return 0; 
	double disc;// 用于记录中间的变量,还有一个作用是判别矩阵是否可逆
	Mat I=Mat::eye(sz,origMat.type());
	Mat D=origMat-I;
	Mat temp;
	disc=1+D.at<double>(0,0);
	if(abs(disc)<1e-12) // 认为小于1e-12为0,这里可以根据个人习惯可运用需要来定
	  return 0; // 矩阵不可逆
	else
	{
		Mat invCk(I-(1/disc)*(I.col(0)*D.row(0)));
		for(long step=1;step<sz.height;step++)
		{
			temp=D.row(step)*invCk;
			disc=1+temp.at<double>(0,step);
			if(abs(disc)<1e-12) // 同上
			  return 0;// 同前意义
			else
			  invCk=invCk-(invCk.col(step)*temp)/disc;
		}
		invert=invCk.clone(); // 拷贝到目标Matrix中;
	}
	return 1;
}

其中Mat就是OpenCV中的矩阵类,不知道的读者就理解为封装了矩阵操作的类就好。

你可能会问这是啥!我说这就是矩阵求逆的算法,详细证明如下(其实它就是一个迭代过程,也就是下面证明的结论部分——两句话):

这里以截图的形式给出(之前是写在word中的),方便阅读:

这里先给出矩阵的求逆引理,及其推论:

《C++矩阵求逆及求数据动态增长的最小二乘中的(X‘X)^-1(一种容易实现的方式,包含详细的证明)》《C++矩阵求逆及求数据动态增长的最小二乘中的(X‘X)^-1(一种容易实现的方式,包含详细的证明)》

《C++矩阵求逆及求数据动态增长的最小二乘中的(X‘X)^-1(一种容易实现的方式,包含详细的证明)》

关于矩阵引理的一个简单的证明,帮我把他放在后面,或者你也可以在矩阵论的相关书籍中找到他的证明,下面将要给出的就是上面的代码的数学支持(注:I为单位阵):

《C++矩阵求逆及求数据动态增长的最小二乘中的(X‘X)^-1(一种容易实现的方式,包含详细的证明)》《C++矩阵求逆及求数据动态增长的最小二乘中的(X‘X)^-1(一种容易实现的方式,包含详细的证明)》

最后的迭代公式就是上面的程序的核心,这里就得到了矩阵的逆,此外这里利用这个思想我们还可以对最小二乘中矩阵运算有所改进:

《C++矩阵求逆及求数据动态增长的最小二乘中的(X‘X)^-1(一种容易实现的方式,包含详细的证明)》

也就是说我们不需要重新计算,只需要再一次的迭代的就可以得到新的矩阵了,时间复杂度为O(n*n).

到这里我已经证明的结果就全都放在这里了,下面是我的一些想法,也是写这篇文章的另一个目的:与诸位探讨一下可能的性能提升以及其他一些问题。

首先,可以看到的是无论是一般的矩阵求逆还是最小二乘中,每次迭代都会涉及到一个矩阵与一个向量的乘积也即是上面的《C++矩阵求逆及求数据动态增长的最小二乘中的(X‘X)^-1(一种容易实现的方式,包含详细的证明)》,同时还有就是一个列向量与一个行向量相乘得到一个矩阵的情形,两者的时间复杂度都为O(n*n),但是两者的计算过程都很简单,尤其是后者,所以一个问题就是可否利用加速的矩阵乘法对上述算法进行加速?

其次,如果上面的问题解决了那么,对矩阵的乘法的算法也会有相应的影响,因为:

《C++矩阵求逆及求数据动态增长的最小二乘中的(X‘X)^-1(一种容易实现的方式,包含详细的证明)》

对应于第一种情况。而:

《C++矩阵求逆及求数据动态增长的最小二乘中的(X‘X)^-1(一种容易实现的方式,包含详细的证明)》

对应第二种情况,所以无论哪一种情况解决了,都会对矩阵的乘法有提高,当然这也是我一直纠结的问题,而且从这里可以看出,上面的算法的时间复杂度其实与矩阵乘法的时间复杂度相当。

最后还有一点猜想,学过矩阵论的人都知道,其实矩阵的逆与矩阵的特征值是有很大的关系的,所以这里猜想这对矩阵特征值的计算会不会有帮助,当然这里说了只是猜想,依据不足仅供参考。

好了!写到这里该总结了:这里的结果其实有点像是习题的结果,完全是自己一个人在下面慢慢弄的,当然我也不知道,是否有其他人证明过相关的结论,因为我目前尚未见到,所以才写了这篇文章,倘若已经有人证明过得话,那就当我是孤陋寡闻,班门弄斧了。不过这里我保证是我原创的。

附录:矩阵引理+测试代码:

《C++矩阵求逆及求数据动态增长的最小二乘中的(X‘X)^-1(一种容易实现的方式,包含详细的证明)》

《C++矩阵求逆及求数据动态增长的最小二乘中的(X‘X)^-1(一种容易实现的方式,包含详细的证明)》

测试代码:

#ifndef    INVMATRIX_CPP
#define    INVMATRIX_CPP

#ifndef    INCLUDE_OPENCV2_CORE
#define    INCLUDE_OPENCV2_CORE
#include   "opencv2/core/core.hpp"
#endif
#include   <iostream>
#include   <iomanip>
#include   <cmath>

using namespace std;
using namespace cv;

void printMat(Mat & m,int precision)// 打印矩阵到屏幕上,且指明精度
{
	Size sz=m.size();
	int i,j;
	for(i=0;i<sz.height;i++)
	{
		for(j=0;j<sz.width;j++)
		{
			if(abs(m.at<double>(i,j)<1e-12)) // 将绝对值小于e-12 的数视为0,可自行决定
			  cout<<0<<"  ";
			else
			  cout<<setprecision(precision)<<m.at<double>(i,j)<<"  ";
		}
		cout<<endl;
	}
}

bool invMat(Mat origMat, Mat& invert)// 这就是上面的算法
{
	Size sz=origMat.size();
	if(sz.height!=sz.width||origMat.depth()!=CV_64F) // 只处理double型的方阵
	  return 0; 
	double disc;//记录中间量,并作判别之用
	Mat I=Mat::eye(sz,origMat.type());
	Mat D=origMat-I;
	Mat temp;
	disc=1+D.at<double>(0,0);
	if(abs(disc)<1e-12) // 同上
	  return 0; // 矩阵不可逆
	else
	{
		Mat invCk(I-(1/disc)*(I.col(0)*D.row(0)));
		for(long step=1;step<sz.height;step++)
		{
			temp=D.row(step)*invCk;
			disc=1+temp.at<double>(0,step);
			if(abs(disc)<1e-12) // 同上
			  return 0;// 同上
			else
			  invCk=invCk-(invCk.col(step)*temp)/disc;
		}
		invert=invCk.clone(); 
	}
	return 1;
}

int main(int argc, char *argv[])
{
	RNG rng(1234567);
	Mat src(26,26,CV_64FC1),dst; //26是我自己选的,你可以用其他的数字
	rng.fill(src,RNG::UNIFORM,Scalar(-100),Scalar(100));//用 -100~100之间的数随机填充矩阵 
	cout<<" 原始矩阵:"<<endl;
		printMat(src,4);
	if(!invMat(src,dst))
	  cout<<"错误!矩阵不可逆"<<endl;
	else
	{
		cout<<"其逆矩阵为:"<<endl;
			printMat(dst,4);
	}
	Mat re=src*dst;
	cout<<"看看两者的乘积:"<<endl;
		printMat(re,4);
	return 0;
}

#endif    /*INVMATRIX_CPP*/

    原文作者:u010842413
    原文地址: https://blog.csdn.net/u010842413/article/details/21122043
    本文转自网络文章,转载此文章仅为分享知识,如有侵权,请联系博主进行删除。
点赞