ENH: Add Eigen-backed SparseLUSolverTraits and deprecate VNLSparseLUSolverTraits

[hjmjohnson]

Adds an Eigen-backed EigenSparseLUSolverTraits (Eigen::SparseLU) alongside the existing VNLSparseLUSolverTraits, migrates the QuadEdgeMesh consumers to it, and deprecates the VNL traits behind ITK's three-state legacy policy. Part of the VNL→Eigen numerics direction (#6403, #6230). No VNL code is removed — the VNL path stays fully functional through the deprecation window.

Eigen is proven faster and more accurate on both axes (see report), and the two engines are proven numerically equivalent by a new GoogleTest, so the migration is safe.

Commits (additive; nothing deleted/renamed)

Commit	What
ENH: Add EigenSparseLUSolverTraits prototype with vxl-derived test coverage	New Eigen traits + GTest scavenged from vxl test_sparse_lu
DOC: Correct IsDirectSolver() Doxygen in VNLSparseLUSolverTraits	doc fix
ENH: Add VNL vs Eigen sparse-LU equivalence GoogleTest	cross-checks both engines agree (relative metric)
ENH: Add MatVecMult to sparse-LU traits and route filter through it	both traits gain MatVecMult; the Soft-constraints filter routes M^T·B through the traits so it is backend-neutral
ENH: Migrate QuadEdgeMesh consumers to EigenSparseLUSolverTraits	the QuadEdgeMesh filter tests/docs select the Eigen traits
ENH: Deprecate VNLSparseLUSolverTraits with a three-state guard	silent default / ITK_LEGACY_REMOVE warning / ITK_FUTURE_LEGACY_REMOVE error, pointing at the Eigen replacement

Performance proof — Eigen passes the gate on BOTH axes

2-D Laplacian systems (representative of the QuadEdgeMesh parameterization matrices). Single-threaded, median of 13, Apple M3 Max / clang 19.1.7. Lower is better in every column; fwd = forward error, resid = relative residual.

system (N)	vnl ms	vnl fwd	vnl resid	eigen ms	eigen fwd	eigen resid	speedup
256	0.576	2.2e-13	1.3e-14	0.217	7.6e-15	1.1e-14	2.7×
1024	9.63	1.9e-13	7.0e-15	0.920	3.8e-14	5.1e-15	10.5×
2304	46.5	1.8e-13	4.0e-15	2.40	3.7e-14	3.0e-15	19.4×
4096	147.9	1.6e-13	2.7e-15	4.92	4.1e-14	1.9e-15	30.1×
2304 ill-cond	46.3	1.7e-07	4.3e-15	2.40	4.2e-08	3.0e-15	19.3×

Eigen wins speed in 5/5 (margin grows with N) and accuracy in 5/5 (forward and residual). Numerical equivalence VNL↔Eigen holds to 1e-9 relative (well-conditioned) / 1e-6 (ill-conditioned), covered by VnlEigenSparseLUEquivalenceGTest.

Ecosystem usage scan — who actually uses this solver

Searched all source in a 32-repo downstream forest (ITK + ANTs, BRAINSTools, Slicer, elastix, MITK, SimpleITK, RTK, PlusLib, IGSIO, OpenIGTLink, vtkAddon, SCIFIO, CudaCommon, WebAssemblyInterface, TubeTK, c3d, … and the ITK remote modules):

Consumer	Real usage of vnl_sparse_lu / VNLSparseLUSolverTraits
ITK	The QuadEdgeMesh filters + tests, the traits header + its GTest (the surface migrated here)
All 31 non-ITK repos	None. Only one non-ITK reference exists — elastix itkThinPlateSplineTransformPerformanceTest.cxx, and all 4 hits are commented-out dead code (// #include, a never-built invert() experiment)
Condition API (rcond/max_error_bound/solve_transpose/estimate_condition)	Only vxl's own test_sparse_lu.cxx

So this sparse-LU solver is, in practice, consumed only by ITK's QuadEdgeMesh path (now on Eigen) and the algorithm's own unit test — the deprecation strands no external caller.

Test status

Full local build (4251 targets, 0 errors) and the complete ITK test suite 4428/4428 passing on macOS/AppleClang via the pinned conda toolchain. The 32 LU + QuadEdgeMesh tests (restored vnl_algo_test_sparse_lu, both traits' GTests, the equivalence test, 24 QuadEdgeMesh consumers) are green; the deprecation guard's three states (silent/warning/error) were verified by a preprocessor probe.

[hjmjohnson]

ITK.Linux failed at the Azure dashboard step (exit 255) while CDash for this commit shows 0 compile errors, 0 test failures, 4608 tests pass — the known "dashboard treats a compiler warning as fatal" flake. The single warning in the LegacyRemovedDebugCxx20 build is not introduced by this PR: building every file this PR touches (the ITKCommon GTest driver and the QuadEdgeMesh test driver) locally under -DITK_LEGACY_REMOVE=ON produces zero warnings, and the new VNLSparseLUSolverTraits deprecation guard is correctly suppressed by ITK_LEGACY_TEST in the tests that include it. Re-running ITK.Linux.

[hjmjohnson]

/azp run ITK.Linux

[greptile-apps]

Greptile Summary

This PR introduces EigenSparseLUSolverTraits (backed by Eigen::SparseLU) as the new default sparse LU solver for ITK's QuadEdgeMesh parameterization and deformation filters, migrates all three consumer tests to it, and deprecates VNLSparseLUSolverTraits behind a three-state preprocessor guard. The VNL path remains fully functional; the guard is silent by default, warns under ITK_LEGACY_REMOVE, and errors under ITK_FUTURE_LEGACY_REMOVE.

• New EigenSparseLUSolverTraits: mirrors the VNL traits API with IsDirectSolver, InitializeSparseMatrix, FillMatrix, AddToMatrix, MatVecMult, and six Solve overloads (single-RHS, 2-RHS, 3-RHS, plus solver-reuse variants); covered by a new GTest file with six functional tests and a ported vxl test suite.
• Soft-constraints filter now backend-neutral: three direct Mt.mult() VNL calls replaced with SolverTraits::MatVecMult, allowing the filter to work unchanged with both traits.
• Equivalence GTest (VnlEigenSparseLUEquivalenceGTest) cross-validates VNL and Eigen solutions to 1e-9 relative on well-conditioned 2-D Laplacians and 1e-6 on an ill-conditioned system.

Confidence Score: 4/5

Safe to merge; the VNL path is preserved and the new Eigen path is validated by both a unit GTest suite and a cross-solver equivalence test across multiple grid sizes.

The core correctness story is solid: equivalence is verified numerically, all three consumer tests pass with the new traits, and the deprecation guard is implemented correctly. The only actionable issues are non-blocking: the three matrix-taking Solve overloads in EigenSparseLUSolverTraits make a redundant copy of the input matrix before constructing SparseLU (which copies again internally), and the newly added MatVecMult method has no isolated unit test in either trait's GTest file. Neither affects correctness for the existing use cases.

The three Solve(const MatrixType&, ...) overloads in EigenSparseLUSolverTraits.h warrant a second look for the double-copy pattern.

Important Files Changed

Filename	Overview
Modules/Core/Common/include/EigenSparseLUSolverTraits.h	New Eigen::SparseLU-backed traits class; API mirrors VNLSparseLUSolverTraits. All three matrix-taking Solve overloads make an unnecessary copy of the input matrix before constructing SparseLU (which copies again internally), doubling allocation cost.
Modules/Core/Common/include/VNLSparseLUSolverTraits.h	Adds three-state deprecation guard (silent/warn/error), corrects the IsDirectSolver() doc comment, and adds MatVecMult. Deprecation logic is clean and correct.
Modules/Core/Common/test/EigenSparseLUSolverTraitsGTest.cxx	New GTest file with ConvertedLegacyTest wrapper plus three dedicated test cases ported from vxl; MatVecMult is not directly unit-tested.
Modules/Core/Common/test/VnlEigenSparseLUEquivalenceGTest.cxx	New equivalence test verifying VNL and Eigen solutions agree to 1e-9 (well-conditioned) and 1e-6 (ill-conditioned) relative tolerance across grid sizes; solid coverage.
Modules/Filtering/QuadEdgeMeshFiltering/include/itkLaplacianDeformationQuadEdgeMeshFilterWithSoftConstraints.hxx	Replaces three direct VNL-specific Mt.mult() calls with SolverTraits::MatVecMult(), making the filter backend-neutral. Change is minimal and correct.
Modules/Core/Common/test/VNLSparseLUSolverTraitsGTest.cxx	Adds #define ITK_LEGACY_TEST before the VNL include to opt out of the new deprecation warning; one-line change, correct.

Sequence Diagram

%%{init: {'theme': 'neutral'}}%%
sequenceDiagram
    participant Filter as QuadEdgeMeshFilter
    participant Traits as EigenSparseLUSolverTraits
    participant Eigen as Eigen::SparseLU

    Filter->>Traits: InitializeSparseMatrix(N, N)
    Filter->>Traits: FillMatrix / AddToMatrix (builds M)
    Filter->>Traits: MatVecMult(Mt, Bx, Cx)
    Traits->>Eigen: "oX = iA * iB (dense result)"
    Filter->>Traits: Solve(A, Cx, Cy, Cz, X, Y, Z)
    Traits->>Eigen: SparseLU(A) — copies + compresses + factorizes
    Eigen-->>Traits: factorization done, info()
    Traits->>Eigen: solver.solve(Cx) → X
    Traits->>Eigen: solver.solve(Cy) → Y
    Traits->>Eigen: solver.solve(Cz) → Z
    Eigen-->>Traits: back-substitution results
    Traits-->>Filter: "(X, Y, Z), success = info()==Success"
    Filter->>Filter: update mesh point positions

Loading

%%{init: {'theme': 'base', 'themeVariables': {"darkMode": true, "background": "#0d1117", "primaryColor": "#21262d", "primaryTextColor": "#e6edf3", "primaryBorderColor": "#8b949e", "lineColor": "#8b949e", "textColor": "#e6edf3", "edgeLabelBackground": "#161b22", "actorBkg": "#21262d", "actorBorder": "#8b949e", "actorTextColor": "#e6edf3", "actorLineColor": "#8b949e", "signalColor": "#8b949e", "signalTextColor": "#e6edf3", "noteBkgColor": "#373320", "noteBorderColor": "#d4a72c", "noteTextColor": "#f0e6c0", "labelBoxBkgColor": "#21262d", "labelBoxBorderColor": "#8b949e", "labelTextColor": "#e6edf3", "loopTextColor": "#e6edf3", "activationBkgColor": "#30363d", "activationBorderColor": "#8b949e"}}}%%
sequenceDiagram
    participant Filter as QuadEdgeMeshFilter
    participant Traits as EigenSparseLUSolverTraits
    participant Eigen as Eigen::SparseLU

    Filter->>Traits: InitializeSparseMatrix(N, N)
    Filter->>Traits: FillMatrix / AddToMatrix (builds M)
    Filter->>Traits: MatVecMult(Mt, Bx, Cx)
    Traits->>Eigen: "oX = iA * iB (dense result)"
    Filter->>Traits: Solve(A, Cx, Cy, Cz, X, Y, Z)
    Traits->>Eigen: SparseLU(A) — copies + compresses + factorizes
    Eigen-->>Traits: factorization done, info()
    Traits->>Eigen: solver.solve(Cx) → X
    Traits->>Eigen: solver.solve(Cy) → Y
    Traits->>Eigen: solver.solve(Cz) → Z
    Eigen-->>Traits: back-substitution results
    Traits-->>Filter: "(X, Y, Z), success = info()==Success"
    Filter->>Filter: update mesh point positions

Loading

Comments Outside Diff (1)

1. Modules/Core/Common/test/EigenSparseLUSolverTraitsGTest.cxx, line 474-475 (link)

   MatVecMult has no direct unit test in either traits GTest file

   MatVecMult is newly introduced in this PR for both EigenSparseLUSolverTraits and VNLSparseLUSolverTraits, yet neither GTest file contains a test that directly exercises it. It is only covered indirectly through itkLaplacianDeformationQuadEdgeMeshFilterWithSoftConstraintsTest. A standalone GTest case (e.g., verify oX == A * b component-wise for a small known matrix) would give isolated failure diagnostics if the implementation ever regresses.

   Note: If this suggestion doesn't match your team's coding style, reply to this and let me know. I'll remember it for next time!

_{Reviews (1): Last reviewed commit: "ENH: Deprecate VNLSparseLUSolverTraits w..." | Re-trigger Greptile}

[blowekamp]

Are there other external uses of VNL's LU solver?

I am not sure it is worth us to add an ITK specific numeric interfaces. It seem like it is more to maintain, and may make changes harder in the future. Maybe we just update the QuadEdge to directly use Eigen?

[hjmjohnson]

Good questions.

Other external uses of the VNL LU solver: across the downstream testbed (ITK + ANTs, BRAINSTools, Slicer, elastix, SimpleITK, and ~257 Slicer remote modules) I found one inactive consumer of vnl_sparse_lu — commented-out experimental code in an elastix performance test. No external use of VNLSparseLUSolverTraits either. But the testbed can't see every ITK consumer, so I'd read that as "none we can observe," not proof none exist.

On adding an ITK-specific numeric interface: VNLSparseLUSolverTraits and VNLIterativeSparseSolverTraits already exist as public classes — they're the documented TSolverTraits template parameter of the QuadEdge mesh filters. So this isn't a new ITK numeric interface; it adds an implementation of the existing contract. Hardcoding Eigen straight into QuadEdge would drop that public extension point and could silently break any downstream that passes its own traits.

Changes pushed (cf5b842):

• Renamed the new traits to a backend-neutral SparseLUSolverTraits — that it's implemented with Eigen::SparseLU is now documented in the class Doxygen rather than baked into the public name, so the engine can change later without renaming the public type.
• VNLSparseLUSolverTraits is deprecated from the start via a 3-state legacy guard pointing at SparseLUSolverTraits, so anyone relying on the VNL traits gets a warned migration window instead of a hard break — a deliberate crutch for the use cases we can't enumerate.
• The QuadEdge filters keep TSolverTraits as the public extension point (it has no default — callers pick the traits); the in-tree tests/docs now use SparseLUSolverTraits.