<figure><img src="image_1.png" /></figure>

HAL
open science

# Designing and exploiting fast algorithms for univariate polynomial matrices

Vincent Neiger

<figure><img src="image_2.png" /></figure>

## To cite this version:

Vincent Neiger. Designing and exploiting fast algorithms for univariate polynomial matrices. École thématique. JNCF, CIRM, Luminy (Marseille), France. 2024. hal-04992937

HAL Id: hal-04992937

https://hal.science/hal-04992937v1

Submitted on 19 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
designing and exploiting fast algorithms
for univariate polynomial matrices

Journées Nationales de Calcul Formel
Centre International de Rencontre Mathématiques
Marseille Luminy, France, 4 March 2024
<figure><img src="image_4.png" /></figure>

computer algebra

<figure><img src="image_5.png" /></figure>

polynomial matrices

<figure><img src="image_6.png" /></figure>

first algorithms

<figure><img src="image_7.png" /></figure>

exercises
## outline

* computer algebra

* polynomial matrices

* first algorithms

* exercises

* efficient algorithms and software

* for matrices over a field

* for univariate polynomials
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

MAGMA
COMPUTER • ALGEBRA

<figure><img src="image_19.png" /></figure>

MATLAB®

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

Maple®

<figure><img src="image_22.png" /></figure>

WOLFRAM
MATHEMATICA 12

<figure><img src="image_23.png" /></figure>

sXge

<figure><img src="image_24.png" /></figure>

<figure><img src="image_25.png" /></figure>

SymPy

<figure><img src="image_26.png" /></figure>
computer algebra

algorithm design
and software implementations
for exact computations
with mathematical objects
Euclid's GCD -300

computer algebra

algorithm design
and software implementations
for exact computations
with mathematical objects
computer algebra

algorithm design
and software implementations
for exact computations
with mathematical objects
computer algebra

algorithm design
and software implementations
for exact computations
with mathematical objects

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>
computer algebra

algorithm design
and software implementations
for exact computations
with mathematical objects

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

SymPy

<figure><img src="image_32.png" /></figure>
computer algebra

algorithm design
and software implementations
for exact computations
with mathematical objects

<figure><img src="image_33.png" /></figure>

<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>
computer algebra

algorithm design
and software implementations
for exact computations
with mathematical objects

<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /></figure>

<figure><img src="image_41.png" /></figure>

<figure><img src="image_42.png" /></figure>

<figure><img src="image_43.png" /></figure>
computer algebra

algorithm design
and software implementations
for exact computations
with mathematical objects
Euclid's GCD -300

Gaussian elimination 179

Newton's method 1669

LLL '82, NFS '88

computer algebra

algorithm design
and software implementations
for exact computations
with mathematical objects

FFT 1805, '65

Buchberger '76

Strassen '69

Karatsuba '62
<table><caption>Principal Discoveries of Efficient Methods of Computing the DFT</caption><thead><tr><th>Researcher(s)</th><th>Date</th><th>Sequence Lengths</th><th>Number of DFT Values</th><th>Application</th></tr></thead><tbody><tr><td>C. F. Gauss [10]</td><td>1805</td><td>Any composite integer</td><td>All</td><td>Interpolation of orbits of celestial bodies</td></tr><tr><td>F. Carlini [28]</td><td>1828</td><td>12</td><td>—</td><td>Harmonic analysis of barometric pressure</td></tr><tr><td>A. Smith [25]</td><td>1846</td><td>4, 8, 16, 32</td><td>5 or 9</td><td>Correcting deviations in compasses on ships</td></tr><tr><td>J. D. Everett [23]</td><td>1860</td><td>12</td><td>5</td><td>Modeling underground temperature deviations</td></tr><tr><td>C. Runge [7]</td><td>1903</td><td>2<sup>8</sup></td><td>All</td><td>Harmonic analysis of functions</td></tr><tr><td>K. Stumpff [16]</td><td>1939</td><td>2<sup>8</sup>, 3<sup>8</sup></td><td>All</td><td>Harmonic analysis of functions</td></tr><tr><td>Danielson and Lanczos [5]</td><td>1942</td><td>2<sup>16</sup></td><td>All</td><td>X-ray diffraction in crystals</td></tr><tr><td>L. H. Thomas [13]</td><td>1948</td><td>Any integer with relatively prime factors</td><td>All</td><td>Harmonic analysis of functions</td></tr><tr><td>I. J. Good [3]</td><td>1958</td><td>Any integer with relatively prime factors</td><td>All</td><td>Harmonic analysis of functions</td></tr><tr><td>Cooley and Tukey [1]</td><td>1965</td><td>Any composite integer</td><td>All</td><td>Harmonic analysis of functions</td></tr><tr><td>S. Winograd [14]</td><td>1976</td><td>Any integer with relatively prime factors</td><td>All</td><td>Use of complexity theory for harmonic analysis</td></tr></tbody></table>

FFT 1805, '65

Karatsuba '62
<figure><img src="image_44.png" /></figure>
<figure><img src="image_45.png" /></figure>
<figure><img src="image_46.png" /></figure>
error correcting codes

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>

cryptographic protocols

<figure><img src="image_49.png" /></figure>

XXth-XXIst centuries: digital data & interconnected networks

**integrity – confidentiality**

discrete structures: *exact* and *intensive* computations

* matrices of large size, with sparsity or structure

* polynomials and polynomial matrices in one variable

* polynomials in several variables

goal of computer algebra

**fast algorithms: complexity & efficient implementations**
<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /></figure>

XXth-XXIst centuries: digital data & interconnected networks

**integrity – confidentiality**

discrete structures: *exact* and *intensive computations*

* matrices of large size, with sparsity or structure

* polynomials and polynomial matrices in one variable

* polynomials in several variables

goal of computer algebra

**fast algorithms: complexity & efficient implementations**

general methodology: **reductions** to efficient basic operations

* IntMul: integer multiplication

* MatMul: matrix multiplication

* PolMul: univariate polynomial multiplication
efficient algorithms for polynomials, matrices, power series, ...
with coefficients in some base field $K$

* low complexity bound

* low execution time

low memory usage, power consumption, ...

prime field $\mathbb{F}_p = \mathbb{Z}/p\mathbb{Z}$

field extension $\mathbb{F}_p[x]/\langle f(x) \rangle$

rationals $\mathbb{Q}$, number fields, ...
measuring efficiency

efficient algorithms for polynomials, matrices, power series, ...
with coefficients in some **base field** $\mathbb{K}$

► low complexity bound

► low execution time

low memory usage, power consumption, ...

prime field $\mathbb{F}_p = \mathbb{Z}/p\mathbb{Z}$
field extension $\mathbb{F}_p[x]/\langle f(x) \rangle$
rationals $\mathbb{Q}$, number fields, ...

algebraic complexity (upper) bounds
$\rightsquigarrow$ count number of operations in $\mathbb{K}$

<figure><img src="image_53.png" /></figure>

standard complexity model for algebraic computations

<figure><img src="image_54.png" /></figure>

often well correlated to implementation timings (e.g. over $\mathbb{K} = \mathbb{F}_p$)

<figure><img src="image_55.png" /></figure>

ignores coefficient growth (e.g. over $\mathbb{K} = \mathbb{Q}$)
measuring efficiency

efficient algorithms for polynomials, matrices, power series, ...
with coefficients in some base field $K$

▶ low complexity bound

▶ low execution time

low memory usage, power consumption, ...

prime field $F_p = \mathbb{Z}/p\mathbb{Z}$
field extension $F_p[x]/\langle f(x)\rangle$
rationals $\mathbb{Q}$, number fields, ...

practical performance

∼ measure software running time

⚠ strongly influenced by the
quality of the implementation

this talk:

▶ working over $K = F_p$ with word-size prime $p$

▶ Intel Core i7-7600U @ 2.80GHz, no multithreading
$$ M = \begin{bmatrix} 28 & 68 & 75 & 70 \\ 38 & 25 & 75 & 55 \\ 24 & 1 & 56 & 28 \end{bmatrix} \in \mathbb{K}^{3 \times 4} \rightarrow 3 \times 4 \text{ matrix over } \mathbb{K} \text{ (here } \mathbb{F}_{97}) $$

fundamental operations on $m \times m$ matrices:

* addition is “quadratic”: $O(m^2)$ operations in $\mathbb{K}$

* naive multiplication is cubic: $O(m^3)$

[Strassen'69]

breakthrough: **subcubic matrix multiplication**
$$ M = \begin{bmatrix} 28 & 68 & 75 & 70 \\ 38 & 25 & 75 & 55 \\ 24 & 1 & 56 & 28 \end{bmatrix} \in \mathbb{K}^{3 \times 4} \longrightarrow 3 \times 4 \text{ matrix over } \mathbb{K} \text{ (here } \mathbb{F}_{97} \text{)} $$

fundamental operations on $m \times m$ matrices:

► addition is “quadratic”: $O(m^2)$ operations in $\mathbb{K}$

► naive multiplication is cubic: $O(m^3)$

[Strassen'69]

breakthrough: **subcubic matrix multiplication**

► complexity exponent $\omega \approx 2.81$ — i.e. $O(m^\omega)$ complexity

► used in practice for $m \ge$ a few 100s
in NTL, FLINT, fflas-ffpack...

[Coppersmith-Winograd 1990]

► best-known exponent $\omega \approx 2.3719$

[Le Gall'14] [Alman-Williams'20] [Duan-Wu-Zhou'23]

► “galactic” algorithms: strongly impractical as such
<figure><img src="image_56.png" /></figure>

Strassen, in his seminal 1969 paper

“Gaussian Elimination is Not Optimal”

sent a clear message to the scientific community:

Natural, obvious and centuries-old
methods for solving important
computational problems may be
far from the fastest.
<figure><img src="image_57.png" /></figure>

Strassen, in his seminal 1969 paper

"Gaussian Elimination is Not Optimal"

sent a clear message to the scientific
community:

Natural, obvious and centuries-old
methods for solving important
computational problems may be
far from the fastest.

<figure><img src="image_58.png" /></figure>

<figure><img src="image_59.png" /></figure>

<figure><img src="image_60.png" /></figure>

1:21 / 44:52 • Gaussian elimination is not optimal >

<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /></figure>

CC

<figure><img src="image_63.png" /></figure>

<figure><img src="image_64.png" /></figure>

<figure><img src="image_65.png" /></figure>

Abel Prize Laureates - Jeroen Zuiddam - Matrix multiplication and Shannon capacity

<figure><img src="image_66.png" /></figure>

Centrum Wiskunde & Informatica
1.16K subscribers

Subscribe

<figure><img src="image_67.png" /></figure>

<figure><img src="image_68.png" /></figure>

<figure><img src="image_69.png" /></figure>

<figure><img src="image_70.png" /></figure>

Share

<figure><img src="image_71.png" /></figure>

189 views 1 year ago Abel Prize Laureates Lectures (8 April 2022)

Asymptotic spectrum duality in computer science and discrete mathematics: Matrix multiplication and Shannon capacity

Abstracti ...more

0 Comments

Sort by
<figure><img src="image_72.png" /></figure>
reductions: a new hope

<figure><img src="image_73.png" /></figure>

<figure><img src="image_74.png" /></figure>
reductions: a new hope

<figure><img src="image_75.png" /></figure>

<figure><img src="image_76.png" /></figure>

take-home messages:

▸ bibliometric indicators measure quantity, and there exist counterexamples to "quantity = quality"

▸ design fast algorithms for the most basic routines

→ MatMul

▸ design efficient reductions to them for other tasks

→ LinSys, Det, Inverse
polynomials: multiplication

$$p = 87x^7 + 74x^6 + 60x^5 + 46x^4 + 16x^3 + 41x^2 + 86x + 69$$

$p \in K[x]_{<8} \rightarrow \text{univariate polynomial in } x \text{ of degree } < 8 \text{ over } K$

fundamental operations on polynomials of degree < d:

* addition and Horner's evaluation are linear: $O(d)$

* naive multiplication is quadratic: $O(d^2)$

$$[\text{Karatsuba'62}] \quad M(d) \in O(d^{1.58})$$

breakthrough: **subquadratic polynomial multiplication**
polynomials: multiplication

$$p = 87x^7 + 74x^6 + 60x^5 + 46x^4 + 16x^3 + 41x^2 + 86x + 69$$

$p \in \mathbb{K}[x]_{<8} \rightarrow \text{univariate polynomial in } x \text{ of degree } < 8 \text{ over } \mathbb{K}$

fundamental operations on polynomials of degree < d:

* addition and Horner's evaluation are linear: $O(d)$

* naive multiplication is quadratic: $O(d^2)$

[Karatsuba'62] $M(d) \in O(d^{1.58})$

breakthrough: subquadratic polynomial multiplication

[Schönhage-Strassen'71] [Nussbaumer'80] [Cantor-Kaltofen'91] $M(d) \in O(d \log(d) \log \log(d))$

breakthrough: quasi-linear polynomial multiplication

research still active, with recent progress by [Harvey-van der Hoeven-Lecerf]

* change of representation by evaluation-interpolation

* used in practice as soon as $d \approx 100$

* FFT techniques using (virtual) roots of unity

note: $M(d) \in O(d \log(d))$
if provided a "good" root of unity
```
3250 void DivRem(zz_pX& q, zz_pX& r, const zz_pX& a,
3251 +--- 6 lines: {-------}
3257
3258 void div(zz_pX& q, const zz_pX& a, const zz_pX& a,
3259 +--- 6 lines: {-------}
3265
3266 void div(zz_pX& q, const zz_pX& a, zz_p b)
3267 +--- 5 lines: {-------}
3272
3273
3274 void rem(zz_pX& r, const zz_pX& a, const zz_pX& a,
3275 +--- 6 lines: {-------}
3281
3282 []
3283
3284 long operator== (const zz_pX& a, long b)
3285 +--- 20 lines: {-------}
3305
3306 long operator== (const zz_pX& a, zz_p b)
3307 +--- 11 lines: {-------}
3318
3319 void power(zz_pX& x, const zz_pX& a, long e)
3320 +--- 40 lines: {-------}
3360
3361 void reverse(zz_pX& x, const zz_pX& a, long hi)
3362 +--- 13 lines: {-------}
3375
3376 NTL_END_IMPL
```

NORMAL > main > software/ntl/src/lzz_pX.cpp

► small prime FFT in NTL:

≈ about **5500** lines of C++

≈ target operation: FFT

(including 1200 lines for vectorized version and 1100 for machine word arithmetic...)

► polynomials in $Z/pZ[x]$:

≈ about **5500** lines as well

≈ target operations include:

► multiplication, truncated inversion, division,

► interpolation, multipoint evaluation,

► XGCD, Berlekamp-Massey, resultant,

► power projection, modular composition, ...
3165 void FFTDiv(zz_pX& q, const zz_pX& a, const zz_pX& b)
3166 {
3167
3168 long n = deg(b);
3169 long m = deg(a);
3170 long k;
3171
3172 +---- 4 lines: if (m < n) {
3173
3174 +---- 6 lines: if (m >= 3*n) {
3183
3184 zz_pX P1, P2, P3;
3185
3186 CopyReverse(P3, b, 0, n);
3187 InvTrunc(P2, P3, m-n+1);
3188 CopyReverse(P1, P2, 0, m-n);
3189
3190 k = NextPowerOfTwo(2*(m-n)+1);
3191
3192 fftRep R1(INIT_SIZE, k), R2(INIT_SIZE, k);
3193
3194
3195
3196
3197
3198 }

small prime FFT in NTL:
about 5500 lines of C++
target operation: FFT
(including 1200 lines for vectorized version
and 1100 for machine word arithmetic...)

► small prime FFT in NTL:

about 5500 lines of C++

target operation: FFT

► polynomials in Z/pZ[x]:

about 5500 lines as well

target operations include:

multiplication, truncated inversion, division,

interpolation, multipoint evaluation,

XGCD, Berlekamp-Massey, resultant,

power projection, modular composition, ...

► reductions are often

concise and readable

close to the pseudocode
concentrate efforts on: basic routines + good reductions

3165 void FFTDiv(zp_x& q, const zz_pX& a, const zz_pX& b)
3166 {
3167
3168 long n = deg(b);
3169 long m = deg(a);
3170 long k;
3171
3172 +---- 4 lines: if (m < n) {
3173
3174 +---- 6 lines: if (m >= 3*n) {
3175
3176
3177
3178
3179
3180
3181
3182
3183
3184
3185
3186
3187
3188
3189
3190
3191
3192
3193
3194
3195
3196
3197
3198
}

zz_pX P1, P2, P3;
CopyReverse(P3, b, 0, n);
InvTrunc(P2, P3, m-n+1);
CopyReverse(P1, P2, 0, m-n);
k = NextPowerOfTwo(2*(m-n)+1);
fftRep R1(INIT_SIZE, k), R2(INIT_SIZE, k);

TofftRep(R1, P1, k);
TofftRep(R2, a, k, n, m);
mul(R1, R1, R2);
FromfftRep(q, R1, m-n, 2*(m-n));

▶ m ← deg(A) and n ← deg(B)

▶ if m < n, return (0, A)

▶ set reversals $\tilde{A} \leftarrow x^m A(1/x)$ and $\tilde{B} \leftarrow x^n B(1/x)$

▶ find $\tilde{Q} = \tilde{A}/\tilde{B}$ mod $x^{m-n+1}$ by power series inversion and product

▶ reverse $\tilde{Q}$ to obtain $\tilde{Q}$
reductions of most problems to matrix multiplication

<figure><img src="image_77.png" /></figure>

PLUQ = Gaussian elimination

TRSM = triangular solving

$$ \left. \begin{array}{l} \text{LinSys} \\ \text{Det} \\ \text{Rank} \\ \text{PLUQ} \\ \text{TRSM} \\ \text{Inverse} \end{array} \right\} = O(\text{MatMul}) $$

not closed:

open:
reductions of most problems to matrix multiplication

<figure><img src="image_78.png" /></figure>

PLUQ = Gaussian elimination

TRSM = triangular solving

$$ \left. \begin{array}{l} \text{LinSys} \\ \text{Det} \\ \text{Rank} \\ \text{PLUQ} \\ \text{TRSM} \\ \text{Inverse} \end{array} \right\} = O(\text{MatMul}) $$

$$ \left. \begin{array}{l} \text{MinPoly} \\ \text{CharPoly} \end{array} \right\} = O(\text{MatMul}) $$

exploiting non-naive PolMul

not closed:

open:
reductions of most problems to matrix multiplication

<figure><img src="image_79.png" /></figure>

PLUQ = Gaussian elimination

TRSM = triangular solving

$$ \left. \begin{array}{l} \text{LinSys} \\ \text{Det} \\ \text{Rank} \\ \text{PLUQ} \\ \text{TRSM} \\ \text{Inverse} \end{array} \right\} = O(\text{MatMul}) $$

$$ \left. \begin{array}{l} \text{MinPoly} \\ \text{CharPoly} \end{array} \right\} = O(\text{MatMul}) $$

exploiting non-naive PolMul

not closed: is Frobenius normal form in O(MatMul)?
open:
reductions of most problems to matrix multiplication

<figure><img src="image_80.png" /></figure>

PLUQ = Gaussian elimination

TRSM = triangular solving

$$ \left. \begin{array}{l} \text{LinSys} \\ \text{Det} \\ \text{Rank} \\ \text{PLUQ} \\ \text{TRSM} \\ \text{Inverse} \end{array} \right\} = O(\text{MatMul}) $$

$$ \left. \begin{array}{l} \text{MinPoly} \\ \text{CharPoly} \end{array} \right\} = O(\text{MatMul}) $$

exploiting non-naive PolMul

not closed: is Frobenius normal form in $O(\text{MatMul})$?
open:
reductions of most problems to matrix multiplication

<figure><img src="image_81.png" /></figure>

PLUQ = Gaussian elimination

TRSM = triangular solving

$$
\left.
\begin{array}{l}
\text{LinSys} \\
\text{Det} \\
\text{Rank} \\
\text{PLUQ} \\
\text{TRSM} \\
\text{Inverse}
\end{array}
\right\} = O(\text{MatMul})
$$

$$
\left.
\begin{array}{l}
\text{MinPoly} \\
\text{CharPoly}
\end{array}
\right\} = O(\text{MatMul})
$$

exploiting non-naive PolMul

not closed: is Frobenius normal form in O(MatMul)?

open: is linear system solving as hard as multiplication?
bonus: some notes/references

[Jeannerod-Pernet-Storjohann 2013] doi.org/10.1016/j.jsc.2013.04.004

► explicit reductions between inversion & MatMul & Gaussian elimination / echelonization

► constants in the $O(\cdot)$ complexities when using classical matrix multiplication ($\omega = 3$) or Strassen's multiplication

“not closed”: it is open, but

► there is a randomized algorithm for Frobenius form computation which has complexity $O(\mathrm{MatMul})$

[Pernet-Storjohann 2007] http://www.cs.uwaterloo.ca/~astorjoh/cpoly.pdf

► recent developments give new insight concerning core operations typically used in Frobenius form algorithms

charpoly in $O(\mathrm{MatMul})$: [Neiger-Pernet 2021] doi.org/10.1016/S0885-064X(22)00005-X
Krylov iterates in $O(\mathrm{MatMul})$: [Neiger-Pernet-Villard 2024] hal.science/hal-04445355
## most problems have quasi-linear complexity

thanks to reductions to PolMul — did we mention the importance of good reductions?

* addition $f + g$, multiplication $f * g$

* division with remainder $f = qg + r$

* truncated inverse $f^{-1} \mod x^d$

* extended GCD $fu + gv = gcd(f, g)$

* multipoint eval. $f \mapsto f(\alpha_1), \dots, f(\alpha_d)$

* interpolation $f(\alpha_1), \dots, f(\alpha_d) \mapsto f$

* Padé approximation $f = [p]_q \mod x^d$

* minpoly of linearly recurrent sequence

<figure><img src="image_82.png" /></figure>

open:

open:

open:

open:
## most problems have quasi-linear complexity

thanks to reductions to PolMul — did we mention the importance of good reductions?

### $O(M(d))$

* addition $f + g$, multiplication $f * g$

* division with remainder $f = qg + r$

* truncated inverse $f^{-1} \mod x^d$

* extended GCD $fu + gv = gcd(f, g)$

* multipoint eval. $f \mapsto f(\alpha_1), \dots, f(\alpha_d)$

* interpolation $f(\alpha_1), \dots, f(\alpha_d) \mapsto f$

* Padé approximation $f = [p]_q \mod x^d$

* minpoly of linearly recurrent sequence

<figure><img src="image_83.png" /></figure>

open: polynomial multiplication in $O(d \log(d))$?

open:

open:

open:
# polynomials: main computational problems

## most problems have quasi-linear complexity

thanks to reductions to PolMul — did we mention the importance of good reductions?

$O(M(d))$

* addition $f + g$, multiplication $f * g$

* division with remainder $f = qg + r$

* truncated inverse $f^{-1} \mod x^d$

* extended GCD $fu + gv = gcd(f, g)$

$O(M(d) \log(d))$

* multipoint eval. $f \mapsto f(\alpha_1), \dots, f(\alpha_d)$

* interpolation $f(\alpha_1), \dots, f(\alpha_d) \mapsto f$

* Padé approximation $f = [p]_q \mod x^d$

* minpoly of linearly recurrent sequence

<figure><img src="image_84.png" /></figure>

* open: polynomial multiplication in $O(d \log(d))$?

* open: interpolation and multipoint eval. in $O(\text{PolMul})$?

* open: XGCD and friends in $O(\text{PolMul})$?

* open:
# polynomials: main computational problems

## most problems have quasi-linear complexity

thanks to reductions to PolMul — did we mention the importance of good reductions?

O(M(d))

* addition $f + g$, multiplication $f * g$

* division with remainder $f = qg + r$

* truncated inverse $f^{-1} \mod x^d$

* extended GCD $fu + gv = gcd(f, g)$

O(M(d) log(d))

* multipoint eval. $f \mapsto f(\alpha_1), \dots, f(\alpha_d)$

* interpolation $f(\alpha_1), \dots, f(\alpha_d) \mapsto f$

* Padé approximation $f = [p]_q \mod x^d$

* minpoly of linearly recurrent sequence

<figure><img src="image_85.png" /></figure>

<figure><img src="image_86.png" /></figure>

<figure><img src="image_87.png" /></figure>

<figure><img src="image_88.png" /></figure>

* open: polynomial multiplication in O(d log(d))?

* open: interpolation and multipoint eval. in O(PolMul)?

* open: XGCD and friends in O(PolMul)?

* open: modular composition $f(g) \mod h$ closer to O(PolMul)?
bonus: some notes/references

polynomial multiplication in O(d log(d))?

► remains open over an arbitrary field, concerning algebraic complexity

► solved when the field possesses suitable roots of unity for FFT

► method of choice in practice (using several primes and CRT if needed) when working over prime finite fields $Z/pZ$

► recent progress in the bit complexity model
[Harvey-van der Hoeven 2019] https://doi.org/10.1016/j.jco.2019.03.004
[Harvey-van der Hoeven 2022] https://doi.org/10.1145/3505584

interpolation and multipoint evaluation in O(PolMul)?

► remains open for an arbitrary set of points, with no assumption, but:

► by design, solved for FFT points, when available

► more generally, solved for points forming a geometric sequence
[Bostan-Schost 2005] https://doi.org/10.1016/j.jco.2004.09.009

► in many applications of interpolation/evaluation, one can choose the points, in which case $O(\text{PolMul})$ is feasible
<figure><img src="image_89.png" /></figure>

high-performance exact linear algebra

LinBox – fflas-ffpack

C/C++

high-performance polynomials (and more)

NTL & FLINT

C/C++

```c
// order that remains to be dealt with
VecLong ren_order();

// Indices of columns/orders that remain to be dealt with
VecLong ren_index(cdrn);

// std::ita(ren_index.begin(), ren_index.end(), 0);

// all along the algorithm, shift shifted row degrees of approximate matrix
// (Initially, input shift = shifted row degree of the identity matrix)

while (not ren_order.empty())
{
    // ** Invariant:
    // - appbas is a shift-orderd weak Popov approxontant basis for
    //   (pnat,reached_order) where doneorder is the tuple such that
    //   -->reached_order[i] = ren_order[i] = order[i] for j appearing in
    //   residual
    // - shift = the "input shift"-row degree of appbas
    // - residual = submatrx of columns (appbas + pnat)[i::] for all
```

matrices

software

polynomials

```c
sage: M = GF(7);
sage: A = Matrix(M, 3, [3, 2, 1, 2x, 1, 2x, 1, 2x, 2])
sage: A = hermite form
1
0
sage: A hermite form(transformation=True)
[ 1   1   2x   [1,0]
  0   0   7x+2, 2] [6,1]
sage: A = matrix(M, 3, [3, 2, 1, 2x, 2, 4x])
sage: A = hermite form(transformation=True, include_zero_rows=False)
[ 1   1   2x   [0,4]
  0   0   8, 5] [5,1]
sage: U = A hermite form(transformation=True, include_zero_rows=True); U
[ 1   1   2x   [0,4]
  0   0   8, 5] [5,1]
sage: U = A^A = M
True
sage: U, U = A hermite form(transformation=True, include_zero_rows=False)
sage: U = A
[ 1   1   2x   [0,4]
  0   0   8, 5] [5,1]
sage: U = A^A = M
True
```

See also: `is_hermite()`.

`is_hermite(row_wise=True, lower_echelon=False, include_zeroTrue) = True`

Return a boolean indicating whether this is the matrix is in Hermite form.

```c
std::distance(ren_order.begin(), std::maxElement(ren_order).end());
long deg = order[ren_index[i]] - ren_order[i];

// record the coefficients of degree deg of the column j of residual
// also keep track of which of these are nonzero,
// and among the nonzero ones, which is the first with smallest shift
VecZZP const_residual;
const residual[i] = order[i] for j appearing in ren_index
VecLong indices_nonzero;
long plv = -1;

for (long i = 0; i < rdn; ++i)
{
    // const_residual[i] = coeff(residual[i][], deg);
    // (const_residual[i] == 0)
    {
        // indices_nonzero.push_back(i);
        // (plv[i] << shift[i] & shift[plv])
        plv = i;
    }

    // (if indices_nonzero is empty, const_residual is already zero, there
    // is not indices_nonzero.empty())
    // ... update all rows of appbas and residual in indices nonzero rows
```

`src/nat_t2pApproximate.cpp`
open-source mathematics software system

<figure><img src="image_90.png" /></figure>

Python/Cython

high-performance exact linear algebra

LinBox – fflas-ffpack C/C++

high-performance polynomials (and more)

NTL & FLINT C/C++

▶ choice of algorithms

▶ data structures and storage

▶ cache efficiency

▶ SIMD vectorization instructions

▶ multithreading, GPU programming

matrices

software

polynomials

what you can compute in about 1 second
with fflas-ffpack
with NTL

▶ PLUQ m = 3800 1.00s

▶ LinSys m = 3800 1.00s

▶ MatMul m = 3000 0.97s

▶ Inverse m = 2800 1.01s

▶ CharPoly m = 2000 1.09s

is_hermiterow_wise=True, lower_echelon=False, include_zero vectors=True)

Return a boolean indicating whether this matrix is in Hermitian form.
open-source mathematics software system

<figure><img src="image_91.png" /></figure>

Python/Cython

high-performance exact linear algebra

LinBox – fflas-ffpack C/C++

high-performance polynomials (and more)

NTL & FLINT C/C++

▶ choice of algorithms

▶ data structures and storage

▶ cache efficiency

▶ SIMD vectorization instructions

▶ multithreading, GPU programming

matrices

software

polynomials

what you can compute in about 1 second
with fflas-ffpack with NTL

▶ PLUQ m = 3800 1.00s

▶ LinSys m = 3800 1.00s

▶ MatMul m = 3000 0.97s

▶ Inverse m = 2800 1.01s

▶ CharPoly m = 2000 1.09s

▶ PolMul d = 7 × 10⁶ 1.03s

▶ Division d = 4 × 10⁶ 0.96s

▶ XGCD d = 2 × 10⁵ 0.99s

▶ MinPoly d = 2 × 10⁵ 1.10s

▶ MPEval d = 1 × 10⁴ 1.01s

is_hermiterow_wise=True, lower_echelon=False, include_zeroTrue

Return a boolean indicating whether this matrix is in Hermite form.
open-source mathematics software system

<figure><img src="image_92.png" /></figure>

Python/Cython

high-performance exact linear algebra

LinBox – fflas-ffpack C/C++

high-performance polynomials (and more)

NTL & FLINT C/C++

▶ choice of algorithms

▶ data structures and storage

▶ cache efficiency

▶ SIMD vectorization instructions

▶ multithreading, GPU programming

matrices

software

polynomials

what you can compute in about 1 second
with fflas-ffpack with NTL

▶ PLUQ m = 3800 1.00s

▶ LinSys m = 3800 1.00s

▶ MatMul m = 3000 0.97s

▶ Inverse m = 2800 1.01s

▶ CharPoly m = 2000 1.09s

▶ PolMul d = 7 × 10⁶ 1.03s

▶ Division d = 4 × 10⁶ 0.96s

▶ XGCD d = 2 × 10⁵ 0.99s

▶ MinPoly d = 2 × 10⁵ 1.10s

▶ MPEval d = 1 × 10⁴ 1.01s

is_hermiterow_wise=True, lower_echelon=False, include_zero vectors=True

Return a boolean indicating whether this matrix is in Hermite form.
matrix exponentiation

input: matrix $A \in \mathbb{K}^{m \times m}$,
integer $k > 0$

output: $A^k$
matrix exponentiation

input: matrix $A \in \mathbb{K}^{m \times m}$,
integer $k > 0$

output: $A^k$

▶ repeated squaring: $O(m^\omega \log(k))$
matrix exponentiation

input: matrix $A \in \mathbb{K}^{m \times m}$,
integer $k > 0$

output: $A^k$

▶ repeated squaring: $O(m^\omega \log(k))$

▶ using Frobenius form:
$O(m^\omega \log(m) \log \log(m))$ if $\log(k) \in O(m)$
[Giesbrecht 1995] [Storjohann 2001]

▶ improvement with polynomial matrices:
$O(m^\omega \log \log(m)^2)$ if $\log(k) \in O(m)$
[Giesbrecht 1995] [Neiger-Pernet-Villard 2024]
matrix exponentiation

input: matrix $A \in \mathbb{K}^{m \times m}$,
integer $k > 0$

output: $A^k$

can we reach $O(m^\omega)$?

▶ repeated squaring: $O(m^\omega \log(k))$

▶ using Frobenius form:
$O(m^\omega \log(m) \log \log(m))$ if $\log(k) \in O(m)$
[Giesbrecht 1995] [Storjohann 2001]

▶ improvement with polynomial matrices:
$O(m^\omega \log \log(m)^2)$ if $\log(k) \in O(m)$
[Giesbrecht 1995] [Neiger-Pernet-Villard 2024]
constant matrices accelerated by polynomial matrices

matrix exponentiation

input: matrix $A \in \mathbb{K}^{m \times m}$,
integer $k > 0$

output: $A^k$

can we reach $O(m^\omega)$?

▶ repeated squaring: $O(m^\omega \log(k))$

▶ using Frobenius form:
$O(m^\omega \log(m) \log \log(m))$ if $\log(k) \in O(m)$
[Giesbrecht 1995] [Storjohann 2001]

▶ improvement with polynomial matrices:
$O(m^\omega \log \log(m)^2)$ if $\log(k) \in O(m)$
[Giesbrecht 1995] [Neiger-Pernet-Villard 2024]

Krylov iterates

input: matrix $A \in \mathbb{K}^{m \times m}$,
vector $v \in \mathbb{K}^{m \times 1}$

output: $v, Av, \dots, A^{m-1}v$
constant matrices accelerated by polynomial matrices

matrix exponentiation

input: matrix $A \in \mathbb{K}^{m \times m}$,
integer $k > 0$

output: $A^k$

can we reach $O(m^\omega)$?

▶ repeated squaring: $O(m^\omega \log(k))$

▶ using Frobenius form:
$O(m^\omega \log(m) \log \log(m))$ if $\log(k) \in O(m)$
[Giesbrecht 1995] [Storjohann 2001]

▶ improvement with polynomial matrices:
$O(m^\omega \log \log(m)^2)$ if $\log(k) \in O(m)$
[Giesbrecht 1995] [Neiger-Pernet-Villard 2024]

Krylov iterates

input: matrix $A \in \mathbb{K}^{m \times m}$,
vector $v \in \mathbb{K}^{m \times 1}$

output: $v, Av, \dots, A^{m-1}v$

▶ repeated matrix-vector products: $O(m^3)$

▶ via repeated squaring: $O(m^\omega \log(m))$ [Keller-Gehrig 1985]
constant matrices accelerated by polynomial matrices

matrix exponentiation

input: matrix $A \in \mathbb{K}^{m \times m}$,
integer $k > 0$

output: $A^k$

can we reach $O(m^\omega)$?

▶ repeated squaring: $O(m^\omega \log(k))$

▶ using Frobenius form:
$O(m^\omega \log(m) \log \log(m))$ if $\log(k) \in O(m)$
[Giesbrecht 1995] [Storjohann 2001]

▶ improvement with polynomial matrices:
$O(m^\omega \log \log(m)^2)$ if $\log(k) \in O(m)$
[Giesbrecht 1995] [Neiger-Pernet-Villard 2024]

Krylov iterates

input: matrix $A \in \mathbb{K}^{m \times m}$,
vector $v \in \mathbb{K}^{m \times 1}$

output: $v, Av, \dots, A^{m-1}v$

▶ repeated matrix-vector products: $O(m^3)$

▶ via repeated squaring: $O(m^\omega \log(m))$
[Keller-Gehrig 1985]

▶ with polynomial matrices: $O(m^\omega)$
[Zhou-Labahn-Storjohann 2012] [Neiger-Pernet 2021]
constant matrices accelerated by polynomial matrices

matrix exponentiation

input: matrix $A \in \mathbb{K}^{m \times m}$,
integer $k > 0$

output: $A^k$

can we reach $O(m^\omega)$?

▶ repeated squaring: $O(m^\omega \log(k))$

▶ using Frobenius form:
$O(m^\omega \log(m) \log \log(m))$ if $\log(k) \in O(m)$
[Giesbrecht 1995] [Storjohann 2001]

▶ improvement with polynomial matrices:
$O(m^\omega \log \log(m)^2)$ if $\log(k) \in O(m)$
[Giesbrecht 1995] [Neiger-Pernet-Villard 2024]

Krylov iterates

input: matrix $A \in \mathbb{K}^{m \times m}$,
vector $v \in \mathbb{K}^{m \times 1}$

output: $v, Av, \dots, A^{m-1}v$

we do reach $O(m^\omega)!$

... can we do better?

▶ repeated matrix-vector products: $O(m^3)$

▶ via repeated squaring: $O(m^\omega \log(m))$
[Keller-Gehrig 1985]

▶ with polynomial matrices: $O(m^\omega)$
[Zhou-Labahn-Storjohann 2012] [Neiger-Pernet 2021]
open-source mathematics software system

<figure><img src="image_93.png" /></figure>

Python / Cython

The Heaviside form is, e.g., the p(x) polynomial is analytic.

INPUT:

**goals: complete, robust, available**
(more than 60k downloads per month)

OUTPUT:

```c
// order that remains to be dealt with
VecLong rem_order;

// Indices of columns/orders that remain to be dealt with
VecLong rem_index[cdn];

// std::ita(rem_index.begin(), rem_index.end(), 0);

// all along the algorithm, shift - shifted row degrees of approximate matrix
// (Initially, input shift = shifted row degree of the identity matrix)
while (not rem_order.empty())
{
    // ** Invariant:
    // - appbas is a shift-orderd weak Popov approxontant basis for
    //   (pnat, reached_order) where doneorder is the tuple such that
    //   -->reached_order[i] = rem_order[i] = order[i] for j appearing in
    //   -->reached_order[i] = order[i] for j not appearing in rem_index
    // - shift == the "input shift" - row degree of appbas
    // - residual = submatrx of columns (appbas + pnat)[i::] for all
```

software development for polynomial matrices

```
sage: M = GF(7)
sage: A = matrix(M, 2, 3, [x, 1, 2*x, x; 1, 2, 2])
sage: A.is_hermitic form
True
0  1  2*x  [1:0]
0  0  7*x+2  [6:1]
sage: A.is_hermitic form(transformation=True)
[ 1  1  2*x  [1:0]
 0  0  7*x+2  [6:1]
True
```

```
sage: A = matrix(M, 2, 3, [x, 1, 2*x, 2, 2, 4*x])
sage: A.is_hermitic form(transformation=True, include_zero_rows=False)
[ 1  1  2*x  [0:4]
 0  0  7  [5:1]
True
```

```
sage: H, U = A.is_hermitic form(transformation=True, include_zero_rows=True); U, H
[ 1  1  2*x  [0:4]
 0  0  7  [5:1]
True
```

```
sage: U, A == H
True
```

```
sage: H, U, A = A.is_hermitic form(transformation=True, include_zero_rows=False);
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

```
sage: U, A == H
True
```

See also: `is_hermitic()`.

`is_hermitic(row_wise=True, lower_echelon=False, include_zeros=True)`

Return a boolean indicating whether this matrix is in Hermite form.

```c
std::distance(rem_order.begin(), std::maxElement(rem_order).begin());
long long deg = order[rem_index[i]] - rem_index[i];
// record the coefficients of degree deg of the column j of residual
// also keep track of which of these are nonzero,
// and among the nonzero ones, which is the first with smallest shift
VecZZP const_residual;
const residual, setLength(residual);
VecLong indices_nonzero;
long plv = 1;
for (long i = 0; i < rdn; ++i)
{
    // const_residual[i] = coeff(residual[i], deg);
    // (const_residual[i] != 0)
    {
        // indices_nonzero.push_back(i);
        // (plv[0] << shift[i] & shift[plv])
        plv = i;
    }
}
// (if indices_nonzero is empty, const_residual is already zero, there
// is not indices_nonzero.empty())
// update all rows of appbas and residual in indices_nonzero res
```

src/nat_t2z pk_aproximate.cpp

17
open-source mathematics software system

<figure><img src="image_94.png" /></figure>

Python/Cython

**goals: complete, robust, available**

(more than 60k downloads per month)

high-performance exact linear algebra

**LinBox – fflas-ffpack** C/C++

**goal: optimized basic operations**

algorithms, vectorization, multithreading

# software development for polynomial matrices

## Polynomial Matrix Library C/C++

533 files, 72k lines of code, including 21k lines of comments

https://github.com/vneiger/pml

[Hyun-Neiger-Schost '19]

* most tools are based on NTL

* work-in-progress version based on FLINT

* welcome comments, suggestions, contributions
  "hey, this doesn't work!"

"yo, plans for implementing this?"

"how to compute this determinant with PML?"

wide range of algorithms

efficiency = state of the art

system solving, determinant, kernel,
Hermite-Padé approximation,
high-order lifting, basis reduction...

is_hermite:row_wise=True, lower_echelon=False, include_zero

Return a boolean indicating whether this matrix is in Hermite form.

C/C++ library for computing approximate eigenvalues and
computer algebra

polynomial matrices

first algorithms

exercises

efficient algorithms and software

for matrices over a field

for univariate polynomials
computer algebra

efficient algorithms and software

for matrices over a field

for univariate polynomials

polynomial matrices

basic definitions and properties

use in various situations

seen as matrices / seen as polynomials

first algorithms

exercises
$\mathbb{K}[x]^{m \times n}$ is set of $m \times n$ matrices over $\mathbb{K}[x]$

called polynomial matrices in what follows

$$ \begin{bmatrix} 3x+4 & x^3+4x+1 & 4x^2+3 \\ 5 & 5x^2+3x+1 & 5x+3 \\ 3x^3+x^2+5x+3 & 6x+5 & 2x+1 \end{bmatrix} \in \mathbb{K}[x]^{3 \times 3} $$

► structure: matrices over $\mathbb{K}[x] \leftrightarrow$ free modules over $\mathbb{K}[x]$

similarly to: matrices over $\mathbb{K} \leftrightarrow$ vector spaces over $\mathbb{K}$

► basic operations: addition and multiplication
defined as usual (multiplication requires compatible dimensions)

► $\mathbb{K}[x]$ is not a field
$\mathbb{K}[x]^{m \times n} = \text{set of } m \times n \text{ matrices over } \mathbb{K}[x]$

called polynomial matrices in what follows

$$ \begin{bmatrix} 3x + 4 & x^3 + 4x + 1 & 4x^2 + 3 \\ 5 & 5x^2 + 3x + 1 & 5x + 3 \\ 3x^3 + x^2 + 5x + 3 & 6x + 5 & 2x + 1 \end{bmatrix} \in \mathbb{K}[x]^{3 \times 3} $$

► structure: matrices over $\mathbb{K}[x] \leftrightarrow$ free modules over $\mathbb{K}[x]$

similarly to: matrices over $\mathbb{K} \leftrightarrow$ vector spaces over $\mathbb{K}$

► basic operations: addition and multiplication
defined as usual (multiplication requires compatible dimensions)

► $\mathbb{K}[x]$ is not a field

↔ algorithms may work in $\mathbb{K}(x)^{m \times n}$, but be careful with “degree explosion”!
**large matrices with small degrees:**

characteristic polynomial $\det(x\mathbf{I}_m - \mathbf{M}) \in \mathbb{K}[x]$ of a matrix $\mathbf{M} \in \mathbb{K}^{m \times m}$
$\rightsquigarrow$ determinant of polynomial matrix $x\mathbf{I}_m - \mathbf{M} \in \mathbb{K}[x]^{m \times m}$

▶ fastest known algorithm uses this viewpoint [N.-Pernet, 2021]

▶ gradually transforms $x\mathbf{I}_m - \mathbf{M}$ to smaller matrices with larger degrees
**large matrices with small degrees:**

characteristic polynomial $det(xI_m - M) \in K[x]$ of a matrix $M \in K^{m \times m}$
$\rightsquigarrow$ determinant of polynomial matrix $xI_m - M \in K[x]^{m \times m}$

▶ fastest known algorithm uses this viewpoint [N.-Pernet, 2021]

▶ gradually transforms $xI_m - M$ to smaller matrices with larger degrees

**small matrices with large degree:**

extended GCD $uf + vg = \gcd(f, g)$ for polynomials $f, g \in K[x] \le d$
$\rightsquigarrow$ corresponds to a polynomial matrix transformation

$$ \begin{bmatrix} u & v \\ \tilde{g} & \tilde{f} \end{bmatrix} \begin{bmatrix} f \\ g \end{bmatrix} = \begin{bmatrix} \gcd(f, g) \\ 0 \end{bmatrix} $$

with the leftmost (polynomial) matrix of determinant in $K \setminus \{0\}$

▶ fastest known “half-gcd” algorithms use this viewpoint
[Knuth, 1970] [Schönhage, 1971] [Brent-Gustavson-Yun, 1980]
operations on sparse matrices

▶ solving sparse linear systems over $K$

▶ computing the minimal polynomial / Frobenius form

▶ introducing parallelism in these computations

[Wiedemann 1986]

[Coppersmith 1993]

[Villard 1997]

example of sparse matrix in $K^{m \times m}$
typical case: $O(m)$ nonzero entries

<figure><img src="image_95.png" /></figure>

uses polynomial matrix generator
of linearly recurrent matrix sequence
operations on structured matrices

▶ matrix-vector multiplication

▶ linear system solving

▶ nullspace computation

[Kailath-Kung-Morf 1979]

[Bostan et al. 2017]

example of Hankel matrix

~~ block-Hankel matrices

~~ Hankel-like matrices

<figure><img src="image_96.png" /></figure>

uses polynomial matrix multiplication and
matrix-Padé approximation / matrix-GCD
bivariate interpolation and multipoint evaluation

problem: given points $(\alpha_1, \beta_1), \dots, (\alpha_n, \beta_n)$ in $\mathbb{K}^2$,

▶ given $p(x, y)$, compute $p(\alpha_i, \beta_i)$ for $1 \le i \le n$

▶ find $p(x, y)$ of small degree such that $p(\alpha_i, \beta_i) = 0$

[Nüsken-Ziegler 2004]

[Beckermann 1992] [van Barel-Bultheel 1992]

[Marinari-Möller-Mora 1993]

bivariate interpolation = main step
in Reed-Solomon list-decoding
(univariate interpolation with errors)

[Guruswami-Sudan 1999]

[Kötter-Vardy 2003]

<figure><img src="image_97.png" /></figure>

uses polynomial matrix multiplication and
matrix rational reconstruction / algebraic approximants
seen as matrices over $K(x)$

linear algebra viewpoint:

matrices in $K[x]^{m×n}$ are also in $K(x)^{m×n}$

(and $K(x)$ is a field)

⇒ usual definition of addition, multiplication, determinant
these do not involve division anyway (... in algorithms?)

⇒ usual definition of rank
coincides with rank of free module

⇒ usual definition of inverse
with inverse over $K(x)$
seen as matrices over $K(x)$

linear algebra viewpoint:

matrices in $K[x]^{m×n}$ are also in $K(x)^{m×n}$

(and $K(x)$ is a field)

⇒ usual definition of addition, multiplication, determinant
these do not involve division anyway (... in algorithms?)

⇒ usual definition of rank
coincides with rank of free module

⇒ usual definition of inverse
with inverse over $K(x)$

inverse is over $K[x] \Leftrightarrow \det(A) \in K \setminus \{0\}$

def.: A is **unimodular**
# polynomial matrices

seen as matrices over $K(x)$

linear algebra viewpoint:

matrices in $K[x]^{m×n}$ are also in $K(x)^{m×n}$

(and $K(x)$ is a field)

⇒ usual definition of addition, multiplication, determinant
these do not involve division anyway (... in algorithms?)

⇒ usual definition of rank
coincides with rank of free module

⇒ usual definition of inverse
with inverse over $K(x)$

≈ algorithms may work in $K(x)^{m×n}$, but be careful with “degree explosion”!

exercise: Gaussian elimination is exponential-time
seen as matrices over $K(x)$

linear algebra viewpoint:

matrices in $K[x]^{m×n}$ are also in $K(x)^{m×n}$

(and $K(x)$ is a field)

* viewpoint useful for definitions and properties

* viewpoint hardly usable for algorithms:
  ignores degree growth + too coarse cost bounds

* cost of naive addition in $K[x]^{m×n}$ ———> $O(mn)$ additions in $K(x)$

* cost of naive multiplication in $K[x]^{m×m}$ ———> $O(m^3)$ ops in $K(x)$
seen as matrices over $K(x)$

linear algebra viewpoint:

matrices in $K[x]^{m×n}$ are also in $K(x)^{m×n}$

(and $K(x)$ is a field)

* viewpoint useful for definitions and properties

* viewpoint **hardly usable for algorithms**:
ignores degree growth + too coarse cost bounds

* cost of naive addition in $K[x]^{m×n}$ ———→ O(mn) additions in $K(x)$

* cost of naive multiplication in $K[x]^{m×m}$ ———→ O(m³) ops in $K(x)$

for algorithms&complexity, considering the degrees of entries is essential
seen as polynomials over $K^{m×n}$

polynomial viewpoint:

$K[x]^{m×n}$ is isomorphic to $K^{m×n}[x]$

$$A = \begin{bmatrix} 3x + 4 & x^3 + 4x + 1 & 4x^2 + 3 \\ 5 & 5x^2 + 3x + 1 & 5x + 3 \\ 3x^3 + x^2 + 5x + 3 & 6x + 5 & 2x + 1 \end{bmatrix}$$

$$= \begin{bmatrix} 4 & 1 & 3 \\ 5 & 1 & 3 \\ 3 & 5 & 1 \end{bmatrix} + \begin{bmatrix} 3 & 4 & 0 \\ 0 & 3 & 5 \\ 5 & 6 & 2 \end{bmatrix} x + \begin{bmatrix} 0 & 0 & 4 \\ 0 & 5 & 0 \\ 1 & 0 & 0 \end{bmatrix} x^2 + \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 3 & 0 & 0 \end{bmatrix} x^3$$
seen as polynomials over $K^{m×n}$

polynomial viewpoint:

$K[x]^{m×n}$ is isomorphic to $K^{m×n}[x]$

$$ 
\begin{aligned}
A &= \begin{bmatrix}
3x + 4 & x^3 + 4x + 1 & 4x^2 + 3 \\
5 & 5x^2 + 3x + 1 & 5x + 3 \\
3x^3 + x^2 + 5x + 3 & 6x + 5 & 2x + 1
\end{bmatrix} \\
&= \begin{bmatrix}
4 & 1 & 3 \\
5 & 1 & 3 \\
3 & 5 & 1
\end{bmatrix} + \begin{bmatrix}
3 & 4 & 0 \\
0 & 3 & 5 \\
5 & 6 & 2
\end{bmatrix} x + \begin{bmatrix}
0 & 0 & 4 \\
0 & 5 & 0 \\
1 & 0 & 0
\end{bmatrix} x^2 + \begin{bmatrix}
0 & 1 & 0 \\
0 & 0 & 0 \\
3 & 0 & 0
\end{bmatrix} x^3
\end{aligned}
$$

A has degree 3; in general, $\deg(AB) \le \deg(A) + \deg(B)$
e.g. $\deg(A^2) = 6$, and $\deg(A^3) = 8$, and $\deg(A^4) = 11$
seen as polynomials over $K^{m×n}$

polynomial viewpoint:

$K[x]^{m×n}$ is isomorphic to $K^{m×n}[x]$

$$ A = \begin{bmatrix} 3x + 4 & x^3 + 4x + 1 & 4x^2 + 3 \\ 5 & 5x^2 + 3x + 1 & 5x + 3 \\ 3x^3 + x^2 + 5x + 3 & 6x + 5 & 2x + 1 \end{bmatrix} $$ $$ = \begin{bmatrix} 4 & 1 & 3 \\ 5 & 1 & 3 \\ 3 & 5 & 1 \end{bmatrix} + \begin{bmatrix} 3 & 4 & 0 \\ 0 & 3 & 5 \\ 5 & 6 & 2 \end{bmatrix} x + \begin{bmatrix} 0 & 0 & 4 \\ 0 & 5 & 0 \\ 1 & 0 & 0 \end{bmatrix} x^2 + \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 3 & 0 & 0 \end{bmatrix} x^3 $$

degree growth enhances computational aspects

example: computing the N-th power $A^N$
seen as polynomials over $K^{m×n}$

polynomial viewpoint:

$K[x]^{m×n}$ is isomorphic to $K^{m×n}[x]$

$$ 
\begin{aligned}
\mathbf{A} &= \begin{bmatrix}
3x + 4 & x^3 + 4x + 1 & 4x^2 + 3 \\
5 & 5x^2 + 3x + 1 & 5x + 3 \\
3x^3 + x^2 + 5x + 3 & 6x + 5 & 2x + 1
\end{bmatrix} \\
&= \begin{bmatrix}
4 & 1 & 3 \\
5 & 1 & 3 \\
3 & 5 & 1
\end{bmatrix} + \begin{bmatrix}
3 & 4 & 0 \\
0 & 3 & 5 \\
5 & 6 & 2
\end{bmatrix} \times \begin{bmatrix}
0 & 0 & 4 \\
0 & 5 & 0 \\
1 & 0 & 0
\end{bmatrix} x^2 + \begin{bmatrix}
0 & 1 & 0 \\
0 & 0 & 0 \\
3 & 0 & 0
\end{bmatrix} x^3
\end{aligned}
$$

degree growth enhances computational aspects

example: computing the N-th power $\mathbf{A}^N$

repeated squaring:

$$ 
\left\{ 
\begin{array}{ll}
\mathbf{A} \times \mathbf{A} & (\deg = 3) \\
\mathbf{A}^2 \times \mathbf{A}^2 & (\deg \le 6) \\
\multicolumn{2}{c}{\vdots} \\
\mathbf{A}^{\frac{N}{2}} \times \mathbf{A}^{\frac{N}{2}} & (\deg \le \frac{3N}{4}) \\
\mathbf{A}^{\frac{N}{2}} \times \mathbf{A}^{\frac{N}{2}} & (\deg \le \frac{3N}{2})
\end{array}
\right. 
$$

find small recurrence + unroll it:

[Fajollet-Salvy 1997; Bostan-Neiger-Yurkevich 2023]

$O(N)$ operations in $\mathbb{K}$

► faster than multiplying $\mathbf{A}^{\frac{N}{2}} \times \mathbf{A}^{\frac{N}{2}}$

► does not require FFT

► prototype: $N = 2^{20} \rightsquigarrow 1.6s$ vs. $11.5s$
seen as polynomials over $K^{m×n}$

polynomial viewpoint:

$K[x]^{m×n}$ is isomorphic to $K^{m×n}[x]$

$$ A = \begin{bmatrix} 3x + 4 & x^3 + 4x + 1 & 4x^2 + 3 \\ 5 & 5x^2 + 3x + 1 & 5x + 3 \\ 3x^3 + x^2 + 5x + 3 & 6x + 5 & 2x + 1 \end{bmatrix} \\ = \begin{bmatrix} 4 & 1 & 3 \\ 5 & 1 & 3 \\ 3 & 5 & 1 \end{bmatrix} + \begin{bmatrix} 3 & 4 & 0 \\ 0 & 3 & 5 \\ 5 & 6 & 2 \end{bmatrix} x + \begin{bmatrix} 0 & 0 & 4 \\ 0 & 5 & 0 \\ 1 & 0 & 0 \end{bmatrix} x^2 + \begin{bmatrix} 0 & 1 & 0 \\ 0 & 0 & 0 \\ 3 & 0 & 0 \end{bmatrix} x^3 $$

* natural notion of degree of a polynomial matrix

* addition of $A, B \in K[x]^{m×n}$ is in $O(mnd)$ operations in $K$ where $d = \min(\deg(A), \deg(B))$

* some other polynomial operations available:

* truncation $A \text{rem } x^N$, shift $x^d A$, evaluation $A(\alpha)$ for $\alpha \in K$
seen as polynomials over $K^{m×n}$

polynomial viewpoint:

$K[x]^{m×n}$ is isomorphic to $K^{m×n}[x]$

when $m = n$, $K^{m×m}$ is a (non-commutative) ring

derived from univariate polynomial algorithms:

* **multiplication** in $K[x]^{m×m}$ seen as a product of polynomials
  complexity? $O(m^ωM(d))$ is tempting... and true for best known $M(d)$

* **truncated inversion** via power series & Newton iteration
  condition for invertibility? complexity?

* **fast Euclidean division with remainder**
  conditions for feasibility? complexity?
seen as polynomials over $K^{m×n}$

polynomial viewpoint:

$K[x]^{m×n}$ is isomorphic to $K^{m×n}[x]$

algorithmically fruitful viewpoint, with some limitations

ignores heterogeneous degrees of matrix entries

consider $\mathbf{A} = \begin{bmatrix} f(x) & a_{01} & \cdots \\ a_{10} & a_{11} & \cdots \\ \vdots & \vdots & \ddots \end{bmatrix} \in K[x]^{m×m}$,

$f(x)$ of degree $d$, other entries in $K$

► data structure: $d+1$ matrices in $K^{m×m}$

► size of representation: $m^2(d+1) \rightarrow m^2 + d$?

► adding two such matrices: $O(m^2(d+1)) \rightarrow m^2 + d$?
computer algebra

efficient algorithms and software

for matrices over a field

for univariate polynomials

polynomial matrices

basic definitions and properties

use in various situations

seen as matrices / seen as polynomials

first algorithms

exercises
## outline

### computer algebra

* efficient algorithms and software

* for matrices over a field

* for univariate polynomials

### polynomial matrices

* basic definitions and properties

* use in various situations

* seen as matrices / seen as polynomials

### first algorithms

* exploiting evaluation-interpolation

* extending algorithms for polynomials

* partial linearization techniques

### exercises
fast multiplication

naive multiplication: $O(m^3d^2)$ operations in $K$ $O(m^\omega M(d))$?
# On fast multiplication of polynomials over arbitrary algebras

David G. Cantor¹ and Erich Kaltofen²*

¹ Department of Mathematics, University of California, Los Angeles, CA 90024-1555, USA

² Department of Computer Science, Rensselaer Polytechnic Institute, Troy, NY 12180-3590, USA

Received January 22, 1988 / May 10, 1991

## 1 Introduction

In this paper we generalize the well-known Schönhage-Strassen algorithm for multiplying large integers to an algorithm for multiplying polynomials with coefficients from an arbitrary, not necessarily commutative, not necessarily associative, algebra $\mathcal{A}$. Our main result is an algorithm to multiply polynomials of degree $<n$ in $O(n \log n)$ algebra multiplications and $O(n \log n \log \log n)$ algebra additions/subtractions (we count a subtraction as an addition). The constant implied by the “O” does not depend upon the algebra $\mathcal{A}$. The parallel complexity of our algorithm, i.e., the depth of the corresponding arithmetic circuit, is
naive multiplication: $O(m^3d^2)$ operations in $\mathbb{K}$ $O(m^\omega M(d))$?

# On fast multiplication of polynomials over arbitrary algebras

David G. Cantor¹ and Erich Kaltofen²*

¹Department of Mathematics, University of California, Los Angeles, CA 90024-1555, USA

²Department of Computer Science, Rensselaer Polytechnic Institute, Troy, NY 12180-3590, USA

Received January 22, 1988 / May 10, 1991

## 1 Introduction

multiplication in $\mathbb{K}^{m \times m}[x]$ with degree $\le$ d:

* $O(d \log(d))$ multiplications in $\mathbb{K}^{m \times m}$

* $O(d \log(d) \log(d))$ additions in $\mathbb{K}^{m \times m}$

$MM(m, d) \in O(m^\omega d \log(d) + m^2 d \log(d) \log(d))$

In this paper we generalize the well-known Schönhage-Strassen algorithm for multiplying large integers to an algorithm for multiplying polynomials with coefficients from an arbitrary, not necessarily commutative, not necessarily associative, algebra $\mathcal{A}$. Our main result is an algorithm to multiply polynomials of degree $<n$ in $O(n \log n)$ algebra multiplications and $O(n \log n \loglog n)$ algebra additions/subtractions (we count a subtraction as an addition). The constant implied by the “O” does not depend upon the algebra $\mathcal{A}$. The parallel complexity of our algorithm, i.e., the depth of the corresponding arithmetic circuit, is
exploiting evaluation-interpolation

**exercise: multiplication, determinant, inversion**

1. adapting the evaluation-interpolation paradigm to
matrices in $K[x]^{m×m}$,

▶ give an explicit **multiplication** algorithm

▶ give a **determinant** algorithm

▶ give an **inversion** algorithm
computing the inverse over the fractions $K(x)$

2. for each of these algorithms,

▶ give a required lower bound on the **cardinality** of $K$

▶ state and prove an upper bound on the **complexity**

*hint*: use **known degree bounds on the output**
**exercise: multiplication, determinant, inversion**

1. adapting the evaluation-interpolation paradigm to
matrices in $K[x]^{m×m}$,

▶ give an explicit **multiplication** algorithm

▶ give a **determinant** algorithm

▶ give an **inversion** algorithm
computing the inverse over the fractions $K(x)$

2. for each of these algorithms,

▶ give a required lower bound on the **cardinality** of $K$

▶ state and prove an upper bound on the **complexity**

**multiplication:** for large enough $K$,
$MM(m, d) \in O(m^\omega d + m^2M(d))$ [Bostan-Schost 2005]

≈ better than $m^\omega M(d)$
exploiting evaluation-interpolation

**exercise: multiplication, determinant, inversion**

1. adapting the evaluation-interpolation paradigm to matrices in $K[x]^{m×m}$,

* give an explicit **multiplication** algorithm

* give a **determinant** algorithm

* give an **inversion** algorithm

computing the inverse over the fractions $K(x)$

2. for each of these algorithms,

* give a required lower bound on the **cardinality** of $K$

* state and prove an upper bound on the **complexity**

<table><thead><tr><th></th><th>evaluation-interpolation, large <i>K</i></th><th>best known, unconditional</th></tr></thead><tbody><tr><th scope="row">determinant</th><td>O~(m<sup>ω+1</sup>d)</td><td>O~(m<sup>ω</sup>d)</td></tr><tr><th scope="row">inversion</th><td>O~(m<sup>ω+1</sup>d)</td><td>O~(m<sup>3</sup>d)</td></tr><tr><td></td><td>reductions to PolMul&amp;MatMul</td><td>reductions to PolMatMul</td></tr></tbody></table>
**Entrée** Un entier $N > 0$, F mod $X^N$ une série tronquée.

**Sortie** $F^{-1}$ mod $X^N$.

Si $N = 1$, alors renvoyer $f_0^{-1}$, où $f_0 = F(0)$.

Sinon :

1. Calculer récursivement l'inverse $G$ de $F$ mod $X^{[N/2]}$.

2. Renvoyer $G + (1 - GF)G$ mod $X^N$.

Algorithme 3.2 – Inverse de série par itération de Newton.

Convergence quadratique pour l'inverse d'une série formelle

**Lemme 3.2** Soient $A$ un anneau non nécessairement commutatif, $F \in A[[X]]$ une série formelle de terme constant inversible et $G$ une série telle que $G - F^{-1} = O(X^n)$ ($n \ge 1$). Alors la série

$$ \mathcal{N}(G) = G + (1 - GF)G \quad (3.2) $$

vérifie $\mathcal{N}(G) - F^{-1} = O(X^{2n})$.
truncated inversion — results

consider a (square) polynomial matrix $A \in K[x]^{m \times m}$

► A is invertible as a power series
  ⇔ its constant term $A(0) \in K^{m \times m}$ is invertible

► if A is invertible as a power series, computing $A^{-1} \mod x^N$ costs $O(MM(m, N))$ operations in $K$

► no additional log: $MM(m, \frac{N}{2}) + MM(m, \frac{N}{4}) + MM(m, \frac{N}{8}) + \dots$

► excellent reduction to PolMatMul!

► timings with the Polynomial Matrix Library:

<table><thead><tr><td>m</td><td>d</td><td>PolMatMul</td><td>TruncInv</td></tr></thead><tbody><tr><td>10</td><td>20000</td><td>0.203</td><td>0.551</td></tr><tr><td>20</td><td>5000</td><td>0.225</td><td>0.639</td></tr><tr><td>40</td><td>2500</td><td>0.528</td><td>1.424</td></tr><tr><td>80</td><td>1250</td><td>1.227</td><td>3.653</td></tr></tbody></table>
division with remainder

**problem:**

given $\mathbf{A}, \mathbf{B} \in \mathbb{K}^{m \times m}[x]$,

compute $\mathbf{Q}, \mathbf{R} \in \mathbb{K}^{m \times m}[x]$ such that

$$ \mathbf{A} = \mathbf{BQ} + \mathbf{R} \quad \text{and} \quad \deg(\mathbf{R}) < \deg(\mathbf{B}) $$

... are we not missing an assumption?
division with remainder

problem:

given $\mathbf{A}, \mathbf{B} \in \mathbb{K}^{m \times m}[x]$, compute $\mathbf{Q}, \mathbf{R} \in \mathbb{K}^{m \times m}[x]$ such that $\mathbf{A} = \mathbf{BQ} + \mathbf{R}$ and $\deg(\mathbf{R}) < \deg(\mathbf{B})$

... are we not missing an assumption?

rule 1: dividing by zero is generally a bad idea

rule 2: if you think you need to divide by zero, refer to rule 1

rule 3: neglecting to check that something is not zero does not make it nonzero etc. etc.
division with remainder

**problem:**

given $\mathbf{A}, \mathbf{B} \in \mathbb{K}^{m \times m}[x]$,

compute $\mathbf{Q}, \mathbf{R} \in \mathbb{K}^{m \times m}[x]$ such that

$$ \mathbf{A} = \mathbf{BQ} + \mathbf{R} \quad \text{and} \quad \deg(\mathbf{R}) < \deg(\mathbf{B}) $$

... are we not missing an assumption?

for a polynomial $p \in \mathcal{A}[x]$, over some ring $\mathcal{A}$, division by $p$ is feasible

► if $p$ is monic (leading coefficient $1_\mathcal{A}$)

► and more generally if the leading coefficient of $p$ is invertible in $\mathcal{A}$

assumption: the leading coefficient of $\mathbf{B}$ is invertible in $\mathbb{K}^{m \times m}$

recall $\mathbf{B} = \mathbf{B}_0 + \mathbf{B}_1 x + \dots + \mathbf{B}_d x^d$ with $\mathbf{B}_i \in \mathbb{K}^{m \times m}$
division with remainder

**problem:**

given $\mathbf{A}, \mathbf{B} \in \mathbb{K}^{m \times m}[x]$ with $\text{lc}(\mathbf{B})$ invertible,
compute $\mathbf{Q}, \mathbf{R} \in \mathbb{K}^{m \times m}[x]$ such that
$\mathbf{A} = \mathbf{B}\mathbf{Q} + \mathbf{R}$ and $\text{deg}(\mathbf{R}) < \text{deg}(\mathbf{B})$

► under this assumption, the usual fast Euclidean algorithm works

► recall:

1. reverse the equation,

2. compute quotient by truncated inverse multiplication
   $\tilde{\mathbf{Q}} = \mathbf{B}^{-1}\tilde{\mathbf{A}} \bmod \chi^{d_A - d_B + 1}$

3. deduce remainder

► complexity is $O(\underbrace{(MM(m, d_A - d_B))}_{\text{find } \mathbf{Q}} + \underbrace{MM(m, d_B))}_{\text{find } \mathbf{R}})$
division with remainder

problem:

given $A, B \in \mathbb{K}^{m \times m}[x]$ with $\text{lc}(B)$ invertible,

compute $Q, R \in \mathbb{K}^{m \times m}[x]$ such that

$$A = BQ + R \quad \text{and} \quad \text{deg}(R) < \text{deg}(B)$$

<table><thead><tr><th>m</th><th>d<sub>A</sub></th><th>d<sub>B</sub></th><th>PolMatMul<br>in deg d<sub>B</sub></th><th>TruncInv<br>in deg d<sub>B</sub></th><th>QuoRem</th></tr></thead><tbody><tr><td>10</td><td>40000</td><td>20000</td><td>0.203</td><td>0.551</td><td>1.873</td></tr><tr><td>20</td><td>10000</td><td>5000</td><td>0.225</td><td>0.639</td><td>2.164</td></tr><tr><td>40</td><td>5000</td><td>2500</td><td>0.528</td><td>1.424</td><td>6.468</td></tr><tr><td>80</td><td>2500</td><td>1250</td><td>1.227</td><td>3.653</td><td>15.59</td></tr></tbody></table>
division with remainder

problem:

given $A, B \in \mathbb{K}^{m \times m}[x]$ with $\text{lc}(B)$ invertible,
compute $Q, R \in \mathbb{K}^{m \times m}[x]$ such that
$A = BQ + R$ and $\text{deg}(R) < \text{deg}(B)$

```R
44 // step 1: reverse input matrices
45 row(Brev, B, rdegB);
46 row(buf, A, rdegA);
47
48 // step 2: compute quotient
49 // Qrev = Brev^{-1} R mod X^{d+1}
50 solve_series(Qrev, Brev, buf, d+1);
51 reverse(0, Qrev, d);
52
53 // step 3: deduce remainder
54 // R = A - B*Q
55 multiply(buf, B, Q);
56 sub(R, A, buf);
```

▸ an efficient reduction, again

▸ rdegA vs. degree d?

▸ row娱 reverse vs. reverse ?

▸ refinement of matrix degree:
**row- or column-wise degrees**
≈ improves applicability & complexity

e.g. division by $B = \text{diag}(x^{d_1}, \dots, x^{d_m})$
row degree of a polynomial matrix
= the list of the maximum degree in each of its rows

for $\mathbf{A} = (a_{i,j}) \in \mathbb{K}[x]^{m \times n}$,

$$rdeg(\mathbf{A}) = (rdeg(\mathbf{A}_{1,\*}), \dots, rdeg(\mathbf{A}_{m,\*}))$$

$$= \left( \max_{1 \le j \le n} \deg(\mathbf{A}_{1,j}), \dots, \max_{1 \le j \le n} \deg(\mathbf{A}_{m,j}) \right) \in \mathbb{Z}^m$$
**row degree** of a polynomial matrix

= the list of the maximum degree in each of its rows

**column degree** of a polynomial matrix

= the list of the maximum degree in each of its columns
**row degree** of a polynomial matrix

= the list of the maximum degree in each of its rows

**column degree** of a polynomial matrix

= the list of the maximum degree in each of its columns

sum of degrees of all entries $\le \frac{n \times \text{sum of row degrees}}{m \times \text{sum of column degrees}} \le mn \times \text{global degree}$

with notation:

$$ \sum_{i,j} \deg(a_{ij}) \le \frac{n|r\deg(\mathbf{A})|}{m|c\deg(\mathbf{A})|} \le mn \deg(\mathbf{A}) $$
**row degree** of a polynomial matrix
= the list of the maximum degree in each of its rows

**column degree** of a polynomial matrix
= the list of the maximum degree in each of its columns

sum of degrees of all entries $\le$ $n \times$ sum of row degrees
$\le$ $m \times$ sum of column degrees
$\le$ $mn \times$ global degree

with notation:

$$ \sum_{i,j} \deg(a_{ij}) \le \frac{n|r\deg(\mathbf{A})|}{m|c\deg(\mathbf{A})|} \le mn \deg(\mathbf{A}) $$

consider $\mathbf{A}$ with degree matrix

$$ \begin{pmatrix} 100 & 5 & 20 & 1 \\ 100 & 5 & 20 & 1 \\ 100 & 5 & 20 & 1 \\ 100 & 5 & 20 & 1 \end{pmatrix} $$

determinant of $\mathbf{A}$: degree $\le 126$

$$ \rightsquigarrow \text{ better than naive bound } 4 \deg(\mathbf{A}) = 400 $$
**row degree** of a polynomial matrix
= the list of the maximum degree in each of its rows

**column degree** of a polynomial matrix
= the list of the maximum degree in each of its columns

sum of degrees of all entries $\le \frac{n \times \text{sum of row degrees}}{m \times \text{sum of column degrees}} \le mn \times \text{global degree}$

with notation:

$$ \sum_{i,j} \deg(a_{ij}) \le \frac{n|\mathrm{rdeg}(A)|}{m|\mathrm{cdeg}(A)|} \le mn \deg(A) $$

more general division with remainder:

▶ take for $\mathrm{lc}(\mathbf{B})$ **row-wise** leading coefficients

▶ if $\mathrm{lc}(\mathbf{B})$ is invertible, division by $\mathbf{B}$ is feasible

▶ with **row-wise** degree bounds on remainder

$$ \begin{bmatrix} 4x^3 + 2x + 2 & 6x^3 + 2x^2 + 5 \\ 4x^2 + 2 & 3x^3 + x + 3 \end{bmatrix} = \begin{bmatrix} x & 0 \\ 0 & 2x^2 \end{bmatrix} \mathbf{Q} + \begin{bmatrix} 2 & 5 \\ 2 & x + 3 \end{bmatrix} $$
reduce unbalanced degrees to some average degree

where degree means row degree, column degree, or related refined measures

[Storjohann 2006] [Zhou-Labahn 2012] [Jeannerod-Neiger-Villard 2020]

**typical properties:**

from a matrix $A \in K[x]^{m \times m}$ with $D = |rdeg(A)| \ll m \deg(A)$ construct a matrix $\bar{A} \in K[x]^{m' \times m'}$ with

► a slight increase of matrix dimension: $m \le \text{m}' \le 2m$

► a strong decrease of matrix degree: $\deg(\bar{A}) \le 2 \frac{D}{m}$

► preservation of the features targeted by our computations

**examples:**

► product $AB$ easily deduced from product $\bar{A}\bar{B}$

► preservation of the determinant $\det(A) = \det(\bar{A})$

► inverse of $\bar{A}$ contains inverse of $A$ as submatrix

► ...
reduce unbalanced degrees to some average degree

**basic illustration:**

► let $\mathbf{A} \in \mathbb{K}[x]^{m \times m}$ of degree < d,

► let $\mathbf{u} \in \mathbb{K}[x]^{m \times 1}$ of degree < md,

then the matrix-vector product $\mathbf{Au}$ can be computed in
$\mathbf{MM}(\mathbf{m}, \mathbf{d}) + \mathbf{O}(\mathbf{m}^2\mathbf{d})$ operations in $\mathbb{K}$

what would be the cost of the "naive" multiplication? $\rightsquigarrow \mathbf{O}(\mathbf{m}^2\mathbf{M}(\mathbf{m}\mathbf{d}))$

**algorithm:**

[Lecerf 2001 (in communication + software)]
reduce unbalanced degrees to some average degree

**basic illustration:**

► let $\mathbf{A} \in \mathbb{K}[x]^{m \times m}$ of degree $< d$,

► let $\mathbf{u} \in \mathbb{K}[x]^{m \times 1}$ of degree $< md$,

then the matrix-vector product $\mathbf{Au}$ can be computed in $\mathbf{MM}(m, d) + \mathbf{O}(m^2d)$ operations in $\mathbb{K}$

what would be the cost of the “naive” multiplication? $\rightsquigarrow \mathbf{O}(m^2M(md))$

**algorithm:**

[Lecerf 2001 (in communication + software)]

$$ \left[ \mathbf{A} \right] \left[ \mathbf{u} \right] = \left[ \mathbf{A} \right] \left[ \bar{\mathbf{U}} \right] \left[ \begin{matrix} 1 \\ \chi^d \\ \chi^{2d} \\ \vdots \end{matrix} \right] $$

where the columns of $\bar{\mathbf{U}} \in \mathbb{K}[x]^{m \times m}$ form the $\chi^d$-adic expansion of $\mathbf{u}$
$\Rightarrow$ here $\deg(\bar{\mathbf{U}}) < d$
reduce unbalanced degrees to some average degree

**basic illustration:**

► let $A \in K[x]^{m \times m}$ of degree < d,

► let $u \in K[x]^{m \times 1}$ of degree < md,

then the matrix-vector product Au can be computed in $MM(m, d) + O(m^2d)$ operations in $K$

what would be the cost of the “naive” multiplication? $\rightsquigarrow O(m^2M(md))$

**algorithm:**

[Lecerf 2001 (in communication + software)]

<table><thead><tr><th>m</th><th>d</th><th>md</th><th>via PolMatMul</th><th>matrix-vector</th></tr></thead><tbody><tr><td>10</td><td>20000</td><td>200000</td><td>0.203</td><td>0.368</td></tr><tr><td>20</td><td>5000</td><td>100000</td><td>0.225</td><td>0.683</td></tr><tr><td>40</td><td>2500</td><td>100000</td><td>0.528</td><td>2.481</td></tr><tr><td>80</td><td>1250</td><td>100000</td><td>1.227</td><td>9.592</td></tr></tbody></table>
## outline

### computer algebra

* efficient algorithms and software

* for matrices over a field

* for univariate polynomials

### polynomial matrices

* basic definitions and properties

* use in various situations

* seen as matrices / seen as polynomials

### first algorithms

* exploiting evaluation-interpolation

* extending algorithms for polynomials

* partial linearization techniques

### exercises
computer algebra

efficient algorithms and software

for matrices over a field

for univariate polynomials

polynomial matrices

basic definitions and properties

use in various situations

seen as matrices / seen as polynomials

first algorithms

exploiting evaluation-interpolation

extending algorithms for polynomials

partial linearization techniques

exercises

evaluation-interpolation-based algorithms

Krylov iterates via repeated squaring

Krylov iterates in MatMul time
evaluation-interpolation-based algorithms

**exercise: multiplication, determinant, inversion**

1. adapting the evaluation-interpolation paradigm to matrices in $K[x]^{m×m}$,

* give an explicit **multiplication** algorithm

* give a **determinant** algorithm

* give an **inversion** algorithm
computing the inverse over the fractions $K(x)$

2. for each of these algorithms,

* give a required lower bound on the **cardinality** of $K$

* state and prove an upper bound on the **complexity**

*hint*: use **known degree bounds on the output**
exercises

evaluation-interpolation: multiplication

given **A** and **B** in $K[x]^{m×m}$ of degree $≤ d$,
we know that $C = AB$ has degree at most 2d, so:

1. **pick points:** pairwise distinct $α_1, ..., α_{2d+1} ∈ K$ \quad $Card(K) ≥ 2d + 1$

2. **evaluate:** $A(α_i)$ and $B(α_i)$, for $i = 1, ..., 2d + 1$ \quad $O(m^2M(d)log(d))$

3. **multiply:** $A(α_i)B(α_i)$, for $i = 1, ..., 2d + 1$ \quad $O(m^ωd)$

4. **interpolate:** find $C$ in $K[x]^{m×m}$ of degree $≤ 2d$ such that $C(α_i) = A(α_i)B(α_i)$, for $i = 1, ..., 2d + 1$ \quad $O(m^2M(d)log(d))$

5. **return** $C$

excellent algorithm:

1. linear in $d$ in the term $m^ωd$ (recall Cantor-Kaltofen: $m^ωdlog(d)$)

2. exponent $ω$ of matrix multiplication

3. the $m^2M(d)log(d)$ term can be improved via points in geometric sequence

4. downside: restriction on $K$ (large degrees + small finite fields does happen)
given $\mathbf{A}$ in $K[x]^{m \times m}$ of degree $\le d$,
we know that $\Delta = \det(\mathbf{A})$ has degree at most $md$, so:

1. pick points: pairwise distinct $\alpha_1, \dots, \alpha_{md+1} \in K$

2. evaluate: $\mathbf{A}(\alpha_i)$ for $i = 1, \dots, md + 1$

3. determinant: $\beta_i = \det(\mathbf{A}(\alpha_i))$, for $i = 1, \dots, md + 1$

4. interpolate: find $\Delta$ in $K[x]$ of degree $\le md$ such that
$\Delta(\alpha_i) = \beta_i$, for $i = 1, \dots, md + 1$

5. return $\Delta$

$$ \mathrm{Card}(\mathbb{K}) \geq md + 1 $$

$$ \mathrm{O}(m^3 M(d) \log(d)) $$

$$ \mathrm{O}(m^{\omega+1} d) $$

$$ \mathrm{O}(M.md) \log.md) $$

* quasi-linear in degree $d$: fast for large $d$, small $m$

* exponent $> 3$ on matrix dimension $m$: slow for large $m$

* best known today: $O(m^\omega d)$
evaluation-interpolation: inversion

given $\mathbf{A}$ in $K[x]^{m \times m}$ of degree $\le d$,
we know that $\mathbf{C} = \mathbf{A}^{-1} = \frac{1}{\Delta} \mathbf{U}$ with
$\deg(\Delta) \le \text{md}$ and $\deg(\mathbf{U}) \le (m-1)d$, so:

0. set $n = (2m-1)d+1$ \quad $n = \Theta(
\mathrm{md})$

1. pick points: pairwise distinct $\alpha_1, \dots, \alpha_n \in K$ \quad $\mathrm{Card}(K) \ge (2m-1)d+1$

2. evaluate: $\mathbf{A}(\alpha_i)$, for $i = 1, \dots, n$ \quad $\mathbf{O}(m^3 M(d) \log(d))$

3. invert: $\mathbf{A}(\alpha_i)^{-1}$, for $i = 1, \dots, n$ \quad $\mathbf{O}(m^{\omega+1} d)$

4. interpolate: using Cauchy interpolation find $\mathbf{C}$ in $K(X)^{m \times m}$ with all numerators of degree $\le (m-1)d$ and all denominators of degree $\le \text{md}$ such that $\mathbf{C}(\alpha_i) = \mathbf{A}(\alpha_i)^{-1}$, for $i = 1, \dots, n$ \quad $\mathbf{O}(m^2 M(\text{md}) \log(\text{md}))$

5. return $\mathbf{C}$

. quasi-linear in degree d: fast for large d, small m

. exponent > 3 on dimension m but recall size of A-1 is typically Θ(m³d)

. best known today: $O'(m^3d)$, and even $O'(m^\omega d)$ for factorized form

. note: one could compute det($\mathbf{A}$) to avoid Cauchy interpolation
**problem (Krylov iterates):**

*input*: matrix $A \in \mathbb{K}^{m \times m}$,
vector $v \in \mathbb{K}^{m \times 1}$
integer $d > 0$

*output*: $v, Av, \dots, A^{d-1}v$

**kernel black box:**

given a matrix $\mathbf{F} \in \mathbb{K}[x]^{m \times (m+1)}$ of rank $m$ and degree $\le 1$, one can compute a nonzero element of degree $\le m$ in the right kernel of $\mathbf{F}$ using $O(m^\omega)$ operations in $\mathbb{K}$

[refined analysis of Algo.1 in Zhou-Labahn-Storjohann 2012]

1. give an algorithm which costs $O(m^\omega \log(d) + m^{\omega-1}d)$ operations in $\mathbb{K}$, based on repeated squaring

2. prove that the generating series of $(A^k v)_{k \ge 0}$ rewrites as a fraction of polynomial matrices:

$$ \sum_{k \ge 0} A^k v x^k = (I - xA)^{-1}v $$

3. using the kernel black box, give a complexity bound for finding $\lambda \in \mathbb{K}[x]$ and $\mathbf{u} \in \mathbb{K}[x]^{m \times 1}$, both of degree $\le m$, such that

$$ \sum_{k \ge 0} A^k v x^k = \mathbf{u}/\lambda $$

4. show that $(A^k v)_{0 \le k < d}$ can be computed in $O(m^\omega + mM(d))$
**problem (Krylov iterates):**

*input*: matrix $A \in \mathbb{K}^{m \times m}$,
vector $v \in \mathbb{K}^{m \times 1}$
integer $d > 0$

*output*: $v, Av, \dots, A^{d-1}v$

**kernel black box:**

given a matrix $\mathbf{F} \in \mathbb{K}[\mathbf{x}]^{m \times (m+1)}$ of rank $m$
and degree $\le 1$, one can compute a nonzero
element of degree $\le m$ in the right kernel of
$\mathbf{F}$ using $O(m^\omega)$ operations in $\mathbb{K}$

[refined analysis of Algo.1 in Zhou-Labahn-Storjohann 2012]

1. give an algorithm which costs $O(m^\omega \log(d) + m^{\omega-1}d)$ operations in $\mathbb{K}$, based on repeated squaring

for simplicity, take $d$ a power of 2

first compute $A^2, A^4, \dots, A^{d/2}$, cost $O(m^\omega \log(d))$

from $v$, compute $Av$

from $[v \ Av]$, compute $A^2[v \ Av] = [A^2v \ A^3v]$

from $[v \ Av \ A^2v \ A^3v]$, compute $A^4[v \ Av \ A^2v \ A^3v] = [A^4v \ A^5v \ A^6v \ A^7v]$

etc...

from $[A^k v]_{0 \le k < d/2}$, compute $A^{d/2} [A^k v]_{0 \le k < d/2} = [A^k v]_{d/2 \le k < d}$
**problem (Krylov iterates):**

*input*: matrix $A \in \mathbb{K}^{m \times m}$,
vector $v \in \mathbb{K}^{m \times 1}$
integer $d > 0$

*output*: $v, Av, \dots, A^{d-1}v$

**kernel black box:**

given a matrix $\mathbf{F} \in \mathbb{K}[\mathbf{x}]^{m \times (m+1)}$ of rank $m$
and degree $\le 1$, one can compute a nonzero
element of degree $\le m$ in the right kernel of
$\mathbf{F}$ using $O(m^\omega)$ operations in $\mathbb{K}$

[refined analysis of Algo.1 in Zhou-Labahn-Storjohann 2012]

1. give an algorithm which costs $O(m^\omega \log(d) + m^{\omega-1}d)$ operations in $\mathbb{K}$, based on repeated squaring

the first min($\log(d)$, $\log(m)$) products involve matrices of dimensions $m$ or less, hence a
total cost bounded by $O(m^\omega \log(d))$

the remaining products (if any) involve a lefthand operand of dimensions $m \times m$ and a
righthand one of dimensions $m \times 2^k$, where $k$ goes from about $\log_2(m)$ to $\log_2(d)$

$\rightsquigarrow$ for a given $k$, the product costs $O(m^{\omega-1}2^k)$

$\rightsquigarrow$ summing this over all $k$, with $\sum_{k \le \log_2(d)} 2^k \in O(d)$, gives $O(m^{\omega-1}d)$
**problem (Krylov iterates):**

*input*: matrix $A \in \mathbb{K}^{m \times m}$,
vector $v \in \mathbb{K}^{m \times 1}$
integer $d > 0$

*output*: $v, Av, \dots, A^{d-1}v$

**kernel black box:**

given a matrix $\mathbf{F} \in \mathbb{K}[x]^{m \times (m+1)}$ of rank $m$
and degree $\le 1$, one can compute a nonzero
element of degree $\le m$ in the right kernel of
$\mathbf{F}$ using $O(m^\omega)$ operations in $\mathbb{K}$

[refined analysis of Algo.1 in Zhou-Labahn-Storjohann 2012]

2. prove that the generating series of $(A^k v)_{k \ge 0}$ rewrites as a
fraction of polynomial matrices:

$$
\sum_{k \ge 0} A^k v x^k = (I - xA)^{-1} v
$$

multiply the left-hand side by $I - xA$, this yields $v$
**problem (Krylov iterates):**

*input*: matrix $A \in \mathbb{K}^{m \times m}$,
vector $v \in \mathbb{K}^{m \times 1}$
integer $d > 0$

*output*: $v, Av, \dots, A^{d-1}v$

**kernel black box:**

given a matrix $\mathbf{F} \in \mathbb{K}[x]^{m \times (m+1)}$ of rank $m$
and degree $\le 1$, one can compute a nonzero
element of degree $\le m$ in the right kernel of
$\mathbf{F}$ using $O(m^\omega)$ operations in $\mathbb{K}$

[refined analysis of Algo.1 in Zhou-Labahn-Storjohann 2012]

3. using the kernel black box, give a complexity bound for finding
$\lambda \in \mathbb{K}[x]$ and $\mathbf{u} \in \mathbb{K}[x]^{m \times 1}$, both of degree $\le m$, such that

$$ \sum_{k \ge 0} \lambda^k v \mathbf{x}^k = \mathbf{u} / \lambda $$

. consider $\mathbf{F} = [\mathbf{I} - x\mathbf{A} - v]$; this matrix has degree $\le 1$ and rank $m$ (its leftmost $m \times m$ submatrix is nonsingular)

. so, in $O(m^\omega)$, we can compute a nonzero element of degree $\le m$ in its right kernel

. this element can be written $[\frac{u}{\lambda}]$, and $\mathbf{F}[\frac{u}{\lambda}] = 0$ rewrites as $(\mathbf{I} - x\mathbf{A})\mathbf{u} = v\lambda$

. observe that $\lambda$ cannot be zero (otherwise, $\mathbf{u}$ would be a nonzero vector in the right kernel of $\mathbf{I} - x\mathbf{A}$, which is not possible)

. thus $(\mathbf{I} - x\mathbf{A})^{-1}v = \frac{1}{\lambda}\mathbf{u}$
**problem (Krylov iterates):**

*input*: matrix $A \in \mathbb{K}^{m \times m}$,
vector $v \in \mathbb{K}^{m \times 1}$
integer $d > 0$

*output*: $v, Av, \dots, A^{d-1}v$

**kernel black box:**

given a matrix $\mathbf{F} \in \mathbb{K}[x]^{m \times (m+1)}$ of rank $m$ and degree $\le 1$, one can compute a nonzero element of degree $\le m$ in the right kernel of $\mathbf{F}$ using $O(m^\omega)$ operations in $\mathbb{K}$

[refined analysis of Algo.1 in Zhou-Labahn-Storjohann 2012]

4. show that $(A^k v)_0 \le k < d$ can be computed in $O(m^\omega + mM(d))$

these $d$ vectors are the first $d$ terms of the series $\sum_{k \ge 0} A^k v x^k$

we have seen that this series is equal to $\frac{1}{\lambda}u$ (with $u$ and $\lambda$ found in $O(m^\omega)$) $\rightsquigarrow$ it suffices to expand $u/\lambda$ as a power series in precision $d$

since $u$ is a vector of $m$ entries, this costs $O(mM(d))$
summary

computer algebra

efficient algorithms and software

for matrices over a field

for univariate polynomials

polynomial matrices

basic definitions and properties

use in various situations

seen as matrices / seen as polynomials

first algorithms

exploiting evaluation-interpolation

extending algorithms for polynomials

partial linearization techniques

exercises

evaluation-interpolation-based algorithms

Krylov iterates via repeated squaring

Krylov iterates in MatMul time