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