matlab如何实现LDLT分解（改进的平方根法）

matlab如何实现LDLT分解（改进的平方根法）

• 前言
• 一、思路说明
• 二、LDLT分解的函数（.m文件后缀）
• 总结

  ---

  前言

  提示：本篇博客仅提供函数代码，不进行数学推导！需要了解数学推导的uu请移步其他博客！

  笔者这学期选修了数值计算这门课程。一次实验报告中需要自行编写LDLT分解的函数，并在脚本文件中调用此函数，返回求得的x、L、D。

  可能由于这个方法较为冷门？（瞎猜的，别cue我）网上相关的资料很少，能找到的、不用开vip就能看的文章只有个位数。。。因此能参考的内容实在太少。

  求人不如求己！ 自！己！来！

  所以笔者的灵感来自于老师已经给出的Cholesky分解、LU分解的函数文件上进行修改，最终完成了LDLT方法的求解函数。（事实上，老师布置这道题的初衷就是希望我们能在他给出的函数文件上进行修改，进而得到LDLT分解的函数）

  秉着乐于助人的原则（always），遂把此代码分享在各个平台。愿在世界各地因同样问题正在挠头的你们，能看到这篇博客吧。

  ps：笔者是刚入门matlab和数值计算的菜鸟，若以下代码有更加简略的写法，欢迎评论区留言或私信我~

  ---

  提示：以下是本篇文章正文内容

  一、思路说明

  1.假设我们要使用LDLT分解解如下一个方程组：

  {   4 x 1 − 2 x 2 − 4 x 3 = 10   − 2 x 1 + 17 x 2 + 10 x 3 = 3   − 4 x 1 + 10 x 2 + 9 x 3 = − 7 \begin{cases} \ 4x_1- 2x_2-4x_3=10\\ \ -2x_1+17x_2+10x_3=3\\ \ -4x_1+10x_2+9x_3=-7 \end{cases} ⎩ ⎨ ⎧​ 4x1​−2x2​−4x3​=10 −2x1​+17x2​+10x3​=3 −4x1​+10x2​+9x3​=−7​

  2.那么系数矩阵为A，A是如下的一个三阶对称矩阵：

  [ 4 − 2 − 4 − 2 17 10 − 4 10 9 ] \left[ \begin{matrix} 4 & -2 & -4 \\ -2 & 17 & 10 \\ -4 & 10 & 9 \end{matrix} \right] ​4−2−4​−21710​−4109​ ​经过初等变换后，A变为上三角矩阵DU。我们仅取DU的对角线元素，并赋值给全零阵D。

  则DU如下：

  [ 4 − 2 − 4 0 16 8 0 0 1 ] \left[ \begin{matrix} 4 & -2 & -4 \\ 0 & 16 & 8 \\ 0 & 0 & 1 \end{matrix} \right] ​400​−2160​−481​ ​

  所以根据上面的规则，矩阵D如下：

  [ 4 0 0 0 16 0 0 0 1 ] \left[ \begin{matrix} 4 & 0 & 0 \\ 0 & 16 & 0 \\ 0 & 0 & 1 \end{matrix} \right] ​400​0160​001​ ​

  至于D为什么长这个样子，大家去看LDLT的数学推导就好，此处不再展开。

  求得LDLT的D之后，接下来求L。这里有着繁琐的数学推导，此处不再展开*2。

  二、LDLT分解的函数（.m文件后缀）

  废话少说，直接上LDLT分解的函数代码：

  %LDL分解MATLAB程序function[x,L,D]=mldl(A,b)%用途：用平方根法解线性方程组Ax=b,A=LL'%格式：[x,L]=mchol(A,b),A为系数矩阵,b为右端向量%返回：x-解向量,L-下三角阵,D-对角阵n=length(b);L=zeros(n,n);L(logical(eye(size(L)))) = repmat(1,1);  %L主对角线全是1DU = diag(basis_change(A)); %DU是把A经过初等行变换之后的上三角矩阵D = zeros(3);D(logical(eye(size(D))))=DU;  %把D对角线上的元素更改为DU的对角线上的元素for k=2:n    L(k:n,1) = A(k:n,1)/D(1,1);    L(k+1:n,k)=[A(k+1:n,k)-L(k+1:n,1:k-1)*D(k-1,k-1)*(L(k,1:k-1)')]/D(2,2);   end%解下三角方程组Ly=by=zeros(n,1);y(1)=b(1);for k=2:n    y(k)=b(k)-L(k,1:k-1)*y(1:k-1);end%解上三角方程组DL'x=yx=zeros(n,1);x(n)=y(n)/(D(n,n)*L(n,n)');DLT = D*L';for k=n-1:-1:1    x(k)=(y(k)-DLT(k,k+1:n)*x(k+1:n))/DLT(k,k);end

  最后求解x的时候，实际上还是偷懒用的LU分解的方法，因为我实在写不出x的推导式。。。

  代码中的basis_change函数是我搜到的能够使用初等变换将矩阵变为上三角矩阵的代码，文章链接放在这里，把函数名修改成basis_change即可：

  matlab中将任意矩阵转换成上三角矩阵的源码

  最后在命令行或脚本文件中调用函数求解即可：

  clc,clear% 定义系数矩阵A和右端项bA = [4 -2 -4;-2 17 10;-4 10 9]b = [10 3 -7]';[x,L,D] = mldl(A,b)

  即可返回求得的x、L和D。

  总结

  matlab官方的内置函数ldl可以进行LDLT分解，但是笔者用主流方法试图查看ldl函数的源代码，都失败了，所以就自行编了一个LDLT函数代码。以上就是这篇博客的全部内容

  其实这篇博客的内容不难，自行推导之后代码很简单。但是对于新手还是一个坎吧~

  写这篇博客的时候，还没有到提交实验报告的时间，或许有人看到后会直接copy也不做修改吧。

  上面的代码可优化的空间很大，但笔者认为目前够用了，就不再修改。

  无所谓，写blog的意义就是给人看的。

  能帮到你们，我很开心。