diff --git a/expui/SvdSignChoice.cc b/expui/SvdSignChoice.cc index bb511e71..14a7b63e 100644 --- a/expui/SvdSignChoice.cc +++ b/expui/SvdSignChoice.cc @@ -109,4 +109,173 @@ namespace MSSA { // Done } + + void SvdSignChoice + (const Eigen::MatrixXd& X, + Eigen::MatrixXd& U, const Eigen::VectorXd& S, const Eigen::MatrixXd& V) + { + // SDV dimensions + int I = U.rows(); + int J = V.rows(); + int K = S.size(); + + // Dimensions + // ---------- + // X is a [I x J] data matrix + // U is a [I x K] matrix + // S is a [K x 1] vector (diagonal matrix) + // V is a [J x K] matrix + + // Sanity checks + if (U.cols() != K) + throw std::invalid_argument("SvdSignChoice: U has wrong dimensions"); + + if (V.cols() != K) + throw std::invalid_argument("SvdSignChoice: V has wrong dimensions"); + + if (X.rows() != I || X.cols() != J) + throw std::invalid_argument("SvdSignChoice: X dimensions do not match SVD input"); + + // Sign determination loop + // + Eigen::VectorXd sL(K), sR(K); + sL.setZero(); + sR.setZero(); + + // Working, non-const instance + auto S1 = S; + + // Get projections from left and right singular vectors onto data + // matrix + // + for (int k=0; k int + { + return (0.0 < val) - (val < 0.0); + }; + + // Determine and apply the sign correction + // + for (int k=0; k int + { + return (0.0 < val) - (val < 0.0); + }; + + // Determine and apply the sign correction + // + for (int k=0; k svd(cov, Eigen::ComputeThinU | Eigen::ComputeThinV); S = svd.singularValues(); U = svd.matrixU(); - if (useSignChoice) { - auto V = svd.matrixV(); - SvdSignChoice(cov, U, S, V); - } + if (useSignChoice) SvdSignChoice(cov, U, S, svd.matrixV()); } } else if (params["BDCSVD"]) { // -->Using BDC @@ -322,20 +322,14 @@ namespace MSSA { svd(YY, Eigen::ComputeThinU | Eigen::ComputeThinV); S = svd.singularValues(); U = svd.matrixV(); - if (useSignChoice) { - auto V = svd.matrixU(); - SvdSignChoice(YY, V, S, U); - } + if (useSignChoice) SvdSignChoice(YY, svd.matrixU(), S, U); } else { // Covariance matrix Eigen::BDCSVD svd(cov, Eigen::ComputeThinU | Eigen::ComputeThinV); S = svd.singularValues(); U = svd.matrixU(); - if (useSignChoice) { - auto V = svd.matrixV(); - SvdSignChoice(cov, U, S, V); - } + if (useSignChoice) SvdSignChoice(cov, U, S, svd.matrixV()); } } else { // -->Use Random approximation algorithm from Halko, Martinsson, @@ -345,28 +339,19 @@ namespace MSSA { RedSVD::RedSVD svd(YY, srank); S = svd.singularValues(); U = svd.matrixV(); - if (useSignChoice) { - auto V = svd.matrixU(); - SvdSignChoice(YY, V, S, U); - } + if (useSignChoice) SvdSignChoice(YY, svd.matrixU(), S, U); } else { // Covariance matrix if (params["RedSym"]) { RedSVD::RedSymEigen eigen(cov, srank); S = eigen.eigenvalues().reverse(); U = eigen.eigenvectors().rowwise().reverse(); - if (useSignChoice) { - Eigen::MatrixXd V = U.transpose(); - SvdSignChoice(cov, U, S, V); - } + if (useSignChoice) SvdSignChoice(cov, U, S, U.transpose()); } else { RedSVD::RedSVD svd(cov, srank); S = svd.singularValues(); U = svd.matrixU(); - if (useSignChoice) { - Eigen::MatrixXd V = U.transpose(); - SvdSignChoice(cov, U, S, V); - } + if (useSignChoice) SvdSignChoice(cov, U, S, U.transpose()); } } } @@ -1840,7 +1825,7 @@ namespace MSSA { } - expMSSA::expMSSA(const mssaConfig& config, int nW, int nPC, const std::string flags) : numW(nW), npc(nPC), trajectory(true), useSignChoice(false) + expMSSA::expMSSA(const mssaConfig& config, int nW, int nPC, const std::string flags) : numW(nW), npc(nPC), trajectory(true), useSignChoice(true) { // Parse the YAML string //