Add minimal relation bases for univariate polynomial matrices

[vneiger]

Fixes #23645

This covers:

• the general problem of minimal relation bases (modular matrix equations), see Approximant bases and Relation bases for polynomial matrices #23645 .
• the particular case of minimal interpolant bases (also known as M-Padé approximation), with specific code for better performance.

The focus of this PR is on providing methods that are as versatile/general as possible (while keeping performance at a reasonable level).

With regards to #23645 , note that approximant bases had already been integrated in SageMath some years ago. Thus with this PR all types of relation bases mentioned in that issue are now covered.

[vneiger]

Still to do: add interpolant bases, and the corresponding documentation and tests.

For reference: this is based directly on a minimal kernel basis computation, with suitable degree constraints to ensure that the sought relation basis appears as a submatrix of this kernel basis.

[github-actions]

Documentation preview for this PR (built with commit 26316fc; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

[vneiger]

This is ready for review.

Note that approximant and interpolation bases could call relation bases with an input of some specific form. Yet having dedicated implementations does make a non-negligible difference in performance, see the timings below (ratios in columns rel/app and int/app):

shift leading to reduced form (weak Popov):
m	n	d	appbas	relbas	rel/app	intbas	relbas	rel/int
6	1	200	5.2e-03	8.7e-03	1.7e+00	7.0e-03	9.9e-03	1.4e+00
6	3	200	1.6e-02	5.4e-02	3.3e+00	2.7e-02	8.1e-02	3.0e+00
6	5	200	2.4e-02	1.5e-01	6.4e+00	5.7e-02	2.7e-01	4.8e+00
14	1	20	1.5e-03	2.2e-03	1.5e+00	1.7e-03	2.3e-03	1.4e+00
14	7	20	9.1e-03	4.7e-02	5.1e+00	1.4e-02	5.7e-02	4.0e+00
14	13	20	9.3e-03	1.7e-01	1.8e+01	3.2e-02	2.4e-01	7.7e+00

shift leading to lower triangular basis (close to HNF):
m	n	d	appbas	relbas	rel/app	intbas	relbas	rel/int
6	1	200	5.4e-03	1.6e-01	2.9e+01	8.2e-03	2.1e-01	2.5e+01
6	3	200	1.7e-02	5.9e-01	3.5e+01	3.7e-02	1.8e+00	4.8e+01
6	5	200	2.2e-02	1.2e+00	5.2e+01	8.5e-02	5.6e+00	6.5e+01
14	1	20	1.6e-03	1.0e-01	6.4e+01	1.6e-03	1.1e-01	7.1e+01
14	7	20	8.7e-03	9.5e-01	1.1e+02	1.6e-02	2.5e+00	1.5e+02
14	13	20	9.0e-03	2.3e+00	2.5e+02	4.3e-02	7.2e+00	1.7e+02

[HugoPasse]

I saw in the reference manual that the methods minimal_kernel_basis and minimal_approximant_basis both have their counterparts is_minimal_kernel_basis and is_minimal_approximant_basis. Do you plan to impement equivalent methods for the suggested features minimal_relation_basis and minimal_interpolation_basis?

[vneiger]

I saw in the reference manual that the methods minimal_kernel_basis and minimal_approximant_basis both have their counterparts is_minimal_kernel_basis and is_minimal_approximant_basis. Do you plan to impement equivalent methods for the suggested features minimal_relation_basis and minimal_interpolation_basis?

Thanks for the suggestion. Here are some thoughts on this question.

• Concerning a relation basis P in general, say for a matrix F modulo a matrix M. Testing that P has the right weak Popov form is very cheap. Testing that P consists of relations is already less cheap (one matrix multiplication P F and one matrix division with remainder PF rem M; not clear to me how to do it faster deterministically). And finally, the worst part: for verifying that P indeed generates the module of relations, I'm not sure to see any (deterministic) method that would be faster than basically recomputing P. Well, any method that provides the degree of determinant of this module would directly yield a test, but same thing, not clear to me how to do this much more efficiently than by computing P. With all of this in mind, it seems that such a method would not bring much.

• Concerning an interpolant basis P, say for a matrix F at points $a_{i,j}, 0 \le i < d_j, 0 \le j < n$. Testing weak Popov form remains easy. Testing that P consists of interpolants is somewhat easier than above: the matrix division with remainder is replaced by usual polynomial quo_rem's (or, one might use multipoint evaluation but that is not yet fast in SageMath, I think). However here, the generation test could become easier; at least that is the case when all points are pairwise distinct. I will think about it, to see if there is an "easy" test for more general points, in which case it could make sense to add is_minimal_interpolant_basis.

[vneiger]

here, the generation test could become easier; at least that is the case when all points are pairwise distinct. I will think about it, to see if there is an "easy" test for more general points, in which case it could make sense to add is_minimal_interpolant_basis.

Concerning this last point, to complete my answer: the method minimal_interpolant_basis allows for possibly repeated points, and in that case, I do not see a way to check deterministically that the computed basis P generates the interpolation module, without doing computations that (in the current state of SageMath at least) take about as much time as computing P itself via minimal_interpolant_basis.

So, to answer the initial question, I do not see a strong point in favor of adding the is_... methods for the moment.

[HugoPasse]

I did some tests and found a small issue. If you take $M \in \mathbb{K}[x]^{n \times n}$ non-singular and $F \in \mathbb{K}[x]^{0 \times n}$, one would expect the call F.minimal_relation_basis(M) to return the empty matrix (this is for instance the case when calling F.minimal_kernel_basis()). Currently, an error is raised. It is caused by a call to _right_quo_rem_reduced() that also calls column_degrees(), which raises a ValueError on the empty matrix

Steps to reproduce:

pring.<x> = GF(9001)[]
M = random_matrix(pring,5,5) #This matrix has full rank with high probability
F = random_matrix(pring,0,5)
F.minimal_relation_basis(M) #This yields an error

[vneiger]

Thanks for the comment. I have done some experiments and this is a more general problem among the available methods for these matrices. This would bring this PR too far from its original goal, so I propose to leave this as another task, cf issue #39587 created for this.

[HugoPasse]

LGTM