用豪斯霍尔德（Householder）变换进行矩阵的QR分解，及其Matlab和OpenCV实现

1、豪斯霍尔德变换
一般地，对给定的$m$维向量$a$，考虑分块 $a=\left[\begin{array}{c}{a}_1\\ a_2\end{array}\right]$，其中$a_1$是$(k-1)$维向量，$1\le k<m$。如果豪斯霍尔德向量为$v=\left[\begin{array}{c}0\\ a_2\end{array}\right]-\alpha e_k$，其中$\alpha =-sign(a_k)||a_2|{|}_2$，$e_k$是单位矩阵的第$k$列，则由此形成的豪斯霍尔德变换可消去$a$的后$m-k$个分量。对$k=1,2,\cdots ,n$用这种方法可定义一系列豪斯霍尔德变换，以消去$m\times n$阶矩阵$A$的副对角线及以下的元素，对$A$的列按从左到右的顺序依次作上述变换，就可将矩阵约化为上三角阵。每个豪斯霍尔德变换只对$A$的未被约化的部分起作用，并不影响前面已约化好的部分，因而在依次变换过程中，已得到的零元素能得以保持。对任意向量$u$，经过豪斯霍尔德变换$H$后，为 $Hu=\left(I-2\frac{v{v}^T}{v^{T}v}\right)u=u-\left(2\frac{v^{T}u}{v^{T}v}\right)v$，这个变换的计算量比通常的矩阵向量相乘要少得多，且只要知道向量$v$即可，不必具体算出矩阵$H$。
2、伪代码

3、Matlab实现

%% Matlab 2018b——用豪斯霍尔德变换进行矩阵的QR分解:A→R（上三角阵）
A=[1,0,0;0,1,0;0,0,1;-1,1,0;-1,0,1;0,-1,1];
% b=[1237,1941,2417,711,1177,475];
% A=[A,b'];
[m,n]=size(A);
alpha=zeros(1,n);
beta=alpha;
gama=alpha;
v=zeros(m,n);
E=eye(m,n);
for i=1:n
   alpha(i)=-1*sign(A(i,i))*sqrt(sum(A(i:m,i).*A(i:m,i)));
   v(:,i)=[zeros(i-1,1);A(i:m,i)]-alpha(i).*E(:,i);
   beta(i)=v(:,i)'*v(:,i);
   if beta(i)==0
      continue; 
   end
   for j=i:n
      gama(j)=v(:,i)'*A(:,j);
      A(:,j)=A(:,j)-(2*gama(j)/beta(i)).*v(:,i);  
   end
end
R=A;

4、OpenCV实现

// OpenCV 4.5.0——用豪斯霍尔德变换进行矩阵的QR分解：A→R（上三角阵）
#include <iostream>
#include <opencv2/opencv.hpp>

using namespace std;
using namespace cv;

int main()
{
	Mat A = (Mat_<float>(6, 3) << 1, 0, 0, 0, 1, 0, 0, 0, 1, -1, 1, 0, -1, 0, 1, 0, -1, 1);
	Size n = A.size();
	Mat alpha = Mat::zeros(1, n.width, CV_32FC1);
	Mat V = Mat::zeros(n, CV_32FC1);
	Mat E = Mat::eye(n, CV_32FC1);
	Mat beta = alpha.clone(), gama = alpha.clone();
	Mat temp, temp1;	// 过程暂存变量
	for (int i = 0; i < n.width; i++)
	{		
		temp = A(Range(i, n.height), Range(i, i + 1)).clone();
		temp1 = temp.mul(temp);		// 各项求平方
		reduce(temp1, temp1, 0, REDUCE_SUM);	// 平方和
		alpha.at<float>(i) = -1 * A.at<float>(i, i) / fabsf(A.at<float>(i, i)) * sqrt(temp1.at<float>(0));
		vconcat(Mat::zeros(i, 1, CV_32FC1), temp, temp);
		V.col(i) = temp - alpha.at<float>(i)*E.col(i);
		temp1 = V.col(i).t()*V.col(i);
		beta.at<float>(i) = temp1.at<float>(0);
		if (beta.at<float>(i) == 0)
			continue;
		for (int j = i; j < n.width; j++)
		{
			temp = V.col(i).t()*A.col(j);
			gama.at<float>(j) = temp.at<float>(0);
			A.col(j) = A.col(j) - (2 * gama.at<float>(j) / beta.at<float>(i))*V.col(i);
		}
	}
	cout << A << endl;
	getchar();
	return 0;
}

参考文献：《科学计算导论》（第2版）——Michael T. Health 著 张威等 译