tests: establish systematic tolerance strategy for numerical operations

[ultimatile]

Summary

Numerical tolerance in tests is currently ad-hoc: most tests use the base eps from the fixture, while 9 linear algebra tests use eps * 100 without documented rationale. A systematic strategy is needed, but the appropriate tolerances depend on factors outside the conformance test's control (BLAS implementation, algorithmic choices, problem conditioning).

Current state

• Base epsilon: hardcoded 1e-10 (see fixture: epsilon should be type-aware based on real_t precision #29 for type-awareness fix)
• 9 tests use eps * 100: QR, LQ, SVD, trunc_SVD, exp (in linear_algebra.h)
• 98 tests use bare eps
• No documentation of why specific multipliers were chosen

Problem

The right tolerance depends on:

1. Operation category — exact operations (copy, reshape) vs. numerically sensitive decompositions (SVD, eig)
2. Problem size — larger matrices accumulate more floating-point error
3. BLAS/LAPACK implementation — different backends (OpenBLAS, MKL, cuSOLVER, etc.) use different algorithms with different numerical properties
4. Conditioning — the test matrices' condition numbers affect result accuracy

A conformance test should be tight enough to catch bugs but loose enough to not reject correct implementations that differ only in rounding.

Proposed approach

Step 1: Categorize operations by expected error scaling

Category	Examples	Expected error	Notes
Exact	copy, reshape, transpose, fill, zeros, eye	0	No floating-point arithmetic
Single arithmetic	scale, set_elem/get_elem	O(eps)	One multiply/assignment
Multi-element reduction	norm, trace, linear_combine	O(sqrt(N) * eps)	Summation error
Factorization	QR, LQ, SVD, eig, eigh	O(N * eps * cond)	BLAS/LAPACK dependent
Iterative/composite	exp, inverse	O(N^k * eps)	Algorithm dependent

Step 2: Design tolerance as a function, not a constant

Rather than hardcoding multipliers, provide a tolerance helper that takes the operation category and problem size:

enum class tol_category { exact, elementwise, reduction, factorization, iterative };

template <typename TenT>
auto tolerance(tci_test_fixture<TenT>& fix, tol_category cat, std::size_t N = 1) 
    -> tci::real_t<TenT>;

Step 3: Allow backend-specific override

Backends should be able to adjust tolerance multipliers if their BLAS implementation has known characteristics. This could be via:

• Specializing the fixture
• Compile-time defines (e.g., TCICT_TOL_FACTORIZATION_MULT=1000)
• Or a configuration mechanism TBD

Open questions

• Should we run a reference benchmark across multiple BLAS implementations (OpenBLAS, MKL, Apple Accelerate) to empirically determine reasonable multiplier ranges?
• Should tolerance be problem-size-dependent (tighter for small test matrices, looser for large)?
• How to handle the case where a backend's algorithm is fundamentally different (e.g., randomized SVD vs. divide-and-conquer)?
• Is the current eps * 100 for factorizations too tight or too loose? Need empirical data.

Non-goals

This issue is NOT about changing the current hardcoded values immediately. It's about designing a principled framework so that tolerance decisions are explicit, documented, and adjustable.

Related

• fixture: epsilon should be type-aware based on real_t precision #29 (type-aware epsilon — prerequisite)

[ultimatile]

Implementation plan (post-research, post-codex-review)

Goal

Introduce a tolerance categorization framework so that future empirical calibration of numerical-test tolerances is principled and adjustable. Non-goal: changing current hardcoded values immediately. This PR introduces the framework and migrates the 9 existing eps * 100 callsites; the framework preserves their numerical semantics exactly.

Files

1. include/tcict/fixture.h — additions only.
2. include/tcict/tests/linear_algebra.h — migrate 9 eps * 100 callsites; remove 4 now-unused auto eps declarations.

include/tcict/fixture.h additions

namespace tcict {

enum class tol_category {
  exact,          // No FP arithmetic (copy, reshape, fill, zeros, eye)
  elementwise,    // O(eps): single multiply/assignment (scale, set/get_elem)
  reduction,      // O(sqrt(N) * eps): summation (norm, trace, linear_combine)
  factorization,  // O(N * eps * cond): BLAS/LAPACK (QR, LQ, SVD, eig, eigh)
  iterative,      // O(N^k * eps): algorithm-dependent (exp, inverse)
};

template <typename TenT>
auto tolerance(tci_test_fixture<TenT>& fix, tol_category cat, std::size_t /*N*/ = 1)
    -> tci::real_t<TenT> {
  using real_type = tci::real_t<TenT>;
  const auto eps = fix.epsilon();
  switch (cat) {
    case tol_category::exact:         return real_type{0};
    case tol_category::elementwise:   return eps;
    case tol_category::reduction:     return eps;
    case tol_category::factorization: return eps * real_type{100};
    case tol_category::iterative:     return eps * real_type{100};
  }
  return eps;  // unreachable; switch is exhaustive over all enum values
}

}  // namespace tcict

Why a free function (not a fixture method)

A method on the primary tci_test_fixture<TenT> template would not be inherited by explicit specializations. The fixture's documented extension point is "Backends may specialize this for custom context creation (e.g. GPU) or epsilon" — a backend specialization that overrides epsilon() would silently lack a tolerance() method, breaking every callsite. A free function template that takes the fixture and calls fix.epsilon() works against any fixture (primary or specialized) because epsilon() is what backends are already documented to provide. Backends that want to override the tolerance computation itself can specialize the free function template.

N parameter

The signature accepts N for forward compatibility with future N-scaling of factorization / iterative tolerances, but v1 ignores it. Implementing N-scaling without empirical residuals across BLAS implementations would risk miscalibration; see Inconclusive / Deferred items.

Callsite migration in linear_algebra.h

Tests are inside namespace tcict::tests, so unqualified tolerance(fix, tol_category::factorization) resolves through the enclosing namespace.

Test	Sites	New category
test_qr	2 (Q*R reconstruction, Q†Q identity)	factorization
test_lq	2 (L*Q reconstruction, QQ† identity)	factorization
test_svd_reconstruction	1	factorization
test_exp_anti_hermitian	4 (4 element comparisons against cos/sin)	iterative

Total: 9 sites — matches the issue's listing.

Now-unused declarations to delete

After migration, eps becomes unused in 4 functions because every site that referenced it is the migrated site. Delete:

• test_qr — auto eps = fix.epsilon();
• test_lq — same
• test_exp_anti_hermitian — same
• test_svd_reconstruction — same

Impact list

• The 9 callsites in linear_algebra.h. No callers outside the test suite.
• Backends that specialize tci_test_fixture<T>: no change required. The free function template instantiates against any specialization that exposes epsilon(), which is the existing documented extension point.

Test plan

• Existing conformance tests must continue to pass with no change in numerical pass/fail behavior. The helper preserves exact eps * 100 semantics for both factorization and iterative.
• No new test cases (this is API refactoring).
• Verification: make integration-check before and after; expected: identical results.

Implementation guards

• Helper return type must equal tci::real_t<TenT> (matches epsilon()).
• switch is exhaustive over the 5 enum values; the trailing return eps is a defensive fallback for compilers that do not realize the switch is exhaustive (no default: is used so future additions to the enum will be caught by -Wswitch when enabled).
• tol_category::exact returns real_type{0} (literal zero), not epsilon().

Derivations (numerical-analysis grounding for the category table)

Standard backward-error scaling laws from Higham, Accuracy and Stability of Numerical Algorithms (2nd ed., 2002). Reproduced here for the plan reader; not novel.

Category	Scaling	Reference
exact	0	no FP arithmetic
elementwise	O(eps)	one multiply / assignment introduces ≤ 1 ulp
reduction	O(sqrt(N) * eps)	pairwise / Kahan summation; Higham §4
factorization	O(N * eps * cond)	QR / SVD / eig backward stability; Higham §19, §20
iterative	O(N^k * eps)	algorithm-dependent (matrix exp via scaling-and-squaring); Higham §10

This PR hardcodes per-category constants matching the current code (eps for the first three, eps*100 for the last two). N-scaling and condition-number scaling are deferred until empirical residuals are available.

Inconclusive / Deferred items

• deferred — Empirical calibration of multipliers
  reason: Open questions in the issue body — comparing residuals across BLAS implementations, deciding whether eps * 100 is too tight or too loose, choosing N-scaling formulas — all require running the conformance tests against real backends with various floating-point types and matrix sizes. The cytnx integration check (Add scripted local integration verification with cytnx backend #38) just landed; the cytnx side has not yet registered float / complex64 conformance instances. Out of scope per the issue's non-goals ("NOT about changing the current hardcoded values immediately").
  resolution-point: After the cytnx-side float / complex64 register PR runs successfully against make integration-check, gather residuals across BLAS implementations and revisit the multipliers (and the N-scaling formula).

• deferred — N-scaling implementation
  reason: The helper signature accepts N but ignores it. Implementing N-scaling without empirical data risks miscalibration.
  resolution-point: Tied to the empirical observation cycle above.

• deferred — bare eps callsite migration
  reason: 308 textual references to eps exist beyond the 9 eps * 100 sites (most are auto eps = fix.epsilon(); declarations and direct comparisons). Per-callsite category classification is a separate, larger refactor. Current bare-eps semantics are equivalent to tolerance(fix, tol_category::elementwise), so leaving them unchanged preserves current behavior without harm.
  resolution-point: Optional follow-up issue when category differentiation is empirically motivated for non-elementwise sites.

• deferred — Backend override beyond fixture specialization
  reason: The issue listed compile-time defines and config files as TBD alternative override mechanisms. The free-function form already supports backend override via free-function specialization, which mirrors the existing fixture-specialization pattern. Adding additional override mechanisms now would expand the API surface without empirical motivation.
  resolution-point: When a backend author requests a non-specialization override path.

Conflicts / dependencies

• Depends on fixture: epsilon should be type-aware based on real_t precision #29 (type-aware epsilon) — already merged.
• Independent of Add scripted local integration verification with cytnx backend #38 (just merged) — orthogonal change.

[ultimatile]

Status update

Landed

• Framework introduction — feat(tolerance): add tol_category framework for backward-error scaling #40 (commit a04a17d) added namespace-scope tol_category enum (exact / elementwise / reduction / factorization / iterative) and the free function tolerance(fix, cat, N=1) to include/tcict/fixture.h.
• Callsite migration (factorization / iterative tier) — the same PR migrated the 9 historical eps * 100 sites in include/tcict/tests/linear_algebra.h (test_qr ×2, test_lq ×2, test_svd_reconstruction ×1, test_exp_anti_hermitian ×4) to tolerance(fix, tol_category::{factorization|iterative}). The framework preserves the existing numerical semantics exactly: both categories still return eps * 100 in this version.

Empirical observation now flowing

The plan deferred per-multiplier calibration to the resolution-point "after the cytnx-side float / complex64 register PR runs successfully against integration-check". That has now happened: cytnx PR r-ccs-cms/tensor-computing-interface-backend-cytnx#12 merged on 2026-05-09 with cytnx_float / cytnx_complex64 instances registered through the bulk macros, expanding the suite from 243 cases (2 types) to 460 (4 types) on a real backend.

The single-precision SVD / trunc_svd / eig failures observed there are tracked at r-ccs-cms/tensor-computing-interface-backend-cytnx#10 — that issue is the empirical input for the calibration follow-up and explicitly notes the dual-direction question (cytnx-side numerical-stability fix vs. tcict-side multiplier calibration).

Still open (carried forward from the plan's Inconclusive / Deferred items)

Item	Status	Resolution-point
Empirical calibration of multipliers across BLAS impls	unlocked: data now collecting via cytnx #10	once enough residuals are gathered to rule for / against eps * 100 for factorization and iterative
N-scaling implementation	still deferred	tied to the calibration question above
Migration of remaining bare-eps callsites (the 308 references the plan deferred)	unchanged	optional follow-up if category differentiation is empirically motivated for non-elementwise sites
Backend override beyond fixture specialization	unchanged (no demand surfaced)	when a backend author requests a non-fixture override path

The original issue body remains the design statement; this comment is the journal entry for what landed in #40 and what unblocked / remained after.