PaperlayerView exact revision on arXiv

Source-linked AI summary

Fast and Stable Arbitrary Precision Approximate Matrix Multiplication

Ioana Dumitriu, Olga Holtz, Abraham Hsu, Nicholas Liskij, Yoav Moran, Ryan Schneider, Oded Schwartz, Shuncheng Yuan

arXiv:2610.13905v1math.NA

TL;DR

APA algorithms can use fewer scalar products than exact multiplication, but their numerical properties are less understood and exactification requires repeated evaluations. The paper develops an APA search procedure and analyzes alternative-basis and exactification schemes, finding a rank-31 APA algorithm for ⟨2, 4, 5⟩ multiplication.

  • Problem

    APA algorithms may use fewer scalar multiplications, but their numerical properties are less understood and exactification requires repeated evaluations.

    Source: §1 · pp. 2–3

  • Method

    The paper searches for APA routines by transforming exact algorithms and deleting a rank-one term, and extends alternative-basis and stability analyses to APA schemes.

    Source: §1.1 · p. 4, §1 · p. 4

  • Results

    The search found a rank-31 APA algorithm for ⟨2, 4, 5⟩ multiplication, lowering that size's border-rank upper bound to 31.

    Source: §2.2.2 · p. 13

  • Takeaways & Limitations

    The paper presents APA search and numerical analysis as ways to renew interest in the theory and implementation of APA matrix multiplication.

    Source: §6 · p. 31

    The stability analysis uses the maximum-entry norm; extending the APA results to arbitrary norms is left for future work.

    Source: §4 · p. 20

Abstract

from arXiv · show

We present several new and surprising results on the theory and practice of Arbitrary Precision Approximate (APA) matrix multiplication, which refers to algorithms that output merely approximations to matrix multiplication even in exact arithmetic. While APA algorithms have been around for nearly half a century, and despite the fact that they have potentially lower (asymptotic) cost than many fast exact matrix multiplication schemes, they have not been studied as closely, to the extent that relatively little is known about their numerical properties. In this paper, we shine a spotlight on these algorithms by presenting three important results. First, we provide a new search algorithm for APA routines that is easy to use, automatically generates encoding matrices with functional entries, and can match or beat the best known APA algorithms for several sizes. In particular, we improve the upper bound on the border rank of $\langle 2,4,5 \rangle$ matrix multiplication. Second, we combine this insight with alternative bases to minimize the number of additions in an APA algorithm, thereby improving its leading coefficient. Finally, we prove the same kind of stability that characterizes fast exact matrix multiplication for an algorithm of Bini [Calcolo 1980] that computes an exact product of two matrices using a modified APA scheme (with or without alternative bases).

1 Introduction

APA algorithms can reduce scalar multiplication counts, but exact use requires repeated evaluations and their search and numerical properties have been less studied. The paper introduces a search method, adapts alternative bases, and analyzes stability for an APA-based exactification procedure.

1 Introduction

  • APA routines can use fewer scalar products than exact algorithms, but exactification requires repeated evaluations and adds a logarithmic complexity factor.The introduction also identifies poorly documented search methods and open stability questions as barriers to practical use.Source: §1 · pp. 2, 4

1.1 Our Contributions

  • A symbolic De-Groote transformation followed by deleting a rank-one term yields an APA routine when the deleted factors' product is nilpotent.The search uses the flip graph to find suitable exact algorithms and produces routines with functional entries.Source: §1.1 · p. 4, §1 · p. 4
  • The search produces APA routines for multiple matrix sizes and improves a previously best-known rank.The introduction describes the search as generating decompositions with functional entries.Source: §1.1 · p. 4, §1 · p. 4
  • Alternative bases sparsify APA algorithms and improve their leading arithmetic-complexity coefficient.Table 2 reports coefficient improvements for several schemes obtained with this approach.Source: §1.1 · p. 5, Table 2 · p. 6
  • The Bini approximate-to-exact procedure has a weak forward stability bound under mild assumptions, including fast basis transformations.The paper presents these as its first stability bounds for APA matrix multiplication.Source: §1.1 · p. 5

1.2 Notation and Conventions

  • Encoding matrices U, V, and W record input and output roles; APA encodings depend on ε, and sparsity can indicate practical efficiency.For implementation, the number of nonzero remainder terms—not merely the largest exponent—determines how many evaluations recover the exact product.Source: §1.2 · pp. 7–8

2 A New Search Algorithm for APA Schemes

The paper turns exact matrix-multiplication decompositions into lower-rank APA schemes by using symbolic basis changes to identify removable terms. A flip-graph search finds suitable decompositions, including a rank-31 APA algorithm for ⟨2, 4, 5⟩.

2 A New Search Algorithm for APA Schemes

  • A symbolic De-Groote transformation followed by deletion of O(ε) terms converts an exact decomposition into an APA scheme that can approach arbitrary accuracy as ε → 0.The method starts from an exact algorithm and requires suitable transformations that make deleted terms vanish with ε.Source: §2 · p. 8

2.1 Sufficient Conditions for Rank Improvements

  • If the transformed factors of l terms have entries with zero constant term, those terms can be removed, yielding an APA algorithm of rank r − l and a border-rank upper bound of r − l.A sufficient condition for removing one term is that its product UVW is nilpotent, which guarantees suitable symbolic basis transformations.Source: §2.1 · p. 9

2.3 Generalization via Null Cone Membership Testing

  • Multiple terms can be removed together only when the same symbolic basis transformations make them all O(ε); the generalization via null-cone membership testing produced no further border-rank improvements in practice.The authors describe their search approach as less general than null-cone methods, while reporting no additional improvements from that broader route.Source: §2.1 · p. 10, §2.3 · p. 13

2.1.1 ⟨2, 2, 3⟩Example

  • In the ⟨2, 2, 3⟩ example, the sixth term has a nilpotent product and is deleted, producing a ⟨2, 2, 3; 10⟩ APA algorithm.The example’s transformed tensor product is O(ε) even though not every transformed factor has entries in εC[ε].Source: §2.1.1 · pp. 11–12

2.2 Implementation Details and Results

  • Nilpotent terms were rare in the initial repository search, so the authors generated more candidates with flip-graph searches and checked them for the removal conditions.They also note that other exact-algorithm search methods could replace the flip graph.Source: §2.2 · p. 12

2.2.1 Flip Graph Search for Exact Decompositions

  • The search uses flip-graph walks to find exact decompositions meeting the term-removal conditions, then lifts them to characteristic zero before constructing symbolic basis transformations.If lifting fails, the algorithm returns to the search; the flip-graph approach was chosen for ease of use and minimal computing resources.Source: §2.2.1 · pp. 12–13, §2.2 · p. 12

2.2.2 Additional Results

  • Rank 31 is the first APA algorithm found for ⟨2, 4, 5⟩, lowering the border-rank upper bound for that size to 31.The search ran for about a week using 23 threads, with promising decompositions selected as seeds for restarted searches.Source: §2.2.2 · p. 13

3 Alternative Basis APA Matrix Multiplication

This section develops alternative-basis and term-deletion methods for sparsifying APA encoding matrices, with stated guarantees on sparsity and weight. It also describes their application to batched Bini approximate-to-exact computation and evaluates recursive leading coefficients and sparsification runtimes.

3 Alternative Basis APA Matrix Multiplication

  • APA basis transformations reduce encoding-matrix nonzeros to lower the leading coefficients in arithmetic and IO complexity.The paper presents a greedy construction of these symbolic transformations that depends on ε.Source: §3 · p. 13

3.1 Background and Motivation

  • Alternative-basis APA multiplication factors each encoding matrix into a basis transformation and a sparser matrix, then applies the bilinear form to transformed inputs when the transformations are cheap to compute.The transformations can be applied recursively to matrices of arbitrary size.Source: §3.1 · p. 14

3.2 Symbolic Sparsification

  • The symbolic greedy algorithm builds basis transformations by adding linearly independent vectors that minimize the weight of their product with an encoding matrix.It tests increasing weight values, promoting optimal sparsity; the weight measure can penalize variable entries through a tunable parameter.Source: §3.2 · p. 15

3.2.1 Implementation Details

  • Under Theorem 3.2’s conditions, admissible rational-function solutions have Laurent-polynomial solutions, making the symbolic search well-posed without sacrificing sparsity.The theorem also bounds the exponents in those solutions, which can narrow the solver’s search range.Source: §3.2.1 · pp. 16–18

3.2.2 Optimality

  • Symbolic sparsification is guaranteed to achieve optimal sparsity under the stated condition on Next-Vector, and supports alternative definitions of the weight measure.The guarantee assumes Next-Vector returns an admissible vector whenever one exists; Theorem 3.2’s first condition for the encoding matrix ensures this availability.Source: §3.2.2 · p. 18

3.2.3 Sparsifying Algorithms Obtained by Term Deletion

  • For APA algorithms whose parameter dependence comes from a De-Groote transformation of an exact decomposition, alternative bases can make all three encoding matrices independent of ε without sacrificing measured sparsity.Source: §3.2.3 · p. 19
  • Under invertible symbolic alternative-basis transformations with Laurent-polynomial entries, equivalent transformations make the encodings constant; for p ≥ 1, their weights do not increase.The method enables applying basis transformations in parallel before batched evaluation of the bilinear algorithm, followed by Bini’s approximate-to-exact procedure.Source: §3.2.3 · p. 19

3.3 Results

  • Leading coefficients were computed from symbolic base-case operation counts and minimized over recursive compositions and permutations using the framework of.For APA weight calculations, p = 2; the reported coefficient is the minimum returned by the recursion-order optimization.Source: §3.3 · p. 19
  • Most encoding matrices were sparsified in under an hour, while three took more than a day: two took two days and one took three.The experiments ran on a laptop with an Intel Core Ultra 9; entry complexity appeared to affect runtime more than matrix size.Source: §3.3 · p. 19

4 APA Matrix Multiplication is Stable

This section analyzes stability in APA matrix multiplication, proving weak forward stability for an approximate-to-exact alternative-basis procedure and extending the analysis to Algorithm 5 by separating truncation and floating-point errors.

4 APA Matrix Multiplication is Stable

  • The stability bounds use the maximum-entry norm; the paper leaves APA bounds for arbitrary norms and scaling techniques to future work.The cited arbitrary-norm results are described as potentially tighter than bounds obtained through norm equivalence.Source: §4 · p. 20

4.1 Preliminaries

  • The analysis builds on a stability bound for a single recursive alternative-basis APA evaluation, extending Schwartz et al.’s framework by allowing the encoding matrices and basis transformations to depend on ε.The analysis defines stability factors and a prefactor for this setting.Source: §4.1 · pp. 21–22

4.2 Main Error Bound

  • Theorem 4.6 establishes a weak forward error bound for Algorithm 4 under assumptions on the inputs, evaluation points, and floating-point subroutines.The conditions include normalized inputs and evaluation points, accurate computation of the Vandermonde solution and encoding entries, and no overflow or underflow.Source: §4.2 · pp. 22–23

4.3 Parameter Selection for Improved Stability

  • Roots of unity simplify the stability prefactor to 1/p, while rounded roots can support the stability bound but require complex encoding matrices that may increase arithmetic cost.The stability of the approximate-to-exact procedure is characterized by the behavior of the APA stability factor and prefactor near the unit circle.Source: §4.3 · pp. 24–25

4.4 Extensions

  • The truncated procedure uses r < p evaluations and cancels the first r − 1 remainder terms, trading a nonzero remainder and tighter parameter spacing for improved efficiency, especially when r = O(1).When all evaluation parameters satisfy |εj| < 1, the output can be expected to be close to the exact product because the remaining terms involve powers of ε of order r or higher.Source: §4.4 · p. 25
  • With p-th roots of unity as evaluation parameters, Algorithm 2 outputs can be combined using Vandermonde-derived coefficients whose magnitudes are 1/p.Source: §4.4 · pp. 25–26
  • For Algorithm 5, the error separates into truncation loss and floating-point loss, and appropriate choices of r and evaluation parameters can make the latter arbitrarily small.The analysis uses Theorem 4.6 to bound the floating-point implementation even though Algorithm 5's exact-arithmetic output is not the exact product.Source: §4.4 · p. 26
  • Under the stated parameter choice, bound (4.15) essentially implies weak forward stability; scaling the roots as epsilon_j = beta * e^(2pi i j/r), for 0 < beta < 1, further simplifies bound (4.16).Source: §4.4 · p. 26

5 Numerical Examples

The numerical experiments assess stability-bound predictions, alternative-basis behavior, and truncation, finding that the bounds describe practical stability while truncation trades some accuracy for speed.

5.4 Effectiveness of Truncation

  • Truncation can significantly improve runtime while trading mild accuracy loss for speed.When β is close to one and only a modestly accurate product is needed, evaluations can be reduced by a constant factor.Source: §5 · p. 27, §5.4 · p. 31

5.1 Experimental Setup

  • The experiments use recursive APA implementations, measure mean relative forward error over ten random input draws, and compare against high-precision symbolic products.The algorithms run in binary64; errors use the maximum-entry norm.Source: §5.1 · p. 27

5.2 Verifying Theorem 4.6

  • The first stability experiment compares a searched ⟨3, 3, 3; 21⟩ routine with Smirnov’s and Schönhage’s algorithms across three input distributions, using classical multiplication as a baseline.The “worse-case” inputs are constructed to increase error for given encoding matrices; the authors do not claim they are actual worst cases.Source: §5.2 · pp. 27–28
  • Theorem 4.6 and 4.7 bounds gauge relative stability, describe measured errors, and rank the three algorithms by expected accuracy.Figure 3 compares measured and theoretical accuracy across input distributions.Source: §5 · p. 27, Figure 3 · p. 28, §5.2 · p. 28

5.3 Stability with Alternative Bases

  • Alternative bases show a larger predicted stability loss but essentially no difference in practical performance.Table 4 compares errors with alternative bases against the standard basis; a ratio of one means equal accuracy.Source: §5 · p. 27, §5.3 · p. 30, Table 4 · p. 31

5.4 Effectiveness of Truncation

  • For the ⟨2, 2, 3; 10⟩ routine, recursive remainder degrees are 3, 6, 9, 12, 15, 18, 21, and 24, so r=25 is not required for exactness.Exactness also occurs at lower r values when r divides none of the degrees present.Source: §5.4 · p. 31

6 Conclusion

The paper presents a search algorithm for APA schemes and numerical analysis of their practical performance. High-performance implementations and stability proofs for other approximate matrix multiplication approaches remain future work.

6 Conclusion

  • The paper contributes a new search algorithm for APA schemes and numerical analysis aimed at understanding their practical performance.The stated goal is to renew interest in APA matrix multiplication from theory and implementation perspectives.Source: §6 · p. 31
  • High-performance APA implementations competitive with fast exact multiplication remain a future-work goal; sparsification and truncation have not yet been tested in that setting.The authors also identify possible applications where lower precision is required.Source: §6 · p. 31
  • Numerical stability remains open for approximate matrix multiplication algorithms derived from floating-point least-squares tensor-decomposition solutions.The passage cites existing examples of such algorithms [35, 36, 56].Source: §6 · p. 32

A Technical Results for Section 4

The appendix collects technical results used to prove the paper’s main stability bound.

A Technical Results for Section 4

A.1 Vandermonde Bounds

This appendix develops floating-point evaluation and perturbation bounds for quantities associated with a Vandermonde system.

A.1 Vandermonde Bounds

  • The appendix treats a Vandermonde system formed from distinct, nonzero points ε1 through εp.Source: §A.1 · p. 36
  • Under stated mild assumptions, the αj values can be computed accurately in relative terms by directly evaluating the expression in floating-point arithmetic.The listed assumptions include normalized floating-point points, controlled evaluation without over- or underflow, and a condition involving p and unit roundoff.Source: §A.1 · p. 36
  • The appendix also gives an exact-arithmetic perturbation bound for the solutions at original and relatively perturbed evaluation points, subject to a stated condition.The passage defines δ as the size of the relative perturbation but does not include the resulting bound itself.Source: §A.1 · pp. 36–37

A.2 Perturbation Bound for Bilinear Algorithms

This appendix proves a perturbation bound for recursively evaluated alternative-basis bilinear algorithms in exact arithmetic. The argument treats entrywise perturbations to the encoding matrices and basis transformations and proceeds by induction on recursion depth.

A.2 Perturbation Bound for Bilinear Algorithms

  • The section proves a perturbation bound for recursively evaluated alternative-basis bilinear algorithms in exact arithmetic.The proof proceeds by induction on recursion depth, with matrix multiplication as the base case.Source: §A.2 · pp. 38, 41
  • The bound concerns entrywise relative perturbations of size at most δ to the encoding matrices and basis transformations, with recursion defaulting to matrix multiplication after L steps.The input matrices are partitioned into blocks, assuming m0, k0, and n0 divide m, k, and n.Source: §A.2 · p. 38
  • The proof relates recursive block evaluation to applying the basis transformations to the inputs first, then evaluating the bilinear algorithm and applying the output transformation.The equivalence is stated in exact arithmetic and extended recursively to the other encoding matrices and basis transformations.Source: §A.2 · p. 39
  • For a hybrid recursive algorithm, the stability factor associated with its top-level perturbed matrices is at most (1 + δ)^6E.The section identifies the maximum-entry norm and E as the stability factor for the encoding matrices and basis transformations.Source: §A.2 · pp. 38, 40

B Encoding Matrices for ⟨2, 4, 5; 31⟩

The appendix gives an APA scheme for multiplying a 2 × 4 matrix by a 4 × 5 matrix, establishing 31 as an upper bound on the border rank. Its approximation has residual terms from powers ε^3 through ε^8.

B Encoding Matrices for ⟨2, 4, 5; 31⟩

B.1 Standard Basis Encoding Matrices

This subsection presents the encoding matrices for the scheme in the standard basis.

B.1 Standard Basis Encoding Matrices

  • The subsection presents the scheme’s encoding matrices in the standard basis.Source: §B.1 · pp. 41–42

B.2 Residual Matrices

This subsection provides the residual matrices associated with the scheme.

B.2 Residual Matrices

  • The subsection gives the remainder matrices R3, R4, R5, R6, R7, and R8.Source: §B.2 · p. 43
  • The displayed residual entries include combinations of input entries such as a23, a24, b13, b23, b34, and b45.Source: §B.2 · p. 43
  • Other displayed residual entries involve a23, a24, a13, a14, b24, b25, and b14.Source: §B.2 · p. 43

B.3 Sparsified Encoding Matrices

The section presents basis transformations and constant encoding matrices obtained by sparsifying the algorithms, using the specialization from Section 3.2.3.

B.3 Sparsified Encoding Matrices

  • Sparsifying the algorithms yields basis transformations and constant encoding matrices, with the specialization of Section 3.2.3 applying.The section introduces the matrices that follow.Source: §B.3 · p. 43
Loading 2610.13905v1…