3D数学之矩阵的各种求逆-白红宇

3D数学之矩阵的各种求逆

阅读量：4041 次

发布时间：2019-05-24

本文共 5885 字，大约阅读时间需要 19 分钟。

经过三天的准备终于把矩阵的各种求逆方法以及代码完成了。心里有点小激动，come on，来吧，点燃你的心中的那团火，跟着游戏音乐的律动一起跟我走入神秘的3D世界。

下面介绍三种方法：

1.用伴随矩阵求逆

2.用高斯－约当消元法求逆

3.用LU分解求逆（听别人说效率最高）

在介绍方法之前，先阐述几个概念

第一种方法：伴随矩阵

为了求A的逆，我们需要求A*和det(A)

为了求A*，我们需要求C

所以我们需要一个求代数余子式的函数，从而实现一个求代数余子式矩阵的函数

由于求代数余子式中需要求余子阵以及其行列式的值。所以我们需要一个函数实现求余子阵，一个求矩阵的行列式的函数。

为了求det(A),我们也需要一个求矩阵行列式的函数。而行列式的求法又牵涉到代数余子式的求法，所以求行列式的函数与求代数余子式的函数相互嵌套调用，最终出口在求行列式的函数里，即当矩阵只有一个元素的时候。

下面提供实现代码：

matrix的定义如下：

#ifndef MATRIX_H#define MATRIX_H#include 
   
    namespace Seamanj{	class Matrix	{	protected:		int Nrows,Ncols;		double *data;	public:		Matrix(int rows = 0,int cols = 0,const double *M = NULL);		Matrix(Matrix& m)//深度复制		{			Nrows=m.Nrows;			Ncols=m.Ncols;			data=new double[Nrows*Ncols];			for(int i=0;i

求余子阵函数的实现

Matrix Matrix::getSubmatrixByDeleting_i_row_j_col(int i,int j)	{		if( i < 0 || j < 0 || i >= Nrows || j >= Ncols )		{			cerr << "getSubmatrixByDeleting_i_row_j_col error!"<
   
    = j)					sm(ii,jj)=(*this)(ii,jj+1);				else if(ii >= i && jj < j)					sm(ii,jj)=(*this)(ii+1,jj);				else if(ii >= i && jj >= j)					sm(ii,jj)=(*this)(ii+1,jj+1);			}		return sm;	}

求代数余子式函数的实现

double Matrix::getCofactor_i_j(int i,int j)	{		if( i < 0 || j < 0 || i >= Nrows || j >= Ncols || Nrows != Ncols)		{			cerr << "getCofactor_i_j error"<

求代数余子式矩阵的实现代码

Matrix Matrix::getCofactorMatrix()	{		Matrix cm(Nrows,Ncols);		for(int i=0;i

求矩阵行列式函数的实现代码

double Matrix::getDetermination()	{		if( Nrows < 1 || Ncols < 1|| Nrows != Ncols)		{			cerr << "getDetermination error"<

最终还需要一个求转置函数，将代数余子式矩阵转置为伴随矩阵

Matrix Matrix::getTranspose()	{		if( Nrows < 1 || Ncols < 1)		{			cerr << "getTranspose error"<

这样用伴随矩阵求逆的过程可以表示如下：

Matrix Matrix::getInverseByAdjugate()	{		if(Nrows != Ncols || Nrows < 1)		{			cerr << "getInverseByAdjugate error"<
   
    -1e-6)		{			cerr<<"This matrix isn't invertible"<

第二种方法：高斯－约当消元法

来个小插曲，为什么叫高斯－约当消元法（Gauss_Jordan_Elimination）呢？

那是因为高斯在消元的时候只想到了往下消元，没有往上消，而约当却想到了往上消，所以叫Gauss_Jordan_Elimination。

对应的代码如下：

Matrix Matrix::getInverseByGauss_Jordan_Elimination()	{		if(Nrows != Ncols || Nrows < 1)		{			cerr << "getgetInverseByGauss_Jordan_Elimination error"<
   
    max)				{					max=abs(m(i,j));					maxrow=i;				}			}			if(maxrow == -1)			{				cerr<<"This matrix isn't invertible"<
    
     =j;k--)//与主元同行的元素（从主元到行末，因为主元前面是0元素，所以没有必要除了）除以主元，注意是从后往前，这样可以保持主元除自己之前保持不变			{				m(j,k)/=m(j,j);			}			for(int k=m.Ncols-1;k>=j;k--)//对每行从后往前进行消元			{				for(int l=0;l

第三种利用LU分解求逆

LU分解函数的实现

Matrix LU_Decompose(Matrix& m, int *index,double* detSign)	{		double *rowNormalizer=new double[m.Nrows];		Matrix result(m);		double exchangeParity = 1.0;		//calulate a normalizer for each row		for(int i=0;i
   
    maxvalue)					maxvalue=abs(result(i,j));			}			if(maxvalue < 1e-6 && maxvalue > -1e-6)			{				cerr<<"This matrix isn't invertible"<
    
     maxvalue)				{					maxvalue=sum;					pivotRow=i;				}			}			if(pivotRow == -1)			{				cerr<<"This matrix isn't invertible"<

LU回带函数的实现

void LU_back_substitution(Matrix& m, int *index, double col[],double result_col[])	{		for(int i=0;i
   
    =0;i--)		{			double sum=result_col[i];			for(int j=i+1;j

最终通过LU分解求逆的函数如下：

Matrix Matrix::getInverseByLU_Decompose()	{		if(Nrows != Ncols || Nrows < 1)		{			cerr << "getInverseByLU_Decompose error"<

下列通过这三种方法进行验证上面那个高斯消元法介绍的矩阵的求逆

void main(){	double m_try[9] = {1,2,3,2,4,5,3,5,6};	Seamanj::Matrix myMatrix(3,3,m_try);	cout << "Origin Matrix:" << endl << myMatrix << endl ;	myMatrix=myMatrix.getInverseByAdjugate();	cout << "After inversing ByAdjugate:"<
   
    <

结果如下：

最后给出该程序的头文件和源文件

Matrix.h头文件内容：

#ifndef MATRIX_H#define MATRIX_H#include 
   
    namespace Seamanj{	class Matrix	{	protected:		int Nrows,Ncols;		double *data;	public:		Matrix(int rows = 0,int cols = 0,const double *M = NULL);		Matrix(Matrix& m)//深度复制		{			Nrows=m.Nrows;			Ncols=m.Ncols;			data=new double[Nrows*Ncols];			for(int i=0;i

Matrix.cpp源文件内容：

#include "Matrix.h"#include 
   
    #include 
    
     using namespace std;namespace Seamanj{	Matrix::Matrix(int rows,int cols,const double *M)	{		if(rows < 0 || cols < 0 || ((rows == 0 || cols == 0) && cols != rows))//注意rows和cols可以同时为0		{			cerr << "Matrix size impossible (negative or zero rows or cols)" << endl;			exit(1);		}		Nrows=rows;		Ncols=cols;		data=new double[Nrows*Ncols];		if(!data)		{			cerr << "not enough memory to allocate " << rows << " x " << cols << " Matrix" << endl;			exit(1);		}		if(M)		{			for(int i=0,k=0;i
     
      = Nrows || j < 0 || j >= Ncols)		{			cerr << " Invalid entry for Matrix " << endl;			exit(0);		}		return data[i*Ncols+j];	}	Matrix Matrix::getInverseByGauss_Jordan_Elimination()	{		if(Nrows != Ncols || Nrows < 1)		{			cerr << "getgetInverseByGauss_Jordan_Elimination error"<
      
       max)				{					max=abs(m(i,j));					maxrow=i;				}			}			if(maxrow == -1)			{				cerr<<"This matrix isn't invertible"<
       
        =j;k--)//与主元同行的元素（从主元到行末，因为主元前面是0元素，所以没有必要除了）除以主元，注意是从后往前，这样可以保持主元除自己之前保持不变			{				m(j,k)/=m(j,j);			}			for(int k=m.Ncols-1;k>=j;k--)//对每行从后往前进行消元			{				for(int l=0;l
        
         -1e-6) { cerr<<"This matrix isn't invertible"<
         
          = Nrows || j >= Ncols || Nrows != Ncols) { cerr << "getCofactor_i_j error"<
          
           = Nrows || j >= Ncols ) { cerr << "getSubmatrixByDeleting_i_row_j_col error!"<
           
            = j) sm(ii,jj)=(*this)(ii,jj+1); else if(ii >= i && jj < j) sm(ii,jj)=(*this)(ii+1,jj); else if(ii >= i && jj >= j) sm(ii,jj)=(*this)(ii+1,jj+1); } return sm; } double Matrix::getDetermination() { if( Nrows < 1 || Ncols < 1|| Nrows != Ncols) { cerr << "getDetermination error"<
            
             maxvalue) maxvalue=abs(result(i,j)); } if(maxvalue < 1e-6 && maxvalue > -1e-6) { cerr<<"This matrix isn't invertible"<
             
              maxvalue) { maxvalue=sum; pivotRow=i; } } if(pivotRow == -1) { cerr<<"This matrix isn't invertible"<
              
               =0;i--) { double sum=result_col[i]; for(int j=i+1;j
               
                -1e-6)?0.0:m(i,j))<<'\t'; os<