《矩阵论》第四章——矩阵分解代码

在矩阵论的学习过程中，编写了课本中的矩阵分解的Matlab代码，以便于加深印象。

高斯消去法

clc;clear;close all;
A = [10,-7,0;-3,2,6;5,-1,5];%系数矩阵
b = [7,4,6];%常数向量
n = length(b);
x = zeros(n,1);
%消去过程
for i = 1:n-1 %该行的主对角元素之前的元素置零，可以看作是基准行，在每一次消去过程中该行不变
    for j = i+1:n  %进行操作的行
        c = (-A(j,i))/A(i,i);   %倍乘因子
        for k = 1:n %进行行操作
            A(j,k) = A(j,k)+c*A(i,k);   %对矩阵进行操作
        end
        b(j) = b(j)+c*b(i);
    end
end
%回代过程
for i = n:-1:1
    x(i) = b(i);
    for j = n:-1:i+1
        x(i) = x(i)-A(i,j)*x(j);
    end
    x(i) = x(i)/A(i,i);
end

高斯列主元消去法

clc;clear;close all;
A = [2,1,2;5,-1,1;1,-3,-4];%系数矩阵
b = [5,8,-4];%常数向量
n = length(b);
x = zeros(n,1);%解向量
tempM = zeros(1,n);%交换的中间向量
tempV = 0;
%消元过程
for i = 1:n-1
    [r,c] = find(abs(A(:,i)) == max(abs(A(:,i))));
    temp = A(i,:);
    A(i,:) = A(r,:);
    A(r,:) = temp;
    tempV = b(i);
    b(i) = b(r);
    b(r) = tempV;
    for j = i+1:n  %进行操作的行
        c = (-A(j,i))/A(i,i);   %倍乘因子
        for k = 1:n %进行行操作
            A(j,k) = A(j,k)+c*A(i,k);   %对矩阵进行操作
        end
        b(j) = b(j)+c*b(i);
    end
end
%回代过程
for i = n:-1:1
    x(i) = b(i);
    for j = n:-1:i+1
        x(i) = x(i)-A(i,j)*x(j);
    end
    x(i) = x(i)/A(i,i);
end

三角分解

clc;clear;close all;
A = [6,2,1,-1;2,4,1,0;1,1,4,-1;-1,0,-1,3];%系数矩阵
b = [6,-1,5,-1];%常数向量
n = length(b);
x = zeros(n,1);
L = zeros(n,n);
x = zeros(n,1);
y = zeros(n,1);
%LU分解
for i = 1:n
    L_1 = eye(n);
    for j = i:n
        L(j,i) = A(j,i)/A(i,i);
        if j > i
            L_1(j,i) = -L(j,i);
        end
    end
    A = L_1*A;
end
U = A;
%回代Ly=b
for i = 1:n
    y(i) = b(i);
    if i == 1
        y(i) = y(i)/L(i,i);
        continue
    end
    for j = 1:i-1
        y(i) = y(i)-L(i,j)*y(j);
    end
    y(i) = y(i)/L(i,i);
end
%回代Ux=y
for i = n:-1:1
    x(i) = y(i);
    for j = n:-1:i+1
        x(i) = x(i)-U(i,j)*x(j);
    end
    x(i) = x(i)/U(i,i);
end

Crout分解（Doolittle分解与Crout分解完全一致，就是LU分解顺序变成UL）

clc;clear;close all;
A = [6,2,1,-1;2,4,1,0;1,1,4,-1;-1,0,-1,3];%系数矩阵
b = [6,-1,5,-1];%常数向量
n = length(b);
L = zeros(n,n);
U = eye(n,n);
L(:,1) = A(:,1);
U(1,2:n) = A(1,2:n)/L(1,1);
%Crout分解
for k = 1:n
    for i = k:n
        L(i,k) = A(i,k);
        for m = 1:k-1
            L(i,k) = L(i,k)-L(i,m)*U(m,k);
        end
    end
    for j = k+1:n
        U(k,j) = A(k,j);
        for m = 1:k-1
            U(k,j) = U(k,j)-L(k,m)*U(m,j);
        end
        U(k,j) = U(k,j)/L(k,k);
    end
end
%回代Ly=b
for i = 1:n
    y(i) = b(i);
    if i == 1
        y(i) = y(i)/L(i,i);
        continue
    end
    for j = 1:i-1
        y(i) = y(i)-L(i,j)*y(j);
    end
    y(i) = y(i)/L(i,i);
end
%回代Ux=y
for i = n:-1:1
    x(i) = y(i);
    for j = n:-1:i+1
        x(i) = x(i)-U(i,j)*x(j);
    end
    x(i) = x(i)/U(i,i);
end

CHolesky分解

clc;clear;close all;
A = [6,2,1,-1;2,4,1,0;1,1,4,-1;-1,0,-1,3];%系数矩阵
n = length(b);
G = zeros(n,n);
%Cholesky分解
G(1,1) = A(1,1)^0.5;
G(2:n,1) = A(2:n,1)/G(1,1);
for i = 2:n
    for j = 2:n
        if i == j
            G(i,j) = A(i,i);
            for k = 1:i-1
                G(i,j) = G(i,j)-G(i,k)^2;
            end
            G(i,j) = G(i,j)^0.5;
        elseif i > j
            G(i,j) = A(i,j);
            for k = 1:j-1
                G(i,j) = G(i,j)-G(i,k)*G(j,k);
            end
            G(i,j) = G(i,j)/G(j,j);
        else 
            G(i,j) = 0;
        end
    end
end

矩阵的QR分解

clc;clear;close all;
A = [1,2,2;2,1,2;1,2,1];
n = length(A(:,1));
B = zeros(n,n);
C = eye(n);
Q = zeros(n,n);
R = zeros(n,n);
B(:,1) = A(:,1);
%Schmidt正交化
for i = 2:n
    B(:,i) = A(:,i);
    for j = 1:i-1
        B(:,i) = B(:,i)-(B(:,j)'*A(:,i))/(B(:,j)'*B(:,j))*B(:,j);
        C(j,i) = (A(:,i)'*B(:,j))/(B(:,j)'*B(:,j));
    end
end
%计算QR矩阵
for i = 1:n
    q = norm(B(:,i));
    Q(:,i) = B(:,i)/q;
    R(i,i) = q;
end
R = R*C;