Skip to content

Commit

Permalink
Turn on unique sign algorithm by default
Browse files Browse the repository at this point in the history
  • Loading branch information
Martin D. Weinberg committed Feb 14, 2025
1 parent 8661d07 commit f19fcd4
Show file tree
Hide file tree
Showing 2 changed files with 184 additions and 30 deletions.
169 changes: 169 additions & 0 deletions expui/SvdSignChoice.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<K; k++) {
// Remove all but target dimension for numerical stability
S1(k) = 0.0;
auto Y = X - U * S1.asDiagonal() * V.transpose();

// Restore the value
S1(k) = S(k);

// d_j = U^T_k * Y_j
Eigen::VectorXd dL = Y.transpose() * U.col(k);

// sum of sgn(dL_j)*dL_j^2
sL(k) += dL.dot(dL.cwiseAbs());

// d_i = V^T_k * (Y^T)_i
Eigen::VectorXd dR = Y * V.col(k);

// sum of sgn(dR_i)*dR_i^2
sR(k) += dR.dot(dR.cwiseAbs());
}

auto sgn = [](double val) -> int
{
return (0.0 < val) - (val < 0.0);
};

// 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));
}

// Done
}

void SvdSignChoice
(const Eigen::MatrixXd& X,
const Eigen::MatrixXd& U, const Eigen::VectorXd& S, 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<K; k++) {
// Remove all but target dimension for numerical stability
S1(k) = 0.0;
auto Y = X - U * S1.asDiagonal() * V.transpose();

// Restore the value
S1(k) = S(k);

// d_j = U^T_k * Y_j
Eigen::VectorXd dL = Y.transpose() * U.col(k);

// sum of sgn(dL_j)*dL_j^2
sL(k) += dL.dot(dL.cwiseAbs());

// d_i = V^T_k * (Y^T)_i
Eigen::VectorXd dR = Y * V.col(k);

// sum of sgn(dR_i)*dR_i^2
sR(k) += dR.dot(dR.cwiseAbs());
}

auto sgn = [](double val) -> int
{
return (0.0 < val) - (val < 0.0);
};

// 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
V.col(k) *= sgn(sR(k));
}

// Done
}

}
45 changes: 15 additions & 30 deletions expui/expMSSA.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -55,9 +55,15 @@

namespace MSSA {

// Update sign of left vectors U only
void SvdSignChoice
(const Eigen::MatrixXd& X,
Eigen::MatrixXd& U, const Eigen::VectorXd& S, Eigen::MatrixXd& V);
Eigen::MatrixXd& U, const Eigen::VectorXd& S, const Eigen::MatrixXd& V);

// Update sign of right vectors V only
void SvdSignChoice
(const Eigen::MatrixXd& X,
const Eigen::MatrixXd& U, const Eigen::VectorXd& S, Eigen::MatrixXd& V);

Eigen::MatrixXd expMSSA::wCorrKey(const Key& key)
{
Expand Down Expand Up @@ -299,20 +305,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::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);
}
if (useSignChoice) SvdSignChoice(cov, U, S, svd.matrixV());
}
} else if (params["BDCSVD"]) {
// -->Using BDC
Expand All @@ -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<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);
}
if (useSignChoice) SvdSignChoice(cov, U, S, svd.matrixV());
}
} else {
// -->Use Random approximation algorithm from Halko, Martinsson,
Expand All @@ -345,28 +339,19 @@ namespace MSSA {
RedSVD::RedSVD<Eigen::MatrixXd> 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::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);
}
if (useSignChoice) SvdSignChoice(cov, U, S, U.transpose());
} 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);
}
if (useSignChoice) SvdSignChoice(cov, U, S, U.transpose());
}
}
}
Expand Down Expand Up @@ -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
//
Expand Down

0 comments on commit f19fcd4

Please sign in to comment.