Trac #28143: Minimal kernel bases

New functionalities:

    * computation of shifted minimal kernel bases (e.g. via approximant
bases at large order and/or via Zhou-Labahn-Storjohann 2012),
    * verification that a matrix is a shifted minimal kernel basis.

This should be done in a general context:

    * accepting non-uniform shifts,
    * allowing to work row-wise or column-wise,
    * offering the possibility of obtaining the canonical basis (that
is, the one in shifted Popov form).

URL: https://trac.sagemath.org/28143
Reported by: vneiger
Ticket author(s): Vincent Neiger
Reviewer(s): Romain Lebreton, Pascal Giorgi, Johan Rosenkilde, Seung Gyu
Hyun

src/sage/matrix/matrix_polynomial_dense.pyx

@@ -16,6 +16,8 @@ AUTHORS:
 - Vincent Neiger (2018-09-29): added functions for computing and for verifying
   minimal approximant bases
 
+- Vincent Neiger (2020-04-01): added functions for computing and for verifying
+  minimal kernel bases
 """
 # ****************************************************************************
 #       Copyright (C) 2016 Kwankyu Lee <ekwankyu@gmail.com>
@@ -1612,12 +1614,10 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense):
         if not self.is_square():
             return False
         # check nonsingular and shifts-(ordered weak) Popov form
-        if normal_form:
-            if not self.is_popov(shifts, row_wise, False, False):
-                return False
-        else:
-            if not self.is_weak_popov(shifts, row_wise, True, False):
-                return False
+        if normal_form and (not self.is_popov(shifts, row_wise, False, False)):
+            return False
+        if (not normal_form) and (not self.is_weak_popov(shifts, row_wise, True, False)):
+            return False
 
         # check that self is a basis of the set of approximants
         if row_wise:
@@ -1672,7 +1672,6 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense):
 
         return True
 
-
     def minimal_approximant_basis(self,
             order,
             shifts=None,
@@ -1714,7 +1713,7 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense):
           ``None`` is interpreted as ``shifts=[0,...,0]``.
 
         - ``row_wise`` -- (optional, default: ``True``) boolean, if ``True``
-          then the output basis considered row-wise and operates on the left
+          then the output basis is considered row-wise and operates on the left
           of ``self``; otherwise it is column-wise and operates on the right
           of ``self``.
 
@@ -1992,3 +1991,250 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense):
             else:
                 rest_order[j] -= 1
         return appbas, rdeg
+
+    def is_minimal_kernel_basis(self,
+            pmat,
+            shifts=None,
+            row_wise=True,
+            normal_form=False):
+        r"""
+        Return ``True`` if and only if this matrix is a left kernel basis in
+        ``shifts``-ordered weak Popov form for the polynomial matrix ``pmat``.
+
+        If ``normal_form`` is ``True``, then the kernel basis must furthermore
+        be in ``shifts``-Popov form. An error is raised if the input dimensions
+        are not sound.
+
+        INPUT:
+
+        - ``pmat`` -- a polynomial matrix.
+
+        - ``shifts`` -- (optional, default: ``None``) list of integers;
+          ``None`` is interpreted as ``shifts=[0,...,0]``.
+
+        - ``row_wise`` -- (optional, default: ``True``) boolean, if ``True``
+          then the basis is considered row-wise and operates on the left of
+          ``pmat``; otherwise it is column-wise and operates on the right of
+          ``pmat``.
+
+        - ``normal_form`` -- (optional, default: ``False``) boolean, if
+          ``True`` then checks for a basis in ``shifts``-Popov form.
+
+        OUTPUT: a boolean.
+
+        ALGORITHM:
+
+        Verification that the matrix has full rank and is in shifted weak
+        Popov form is done via :meth:`is_weak_popov`; verification that the
+        matrix is a left kernel basis is done by checking that the rank is
+        correct, that the product is indeed zero, and that the matrix is
+        saturated, i.e. it has unimodular column bases (see Lemma 6.10 of
+        https://arxiv.org/pdf/1807.01272.pdf for details).
+
+        EXAMPLES::
+
+            sage: pR.<x> = GF(97)[]
+            sage: pmat = Matrix(pR, [[1],[x],[x**2]])
+
+            sage: kerbas = Matrix(pR, [[x,-1,0],[0,x,-1]])
+            sage: kerbas.is_minimal_kernel_basis(pmat)
+            True
+
+        A matrix in Popov form which has the right rank, all rows in the
+        kernel, but does not generate the kernel::
+
+            sage: kerbas = Matrix(pR, [[x**2,0,-1],[0,x,-1]])
+            sage: kerbas.is_minimal_kernel_basis(pmat)
+            False
+
+        Shifts and right kernel bases are supported (with ``row_wise``), and one can test whether the kernel basis is normalized in shifted-Popov form (with ``normal_form``)::
+
+            sage: kerbas = Matrix(pR, [[-x,-x**2],[1,0],[0,1]])
+            sage: kerbas.is_minimal_kernel_basis(pmat.transpose(),row_wise=False,normal_form=True,shifts=[0,1,2])
+            True
+        """
+        m = pmat.nrows()
+        n = pmat.ncols()
+
+        # set default shifts / check shifts dimension
+        if shifts is None:
+            shifts = [0] * m if row_wise else [0] * n
+        elif row_wise and len(shifts) != m:
+            raise ValueError('shifts length should be the row dimension of' \
+                                                      + ' the input matrix')
+        elif (not row_wise) and len(shifts) != n:
+            raise ValueError('shifts length should be the column dimension' \
+                                                   + ' of the input matrix')
+
+        # raise an error if self does not have the right dimension
+        if row_wise and self.ncols() != m:
+            raise ValueError("column dimension should be the row dimension" \
+                                                    + " of the input matrix")
+        elif (not row_wise) and self.nrows() != n:
+            raise ValueError("row dimension should be the column dimension" \
+                                                    + " of the input matrix")
+
+        # check full rank and shifts-(ordered weak) Popov form
+        if normal_form and (not self.is_popov(shifts, row_wise, False, False)):
+            return False
+        if (not normal_form) and (not self.is_weak_popov(shifts, row_wise, True, False)):
+            return False
+
+        # check that self consists of kernel vectors
+        if row_wise and self * pmat != 0:
+            return False
+        if (not row_wise) and pmat * self != 0:
+            return False
+
+        # check self.rank() is right (the above weak Popov test ensures self
+        # has full row rank if row wise, and full column rank if column wise)
+        rk = pmat.rank()
+        if row_wise and self.nrows()!=m-rk:
+            return False
+        if (not row_wise) and self.ncols()!=n-rk:
+            return False
+
+        # final check: self is row saturated (assuming row wise),
+        # since self has full rank this is equivalent to the fact that its
+        # columns generate the identity matrix of size m; in particular the
+        # column Popov and column Hermite forms of self are the identity
+        if row_wise:
+            hnf = self.transpose().hermite_form(include_zero_rows=False)
+            # TODO replace by popov_form (likely more efficient) once it is written
+        else:
+            hnf = self.hermite_form(include_zero_rows=False)
+            # TODO replace by popov_form (likely more efficient) once it is written
+        if hnf != 1:
+            return False
+
+        return True
+
+    def minimal_kernel_basis(self,
+            shifts=None,
+            row_wise=True,
+            normal_form=False):
+        r"""
+        Return a left kernel basis in ``shifts``-ordered weak Popov form for
+        this polynomial matrix.
+
+        Assuming we work row-wise, if `F` is an `m \times n` polynomial matrix,
+        then a left kernel basis for `F` is a polynomial matrix whose rows form
+        a basis of the left kernel of `F`, which is the module of polynomial
+        vectors `p` of size `m` such that `p F` is zero.
+
+        If ``normal_form`` is ``True``, then the output basis `P` is
+        furthermore in ``shifts``-Popov form. By default, `P` is considered
+        row-wise, that is, its rows are left kernel vectors for ``self``; if
+        ``row_wise`` is ``False`` then its columns are right kernel vectors for
+        ``self``.
+
+        An error is raised if the input dimensions are not sound: if working
+        row-wise (resp. column-wise), the length of ``shifts`` must be the
+        number of rows (resp. columns) of ``self``.
+
+        INPUT:
+
+        - ``shifts`` -- (optional, default: ``None``) list of integers;
+          ``None`` is interpreted as ``shifts=[0,...,0]``.
+
+        - ``row_wise`` -- (optional, default: ``True``) boolean, if ``True``
+          then the output basis considered row-wise and operates on the left
+          of ``self``; otherwise it is column-wise and operates on the right
+          of ``self``.
+
+        - ``normal_form`` -- (optional, default: ``False``) boolean, if
+          ``True`` then the output basis is in ``shifts``-Popov form.
+
+        OUTPUT: a polynomial matrix.
+
+        ALGORITHM: uses minimal approximant basis computation at a
+        sufficiently large order so that the approximant basis contains
+        a kernel basis as a submatrix.
+
+        EXAMPLES::
+
+            sage: pR.<x> = GF(7)[]
+            sage: pmat = Matrix([[(x+1)*(x+3)],[(x+1)*(x+3)+1]])
+            sage: pmat.minimal_kernel_basis()
+            [6*x^2 + 3*x + 3   x^2 + 4*x + 3]
+
+            sage: pmat = Matrix([[(x+1)*(x+3)],[(x+1)*(x+4)]])
+            sage: pmat.minimal_kernel_basis()
+            [6*x + 3   x + 3]
+
+            sage: pmat.minimal_kernel_basis(row_wise=False)
+            []
+
+            sage: pmat = Matrix(pR, [[1,x,x**2]])
+            sage: pmat.minimal_kernel_basis(row_wise=False,normal_form=True)
+            [x 0]
+            [6 x]
+            [0 6]
+
+            sage: pmat.minimal_kernel_basis(row_wise=False,normal_form=True,shifts=[0,1,2])
+            [  6*x 6*x^2]
+            [    1     0]
+            [    0     1]
+        """
+        m = self.nrows()
+        n = self.ncols()
+        d = self.degree()
+
+        # set default shifts / check shifts dimension
+        if shifts is None:
+            shifts = [0] * m if row_wise else [0] * n
+        elif row_wise and len(shifts) != m:
+            raise ValueError('shifts length should be the row dimension')
+        elif (not row_wise) and len(shifts) != n:
+            raise ValueError('shifts length should be the column dimension')
+
+        # compute kernel basis
+        if row_wise:
+            if d is -1: # matrix is zero
+                from sage.matrix.constructor import Matrix
+                return Matrix.identity(self.base_ring(), m, m)
+
+            if m <= n and self(0).rank() == m: # early exit: kernel is empty
+                from sage.matrix.constructor import Matrix
+                return Matrix(self.base_ring(), 0, m)
+
+            # degree bounds on the kernel basis
+            degree_bound = min(m,n)*d+max(shifts)
+            degree_bounds = [degree_bound - shifts[i] for i in range(m)]
+
+            # orders for approximation
+            orders = self.column_degrees(degree_bounds)
+            for i in range(n): orders[i] = orders[i]+1
+
+            # compute approximant basis and retrieve kernel rows
+            P = self.minimal_approximant_basis(orders,shifts,True,normal_form)
+            row_indices = []
+            for i in range(m):
+                if P[i,i].degree() + shifts[i] <= degree_bound:
+                    row_indices.append(i)
+            return P[row_indices,:]
+
+        else:
+            if d is -1: # matrix is zero
+                from sage.matrix.constructor import Matrix
+                return Matrix.identity(self.base_ring(), n, n)
+
+            if n <= m and self(0).rank() == n: # early exit: kernel is empty
+                from sage.matrix.constructor import Matrix
+                return Matrix(self.base_ring(), n, 0)
+
+            # degree bounds on the kernel basis
+            degree_bound = min(m,n)*d+max(shifts)
+            degree_bounds = [degree_bound - shifts[i] for i in range(n)]
+
+            # orders for approximation
+            orders = self.row_degrees(degree_bounds)
+            for i in range(m): orders[i] = orders[i]+1
+
+            # compute approximant basis and retrieve kernel columns
+            P = self.minimal_approximant_basis(orders,shifts,False,normal_form)
+            column_indices = []
+            for j in range(n):
+                if P[j,j].degree() + shifts[j] <= degree_bound:
+                    column_indices.append(j)
+            return P[:,column_indices]