矩阵QR 分解的MATLAB 与C++实现
⼀：矩阵QR 分解
矩阵的QR 分解⽬的是将⼀个列满秩矩阵A 分解成A =QR 的形式，我们这⾥暂时讨论A 为⽅阵的情况。其中Q 为正交矩阵；R 为正线（主对⾓线元素为正）上三⾓矩阵，且分解是唯⼀的。
⽐如A =1
printf输出格式matlab
22
2
121
2
1
，我们最终要分解成如下形式：A =Q ⋅R =1
√6
1
√31
√22
√6
−1
√
3
1
√61
√3
−1
√
2⋅
√6√6
7√66
√3
√3
3
√2
2
现在主要的问题是如何由矩阵A 计算得到矩阵Q 和R 呢？我们将在下⾯讨论。
1.1 QR 分解原理
在线性代数或矩阵理论中，我们肯定都学过斯密特正交化（Gram-Schmidt Orthogonalization ），正交化过程即将欧⽒空间的任⼀基化为标准正交基，构造出的标准正交基正好构成了我们想要的Q 矩阵，⽽R 矩阵由正交化过程的公式倒推即可得到。⾸先假设初始⽅阵为A ，→x i 、→y i 、→
z i 都为列向量。我们学过斯密特正交化的步骤如下：
A =→x 1→x 2→x 3
正交化→→y 1→y 2→y 3单位化→→z 1→z 2→
z 3=Q 再具体⼀点（为了好写，之后的→x i 、→y i 、→
z i 都不加箭头了，默认为列向量）：
y k =x k −k −1∑
i =1
(x
k ,y i )
(y i ,y i )
y i =x k −k −1
∑
i =1
(x k ,y i )||y i ||2
y i =x k −k −1
∑
i =1(x k ,z i )z i
z k =
y k
||y k ||
,
k =1...n
Q =z 1⋯
z n
R =
||y 1||(x 2,z 1)⋯(x n ,z 1)||y 2||
⋯(x n ,z 2)⋱
⋮0
||y n ||
由上述公式写出计算Q 和R 的伪代码为：
[]
[][]
[][][]
[]
[
]
for k=1:n
R kk=||A:k||
Q:k=A:k/R kk
for i=k+1:n
R ki=A′:i∗Q:k
A:i=A:i−R ki.∗Q:k
end
end
注：A:k表⽰A的第k列向量。
可以看出其实矩阵的QR分解的步骤并不多，就是不断地循环进⾏A的正交化、标准化、求Q、求R这⼏步。⼆：矩阵QR分解的MATLAB实现
clc, clear all, close all
% 矩阵的QR分解
A = [1 2 2;2 1 2;1 2 1] % 考虑⾮奇异⽅阵
[m,n] = size(A);
Q = zeros(n,n);
X = zeros(n,1);
R = zeros(n);
for k = 1 : n
R(k,k) = norm(A(:,k)); % 计算R的对⾓线元素
Q(:,k) = A(:,k) / R(k,k); % A已正交化，现在做标准化，得到正交矩阵Q
for i = k + 1 : n
R(k,i) = A(:,i)' * Q(:,k); % 计算R的上三⾓部分
A(:,i) = A(:,i) - R(k,i) .* Q(:,k); % 更新矩阵A，斯密特正交公式
end
end
Q
R
三：矩阵QR分解的C++实现
#include <iostream>
#include <vector>
using namespace std;
int main() /* 矩阵A的QR分解*/
{
vector<vector<double>> a = { {1,2,2},{2,1,2},{1,2,1} };
int n = a.size();
vector<vector<double>> q(n, vector<double>(n));
vector<vector<double>> r(n, vector<double>(n));
cout << "A:" << endl; //输出矩阵A
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
printf("%.4f ", a[i][j]);
}
cout << endl;
}
for (int k = 0; k < n; k++)
{
double MOD = 0;
for (int i = 0; i < n; i++)
{
MOD += a[i][k] * a[i][k];
}
r[k][k] = sqrt(MOD); // 计算A第k列的模长，由公式(4)等于R的对⾓线元素||A:k||
for (int i = 0; i < n; i++)
{
q[i][k] = a[i][k] / r[k][k]; // 由公式(2)，A第k列标准化之后成为Q的第k列
}
for (int i = k + 1; i < n; i++)
{
for (int j = 0; j < n; j++)
{
r[k][i] += a[j][i] * q[j][k]; // 由公式(4)，计算R的上三⾓部分
}
for (int j = 0; j < n; j++)
{
a[j][i] -= r[k][i] * q[j][k]; // 由公式(1)，计算更新A的每⼀列
}
}
}
cout << endl;
cout << "Q:" << endl; //输出矩阵Q
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
printf("%.4f ", q[i][j]);
}
cout << endl;
}
cout << endl;
cout << "R:" << endl; //输出矩阵R
for (int i = 0; i < n; i++)
{
for (int j = 0; j < n; j++)
{
printf("%.4f ", r[i][j]);
}
cout << endl;
}
return 0;
}
四：结果对⽐
由下图可以看到，由MATLAB和C++计算出的Q和R 矩阵完全相同。Processing math: 100%