电子工程数学方法-Matlab Doolittle Court Householder Givens

1.基于Givens变换QR分解Matlab代码

function [Q,R]=qrgv(A)
% 基于Givens变换，将方阵A分解为A=QR，其中Q为正交矩阵，R为上三角阵
%
% 参数说明
% A：需要进行QR分解的方阵
% Q：分解得到的正交矩阵
% R：分解得到的上三角阵
%
% 实例说明
% A=[-12 3 3;3 1 -2;3 -2 7];
% [Q,R]=qr(A) % 调用MATLAB自带的QR分解函数进行验证
% [q,r]=qrgv(A) % 调用本函数进行QR分解
% q*r-A % 验证 A=QR
% q'*q % 验证q的正交性
% norm(q) % 验证q的标准化，即二范数等于1
%
% 线性代数基础知识
% 1.B=P*A*inv(P)，称A与B相似，相似矩阵具有相同的特征值
% 2.Q*Q'=I，称Q为正交矩阵，正交矩阵的乘积仍为正交矩阵
%
% by dynamic of Matlab技术论坛
% see also http://www.matlabsky.com
% contact me matlabsky@gmail.com
% 2010-01-17 22:51:18
%
n=size(A,1);
R=A;
Q=eye(n);
for i=1:n-1
    for j=2:n-i+1
        x=R(i:n,i);
        rt=givens(x,1,j);
        r=blkdiag(eye(i-1),rt);
        Q=Q*r';
        R=r*R;
    end
end

function [R,y]=givens(x,i,j)
% 求解标准正交的Given变换矩阵R，使用Rx=y，其中y(j)=0，y(i)=sqrt(x(i)^2+x(j)^2)
%
% 参数说明
% x：需要进行Givens变换的列向量
% i：变为sqrt(x(i)^2+x(j)^2)的元素下标
% j：变为0的元素的下标
% R：Givens变换矩阵
% y：Givens变换结果
%
% 实例说明
% x=[1 3 5 9 6]'; % 将3等效到9上
% [R,y]=givens(x,4,2) % 注意3的下标为2，9的下标为4
% R*x-y % 验证Rx=y
% R'*R % 验证正交性
% norm(R) % 验证标准性，就是范数为1
%
% 关于Givens变换说明
% 1.Givens矩阵是标准正交矩阵，也叫平面旋转矩阵，它是通过坐标旋转的原理将元素j的数值等效到元素i上
% 2.Givens变换每次只能将一个元素变为0，而Householder变换则一次可以将任意个元素变为0
% 3.Givens变换常用于将矩阵A变为对角阵
%
xi=x(i);
xj=x(j);
r=sqrt(xi^2+xj^2);
cost=xi/r;
sint=xj/r;
R=eye(length(x));
R(i,i)=cost;
R(i,j)=sint;
R(j,i)=-sint;
R(j,j)=cost;
y=x(:);
y([i,j])=[r,0];

但是我跑不出来？不知道为什么

原作者：【原创】基于Givens变换QR分解Matlab代码|MATLAB 数学统计与优化|MATLAB技术论坛 - Powered by Discuz! (matlabsky.com)

2. Doolittle分解

function [L,U] = myLU(A)
    n=length(A);
    L=eye(n,n);
    U=zeros(n,n);
    U(1,1:n)=A(1,1:n);
    L(2:n,1)=A(2:n,1)./U(1,1);

    for r=2:n
        i=r:n;
        U(r,i)=A(r,i)-L(r,1:r-1)*U(1:r-1,i);
        i=r+1:n;
        L(i,r)=(A(i,r)-L(i,1:r-1)*U(1:r-1,r))/U(r,r);
    end
end

 

原作者：Doolittle分解（matlab代码）_I_Like_Algorithms的博客-CSDN博客

 

3.Householder QR分解

 

A=[];
[M,N]=size(A);
%获得矩阵维数
A1=A;
H1=zeros(M,M);
for j=1:M
 H1(j,j)=1;
end 
%k表示对所有的列
for k=1:N
 %设置H矩阵初值，这里设置为单位矩阵
  H0=zeros(M,M);
  %设置为单位矩阵
  for i=1:M
  H0(i,i)=1;
  end 
  s=0;
  for i=k:M
   s=s+A1(i,k)*A1(i,k);
  end 
%算的向量的2范数  
  s=sqrt(s);
  u=zeros(N,1);
%---------------------------------
% 这段甚为重要，关系到数值稳定性问题
% 目的是使得u的2范数尽可能大
% 原则是，如果首元素大于零，则u的第一个元素是正+正
%         如果首元素小于零，则u的第一个元素是负-负
  if (A1(k,k)>=0) 
  u(k)=A1(k,k)+s;
  else 
  u(k)=A1(k,k)-s;
  end 
%-------------------------------  
  for i=k+1:M
  u(i)=A1(i,k);
  end   
  du=0;
  for i=k:M
  %求的单位u 长度平方
  du=du+u(i)*u(i);
  end 
  %计算得到大的H矩阵
  for i=k:M
    for j=k:M
        H0(i,j)=-2*u(i)*u(j)/du;            
        if i==j 
        H0(i,j)=1+H0(i,j);
        end
    end 
  end 
%千万要注意矩阵相乘的次序
%先更新矩阵
A2=H0*A1;
A1=A2;
%H1初值为单位矩阵，后逐步更新
H1=H1*H0;
%更新变换后的矩阵
end 
Q=H1;
R=A1;

算出来的结果 符号可能不一样 

原作者：

Householder Transformation 进行QR分解Matlab代码__VioletHan_的博客-CSDN博客_householder transformations matlab

-----------另一种

function [Q,R]=qrhs(A)
% 基于Householder变换，将方阵A分解为A=QR，其中Q为正交矩阵，R为上三角阵
%
% 参数说明
% A：需要进行QR分解的方阵
% Q：分解得到的正交矩阵
% R：分解得到的上三角阵
%
% 实例说明
% A=[-12 3 3;3 1 -2;3 -2 7];
% [Q,R]=qr(A) % 调用MATLAB自带的QR分解函数进行验证
% [q,r]=qrhs(A) % 调用本函数进行QR分解
% q*r-A % 验证 A=QR
% q'*q % 验证q的正交性
% norm(q) % 验证q的标准化，即二范数等于1
% 
% 线性代数基础知识
% 1.B=P*A*inv(P)，称A与B相似，相似矩阵具有相同的特征值
% 2.Q*Q'=I，称Q为正交矩阵，正交矩阵的乘积仍为正交矩阵
% 
% 注意：我们也可以基于Givens变换，对方阵A进行QR分解，但是相对繁琐些，参见http://www.matlabsky.com/thread-4850-1-1.html
%
% by dynamic of Matlab技术论坛
% see also http://www.matlabsky.com
% contact me matlabsky@gmail.com
% 2010-01-17 22:49:51
%
n=size(A,1);
R=A;
Q=eye(n);
for i=1:n-1
    x=R(i:n,i);
    y=[1;zeros(n-i,1)];
    Ht=householder(x,y);
    H=blkdiag(eye(i-1),Ht);
    Q=Q*H;
    R=H*R;
end

function [H,rho]=householder(x,y)
% 求解正交对称的Householder矩阵H，使Hx=rho*y，其中rho=-sign(x(1))*||x||/||y||
%
% 参数说明
% x：列向量
% y：列向量，x和y必须具有相同的维数
%
% 实例说明
% x=[3 5 6 8]';
% y=[1 0 0 0]';
% [H,rho]=householder(x,y);
% H*x-rho*y % 验证Hx=rho*y
% H'*H % 验证正交
%
% 关于HouseHolder变换
% 1.H=I-2vv'，其中||v||=1，我们称H为反射(HouseHolder)矩阵，易证H对称且正交
% 2.如果||x||=||y||，那么存在HouseHolder矩阵H=I-2vv，其中v=±(x-y)/||x-y||，使Hx=y
% 3.如果||x||≠||y||，那么存在HouseHolder矩阵H，使Hx=rho*y，其中rho=-sign(x(1))*||x||/||y||
% 4.Householder矩阵，常用于将一个矩阵A通过正交变换对角阵
%
x=x(:);
y=y(:);
if length(x)~=length(y)
    error('The Column Vectors X and Y Must Have The Same Length!');
end
rho=-sign(x(1))*norm(x)/norm(y);
y=rho*y;
v=(x-y)/norm(x-y);
I=eye(length(x));
H=I-2*v*v';

    

同样跑不出来不知道为什么

 原作者：[原创]基于Householder变换QR分解Matlab代码|MATLAB 数学统计与优化|MATLAB技术论坛 - Powered by Discuz! (matlabsky.com)

4.Court分解

仅适用于3*3

 傻瓜代码 无任何参考价值 但是可以减少手算

A=[];
L=zeros(3);
U=eye(3,3);
L(1)=A(1);
L(2)=A(2);
L(3)=A(3);
U(4)=A(4)/L(1);
U(7)=A(7)/L(1);
L(5)=A(5)-L(2)*U(4);
L(6)=A(6)-L(3)*U(4);
U(8)=(A(8)-L(2)*U(7))/L(5);
L(9)=A(9)-L(3)*U(7)-L(6)*U(8);

 

 

收集整理 仅供个人学习使用