One: matrix QR decomposition

The QR decomposition of A matrix is to decompose A full rank matrix AAA into the form A=QRA=QRA=QR. Let’s discuss the case where AAA is A square matrix for the moment. QQQ is the orthogonal matrix; RRR is the upper triangular matrix of the positive line (the principal diagonal element is positive) and the decomposition is unique.

For example, A=[122212121]A= begin{bmatrix} 1&2&2\2&1&2\1&2&1 \\ end{bmatrix}A=⎣⎢, 121212221 nation ⎥⎤.


[ 1 6 1 3 1 2 2 6 – 1 3 0 1 6 1 3 – 1 2 ] ⋅ [ 6 6 7 6 6 0 3 3 3 0 0 2 2 ] \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}

Now the main question is how do you compute the matrices QQQ and RRR from matrix AAA? We will discuss this below.

1.1 QR decomposition principle

In linear algebra or matrix theory, we must have learned Schmidt Orthogonalization, which converts any basis of a Euclidean space to an orthonormal basis, producing an orthonormal basis that forms the desired QQQ matrix, The RRR matrix can be obtained by the inverse of the orthogonalization formula.

First, assume that the initial square matrix is AAA, xi⃗\vec{x_i}xi, Yi ⃗\vec{y_i}yi, zi⃗\vec{z_i}zi are all column vectors. We learned the steps of Schmidt’s orthogonalization as follows:


A = [ x 1 ⃗ x 2 ⃗ x 3 ⃗ ] – orthogonalization [ y 1 ⃗ y 2 ⃗ y 3 ⃗ ] – Basis, [ z 1 ⃗ z 2 ⃗ z 3 ⃗ ] = Q A = \ begin {bmatrix} \ vec {x_1} & \ vec {x_2} & \ vec {x_3} {bmatrix} \ \ end overset {orthogonalization} {\ underset {} {\ to}} \ begin {bmatrix} & \ vec \ vec {y_1} {y_2} & \ vec {y_3} {bmatrix} \ \ end overset {unit change} {\ underset {} {\ to}} \ begin {bmatrix} \ vec {z_1} & \ vec & {z_2} \vec{z_3} \end{bmatrix} = Q

To be more specific, xi⃗\vec{x_i}xi, yi⃗\vec{y_i}yi, zi⃗\vec{z_i}zi do not include arrows, and default to column vectors:


y k = x k – ∑ i = 1 k – 1 ( x k . y i ) ( y i . y i ) y i = x k – ∑ i = 1 k – 1 ( x k . y i ) ∣ ∣ y i ∣ ∣ 2 y i = x k – ∑ i = 1 k – 1 ( x k . z i ) z i (1) 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 = y k ∣ ∣ y k ∣ ∣ . k = 1… n (2) z_k = \frac{y_k}{||y_k||} ,k=1… n \tag{2}

Q = [ z 1 … z n ] (3) Q = \begin{bmatrix} z_1 & \cdots & z_n \tag{3} \end{bmatrix}

R = [ ∣ ∣ y 1 ∣ ∣ ( x 2 . z 1 ) … ( x n . z 1 ) ∣ ∣ y 2 ∣ ∣ … ( x n . z 2 ) ⋱ ⋮ 0 ∣ ∣ y n ∣ ∣ ] (4) 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}

The pseudocode of QQQ and RRR calculated by the above formula is:


f o r k = 1 : n R k k = ∣ ∣ A : k ∣ ∣ Q : k = A : k / R k k f o r i = k + 1 : n R k i = A : i ‘ ∗ Q : k A : i = A : i – R k i . ∗ Q : k e n d e n d \begin{aligned} & for \quad k=1:n \\ & \qquad R_{kk}=||A_{:k}|| \\ & \qquad Q_{:k}=A_{:k} / R_{kk} \\ & \qquad for \quad i = k + 1 : n \\ & \qquad \qquad R_{ki} = A_{:i}’ * Q_{:k} \\ & \qquad \qquad A_{:i} = A_{:i} – R_{ki} .* Q_{:k} \\ & \qquad end \\ & end \\ \end{aligned}

Note: A:kA_{:k}A:k represents the KKK column vector of AAA.

It can be seen that there are not many steps of QR decomposition of the matrix, that is, the constant cycle of AAA orthogenization, standardization, QQQ, RRR these steps.

Two: matrix QR decomposition MATLAB implementation

clc, clear all, close all

QR decomposition of % matrix
A = [1 2 2;2 1 2;1 2 1] Consider a nonsingular square matrix
[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)); % computes the diagonal element of R
    Q(:,k) = A(:,k) / R(k,k); % A has been orthogonalized, and now if I normalize it, I get the orthogonal matrix Q
    for i = k + 1 : n
        R(k,i) = A(:,i)' * Q(:,k); % computes the upper triangle of R
        A(:,i) = A(:,i) - R(k,i) .* Q(:,k); Update matrix A, Schmidt orthogonal formula
    end
end
Q
R
Copy the code

Three: matrix QR decomposition of C++ implementation

#include <iostream>
#include <vector>
using namespace std;

int main(a) /* QR factorization of matrix A */
{
    vector<vector<double>> a = { {1.2.2}, {2.1.2}, {1.2.1}};int n = a.size(a); vector<vector<double>> q(n, vector<double>(n));
    vector<vector<double>> r(n, vector<double>(n));

    cout << "A:" << endl; // Output matrix 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); / / calculate the first k A column length, by the formula (4) is equal to the diagonal elements of R | | A: k | |
        for (int i = 0; i < n; i++)
        {
            q[i][k] = a[i][k] / r[k][k]; // From formula (2), the KTH column of A becomes the KTH column of Q
        }

        for (int i = k + 1; i < n; i++)
        {
            for (int j = 0; j < n; j++)
            {
                r[k][i] += a[j][i] * q[j][k]; // From formula (4), calculate the upper triangle of R
            }
            for (int j = 0; j < n; j++)
            {
                a[j][i] -= r[k][i] * q[j][k]; // From formula (1), calculate and update each column of A
            }
        }
    }

    cout << endl;
    cout << "Q:" << endl; // Output matrix 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; // output matrix R
    for (int i = 0; i < n; i++)
    {
        for (int j = 0; j < n; j++)
        {
            printf("%.4f ", r[i][j]);
        }
        cout << endl;
    }

    return 0;
}
Copy the code

Four: comparison of results

As can be seen from the following figure, the QQQ and RRR matrices calculated by MATLAB and C++ are exactly the same.