Disambiguate the sign (phase) of the mSSA PCs

expui/CMakeLists.txt

@@ -38,7 +38,7 @@ endif()
 set(expui_SOURCES BasisFactory.cc BiorthBasis.cc FieldBasis.cc
   CoefContainer.cc CoefStruct.cc FieldGenerator.cc expMSSA.cc
   Coefficients.cc KMeans.cc Centering.cc ParticleIterator.cc
-  Koopman.cc BiorthBess.cc)
+  Koopman.cc BiorthBess.cc SvdSignChoice.cc)
 add_library(expui ${expui_SOURCES})
 set_target_properties(expui PROPERTIES OUTPUT_NAME expui)
 target_include_directories(expui PUBLIC ${common_INCLUDE})

expui/SvdSignChoice.cc

@@ -0,0 +1,110 @@
+#include <Eigen/Dense>
+
+namespace MSSA {
+
+  // SVD with corrected signs
+  //
+  // Input
+  // -----
+  // X is data matrix (constant)
+  // U is the matrix of left singular vectors
+  // S is the array of singular values (constant)
+  // V is the matrix of right singular vectors
+  //
+  // Output
+  // ------
+  // U, V returned corrected, disambiguated signs
+  //
+  // Reference
+  // ---------
+  // Bro, R., Acar, E., & Kolda, T. G. (2008). Resolving the sign
+  // ambiguity in the singular value decomposition.  Journal of
+  // Chemometrics: A Journal of the Chemometrics Society, 22(2),
+  // 135-140.
+  //
+  // URL:
+  // https://prod-ng.sandia.gov/techlib-noauth/access-control.cgi/2007/076422.pdf
+  //
+  void SvdSignChoice
+  (const Eigen::MatrixXd& X,
+   Eigen::MatrixXd& U, const Eigen::VectorXd& S, Eigen::MatrixXd& V)
+  {
+    // SDV dimensions
+    int I = U.rows();
+    int J = V.cols();
+    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.rows() != 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();
+
+    auto S1 = S;
+
+    auto sgn = [](double val) -> int
+    {
+      return (0.0 < val) - (val < 0.0);
+    };
+
+    // Get projects from left and right singular vectors onto data matrix
+    //
+    for (int k=0; k<K; k++) {
+
+      // Remove all but target dimension for numerical stability
+      S1(k) = 0.0;
+      auto Y = X - U * S1.asDiagonal() * V.transpose();
+      S1(k) = S(k);
+
+      // d_j = U^T_k * Y_j
+      for (int j=0; j<Y.cols(); j++) {
+	double d = U.col(k).dot(Y.col(j));
+	sL(k) += sgn(d) * d*d;
+      }
+
+      // d_i = V^T_k * (Y^T)_i
+      for (int i=0; i<Y.rows(); i++) {
+	double d = V.col(k).dot(Y.row(i));
+	sR(k) += sgn(d) * d*d;
+      }
+    }
+
+
+    // Determine and apply the sign correction
+    //
+    for (int k=0; k<K; k++) {
+      // If signs are opposite, flip the one with the smaller absolute
+      // value
+      //
+      if (sL(k)*sR(k) < 0.0) {
+	if (std::abs(sL(k)) < std::abs(sR(k)))
+	  sL(k) = -sL(k);
+	else
+	  sR(k) = -sR(k);
+      }
+
+      // Apply
+      U.col(k) *= sgn(sL(k));
+      V.col(k) *= sgn(sR(k));
+    }
+
+    // Done
+  }
+}

expui/expMSSA.cc

@@ -55,6 +55,12 @@
 
 namespace MSSA {
 
+  const bool useSignChoice = true;
+
+  void SvdSignChoice
+  (const Eigen::MatrixXd& X,
+   Eigen::MatrixXd& U, const Eigen::VectorXd& S, Eigen::MatrixXd& V);
+
   Eigen::MatrixXd expMSSA::wCorrKey(const Key& key)
   {
     if (RC.find(key)==RC.end()) {
@@ -289,52 +295,80 @@ namespace MSSA {
     //
     if (params["Jacobi"]) {
       // -->Using Jacobi
-      if (trajectory) {	// Trajectory matrix
+      if (trajectory) {		// Trajectory matrix
 	auto YY = Y/Scale;
 	Eigen::JacobiSVD<Eigen::MatrixXd>
 	  svd(YY, Eigen::ComputeThinU | Eigen::ComputeThinV);
 	S = svd.singularValues();
 	U = svd.matrixV();
+	if (useSignChoice) {
+	  auto V = svd.matrixU();
+	  SvdSignChoice(YY, V, S, U);
+	}
       }
       else {			// Covariance matrix
 	Eigen::JacobiSVD<Eigen::MatrixXd>
 	  svd(cov, Eigen::ComputeThinU | Eigen::ComputeThinV);
 	S = svd.singularValues();
 	U = svd.matrixU();
+	if (useSignChoice) {
+	  auto V = svd.matrixV();
+	  SvdSignChoice(cov, U, S, V);
+	}
       }
     } else if (params["BDCSVD"]) {
       // -->Using BDC
-      if (trajectory) {	// Trajectory matrix
+      if (trajectory) {		// Trajectory matrix
 	auto YY = Y/Scale;
 	Eigen::BDCSVD<Eigen::MatrixXd>
 	  svd(YY, Eigen::ComputeThinU | Eigen::ComputeThinV);
 	S = svd.singularValues();
 	U = svd.matrixV();
+	if (useSignChoice) {
+	  auto V = svd.matrixU();
+	  SvdSignChoice(YY, V, S, U);
+	}
       }
       else {			// Covariance matrix
 	Eigen::BDCSVD<Eigen::MatrixXd>
 	  svd(cov, Eigen::ComputeThinU | Eigen::ComputeThinV);
 	S = svd.singularValues();
 	U = svd.matrixU();
+	if (useSignChoice) {
+	  auto V = svd.matrixV();
+	  SvdSignChoice(cov, U, S, V);
+	}
       }
     } else {
       // -->Use Random approximation algorithm from Halko, Martinsson,
       //    and Tropp
-      if (trajectory) {	// Trajectory matrix
+      if (trajectory) {		// Trajectory matrix
 	auto YY = Y/Scale;
 	RedSVD::RedSVD<Eigen::MatrixXd> svd(YY, srank);
 	S = svd.singularValues();
 	U = svd.matrixV();
+	if (useSignChoice) {
+	  auto V = svd.matrixU();
+	  SvdSignChoice(YY, V, S, U);
+	}
       }
       else {			// Covariance matrix
 	if (params["RedSym"]) {
 	  RedSVD::RedSymEigen<Eigen::MatrixXd> eigen(cov, srank);
 	  S = eigen.eigenvalues().reverse();
 	  U = eigen.eigenvectors().rowwise().reverse();
+	  if (useSignChoice) {
+	    Eigen::MatrixXd V = U.transpose();
+	    SvdSignChoice(cov, U, S, V);
+	  }
 	} else {
 	  RedSVD::RedSVD<Eigen::MatrixXd> svd(cov, srank);
 	  S = svd.singularValues();
 	  U = svd.matrixU();
+	  if (useSignChoice) {
+	    Eigen::MatrixXd V = U.transpose();
+	    SvdSignChoice(cov, U, S, V);
+	  }
 	}
       }
     }