对称矩阵的RtDR分解（LDLt分解）C代码

最新推荐文章于 2024-06-11 16:14:15 发布

yu_xianguo

最新推荐文章于 2024-06-11 16:14:15 发布

阅读量3.7k

点赞数 5

分类专栏：编程文章标签： C代码数值计算

本文链接：https://blog.csdn.net/YU_Xianguo/article/details/84665547

版权

博主因求解线性方程组和矩阵特征值问题，发现OpenCV和Eigen库效率不理想，于是自学数值分析并编写C代码实现对称矩阵的RtDR分解。此方法适用于N<200的矩阵，分解结果D和R直接覆盖输入矩阵A的内存。文章提供Matlab和C代码示例，并讨论了矩阵存储优化、错误处理及大型矩阵的处理策略。

摘要由CSDN通过智能技术生成

最近碰到求解线性方程组以及求矩阵的特征值等问题，OpenCV自带的算法实在是太慢了，另外我还试了Eigen库，比OpenCv虽然快了一倍，但是比Matlab还是慢了一个量级不止。。。因此我决定自己编写几个程序以满足我的特定需要。这篇博文将给出一个对称矩阵的RtDR分解方法（书里面一般都是LDLt分解，我直接求得是转置，即R=L'）。

数值计算不同于基本数学理论，《线性代数》以及《高等数学基础》是工科的数学课程，里面介绍了很多一般线性代数问题的解法，但是那些只是在理论上可行，并没有考虑计算机存储的舍入误差的影响，如果照着课本上的思路来实现算法，则在计算效率和计算结果的精度上都与Matlab相差甚远。为此，我花了一些时间学习了两本数值分析的书，其中《Matrix Computations》(4th Edition)是非常经典的书了，写得也很赞。根据书中介绍的方法写了一些函数，在处理中小型矩阵（e.g.N在200以内）时，具有较快的速度和较高的精度。

本文将给出这样一个函数：输入一个实对称矩阵A，计算它的RtDR分解：A=R'*D*R，其中R是单位上三角矩阵，D是对角矩阵。该方法要求A的所有顺序主子式都是非奇异的。

这是Matlab代码 version#1：

function  [L,D] = LDLtDecomp(A)
% A = L*D*L'
% A is symmetric
% L is unit lower triangular
% D is diagonal
n = size(A,1);
v = zeros(n,1);
A(2:n,1) = A(2:n,1) / A(1,1);
for i = 2:n
    v(1:i-1) = A(i,1:i-1) .* diag(A(1:i-1,1:i-1))';
    A(i,i) = A(i,i) - A(i,1:i-1) * v(1:i-1);
    A(i+1:n,i) = (A(i+1:n,i) - A(i+1:n,1:i-1) * v(1:i-1)) / A(i,i);
end
D = diag(diag(A));
L = A;
for i = 1:n
    L(i,i