How to Learn the Paterson–Stockmeyer Algorithm: Baby Steps, Giant Steps, Matrix Polynomials and Multiplication-Minimising Evaluation

What if addition is cheap, multiplying by a scalar is cheap, but multiplying two matrices is extremely expensive? Horner’s method is elegant, but it does not minimise that expensive operation. Paterson–Stockmeyer reorganises a polynomial so that a degree-n expression can be evaluated with only O(√n) nonscalar multiplications.

This Learning Hall article develops the method from ordinary polynomial evaluation to baby steps, giant steps, matrix powers, multiplication counting, block-size choice, memory–compute trade-offs, numerical concerns and current mixed-precision work. The key lesson is not simply one formula: it is how changing the cost model can change the best algorithm.

Quick Read

We want to evaluate p(A)=c0+c1A+…+cnA^n when multiplying two A-like objects is expensive.
Horner’s method uses about n nonscalar multiplications.
Paterson–Stockmeyer chooses a block size s near √n.
Precompute the baby powers I,A,A²,…,A^(s−1).
Let B=A^s and rewrite the polynomial in blocks q0(A)+Bq1(A)+B²q2(A)+….
Evaluate the outer polynomial in B, usually by Horner’s rule.
The expensive multiplication count becomes roughly the cost of building baby powers plus the number of giant-step blocks: O(√n).
The method is especially important for matrix polynomials and matrix-function algorithms.
Professional implementations must trade multiplications against additions, storage, precision, coefficient structure and hardware.

1. Beginner Level: Start With Horner’s Rule

For a scalar polynomial:

p(x) = c0 + c1x + c2x² + c3x³

Horner rewrites it as:

p(x) = c0 + x(c1 + x(c2 + xc3))

This needs one multiplication by x for each degree. For ordinary floating-point scalars, that is often excellent. But when x is a large matrix A, every multiplication by A is a matrix–matrix multiplication and may dominate the whole computation.

2. The Cost Model Comes First

Paterson and Stockmeyer distinguish nonscalar multiplications from cheap scalar operations. In a matrix polynomial, multiplying a matrix by a scalar coefficient is much cheaper than multiplying two dense matrices of the same size.

So the optimisation target becomes:

minimise the number of matrix × matrix multiplications

That single change makes a different evaluation scheme rational.

3. Break the Exponents Into Blocks

Choose a positive block size s. Every exponent k can be written as:

k = js + r,    0 ≤ r < s

Therefore:

A^k = (A^s)^j A^r

This is the baby-step/giant-step idea. The small remainder powers A^r are baby steps. Repeated powers of A^s move through the polynomial in giant blocks.

4. Group the Polynomial

Define B=A^s. Group coefficients by their quotient j:

q0(A) = c0 + c1A + ... + c(s-1)A^(s-1)
q1(A) = cs + c(s+1)A + ...
q2(A) = c(2s) + ...

Then:

p(A) = q0(A) + B q1(A) + B² q2(A) + ...

Every qj uses only the already computed baby powers. No new matrix–matrix multiplication is needed to form a linear combination of those powers.

5. Build the Baby Powers

Compute:

A² = A·A
A³ = A²·A
...
A^s = A^(s-1)·A

This costs roughly s−1 nonscalar multiplications. The last result A^s becomes B, the giant-step base.

6. Evaluate the Giant-Step Polynomial

If there are m coefficient blocks, evaluate the outer polynomial in B using Horner:

R = q(m-1)
for j from m-2 down to 0:
    R = qj + B·R
return R

This adds about m−1 expensive multiplications. The total is therefore approximately:

(s - 1) + (m - 1)

with m about n/s. Balancing s and n/s gives s near √n.

7. Why √n Appears

The rough expensive-operation count is:

cost(s) ≈ s + n/s

If s is too small, there are too many giant blocks. If s is too large, baby-power construction is expensive. The two terms balance near s=√n, producing about 2√n expensive multiplications rather than n.

Paterson and Stockmeyer proved O(√n) nonscalar multiplication algorithms and a matching order-of-growth lower bound in their model. The square root is therefore structural, not an arbitrary tuning rule.

8. A Worked Degree-8 Example

Take degree 8 and choose s=3. Compute I,A,A² and B=A³. Group:

q0 = c0 I + c1 A + c2 A²
q1 = c3 I + c4 A + c5 A²
q2 = c6 I + c7 A + c8 A²

p(A) = q0 + B(q1 + B q2)

We spend two multiplications to reach A³ and then two giant multiplications by B: four matrix–matrix multiplications. Straight Horner would use eight.

9. High-Level Pseudocode

paterson_stockmeyer(A, c[0..n]):
    s = choose_block_size(n)

    powers[0] = I
    powers[1] = A
    for r = 2..s:
        powers[r] = powers[r-1] * A

    B = powers[s]

    blocks = []
    for each coefficient block j:
        q = zero_matrix
        for r = 0..s-1:
            k = j*s + r
            if k <= n:
                q += c[k] * powers[r]
        blocks.append(q)

    R = blocks[last]
    for j from last-1 down to 0:
        R = blocks[j] + B * R

    return R

A production implementation can avoid storing every qj simultaneously and can exploit zero coefficients, symmetry, sparsity, fused kernels or precision changes. The pseudocode isolates the mathematical decomposition first.

10. Why Matrix Polynomials Are the Natural Application

Many numerical algorithms approximate matrix functions—such as exp(A), log(A), cos(A) or related φ-functions—using polynomials or rational approximants. Once a polynomial approximation is chosen, evaluating it efficiently can dominate cost.

Recent numerical linear algebra still uses and extends Paterson–Stockmeyer. A 2025 SIAM paper developed a mixed-precision variant for matrix polynomials, showing that this 1970s cost-model idea remains relevant in modern high-precision computing.

11. Multiplication Count Is Not the Whole Runtime

Baby powers consume memory.
Forming many qj blocks requires additions and scalar multiplications.
Dense matrix multiplication may be accelerator-friendly, while many small operations may be bandwidth-bound.
Sparse matrices can become dense under powers.
Coefficient patterns may make some blocks cheap or zero.
Parallel hardware may prefer a different evaluation tree even if it uses slightly more multiplications.

The theoretical cost model identifies a valuable resource, but professional engineering measures the actual machine.

12. Numerical Accuracy and Stability

Algebraically equivalent evaluation schemes need not have identical floating-point behaviour. Rounding errors enter through matrix products and sums, and poorly scaled coefficients can amplify cancellation.

Therefore a production decision should consider backward error, conditioning of the target function, scaling, coefficient magnitudes and the arithmetic precision. Higham’s work on matrix functions and more recent mixed-precision Paterson–Stockmeyer research are useful guides for this layer.

13. Choosing the Block Size in Practice

√n is the first estimate, not always the final engineering choice. Try nearby integer values and count the exact number of required powers and giant blocks for the actual degree. If some highest coefficients vanish, the effective degree changes. If memory is tight, a smaller power table may be preferable even with an extra multiplication.

For fixed polynomial families used repeatedly, block size can be selected offline and benchmarked rather than recomputed on every call.

14. Implementation Failure Modes

Counting scalar multiplications as equally expensive: that erases the cost model the method exploits.
Off-by-one block boundaries: degree n means coefficients c0 through cn.
Forgetting I=A⁰: every block needs its constant coefficient.
Computing B twice: A^s should usually be the final baby power and giant base.
Using scalar-only tests: scalars commute with everything and may hide matrix implementation errors.
Assuming sparsity survives powers: fill-in can destroy memory advantages.
Ignoring overflow or scaling: large intermediate matrix norms can damage accuracy.
Picking s=floor(√n) without counting: adjacent block sizes may have different exact costs.

15. How to Test It

Compare against direct power-sum evaluation on small matrices.
Compare against Horner and verify numerical agreement within a justified tolerance.
Instrument matrix–matrix multiplication calls and confirm the predicted count.
Test diagonal matrices, random dense matrices, nilpotent matrices and matrices with large norms.
Use exact rational or symbolic matrices for tiny cases to separate algebraic correctness from floating-point error.
Benchmark several block sizes around √n.
Track peak memory as well as runtime.

16. Beginner-to-Professional Learning Ladder

Beginner: compare ordinary power-sum evaluation with Horner on scalar polynomials.
Intermediate: regroup a degree-8 or degree-15 polynomial by hand using s=3 or s=4.
Advanced: implement the method for matrices and instrument every expensive multiplication.
Professional: model arithmetic intensity, memory, sparsity, precision and stability; benchmark block sizes; and integrate the evaluator into a real matrix-function or approximation pipeline.

17. How to Learn the Decomposition

A useful sequence is Predict–Run–Investigate–Modify–Make. Predict how many multiplications Horner needs. Run a supplied degree-8 Paterson–Stockmeyer grouping. Investigate where every exponent went. Modify the block size and recount. Finally make the grouping automatically for arbitrary n. Subgoal labels—choose s, build powers, form blocks, giant-step Horner, verify count—reduce cognitive load without hiding the mathematics.

18. Practice Problems

Evaluate a degree-15 polynomial symbolically with block size 4 and count nonscalar multiplications.
For n=100, compare s=9,10 and11 using an exact block-count formula.
Implement Horner and Paterson–Stockmeyer with a multiplication counter rather than a timer.
Test both methods on a diagonal matrix and then on a random dense matrix.
Construct a sparse matrix whose square is much denser and discuss how that changes the practical cost.
Use a polynomial with many zero coefficients and redesign the evaluation plan.
Investigate whether a mixed-precision implementation can reduce cost while meeting a specified final error tolerance.

19. Sources and Further Reading

Final idea: Paterson–Stockmeyer is a masterclass in cost-aware reformulation. The polynomial has not changed; the algebraic answer has not changed. What changes is the route through the computation. Once matrix multiplication is identified as the scarce operation, exponents are reorganised into baby and giant steps so expensive work is reused rather than repeated.

Small Group Tutorials