网上大部分C#实现矩阵求逆都比较复杂,现在在这里分享一种很好理解的矩阵求逆方法,而且可以适用于任何形式的可逆矩阵求逆,但是肯定运行效率不如其它的算法,正所谓鱼和熊掌不可兼得。
我们采用的是通过单位矩阵变换的这种方法来实现的,话不多说,下面解释实现原理。
将需要变化的矩阵与单位矩阵拼在一起形成增广矩阵。
A为需要求逆的矩阵,E为单位矩阵。
如图
然后我们经过初等行列式变换,将增广矩阵左半部分变为单位矩阵,那么右半部分我们就可以得到A的逆矩阵,具体原理可以自行百度,这里就不再阐述原理了。
虽然我们通过手算可以轻易的进行初等行变换,知道通过调换两行来让计算更加简便,知道如何才能让除对角线元素都变为0,但是计算机想要实现我们灵活的思想似乎不太容易,所以此时我们想是否可以通过一种通用的行列式变换来达到目的。
我的想法是这样的,先让需要求逆的矩阵的左下部分变为零,之后再让右上部分变为零,这样似乎就可以实现我们的求逆了。
如图,按照以下逻辑进行变换,更高阶数的矩阵也可以使用同样的逻辑进行变换。
文章来源地址https://uudwc.com/A/Mx3vn
此时我们发现右上角也可以用同样的逻辑进行变换,只不过用一个与之前相反的循环就可以了。
代码实现
//input_matrix应该为一个二维数组,如何创建二维数组这里不再阐述,请先将input_matrix变为二维数组
int row = input_matrix.GetLength(0);
int col = input_matrix.GetLength(1); //获取数组长度
double[,] temp_matrix = new double[1, col]; //创建一个临时数组,用以承接变量
//消除左下元素
for (int j = 0; j < col; j++)
{
for (int i = 0; i < row; i++) //两个循环来遍历整个数组
{
if (i == j) //为了消除左下半元素,做一个元素位置的过滤,右上半元素不看
{
for (int item_i = i + 1; item_i < row; item_i++) //循环行
{
if (item_i == row)
{
continue; //边界判断,如果到达矩阵行边界,便跳出循环
}
else if (input_matrix[item_i, j] != 0) //如果下一行元素不为零,就开始进行行变换
{
double temp_value = input_matrix[item_i, j] / input_matrix[i, j]; //得到变换的倍数
for (int item_j = j; item_j < col; item_j++) //对同一行的每一列进行变换
{
input_matrix[item_i, item_j] -= temp_value * input_matrix[i, item_j]; //对每一列乘以倍数以便进行行变换
}
else
{
//没有意义,单纯为了if语句的完整
}
}
}
}
}
//消除右上元素
for (int j = col - 1; j >= 0; j--)
{
for (int i = row - 1; i >= 0; i--) //两个循环遍历数组
{
if (i == j) //为了消除右上元素,做一个元素过滤,左下角元素不看
{
for (int item_i = i - 1; item_i >= 0; item_i--) //开始从右下角变换
{
if (item_i == -1)
{
continue; //边界判断,到达矩阵行边界,跳出
}
else if (input_matrix[item_i, j] != 0)
{
double temp_value = input_matrix[item_i, j] / input_matrix[i, j]; //得到变换的倍数
for (int item_j = j; item_j >= 0; item_j--) //对同一行的非零列进行变换,
{
input_matrix[item_i, item_j] -= temp_value * input_matrix[i, item_j]; //对每一列乘以倍数以便进行行变换
}
else
{
//无意义,只是为了让if语句完整
}
}
}
}
}
经过我们的初等行变换,我们已经将除对角线的元素全部变为0了,此时只需要将对角线元素都变为1即可
//开始单位化输入的矩阵
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++) //遍历数组
{
if (i == j) //对对角线元素进行单位化
{
double temp_value = 1 / input_matrix[i, j];
input_matrix[i, j] *= temp_value; //单位化成功
我们此时基本上完成了将矩阵变为单位矩阵了,但是还有一个问题,如果对角线的头一个元素为0或者对角线中间的某个元素待相加时为0,那么这个算法就寄了。
比如说这个矩阵,当程序用第一行加到第二行时,由于(1,1)为0,后续的计算都会异常,那么如何改变这种状况呢?我们想到可以把对角线带有0的行全部换到最底部的行,这样在相加的时候就不会出现以上的情况,除非你的矩阵本身就不可逆。
//将对角线中的0值都放到末尾
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++) //遍历数组
{
if (i == j) //得到对角元素
{
if (i == row - 1)
{
continue; //检测边界,到了边界便不再执行以下语句
}
else if (input_matrix[i, j] == 0) //找到对角矩阵为零的元素
{
for (int temp_j = 0; temp_j < col; temp_j++)
{
temp_matrix[0, temp_j] = input_matrix[i, temp_j];
}
for (int temp_j = 0; temp_j < col; temp_j++)
{
input_matrix[i, temp_j] = input_matrix[i + 1, temp_j];
}
for (int temp_j = 0; temp_j < col; temp_j++)
{
input_matrix[i + 1, temp_j] = temp_matrix[0, temp_j];
}
//3个for循环的功能是将对角线为零的元素与下一列对换
}
else
{
//无意义,只是为了if语句完整
}
}
}
}
文章来源:https://uudwc.com/A/Mx3vn
此时我们已经基本完成了程序的部分,此时我们只需要创建一个与可逆矩阵同行列的单位矩阵跟着变化即可。
以下为完整代码
public static double[,] matrix_inv(double[,] input_matrix)
{
int row = input_matrix.GetLength(0);
int col = input_matrix.GetLength(1); //获取数组长度
double[,] unit_matrix = new double[row, col]; //定义一个数组,作为单位矩阵
double[,] temp_matrix = new double[1, col]; //创建一个临时数组,用于传递参数
double[,] temp_unit_matrix = new double[1, col]; //创建一个临时单位数组,用于传递参数
for (int i = 0; i < row; i++) //创建单位矩阵
{
for (int j = 0; j < col; j++)
{
if (i == j)
{
unit_matrix[i, j] = 1;
}
else
{
unit_matrix[i, j] = 0;
}
}
} //得到单位矩阵
//因本程序逆矩阵的算法需要对角线的元素不为零,故要对对角线为零的元素进行行变换
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++) //遍历数组
{
if (i == j) //得到对角元素
{
if (i == row - 1)
{
continue; //检测边界,到了边界便不再执行以下语句
}
else if (input_matrix[i, j] == 0) //找到对角矩阵为零的元素
{
for (int temp_j = 0; temp_j < col; temp_j++)
{
temp_matrix[0, temp_j] = input_matrix[i, temp_j];
temp_unit_matrix[0, temp_j] = unit_matrix[i, temp_j];
}
for (int temp_j = 0; temp_j < col; temp_j++)
{
input_matrix[i, temp_j] = input_matrix[i + 1, temp_j];
unit_matrix[i, temp_j] = unit_matrix[i + 1, temp_j];
}
for (int temp_j = 0; temp_j < col; temp_j++)
{
input_matrix[i + 1, temp_j] = temp_matrix[0, temp_j];
unit_matrix[i + 1, temp_j] = temp_unit_matrix[0, temp_j];
}
//3个for循环的功能是将对角线为零的元素与下一列对换
//当然单位矩阵也得跟着变换
}
else
{
//无意义,只是为了if语句完整
}
}
}
}
//消除左下元素
for (int j = 0; j < col; j++)
{
for (int i = 0; i < row; i++) //两个循环来遍历整个数组
{
if (i == j) //为了消除左下半元素,做一个元素位置的过滤,右上半元素不看
{
for (int item_i = i + 1; item_i < row; item_i++) //循环行
{
if (item_i == row)
{
continue; //边界判断,如果到达矩阵行边界,便跳出循环
}
else if (input_matrix[item_i, j] != 0) //如果下一行元素不为零,就开始进行行变换
{
double temp_value = input_matrix[item_i, j] / input_matrix[i, j]; //得到变换的倍数
for (int item_j = j; item_j < col; item_j++) //对同一行的每一列进行变换
{
input_matrix[item_i, item_j] -= temp_value * input_matrix[i, item_j]; //对每一列乘以倍数以便进行行变换
}
for (int unit_j = 0; unit_j < col; unit_j++)
{
unit_matrix[item_i, unit_j] -= temp_value * unit_matrix[i, unit_j]; //单位矩阵跟着变化
}
}
else
{
//没有意义,单纯为了if语句的完整
}
}
}
}
}
//消除右上元素
for (int j = col - 1; j >= 0; j--)
{
for (int i = row - 1; i >= 0; i--) //两个循环遍历数组
{
if (i == j) //为了消除右上元素,做一个元素过滤,左下角元素不看
{
for (int item_i = i - 1; item_i >= 0; item_i--) //开始从右下角变换
{
if (item_i == -1)
{
continue; //边界判断,到达矩阵行边界,跳出
}
else if (input_matrix[item_i, j] != 0)
{
double temp_value = input_matrix[item_i, j] / input_matrix[i, j]; //得到变换的倍数
for (int item_j = j; item_j >= 0; item_j--) //对同一行的非零列进行变换,
{
input_matrix[item_i, item_j] -= temp_value * input_matrix[i, item_j]; //对每一列乘以倍数以便进行行变换
}
for (int unit_j = 0; unit_j < col; unit_j++)
{
unit_matrix[item_i, unit_j] -= temp_value * unit_matrix[i, unit_j]; //单位矩阵跟着变化
}
}
else
{
//无意义,只是为了让if语句完整
}
}
}
}
}
//开始单位化输入的矩阵
for (int i = 0; i < row; i++)
{
for (int j = 0; j < col; j++) //遍历数组
{
if (i == j) //对对角线元素进行单位化
{
double temp_value = 1 / input_matrix[i, j];
input_matrix[i, j] *= temp_value; //单位化成功
for (int unit_j = 0; unit_j < col; unit_j++)
{
unit_matrix[i, unit_j] *= temp_value; //单位矩阵要乘以单位化的系数
}
}
}
}
return unit_matrix;
}
本人能力有限,欢迎各位大佬给出更优的建议。