diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index c70ea5b9b4e..7348c332d79 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/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 @@ -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. = 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. = 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]