Back to Documentations

Signature Description
struct  PCAParams  {

    normalization_type  norm_type { normalization_type::z_score };

    // If populated (set above zero), number of top eigen values to keep.
    //
    long                num_comp_to_keep { 0 };

    // If populated (num_comp_is 0), percentage of eigenvalues to keep. 0.9 means 90%.
    //
    double              pct_comp_to_keep { 0.9 };
};
A structure containing the parameteres to pca_by_eigen() call

Signature Description Parameters
template<typename T>
Matrix<T, matrix_orient::column_major>
pca_by_eigen(std::vector<const char *> &&col_names,
             const PCAParams params = { }) const;
This uses Eigenspace evaluation to calculate Principal Component Analysis (PCA). It returns a matrix whose columns are the reduced dimensions with most significant information.

PCA is a dimensionality reduction method that is often used to reduce the dimensionality of large data sets, by transforming a large set of variables into a smaller one that still contains most of the information in the large set.
Reducing the number of variables of a data set naturally comes at the expense of accuracy, but the trick in dimensionality reduction is to trade a little accuracy for simplicity. Because smaller data sets are easier to explore and visualize, and thus make analyzing data points much easier and faster for machine learning algorithms without extraneous variables to process.
The final projection of data is done with the original un-normalized data. This is implementing a whitening or reconstruction step where you intentionally want scores in original units.

This works with both scalar and multidimensional (MD), vectors and arrays, data. In case of MD data, the data is flattened row-wise into a scalar matrix.
T: Type of the named columns
col_names: Vector of column names
params: Parameters necessary for for this operation
template<typename T>
std::tuple<Matrix<T, matrix_orient::column_major>, // U
           Matrix<T, matrix_orient::column_major>, // Σ
           Matrix<T, matrix_orient::column_major>> // V
compact_svd(std::vector<const char *> &&col_names,
            normalization_type norm_type =
                normalization_type::z_score) const;
This calculates Singular Value Decomposition (SVD). Optionaly it may normalize the original matrix first.
In linear algebra, SVD is a factorization of a real or complex matrix into a rotation, followed by a rescaling followed by another rotation. It generalizes the eigen-decomposition of a square normal matrix with an orthonormal eigenbasis to any ⁠mXn matrix.
It returns the 3 metrices U, Σ, and V inside a std::tuple.

U contains the left singular vectors of the original matrix, meaning its columns are orthonormal vectors that span the row space of the matrix.
Σ is a diagonal matrix that contains sqrt of eigenvalues of the original matrix, arranged in descending order.
V contains the right singular vectors of the original matrix, represented as its columns.
Original matrix (A) = U * Σ * VT

This works with both scalar and multidimensional (MD), vectors and arrays, data. In case of MD data, the data is flattened row-wise into a scalar matrix.
T: Type of the named columns
col_names: Vector of column names
norm_type: Type of normalization applied to raw data first
static void test_pca_by_eigen()  {

    std::cout << "\nTesting pca_by_eigen( ) ..." << std::endl;

    StrDataFrame    df;

    try  {
        df.read("IBM.csv", io_format::csv2);
    }
    catch (const DataFrameError &ex)  {
        std::cout << ex.what() << std::endl;
        ::exit(-1);
    }

    const auto  pca_mat = df.pca_by_eigen<double>({ "IBM_Close", "IBM_Open", "IBM_High", "IBM_Low" });

    // Dimensions were reduced to 1 containing at least 90% of the information.
    // This makes sense, since these 4 columns are highly correlated.
    //
    assert(pca_mat.cols() == 1);
    assert(pca_mat.rows() == 5031);
    assert(std::fabs(pca_mat(0, 0) - 197.063) < 0.001);
    assert(std::fabs(pca_mat(1, 0) - 200.875) < 0.001);
    assert(std::fabs(pca_mat(491, 0) - 149.02) < 0.01);
    assert(std::fabs(pca_mat(1348, 0) - 166.44) < 0.01);
    assert(std::fabs(pca_mat(2677, 0) - 333.405) < 0.001);
    assert(std::fabs(pca_mat(5029, 0) - 216.175) < 0.001);
    assert(std::fabs(pca_mat(5030, 0) - 219.555) < 0.001);

    const auto  pca_mat2 = df.pca_by_eigen<double>({ "IBM_Close", "IBM_Open", "IBM_High", "IBM_Low" }, { .num_comp_to_keep = 3 });

    // 3 most significant dimensions are kept.
    // As you can see the first column is unchanged and clearly contains
    // almost all of the information.
    //
    assert(pca_mat2.cols() == 3);
    assert(pca_mat2.rows() == 5031);

    assert(std::fabs(pca_mat2(0, 0) - 197.063) < 0.001);
    assert(std::fabs(pca_mat2(0, 1) - -0.0951913) < 0.001);
    assert(std::fabs(pca_mat2(0, 2) - 1.85473) < 0.001);

    assert(std::fabs(pca_mat2(1, 0) - 200.875) < 0.001);
    assert(std::fabs(pca_mat2(1, 1) - -2.08604) < 0.001);
    assert(std::fabs(pca_mat2(1, 2) - 2.68895) < 0.001);

    assert(std::fabs(pca_mat2(491, 0) - 149.02) < 0.01);
    assert(std::fabs(pca_mat2(491, 1) - -1.34957) < 0.01);
    assert(std::fabs(pca_mat2(491, 2) - 2.09026) < 0.01);

    assert(std::fabs(pca_mat2(1348, 0) - 166.44) < 0.01);
    assert(std::fabs(pca_mat2(1348, 1) - 0.0354559) < 0.01);
    assert(std::fabs(pca_mat2(1348, 2) - 0.41972) < 0.01);

    assert(std::fabs(pca_mat2(2677, 0) - 333.405) < 0.001);
    assert(std::fabs(pca_mat2(2677, 1) - -1.33686) < 0.001);
    assert(std::fabs(pca_mat2(2677, 2) - 2.13684) < 0.001);

    assert(std::fabs(pca_mat2(5029, 0) - 216.175) < 0.001);
    assert(std::fabs(pca_mat2(5029, 1) - -1.18141) < 0.001);
    assert(std::fabs(pca_mat2(5029, 2) - 2.18029) < 0.001);

    assert(std::fabs(pca_mat2(5030, 0) - 219.555) < 0.001);
    assert(std::fabs(pca_mat2(5030, 1) - -2.66858) < 0.001);
    assert(std::fabs(pca_mat2(5030, 2) - 2.85412) < 0.001);

    // Now multidimensional data
    //
    using ary_col_t2 = std::array<double, 2>;
    using ary_col_t3 = std::array<double, 3>;
    using vec_col_t = std::vector<double>;

    // 1) Variance is almost entirely in the x-component
    //
    std::vector<ary_col_t2> ary_xvar_1  { { 1.0, 0.1 }, { 4.0, 0.2 }, { 7.0, 0.1 }, { 2.0, 0.0 }, { 3.0, 0.2 }, { 1.0, 0.1 } };
    std::vector<ary_col_t2> ary_xvar_2  { { 2.0, 0.2 }, { 5.0, 0.0 }, { 8.0, 0.3 }, { 4.0, 0.1 }, { 6.0, 0.1 }, { 3.0, 0.0 } };
    std::vector<ary_col_t2> ary_xvar_3  { { 3.0, -0.1 }, { 6.0, 0.1 }, { 9.0, -0.2 }, { 6.0, 0.2 }, { 9.0, -0.1 }, { 5.0, 0.1 } };

    std::vector<vec_col_t>  vec_xvar_1  { { 1.0, 0.1 }, { 4.0, 0.2 }, { 7.0, 0.1 }, { 2.0, 0.0 }, { 3.0, 0.2 }, { 1.0, 0.1 } };
    std::vector<vec_col_t>  vec_xvar_2  { { 2.0, 0.2 }, { 5.0, 0.0 }, { 8.0, 0.3 }, { 4.0, 0.1 }, { 6.0, 0.1 }, { 3.0, 0.0 } };
    std::vector<vec_col_t>  vec_xvar_3  { { 3.0, -0.1 }, { 6.0, 0.1 }, { 9.0, -0.2 }, { 6.0, 0.2 }, { 9.0, -0.1 }, { 5.0, 0.1 } };

    df.load_column<ary_col_t2>("ARY XVAR 1", std::move(ary_xvar_1), nan_policy::dont_pad_with_nans);
    df.load_column<ary_col_t2>("ARY XVAR 2", std::move(ary_xvar_2), nan_policy::dont_pad_with_nans);
    df.load_column<ary_col_t2>("ARY XVAR 3", std::move(ary_xvar_3), nan_policy::dont_pad_with_nans);

    df.load_column<vec_col_t>("VEC XVAR 1", std::move(vec_xvar_1), nan_policy::dont_pad_with_nans);
    df.load_column<vec_col_t>("VEC XVAR 2", std::move(vec_xvar_2), nan_policy::dont_pad_with_nans);
    df.load_column<vec_col_t>("VEC XVAR 3", std::move(vec_xvar_3), nan_policy::dont_pad_with_nans);

    const auto  ary_xvar_res = df.pca_by_eigen<ary_col_t2>({ "ARY XVAR 1", "ARY XVAR 2", "ARY XVAR 3" }, { .num_comp_to_keep = 1 });
    const auto  vec_xvar_res = df.pca_by_eigen<vec_col_t>({ "VEC XVAR 1", "VEC XVAR 2", "VEC XVAR 3" }, { .num_comp_to_keep = 1 });

    assert(ary_xvar_res.cols() == 1);
    assert(ary_xvar_res.rows() == 6);
    assert(vec_xvar_res.cols() == 1);
    assert(vec_xvar_res.rows() == 6);

    assert(std::abs(ary_xvar_res(0, 0) - -2.98144) < 0.00001);
    assert(std::abs(ary_xvar_res(3, 0) - -5.6686) < 0.0001);
    assert(std::abs(ary_xvar_res(5, 0) - -4.2447) < 0.0001);
    assert(std::abs(vec_xvar_res(0, 0) - -2.98144) < 0.00001);
    assert(std::abs(vec_xvar_res(3, 0) - -5.6686) < 0.0001);
    assert(std::abs(vec_xvar_res(5, 0) - -4.2447) < 0.0001);

    // 2) x and y components carry equal but orthogonal variance
    //
    std::vector<ary_col_t2> ary_xyvar_1  { { 1.0, 4.0 }, { 2.0, -4.0 }, { 3.0, 4.0 }, { 4.0, -4.0 }, { 1.5, 4.0 }, { 2.5, -4.0 }, { 3.5, 4.0 }, { 4.5, -4.0 } };
    std::vector<ary_col_t2> ary_xyvar_2  { { 2.0, 3.0 }, { 4.0, -3.0 }, { 6.0, 3.0 }, { 8.0, -3.0 }, { 3.0, 3.0 }, { 5.0, -3.0 }, { 7.0, 3.0 }, { 9.0, -3.0 } };
    std::vector<ary_col_t2> ary_xyvar_3  { { 3.0, 2.0 }, { 6.0, -2.0 }, { 9.0, 2.0 }, { 12.0, -2.0 }, { 4.5, 2.0 }, { 7.5, -2.0 }, { 10.5, 2.0 }, { 13.5, -2.0 } };
    std::vector<ary_col_t2> ary_xyvar_4  { { 4.0, 1.0 }, { 8.0, -1.0 }, { 12.0, 1.0 }, { 16.0, -1.0 }, { 6.0, 1.0 }, { 10.0, -1.0 }, { 14.0, 1.0 }, { 18.0, -1.0 } };

    df.load_column<ary_col_t2>("ARY XYVAR 1", std::move(ary_xyvar_1), nan_policy::dont_pad_with_nans);
    df.load_column<ary_col_t2>("ARY XYVAR 2", std::move(ary_xyvar_2), nan_policy::dont_pad_with_nans);
    df.load_column<ary_col_t2>("ARY XYVAR 3", std::move(ary_xyvar_3), nan_policy::dont_pad_with_nans);
    df.load_column<ary_col_t2>("ARY XYVAR 4", std::move(ary_xyvar_4), nan_policy::dont_pad_with_nans);

    const auto  ary_xyvar_res = df.pca_by_eigen<ary_col_t2>({ "ARY XYVAR 1", "ARY XYVAR 2", "ARY XYVAR 3", "ARY XYVAR 4" }, { .pct_comp_to_keep = 0.95 });

    assert(ary_xyvar_res.cols() == 2);
    assert(ary_xyvar_res.rows() == 8);

    assert(std::abs(ary_xyvar_res(0, 0) - 0.0) < 0.00000000001);
    assert(std::abs(ary_xyvar_res(0, 1) - -7.07107) < 0.0001);
    assert(std::abs(ary_xyvar_res(3, 1) - -10.6066) < 0.0001);
    assert(std::abs(ary_xyvar_res(5, 0) - -12.3744) < 0.0001);
    assert(std::abs(ary_xyvar_res(7, 0) - -19.4454) < 0.0001);
    assert(std::abs(ary_xyvar_res(7, 1) - -12.3744) < 0.0001);

    // 3) Redundant component (rank-deficient). cOmponent 2 is always x + y
    //
    std::vector<ary_col_t3> ary_red_1  { { 1.0, 2.0, 3.0 }, { 2.0, 3.0, 5.0 }, { 3.0, 1.0, 4.0 }, { 4.0, 2.0, 6.0 }, { 5.0, 3.0, 8.0 }, { 6.0, 1.0, 7.0 } };
    std::vector<ary_col_t3> ary_red_2  { { 2.0, 1.0, 3.0 }, { 3.0, 2.0, 5.0 }, { 4.0, 3.0, 7.0 }, { 5.0, 4.0, 9.0 }, { 6.0, 2.0, 8.0 }, { 7.0, 3.0, 10.0 } };
    std::vector<ary_col_t3> ary_red_3  { { 3.0, 3.0, 6.0 }, { 4.0, 4.0, 8.0 }, { 5.0, 2.0, 7.0 }, { 6.0, 1.0, 7.0 }, { 7.0, 3.0, 10.0 }, { 8.0, 2.0, 10.0 } };
    std::vector<ary_col_t3> ary_red_4  { { 4.0, 2.0, 6.0 }, { 5.0, 3.0, 8.0 }, { 6.0, 1.0, 7.0 }, { 7.0, 2.0, 9.0 }, { 8.0, 4.0, 12.0 }, { 9.0, 1.0, 10.0 } };

    df.load_column<ary_col_t3>("ARY RED 1", std::move(ary_red_1), nan_policy::dont_pad_with_nans);
    df.load_column<ary_col_t3>("ARY RED 2", std::move(ary_red_2), nan_policy::dont_pad_with_nans);
    df.load_column<ary_col_t3>("ARY RED 3", std::move(ary_red_3), nan_policy::dont_pad_with_nans);
    df.load_column<ary_col_t3>("ARY RED 4", std::move(ary_red_4), nan_policy::dont_pad_with_nans);

    const auto  ary_red_res = df.pca_by_eigen<ary_col_t3>({ "ARY RED 1", "ARY RED 2", "ARY RED 3", "ARY RED 4" }, { .num_comp_to_keep = 4 });

    assert(ary_red_res.cols() == 4);
    assert(ary_red_res.rows() == 6);
    assert(std::abs(ary_red_res(0, 0) - 8.77192) < 0.00001);
    assert(std::abs(ary_red_res(0, 2) - 1.20476) < 0.00001);
    assert(std::abs(ary_red_res(3, 1) - 5.40356) < 0.00001);
    assert(std::abs(ary_red_res(3, 2) - -0.937747) < 0.000001);
    assert(std::abs(ary_red_res(5, 0) - 22.7284) < 0.0001);
    assert(std::abs(ary_red_res(5, 1) - 6.16921) < 0.00001);
    assert(std::abs(ary_red_res(5, 3) - 3.43023) < 0.00001);
}

// ----------------------------------------------------------------------------

static void test_compact_svd()  {

    std::cout << "\nTesting compact_svd( ) ..." << std::endl;

    StrDataFrame    df;

    try  {
        df.read("IBM.csv", io_format::csv2);
    }
    catch (const DataFrameError &ex)  {
        std::cout << ex.what() << std::endl;
        ::exit(-1);
    }

    const auto  [U, S, V] = df.compact_svd<double>({ "IBM_Close", "IBM_Open", "IBM_High", "IBM_Low" });

    assert(U.rows() == 5031);
    // Compact version has the same column # as the original matrix
    //
    assert(U.cols() == 4);
    assert(std::fabs(U(0, 0) - -0.0115747) < 0.000001);
    assert(std::fabs(U(2, 3) - -0.0110622) < 0.000001);
    assert(std::fabs(U(4040, 2) - -0.0147074) < 0.000001);
    assert(std::fabs(U(4994, 1) - 0.0194639) < 0.000001);
    assert(std::fabs(U(5030, 3) - -0.000878688) < 0.000001);

    // In compact version zero rows at the end are omitted
    //
    assert(S.rows() == 4);
    assert(S.cols() == 4);
    assert(std::fabs(S(0, 0) - 141.821) < 0.001);
    assert(std::fabs(S(1, 1) - 1.91734) < 0.00001);
    assert(std::fabs(S(2, 2) - 1.62214) < 0.00001);
    assert(std::fabs(S(3, 3) - 0.73194) < 0.00001);
    assert(S(0, 2) == 0.0);
    assert(S(1, 2) == 0.0);
    assert(S(3, 0) == 0.0);

    assert(V.rows() == 4);
    assert(V.cols() == 4);
    assert(std::fabs(V(0, 0) - 0.499988) < 0.000001);
    assert(std::fabs(V(0, 2) - 0.003710) < 0.000001);
    assert(std::fabs(V(2, 2) - 0.700869) < 0.000001);
    assert(std::fabs(V(3, 1) - -0.00079) < 0.000001);
    assert(std::fabs(V(3, 3) - 0.491216) < 0.000001);

    // Now multidimensional data
    //
    constexpr std::size_t   dim { 3 };

    using ary_col_t = std::array<double, dim>;
    using vec_col_t = std::vector<double>;

    std::vector<ary_col_t>  md_ary_col1  {
        {3,1,3}, {3,4,3}, {6,3,2}, {4,4,7}, {4,1,5}, {4,6,1}, {8,2,6}, {7,8,4}, {8,3,1}, {7,4,6}, {5,4,3}, {8,6,4}, {4,5,3}, {2,4,2}, {5,5,4}, {1,4,4},
        {2,3,1}, {8,1,7}, {1,1,2}, {5,8,1}, {1,4,5}, {7,8,2}, {8,8,3}, {7,3,6}, {5,6,5}, {4,4,1}, {2,1,2}, {1,6,5}, {1,6,7}, {6,4,4},
    };
    std::vector<ary_col_t>  md_ary_col2  {
        {1,7,3}, {5,7,3}, {4,5,3}, {6,4,2}, {3,7,3}, {4,3,1}, {8,8,6}, {4,8,7}, {3,6,7}, {1,6,8}, {2,5,4}, {2,2,2}, {6,2,8}, {3,2,8}, {7,5,4}, {4,2,1},
        {3,3,3}, {8,3,2}, {4,8,5}, {6,2,6}, {8,8,2}, {7,3,6}, {1,7,7}, {7,8,8}, {7,5,7}, {1,3,4}, {5,4,1}, {5,1,8}, {2,3,2}, {5,6,2},
    };
    std::vector<ary_col_t>  md_ary_col3  {
        {5,2,5}, {7,7,6}, {1,4,6}, {3,1,5}, {5,6,8}, {6,2,1}, {3,4,3}, {7,4,1}, {6,7,5}, {2,6,2}, {1,6,6}, {1,3,8}, {5,4,7}, {4,5,1}, {7,6,6}, {5,6,2},
        {2,3,5}, {5,8,2}, {8,1,4}, {6,3,7}, {1,1,8}, {7,2,7}, {6,3,2}, {1,6,7}, {4,1,2}, {3,4,5}, {7,1,6}, {5,8,5}, {8,8,7}, {8,8,8},
    };

    std::vector<vec_col_t>  md_vec_col1  {
        {3,1,3}, {3,4,3}, {6,3,2}, {4,4,7}, {4,1,5}, {4,6,1}, {8,2,6}, {7,8,4}, {8,3,1}, {7,4,6}, {5,4,3}, {8,6,4}, {4,5,3}, {2,4,2}, {5,5,4}, {1,4,4},
        {2,3,1}, {8,1,7}, {1,1,2}, {5,8,1}, {1,4,5}, {7,8,2}, {8,8,3}, {7,3,6}, {5,6,5}, {4,4,1}, {2,1,2}, {1,6,5}, {1,6,7}, {6,4,4},
    };
    std::vector<vec_col_t>  md_vec_col2  {
        {1,7,3}, {5,7,3}, {4,5,3}, {6,4,2}, {3,7,3}, {4,3,1}, {8,8,6}, {4,8,7}, {3,6,7}, {1,6,8}, {2,5,4}, {2,2,2}, {6,2,8}, {3,2,8}, {7,5,4}, {4,2,1},
        {3,3,3}, {8,3,2}, {4,8,5}, {6,2,6}, {8,8,2}, {7,3,6}, {1,7,7}, {7,8,8}, {7,5,7}, {1,3,4}, {5,4,1}, {5,1,8}, {2,3,2}, {5,6,2},
    };
    std::vector<vec_col_t>  md_vec_col3  {
        {5,2,5}, {7,7,6}, {1,4,6}, {3,1,5}, {5,6,8}, {6,2,1}, {3,4,3}, {7,4,1}, {6,7,5}, {2,6,2}, {1,6,6}, {1,3,8}, {5,4,7}, {4,5,1}, {7,6,6}, {5,6,2},
        {2,3,5}, {5,8,2}, {8,1,4}, {6,3,7}, {1,1,8}, {7,2,7}, {6,3,2}, {1,6,7}, {4,1,2}, {3,4,5}, {7,1,6}, {5,8,5}, {8,8,7}, {8,8,8},
    };

    df.load_column<ary_col_t>("ARY COL 1", std::move(md_ary_col1), nan_policy::dont_pad_with_nans);
    df.load_column<ary_col_t>("ARY COL 2", std::move(md_ary_col2), nan_policy::dont_pad_with_nans);
    df.load_column<ary_col_t>("ARY COL 3", std::move(md_ary_col3), nan_policy::dont_pad_with_nans);

    df.load_column<vec_col_t>("VEC COL 1", std::move(md_vec_col1), nan_policy::dont_pad_with_nans);
    df.load_column<vec_col_t>("VEC COL 2", std::move(md_vec_col2), nan_policy::dont_pad_with_nans);
    df.load_column<vec_col_t>("VEC COL 3", std::move(md_vec_col3), nan_policy::dont_pad_with_nans);

    const auto  [ary_U, ary_S, ary_V] = df.compact_svd<ary_col_t>({ "ARY COL 1", "ARY COL 2", "ARY COL 3" });
    const auto  [vec_U, vec_S, vec_V] = df.compact_svd<vec_col_t>({ "VEC COL 1", "VEC COL 2", "VEC COL 3" });

    assert(ary_U.rows() == 30);
    assert(ary_U.cols() == 9);
    assert(vec_U.rows() == 30);
    assert(vec_U.cols() == 9);
    assert(std::fabs(ary_U(0, 0) - 0.18325) < 0.00001);
    assert(std::fabs(ary_U(0, 5) - 0.063227) < 0.000001);
    assert(std::fabs(ary_U(5, 2) - 0.094879) < 0.000001);
    assert(std::fabs(ary_U(5, 6) - 0.338979) < 0.000001);
    assert(std::fabs(ary_U(29, 0) - 0.006552) < 0.000001);
    assert(std::fabs(ary_U(29, 4) - -0.134926) < 0.000001);
    assert(std::fabs(ary_U(29, 8) - 0.075926) < 0.000001);
    assert(std::fabs(vec_U(0, 0) - 0.18325) < 0.00001);
    assert(std::fabs(vec_U(0, 5) - 0.063227) < 0.000001);
    assert(std::fabs(vec_U(5, 2) - 0.094879) < 0.000001);
    assert(std::fabs(vec_U(5, 6) - 0.338979) < 0.000001);
    assert(std::fabs(vec_U(29, 0) - 0.006552) < 0.000001);
    assert(std::fabs(vec_U(29, 4) - -0.134926) < 0.000001);
    assert(std::fabs(vec_U(29, 8) - 0.075926) < 0.000001);

    assert(ary_S.rows() == 9);
    assert(ary_S.cols() == 9);
    assert(vec_S.rows() == 9);
    assert(vec_S.cols() == 9);
    assert(std::fabs(ary_S(0, 0) - 7.08301) < 0.00001);
    assert(std::fabs(ary_S(1, 1) - 6.97274) < 0.00001);
    assert(std::fabs(ary_S(4, 4) - 5.57452) < 0.00001);
    assert(std::fabs(ary_S(8, 8) - 2.63947) < 0.00001);
    assert(std::fabs(ary_S(0, 1) - 0.0) < 0.00000001);
    assert(std::fabs(ary_S(4, 5) - 0.0) < 0.00000001);
    assert(std::fabs(ary_S(8, 0) - 0.0) < 0.00000001);
    assert(std::fabs(ary_S(8, 5) - 0.0) < 0.00000001);
    assert(std::fabs(vec_S(0, 0) - 7.08301) < 0.00001);
    assert(std::fabs(vec_S(1, 1) - 6.97274) < 0.00001);
    assert(std::fabs(vec_S(4, 4) - 5.57452) < 0.00001);
    assert(std::fabs(vec_S(8, 8) - 2.63947) < 0.00001);
    assert(std::fabs(vec_S(0, 1) - 0.0) < 0.00000001);
    assert(std::fabs(vec_S(4, 5) - 0.0) < 0.00000001);
    assert(std::fabs(vec_S(8, 0) - 0.0) < 0.00000001);
    assert(std::fabs(vec_S(8, 5) - 0.0) < 0.00000001);

    assert(ary_V.rows() == 9);
    assert(ary_V.cols() == 9);
    assert(vec_V.rows() == 9);
    assert(vec_V.cols() == 9);
    assert(std::fabs(ary_V(0, 0) - -0.57947) < 0.00001);
    assert(std::fabs(ary_V(0, 5) - 0.32352) < 0.00001);
    assert(std::fabs(ary_V(6, 3) - -0.03933) < 0.00001);
    assert(std::fabs(ary_V(6, 8) - -0.34836) < 0.00001);
    assert(std::fabs(ary_V(8, 1) - -0.314893) < 0.000001);
    assert(std::fabs(ary_V(8, 6) - -0.200483) < 0.000001);
    assert(std::fabs(vec_V(0, 0) - -0.57947) < 0.00001);
    assert(std::fabs(vec_V(0, 5) - 0.32352) < 0.00001);
    assert(std::fabs(vec_V(6, 3) - -0.03933) < 0.00001);
    assert(std::fabs(vec_V(6, 8) - -0.34836) < 0.00001);
    assert(std::fabs(vec_V(8, 1) - -0.314893) < 0.000001);
    assert(std::fabs(vec_V(8, 6) - -0.200483) < 0.000001);
}

C++ DataFrame