Replace VNL eigensystem with Eigen SVD for PCA computation

Switch ParticleShapeStatistics::compute_modes() from VNL's
vnl_symmetric_eigensystem to Eigen's BDCSVD algorithm. This provides:

- Better numerical stability by computing SVD directly on the centered
  data matrix rather than forming the covariance matrix A^T*A first
- Reduced dependency on legacy VNL library in favor of Eigen
- Cleaner code using Eigen's built-in vector normalization

The SVD singular values squared equal the eigenvalues, and the left
singular vectors (U) correspond to the eigenvectors in particle space.
Results are reversed to maintain ascending eigenvalue order for backward
compatibility.

Also add profiling instrumentation to compute_modes() and Studio main.

Author: akenmorris

Libs/Common/Profiling.cpp

@@ -157,6 +157,8 @@ void Profiler::finalize() {
     write_trace_file();
   }
 }
+
+//---------------------------------------------------------------------------
 void Profiler::write_profile_report() {
   // Calculate total time
   double total_time_ms = elapsed_timer_.elapsed();

Libs/Particles/ParticleShapeStatistics.cpp

@@ -7,6 +7,7 @@
 #include <vnl/algo/vnl_symmetric_eigensystem.h>
 
 #include "ExternalLibs/tinyxml/tinyxml.h"
+#include "Profiling.h"
 #include "ShapeEvaluation.h"
 
 namespace shapeworks {
@@ -577,49 +578,48 @@ int ParticleShapeStatistics::do_pca(std::shared_ptr<Project> project) {
 }
 
 //---------------------------------------------------------------------------
+
 int ParticleShapeStatistics::compute_modes() {
+  TIME_SCOPE("ParticleShapeStatistics::compute_modes");
   SW_DEBUG("computing PCA modes");
-  Eigen::MatrixXd A = points_minus_mean_.transpose() * points_minus_mean_ * (1.0 / ((double)(num_samples_ - 1)));
 
-  auto vnlA = vnl_matrix<double>(A.data(), A.rows(), A.cols());
-  vnl_symmetric_eigensystem<double> symEigen(vnlA);
+  double scale = std::sqrt(num_samples_ - 1);
+  Eigen::BDCSVD<Eigen::MatrixXd> svd(points_minus_mean_ / scale, Eigen::ComputeThinU);
 
-  Eigen::MatrixXd eigenSymEigenV =
-      Eigen::Map<Eigen::MatrixXd>(symEigen.V.transpose().data_block(), symEigen.V.rows(), symEigen.V.cols());
-  Eigen::VectorXd eigenSymEigenD = Eigen::Map<Eigen::VectorXd>(symEigen.D.data_block(), symEigen.D.rows(), 1);
+  Eigen::VectorXd eigenvalues_eigen = svd.singularValues().array().square();
 
-  eigenvectors_ = points_minus_mean_ * eigenSymEigenV;
-  eigenvalues_.resize(num_samples_);
+  // U * S gives us the eigenvectors in particle space (like points_minus_mean_ * V did before)
+  // But we want them normalized, so just use U (which is already orthonormal)
+  eigenvectors_ = svd.matrixU() * scale;  // Scale back to match original magnitudes
 
-  // normalize the eigenvectors
-  for (unsigned int i = 0; i < num_samples_; i++) {
-    double total = 0.0f;
-    for (unsigned int j = 0; j < num_dimensions_; j++) {
-      total += eigenvectors_(j, i) * eigenvectors_(j, i);
-    }
-    total = sqrt(total);
+  int num_modes = eigenvalues_eigen.size();
 
-    for (unsigned int j = 0; j < num_dimensions_; j++) {
-      eigenvectors_(j, i) = eigenvectors_(j, i) / (total + 1.0e-15);
-    }
+  // normalize the eigenvectors
+  for (int i = 0; i < eigenvectors_.cols(); i++) {
+    eigenvectors_.col(i).normalize();
+  }
 
-    eigenvalues_[i] = eigenSymEigenD(i);
+  // SVD returns values in descending order, but we need ascending order for backward compatibility
+  // Reverse both eigenvalues and eigenvector columns
+  eigenvalues_.resize(num_modes);
+  for (int i = 0; i < num_modes; i++) {
+    eigenvalues_[i] = eigenvalues_eigen[num_modes - 1 - i];
   }
+  eigenvectors_ = eigenvectors_.rowwise().reverse().eval();
 
-  float sum = 0.0;
-  for (unsigned int n = 0; n < num_samples_; n++) {
-    sum += eigenvalues_[(num_samples_ - 1) - n];
+  // eigenvalues are now in ascending order (smallest to largest)
+  double sum = 0.0;
+  for (int n = 0; n < num_modes; n++) {
+    sum += eigenvalues_[n];
   }
 
-  float sum2 = 0.0;
-  bool found = false;
-  for (unsigned int n = 0; n < num_samples_; n++) {
-    sum2 += eigenvalues_[(num_samples_ - 1) - n];
+  // percent_variance_by_mode_ accumulates from largest eigenvalue (at end) to smallest
+  // to match the old behavior
+  percent_variance_by_mode_.clear();
+  double sum2 = 0.0;
+  for (int n = num_modes - 1; n >= 0; n--) {
+    sum2 += eigenvalues_[n];
     percent_variance_by_mode_.push_back(sum2 / sum);
-
-    if ((sum2 / sum) >= 0.95 && found == false) {
-      found = true;
-    }
   }
 
   return 0;

Studio/main.cpp

@@ -10,6 +10,8 @@
 #include <QSurfaceFormat>
 #include <iostream>
 
+#include "Profiling.h"
+
 // itk
 #include <itkMacro.h>
 
@@ -77,8 +79,8 @@ static void new_log() {
 int main(int argc, char** argv) {
   // tbb::task_scheduler_init init(1);
 
+  TIME_SCOPE("ShapeWorksStudio");
   try {
-
 #ifdef Q_OS_MACOS
     // Prevent cursor crashes on Apple Silicon
     qputenv("QT_MAC_DISABLE_NATIVE_CURSORS", "1");
@@ -148,7 +150,9 @@ int main(int argc, char** argv) {
       app.file_open_callback_(app.stored_filename_);
     }
 
-    return app.exec();
+    auto rc = app.exec();
+    TIME_FINALIZE();
+    return rc;
   } catch (itk::ExceptionObject& excep) {
     std::cerr << excep << std::endl;
   } catch (std::exception& e) {