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));
}

C++ 使用 householder 將矩陣分解成 QR 矩陣

使用 householdler 演算, 可以將任意矩陣分解成 Q R 兩矩陣相乘. 原理是先從矩陣擷取出行向量, 算出單位向量 ûₖ, 再計算對應的屋頂保持矩陣(Householder matrix) Hₖ, 該矩陣的特性是: 1. 兩個矩陣相乘會變成單位矩陣 HₖHₖ = I, 也就是說原矩陣就是反矩陣 Hₖ⁻ = Hₖ, 就好像是一面鏡子可分拆兩邊的特性.  2. 當一個矩陣左邊乘上 H 後, 行向量所在矩陣, 正對角以下元素全部清為 0, 矩陣屋頂(正對角)上的元素則保持不受影響, 也就是說乘一次 H 就清一行 0, 要注意的是 Hₖ 會隨著矩陣的更動而變化, 根據上述兩特性, 任一矩陣 A 就能分解成 Q 及 R 矩陣(A = QR):
    A = IA = I(HₖHₖ)A
    = I(H₁H₁)A = (IH₁) * (H₁A)
    = (IH₁)(H₂H₂)(H₁A) = (IH₁H₂) * (H₂H₁A)
    = (IH₁H₂)(H₃H₃)(H₂H₁A = (IH₁H₂H₃) * (H₃H₂H₁A)
    =  ...
    = (IH₁H₂H₃ ... Hₙ) * (Hₙ ... H₃H₂H₁A)
右式將矩陣 A 不斷的左乘 H, 隨著對角線以下, 一行一行清為零, 自然就變成上三角矩陣 R(又稱右矩陣) , 左式則將單位矩陣 I 不斷右乘 H 的結果, 最後演變成 Q 矩陣(它是正交矩陣, 因為 H 本身就是正交矩陣)
           Hₖ = I - 2*ûₖ*ûₖᵗ , 其中 ûₖ 為 Householder 單位向量
           R₁ = H₁A  = (I - 2*û₁*û₁ᵗ)*A   = A - 2*û₁*û₁ᵗ*A
           Q₁ = IH₁   = I(I - 2*Q*û₁*û₁ᵗ) =  I - 2*û₁*û₁ᵗ
           Rₖ <= HₖRₖ₋₁  = (I - 2*ûₖ*ûₖᵗ)*Rₖ₋₁    = Rₖ₋₁   -  2*ûₖ*ûₖᵗ*Rₖ₋₁
           Qₖ <= Qₖ₋₁Hₖ = Qₖ₋₁(I - 2*Q*ûₖ*ûₖᵗ) = Qₖ₋₁  -  2*Qₖ₋₁*ûₖ*ûₖᵗ
隨著 uₖ => Hₖ => Rₖ₊₁ => uₖ₊₁ => Hₖ₊₁ 不斷推陳出新, 利用疊代法就能算出最終的 Q, R:
           R -=  2*ûₖ*ûₖᵗ*R
           Q -= 2*Q*ûₖ*ûₖᵗ

#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;
}
F_vector householdvector(F_matrix &R, int i) {
    int m = R.size();    
    int mi = m - i;
    F_vector ut(mi);
    ut[0] = R[i][i];
    double s = 0;
    for (int jk = 1, j = i + 1; jk < mi ; jk ++) {
        double r = R[j ++][i];
        s += r * r;
        ut[jk] = r;
    }
    ut[0] += sqrt(s + ut[0]*ut[0]);
    s = sqrt(s + ut[0]*ut[0]);
    if (s < tol) return F_vector(0);
    for (int jk = 0; jk < mi ; jk ++) {
        ut[jk] /= s;
    }
    return ut;
}
int householdQRD(F_matrix &xx, F_matrix &qq, F_matrix &rr) {  // xx => QR
    int rank = 0;
    int m = xx.size();
    int n = xx[0].size();
    int maxSize = m;
    if (maxSize < n) maxSize = n;
    F_matrix Q = eye(m);
    F_matrix R = xx;
    F_matrix tempH(maxSize, F_vector(maxSize));
    F_vector uu(maxSize);
    for (int i = 0; i < n ; i ++) { // n columns
        F_vector ut = householdvector(R, i);
        int mi = ut.size();
        if (mi == 0) break;
        rank ++;
        for (int k = i, ki = 0; k < n; k ++, ki ++) {
            double sum = 0;
            for (int kj = 0, j = i; kj < mi; kj ++, j++) {
                sum +=  ut[kj] * R[j][k];      // uᵗ•Rj
            }
            uu[ki] = sum + sum;             // w = 2uᵗ*R, which is n elements of row vector
        }
        for (int ki = 0; ki < n - i; ki ++) {
            for (int kj = 0; kj < mi; kj ++) {
                tempH[kj][ki] = ut[kj] * uu[ki];// tempH = u*w = 2u*(uᵗ*R), which is a miXni matrix
            }
        }
        for (int k = i,  ki = 0; k < n; k ++, ki ++) {
            for(int j = i, kj = 0; j < m; j ++, kj ++) {
                R[j][k] -= tempH[kj][ki];         // R -= tempH = 2u*(uᵗ*R)  
            }
        }
        for (int j = 0; j < m; j ++) {
            double sum = 0;
            for (int ki = 0, k = i; ki < mi; ki ++) {
                sum += Q[j][k ++] * ut[ki];        // Qj•u
            }
            uu[j] = sum + sum;               // v = 2Q*u; It is a column vector with m elements now.
        }
        for (int j = 0; j < m ; j ++) {                    
            for (int ki = 0; ki < mi ; ki ++) {
                tempH[j][ki] = uu[j] * ut[ki];// tempH = v*uᵗ = 2Q*u*uᵗ; which is a mXm matrix
            }
        }
        for(int j = 0; j < m ; j++) {
            for (int ki = 0, k = i; ki < mi; ki ++, k ++) {
                Q[j][k] -= tempH[j][ki];      // Q -= tempH = (2Q*u)*uᵗ
            }
        }
    }
    swap(Q, qq);
    swap(R, rr);
    return rank;
}
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[] = {
          4,  2, 3, 4,
        17, 8, 9, 13
    };
    F_matrix B = arrayclone(a, 4, 2);
    F_matrix Q, R;
    matrixShow(B);
    householdQRD(B, Q, R);
    matrixShow(Q);
    matrixShow(R);
    matrixShow(product(Q, R));

    F_matrix C = transpose(B);
    matrixShow(C);
    householdQRD(C, Q, R);
    matrixShow(Q);
    matrixShow(R);
    matrixShow(product(Q, R));
}

linux 用 c++ 玩線性代數

在linux mint 下 c++ 使用 eigen3 作矩陣運算, 用 matplotlib-cpp 來繪圖, 可謂天作之合, 兩者都是標準純 C++ 標頭檔(header file), 所有程序都定義在類型內(class), 只要在程式裏面include 標頭檔, 並不需要聯結任何程式庫, 但編譯 eigen3 所寫的程式時,會很耗時, 要解決編譯的問題參考 https://stackoverflow.com/questions/67529835/eigen3-take-a-long-time-to-compile-and-very-slow-when-debug
就是把要用的程式先寫成 c++ 程式碼(.cpp), 編譯成獨立的模組, 作成程式庫, 後續透過聯結該模組就可. 首先上官網下載原始檔: https://eigen.tuxfamily.org/index.php?title=Main_Page
解壓縮到專案目錄內,只需用到目錄底下的 Eigen 目錄, 用符號聯結就可以了, 並不需要特別安裝
     ln   -sf   /home/mint/Downloads/eigen-3.4.0/Eigen   .
例如我只要用到 JacobiSVD, 但又不想用裡面的其他東西, 因習慣用 c++ 標準 <vector>, 只好先轉成 eigen3 的 matrix 才能運算, 結束時用 <vector> 傳回資料:
// jacobi.cpp
#include "Eigen/SVD"
#include <vector>
using namespace Eigen;
using namespace std;
using F_vector = vector<double>;
using F_matrix = vector<F_vector>;
extern "C" { // to export function
    F_vector jacobiSVD(F_matrix &x, F_matrix &U, F_matrix &V)  {
        int m = x.size();
        int n = x[0].size();
        MatrixXf temp (m, n);
        for (int i = 0; i < m; i++) {
            for (int j = 0; j < n; j++) { temp(i, j) = x[i][j]; }// copy x to temp
        }       
        JacobiSVD<MatrixXf> svd(temp, ComputeThinU | ComputeThinV );
        MatrixXf U2 = svd.matrixU();
        MatrixXf V2 = svd.matrixV();
        MatrixXf S2 = svd.singularValues();
        int ss = S2.rows();
        F_vector S(ss);
        U.resize(U2.rows(), F_vector(ss));
        V.resize(V2.rows(), F_vector(ss));
        for (int i = 0; i < ss; i ++)     { S[i]    = S2(i, 0); }// copy from S2
        for (int i = 0; i < U2.rows() ; i ++) {
            for (int j = 0; j < ss; j ++) { U[i][j] = U2(i, j); }// copy from U2
        }  
        for (int i = 0; i < V2.rows() ; i ++) {
            for (int j = 0; j < ss; j ++) { V[i][j] = V2(i, j); }// copy from V2
        }
        return S;// singular values in a vector
    }
}
直接用 g++ 編譯成模組檔 .o
    g++   jacobi.cpp   -c   -o   jacobi.o
後續只要聯結該模組無須重新編譯, 也完全不用再 include  任何 eigen3 的檔案, 解決編譯緩慢的問題, 但主程式必須要用 extern "C" 方式呼叫該模組內的函式:
#include <random>
#include <vector>
#include "matplotlibcpp.h"
#define tol ((double)1e-12)
namespace plt = matplotlibcpp;
using namespace std;
using F_vector = vector<double>;
using F_matrix = vector<F_vector>;
extern "C" {
    F_vector jacobiSVD(F_matrix &x, F_matrix &U, F_matrix &V);
}
double tpower(double t, int n) { // f(t, n) = tⁿ
    if (n == 0) return 1.0;
    if (n == 1) return t;
    double temp = t;
    for (int i = 1; i < n; i ++) temp *= t;
    return temp;
}
double absval(double t) { // f(t) = |t|
    return t < 0 ? -t : t;
}
struct svdSolver {
    F_matrix   U;
    F_vector   S;
    F_matrix   V;
    static svdSolver jacobiMethod(F_matrix &x) {
        F_matrix U;
        F_matrix V;
        F_vector S = jacobiSVD(x, U, V);// call external library  
        return svdSolver {.U = U, .S = S, .V = V}; // wrap U,S,V
    }
    static F_vector polyfit(F_vector y, int parameters = 2) {
        int size = y.size();
        F_matrix At(size, F_vector(parameters));
        for (int j = 0; j < parameters; j ++) {// Aᵗ xᵗ = yᵗ
            for (int i = 0; i < size ; i ++) {
                At[i][j] = tpower(i, j);// iᶨ
            }
        }
        svdSolver svd = jacobiMethod(At);
        svd.show(svd.S);// singular value
        svd.show(svd.V);// all of the right singular vectors
        return svd.svdfit(y);
    }
    void show(F_vector v) {
        for (int i = 0; i < v.size(); i ++) { printf("%f\t", v[i]); }
        printf("\n");
    }
    void show(F_matrix x) {
        for (int j = 0; j < x.size(); j ++) { show(x[j]); }
        printf("\n");
    }
    F_vector svdfit(F_vector &y) {// x̂ᵗ = A⁺ŷᵗ = (UΣVᵗ)⁻ŷᵗ = VΣ⁻Uᵗŷᵗ= V(Uᵗŷᵗ/Σ)
        int udim = U.size();
        if (y.size() != udim) return F_vector(0);
        int ss = S.size();
        F_vector uyisTemp(ss); // uyisTemp[j]ᵗ = Σ⁻(Ujᵗ•ŷᵗ) = Ujᵗ•ŷᵗ / Σ
        for (int j = 0; j < ss; j ++) {
            double sum = 0.0;// inner product
            if (absval(S[j]) > tol) {
                for (int i = 0; i < udim; i ++) { //  Ujᵗ•ŷᵗ             
                    sum += U[i][j] * y[i];
                }
                sum /= S[j];// (Ujᵗ • ŷᵗ)/Σ
            }
            uyisTemp[j] = sum;
        }
        int vdim = V.size();
        F_vector x(vdim);
        for (int i = 0; i < vdim; i ++) { // Vi • uyisTempᵗ
            double sum = 0.0;
            for(int j = 0; j < ss; j++) {
                sum += V[i][j] * uyisTemp[j];
            }
            x[i] = sum;
        }
        return x;
    }
};
int main() {  
    int fs = 512; // number of samples
    int np = 6;    // number of parameters to fit
    vector<double> xn(fs), yn(fs), gt(fs), cx(np);
    srand(time(nullptr));
    for (int i = 0 ; i < cx.size() ; i++) { cx[i] = rand(); }
    for (int t = 0 ; t < fs ; t ++)  {
        double st = 0;
        for (int i = 0; i < cx.size(); i ++) {
            st += tpower(t, i); // tᶤ
        }
        double err = (double)rand()/RAND_MAX - 0.5;
        double noise = st * err; // 50 % random noise
        yn[t] = st + noise;
        xn[t] = t;
    }
 
    vector<double> p = svdSolver::polyfit(yn, cx.size());
    for (int t = 0 ; t < fs ; t ++) {
        double yt = 0;
        for (int i = 0; i < p.size(); i ++) {
            yt += p[i] * tpower(t, i);
        }
        gt[t] = yt;
    }
    plt::plot(xn, yn, "g:x", xn, gt, "r-"); // x-y plot, red solid line.
    plt::show();
}
程式用到 c++ matplotlib, 參考:  https://masontseng.blogspot.com/2023/04/linum-mint-c-matplotlib.html  下載  matplotlibcpp.h
接著將主程式聯結模組成功後, 就能執行:
      g++  plot.c   jacobi.o   -std=c++11   -I/usr/include/python3.8    -lpython3.8  &&   ./a.out

矩陣向量一些概念

 用內積來表示兩向量 r̂, n̂ 之間的正投影長度 r̂•n̂ =|r̂||n̂|cosθ, 假設 r̂, n̂ 是 k 維空間的向量：
        r̂•n̂ = r₁*n₁ + r₂*n₂ + ... + rₖ*nₖ =  Σ(rₖ * nₖ) = Σ(nₖ*rₖ) = n̂•r̂
        |r̂| = √(r₁² + r₂² + ..., rₖ²)
        |n̂| = √(n₁² + n₂² + ..., nₖ²)
        cosθ = r̂•n̂/(|r̂||n̂|)

    想像空間是由一群向量 r̂₁, r̂₂, ... , r̂ₖ 等獨立基底向量 [r̂], 以原點 ô 為中心, 所展示的 k 維輻射空間,   描述任一向量 n̂ 是 r̂ 的線性組合, 用內積代表的座標, 就是在基底向量上的投影分量:
        n̂ = Σ nₖr̂ₖ = Σ(n̂•r̂ₖ)r̂ₖ = Σ(r̂ₖ•n̂)r̂ₖ
            n₁ = r̂₁•n̂
            n₂ = r̂₂•n̂
            .........
            nₖ = r̂ₖ•n̂

    底下符號 M 用來代表一個 nxn 方陣, 行向量彼此獨立, 而且列向量都是獨立向量
    底下符號 λ 用來代表 eigenvalue, 它是一個純量, 用一個數字來量化
    底下符號 ê 用來代表 eigenvector 它是一個向量, 用一組數字來量化(e₁, e₂,...,eₖ), 向量的元素個數 k, 用來描述該空間的維度, 用 [ê] 代表 eigenvector 行向量形成的矩陣 V
    底下符號 Λ 用來代表對角線矩陣, 除對角線上有數值外, 其餘位置都為 0
    底下符號 [λ] 用來代表由 eigenvalue 形成的對角線矩陣 Λ

    將 M 分解為由 eigenvalue 所組成對角線矩陣 Λ, 及 eigenvector 行向量 ê 所組成矩陣 V = [ê] 時: [V, D] = eig(M)
        M = V Λ V⁻ = [ê] [λ] [ê]⁻
            V = [ê₁, ê₂, ê₃, ..., êₖ]
            V⁻ 是 V 的反矩陣 = [ê]⁻
            Λ = D = [λ₁,0,...,0; 0,λ₂,0, ...,0; ... ; 0,0,0,...,λₖ] = [λ]
    就可以看出 M 是將向量先經反矩陣 V⁻ 轉換, 用 eigenvalue 伸縮座標軸, 投影在 eigenvector 上:
        M x̂ᵗ = (V Λ V⁻)x̂ᵗ
    eigenvalue λ , eigenvector ê 與矩陣 Ｍ 的關係式: M•ê = λ•ê
    當 Mᵢⱼ= Mⱼᵢ 時, M 就是一個對稱矩陣, 其 eigenvector 必定互相垂直, 由 eigenvector 組成矩陣的轉置矩陣就是反矩陣  
        ∵ Mᵢⱼ= Mⱼᵢ
        ∴ V⁻ = Vᵗ
        ∴ M  = V Λ Vᵗ
        M x̂ᵗ = (V Λ Vᵗ)x̂ᵗ = V Λ Vᵗ x̂ᵗ
   方陣轉置(將行列互換)後 eigenvalue 不會變
    任何一個矩陣轉置後, 再與原矩陣相乘就會形成一個對稱方陣:
        T = AᵗA = AₙxₘAₘxₙ= Tₙxₙ
    如果相乘的順序反過來也會形成另一個對稱方陣:
        N = AAᵗ =AₘxₙAₙxₘ= Nₘxₘ
    當矩陣不是方陣時, 任何一個 ₘxₙ的矩陣 A 可以透過 svd 分解成 3 個矩陣: [U, S, V] = svd(A)
    其中 S 是 singular value [σ] 形成的矩陣, V 是ＡᵗA 的 eigenvector 形成的矩陣 [v̂],  但 U 是 AAᵗ 的 eigenvector 形成的矩陣 [û], 要注意的是  AᵗA 與 AAᵗ 乘積是兩個不同的矩陣

              A = U S Vᵗ
        Aₘxₙ = Uₘxₘ • Sₘxₙ • Vᵗₙxₙ
        U 是由 AAᵗ 轉置矩的 eigenvector 行向量所組成 ₘxₘ 方陣 [û], m 是目標空間維度
            U = [û₁; û₂; û₃; ...; ûₘxₘ]ᵗ = [û₁, û₂, û₃, ..., ûₘxₘ] = [û]
        Uᵗ 是 U 的轉置矩陣也是 U 的反矩陣
            Uᵗ = [û₁, û₂, û₃, ..., ûₘ]⁻ = [û₁; û₂; û₃; ...; ûₘ] = [û]⁻
        S 對角線由 singular value σ 所組成, 其餘為 0, 組成一個 ₘxₙ 矩陣
            S =  [σ]ₘxₙ = [σ]
        V 是由 AᵗA 的 eigenvector 行向量所組成 ₙxₙ 方陣 [v̂], n 是來源空間維度
            V  = [v̂₁; v̂₂; v̂₃; ...; v̂ₙ]ᵗ = [v̂₁, v̂₂, v̂₃, ..., v̂ₙ] = [v̂]
        Vᵗ 是 V 的轉置矩陣也是 V 的反矩陣, n 是來源空間維度
            Vᵗ = [v̂₁, v̂₂, v̂₃, ..., v̂ₙ]⁻ = [v̂₁; v̂₂; v̂₃; ...; v̂ₙ] = [v̂]⁻

  備註: 只有方陣才能計算出 eigenvalue 與 eigenvector, 非方陣由 SVD 算出的是 singular value 矩陣及左右伴隨矩陣, 伴隨矩陣分別由 right singular vector及 left singular vector 行向量所組成的矩陣 V=[v̂], U=[û]
    1.矩陣 A = Aₘxₙ = Uₘxₘ•Sₘxₙ•Vᵗₙxₙ 若 singular value 只有 k 組 (m!=n, k=m 或是  k=n)時,  改用 Uₘxₖ•Sₖxₖ•Vᵗₖxₙ, 也就是說讓 S 組成稀疏對角線 ₖxₖ 方陣 S=Λ=[σ], U 變成瘦型矩陣 ₘxₖ, V 也順勢成瘦型矩陣 ₙxₖ, 較省記憶空間
    2.矩陣 AᵗA 的 eigenvalue 是 A singular value 的二次方:
        T = AᵗA = (U S Vᵗ)ᵗ U S Vᵗ = VSUᵗUSVᵗ = VS²Vᵗ = VΛVᵗ
        因此 A 的 singular value 平方就是 AᵗA eigenvalue, 反過來說 A 的 singular value 是 AᵗA eigenvalue 開根號
    3. AAᵗ 的 eigenvalue 也是  A的 singular value 二次方, 反過來說 A 的 singular value 是 AAᵗ eigenvalue 開根號
        AAᵗ= U S Vᵗ (U S Vᵗ)ᵗ = USVᵗ VSUᵗ= US²Uᵗ = UΛUᵗ
    4.矩陣 A 的 SVD 算法： Aₘxₙ = Uₘxₘ • Sₘxₙ • Vᵗₙxₙ
        a. 計算 T = AᵗA 用來找尋 A 的 singular value 及左右伴隨矩陣 UV
        b. Tv̂ = σ²v̂,透過 QR decoposition 解出 T 的 eigenvalue 矩陣 [σ²]ₘxₙ 及 eigenvector 所組成的矩陣 Vₙxₙ= [v̂]
        c. 將上述 eigenvalue 開根號得出 σ,  S = Λ = [σ] ₘxₙ, 將 σ 取倒數放入矩陣的對角線形成反矩陣,  注意反矩陣行列要交換才能與後續 Ｖ相乘
           S⁻ = Λ⁻ = [1/σ] ₙxₘ,
        d. U = U (Λ (Vᵗ V) Λ⁻) = (UΛVᵗ) V Λ⁻ = A V Λ⁻
           ＡV 先相乘, 最後再乘步驟 c. 的 singular value 反矩陣 Λ⁻, 就能算出 U = A • V • Λ⁻

c++ 反矩陣解矩陣方程式

#include <math.h>

#include <stdio.h>

#include <vector>

using namespace std;

using F_matrix = vector<vector<double>>;// matrix store raw vectors

void vectorShow(vector<double> v) {

for (int i = 0; i < v.size(); i ++) { printf("%10.4f", v[i]); }

printf("\n");

}

void matrixShow(F_matrix matrix) {

for (int j = 0; j < matrix.size(); j ++) { vectorShow(matrix[j]); }

printf("\n");

}

F_matrix arrayClone(double *array, int w, int h) {// copy array elements one by one

F_matrix temp(h);

for (int i = 0; i < h; i++) {

temp[i].resize(w);

for (int n = 0; n < w; n++) { temp[i][n]= *array ++; }

}

return temp;

}

F_matrix subMatrix(F_matrix &mat, int size, int xcol, int xrow = 0) {// remove xrow and xcol

int nn = mat.size();

if (size <= 0 || size > nn) return F_matrix(0);

F_matrix temp(size, vector<double>(size));

for (int j = 0, xj = 0; j < nn; j ++) {// row scan, skip 1st row

if (j == xrow) continue; // skip row xrow

for (int i = 0, xi = 0; i < nn; i++) { // column scan

if (i == xcol) continue;// skip column xcol

temp[xj][xi ++] = mat[j][i];// copy element one by one,

}

xj ++;

}

return temp;

}

double determinant(F_matrix &x, int s) {

if (s == 0) return 0.0;

int n = x.size();

if (n != x[0].size()) return 0;

if (n == 1) return x[0][0];

if (n == 2) return x[0][0] * x[1][1] - x[0][1] * x[1][0];

int det = 0;

bool positive = true;

for (int k = 0; k < n; k++) {

F_matrix subXk = subMatrix(x, n - 1, k);// remove column k, row 0

double temp = x[0][k] * determinant(subXk, n - 1);// recursive

det += positive ? temp : -temp;

positive = ! positive; // toggle the sign

}

return det;

}

double determinant(F_matrix &x) { return determinant(x, x.size()); }

F_matrix inverse(F_matrix &x) {// cramer's rule [x]⁻ = adj([x])ᵗ/det([x])

double det = determinant(x);

if (fabs(det) < 1e-16) return F_matrix(0);

int xdim = x.size();

F_matrix inv(xdim, vector<double>(xdim));

int k = xdim - 1; // to find cofactor

for (int j = 0; j < xdim; j ++) {

for (int i = 0; i < xdim; i ++) {

F_matrix xijSub = subMatrix(x, k, i, j);// remove col i and row j

double cji = determinant(xijSub); // cofactor of x[j][i]

if ((j + i) % 2) cji = - cji; // sign if odd

inv[i][j] = cji / det;// transpose and divdide determinant at the same time

}

}

return inv;

}

vector<double> xSolve(F_matrix &A, vector<double> b) { // x̂ᵗ = A⁻¹b̂ᵗ

int n = A.size();

vector<double> ans(A.size());

if (n == b.size()) {

F_matrix invA = inverse(A);

if (invA.size() > 0) {

for(int j = 0; j < n; j++) {

double temp = 0;

for(int i=0; i < n; i++) { // Aj⁻¹ • b̂ᵗ

temp += invA[j][i] * b[i];

}

ans[j] = temp;

}

}

}

return ans;

}

int main(int argc, char **argv) {

double a[] = { // Axᵗ = bᵗ, xᵗ = [w,y,z,u]ᵗ, bᵗ = [1,3,4,5]ᵗ, A:

1 , 2, 3, 4, // 1w + 2y + 3z + 4u = 1

4 , 5, 6, 2, // 4w + 5y + 6z + 2u = 3

17, 8, 9, 3, // 17w + 8y + 9z + 3u = 4

14, 7, 5, 6 // 14w + 7y + 5z + 6u = 5

};

vector<double> b{1, 3, 4, 5};

F_matrix A = arrayClone(a, 4, 4);

printf("\nSolve equation: Ax̂ᵗ = b̂ᵗ\n");

matrixShow(A);

vectorShow(b);

printf("\nAnswer is x̂ᵗ = A⁻¹b̂ᵗ:\n");

vector<double> x = xSolve(A, b);

vectorShow(x);

}