Speed up powering a matrix by a large exponent when its minimal polynomial is small - #6495
Speed up powering a matrix by a large exponent when its minimal polynomial is small#6495fingolfin wants to merge 3 commits into
Conversation
…omial is small Powering a matrix by a large exponent reduces x^n modulo the characteristic polynomial, after a base change making the matrix block companion. If the minimal polynomial has degree e rather than d, both the Log2(n) polynomial multiplications and the final evaluation work with degree e, which is worth a factor of 2 to 60. Computing the minimal polynomial of a dense matrix costs about as much as the whole characteristic polynomial method, so probe first with a single spun random vector, abandoned once its order polynomial grows past the bound where the minimal polynomial stops being worth using. That order polynomial divides the minimal polynomial, so abandoning settles the question, and the probe does not show up in the timings. The break even points were determined over a grid of dimensions from 60 to 2000, representations GF(2), GF(5), GF(251), GF(257) and GF(3^10), minimal polynomial degrees from d/30 up to d, and exponents from 2^20 to 2^1000. They depend mostly on the representation: 6*e <= d for compressed matrices, 2*e <= d for the rest, where matrix multiplication is far more expensive than polynomial arithmetic. The choice is within a factor of 1.07 of the best of the three methods everywhere measured, and at least 1.5 times faster than before on a third of the cells. Matrix objects without access to their rows, such as those in IsGenericMatrixRep, now work here at all: the helper evaluating a polynomial at a matrix tested mutability by looking at the first row. With the preceding commit they reach the generic methods for the characteristic and minimal polynomial again, so powering them uses the polynomial method rather than failing with "row access unsupported". For this, Matrix_OrderPolynomialInner gained an optional degree bound, and the helpers of POW_MAT_INT moved out of its body to be shared. This commit was prepared with assistance from the AI tool Claude Code (benchmarking, analysis and drafting of the implementation and tests). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
fingolfin
added
release notes: use title
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## master #6495 +/- ##
==========================================
- Coverage 79.03% 78.97% -0.07%
==========================================
Files 685 685
Lines 293805 294269 +464
Branches 8664 8646 -18
==========================================
+ Hits 232220 232403 +183
- Misses 59786 60063 +277
- Partials 1799 1803 +4 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
|
@frankluebeck @ThomasBreuer any comments? |
What is d here? |
|
no comments from my side |
|
I'll not have time to look at this code for 4 weeks from now. But here are some quick comments: I don't understand most of the initial comment above. What is compared and benchmarked there? The task is to compute x^n for a d x d -matrix x over a finite field GF(q) for large exponents n. One idea is to conjugate x to sparse block triangular form with companion matrices in the diagonal blocks. Of course, when the minimal polynomial of x is known and of smaller degree than d one can do the last trick with the minimal polynomial instead. But is it worth to spent any effort to find the minimal polynomial? If the matrices in applications are close to random matrices I would say "no": by results of Neumann and Praeger on "Cyclic matrices" (meaning that characteristic and minimal polynomial are the same), the proportion of matrices which have minimal polynomial of degree < d is about 1/q^3. If I understand the initial comment correctly, a random vector is spinned up to find a factor of the minimal polynomial. But a random matrix is cyclic with high probablity and in that case the cyclic subspaces of most random vectors are the full space (and the sparse conjugate of x and its characteristic polynomial have been computed), (see also Neumann-Praeger for more precise statements. Do these comments fit with this pull request? |
…omial POW_MAT_INT probed for a small minimal polynomial and, when the probe succeeded, called MinimalPolynomial. That call dominated everything else: on a 1000 x 1000 permutation matrix over GF(2) with minimal polynomial of degree 4 it took 2.9 s, while reducing modulo the characteristic polynomial -- the method it was supposed to beat -- took 0.29 s. For d = 504 the minimal polynomial method was between 1.15 and 5.4 times slower over GF(2), and up to 1.4 times slower over GF(251). Any annihilating polynomial serves as the modulus, not just the minimal one, so use the probe's own order polynomial once Horner confirms that it kills the matrix, at a cost of e-1 matrix multiplications. Do that at the block companion conjugate rather than at the matrix itself unless e is tiny; the base change is needed by the fallback anyway, so a failed verification wastes nothing. For d = 504 and e from 2 to d/6 this is now faster than the characteristic polynomial method on every measured cell: 1.4 to 115 times over GF(2), 1.8 to 9 times over GF(251). On random matrices, which are cyclic with probability about 1-1/q^3 and hence never take this path, the probe costs 0.5% to 2% of the total. The break even points 6*e <= d and 2*e <= d are inherited from the previous cost model and are now conservative: over GF(2) and GF(251) the order polynomial method still wins at e = d/2. Prepared with assistance from the AI tool Claude Code (benchmarking, analysis and drafting of the implementation and tests). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
POW_MAT_INT spun a random vector to guess the minimal polynomial, threw
that spinning away, and then spun standard basis vectors from scratch to
build the base change, after which CharacteristicPolynomial spun a third
time. Frank Luebeck pointed out that for a cyclic matrix the first
spinning already is the base change.
So spin the random vector as the first Krylov chain of the base change
itself. For a cyclic matrix -- which a random matrix is with probability
about 1-1/q^3 -- that one chain spans everything, and its order
polynomial is the characteristic polynomial; in general the product of
the chains' order polynomials is. Either way the separate
CharacteristicPolynomial call is gone, and the guess costs nothing even
when it fails.
Since nothing is staked on the guess any more, the break even point need
no longer be fitted: verifying a candidate costs e-1 of the matrix
multiplications that the evaluation then saves d-e of, so it is worth
trying as soon as 2*e <= d. That replaces the measured 6*e <= d for
compressed matrices, widening their range from d/6 to d/2.
Timings for n = 2^500+1, milliseconds, best of 5, against the same
matrices; e = d is the random matrix case:
e=d/12 e=d/6 e=d/4 e=d/2 e=d
GF(2) d=504 133 107 106 108 109 before
36 38 49 78 82 after
GF(251) d=504 848 847 852 867 888 before
409 460 508 684 855 after
GF(2) d=1008 596 627 569 579 589 before
140 211 270 483 491 after
GF(251) d=1008 5300 5374 5441 5564 5695 before
3028 3331 3661 4756 5556 after
Matrix_OrderPolynomialInner gained an optional argument handing back the
Krylov chain it computes anyway, which lets POW_MAT_INT_ADDB -- a second
implementation of the same spinning -- be dropped.
Prepared with assistance from the AI tool Claude Code (benchmarking,
analysis and drafting of the implementation).
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
4 participants
Powering a matrix by a large exponent reduces
x^nmodulo the characteristic polynomial. Doing it modulo the minimal polynomial instead, when that has degreeerather thand, is worth a factor of 2 to 60.The minimal polynomial of a dense matrix costs about as much as the whole method, so a single spun random vector probes for it first, abandoned once its order polynomial passes the break even point. Depending on
e, the reduced polynomial is then evaluated at the matrix itself or at its block companion form.Break even points measured over dimensions 60 to 2000,
GF(2),GF(5),GF(251),GF(257),GF(3^10), minimal polynomial degreesd/30tod, and exponents2^20to2^1000:6*e <= dfor compressed matrices,2*e <= dotherwise, where matrix multiplication is far more expensive than polynomial arithmetic. The choice is within a factor of 1.07 of the best of the three methods on every measured cell, and at least 1.5 times faster than before on a third of them.This also fixes
m^nfor matrix objects without access to their rows, such as those inIsGenericMatrixRep, which failed withrow access unsupported.AI disclosure: benchmarking, analysis and drafting were done with the AI tool Claude Code; the commit carries a corresponding
Co-authored-by:line.