矩阵的各类分解

最新推荐文章于 2024-03-12 14:43:51 发布

papaofdoudou

最新推荐文章于 2024-03-12 14:43:51 发布

阅读量1k

点赞数 7

分类专栏：工程数学

本文链接：https://blog.csdn.net/tugouxp/article/details/112553279

版权

数学同时被 2 个专栏收录

147 篇文章 19 订阅

订阅专栏

工程

52 篇文章 1 订阅

订阅专栏

1.LU分解

在线性代数中已经证明，如果方阵是非奇异的，即的行列式不为0，LU分解总是存在的，它的形式为：

$\begin{bmatrix} a_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ \vdots & \vdots & \ddots & \vdots \\ a_{n1}& a_{n2}& \cdots & a_{nn} \end{bmatrix}=\begin{bmatrix} 1& 0& \cdots & 0\\ l_{21}& 1& \cdots & 0\\ \vdots & \vdots & \ddots & \vdots \\ l_{n1}& l_{n2}& \cdots & 1 \end{bmatrix}\begin{bmatrix} u_{11}& u_{12}& \cdots & u_{1n}\\ 0& u_{22}& \cdots & u_{2n}\\ \vdots & \vdots & \ddots & \vdots \\ 0& 0& \cdots & u_{nn} \end{bmatrix}$

例如，对于四阶矩阵来说，可以分解为:

$\begin{bmatrix} a_{11} & a_{12} & a_{13} & a_{14} \\ a_{21} & a_{22} & a_{23} &a_{24} \\ a_{31} & a_{32} & a_{33} & a_{34}\\ a_{41} &a_{42} &a_{43} &a_{44} \end{bmatrix} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ l_{21} & 1 & 0 &0 \\ l_{31} & l_{32} & 1} & 0\\ l_{41} &l_{42} &l_{43} &1 \end{bmatrix}\begin{bmatrix} u_{11} & u_{12} & u_{13} & u_{14} \\ 0 & u_{22} & u_{23} &u_{24} \\0 & 0 & u_{33} & u_{34}\\ 0 &0 &0 &u_{44} \end{bmatrix}=\begin{bmatrix} u_{11} & u_{12} & u_{13} & u_{14}\\ l_{21}u_{11}& l_{21}u_{12}+u_{22}&l_{21}u_{13}+u_{23} & l_{21}u_{14} +u_{24} \\ l_{31}u_{11}&l_{31}u_{12} + l_{32}u_{22} &l_{31}u_{13}+l_{32}u_{23}+u_{33} &l_{31}u_{14}+l_{32}u_{24}+u_{34}\\ l_{41}u_{11}&l_{41}u_{12}+l_{42}u_{22} &l_{41}u_{13}+l_{42}u_{23}+l_{43}u_{33} &l_{41}u_{14}+l_{42}u_{24}+l_{43}u_{34} + u_{44} \end{bmatrix}$

寻找规律：

第一行没有依赖因素，因此可以直接得出，第一列仅仅依赖第一行，所以在得到第一行之后，第一列也可以计算求得.
同理，第二主行也是这样，可以先计算主行，在计算主列，也就是按照下图中的，先实线，再虚线的顺序计算.
仔细观察，实际上只要先计算出U矩阵主对角线元素 $u_{ii}$ ，至于先计算行还是列，无关紧要，因为除了对应行，列的未知lu元素，不在依赖其他参变量.
每个原矩阵中的元素 $a_{ij}$ 仅仅使用一次，计算出对应的 $u_{ij}$ 后，其信息已经扩散到整个计算结构中，可以舍弃。后续计算不再依赖，对计算机程序计算来说，好处是为算法的内存优化，提供的抓手。

按照上面的发现的规律，尝试推导计算公式如下:

按照U的行线性组合

$Row_1(A)=1*Row_{1}(U) +0*Row_2(U)+\cdots+0*Row_n(U) = Row_1(U)$

得到:

$\\Row_1(A) = Row_1(U)\Rightarrow u_{1j} = a_{1j}, \qquad j=1,2,\cdots,n$

另外,按照L的列的线性组合

$\\Colum_1(A)=u_{11}*Colum_{1}(U) +0*Colum_2(U)+\cdots+0*Colum_n(U) = u_{11}Colum_1(U)\Rightarrow l_{i1}=\frac{a_{i1}}{u_{11}},\qquad i=2,3,\cdots,n$

所以得到:

$\left\{\begin{matrix} u_{1j} = a_{1j}, \qquad j=1,2,\cdots,n\\ \\ l_{i1}=\frac{a_{i1}}{u_{11}},\qquad i=2,3,\cdots,n\end{matrix}\right.$

对于内层矩阵，由于 i>j 时， $u_{ij}==0$ ,所以，为求 $u_{ij}, i>1,j\geqslant i$ 时候的值，可以计算得到:

$\\a_{ij} = l_{i1}u_{1j}+l_{i2}u_{2j}+\cdots+l_{i,i-1}u_{i-1, j}+1\cdot u_{ij} \quad u_{1j},u_{2j},\cdots,u_{i-1, j}, l_{i1}, l_{i2},\cdots, l_{i, i-1} \ are \ known.$

所以:

$\\u_{ij} = a_{ij}-\sum_{s=1}^{i-1}l_{is}u_{sj}, \qquad i > 1, j \geqslant i$

由于 i< j 时， $l_{ij}==0$ ,所以，为求 $l_{ij}, j>1,i> j$ 时候的值，可以计算得到:

$\\a_{ij} = l_{i1}u_{1j}+l_{i2}u_{2j}+\cdots+l_{i,j-1}u_{j-1, j}+l_{ij}u_{jj} \quad u_{1j},u_{2j},\cdots,u_{j, j},l_{i1},l_{i2},\cdots,l_{i,j-1} \ are \ known.$

$\\l_{ij} = (a_{ij}-\sum_{s=1}^{j-1}l_{is}u_{sj})/u_{jj}, \qquad j > 1, i > j$

所以，得到:

$\left\{\begin{matrix} u_{1j} = a_{1j}, \qquad j=1,2,\cdots,n\\ \\ l_{i1}=\frac{a_{i1}}{u_{11}},\qquad i=2,3,\cdots,n\end{matrix}\right.$

$\\ \left\{\begin{matrix} u_{ij} = a_{ij}-\sum_{s=1}^{i-1}l_{is}u_{sj}, \qquad \qquad i > 1, j \geqslant i \\ \\ \\ l_{ij} = (a_{ij}-\sum_{s=1}^{j-1}l_{is}u_{sj})/u_{jj}, \qquad j > 1, i > j \end{matrix}\right.$

根据上述两个公式，在octave上编码实现为：

% find the LU factorization of the matrix
% input:
%   a: the matrix need to be factorize
%   n: the number of the low or column in the matrix
% output:
%   no output
function LU(a,n)
    m=eye(n);    %genernate the unit matrix.
    for j = 1 : n-1
        if abs(a(j,j))<eps; 
            error('zero pivot encountered');    % when the zero pivot happens,end the process
        end
        for i = j+1 : n
            mult = a(i,j)/a(j,j);
            m(i,j) = mult;  
            for k = j:n
                a(i,k) = a(i,k) - mult*a(j,k);
            end
        end
    end
disp('  L=');  disp(m);
disp('  U=');  disp(a);
disp('  LU='); disp(m*a);          % to check if the result is right

运行结果为: