From 9b7e0fac6fb88a2154697a3afd2bf11c101af5f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Fran=C3=A7ois-David=20Collin?= Date: Thu, 4 Jun 2020 17:15:03 +0200 Subject: [PATCH] Fix lda with constant within class variables --- src/lda-eigen.hpp | 22 ++++++++++------- src/matutils.hpp | 25 ++++++++++++++++--- test/lda-eigen-test.cpp | 53 +++++++++++++++++++++++++++++++++++++++++ 3 files changed, 89 insertions(+), 11 deletions(-) diff --git a/src/lda-eigen.hpp b/src/lda-eigen.hpp index 2557d3b..8a05bf4 100644 --- a/src/lda-eigen.hpp +++ b/src/lda-eigen.hpp @@ -25,13 +25,13 @@ typedef Matrix VectorXs; * @param Ld computes matrix whoses columns are Ld vectors */ template -void lda(const MatrixBase &x, +std::vector lda(const MatrixBase &x, const VectorXs &y, MatrixXd& Ld) { auto n = x.rows(); - auto p = x.cols(); auto K = y.maxCoeff() + 1; + auto p = x.cols(); // M = Centroids MatrixXd M = MatrixXd::Zero(K, p); @@ -59,7 +59,7 @@ void lda(const MatrixBase &x, MatrixXd D(n, p); // W = Within-class covariance matrix - MatrixXd W = MatrixXd::Zero(p, p); + MatrixXd Wraw = MatrixXd::Zero(p, p); size_t slicepos = 0; for (auto c = 0; c < K; c++) { @@ -67,10 +67,15 @@ void lda(const MatrixBase &x, MatrixXd Dc = x.block(slicepos, 0, d[c], p); // Centering Dc.rowwise() -= M.row(c); - W += Dc.transpose() * Dc; + Wraw += Dc.transpose() * Dc; slicepos += d[c]; } - W /= static_cast(n - K); + Wraw /= static_cast(n - K); + auto validvars = std::vector(); + for(auto i = 0; i< x.cols(); i++) { + if (std::abs(Wraw(i,i)) >= 1.0e-8) validvars.push_back(i); + } + auto W = Wraw(validvars,validvars); // [4,4]((0.265008,0.0927211,0.167514,0.0384014),(0.0927211,0.115388,0.0552435,0.0327102),(0.167514,0.0552435,0.185188,0.0426653),(0.0384014,0.0327102,0.0426653,0.0418816)) // Calculate pseudo-inverse square root for W with SVD @@ -84,9 +89,9 @@ void lda(const MatrixBase &x, // } //auto W12 = svd.matrixU() * svd.eigenvalues().array().inverse().sqrt().matrix().asDiagonal() * svd.matrixU().transpose(); // auto W12 = svd.operatorInverseSqrt(); - - MatrixXd Mstar = M * W12; - VectorXd mstar = VectorXd::Zero(p); + MatrixXd Mvalid = M(all,validvars); + MatrixXd Mstar = Mvalid * W12; + VectorXd mstar = VectorXd::Zero(validvars.size()); for (auto c = 0; c < K; c++) mstar += Mstar.row(c); mstar /= static_cast(K); @@ -98,4 +103,5 @@ void lda(const MatrixBase &x, JacobiSVD svd2(Bstar,ComputeThinU); MatrixXd Vl = W12 * svd2.matrixU(); Ld = Vl.leftCols(K-1); + return validvars; } \ No newline at end of file diff --git a/src/matutils.hpp b/src/matutils.hpp index b971075..a16c4e9 100644 --- a/src/matutils.hpp +++ b/src/matutils.hpp @@ -75,6 +75,13 @@ void addLinearComb(Eigen::Ref +void addLinearComb(Eigen::Ref> ref, MatrixBase const & X_,const MatrixBase& M, std::vector& validvars) { + auto ncols = X_.cols(); + auto& X = constCastAddColsMatrix(X_,M.cols()); + X.block(0,ncols,X.rows(),M.cols()) = ref(all,validvars).eval() * M; +} + template void addLinearComb(MatrixBase const & X_, const MatrixBase& M) { auto ncols = X_.cols(); @@ -82,14 +89,26 @@ void addLinearComb(MatrixBase const & X_, const MatrixBase +void addLinearComb(MatrixBase const & X_, const MatrixBase& M, std::vector& validvars) { + auto ncols = X_.cols(); + auto& X = constCastAddColsMatrix(X_,M.cols()); + X.block(0,ncols,X.rows(),M.cols()) = X.block(0,0,X.rows(),ncols)(all,validvars).eval() * M; +} + + template void addLda(Reftable& rf, MatrixXd& data, MatrixBase const &statobs) { VectorXs scen(rf.nrec); for(auto i = 0; i < rf.nrec; i++) scen(i) = static_cast(rf.scenarios[i]) - 1; MatrixXd Ld; - lda(rf.stats, scen, Ld); - addLinearComb(rf.stats,data,Ld); - addLinearComb(statobs,Ld); + std::vector validvars = lda(rf.stats, scen, Ld); + addLinearComb(rf.stats,data,Ld,validvars); + addLinearComb(statobs,Ld,validvars); + for(auto i = 0; i < rf.stats.cols(); i++) { + if (std::find(std::begin(validvars),std::end(validvars),i) == std::end(validvars)) + std::cout << "LDA Warning : " << rf.stats_names[i] << " is constant within class, removed." << std::endl; + } for(auto i = 0; i < rf.nrecscen.size() - 1; i++) { rf.stats_names.push_back("LDA" + std::to_string(i+1)); } diff --git a/test/lda-eigen-test.cpp b/test/lda-eigen-test.cpp index aaa5cd2..1b34c2b 100644 --- a/test/lda-eigen-test.cpp +++ b/test/lda-eigen-test.cpp @@ -213,3 +213,56 @@ TEST_CASE("LDA with Eigen") CHECK(std::min(LdPlus.col(c).lpNorm(),LdMinus.col(c).lpNorm()) == Approx(0.0).margin(1e-13)); } } + +TEST_CASE("LDA with Eigen and constant variable") +{ + const size_t K = 3; + const size_t p = 4; + const size_t n = (S.size() + C.size() + V.size()) / p; + Matrix y(n); + MatrixXd x(n,p+1); + // S + for (auto i = 0; i < (S.size()) / p; i++) + { + y[i] = 0; + for (auto j = 0; j < p; j++) + x(i, j) = S[i * p + j]; + } + for (auto i = 0; i < (C.size()) / p; i++) + { + auto ii = (S.size() / p) + i; + y[ii] = 1; + for (auto j = 0; j < p; j++) + x(ii, j) = C[i * p + j]; + } + for (auto i = 0; i < (V.size()) / p; i++) + { + auto ii = (S.size() + C.size()) / p + i; + y[ii] = 2; + for (auto j = 0; j < p; j++) + x(ii, j) = V[i * p + j]; + } + for(auto i = 0; i < n; i++) { + x(i,p) = x(i,2); + x(i,2) = y[i]; + } + MatrixXd Ld; + lda(x,y,Ld); + + MatrixXd LdMass(p,K-1); + LdMass << + // 0.8293776, 0.02410215, + // 1.5344731, 2.16452123, + // -2.2012117, -0.93192121, + // -2.8104603, 2.83918785; + -0.82937764226600674, 0.024102148876954166, + -1.53447306770001091, 2.164521234658435489, + 2.81046030884310172, 2.839187852982734128, + 2.20121165556177356, -0.931921210029371894; + + auto LdPlus = Ld + LdMass; + auto LdMinus = Ld - LdMass; + for(auto c = 0; c < K - 1; c++) { + CHECK(std::min(LdPlus.col(c).lpNorm(),LdMinus.col(c).lpNorm()) == Approx(0.0).margin(1e-13)); + } +}