Skip to content

Fix #2495 - Fix PPCA noise variance dilution and null-space leakage in MorphologicalDeviationScore #2496

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
akenmorris merged 2 commits into from
Feb 19, 2026
+59 −104
Merged

Fix #2495 - Fix PPCA noise variance dilution and null-space leakage in MorphologicalDeviationScore #2496

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
42af51f
Fix #2495 - Fix PPCA noise variance dilution and null-space leakage i…
akenmorris Feb 18, 2026
90142b6
Add back comment
akenmorris Feb 18, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
134 changes: 52 additions & 82 deletions Libs/Optimize/Function/EarlyStop/MorphologicalDeviationScore.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
#include "MorphologicalDeviationScore.h"

#include <Logging.h>
namespace shapeworks {

Expand All @@ -8,11 +9,9 @@ MorphologicalDeviationScore::MorphologicalDeviationScore() = default;
bool MorphologicalDeviationScore::SetControlShapes(const Eigen::MatrixXd& X) {
try {
bool ppca_status = FitPPCA(X);
precision_matrix_ = ComputePrecisionMatrix();
is_fitted_ = ppca_status;
} catch (std::exception& e) {
SW_ERROR("Exception in setting control shapes for early stopping {}",
e.what());
SW_ERROR("Exception in setting control shapes for early stopping {}", e.what());
return false;
}
return is_fitted_;
Expand All @@ -27,108 +26,79 @@ bool MorphologicalDeviationScore::FitPPCA(const Eigen::MatrixXd& X) {
Eigen::MatrixXd X_c = X.rowwise() - mean_; // (n x d)

try {
Eigen::JacobiSVD<Eigen::MatrixXd> svd(
X_c, Eigen::ComputeThinU | Eigen::ComputeThinV);
Eigen::JacobiSVD<Eigen::MatrixXd> svd(X_c, Eigen::ComputeThinU | Eigen::ComputeThinV);
Eigen::VectorXd s = svd.singularValues(); // (r,)
Eigen::MatrixXd V = svd.matrixV(); // (d x r)
Eigen::VectorXd eigvals = (s.array().square()) / (n - 1);

// save_matrix(X, std::string(DEBUG_FILES_PTH) + "/X.npy");
// save_matrix(X_c, std::string(DEBUG_FILES_PTH) + "/X_c.npy");
// save_matrix(V, std::string(DEBUG_FILES_PTH) + "/V.npy");
// save_vector(eigvals, std::string(DEBUG_FILES_PTH) + "/eigvals.npy");

Eigen::VectorXd cumsum = eigvals;
for (int i = 1; i < eigvals.size(); ++i) cumsum[i] += cumsum[i - 1];
Eigen::VectorXd cumvar = cumsum / eigvals.sum();
n_components_ = 0;
while (n_components_ < cumvar.size() && cumvar[n_components_] <= retained_variance_ratio_)
++n_components_;
int q = n_components_;
principal_components_variance_ = eigvals.head(q); // (q,)
components_ = V.leftCols(q); // (d x q)
// save_matrix(components_, std::string(DEBUG_FILES_PTH) + "/components.npy");
// save_vector(principal_components_variance_,
// std::string(DEBUG_FILES_PTH) + "/eigvals.npy");
noise_variance_ = eigvals.tail(eigvals.size() - q).sum() / double(d - q);
return true;
int q = 0;
while (q < cumvar.size() && cumvar[q] <= retained_variance_ratio_) ++q;

} catch (std::exception& e) {
SW_ERROR(
"SVD computation failed for MorphologicalDeviationScore: {}",
e.what());
return false;
}
}
if (q == 0) {
SW_ERROR("PPCA: no components retained at {}% variance", retained_variance_ratio_ * 100);
return false;
}

//---------------------------------------------------------------------------
Eigen::MatrixXd MorphologicalDeviationScore::ComputeCovarianceMatrix() {
const int d = components_.rows();
const int q = components_.cols();
// Determine effective rank: trim structurally zero eigenvalues
const double eps = 1e-10;
int rank = eigvals.size();
while (rank > q && eigvals[rank - 1] < eps) {
--rank;
}

if (q == 0) {
return noise_variance_ * Eigen::MatrixXd::Identity(d, d);
}
if (rank <= q) {
SW_ERROR("PPCA: effective rank ({}) <= retained components ({}), cannot estimate noise variance", rank, q);
return false;
}

Eigen::MatrixXd diag_lambda = principal_components_variance_.asDiagonal(); // (q x q)
Eigen::MatrixXd cov =
components_ *
(diag_lambda - noise_variance_ * Eigen::MatrixXd::Identity(q, q)) *
components_.transpose() +
noise_variance_ * Eigen::MatrixXd::Identity(d, d); // (d x d)
return cov;
}
// Noise variance: average of excluded-but-real eigenvalues
// Uses (rank - q) denominator instead of (d - q) to avoid dilution
// by structurally zero dimensions beyond the data rank
noise_variance_ = eigvals.segment(q, rank - q).sum() / double(rank - q);

//---------------------------------------------------------------------------
Eigen::MatrixXd MorphologicalDeviationScore::ComputePrecisionMatrix() {
try {
const int d = components_.rows();
const int q = components_.cols();
// Store full rank-dimensional basis for projection-based scoring
all_components_ = V.leftCols(rank); // (d x rank)

Eigen::MatrixXd A_inv = (1 / noise_variance_) * Eigen::MatrixXd::Identity(d, d);
// Build per-dimension precision weights
precision_weights_.resize(rank);
for (int i = 0; i < q; ++i) {
precision_weights_[i] = 1.0 / eigvals[i];
}
for (int i = q; i < rank; ++i) {
precision_weights_[i] = 1.0 / noise_variance_;
}

Eigen::MatrixXd U = components_;
Eigen::MatrixXd C_inv = principal_components_variance_.cwiseInverse().asDiagonal();
Eigen::MatrixXd V = components_.transpose();
// Woodbury matrix identity
// (A + U C V)^(-1)
// = A^(-1)
// - A^(-1) * U * (C^(-1) + V * A^(-1) * U)^(-1) * V * A^(-1)
SW_LOG("PPCA fit: n={}, d={}, rank={}, q={}, noise_var={:.6f}", n, d, rank, q, noise_variance_);

Eigen::MatrixXd precision = A_inv - A_inv * U * (C_inv + V * A_inv * U).inverse() * V * A_inv;
// save_matrix(precision, std::string(DEBUG_FILES_PTH)+"/precision.npy");
return true;

return precision;
} catch (std::exception& e) {
SW_ERROR(
"Failed to compute precision matrix in the early stopping score function: {}",
e.what());
return Eigen::MatrixXd();
SW_ERROR("SVD computation failed for MorphologicalDeviationScore: {}", e.what());
return false;
}
}

//---------------------------------------------------------------------------
Eigen::VectorXd MorphologicalDeviationScore::GetMorphoDevScore(
const Eigen::MatrixXd& X) {
Eigen::VectorXd MorphologicalDeviationScore::GetMorphoDevScore(const Eigen::MatrixXd& X) {
try {

if (!is_fitted_) {
throw std::runtime_error(
"PPCA model is not fitted on control shapes.");
throw std::runtime_error("PPCA model is not fitted on control shapes.");
}
const int n = X.rows();
// Eigen::MatrixXd X_c = X.rowwise() - mean_;
// Eigen::VectorXd dist(n);
// for (int i = 0; i < n; ++i) {
// Eigen::RowVectorXd xi = X_c.row(i);
// dist(i) = std::sqrt(
// std::max(0.0, (xi * precision_matrix_ * xi.transpose())(0)));
// }

Eigen::MatrixXd X_bar = X.rowwise() - mean_;
Eigen::MatrixXd Y = X_bar * precision_matrix_; // (n_samples × n_features)
Eigen::VectorXd sq_mahal = (X_bar.array() * Y.array()).rowwise().sum(); // quadratic forms
Eigen::VectorXd mahalanobis = sq_mahal.array().sqrt(); // final distances

// Project into the rank-dimensional subspace, eliminating
// structurally zero dimensions that cannot be estimated from data
Eigen::MatrixXd X_bar = X.rowwise() - mean_; // (n x d)
Eigen::MatrixXd Z = X_bar * all_components_; // (n x rank)

// Apply per-dimension precision weights: 1/lambda_i for retained PCs,
// 1/noise_variance for excluded-but-real PCs
Eigen::MatrixXd Z_weighted = Z.array().rowwise() * precision_weights_.transpose().array(); // (n x rank)

Eigen::VectorXd sq_mahal = (Z.array() * Z_weighted.array()).rowwise().sum(); // (n,)
Eigen::VectorXd mahalanobis = sq_mahal.array().max(0.0).sqrt();
return mahalanobis;

} catch (std::exception& e) {
Expand All @@ -140,4 +110,4 @@ Eigen::VectorXd MorphologicalDeviationScore::GetMorphoDevScore(
}
}

} // namespace shapeworks
} // namespace shapeworks
29 changes: 7 additions & 22 deletions Libs/Optimize/Function/EarlyStop/MorphologicalDeviationScore.h
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,6 @@
#pragma once

#include <Eigen/Dense>
// #include "cnpy.h"

namespace shapeworks {
class MorphologicalDeviationScore {
Expand All @@ -18,28 +17,14 @@ class MorphologicalDeviationScore {
bool is_fitted_ = false;

// Fitted model parameters
Eigen::RowVectorXd mean_; // (1 x d)
Eigen::MatrixXd components_; // (d x q)
Eigen::VectorXd principal_components_variance_; // (q,)
int n_components_ = 0;
Eigen::RowVectorXd mean_; // (1 x d)
double noise_variance_ = 0.0;
double retained_variance_ratio_ = 0.95;
// Derived matrices
Eigen::MatrixXd precision_matrix_; // (d x d)
// Helper functions

// Full-rank basis and per-dimension precision weights for scoring
Eigen::MatrixXd all_components_; // (d x rank) — all non-zero eigenvectors
Eigen::VectorXd precision_weights_; // (rank,) — 1/lambda_i or 1/noise_variance

bool FitPPCA(const Eigen::MatrixXd& X);
Eigen::MatrixXd ComputeCovarianceMatrix();
Eigen::MatrixXd ComputePrecisionMatrix();
// inline static void save_vector(const Eigen::VectorXd& v,
// const std::string& fname) {
// cnpy::npy_save(fname, v.data(), {(size_t)v.size()}, "w");
// };
// inline static void save_matrix(const Eigen::MatrixXd m,
// const std::string& fname) {
// Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
// X_rm = m;
// cnpy::npy_save(fname, X_rm.data(), {(size_t)m.rows(), (size_t)m.cols()},
// "w");
// };
};
} // namespace shapeworks
} // namespace shapeworks
Loading