矩阵的QR分解求解线性方程组

矩阵的QR分解求解线性方程组

一.QR分解概念

QR分解指的把矩阵分解成一个正交矩阵与一个上三角矩阵的积。
正交矩阵:如果n阶方阵A满足 A A T = A T A = E ( 即 A − 1 = A T ) {AA^T=A^TA=E(即A^{-1}=A^T)} AAT=ATA=E(即A−1=AT),那么称A为正交矩阵
上三角矩阵:
矩阵的QR分解求解线性方程组
A = Q R {A=QR} A=QR,Q是正交矩阵,R是上三角矩阵。

二.使用HouseHolder变换来实现QR分解

1.共轭转置

对于一个矩阵 A {A} A,它的共轭转置我们记为 A ∗ {A^*} A∗,他们之间满足如下关系:
( A ∗ ) i , j = A j , i ‾ (A^*)_{i,j}=\overline{A_{j,i}} (A∗)i,j​=Aj,i​​,当然也有其他等价形式:
A ∗ = A T ‾ = ( A ‾ ) T {A^*}=\overline{A^T}=({\overline A})^T A∗=AT=(A)T,那么我们这里给出一个例子:
A = [ 3 + i 5 2 − 2 i i ] A= \left[ \begin{matrix} 3+i& 5\\ 2-2i & i \end{matrix} \right] A=[3+i2−2i​5i​]
A ∗ = [ 3 − i 2 + 2 i 5 − i ] A^*= \left[ \begin{matrix} 3-i& 2+2i\\ 5& -i \end{matrix} \right] A∗=[3−i5​2+2i−i​]

2.HouseHolder变换实现QR分解

对于任意的一个模长为1的列向量w,我们构造豪斯霍尔德矩阵:
H = I − 2 w w T {H=I-2ww^T} H=I−2wwT,其中 I I I为单位列向量。这里 w T {w^T} wT没使用共轭转置,是因为目前我们先只考虑实数的形式。从这个形式我们可以得到下面一般性的结论"

H T = H {H^T=H} HT=H
H ∗ H = I {H*H=I} H∗H=I

这说明 H {H} H是一个对称正交阵。实际上H也是一个反射矩阵。下面我们来说说为什么?
首先介绍超平面,指和 w {w} w垂直(即乘积为0)的所有向量 x x x组成的平面。
记作:

S = { x ∣ w T x = 0 , ∀ x ϵ R n } S=\{x|w^Tx=0, \forall x\epsilon R^n \} S={x∣wTx=0,∀xϵRn}

基于这样的定义,我们很容易明白,对于一个任意的向量 z ϵ R n z\epsilon R^n zϵRn,向S与w张成的这个空间上一定可以进行唯一的分解。
记作 z = x + y , x z=x+y,x z=x+y,x在超平面 S S S上,而 y = α w , α y=\alpha w,\alpha y=αw,α是一个实数。
那么有:

H z = x + ( I − 2 w w T ) ∗ y = x + y − 2 α w ( w T w ) = x − y Hz=x+(I-2ww^T)*y=x+y-2\alpha w(w^Tw)=x-y Hz=x+(I−2wwT)∗y=x+y−2αw(wTw)=x−y
这就很明显了, H H H相当于将 z z z按照S平面进行了依次反射对称。这就是H为反射矩阵的由来。

这里我们介绍一个定理:

对于任意的模相等的向量 x , y x,y x,y,我们都可以找到一个反射矩阵 H H H,使得 H x = y , H = I − 2 w w T , w = x − y ∣ ∣ x − y ∣ ∣ Hx=y,H=I-2ww^T,w=\frac{x-y}{||x-y||} Hx=y,H=I−2wwT,w=∣∣x−y∣∣x−y​
这个定理我们的用法是这样的,给定一个向量 x x x,它的模长是 l l l,那么我们可以找到 H H H使得 H x = [ ± l , 0 , 0 , . . . , 0 ] T Hx=[\pm l,0,0,...,0]^T Hx=[±l,0,0,...,0]T

矩阵的QR分解求解线性方程组

3.QR分解解线性方程组实现代码

//
int DSPF_dp_qrd_solver(const int Nrows,const int Ncols,double *restrict Q,double *restrict R,double *restrict b,double *restrict y,double *restrict x) {

  short row,col,loop_cnt;
  double sum;
  double xx,yy;

  _nassert(Nrows>0);
  _nassert(Ncols>0);
  _nassert((int)Q % 8 == 0);
  _nassert((int)R % 8 == 0);
  _nassert((int)b % 8 == 0);
  _nassert((int)y % 8 == 0);
  _nassert((int)x % 8 == 0);

  /* generate y=Q'*b 见第一张图*/
#pragma MUST_ITERATE(1,,)
  for (row=0;row<Nrows;row++) {
	sum=0.0;
#pragma MUST_ITERATE(1,,)
	for (col=0;col<Nrows;col++) {
	  sum+=Q[row+col*Nrows]*b[col];
	} 
	y[row]=sum;
  }

  /* use backward substitution to solve x=inv(R)*y 见第二张图*/
  if (Nrows>=Ncols) loop_cnt=Ncols;
  else              loop_cnt=Nrows;
  memset(x,0,sizeof(double)*Ncols);
#pragma MUST_ITERATE(1,,)
  for (row=loop_cnt-1;row>=0;row--) {
	sum=0.0;
	xx=R[row+row*Ncols];
	yy=_rcpdp(xx);//yy是xx的倒数
	yy=yy*(2.0-xx*yy);
	yy=yy*(2.0-xx*yy);
	yy=yy*(2.0-xx*yy);
#pragma MUST_ITERATE(1,,)
	//模拟高斯消元的方法 
	for (col=row+1;col<loop_cnt;col++) {
	  sum+=R[col+row*Ncols]*x[col];//如果有无穷多组解时最下面的几个就赋值为0
	}
	x[row]=(y[row]-sum)*yy;
  }
  return 0;
}

矩阵的QR分解求解线性方程组
矩阵的QR分解求解线性方程组

上一篇:【Ybtoj 第5章 例题1】走迷宫【广搜】


下一篇:HTML5 拖放事件与dataTransfer对象 ---- 2021/07/19