MATLAB 矩阵分解算法大全

(1)LU 分解法程序：function x=solvebyLU(A,b)
% 该函数利用LU分解法求线性方程组Ax=b的解
flag=isexist(A,b); %调用第一小节中的isexist函数判断方程组解的情况if flag==0
disp('该方程组无解！');
x=[];
return;
else
r=rank(A);
[m,n]=size(A);
[L,U,P]=lu(A);
y(1)=b(1);
if m>1
for i=2:m
y(i)=b(i)-L(i,1:i-1)*y(1:i-1)';
end
end
y=y'; % 解Ux=y得原方程组的一个特解
x0(r)=y(r)/U(r,r);
if r>1
for i=r-1:-1:1
x0(i)=(y(i)-U(i,i+1:r)*x0(i+1:r)')/U(i,i);
end
end
x0=x0';
if flag==1 %若方程组有唯一解
x=x0;
return;
else %若方程组有无穷多解
format rat;
Z=null(A,'r'); %求出对应齐次方程组的基础解系
[mZ,nZ]=size(Z);
x0(r+1:n)=0;
for i=1:nZ
t=sym(char([107 48+i]));
k(i)=t; %取k=[k1,k2...,];
end
x=x0;
for i=1:nZ
x=x+k(i)*Z(:,i); %将方程组的通解表示为特解加对应齐次通解形式
end
end
end
(2)矩阵的QR分解法（c语言）：
void QR(double a[N][N],double q[N][N],double r1[N][N],int n) /*QR分解*/ {
int i,j,k,r,m;
double temp,sum,dr,cr,hr;
double ur[N],pr[N],wr[N];
double q1[N][N],emp[N][N];
for(i=1;i<n+1;i++)//将a放入temp中
for(j=1;j<n+1;j++)
{
emp[i][j]=a[i][j];
};
for(i=1;i<n+1;i++)//定义单位矩阵
for(j=1;j<n+1;j++)
{
if(i==j)q[i][j]=1;
else q[i][j]=0;
};
for(r=1;r<n;r++)
{
temp=0;
for(k=r+1;k<n+1;k++)
temp+=fabs(a[k][r]);
if(temp>=ZERO)
{
sum=0;
for(k=r;k<n+1;k++)
sum+=a[k][r]*a[k][r];
dr=sqrt(sum);
if(a[r][r]>ZERO)m=-1;
else m=1;
cr=m*dr;
hr=cr*(cr-a[r][r]);
for(i=1;i<n+1;i++)//定义ur
{
if(i<r)ur[i]=0;
if(i==r)
ur[i]=a[r][r]-cr;
if(i>r)
ur[i]=a[i][r];
};
for(i=1;i<n+1;i++)//定义wr
{
for(j=1;j<n+1;j++)
sum+=q[i][j]*ur[j];
wr[i]=sum;
};
for(i=1;i<n+1;i++)//定义qr
for(j=1;j<n+1;j++)
{
q1[i][j]=q[i][j]-wr[i]*ur[j]/hr;
};
for(i=1;i<n+1;i++)//定义qr+1
for(j=1;j<n+1;j++)
{
q[i][j]=q1[i][j];
};
for(i=1;i<n+1;i++)//定义pr
{
sum=0;
for(j=1;j<n+1;j++)
sum+=a[j][i]*ur[j];
pr[i]=sum/hr;
};
for(i=1;i<n+1;i++)
for(j=1;j<n+1;j++)
{
a[i][j]=a[i][j]-ur[i]*pr[j];
};
};
};
for(i=1;i<n+1;i++)
for(j=1;j<n+1;j++)
{
if(fabs(a[i][j])<ZERO)a[i][j]=0;
};
for(i=1;i<n+1;i++)
for(j=1;j<n+1;j++)
{
r1[i][j]=a[i][j];
};
for(i=1;i<n+1;i++)//将a取出
for(j=1;j<n+1;j++)
{
a[i][j]=emp[i][j];
}
(3)对称正定矩阵的Cholesky分解:
function [X]=m_chol(A,b)
[N,N]=size(A);
X=zeros(N,1);
Y=zeros(N,1);
for i=1:N
A(i,i)=sqrt(A(i,i)-A(i,1:i-1)*A(i,1:i-1)');
if A(i,i)==0
'A is singular. no unique solution'
break;
end;
for j=i+1:N
A(j,i)=(A(j,i)-A(j,1:i-1)*A(i,1:i-1)')/A(i,i);
end
end
A
b
%前代法
for j=1:N
Y(j)=(b(j)-A(j,1:j-1)*Y(1:j-1))/A(j,j);
end
Y;
A=A'
for k=N:-1:1
X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);
end。