Add functions for computation of krylov rank profile and linear interpolation basis

[Biffo89]

This PR implements the following functions, adds aliases row_rank_profile for pivot_rows and column_rank_profile for pivots, and caches the row rank profile that is calculated as part of linbox_echelonize in matrix_modn_dense_template.pxi.

Given a $\mathbb{K}[x]$-module $V$ with multiplication matrix $J \in \mathbb{K}^{\sigma\times\sigma}$, a set of $m$ elements of $V$ expressed in matrix form $E \in \mathbb{K}^{m\times\sigma}$, a priority shift $s = (s_1,\dots,s_m) \in \mathbb{Z}^m$, and an upper bound $d$ on the degree of the minimal polynomial of $J$, krylov_rank_profile finds the lexicographical first maximal linearly independent set of rows from the matrix, with rows ordered according to the priority $P(E_{c,*} J^d) = s_c + d$:

[   E  ]
[  EJ  ]
[  ... ]
[ EJ^δ ]

It returns the positions of these rows in the permuted matrix, the submatrix formed by these rows, and the lexicographical first maximal linearly independent set of columns of the submatrix.

Given the same arguments and a variable name, linear_interpolation_basis calculates an $s$-minimal interpolation basis for $(E,J)$, i.e. an $s$-reduced basis for the space of polynomials satisfying $p_1 E_{1,*}+\dots+p_mE_{m,*}=0$.

📝 Checklist

• The title is concise and informative.
• The description explains in detail what this PR is about.
• I have linked a relevant issue or discussion.
• I have created tests covering the changes.
• I have updated the documentation and checked the documentation preview.

[Biffo89]

Some benchmarks are provided here (done in a branch with #40423 and #40435 merged):

GF(2)

matrix multiplication
+------+---------+
|    n | time    |
+======+=========+
|  128 | 31.4 μs |
+------+---------+
|  256 | 40 μs   |
+------+---------+
|  512 | 181 μs  |
+------+---------+
| 1024 | 980 μs  |
+------+---------+
| 2048 | 10.8 ms |
+------+---------+
| 4096 | 73.5 ms |
+------+---------+

uniform profile
+-------------+--------+--------+---------+--------+--------+
|   *   sigma | 256    | 512    | 1024    | 2048   | 4096   |
|       *     |        |        |         |        |        |
|   m       * |        |        |         |        |        |
+=============+========+========+=========+========+========+
|           1 |        |        | 121 ms  | 209 ms | 1.36 s |
+-------------+--------+--------+---------+--------+--------+
|           4 |        |        | 92.2 ms | 478 ms | 1.27 s |
+-------------+--------+--------+---------+--------+--------+
|          16 |        |        | 113 ms  | 243 ms | 1.67 s |
+-------------+--------+--------+---------+--------+--------+
|         256 | 115 ms | 236 ms | 822 ms  |        |        |
+-------------+--------+--------+---------+--------+--------+
|         512 | 207 ms | 459 ms | 1.18 s  |        |        |
+-------------+--------+--------+---------+--------+--------+
|        1024 | 404 ms | 1.13 s | 2.75 s  |        |        |
+-------------+--------+--------+---------+--------+--------+

hermite profile
+-------------+--------+--------+---------+--------+--------+
|   *   sigma | 256    | 512    | 1024    | 2048   | 4096   |
|       *     |        |        |         |        |        |
|   m       * |        |        |         |        |        |
+=============+========+========+=========+========+========+
|           1 |        |        | 53.7 ms | 451 ms | 1.28 s |
+-------------+--------+--------+---------+--------+--------+
|           4 |        |        | 92.2 ms | 364 ms | 2.27 s |
+-------------+--------+--------+---------+--------+--------+
|          16 |        |        | 335 ms  | 488 ms | 3.25 s |
+-------------+--------+--------+---------+--------+--------+
|         256 | 309 ms | 464 ms | 789 ms  |        |        |
+-------------+--------+--------+---------+--------+--------+
|         512 | 400 ms | 898 ms | 1.78 s  |        |        |
+-------------+--------+--------+---------+--------+--------+
|        1024 | 515 ms | 1.41 s | 2.76 s  |        |        |
+-------------+--------+--------+---------+--------+--------+

uniform basis
+-------------+--------+--------+--------+--------+--------+--------+
|   *   sigma | 128    | 256    | 512    | 1024   | 2048   | 4096   |
|       *     |        |        |        |        |        |        |
|   m       * |        |        |        |        |        |        |
+=============+========+========+========+========+========+========+
|           1 |        |        |        | 146 ms | 249 ms | 1.7 s  |
+-------------+--------+--------+--------+--------+--------+--------+
|           4 |        |        |        | 158 ms | 263 ms | 1.65 s |
+-------------+--------+--------+--------+--------+--------+--------+
|          16 |        |        |        | 253 ms | 628 ms | 2 s    |
+-------------+--------+--------+--------+--------+--------+--------+
|         128 | 177 ms | 225 ms | 779 ms |        |        |        |
+-------------+--------+--------+--------+--------+--------+--------+
|         256 | 461 ms | 686 ms | 1.23 s |        |        |        |
+-------------+--------+--------+--------+--------+--------+--------+
|         512 | 1.91 s | 2.88 s | 4.22 s |        |        |        |
+-------------+--------+--------+--------+--------+--------+--------+

hermite basis
+-------------+--------+--------+--------+--------+--------+--------+
|   *   sigma | 128    | 256    | 512    | 1024   | 2048   | 4096   |
|       *     |        |        |        |        |        |        |
|   m       * |        |        |        |        |        |        |
+=============+========+========+========+========+========+========+
|           1 |        |        |        | 143 ms | 247 ms | 1.7 s  |
+-------------+--------+--------+--------+--------+--------+--------+
|           4 |        |        |        | 264 ms | 380 ms | 2.57 s |
+-------------+--------+--------+--------+--------+--------+--------+
|          16 |        |        |        | 180 ms | 614 ms | 3.38 s |
+-------------+--------+--------+--------+--------+--------+--------+
|         128 | 259 ms | 185 ms | 745 ms |        |        |        |
+-------------+--------+--------+--------+--------+--------+--------+
|         256 | 332 ms | 700 ms | 762 ms |        |        |        |
+-------------+--------+--------+--------+--------+--------+--------+
|         512 | 1.68 s | 2.05 s | 3.26 s |        |        |        |
+-------------+--------+--------+--------+--------+--------+--------+

GF(243)

matrix multiplication
+-----+---------+
|   n | time    |
+=====+=========+
|  16 | 1.77 ms |
+-----+---------+
|  32 | 14 ms   |
+-----+---------+
|  64 | 111 ms  |
+-----+---------+

uniform profile
+-------------+---------+---------+--------+
|   *   sigma | 16      | 32      | 64     |
|       *     |         |         |        |
|   m       * |         |         |        |
+=============+=========+=========+========+
|           1 | 17.1 ms | 119 ms  | 485 ms |
+-------------+---------+---------+--------+
|           4 | 11.9 ms | 86.7 ms | 387 ms |
+-------------+---------+---------+--------+
|          16 | 6.35 ms | 26.5 ms | 450 ms |
+-------------+---------+---------+--------+
|          32 | 3.83 ms | 16.4 ms | 395 ms |
+-------------+---------+---------+--------+
|          64 | 10.6 ms | 40.8 ms | 226 ms |
+-------------+---------+---------+--------+

hermite profile
+-------------+---------+---------+--------+
|   *   sigma | 16      | 32      | 64     |
|       *     |         |         |        |
|   m       * |         |         |        |
+=============+=========+=========+========+
|           1 | 16.9 ms | 118 ms  | 999 ms |
+-------------+---------+---------+--------+
|           4 | 23 ms   | 162 ms  | 777 ms |
+-------------+---------+---------+--------+
|          16 | 28.5 ms | 97.2 ms | 1.65 s |
+-------------+---------+---------+--------+
|          32 | 14.6 ms | 220 ms  | 1.14 s |
+-------------+---------+---------+--------+
|          64 | 33.4 ms | 226 ms  | 1.87 s |
+-------------+---------+---------+--------+

uniform basis
+-------------+---------+--------+--------+
|   *   sigma | 16      | 32     | 64     |
|       *     |         |        |        |
|   m       * |         |        |        |
+=============+=========+========+========+
|           1 | 21.5 ms | 111 ms | 975 ms |
+-------------+---------+--------+--------+
|           4 | 25.6 ms | 117 ms | 883 ms |
+-------------+---------+--------+--------+
|          16 | 36 ms   | 110 ms | 701 ms |
+-------------+---------+--------+--------+
|          32 | 50.8 ms | 110 ms | 318 ms |
+-------------+---------+--------+--------+
|          64 | 66.3 ms | 139 ms | 554 ms |
+-------------+---------+--------+--------+

hermite basis
+-------------+---------+---------+--------+
|   *   sigma | 16      | 32      | 64     |
|       *     |         |         |        |
|   m       * |         |         |        |
+=============+=========+=========+========+
|           1 | 20.4 ms | 58.1 ms | 1.08 s |
+-------------+---------+---------+--------+
|           4 | 36.3 ms | 191 ms  | 1.3 s  |
+-------------+---------+---------+--------+
|          16 | 59.3 ms | 256 ms  | 1.01 s |
+-------------+---------+---------+--------+
|          32 | 72.2 ms | 287 ms  | 1.98 s |
+-------------+---------+---------+--------+
|          64 | 90.3 ms | 154 ms  | 1.44 s |
+-------------+---------+---------+--------+

GF(256)

matrix multiplication
+------+---------+
|    n | time    |
+======+=========+
|   64 | 558 μs  |
+------+---------+
|  128 | 1.52 ms |
+------+---------+
|  256 | 5.45 ms |
+------+---------+
|  512 | 16.2 ms |
+------+---------+
| 1024 | 36.9 ms |
+------+---------+

uniform profile
+-------------+---------+---------+--------+--------+
|   *   sigma | 128     | 256     | 512    | 1024   |
|       *     |         |         |        |        |
|   m       * |         |         |        |        |
+=============+=========+=========+========+========+
|           1 |         | 61.9 ms | 515 ms | 1.99 s |
+-------------+---------+---------+--------+--------+
|           4 |         | 56.7 ms | 471 ms | 2.09 s |
+-------------+---------+---------+--------+--------+
|          16 |         | 104 ms  | 421 ms | 1.84 s |
+-------------+---------+---------+--------+--------+
|         128 | 77.6 ms | 179 ms  | 544 ms |        |
+-------------+---------+---------+--------+--------+
|         256 | 126 ms  | 172 ms  | 751 ms |        |
+-------------+---------+---------+--------+--------+
|         512 | 236 ms  | 499 ms  | 1.07 s |        |
+-------------+---------+---------+--------+--------+

hermite profile
+-------------+--------+---------+--------+--------+
|   *   sigma | 128    | 256     | 512    | 1024   |
|       *     |        |         |        |        |
|   m       * |        |         |        |        |
+=============+========+=========+========+========+
|           1 |        | 59.9 ms | 519 ms | 2.37 s |
+-------------+--------+---------+--------+--------+
|           4 |        | 83.2 ms | 797 ms | 3.62 s |
+-------------+--------+---------+--------+--------+
|          16 |        | 220 ms  | 1.09 s | 3.87 s |
+-------------+--------+---------+--------+--------+
|         128 | 145 ms | 390 ms  | 1.72 s |        |
+-------------+--------+---------+--------+--------+
|         256 | 133 ms | 537 ms  | 2.18 s |        |
+-------------+--------+---------+--------+--------+
|         512 | 343 ms | 396 ms  | 2.76 s |        |
+-------------+--------+---------+--------+--------+

uniform basis
+-------------+---------+---------+--------+--------+
|   *   sigma | 64      | 128     | 256    | 512    |
|       *     |         |         |        |        |
|   m       * |         |         |        |        |
+=============+=========+=========+========+========+
|           1 |         | 34.9 ms | 208 ms | 750 ms |
+-------------+---------+---------+--------+--------+
|           4 |         | 98.3 ms | 218 ms | 720 ms |
+-------------+---------+---------+--------+--------+
|          16 |         | 111 ms  | 215 ms | 319 ms |
+-------------+---------+---------+--------+--------+
|          64 | 93.4 ms | 71.5 ms | 293 ms |        |
+-------------+---------+---------+--------+--------+
|         128 | 134 ms  | 101 ms  | 457 ms |        |
+-------------+---------+---------+--------+--------+
|         256 | 273 ms  | 408 ms  | 635 ms |        |
+-------------+---------+---------+--------+--------+

hermite basis
+-------------+--------+---------+--------+--------+
|   *   sigma | 64     | 128     | 256    | 512    |
|       *     |        |         |        |        |
|   m       * |        |         |        |        |
+=============+========+=========+========+========+
|           1 |        | 72.8 ms | 227 ms | 360 ms |
+-------------+--------+---------+--------+--------+
|           4 |        | 130 ms  | 335 ms | 1.08 s |
+-------------+--------+---------+--------+--------+
|          16 |        | 200 ms  | 585 ms | 2.47 s |
+-------------+--------+---------+--------+--------+
|          64 | 167 ms | 184 ms  | 805 ms |        |
+-------------+--------+---------+--------+--------+
|         128 | 121 ms | 686 ms  | 2.85 s |        |
+-------------+--------+---------+--------+--------+
|         256 | 529 ms | 1.28 s  | 3.95 s |        |
+-------------+--------+---------+--------+--------+

GF(257)

matrix multiplication
+------+---------+
|    n | time    |
+======+=========+
|  128 | 532 μs  |
+------+---------+
|  256 | 2.39 ms |
+------+---------+
|  512 | 12.4 ms |
+------+---------+
| 1024 | 73.4 ms |
+------+---------+

uniform profile
+-------------+---------+---------+--------+--------+
|   *   sigma | 128     | 256     | 512    | 1024   |
|       *     |         |         |        |        |
|   m       * |         |         |        |        |
+=============+=========+=========+========+========+
|           1 |         | 70 ms   | 308 ms | 1.62 s |
+-------------+---------+---------+--------+--------+
|           4 |         | 72.8 ms | 296 ms | 1.55 s |
+-------------+---------+---------+--------+--------+
|          16 |         | 78 ms   | 286 ms | 1.44 s |
+-------------+---------+---------+--------+--------+
|         128 | 71.3 ms | 169 ms  | 459 ms |        |
+-------------+---------+---------+--------+--------+
|         256 | 127 ms  | 288 ms  | 684 ms |        |
+-------------+---------+---------+--------+--------+
|         512 | 232 ms  | 506 ms  | 857 ms |        |
+-------------+---------+---------+--------+--------+

hermite profile
+-------------+--------+---------+--------+--------+
|   *   sigma | 128    | 256     | 512    | 1024   |
|       *     |        |         |        |        |
|   m       * |        |         |        |        |
+=============+========+=========+========+========+
|           1 |        | 70.6 ms | 302 ms | 1.73 s |
+-------------+--------+---------+--------+--------+
|           4 |        | 114 ms  | 477 ms | 2.48 s |
+-------------+--------+---------+--------+--------+
|          16 |        | 160 ms  | 641 ms | 3.31 s |
+-------------+--------+---------+--------+--------+
|         128 | 119 ms | 325 ms  | 1.09 s |        |
+-------------+--------+---------+--------+--------+
|         256 | 188 ms | 472 ms  | 1.01 s |        |
+-------------+--------+---------+--------+--------+
|         512 | 322 ms | 757 ms  | 2.01 s |        |
+-------------+--------+---------+--------+--------+

uniform basis
+-------------+--------+---------+--------+--------+
|   *   sigma | 128    | 256     | 512    | 1024   |
|       *     |        |         |        |        |
|   m       * |        |         |        |        |
+=============+========+=========+========+========+
|           1 |        | 97.1 ms | 200 ms | 2.19 s |
+-------------+--------+---------+--------+--------+
|           4 |        | 94 ms   | 390 ms | 2.08 s |
+-------------+--------+---------+--------+--------+
|          16 |        | 108 ms  | 399 ms | 2.04 s |
+-------------+--------+---------+--------+--------+
|         128 | 337 ms | 480 ms  | 874 ms |        |
+-------------+--------+---------+--------+--------+
|         256 | 343 ms | 1.32 s  | 2.11 s |        |
+-------------+--------+---------+--------+--------+
|         512 | 1.98 s | 2.24 s  | 4.72 s |        |
+-------------+--------+---------+--------+--------+

hermite basis
+-------------+--------+--------+--------+--------+
|   *   sigma | 128    | 256    | 512    | 1024   |
|       *     |        |        |        |        |
|   m       * |        |        |        |        |
+=============+========+========+========+========+
|           1 |        | 94 ms  | 206 ms | 2.18 s |
+-------------+--------+--------+--------+--------+
|           4 |        | 135 ms | 557 ms | 3.04 s |
+-------------+--------+--------+--------+--------+
|          16 |        | 195 ms | 739 ms | 3.88 s |
+-------------+--------+--------+--------+--------+
|         128 | 191 ms | 432 ms | 1.24 s |        |
+-------------+--------+--------+--------+--------+
|         256 | 414 ms | 408 ms | 1.85 s |        |
+-------------+--------+--------+--------+--------+
|         512 | 1.2 s  | 1.71 s | 3.29 s |        |
+-------------+--------+--------+--------+--------+

GF(59049)

matrix multiplication
+-----+---------+
|   n | time    |
+=====+=========+
|  16 | 3.37 ms |
+-----+---------+
|  32 | 34.2 ms |
+-----+---------+
|  64 | 283 ms  |
+-----+---------+

uniform profile
+-------------+---------+---------+--------+
|   *   sigma | 16      | 32      | 64     |
|       *     |         |         |        |
|   m       * |         |         |        |
+=============+=========+=========+========+
|           1 | 37.5 ms | 286 ms  | 2.82 s |
+-------------+---------+---------+--------+
|           4 | 34 ms   | 239 ms  | 2.55 s |
+-------------+---------+---------+--------+
|          16 | 8.33 ms | 119 ms  | 1.03 s |
+-------------+---------+---------+--------+
|          32 | 12.9 ms | 87.3 ms | 757 ms |
+-------------+---------+---------+--------+
|          64 | 21.5 ms | 94.8 ms | 542 ms |
+-------------+---------+---------+--------+

hermite profile
+-------------+---------+--------+--------+
|   *   sigma | 16      | 32     | 64     |
|       *     |         |        |        |
|   m       * |         |        |        |
+=============+=========+========+========+
|           1 | 49.8 ms | 290 ms | 3.23 s |
+-------------+---------+--------+--------+
|           4 | 67.9 ms | 495 ms | 2.74 s |
+-------------+---------+--------+--------+
|          16 | 43 ms   | 406 ms | 3.38 s |
+-------------+---------+--------+--------+
|          32 | 48.9 ms | 421 ms | 3.49 s |
+-------------+---------+--------+--------+
|          64 | 45.5 ms | 431 ms | 4.1 s  |
+-------------+---------+--------+--------+

uniform basis
+-------------+---------+--------+--------+
|   *   sigma | 16      | 32     | 64     |
|       *     |         |        |        |
|   m       * |         |        |        |
+=============+=========+========+========+
|           1 | 27.2 ms | 241 ms | 2.23 s |
+-------------+---------+--------+--------+
|           4 | 52.6 ms | 295 ms | 2.24 s |
+-------------+---------+--------+--------+
|          16 | 101 ms  | 357 ms | 2.24 s |
+-------------+---------+--------+--------+
|          32 | 180 ms  | 488 ms | 2.37 s |
+-------------+---------+--------+--------+
|          64 | 336 ms  | 809 ms | 2.77 s |
+-------------+---------+--------+--------+

hermite basis
+-------------+---------+----------+--------+
|   *   sigma | 16      | 32       | 64     |
|       *     |         |          |        |
|   m       * |         |          |        |
+=============+=========+==========+========+
|           1 | 26.8 ms | 268 ms   | 2.07 s |
+-------------+---------+----------+--------+
|           4 | 78.6 ms | 477 ms   | 4 s    |
+-------------+---------+----------+--------+
|          16 | 146 ms  | 775 ms   | 5.56 s |
+-------------+---------+----------+--------+
|          32 | 224 ms  | 1e+03 ms | 6.65 s |
+-------------+---------+----------+--------+
|          64 | 434 ms  | 1.32 s   | 6.91 s |
+-------------+---------+----------+--------+

GF(65536)

matrix multiplication
+-----+---------+
|   n | time    |
+=====+=========+
|  16 | 13.4 ms |
+-----+---------+
|  32 | 37.1 ms |
+-----+---------+
|  64 | 129 ms  |
+-----+---------+

uniform profile
+-------------+---------+--------+--------+
|   *   sigma | 16      | 32     | 64     |
|       *     |         |        |        |
|   m       * |         |        |        |
+=============+=========+========+========+
|           1 | 140 ms  | 426 ms | 1.91 s |
+-------------+---------+--------+--------+
|           4 | 93.1 ms | 355 ms | 1.45 s |
+-------------+---------+--------+--------+
|          16 | 51.5 ms | 222 ms | 930 ms |
+-------------+---------+--------+--------+
|          32 | 57.9 ms | 153 ms | 700 ms |
+-------------+---------+--------+--------+
|          64 | 70.3 ms | 163 ms | 471 ms |
+-------------+---------+--------+--------+

hermite profile
+-------------+--------+--------+--------+
|   *   sigma | 16     | 32     | 64     |
|       *     |        |        |        |
|   m       * |        |        |        |
+=============+========+========+========+
|           1 | 114 ms | 455 ms | 1.89 s |
+-------------+--------+--------+--------+
|           4 | 187 ms | 642 ms | 2.25 s |
+-------------+--------+--------+--------+
|          16 | 213 ms | 689 ms | 2.44 s |
+-------------+--------+--------+--------+
|          32 | 230 ms | 738 ms | 2.64 s |
+-------------+--------+--------+--------+
|          64 | 233 ms | 743 ms | 2.79 s |
+-------------+--------+--------+--------+

uniform basis
+-------------+--------+--------+--------+
|   *   sigma | 16     | 32     | 64     |
|       *     |        |        |        |
|   m       * |        |        |        |
+=============+========+========+========+
|           1 | 147 ms | 462 ms | 1.62 s |
+-------------+--------+--------+--------+
|           4 | 141 ms | 408 ms | 1.51 s |
+-------------+--------+--------+--------+
|          16 | 194 ms | 548 ms | 1.77 s |
+-------------+--------+--------+--------+
|          32 | 311 ms | 652 ms | 1.96 s |
+-------------+--------+--------+--------+
|          64 | 499 ms | 1.16 s | 2.84 s |
+-------------+--------+--------+--------+

hermite basis
+-------------+--------+--------+--------+
|   *   sigma | 16     | 32     | 64     |
|       *     |        |        |        |
|   m       * |        |        |        |
+=============+========+========+========+
|           1 | 147 ms | 490 ms | 1.62 s |
+-------------+--------+--------+--------+
|           4 | 216 ms | 637 ms | 2.26 s |
+-------------+--------+--------+--------+
|          16 | 340 ms | 973 ms | 3.19 s |
+-------------+--------+--------+--------+
|          32 | 442 ms | 1.19 s | 3.84 s |
+-------------+--------+--------+--------+
|          64 | 801 ms | 1.76 s | 5.51 s |
+-------------+--------+--------+--------+

GF(65537)

matrix multiplication
+------+---------+
|    n | time    |
+======+=========+
|   64 | 203 μs  |
+------+---------+
|  128 | 689 μs  |
+------+---------+
|  256 | 3.48 ms |
+------+---------+
|  512 | 19.6 ms |
+------+---------+
| 1024 | 123 ms  |
+------+---------+

uniform profile
+-------------+---------+---------+--------+--------+
|   *   sigma | 128     | 256     | 512    | 1024   |
|       *     |         |         |        |        |
|   m       * |         |         |        |        |
+=============+=========+=========+========+========+
|           1 |         | 86.1 ms | 547 ms | 2.65 s |
+-------------+---------+---------+--------+--------+
|           4 |         | 104 ms  | 408 ms | 2.36 s |
+-------------+---------+---------+--------+--------+
|          16 |         | 92.4 ms | 383 ms | 2.06 s |
+-------------+---------+---------+--------+--------+
|         128 | 84.6 ms | 191 ms  | 532 ms |        |
+-------------+---------+---------+--------+--------+
|         256 | 136 ms  | 315 ms  | 777 ms |        |
+-------------+---------+---------+--------+--------+
|         512 | 248 ms  | 556 ms  | 1.26 s |        |
+-------------+---------+---------+--------+--------+

hermite profile
+-------------+--------+---------+--------+--------+
|   *   sigma | 128    | 256     | 512    | 1024   |
|       *     |        |         |        |        |
|   m       * |        |         |        |        |
+=============+========+=========+========+========+
|           1 |        | 98.6 ms | 478 ms | 2.64 s |
+-------------+--------+---------+--------+--------+
|           4 |        | 147 ms  | 646 ms | 3.51 s |
+-------------+--------+---------+--------+--------+
|          16 |        | 193 ms  | 810 ms | 4.57 s |
+-------------+--------+---------+--------+--------+
|         128 | 135 ms | 378 ms  | 1.34 s |        |
+-------------+--------+---------+--------+--------+
|         256 | 209 ms | 571 ms  | 1.7 s  |        |
+-------------+--------+---------+--------+--------+
|         512 | 355 ms | 847 ms  | 2.33 s |        |
+-------------+--------+---------+--------+--------+

uniform basis
+-------------+---------+--------+--------+--------+--------+
|   *   sigma | 64      | 128    | 256    | 512    | 1024   |
|       *     |         |        |        |        |        |
|   m       * |         |        |        |        |        |
+=============+=========+========+========+========+========+
|           1 |         |        | 150 ms | 672 ms | 3.63 s |
+-------------+---------+--------+--------+--------+--------+
|           4 |         |        | 149 ms | 673 ms | 3.35 s |
+-------------+---------+--------+--------+--------+--------+
|          16 |         |        | 150 ms | 612 ms | 2.91 s |
+-------------+---------+--------+--------+--------+--------+
|          64 | 96.4 ms | 153 ms | 282 ms |        |        |
+-------------+---------+--------+--------+--------+--------+
|         128 | 207 ms  | 370 ms | 561 ms |        |        |
+-------------+---------+--------+--------+--------+--------+
|         256 | 496 ms  | 801 ms | 1.46 s |        |        |
+-------------+---------+--------+--------+--------+--------+

hermite basis
+-------------+---------+--------+--------+--------+--------+
|   *   sigma | 64      | 128    | 256    | 512    | 1024   |
|       *     |         |        |        |        |        |
|   m       * |         |        |        |        |        |
+=============+=========+========+========+========+========+
|           1 |         |        | 151 ms | 667 ms | 3.61 s |
+-------------+---------+--------+--------+--------+--------+
|           4 |         |        | 198 ms | 855 ms | 4.25 s |
+-------------+---------+--------+--------+--------+--------+
|          16 |         |        | 238 ms | 982 ms | 5.22 s |
+-------------+---------+--------+--------+--------+--------+
|          64 | 64.3 ms | 132 ms | 364 ms |        |        |
+-------------+---------+--------+--------+--------+--------+
|         128 | 128 ms  | 224 ms | 520 ms |        |        |
+-------------+---------+--------+--------+--------+--------+
|         256 | 332 ms  | 468 ms | 878 ms |        |        |
+-------------+---------+--------+--------+--------+--------+

GF(67108879)

matrix multiplication
+-----+---------+
|   n | time    |
+=====+=========+
|  64 | 447 μs  |
+-----+---------+
| 128 | 3.1 ms  |
+-----+---------+
| 256 | 27.4 ms |
+-----+---------+

uniform profile
+-------------+---------+---------+--------+
|   *   sigma | 64      | 128     | 256    |
|       *     |         |         |        |
|   m       * |         |         |        |
+=============+=========+=========+========+
|           1 | 12.5 ms | 50.2 ms | 384 ms |
+-------------+---------+---------+--------+
|           4 | 11.7 ms | 45 ms   | 333 ms |
+-------------+---------+---------+--------+
|          16 | 12.5 ms | 42.8 ms | 287 ms |
+-------------+---------+---------+--------+
|          64 | 20.5 ms | 55.2 ms | 275 ms |
+-------------+---------+---------+--------+
|         128 | 37 ms   | 78.3 ms | 316 ms |
+-------------+---------+---------+--------+
|         256 | 61.2 ms | 135 ms  | 394 ms |
+-------------+---------+---------+--------+

hermite profile
+-------------+---------+---------+--------+
|   *   sigma | 64      | 128     | 256    |
|       *     |         |         |        |
|   m       * |         |         |        |
+=============+=========+=========+========+
|           1 | 12.4 ms | 50.3 ms | 383 ms |
+-------------+---------+---------+--------+
|           4 | 18.4 ms | 73.9 ms | 632 ms |
+-------------+---------+---------+--------+
|          16 | 26.2 ms | 102 ms  | 812 ms |
+-------------+---------+---------+--------+
|          64 | 42 ms   | 145 ms  | 1.05 s |
+-------------+---------+---------+--------+
|         128 | 59.5 ms | 193 ms  | 1.2 s  |
+-------------+---------+---------+--------+
|         256 | 92.2 ms | 263 ms  | 1.44 s |
+-------------+---------+---------+--------+

uniform basis
+-------------+---------+---------+--------+
|   *   sigma | 64      | 128     | 256    |
|       *     |         |         |        |
|   m       * |         |         |        |
+=============+=========+=========+========+
|           1 | 16.8 ms | 66.5 ms | 500 ms |
+-------------+---------+---------+--------+
|           4 | 16.9 ms | 63 ms   | 440 ms |
+-------------+---------+---------+--------+
|          16 | 23.6 ms | 68.8 ms | 410 ms |
+-------------+---------+---------+--------+
|          64 | 96.6 ms | 161 ms  | 510 ms |
+-------------+---------+---------+--------+
|         128 | 209 ms  | 383 ms  | 792 ms |
+-------------+---------+---------+--------+
|         256 | 496 ms  | 823 ms  | 1 s    |
+-------------+---------+---------+--------+

hermite basis
+-------------+---------+---------+--------+
|   *   sigma | 64      | 128     | 256    |
|       *     |         |         |        |
|   m       * |         |         |        |
+=============+=========+=========+========+
|           1 | 16.8 ms | 66.3 ms | 502 ms |
+-------------+---------+---------+--------+
|           4 | 23.3 ms | 91.1 ms | 686 ms |
+-------------+---------+---------+--------+
|          16 | 34.4 ms | 124 ms  | 898 ms |
+-------------+---------+---------+--------+
|          64 | 70.5 ms | 197 ms  | 1.26 s |
+-------------+---------+---------+--------+
|         128 | 134 ms  | 296 ms  | 1.5 s  |
+-------------+---------+---------+--------+
|         256 | 329 ms  | 575 ms  | 1.93 s |
+-------------+---------+---------+--------+

[xcaruso]

Just a few comments on the doctests to start with.

[xcaruso]

Is it really a method we want to expose to the user? (It's really a question, I don't know.)
If you choose to keep it, you probably need to explain briefly what it is and why it is important (with references).

[vneiger]

It gives some information (monomial basis, vector space dimension...) but this information can also be found from the output of the main function (currently called linear interpolation basis). And calling that main function should rarely be slower than calling this one, so my vote would go against exposing it.

[vneiger]

Also, I would suggest permuting the order of arguments: J, shift, degree, output_pairs.

[vneiger]

It took 1 hour to change my mind :-) . This function is computing not only the row rank profile but also the corresponding rows of the Krylov matrix (output pivot_matrix). As such, it is computing what is sometimes referred to as a Krylov basis for the row vectors of self and the multiplication by J, with order of vectors defined by the shift. This is instrumental in computations for various matrix transformations and properties (see e.g. Chapter 9 of this PhD thesis and this paper ). We could keep this available to the user, but name it more explicitly such as krylov_basis, and possibly slightly rework the interface (input/output/options).

[xcaruso]

Sounds like a good solution!

[vneiger]

"striped Krylov matrix" is not so classical. Suggestion:
"the striped Krylov matrix built from the rows of self and the matrix J, with rows ordered according to a priority defined by shift. See [Beckermann and Labahn, 2000, https://doi.org/10.1137/S0895479897326912] and [Jeannerod, Neiger, Schost, Villard, 2017, https://doi.org/10.1016/j.jsc.2016.11.015]."

Note that we might add these to the SageMath list of references that go into the documentation, but I don't know if that is necessary here (assuming this becomes an auxiliary function hidden from the user).

[vneiger]

Suggestion: using (here and everywhere) a more transparent name for this matrix J. E.g. M or mulmat, shorthand for "multiplication matrix".

[xcaruso]

I vote for M.

[vneiger]

Indicate that the number of integers should equal the number of rows of self.

[vneiger]

Indicate that this is optional, and if not provided, a suitable maximum degree will be taken, namely the number of columns c of J (higher powers of J are known not to be needed since there exists an annihilating polynomial of degree c for J).

[vneiger]

Add some note about the input output_pairs. Since in the main function this optional argument is set to True, maybe this could be the default? Or, maybe better (?), the optional argument could be output_row_indices, set to False by default (which would correspond to the current output_pairs=True), and when set to True this would correspond to the current output_pairs=False.

[Biffo89]

I've removed the modifications to matrix_modn_dense_template.pxi. I played around with modifying the code and although it can almost be made to work, it looks like reconstructing the row rank profile of the matrix from the outputs of P and Q without rewriting part of an RREF algorithm is near impossible.

P returned from linbox_echelonize seems to be list of indices of the swapped rows at the time they were swapped. Unfortunately, the swapping can disturb the lexicographic order. I ran a test to see if a function mapping P and Q to the desired row rank profile existed, and there is the following inconsistency in GF(3):

[1 0 1 0]
[1 0 0 0]
[1 0 0 0]
[0 1 0 0]

[1 0 1 0]
[1 0 1 0]
[1 0 0 0]
[0 1 0 0]

Both return P = [0, 3, 2] and Q = [0, 1, 2], while the row rank profile for the first is [0,1,3] and [0,2,3] for the second. For illustration the value of P seems to be calculated as follows for the first matrix:

[1 0 1 0]
[1 0 0 0]
[1 0 0 0]
[0 1 0 0]

The first row with a non-zero entry in the first column is index 0. So P = [0], and we update the matrix to:

[1 0 1 0]
[0 0 2 0]
[0 0 2 0]
[0 1 0 0]

The first remaining row with a non-zero entry in the second column is index 3. So P = [0, 3], and we update the matrix to:

[1 0 1 0]
[0 1 0 0]
[0 0 2 0]
[0 0 2 0]

The first remaining row with a non-zero entry in the third column is index 2. So P = [0, 3, 2], and we update the matrix to:

[1 0 0 0]
[0 1 0 0]
[0 0 1 0]
[0 0 0 0]

[vneiger]

Here this reference is maybe not the most relevant, however it could be cited for the algorithm itself, or for the specific definition of striped Krylov matrix with a priority given by shift (the number of the relevant section should be indicated). See this about the ALGORITHM block and also about adding references to the master bibliography file.

For, more generally, the relevant notion of basis of such modules I would rather cite Kailath's book Linear Systems, Section 6.7 (the book is already in Sage's biblio as [Kai1980]) and also Section 16.6 of the book Algebraic Complexity Theory (which is not yet in Sage's biblio).

[vneiger]

Same remark as above concerning notation (changing J into e.g. M or mulmat).

[vneiger]

This name comes from the reference that was used for the algorithm, where interpolation came naturally from its main contribution (which is on a particular case of this problem). Here, in the general case, it is not very clear why this should be seen as an interpolation problem. The method name could be for example krylov_kernel_basis (open to discussion!).

[vneiger]

I've removed the modifications to matrix_modn_dense_template.pxi. I played around with modifying the code and although it can almost be made to work, it looks like reconstructing the row rank profile of the matrix from the outputs of P and Q without rewriting part of an RREF algorithm is near impossible.
[...]

In terms of efficiency, since this has been dropped (it should be if it is not correct, thanks for checking carefully), I assume this means that an additional RREF is computed on the transpose?

LinBox/FFLAS-FFPACK does incorporate echelonization procedures which preserve the row rank profile as well as the column rank profile. We could try to include this in Sage. This would avoid wasting a factor of 2 by calling the echelonization twice when both are needed. Would you mind opening an issue for this? (and you could link to your above message or copy it there, to show the observation that a trivial deduction from P does not always work)

[xcaruso]

I'm not 100% sure but I think that indentation is still incorrect.
The classical format is

A tuple:
- ``row_profile``: list of the first r row indices in the striped
  Krylov matrix ``K`` corresponding to independants rows.
  [...]
- ``pivot_matrix``: the submatrix of ``K`` given by ``row_profile``
- ``column_profile``: the first r independent column indices in
  ``pivot_matrix``

[xcaruso]

Maybe, it could be a good idea to have a method krylov_matrix for constructing this matrix.
I think it makes sense if we also have krylov_basis and krylov_kernel_basis.

By the way, what "striped" means here?

[vneiger]

After inspection (web search), "striped" is used when describing matrices with particular structures, but usually that revolve around constant (anti)diagonals including some that are zero. Which is not really the case here.

Here the term was taken from a less common occurrence: it was used in Beckermann & Labahn (SIAM, 2000) to describe the matrices we have here, and just in a few other documents that mention this B&L's work. "Block Krylov" would probably have been a more classical terminology (maybe this was not chosen due to confusion with the existing block-Krylov / block-Wiedemann method?).

In any case, I think we could avoid this unclear terminology and just talk about a Krylov matrix (it is just a usual Krylov matrix, with several initial vectors; and the presence of several vectors gives rise to several row orderings that we describe through the shift).

[vneiger]

+1 for a method krylov_matrix which builds the matrix. It should take some bound as input: it could be the maximum exponent of the multiplication matrix that is used, this for each of the input rows of self.