From c2d6803400ed6ea514b4a0649e665b8cd557409c Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Wed, 10 Jul 2019 11:26:30 +0200 Subject: [PATCH 01/13] working on kernel basis via approx at large order --- src/sage/matrix/matrix_polynomial_dense.pyx | 164 ++++++++++++++++++++ 1 file changed, 164 insertions(+) diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index 8a6f3089fb1..30b3bb98d4e 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 +- Romain Lebreton and Vincent Neiger (2019-??-??): added functions for + computing and for verifying minimal kernel bases """ # **************************************************************************** # Copyright (C) 2016 Kwankyu Lee @@ -1995,3 +1997,165 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): else: rest_order[j] -= 1 return appbas, rdeg + + + 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: + + TODO (approximation large order + ZLS 12 ?). + + EXAMPLES:: + + sage: pR. = GF(7)[] + """ + m = self.nrows() + n = self.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') + elif (not row_wise) and len(shifts) != n: + raise ValueError('shifts length should be the column dimension') + + # compute approximant basis + # if required, normalize it into shifted Popov form + if row_wise: + P,rdeg = self._approximant_basis_iterative(order, shifts) + if normal_form: + # compute the list "- pivot degree" + # (since weak Popov, pivot degree is rdeg-shifts entrywise) + # Note: -deg(P[i,i]) = shifts[i] - rdeg[i] + degree_shifts = [shifts[i] - rdeg[i] for i in range(m)] + # compute approximant basis with that list as shifts + P,rdeg = self._approximant_basis_iterative(order, + degree_shifts) + # left-multiply by inverse of leading matrix + lmat = P.leading_matrix(shifts=degree_shifts) + P = lmat.inverse() * P + else: + P,rdeg = self.transpose()._approximant_basis_iterative(order, + shifts) + if normal_form: + # compute the list "- pivot degree" + # (since weak Popov, pivot degree is rdeg-shifts entrywise) + degree_shifts = [shifts[i] - rdeg[i] for i in range(n)] + # compute approximant basis with that list as shifts + P,rdeg = self.transpose()._approximant_basis_iterative( \ + order, degree_shifts) + P = P.transpose() + # right-multiply by inverse of leading matrix + lmat = P.leading_matrix(shifts=degree_shifts,row_wise=False) + P = P * lmat.inverse() + else: + P = P.transpose() + + return P + + def _kernel_basis_via_approximant(self, shifts, degree_bounds=None): + r""" + Return a ``shifts``-ordered weak Popov kernel basis for this polynomial + matrix (see :meth:`minimal_kernel_basis` for definitions). The output + basis is considered row-wise, that is, its rows are in the left kernel + of ``self``. + + If ``degree_bounds`` is provided, it must be a list of integers such + that, if ``d`` is the list of column degrees of ``self`` shifted by + ``degree_bounds`` to which we add ``1`` entry-wise, then a + ``shifts``-minimal approximant basis of ``self`` at order ``d`` + contains a ``shifts``-minimal kernel basis of ``self``. If this + argument is not provided, then general bounds are used; namely, we take + ``degree_bounds`` as the list $(c+a,\ldots,c+a)$ where $c$ is the sum + of column degrees of ``self`` and $a$ is the amplitude of ``shifts`` + (i.e. the difference between its maximum and its minimum). + + The input dimensions are supposed to be sound: the length of ``shifts`` + must be the number of rows of ``self``. + + INPUT: + + - ``shifts`` -- a list of integers. + + - ``degree_bounds`` -- (optional, default: ``None``) list of integers, + which must be either ``None`` or a list of bounds on the column + degrees of any `shifts`-minimal kernel basis of ``self``. + + OUTPUT: + + - a polynomial matrix (the kernel basis ``P``). + + - a list of integers (the shifts-row degree of ``P``). + + - a list of integers (the shifts-leading positions of ``P``). + + ALGORITHM: + + This computes a `shifts`-minimal approximant basis at sufficiently + large order so that a subset of the rows of this basis form a kernel + basis. This subset is identified using degree properties. + + EXAMPLES:: + + sage: pR. = GF(7)[] + """ + m = self.nrows() + n = self.ncols() + d = self.degree() + + if d is -1: # matrix is zero + return Matrix.identity(self.base_ring(), m, m), shifts, [i for i in range(m)] + + if m <= n and self(0).rank() == m: # early exit: kernel is empty + return Matrix.zero(self.base_ring(), 0, m), [], [] + + if degree_bounds==None: + # list of column degrees with zero columns discarded + cdeg = self.column_degrees() + cdeg = [cdeg[i] for i in range(n) if cdeg[i] >= 0] + # take general bound: bound on pivot degree + amplitude of shift + degree_bound = sum(cdeg)+max(shifts)-min(shifts) + degree_bounds = [degree_bound]*m + + orders = self.column_degrees(degree_bounds) + for i in range(n): orders[i] = orders[i]+1 + + P = self.minimal_approximant_basis(orders,shifts) + return P From 6d6cf1ff769cd362b9175e30c180309824a8bbf2 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Mon, 13 Jan 2020 12:56:11 +0100 Subject: [PATCH 02/13] minor diff --- src/sage/matrix/matrix_polynomial_dense.pyx | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index 30b3bb98d4e..4016cb966c5 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -2059,28 +2059,25 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): # compute approximant basis # if required, normalize it into shifted Popov form if row_wise: - P,rdeg = self._approximant_basis_iterative(order, shifts) + P,rdeg = self._approximant_basis_iterative(shifts) if normal_form: # compute the list "- pivot degree" # (since weak Popov, pivot degree is rdeg-shifts entrywise) # Note: -deg(P[i,i]) = shifts[i] - rdeg[i] degree_shifts = [shifts[i] - rdeg[i] for i in range(m)] # compute approximant basis with that list as shifts - P,rdeg = self._approximant_basis_iterative(order, - degree_shifts) + P,rdeg = self._approximant_basis_iterative(degree_shifts) # left-multiply by inverse of leading matrix lmat = P.leading_matrix(shifts=degree_shifts) P = lmat.inverse() * P else: - P,rdeg = self.transpose()._approximant_basis_iterative(order, - shifts) + P,rdeg = self.transpose()._approximant_basis_iterative(shifts) if normal_form: # compute the list "- pivot degree" # (since weak Popov, pivot degree is rdeg-shifts entrywise) degree_shifts = [shifts[i] - rdeg[i] for i in range(n)] # compute approximant basis with that list as shifts - P,rdeg = self.transpose()._approximant_basis_iterative( \ - order, degree_shifts) + P,rdeg = self.transpose()._approximant_basis_iterative(degree_shifts) P = P.transpose() # right-multiply by inverse of leading matrix lmat = P.leading_matrix(shifts=degree_shifts,row_wise=False) From c04d0a4d0af40927f6876adf59031ffb69d0e2c6 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Mon, 13 Jan 2020 15:56:35 +0100 Subject: [PATCH 03/13] finishing kernel via approximation; inserting two simple examples --- src/sage/matrix/matrix_polynomial_dense.pyx | 86 +++++++-------------- 1 file changed, 27 insertions(+), 59 deletions(-) diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index 4016cb966c5..9cd1140c90f 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -2044,6 +2044,16 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): EXAMPLES:: sage: pR. = GF(7)[] + sage: F = Matrix([[(x+1)*(x+3)],[(x+1)*(x+3)+1]]) + sage: F.minimal_kernel_basis() + [6*x^2 + 3*x + 3 x^2 + 4*x + 3] + + sage: F = Matrix([[(x+1)*(x+3)],[(x+1)*(x+4)]]) + sage: F.minimal_kernel_basis() + [6*x + 3 x + 3] + + sage: F.minimal_kernel_basis(row_wise=False) + [6*x^2 + 3*x + 3 x^2 + 4*x + 3] """ m = self.nrows() n = self.ncols() @@ -2056,54 +2066,22 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): elif (not row_wise) and len(shifts) != n: raise ValueError('shifts length should be the column dimension') - # compute approximant basis - # if required, normalize it into shifted Popov form + # compute kernel basis if row_wise: - P,rdeg = self._approximant_basis_iterative(shifts) - if normal_form: - # compute the list "- pivot degree" - # (since weak Popov, pivot degree is rdeg-shifts entrywise) - # Note: -deg(P[i,i]) = shifts[i] - rdeg[i] - degree_shifts = [shifts[i] - rdeg[i] for i in range(m)] - # compute approximant basis with that list as shifts - P,rdeg = self._approximant_basis_iterative(degree_shifts) - # left-multiply by inverse of leading matrix - lmat = P.leading_matrix(shifts=degree_shifts) - P = lmat.inverse() * P + P = self._kernel_basis_via_approximant(shifts,normal_form) else: - P,rdeg = self.transpose()._approximant_basis_iterative(shifts) - if normal_form: - # compute the list "- pivot degree" - # (since weak Popov, pivot degree is rdeg-shifts entrywise) - degree_shifts = [shifts[i] - rdeg[i] for i in range(n)] - # compute approximant basis with that list as shifts - P,rdeg = self.transpose()._approximant_basis_iterative(degree_shifts) - P = P.transpose() - # right-multiply by inverse of leading matrix - lmat = P.leading_matrix(shifts=degree_shifts,row_wise=False) - P = P * lmat.inverse() - else: - P = P.transpose() + P = self.transpose()._kernel_basis_via_approximant(shifts,normal_form) + print P return P - def _kernel_basis_via_approximant(self, shifts, degree_bounds=None): + def _kernel_basis_via_approximant(self, shifts, normal_form=False): r""" Return a ``shifts``-ordered weak Popov kernel basis for this polynomial matrix (see :meth:`minimal_kernel_basis` for definitions). The output basis is considered row-wise, that is, its rows are in the left kernel of ``self``. - If ``degree_bounds`` is provided, it must be a list of integers such - that, if ``d`` is the list of column degrees of ``self`` shifted by - ``degree_bounds`` to which we add ``1`` entry-wise, then a - ``shifts``-minimal approximant basis of ``self`` at order ``d`` - contains a ``shifts``-minimal kernel basis of ``self``. If this - argument is not provided, then general bounds are used; namely, we take - ``degree_bounds`` as the list $(c+a,\ldots,c+a)$ where $c$ is the sum - of column degrees of ``self`` and $a$ is the amplitude of ``shifts`` - (i.e. the difference between its maximum and its minimum). - The input dimensions are supposed to be sound: the length of ``shifts`` must be the number of rows of ``self``. @@ -2111,48 +2089,38 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): - ``shifts`` -- a list of integers. - - ``degree_bounds`` -- (optional, default: ``None``) list of integers, - which must be either ``None`` or a list of bounds on the column - degrees of any `shifts`-minimal kernel basis of ``self``. + - ``normal_form`` -- (optional, default: ``False``) boolean, if + ``True`` then the output basis is in ``shifts``-Popov form. OUTPUT: - a polynomial matrix (the kernel basis ``P``). - - a list of integers (the shifts-row degree of ``P``). - - - a list of integers (the shifts-leading positions of ``P``). - ALGORITHM: This computes a `shifts`-minimal approximant basis at sufficiently large order so that a subset of the rows of this basis form a kernel basis. This subset is identified using degree properties. - - EXAMPLES:: - - sage: pR. = GF(7)[] """ m = self.nrows() n = self.ncols() d = self.degree() if d is -1: # matrix is zero - return Matrix.identity(self.base_ring(), m, m), shifts, [i for i in range(m)] + return Matrix.identity(self.base_ring(), m, m) if m <= n and self(0).rank() == m: # early exit: kernel is empty - return Matrix.zero(self.base_ring(), 0, m), [], [] + return Matrix(self.base_ring(), 0, m) - if degree_bounds==None: - # list of column degrees with zero columns discarded - cdeg = self.column_degrees() - cdeg = [cdeg[i] for i in range(n) if cdeg[i] >= 0] - # take general bound: bound on pivot degree + amplitude of shift - degree_bound = sum(cdeg)+max(shifts)-min(shifts) - degree_bounds = [degree_bound]*m + degree_bound = min(m,n)*d+max(shifts) + degree_bounds = [degree_bound - shifts[i] for i in range(m)] orders = self.column_degrees(degree_bounds) for i in range(n): orders[i] = orders[i]+1 - P = self.minimal_approximant_basis(orders,shifts) - return P + 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,:] From edddf97c6ec1b50727546d7e1c7e8f1aef20adc0 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Tue, 14 Jan 2020 10:35:30 +0100 Subject: [PATCH 04/13] simplify --- src/sage/matrix/matrix_polynomial_dense.pyx | 83 +++++++++------------ 1 file changed, 37 insertions(+), 46 deletions(-) diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index 9cd1140c90f..4517cf82f2a 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -2057,6 +2057,7 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): """ m = self.nrows() n = self.ncols() + d = self.degree() # set default shifts / check shifts dimension if shifts is None: @@ -2068,59 +2069,49 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): # compute kernel basis if row_wise: - P = self._kernel_basis_via_approximant(shifts,normal_form) - else: - P = self.transpose()._kernel_basis_via_approximant(shifts,normal_form) - print P - - return P - - def _kernel_basis_via_approximant(self, shifts, normal_form=False): - r""" - Return a ``shifts``-ordered weak Popov kernel basis for this polynomial - matrix (see :meth:`minimal_kernel_basis` for definitions). The output - basis is considered row-wise, that is, its rows are in the left kernel - of ``self``. + if d is -1: # matrix is zero + return Matrix.identity(self.base_ring(), m, m) - The input dimensions are supposed to be sound: the length of ``shifts`` - must be the number of rows of ``self``. + if m <= n and self(0).rank() == m: # early exit: kernel is empty + return Matrix(self.base_ring(), 0, m) - INPUT: - - - ``shifts`` -- a list of integers. + # 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)] - - ``normal_form`` -- (optional, default: ``False``) boolean, if - ``True`` then the output basis is in ``shifts``-Popov form. + # orders for approximation + orders = self.column_degrees(degree_bounds) + for i in range(n): orders[i] = orders[i]+1 - OUTPUT: - - - a polynomial matrix (the kernel basis ``P``). - - ALGORITHM: + # 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,:] - This computes a `shifts`-minimal approximant basis at sufficiently - large order so that a subset of the rows of this basis form a kernel - basis. This subset is identified using degree properties. - """ - m = self.nrows() - n = self.ncols() - d = self.degree() + else: + if d is -1: # matrix is zero + return Matrix.identity(self.base_ring(), n, n) - if d is -1: # matrix is zero - return Matrix.identity(self.base_ring(), m, m) + if n <= m and self(0).rank() == n: # early exit: kernel is empty + return Matrix(self.base_ring(), n, 0) - if m <= n and self(0).rank() == m: # early exit: kernel is empty - 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(n)] - degree_bound = min(m,n)*d+max(shifts) - degree_bounds = [degree_bound - shifts[i] for i in range(m)] + # orders for approximation + orders = self.row_degrees(degree_bounds) + for i in range(m): orders[i] = orders[i]+1 - orders = self.column_degrees(degree_bounds) - for i in range(n): 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] - 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,:] + return P From 94d6e0a5924ccc3546c9f41196b0346443e46a93 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Wed, 15 Jan 2020 16:05:38 +0100 Subject: [PATCH 05/13] ok --- src/sage/matrix/matrix_polynomial_dense.pyx | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index 4517cf82f2a..5848e2604a3 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -2096,6 +2096,7 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): return Matrix.identity(self.base_ring(), n, n) if n <= m and self(0).rank() == n: # early exit: kernel is empty + print self.base_ring() return Matrix(self.base_ring(), n, 0) # degree bounds on the kernel basis @@ -2113,5 +2114,3 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): if P[j,j].degree() + shifts[j] <= degree_bound: column_indices.append(j) return P[:,column_indices] - - return P From 132afb663b0816a30070a94deb2b173a7d3dfedc Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Wed, 1 Apr 2020 17:23:21 +0200 Subject: [PATCH 06/13] starting is-kernel-basis: doc text --- src/sage/matrix/matrix_polynomial_dense.pyx | 48 +++++++++++++++++++-- 1 file changed, 45 insertions(+), 3 deletions(-) diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index 5848e2604a3..c460af280e2 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -16,8 +16,8 @@ AUTHORS: - Vincent Neiger (2018-09-29): added functions for computing and for verifying minimal approximant bases -- Romain Lebreton and Vincent Neiger (2019-??-??): added functions for - computing and for verifying minimal kernel bases +- Vincent Neiger (2020-04-01): added functions for computing and for verifying + minimal kernel bases """ # **************************************************************************** # Copyright (C) 2016 Kwankyu Lee @@ -1719,7 +1719,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``. @@ -1998,6 +1998,48 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): 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:: + """ + return True def minimal_kernel_basis(self, shifts=None, From a9b6abb3ae74bc7f55b3dafc0ff6a162301cb841 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Wed, 1 Apr 2020 17:33:54 +0200 Subject: [PATCH 07/13] is-kernel-basis: some tests --- src/sage/matrix/matrix_polynomial_dense.pyx | 40 +++++++++++++++++++-- 1 file changed, 38 insertions(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index c460af280e2..492cbad1f2c 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -2030,7 +2030,7 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): ALGORITHM: - Verification that the matrix (has full rank and) is in shifted weak + 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 @@ -2039,7 +2039,43 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): EXAMPLES:: """ - return 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: + 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 + + # check that self consists of kernel vectors + if row_wise and self * pmat != 0: + return False + elif (not row_wise) and pmat * self != 0: + return False + + # check rank TODO + # check saturated TODO def minimal_kernel_basis(self, shifts=None, From 19958708c110464b62d8d1dc8b321c663d25eccb Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Wed, 1 Apr 2020 20:46:58 +0200 Subject: [PATCH 08/13] is-kernel-basis: is in kernel + rank is right --- src/sage/matrix/matrix_polynomial_dense.pyx | 33 +++++++++++---------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index 492cbad1f2c..e6ade6e3245 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -1617,12 +1617,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: @@ -2061,21 +2059,26 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): + " of the input matrix") # check full rank 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 consists of kernel vectors if row_wise and self * pmat != 0: return False - elif (not row_wise) and pmat * self != 0: + 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 - # check rank TODO - # check saturated TODO + # check saturated def minimal_kernel_basis(self, shifts=None, From c674d69803d9536087c43b0245b8612eb3ae723c Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Thu, 2 Apr 2020 00:31:02 +0200 Subject: [PATCH 09/13] is-kernel-basis: is saturated --- src/sage/matrix/matrix_polynomial_dense.pyx | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index e6ade6e3245..ec318e790ac 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -1675,7 +1675,6 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): return True - def minimal_approximant_basis(self, order, shifts=None, @@ -2078,7 +2077,20 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): if (not row_wise) and self.ncols()==n-rk: return False - # check saturated + # 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, From dc440148ada9fba38f917bdf4addeb01f9a25c85 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Thu, 2 Apr 2020 01:40:35 +0200 Subject: [PATCH 10/13] fix equality testing --- src/sage/matrix/matrix_polynomial_dense.pyx | 24 +++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index ec318e790ac..85c577fb5cd 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -2035,6 +2035,11 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): 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 """ m = pmat.nrows() n = pmat.ncols() @@ -2063,20 +2068,26 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): if (not normal_form) and (not self.is_weak_popov(shifts, row_wise, True, False)): return False + #print "full rank / owP --> ok" + # 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 + #print "in kernel --> ok" + # 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: + if row_wise and self.nrows()!=m-rk: return False - if (not row_wise) and self.ncols()==n-rk: + if (not row_wise) and self.ncols()!=n-rk: return False + #print "good rank --> ok" + # 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 @@ -2090,6 +2101,8 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): if hnf != 1: return False + #print "saturated --> ok" + return True def minimal_kernel_basis(self, @@ -2130,9 +2143,9 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): OUTPUT: a polynomial matrix. - ALGORITHM: - - TODO (approximation large order + ZLS 12 ?). + ALGORITHM: uses minimal approximant basis computation at a + sufficiently large order so that the approximant basis contains + a kernel basis as a submatrix. EXAMPLES:: @@ -2189,7 +2202,6 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): return Matrix.identity(self.base_ring(), n, n) if n <= m and self(0).rank() == n: # early exit: kernel is empty - print self.base_ring() return Matrix(self.base_ring(), n, 0) # degree bounds on the kernel basis From 3d1097690e7004fcefd2b2f54e4993d2d108cdd9 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Thu, 2 Apr 2020 01:55:37 +0200 Subject: [PATCH 11/13] some examples in doc --- src/sage/matrix/matrix_polynomial_dense.pyx | 32 ++++++++++++--------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index 85c577fb5cd..20483cd7446 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -2037,9 +2037,17 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): 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 """ m = pmat.nrows() n = pmat.ncols() @@ -2068,16 +2076,12 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): if (not normal_form) and (not self.is_weak_popov(shifts, row_wise, True, False)): return False - #print "full rank / owP --> ok" - # 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 - #print "in kernel --> ok" - # 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() @@ -2086,8 +2090,6 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): if (not row_wise) and self.ncols()!=n-rk: return False - #print "good rank --> ok" - # 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 @@ -2101,8 +2103,6 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): if hnf != 1: return False - #print "saturated --> ok" - return True def minimal_kernel_basis(self, @@ -2150,16 +2150,22 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): EXAMPLES:: sage: pR. = GF(7)[] - sage: F = Matrix([[(x+1)*(x+3)],[(x+1)*(x+3)+1]]) - sage: F.minimal_kernel_basis() + 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: F = Matrix([[(x+1)*(x+3)],[(x+1)*(x+4)]]) - sage: F.minimal_kernel_basis() + sage: pmat = Matrix([[(x+1)*(x+3)],[(x+1)*(x+4)]]) + sage: pmat.minimal_kernel_basis() [6*x + 3 x + 3] - sage: F.minimal_kernel_basis(row_wise=False) + sage: pmat.minimal_kernel_basis(row_wise=False) [6*x^2 + 3*x + 3 x^2 + 4*x + 3] + + sage: pmat = Matrix(pR, [[1,x,x**2]]) + sage: pmat.minimal_kernel_basis(row_wise=False) + [x 0] + [6 x] + [0 6] """ m = self.nrows() n = self.ncols() From 32bfc3b14d89a1446819f875a57d5ec110949086 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Fri, 3 Apr 2020 13:43:22 +0200 Subject: [PATCH 12/13] fix small bug + add some examples --- src/sage/matrix/matrix_polynomial_dense.pyx | 19 +++++++++++++++++-- 1 file changed, 17 insertions(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index 20483cd7446..e11cea5e47e 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -2048,6 +2048,12 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): 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`); 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() @@ -2159,13 +2165,18 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): [6*x + 3 x + 3] sage: pmat.minimal_kernel_basis(row_wise=False) - [6*x^2 + 3*x + 3 x^2 + 4*x + 3] + [] sage: pmat = Matrix(pR, [[1,x,x**2]]) - sage: pmat.minimal_kernel_basis(row_wise=False) + 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() @@ -2182,9 +2193,11 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): # 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 @@ -2205,9 +2218,11 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): 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 From 0c9bca03be079e1a904d33d3e1587375a3a50751 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Fri, 3 Apr 2020 13:54:03 +0200 Subject: [PATCH 13/13] minor doc fixes --- src/sage/matrix/matrix_polynomial_dense.pyx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index e11cea5e47e..24eb3e299e5 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -2035,6 +2035,7 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): https://arxiv.org/pdf/1807.01272.pdf for details). EXAMPLES:: + sage: pR. = GF(97)[] sage: pmat = Matrix(pR, [[1],[x],[x**2]]) @@ -2049,7 +2050,7 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): sage: kerbas.is_minimal_kernel_basis(pmat) False - Shifts and right kernel bases are supported (with `row_wise`); one can test whether the kernel basis is normalized in shifted-Popov form (with `normal_form`):: + 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])