本章介绍了LU分解法，以及如何利用LU分解法求逆、行列式，针对每个公式、原理、代码进行了详细介绍，希望可以给大家带来帮助。
LU分解法与高斯法求逆一样，可以进行较高维数的矩阵运算（可计算万维及以上，但是精度不能保证，并且占有内存大，高维矩阵需要进行分块运算）
更改：
- 之前的代码使用_msize函数求矩阵维数，然而_msize函数只适用于windows，提供了一种可移植度高的版本。
- 之前在计算L、U的逆矩阵时，使用的是高斯求逆法，没有充分利用L、U矩阵的特性，已更改。
- 修改了选主元时的bug。
目录
LU分解法
概念
将系数矩阵A转变成等价两个矩阵L和U的乘积 ，其中L和U分别是单位下三角矩阵和上三角矩阵。当A的所有顺序主子式都不为0时，矩阵A可以唯一地分解为A=LU。其中L是下三角矩阵（主对角线元素为1），U是上三角矩阵。
于是，对矩阵A求逆就变成了：
因为LU分别为下三角矩阵和上三角矩阵，再进行高斯变换求逆矩阵时，浮点运算的复杂度将减少一半。
对矩阵A求其行列式的值变成：
因为L的主对角元素全1，故 A的行列式的值等于U主对角线元素的乘积。
确定L、U矩阵
因为矩阵L的主对角元素定死为1，因此可以通过矩阵乘法原理，逐行、逐列的逆推出矩阵L和U的元素
方法如下：
- 确定U的第一行
- 确定L的第一列
- 确定U的第二行
- 确定L的第二列
- ...
- 确定U的第n行
- 确定L的第n列
确定U的第一行：
- L的第一行与U的第一列对应元素相乘再相加，结果等于
，即：
- L的第一行与U的第二列对应元素相乘再相加，结果等于
，即：
- 不难看出，U的第一行与A的第一行元素相同：
确定L的第一列：
- L每列的首元素都为1，从第二个元素开始判断
- L的第二行与U的第一列对应元素相乘再相加，结果等于
，即：
- L的第三行与U的第一列对应元素相乘再相加，结果等于
，即：
- 故，L的第一列为：
仿照上述步骤，即可求出L、U，整理公式如下：
在进行代码设计时，求和部分可用for循环单独计算，存储到变量a中，变量a需初始化为0，当i=0时，不满足for循环条件，a=0，即可推出上述U的第一行和L的第一列。
LU分解法的推导与验证请参考相关数值分析教材，推荐参考博客矩阵的直接LU分解法
求L、U的逆矩阵
由于L、U矩阵自身的特性，这时我们可以通过L、U矩阵直接逆推出其逆矩阵，无需再用高斯求逆，逆推思想与上述思想类型，这里直接给出下三角矩阵和上三角矩阵求逆公式：
对于下三角矩阵：
对于上三角矩阵：
需要注意：上述公式是基于程序设计的思维，也就是我们对L、U矩阵本身进行改动，使其变成各自的逆矩阵，这样就无需再开辟新的内存空间。
- 计算下三角矩阵的逆时，我们是按行，从第1行到第n行开始计算，第一行第一列，第二行第一列，第二行第二列，以此类推。
- 计算上三角矩阵的逆时，我们是按列，从第n列到第1列开始计算，第n列第n行，第n列第n-1行，以此类推。
- 二者的计算次序正好相反，该点在代码设计中也有所体现。
LU分解法的意义
- LU具有承袭性，这是LU的优点。
- LU只适用于解所有顺序主子式都大于0的，通用性欠缺，这是LU的缺点。
- LU法不保证具有数值稳定性，这是LU的缺点。（Gauss法可以用选取列主元技巧保证数值稳定性）
集合LU与Gauss优点，同时规避掉这些缺点的，是LUP分解法。
作者：寨森Lambda-CDM
我个人的理解是：计算机在处理超维矩阵时（例如一万维），会采用分块矩阵的方式进行求解，通过先分块再LU，将矩阵分成一块块下三角矩阵和上三角矩阵
存储起来，在后续其他计算中，直接调用下三角矩阵和上三角矩阵计算的效率更高。（该部分我也在学习，欢迎大家讨论）
程序设计
该部分程序实际上是LUP分解法，增加了选择主元的过程，因为在分解LU以及高斯消元求逆时，如果主元的元素是0，那么计算机将无法计算，所以在分解前需先选择主元。
在矩阵设计时，最开始使用的是_msize函数通过内存计算矩阵的维数，而在后来使用Linux和Mac系统后，发现不支持_msize函数，同时又不想每次计算时都输入矩阵的维数，因此设计一个结构体来存储矩阵。（_msize函数只支持windows，我的其他博客有相关介绍：判断矩阵维数）
整体流程如下：
- 检查是否符合计算要求
- 保护原矩阵，进行临时矩阵的拷贝
- 选择主元
- 构建L、U矩阵
- 计算L、U矩阵的逆（如果求行列式，无需求逆，直接求U矩阵主对角线的乘积）
- 利用矩阵乘法求得最终结果
- 释放内存
LUP求逆
1）结构体定义与头文件
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <time.h>
#include <math.h>
typedef struct Matrix
{
int row;
int col;
double **data;
} Matrix;2）辅助函数
// 创建矩阵
Matrix MakeMatrix(int row, int col)
{
int i = 0;
Matrix arr = {0};
arr.row = row;
arr.col = col;
arr.data = (double **)malloc(sizeof(double *) * arr.row);
if (arr.data == NULL)
exit(-1);
for (i = 0; i < arr.row; i++)
{
arr.data[i] = (double *)malloc(sizeof(double) * arr.col);
memset(arr.data[i], 0, sizeof(double) * arr.col);
}
return arr;
}
// 释放矩阵内存
void free_Matrix(Matrix src)
{
int i;
for (i = 0; i < src.row; i++)
{
free(src.data[i]);
}
free(src.data);
}
// 打印矩阵
void print_Matrix(Matrix arr)
{
int i, j;
putchar('\n');
for (i = 0; i < arr.row; i++)
{
for (j = 0; j < arr.col; j++)
{
printf("%10.4lf ", arr.data[i][j]);
}
putchar('\n');
}
}
// 生成随机矩阵
void rand_Matrix(Matrix arr)
{
int i, j;
double num, res, res1, res2 = 0;
srand((unsigned int)time(NULL));
for (i = 0; i < arr.row; i++)
{
for (j = 0; j < arr.col; j++)
{
//随机生成整数
/* num = pow(-1, rand() % 2 + 1) * (rand() % 20);
arr.data[i][j] = num; */
//4位浮点数
res1 = res2 = 0;
while (res1 == 0)
{
res1 = pow(-1, rand() % 2 + 1) * (rand() % 101);
res2 = pow(-1, rand() % 2 + 1) * (rand() % (int)(1e4));
}
res = res1 + (res2 * 1e-4);
arr.data[i][j] = res;
}
}
}上述辅助函数需要注意两点：
- 释放内存时，先便利释放每一行的内存，最终将整个data内存释放掉。
- 创建随机矩阵时，srand函数用于修饰rand，以time为标准生成随机数。
- 这里创建浮点数的思想为：1）创建一个范围在[-100,100]内的整数。2）创建一个范围在[0,1000)内的整数，然后对其除以10的-4次方，使其变成小数部分。3）将二者相加，即得到了小数点后4位的随机数。（创建随机矩阵时，笔者之前使用的是将两个整数做除法，但该种做法无法与matlab结果验证）
3）矩阵乘法
矩阵乘法较为简单，笔者在C语言矩阵乘法中介绍了其原理，以及多种程序设计思路，这里不再赘述。
Matrix Matrix_Mul(const Matrix left, const Matrix right)
{
if (left.col != right.row)
exit(-1); // 判断左列是否等于右行
int i, j, k;
Matrix res = MakeMatrix(left.row, left.row);
for (i = 0; i < left.row; i++)
{
for (j = 0; j < right.col; j++)
{
for (k = 0; k < left.col; k++)
{
res.data[i][j] += left.data[i][k] * right.data[k][j];
}
}
}
return res;
}4） LU求逆
Matrix Matrix_LU_Inv(Matrix src)
{
assert(src.row == src.col);
int i, j, k, principal;
Matrix L, U, tmp, res;
double Lsum, Usum, p, Max;
L = MakeMatrix(src.row, src.col);
U = MakeMatrix(src.row, src.col);
tmp = MakeMatrix(src.row, src.col);
for (i = 0; i < src.row; i++)
{
memcpy(tmp.data[i], src.data[i], sizeof(src.data[0][0]) * src.col);
}
// 初始化斜对角线
for (i = 0; i < L.row; i++)
{
for (j = 0; j < L.col; j++)
{
if (i == j)
L.data[i][j] = 1; // L主对角线为1
}
}
// 选主元
for (j = 0; j < tmp.col; j++)
{
// 默认当前主元最大
principal = j;
Max = fabs(tmp.data[principal][j]); // 用绝对值比较
// 只遍历当前主元下方的元素
for (i = j + 1; i < tmp.row; i++)
{
if (fabs(tmp.data[i][j]) > Max)
{
principal = i;
Max = fabs(tmp.data[i][j]);
}
}
if (j != principal)
{
for (k = 0; k < tmp.col; k++)
{
p = tmp.data[principal][k];
tmp.data[principal][k] = tmp.data[j][k]; // 交换两行
tmp.data[j][k] = p;
}
}
}
// 计算L、U
for (i = 0; i < tmp.row; i++)
{
for (j = 0; j < tmp.col; j++)
{
if (i < 1 && j < 1)
{
U.data[i][j] = tmp.data[i][j];
}
else if (i <= j)
{
Usum = 0;
for (k = 0; k < i; k++)
{
Usum += L.data[i][k] * U.data[k][j]; // 计算求和部分
}
U.data[i][j] = tmp.data[i][j] - Usum;
}
else
{
Lsum = 0;
for (k = 0; k < j; k++)
{
Lsum += L.data[i][k] * U.data[k][j]; // 计算求和部分
}
L.data[i][j] = (tmp.data[i][j] - Lsum) / U.data[j][j];
}
}
}
free_Matrix(tmp);
//L求逆
for (i = 0; i < L.row; i++)
{
for (j = 0; j <= i; j++)
{
if (i == j)
{
L.data[i][j] = 1.0 / L.data[i][j];
}
else
{
Lsum = 0;
for (k = j; k < i; k++)
{
Lsum += L.data[i][k] * L.data[k][j];
}
L.data[i][j] = -1.0 * Lsum / L.data[i][i];
}
}
}
// U求逆
for (j = U.col - 1; j >= 0; j--)
{
for (i = j; i >= 0; i--)
{
if (i == j)
{
U.data[i][j] = 1.0 / U.data[i][j];
}
else
{
Usum = 0;
for (k = j; k >= i + 1; k--)
{
Usum += U.data[i][k] * U.data[k][j];
}
U.data[i][j] = -1.0 * Usum / U.data[i][i];
}
}
}
res = Matrix_Mul(U, L); // 乘法
free_Matrix(L), free_Matrix(U);
return res;
}测试
void LU_inv_test()
{
Matrix arr = MakeMatrix(3, 3);
rand_Matrix(arr);
print_Matrix(arr);
Matrix res = Matrix_LU_Inv(arr);
print_Matrix(res);
Matrix res2 = Matrix_Guass_inver(arr);
print_Matrix(res2);
}
int main()
{
LU_inv_test();
return 0;
}整数矩阵计算：
C语言计算结果：
Matlab计算结果：
显然，二者计算的结果是等价的。
浮点矩阵计算：
C语言计算结果：
Matlab计算结果：
显然，二者也是等价的。有兴趣的同学可以尝试创建一万维的矩阵进行运算，该方法也是能够计算的，当然会牺牲一些内存。
- 矩阵求逆时只涉及到简单的加减乘除运算，运算结果与matlab结果并无差异。而涉及到复杂的浮点数的运算，C语言的计算精度稍显不足。例如：C语言中sqrt开方的处理，与Fortran中的处理就略有差别。对于精度较高的运算，笔者认为使用Fortran是个不错的选择，或者使用各大机构发布的数值运算库（例如Intel OneAPI的mkl库）。
- 对于极小或者极大的数，可以考虑加入缩放因子来提高精度。在计算LU矩阵前对原矩阵进行放大，每个元素乘以缩放因子beta，得到逆矩阵同样乘以beta即可。（本文算例不需要加入缩放因子，故所以上文就没有提及)
LUP求行列式
行列式的求法相对简单，依据前文的公式，最后只需要对U矩阵对主对角线求乘积即可。这里给出代码，各位同学请自行测试。
#include "matrix.h"
double Matrix_LU_Det(Matrix src)
{
assert(src.row == src.col);
int i, j, k, principal;
Matrix L, U, tmp, res;
double Lsum, Usum, p, Max;
L = MakeMatrix(src.row, src.col);
U = MakeMatrix(src.row, src.col);
tmp = MakeMatrix(src.row, src.col);
for (i = 0; i < src.row; i++)
{
memcpy(tmp.data[i], src.data[i], sizeof(src.data[0][0]) * src.col);
}
// 初始化斜对角线
for (i = 0; i < L.row; i++)
{
for (j = 0; j < L.col; j++)
{
if (i == j)
L.data[i][j] = 1; // L主对角线为1
}
}
// 选主元
for (j = 0; j < tmp.col; j++)
{
// 默认当前主元最大
principal = j;
Max = fabs(tmp.data[principal][j]); // 用绝对值比较
// 只遍历当前主元下方的元素
for (i = j + 1; i < tmp.row; i++)
{
if (fabs(tmp.data[i][j]) > Max)
{
principal = i;
Max = fabs(tmp.data[i][j]);
}
}
if (j != principal)
{
for (k = 0; k < tmp.col; k++)
{
p = tmp.data[principal][k];
tmp.data[principal][k] = tmp.data[j][k]; // 交换两行
tmp.data[j][k] = p;
}
}
}
// 计算L、U
for (i = 0; i < tmp.row; i++)
{
for (j = 0; j < tmp.col; j++)
{
if (i < 1 && j < 1)
{
U.data[i][j] = tmp.data[i][j];
}
else if (i <= j)
{
Usum = 0;
for (k = 0; k < i; k++)
{
Usum += L.data[i][k] * U.data[k][j]; // 计算求和部分
}
U.data[i][j] = tmp.data[i][j] - Usum;
}
else
{
Lsum = 0;
for (k = 0; k < j; k++)
{
Lsum += L.data[i][k] * U.data[k][j]; // 计算求和部分
}
L.data[i][j] = (tmp.data[i][j] - Lsum) / U.data[j][j];
}
}
}
free_Matrix(tmp);
p = 1;
for (i = 0; i < U.row; i++)
{
for (j = 0; j < U.col; j++)
{
if (i == j)
p *= U.data[i][j];
}
}
free_Matrix(L), free_Matrix(U);
return p;
}结语：
通过上述的实验，我们发现单纯C语言设计矩阵运算时有些劣势，代码有些冗余而且不易阅读，这些问题将在C++中得到解决，以利用C++将上述结构体与函数统一归为类，然后利用函数重载，优化代码的可阅读性。
7766
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
