矩阵的各类分解

原创

豆豆爹 2022-12-04 00:12:33 博主文章分类：工程 ©著作权

©著作权归作者所有：来自51CTO博客作者豆豆爹的原创作品，请联系作者获取转载授权，否则将追究法律责任

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矩阵主对角线元素，至于先计算行还是列，无关紧要，因为除了对应行，列的未知lu元素，不在依赖其他参变量.
每个原矩阵中的元素仅仅使用一次，计算出对应的后，其信息已经扩散到整个计算结构中，可以舍弃。后续计算不再依赖，对计算机程序计算来说，好处是为算法的内存优化，提供的抓手。

矩阵的各类分解_计算机程序_03

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

按照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

时，

矩阵的各类分解_线性代数_09

,所以，为求

矩阵的各类分解_线性代数_10

时候的值，可以计算得到:

$\\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$

由于

矩阵的各类分解_内存优化_13

时，

$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