矩陣QR分解的MATLAB與C++實現

一：矩陣QR分解

矩陣的QR分解目的是將一個列滿秩矩陣\(A\)分解成\(A=QR\)的形式，我們這里暫時討論\(A\)為方陣的情況。其中\(Q\)為正交矩陣；\(R\)為正線（主對角線元素為正）上三角矩陣，且分解是唯一的。

比如\(A= \begin{bmatrix} 1 & 2 & 2 \\ 2 & 1 & 2 \\ 1 & 2 & 1 \\ \end{bmatrix}\)，我們最終要分解成如下形式：

\[A=Q \cdot R = \begin{bmatrix} \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{3}} & \frac{1}{\sqrt{2}} \\ \frac{2}{\sqrt{6}} & -\frac{1}{\sqrt{3}} & 0 \\ \frac{1}{\sqrt{6}} & \frac{1}{\sqrt{3}} & -\frac{1}{\sqrt{2}} \\ \end{bmatrix} \cdot \begin{bmatrix} \sqrt{6} & \sqrt{6} & \frac{7\sqrt{6}}{6} \\ 0 & \sqrt{3} & \frac{\sqrt{3}}{3} \\ 0 & 0 & \frac{\sqrt{2}}{2} \\ \end{bmatrix} \]

現在主要的問題是如何由矩陣\(A\)計算得到矩陣\(Q\)和\(R\)呢？我們將在下面討論。

1.1 QR分解原理

在線性代數或矩陣理論中，我們肯定都學過斯密特正交化（Gram-Schmidt Orthogonalization）,正交化過程即將歐氏空間的任一基化為標准正交基，構造出的標准正交基正好構成了我們想要的\(Q\)矩陣，而\(R\)矩陣由正交化過程的公式倒推即可得到。

首先假設初始方陣為\(A\)，\(\vec{x_i}\)、\(\vec{y_i}\)、\(\vec{z_i}\)都為列向量。我們學過斯密特正交化的步驟如下：

\[A=\begin{bmatrix} \vec{x_1} & \vec{x_2} & \vec{x_3} \end{bmatrix} \overset{正交化}{\underset{}{\to}} \begin{bmatrix} \vec{y_1} & \vec{y_2} & \vec{y_3} \end{bmatrix} \overset{單位化}{\underset{}{\to}} \begin{bmatrix} \vec{z_1} & \vec{z_2} & \vec{z_3} \end{bmatrix} = Q \]

再具體一點（為了好寫，之后的\(\vec{x_i}\)、\(\vec{y_i}\)、\(\vec{z_i}\)都不加箭頭了，默認為列向量）：

\[y_k = x_k - \sum_{i=1}^{k-1} \frac{(x_k,y_i)}{(y_i,y_i)}y_i = x_k - \sum_{i=1}^{k-1} \frac{(x_k,y_i)}{||y_i||^2}y_i = x_k - \sum_{i=1}^{k-1} (x_k,z_i)z_i \tag{1} \]

\[z_k = \frac{y_k}{||y_k||} ,k=1...n \tag{2} \]

\[Q = \begin{bmatrix} z_1 & \cdots & z_n \tag{3} \end{bmatrix} \]

\[R= \begin{bmatrix} ||y_1|| & (x_2,z_1) & \cdots & (x_n,z_1) \\ & ||y_2|| & \cdots & (x_n,z_2) \\ & & \ddots & \vdots\\ \mathsf 0 & & &||y_n|| \end{bmatrix} \tag{4} \]

由上述公式寫出計算\(Q\)和\(R\)的偽代碼為：

\[\begin{align} & for \quad k=1:n \notag\\ & \qquad R_{kk}=||A_{:k}|| \notag\\ & \qquad Q_{:k}=A_{:k} / R_{kk} \notag\\ & \qquad for \quad i = k + 1 : n \notag\\ & \qquad \qquad R_{ki} = A_{:i}' * Q_{:k} \notag\\ & \qquad \qquad A_{:i} = A_{:i} - R_{ki} .* Q_{:k} \notag\\ & \qquad end \notag\\ & end \notag\\ \end{align} \]

注：\(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\)矩陣完全相同。