数值分析中的几个三角分解

这个系列是学习《数值分析》时，将所学内容进行简单记录，以及自己用matlab的一些实现过程。

这次是这周学的三角分解，matlab有对应的函数，在这里主要记录一下大概的分解结果和方法，最后自己用代码实现了一下。具体的求解方法和直接三角分解的方法一致，都是写出LU的形式，然后展开LU并与A对应，求取L和U的元素。

直接三角分解法

如果A是n阶方阵，并且A的前n阶顺序主子式均不为0，则A存在唯一的分解A=LU
$$A=\left[\begin{array}{ccccc}1& & & & \\ l_{21}& 1& & & \\ l_{31}& l_{32}& 1& & \\ ⋮& ⋮& ⋮& \ddots & \\ l_{n1}& l_{n2}& l_{n3}& \cdots & 1\end{array}\right]\left[\begin{array}{ccccc}{u}_{11}& u_{12}& u_{13}& \cdots & u_{1n}\\ & u_{22}& u_{23}& \cdots & u_{2n}\\ & & u_{33}& \cdots & u_{3n}\\ & & & \ddots & ⋮\\ & & & & u_{nn}\end{array}\right]$$
前边的下三角矩阵是L，后边的上三角矩阵是U。

求解时，先求U的第一行，用L的第一列与U做乘法，得到$A(1,:)=L(1,:)U=U(1,:)$，即U的第一行与A的第一行相等。
然后求L的第一列，用L与U的第一列做乘法，即$A(:,1)=L\ast U(:,1)=L(:,1)\ast u_{11}$，所以，U的第一列就是$A(:,1)/u_{11}$，第一个元素是1，也可以不用计算，$A(2:end,1)\ast U(:,1)$。

对于L和U中间行列的计算，也是同理，假设前L的前r-1列和U的前r行已经得到，现在求取L的第r列和U的第r行。
求U的第r行第j列的元素（j>=r）
$$\begin{gathered} a_{r,j}=L(r,:)\ast U(:,j) \\ =\left[\begin{array}{ccccccc}{l}_{r,1}& \cdots & l_{r,r-1}& 1& 0& \cdots & 0\end{array}\right]\ast {\left[\begin{array}{ccccccc}{u}_{1,j}& \cdots & u_{r-1,j}& u_{r,j}& 0& \cdots & 0\end{array}\right]}^T \\ =\sum _{m=1}^{r-1}(l_{r,m}{u}_{m,j})+u(r,j) \end{gathered}$$
所以
$$u_{r,j}=a_{r,j}-\sum _{m=1}^{r-1}(l_{r,m}{u}_{m,j})$$

对于L的第i行也是类似，求L的第r列的第i个元素，r<i<=n。

$$\begin{gathered} a_{i,r}=L(i,:)\ast U(:,r) \\ =\left[\begin{array}{cccccccc}{l}_{i,1}& \cdots & l_{i,r}& \cdots & l_{i,i}& 0& \cdots & 0\end{array}\right]\ast \\ {\left[\begin{array}{cccccc}{u}_{1,r}& \cdots & u_{r,r}& 0& \cdots & 0\end{array}\right]}^T \\ =\sum _{m=1}^{r-1}(l_{i,m}\ast u_{m,r})+l_{i,r}\ast u_{r,r} \end{gathered}$$
所以$$l_{i,r}=\frac{a_{i,r}-{\sum }_{m=1}^{r-1}(l_{i,m}\ast u_{m,r})}{u_{r,r}}$$
下边是一张图比较清晰地说明这个步骤：

比如求U的第r行第j列时，用这一行左边黄色部分（属于L的前r-1列）和第r行第j列上方的黄色部分（属于U的前r-1行）作乘积，然后用A(r,j)减去这个乘积就是$u_{r,j}$。
求L的第i行第r列时，用这一列上方的蓝色部分（属于U的前r-1行）与对应列的蓝色部分（属于L的前r-1列）做乘积，假设值是sum，然后$$l_{i,r}=\frac{a_{i,r}-sum}{u_{r,r}}$$

自己用matlab实现的代码：

%对应于高斯消元法的LU分解，结果都存在A中
function A = LU(A)
%U的第一行就是A的第一行
%一个循环处理L的第i列（i+1:n,i），U的第i+1行（i+1,i+1:n）
[rows,cols] = size(A);
if rows~=cols
    error("rows != cols");
end

for i=1:rows-1
    if i==1
        A(i+1:end,i)=A(i+1:end,i)/A(i,i);%第一列
        A(i+1,i+1:end) = A(i+1,i+1:end)- A(i+1,1:i)*A(1:i,i+1:end);
        continue
    end
    
    %列
    A(i+1:end,i) = (A(i+1:end,i) - A(i+1:end,1:i-1)*A(1:i-1,i))/A(i,i); 
    A(i+1,i+1:end) = A(i+1,i+1:end)- A(i+1,1:i)*A(1:i,i+1:end);
end
end

这里直接将LR合并写到了A中，方便和手动计算题目作比较。

列主元的三角分解

列主元的三角分解和之前列主元法一样，只不过是把消元过程的系数存储到A的下半部分，并且在做行变换时，对当前A的所有列都做行变换（即L也参与行变换）。列主元法的三角分解可以写为：
$$PA=LU$$
从这里也可以看出为什么要对L进行行变换，普通的三角分解假设是$A=lU$，则列主元法的三角分解为$PA=PlU=LU$。
下边是实现的代码，在上一篇中列主元法的基础上稍作修改。

%对应于列主元高斯消元法的LU分解，结果都存在A中
%和普通的相比，就是多了一部初等行交换。
%PA=LU，P是多次初等行变换叠加的矩阵
%相对于LU分解，这里的L也要行交换，多了次找列主元
%L在A的下三角部分，U在A的对角线及以上
function [P,A] = lzy_LU(A)
%U的第一行就是A的第一行
%一个循环处理L的第i列（i+1:n,i），U的第i+1行（i+1,i+1:n）
[rows,cols] = size(A);
if rows~=cols
    error("rows != cols");
end
P=eye(rows);

for i=1:rows-1
    %第i行，第i列
    [M,idx]=max(abs(A(i:end,i)));
    idx=idx+i-1;
    %注意这里是全部交换
    A([i,idx],:) = A([idx,i],:);
    P([i,idx],:)=P([idx,i],:);
    
    A(i+1:end,i) = A(i+1:end,i)/A(i,i);%L的第i列
    A(i+1:end,i+1:end) = A(i+1:end,i+1:end)-A(i+1:end,i)*A(i,i+1:end); %更新U及A

   
end

end

三对角矩阵的三角分解（追赶法）

如果矩阵A的形式为
$$\left[\begin{array}{ccccc}{b}_1& c_1& & & \\ a_2& b_2& c_2& & \\ & \ddots & \ddots & \ddots & \\ & & a_{n-1}& b_{n-1}& c_{n-1}\\ & & & a_n& b_n\end{array}\right]$$
并且满足$|b_1|>|c_1|>0,|b_i|\ge |a_i|+|c_i|,|b_n|>|a_n|>0$，则A可分解为：
$$LU=\left[\begin{array}{ccccc}1& & & & \\ l_2& 1& & & \\ & \ddots & \ddots & & \\ & & l_{n-1}& 1& \\ & & & l_n& 1\end{array}\right]\left[\begin{array}{ccccc}{u}_1& c_1& & & \\ & u_2& c_2& & \\ & & \ddots & \ddots & \\ & & & u_{n-1}& c_{n-1}\\ & & & & 1u_n\end{array}\right]$$
根据这个形式可以仿照上边计算LU矩阵，并能简化。
$u_1=b_1,l_2=a_2/u_1,u_2=b_2-c_1\ast l_2$
$l_k=a_k/u_{k-1},u_k=b_k-l_k\ast c_{k-1}$
这里计算$l_k$时，$l_k$左边全为0，所以$a_k$不用减左边向量和当前列上边向量的乘积。计算$u_k$时，$u_k$左边为$l_k$，上边为$c_{k-1}$，内积也只有这两个值相乘。

对称正定矩阵的三角分解

后边几个晚上再写。