数值线性代数课程设计

《数值线性代数》课程设计

题目

学生

指导教师

徐州师范大学

课程设计任务书

院班学生

课程设计课题：

矩阵的QR 分解

课程设计任务要求（包括课题来源、类型、目的和意义、基本要求、参考资料等）：

✧ 来源与意义:

本课题来源于教材第三章第三节——正交化算法，目的是通过QR 分解求解方程Ax=b ，以求解最小二乘问题。

✧ 基本要求: 1. 要求自编程序；

2. 掌握编程思想，学会一门编程语言；

3. 报告要有较强的理论分析；有较强说服力的数据表或图像；

4. 对结果进行分析；

5. 给出相应结论。

✧ 参考资料:

《数值线性代数》：徐树芳 高立 张平文

指导教师签字：

一、 问题提出：

在探究最小二乘问题时，我们得出结论，求解最小二乘问题等价于求

2min ||||T T Q Ax Q b 的问题，其中Q 为正交矩阵，因此，我们希望通过选取适当的正交矩阵Q ，使原问题转化为较易求解的最小二乘问题。这里，我们将讨论对任意A 的QR 分解的方法。

二、理论基础

1、QR 分解定理：设()m n A R m n ⨯∈≥，则A 有QR 分解：

0R A Q ⎡⎤=⎢⎥⎣⎦

，

其中m m Q R ⨯∈是正交矩阵，n n R R ⨯∈是具有非负对角元的上三角阵；而且当m n =且A 非奇异时，上述分解还是唯一的。

3）说明：在这个定理的基础上，我们就可以放心的研究任意矩阵的QR 分解，因为它总是存在的。

2、Householder 变换：

1) 定义：设n R ω∈满足21ω=，定义n n H R ⨯∈为

2T H I ωω=-，

则称H 为Householder 变换。

2）作用：H x ⨯可以将12

n x x x x ⎡⎤

⎢⎥

⎢⎥=⎢⎥

⎢⎥

⎣⎦ 化成(1)100x x ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦

的形式。

3）说明：利用Householder 变换，我们可以将A 的对角线下元素逐列的消为0。

3、Givens 变换

1）定义：形如G=11i c s

j s c i j

⎛⎫

⎪

⎪ ⎪ ⎪ ⎪

⎪- ⎪

⎪ ⎪⎝⎭

的矩阵称为Givens 矩阵。

2）作用：G x ⨯可以将12

n x x x x ⎡⎤

⎢⎥

⎢⎥=⎢⎥

⎢⎥

⎣⎦

化成0

0i x j ⎡

⎤⎢⎥⎢⎥⎢⎥=⎢⎥

⎢⎥⎢⎥

⎢⎥

⎢⎥⎣⎦

的形式。

3）说明：利用Givens 变换，我们可以将A 的对角线下元素逐个消为0。

三、实验内容

模型一：利用householder变换对矩阵A进行QR分解。

1.用MATLAB建立m文件house.m,程序如下：function [v,b]=house(x)

n=length(x);

v=zeros(n,1);

eta=norm(x,inf);

x=x/eta;

v(1)=1;

v(2:n)=x(2:n);

sig=x(2:n)'*x(2:n);

if (sig==0)

b==0;

else

a=(x(1)*x(1)+sig)^(1/2);

if (x(1)<=0)

v(1)=x(1)-a;

else

v(1)=-sig/(x(1)+a);

end

b=2*v(1)*v(1)/(sig+v(1)*v(1));

v=v/v(1);

end

2.用MATLAB建立m文件：qr1.m

A=input('A=');

[m,n]=size(A);

for j=1:n

if j<=m

[v,b]=house(A(j:m,j))

A(j:m,j:n)=(eye(m-j+1)-b*v*v')*A(j:m,j:n); d(j)=b;

A(j+1:m,j)=v(2:m-j+1);

end

end

R=A(1:n,1:n);

R=triu(R);

Q=eye(m);

for j=1:n

A(j,j)=1;

end

for j=1:n

H=eye(m);

H(j:m,j:m)=eye(m-j+1)-d(j)*A(j:m,j)*A(j:m,j)';

Q=Q*H;

end

Q

R

Q*Q' %这里是为了检验Q是否为正交矩阵

3.结果：利用Householder变换，我们可以将A的对角线下元素逐列的消为0，从而得到R，每一步中的Householder变换矩阵H的乘积就是我们所求的Q。（程序运行的结果数据见附录一）。

模型二：利用givens变换对矩阵A进行QR分解。

1.用MATLAB建立m文件givens.m,程序如下：

function[c,s]=givens (a,b)

if (b==0)

c=1;

s=0;

else if(abs(b)>abs(a))

r=a/b;

s=1/sqrt(1+r*r);

c=s*r;

else

r=b/a;

c=1/sqrt(1+r*r);

s=c*r;

end

end

2. 用MATLAB建立m文件：qr2.m

A=input('A=');

[m,n]=size(A);

Q=eye(m);

for j=1:n

for i=j+1:m

[c,s]=givens(A(j,j),A(i,j));