matlab实现矩阵节约存储的LR(Doolittle )分解以及LDR分解

• LR 分解 (Doolittle 分解)
  $$A=LR$$ 其中，L 是单位下三角矩阵, R 是上三角矩阵
  设 A 有 LR 分解：
  比较矩阵两边元素，由第一行元素相等和第一列元素分别相等，可以得到：
  
  按照第 1 行 + 第 1 列 → 第 2 行 + 第 2 列 →…的顺序类推可以得到：
  对于 $i=1,...,n$，有
  
• LR 分解matlab程序

function[L_matrix,R_matrix] = LR_1(A_matrix)
[row_A,col_A] = size(A_matrix);

% 判断是否是非奇异方阵
if row_A ~= col_A
    error('矩阵非方阵，无法进行LR分解');
end
if row_A ~= rank(A_matrix)
    error('矩阵奇异，无法进行LR分解');
end

% 初始化
L_matrix = zeros(size(A_matrix));
R_matrix = zeros(size(A_matrix));

for j = 1:col_A
    R_matrix(1,j) = A_matrix(1,j);
end
for i = 1:col_A
    L_matrix(i,i) = 1;
end
for i = 2:row_A
    L_matrix(i,1) = A_matrix(i,1) / A_matrix(1,1);
end

for i = 1:row_A
    for j = i:row_A
        sum = 0;
        for k = 1:i-1
            sum = sum + L_matrix(i,k) * R_matrix(k,j);
        end
        R_matrix(i,j) = A_matrix(i,j) - sum;
        sum_1=0;
        for k = 1:i-1
            sum_1 = sum_1 + L_matrix(j,k) * R_matrix(k,i);
        end
        L_matrix(j,i) = (A_matrix(j,i) - sum_1) / R_matrix(i,i);
    end
end

• 节约储存的LR分解
  将 A 进行 LR 分解后，不再储存原来的 A ，直接将分解的 L 和 R 按照
  
  这样的方式来存放，存放于矩阵 A 中，得到算法：
  对于 $i=1,...,n$，有
  

function[A_matrix] = LR_2(A_matrix)
% 节省存储空间的 LR 分解算法
row_A = size(A_matrix,1);
for i = 1:row_A
    T = det(A_matrix(1:i,1:i));
    if T == 0
        error('矩阵的顺序主子式为零，无法进行LR分解');
    else
        for j = i:row_A
            sum = 0;
            for k = 1:i-1
                sum = sum+A_matrix(i,k) * A_matrix(k,j);
            end
            A_matrix(i,j) = A_matrix(i,j) - sum;
        end
        for j = i + 1:row_A
            sum_1 = 0;
            for k = 1:i-1
                sum_1 = sum_1 + A_matrix(j,k) * A_matrix(k,i);
            end
            A_matrix(j,i) = (A_matrix(j,i) - sum_1) / A_matrix(i,i);
        end
    end
end

• LDR分解
  $$A=LDR$$ 其中，L 是单位下三角矩阵, D 是对角矩阵, R 是单位上三角矩阵
  可以直接在 LR 分解的基础上进行， $$A=LR=LD{D}^{-1}R=LD\tilde{R}$$其中， D 即为 LR 分解中的 R 的对角线元素组成的矩阵， $$D=diag(r_{11},r_{22},...r_{nn})$$此时 $\tilde{R}$ 为单位上三角矩阵
• LDR分解matlab代码

function[L_matrix,D_matrix,R_matrix] = LDR_1(A_matrix)
[row_A,col_A] = size(A_matrix);

% 判断是否是非奇异方阵
if row_A ~= col_A
    error('矩阵非方阵，无法进行LR分解');
end
if row_A ~= rank(A_matrix)
    error('矩阵奇异，无法进行LR分解');
end

% 初始化
L_matrix = zeros(size(A_matrix));
D_matrix = zeros(size(A_matrix));
R_matrix = zeros(size(A_matrix));

% 先进行LR分解
for j = 1:col_A
    R_matrix(1,j) = A_matrix(1,j);
end
for i = 1:col_A
    L_matrix(i,i) = 1;
end
for i = 2:row_A
    L_matrix(i,1) = A_matrix(i,1) / A_matrix(1,1);
end

for i = 1:row_A
    for j = i:row_A
        sum = 0;
        for k = 1:i-1
            sum = sum + L_matrix(i,k) * R_matrix(k,j);
        end
        R_matrix(i,j) = A_matrix(i,j) - sum;
        sum_1=0;
        for k = 1:i-1
            sum_1 = sum_1 + L_matrix(j,k) * R_matrix(k,i);
        end
        L_matrix(j,i) = (A_matrix(j,i) - sum_1) / R_matrix(i,i);
    end
end
% 进行完LR分解之后，将R分解为D*R
for i = 1:row_A
    D_matrix(i,i) = R_matrix(i,i);
end
for i = 1:row_A
    for j = 1:row_A
        if i <= j
            R_matrix(i,j) = R_matrix(i,j) / D_matrix(i,i);
        end
    end
end

• 节约储存的LDR分解matlab代码
  同样可以得到节约存储的LDR分解，将 L 放入原来 A 矩阵的下三角， D 放入对角线，R 放入上三角

function[A_matrix]=LDR_2(A_matrix)
% 节省存储空间的 LDR 分解算法
row_A = size(A_matrix,1);

% 先进行节省存储空间的LR分解
for i = 1:row_A
    T = det(A_matrix(1:i,1:i));
    if T == 0
        error('矩阵的顺序主子式为零，无法进行LR分解');
    else
        for j = i:row_A
            sum = 0;
            for k = 1:i-1
                sum = sum+A_matrix(i,k) * A_matrix(k,j);
            end
            A_matrix(i,j) = A_matrix(i,j) - sum;
        end
        for j = i + 1:row_A
            sum_1 = 0;
            for k = 1:i-1
                sum_1 = sum_1 + A_matrix(j,k) * A_matrix(k,i);
            end
            A_matrix(j,i) = (A_matrix(j,i)-sum_1) / A_matrix(i,i);
        end
    end
end

% 进行完LR分解之后，将R分解为D*R，用节省存储空间的方式储存
for i = 1:row_A
    for j = 1:row_A
        if i < j
            A_matrix(i,j) = A_matrix(i,j) / A_matrix(i,i);
            j=j+1;
        else
            A_matrix(i,j) = A_matrix(i,j);
        end
    end
end

矩阵的分解还有很多很多，这里空白太小,我写不下了 0.0