C++ 使用 Givens 旋轉矩陣將矩陣分解 QR

將矩陣分解成 QR 兩矩陣相乘, 可以利用一個類似旋轉矩陣,  一次清除一個位置, 這種矩陣稱為 Givens 旋轉矩陣 Gⱼ(θ), 有別於用屋頂保持矩陣 H(û) (Householder matrix)一次清除一行,  Gⱼ(θ) 是在單位矩陣對角線上, I(j,j) 的位置附近附加 2x2 旋轉角 θ 的元素 (cosθ, sinθ), 用來將位置 (j,i) 的元素清為 0, 它有旋轉矩陣的特性: Gⱼ 轉置矩陣 Gⱼᵗ 就是反矩陣 Gⱼ⁻ = Gⱼᵗ, 也就是說 GⱼᵗGⱼ = I, 因此只要將對角線以下的元素一個一個清為零, 自然就化解出 QR 矩陣:
           A = IA = I(GⱼᵗGⱼ)A  = (IGⱼᵗ ...) * (... Gⱼ)A
               = Q * R
要注意的是 Gⱼ(θ) 與矩陣的更動及要清 0 位置 (j, i) 的旋轉角 θ 有關, 需要動態更新,上述右式隨著矩陣 A 不斷左乘 Gⱼ 矩陣, 對角線下的元素一一清為  0, 最終自然浮現出上三角矩陣 R, 左式伴隨著將單位矩陣不斷右乘 Gⱼᵗ 矩陣, 最終演化成 Q 矩陣.  因此可以用下面疊代法,算出最終的 Q, R:
           初始化 Q =  I, R = A, Gⱼ(θ)=[cosθ, sinθ ; -sinθ, cosθ]
           R =  Gⱼ(θ) * R
           Q = Q * Gⱼᵗ(θ)
其中 Gⱼ(θ) 附加的 2x2 矩陣 Gⱼ(θ)=[cosθ, sinθ; -sinθ, cosθ], Gᵗ(θ)= [cosθ, -sinθ; sinθ, cosθ].
Givens 運算適合用在下三角位置上 0 較多的矩陣. 等於 0 時就不用計算. 實際上A[j][i] = 0 時,  Gⱼ(θ) 的旋轉角正好是 90ᵒ.  而 Gⱼ(90ᵒ) 等於單位矩陣 I, 運算就是多餘的. 當矩陣與 Givens 矩陣相乘時,  其實只會影響到附近 2x2  整行與整列計算, 有很多乘 0 的部份其實可以省略, 底下只闡述實作,  並未作最佳化:
#include <stdio.h>
#include <math.h>
#include <vector>
#define tol ((double)1e-12)
using namespace std;
using F_vector = vector<double>;
using F_matrix = vector<F_vector>;
void vectorShow(F_vector v) {
    for (int i = 0; i < v.size(); i ++) { printf("%-10.4f   \t", v[i]); }
    printf("\n");
}
void matrixShow(F_matrix matrix) {
    for (int j = 0; j < matrix.size(); j ++) {
        printf("<%-2d:>\t", j);
        vectorShow(matrix[j]);
    }
    printf("\n");
}

F_matrix transpose(F_matrix &matrix) { //  矩陣轉置
    int row = matrix.size();
    int col = matrix[0].size();
    F_matrix temp(col, F_vector(row));
    for (int i = 0; i < row; i++) {
        for (int j = 0; j < col; j ++) {  
            temp[j][i] = matrix[i][j];
        }
    }
    return temp;
}
F_matrix eye(int n) {// 單位矩陣
    F_matrix temp(n, F_vector(n));
    for (int i = 0; i < n; i ++) {
        for (int j = 0; j < n; j++) {
            temp[i][j] = (i == j);
        }
    }
    return temp;
}
F_matrix product(F_matrix q, F_matrix r) { // 矩陣相乘
    int dim = q[0].size();
    if (dim != r.size()) return F_matrix(0);
    int m = q.size();
    int n = r[0].size();
    F_matrix x(m, F_vector(n));
    for (int i = 0; i < m; i ++) {
        for (int j = 0; j < n; j ++) {
            double sum = 0;
            for (int k = 0; k < dim; k++) {
                sum += q[i][k] * r[k][j];
            }
            x[i][j] = sum;
        }
    }
    return x;
}

void givensQRD(F_matrix &A, F_matrix &qq, F_matrix &rr) {    // A = QR
    int m = A.size();
    int n = A[0].size();
    F_matrix Q = eye(m);
    F_matrix R = A;
    double c = 1;// dummy, 預設 θ = 90ᵒ, c = cos(πθ/180) = 1
    double s = 0;// dummy, 預設 θ = 90ᵒ, s = sin(πθ/180) = 0
    for (int i = 0; i < n; i ++) {// 從第一行開始掃描
        for (int j = m - 1; j > i; j --) {// 行向量 i 從最後一個元素開始往前直到對角線, 一一清為 0
                double a = R[j][i]; // 目標行向量 i,  位置 j 元素準備清為 0, 先計算出旋轉矩陣 Gⱼᵢ(θ)
                if (fabs(a) < tol) continue; // 當 a != 0 才要清除
                double b = R[j-1][i];// 行向量 i 只剩位置 j-1 的元素值會受影響
                if (fabs(a) > fabs(b)) { // 找出  Gⱼ(θ) 旋轉, 只要推導出  cosθ, sinθ, 不用算出實際 θ 值
                    double t = b / a;
                    s = 1 / sqrt(1 + t * t);
                    c = s * t;
                } else {
                    double t = a / b;
                    c = 1 / sqrt(1 + t * t);
                    s = c * t;
                }
                F_matrix G = eye(m);// 初始化旋轉矩陣 Gⱼ(θ) = I,  旋轉角 θ = 90ᵒ
                G[j][j] = c;// 更新旋轉矩陣,準對角線上擺 cosθ, 目標清除行向量 i 位置 j 元素 a => 0
                G[j-1][j-1] = c;// 更新旋轉矩陣, 準對角線左上擺 cosθ
                G[j][j-1] = -s;// 更新旋轉矩陣,準對角, 正左擺 -sinθ
                G[j-1][j] = s;// 更新旋轉矩陣,準對角, 正上擺 sinθ
                R = product(G, R); // 乘完後 R[j[i] => 0
                G[j][j-1] = s;// 矩陣轉置 Gⱼ => Gⱼᵗ
                G[j-1][j] = -s;// GⱼᵗGⱼ = I
                Q = product(Q, G);// A = QGⱼᵗ GⱼR,            
        }
    }
    swap(Q, qq);
    swap(R, rr);
}
F_matrix arrayclone(double *array, int w, int h) {
    F_matrix temp(h, F_vector(w));
    for (int i = 0; i < h; i++) {
        for (int j = 0; j < w; j++) { temp[i][j]= *array ++; }
    }
    return temp;
}
int main(int argc, char **argv) {
    double a[] = {
        1,  2, 3,  4, 5,
        6,  7, 8,  9,10,
        17, 8, 7, 12, 4,
        17, 8, 1, 11, 5,
         7, 2, 5,  3, 1
    };
    F_matrix B = arrayclone(a, 5, 5);
    F_matrix C = transpose(B);
    F_matrix Q, R;
 
    matrixShow(B);
    givensQRD(B, Q, R);
    matrixShow(R);
    matrixShow(Q);
    matrixShow(product(Q,R));
 
    matrixShow(C);
    givensQRD(C, Q, R);
    matrixShow(R);
    matrixShow(Q);
    matrixShow(product(Q,R));
}