Add functions for computation of krylov basis and krylov kernel 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.

[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)

[Biffo89]

I've made many of the changes in comments, but realised my original naive implementation during testing had performance problems and with the necessary fixes should be the implementation for smaller cases. I'm currently working on mapping the parameter criteria as the relation is not trivial.

[xcaruso]

Some additional comments; however, my review is not finished (I will continue later).

[xcaruso]

I'm not sure we want to change the base ring of self.
I think I would prefer:

    M = M.change_ring(E.base_ring())

but I'm open to the discussion if you think that it is not a good idea.

[Biffo89]

Are we happy that in the case E has a smaller base ring than M, for example E is over GF(7) and M is over GF(49) that the operation should fail (given M actually has elements not in GF(7))? Currently, E would be coerced to GF(49) in this case.

[xcaruso]

OK.
I think you cannot use it here because you have extra arguments but let me mention the decorator @coerce_binop which basically does what you are doing.

[xcaruso]

We should update these lines, I think.

[xcaruso]

Same remark as above.

[xcaruso]

Some of them should be EXAMPLES.

[Biffo89]

I have one uncertainty about what expected behaviour should be. Should the degree bounds impact the returned row indices that refer to the Krylov matrix? Currently placing different degree bounds for each row changes the returned row indices for certain inputs, such as in the following:

sage: E
[27 30 37 81 63 68 94 44 35 25]
[ 0  0  0  0  0  0  0  0  0  1]
sage: M
[0 1 0 0 0 0 0 0 0 0]
[0 0 1 0 0 0 0 0 0 0]
[0 0 0 1 0 0 0 0 0 0]
[0 0 0 0 1 0 0 0 0 0]
[0 0 0 0 0 1 0 0 0 0]
[0 0 0 0 0 0 1 0 0 0]
[0 0 0 0 0 0 0 1 0 0]
[0 0 0 0 0 0 0 0 1 0]
[0 0 0 0 0 0 0 0 0 1]
[0 0 0 0 0 0 0 0 0 0]
sage: E.krylov_basis(M,shifts=[5,0])
(
[ 0  0  0  0  0  0  0  0  0  1]
[27 30 37 81 63 68 94 44 35 25]
[ 0 27 30 37 81 63 68 94 44 35]
[ 0  0 27 30 37 81 63 68 94 44]
[ 0  0  0 27 30 37 81 63 68 94]
[ 0  0  0  0 27 30 37 81 63 68]
[ 0  0  0  0  0 27 30 37 81 63]
[ 0  0  0  0  0  0 27 30 37 81]
[ 0  0  0  0  0  0  0 27 30 37]
[ 0  0  0  0  0  0  0  0 27 30], ((1, 0, 0), (0, 0, 5), (0, 1, 7), (0, 2, 9), (0, 3, 11), (0, 4, 13), (0, 5, 15), (0, 6, 17), (0, 7, 18), (0, 8, 19))
)
sage: E.krylov_basis(M,shifts=[5,0],degree=[10,0])
(
[ 0  0  0  0  0  0  0  0  0  1]
[27 30 37 81 63 68 94 44 35 25]
[ 0 27 30 37 81 63 68 94 44 35]
[ 0  0 27 30 37 81 63 68 94 44]
[ 0  0  0 27 30 37 81 63 68 94]
[ 0  0  0  0 27 30 37 81 63 68]
[ 0  0  0  0  0 27 30 37 81 63]
[ 0  0  0  0  0  0 27 30 37 81]
[ 0  0  0  0  0  0  0 27 30 37]
[ 0  0  0  0  0  0  0  0 27 30], ((1, 0, 0), (0, 0, 1), (0, 1, 2), (0, 2, 3), (0, 3, 4), (0, 4, 5), (0, 5, 6), (0, 6, 7), (0, 7, 8), (0, 8, 9))
)

This is because the naive method computes a krylov matrix only up to the degrees specified. (The elimination method is adjusted to return degrees consistent with this.)

sage: print(E.krylov_matrix(M,shifts=[5,0]))
[ 0  0  0  0  0  0  0  0  0  1]
[ 0  0  0  0  0  0  0  0  0  0]
[ 0  0  0  0  0  0  0  0  0  0]
[ 0  0  0  0  0  0  0  0  0  0]
[ 0  0  0  0  0  0  0  0  0  0]
[27 30 37 81 63 68 94 44 35 25]
[ 0  0  0  0  0  0  0  0  0  0]
[ 0 27 30 37 81 63 68 94 44 35]
[ 0  0  0  0  0  0  0  0  0  0]
[ 0  0 27 30 37 81 63 68 94 44]
[ 0  0  0  0  0  0  0  0  0  0]
[ 0  0  0 27 30 37 81 63 68 94]
[ 0  0  0  0  0  0  0  0  0  0]
[ 0  0  0  0 27 30 37 81 63 68]
[ 0  0  0  0  0  0  0  0  0  0]
[ 0  0  0  0  0 27 30 37 81 63]
[ 0  0  0  0  0  0  0  0  0  0]
[ 0  0  0  0  0  0 27 30 37 81]
[ 0  0  0  0  0  0  0 27 30 37]
[ 0  0  0  0  0  0  0  0 27 30]
[ 0  0  0  0  0  0  0  0  0 27]
[ 0  0  0  0  0  0  0  0  0  0]
sage: print(E.krylov_matrix(M,shifts=[5,0],degree=[10,0]))
[ 0  0  0  0  0  0  0  0  0  1]
[27 30 37 81 63 68 94 44 35 25]
[ 0 27 30 37 81 63 68 94 44 35]
[ 0  0 27 30 37 81 63 68 94 44]
[ 0  0  0 27 30 37 81 63 68 94]
[ 0  0  0  0 27 30 37 81 63 68]
[ 0  0  0  0  0 27 30 37 81 63]
[ 0  0  0  0  0  0 27 30 37 81]
[ 0  0  0  0  0  0  0 27 30 37]
[ 0  0  0  0  0  0  0  0 27 30]
[ 0  0  0  0  0  0  0  0  0 27]
[ 0  0  0  0  0  0  0  0  0  0]

It doesn't sit completely right with me, having an argument that is provided for performance, affect the value returned. I was considering adjusting the returned row indices.

However, even choosing a uniform bound causes inconsistencies:

sage: E
[31 50 95 22  3 75 65 57 15 26]
[37 10 73 28 62 29  3 89 34 15]
[ 0  0  0  0  0  0  0  0  0  1]
sage: M
[0 1 0 0 0 0 0 0 0 0]
[0 0 1 0 0 0 0 0 0 0]
[0 0 0 1 0 0 0 0 0 0]
[0 0 0 0 1 0 0 0 0 0]
[0 0 0 0 0 1 0 0 0 0]
[0 0 0 0 0 0 1 0 0 0]
[0 0 0 0 0 0 0 1 0 0]
[0 0 0 0 0 0 0 0 1 0]
[0 0 0 0 0 0 0 0 0 1]
[0 0 0 0 0 0 0 0 0 0]
sage: E.krylov_basis(M,shifts=[5,5,0])
(
[ 0  0  0  0  0  0  0  0  0  1]
[31 50 95 22  3 75 65 57 15 26]
[37 10 73 28 62 29  3 89 34 15]
[ 0 31 50 95 22  3 75 65 57 15]
[ 0 37 10 73 28 62 29  3 89 34]
[ 0  0 31 50 95 22  3 75 65 57]
[ 0  0 37 10 73 28 62 29  3 89]
[ 0  0  0 31 50 95 22  3 75 65]
[ 0  0  0 37 10 73 28 62 29  3]
[ 0  0  0  0 31 50 95 22  3 75], ((2, 0, 0), (0, 0, 5), (1, 0, 6), (0, 1, 8), (1, 1, 9), (0, 2, 11), (1, 2, 12), (0, 3, 14), (1, 3, 15), (0, 4, 17))
)
sage: E.krylov_basis(M,shifts=[5,5,0],degree=4)
(
[ 0  0  0  0  0  0  0  0  0  1]
[31 50 95 22  3 75 65 57 15 26]
[37 10 73 28 62 29  3 89 34 15]
[ 0 31 50 95 22  3 75 65 57 15]
[ 0 37 10 73 28 62 29  3 89 34]
[ 0  0 31 50 95 22  3 75 65 57]
[ 0  0 37 10 73 28 62 29  3 89]
[ 0  0  0 31 50 95 22  3 75 65]
[ 0  0  0 37 10 73 28 62 29  3]
[ 0  0  0  0 31 50 95 22  3 75], ((2, 0, 0), (0, 0, 5), (1, 0, 6), (0, 1, 7), (1, 1, 8), (0, 2, 9), (1, 2, 10), (0, 3, 11), (1, 3, 12), (0, 4, 13))
)

An alternative solution would be to give the indices as they would appear in an "infinite" Krylov matrix, i.e. where all degrees of M appear. This is inconsistent with the examples given in https://arxiv.org/abs/1512.03503 however.

[xcaruso]

IMO, the outputted index row should be the index of the row as it appears in the output of krylov_matrix with the same parameters.

[xcaruso]

Some further comments.

[xcaruso]

Is it important that degree is a vector and not just a list?
By the way, don't we need to use the plural for degrees (similarly to shifts)?

[Biffo89]

The main advantage to this is the automatic type conversion. If degree contains 1.0 or 10/5 then ZZ**E.nrows() will automatically check if conversion is possible. A secondary advantage is that this automatically checks that the length of the list, tuple, or vector is correct.

I'll change degree to degrees.

[xcaruso]

I just noticed that there exists a method _check_shift_dimension.
Probably, it is better to use it for shifts and have a similar methods for degrees (it will avoid code duplication). This method currently does not check that the entries are integers but maybe we can add it.

[xcaruso]

Maybe we should say somewhere that if the bound is not correct, then the output may be incorrect as well.

[xcaruso]

This probably should be fixed in the __init__/_element_constructor_ method of the relevant class.
But that's for another PR.

[vneiger]

IMO, the outputted index row should be the index of the row as it appears in the output of krylov_matrix with the same parameters.

Yes, I agree with this. @Biffo89 can you confirm this is the behaviour currently? (it seems to be, from your examples above) If yes, then I would not change anything regarding this point.

It doesn't sit completely right with me, having an argument that is provided for performance, affect the value returned. I was considering adjusting the returned row indices.

We do not have to see this as some inconsistency, and we do not have to see the degree bounds as being here for performance only. To expand on xcaruso's comment: the degree bounds w.r.t each row of E allow us to define the Krylov matrix. It is entirely defined from E, M, the shifts, and the degrees. If you change the degrees, the matrix changes; this method krylov_matrix has interest independently from krylov_basis, and for this one, the degrees are not given only for performance, they define the output. Since krylov_basis is about selecting suitable rows from this Krylov matrix, one would naturally expect the row indices to correspond to the rows in the Krylov matrix as returned by krylov_matrix with the same parameters. These indices may indeed change depending on the degrees (and the shifts) provided as input. If I'm not mistaken, however, the output pairs indicating the rows of E and the power of M for each of these basis rows should be the same whatever the input degrees are (assuming they have been given sufficiently large for correctness).

[xcaruso]

Again, I think that those lines are not needed for a private method.

[xcaruso]

OK.
I think you cannot use it here because you have extra arguments but let me mention the decorator @coerce_binop which basically does what you are doing.

[xcaruso]

What's the difference between row_coordinates and row_profile?

[Biffo89]

row_coordinates[_][:2] is a sorted list of the row-degree pairs (c,d) in K. row_profile is initially the independent rows in K. Line 19523 uses row_coordinates to update row_profile to add the corresponding row in self (row_coordinates[_][0]) and degree of M (row_coordinates[_][1]).

[xcaruso]

Sorry, my comment was very very badly formulated.
I was in fact wondering if you cannot simply write

    row_profile = tuple([row_coordinates[i] for i in row_profile])

[xcaruso]

This line is very long. For readibility, I suggest to cut it:

Suggested change

	R = matrix.block([[exhausted],[R],[(R*M_L).matrix_from_rows([i for i in range(len(row_profile_self)) if row_profile_self[i][1] + L <= degrees[row_profile_self[i][0]]])]], subdivide=False)
	rows = [i for i, x in enumerate(row_profile_self) if x[1] + L <= degrees[x[0]]]
	S = (R*M_L).matrix_from_rows(rows)
	R = matrix.block([[exhausted],[R],[S]], subdivide=False)

[xcaruso]

Suggested change

	k_coordinates = self._krylov_row_coordinates(shifts, degrees, k)
	k_perm = Permutation([x[2] + 1 for x in k_coordinates])
	k_rows = [x[2] for x in self._krylov_row_coordinates(shifts, degrees, k)]
	k_perm = Permutation([x + 1 for x in k_rows])

so that, instead of writing k_perm(i+1)-1 everywhere afterwards, we can just write k_rows[i].

[xcaruso]

Or maybe _krylov_row_coordinates could have an option for just returning the row index.

[xcaruso]

TypeError or NotImplementedError?