MiniTensor: modernize the general log() and exp() algorithms - #15515
Merged
Conversation
log(), dimensions 3 and below: hybrid of two methods, chosen by the distance from the identity. - Within the Pade trust region (norm_1(A - I) <= 0.72) the logarithm is a direct Gauss-Legendre partial-fraction Pade evaluation on A - I, which is computed exactly. This preserves the relative accuracy of small logarithms, which any similarity-based method loses to the ~eps*norm(A) rounding of the transformation itself. - Away from the identity, a Schur-based inverse scaling and squaring in the spirit of Al-Mohy and Higham, Improved Inverse Scaling and Squaring Algorithms for the Matrix Logarithm, SIAM J. Sci. Comput. 34(4), 2012: real Schur decomposition (Givens Hessenberg reduction plus Francis double-shift QR with standardized 2x2 blocks and the classical difference form of the shift vector), exact quasi-triangular square roots in place of Denman-Beavers iterations, Pade evaluation on the quasi-triangular factor, and exact recomputation of the diagonal blocks of the logarithm from the eigenvalues. Inputs with eigenvalues on the closed negative real axis fail with a clean error instead of failing to converge. Larger dimensions keep the previous transformation-free path. Measured on 500-sample suites per parameter point against an 80-bit eigendecomposition reference: identical accuracy and speed near the identity, 1.8-4.3x more accurate and 1.4-4.7x faster elsewhere (condition numbers 10 to 1e14), no failures. exp(): parameter selection of Al-Mohy and Higham, A New Scaling and Squaring Algorithm for the Matrix Exponential, SIAM J. Matrix Anal. Appl. 31(3), 2009. Their theta thresholds make the norm-based order choice exact short-circuit evaluation of the full algorithm whenever norm_1(A) <= theta_13 = 4.25, so the norm(A^p)^(1/p) overscaling analysis and the rounding correction term run only above that threshold, where they prevent overscaling of nonnormal arguments. The correction term exits early through a conservative bound and computes its matrix powers by binary powering. Median cost is unchanged; the extreme-norm cases pay up to ~1.9x for up to 1.8x accuracy gains at strong stretches. Also fixes a decimal-point typo (9,4e-0) in the scaling_squaring_theta table that silently shifted the entries above order 16. New tests: nonsymmetric logarithm with a complex-conjugate eigenvalue pair recovered through log(exp(L)), the analytic logarithm of a scaled 2x2 rotation, and agreement of exp() with the Taylor series on a nonsymmetric argument. Signed-off-by: Alejandro Mota <amota@sandia.gov>
Reproducing CI builds locallyCI builds can be reproduced locally in a container. Requirements: Python 3.9+, Podman, CMake Steps:
|
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
@trilinos/minitensor
Motivation
The general
log()used a transformation-free inverse scaling and squaring whose Denman–Villiers square-root iterations and determinant-based scaling are both expensive and fragile for ill-conditioned input (see #15513), andexp()used the 2005-era parameter selection with rounded threshold constants. This PR brings both up to the current published algorithms, specialized to the small dimensions MiniTensor serves.Changes
log(), dimensions <= 3 (larger dimensions keep the previous path):
norm_1(A - I) <= 0.72): direct Gauss–Legendre partial-fraction Padé evaluation on the exactly computed A − I. This preserves the relative accuracy of small logarithms, which any similarity-based method loses to the ~eps·norm(A) rounding of the transformation itself.exp(): parameter selection of Al-Mohy & Higham, SIAM J. Matrix Anal. Appl. 31(3), 2009. Their θ thresholds make the norm-based order choice provably exact short-circuit evaluation of the full algorithm whenever
norm_1(A) <= θ13 = 4.25; the norm(A^p)^(1/p) overscaling analysis and the rounding-correction term (with a conservative early exit and binary powering) run only above that threshold. Also fixes a decimal-point typo (9,4e-0) in thescaling_squaring_thetatable that silently shifted entries above order 16.Related Issues
Stacked on #15513 and includes its commit; this PR will be rebased to a single commit once #15513 merges. Follows #15514.
Testing