原创，注释全，可以参考也可以拿来直接用。内含五大分解，计算矩阵行列式等matlab代码
本人使用的是MATLAB R2022b，如果你的MATLAB版本较低可能会报错。
老师给出的要求：
一个综合程序，根据选择参数的不同，实现不同的矩阵分解；
在此基础上，实现Ax=b方程组的求解，以及计算A的行列式；
Faciorization.m根据选择参数的不同，实现不同的矩阵分解。LU.m, GramSchmidt.m, Householder.m, Givens.m, URV.m分别为LU分解，Gram-Schmidt，Householder Reduction，Givens Reduction，URV相应的程序。Axb.m为求解方程Ax=b的程序。DET.m为求解方程行列式的程序。main.m为最终执行文件，具体用法均有注释。
LU.m
function[l,u,p] = LU(a)
% 示例：输入矩阵a=
% 1 2 -3 4
% 4 8 12 -8
% 2 3 2 1
% 3 -1 1 -4
[m,n]=size(a); % 获取a的行列数
if m~=n
error('矩阵非方阵！')
end
p = eye(n); % 初始化P L U矩阵
l = zeros(n,n);
u = a;
for j = 1 : n-1
% 示例：第二轮循环开始时，j=2
% u = 4 8 12 −8
% 0 0 −6 6
% 0 −1 −4 5
% 0 5 10 −10
[~, index] = max(abs(u(j:n,j))); % 寻找当前枢轴下最大主元的相对位置索引
% u(2:4,2)为[0 -1 5]'，max得到index=3
index = index + j - 1; % index表示最大主元在矩阵中的行数索引
% 第二轮循环中，5所在行数为4
p([j,index],:) = p([index,j],:); % P矩阵进行相同的行变换
% p([2 4],:) = p([4 2],:) 交换p矩阵R2与R4
u([j,index],:) = u([index,j],:); % 行变换使当前列的最大主元到枢轴位置
% u([2 4],:) = u([4 2],:) 交换u矩阵R2与R4
l([j,index],:) = l([index,j],:);
% u = 4 8 12 −8
% 0 5 10 −10
% 0 −1 −4 5
% 0 0 −6 6
for i = j+1 : n % j=2,i=3时
multi = u(i,j)/u(j,j); % multi = u(3,2)/u(2,2) = -1/5
l(i,j) = multi; % l矩阵的(3,2)位置写入-1/5
for k = j:n % 改变第i行元素，j列之前元素已经为0
u(i,k) = u(i,k) - multi * u(j,k);
end
% u = 4 8 12 −8
% 0 5 10 −10
% 0 0 −6 6
% 0 0 −2 3
end
end
% 将l对角线置为1
for k = 1:n
l(k,k)=1;
end
endGramSchmidt.m
function[q,r] = GramSchmidt(a)
% a= 0 −20 −14
% 3 27 −4
% 4 11 −2
[m,n]=size(a);
if rank(a) ~= n
error('!!!!!!!!!!!!矩阵列不满秩，无法分解!!!!!!!!!!!');
end
q = zeros(m,n); % 初始化Q R矩阵
r = zeros(n,n);
r(1,1) = norm(a(:,1)); % r(1,1) = norm([0,3,4]') = 5
q(:,1) = a(:,1)/r(1,1); % q1=a1/5
% q= 0 0 0 r = 5 0 0
% 0.6 0 0 0 0 0
% 0.8 0 0 0 0 0
for j = 2:n % 从第2列到第n列
for i = 1:j-1 % j=2时，从第2列第1行到第2列第1行
r(i,j) = q(:,i)' * a(:,j); % r(i,j) = <qi|aj>
% r(1,2) = -20*0+27*0.6+11*0.8 = 25
end
% 示例：j=2时结果
% q= 0 0 0 r = 5 25 0
% 0.6 0 0 0 0 0
% 0.8 0 0 0 0 0
qq = zeros(m,1);
for i = 1: j-1
qq = q(:,i) * r(i,j) + qq; % qq为qi<qi|aj>对i求和
end
r(j,j) = norm(a(:,j)-qq); % 对角线上元素||a2-q1<q1|a2>||
% [-20,12,-9]'的模长，25
q(:,j) = (a(:,j) - qq) / r(j,j); % [-20,12,-9]' * 1/25 作为q的第2列
% 示例：j=2时结果
% q= 0 -0.8 0 r = 5 25 0
% 0.6 0.48 0 0 25 0
% 0.8 -0.36 0 0 0 0
end
r(abs(r)<1e-12) = 0; % 置0
endHouseholder.m
function[q,r] = Householder(a)
% 考察教材QR factorization的定义，Q为mxn维，R为nxn上三角矩阵
% 为方便计算URV分解，Householder结果为full QR factorization
% Q为mxm维，R为mxn维
% a= 0 −20 −14
% 3 27 −4
% 4 11 −2
[m,n]=size(a);
q = eye(m);
for k=1:m-1 % 当k=2时第2轮循环
ak = a(k:m,k:n); % a = 5 25 −4 ak = 0 -10
% 0 0 −10 -25 -10
% 0 -25 −10
mk = m-k+1; % mk是ak的行数
e = zeros(mk, 1);
e(1,1) = 1; % e=[1,0]'
u = ak(:,1)-norm(ak(:,1))*e; % u=25*[-1,-1]'
R = eye(mk)-2*(u*u')/(u'*u);
I = eye(m);
I(k:m,k:m) = R; % R = 1 0
R = I; % 0 R
a = R*a; % …R3*R2*R1*A=T=r
q = R*q; % q= R3*R2*R1
end
r=a;
q=q'; % q= (R3*R2*R1)'
r(abs(r)<1e-12) = 0; % 置0
endGivens.m
function[q,r] = Givens(a)
% a= 0 −20 −14
% 3 27 −4
% 4 11 −2
[m,n]=size(a);
q = eye(m); % 初始化Q矩阵
for i = 1:n
if a(i,i)==0 % 当前位置元素为0时
p = eye(m);
[pivot, index] = max(abs(a(i:m,i))); % 寻找当前枢轴下最大主元的相对位置
index = index + i - 1; % index表示最大主元在矩阵中的行数索引
p([i,index],:) = p([index,i],:); % 将最大主元换到上面
a = p*a;
q = p*q;
end
for j = i+1:m
c = a(i,i)/(a(i,i)^2 + a(j,i)^2)^0.5;
s = a(j,i)/(a(i,i)^2 + a(j,i)^2)^0.5;
p = eye(m);
% p = 1 0 0
% 0 1 0
% 0 0 1
p([i,j],[i,j]) = [c,s;-s,c];
% p = c s 0
% -s c 0
% 0 0 1
a = p*a;
q = p*q;
end
if i == m % 当m<n时防止数组越界
break;
end
end
r = a;
q = q';
r(abs(r)<1e-12) = 0; % 置0
endURV.m
function[U,R,V] = URV(A)
[m,n] = size(A);
[Q1,R1] = Householder(A); % PA = (B 0)' 其中P=Q1',[B 0]'=R1
r = rank(R1);
B = R1(1:r, :);
[Q2,R2] = Householder(B'); % QB'=(T 0)' -> B=(T'|0)Q, Q=Q2'
r = rank(R2);
T = R2(1:r, 1:r);
U = Q1;
R = zeros(m,n);
R(1:r,1:r) = T';
V = Q2;
endDET.m
function[x] = DET(A)
% det(A)=det(P)*det(L)*det(U)
% L为主对角线全为1的下三角矩阵，det(L)=1
[m,n] = size(A);
if m ~= n
error('矩阵非方阵，无法求解！！');
end
% 进行LU分解
[~,U,P] = Factorization(A,0);
% 计算det(U)
% U为上三角矩阵,det(U)为主对角线元素之积
x = 1;
for i = 1:n
x = x*U(i,i);
end
% 计算det(P)
% P为单位矩阵经过行变换，det(P)= 1或-1
[kk,jj] = find(P==1); % 寻找P的元素为1的坐标，kk为所在行，jj为所在列
kk = kk';
% 计算kk的逆序数
count = 0;
for i = 1:n
for j = i:n
if kk(i) > kk(j)
count = count + 1;
end
end
end
if mod(count,2)==1 % 逆序数为奇数时
x = -x;
end
endAxb.m
% 利用LU分解求Ax=b
function[x] = Axb(A,b)
[m,n] = size(A);
if m ~= n
error('矩阵非方阵，无法求解！！');
end
if m ~= size(b,1)
error('A与b行数不同，无法求解！！');
end
x = zeros(n,1);
y = zeros(n,1);
% 对A进行LU分解
[L,U,P] = Factorization(A,0);
b = P*b;
% 求Ly=b
% 递推公式yi=bi-(求和yj*Lij)，求和j从1到i-1
for i = 1:n
y(i,1) = b(i,1);
for j = 1:i-1
y(i,1) = y(i,1)-y(j,1)*L(i,j);
end
end
% 求解Ux=y
% 递推公式xi=(yi-(求和yj*uij))/uii，求和j从i+1到n
for i = n:-1:1 % 从n到1
x(i,1) = y(i,1);
for j = i+1:n
x(i,1) = x(i,1)-x(j,1)*U(i,j);
end
x(i,1) = x(i,1)/U(i,i);
end
endFactorization.m
% method为选取的分解方法
% method=0 LU分解 [L,U,P] = Factorization(a, 0)
% method=1 GramSchmidt [Q,R] = Factorization(a, 1)
% method=2 Householder Reduction [Q,R] = Factorization(a, 2)
% method=3 Givens Reduction [Q,R] = Factorization(a, 3)
% method=4 URV [U,R,V] = Factorization(a, 4)
% 考察教材QR factorization的定义，Q为mxn维，R为nxn上三角矩阵
% GramSchmidt下的QR分解采取这种形式
% 为方便进行URV分解，Householder和Givens采取full QR factorization
% 即Q为mxm维，R为mxn维
function[a1,a2,a3] = Factorization(a, method)
if method==0
[a1,a2,a3] = LU(a);
elseif method==1
[a1,a2] = GramSchmidt(a);
elseif method==2
[a1,a2] = Householder(a);
elseif method==3
[a1,a2] = Givens(a);
elseif method==4
[a1,a2,a3] = URV(a);
else
fprintf('输入参数有误！')
end
endmain.m
clear all
clc
%% Faciorization函数实现不同的矩阵分解
% method为选取的分解方法
% method=0 LU分解 [L,U,P] = Factorization(a, 0)
% method=1 GramSchmidt [Q,R] = Factorization(a, 1)
% method=2 Householder Reduction [Q,R] = Factorization(a, 2)
% method=3 Givens Reduction [Q,R] = Factorization(a, 3)
% method=4 URV [U,R,V] = Factorization(a, 4)
% B = [0,-20,-40;3,27,-4;4,11,-2];
B = [1,2,4,17;3,6,-12,3;2,3,-3,2;0,2,-2,6];
[L,U,P] = Factorization(B, 0);
[Q1,R1] = Factorization(B, 1);
[Q2,R2] = Factorization(B, 2);
[Q3,R3] = Factorization(B, 3);
[u,r,v] = Factorization(B, 4);
% L,U,P
% Q1,R1
% Q2,R2
% Q3,R3
% u,r,v
%% Axb函数实现求解方程Ax=b
A = [1,2,4,17;3,6,-12,3;2,3,-3,2;0,2,-2,6];
b = [17;3;3;4];
x1 = Axb(A,b);
disp(x1);
% 输出Ax=b的解x
% 2
% -1
% 0
% 1
%% DET函数实现求矩阵行列式
A = [1,2,4,17;3,6,-12,3;2,3,-3,2;0,2,-2,6];
x2 = DET(A);
disp(x2);