From e13a26b4e2293fc28a407865abfebfba13cc70a6 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Fri, 13 Jun 2025 08:13:42 +0200 Subject: [PATCH 01/28] implementation --- src/sage/matrix/matrix2.pyx | 343 ++++++++++++++++++++ src/sage/matrix/matrix_polynomial_dense.pyx | 41 +++ 2 files changed, 384 insertions(+) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index a5cfab6aa02..15cdcbf7b90 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -18939,6 +18939,347 @@ cdef class Matrix(Matrix1): U.rescale_col(n - 1, -1) return U + def polynomial_compression(self,degree,variable): + """ + Returns the corresponding polynomial matrix where the coefficient matrix of degree i is the ith block of ``self``. + + INPUT: + + - ``degree`` -- integer: the number of blocks in the expanded matrix. If ``self`` does not have zero-blocks, + this corresponds to the degree of the resulting matrix. + - ``variable`` -- a symbolic variable (e.g., ``x``) to use in the polynomial. + + OUTPUT: + + - a polynomial matrix with the same number of rows and self.ncols()//(degree+1) columns + + EXAMPLES: + + The first 3 x 3 block of N corresponds to the + + sage: R. = GF(97)[] + sage: S. = ZZ[] + sage: M = matrix([[x^2+36*x,31*x,0],[3*x+13,x+57,0],[96,96,1]]) + sage: M + [x^2 + 36*x 31*x 0] + [ 3*x + 13 x + 57 0] + [ 96 96 1] + sage: N = M.expansion() + sage: N + [ 0 0 0 36 31 0 1 0 0] + [13 57 0 3 1 0 0 0 0] + [96 96 1 0 0 0 0 0 0] + sage: N.polynomial_compression(2,x) + [x^2 + 36*x 31*x 0] + [ 3*x + 13 x + 57 0] + [ 96 96 1] + sage: N.polynomial_compression(2,y) + [y^2 + 36*y 31*y 0] + [ 3*y + 13 y + 57 0] + [ 96 96 1] + """ + if self.ncols() % (degree+1) != 0: + raise ValueError('The column number must be divisible by degree+1.') + coeff_matrices = [self[:,i*(self.ncols()//(degree+1)):(i+1)*(self.ncols()//(degree+1))] for i in range(degree+1)] + return sum([coeff_matrices[i]*variable**i for i in range(degree+1)]) + + def striped_krylov_matrix(self, J, degree, shift=None): + r""" + Return the block matrix of the following form, with rows permuted according to shift. + The following uses `E` to refer to ``self``, and `d` to refer to ``degree``. + + [ ] + [ E ] + [ ] + [-----] + [ ] + [ EJ ] + [ ] + [-----] + [ . ] + [ . ] + [ . ] + [-----] + [ ] + [ EJ^d] + [ ] + + INPUT: + + - ``J`` -- a square matrix of size equal to the number of columns of ``self``. + - ``degree`` -- integer, the maximum exponent for the Krylov matrix. + - ``shift`` -- list of integers (optional), row priority shifts. If ``None``, + defaults to all zero. + + OUTPUT: + + - A matrix with block rows [E, EJ, EJ^2, ..., EJ^d], row-permuted by shift. + + EXAMPLES: + + sage: R. = GF(97)[] + sage: E = matrix([[27,49,29],[50,58,0],[77,10,29]]) + sage: E + [27 49 29] + [50 58 0] + [77 10 29] + sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) + sage: J + [0 1 0] + [0 0 1] + [0 0 0] + sage: E.striped_krylov_matrix(J,3) + [27 49 29] + [50 58 0] + [77 10 29] + [ 0 27 49] + [ 0 50 58] + [ 0 77 10] + [ 0 0 27] + [ 0 0 50] + [ 0 0 77] + [ 0 0 0] + [ 0 0 0] + [ 0 0 0] + sage: E.striped_krylov_matrix(J,3,[0,3,6]) + [27 49 29] + [ 0 27 49] + [ 0 0 27] + [ 0 0 0] + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [ 0 0 0] + [77 10 29] + [ 0 77 10] + [ 0 0 77] + [ 0 0 0] + sage: E.striped_krylov_matrix(J,3,[3,0,2]) + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [77 10 29] + [27 49 29] + [ 0 0 0] + [ 0 77 10] + [ 0 27 49] + [ 0 0 77] + [ 0 0 27] + [ 0 0 0] + [ 0 0 0] + + """ + from sage.combinat.permutation import Permutation + from sage.matrix.constructor import matrix # for matrix.block + + m = self.nrows() + + # if no shift, this is equivalent to 0s + if shift is None: + shift = [0]*self.nrows() + + # size checks + if len(shift) != m: + raise ValueError(f"The shift should have the same number of elements as the rows of the matrix ({m}), but had {len(shift)}.") + if not J.is_square() or J.nrows() != self.ncols(): + raise ValueError(f"The matrix J should be a square matrix and match the number of columns of self ({self.ncols()}), but is of dimension {J.nrows()} x {J.ncols()}.") + + # define priority and indexing functions + # index: [0,1,...,m-1,m,...,2m-1,...,m*degree-1] -> [(0,0),(1,0),...,(m-1,0),(0,1),...,(m-1,1),...,(m-1,degree)] + priority = lambda c,d : shift[c] + d + index = lambda i : (i%m,i//m) + + # deduce permutation by sorting priorities + # rows should be ordered by ascending priority, then by associated row number in E (i%m) + priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) + priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) + + # build all blocks for the striped matrix + blocks = [] + + # for i from 0 to degree (inclusive), add E*J^i + J_i = matrix.identity(self.ncols()) + for i in range(degree+1): + blocks.append([self*J_i]) + if i < degree: + J_i *= J + + # return block matrix permuted according to shift + krylov = matrix.block(blocks,subdivide=False) + krylov.permute_rows(priority_permutation) + return krylov + + def naive_krylov_rank_profile(self, J, degree, shift=None): + r""" + Compute the rank profile (row and column) of the striped Krylov matrix + built from ``self`` and matrix `J`. + + INPUT: + + - ``J`` -- square matrix used in the Krylov construction. + - ``degree`` -- maximum power of J in Krylov matrix. + - ``shift`` -- list of integers (optional): priority row shift. + + OUTPUT: + + - A tuple (row_profile, pivot_matrix, column_profile): + * row_profile: list of the first r row indices in the striped Krylov matrix `K` corresponding to independent rows + * pivot: the submatrix of `K` given by those rows + * column_profile: list of the first r independent column indices in ``pivot_matrix`` + where r is the rank of `K`. + + EXAMPLES: + + sage: R. = GF(97)[] + sage: E = matrix([[27,49,29],[50,58,0],[77,10,29]]) + sage: E + [27 49 29] + [50 58 0] + [77 10 29] + sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) + sage: J + [0 1 0] + [0 0 1] + [0 0 0] + sage: E.striped_krylov_matrix(J,3) + [27 49 29] + [50 58 0] + [77 10 29] + [ 0 27 49] + [ 0 50 58] + [ 0 77 10] + [ 0 0 27] + [ 0 0 50] + [ 0 0 77] + [ 0 0 0] + [ 0 0 0] + [ 0 0 0] + sage: E.naive_krylov_rank_profile(J,3) + ( + [27 49 29] + [50 58 0] + (0, 1, 3), [ 0 27 49], (0, 1, 2) + ) + sage: E.striped_krylov_matrix(J,3,[0,3,6]) + [27 49 29] + [ 0 27 49] + [ 0 0 27] + [ 0 0 0] + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [ 0 0 0] + [77 10 29] + [ 0 77 10] + [ 0 0 77] + [ 0 0 0] + sage: E.naive_krylov_rank_profile(J,3,[0,3,6]) + ( + [27 49 29] + [ 0 27 49] + (0, 1, 2), [ 0 0 27], (0, 1, 2) + ) + sage: E.striped_krylov_matrix(J,3,[3,0,2]) + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [77 10 29] + [27 49 29] + [ 0 0 0] + [ 0 77 10] + [ 0 27 49] + [ 0 0 77] + [ 0 0 27] + [ 0 0 0] + [ 0 0 0] + sage: E.naive_krylov_rank_profile(J,3,[3,0,2]) + ( + [50 58 0] + [ 0 50 58] + (0, 1, 2), [ 0 0 50], (0, 1, 2) + ) + """ + from sage.matrix.constructor import matrix + + K = self.striped_krylov_matrix(J,degree,shift) + + # calculate first independent rows of K + row_profile = K.transpose().pivots() + + # construct submatrix + pivot = K.matrix_from_rows(row_profile) + + # calculate first independent columns of pivot + col_profile = pivot.pivots() + + return row_profile,pivot,col_profile + + def linear_interpolation_basis(self, J, degree, shift=None): + r""" + Construct a linear interpolant basis for (``self``,`J`) in `s`-Popov form. + + INPUT: + + - ``J`` -- square matrix for Krylov construction. + - ``degree`` -- upper bound on degree of minpoly(`J`), power of + - ``shift`` -- list of integers (optional): controls priority of rows. + + OUTPUT: + + - Matrix whose rows form an interpolation basis. + + """ + from sage.combinat.permutation import Permutation + from sage.matrix.constructor import matrix + + m = self.nrows() + + # if no shift, this is equivalent to 0s + if shift is None: + shift = [0]*self.nrows() + + # calculate shift priority function and permutation + priority = lambda c,d : shift[c] + d + index = lambda i : (i%m,i//m) + index_inv = lambda c,d : c + d*m + + priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) + priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) + + # maps row c of EJ^d to position i in striped Krylov matrix + # +/- 1 as permutations are 1-indexed, matrices are 0-indexed + phi = lambda c,d : priority_permutation.inverse()(index_inv(c,d) + 1) - 1 + phi_inv = lambda i : index(priority_permutation(i + 1) - 1) + + # calculate krylov profile + row_profile, pivot, col_profile = self.naive_krylov_rank_profile(J,degree,shift) + + # (c_k, d_k) = phi^-1 (row_i) + c, d = zip(*(phi_inv(i) for i in row_profile)) + + # degree_c = 0 or max d[k] such that c[k] = c + degree_c = [0]*m + for k in range(len(row_profile)): + if d[k] >= degree_c[c[k]]: + degree_c[c[k]] = d[k] + 1 + + # compute striped Krylov matrix + K = self.striped_krylov_matrix(J,degree,shift) + + # compute submatrix of target with rows phi(i,degree_c[i]) + target = K.matrix_from_rows([phi(i,degree_c[i]) for i in range(m)]) + + # compute submatrices of pivot and target with col_profile + C = pivot.matrix_from_columns(col_profile) + D = target.matrix_from_columns(col_profile) + + relation = D*C.inverse() + + # linear interpolation basis in shifted Popov form + uncompressed_basis = matrix.block([[-relation,matrix.identity(m)]],subdivide=False) + + return uncompressed_basis + # a limited number of access-only properties are provided for matrices @property def T(self): @@ -19562,6 +19903,8 @@ def _matrix_power_symbolic(A, n): return P * M * Pinv + + class NotFullRankError(ValueError): """ diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index db6ff8f4ccf..4e147c2f82e 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -5337,3 +5337,44 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): return False return True + + def expansion(self, degree=None): + r""" + Return the horizontal expansion of the matrix, + where the ith block corresponds to the coefficients for x^i. + + INPUT: + + - ``degree`` -- integer, the degree bound on the matrix; + if ``None``, defaults to degree of matrix. + + OUTPUT: matrix + + EXAMPLES: + sage: R. = GF(97)[] + sage: M = matrix([[x^2+36*x,31*x,0],[3*x+13,x+57,0],[96,96,1]]) + sage: M + [x^2 + 36*x 31*x 0] + [ 3*x + 13 x + 57 0] + [ 96 96 1] + sage: M.expansion() + [ 0 0 0 36 31 0 1 0 0] + [13 57 0 3 1 0 0 0 0] + [96 96 1 0 0 0 0 0 0] + sage: M.expansion(1) + [ 0 0 0 36 31 0] + [13 57 0 3 1 0] + [96 96 1 0 0 0] + sage: M.expansion(3) + [ 0 0 0 36 31 0 1 0 0 0 0 0] + [13 57 0 3 1 0 0 0 0 0 0 0] + [96 96 1 0 0 0 0 0 0 0 0 0] + """ + from sage.matrix.constructor import matrix # for matrix.block + + if degree is None: + degree = self.degree() + + # compute each coefficient matrix + coeff_matrices = [[self.coefficient_matrix(i) for i in range(degree+1)]] + return matrix.block(coeff_matrices,subdivide=False) \ No newline at end of file From f69c4026e21ec3c670c4f9abc4657ba78aff7d9e Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Mon, 16 Jun 2025 13:03:35 +0200 Subject: [PATCH 02/28] compression for linear interpolation basis --- src/sage/matrix/matrix2.pyx | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 15cdcbf7b90..e3289fb0e1c 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19214,7 +19214,7 @@ cdef class Matrix(Matrix1): return row_profile,pivot,col_profile - def linear_interpolation_basis(self, J, degree, shift=None): + def linear_interpolation_basis(self, J, degree, variable, shift=None): r""" Construct a linear interpolant basis for (``self``,`J`) in `s`-Popov form. @@ -19278,7 +19278,14 @@ cdef class Matrix(Matrix1): # linear interpolation basis in shifted Popov form uncompressed_basis = matrix.block([[-relation,matrix.identity(m)]],subdivide=False) - return uncompressed_basis + basis_rows = [[0]*m for i in range(m)] + for col in range(relation.ncols()): + for row in range(m): + basis_rows[row][c[col]] += uncompressed_basis[row][col] * variable**d[col] + for i in range(m): + basis_rows[i][i] += variable**degree_c[i] + + return matrix(basis_rows) # a limited number of access-only properties are provided for matrices @property From a415dccedfba2cab506e9a5f548e11cfa266a9a1 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Mon, 16 Jun 2025 17:11:53 +0200 Subject: [PATCH 03/28] fast krylov almost working, problem with indices being too large for lambda --- src/sage/matrix/matrix2.pyx | 139 ++++++++++++++++++++++++++++++++++++ 1 file changed, 139 insertions(+) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index e3289fb0e1c..a7a88491811 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19214,6 +19214,144 @@ cdef class Matrix(Matrix1): return row_profile,pivot,col_profile + def krylov_rank_profile(self, J, degree, shift=None): + r""" + Compute the rank profile (row and column) of the striped Krylov matrix + built from ``self`` and matrix `J`. + + INPUT: + + - ``J`` -- square matrix used in the Krylov construction. + - ``degree`` -- maximum power of J in Krylov matrix. + - ``shift`` -- list of integers (optional): priority row shift. + + OUTPUT: + + - A tuple (row_profile, pivot_matrix, column_profile): + * row_profile: list of the first r row indices in the striped Krylov matrix `K` corresponding to independent rows + * pivot: the submatrix of `K` given by those rows + * column_profile: list of the first r independent column indices in ``pivot_matrix`` + where r is the rank of `K`. + + EXAMPLES: + + sage: R. = GF(97)[] + sage: E = matrix([[27,49,29],[50,58,0],[77,10,29]]) + sage: E + [27 49 29] + [50 58 0] + [77 10 29] + sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) + sage: J + [0 1 0] + [0 0 1] + [0 0 0] + sage: E.striped_krylov_matrix(J,3) + [27 49 29] + [50 58 0] + [77 10 29] + [ 0 27 49] + [ 0 50 58] + [ 0 77 10] + [ 0 0 27] + [ 0 0 50] + [ 0 0 77] + [ 0 0 0] + [ 0 0 0] + [ 0 0 0] + sage: E.naive_krylov_rank_profile(J,3) + ( + [27 49 29] + [50 58 0] + (0, 1, 3), [ 0 27 49], (0, 1, 2) + ) + sage: E.striped_krylov_matrix(J,3,[0,3,6]) + [27 49 29] + [ 0 27 49] + [ 0 0 27] + [ 0 0 0] + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [ 0 0 0] + [77 10 29] + [ 0 77 10] + [ 0 0 77] + [ 0 0 0] + sage: E.naive_krylov_rank_profile(J,3,[0,3,6]) + ( + [27 49 29] + [ 0 27 49] + (0, 1, 2), [ 0 0 27], (0, 1, 2) + ) + sage: E.striped_krylov_matrix(J,3,[3,0,2]) + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [77 10 29] + [27 49 29] + [ 0 0 0] + [ 0 77 10] + [ 0 27 49] + [ 0 0 77] + [ 0 0 27] + [ 0 0 0] + [ 0 0 0] + sage: E.naive_krylov_rank_profile(J,3,[3,0,2]) + ( + [50 58 0] + [ 0 50 58] + (0, 1, 2), [ 0 0 50], (0, 1, 2) + ) + """ + from sage.matrix.constructor import matrix + import math + from sage.combinat.permutation import Permutation + + m = self.nrows() + + # calculate shift priority function and permutation + priority = lambda c,d : shift[c] + d + index = lambda i : (i%m,i//m) + index_inv = lambda c,d : c + d*m + + priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) + priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) + + # maps row c of EJ^d to position i in striped Krylov matrix + # +/- 1 as permutations are 1-indexed, matrices are 0-indexed + phi = lambda c,d : priority_permutation.inverse()(index_inv(c,d) + 1) - 1 + phi_inv = lambda i : index(priority_permutation(i + 1) - 1) + + row_profile_self = self.transpose().pivots() + r = len(row_profile_self) + row_profile_K = sorted([phi(row,0) for row in row_profile_self]) + + M = matrix([self.row(index_inv(*phi_inv(i))) for i in row_profile_K]) + + J_L = None + for l in range(math.ceil(math.log(degree,2)) + 1): + c,d = zip(*(phi_inv(i) for i in row_profile_K)) + L = pow(2,l) + row_extension = sorted([phi(c[i],d[i] + L) for i in range(r)]) + + k = sorted(row_profile_K + row_extension) + + if J_L is None: + J_L = J + else: + J_L = J_L * J_L + + M = matrix.block([[M],[M*J_L]],subdivide=False) + + row_profile_M = M.transpose().pivots() + M = matrix([M.row(i) for i in row_profile_M]) + + row_profile_K = [k[i] for i in row_profile_M] + r = len(row_profile_K) + + print(c,d,L,row_extension,k) + def linear_interpolation_basis(self, J, degree, variable, shift=None): r""" Construct a linear interpolant basis for (``self``,`J`) in `s`-Popov form. @@ -19278,6 +19416,7 @@ cdef class Matrix(Matrix1): # linear interpolation basis in shifted Popov form uncompressed_basis = matrix.block([[-relation,matrix.identity(m)]],subdivide=False) + # compression of basis into polynomial form basis_rows = [[0]*m for i in range(m)] for col in range(relation.ncols()): for row in range(m): From 59663bcb43a24ccafe24e1952ef9a222bd11fc1c Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Wed, 18 Jun 2025 12:54:51 +0200 Subject: [PATCH 04/28] krylov_rank_profile working now --- src/sage/matrix/matrix2.pyx | 41 ++++++++++++++++++++++++------------- 1 file changed, 27 insertions(+), 14 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index a7a88491811..1291dbe3cfb 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19259,7 +19259,7 @@ cdef class Matrix(Matrix1): [ 0 0 0] [ 0 0 0] [ 0 0 0] - sage: E.naive_krylov_rank_profile(J,3) + sage: E.krylov_rank_profile(J,3) ( [27 49 29] [50 58 0] @@ -19278,7 +19278,7 @@ cdef class Matrix(Matrix1): [ 0 77 10] [ 0 0 77] [ 0 0 0] - sage: E.naive_krylov_rank_profile(J,3,[0,3,6]) + sage: E.krylov_rank_profile(J,3,[0,3,6]) ( [27 49 29] [ 0 27 49] @@ -19297,7 +19297,7 @@ cdef class Matrix(Matrix1): [ 0 0 27] [ 0 0 0] [ 0 0 0] - sage: E.naive_krylov_rank_profile(J,3,[3,0,2]) + sage: E.krylov_rank_profile(J,3,[3,0,2]) ( [50 58 0] [ 0 50 58] @@ -19318,39 +19318,52 @@ cdef class Matrix(Matrix1): priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) - # maps row c of EJ^d to position i in striped Krylov matrix + # maps row c of self*J^d to position i in striped Krylov matrix # +/- 1 as permutations are 1-indexed, matrices are 0-indexed phi = lambda c,d : priority_permutation.inverse()(index_inv(c,d) + 1) - 1 phi_inv = lambda i : index(priority_permutation(i + 1) - 1) - row_profile_self = self.transpose().pivots() - r = len(row_profile_self) - row_profile_K = sorted([phi(row,0) for row in row_profile_self]) + # calculate row profile of self, with shift applied + self_permutation = Permutation([pair[1]+1 for pair in sorted([[phi(i,0),i] for i in range(m)])]) + row_profile_self_permuted = self.with_permuted_rows(self_permutation).transpose().pivots() + row_profile_self = [self_permutation(i+1)-1 for i in row_profile_self_permuted] + + # base row_profile_K using row_profile_self + row_profile_K = sorted([phi(i,0) for i in row_profile_self]) + r = len(row_profile_K) M = matrix([self.row(index_inv(*phi_inv(i))) for i in row_profile_K]) J_L = None for l in range(math.ceil(math.log(degree,2)) + 1): + # row and degree of profile within unshifted expansion c,d = zip(*(phi_inv(i) for i in row_profile_K)) L = pow(2,l) - row_extension = sorted([phi(c[i],d[i] + L) for i in range(r)]) + # adding 2^l to each degree + row_extension = sorted([phi(c[i],d[i] + L) for i in range(r) if d[i] + L <= degree]) - k = sorted(row_profile_K + row_extension) + # concatenate two sequences (order preserved) + k = row_profile_K + row_extension + # calculate sorting permutation + k_perm = Permutation([x[1] for x in sorted([k[i],i+1] for i in range(len(k)))]) + # fast calculation of rows formed by indices in k if J_L is None: J_L = J else: J_L = J_L * J_L - M = matrix.block([[M],[M*J_L]],subdivide=False) + # sort rows of M, find profile, translate to k (indices of full krylov matrix) + M.permute_rows(k_perm) row_profile_M = M.transpose().pivots() - M = matrix([M.row(i) for i in row_profile_M]) - - row_profile_K = [k[i] for i in row_profile_M] r = len(row_profile_K) + row_profile_K = [k[k_perm(i+1)-1] for i in row_profile_M] - print(c,d,L,row_extension,k) + # calculate new M for return value or next loop + M = matrix([M.row(i) for i in row_profile_M]) + col_profile = M.pivots() + return tuple(row_profile_K), M, col_profile def linear_interpolation_basis(self, J, degree, variable, shift=None): r""" From 21882b5973a621c84d8a03fcbff0cf29c6de3786 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Wed, 18 Jun 2025 14:28:52 +0200 Subject: [PATCH 05/28] optimised target calculation --- src/sage/matrix/matrix2.pyx | 25 +++++++++++++++++++++++-- 1 file changed, 23 insertions(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 1291dbe3cfb..aefc35e7da7 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19409,16 +19409,37 @@ cdef class Matrix(Matrix1): c, d = zip(*(phi_inv(i) for i in row_profile)) # degree_c = 0 or max d[k] such that c[k] = c + inv_c = [None]*m degree_c = [0]*m for k in range(len(row_profile)): if d[k] >= degree_c[c[k]]: degree_c[c[k]] = d[k] + 1 + inv_c[c[k]] = k + ################# + # OPTIMISATION + ################# + T_present = matrix([pivot.row(inv_c[i]) for i in range(m) if inv_c[i] is not None]) * J + T_absent = matrix([self.row(i) for i in range(m) if inv_c[i] is None]) + T_rows = [] + idx_p = 0 + idx_a = 0 + for i in range(m): + if inv_c[i] is None: + T_rows.append(T_absent[idx_a]) + idx_a += 1 + else: + T_rows.append(T_present[idx_p]) + idx_p += 1 + target = matrix(T_rows) + ################# + # ELSE # compute striped Krylov matrix - K = self.striped_krylov_matrix(J,degree,shift) + #K = self.striped_krylov_matrix(J,degree,shift) # compute submatrix of target with rows phi(i,degree_c[i]) - target = K.matrix_from_rows([phi(i,degree_c[i]) for i in range(m)]) + #target = K.matrix_from_rows([phi(i,degree_c[i]) for i in range(m)]) + ################# # compute submatrices of pivot and target with col_profile C = pivot.matrix_from_columns(col_profile) From ce717fedbee83298727650d0aa35e0b3be51812c Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Thu, 19 Jun 2025 03:39:43 +0200 Subject: [PATCH 06/28] combine row and column submatrix operations in basis calc --- src/sage/matrix/matrix2.pyx | 39 +++++++++++++------------------------ 1 file changed, 14 insertions(+), 25 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index aefc35e7da7..e5e4902bdfc 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19365,7 +19365,7 @@ cdef class Matrix(Matrix1): col_profile = M.pivots() return tuple(row_profile_K), M, col_profile - def linear_interpolation_basis(self, J, degree, variable, shift=None): + def linear_interpolation_basis(self, J, degree, variable, shift=None, opt=False): r""" Construct a linear interpolant basis for (``self``,`J`) in `s`-Popov form. @@ -19416,34 +19416,23 @@ cdef class Matrix(Matrix1): degree_c[c[k]] = d[k] + 1 inv_c[c[k]] = k - ################# - # OPTIMISATION - ################# - T_present = matrix([pivot.row(inv_c[i]) for i in range(m) if inv_c[i] is not None]) * J - T_absent = matrix([self.row(i) for i in range(m) if inv_c[i] is None]) - T_rows = [] + T_absent_indices = [i for i in range(m) if inv_c[i] is None] + T_present_indices = [inv_c[i] for i in range(m) if inv_c[i] is not None] + D_absent = self.matrix_from_rows_and_columns(T_absent_indices, col_profile) + D_present = (pivot.matrix_from_rows(T_present_indices)*J).matrix_from_columns(col_profile) + D_rows = [] idx_p = 0 idx_a = 0 for i in range(m): - if inv_c[i] is None: - T_rows.append(T_absent[idx_a]) - idx_a += 1 - else: - T_rows.append(T_present[idx_p]) - idx_p += 1 - target = matrix(T_rows) - ################# - # ELSE - # compute striped Krylov matrix - #K = self.striped_krylov_matrix(J,degree,shift) - - # compute submatrix of target with rows phi(i,degree_c[i]) - #target = K.matrix_from_rows([phi(i,degree_c[i]) for i in range(m)]) - ################# - - # compute submatrices of pivot and target with col_profile + if inv_c[i] is None: + D_rows.append(D_absent[idx_a]) + idx_a += 1 + else: + D_rows.append(D_present[idx_p]) + idx_p += 1 + C = pivot.matrix_from_columns(col_profile) - D = target.matrix_from_columns(col_profile) + D = matrix(D_rows) relation = D*C.inverse() From c4ae3ea14abcd9347c65f67c26cc17844f3655cf Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Thu, 19 Jun 2025 10:34:54 +0200 Subject: [PATCH 07/28] small bug fixes and cases with 0 rows --- src/sage/matrix/matrix2.pyx | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index e5e4902bdfc..4edc907ed38 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19310,6 +19310,9 @@ cdef class Matrix(Matrix1): m = self.nrows() + if m == 0: + return (),self,() + # calculate shift priority function and permutation priority = lambda c,d : shift[c] + d index = lambda i : (i%m,i//m) @@ -19332,7 +19335,10 @@ cdef class Matrix(Matrix1): row_profile_K = sorted([phi(i,0) for i in row_profile_self]) r = len(row_profile_K) - M = matrix([self.row(index_inv(*phi_inv(i))) for i in row_profile_K]) + if r == 0: + return (),self.matrix_from_rows([]),() + + M = self.matrix_from_rows([index_inv(*phi_inv(i)) for i in row_profile_K]) J_L = None for l in range(math.ceil(math.log(degree,2)) + 1): @@ -19357,7 +19363,7 @@ cdef class Matrix(Matrix1): # sort rows of M, find profile, translate to k (indices of full krylov matrix) M.permute_rows(k_perm) row_profile_M = M.transpose().pivots() - r = len(row_profile_K) + r = len(row_profile_M) row_profile_K = [k[k_perm(i+1)-1] for i in row_profile_M] # calculate new M for return value or next loop @@ -19365,7 +19371,7 @@ cdef class Matrix(Matrix1): col_profile = M.pivots() return tuple(row_profile_K), M, col_profile - def linear_interpolation_basis(self, J, degree, variable, shift=None, opt=False): + def linear_interpolation_basis(self, J, degree, variable, shift=None): r""" Construct a linear interpolant basis for (``self``,`J`) in `s`-Popov form. @@ -19403,7 +19409,10 @@ cdef class Matrix(Matrix1): phi_inv = lambda i : index(priority_permutation(i + 1) - 1) # calculate krylov profile - row_profile, pivot, col_profile = self.naive_krylov_rank_profile(J,degree,shift) + row_profile, pivot, col_profile = self.krylov_rank_profile(J,degree,shift) + + if len(row_profile) == 0: + return matrix.identity(self.base_ring(),m) # (c_k, d_k) = phi^-1 (row_i) c, d = zip(*(phi_inv(i) for i in row_profile)) From 308333f230e4577a2e7c88389251c8aaee32b75e Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Thu, 19 Jun 2025 14:18:50 +0200 Subject: [PATCH 08/28] performance fix for lambda - calculate inverse once --- src/sage/matrix/matrix2.pyx | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 4edc907ed38..de8bbd3f3e1 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19306,8 +19306,14 @@ cdef class Matrix(Matrix1): """ from sage.matrix.constructor import matrix import math + import time from sage.combinat.permutation import Permutation + #DEBUG + time_idx = 0 + start = time.time() + #DEBUG + m = self.nrows() if m == 0: @@ -19320,14 +19326,16 @@ cdef class Matrix(Matrix1): priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) + priority_permutation_inv = priority_permutation.inverse() # maps row c of self*J^d to position i in striped Krylov matrix # +/- 1 as permutations are 1-indexed, matrices are 0-indexed - phi = lambda c,d : priority_permutation.inverse()(index_inv(c,d) + 1) - 1 + phi = lambda c,d : priority_permutation_inv(index_inv(c,d) + 1) - 1 phi_inv = lambda i : index(priority_permutation(i + 1) - 1) # calculate row profile of self, with shift applied self_permutation = Permutation([pair[1]+1 for pair in sorted([[phi(i,0),i] for i in range(m)])]) + row_profile_self_permuted = self.with_permuted_rows(self_permutation).transpose().pivots() row_profile_self = [self_permutation(i+1)-1 for i in row_profile_self_permuted] @@ -19402,10 +19410,11 @@ cdef class Matrix(Matrix1): priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) + priority_permutation_inv = priority_permutation.inverse() # maps row c of EJ^d to position i in striped Krylov matrix # +/- 1 as permutations are 1-indexed, matrices are 0-indexed - phi = lambda c,d : priority_permutation.inverse()(index_inv(c,d) + 1) - 1 + phi = lambda c,d : priority_permutation_inv(index_inv(c,d) + 1) - 1 phi_inv = lambda i : index(priority_permutation(i + 1) - 1) # calculate krylov profile From 4126fec1cfaec559d100bc933ccce2c2764d69bf Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Fri, 20 Jun 2025 13:19:23 +0200 Subject: [PATCH 09/28] remove debug code, tidy up variable handling for polynomial ring --- src/sage/matrix/matrix2.pyx | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index de8bbd3f3e1..2852e8084ab 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -18939,7 +18939,7 @@ cdef class Matrix(Matrix1): U.rescale_col(n - 1, -1) return U - def polynomial_compression(self,degree,variable): + def polynomial_compression(self,degree,var_name): """ Returns the corresponding polynomial matrix where the coefficient matrix of degree i is the ith block of ``self``. @@ -18978,6 +18978,11 @@ cdef class Matrix(Matrix1): [ 3*y + 13 y + 57 0] [ 96 96 1] """ + from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing + + poly_ring = PolynomialRing(self.base_ring(),var_name) + variable = poly_ring.gen() + if self.ncols() % (degree+1) != 0: raise ValueError('The column number must be divisible by degree+1.') coeff_matrices = [self[:,i*(self.ncols()//(degree+1)):(i+1)*(self.ncols()//(degree+1))] for i in range(degree+1)] @@ -19306,14 +19311,8 @@ cdef class Matrix(Matrix1): """ from sage.matrix.constructor import matrix import math - import time from sage.combinat.permutation import Permutation - #DEBUG - time_idx = 0 - start = time.time() - #DEBUG - m = self.nrows() if m == 0: @@ -19379,7 +19378,7 @@ cdef class Matrix(Matrix1): col_profile = M.pivots() return tuple(row_profile_K), M, col_profile - def linear_interpolation_basis(self, J, degree, variable, shift=None): + def linear_interpolation_basis(self, J, degree, var_name, shift=None): r""" Construct a linear interpolant basis for (``self``,`J`) in `s`-Popov form. @@ -19396,6 +19395,9 @@ cdef class Matrix(Matrix1): """ from sage.combinat.permutation import Permutation from sage.matrix.constructor import matrix + from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing + + poly_ring = PolynomialRing(self.base_ring(),var_name) m = self.nrows() @@ -19421,7 +19423,7 @@ cdef class Matrix(Matrix1): row_profile, pivot, col_profile = self.krylov_rank_profile(J,degree,shift) if len(row_profile) == 0: - return matrix.identity(self.base_ring(),m) + return matrix.identity(poly_ring,m) # (c_k, d_k) = phi^-1 (row_i) c, d = zip(*(phi_inv(i) for i in row_profile)) @@ -19457,6 +19459,9 @@ cdef class Matrix(Matrix1): # linear interpolation basis in shifted Popov form uncompressed_basis = matrix.block([[-relation,matrix.identity(m)]],subdivide=False) + # construct variable + variable = poly_ring.gen() + # compression of basis into polynomial form basis_rows = [[0]*m for i in range(m)] for col in range(relation.ncols()): From 346aaab776aedf5bbd83fd9452ff774239ad4172 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Tue, 1 Jul 2025 21:06:47 +0200 Subject: [PATCH 10/28] add caching for row_pivots in calculation of col_pivots for matrix_modn_dense, add aliases for pivots and pivot_rows, optimisation for krylov_rank_profile by minimising index conversion and full permutation calculation --- src/sage/matrix/matrix2.pyx | 178 +++++++++++++++++- .../matrix/matrix_modn_dense_template.pxi | 8 +- 2 files changed, 177 insertions(+), 9 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 2852e8084ab..20ef0354464 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -1027,6 +1027,22 @@ cdef class Matrix(Matrix1): v = self.transpose().pivots() self.cache('pivot_rows', v) return v + + def row_rank_profile(self): + """ + Alias for pivot_rows(self). Return the pivot row positions for this matrix, which are a topmost + subset of the rows that span the row space and are linearly + independent. + """ + return self.pivot_rows() + + def col_rank_profile(self): + """ + Alias for pivots(self). Return the pivot column positions for this matrix, which are a leftmost + subset of the columns that span the column space and are linearly + independent. + """ + return self.pivots() def _solve_right_general(self, B, check=True): r""" @@ -19209,13 +19225,13 @@ cdef class Matrix(Matrix1): K = self.striped_krylov_matrix(J,degree,shift) # calculate first independent rows of K - row_profile = K.transpose().pivots() + row_profile = K.row_rank_profile() # construct submatrix pivot = K.matrix_from_rows(row_profile) # calculate first independent columns of pivot - col_profile = pivot.pivots() + col_profile = pivot.col_rank_profile() return row_profile,pivot,col_profile @@ -19335,7 +19351,7 @@ cdef class Matrix(Matrix1): # calculate row profile of self, with shift applied self_permutation = Permutation([pair[1]+1 for pair in sorted([[phi(i,0),i] for i in range(m)])]) - row_profile_self_permuted = self.with_permuted_rows(self_permutation).transpose().pivots() + row_profile_self_permuted = self.with_permuted_rows(self_permutation).row_rank_profile() row_profile_self = [self_permutation(i+1)-1 for i in row_profile_self_permuted] # base row_profile_K using row_profile_self @@ -19352,11 +19368,13 @@ cdef class Matrix(Matrix1): # row and degree of profile within unshifted expansion c,d = zip(*(phi_inv(i) for i in row_profile_K)) L = pow(2,l) + # adding 2^l to each degree row_extension = sorted([phi(c[i],d[i] + L) for i in range(r) if d[i] + L <= degree]) # concatenate two sequences (order preserved) k = row_profile_K + row_extension + # calculate sorting permutation k_perm = Permutation([x[1] for x in sorted([k[i],i+1] for i in range(len(k)))]) @@ -19365,17 +19383,165 @@ cdef class Matrix(Matrix1): J_L = J else: J_L = J_L * J_L - M = matrix.block([[M],[M*J_L]],subdivide=False) + M = M.stack(M*J_L,subdivide=False) # sort rows of M, find profile, translate to k (indices of full krylov matrix) M.permute_rows(k_perm) - row_profile_M = M.transpose().pivots() + row_profile_M = M.row_rank_profile() r = len(row_profile_M) row_profile_K = [k[k_perm(i+1)-1] for i in row_profile_M] # calculate new M for return value or next loop M = matrix([M.row(i) for i in row_profile_M]) - col_profile = M.pivots() + col_profile = M.col_rank_profile() + return tuple(row_profile_K), M, col_profile + + def krylov_rank_profile_fast_perm(self, J, degree, shift=None): + r""" + Compute the rank profile (row and column) of the striped Krylov matrix + built from ``self`` and matrix `J`. + + INPUT: + + - ``J`` -- square matrix used in the Krylov construction. + - ``degree`` -- maximum power of J in Krylov matrix. + - ``shift`` -- list of integers (optional): priority row shift. + + OUTPUT: + + - A tuple (row_profile, pivot_matrix, column_profile): + * row_profile: list of the first r row indices in the striped Krylov matrix `K` corresponding to independent rows + * pivot: the submatrix of `K` given by those rows + * column_profile: list of the first r independent column indices in ``pivot_matrix`` + where r is the rank of `K`. + + EXAMPLES: + + sage: R. = GF(97)[] + sage: E = matrix([[27,49,29],[50,58,0],[77,10,29]]) + sage: E + [27 49 29] + [50 58 0] + [77 10 29] + sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) + sage: J + [0 1 0] + [0 0 1] + [0 0 0] + sage: E.striped_krylov_matrix(J,3) + [27 49 29] + [50 58 0] + [77 10 29] + [ 0 27 49] + [ 0 50 58] + [ 0 77 10] + [ 0 0 27] + [ 0 0 50] + [ 0 0 77] + [ 0 0 0] + [ 0 0 0] + [ 0 0 0] + sage: E.krylov_rank_profile(J,3) + ( + [27 49 29] + [50 58 0] + (0, 1, 3), [ 0 27 49], (0, 1, 2) + ) + sage: E.striped_krylov_matrix(J,3,[0,3,6]) + [27 49 29] + [ 0 27 49] + [ 0 0 27] + [ 0 0 0] + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [ 0 0 0] + [77 10 29] + [ 0 77 10] + [ 0 0 77] + [ 0 0 0] + sage: E.krylov_rank_profile(J,3,[0,3,6]) + ( + [27 49 29] + [ 0 27 49] + (0, 1, 2), [ 0 0 27], (0, 1, 2) + ) + sage: E.striped_krylov_matrix(J,3,[3,0,2]) + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [77 10 29] + [27 49 29] + [ 0 0 0] + [ 0 77 10] + [ 0 27 49] + [ 0 0 77] + [ 0 0 27] + [ 0 0 0] + [ 0 0 0] + sage: E.krylov_rank_profile(J,3,[3,0,2]) + ( + [50 58 0] + [ 0 50 58] + (0, 1, 2), [ 0 0 50], (0, 1, 2) + ) + """ + from sage.matrix.constructor import matrix + import math + from sage.combinat.permutation import Permutation + + m = self.nrows() + + if m == 0: + return (),self,() + + # calculate row profile of self, with shift applied + self_permutation = Permutation(sorted([i+1 for i in range(m)],key=lambda x:(shift[x-1],x-1))) + + self_permuted = self.with_permuted_rows(self_permutation) + row_profile_self_permuted = self_permuted.row_rank_profile() + row_profile_self = [(self_permutation(i+1)-1,0) for i in row_profile_self_permuted] + + r = len(row_profile_self) + + if r == 0: + return (),self.matrix_from_rows([]),() + + M = self_permuted.matrix_from_rows(row_profile_self_permuted) + + J_L = None + + for l in range(math.ceil(math.log(degree,2)) + 1): + L = pow(2,l) + # adding 2^l to each degree + row_extension = [(x[0],x[1] + L) for x in row_profile_self if x[1] + L <= degree] + + # concatenate two sequences (order preserved) + k = row_profile_self + row_extension + + # calculate sorting permutation, sort k by (shift[c]+d, c) + k_perm = Permutation(sorted([i+1 for i in range(len(k))],key=lambda x: (shift[k[x-1][0]] + k[x-1][1],k[x-1][0]))) + + # fast calculation of rows formed by indices in k + if J_L is None: + J_L = J + else: + J_L = J_L * J_L + M = M.stack(M*J_L,subdivide=False) + + # sort rows of M, find profile, translate to k (indices of full krylov matrix) + M.permute_rows(k_perm) + row_profile_M = M.row_rank_profile() + r = len(row_profile_M) + row_profile_self = sorted([k[k_perm(i+1)-1] for i in row_profile_M],key=lambda x: (shift[x[0]]+x[1],x[0])) + + # calculate new M for return value or next loop + M = matrix([M.row(i) for i in row_profile_M]) + # convert c,d to actual position in striped Krylov matrix + phi = lambda c,d : sum(min(max(shift[c] - shift[i] + d + (i < c and shift[i] <= shift[c] + degree),0),degree+1) for i in range(len(shift))) + row_profile_K = [phi(*row) for row in row_profile_self] + + col_profile = M.col_rank_profile() return tuple(row_profile_K), M, col_profile def linear_interpolation_basis(self, J, degree, var_name, shift=None): diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi index 8e9e7bf0f89..63fac0af4b0 100644 --- a/src/sage/matrix/matrix_modn_dense_template.pxi +++ b/src/sage/matrix/matrix_modn_dense_template.pxi @@ -204,12 +204,13 @@ cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nr applyP(F[0], FflasRight, FflasNoTrans, nrows, 0, r, entries, ncols, Q) - cdef list pivots = [int(Q[i]) for i in range(r)] + cdef list column_pivots = [int(Q[i]) for i in range(r)] + cdef list row_pivots = sorted([int(P[i]) for i in range(r)]) sig_free(P) sig_free(Q) del F - return r, pivots + return r, column_pivots, row_pivots cdef inline linbox_echelonize_efd(celement modulus, celement* entries, Py_ssize_t nrows, Py_ssize_t ncols): # See trac #13878: This is to avoid sending invalid data to linbox, @@ -1788,8 +1789,9 @@ cdef class Matrix_modn_dense_template(Matrix_dense): r, pivots = linbox_echelonize_efd(self.p, self._entries, self._nrows, self._ncols) else: - r, pivots = linbox_echelonize(self.p, self._entries, + r, pivots, row_pivots = linbox_echelonize(self.p, self._entries, self._nrows, self._ncols) + self.cache('pivot_rows',tuple(row_pivots)) verbose('done with echelonize', t) self.cache('in_echelon_form', True) self.cache('rank', r) From a41e6ec5d08b60428c3c332387885440fba64bfe Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Wed, 2 Jul 2025 14:01:47 +0200 Subject: [PATCH 11/28] make optimisation fast perm profile main implementation --- src/sage/matrix/matrix2.pyx | 161 ------------------------------------ 1 file changed, 161 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 20ef0354464..779cf473084 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19331,167 +19331,6 @@ cdef class Matrix(Matrix1): m = self.nrows() - if m == 0: - return (),self,() - - # calculate shift priority function and permutation - priority = lambda c,d : shift[c] + d - index = lambda i : (i%m,i//m) - index_inv = lambda c,d : c + d*m - - priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) - priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) - priority_permutation_inv = priority_permutation.inverse() - - # maps row c of self*J^d to position i in striped Krylov matrix - # +/- 1 as permutations are 1-indexed, matrices are 0-indexed - phi = lambda c,d : priority_permutation_inv(index_inv(c,d) + 1) - 1 - phi_inv = lambda i : index(priority_permutation(i + 1) - 1) - - # calculate row profile of self, with shift applied - self_permutation = Permutation([pair[1]+1 for pair in sorted([[phi(i,0),i] for i in range(m)])]) - - row_profile_self_permuted = self.with_permuted_rows(self_permutation).row_rank_profile() - row_profile_self = [self_permutation(i+1)-1 for i in row_profile_self_permuted] - - # base row_profile_K using row_profile_self - row_profile_K = sorted([phi(i,0) for i in row_profile_self]) - r = len(row_profile_K) - - if r == 0: - return (),self.matrix_from_rows([]),() - - M = self.matrix_from_rows([index_inv(*phi_inv(i)) for i in row_profile_K]) - - J_L = None - for l in range(math.ceil(math.log(degree,2)) + 1): - # row and degree of profile within unshifted expansion - c,d = zip(*(phi_inv(i) for i in row_profile_K)) - L = pow(2,l) - - # adding 2^l to each degree - row_extension = sorted([phi(c[i],d[i] + L) for i in range(r) if d[i] + L <= degree]) - - # concatenate two sequences (order preserved) - k = row_profile_K + row_extension - - # calculate sorting permutation - k_perm = Permutation([x[1] for x in sorted([k[i],i+1] for i in range(len(k)))]) - - # fast calculation of rows formed by indices in k - if J_L is None: - J_L = J - else: - J_L = J_L * J_L - M = M.stack(M*J_L,subdivide=False) - - # sort rows of M, find profile, translate to k (indices of full krylov matrix) - M.permute_rows(k_perm) - row_profile_M = M.row_rank_profile() - r = len(row_profile_M) - row_profile_K = [k[k_perm(i+1)-1] for i in row_profile_M] - - # calculate new M for return value or next loop - M = matrix([M.row(i) for i in row_profile_M]) - col_profile = M.col_rank_profile() - return tuple(row_profile_K), M, col_profile - - def krylov_rank_profile_fast_perm(self, J, degree, shift=None): - r""" - Compute the rank profile (row and column) of the striped Krylov matrix - built from ``self`` and matrix `J`. - - INPUT: - - - ``J`` -- square matrix used in the Krylov construction. - - ``degree`` -- maximum power of J in Krylov matrix. - - ``shift`` -- list of integers (optional): priority row shift. - - OUTPUT: - - - A tuple (row_profile, pivot_matrix, column_profile): - * row_profile: list of the first r row indices in the striped Krylov matrix `K` corresponding to independent rows - * pivot: the submatrix of `K` given by those rows - * column_profile: list of the first r independent column indices in ``pivot_matrix`` - where r is the rank of `K`. - - EXAMPLES: - - sage: R. = GF(97)[] - sage: E = matrix([[27,49,29],[50,58,0],[77,10,29]]) - sage: E - [27 49 29] - [50 58 0] - [77 10 29] - sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) - sage: J - [0 1 0] - [0 0 1] - [0 0 0] - sage: E.striped_krylov_matrix(J,3) - [27 49 29] - [50 58 0] - [77 10 29] - [ 0 27 49] - [ 0 50 58] - [ 0 77 10] - [ 0 0 27] - [ 0 0 50] - [ 0 0 77] - [ 0 0 0] - [ 0 0 0] - [ 0 0 0] - sage: E.krylov_rank_profile(J,3) - ( - [27 49 29] - [50 58 0] - (0, 1, 3), [ 0 27 49], (0, 1, 2) - ) - sage: E.striped_krylov_matrix(J,3,[0,3,6]) - [27 49 29] - [ 0 27 49] - [ 0 0 27] - [ 0 0 0] - [50 58 0] - [ 0 50 58] - [ 0 0 50] - [ 0 0 0] - [77 10 29] - [ 0 77 10] - [ 0 0 77] - [ 0 0 0] - sage: E.krylov_rank_profile(J,3,[0,3,6]) - ( - [27 49 29] - [ 0 27 49] - (0, 1, 2), [ 0 0 27], (0, 1, 2) - ) - sage: E.striped_krylov_matrix(J,3,[3,0,2]) - [50 58 0] - [ 0 50 58] - [ 0 0 50] - [77 10 29] - [27 49 29] - [ 0 0 0] - [ 0 77 10] - [ 0 27 49] - [ 0 0 77] - [ 0 0 27] - [ 0 0 0] - [ 0 0 0] - sage: E.krylov_rank_profile(J,3,[3,0,2]) - ( - [50 58 0] - [ 0 50 58] - (0, 1, 2), [ 0 0 50], (0, 1, 2) - ) - """ - from sage.matrix.constructor import matrix - import math - from sage.combinat.permutation import Permutation - - m = self.nrows() - if m == 0: return (),self,() From 3cc05bd13300c42819d92f224ee49b6b5ac42896 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Thu, 3 Jul 2025 12:27:59 +0200 Subject: [PATCH 12/28] optimised version of basis calculation without permutation --- src/sage/matrix/matrix2.pyx | 104 +++++++++++++++++++++++++++++++++++- 1 file changed, 102 insertions(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 779cf473084..df4f389eedf 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19235,7 +19235,7 @@ cdef class Matrix(Matrix1): return row_profile,pivot,col_profile - def krylov_rank_profile(self, J, degree, shift=None): + def krylov_rank_profile(self, J, degree, shift=None, output_pairs=False): r""" Compute the rank profile (row and column) of the striped Krylov matrix built from ``self`` and matrix `J`. @@ -19376,11 +19376,16 @@ cdef class Matrix(Matrix1): # calculate new M for return value or next loop M = matrix([M.row(i) for i in row_profile_M]) + + col_profile = M.col_rank_profile() + + if output_pairs: + return tuple(row_profile_self), M, col_profile + # convert c,d to actual position in striped Krylov matrix phi = lambda c,d : sum(min(max(shift[c] - shift[i] + d + (i < c and shift[i] <= shift[c] + degree),0),degree+1) for i in range(len(shift))) row_profile_K = [phi(*row) for row in row_profile_self] - col_profile = M.col_rank_profile() return tuple(row_profile_K), M, col_profile def linear_interpolation_basis(self, J, degree, var_name, shift=None): @@ -19476,6 +19481,101 @@ cdef class Matrix(Matrix1): basis_rows[i][i] += variable**degree_c[i] return matrix(basis_rows) + + def linear_interpolation_basis_fast_perm(self, J, degree, var_name, shift=None): + r""" + Construct a linear interpolant basis for (``self``,`J`) in `s`-Popov form. + + INPUT: + + - ``J`` -- square matrix for Krylov construction. + - ``degree`` -- upper bound on degree of minpoly(`J`), power of + - ``shift`` -- list of integers (optional): controls priority of rows. + + OUTPUT: + + - Matrix whose rows form an interpolation basis. + + """ + from sage.combinat.permutation import Permutation + from sage.matrix.constructor import matrix + from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing + + poly_ring = PolynomialRing(self.base_ring(),var_name) + + m = self.nrows() + + # if no shift, this is equivalent to 0s + if shift is None: + shift = [0]*self.nrows() + + # # calculate shift priority function and permutation + # priority = lambda c,d : shift[c] + d + # index = lambda i : (i%m,i//m) + # index_inv = lambda c,d : c + d*m + + # priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) + # priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) + # priority_permutation_inv = priority_permutation.inverse() + + # # maps row c of EJ^d to position i in striped Krylov matrix + # # +/- 1 as permutations are 1-indexed, matrices are 0-indexed + # phi = lambda c,d : priority_permutation_inv(index_inv(c,d) + 1) - 1 + # phi_inv = lambda i : index(priority_permutation(i + 1) - 1) + + # calculate krylov profile + row_profile, pivot, col_profile = self.krylov_rank_profile(J,degree,shift,output_pairs=True) + + if len(row_profile) == 0: + return matrix.identity(poly_ring,m) + + # (c_k, d_k) = phi^-1 (row_i) + #c, d = zip(*(phi_inv(i) for i in row_profile)) + c, d = zip(*(row for row in row_profile)) + + # degree_c = 0 or max d[k] such that c[k] = c + inv_c = [None]*m + degree_c = [0]*m + for k in range(len(row_profile)): + if d[k] >= degree_c[c[k]]: + degree_c[c[k]] = d[k] + 1 + inv_c[c[k]] = k + + T_absent_indices = [i for i in range(m) if inv_c[i] is None] + T_present_indices = [inv_c[i] for i in range(m) if inv_c[i] is not None] + D_absent = self.matrix_from_rows_and_columns(T_absent_indices, col_profile) + D_present = (pivot.matrix_from_rows(T_present_indices)*J).matrix_from_columns(col_profile) + D_rows = [] + idx_p = 0 + idx_a = 0 + for i in range(m): + if inv_c[i] is None: + D_rows.append(D_absent[idx_a]) + idx_a += 1 + else: + D_rows.append(D_present[idx_p]) + idx_p += 1 + + C = pivot.matrix_from_columns(col_profile) + D = matrix(D_rows) + + relation = D*C.inverse() + + # linear interpolation basis in shifted Popov form + uncompressed_basis = matrix.block([[-relation,matrix.identity(m)]],subdivide=False) + + # construct variable + variable = poly_ring.gen() + + # compression of basis into polynomial form + basis_rows = [[0]*m for i in range(m)] + for col in range(relation.ncols()): + for row in range(m): + basis_rows[row][c[col]] += uncompressed_basis[row][col] * variable**d[col] + for i in range(m): + basis_rows[i][i] += variable**degree_c[i] + + return matrix(basis_rows) # a limited number of access-only properties are provided for matrices @property From 59dea5739548c3749dc219417229368667b1e203 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Thu, 3 Jul 2025 12:31:48 +0200 Subject: [PATCH 13/28] make optimisation final --- src/sage/matrix/matrix2.pyx | 96 +------------------------------------ 1 file changed, 1 insertion(+), 95 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index df4f389eedf..125fdeb1af3 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19387,102 +19387,8 @@ cdef class Matrix(Matrix1): row_profile_K = [phi(*row) for row in row_profile_self] return tuple(row_profile_K), M, col_profile - - def linear_interpolation_basis(self, J, degree, var_name, shift=None): - r""" - Construct a linear interpolant basis for (``self``,`J`) in `s`-Popov form. - - INPUT: - - - ``J`` -- square matrix for Krylov construction. - - ``degree`` -- upper bound on degree of minpoly(`J`), power of - - ``shift`` -- list of integers (optional): controls priority of rows. - OUTPUT: - - - Matrix whose rows form an interpolation basis. - - """ - from sage.combinat.permutation import Permutation - from sage.matrix.constructor import matrix - from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing - - poly_ring = PolynomialRing(self.base_ring(),var_name) - - m = self.nrows() - - # if no shift, this is equivalent to 0s - if shift is None: - shift = [0]*self.nrows() - - # calculate shift priority function and permutation - priority = lambda c,d : shift[c] + d - index = lambda i : (i%m,i//m) - index_inv = lambda c,d : c + d*m - - priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) - priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) - priority_permutation_inv = priority_permutation.inverse() - - # maps row c of EJ^d to position i in striped Krylov matrix - # +/- 1 as permutations are 1-indexed, matrices are 0-indexed - phi = lambda c,d : priority_permutation_inv(index_inv(c,d) + 1) - 1 - phi_inv = lambda i : index(priority_permutation(i + 1) - 1) - - # calculate krylov profile - row_profile, pivot, col_profile = self.krylov_rank_profile(J,degree,shift) - - if len(row_profile) == 0: - return matrix.identity(poly_ring,m) - - # (c_k, d_k) = phi^-1 (row_i) - c, d = zip(*(phi_inv(i) for i in row_profile)) - - # degree_c = 0 or max d[k] such that c[k] = c - inv_c = [None]*m - degree_c = [0]*m - for k in range(len(row_profile)): - if d[k] >= degree_c[c[k]]: - degree_c[c[k]] = d[k] + 1 - inv_c[c[k]] = k - - T_absent_indices = [i for i in range(m) if inv_c[i] is None] - T_present_indices = [inv_c[i] for i in range(m) if inv_c[i] is not None] - D_absent = self.matrix_from_rows_and_columns(T_absent_indices, col_profile) - D_present = (pivot.matrix_from_rows(T_present_indices)*J).matrix_from_columns(col_profile) - D_rows = [] - idx_p = 0 - idx_a = 0 - for i in range(m): - if inv_c[i] is None: - D_rows.append(D_absent[idx_a]) - idx_a += 1 - else: - D_rows.append(D_present[idx_p]) - idx_p += 1 - - C = pivot.matrix_from_columns(col_profile) - D = matrix(D_rows) - - relation = D*C.inverse() - - # linear interpolation basis in shifted Popov form - uncompressed_basis = matrix.block([[-relation,matrix.identity(m)]],subdivide=False) - - # construct variable - variable = poly_ring.gen() - - # compression of basis into polynomial form - basis_rows = [[0]*m for i in range(m)] - for col in range(relation.ncols()): - for row in range(m): - basis_rows[row][c[col]] += uncompressed_basis[row][col] * variable**d[col] - for i in range(m): - basis_rows[i][i] += variable**degree_c[i] - - return matrix(basis_rows) - - def linear_interpolation_basis_fast_perm(self, J, degree, var_name, shift=None): + def linear_interpolation_basis(self, J, degree, var_name, shift=None): r""" Construct a linear interpolant basis for (``self``,`J`) in `s`-Popov form. From ae04c03f4cd4cdba5c60d776f654385b8e097b07 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Fri, 4 Jul 2025 13:12:14 +0200 Subject: [PATCH 14/28] optimise matrix creation for basis, fix echelon_form bug --- src/sage/matrix/matrix2.pyx | 49 ++++++++++--------- .../matrix/matrix_modn_dense_template.pxi | 2 +- 2 files changed, 27 insertions(+), 24 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 125fdeb1af3..a6ef7f307ce 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19415,20 +19415,6 @@ cdef class Matrix(Matrix1): if shift is None: shift = [0]*self.nrows() - # # calculate shift priority function and permutation - # priority = lambda c,d : shift[c] + d - # index = lambda i : (i%m,i//m) - # index_inv = lambda c,d : c + d*m - - # priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) - # priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) - # priority_permutation_inv = priority_permutation.inverse() - - # # maps row c of EJ^d to position i in striped Krylov matrix - # # +/- 1 as permutations are 1-indexed, matrices are 0-indexed - # phi = lambda c,d : priority_permutation_inv(index_inv(c,d) + 1) - 1 - # phi_inv = lambda i : index(priority_permutation(i + 1) - 1) - # calculate krylov profile row_profile, pivot, col_profile = self.krylov_rank_profile(J,degree,shift,output_pairs=True) @@ -19467,19 +19453,36 @@ cdef class Matrix(Matrix1): relation = D*C.inverse() - # linear interpolation basis in shifted Popov form - uncompressed_basis = matrix.block([[-relation,matrix.identity(m)]],subdivide=False) + coeffs_map = {} + for i in range(m): + coeffs_map[i] = {} + for j in range(m): + coeffs_map[i][j] = {} - # construct variable - variable = poly_ring.gen() + for i in range(m): + coeffs_map[i][i][degree_c[i]] = poly_ring.base_ring().one() - # compression of basis into polynomial form - basis_rows = [[0]*m for i in range(m)] for col in range(relation.ncols()): for row in range(m): - basis_rows[row][c[col]] += uncompressed_basis[row][col] * variable**d[col] - for i in range(m): - basis_rows[i][i] += variable**degree_c[i] + coeffs_map[row][c[col]][d[col]] = coeffs_map[row][c[col]].get(d[col],poly_ring.base_ring().zero()) - relation[row][col] + + basis_rows = [[None for _ in range(m)] for _ in range(m)] + + monomial_cache = {} + for row in range(m): + for col in range(m): + if not coeffs_map[row][col]: + basis_rows[row][col] = poly_ring.zero() + elif len(coeffs_map[row][col]) == 1: + deg, coeff = coeffs_map[row][col].popitem() + if coeff not in monomial_cache: + monomial_cache[coeff] = poly_ring(coeff) + basis_rows[row][col] = monomial_cache[coeff].shift(deg) + else: + coeffs = [poly_ring.base_ring().zero()] * (max(coeffs_map[row][col].keys()) + 1) + for deg, coeff in coeffs_map[row][col].items(): + coeffs[deg] = coeff + basis_rows[row][col] = poly_ring(coeffs) return matrix(basis_rows) diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi index 63fac0af4b0..4655b3be180 100644 --- a/src/sage/matrix/matrix_modn_dense_template.pxi +++ b/src/sage/matrix/matrix_modn_dense_template.pxi @@ -176,7 +176,7 @@ cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nr Return the reduced row echelon form of this matrix. """ if linbox_is_zero(modulus, entries, nrows, ncols): - return 0, [] + return 0, [], [] cdef Py_ssize_t i, j cdef ModField *F = new ModField(modulus) From 796155e249576d91cfa18da43cb6e1de75076e53 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Tue, 8 Jul 2025 15:52:52 +0200 Subject: [PATCH 15/28] use matrix_from_rows in krylov_rank_profile, optimisation for polynomial creation in basis calculation --- src/sage/matrix/matrix2.pyx | 104 ++++++++++++++++++++++++------------ 1 file changed, 69 insertions(+), 35 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index a6ef7f307ce..5f8ad5fb952 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19349,7 +19349,7 @@ cdef class Matrix(Matrix1): M = self_permuted.matrix_from_rows(row_profile_self_permuted) J_L = None - + for l in range(math.ceil(math.log(degree,2)) + 1): L = pow(2,l) # adding 2^l to each degree @@ -19366,16 +19366,19 @@ cdef class Matrix(Matrix1): J_L = J else: J_L = J_L * J_L + M = M.stack(M*J_L,subdivide=False) # sort rows of M, find profile, translate to k (indices of full krylov matrix) M.permute_rows(k_perm) + row_profile_M = M.row_rank_profile() + r = len(row_profile_M) row_profile_self = sorted([k[k_perm(i+1)-1] for i in row_profile_M],key=lambda x: (shift[x[0]]+x[1],x[0])) # calculate new M for return value or next loop - M = matrix([M.row(i) for i in row_profile_M]) + M = M.matrix_from_rows(row_profile_M) col_profile = M.col_rank_profile() @@ -19383,12 +19386,12 @@ cdef class Matrix(Matrix1): return tuple(row_profile_self), M, col_profile # convert c,d to actual position in striped Krylov matrix - phi = lambda c,d : sum(min(max(shift[c] - shift[i] + d + (i < c and shift[i] <= shift[c] + degree),0),degree+1) for i in range(len(shift))) + phi = lambda c,d : sum(min(max(shift[c] - shift[i] + d + (i < c and shift[i] <= shift[c] + degree),0),degree+1) for i in range(m)) row_profile_K = [phi(*row) for row in row_profile_self] return tuple(row_profile_K), M, col_profile - def linear_interpolation_basis(self, J, degree, var_name, shift=None): + def linear_interpolation_basis(self, J, degree, var_name, shift=None, polynomial_output=True): r""" Construct a linear interpolant basis for (``self``,`J`) in `s`-Popov form. @@ -19449,42 +19452,73 @@ cdef class Matrix(Matrix1): idx_p += 1 C = pivot.matrix_from_columns(col_profile) + D = matrix(D_rows) relation = D*C.inverse() - coeffs_map = {} - for i in range(m): - coeffs_map[i] = {} - for j in range(m): - coeffs_map[i][j] = {} - - for i in range(m): - coeffs_map[i][i][degree_c[i]] = poly_ring.base_ring().one() - - for col in range(relation.ncols()): - for row in range(m): - coeffs_map[row][c[col]][d[col]] = coeffs_map[row][c[col]].get(d[col],poly_ring.base_ring().zero()) - relation[row][col] - - basis_rows = [[None for _ in range(m)] for _ in range(m)] - - monomial_cache = {} - for row in range(m): - for col in range(m): - if not coeffs_map[row][col]: - basis_rows[row][col] = poly_ring.zero() - elif len(coeffs_map[row][col]) == 1: - deg, coeff = coeffs_map[row][col].popitem() - if coeff not in monomial_cache: - monomial_cache[coeff] = poly_ring(coeff) - basis_rows[row][col] = monomial_cache[coeff].shift(deg) - else: - coeffs = [poly_ring.base_ring().zero()] * (max(coeffs_map[row][col].keys()) + 1) - for deg, coeff in coeffs_map[row][col].items(): - coeffs[deg] = coeff - basis_rows[row][col] = poly_ring(coeffs) + if polynomial_output: + # construct polynomial matrix - potentially costly at polynomial construction stage for extension fields + coeffs_map = [[{} for _ in range(m)] for _ in range(m)] + # add identity part of basis + for i in range(m): + coeffs_map[i][i][degree_c[i]] = poly_ring.base_ring().one() + # add relation part of basis + for col in range(relation.ncols()): + for row in range(m): + coeffs_map[row][c[col]][d[col]] = coeffs_map[row][c[col]].get(d[col],poly_ring.base_ring().zero()) - relation[row,col] + # construct rows + basis_rows = [[None for _ in range(m)] for _ in range(m)] + if poly_ring.base_ring().degree() <= 1 or m*m > poly_ring.base_ring().order(): + # if base field is large, don't bother caching monomials + # or if polynomial construction is efficient (order = 1) + for row in range(m): + for col in range(m): + if not coeffs_map[row][col]: + basis_rows[row][col] = poly_ring.base_ring().zero() + else: + # construct a small polynomial (less costly), then shift + mindeg = min(coeffs_map[row][col].keys()) + maxdeg = max(coeffs_map[row][col].keys()) + entries = [poly_ring.base_ring().zero()] * (maxdeg-mindeg+1) + for deg, coeff in coeffs_map[row][col].items(): + entries[deg-mindeg] = coeff + basis_rows[row][col] = poly_ring(entries).shift(mindeg) + else: + # if base field is small, cache monomials to minimise number of constructions + monomial_cache = {} + for row in range(m): + for col in range(m): + if not coeffs_map[row][col]: + basis_rows[row][col] = poly_ring.zero() + elif len(coeffs_map[row][col]) == 1: + # monomial case + deg, coeff = coeffs_map[row][col].popitem() + if coeff not in monomial_cache: + monomial_cache[coeff] = poly_ring(coeff) + basis_rows[row][col] = monomial_cache[coeff].shift(deg) + else: + # dense case + mindeg = min(coeffs_map[row][col].keys()) + maxdeg = max(coeffs_map[row][col].keys()) + entries = [poly_ring.base_ring().zero()] * (maxdeg-mindeg+1) + for deg, coeff in coeffs_map[row][col].items(): + entries[deg-mindeg] = coeff + basis_rows[row][col] = poly_ring(entries).shift(mindeg) + # convert to matrix + output = matrix(poly_ring,basis_rows) + else: + # expected output is expanded matrix + row_width = m * (max(degree_c) + 1) + basis_rows = [[poly_ring.base_ring().zero()]*row_width for i in range(m)] + for i in range(m): + basis_rows[i][i+degree_c[i]*m] = poly_ring.base_ring().one() + for col in range(relation.ncols()): + for row in range(m): + basis_rows[row][c[col]+d[col]*m] -= relation[row,col] + output = matrix(poly_ring.base_ring(),basis_rows) - return matrix(basis_rows) + return output # a limited number of access-only properties are provided for matrices @property From a6de9529bf3a25a259eb0914a7b7d6589c62e00d Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Tue, 8 Jul 2025 16:44:37 +0200 Subject: [PATCH 16/28] add early exit filtering for krylov rank profile --- src/sage/matrix/matrix2.pyx | 180 ++++++++++++++++++++++++++++++++++++ 1 file changed, 180 insertions(+) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 5f8ad5fb952..2d6de67ce9a 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19391,6 +19391,186 @@ cdef class Matrix(Matrix1): return tuple(row_profile_K), M, col_profile + def krylov_rank_profile_early_exit(self, J, degree, shift=None, output_pairs=False): + r""" + Compute the rank profile (row and column) of the striped Krylov matrix + built from ``self`` and matrix `J`. + + INPUT: + + - ``J`` -- square matrix used in the Krylov construction. + - ``degree`` -- maximum power of J in Krylov matrix. + - ``shift`` -- list of integers (optional): priority row shift. + + OUTPUT: + + - A tuple (row_profile, pivot_matrix, column_profile): + * row_profile: list of the first r row indices in the striped Krylov matrix `K` corresponding to independent rows + * pivot: the submatrix of `K` given by those rows + * column_profile: list of the first r independent column indices in ``pivot_matrix`` + where r is the rank of `K`. + + EXAMPLES: + + sage: R. = GF(97)[] + sage: E = matrix([[27,49,29],[50,58,0],[77,10,29]]) + sage: E + [27 49 29] + [50 58 0] + [77 10 29] + sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) + sage: J + [0 1 0] + [0 0 1] + [0 0 0] + sage: E.striped_krylov_matrix(J,3) + [27 49 29] + [50 58 0] + [77 10 29] + [ 0 27 49] + [ 0 50 58] + [ 0 77 10] + [ 0 0 27] + [ 0 0 50] + [ 0 0 77] + [ 0 0 0] + [ 0 0 0] + [ 0 0 0] + sage: E.krylov_rank_profile(J,3) + ( + [27 49 29] + [50 58 0] + (0, 1, 3), [ 0 27 49], (0, 1, 2) + ) + sage: E.striped_krylov_matrix(J,3,[0,3,6]) + [27 49 29] + [ 0 27 49] + [ 0 0 27] + [ 0 0 0] + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [ 0 0 0] + [77 10 29] + [ 0 77 10] + [ 0 0 77] + [ 0 0 0] + sage: E.krylov_rank_profile(J,3,[0,3,6]) + ( + [27 49 29] + [ 0 27 49] + (0, 1, 2), [ 0 0 27], (0, 1, 2) + ) + sage: E.striped_krylov_matrix(J,3,[3,0,2]) + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [77 10 29] + [27 49 29] + [ 0 0 0] + [ 0 77 10] + [ 0 27 49] + [ 0 0 77] + [ 0 0 27] + [ 0 0 0] + [ 0 0 0] + sage: E.krylov_rank_profile(J,3,[3,0,2]) + ( + [50 58 0] + [ 0 50 58] + (0, 1, 2), [ 0 0 50], (0, 1, 2) + ) + """ + from sage.matrix.constructor import matrix + import math + from sage.combinat.permutation import Permutation + + m = self.nrows() + + if m == 0: + return (),self,() + + # calculate row profile of self, with shift applied + self_permutation = Permutation(sorted([i+1 for i in range(m)],key=lambda x:(shift[x-1],x-1))) + + self_permuted = self.with_permuted_rows(self_permutation) + row_profile_self_permuted = self_permuted.row_rank_profile() + row_profile_self = [(self_permutation(i+1)-1,0) for i in row_profile_self_permuted] + + exhausted = matrix.zero(self.base_ring(), 0, self.ncols()) + row_profile_exhausted = [] + excluded_rows = set() + + r = len(row_profile_self) + + if r == 0: + return (),self.matrix_from_rows([]),() + + M = self_permuted.matrix_from_rows(row_profile_self_permuted) + + J_L = None + + for l in range(math.ceil(math.log(degree,2)) + 1): + L = pow(2,l) + # adding 2^l to each degree + row_extension = [(x[0],x[1] + L) for x in row_profile_self if x[1] + L <= degree] + if len(row_extension) == 0: + break + + # concatenate two sequences (order preserved) + k = row_profile_exhausted + row_profile_self + row_extension + + # calculate sorting permutation, sort k by (shift[c]+d, c) + k_perm = Permutation(sorted([i+1 for i in range(len(k))],key=lambda x: (shift[k[x-1][0]] + k[x-1][1],k[x-1][0]))) + + # fast calculation of rows formed by indices in k + if J_L is None: + J_L = J + else: + J_L = J_L * J_L + + M = matrix.block([[exhausted],[M],[M*J_L]],subdivide=False) + + # sort rows of M, find profile, translate to k (indices of full krylov matrix) + M.permute_rows(k_perm) + + row_profile_M = M.row_rank_profile() + r = len(row_profile_M) + + if r == self.ncols(): + tail = list(range(row_profile_M[-1]+1,M.nrows())) + excluded_rows.update(set([k[k_perm(i+1)-1][0] for i in tail if i < len(k)])) + + xmi = [i for i in row_profile_M if k[k_perm(i+1)-1][0] in excluded_rows] + imi = [i for i in row_profile_M if k[k_perm(i+1)-1][0] not in excluded_rows] + + row_profile_exhausted = [k[k_perm(i+1)-1] for i in xmi] + row_profile_self = [k[k_perm(i+1)-1] for i in imi] + + exhausted = M.matrix_from_rows(xmi) + M = M.matrix_from_rows(imi) + else: + row_profile_self = [k[k_perm(i+1)-1] for i in row_profile_M] + # calculate new M for return value or next loop + M = M.matrix_from_rows(row_profile_M) + if exhausted.nrows() != 0: + k = row_profile_exhausted + row_profile_self + M = exhausted.stack(M) + k_perm = Permutation(sorted([i+1 for i in range(len(k))],key=lambda x: (shift[k[x-1][0]] + k[x-1][1],k[x-1][0]))) + + M.permute_rows(k_perm) + row_profile_self = [k[k_perm(i+1)-1] for i in range(len(k))] + col_profile = M.col_rank_profile() + + if output_pairs: + return tuple(row_profile_self), M, col_profile + + # convert c,d to actual position in striped Krylov matrix + phi = lambda c,d : sum(min(max(shift[c] - shift[i] + d + (i < c and shift[i] <= shift[c] + degree),0),degree+1) for i in range(m)) + row_profile_K = [phi(*row) for row in row_profile_self] + + return tuple(row_profile_K), M, col_profile + def linear_interpolation_basis(self, J, degree, var_name, shift=None, polynomial_output=True): r""" Construct a linear interpolant basis for (``self``,`J`) in `s`-Popov form. From a010078eef0040fcd1a9116c1ef55f766c09d83a Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Wed, 9 Jul 2025 13:55:09 +0200 Subject: [PATCH 17/28] mitigate performance loss on hermite shift --- src/sage/matrix/matrix2.pyx | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 2d6de67ce9a..4df0bc16bd8 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19537,7 +19537,7 @@ cdef class Matrix(Matrix1): row_profile_M = M.row_rank_profile() r = len(row_profile_M) - if r == self.ncols(): + if r == self.ncols() and l < math.ceil(math.log(degree,2)): tail = list(range(row_profile_M[-1]+1,M.nrows())) excluded_rows.update(set([k[k_perm(i+1)-1][0] for i in tail if i < len(k)])) @@ -19550,10 +19550,16 @@ cdef class Matrix(Matrix1): exhausted = M.matrix_from_rows(xmi) M = M.matrix_from_rows(imi) else: + if l == math.ceil(math.log(degree,2)): + row_profile_exhausted = [] + exhausted = matrix.zero(self.base_ring(), 0, self.ncols()) row_profile_self = [k[k_perm(i+1)-1] for i in row_profile_M] # calculate new M for return value or next loop M = M.matrix_from_rows(row_profile_M) - if exhausted.nrows() != 0: + if M.nrows() == 0: + M = exhausted + row_profile_self = row_profile_exhausted + elif exhausted.nrows() != 0: k = row_profile_exhausted + row_profile_self M = exhausted.stack(M) k_perm = Permutation(sorted([i+1 for i in range(len(k))],key=lambda x: (shift[k[x-1][0]] + k[x-1][1],k[x-1][0]))) From b6df7d32c3b1fa72a49fb3154fdef1a327d40247 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Thu, 17 Jul 2025 20:07:08 +0200 Subject: [PATCH 18/28] make early exit permanent --- src/sage/matrix/matrix2.pyx | 158 +----------------------------------- 1 file changed, 1 insertion(+), 157 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 13303ce0fbd..6dbc119ed8e 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19319,164 +19319,8 @@ cdef class Matrix(Matrix1): col_profile = pivot.col_rank_profile() return row_profile,pivot,col_profile - - def krylov_rank_profile(self, J, degree, shift=None, output_pairs=False): - r""" - Compute the rank profile (row and column) of the striped Krylov matrix - built from ``self`` and matrix `J`. - - INPUT: - - - ``J`` -- square matrix used in the Krylov construction. - - ``degree`` -- maximum power of J in Krylov matrix. - - ``shift`` -- list of integers (optional): priority row shift. - - OUTPUT: - - - A tuple (row_profile, pivot_matrix, column_profile): - * row_profile: list of the first r row indices in the striped Krylov matrix `K` corresponding to independent rows - * pivot: the submatrix of `K` given by those rows - * column_profile: list of the first r independent column indices in ``pivot_matrix`` - where r is the rank of `K`. - - EXAMPLES: - - sage: R. = GF(97)[] - sage: E = matrix([[27,49,29],[50,58,0],[77,10,29]]) - sage: E - [27 49 29] - [50 58 0] - [77 10 29] - sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) - sage: J - [0 1 0] - [0 0 1] - [0 0 0] - sage: E.striped_krylov_matrix(J,3) - [27 49 29] - [50 58 0] - [77 10 29] - [ 0 27 49] - [ 0 50 58] - [ 0 77 10] - [ 0 0 27] - [ 0 0 50] - [ 0 0 77] - [ 0 0 0] - [ 0 0 0] - [ 0 0 0] - sage: E.krylov_rank_profile(J,3) - ( - [27 49 29] - [50 58 0] - (0, 1, 3), [ 0 27 49], (0, 1, 2) - ) - sage: E.striped_krylov_matrix(J,3,[0,3,6]) - [27 49 29] - [ 0 27 49] - [ 0 0 27] - [ 0 0 0] - [50 58 0] - [ 0 50 58] - [ 0 0 50] - [ 0 0 0] - [77 10 29] - [ 0 77 10] - [ 0 0 77] - [ 0 0 0] - sage: E.krylov_rank_profile(J,3,[0,3,6]) - ( - [27 49 29] - [ 0 27 49] - (0, 1, 2), [ 0 0 27], (0, 1, 2) - ) - sage: E.striped_krylov_matrix(J,3,[3,0,2]) - [50 58 0] - [ 0 50 58] - [ 0 0 50] - [77 10 29] - [27 49 29] - [ 0 0 0] - [ 0 77 10] - [ 0 27 49] - [ 0 0 77] - [ 0 0 27] - [ 0 0 0] - [ 0 0 0] - sage: E.krylov_rank_profile(J,3,[3,0,2]) - ( - [50 58 0] - [ 0 50 58] - (0, 1, 2), [ 0 0 50], (0, 1, 2) - ) - """ - from sage.matrix.constructor import matrix - import math - from sage.combinat.permutation import Permutation - - m = self.nrows() - - if m == 0: - return (),self,() - - # calculate row profile of self, with shift applied - self_permutation = Permutation(sorted([i+1 for i in range(m)],key=lambda x:(shift[x-1],x-1))) - - self_permuted = self.with_permuted_rows(self_permutation) - row_profile_self_permuted = self_permuted.row_rank_profile() - row_profile_self = [(self_permutation(i+1)-1,0) for i in row_profile_self_permuted] - - r = len(row_profile_self) - - if r == 0: - return (),self.matrix_from_rows([]),() - - M = self_permuted.matrix_from_rows(row_profile_self_permuted) - - J_L = None - - for l in range(math.ceil(math.log(degree,2)) + 1): - L = pow(2,l) - # adding 2^l to each degree - row_extension = [(x[0],x[1] + L) for x in row_profile_self if x[1] + L <= degree] - - # concatenate two sequences (order preserved) - k = row_profile_self + row_extension - - # calculate sorting permutation, sort k by (shift[c]+d, c) - k_perm = Permutation(sorted([i+1 for i in range(len(k))],key=lambda x: (shift[k[x-1][0]] + k[x-1][1],k[x-1][0]))) - - # fast calculation of rows formed by indices in k - if J_L is None: - J_L = J - else: - J_L = J_L * J_L - - M = M.stack(M*J_L,subdivide=False) - - # sort rows of M, find profile, translate to k (indices of full krylov matrix) - M.permute_rows(k_perm) - - row_profile_M = M.row_rank_profile() - - r = len(row_profile_M) - row_profile_self = sorted([k[k_perm(i+1)-1] for i in row_profile_M],key=lambda x: (shift[x[0]]+x[1],x[0])) - - # calculate new M for return value or next loop - M = M.matrix_from_rows(row_profile_M) - - col_profile = M.col_rank_profile() - - if output_pairs: - return tuple(row_profile_self), M, col_profile - - # convert c,d to actual position in striped Krylov matrix - phi = lambda c,d : sum(min(max(shift[c] - shift[i] + d + (i < c and shift[i] <= shift[c] + degree),0),degree+1) for i in range(m)) - row_profile_K = [phi(*row) for row in row_profile_self] - - return tuple(row_profile_K), M, col_profile - def krylov_rank_profile_early_exit(self, J, degree, shift=None, output_pairs=False): + def krylov_rank_profile(self, J, degree, shift=None, output_pairs=False): r""" Compute the rank profile (row and column) of the striped Krylov matrix built from ``self`` and matrix `J`. From 661db965d4dbd803a7df1bbe59ff5a4cd115cd8a Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Mon, 21 Jul 2025 12:18:46 +0200 Subject: [PATCH 19/28] fix boundary conditions for polynomial conversion --- src/sage/matrix/matrix2.pyx | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 6dbc119ed8e..64c5431809a 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19584,7 +19584,7 @@ cdef class Matrix(Matrix1): coeffs_map[row][c[col]][d[col]] = coeffs_map[row][c[col]].get(d[col],poly_ring.base_ring().zero()) - relation[row,col] # construct rows basis_rows = [[None for _ in range(m)] for _ in range(m)] - if poly_ring.base_ring().degree() <= 1 or m*m > poly_ring.base_ring().order(): + if poly_ring.base_ring().degree() <= 1 or m*m*self.ncols() < poly_ring.base_ring().order(): # if base field is large, don't bother caching monomials # or if polynomial construction is efficient (order = 1) for row in range(m): @@ -19613,13 +19613,19 @@ cdef class Matrix(Matrix1): monomial_cache[coeff] = poly_ring(coeff) basis_rows[row][col] = monomial_cache[coeff].shift(deg) else: - # dense case + # general case mindeg = min(coeffs_map[row][col].keys()) - maxdeg = max(coeffs_map[row][col].keys()) - entries = [poly_ring.base_ring().zero()] * (maxdeg-mindeg+1) + # maxdeg = max(coeffs_map[row][col].keys()) + poly = poly_ring.zero() for deg, coeff in coeffs_map[row][col].items(): - entries[deg-mindeg] = coeff - basis_rows[row][col] = poly_ring(entries).shift(mindeg) + if coeff not in monomial_cache: + monomial_cache[coeff] = poly_ring(coeff) + poly += monomial_cache[coeff].shift(deg-mindeg) + basis_rows[row][col] = poly.shift(mindeg) + # entries = [poly_ring.base_ring().zero()] * (maxdeg-mindeg+1) + # for deg, coeff in coeffs_map[row][col].items(): + # entries[deg-mindeg] = coeff + # basis_rows[row][col] = poly_ring(entries).shift(mindeg) # convert to matrix output = matrix(poly_ring,basis_rows) else: From 01f9301c8f0ab9d78e57528e61621703e71b68fa Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Mon, 28 Jul 2025 15:41:19 +0200 Subject: [PATCH 20/28] tidy up code --- src/sage/matrix/matrix2.pyx | 346 ++++++++++++++++-------------------- 1 file changed, 152 insertions(+), 194 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 64c5431809a..fc101b3000a 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19040,60 +19040,11 @@ cdef class Matrix(Matrix1): U.rescale_col(n - 1, -1) return U - def polynomial_compression(self,degree,var_name): - """ - Returns the corresponding polynomial matrix where the coefficient matrix of degree i is the ith block of ``self``. - - INPUT: - - - ``degree`` -- integer: the number of blocks in the expanded matrix. If ``self`` does not have zero-blocks, - this corresponds to the degree of the resulting matrix. - - ``variable`` -- a symbolic variable (e.g., ``x``) to use in the polynomial. - - OUTPUT: - - - a polynomial matrix with the same number of rows and self.ncols()//(degree+1) columns - - EXAMPLES: - - The first 3 x 3 block of N corresponds to the - - sage: R. = GF(97)[] - sage: S. = ZZ[] - sage: M = matrix([[x^2+36*x,31*x,0],[3*x+13,x+57,0],[96,96,1]]) - sage: M - [x^2 + 36*x 31*x 0] - [ 3*x + 13 x + 57 0] - [ 96 96 1] - sage: N = M.expansion() - sage: N - [ 0 0 0 36 31 0 1 0 0] - [13 57 0 3 1 0 0 0 0] - [96 96 1 0 0 0 0 0 0] - sage: N.polynomial_compression(2,x) - [x^2 + 36*x 31*x 0] - [ 3*x + 13 x + 57 0] - [ 96 96 1] - sage: N.polynomial_compression(2,y) - [y^2 + 36*y 31*y 0] - [ 3*y + 13 y + 57 0] - [ 96 96 1] - """ - from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing - - poly_ring = PolynomialRing(self.base_ring(),var_name) - variable = poly_ring.gen() - - if self.ncols() % (degree+1) != 0: - raise ValueError('The column number must be divisible by degree+1.') - coeff_matrices = [self[:,i*(self.ncols()//(degree+1)):(i+1)*(self.ncols()//(degree+1))] for i in range(degree+1)] - return sum([coeff_matrices[i]*variable**i for i in range(degree+1)]) - def striped_krylov_matrix(self, J, degree, shift=None): r""" Return the block matrix of the following form, with rows permuted according to shift. The following uses `E` to refer to ``self``, and `d` to refer to ``degree``. - + [ ] [ E ] [ ] @@ -19109,20 +19060,20 @@ cdef class Matrix(Matrix1): [ ] [ EJ^d] [ ] - + INPUT: - + - ``J`` -- a square matrix of size equal to the number of columns of ``self``. - ``degree`` -- integer, the maximum exponent for the Krylov matrix. - ``shift`` -- list of integers (optional), row priority shifts. If ``None``, defaults to all zero. - + OUTPUT: - A matrix with block rows [E, EJ, EJ^2, ..., EJ^d], row-permuted by shift. - + EXAMPLES: - + sage: R. = GF(97)[] sage: E = matrix([[27,49,29],[50,58,0],[77,10,29]]) sage: E @@ -19173,48 +19124,48 @@ cdef class Matrix(Matrix1): [ 0 0 27] [ 0 0 0] [ 0 0 0] - + """ from sage.combinat.permutation import Permutation from sage.matrix.constructor import matrix # for matrix.block - + m = self.nrows() - + # if no shift, this is equivalent to 0s if shift is None: shift = [0]*self.nrows() - + # size checks if len(shift) != m: raise ValueError(f"The shift should have the same number of elements as the rows of the matrix ({m}), but had {len(shift)}.") if not J.is_square() or J.nrows() != self.ncols(): raise ValueError(f"The matrix J should be a square matrix and match the number of columns of self ({self.ncols()}), but is of dimension {J.nrows()} x {J.ncols()}.") - + # define priority and indexing functions # index: [0,1,...,m-1,m,...,2m-1,...,m*degree-1] -> [(0,0),(1,0),...,(m-1,0),(0,1),...,(m-1,1),...,(m-1,degree)] priority = lambda c,d : shift[c] + d index = lambda i : (i%m,i//m) - + # deduce permutation by sorting priorities # rows should be ordered by ascending priority, then by associated row number in E (i%m) priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) - + # build all blocks for the striped matrix blocks = [] - + # for i from 0 to degree (inclusive), add E*J^i J_i = matrix.identity(self.ncols()) for i in range(degree+1): blocks.append([self*J_i]) if i < degree: J_i *= J - + # return block matrix permuted according to shift krylov = matrix.block(blocks,subdivide=False) krylov.permute_rows(priority_permutation) return krylov - + def naive_krylov_rank_profile(self, J, degree, shift=None): r""" Compute the rank profile (row and column) of the striped Krylov matrix @@ -19233,9 +19184,9 @@ cdef class Matrix(Matrix1): * pivot: the submatrix of `K` given by those rows * column_profile: list of the first r independent column indices in ``pivot_matrix`` where r is the rank of `K`. - + EXAMPLES: - + sage: R. = GF(97)[] sage: E = matrix([[27,49,29],[50,58,0],[77,10,29]]) sage: E @@ -19306,21 +19257,21 @@ cdef class Matrix(Matrix1): ) """ from sage.matrix.constructor import matrix - + K = self.striped_krylov_matrix(J,degree,shift) - + # calculate first independent rows of K row_profile = K.row_rank_profile() - + # construct submatrix pivot = K.matrix_from_rows(row_profile) - + # calculate first independent columns of pivot col_profile = pivot.col_rank_profile() - + return row_profile,pivot,col_profile - def krylov_rank_profile(self, J, degree, shift=None, output_pairs=False): + def krylov_rank_profile(self, J, degree=None, shift=None, output_pairs=False): r""" Compute the rank profile (row and column) of the striped Krylov matrix built from ``self`` and matrix `J`. @@ -19339,10 +19290,10 @@ cdef class Matrix(Matrix1): * column_profile: list of the first r independent column indices in ``pivot_matrix`` where r is the rank of `K`. - EXAMPLES: + TESTS:: - sage: R. = GF(97)[] - sage: E = matrix([[27,49,29],[50,58,0],[77,10,29]]) + sage: R = GF(97) + sage: E = matrix(R,[[27,49,29],[50,58,0],[77,10,29]]) sage: E [27 49 29] [50 58 0] @@ -19352,70 +19303,71 @@ cdef class Matrix(Matrix1): [0 1 0] [0 0 1] [0 0 0] - sage: E.striped_krylov_matrix(J,3) - [27 49 29] - [50 58 0] - [77 10 29] - [ 0 27 49] - [ 0 50 58] - [ 0 77 10] - [ 0 0 27] - [ 0 0 50] - [ 0 0 77] - [ 0 0 0] - [ 0 0 0] - [ 0 0 0] sage: E.krylov_rank_profile(J,3) ( [27 49 29] [50 58 0] (0, 1, 3), [ 0 27 49], (0, 1, 2) ) - sage: E.striped_krylov_matrix(J,3,[0,3,6]) - [27 49 29] - [ 0 27 49] - [ 0 0 27] - [ 0 0 0] - [50 58 0] - [ 0 50 58] - [ 0 0 50] - [ 0 0 0] - [77 10 29] - [ 0 77 10] - [ 0 0 77] - [ 0 0 0] + sage: E.krylov_rank_profile(J,3,None,True) + ( + [27 49 29] + [50 58 0] + ((0, 0), (1, 0), (0, 1)), [ 0 27 49], (0, 1, 2) + ) sage: E.krylov_rank_profile(J,3,[0,3,6]) ( [27 49 29] [ 0 27 49] (0, 1, 2), [ 0 0 27], (0, 1, 2) ) - sage: E.striped_krylov_matrix(J,3,[3,0,2]) - [50 58 0] - [ 0 50 58] - [ 0 0 50] - [77 10 29] - [27 49 29] - [ 0 0 0] - [ 0 77 10] - [ 0 27 49] - [ 0 0 77] - [ 0 0 27] - [ 0 0 0] - [ 0 0 0] + sage: E.krylov_rank_profile(J,3,[0,3,6],True) + ( + [27 49 29] + [ 0 27 49] + ((0, 0), (0, 1), (0, 2)), [ 0 0 27], (0, 1, 2) + ) sage: E.krylov_rank_profile(J,3,[3,0,2]) ( [50 58 0] [ 0 50 58] (0, 1, 2), [ 0 0 50], (0, 1, 2) ) + sage: E.krylov_rank_profile(J,3,[3,0,2],True) + ( + [50 58 0] + [ 0 50 58] + ((1, 0), (1, 1), (1, 2)), [ 0 0 50], (0, 1, 2) + ) """ from sage.matrix.constructor import matrix import math from sage.combinat.permutation import Permutation - + + if not isinstance(J, Matrix): + raise TypeError("krylov_rank_profile: J is not a matrix") + if J.nrows() != self.ncols() or J.ncols() != self.ncols(): + raise ValueError("krylov_rank_profile: matrix J does not have correct dimensions.") + if J.base_ring() != self.base_ring(): + raise ValueError("krylov_rank_profile: matrix J does not have same base ring as E.") + + if degree is None: + degree = self.ncols() + if not isinstance(degree, (int, sage.rings.integer.Integer)): + raise TypeError("krylov_rank_profile: degree is not an integer.") + if degree < 0: + raise ValueError("krylov_rank_profile: degree bound cannot be negative.") + + if shift is None or shift == 0: + shift = (ZZ**self.nrows()).zero() + if isinstance(shift, (list, tuple)): + shift = (ZZ**self.nrows())(shift) + if shift not in ZZ**self.nrows(): + raise ValueError(f"krylov_rank_profile: shift is not a Z-vector of length {self.nrows()}.") + m = self.nrows() - + sigma = self.ncols() + if m == 0: return (),self,() @@ -19426,7 +19378,7 @@ cdef class Matrix(Matrix1): row_profile_self_permuted = self_permuted.row_rank_profile() row_profile_self = [(self_permutation(i+1)-1,0) for i in row_profile_self_permuted] - exhausted = matrix.zero(self.base_ring(), 0, self.ncols()) + exhausted = matrix.zero(self.base_ring(), 0, sigma) row_profile_exhausted = [] excluded_rows = set() @@ -19466,7 +19418,7 @@ cdef class Matrix(Matrix1): row_profile_M = M.row_rank_profile() r = len(row_profile_M) - if r == self.ncols() and l < math.ceil(math.log(degree,2)): + if r == sigma and l < math.ceil(math.log(degree,2)): tail = list(range(row_profile_M[-1]+1,M.nrows())) excluded_rows.update(set([k[k_perm(i+1)-1][0] for i in tail if i < len(k)])) @@ -19481,7 +19433,7 @@ cdef class Matrix(Matrix1): else: if l == math.ceil(math.log(degree,2)): row_profile_exhausted = [] - exhausted = matrix.zero(self.base_ring(), 0, self.ncols()) + exhausted = matrix.zero(self.base_ring(), 0, sigma) row_profile_self = [k[k_perm(i+1)-1] for i in row_profile_M] # calculate new M for return value or next loop M = M.matrix_from_rows(row_profile_M) @@ -19506,13 +19458,14 @@ cdef class Matrix(Matrix1): return tuple(row_profile_K), M, col_profile - def linear_interpolation_basis(self, J, degree, var_name, shift=None, polynomial_output=True): + def linear_interpolation_basis(self, J, var_name, degree=None, shift=None): r""" Construct a linear interpolant basis for (``self``,`J`) in `s`-Popov form. INPUT: - ``J`` -- square matrix for Krylov construction. + - ``var_name`` -- variable name for the returned polynomial matrix - ``degree`` -- upper bound on degree of minpoly(`J`), power of - ``shift`` -- list of integers (optional): controls priority of rows. @@ -19524,18 +19477,40 @@ cdef class Matrix(Matrix1): from sage.combinat.permutation import Permutation from sage.matrix.constructor import matrix from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing - - poly_ring = PolynomialRing(self.base_ring(),var_name) - + + # INPUT VALIDATION + if not isinstance(J, Matrix): + raise TypeError("krylov_rank_profile: J is not a matrix") + if J.nrows() != self.ncols() or J.ncols() != self.ncols(): + raise ValueError("krylov_rank_profile: matrix J does not have correct dimensions.") + if J.base_ring() != self.base_ring(): + raise ValueError("krylov_rank_profile: matrix J does not have same base ring as E.") + + if not isinstance(var_name, str): + raise TypeError("var_name is not a string") + + if degree is None: + degree = self.ncols() + if not isinstance(degree, (int, sage.rings.integer.Integer)): + raise TypeError("krylov_rank_profile: degree is not an integer.") + if degree < 0: + raise ValueError("krylov_rank_profile: degree bound cannot be negative.") + + if shift is None or shift == 0: + shift = (ZZ**self.nrows()).zero() + if isinstance(shift, (list, tuple)): + shift = (ZZ**self.nrows())(shift) + if shift not in ZZ**self.nrows(): + raise ValueError(f"krylov_rank_profile: shift is not a Z-vector of length {self.nrows()}.") + m = self.nrows() + sigma = self.ncols() - # if no shift, this is equivalent to 0s - if shift is None: - shift = [0]*self.nrows() + poly_ring = PolynomialRing(self.base_ring(),var_name) # calculate krylov profile row_profile, pivot, col_profile = self.krylov_rank_profile(J,degree,shift,output_pairs=True) - + if len(row_profile) == 0: return matrix.identity(poly_ring,m) @@ -19572,74 +19547,57 @@ cdef class Matrix(Matrix1): relation = D*C.inverse() - if polynomial_output: - # construct polynomial matrix - potentially costly at polynomial construction stage for extension fields - coeffs_map = [[{} for _ in range(m)] for _ in range(m)] - # add identity part of basis - for i in range(m): - coeffs_map[i][i][degree_c[i]] = poly_ring.base_ring().one() - # add relation part of basis - for col in range(relation.ncols()): - for row in range(m): - coeffs_map[row][c[col]][d[col]] = coeffs_map[row][c[col]].get(d[col],poly_ring.base_ring().zero()) - relation[row,col] - # construct rows - basis_rows = [[None for _ in range(m)] for _ in range(m)] - if poly_ring.base_ring().degree() <= 1 or m*m*self.ncols() < poly_ring.base_ring().order(): - # if base field is large, don't bother caching monomials - # or if polynomial construction is efficient (order = 1) - for row in range(m): - for col in range(m): - if not coeffs_map[row][col]: - basis_rows[row][col] = poly_ring.base_ring().zero() - else: - # construct a small polynomial (less costly), then shift - mindeg = min(coeffs_map[row][col].keys()) - maxdeg = max(coeffs_map[row][col].keys()) - entries = [poly_ring.base_ring().zero()] * (maxdeg-mindeg+1) - for deg, coeff in coeffs_map[row][col].items(): - entries[deg-mindeg] = coeff - basis_rows[row][col] = poly_ring(entries).shift(mindeg) - else: - # if base field is small, cache monomials to minimise number of constructions - monomial_cache = {} - for row in range(m): - for col in range(m): - if not coeffs_map[row][col]: - basis_rows[row][col] = poly_ring.zero() - elif len(coeffs_map[row][col]) == 1: - # monomial case - deg, coeff = coeffs_map[row][col].popitem() + # construct polynomial matrix - potentially costly at polynomial construction stage for extension fields + coeffs_map = [[{} for _ in range(m)] for _ in range(m)] + # add identity part of basis + for i in range(m): + coeffs_map[i][i][degree_c[i]] = self.base_ring().one() + # add relation part of basis + for col in range(relation.ncols()): + for row in range(m): + coeffs_map[row][c[col]][d[col]] = coeffs_map[row][c[col]].get(d[col], self.base_ring().zero()) - relation[row,col] + # construct rows + basis_rows = [[None] * m for _ in range(m)] + if self.base_ring().degree() <= 1 or m*m*sigma < self.base_ring().order(): + # if base field is large, don't bother caching monomials + # or if polynomial construction is efficient (order = 1) + for row in range(m): + for col in range(m): + if not coeffs_map[row][col]: + basis_rows[row][col] = self.base_ring().zero() + else: + # construct a small polynomial (less costly), then shift + mindeg = min(coeffs_map[row][col].keys()) + maxdeg = max(coeffs_map[row][col].keys()) + entries = [self.base_ring().zero()] * (maxdeg-mindeg+1) + for deg, coeff in coeffs_map[row][col].items(): + entries[deg-mindeg] = coeff + basis_rows[row][col] = poly_ring(entries).shift(mindeg) + return matrix(poly_ring, basis_rows) + else: + # if base field is small, cache monomials to minimise number of constructions + monomial_cache = {} + for row in range(m): + for col in range(m): + if not coeffs_map[row][col]: + basis_rows[row][col] = poly_ring.zero() + elif len(coeffs_map[row][col]) == 1: + # monomial case + deg, coeff = coeffs_map[row][col].popitem() + if coeff not in monomial_cache: + monomial_cache[coeff] = poly_ring(coeff) + basis_rows[row][col] = monomial_cache[coeff].shift(deg) + else: + # general case + mindeg = min(coeffs_map[row][col].keys()) + poly = poly_ring.zero() + for deg, coeff in coeffs_map[row][col].items(): if coeff not in monomial_cache: monomial_cache[coeff] = poly_ring(coeff) - basis_rows[row][col] = monomial_cache[coeff].shift(deg) - else: - # general case - mindeg = min(coeffs_map[row][col].keys()) - # maxdeg = max(coeffs_map[row][col].keys()) - poly = poly_ring.zero() - for deg, coeff in coeffs_map[row][col].items(): - if coeff not in monomial_cache: - monomial_cache[coeff] = poly_ring(coeff) - poly += monomial_cache[coeff].shift(deg-mindeg) - basis_rows[row][col] = poly.shift(mindeg) - # entries = [poly_ring.base_ring().zero()] * (maxdeg-mindeg+1) - # for deg, coeff in coeffs_map[row][col].items(): - # entries[deg-mindeg] = coeff - # basis_rows[row][col] = poly_ring(entries).shift(mindeg) - # convert to matrix - output = matrix(poly_ring,basis_rows) - else: - # expected output is expanded matrix - row_width = m * (max(degree_c) + 1) - basis_rows = [[poly_ring.base_ring().zero()]*row_width for i in range(m)] - for i in range(m): - basis_rows[i][i+degree_c[i]*m] = poly_ring.base_ring().one() - for col in range(relation.ncols()): - for row in range(m): - basis_rows[row][c[col]+d[col]*m] -= relation[row,col] - output = matrix(poly_ring.base_ring(),basis_rows) - - return output + poly += monomial_cache[coeff].shift(deg-mindeg) + basis_rows[row][col] = poly.shift(mindeg) + # convert to matrix + return matrix(poly_ring,basis_rows) # a limited number of access-only properties are provided for matrices @property From 926420550c18596b4166a2abb474793d16570b9e Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Tue, 29 Jul 2025 15:55:10 +0200 Subject: [PATCH 21/28] update doc, add tests, remove whitespace --- src/sage/matrix/matrix2.pyx | 257 +++++++++++++++++++++++++++++------- 1 file changed, 211 insertions(+), 46 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index fc101b3000a..3a82627f472 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19285,13 +19285,20 @@ cdef class Matrix(Matrix1): OUTPUT: - A tuple (row_profile, pivot_matrix, column_profile): - * row_profile: list of the first r row indices in the striped Krylov matrix `K` corresponding to independent rows - * pivot: the submatrix of `K` given by those rows - * column_profile: list of the first r independent column indices in ``pivot_matrix`` + * ``row_profile``: list of the first r row indices in the striped + Krylov matrix ``K`` corresponding to + independent rows. If ``output_pairs`` is True, + then the output is the list of pairs (c_i, d_i) + where c_i is the corresponding row in the + original matrix and d_i the corresponding power + of ``J``. + * ``pivot_matrix``: the submatrix of ``K`` given by ``row_profile`` + * ``column_profile``: the first r independent column indices in + ``pivot_matrix`` where r is the rank of `K`. - + TESTS:: - + sage: R = GF(97) sage: E = matrix(R,[[27,49,29],[50,58,0],[77,10,29]]) sage: E @@ -19303,37 +19310,124 @@ cdef class Matrix(Matrix1): [0 1 0] [0 0 1] [0 0 0] - sage: E.krylov_rank_profile(J,3) + sage: degree = 3 + sage: striped_krylov_matrix = matrix.block([[E*J^i] for i in range(degree+1)],subdivide=False) + sage: striped_krylov_matrix + [27 49 29] + [50 58 0] + [77 10 29] + [ 0 27 49] + [ 0 50 58] + [ 0 77 10] + [ 0 0 27] + [ 0 0 50] + [ 0 0 77] + [ 0 0 0] + [ 0 0 0] + [ 0 0 0] + sage: E.krylov_rank_profile(J) ( [27 49 29] [50 58 0] (0, 1, 3), [ 0 27 49], (0, 1, 2) ) - sage: E.krylov_rank_profile(J,3,None,True) + sage: E.krylov_rank_profile(J,output_pairs=True) ( [27 49 29] [50 58 0] ((0, 0), (1, 0), (0, 1)), [ 0 27 49], (0, 1, 2) ) - sage: E.krylov_rank_profile(J,3,[0,3,6]) + + sage: rows = [(i+j*E.nrows()+1,i,j) for j in range(degree+1) for i in range(E.nrows())] + sage: rows + [(1, 0, 0), + (2, 1, 0), + (3, 2, 0), + (4, 0, 1), + (5, 1, 1), + (6, 2, 1), + (7, 0, 2), + (8, 1, 2), + (9, 2, 2), + (10, 0, 3), + (11, 1, 3), + (12, 2, 3)] + sage: shift = [0,3,6] + sage: rows.sort(key=lambda x: (x[2] + shift[x[1]],x[1])) + sage: rows + [(1, 0, 0), + (4, 0, 1), + (7, 0, 2), + (10, 0, 3), + (2, 1, 0), + (5, 1, 1), + (8, 1, 2), + (11, 1, 3), + (3, 2, 0), + (6, 2, 1), + (9, 2, 2), + (12, 2, 3)] + sage: striped_krylov_matrix.with_permuted_rows(Permutation([x[0] for x in rows])) + [27 49 29] + [ 0 27 49] + [ 0 0 27] + [ 0 0 0] + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [ 0 0 0] + [77 10 29] + [ 0 77 10] + [ 0 0 77] + [ 0 0 0] + sage: E.krylov_rank_profile(J,shift=shift) ( [27 49 29] [ 0 27 49] (0, 1, 2), [ 0 0 27], (0, 1, 2) ) - sage: E.krylov_rank_profile(J,3,[0,3,6],True) + sage: E.krylov_rank_profile(J,shift=shift,output_pairs=True) ( [27 49 29] [ 0 27 49] ((0, 0), (0, 1), (0, 2)), [ 0 0 27], (0, 1, 2) ) - sage: E.krylov_rank_profile(J,3,[3,0,2]) + + sage: shift = [3,0,2] + sage: rows.sort(key=lambda x: (x[2] + shift[x[1]],x[1])) + sage: rows + [(2, 1, 0), + (5, 1, 1), + (8, 1, 2), + (3, 2, 0), + (1, 0, 0), + (11, 1, 3), + (6, 2, 1), + (4, 0, 1), + (9, 2, 2), + (7, 0, 2), + (12, 2, 3), + (10, 0, 3)] + sage: striped_krylov_matrix.with_permuted_rows(Permutation([x[0] for x in rows])) + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [77 10 29] + [27 49 29] + [ 0 0 0] + [ 0 77 10] + [ 0 27 49] + [ 0 0 77] + [ 0 0 27] + [ 0 0 0] + [ 0 0 0] + sage: E.krylov_rank_profile(J,shift=shift) ( [50 58 0] [ 0 50 58] (0, 1, 2), [ 0 0 50], (0, 1, 2) ) - sage: E.krylov_rank_profile(J,3,[3,0,2],True) + sage: E.krylov_rank_profile(J,shift=shift,output_pairs=True) ( [50 58 0] [ 0 50 58] @@ -19370,64 +19464,64 @@ cdef class Matrix(Matrix1): if m == 0: return (),self,() - + # calculate row profile of self, with shift applied self_permutation = Permutation(sorted([i+1 for i in range(m)],key=lambda x:(shift[x-1],x-1))) - + self_permuted = self.with_permuted_rows(self_permutation) row_profile_self_permuted = self_permuted.row_rank_profile() row_profile_self = [(self_permutation(i+1)-1,0) for i in row_profile_self_permuted] - + exhausted = matrix.zero(self.base_ring(), 0, sigma) row_profile_exhausted = [] excluded_rows = set() - + r = len(row_profile_self) - + if r == 0: return (),self.matrix_from_rows([]),() - + M = self_permuted.matrix_from_rows(row_profile_self_permuted) - + J_L = None - + for l in range(math.ceil(math.log(degree,2)) + 1): L = pow(2,l) # adding 2^l to each degree row_extension = [(x[0],x[1] + L) for x in row_profile_self if x[1] + L <= degree] if len(row_extension) == 0: break - + # concatenate two sequences (order preserved) k = row_profile_exhausted + row_profile_self + row_extension - + # calculate sorting permutation, sort k by (shift[c]+d, c) k_perm = Permutation(sorted([i+1 for i in range(len(k))],key=lambda x: (shift[k[x-1][0]] + k[x-1][1],k[x-1][0]))) - + # fast calculation of rows formed by indices in k if J_L is None: J_L = J else: J_L = J_L * J_L - + M = matrix.block([[exhausted],[M],[M*J_L]],subdivide=False) - + # sort rows of M, find profile, translate to k (indices of full krylov matrix) M.permute_rows(k_perm) - + row_profile_M = M.row_rank_profile() r = len(row_profile_M) - + if r == sigma and l < math.ceil(math.log(degree,2)): tail = list(range(row_profile_M[-1]+1,M.nrows())) excluded_rows.update(set([k[k_perm(i+1)-1][0] for i in tail if i < len(k)])) - + xmi = [i for i in row_profile_M if k[k_perm(i+1)-1][0] in excluded_rows] imi = [i for i in row_profile_M if k[k_perm(i+1)-1][0] not in excluded_rows] - + row_profile_exhausted = [k[k_perm(i+1)-1] for i in xmi] row_profile_self = [k[k_perm(i+1)-1] for i in imi] - + exhausted = M.matrix_from_rows(xmi) M = M.matrix_from_rows(imi) else: @@ -19444,35 +19538,108 @@ cdef class Matrix(Matrix1): k = row_profile_exhausted + row_profile_self M = exhausted.stack(M) k_perm = Permutation(sorted([i+1 for i in range(len(k))],key=lambda x: (shift[k[x-1][0]] + k[x-1][1],k[x-1][0]))) - + M.permute_rows(k_perm) row_profile_self = [k[k_perm(i+1)-1] for i in range(len(k))] col_profile = M.col_rank_profile() - + if output_pairs: return tuple(row_profile_self), M, col_profile - + # convert c,d to actual position in striped Krylov matrix phi = lambda c,d : sum(min(max(shift[c] - shift[i] + d + (i < c and shift[i] <= shift[c] + degree),0),degree+1) for i in range(m)) row_profile_K = [phi(*row) for row in row_profile_self] - + return tuple(row_profile_K), M, col_profile def linear_interpolation_basis(self, J, var_name, degree=None, shift=None): r""" - Construct a linear interpolant basis for (``self``,`J`) in `s`-Popov form. + Construct a linear interpolant basis for (``self``,``J``) in ``s``-Popov form. INPUT: - ``J`` -- square matrix for Krylov construction. - ``var_name`` -- variable name for the returned polynomial matrix - ``degree`` -- upper bound on degree of minpoly(`J`), power of - - ``shift`` -- list of integers (optional): controls priority of rows. + - ``shift`` -- list of self.nrows() integers (optional): controls + priority of rows. OUTPUT: - - Matrix whose rows form an interpolation basis. + - A self.ncols() * self.ncols() matrix whose rows form an interpolation + basis. + + TESTS:: + + sage: R = GF(97) + sage: E = matrix(R,[[27,49,29],[50,58,0],[77,10,29]]) + sage: E + [27 49 29] + [50 58 0] + [77 10 29] + sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) + sage: J + [0 1 0] + [0 0 1] + [0 0 0] + sage: degree = E.ncols() + sage: striped_krylov_matrix = matrix.block([[E*J^i] for i in range(degree+1)],subdivide=False) + sage: striped_krylov_matrix + [27 49 29] + [50 58 0] + [77 10 29] + [ 0 27 49] + [ 0 50 58] + [ 0 77 10] + [ 0 0 27] + [ 0 0 50] + [ 0 0 77] + [ 0 0 0] + [ 0 0 0] + [ 0 0 0] + sage: basis = E.linear_interpolation_basis(J,'x') + sage: basis + [x^2 + 40*x + 82 76 0] + [ 3*x + 13 x + 57 0] + [ 96 96 1] + sage: basis.is_popov() + True + sage: basis.degree() <= E.ncols() + True + sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) + True + sage: len(E.krylov_rank_profile(J)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + True + + sage: shift = [0,3,6] + sage: basis = E.linear_interpolation_basis(J,'x',shift=shift) + sage: basis + [ x^3 0 0] + [60*x^2 + 72*x + 70 1 0] + [60*x^2 + 72*x + 69 0 1] + sage: basis.is_popov(shifts=shift) + True + sage: basis.degree() <= E.ncols() + True + sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) + True + sage: len(E.krylov_rank_profile(J,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + True + sage: shift = [3,0,2] + sage: basis = E.linear_interpolation_basis(J,'x',shift=shift) + sage: basis + [ 1 26*x^2 + 49*x + 79 0] + [ 0 x^3 0] + [ 0 26*x^2 + 49*x + 78 1] + sage: basis.is_popov(shifts=shift) + True + sage: basis.degree() <= E.ncols() + True + sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) + True + sage: len(E.krylov_rank_profile(J,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + True """ from sage.combinat.permutation import Permutation from sage.matrix.constructor import matrix @@ -19505,19 +19672,19 @@ cdef class Matrix(Matrix1): m = self.nrows() sigma = self.ncols() - + poly_ring = PolynomialRing(self.base_ring(),var_name) - + # calculate krylov profile row_profile, pivot, col_profile = self.krylov_rank_profile(J,degree,shift,output_pairs=True) if len(row_profile) == 0: return matrix.identity(poly_ring,m) - + # (c_k, d_k) = phi^-1 (row_i) #c, d = zip(*(phi_inv(i) for i in row_profile)) c, d = zip(*(row for row in row_profile)) - + # degree_c = 0 or max d[k] such that c[k] = c inv_c = [None]*m degree_c = [0]*m @@ -19525,7 +19692,7 @@ cdef class Matrix(Matrix1): if d[k] >= degree_c[c[k]]: degree_c[c[k]] = d[k] + 1 inv_c[c[k]] = k - + T_absent_indices = [i for i in range(m) if inv_c[i] is None] T_present_indices = [inv_c[i] for i in range(m) if inv_c[i] is not None] D_absent = self.matrix_from_rows_and_columns(T_absent_indices, col_profile) @@ -19540,13 +19707,11 @@ cdef class Matrix(Matrix1): else: D_rows.append(D_present[idx_p]) idx_p += 1 - + C = pivot.matrix_from_columns(col_profile) - D = matrix(D_rows) - relation = D*C.inverse() - + # construct polynomial matrix - potentially costly at polynomial construction stage for extension fields coeffs_map = [[{} for _ in range(m)] for _ in range(m)] # add identity part of basis @@ -19598,7 +19763,7 @@ cdef class Matrix(Matrix1): basis_rows[row][col] = poly.shift(mindeg) # convert to matrix return matrix(poly_ring,basis_rows) - + # a limited number of access-only properties are provided for matrices @property def T(self): From fa0bf8c31f7b59ba4e3fa04ecd15d9084d9c0c74 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Tue, 29 Jul 2025 18:07:45 +0200 Subject: [PATCH 22/28] remove testing functions from development --- src/sage/matrix/matrix2.pyx | 233 -------------------- src/sage/matrix/matrix_polynomial_dense.pyx | 41 ---- 2 files changed, 274 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 3a82627f472..3d239472321 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19040,237 +19040,6 @@ cdef class Matrix(Matrix1): U.rescale_col(n - 1, -1) return U - def striped_krylov_matrix(self, J, degree, shift=None): - r""" - Return the block matrix of the following form, with rows permuted according to shift. - The following uses `E` to refer to ``self``, and `d` to refer to ``degree``. - - [ ] - [ E ] - [ ] - [-----] - [ ] - [ EJ ] - [ ] - [-----] - [ . ] - [ . ] - [ . ] - [-----] - [ ] - [ EJ^d] - [ ] - - INPUT: - - - ``J`` -- a square matrix of size equal to the number of columns of ``self``. - - ``degree`` -- integer, the maximum exponent for the Krylov matrix. - - ``shift`` -- list of integers (optional), row priority shifts. If ``None``, - defaults to all zero. - - OUTPUT: - - - A matrix with block rows [E, EJ, EJ^2, ..., EJ^d], row-permuted by shift. - - EXAMPLES: - - sage: R. = GF(97)[] - sage: E = matrix([[27,49,29],[50,58,0],[77,10,29]]) - sage: E - [27 49 29] - [50 58 0] - [77 10 29] - sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) - sage: J - [0 1 0] - [0 0 1] - [0 0 0] - sage: E.striped_krylov_matrix(J,3) - [27 49 29] - [50 58 0] - [77 10 29] - [ 0 27 49] - [ 0 50 58] - [ 0 77 10] - [ 0 0 27] - [ 0 0 50] - [ 0 0 77] - [ 0 0 0] - [ 0 0 0] - [ 0 0 0] - sage: E.striped_krylov_matrix(J,3,[0,3,6]) - [27 49 29] - [ 0 27 49] - [ 0 0 27] - [ 0 0 0] - [50 58 0] - [ 0 50 58] - [ 0 0 50] - [ 0 0 0] - [77 10 29] - [ 0 77 10] - [ 0 0 77] - [ 0 0 0] - sage: E.striped_krylov_matrix(J,3,[3,0,2]) - [50 58 0] - [ 0 50 58] - [ 0 0 50] - [77 10 29] - [27 49 29] - [ 0 0 0] - [ 0 77 10] - [ 0 27 49] - [ 0 0 77] - [ 0 0 27] - [ 0 0 0] - [ 0 0 0] - - """ - from sage.combinat.permutation import Permutation - from sage.matrix.constructor import matrix # for matrix.block - - m = self.nrows() - - # if no shift, this is equivalent to 0s - if shift is None: - shift = [0]*self.nrows() - - # size checks - if len(shift) != m: - raise ValueError(f"The shift should have the same number of elements as the rows of the matrix ({m}), but had {len(shift)}.") - if not J.is_square() or J.nrows() != self.ncols(): - raise ValueError(f"The matrix J should be a square matrix and match the number of columns of self ({self.ncols()}), but is of dimension {J.nrows()} x {J.ncols()}.") - - # define priority and indexing functions - # index: [0,1,...,m-1,m,...,2m-1,...,m*degree-1] -> [(0,0),(1,0),...,(m-1,0),(0,1),...,(m-1,1),...,(m-1,degree)] - priority = lambda c,d : shift[c] + d - index = lambda i : (i%m,i//m) - - # deduce permutation by sorting priorities - # rows should be ordered by ascending priority, then by associated row number in E (i%m) - priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) - priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) - - # build all blocks for the striped matrix - blocks = [] - - # for i from 0 to degree (inclusive), add E*J^i - J_i = matrix.identity(self.ncols()) - for i in range(degree+1): - blocks.append([self*J_i]) - if i < degree: - J_i *= J - - # return block matrix permuted according to shift - krylov = matrix.block(blocks,subdivide=False) - krylov.permute_rows(priority_permutation) - return krylov - - def naive_krylov_rank_profile(self, J, degree, shift=None): - r""" - Compute the rank profile (row and column) of the striped Krylov matrix - built from ``self`` and matrix `J`. - - INPUT: - - - ``J`` -- square matrix used in the Krylov construction. - - ``degree`` -- maximum power of J in Krylov matrix. - - ``shift`` -- list of integers (optional): priority row shift. - - OUTPUT: - - - A tuple (row_profile, pivot_matrix, column_profile): - * row_profile: list of the first r row indices in the striped Krylov matrix `K` corresponding to independent rows - * pivot: the submatrix of `K` given by those rows - * column_profile: list of the first r independent column indices in ``pivot_matrix`` - where r is the rank of `K`. - - EXAMPLES: - - sage: R. = GF(97)[] - sage: E = matrix([[27,49,29],[50,58,0],[77,10,29]]) - sage: E - [27 49 29] - [50 58 0] - [77 10 29] - sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) - sage: J - [0 1 0] - [0 0 1] - [0 0 0] - sage: E.striped_krylov_matrix(J,3) - [27 49 29] - [50 58 0] - [77 10 29] - [ 0 27 49] - [ 0 50 58] - [ 0 77 10] - [ 0 0 27] - [ 0 0 50] - [ 0 0 77] - [ 0 0 0] - [ 0 0 0] - [ 0 0 0] - sage: E.naive_krylov_rank_profile(J,3) - ( - [27 49 29] - [50 58 0] - (0, 1, 3), [ 0 27 49], (0, 1, 2) - ) - sage: E.striped_krylov_matrix(J,3,[0,3,6]) - [27 49 29] - [ 0 27 49] - [ 0 0 27] - [ 0 0 0] - [50 58 0] - [ 0 50 58] - [ 0 0 50] - [ 0 0 0] - [77 10 29] - [ 0 77 10] - [ 0 0 77] - [ 0 0 0] - sage: E.naive_krylov_rank_profile(J,3,[0,3,6]) - ( - [27 49 29] - [ 0 27 49] - (0, 1, 2), [ 0 0 27], (0, 1, 2) - ) - sage: E.striped_krylov_matrix(J,3,[3,0,2]) - [50 58 0] - [ 0 50 58] - [ 0 0 50] - [77 10 29] - [27 49 29] - [ 0 0 0] - [ 0 77 10] - [ 0 27 49] - [ 0 0 77] - [ 0 0 27] - [ 0 0 0] - [ 0 0 0] - sage: E.naive_krylov_rank_profile(J,3,[3,0,2]) - ( - [50 58 0] - [ 0 50 58] - (0, 1, 2), [ 0 0 50], (0, 1, 2) - ) - """ - from sage.matrix.constructor import matrix - - K = self.striped_krylov_matrix(J,degree,shift) - - # calculate first independent rows of K - row_profile = K.row_rank_profile() - - # construct submatrix - pivot = K.matrix_from_rows(row_profile) - - # calculate first independent columns of pivot - col_profile = pivot.col_rank_profile() - - return row_profile,pivot,col_profile - def krylov_rank_profile(self, J, degree=None, shift=None, output_pairs=False): r""" Compute the rank profile (row and column) of the striped Krylov matrix @@ -20387,8 +20156,6 @@ def _matrix_power_symbolic(A, n): return P * M * Pinv - - class NotFullRankError(ValueError): """ diff --git a/src/sage/matrix/matrix_polynomial_dense.pyx b/src/sage/matrix/matrix_polynomial_dense.pyx index 4e147c2f82e..db6ff8f4ccf 100644 --- a/src/sage/matrix/matrix_polynomial_dense.pyx +++ b/src/sage/matrix/matrix_polynomial_dense.pyx @@ -5337,44 +5337,3 @@ cdef class Matrix_polynomial_dense(Matrix_generic_dense): return False return True - - def expansion(self, degree=None): - r""" - Return the horizontal expansion of the matrix, - where the ith block corresponds to the coefficients for x^i. - - INPUT: - - - ``degree`` -- integer, the degree bound on the matrix; - if ``None``, defaults to degree of matrix. - - OUTPUT: matrix - - EXAMPLES: - sage: R. = GF(97)[] - sage: M = matrix([[x^2+36*x,31*x,0],[3*x+13,x+57,0],[96,96,1]]) - sage: M - [x^2 + 36*x 31*x 0] - [ 3*x + 13 x + 57 0] - [ 96 96 1] - sage: M.expansion() - [ 0 0 0 36 31 0 1 0 0] - [13 57 0 3 1 0 0 0 0] - [96 96 1 0 0 0 0 0 0] - sage: M.expansion(1) - [ 0 0 0 36 31 0] - [13 57 0 3 1 0] - [96 96 1 0 0 0] - sage: M.expansion(3) - [ 0 0 0 36 31 0 1 0 0 0 0 0] - [13 57 0 3 1 0 0 0 0 0 0 0] - [96 96 1 0 0 0 0 0 0 0 0 0] - """ - from sage.matrix.constructor import matrix # for matrix.block - - if degree is None: - degree = self.degree() - - # compute each coefficient matrix - coeff_matrices = [[self.coefficient_matrix(i) for i in range(degree+1)]] - return matrix.block(coeff_matrices,subdivide=False) \ No newline at end of file From e55bb367963452232f0dc80f2ffaa7a7b0b7ab14 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Thu, 31 Jul 2025 18:43:45 +0200 Subject: [PATCH 23/28] updates --- src/sage/matrix/matrix2.pyx | 442 +++++++++--------- .../matrix/matrix_modn_dense_template.pxi | 8 +- 2 files changed, 221 insertions(+), 229 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 3d239472321..adbd9117485 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -1065,22 +1065,6 @@ cdef class Matrix(Matrix1): v = self.transpose().pivots() self.cache('pivot_rows', v) return v - - def row_rank_profile(self): - """ - Alias for pivot_rows(self). Return the pivot row positions for this matrix, which are a topmost - subset of the rows that span the row space and are linearly - independent. - """ - return self.pivot_rows() - - def col_rank_profile(self): - """ - Alias for pivots(self). Return the pivot column positions for this matrix, which are a leftmost - subset of the columns that span the column space and are linearly - independent. - """ - return self.pivots() def _solve_right_general(self, B, check=True): r""" @@ -19053,155 +19037,155 @@ cdef class Matrix(Matrix1): OUTPUT: - - A tuple (row_profile, pivot_matrix, column_profile): - * ``row_profile``: list of the first r row indices in the striped + A tuple (row_profile, pivot_matrix, column_profile): + - ``row_profile``: list of the first r row indices in the striped Krylov matrix ``K`` corresponding to independent rows. If ``output_pairs`` is True, then the output is the list of pairs (c_i, d_i) where c_i is the corresponding row in the original matrix and d_i the corresponding power of ``J``. - * ``pivot_matrix``: the submatrix of ``K`` given by ``row_profile`` - * ``column_profile``: the first r independent column indices in + - ``pivot_matrix``: the submatrix of ``K`` given by ``row_profile`` + - ``column_profile``: the first r independent column indices in ``pivot_matrix`` where r is the rank of `K`. TESTS:: - sage: R = GF(97) - sage: E = matrix(R,[[27,49,29],[50,58,0],[77,10,29]]) - sage: E - [27 49 29] - [50 58 0] - [77 10 29] - sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) - sage: J - [0 1 0] - [0 0 1] - [0 0 0] - sage: degree = 3 - sage: striped_krylov_matrix = matrix.block([[E*J^i] for i in range(degree+1)],subdivide=False) - sage: striped_krylov_matrix - [27 49 29] - [50 58 0] - [77 10 29] - [ 0 27 49] - [ 0 50 58] - [ 0 77 10] - [ 0 0 27] - [ 0 0 50] - [ 0 0 77] - [ 0 0 0] - [ 0 0 0] - [ 0 0 0] - sage: E.krylov_rank_profile(J) - ( - [27 49 29] - [50 58 0] - (0, 1, 3), [ 0 27 49], (0, 1, 2) - ) - sage: E.krylov_rank_profile(J,output_pairs=True) - ( - [27 49 29] - [50 58 0] - ((0, 0), (1, 0), (0, 1)), [ 0 27 49], (0, 1, 2) - ) + sage: R = GF(97) + sage: E = matrix(R,[[27,49,29],[50,58,0],[77,10,29]]) + sage: E + [27 49 29] + [50 58 0] + [77 10 29] + sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) + sage: J + [0 1 0] + [0 0 1] + [0 0 0] + sage: degree = 3 + sage: striped_krylov_matrix = matrix.block([[E*J^i] for i in range(degree+1)],subdivide=False) + sage: striped_krylov_matrix + [27 49 29] + [50 58 0] + [77 10 29] + [ 0 27 49] + [ 0 50 58] + [ 0 77 10] + [ 0 0 27] + [ 0 0 50] + [ 0 0 77] + [ 0 0 0] + [ 0 0 0] + [ 0 0 0] + sage: E.krylov_rank_profile(J) + ( + [27 49 29] + [50 58 0] + (0, 1, 3), [ 0 27 49], (0, 1, 2) + ) + sage: E.krylov_rank_profile(J,output_pairs=True) + ( + [27 49 29] + [50 58 0] + ((0, 0), (1, 0), (0, 1)), [ 0 27 49], (0, 1, 2) + ) - sage: rows = [(i+j*E.nrows()+1,i,j) for j in range(degree+1) for i in range(E.nrows())] - sage: rows - [(1, 0, 0), - (2, 1, 0), - (3, 2, 0), - (4, 0, 1), - (5, 1, 1), - (6, 2, 1), - (7, 0, 2), - (8, 1, 2), - (9, 2, 2), - (10, 0, 3), - (11, 1, 3), - (12, 2, 3)] - sage: shift = [0,3,6] - sage: rows.sort(key=lambda x: (x[2] + shift[x[1]],x[1])) - sage: rows - [(1, 0, 0), - (4, 0, 1), - (7, 0, 2), - (10, 0, 3), - (2, 1, 0), - (5, 1, 1), - (8, 1, 2), - (11, 1, 3), - (3, 2, 0), - (6, 2, 1), - (9, 2, 2), - (12, 2, 3)] - sage: striped_krylov_matrix.with_permuted_rows(Permutation([x[0] for x in rows])) - [27 49 29] - [ 0 27 49] - [ 0 0 27] - [ 0 0 0] - [50 58 0] - [ 0 50 58] - [ 0 0 50] - [ 0 0 0] - [77 10 29] - [ 0 77 10] - [ 0 0 77] - [ 0 0 0] - sage: E.krylov_rank_profile(J,shift=shift) - ( - [27 49 29] - [ 0 27 49] - (0, 1, 2), [ 0 0 27], (0, 1, 2) - ) - sage: E.krylov_rank_profile(J,shift=shift,output_pairs=True) - ( - [27 49 29] - [ 0 27 49] - ((0, 0), (0, 1), (0, 2)), [ 0 0 27], (0, 1, 2) - ) + sage: rows = [(i+j*E.nrows()+1,i,j) for j in range(degree+1) for i in range(E.nrows())] + sage: rows + [(1, 0, 0), + (2, 1, 0), + (3, 2, 0), + (4, 0, 1), + (5, 1, 1), + (6, 2, 1), + (7, 0, 2), + (8, 1, 2), + (9, 2, 2), + (10, 0, 3), + (11, 1, 3), + (12, 2, 3)] + sage: shift = [0,3,6] + sage: rows.sort(key=lambda x: (x[2] + shift[x[1]],x[1])) + sage: rows + [(1, 0, 0), + (4, 0, 1), + (7, 0, 2), + (10, 0, 3), + (2, 1, 0), + (5, 1, 1), + (8, 1, 2), + (11, 1, 3), + (3, 2, 0), + (6, 2, 1), + (9, 2, 2), + (12, 2, 3)] + sage: striped_krylov_matrix.with_permuted_rows(Permutation([x[0] for x in rows])) + [27 49 29] + [ 0 27 49] + [ 0 0 27] + [ 0 0 0] + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [ 0 0 0] + [77 10 29] + [ 0 77 10] + [ 0 0 77] + [ 0 0 0] + sage: E.krylov_rank_profile(J,shift=shift) + ( + [27 49 29] + [ 0 27 49] + (0, 1, 2), [ 0 0 27], (0, 1, 2) + ) + sage: E.krylov_rank_profile(J,shift=shift,output_pairs=True) + ( + [27 49 29] + [ 0 27 49] + ((0, 0), (0, 1), (0, 2)), [ 0 0 27], (0, 1, 2) + ) - sage: shift = [3,0,2] - sage: rows.sort(key=lambda x: (x[2] + shift[x[1]],x[1])) - sage: rows - [(2, 1, 0), - (5, 1, 1), - (8, 1, 2), - (3, 2, 0), - (1, 0, 0), - (11, 1, 3), - (6, 2, 1), - (4, 0, 1), - (9, 2, 2), - (7, 0, 2), - (12, 2, 3), - (10, 0, 3)] - sage: striped_krylov_matrix.with_permuted_rows(Permutation([x[0] for x in rows])) - [50 58 0] - [ 0 50 58] - [ 0 0 50] - [77 10 29] - [27 49 29] - [ 0 0 0] - [ 0 77 10] - [ 0 27 49] - [ 0 0 77] - [ 0 0 27] - [ 0 0 0] - [ 0 0 0] - sage: E.krylov_rank_profile(J,shift=shift) - ( - [50 58 0] - [ 0 50 58] - (0, 1, 2), [ 0 0 50], (0, 1, 2) - ) - sage: E.krylov_rank_profile(J,shift=shift,output_pairs=True) - ( - [50 58 0] - [ 0 50 58] - ((1, 0), (1, 1), (1, 2)), [ 0 0 50], (0, 1, 2) - ) + sage: shift = [3,0,2] + sage: rows.sort(key=lambda x: (x[2] + shift[x[1]],x[1])) + sage: rows + [(2, 1, 0), + (5, 1, 1), + (8, 1, 2), + (3, 2, 0), + (1, 0, 0), + (11, 1, 3), + (6, 2, 1), + (4, 0, 1), + (9, 2, 2), + (7, 0, 2), + (12, 2, 3), + (10, 0, 3)] + sage: striped_krylov_matrix.with_permuted_rows(Permutation([x[0] for x in rows])) + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [77 10 29] + [27 49 29] + [ 0 0 0] + [ 0 77 10] + [ 0 27 49] + [ 0 0 77] + [ 0 0 27] + [ 0 0 0] + [ 0 0 0] + sage: E.krylov_rank_profile(J,shift=shift) + ( + [50 58 0] + [ 0 50 58] + (0, 1, 2), [ 0 0 50], (0, 1, 2) + ) + sage: E.krylov_rank_profile(J,shift=shift,output_pairs=True) + ( + [50 58 0] + [ 0 50 58] + ((1, 0), (1, 1), (1, 2)), [ 0 0 50], (0, 1, 2) + ) """ from sage.matrix.constructor import matrix import math @@ -19323,7 +19307,19 @@ cdef class Matrix(Matrix1): def linear_interpolation_basis(self, J, var_name, degree=None, shift=None): r""" - Construct a linear interpolant basis for (``self``,``J``) in ``s``-Popov form. + Return a shifted minimal linear interpolation basis for + (``self``,``J``) in ``s``-Popov form with respect to a shift ``s``. + + Given a K[x]-module V defined by multiplication matrix J, i.e. elements + are vectors of length n in K, scalars are polynomials in K[x], and + scalar multiplication is defined by p v := v p(J), ``self`` represents + ``self.nrows()`` elements e_1, ..., e_m in V: + + A linear interpolation basis is a set of tuples of length + ``self.nrows()`` in K[x] that form a basis for solutions to the + equation ``p_1 e_1 + ... + p_n e_n = 0``. + + See https://arxiv.org/abs/1512.03503 INPUT: @@ -19340,75 +19336,75 @@ cdef class Matrix(Matrix1): TESTS:: - sage: R = GF(97) - sage: E = matrix(R,[[27,49,29],[50,58,0],[77,10,29]]) - sage: E - [27 49 29] - [50 58 0] - [77 10 29] - sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) - sage: J - [0 1 0] - [0 0 1] - [0 0 0] - sage: degree = E.ncols() - sage: striped_krylov_matrix = matrix.block([[E*J^i] for i in range(degree+1)],subdivide=False) - sage: striped_krylov_matrix - [27 49 29] - [50 58 0] - [77 10 29] - [ 0 27 49] - [ 0 50 58] - [ 0 77 10] - [ 0 0 27] - [ 0 0 50] - [ 0 0 77] - [ 0 0 0] - [ 0 0 0] - [ 0 0 0] - sage: basis = E.linear_interpolation_basis(J,'x') - sage: basis - [x^2 + 40*x + 82 76 0] - [ 3*x + 13 x + 57 0] - [ 96 96 1] - sage: basis.is_popov() - True - sage: basis.degree() <= E.ncols() - True - sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) - True - sage: len(E.krylov_rank_profile(J)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) - True - - sage: shift = [0,3,6] - sage: basis = E.linear_interpolation_basis(J,'x',shift=shift) - sage: basis - [ x^3 0 0] - [60*x^2 + 72*x + 70 1 0] - [60*x^2 + 72*x + 69 0 1] - sage: basis.is_popov(shifts=shift) - True - sage: basis.degree() <= E.ncols() - True - sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) - True - sage: len(E.krylov_rank_profile(J,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) - True - - sage: shift = [3,0,2] - sage: basis = E.linear_interpolation_basis(J,'x',shift=shift) - sage: basis - [ 1 26*x^2 + 49*x + 79 0] - [ 0 x^3 0] - [ 0 26*x^2 + 49*x + 78 1] - sage: basis.is_popov(shifts=shift) - True - sage: basis.degree() <= E.ncols() - True - sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) - True - sage: len(E.krylov_rank_profile(J,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) - True + sage: R = GF(97) + sage: E = matrix(R,[[27,49,29],[50,58,0],[77,10,29]]) + sage: E + [27 49 29] + [50 58 0] + [77 10 29] + sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) + sage: J + [0 1 0] + [0 0 1] + [0 0 0] + sage: degree = E.ncols() + sage: striped_krylov_matrix = matrix.block([[E*J^i] for i in range(degree+1)],subdivide=False) + sage: striped_krylov_matrix + [27 49 29] + [50 58 0] + [77 10 29] + [ 0 27 49] + [ 0 50 58] + [ 0 77 10] + [ 0 0 27] + [ 0 0 50] + [ 0 0 77] + [ 0 0 0] + [ 0 0 0] + [ 0 0 0] + sage: basis = E.linear_interpolation_basis(J,'x') + sage: basis + [x^2 + 40*x + 82 76 0] + [ 3*x + 13 x + 57 0] + [ 96 96 1] + sage: basis.is_popov() + True + sage: basis.degree() <= E.ncols() + True + sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) + True + sage: len(E.krylov_rank_profile(J)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + True + + sage: shift = [0,3,6] + sage: basis = E.linear_interpolation_basis(J,'x',shift=shift) + sage: basis + [ x^3 0 0] + [60*x^2 + 72*x + 70 1 0] + [60*x^2 + 72*x + 69 0 1] + sage: basis.is_popov(shifts=shift) + True + sage: basis.degree() <= E.ncols() + True + sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) + True + sage: len(E.krylov_rank_profile(J,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + True + + sage: shift = [3,0,2] + sage: basis = E.linear_interpolation_basis(J,'x',shift=shift) + sage: basis + [ 1 26*x^2 + 49*x + 79 0] + [ 0 x^3 0] + [ 0 26*x^2 + 49*x + 78 1] + sage: basis.is_popov(shifts=shift) + True + sage: basis.degree() <= E.ncols() + True + sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) + True + sage: len(E.krylov_rank_profile(J,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + True """ from sage.combinat.permutation import Permutation from sage.matrix.constructor import matrix @@ -19450,8 +19446,6 @@ cdef class Matrix(Matrix1): if len(row_profile) == 0: return matrix.identity(poly_ring,m) - # (c_k, d_k) = phi^-1 (row_i) - #c, d = zip(*(phi_inv(i) for i in row_profile)) c, d = zip(*(row for row in row_profile)) # degree_c = 0 or max d[k] such that c[k] = c diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi index 4655b3be180..2a304da51d8 100644 --- a/src/sage/matrix/matrix_modn_dense_template.pxi +++ b/src/sage/matrix/matrix_modn_dense_template.pxi @@ -204,13 +204,12 @@ cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nr applyP(F[0], FflasRight, FflasNoTrans, nrows, 0, r, entries, ncols, Q) - cdef list column_pivots = [int(Q[i]) for i in range(r)] - cdef list row_pivots = sorted([int(P[i]) for i in range(r)]) + cdef list pivots = [int(Q[i]) for i in range(r)] sig_free(P) sig_free(Q) del F - return r, column_pivots, row_pivots + return r, pivots cdef inline linbox_echelonize_efd(celement modulus, celement* entries, Py_ssize_t nrows, Py_ssize_t ncols): # See trac #13878: This is to avoid sending invalid data to linbox, @@ -1789,9 +1788,8 @@ cdef class Matrix_modn_dense_template(Matrix_dense): r, pivots = linbox_echelonize_efd(self.p, self._entries, self._nrows, self._ncols) else: - r, pivots, row_pivots = linbox_echelonize(self.p, self._entries, + r, pivots = linbox_echelonize(self.p, self._entries, self._nrows, self._ncols) - self.cache('pivot_rows',tuple(row_pivots)) verbose('done with echelonize', t) self.cache('in_echelon_form', True) self.cache('rank', r) From 4f4d6e811601d8209b4e5b96c2a8d31cb5ccf973 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Thu, 31 Jul 2025 22:31:14 +0200 Subject: [PATCH 24/28] changes --- src/sage/matrix/matrix2.pyx | 6 +++--- src/sage/matrix/matrix_modn_dense_template.pxi | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 2e580141736..53aaa80572c 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19301,7 +19301,7 @@ cdef class Matrix(Matrix1): self_permutation = Permutation(sorted([i+1 for i in range(m)],key=lambda x:(shift[x-1],x-1))) self_permuted = self.with_permuted_rows(self_permutation) - row_profile_self_permuted = self_permuted.row_rank_profile() + row_profile_self_permuted = self_permuted.pivot_rows() row_profile_self = [(self_permutation(i+1)-1,0) for i in row_profile_self_permuted] exhausted = matrix.zero(self.base_ring(), 0, sigma) @@ -19341,7 +19341,7 @@ cdef class Matrix(Matrix1): # sort rows of M, find profile, translate to k (indices of full krylov matrix) M.permute_rows(k_perm) - row_profile_M = M.row_rank_profile() + row_profile_M = M.pivot_rows() r = len(row_profile_M) if r == sigma and l < math.ceil(math.log(degree,2)): @@ -19373,7 +19373,7 @@ cdef class Matrix(Matrix1): M.permute_rows(k_perm) row_profile_self = [k[k_perm(i+1)-1] for i in range(len(k))] - col_profile = M.col_rank_profile() + col_profile = M.pivots() if output_pairs: return tuple(row_profile_self), M, col_profile diff --git a/src/sage/matrix/matrix_modn_dense_template.pxi b/src/sage/matrix/matrix_modn_dense_template.pxi index 2a304da51d8..8e9e7bf0f89 100644 --- a/src/sage/matrix/matrix_modn_dense_template.pxi +++ b/src/sage/matrix/matrix_modn_dense_template.pxi @@ -176,7 +176,7 @@ cdef inline linbox_echelonize(celement modulus, celement* entries, Py_ssize_t nr Return the reduced row echelon form of this matrix. """ if linbox_is_zero(modulus, entries, nrows, ncols): - return 0, [], [] + return 0, [] cdef Py_ssize_t i, j cdef ModField *F = new ModField(modulus) From 1a89eee290791f56c9fc27ba36f10df821cce87e Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Mon, 4 Aug 2025 16:14:55 +0200 Subject: [PATCH 25/28] readd naive methods, changes from review --- src/sage/matrix/matrix2.pyx | 517 ++++++++++++++++++++++++++++++------ 1 file changed, 432 insertions(+), 85 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 53aaa80572c..d59eefea214 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19103,31 +19103,378 @@ cdef class Matrix(Matrix1): U.rescale_col(n - 1, -1) return U - def krylov_rank_profile(self, J, degree=None, shift=None, output_pairs=False): + def striped_krylov_matrix(self, M, shift=None, degree=None): + r""" + Return the block matrix of the following form, with rows permuted according to shift. + The following uses `E` to refer to ``self``, and `d` to refer to ``degree``. + + [ ] + [ E ] + [ ] + [-----] + [ ] + [ EM ] + [ ] + [-----] + [ . ] + [ . ] + [ . ] + [-----] + [ ] + [ EM^d] + [ ] + + INPUT: + + - ``M`` -- a square matrix of size equal to the number of columns of ``self``. + - ``degree`` -- integer, the maximum exponent for the Krylov matrix. + - ``shift`` -- list of integers (optional), row priority shifts. If ``None``, + defaults to all zero. + + OUTPUT: + + - A matrix with block rows [E, EM, EM^2, ..., EM^d], row-permuted by shift. + + EXAMPLES:: + + sage: R. = GF(97) + sage: E = matrix(R,[[27,49,29],[50,58,0],[77,10,29]]) + sage: E + [27 49 29] + [50 58 0] + [77 10 29] + sage: M = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) + sage: M + [0 1 0] + [0 0 1] + [0 0 0] + sage: E.striped_krylov_matrix(M) + [27 49 29] + [50 58 0] + [77 10 29] + [ 0 27 49] + [ 0 50 58] + [ 0 77 10] + [ 0 0 27] + [ 0 0 50] + [ 0 0 77] + [ 0 0 0] + [ 0 0 0] + [ 0 0 0] + sage: E.striped_krylov_matrix(M,[0,3,6]) + [27 49 29] + [ 0 27 49] + [ 0 0 27] + [ 0 0 0] + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [ 0 0 0] + [77 10 29] + [ 0 77 10] + [ 0 0 77] + [ 0 0 0] + sage: E.striped_krylov_matrix(M,[3,0,2]) + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [77 10 29] + [27 49 29] + [ 0 0 0] + [ 0 77 10] + [ 0 27 49] + [ 0 0 77] + [ 0 0 27] + [ 0 0 0] + [ 0 0 0] + + """ + from sage.combinat.permutation import Permutation + import math + from sage.matrix.constructor import matrix # for matrix.block + + if not isinstance(M, Matrix): + raise TypeError("striped_krylov_matrix: M is not a matrix") + if M.nrows() != self.ncols() or M.ncols() != self.ncols(): + raise ValueError("striped_krylov_matrix: matrix M does not have correct dimensions.") + if M.base_ring() != self.base_ring(): + raise ValueError("striped_krylov_matrix: matrix M does not have same base ring as E.") + + if degree is None: + degree = self.ncols() + if not isinstance(degree, (int, sage.rings.integer.Integer)): + raise TypeError("striped_krylov_matrix: degree is not an integer.") + if degree < 0: + raise ValueError("striped_krylov_matrix: degree bound cannot be negative.") + + if shift is None or shift == 0: + shift = (ZZ**self.nrows()).zero() + if isinstance(shift, (list, tuple)): + shift = (ZZ**self.nrows())(shift) + if shift not in ZZ**self.nrows(): + raise ValueError(f"striped_krylov_matrix: shift is not an integer vector of length {self.nrows()}.") + + m = self.nrows() + + # define priority and indexing functions + # index: [0,1,...,m-1,m,...,2m-1,...,m*degree-1] -> [(0,0),(1,0),...,(m-1,0),(0,1),...,(m-1,1),...,(m-1,degree)] + priority = lambda c,d : shift[c] + d + index = lambda i : (i%m,i//m) + + # deduce permutation by sorting priorities + # rows should be ordered by ascending priority, then by associated row number in E (i%m) + priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) + priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) + + # build all blocks for the striped matrix + blocks = [self] + + # for i from 0 to degree (inclusive), add E*M^i + M_i = None + for i in range(math.ceil(math.log(degree + 1, 2))): + if M_i is None: + M_i = M + else: + M_i = M_i * M_i + blocks = [matrix.block([[blocks[0]],[blocks[1]]],subdivide=False)] + blocks.append(blocks[0] * M_i) + + if (degree + 1) & degree != 0: + blocks[1] = blocks[1].matrix_from_rows(range((degree + 1) * m - blocks[1].nrows())) + + # return block matrix permuted according to shift + krylov = matrix.block([[blocks[0]],[blocks[1]]],subdivide=False) + krylov.permute_rows(priority_permutation) + return krylov + + def naive_krylov_profile(self, M, shift=None, degree=None, output_pairs=True): r""" Compute the rank profile (row and column) of the striped Krylov matrix - built from ``self`` and matrix `J`. + built from ``self`` and matrix ``M``. INPUT: - - ``J`` -- square matrix used in the Krylov construction. - - ``degree`` -- maximum power of J in Krylov matrix. + - ``M`` -- square matrix used in the Krylov construction. + - ``degree`` -- maximum power of ``M`` in Krylov matrix. - ``shift`` -- list of integers (optional): priority row shift. OUTPUT: + - A tuple (row_profile, pivot_matrix, column_profile): + * row_profile: list of the first r row indices in the striped Krylov matrix ``K`` corresponding to independent rows + * pivot: the submatrix of ``K`` given by those rows + * column_profile: list of the first r independent column indices in ``pivot_matrix`` + where r is the rank of ``K``. + + TESTS:: + + sage: R = GF(97) + sage: E = matrix(R,[[27,49,29],[50,58,0],[77,10,29]]) + sage: E + [27 49 29] + [50 58 0] + [77 10 29] + sage: M = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) + sage: M + [0 1 0] + [0 0 1] + [0 0 0] + sage: degree = 3 + sage: striped_krylov_matrix = matrix.block([[E*M^i] for i in range(degree+1)],subdivide=False) + sage: striped_krylov_matrix + [27 49 29] + [50 58 0] + [77 10 29] + [ 0 27 49] + [ 0 50 58] + [ 0 77 10] + [ 0 0 27] + [ 0 0 50] + [ 0 0 77] + [ 0 0 0] + [ 0 0 0] + [ 0 0 0] + sage: E.krylov_profile(M, output_pairs=False) + ( + [27 49 29] + [50 58 0] + (0, 1, 3), [ 0 27 49], (0, 1, 2) + ) + sage: E.krylov_profile(M) + ( + [27 49 29] + [50 58 0] + ((0, 0), (1, 0), (0, 1)), [ 0 27 49], (0, 1, 2) + ) + + sage: rows = [(i+j*E.nrows()+1,i,j) for j in range(degree+1) for i in range(E.nrows())] + sage: rows + [(1, 0, 0), + (2, 1, 0), + (3, 2, 0), + (4, 0, 1), + (5, 1, 1), + (6, 2, 1), + (7, 0, 2), + (8, 1, 2), + (9, 2, 2), + (10, 0, 3), + (11, 1, 3), + (12, 2, 3)] + sage: shift = [0,3,6] + sage: rows.sort(key=lambda x: (x[2] + shift[x[1]],x[1])) + sage: rows + [(1, 0, 0), + (4, 0, 1), + (7, 0, 2), + (10, 0, 3), + (2, 1, 0), + (5, 1, 1), + (8, 1, 2), + (11, 1, 3), + (3, 2, 0), + (6, 2, 1), + (9, 2, 2), + (12, 2, 3)] + sage: striped_krylov_matrix.with_permuted_rows(Permutation([x[0] for x in rows])) + [27 49 29] + [ 0 27 49] + [ 0 0 27] + [ 0 0 0] + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [ 0 0 0] + [77 10 29] + [ 0 77 10] + [ 0 0 77] + [ 0 0 0] + sage: E.krylov_profile(M, shift=shift, output_pairs=False) + ( + [27 49 29] + [ 0 27 49] + (0, 1, 2), [ 0 0 27], (0, 1, 2) + ) + sage: E.krylov_profile(M, shift=shift) + ( + [27 49 29] + [ 0 27 49] + ((0, 0), (0, 1), (0, 2)), [ 0 0 27], (0, 1, 2) + ) + + sage: shift = [3,0,2] + sage: rows.sort(key=lambda x: (x[2] + shift[x[1]],x[1])) + sage: rows + [(2, 1, 0), + (5, 1, 1), + (8, 1, 2), + (3, 2, 0), + (1, 0, 0), + (11, 1, 3), + (6, 2, 1), + (4, 0, 1), + (9, 2, 2), + (7, 0, 2), + (12, 2, 3), + (10, 0, 3)] + sage: striped_krylov_matrix.with_permuted_rows(Permutation([x[0] for x in rows])) + [50 58 0] + [ 0 50 58] + [ 0 0 50] + [77 10 29] + [27 49 29] + [ 0 0 0] + [ 0 77 10] + [ 0 27 49] + [ 0 0 77] + [ 0 0 27] + [ 0 0 0] + [ 0 0 0] + sage: E.krylov_profile(M,shift=shift,output_pairs=False) + ( + [50 58 0] + [ 0 50 58] + (0, 1, 2), [ 0 0 50], (0, 1, 2) + ) + sage: E.krylov_profile(M,shift=shift) + ( + [50 58 0] + [ 0 50 58] + ((1, 0), (1, 1), (1, 2)), [ 0 0 50], (0, 1, 2) + ) + """ + from sage.matrix.constructor import matrix + + if not isinstance(M, Matrix): + raise TypeError("naive_krylov_profile: M is not a matrix") + if M.nrows() != self.ncols() or M.ncols() != self.ncols(): + raise ValueError("naive_krylov_profile: matrix M does not have correct dimensions.") + if M.base_ring() != self.base_ring(): + raise ValueError("naive_krylov_profile: matrix M does not have same base ring as E.") + + if degree is None: + degree = self.ncols() + if not isinstance(degree, (int, sage.rings.integer.Integer)): + raise TypeError("naive_krylov_profile: degree is not an integer.") + if degree < 0: + raise ValueError("naive_krylov_profile: degree bound cannot be negative.") + + if shift is None or shift == 0: + shift = (ZZ**self.nrows()).zero() + if isinstance(shift, (list, tuple)): + shift = (ZZ**self.nrows())(shift) + if shift not in ZZ**self.nrows(): + raise ValueError(f"naive_krylov_profile: shift is not an integer vector of length {self.nrows()}.") + + m = self.nrows() + + K = self.striped_krylov_matrix(M,shift,degree) + + # define priority and indexing functions + # index: [0,1,...,m-1,m,...,2m-1,...,m*degree-1] -> [(0,0),(1,0),...,(m-1,0),(0,1),...,(m-1,1),...,(m-1,degree)] + priority = lambda c,d : shift[c] + d + index = lambda i : (i%m,i//m) + + # deduce permutation by sorting priorities + # rows should be ordered by ascending priority, then by associated row number in E (i%m) + priority_pairs = sorted([[priority(*index(i)),index(i)] for i in range(m*(degree+1))]) + + # calculate first independent rows of K + row_profile = K.pivot_rows() + + # construct submatrix + pivot = K.matrix_from_rows(row_profile) + + # calculate first independent columns of pivot + col_profile = pivot.pivots() + + row_profile = tuple([priority_pairs[i][1] for i in row_profile]) + + return row_profile,pivot,col_profile + + def krylov_profile(self, M, shift=None, degree=None, output_pairs=True): + r""" + Compute the rank profile (row and column) of the block Krylov matrix + built from ``self`` and matrix ``M``. + + INPUT: + + - ``M`` -- square matrix used in the Krylov construction. + - ``degree`` -- maximum power of ``M`` in Krylov matrix. If None, a suitable upper bound of ``self.ncols()`` is default. + - ``shift`` -- list of ``self.nrows()`` integers (optional): priority + row shift. + + OUTPUT: + A tuple (row_profile, pivot_matrix, column_profile): - - ``row_profile``: list of the first r row indices in the striped - Krylov matrix ``K`` corresponding to - independent rows. If ``output_pairs`` is True, - then the output is the list of pairs (c_i, d_i) - where c_i is the corresponding row in the - original matrix and d_i the corresponding power - of ``J``. - - ``pivot_matrix``: the submatrix of ``K`` given by ``row_profile`` - - ``column_profile``: the first r independent column indices in - ``pivot_matrix`` - where r is the rank of `K`. + - ``row_profile``: list of the first r row indices in the striped + Krylov matrix ``K`` corresponding to independent rows. If + ``output_pairs`` is True, then the output is the list of pairs + (c_i, d_i) where c_i is the corresponding row in the original matrix + and d_i the corresponding power of ``M``. + - ``pivot_matrix``: the submatrix of ``K`` given by ``row_profile`` + - ``column_profile``: the first r independent column indices in + ``pivot_matrix`` where r is the rank of `K`. TESTS:: @@ -19137,13 +19484,13 @@ cdef class Matrix(Matrix1): [27 49 29] [50 58 0] [77 10 29] - sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) - sage: J + sage: M = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) + sage: M [0 1 0] [0 0 1] [0 0 0] sage: degree = 3 - sage: striped_krylov_matrix = matrix.block([[E*J^i] for i in range(degree+1)],subdivide=False) + sage: striped_krylov_matrix = matrix.block([[E*M^i] for i in range(degree+1)],subdivide=False) sage: striped_krylov_matrix [27 49 29] [50 58 0] @@ -19157,13 +19504,13 @@ cdef class Matrix(Matrix1): [ 0 0 0] [ 0 0 0] [ 0 0 0] - sage: E.krylov_rank_profile(J) + sage: E.krylov_profile(M, output_pairs=False) ( [27 49 29] [50 58 0] (0, 1, 3), [ 0 27 49], (0, 1, 2) ) - sage: E.krylov_rank_profile(J,output_pairs=True) + sage: E.krylov_profile(M) ( [27 49 29] [50 58 0] @@ -19212,13 +19559,13 @@ cdef class Matrix(Matrix1): [ 0 77 10] [ 0 0 77] [ 0 0 0] - sage: E.krylov_rank_profile(J,shift=shift) + sage: E.krylov_profile(M, shift=shift, output_pairs=False) ( [27 49 29] [ 0 27 49] (0, 1, 2), [ 0 0 27], (0, 1, 2) ) - sage: E.krylov_rank_profile(J,shift=shift,output_pairs=True) + sage: E.krylov_profile(M, shift=shift) ( [27 49 29] [ 0 27 49] @@ -19253,13 +19600,13 @@ cdef class Matrix(Matrix1): [ 0 0 27] [ 0 0 0] [ 0 0 0] - sage: E.krylov_rank_profile(J,shift=shift) + sage: E.krylov_profile(M,shift=shift,output_pairs=False) ( [50 58 0] [ 0 50 58] (0, 1, 2), [ 0 0 50], (0, 1, 2) ) - sage: E.krylov_rank_profile(J,shift=shift,output_pairs=True) + sage: E.krylov_profile(M,shift=shift) ( [50 58 0] [ 0 50 58] @@ -19270,26 +19617,26 @@ cdef class Matrix(Matrix1): import math from sage.combinat.permutation import Permutation - if not isinstance(J, Matrix): - raise TypeError("krylov_rank_profile: J is not a matrix") - if J.nrows() != self.ncols() or J.ncols() != self.ncols(): - raise ValueError("krylov_rank_profile: matrix J does not have correct dimensions.") - if J.base_ring() != self.base_ring(): - raise ValueError("krylov_rank_profile: matrix J does not have same base ring as E.") + if not isinstance(M, Matrix): + raise TypeError("krylov_profile: M is not a matrix") + if M.nrows() != self.ncols() or M.ncols() != self.ncols(): + raise ValueError("krylov_profile: matrix M does not have correct dimensions.") + if M.base_ring() != self.base_ring(): + raise ValueError("krylov_profile: matrix M does not have same base ring as E.") if degree is None: degree = self.ncols() if not isinstance(degree, (int, sage.rings.integer.Integer)): - raise TypeError("krylov_rank_profile: degree is not an integer.") + raise TypeError("krylov_profile: degree is not an integer.") if degree < 0: - raise ValueError("krylov_rank_profile: degree bound cannot be negative.") + raise ValueError("krylov_profile: degree bound cannot be negative.") if shift is None or shift == 0: shift = (ZZ**self.nrows()).zero() if isinstance(shift, (list, tuple)): shift = (ZZ**self.nrows())(shift) if shift not in ZZ**self.nrows(): - raise ValueError(f"krylov_rank_profile: shift is not a Z-vector of length {self.nrows()}.") + raise ValueError(f"krylov_profile: shift is not a Z-vector of length {self.nrows()}.") m = self.nrows() sigma = self.ncols() @@ -19313,11 +19660,11 @@ cdef class Matrix(Matrix1): if r == 0: return (),self.matrix_from_rows([]),() - M = self_permuted.matrix_from_rows(row_profile_self_permuted) + R = self_permuted.matrix_from_rows(row_profile_self_permuted) - J_L = None + M_L = None - for l in range(math.ceil(math.log(degree,2)) + 1): + for l in range(math.ceil(math.log(degree + 1, 2))): L = pow(2,l) # adding 2^l to each degree row_extension = [(x[0],x[1] + L) for x in row_profile_self if x[1] + L <= degree] @@ -19331,67 +19678,67 @@ cdef class Matrix(Matrix1): k_perm = Permutation(sorted([i+1 for i in range(len(k))],key=lambda x: (shift[k[x-1][0]] + k[x-1][1],k[x-1][0]))) # fast calculation of rows formed by indices in k - if J_L is None: - J_L = J + if M_L is None: + M_L = M else: - J_L = J_L * J_L + M_L = M_L * M_L - M = matrix.block([[exhausted],[M],[M*J_L]],subdivide=False) + R = matrix.block([[exhausted],[R],[R*M_L]],subdivide=False) - # sort rows of M, find profile, translate to k (indices of full krylov matrix) - M.permute_rows(k_perm) + # sort rows of R, find profile, translate to k (indices of full krylov matrix) + R.permute_rows(k_perm) - row_profile_M = M.pivot_rows() - r = len(row_profile_M) + row_profile_R = R.pivot_rows() + r = len(row_profile_R) if r == sigma and l < math.ceil(math.log(degree,2)): - tail = list(range(row_profile_M[-1]+1,M.nrows())) + tail = list(range(row_profile_R[-1]+1,R.nrows())) excluded_rows.update(set([k[k_perm(i+1)-1][0] for i in tail if i < len(k)])) - xmi = [i for i in row_profile_M if k[k_perm(i+1)-1][0] in excluded_rows] - imi = [i for i in row_profile_M if k[k_perm(i+1)-1][0] not in excluded_rows] + xmi = [i for i in row_profile_R if k[k_perm(i+1)-1][0] in excluded_rows] + imi = [i for i in row_profile_R if k[k_perm(i+1)-1][0] not in excluded_rows] row_profile_exhausted = [k[k_perm(i+1)-1] for i in xmi] row_profile_self = [k[k_perm(i+1)-1] for i in imi] - exhausted = M.matrix_from_rows(xmi) - M = M.matrix_from_rows(imi) + exhausted = R.matrix_from_rows(xmi) + R = R.matrix_from_rows(imi) else: if l == math.ceil(math.log(degree,2)): row_profile_exhausted = [] exhausted = matrix.zero(self.base_ring(), 0, sigma) - row_profile_self = [k[k_perm(i+1)-1] for i in row_profile_M] - # calculate new M for return value or next loop - M = M.matrix_from_rows(row_profile_M) - if M.nrows() == 0: - M = exhausted + row_profile_self = [k[k_perm(i+1)-1] for i in row_profile_R] + # calculate new R for return value or next loop + R = R.matrix_from_rows(row_profile_R) + if R.nrows() == 0: + R = exhausted row_profile_self = row_profile_exhausted elif exhausted.nrows() != 0: k = row_profile_exhausted + row_profile_self - M = exhausted.stack(M) + R = exhausted.stack(R) k_perm = Permutation(sorted([i+1 for i in range(len(k))],key=lambda x: (shift[k[x-1][0]] + k[x-1][1],k[x-1][0]))) - M.permute_rows(k_perm) + R.permute_rows(k_perm) row_profile_self = [k[k_perm(i+1)-1] for i in range(len(k))] - col_profile = M.pivots() + col_profile = R.pivots() if output_pairs: - return tuple(row_profile_self), M, col_profile + return tuple(row_profile_self), R, col_profile # convert c,d to actual position in striped Krylov matrix phi = lambda c,d : sum(min(max(shift[c] - shift[i] + d + (i < c and shift[i] <= shift[c] + degree),0),degree+1) for i in range(m)) row_profile_K = [phi(*row) for row in row_profile_self] - return tuple(row_profile_K), M, col_profile + return tuple(row_profile_K), R, col_profile - def linear_interpolation_basis(self, J, var_name, degree=None, shift=None): + def linear_interpolation_basis(self, M, var_name, degree=None, shift=None): r""" Return a shifted minimal linear interpolation basis for - (``self``,``J``) in ``s``-Popov form with respect to a shift ``s``. + (``self``,``M``) in ``s``-Popov form with respect to a shift ``s``. - Given a K[x]-module V defined by multiplication matrix J, i.e. elements + Given a K[x]-module V defined by multiplication matrix ``M``, i.e. elements are vectors of length n in K, scalars are polynomials in K[x], and - scalar multiplication is defined by p v := v p(J), ``self`` represents + scalar multiplication is defined by p v := v p(``M``), ``self`` represents ``self.nrows()`` elements e_1, ..., e_m in V: A linear interpolation basis is a set of tuples of length @@ -19402,9 +19749,9 @@ cdef class Matrix(Matrix1): INPUT: - - ``J`` -- square matrix for Krylov construction. + - ``M`` -- square matrix for Krylov construction. - ``var_name`` -- variable name for the returned polynomial matrix - - ``degree`` -- upper bound on degree of minpoly(`J`), power of + - ``degree`` -- upper bound on degree of minpoly(`M`). If None, a suitable upper bound of ``self.ncols()`` is default. - ``shift`` -- list of self.nrows() integers (optional): controls priority of rows. @@ -19421,13 +19768,13 @@ cdef class Matrix(Matrix1): [27 49 29] [50 58 0] [77 10 29] - sage: J = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) - sage: J + sage: M = matrix(R,[[0,1,0],[0,0,1],[0,0,0]]) + sage: M [0 1 0] [0 0 1] [0 0 0] sage: degree = E.ncols() - sage: striped_krylov_matrix = matrix.block([[E*J^i] for i in range(degree+1)],subdivide=False) + sage: striped_krylov_matrix = matrix.block([[E*M^i] for i in range(degree+1)],subdivide=False) sage: striped_krylov_matrix [27 49 29] [50 58 0] @@ -19441,7 +19788,7 @@ cdef class Matrix(Matrix1): [ 0 0 0] [ 0 0 0] [ 0 0 0] - sage: basis = E.linear_interpolation_basis(J,'x') + sage: basis = E.linear_interpolation_basis(M,'x') sage: basis [x^2 + 40*x + 82 76 0] [ 3*x + 13 x + 57 0] @@ -19452,11 +19799,11 @@ cdef class Matrix(Matrix1): True sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) True - sage: len(E.krylov_rank_profile(J)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + sage: len(E.krylov_profile(M)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) True sage: shift = [0,3,6] - sage: basis = E.linear_interpolation_basis(J,'x',shift=shift) + sage: basis = E.linear_interpolation_basis(M,'x',shift=shift) sage: basis [ x^3 0 0] [60*x^2 + 72*x + 70 1 0] @@ -19467,11 +19814,11 @@ cdef class Matrix(Matrix1): True sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) True - sage: len(E.krylov_rank_profile(J,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + sage: len(E.krylov_profile(M,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) True sage: shift = [3,0,2] - sage: basis = E.linear_interpolation_basis(J,'x',shift=shift) + sage: basis = E.linear_interpolation_basis(M,'x',shift=shift) sage: basis [ 1 26*x^2 + 49*x + 79 0] [ 0 x^3 0] @@ -19482,7 +19829,7 @@ cdef class Matrix(Matrix1): True sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) True - sage: len(E.krylov_rank_profile(J,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + sage: len(E.krylov_profile(M,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) True """ from sage.combinat.permutation import Permutation @@ -19490,12 +19837,12 @@ cdef class Matrix(Matrix1): from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing # INPUT VALIDATION - if not isinstance(J, Matrix): - raise TypeError("krylov_rank_profile: J is not a matrix") - if J.nrows() != self.ncols() or J.ncols() != self.ncols(): - raise ValueError("krylov_rank_profile: matrix J does not have correct dimensions.") - if J.base_ring() != self.base_ring(): - raise ValueError("krylov_rank_profile: matrix J does not have same base ring as E.") + if not isinstance(M, Matrix): + raise TypeError("linear_interpolation_basis: M is not a matrix") + if M.nrows() != self.ncols() or M.ncols() != self.ncols(): + raise ValueError("linear_interpolation_basis: matrix M does not have correct dimensions.") + if M.base_ring() != self.base_ring(): + raise ValueError("linear_interpolation_basis: matrix M does not have same base ring as E.") if not isinstance(var_name, str): raise TypeError("var_name is not a string") @@ -19503,16 +19850,16 @@ cdef class Matrix(Matrix1): if degree is None: degree = self.ncols() if not isinstance(degree, (int, sage.rings.integer.Integer)): - raise TypeError("krylov_rank_profile: degree is not an integer.") + raise TypeError("linear_interpolation_basis: degree is not an integer.") if degree < 0: - raise ValueError("krylov_rank_profile: degree bound cannot be negative.") + raise ValueError("linear_interpolation_basis: degree bound cannot be negative.") if shift is None or shift == 0: shift = (ZZ**self.nrows()).zero() if isinstance(shift, (list, tuple)): shift = (ZZ**self.nrows())(shift) if shift not in ZZ**self.nrows(): - raise ValueError(f"krylov_rank_profile: shift is not a Z-vector of length {self.nrows()}.") + raise ValueError(f"linear_interpolation_basis: shift is not a Z-vector of length {self.nrows()}.") m = self.nrows() sigma = self.ncols() @@ -19520,7 +19867,7 @@ cdef class Matrix(Matrix1): poly_ring = PolynomialRing(self.base_ring(),var_name) # calculate krylov profile - row_profile, pivot, col_profile = self.krylov_rank_profile(J,degree,shift,output_pairs=True) + row_profile, pivot, col_profile = self.krylov_profile(M, shift, degree) if len(row_profile) == 0: return matrix.identity(poly_ring,m) @@ -19538,7 +19885,7 @@ cdef class Matrix(Matrix1): T_absent_indices = [i for i in range(m) if inv_c[i] is None] T_present_indices = [inv_c[i] for i in range(m) if inv_c[i] is not None] D_absent = self.matrix_from_rows_and_columns(T_absent_indices, col_profile) - D_present = (pivot.matrix_from_rows(T_present_indices)*J).matrix_from_columns(col_profile) + D_present = (pivot.matrix_from_rows(T_present_indices)*M).matrix_from_columns(col_profile) D_rows = [] idx_p = 0 idx_a = 0 From 6b0a721a5f116a38eeefbfa19cedc5d3fcab43d4 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Wed, 6 Aug 2025 16:51:00 +0200 Subject: [PATCH 26/28] implement coefficient output, various doc changes --- src/sage/matrix/matrix2.pyx | 354 ++++++++++++++++++++++-------------- 1 file changed, 215 insertions(+), 139 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index b5105422411..44965f81c4d 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19111,9 +19111,10 @@ cdef class Matrix(Matrix1): U.rescale_col(n - 1, -1) return U - def striped_krylov_matrix(self, M, shift=None, degree=None): + def krylov_matrix(self, M, shift=None, degree=None): r""" - Return the block matrix of the following form, with rows permuted according to shift. + Return the block matrix built from the rows of ``self`` and the matrix ``M``, with rows ordered according to a priority defined by ``shift``. See [Beckermann and Labahn, 2000, https://doi.org/10.1137/S0895479897326912] and [Jeannerod, Neiger, Schost, Villard, 2017, https://doi.org/10.1016/j.jsc.2016.11.015]." + The following uses `E` to refer to ``self``, and `d` to refer to ``degree``. [ ] @@ -19136,12 +19137,13 @@ cdef class Matrix(Matrix1): - ``M`` -- a square matrix of size equal to the number of columns of ``self``. - ``degree`` -- integer, the maximum exponent for the Krylov matrix. - - ``shift`` -- list of integers (optional), row priority shifts. If ``None``, - defaults to all zero. + - ``shift`` -- list of ``self.nrows()`` integers (optional), row priority + shifts. If ``None``, defaults to all zeroes. OUTPUT: - - A matrix with block rows [E, EM, EM^2, ..., EM^d], row-permuted by shift. + - A matrix with block rows [E, EM, EM^2, ..., EM^d], row-permuted by + shift. EXAMPLES:: @@ -19156,7 +19158,7 @@ cdef class Matrix(Matrix1): [0 1 0] [0 0 1] [0 0 0] - sage: E.striped_krylov_matrix(M) + sage: E.krylov_matrix(M) [27 49 29] [50 58 0] [77 10 29] @@ -19169,7 +19171,7 @@ cdef class Matrix(Matrix1): [ 0 0 0] [ 0 0 0] [ 0 0 0] - sage: E.striped_krylov_matrix(M,[0,3,6]) + sage: E.krylov_matrix(M,[0,3,6]) [27 49 29] [ 0 27 49] [ 0 0 27] @@ -19182,7 +19184,7 @@ cdef class Matrix(Matrix1): [ 0 77 10] [ 0 0 77] [ 0 0 0] - sage: E.striped_krylov_matrix(M,[3,0,2]) + sage: E.krylov_matrix(M,[3,0,2]) [50 58 0] [ 0 50 58] [ 0 0 50] @@ -19202,25 +19204,25 @@ cdef class Matrix(Matrix1): from sage.matrix.constructor import matrix # for matrix.block if not isinstance(M, Matrix): - raise TypeError("striped_krylov_matrix: M is not a matrix") + raise TypeError("krylov_matrix: M is not a matrix") if M.nrows() != self.ncols() or M.ncols() != self.ncols(): - raise ValueError("striped_krylov_matrix: matrix M does not have correct dimensions.") + raise ValueError("krylov_matrix: matrix M does not have correct dimensions.") if M.base_ring() != self.base_ring(): - raise ValueError("striped_krylov_matrix: matrix M does not have same base ring as E.") + raise ValueError("krylov_matrix: matrix M does not have same base ring as E.") if degree is None: degree = self.ncols() if not isinstance(degree, (int, sage.rings.integer.Integer)): - raise TypeError("striped_krylov_matrix: degree is not an integer.") + raise TypeError("krylov_matrix: degree is not an integer.") if degree < 0: - raise ValueError("striped_krylov_matrix: degree bound cannot be negative.") + raise ValueError("krylov_matrix: degree bound cannot be negative.") if shift is None or shift == 0: shift = (ZZ**self.nrows()).zero() if isinstance(shift, (list, tuple)): shift = (ZZ**self.nrows())(shift) if shift not in ZZ**self.nrows(): - raise ValueError(f"striped_krylov_matrix: shift is not an integer vector of length {self.nrows()}.") + raise ValueError(f"krylov_matrix: shift is not an integer vector of length {self.nrows()}.") m = self.nrows() @@ -19234,10 +19236,10 @@ cdef class Matrix(Matrix1): priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) - # build all blocks for the striped matrix + # store 2 blocks for the Krylov matrix blocks = [self] - # for i from 0 to degree (inclusive), add E*M^i + # for i from 0 to degree (inclusive), merge existing blocks, add E*M^i M_i = None for i in range(math.ceil(math.log(degree + 1, 2))): if M_i is None: @@ -19255,21 +19257,21 @@ cdef class Matrix(Matrix1): krylov.permute_rows(priority_permutation) return krylov - def naive_krylov_profile(self, M, shift=None, degree=None, output_pairs=True): + def _naive_krylov_basis(self, M, shift=None, degree=None, output_pairs=True): r""" - Compute the rank profile (row and column) of the striped Krylov matrix + Compute the rank profile (row and column) of the Krylov matrix built from ``self`` and matrix ``M``. INPUT: - ``M`` -- square matrix used in the Krylov construction. - - ``degree`` -- maximum power of ``M`` in Krylov matrix. - - ``shift`` -- list of integers (optional): priority row shift. + - ``degree`` -- maximum power of ``M`` in Krylov matrix. If None, a suitable upper bound of ``self.ncols()`` is default. + - ``shift`` -- list of ``self.nrows()`` integers (optional): priority row shift. If ``None``, defaults to all zeroes. OUTPUT: - A tuple (row_profile, pivot_matrix, column_profile): - * row_profile: list of the first r row indices in the striped Krylov matrix ``K`` corresponding to independent rows + * row_profile: list of the first r row indices in the Krylov matrix ``K`` corresponding to independent rows * pivot: the submatrix of ``K`` given by those rows * column_profile: list of the first r independent column indices in ``pivot_matrix`` where r is the rank of ``K``. @@ -19288,8 +19290,8 @@ cdef class Matrix(Matrix1): [0 0 1] [0 0 0] sage: degree = 3 - sage: striped_krylov_matrix = matrix.block([[E*M^i] for i in range(degree+1)],subdivide=False) - sage: striped_krylov_matrix + sage: K = E.krylov_matrix(M) + sage: K [27 49 29] [50 58 0] [77 10 29] @@ -19302,13 +19304,13 @@ cdef class Matrix(Matrix1): [ 0 0 0] [ 0 0 0] [ 0 0 0] - sage: E.krylov_profile(M, output_pairs=False) + sage: E.krylov_basis(M, output_pairs=False) ( [27 49 29] [50 58 0] (0, 1, 3), [ 0 27 49], (0, 1, 2) ) - sage: E.krylov_profile(M) + sage: E.krylov_basis(M) ( [27 49 29] [50 58 0] @@ -19344,7 +19346,7 @@ cdef class Matrix(Matrix1): (6, 2, 1), (9, 2, 2), (12, 2, 3)] - sage: striped_krylov_matrix.with_permuted_rows(Permutation([x[0] for x in rows])) + sage: E.krylov_matrix(M,shift=shift) [27 49 29] [ 0 27 49] [ 0 0 27] @@ -19357,13 +19359,13 @@ cdef class Matrix(Matrix1): [ 0 77 10] [ 0 0 77] [ 0 0 0] - sage: E.krylov_profile(M, shift=shift, output_pairs=False) + sage: E.krylov_basis(M, shift=shift, output_pairs=False) ( [27 49 29] [ 0 27 49] (0, 1, 2), [ 0 0 27], (0, 1, 2) ) - sage: E.krylov_profile(M, shift=shift) + sage: E.krylov_basis(M,shift=shift) ( [27 49 29] [ 0 27 49] @@ -19385,7 +19387,7 @@ cdef class Matrix(Matrix1): (7, 0, 2), (12, 2, 3), (10, 0, 3)] - sage: striped_krylov_matrix.with_permuted_rows(Permutation([x[0] for x in rows])) + sage: E.krylov_matrix(M,shift=shift) [50 58 0] [ 0 50 58] [ 0 0 50] @@ -19398,13 +19400,13 @@ cdef class Matrix(Matrix1): [ 0 0 27] [ 0 0 0] [ 0 0 0] - sage: E.krylov_profile(M,shift=shift,output_pairs=False) + sage: E.krylov_basis(M,shift=shift,output_pairs=False) ( [50 58 0] [ 0 50 58] (0, 1, 2), [ 0 0 50], (0, 1, 2) ) - sage: E.krylov_profile(M,shift=shift) + sage: E.krylov_basis(M,shift=shift) ( [50 58 0] [ 0 50 58] @@ -19414,38 +19416,29 @@ cdef class Matrix(Matrix1): from sage.matrix.constructor import matrix if not isinstance(M, Matrix): - raise TypeError("naive_krylov_profile: M is not a matrix") + raise TypeError("naive_krylov_basis: M is not a matrix") if M.nrows() != self.ncols() or M.ncols() != self.ncols(): - raise ValueError("naive_krylov_profile: matrix M does not have correct dimensions.") + raise ValueError("naive_krylov_basis: matrix M does not have correct dimensions.") if M.base_ring() != self.base_ring(): - raise ValueError("naive_krylov_profile: matrix M does not have same base ring as E.") + raise ValueError("naive_krylov_basis: matrix M does not have same base ring as E.") if degree is None: degree = self.ncols() if not isinstance(degree, (int, sage.rings.integer.Integer)): - raise TypeError("naive_krylov_profile: degree is not an integer.") + raise TypeError("naive_krylov_basis: degree is not an integer.") if degree < 0: - raise ValueError("naive_krylov_profile: degree bound cannot be negative.") + raise ValueError("naive_krylov_basis: degree bound cannot be negative.") if shift is None or shift == 0: shift = (ZZ**self.nrows()).zero() if isinstance(shift, (list, tuple)): shift = (ZZ**self.nrows())(shift) if shift not in ZZ**self.nrows(): - raise ValueError(f"naive_krylov_profile: shift is not an integer vector of length {self.nrows()}.") + raise ValueError(f"naive_krylov_basis: shift is not an integer vector of length {self.nrows()}.") m = self.nrows() - K = self.striped_krylov_matrix(M,shift,degree) - - # define priority and indexing functions - # index: [0,1,...,m-1,m,...,2m-1,...,m*degree-1] -> [(0,0),(1,0),...,(m-1,0),(0,1),...,(m-1,1),...,(m-1,degree)] - priority = lambda c,d : shift[c] + d - index = lambda i : (i%m,i//m) - - # deduce permutation by sorting priorities - # rows should be ordered by ascending priority, then by associated row number in E (i%m) - priority_pairs = sorted([[priority(*index(i)),index(i)] for i in range(m*(degree+1))]) + K = self.krylov_matrix(M,shift,degree) # calculate first independent rows of K row_profile = K.pivot_rows() @@ -19456,11 +19449,22 @@ cdef class Matrix(Matrix1): # calculate first independent columns of pivot col_profile = pivot.pivots() - row_profile = tuple([priority_pairs[i][1] for i in row_profile]) + if output_pairs: + + # define priority and indexing functions + # index: [0,1,...,m-1,m,...,2m-1,...,m*degree-1] -> [(0,0),(1,0),...,(m-1,0),(0,1),...,(m-1,1),...,(m-1,degree)] + priority = lambda c,d : shift[c] + d + index = lambda i : (i%m,i//m) + + # deduce permutation by sorting priorities + # rows should be ordered by ascending priority, then by associated row number in E (i%m) + priority_pairs = sorted([[priority(*index(i)),index(i)] for i in range(m*(degree+1))]) + + row_profile = tuple([priority_pairs[i][1] for i in row_profile]) return row_profile,pivot,col_profile - def krylov_profile(self, M, shift=None, degree=None, output_pairs=True): + def _elimination_krylov_basis(self, M, shift=None, degree=None, output_pairs=True): r""" Compute the rank profile (row and column) of the block Krylov matrix built from ``self`` and matrix ``M``. @@ -19470,7 +19474,7 @@ cdef class Matrix(Matrix1): - ``M`` -- square matrix used in the Krylov construction. - ``degree`` -- maximum power of ``M`` in Krylov matrix. If None, a suitable upper bound of ``self.ncols()`` is default. - ``shift`` -- list of ``self.nrows()`` integers (optional): priority - row shift. + row shift. If ``None``, defaults to all zeroes. OUTPUT: @@ -19498,8 +19502,8 @@ cdef class Matrix(Matrix1): [0 0 1] [0 0 0] sage: degree = 3 - sage: striped_krylov_matrix = matrix.block([[E*M^i] for i in range(degree+1)],subdivide=False) - sage: striped_krylov_matrix + sage: K = E.krylov_matrix(M) + sage: K [27 49 29] [50 58 0] [77 10 29] @@ -19512,19 +19516,18 @@ cdef class Matrix(Matrix1): [ 0 0 0] [ 0 0 0] [ 0 0 0] - sage: E.krylov_profile(M, output_pairs=False) + sage: E.krylov_basis(M,output_pairs=False) ( [27 49 29] [50 58 0] (0, 1, 3), [ 0 27 49], (0, 1, 2) ) - sage: E.krylov_profile(M) + sage: E.krylov_basis(M) ( [27 49 29] [50 58 0] ((0, 0), (1, 0), (0, 1)), [ 0 27 49], (0, 1, 2) ) - sage: rows = [(i+j*E.nrows()+1,i,j) for j in range(degree+1) for i in range(E.nrows())] sage: rows [(1, 0, 0), @@ -19554,7 +19557,7 @@ cdef class Matrix(Matrix1): (6, 2, 1), (9, 2, 2), (12, 2, 3)] - sage: striped_krylov_matrix.with_permuted_rows(Permutation([x[0] for x in rows])) + sage: E.krylov_matrix(M,shift=shift) [27 49 29] [ 0 27 49] [ 0 0 27] @@ -19567,19 +19570,18 @@ cdef class Matrix(Matrix1): [ 0 77 10] [ 0 0 77] [ 0 0 0] - sage: E.krylov_profile(M, shift=shift, output_pairs=False) + sage: E.krylov_basis(M,shift=shift,output_pairs=False) ( [27 49 29] [ 0 27 49] (0, 1, 2), [ 0 0 27], (0, 1, 2) ) - sage: E.krylov_profile(M, shift=shift) + sage: E.krylov_basis(M,shift=shift) ( [27 49 29] [ 0 27 49] ((0, 0), (0, 1), (0, 2)), [ 0 0 27], (0, 1, 2) ) - sage: shift = [3,0,2] sage: rows.sort(key=lambda x: (x[2] + shift[x[1]],x[1])) sage: rows @@ -19595,7 +19597,7 @@ cdef class Matrix(Matrix1): (7, 0, 2), (12, 2, 3), (10, 0, 3)] - sage: striped_krylov_matrix.with_permuted_rows(Permutation([x[0] for x in rows])) + sage: E.krylov_matrix(M,shift=shift) [50 58 0] [ 0 50 58] [ 0 0 50] @@ -19608,13 +19610,13 @@ cdef class Matrix(Matrix1): [ 0 0 27] [ 0 0 0] [ 0 0 0] - sage: E.krylov_profile(M,shift=shift,output_pairs=False) + sage: E.krylov_basis(M,shift=shift,output_pairs=False) ( [50 58 0] [ 0 50 58] (0, 1, 2), [ 0 0 50], (0, 1, 2) ) - sage: E.krylov_profile(M,shift=shift) + sage: E.krylov_basis(M,shift=shift) ( [50 58 0] [ 0 50 58] @@ -19626,25 +19628,25 @@ cdef class Matrix(Matrix1): from sage.combinat.permutation import Permutation if not isinstance(M, Matrix): - raise TypeError("krylov_profile: M is not a matrix") + raise TypeError("krylov_basis: M is not a matrix") if M.nrows() != self.ncols() or M.ncols() != self.ncols(): - raise ValueError("krylov_profile: matrix M does not have correct dimensions.") + raise ValueError("krylov_basis: matrix M does not have correct dimensions.") if M.base_ring() != self.base_ring(): - raise ValueError("krylov_profile: matrix M does not have same base ring as E.") + raise ValueError("krylov_basis: matrix M does not have same base ring as E.") if degree is None: degree = self.ncols() if not isinstance(degree, (int, sage.rings.integer.Integer)): - raise TypeError("krylov_profile: degree is not an integer.") + raise TypeError("krylov_basis: degree is not an integer.") if degree < 0: - raise ValueError("krylov_profile: degree bound cannot be negative.") + raise ValueError("krylov_basis: degree bound cannot be negative.") if shift is None or shift == 0: shift = (ZZ**self.nrows()).zero() if isinstance(shift, (list, tuple)): shift = (ZZ**self.nrows())(shift) if shift not in ZZ**self.nrows(): - raise ValueError(f"krylov_profile: shift is not a Z-vector of length {self.nrows()}.") + raise ValueError(f"krylov_basis: shift is not a Z-vector of length {self.nrows()}.") m = self.nrows() sigma = self.ncols() @@ -19739,7 +19741,54 @@ cdef class Matrix(Matrix1): return tuple(row_profile_K), R, col_profile - def linear_interpolation_basis(self, M, var_name, degree=None, shift=None): + def krylov_basis(self, M, shift=None, degree=None, output_pairs=True, algorithm='default'): + r""" + Compute the rank profile (row and column) of the block Krylov matrix + built from ``self`` and matrix ``M``. + + INPUT: + + - ``M`` -- square matrix used in the Krylov construction. + - ``degree`` -- maximum power of ``M`` in Krylov matrix. If None, a suitable upper bound of ``self.ncols()`` is default. + - ``shift`` -- list of ``self.nrows()`` integers (optional): priority + row shift. If ``None``, defaults to all zeroes. + + OUTPUT: + + A tuple (row_profile, pivot_matrix, column_profile): + - ``row_profile``: list of the first r row indices in the striped + Krylov matrix ``K`` corresponding to independent rows. If + ``output_pairs`` is True, then the output is the list of pairs + (c_i, d_i) where c_i is the corresponding row in the original matrix + and d_i the corresponding power of ``M``. + - ``pivot_matrix``: the submatrix of ``K`` given by ``row_profile`` + - ``column_profile``: the first r independent column indices in + ``pivot_matrix`` where r is the rank of `K`. + """ + if algorithm == 'default': + if self.base_ring().order() == 2: + if self.nrows() <= 12: + algorithm = 'naive' + else: + algorithm = 'elimination' + elif self.base_ring().degree() == 1: + if self.nrows() <= 4: + algorithm = 'naive' + else: + algorithm = 'elimination' + else: + if self.nrows() <= 2: + algorithm = 'naive' + else: + algorithm = 'elimination' + if algorithm == 'naive': + return self._naive_krylov_basis(M, shift, degree, output_pairs) + elif algorithm == 'elimination': + return self._elimination_krylov_basis(M, shift, degree, output_pairs) + else: + raise ValueError("algorithm must be one of \"default\", \"naive\" or \"elimination\".") + + def krylov_kernel_basis(self, M, var_name='x', degree=None, shift=None, output_coefficients=False): r""" Return a shifted minimal linear interpolation basis for (``self``,``M``) in ``s``-Popov form with respect to a shift ``s``. @@ -19753,15 +19802,18 @@ cdef class Matrix(Matrix1): ``self.nrows()`` in K[x] that form a basis for solutions to the equation ``p_1 e_1 + ... + p_n e_n = 0``. - See https://arxiv.org/abs/1512.03503 + See [Kai1980] INPUT: - ``M`` -- square matrix for Krylov construction. - ``var_name`` -- variable name for the returned polynomial matrix - - ``degree`` -- upper bound on degree of minpoly(`M`). If None, a suitable upper bound of ``self.ncols()`` is default. + - ``degree`` -- upper bound on degree of minpoly(`M`). If None, a + suitable upper bound of ``self.ncols()`` is default. - ``shift`` -- list of self.nrows() integers (optional): controls - priority of rows. + priority of rows. If ``None``, defaults to all zeroes. + - ``output_coefficients`` -- If True, outputs the matrix as a sparse + block matrix of coefficients. OUTPUT: @@ -19782,8 +19834,8 @@ cdef class Matrix(Matrix1): [0 0 1] [0 0 0] sage: degree = E.ncols() - sage: striped_krylov_matrix = matrix.block([[E*M^i] for i in range(degree+1)],subdivide=False) - sage: striped_krylov_matrix + sage: krylov_matrix = matrix.block([[E*M^i] for i in range(degree+1)],subdivide=False) + sage: krylov_matrix [27 49 29] [50 58 0] [77 10 29] @@ -19796,48 +19848,63 @@ cdef class Matrix(Matrix1): [ 0 0 0] [ 0 0 0] [ 0 0 0] - sage: basis = E.linear_interpolation_basis(M,'x') + sage: basis = E.krylov_kernel_basis(M,'x') sage: basis [x^2 + 40*x + 82 76 0] [ 3*x + 13 x + 57 0] [ 96 96 1] + sage: basis_coeff = E.krylov_kernel_basis(M,'x',output_coefficients=True) + sage: basis_coeff + [82 76 0 40 0 0 1 0 0] + [13 57 0 3 1 0 0 0 0] + [96 96 1 0 0 0 0 0 0] sage: basis.is_popov() True sage: basis.degree() <= E.ncols() True - sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) + sage: basis_coeff.augment(matrix.zero(R,E.nrows())) * krylov_matrix == matrix.zero(R,E.nrows()) True - sage: len(E.krylov_profile(M)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + sage: len(E.krylov_basis(M)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) True sage: shift = [0,3,6] - sage: basis = E.linear_interpolation_basis(M,'x',shift=shift) + sage: basis = E.krylov_kernel_basis(M,'x',shift=shift) sage: basis [ x^3 0 0] [60*x^2 + 72*x + 70 1 0] [60*x^2 + 72*x + 69 0 1] + sage: basis_coeff = E.krylov_kernel_basis(M,'x',shift=shift,output_coefficients=True) + sage: basis_coeff + [ 0 0 0 0 0 0 0 0 0 1 0 0] + [70 1 0 72 0 0 60 0 0 0 0 0] + [69 0 1 72 0 0 60 0 0 0 0 0] sage: basis.is_popov(shifts=shift) True sage: basis.degree() <= E.ncols() True - sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) + sage: basis_coeff * krylov_matrix == matrix.zero(R,E.nrows()) True - sage: len(E.krylov_profile(M,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + sage: len(E.krylov_basis(M,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) True sage: shift = [3,0,2] - sage: basis = E.linear_interpolation_basis(M,'x',shift=shift) + sage: basis = E.krylov_kernel_basis(M,'x',shift=shift) sage: basis [ 1 26*x^2 + 49*x + 79 0] [ 0 x^3 0] [ 0 26*x^2 + 49*x + 78 1] + sage: basis_coeff = E.krylov_kernel_basis(M,'x',shift=shift,output_coefficients=True) + sage: basis_coeff + [ 1 79 0 0 49 0 0 26 0 0 0 0] + [ 0 0 0 0 0 0 0 0 0 0 1 0] + [ 0 78 1 0 49 0 0 26 0 0 0 0] sage: basis.is_popov(shifts=shift) True sage: basis.degree() <= E.ncols() True - sage: matrix.block([[basis.coefficient_matrix(i) for i in range(E.ncols()+1)]],subdivide=False) * striped_krylov_matrix == matrix.zero(R,E.nrows()) + sage: basis_coeff * krylov_matrix == matrix.zero(R,E.nrows()) True - sage: len(E.krylov_profile(M,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + sage: len(E.krylov_basis(M,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) True """ from sage.combinat.permutation import Permutation @@ -19846,11 +19913,11 @@ cdef class Matrix(Matrix1): # INPUT VALIDATION if not isinstance(M, Matrix): - raise TypeError("linear_interpolation_basis: M is not a matrix") + raise TypeError("krylov_kernel_basis: M is not a matrix") if M.nrows() != self.ncols() or M.ncols() != self.ncols(): - raise ValueError("linear_interpolation_basis: matrix M does not have correct dimensions.") + raise ValueError("krylov_kernel_basis: matrix M does not have correct dimensions.") if M.base_ring() != self.base_ring(): - raise ValueError("linear_interpolation_basis: matrix M does not have same base ring as E.") + raise ValueError("krylov_kernel_basis: matrix M does not have same base ring as E.") if not isinstance(var_name, str): raise TypeError("var_name is not a string") @@ -19858,16 +19925,16 @@ cdef class Matrix(Matrix1): if degree is None: degree = self.ncols() if not isinstance(degree, (int, sage.rings.integer.Integer)): - raise TypeError("linear_interpolation_basis: degree is not an integer.") + raise TypeError("krylov_kernel_basis: degree is not an integer.") if degree < 0: - raise ValueError("linear_interpolation_basis: degree bound cannot be negative.") + raise ValueError("krylov_kernel_basis: degree bound cannot be negative.") if shift is None or shift == 0: shift = (ZZ**self.nrows()).zero() if isinstance(shift, (list, tuple)): shift = (ZZ**self.nrows())(shift) if shift not in ZZ**self.nrows(): - raise ValueError(f"linear_interpolation_basis: shift is not a Z-vector of length {self.nrows()}.") + raise ValueError(f"krylov_kernel_basis: shift is not a Z-vector of length {self.nrows()}.") m = self.nrows() sigma = self.ncols() @@ -19875,14 +19942,14 @@ cdef class Matrix(Matrix1): poly_ring = PolynomialRing(self.base_ring(),var_name) # calculate krylov profile - row_profile, pivot, col_profile = self.krylov_profile(M, shift, degree) + row_profile, pivot, col_profile = self.krylov_basis(M, shift, degree) if len(row_profile) == 0: return matrix.identity(poly_ring,m) c, d = zip(*(row for row in row_profile)) - # degree_c = 0 or max d[k] such that c[k] = c + # degree_c = 0 or max d[k] + 1 such that c[k] = c inv_c = [None]*m degree_c = [0]*m for k in range(len(row_profile)): @@ -19909,57 +19976,66 @@ cdef class Matrix(Matrix1): D = matrix(D_rows) relation = D*C.inverse() - # construct polynomial matrix - potentially costly at polynomial construction stage for extension fields - coeffs_map = [[{} for _ in range(m)] for _ in range(m)] - # add identity part of basis - for i in range(m): - coeffs_map[i][i][degree_c[i]] = self.base_ring().one() - # add relation part of basis - for col in range(relation.ncols()): - for row in range(m): - coeffs_map[row][c[col]][d[col]] = coeffs_map[row][c[col]].get(d[col], self.base_ring().zero()) - relation[row,col] - # construct rows - basis_rows = [[None] * m for _ in range(m)] - if self.base_ring().degree() <= 1 or m*m*sigma < self.base_ring().order(): - # if base field is large, don't bother caching monomials - # or if polynomial construction is efficient (order = 1) - for row in range(m): - for col in range(m): - if not coeffs_map[row][col]: - basis_rows[row][col] = self.base_ring().zero() - else: - # construct a small polynomial (less costly), then shift - mindeg = min(coeffs_map[row][col].keys()) - maxdeg = max(coeffs_map[row][col].keys()) - entries = [self.base_ring().zero()] * (maxdeg-mindeg+1) - for deg, coeff in coeffs_map[row][col].items(): - entries[deg-mindeg] = coeff - basis_rows[row][col] = poly_ring(entries).shift(mindeg) - return matrix(poly_ring, basis_rows) + if output_coefficients: + coefficients = matrix.zero(self.base_ring(),m,m*(max(degree_c)+1)) + for i in range(m): + coefficients[i,i+degree_c[i]*m] = self.base_ring().one() + for col in range(relation.ncols()): + for row in range(m): + coefficients[row,c[col]+d[col]*m] -= relation[row,col] + return coefficients else: - # if base field is small, cache monomials to minimise number of constructions - monomial_cache = {} - for row in range(m): - for col in range(m): - if not coeffs_map[row][col]: - basis_rows[row][col] = poly_ring.zero() - elif len(coeffs_map[row][col]) == 1: - # monomial case - deg, coeff = coeffs_map[row][col].popitem() - if coeff not in monomial_cache: - monomial_cache[coeff] = poly_ring(coeff) - basis_rows[row][col] = monomial_cache[coeff].shift(deg) - else: - # general case - mindeg = min(coeffs_map[row][col].keys()) - poly = poly_ring.zero() - for deg, coeff in coeffs_map[row][col].items(): + # construct polynomial matrix - potentially costly at polynomial construction stage for extension fields + coeffs_map = [[{} for _ in range(m)] for _ in range(m)] + # add identity part of basis + for i in range(m): + coeffs_map[i][i][degree_c[i]] = self.base_ring().one() + # add relation part of basis + for col in range(relation.ncols()): + for row in range(m): + coeffs_map[row][c[col]][d[col]] = coeffs_map[row][c[col]].get(d[col], self.base_ring().zero()) - relation[row,col] + # construct rows + basis_rows = [[None] * m for _ in range(m)] + if self.base_ring().degree() <= 1 or m*m*sigma < self.base_ring().order(): + # if base field is large, don't bother caching monomials + # or if polynomial construction is efficient (order = 1) + for row in range(m): + for col in range(m): + if not coeffs_map[row][col]: + basis_rows[row][col] = self.base_ring().zero() + else: + # construct a small polynomial (less costly), then shift + mindeg = min(coeffs_map[row][col].keys()) + maxdeg = max(coeffs_map[row][col].keys()) + entries = [self.base_ring().zero()] * (maxdeg-mindeg+1) + for deg, coeff in coeffs_map[row][col].items(): + entries[deg-mindeg] = coeff + basis_rows[row][col] = poly_ring(entries).shift(mindeg) + return matrix(poly_ring, basis_rows) + else: + # if base field is small, cache monomials to minimise number of constructions + monomial_cache = {} + for row in range(m): + for col in range(m): + if not coeffs_map[row][col]: + basis_rows[row][col] = poly_ring.zero() + elif len(coeffs_map[row][col]) == 1: + # monomial case + deg, coeff = coeffs_map[row][col].popitem() if coeff not in monomial_cache: monomial_cache[coeff] = poly_ring(coeff) - poly += monomial_cache[coeff].shift(deg-mindeg) - basis_rows[row][col] = poly.shift(mindeg) - # convert to matrix - return matrix(poly_ring,basis_rows) + basis_rows[row][col] = monomial_cache[coeff].shift(deg) + else: + # general case + mindeg = min(coeffs_map[row][col].keys()) + poly = poly_ring.zero() + for deg, coeff in coeffs_map[row][col].items(): + if coeff not in monomial_cache: + monomial_cache[coeff] = poly_ring(coeff) + poly += monomial_cache[coeff].shift(deg-mindeg) + basis_rows[row][col] = poly.shift(mindeg) + # convert to matrix + return matrix(poly_ring,basis_rows) # a limited number of access-only properties are provided for matrices @property From b9389ef5003a3023c9922159cdf0ce35d3f693e8 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Mon, 11 Aug 2025 20:04:32 +0200 Subject: [PATCH 27/28] changes from review, add support for multiple degrees (todo leftover) --- src/sage/matrix/matrix2.pyx | 445 +++++++++++++++++++----------------- 1 file changed, 229 insertions(+), 216 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 44965f81c4d..3efe21af167 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19111,11 +19111,19 @@ cdef class Matrix(Matrix1): U.rescale_col(n - 1, -1) return U - def krylov_matrix(self, M, shift=None, degree=None): + def krylov_matrix(self, M, shifts=None, degree=None): r""" - Return the block matrix built from the rows of ``self`` and the matrix ``M``, with rows ordered according to a priority defined by ``shift``. See [Beckermann and Labahn, 2000, https://doi.org/10.1137/S0895479897326912] and [Jeannerod, Neiger, Schost, Villard, 2017, https://doi.org/10.1016/j.jsc.2016.11.015]." - - The following uses `E` to refer to ``self``, and `d` to refer to ``degree``. + Return the block matrix built from the rows of ``self`` and the matrix + ``M``, with rows ordered according to a priority defined by + ``shifts``. See + [Beckermann and Labahn, 2000, https://doi.org/10.1137/S0895479897326912] + and + [Jeannerod, Neiger, Schost, Villard, 2017, https://doi.org/10.1016/j.jsc.2016.11.015]." + + The following example uses `E` to refer to ``self``, and `d` to refer + to ``degree`` (in the case where a single integer is provided). This is + the case where ``shifts`` is `[0,...,0]``; another shift will + correspond to some row permutation of this matrix. [ ] [ E ] @@ -19135,15 +19143,16 @@ cdef class Matrix(Matrix1): INPUT: - - ``M`` -- a square matrix of size equal to the number of columns of ``self``. + - ``M`` -- a square matrix of size equal to the number of columns of + ``self``. - ``degree`` -- integer, the maximum exponent for the Krylov matrix. - - ``shift`` -- list of ``self.nrows()`` integers (optional), row priority - shifts. If ``None``, defaults to all zeroes. + - ``shifts`` -- list of ``self.nrows()`` integers (optional), row + priority shifts. If ``None``, defaults to all zeroes. OUTPUT: - A matrix with block rows [E, EM, EM^2, ..., EM^d], row-permuted by - shift. + ``shifts``. EXAMPLES:: @@ -19197,11 +19206,11 @@ cdef class Matrix(Matrix1): [ 0 0 27] [ 0 0 0] [ 0 0 0] - """ from sage.combinat.permutation import Permutation import math - from sage.matrix.constructor import matrix # for matrix.block + from sage.matrix.constructor import matrix + from sage.modules.free_module_element import vector if not isinstance(M, Matrix): raise TypeError("krylov_matrix: M is not a matrix") @@ -19211,53 +19220,60 @@ cdef class Matrix(Matrix1): raise ValueError("krylov_matrix: matrix M does not have same base ring as E.") if degree is None: - degree = self.ncols() - if not isinstance(degree, (int, sage.rings.integer.Integer)): - raise TypeError("krylov_matrix: degree is not an integer.") - if degree < 0: - raise ValueError("krylov_matrix: degree bound cannot be negative.") - - if shift is None or shift == 0: - shift = (ZZ**self.nrows()).zero() - if isinstance(shift, (list, tuple)): - shift = (ZZ**self.nrows())(shift) - if shift not in ZZ**self.nrows(): - raise ValueError(f"krylov_matrix: shift is not an integer vector of length {self.nrows()}.") + degree = vector(ZZ, [self.ncols()] * self.nrows()) + if isinstance(degree, (int, sage.rings.integer.Integer)): + degree = vector(ZZ, [degree] * self.nrows()) + if isinstance(degree, (list, tuple)): + degree = vector(ZZ, degree) + if degree not in ZZ**self.nrows(): + raise TypeError(f"krylov_matrix: degree must be an None, an integer, or an integer vector of length self.nrows().") + if self.nrows() > 0 and min(degree) < 0: + raise ValueError(f"krylov_matrix: degree must not contain a negative bound.") + + if shifts is None or shifts == 0: + shifts = (ZZ**self.nrows()).zero() + if isinstance(shifts, (list, tuple)): + shifts = (ZZ**self.nrows())(shifts) + if shifts not in ZZ**self.nrows(): + raise ValueError(f"krylov_matrix: shifts is not an integer vector of length {self.nrows()}.") m = self.nrows() # define priority and indexing functions # index: [0,1,...,m-1,m,...,2m-1,...,m*degree-1] -> [(0,0),(1,0),...,(m-1,0),(0,1),...,(m-1,1),...,(m-1,degree)] - priority = lambda c,d : shift[c] + d + priority = lambda c,d : shifts[c] + d index = lambda i : (i%m,i//m) # deduce permutation by sorting priorities # rows should be ordered by ascending priority, then by associated row number in E (i%m) - priority_triplets = sorted([[priority(*index(i)),index(i),i] for i in range(m*(degree+1))]) + priority_triplets = [[priority(*index(i)),index(i)] for i in range(m*(max(degree, default=0)+1)) if index(i)[1] <= degree[index(i)[0]]] + priority_triplets = sorted([[t[0],t[1],i] for i,t in enumerate(priority_triplets)]) priority_permutation = Permutation([t[2]+1 for t in priority_triplets]) # store 2 blocks for the Krylov matrix blocks = [self] + # keep track of indices + indices = [(i,0) for i in range(m)] + # for i from 0 to degree (inclusive), merge existing blocks, add E*M^i M_i = None - for i in range(math.ceil(math.log(degree + 1, 2))): + for i in range(math.ceil(math.log(max(degree, default=0) + 1, 2))): if M_i is None: M_i = M else: M_i = M_i * M_i - blocks = [matrix.block([[blocks[0]],[blocks[1]]],subdivide=False)] - blocks.append(blocks[0] * M_i) + blocks = [blocks[0].stack(blocks[1],subdivide=False)] + new_indices = [(row[0], row[1] + 2**i) for row in indices] + blocks.append((blocks[0].matrix_from_rows([i for i,x in enumerate(new_indices) if x[1] <= degree[x[0]]]) * M_i)) + indices.extend([x for x in new_indices if x[1] <= degree[x[0]]]) - if (degree + 1) & degree != 0: - blocks[1] = blocks[1].matrix_from_rows(range((degree + 1) * m - blocks[1].nrows())) - - # return block matrix permuted according to shift - krylov = matrix.block([[blocks[0]],[blocks[1]]],subdivide=False) + # return block matrix permuted according to shifts + krylov = blocks[0].stack(blocks[1],subdivide=False) if len(blocks) > 1 else blocks[0] krylov.permute_rows(priority_permutation) return krylov - def _naive_krylov_basis(self, M, shift=None, degree=None, output_pairs=True): + def _naive_krylov_basis(self, M, shifts, degree, output_row_indices): r""" Compute the rank profile (row and column) of the Krylov matrix built from ``self`` and matrix ``M``. @@ -19265,16 +19281,19 @@ cdef class Matrix(Matrix1): INPUT: - ``M`` -- square matrix used in the Krylov construction. - - ``degree`` -- maximum power of ``M`` in Krylov matrix. If None, a suitable upper bound of ``self.ncols()`` is default. - - ``shift`` -- list of ``self.nrows()`` integers (optional): priority row shift. If ``None``, defaults to all zeroes. + - ``degree`` -- maximum power of ``M`` in Krylov matrix. If None, a + suitable upper bound of ``self.ncols()`` is default. + - ``shifts`` -- list of ``self.nrows()`` integers (optional): priority + row shifts. If ``None``, defaults to all zeroes. OUTPUT: - - A tuple (row_profile, pivot_matrix, column_profile): - * row_profile: list of the first r row indices in the Krylov matrix ``K`` corresponding to independent rows - * pivot: the submatrix of ``K`` given by those rows - * column_profile: list of the first r independent column indices in ``pivot_matrix`` - where r is the rank of ``K``. + A tuple (row_profile, pivot_matrix, column_profile): + - pivot: the submatrix of ``K`` given by those rows + - row_profile: list of the first r row indices in the Krylov matrix + ``K`` corresponding to independent rows + - column_profile: list of the first r independent column indices in + ``pivot_matrix`` where r is the rank of ``K``. TESTS:: @@ -19304,17 +19323,17 @@ cdef class Matrix(Matrix1): [ 0 0 0] [ 0 0 0] [ 0 0 0] - sage: E.krylov_basis(M, output_pairs=False) + sage: E.krylov_basis(M, output_row_indices=True) ( - [27 49 29] - [50 58 0] - (0, 1, 3), [ 0 27 49], (0, 1, 2) + [27 49 29] + [50 58 0] + [ 0 27 49], (0, 1, 3), (0, 1, 2) ) sage: E.krylov_basis(M) ( - [27 49 29] - [50 58 0] - ((0, 0), (1, 0), (0, 1)), [ 0 27 49], (0, 1, 2) + [27 49 29] + [50 58 0] + [ 0 27 49], ((0, 0), (1, 0), (0, 1)), (0, 1, 2) ) sage: rows = [(i+j*E.nrows()+1,i,j) for j in range(degree+1) for i in range(E.nrows())] @@ -19331,8 +19350,8 @@ cdef class Matrix(Matrix1): (10, 0, 3), (11, 1, 3), (12, 2, 3)] - sage: shift = [0,3,6] - sage: rows.sort(key=lambda x: (x[2] + shift[x[1]],x[1])) + sage: shifts = [0,3,6] + sage: rows.sort(key=lambda x: (x[2] + shifts[x[1]],x[1])) sage: rows [(1, 0, 0), (4, 0, 1), @@ -19346,7 +19365,7 @@ cdef class Matrix(Matrix1): (6, 2, 1), (9, 2, 2), (12, 2, 3)] - sage: E.krylov_matrix(M,shift=shift) + sage: E.krylov_matrix(M,shifts=shifts) [27 49 29] [ 0 27 49] [ 0 0 27] @@ -19359,21 +19378,21 @@ cdef class Matrix(Matrix1): [ 0 77 10] [ 0 0 77] [ 0 0 0] - sage: E.krylov_basis(M, shift=shift, output_pairs=False) + sage: E.krylov_basis(M, shifts=shifts, output_row_indices=True) ( - [27 49 29] - [ 0 27 49] - (0, 1, 2), [ 0 0 27], (0, 1, 2) + [27 49 29] + [ 0 27 49] + [ 0 0 27], (0, 1, 2), (0, 1, 2) ) - sage: E.krylov_basis(M,shift=shift) + sage: E.krylov_basis(M,shifts=shifts) ( - [27 49 29] - [ 0 27 49] - ((0, 0), (0, 1), (0, 2)), [ 0 0 27], (0, 1, 2) + [27 49 29] + [ 0 27 49] + [ 0 0 27], ((0, 0), (0, 1), (0, 2)), (0, 1, 2) ) - sage: shift = [3,0,2] - sage: rows.sort(key=lambda x: (x[2] + shift[x[1]],x[1])) + sage: shifts = [3,0,2] + sage: rows.sort(key=lambda x: (x[2] + shifts[x[1]],x[1])) sage: rows [(2, 1, 0), (5, 1, 1), @@ -19387,7 +19406,7 @@ cdef class Matrix(Matrix1): (7, 0, 2), (12, 2, 3), (10, 0, 3)] - sage: E.krylov_matrix(M,shift=shift) + sage: E.krylov_matrix(M,shifts=shifts) [50 58 0] [ 0 50 58] [ 0 0 50] @@ -19400,45 +19419,24 @@ cdef class Matrix(Matrix1): [ 0 0 27] [ 0 0 0] [ 0 0 0] - sage: E.krylov_basis(M,shift=shift,output_pairs=False) + sage: E.krylov_basis(M,shifts=shifts,output_row_indices=True) ( - [50 58 0] - [ 0 50 58] - (0, 1, 2), [ 0 0 50], (0, 1, 2) + [50 58 0] + [ 0 50 58] + [ 0 0 50], (0, 1, 2), (0, 1, 2) ) - sage: E.krylov_basis(M,shift=shift) + sage: E.krylov_basis(M,shifts=shifts) ( - [50 58 0] - [ 0 50 58] - ((1, 0), (1, 1), (1, 2)), [ 0 0 50], (0, 1, 2) + [50 58 0] + [ 0 50 58] + [ 0 0 50], ((1, 0), (1, 1), (1, 2)), (0, 1, 2) ) """ from sage.matrix.constructor import matrix - if not isinstance(M, Matrix): - raise TypeError("naive_krylov_basis: M is not a matrix") - if M.nrows() != self.ncols() or M.ncols() != self.ncols(): - raise ValueError("naive_krylov_basis: matrix M does not have correct dimensions.") - if M.base_ring() != self.base_ring(): - raise ValueError("naive_krylov_basis: matrix M does not have same base ring as E.") - - if degree is None: - degree = self.ncols() - if not isinstance(degree, (int, sage.rings.integer.Integer)): - raise TypeError("naive_krylov_basis: degree is not an integer.") - if degree < 0: - raise ValueError("naive_krylov_basis: degree bound cannot be negative.") - - if shift is None or shift == 0: - shift = (ZZ**self.nrows()).zero() - if isinstance(shift, (list, tuple)): - shift = (ZZ**self.nrows())(shift) - if shift not in ZZ**self.nrows(): - raise ValueError(f"naive_krylov_basis: shift is not an integer vector of length {self.nrows()}.") - m = self.nrows() - K = self.krylov_matrix(M,shift,degree) + K = self.krylov_matrix(M,shifts,degree) # calculate first independent rows of K row_profile = K.pivot_rows() @@ -19449,22 +19447,22 @@ cdef class Matrix(Matrix1): # calculate first independent columns of pivot col_profile = pivot.pivots() - if output_pairs: + if not output_row_indices: # define priority and indexing functions # index: [0,1,...,m-1,m,...,2m-1,...,m*degree-1] -> [(0,0),(1,0),...,(m-1,0),(0,1),...,(m-1,1),...,(m-1,degree)] - priority = lambda c,d : shift[c] + d + priority = lambda c,d : shifts[c] + d index = lambda i : (i%m,i//m) # deduce permutation by sorting priorities # rows should be ordered by ascending priority, then by associated row number in E (i%m) - priority_pairs = sorted([[priority(*index(i)),index(i)] for i in range(m*(degree+1))]) + priority_pairs = sorted([[priority(*index(i)),index(i)] for i in range(m*(max(degree, default=0)+1)) if index(i)[1] <= degree[index(i)[0]]]) row_profile = tuple([priority_pairs[i][1] for i in row_profile]) - return row_profile,pivot,col_profile + return pivot, row_profile, col_profile - def _elimination_krylov_basis(self, M, shift=None, degree=None, output_pairs=True): + def _elimination_krylov_basis(self, M, shifts, degree, output_row_indices): r""" Compute the rank profile (row and column) of the block Krylov matrix built from ``self`` and matrix ``M``. @@ -19472,19 +19470,21 @@ cdef class Matrix(Matrix1): INPUT: - ``M`` -- square matrix used in the Krylov construction. - - ``degree`` -- maximum power of ``M`` in Krylov matrix. If None, a suitable upper bound of ``self.ncols()`` is default. - - ``shift`` -- list of ``self.nrows()`` integers (optional): priority - row shift. If ``None``, defaults to all zeroes. + - ``degree`` -- maximum power of ``M`` in Krylov matrix. If None, a + suitable upper bound of ``self.ncols()`` is default. + - ``shifts`` -- list of ``self.nrows()`` integers (optional): priority + row shifts. If ``None``, defaults to all zeroes. OUTPUT: A tuple (row_profile, pivot_matrix, column_profile): + - ``pivot_matrix``: the submatrix of ``K`` given by ``row_profile`` - ``row_profile``: list of the first r row indices in the striped Krylov matrix ``K`` corresponding to independent rows. If - ``output_pairs`` is True, then the output is the list of pairs + ``output_row_indices`` is False, then the output is the list of pairs (c_i, d_i) where c_i is the corresponding row in the original matrix - and d_i the corresponding power of ``M``. - - ``pivot_matrix``: the submatrix of ``K`` given by ``row_profile`` + and d_i the corresponding power of ``M``. If True, the output is the + row indices in the Krylov matrix. - ``column_profile``: the first r independent column indices in ``pivot_matrix`` where r is the rank of `K`. @@ -19516,17 +19516,17 @@ cdef class Matrix(Matrix1): [ 0 0 0] [ 0 0 0] [ 0 0 0] - sage: E.krylov_basis(M,output_pairs=False) + sage: E.krylov_basis(M,output_row_indices=True) ( - [27 49 29] - [50 58 0] - (0, 1, 3), [ 0 27 49], (0, 1, 2) + [27 49 29] + [50 58 0] + [ 0 27 49], (0, 1, 3), (0, 1, 2) ) sage: E.krylov_basis(M) ( - [27 49 29] - [50 58 0] - ((0, 0), (1, 0), (0, 1)), [ 0 27 49], (0, 1, 2) + [27 49 29] + [50 58 0] + [ 0 27 49], ((0, 0), (1, 0), (0, 1)), (0, 1, 2) ) sage: rows = [(i+j*E.nrows()+1,i,j) for j in range(degree+1) for i in range(E.nrows())] sage: rows @@ -19542,8 +19542,8 @@ cdef class Matrix(Matrix1): (10, 0, 3), (11, 1, 3), (12, 2, 3)] - sage: shift = [0,3,6] - sage: rows.sort(key=lambda x: (x[2] + shift[x[1]],x[1])) + sage: shifts = [0,3,6] + sage: rows.sort(key=lambda x: (x[2] + shifts[x[1]],x[1])) sage: rows [(1, 0, 0), (4, 0, 1), @@ -19557,7 +19557,7 @@ cdef class Matrix(Matrix1): (6, 2, 1), (9, 2, 2), (12, 2, 3)] - sage: E.krylov_matrix(M,shift=shift) + sage: E.krylov_matrix(M,shifts=shifts) [27 49 29] [ 0 27 49] [ 0 0 27] @@ -19570,20 +19570,20 @@ cdef class Matrix(Matrix1): [ 0 77 10] [ 0 0 77] [ 0 0 0] - sage: E.krylov_basis(M,shift=shift,output_pairs=False) + sage: E.krylov_basis(M,shifts=shifts,output_row_indices=True) ( - [27 49 29] - [ 0 27 49] - (0, 1, 2), [ 0 0 27], (0, 1, 2) + [27 49 29] + [ 0 27 49] + [ 0 0 27], (0, 1, 2), (0, 1, 2) ) - sage: E.krylov_basis(M,shift=shift) + sage: E.krylov_basis(M,shifts=shifts) ( - [27 49 29] - [ 0 27 49] - ((0, 0), (0, 1), (0, 2)), [ 0 0 27], (0, 1, 2) + [27 49 29] + [ 0 27 49] + [ 0 0 27], ((0, 0), (0, 1), (0, 2)), (0, 1, 2) ) - sage: shift = [3,0,2] - sage: rows.sort(key=lambda x: (x[2] + shift[x[1]],x[1])) + sage: shifts = [3,0,2] + sage: rows.sort(key=lambda x: (x[2] + shifts[x[1]],x[1])) sage: rows [(2, 1, 0), (5, 1, 1), @@ -19597,7 +19597,7 @@ cdef class Matrix(Matrix1): (7, 0, 2), (12, 2, 3), (10, 0, 3)] - sage: E.krylov_matrix(M,shift=shift) + sage: E.krylov_matrix(M,shifts=shifts) [50 58 0] [ 0 50 58] [ 0 0 50] @@ -19610,52 +19610,31 @@ cdef class Matrix(Matrix1): [ 0 0 27] [ 0 0 0] [ 0 0 0] - sage: E.krylov_basis(M,shift=shift,output_pairs=False) + sage: E.krylov_basis(M,shifts=shifts,output_row_indices=True) ( - [50 58 0] - [ 0 50 58] - (0, 1, 2), [ 0 0 50], (0, 1, 2) + [50 58 0] + [ 0 50 58] + [ 0 0 50], (0, 1, 2), (0, 1, 2) ) - sage: E.krylov_basis(M,shift=shift) + sage: E.krylov_basis(M,shifts=shifts) ( - [50 58 0] - [ 0 50 58] - ((1, 0), (1, 1), (1, 2)), [ 0 0 50], (0, 1, 2) + [50 58 0] + [ 0 50 58] + [ 0 0 50], ((1, 0), (1, 1), (1, 2)), (0, 1, 2) ) """ from sage.matrix.constructor import matrix import math from sage.combinat.permutation import Permutation - if not isinstance(M, Matrix): - raise TypeError("krylov_basis: M is not a matrix") - if M.nrows() != self.ncols() or M.ncols() != self.ncols(): - raise ValueError("krylov_basis: matrix M does not have correct dimensions.") - if M.base_ring() != self.base_ring(): - raise ValueError("krylov_basis: matrix M does not have same base ring as E.") - - if degree is None: - degree = self.ncols() - if not isinstance(degree, (int, sage.rings.integer.Integer)): - raise TypeError("krylov_basis: degree is not an integer.") - if degree < 0: - raise ValueError("krylov_basis: degree bound cannot be negative.") - - if shift is None or shift == 0: - shift = (ZZ**self.nrows()).zero() - if isinstance(shift, (list, tuple)): - shift = (ZZ**self.nrows())(shift) - if shift not in ZZ**self.nrows(): - raise ValueError(f"krylov_basis: shift is not a Z-vector of length {self.nrows()}.") - m = self.nrows() sigma = self.ncols() if m == 0: - return (),self,() + return self, (), () - # calculate row profile of self, with shift applied - self_permutation = Permutation(sorted([i+1 for i in range(m)],key=lambda x:(shift[x-1],x-1))) + # calculate row profile of self, with shifts applied + self_permutation = Permutation(sorted([i+1 for i in range(m)],key=lambda x:(shifts[x-1],x-1))) self_permuted = self.with_permuted_rows(self_permutation) row_profile_self_permuted = self_permuted.pivot_rows() @@ -19668,24 +19647,24 @@ cdef class Matrix(Matrix1): r = len(row_profile_self) if r == 0: - return (),self.matrix_from_rows([]),() + return self.matrix_from_rows([]),(),() R = self_permuted.matrix_from_rows(row_profile_self_permuted) M_L = None - for l in range(math.ceil(math.log(degree + 1, 2))): + for l in range(math.ceil(math.log(max(degree, default=0) + 1, 2))): L = pow(2,l) # adding 2^l to each degree - row_extension = [(x[0],x[1] + L) for x in row_profile_self if x[1] + L <= degree] + row_extension = [(x[0],x[1] + L) for x in row_profile_self if x[1] + L <= degree[x[0]]] if len(row_extension) == 0: break # concatenate two sequences (order preserved) k = row_profile_exhausted + row_profile_self + row_extension - # calculate sorting permutation, sort k by (shift[c]+d, c) - k_perm = Permutation(sorted([i+1 for i in range(len(k))],key=lambda x: (shift[k[x-1][0]] + k[x-1][1],k[x-1][0]))) + # calculate sorting permutation, sort k by (shifts[c]+d, c) + k_perm = Permutation(sorted([i+1 for i in range(len(k))],key=lambda x: (shifts[k[x-1][0]] + k[x-1][1],k[x-1][0]))) # fast calculation of rows formed by indices in k if M_L is None: @@ -19693,7 +19672,7 @@ cdef class Matrix(Matrix1): else: M_L = M_L * M_L - R = matrix.block([[exhausted],[R],[R*M_L]],subdivide=False) + R = matrix.block([[exhausted],[R],[(R*M_L).matrix_from_rows([i for i in range(len(row_profile_self)) if row_profile_self[i][1] + L <= degree[row_profile_self[i][0]]])]], subdivide=False) # sort rows of R, find profile, translate to k (indices of full krylov matrix) R.permute_rows(k_perm) @@ -19701,7 +19680,7 @@ cdef class Matrix(Matrix1): row_profile_R = R.pivot_rows() r = len(row_profile_R) - if r == sigma and l < math.ceil(math.log(degree,2)): + if r == sigma and l < math.ceil(math.log(max(degree, default=0) + 1, 2)) - 1: tail = list(range(row_profile_R[-1]+1,R.nrows())) excluded_rows.update(set([k[k_perm(i+1)-1][0] for i in tail if i < len(k)])) @@ -19714,7 +19693,7 @@ cdef class Matrix(Matrix1): exhausted = R.matrix_from_rows(xmi) R = R.matrix_from_rows(imi) else: - if l == math.ceil(math.log(degree,2)): + if l == math.ceil(math.log(max(degree, default=0) + 1, 2)) - 1: row_profile_exhausted = [] exhausted = matrix.zero(self.base_ring(), 0, sigma) row_profile_self = [k[k_perm(i+1)-1] for i in row_profile_R] @@ -19726,22 +19705,22 @@ cdef class Matrix(Matrix1): elif exhausted.nrows() != 0: k = row_profile_exhausted + row_profile_self R = exhausted.stack(R) - k_perm = Permutation(sorted([i+1 for i in range(len(k))],key=lambda x: (shift[k[x-1][0]] + k[x-1][1],k[x-1][0]))) + k_perm = Permutation(sorted([i+1 for i in range(len(k))],key=lambda x: (shifts[k[x-1][0]] + k[x-1][1],k[x-1][0]))) R.permute_rows(k_perm) row_profile_self = [k[k_perm(i+1)-1] for i in range(len(k))] col_profile = R.pivots() - if output_pairs: - return tuple(row_profile_self), R, col_profile + if not output_row_indices: + return R, tuple(row_profile_self), col_profile # convert c,d to actual position in striped Krylov matrix - phi = lambda c,d : sum(min(max(shift[c] - shift[i] + d + (i < c and shift[i] <= shift[c] + degree),0),degree+1) for i in range(m)) + phi = lambda c,d : sum(min(max(shifts[c] - shifts[i] + d + (i < c and shifts[i] <= shifts[c] + max(degree)),0),max(degree)+1) for i in range(m)) # TODO: fix row_profile_K = [phi(*row) for row in row_profile_self] - return tuple(row_profile_K), R, col_profile + return R, tuple(row_profile_K), col_profile - def krylov_basis(self, M, shift=None, degree=None, output_pairs=True, algorithm='default'): + def krylov_basis(self, M, shifts=None, degree=None, output_row_indices=False, algorithm=None): r""" Compute the rank profile (row and column) of the block Krylov matrix built from ``self`` and matrix ``M``. @@ -19749,23 +19728,57 @@ cdef class Matrix(Matrix1): INPUT: - ``M`` -- square matrix used in the Krylov construction. - - ``degree`` -- maximum power of ``M`` in Krylov matrix. If None, a suitable upper bound of ``self.ncols()`` is default. - - ``shift`` -- list of ``self.nrows()`` integers (optional): priority - row shift. If ``None``, defaults to all zeroes. + - ``degree`` -- maximum power of ``M`` in Krylov matrix. If None, a + suitable upper bound of ``self.ncols()`` is default. + - ``shifts`` -- list of ``self.nrows()`` integers (optional): priority + row shifts. If ``None``, defaults to all zeroes. + - ``output_row_indices`` -- determines whether ``row_profile`` is + returned as pairs of row indices in ``self`` and degrees of ``M`` + (``False``), or row indices in the Krylov matrix (``True``). + - ``algorithm`` -- either 'naive', 'elimination', or None (let + Sage choose). OUTPUT: A tuple (row_profile, pivot_matrix, column_profile): + - ``pivot_matrix``: the submatrix of ``K`` given by ``row_profile`` - ``row_profile``: list of the first r row indices in the striped Krylov matrix ``K`` corresponding to independent rows. If - ``output_pairs`` is True, then the output is the list of pairs + ``output_row_indices`` is False, then the output is the list of pairs (c_i, d_i) where c_i is the corresponding row in the original matrix - and d_i the corresponding power of ``M``. - - ``pivot_matrix``: the submatrix of ``K`` given by ``row_profile`` + and d_i the corresponding power of ``M``. If True, the output is the + row indices in the Krylov matrix. - ``column_profile``: the first r independent column indices in ``pivot_matrix`` where r is the rank of `K`. """ - if algorithm == 'default': + from sage.modules.free_module_element import vector + + if not isinstance(M, Matrix): + raise TypeError("_naive_krylov_basis: M is not a matrix") + if M.nrows() != self.ncols() or M.ncols() != self.ncols(): + raise ValueError("_naive_krylov_basis: matrix M does not have correct dimensions.") + if M.base_ring() != self.base_ring(): + raise ValueError("naive_krylov_basis: matrix M does not have same base ring as E.") + + if degree is None: + degree = vector(ZZ, [self.ncols()] * self.nrows()) + if isinstance(degree, (int, sage.rings.integer.Integer)): + degree = vector(ZZ, [degree] * self.nrows()) + if isinstance(degree, (list, tuple)): + degree = vector(ZZ, degree) + if degree not in ZZ**self.nrows(): + raise TypeError(f"_naive_krylov_basis: degree must be an None, an integer, or an integer vector of length self.nrows().") + if self.nrows() > 0 and min(degree) < 0: + raise ValueError(f"_naive_krylov_basis: degree must not contain a negative bound.") + + if shifts is None or shifts == 0: + shifts = (ZZ**self.nrows()).zero() + if isinstance(shifts, (list, tuple)): + shifts = (ZZ**self.nrows())(shifts) + if shifts not in ZZ**self.nrows(): + raise ValueError(f"_naive_krylov_basis: shifts is not an integer vector of length {self.nrows()}.") + + if algorithm is None: if self.base_ring().order() == 2: if self.nrows() <= 12: algorithm = 'naive' @@ -19782,38 +19795,38 @@ cdef class Matrix(Matrix1): else: algorithm = 'elimination' if algorithm == 'naive': - return self._naive_krylov_basis(M, shift, degree, output_pairs) + return self._naive_krylov_basis(M, shifts, degree, output_row_indices) elif algorithm == 'elimination': - return self._elimination_krylov_basis(M, shift, degree, output_pairs) + return self._elimination_krylov_basis(M, shifts, degree, output_row_indices) else: - raise ValueError("algorithm must be one of \"default\", \"naive\" or \"elimination\".") + raise ValueError("algorithm must be one of None, \"naive\" or \"elimination\".") - def krylov_kernel_basis(self, M, var_name='x', degree=None, shift=None, output_coefficients=False): + def krylov_kernel_basis(self, M, degree=None, shifts=None, output_coefficients=False, var='x'): r""" - Return a shifted minimal linear interpolation basis for - (``self``,``M``) in ``s``-Popov form with respect to a shift ``s``. + Return a shifted minimal krylov kernel basis for (``self``,``M``) in + ``s``-Popov form with respect to a set of shifts ``s``. - Given a K[x]-module V defined by multiplication matrix ``M``, i.e. elements - are vectors of length n in K, scalars are polynomials in K[x], and - scalar multiplication is defined by p v := v p(``M``), ``self`` represents - ``self.nrows()`` elements e_1, ..., e_m in V: + Given a K[x]-module V defined by multiplication matrix ``M``, i.e. + elements are vectors of length n in K, scalars are polynomials in K[x], + and scalar multiplication is defined by p v := v p(``M``), ``self`` + represents ``self.nrows()`` elements e_1, ..., e_m in V: A linear interpolation basis is a set of tuples of length ``self.nrows()`` in K[x] that form a basis for solutions to the equation ``p_1 e_1 + ... + p_n e_n = 0``. - See [Kai1980] + See [Kai1980]_ INPUT: - ``M`` -- square matrix for Krylov construction. - - ``var_name`` -- variable name for the returned polynomial matrix + - ``shifts`` -- list of self.nrows() integers (optional): controls + priority of rows. If ``None``, defaults to all zeroes. - ``degree`` -- upper bound on degree of minpoly(`M`). If None, a suitable upper bound of ``self.ncols()`` is default. - - ``shift`` -- list of self.nrows() integers (optional): controls - priority of rows. If ``None``, defaults to all zeroes. - ``output_coefficients`` -- If True, outputs the matrix as a sparse block matrix of coefficients. + - ``var`` -- variable name for the returned polynomial matrix OUTPUT: @@ -19834,7 +19847,7 @@ cdef class Matrix(Matrix1): [0 0 1] [0 0 0] sage: degree = E.ncols() - sage: krylov_matrix = matrix.block([[E*M^i] for i in range(degree+1)],subdivide=False) + sage: krylov_matrix = E.krylov_matrix(M) sage: krylov_matrix [27 49 29] [50 58 0] @@ -19848,12 +19861,12 @@ cdef class Matrix(Matrix1): [ 0 0 0] [ 0 0 0] [ 0 0 0] - sage: basis = E.krylov_kernel_basis(M,'x') + sage: basis = E.krylov_kernel_basis(M) sage: basis [x^2 + 40*x + 82 76 0] [ 3*x + 13 x + 57 0] [ 96 96 1] - sage: basis_coeff = E.krylov_kernel_basis(M,'x',output_coefficients=True) + sage: basis_coeff = E.krylov_kernel_basis(M,output_coefficients=True) sage: basis_coeff [82 76 0 40 0 0 1 0 0] [13 57 0 3 1 0 0 0 0] @@ -19864,47 +19877,47 @@ cdef class Matrix(Matrix1): True sage: basis_coeff.augment(matrix.zero(R,E.nrows())) * krylov_matrix == matrix.zero(R,E.nrows()) True - sage: len(E.krylov_basis(M)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + sage: len(E.krylov_basis(M)[1]) == sum(basis[i,i].degree() for i in range(E.nrows())) True - sage: shift = [0,3,6] - sage: basis = E.krylov_kernel_basis(M,'x',shift=shift) + sage: shifts = [0,3,6] + sage: basis = E.krylov_kernel_basis(M,shifts=shifts) sage: basis [ x^3 0 0] [60*x^2 + 72*x + 70 1 0] [60*x^2 + 72*x + 69 0 1] - sage: basis_coeff = E.krylov_kernel_basis(M,'x',shift=shift,output_coefficients=True) + sage: basis_coeff = E.krylov_kernel_basis(M,shifts=shifts,output_coefficients=True) sage: basis_coeff [ 0 0 0 0 0 0 0 0 0 1 0 0] [70 1 0 72 0 0 60 0 0 0 0 0] [69 0 1 72 0 0 60 0 0 0 0 0] - sage: basis.is_popov(shifts=shift) + sage: basis.is_popov(shifts=shifts) True sage: basis.degree() <= E.ncols() True sage: basis_coeff * krylov_matrix == matrix.zero(R,E.nrows()) True - sage: len(E.krylov_basis(M,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + sage: len(E.krylov_basis(M,shifts=shifts)[1]) == sum(basis[i,i].degree() for i in range(E.nrows())) True - sage: shift = [3,0,2] - sage: basis = E.krylov_kernel_basis(M,'x',shift=shift) + sage: shifts = [3,0,2] + sage: basis = E.krylov_kernel_basis(M,shifts=shifts) sage: basis [ 1 26*x^2 + 49*x + 79 0] [ 0 x^3 0] [ 0 26*x^2 + 49*x + 78 1] - sage: basis_coeff = E.krylov_kernel_basis(M,'x',shift=shift,output_coefficients=True) + sage: basis_coeff = E.krylov_kernel_basis(M,shifts=shifts,output_coefficients=True) sage: basis_coeff [ 1 79 0 0 49 0 0 26 0 0 0 0] [ 0 0 0 0 0 0 0 0 0 0 1 0] [ 0 78 1 0 49 0 0 26 0 0 0 0] - sage: basis.is_popov(shifts=shift) + sage: basis.is_popov(shifts=shifts) True sage: basis.degree() <= E.ncols() True sage: basis_coeff * krylov_matrix == matrix.zero(R,E.nrows()) True - sage: len(E.krylov_basis(M,shift=shift)[0]) == sum(basis[i,i].degree() for i in range(E.nrows())) + sage: len(E.krylov_basis(M,shifts=shifts)[1]) == sum(basis[i,i].degree() for i in range(E.nrows())) True """ from sage.combinat.permutation import Permutation @@ -19919,30 +19932,30 @@ cdef class Matrix(Matrix1): if M.base_ring() != self.base_ring(): raise ValueError("krylov_kernel_basis: matrix M does not have same base ring as E.") - if not isinstance(var_name, str): - raise TypeError("var_name is not a string") + if not isinstance(var, str): + raise TypeError("var is not a string") if degree is None: degree = self.ncols() if not isinstance(degree, (int, sage.rings.integer.Integer)): raise TypeError("krylov_kernel_basis: degree is not an integer.") - if degree < 0: + if self.nrows() > 0 and degree < 0: raise ValueError("krylov_kernel_basis: degree bound cannot be negative.") - if shift is None or shift == 0: - shift = (ZZ**self.nrows()).zero() - if isinstance(shift, (list, tuple)): - shift = (ZZ**self.nrows())(shift) - if shift not in ZZ**self.nrows(): - raise ValueError(f"krylov_kernel_basis: shift is not a Z-vector of length {self.nrows()}.") + if shifts is None or shifts == 0: + shifts = (ZZ**self.nrows()).zero() + if isinstance(shifts, (list, tuple)): + shifts = (ZZ**self.nrows())(shifts) + if shifts not in ZZ**self.nrows(): + raise ValueError(f"krylov_kernel_basis: shifts is not a Z-vector of length {self.nrows()}.") m = self.nrows() sigma = self.ncols() - poly_ring = PolynomialRing(self.base_ring(),var_name) + poly_ring = PolynomialRing(self.base_ring(),var) # calculate krylov profile - row_profile, pivot, col_profile = self.krylov_basis(M, shift, degree) + pivot, row_profile, col_profile = self.krylov_basis(M, shifts, degree) if len(row_profile) == 0: return matrix.identity(poly_ring,m) @@ -19977,7 +19990,7 @@ cdef class Matrix(Matrix1): relation = D*C.inverse() if output_coefficients: - coefficients = matrix.zero(self.base_ring(),m,m*(max(degree_c)+1)) + coefficients = matrix.zero(self.base_ring(),m,m*(max(degree_c, default=0)+1), sparse=True) for i in range(m): coefficients[i,i+degree_c[i]*m] = self.base_ring().one() for col in range(relation.ncols()): From 4a1084ec69b492453d59551f564deb583e988119 Mon Sep 17 00:00:00 2001 From: Nicholas Bell Date: Wed, 13 Aug 2025 17:33:20 +0200 Subject: [PATCH 28/28] finish multiple degree bound modifications --- src/sage/matrix/matrix2.pyx | 86 ++++++++++++++++++++++++++----------- 1 file changed, 60 insertions(+), 26 deletions(-) diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx index 3efe21af167..7ce08a27b30 100644 --- a/src/sage/matrix/matrix2.pyx +++ b/src/sage/matrix/matrix2.pyx @@ -19145,9 +19145,11 @@ cdef class Matrix(Matrix1): - ``M`` -- a square matrix of size equal to the number of columns of ``self``. - - ``degree`` -- integer, the maximum exponent for the Krylov matrix. - ``shifts`` -- list of ``self.nrows()`` integers (optional), row priority shifts. If ``None``, defaults to all zeroes. + - ``degree`` -- the maximum power or list of ``self.nrows()`` maximum + powers of ``M`` for each row in ``self`` in the output. If None, + ``self.ncols()`` is chosen by default for all rows. OUTPUT: @@ -19270,21 +19272,30 @@ cdef class Matrix(Matrix1): # return block matrix permuted according to shifts krylov = blocks[0].stack(blocks[1],subdivide=False) if len(blocks) > 1 else blocks[0] - krylov.permute_rows(priority_permutation) - return krylov + return krylov.with_permuted_rows(priority_permutation) def _naive_krylov_basis(self, M, shifts, degree, output_row_indices): r""" Compute the rank profile (row and column) of the Krylov matrix built from ``self`` and matrix ``M``. + .. WARNING:: + + If the degree bounds are not true upper bounds, no guarantees are + made on the value of the output. + INPUT: - ``M`` -- square matrix used in the Krylov construction. - - ``degree`` -- maximum power of ``M`` in Krylov matrix. If None, a - suitable upper bound of ``self.ncols()`` is default. - ``shifts`` -- list of ``self.nrows()`` integers (optional): priority row shifts. If ``None``, defaults to all zeroes. + - ``degree`` -- an upper bound or list of ``self.nrows()`` upper bounds + on the power of ``M`` for each row in ``self`` in the + ``krylov_basis`` of ``M`` in Krylov matrix. If None, a suitable upper + bound of ``self.ncols()`` is default for all rows. + - ``output_row_indices`` -- determines whether the output row indices + are pairs of row indices in ``self`` and degrees of ``M`` (False) or + row indices in the Krylov matrix (True). OUTPUT: @@ -19467,13 +19478,23 @@ cdef class Matrix(Matrix1): Compute the rank profile (row and column) of the block Krylov matrix built from ``self`` and matrix ``M``. + .. WARNING:: + + If the degree bounds are not true upper bounds, no guarantees are + made on the value of the output. + INPUT: - ``M`` -- square matrix used in the Krylov construction. - - ``degree`` -- maximum power of ``M`` in Krylov matrix. If None, a - suitable upper bound of ``self.ncols()`` is default. - ``shifts`` -- list of ``self.nrows()`` integers (optional): priority row shifts. If ``None``, defaults to all zeroes. + - ``degree`` -- an upper bound or list of ``self.nrows()`` upper bounds + on the power of ``M`` for each row in ``self`` in the + ``krylov_basis`` of ``M`` in Krylov matrix. If None, a suitable upper + bound of ``self.ncols()`` is default for all rows. + - ``output_row_indices`` -- determines whether the output row indices + are pairs of row indices in ``self`` and degrees of ``M`` (False) or + row indices in the Krylov matrix (True). OUTPUT: @@ -19654,9 +19675,9 @@ cdef class Matrix(Matrix1): M_L = None for l in range(math.ceil(math.log(max(degree, default=0) + 1, 2))): - L = pow(2,l) + L = pow(2, l) # adding 2^l to each degree - row_extension = [(x[0],x[1] + L) for x in row_profile_self if x[1] + L <= degree[x[0]]] + row_extension = [(x[0], x[1] + L) for x in row_profile_self if x[1] + L <= degree[x[0]]] if len(row_extension) == 0: break @@ -19715,7 +19736,7 @@ cdef class Matrix(Matrix1): return R, tuple(row_profile_self), col_profile # convert c,d to actual position in striped Krylov matrix - phi = lambda c,d : sum(min(max(shifts[c] - shifts[i] + d + (i < c and shifts[i] <= shifts[c] + max(degree)),0),max(degree)+1) for i in range(m)) # TODO: fix + phi = lambda c,d : sum(min(max(shifts[c] - shifts[i] + d + (i < c and shifts[i] <= shifts[c] + d), 0), degree[i] + 1) for i in range(m)) row_profile_K = [phi(*row) for row in row_profile_self] return R, tuple(row_profile_K), col_profile @@ -19725,13 +19746,21 @@ cdef class Matrix(Matrix1): Compute the rank profile (row and column) of the block Krylov matrix built from ``self`` and matrix ``M``. + .. WARNING:: + + If the degree bounds are not true upper bounds, no guarantees are + made on the value of the output, including whether it is consistent + between different algorithms. + INPUT: - ``M`` -- square matrix used in the Krylov construction. - - ``degree`` -- maximum power of ``M`` in Krylov matrix. If None, a - suitable upper bound of ``self.ncols()`` is default. - ``shifts`` -- list of ``self.nrows()`` integers (optional): priority row shifts. If ``None``, defaults to all zeroes. + - ``degree`` -- an upper bound or list of ``self.nrows()`` upper bounds + on the power of ``M`` for each row in ``self`` in the + ``krylov_basis`` of ``M`` in Krylov matrix. If None, a suitable upper + bound of ``self.ncols()`` is default for all rows. - ``output_row_indices`` -- determines whether ``row_profile`` is returned as pairs of row indices in ``self`` and degrees of ``M`` (``False``), or row indices in the Krylov matrix (``True``). @@ -19754,11 +19783,11 @@ cdef class Matrix(Matrix1): from sage.modules.free_module_element import vector if not isinstance(M, Matrix): - raise TypeError("_naive_krylov_basis: M is not a matrix") + raise TypeError("krylov_basis: M is not a matrix") if M.nrows() != self.ncols() or M.ncols() != self.ncols(): - raise ValueError("_naive_krylov_basis: matrix M does not have correct dimensions.") + raise ValueError("krylov_basis: matrix M does not have correct dimensions.") if M.base_ring() != self.base_ring(): - raise ValueError("naive_krylov_basis: matrix M does not have same base ring as E.") + raise ValueError("krylov_basis: matrix M does not have same base ring as E.") if degree is None: degree = vector(ZZ, [self.ncols()] * self.nrows()) @@ -19767,16 +19796,16 @@ cdef class Matrix(Matrix1): if isinstance(degree, (list, tuple)): degree = vector(ZZ, degree) if degree not in ZZ**self.nrows(): - raise TypeError(f"_naive_krylov_basis: degree must be an None, an integer, or an integer vector of length self.nrows().") + raise TypeError(f"krylov_basis: degree must be an None, an integer, or an integer vector of length self.nrows().") if self.nrows() > 0 and min(degree) < 0: - raise ValueError(f"_naive_krylov_basis: degree must not contain a negative bound.") + raise ValueError(f"krylov_basis: degree must not contain a negative bound.") if shifts is None or shifts == 0: shifts = (ZZ**self.nrows()).zero() if isinstance(shifts, (list, tuple)): shifts = (ZZ**self.nrows())(shifts) if shifts not in ZZ**self.nrows(): - raise ValueError(f"_naive_krylov_basis: shifts is not an integer vector of length {self.nrows()}.") + raise ValueError(f"krylov_basis: shifts is not an integer vector of length {self.nrows()}.") if algorithm is None: if self.base_ring().order() == 2: @@ -19801,7 +19830,7 @@ cdef class Matrix(Matrix1): else: raise ValueError("algorithm must be one of None, \"naive\" or \"elimination\".") - def krylov_kernel_basis(self, M, degree=None, shifts=None, output_coefficients=False, var='x'): + def krylov_kernel_basis(self, M, shifts=None, degree=None, output_coefficients=False, var='x'): r""" Return a shifted minimal krylov kernel basis for (``self``,``M``) in ``s``-Popov form with respect to a set of shifts ``s``. @@ -19922,6 +19951,7 @@ cdef class Matrix(Matrix1): """ from sage.combinat.permutation import Permutation from sage.matrix.constructor import matrix + from sage.modules.free_module_element import vector from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing # INPUT VALIDATION @@ -19936,11 +19966,15 @@ cdef class Matrix(Matrix1): raise TypeError("var is not a string") if degree is None: - degree = self.ncols() - if not isinstance(degree, (int, sage.rings.integer.Integer)): - raise TypeError("krylov_kernel_basis: degree is not an integer.") - if self.nrows() > 0 and degree < 0: - raise ValueError("krylov_kernel_basis: degree bound cannot be negative.") + degree = vector(ZZ, [self.ncols()] * self.nrows()) + if isinstance(degree, (int, sage.rings.integer.Integer)): + degree = vector(ZZ, [degree] * self.nrows()) + if isinstance(degree, (list, tuple)): + degree = vector(ZZ, degree) + if degree not in ZZ**self.nrows(): + raise TypeError(f"krylov_kernel_basis: degree must be an None, an integer, or an integer vector of length self.nrows().") + if self.nrows() > 0 and min(degree) < 0: + raise ValueError(f"krylov_kernel_basis: degree must not contain a negative bound.") if shifts is None or shifts == 0: shifts = (ZZ**self.nrows()).zero() @@ -19990,7 +20024,7 @@ cdef class Matrix(Matrix1): relation = D*C.inverse() if output_coefficients: - coefficients = matrix.zero(self.base_ring(),m,m*(max(degree_c, default=0)+1), sparse=True) + coefficients = matrix.zero(self.base_ring(),m,m*(max(degree_c, default=0)+1)) for i in range(m): coefficients[i,i+degree_c[i]*m] = self.base_ring().one() for col in range(relation.ncols()): @@ -20037,7 +20071,7 @@ cdef class Matrix(Matrix1): deg, coeff = coeffs_map[row][col].popitem() if coeff not in monomial_cache: monomial_cache[coeff] = poly_ring(coeff) - basis_rows[row][col] = monomial_cache[coeff].shift(deg) + basis_rows[row][col] = monomial_cache[coeff].shift (deg) else: # general case mindeg = min(coeffs_map[row][col].keys())