Wrong kernel dimension nmod_poly_mat_kernel_zls_approx

[gilvillard]

In nmod_poly_mat_kernel_zls_approx called by nmod_poly_mat_kernel

nmod_poly_mat_kernel(NN, pivind, shift, A, ORD_WEAK_POPOV, COL_UPPER);

See two examples below. Actually, Example 1 is created recursively from Example 2 by the implementation.

After an initial inspection with Example 1
Its seems that the easy degree-based detection if (shift[i] < order + diff_shift) detects only 4 kernel vectors, then, after the early termination test if (nz >= m - n), which is 7-3 >= 4 the algorithm terminates.

Hint (?): order and diff_shift have to be checked.

Example 1. Dimensions nxm = 3x7, the dimension of the kernel is 5, 4 is returned

prime = 2;
m= 3;
n= 7;
nmod_poly_mat_init(A, m, n , prime);
nmod_poly_set_str(nmod_poly_mat_entry(A,0,0),"1 2  1");
nmod_poly_set_str(nmod_poly_mat_entry(A,1,1),"1 2  1");
slong shift[n];
shift[0]=9;
shift[1]=9;
shift[2]=10;
shift[3]=8;
shift[4]=8;
shift[5]=8;
shift[6]=8;

nmod_poly_mat_kernel(NN, pivind, shift, A, ORD_WEAK_POPOV, COL_UPPER);

Example 2. Dimensions nxm = 12x31, the dimension of the kernel is 22, 20 is returned

prime= 2;
m= 12;
n= 31;
nmod_poly_mat_init(A, m, n , prime);
nmod_poly_set_str(nmod_poly_mat_entry(A,0,0),"5 2  0 1 0 0 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,0,11),"8 2  0 0 0 0 0 0 0 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,3,2),"7 2  0 1 0 0 0 0 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,3,8),"1 2  1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,3,15),"4 2  0 0 0 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,3,17),"8 2  1 0 0 1 0 0 1 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,3,23),"7 2  1 0 1 1 1 1 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,4,7),"2 2  0 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,4,17),"8 2  1 0 0 1 1 1 1 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,5,13),"3 2  0 1 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,5,22),"4 2  1 1 1 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,5,24),"6 2  1 0 1 0 0 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,5,29),"1 2  1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,6,1),"6 2  1 0 0 0 1 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,6,5),"5 2  1 1 1 1 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,6,9),"1 2  1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,6,12),"8 2  0 1 0 1 1 1 0 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,7,12),"6 2  1 1 0 1 1 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,7,21),"7 2  0 0 0 0 0 0 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,8,20),"5 2  1 0 0 0 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,8,24),"4 2  1 1 1 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,8,26),"1 2  1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,9,2),"8 2  0 1 0 1 0 1 1 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,9,13),"7 2  0 0 0 0 0 0 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,9,20),"1 2  1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,9,23),"5 2  1 1 0 1 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,9,25),"9 2  0 0 0 1 1 0 0 0 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,10,3),"3 2  1 0 1 ");
nmod_poly_set_str(nmod_poly_mat_entry(A,10,22),"1 2  1 ");

nmod_poly_mat_kernel(NN, pivind, NULL, A, ORD_WEAK_POPOV, COL_UPPER);

Configuration
Apple M2 macOS Tahoe clang version 17.0.0

[vneiger]

Thanks for reporting! Chasing the source of this bug from the first example, it seems to be due to the early exit, which is too optimistically early when the sought kernel basis contains pivot entries that are constants. The test checks for sum of pivot degrees, which obviously does not take good account of degree 0 pivots... it seems to me that the early exit is correct except for possibly missing kernel basis vectors with constant pivot (and even then, only for input matrices that are rank deficient), but let me know if you think the bug is deeper than this.

Just to check that a better handling of zero columns/rows would not suffice to resolve this, here is another failing example. For the left kernel of the 5 x 2 matrix with all entries 1, the algorithm returns nullity 3 instead of 4.

#include <pml/nmod_poly_mat_extra.h>

int main(int argc, char *argv[])
{
    slong prime = 997;
    slong m = 5;
    slong n = 2;
    nmod_poly_mat_t A;
    nmod_poly_mat_init(A, m, n, prime);
    for (slong i = 0; i < m; i++)
        for (slong j = 0; j < n; j++)
            nmod_poly_one(nmod_poly_mat_entry(A, i, j));

    slong shift[] = {9, 10, 8, 8, 8};
    slong pivind[m];

    nmod_poly_mat_t K;
    nmod_poly_mat_init(K, m, m, prime);
    slong nz = nmod_poly_mat_kernel_zls_approx(K, pivind, shift, A);
    flint_printf("nullity: %ld\n", nz);
    nmod_poly_mat_clear(A);
    nmod_poly_mat_clear(K);
}

[gilvillard]

I don't know if it's limited to constant pivot or not, I'll keep working on it too. I'm a little unsure because I tried to add a +1 to one of the tests that precedes the early exit (I was suspicious of FLINT_MAX(0, d)), but on the 12x31 example the nullity only increased by 1 (not 2 as expected). Could also be related only to zero columns/rows obtained during algorithm execution using recursion.

[vneiger]

Here are some arguments, trying to make things clearer with respect to two things: the definition of diff_shift, and the design of some correct early exit.

Let's call s the shift shift before the call to pmbasis (currently around line 350), and let t be the same shift but to which we subtract diff_shift to all entries. The choice of diff_shift is such that the entries of t are at least the row degree of pmat, indeed,
t[i] = s[i] - diff_shift >= s[i] - (s[i] - max(0,rdeg(pmat[i,:]))) >= rdeg(pmat[i,:]))
and in fact, diff_shift is the largest value such that this is true, or in other words, it minimizes t while ensuring t >= rdeg(pmat) (i.e., any larger diff_shift would lead to some entry of t being strictly less than the corresponding entry of rdeg(pmat)).

This pmbasis call computes ker which is an s-ordered weak Popov approximant basis of pmat at order order. It also in-place modifies shift into the s-row degree of ker. Since s is t translated by +diff_shift, then ker is also in t-ordered weak Popov form, and the modified shift contains the t-row degree of ker translated by +diff_shift, that is, after the pmbasis call,
shift[i] = rdeg_s(ker[i,:]) = rdeg_t(ker[i,:]) + diff_shift
from which we can deduce that
shift[i] >= rdeg(ker[i,:] * pmat) + diff_shift

Then, what is said to be an "easy degree-based detection of rows in the kernel" in the code is based on the fact that any i such that shift[i] < order + diff_shift corresponds to a row in the kernel of pmat. Indeed this inequality is equivalent to the fact that the t-row degree of the row ker[i,:] is strictly less than order, and since t >= rdeg(pmat) entry-wise, this means that the row degree of ker[i,:] * pmat is strictly less than order. But at the same time this vector is zero modulo x**order, so it must be actually zero, meaning that ker[i,:] is in the kernel of pmat.

This is in good part a mix of classical arguments concerning kernel bases via approximant bases, but hopefully this (technically) highlights the interplay between diff_shift and the row degree of pmat. Another classical thing is that every such row is actually minimal, in the sense that there is no vector in the kernel of pmat which would have s-pivot at index i and whose s-pivot degree would be strictly less than the s-pivot degree deg(ker[i,i]) of the i-th row of the approximant basis ker.

Concerning the early exit, the quantity sum_pmatdeg is the sum of degrees of rows of pmat that are not indexed by some i that are already known to be an s-pivot index of some kernel row (based on the above detection of kernel rows). So this subset of rows includes all rows of pmat that are not indexed by the s-pivot index of some row in the kernel, and this means that sum_pmatdeg is a (non-strict) upper bound on the sum of s-pivot degrees of an s-weak Popov kernel basis. In particular, if the sum of s-pivot degrees of already found kernel rows is equal to sum_pmatdeg, then the only possibly missing kernel rows must have s-pivot degree equal to zero.

[gilvillard]

Your argument seems convincing to me (with the convention degree zero and not -1 for zero if I understand well).

[vneiger]

Your argument seems convincing to me (with the convention degree zero and not -1 for zero if I understand well).

Yes, exactly. This was a bit unclear in the code, with a mix of $-1$ and sometimes of max(0, degree).

then the only possibly missing kernel rows must have s-pivot degree equal to zero.

Expanding upon this, we can say a bit more: "the corresponding row of pmat must be constant as well".

Through pmbasis, we have found indices $i_0 < \cdots < i_{k-1}$ of rows that are in the kernel of pmat (where $k$ is the partial nullity computed in the code: number of already found kernel rows). Because they are extracted from an s-ordered weak Popov basis, their s-pivots are on the diagonal, so at the column positions $i_0 < \cdots < i_{k-1}$. Let $\delta_0,\ldots,\delta_{k-1}$ be their s-pivot degrees. Let $j_0 < \cdots < j_{m-k-1}$ be the $m-k$ indices that form the complement, so that the union of the $i$'s and the $j$'s gives all row indices ${0,\ldots,m-1}$ of the input (or column indices of the kernel). Finally, let $r_0,\ldots,r_{m-k-1}$ be the degrees of the $m-k$ rows indexed by these $j$'s in the input matrix pmat, taken as 0 when the row is zero.

What we know is that $\delta_0 + \cdots + \delta_{k-1} \le r_0 + \cdots + r_{m-k-1}$. This is possibly a strict inequality, even when we have found the whole kernel. Now suppose this is an equality. If we are missing any row from the kernel, say a row with s-pivot index $\pi$, adding this row to the collection of "already known kernel rows" makes both the left-hand side increase (or stay the same) and the right-hand side decrease (or stay the same). But the inequality must remain valid, so, since we already had equality, both sides must stay the same. What this means is that this kernel row must have s-pivot degree 0, and also, the corresponding row $\pi$ of pmat must have degree 0 (or be zero).

This seems sufficient to make the existing early exit strategy correct without impacting its efficiency: to make sure we can early exit, we just have to add a check that, among the rows not detected to be in the kernel, none has pivot degree $0$ with also the corresponding row of pmat having degree $0$.

[vneiger]

This is fixed by #27 . I will also copy the above notes as a comment in the source file, in case questions arise again about the early exit.

Please reopen if the issue persists on your side. Thanks again for reporting!