几种常见的矩阵分解综合程序_matlab

William国学

于 2022-01-17 16:17:10 发布

阅读量924

点赞数

文章标签：矩阵 matlab 线性代数

本文链接：https://blog.csdn.net/qq_41277628/article/details/122542294

版权

矩阵分析与应用–矩阵分解

参数输入：

待分解矩阵： $A$

求解方法： $m e t h o d$

LU分解（‘LU’）
QR分解-古典施密特正交法（‘QR-古典型施密特’）
QR分解-改进施密特正交法（‘QR-改进型施密特’）
Householder分解（‘HouseholderReduction’）
Givens分解（‘GivensReduction’）
URV分解（‘URV’）
求行列式（‘Det’）
求解向量： $b$

原理说明：

LU分解
本程序使用部分主元法进行LU分解，要求输入矩阵可逆。
分解结果 $P A = L U$ ，其中， $P$ 为初等行变换阵， $L$ 为主对角线元素为1的下三角矩阵， $U$ 为上三角矩阵。

本程序实现了四种QR分解的方法，分别是 Gram-Schmidt、modified Gram-Schmidt、Householder reduction 和 Givens reduction。要求输入矩阵列向量无关，分解所得矩阵 $A = Q R$ , 其中Q为正交方阵，R为上三角矩阵。

QR分解-古典施密特正交法
QR分解-改进施密特正交法
QR分解-Householder分解
QR分解-Givens分解
URV分解
$A=URV^T$ ,其中， $A$ 为 $m * n$ 矩阵， $U$ 和 $V$ 分别为 $m$ 阶和 $n$ 阶正交矩阵。
$U$ 的前 $r$ 列是 $R (A)$ 的标准正交基，后 $m - r$ 列是 $N(A^T)$ 的标准正交基。
$V$ 的前 $r$ 列是 $R(A^T)$ 的标准正交基，后 $m - r$ 列是 $N (A)$ 的标准正交基.
求行列式
1. $d e t (A B) = d e t (A) d e t (B)$ ；
2. 每进行一次基本行变换，行列式的值取负。
3. 三角阵的行列式的值为主对角线元素乘积。
  因此进行LU分解，PA=LU，记录P的行交换次数n，即 $1）^n*det(A)=det(L)det(U)$ .
求解向量
求解向量的有三种情况，即无解、有唯一解、有无穷多解。
本程序在QR分解的基础上：
- 当存在唯一解时，使用解方程法逐步求解上三角阵R的逆，并将其代入方程中进行求解方程组。
- 当方程不可解时返回错误信息
- 当方程有无穷多解时返回基础解系和特解，

输入输出示例：

LU分解

输入：
- A=[0 -20 -14; 3 27 -4; 4 11 -2]
- method = ‘LU’
调用：
- [P,L,U] = matrixAnalysis(A,‘LU’)

QR分解（古典施密特正交法）

输入：
- A=[0 -20 -14; 3 27 -4; 4 11 -2]
- method = ‘QR-古典型施密特’
调用：
- [Q,R] = matrixAnalysis(A,‘QR-古典型施密特’)

QR分解（改进施密特正交法）

输入：
- A=[0 -20 -14; 3 27 -4; 4 11 -2]
- method = ‘QR-古典型施密特’
调用：
- [Q,R] = matrixAnalysis(A,‘QR-改进型施密特’)
输入输出程序截图

Householder分解

输入：
- A=[0 -20 -14; 3 27 -4; 4 11 -2]
- method = ‘HouseholderReduction’
调用：
- [Q,R] = matrixAnalysis(A,‘HouseholderReduction’)

Givens分解

输入：
- A=[0 -20 -14; 3 27 -4; 4 11 -2]
- method = ‘GivensReduction’
调用：
- [Q,R] = matrixAnalysis(A,‘GivensReduction’)

URV分解

输入：
- A = [-4 -2 -4 -2;2 -2 2 1;-4 1 -4 -2];
- method = ‘URV2’
调用：
- [U,V,R] = matrixAnalysis(A,‘URV2’)

求行列式（‘Det’）

输入：
- A=[0 -20 -14; 3 27 -4; 4 11 -2]
- method = ‘det’
调用：
- DetA = matrixAnalysis(A,‘det’)

求解向量： $b$

输入：
- A = [1 2 -3 4;4 8 12 -18;2 3 2 1;-3 -1 1 -4];
- b = [4 6 8 -7]’;
- method = b
调用：
- x = matrixAnalysis(A,b)
输入：
- A = [1 2 -3 4;4 8 12 -18;2 3 2 1;1 2 -3 4];
- b = [4 6 8 -7]’;
- method = b
调用：
- x = matrixAnalysis(A,b)
输入：
- A = [1 2 -3 4;4 8 12 -18;2 3 2 1;1 2 -3 4];
- b = [4 6 8 4]’;
- method = b
调用：
- x = matrixAnalysis(A,b)

程序

% % LU
% % QR（Gram-Schmidt）
%Orthogonal Reduction
% % Householder Reduction 
% % Givens reduction)
% % URV
% %程序实现，要求如下：
% 
% 1、一个综合程序，根据选择参数的不同，实现不同的矩阵分解；在此基础上，实现Ax=b方程组的求解，以及计算A的行列式；
% 
% 2、可以用matlab、Python等编写程序，需附上简单的程序说明，比如参数代表什么意思，输入什么，输出什么等等，附上相应的例子；
% 
% 3、一定是可执行文件，例如 .m文件等,不能是word或者txt文档。附上源代码，不能为直接调用matlab等函数库;
function varargout = matrixAnalysis(A,method)
    A
    [m,n] = size(A);
    if isnumeric(method) 
        b = method
        [~,Q1,R1] =  Hous(A);
        [~,~,R2] =  Hous([A b]);
        R1 = roundn(R1,-10);
        R2 = roundn(R2,-10);
        R1(all(R1==0,2),:) = [];
        R2(all(R2==0,2),:) = [];
        
        if(size(R1,1)~=size(R2,1))
            fprintf("方程无解");
            varargout{1} = [];
            return;
            
        elseif(size(R1,1)==size(R2,1)&&size(R2,1)==n)
            R1contrary = contraryTra(R1);
            fprintf("·······方程有唯一解为：\n");
            x = R1contrary*Q1'*b;
            varargout{1} = x;
        else
%             
            R3= [A b];
            P1 = eye(m);
            for i=1:m-1
                [~,currentMaxIndex] = max(abs(R3(i:m,i)));%寻找主元
                P1([i currentMaxIndex+i-1],:)=P1([currentMaxIndex+i-1 i],:);%主元行变换
                R3([i currentMaxIndex+i-1],:)=R3([currentMaxIndex+i-1 i],:);%主元行变换

                for j=i+1:m
                    k = (R3(j,i)/R3(i,i));
                    R3(j,:) = R3(j,:)-k*R3(i,:);
                    P1(j,:) = P1(j,:)-k*P1(i,:);   
                end
            end
            R3;
            % U
            R3(all(abs(R3)<=1e-10,2),:) = [];
            [m1,~] = size(R3);

            [~,Ind]=max(R3~=0 , [] , 2);%z主元位置
            for i=1:m1
                R3(i,:) = R3(i,:)/R3(i,Ind(i));
                for j=1:i-1
                    k = (R3(j,Ind(i))/R3(i,Ind(i)));
                    R3(j,:) = R3(j,:) - R3(i,:)*k;
                    P1(j,:) = P1(j,:) - P1(i,:)*k;
                end
            end
            R3;
            %U = [A(:,Ind) P1(m1+1:m,:)'];

            % R

            funSol = zeros(n,n-m1);
            a=1:n;
            [c i]=setdiff(a,Ind);
            t=a(sort(i));

            funSol(1:m1,:) = -R3(:,t);
            fprintf("·······方程有无穷多解,基础解系为：\n");
            funSol(m1+1:n,:) = eye(n-length(Ind))
            fprintf("·······方程有无穷多解,其特解为：\n");
            praSol = [R3(:,m+1);zeros(1,n-size(R2,1))]
            x = [funSol praSol]; 
            varargout{1} = x;
        end
        
    else
        switch method
            case "LU"
                [~,P,Q,R] = LU(A);
                varargout{1} = P;
                varargout{2} = Q;
                varargout{3} = R;
            case "QR-古典型施密特"
                [P,Q,R] = smiClass(A);
                varargout{1} = Q;
                varargout{2} = R;
            case "QR-改进型施密特"
                [P,Q,R] = smiImpr(A);
                varargout{1} = Q;
                varargout{2} = R;

            case "GivensReduction"
                [P,Q,R] = Givens(A);
                varargout{1} = Q;
                varargout{2} = R;
            case "HouseholderReduction"    
                [P,Q,R] = Hous(A);
                varargout{1} = Q;
                varargout{2} = R;
            case "URV1"
                [U,V,R] = URV1(A);
                varargout{1} = U;
                varargout{2} = V;
                varargout{3} = R;
            case 'URV2'
                [U,V,R] = URV2(A);
                varargout{1} = U;
                varargout{2} = V;
                varargout{3} = R;
            case 'det'
                [P,Q,R,Det] = ourDet(A);
                varargout{1} = Det;
            
        end
    end
end
    

function [count,P,Q,R] = LU(A)
    [m,n] = size(A);
    if m~=n
        error("······您输入的不是方阵,无法进行LU分解或行列式计算·······");
    end
    count=0;
    R = A;
    P = eye(n);
    Q = eye(n);
    for i=1:n-1
        a =eye(n);%初始化P0矩阵
        b = a;%初始化M矩阵
        [~,currentMaxIndex] = max(abs(R(i:n,i)));%寻找主元
        if((currentMaxIndex+i-1)~=i)
            count=count+1;
        end
        a([i currentMaxIndex+i-1],i:n)=a([currentMaxIndex+i-1 i],i:n);%主元行变换
        %逐步计算 P M
        P0{i} = a;
        R = P0{i}*R;
        M{i} = b;
        for j=i+1:n
            M{i}(j,i) = -(R(j,i)/R(i,i));
        end
        R = M{i}*R;%迭代变换
        P = P0{i}*P;
        Q = P0{i}*Q/(P0{i})*inv(M{i});
    end
%     P = roundn(P,-3);
%     Q = roundn(Q,-3);
%     R = roundn(R,-3);
end

function [P,Q,R] = smiClass(A)
    [m,n] = size(A);
    if m<n
        error("········您输入的矩阵列相关········");
    end
    P = zeros(m);
    for i =1:m
        u = A(:,i);
        for j = 1:i-1
            u = u - P(:,j)'*A(:,i)*P(:,j);
        end
        P(:,i) = u/(sqrt(sum(u.^2)));
    end
    P;
    Q = P';
    R = Q*A;
%     P = roundn(P,-3);
%     Q = roundn(Q,-3);
%     R = roundn(R,-3);
end
function [P,Q,R] = smiImpr(A)
    [m,n] = size(A);
    if m<n
        error("········您输入的矩阵列相关········");
    end
    [m,n] = size(A);
    P = zeros(m);
    U = A;
    for i =1:n   
        P(:,i) = U(:,i)/(sqrt(sum(U(:,i).^2)));
        for j = i+1:n
            U(:,j) = U(:,j) - P(:,i)'*A(:,j)*P(:,i);
        end
    end
    Q = P';
    R = Q*A;
%     P = roundn(P,-3);
%     Q = roundn(Q,-3);
%     R = roundn(R,-3);
end
function [P,Q,R] = Givens(A)
    [m,n] = size(A);
    R=A;
    P = eye(m);
    for i=1:n-1 
        T_i = eye(m);
        for j = i+1:m 
            T_ii = eye(m);
            u = R(:,i);
            var1 = u(i);
            var2 = u(j);
            cos_a = var1/sqrt(var1^2+var2^2);
            sin_a = var2/sqrt(var1^2+var2^2);
            T_ii(i,i) = cos_a;
            T_ii(i,j) = sin_a;
            T_ii(j,j) = cos_a;
            T_ii(j,i) = -sin_a;
            R = T_ii*R;
            T_i = T_ii*T_i;
        end
        P = T_i*P;
    end
    Q = P';
    R;
%     P = roundn(P,-3);
%     Q = roundn(Q,-3);
%     R = roundn(R,-3);
end
function [P,Q,R] = Hous(A)
    [m,n] = size(A);
    t = min(n,m-1);%约简次数
    E = eye(m);
    R = A;
    P = E;
    for i=1:t
        u = R(i:m,i);
        u(1) = u(1)-sqrt(sum(abs(R(i:m,i)).^2));
        if(u(1:m-i+1)==0)
            P_small = eye(m-i+1);
        else
             P_small = eye(m-i+1)-((2/(u(1:m-i+1)'*u(1:m-i+1)))*(u(1:m-i+1)*u(1:m-i+1)'));
        end
        R(i:m,i:n) = P_small*R(i:m,i:n);
        P_i = E;
        P_i(i:m,i:m) = P_small;
        P = P_i*P;
    end
    Q = P'; 
%     P = roundn(P,-3);
%     Q = roundn(Q,-3);
%     R = roundn(R,-3);
end
function [U1,V1,R1] = URV2(A)
    [m,n] = size(A);
    R= A;
    P = eye(m);
    %部分主元法化简
    for i=1:min(m-1,n-1)
        [~,currentMaxIndex] = max(abs(R(i:m,i)));%寻找主元
        P([i currentMaxIndex+i-1],:)=P([currentMaxIndex+i-1 i],:);%主元行变换
        R([i currentMaxIndex+i-1],:)=R([currentMaxIndex+i-1 i],:);%主元行变换
        for j=i+1:m
            k = (R(j,i)/R(i,i));
            R(j,:) = R(j,:)-k*R(i,:);
            P(j,:) = P(j,:)-k*P(i,:);      
        end
    end
    P;
    R;
    %% U
    R(all(abs(R)<=1e-10,2),:) = [];%消除精度影响
    [m1,~] = size(R);%行秩   
    [~,Ind]=max(R~=0 , [] , 2);%z主元位置
    %极简化
    for i=1:m1
        R(i,:) = R(i,:)/R(i,Ind(i));
        for j=1:i-1
            k = R(j,Ind(i))/R(i,Ind(i));
            R(j,:) = R(j,:) - R(i,:)*k;
            P(j,:) = P(j,:) - P(i,:)*k;
        end
    end
    R;
    P;
    %主值空间计算
    U = [A(:,Ind) P(m1+1:m,:)'];%R(A)  N(A^T)
    N_a = zeros(n,n-m1);% N(A)
    a=1:n;
    [c i]=setdiff(a,Ind);
    x=a(sort(i));

    N_a(1:m1,:) = -R(:,x);
    N_a(m1+1:n,:) = eye(n-length(Ind));
    R';%R(A^T)
    V = [R' N_a];

    U1 = smiImpr(U);
    V1 = smiImpr(V);

    R1 = U1'*A*V1;
%     U1 = roundn(U1,-3);
%     V1 = roundn(V1,-3);
%     R1 = roundn(R1,-3);
end

function [U,V,R] = URV1(A)
    [m,n] = size(A);
    [P,~,B] = Hous(A);
    B = roundn(B,-3);
    B(all(B==0,2),:) = [];
    B;
    [Q,~,~] = Hous(B');
    U = P';
    V = Q;
    R = U'*A*V;
    R = roundn(R,-3);
    U*R*V;
%     U = roundn(U,-3);
%     V = roundn(V,-3);
%     R = roundn(R,-3);
end

function [P,L,U,Det] = ourDet(A)
    [m,~] = size(A);
    [count,P,L,U] = LU(A);
    Det = 1;
    for i=1:m
        Det = Det*U(i,i);
    end
    Det = ((-1)^count)*Det;
%     P = roundn(P,-3);
%     L = roundn(L,-3);
%     U = roundn(U,-3);
end

%计算逆矩阵(解方程法）
function returnContrary = contraryTra(matrixA)
    [n ,m] = size(matrixA);
    if n~=m
     fprintf("您输入的不是方阵\n");
     return;
    end
    returnContrary = zeros(n);
    for i = 1:n
        returnContrary(i,i) = 1/ matrixA(i,i);
    end
    for t=1:n-1
            for k = t+1:n
                j = n-k+t+1;
                i = j - t;
                for m = (i+1):j
                    returnContrary(i,j) = (returnContrary(i,j)+returnContrary(m,j)*matrixA(i,m));            
                end
                returnContrary(i,j) = returnContrary(i,j)*(-1/matrixA(i,i));
            end
    end
end