Deterministic computation of the characteristic polynomial in the time of matrix multiplication
DDeterministic computation of the characteristic polynomialin the time of matrix multiplication
Vincent Neiger
Univ. Limoges, CNRS, XLIM, UMR 7252, F-87000 Limoges, France
Cl´ement Pernet
Universit´e Grenoble Alpes, Laboratoire Jean Kuntzmann, CNRS, UMR 5224700 avenue centrale, IMAG - CS 40700, 38058 Grenoble cedex 9, France
Abstract
This paper describes an algorithm which computes the characteristic polynomial of a matrix overa field within the same asymptotic complexity, up to constant factors, as the multiplication of twosquare matrices. Previously, to our knowledge, this was only achieved by resorting to genericityassumptions or randomization techniques, while the best known complexity bound with a generaldeterministic algorithm was obtained by Keller-Gehrig in 1985 and involves logarithmic factors.Our algorithm computes more generally the determinant of a univariate polynomial matrix in re-duced form, and relies on new subroutines for transforming shifted reduced matrices into shiftedweak Popov matrices, and shifted weak Popov matrices into shifted Popov matrices.
Keywords:
Characteristic polynomial, polynomial matrices, determinant, fast linear algebra.
1. Introduction
The last five decades witnessed a constant e ff ort towards computational reductions of linearalgebra problems to matrix multiplication. It has been showed that most classical problems arenot harder than multiplying two square matrices, such as matrix inversion, LU decomposition,nullspace basis computation, linear system solving, rank and determinant computations, etc. [7][25] [8, Chap. 16]. In this context, one major challenge stands out: designing a similar reductionto matrix multiplication for the computation of characteristic polynomials and related objectssuch as minimal polynomials and Frobenius forms. For the characteristic polynomial, significantprogress was achieved by Keller-Gehrig [31], and more recently by Pernet and Storjohann [39]who solved the problem if one allows randomization. This paper closes the problem by providinga deterministic algorithm with the same asymptotic complexity as matrix multiplication.The characteristic polynomial of a square matrix over a field K , say M ∈ K m × m , is defined asdet( x I m − M ). Specific algorithms exist for sparse or structured matrices; here we consider theclassical, dense case. Hereafter the complexity of an algorithm is measured as an upper boundon its arithmetic cost, that is, the number of basic field operations it uses to compute the output. Theorem 1.1.
Let K be a field. Using a subroutine which multiplies two matrices in K m × m inO ( m ω ) field operations for some ω > , the characteristic polynomial of a matrix in K m × m canbe computed deterministically in O ( m ω ) field operations. a r X i v : . [ c s . S C ] O c t utline. The rest of this introduction gives more details about our framework for complexitybounds (Section 1.1), summarizes previous work (Section 1.2), describes our contribution onpolynomial matrix determinant computation (Section 1.3), gives an overview of our approachand of new tools that we designed to avoid logarithmic factors (Sections 1.4 and 1.5), and finallylists a few perspectives (Section 1.6). Section 2 introduces the notation, main definitions, andbasic properties used in this paper. Then Section 3 presents the main algorithm of this paper alongwith a detailed complexity analysis. This algorithm uses two new technical tools described inSections 4 and 5: the transformation of reduced forms into weak Popov forms and of weak Popovforms into Popov forms, in the case of shifted forms.
In this paper, K is any field and we seek upper bounds on the complexity of algorithms whichoperate on objects such as matrices and polynomials over K . We consider the arithmetic cost ofthese algorithms, i.e. the number of basic operations in K that are used to compute the outputfrom some input of a given size. The basic operations are addition, subtraction, multiplication,and inversion in the field, as well as testing whether a given field element is zero.As already highlighted in Theorem 1.1, in this paper we fix any 2 < ω ≤ K m × m using O ( m ω ) operations in K : this algorithm isassumed to be the one used as a black box for all matrix multiplications arising in the algorithmswe design. The current best known cost bounds ensure that any ω > .
373 is feasible [33].In practice, one often considers a cubic algorithm with ω = ω = log (7) [46]. Our results hold with the only assumption that 2 < ω ≤ ω known tobe feasible at the time of writing, it is known that ω − ε is feasible as well for a su ffi ciently small ε >
0; then one might rather rely on this faster subroutine, and apply Keller-Gehrig’s algorithmto obtain the characteristic polynomial in O ( m ω − ε log( m )) operations in K , which is in O ( m ω ).Similarly, we consider a nondecreasing function d (cid:55)→ M ( d ) and an algorithm which multipliestwo polynomials in K [ x ] of degree at most d using at most M ( d ) operations in K ; our algorithmsrely on this subroutine for polynomial multiplication. Here d is any nonnegative real number; itwill often be a fraction D / m of positive integers; we assume that M ( d ) = ≤ d <
1, sothat M ( d ) ≥ d ≥
0. To help deriving complexity upper bounds, we will also consider thefollowing assumptions H sl , H sm , and H ω . H sl : 2 M ( d ) ≤ M (2 d ) for all d ≥ H sm : M ( d d ) ≤ M ( d ) M ( d ) for all d , d ≥ H ω : M ( d ) ∈ O ( d ω − − (cid:15) ) for some (cid:15) > M ( d ) ≥ d forall d ≥
1. The second and last assumptions are commonly made in complexity analyses fordivide and conquer algorithms on polynomial matrices [44, 23]: we refer to [44, Sec. 2] forfurther comments on these assumptions. They are satisfied by the cost bounds of polynomialmultiplication algorithms such as the quasi-linear algorithm of Cantor and Kaltofen [9] and, forsuitable fields K , the quasi-linear algorithm of Harvey and van der Hoeven and Lecerf [24],and most of Toom-Cook subquadratic algorithms [49, 10]. For the latter only H ω might not2e satisfied, depending on ω and on the number of points used. Note that with the currentestimates having ω > . / log(5) ≈ . < ω −
1; thus for such exponents ω all Toom-Cookalgorithms of order 5 or more satisfy all the above assumptions.Following [44, 23], we also define a function d (cid:55)→ M (cid:48) ( d ) related to the cost of divide andconquer methods such as the half-gcd algorithm: M (cid:48) ( d ) = (cid:80) ≤ i ≤(cid:100) log ( d ) (cid:101) i M (2 − i d ) for d ≥
1, and M (cid:48) ( d ) = ≤ d ≤
1. By definition one has M (cid:48) ( d ) ≥ M ( d ) ≥ d ≥
0, and theidentity M (cid:48) (2 d ) = M (cid:48) ( d ) + M (2 d ) for d ≥ M (cid:48) ( d ) is superlinear: 2 M (cid:48) ( d ) ≤ M (cid:48) (2 d )for all d ≥
1. Assuming H sl yields the asymptotic bound M (cid:48) ( d ) ∈ O ( M ( d ) log( d )) (which can bepessimistic: for example M (cid:48) ( d ) ∈ O ( M ( d )) holds in the case of Toom-Cook multiplication); inparticular, H sl and H ω imply M (cid:48) ( d ) ∈ O ( d ω − − ε ) for some ε >
0. Furthermore, if one assumes H sl and H sm , then M (cid:48) ( · ) is submultiplicative as well: M (cid:48) ( d d ) ≤ M (cid:48) ( d ) M (cid:48) ( d ) for all d , d ≥ K [ x ] m × m of degree at most d ≥ O ( m ω M ( d )) operations in K ; this holds when M ( d ) corresponds to oneof the above-mentioned polynomial multiplication algorithms. However this might not hold forother functions M ( d ) or matrix multiplication algorithms that are yet to be discovered. Besides,this bound is slightly worse than the best known ones [9, 24]; for example, Cantor and Kaltofen’salgorithm performs polynomial matrix multiplication in O ( m ω d log( d ) + m d log( d ) log(log( d )))field operations, which is finer than the bound O ( m ω M ( d )) with M ( d ) = Θ ( d log( d ) log(log( d )))in that case. This simplification is frequent in the polynomial matrix literature, and it is madehere for the sake of presentation, to improve the clarity of our main complexity results and of theanalyses that lead to them. Previous algorithms based on linear algebra over K for computing the characteristic polyno-mial of M ∈ K m × m mainly fall in three types of methods. Traces of powers: combining the traces of the first n powers of the input matrix using theNewton identities reveals the coe ffi cients of the characteristic polynomial. Known as theFaddeev-LeVerrier algorithm, it was introduced in [34], refined and rediscovered in [43,16, 18], and used in [11] to prove that the problem is in the NC parallel complexity class. Determinant expansion formula: introduced in [41] and improved in [6], this approach doesnot involve division, and is therefore well suited for computing over integral domains.Later developments in this field include [1, 30].
Krylov methods: based on sequences of iterates of vectors under the application of the matrix:( v , Mv , M v , . . . ). These methods rely on the fact that the first linear dependency betweenthese iterates defines a polynomial which divides the characteristic polynomial. Some al-gorithms construct the Krylov basis explicitly [31, 20, 14], while others can be interpretedas an implicit Krylov iteration with structured vectors [12, 39].Methods based on traces of powers use O ( m ) or O ( m ω + ) field operations, and are mostly com-petitive for their parallel complexity. Methods based on determinant expansions use O ( m ) or O ( m ω + ) field operations and are relevant for division-free algorithms. Lastly, the Krylov meth-ods run in O ( m ) [12, 14] or O ( m ω log m ) [31] field operations with deterministic algorithms, orin O ( m ω ) field operations with the Las Vegas randomized algorithm in [39].Note that the characteristic polynomial of M cannot be computed faster than the determinantof M , since the latter is the constant coe ffi cient of the former. Furthermore, under the model of3omputation trees, the determinant of m × m matrices cannot be computed faster than the productof two m × m matrices [8, Sec. 16.4], a consequence of Baur and Strassen’s theorem [2].Another type of characteristic polynomial algorithms is based on operations on matrices over K [ x ], called polynomial matrices in what follows. Indeed the characteristic polynomial may beobtained by calling a determinant algorithm on the characteristic matrix x I m − M , which is in K [ x ] m × m . Existing algorithms, which accept any matrix in K [ x ] m × m of degree d as input, include • the evaluation-interpolation method, which costs O ( m ω + d ) field operations, assumes thefield K is large enough, and relies on the computation of about md determinants in K m × m ; • the algorithm of Mulders and Storjohann [36] based on weak Popov form computation,which uses O ( m d ) field operations; • retrieving the determinant as the product of the diagonal entries of the Smith form, itselfcomputed by a Las Vegas randomized algorithm in O ( m ω M (cid:48) ( d ) log( m ) ) field operations[44, Prop. 41], assuming M (cid:48) ( d ) ∈ O ( d ω − ); • the algorithm based on unimodular triangularization in [32], which is deterministic anduses O ( m ω d log( d ) a log( m ) b ) field operations for some constants a , b ∈ Z > .In the last two items the cost bound is, up to logarithmic factors, the same as the cost of multiply-ing matrices K [ x ] m × m of degree d by relying on both fast linear algebra over K and fast arithmeticin K [ x ], as showed in [9]. The last two of these cost bounds do involve factors logarithmic in m ,whereas the first two have an exponent on m which exceeds ω .In summary, the fastest characteristic polynomial algorithms either are randomized or have acost a logarithmic factor away from the lower bound. This paper, with Theorem 1.1, bridges thisgap by proposing the first deterministic algorithm with cost O ( m ω ). Our algorithm falls within the category of polynomial matrix determinant computation. Yetunlike the above-listed approaches ours is tailored to a specific family of polynomial matrices,which contains the characteristic matrix x I m − M : the family of row reduced matrices [51, 29].Restricting to such matrices provides us with good control of the degrees in computations; asa typical example, it is easy to predict the degree of a vector-matrix product v ( x I m − M ) byobserving the degrees in v , without actually computing the product. As we explain below, thisdegree control allows us to avoid searches of degree profiles, which would add logarithmic termsto the cost. Although the characteristic matrix has other properties besides row reducedness (ithas degree 1, and is in Popov form [40] hence column reduced), we do not exploit them.When appropriate, the average row degree D / m , where D is the sum of the degrees of the rowsof the matrix, is chosen as a measure of the input degree which refines the matrix degree d usedabove. This gives cost bounds more sensitive to the input degrees and also, most importantly,leverages the fact that even if the algorithm starts from a matrix with uniform degrees such as x I m − M , it may end up handling matrices with unbalanced row degrees in the process. Theorem 1.2.
Assuming H sl , H sm , and H ω , there is an algorithm which takes as input a rowreduced matrix A ∈ K [ x ] m × m and computes its determinant usingO (cid:0) m ω M (cid:48) ( D / m ) (cid:1) ⊆ O (cid:0) m ω M (cid:48) (cid:0) deg( A ) (cid:1)(cid:1) operations in K , where D = deg(det( A )) is equal to the sum of the degrees of the rows of A . A )) is the sum of row degrees is the consequence of row reducedness [29],and the cost bound inclusion follows from deg(det( A )) ≤ m deg( A ). Taking A = x I m − M for M ∈ K m × m , Theorem 1.1 is a direct corollary of Theorem 1.2. The only assumption needed inTheorem 1.1 is ω >
2, since it implies the existence of a polynomial multiplication algorithmsuch that H sl , H sm , and H ω hold, such as Cantor and Kaltofen’s algorithm [9].Previous polynomial matrix determinant algorithms with costs of the order of m ω deg( A ), upto logarithmic factors, have been listed above: a randomized one from [44], and a deterministicone from [32]. To our knowledge, this paper gives the first description of an algorithm achievingsuch a cost involving no factor logarithmic in m . Our approach partially follows the algorithm of[32], but also substantially di ff ers from it in a way that allows us to benefit from the reducednessof A . The cost bound O ( m ω M (cid:48) (deg( A ))) has been obtained before in [21, Sec. 4.2.2] in theparticular case of a “su ffi ciently generic ” matrix A . In that case, both the algorithm of [32] andthe one here coincide and become the algorithm of [21, Sec. 4.2.2]; when A is the characteristicmatrix x I m − M , this also relates to the fast algorithm in [31, Sec. 6] for a generic M . For the sake of presentation, suppose m is a power of 2. Writing A = [ A A A A ] with the A i ’s ofdimensions ( m / × ( m / (cid:34) ∗ ∗ K K (cid:35) (cid:34) A A A A (cid:35) = (cid:34) R ∗ (cid:35) where the entries “ ∗ ” are not computed, B = K A + K A , and R and [ K K ] are computedfrom [ A A ] as a row basis and a kernel basis , respectively (see Section 2.2 for definitions). Thenthe leftmost matrix in the above identity is unimodular [32, Lemma 3.1] and thus, up to a constantfactor, det( A ) can be computed recursively as det( R ) det( B ).A first observation is that neither the kernel basis computation nor the matrix multiplicationgiving B is an obstacle towards a cost which is free of log( m ). (The reason why this is not obviousfor matrix multiplication is that here the product K A + K A involves matrices with possiblyunbalanced degrees, and the current fastest method for this splits the computation into O (log( m ))multiplications of smaller matrices which have balanced degrees [57, Sec. 3.6].) Indeed we showthat, under the above assumptions on M ( · ), the cost of these operations is in O ( m ω M ( D / m )) and O ( m ω M (cid:48) ( D / m )), thus only involving factors logarithmic in D / m . In previous work, either costbounds were hiding logarithmic factors [57] or were derived without assuming H sm and had theform O ( m ω − M ( D )) and O ( m ω − M (cid:48) ( D )) [27], thus resulting in factors logarithmic in D . Provingthis observation is straightforward from the analyses in [57, 27] (see Section 2.5). This is a firstkey towards our main result: the characteristic matrix has D = m , and O ( m ω M (1)) is the same as O ( m ω ) whereas O ( m ω − M ( m )) involves logarithmic factors.However, the computation of the row basis R remains an obstacle which prevents the algo-rithm of [32] from being a candidate for Theorem 1.2. Indeed, among the row basis algorithmswe are aware of, only one has a cost bound which fits into our target up to logarithmic factors: theone of [55]. It relies on three kernel bases computations, and while one of them is similar to thecomputation of [ K K ] and is handled via the algorithm of [57], the two others have di ff erentconstraints on the input and was the subject of a specific algorithm described in [55, Sec. 4]. In Precisely, if the upper triangular, row-wise Hermite normal form of A has diagonal entries (1 , . . . , , λ det( A )), forsome λ ∈ K \ { } making λ det( A ) monic. m . The algorithm has a loop over Θ (log( m ))iterations, each of them calling [54, Algo. 2] for minimal approximant bases with unbalanced in-put. This approximant basis algorithm may spend a logarithmic number of iterations for findingsome degree profile of the output basis, in a way reminiscent of Keller-Gehrig’s algorithm in [31,Sec. 5] which finds the lengths of Krylov sequences (the link between the two situations becomesmore explicit for approximant bases at small orders, see [27, Sec. 7]).Our attempts at accelerating the row basis algorithm of [55] having not succeeded, the algo-rithm in this paper follows an approach which is more direct at first: remove the obstacle. Insteadof computing a row basis R and relying on the identity det( A ) = det( R ) det( B ) (up to a constant),keep the first block row of A : (cid:34) I m / K (cid:35) (cid:34) A A A A (cid:35) = (cid:34) A A (cid:35) (1)and rely on the identity det( A ) = det( A ) det( B ) / det( K ). The nonsingularity of A and K iseasily ensured thanks to the assumption that A is reduced, as discussed in Section 1.5.This leads to an unusual recursion scheme: we are not aware of a similar scheme being usedin the literature on computational linear algebra. The algorithm uses three recursive calls with( m / × ( m /
2) matrices whose determinant has degree at most D / D for the third; our complexity analysis in Section 3.3 shows that such a recursion gives the costin Theorem 1.2. Precisely, if deg(det( A )) ≤ D / K )) ≤ D /
2, yielding the two calls in half the degree; otherwise the algorithmuses inexpensive row and column operations on A to reduce to the case deg(det( A )) ≤ D / A and B thanks to a straightforward gener-alization of [42, Sec. 3], and we describe a new algorithm which handles the more involved caseof K . When outlining the approach of our determinant algorithm via the identity in Eq. (1), weimplicitly assumed that the matrices used as input in recursive calls, i.e. A and K and B , dosatisfy the input requirement of row reducedness: this is not necessarily the case, even if startingfrom a reduced matrix A .Concerning A , one may locate such a reduced submatrix of A and then permute rows andcolumns of A (which only a ff ects the sign of det( A )) to make this submatrix become the leadingprincipal submatrix A . This is a classical operation on reduced matrices which suggests using aform slightly stronger than reducedness called weak Popov form [36] (see Section 2.4). Assumingthat A has this form ensures that its leading principal submatrix A has it as well. This assumptionis acceptable in terms of complexity since one can transform a reduced A into a weak Popov P by means of fast linear algebra in a cost negligible compared to our target [42, Sec. 3]; note that A and P have the same determinant up to an easily found constant (see Algorithm 1).Next, the cases of K and B are strongly linked. First, we will not discuss K but the wholekernel basis [ K K ]. The fastest known algorithm for computing such a basis is that of [57],and for best e ffi ciency it outputs a matrix in shifted reduced form, which is a generalization ofreducedness involving degree weights given by a tuple s ∈ Z m called a shift (see Sections 2.36nd 2.4 for definitions); the non-shifted case is for s = . As in the determinant algorithm of [32],here the shift for [ K K ] is taken as the list of row degrees of A , denoted by s = rdeg( A ); for thecharacteristic matrix one has s = (1 , . . . ,
1) but non-uniform shifts may arise in recursive calls.We want [ K K ] to be not only s -reduced, but in s -weak Popov form: a direct consequence isthat B is in weak Popov form, and is thus suitable input for a recursive call.To obtain [ K K ] we use the kernel basis algorithm of [57] and transform its output into s -weak Popov form. A minor issue is that the fastest known algorithm for such transformationswas written in [42, Sec. 3] for non-shifted forms; yet it easily extends to shifted forms as weshow in Section 4, obtaining the next result. Theorem 1.3.
There is an algorithm R educed T o W eak P opov which takes as input a matrix A ∈ K [ x ] m × n with m ≤ n and a shift s ∈ Z n such that A is in s -reduced form, and returns an s -weakPopov form of A using O ( m ω − nD + m ω − n ) operations in K , where D = | rdeg s ( A ) | − m · min( s ) . Here, following usual notation recalled in Section 2.1, | rdeg s ( A ) | is the sum of the s -degrees ofthe rows of A . This result extends [42, Thm. 13] since for s = the quantity D is the sum of therow degrees of A and in particular D ≤ m deg( A ), leading to the cost bound O ( m ω − n deg( A )).To summarize, at this stage we have outlined how to ensure, without exceeding our targetcost bound, that A and B are valid input for recursive calls, i.e. are in weak Popov form. Havingdet( A ) and det( B ), it remains to find det( K ) and then the sought det( A ) follows. We noted that,to ensure the form of B but also for e ffi ciency reasons, the kernel basis [ K K ] is computedin s -weak Popov form for the shift s = rdeg( A ). This causes the main di ffi culty related to ourmodification of the determinant algorithm of [32]: K is not valid input for a recursive call sinceit is in v -weak Popov form for some shift v , a subtuple of s which is possibly nonzero.A first idea is to extend our approach to the shifted case, allowing recursive calls with sucha v -reduced matrix: this is straightforward but gives an ine ffi cient algorithm. Indeed, along therecursion the shift drifts away from its initial value and becomes arbitrarily large and unbalancedwith respect to the degrees of the input matrices of recursive calls. For example, as mentionedabove the sum of row degrees of the initial non-shifted m × m matrix A is D = deg(det( A )),whereas for the v -shifted ( m / × ( m /
2) matrix K we only have the same bound D instead ofone related to deg(det( K )) itself, which is known to be at most D / D and D /
2, will only grow as the algorithm goes down the tree of recursive calls,meaning that degrees in matrices handled recursively are not su ffi ciently well controlled.Another idea is to compute a -reduced matrix which has the same determinant as K . Find-ing a -reduced form of K within our target cost seems to be a di ffi cult problem. The best knownalgorithms for general -reduction involve log( m ) factors, either explicitly [23] or implicitly [38](in the latter approach one starts by using the above-discussed triangularization procedure of [32]which we are modifying here to avoid log( m ) factors). More specific algorithms exploit the formof K , interpreting the problem as a change of shift from v to ; yet at the time of writing e ffi cientchanges of shifts have only been achieved when the target shift is larger than the origin shift [27,Sec. 5], a fact that o ff ers degree control for the transformation between the two matrices. Anotherpossibility is to compute the so-called v -Popov form P of K , since its transpose P T is -reducedby definition (see Section 2.4), and det( P T ) = det( P ) is det( K ) up to a constant. However thissu ff ers from the same issue, as computing P is essentially the same as changing the shift v intothe nonpositive shift − δ , where δ is the list of diagonal degrees of K [42, 26].To circumvent these issues, we use the property that the transpose K T of a v -reduced matrixis in − d -reduced form where d = rdeg v ( K ). This fact naturally comes up here since det( K ) = K T ), but seems otherwise rarely exploited in polynomial matrix algorithms: in fact we didnot find a previous occurrence of it apart from related degree considerations in [55, Lem. 2.2].Transposing the above two approaches using K T instead of K , we observe that computinga -reduced form of K T is a change of shift from − d to , and computing the − d -Popov form P of K T is essentially a change of shift from − d to − δ . In both cases the target shift is larger thanthe origin shift, implying that the kernel-based change of shift of [27, Sec. 5] involves matricesof well-controlled degrees. Still, this is not enough to make this change of shift e ffi cient as such,the di ffi culty being now that the average row degree of K T may not be small: only its averagecolumn degree, which corresponds to the average row degree of K , is controlled.Our solution uses the second approach, computing the − d -Popov form P , because it o ff ers thea priori knowledge that the column degrees of P are exactly δ . We exploit this degree knowledgeto carry out partial linearization techniques, originally designed for approximant bases [45, 28],which we extend here to kernel bases. These techniques allow us to reduce our problem to akernel basis computation where the matrix entries have uniformly small degrees, implying thatit can be e ffi ciently handled via the minimal approximant basis algorithm PM-B asis from [21].The next result summarizes the new algorithmic tool developed in Section 5 for finding P . Theorem 1.4.
Let s ∈ Z m , let A ∈ K [ x ] m × m be in − s -weak Popov form, let δ ∈ Z m ≥ be the − s -pivot degree of A , and assume that s ≥ δ . There is an algorithm W eak P opov T o P opov whichtakes as input ( A , s ) and computes the − s -Popov form of A by • performing PM-B asis at order less than | s | / m + on an input matrix of row dimension atmost m and column dimension at most m, • multiplying the inverse of a matrix in K m × m by a matrix in K [ x ] m × m of column degree δ , • and performing O ( m ) extra operations in K .Thus, computing the − s -Popov form of A can be done in O ( m ω M (cid:48) ( | s | / m )) operations in K . This theorem is a generalization of [42, Sec. 4] to shifted forms, for shifts − s that satisfy theassumption s ≥ δ . Indeed, if A is -weak Popov, then one recovers [42, Thm. 20] by taking s = (deg( A ) , . . . , deg( A )) in the above theorem. For comparison, the naive generalization of [42,Sec. 4] to shifted forms runs in O ( m ω M (cid:48) (max( s ))), which exceeds our target complexity as soonas max( s ) (cid:29) | s | / m . Hence the use of partial linearization techniques, which were not needed inthe non-shifted case featuring max( s ) = | s | / m = deg( A ).As mentioned above, our Algorithm W eak P opov T o P opov is based on the computation of akernel basis with a priori knowledge of the degree profile of the output. This kernel problem isvery close to the one handled in [55, Sec. 4], except that in this reference one only has upperbounds on the output degrees, implying a certain number—possibly logarithmic in m —of callsto PM-B asis to recover the output and its actual degrees. In the same spirit but in the context ofapproximant bases, [28, Sec. 5] uses partial linearization techniques to reduce an arbitrary inputwith known output degrees to essentially one call to PM-B asis , whereas [54, Algo. 2] assumesweaker output degree information and makes a potentially logarithmic number of calls to PM-B asis . We plan to implement our characteristic polynomial algorithm in the LinBox ecosystem [47,48]. First prototype experiments suggest that, for large finite fields, it could be competitivewith the existing fastest-known implementation, based on the randomized algorithm of [39].The native support for small fields of our algorithm should outperform the algorithm of [39]8hich requires expensive field extensions. Another perspective stems from the remark that ouralgorithm resorts to fast polynomial multiplication (see assumption H ω ), while previous onesdid not [31, 39]: we woud like to understand whether the same cost can be achieved by a purelylinear algebraic approach. Finally, perhaps the most challenging problem related to characteristicpolynomial computation is to compute Frobenius forms deterministically in the time of matrixmultiplication, and more generally computing Smith forms of polynomial matrices with a costfree of factors logarithmic in the matrix dimension.
2. Preliminaries on polynomial matrices
In this section we present the notation as well as basic definitions and properties that will beused throughout the paper.
Tuples of integers will often be manipulated entry-wise. In particular, for tuples s , t ∈ Z n ofthe same length n , we write s + t for their entry-wise sum, and the inequality s ≤ t means thateach entry in s is less than or equal to the corresponding entry in t . The concatenation of tuplesis denoted by ( s , t ). We write | t | for the sum of the entries of t . The tuple of zeros is denoted by = (0 , . . . , m × n matrix A over some ring, we write A i , j for its entry at index ( i , j ). We extendthis to submatrices: given sets I ⊆ { , . . . , m } and J ⊆ { , . . . , n } of row and column indices, wewrite A I , J for the submatrix of A formed by its entries indexed by I × J . Besides, A I , ∗ stands forthe submatrix of A formed by its rows with index in I , and we use the similar notation A ∗ , J . Thetranspose of A is denoted by A T . The identity matrix of size n is denoted by I n , while the n × n matrix with 1 on the antidiagonal and 0 elsewhere is denoted by J n . In particular, when writing s J n for a tuple s = ( s , . . . , s n ) ∈ Z n , we mean the reversed tuple s J n = ( s n , . . . , s ).Now consider A with polynomial entries, i.e. A ∈ K [ x ] m × n . The degree of A is denoted bydeg( A ) and is the largest of the degrees of its entries, or −∞ if A = . The row degree of A isthe tuple rdeg( A ) ∈ ( Z ≥ ∪ {−∞} ) m whose i th entry is max ≤ j ≤ n (deg( A i , j )). More generally, for atuple s = ( s , . . . , s n ) ∈ Z n , the s -row degree of A is the tuple rdeg s ( A ) ∈ ( Z ∪ {−∞} ) m whose i thentry is max ≤ j ≤ n (deg( A i , j ) + s j ). In this context, the tuple s is commonly called a (degree) shift [4]. The (shifted) column degree of A is defined similarly.We write X s for the n × n diagonal matrix diag( x s , . . . , x s n ) which is over the ring K [ x , x − ] ofLaurent polynomials over K . Note that, hereafter, Laurent polynomials will only arise in proofsand explanations, more specifically in considerations about shifted degrees: they never arise inalgorithms, which for the sake of clarity only involve polynomials in K [ x ]. The usefulness ofthis matrix X s will become clear in the definition of leading matrices in the next subsection.The next lemma gives a link between shifted row degrees and shifted column degrees. Wewill mostly use the following particular case of it: the column degree of A is at most d ∈ Z n ≥ (entry-wise) if and only if the − d -row degree of A is nonpositive. Lemma 2.1 ([55, Lemma 2.2]) . Let A be a matrix in K [ x ] m × n , d be a tuple in Z n , and t be atuple in Z m . Then, cdeg t ( A ) ≤ d if and only if rdeg − d ( A ) ≤ − t . .2. Bases of modules, kernel bases and approximant bases We recall that any K [ x ]-submodule M of K [ x ] n is free, and admits a basis formed by r elements of K [ x ] n , where r ≤ n is called the rank of M [see e.g. 15]. Such a basis can thus berepresented as an r × n matrix B over K [ x ] whose rows are the basis elements; this basis matrix B has rank r .For a matrix A ∈ K [ x ] m × n , its row space is the K [ x ]-submodule { pA , p ∈ K [ x ] × m } of K [ x ] × n ,that is, the set of all K [ x ]-linear combinations of its rows. If B ∈ K [ x ] r × n is a basis of this rowspace, then B is said to be a row basis of A ; in particular, r is the rank of B and of A .The left kernel of A , denoted by K ( A ), is the K [ x ]-module { p ∈ K [ x ] × m | pA = } . A matrix K ∈ K [ x ] k × m is a left kernel basis of A if its rows form a basis of K ( A ), in which case k = m − r .Similarly, a right kernel basis of A is a matrix K ∈ K [ x ] n × ( n − r ) whose columns form a basis ofthe right kernel of A .Given positive integers γ = ( γ , . . . , γ n ) ∈ Z n > and a matrix F ∈ K [ x ] m × n , the set of approxi-mants for F at order γ [see e.g. 50, 3] is the K [ x ]-submodule of K [ x ] × m defined as A γ ( F ) = { p ∈ K [ x ] × m | pF = mod X γ } . The identity pF = mod X γ means that pF ∗ , j = x γ j for 1 ≤ j ≤ n . Since all m rows ofthe matrix x max( γ ) I m are in A γ ( F ), this module has rank m . We will often compute with polynomial matrices that have a special form, called the (shifted)reduced form . It corresponds to a type of minimality of the degrees of such matrices, and alsoprovides good control of these degrees during computations as illustrated by the predictabledegree property [17] [29, Thm. 6.3-13] which we recall below. In this section, we introduce thenotion of row reducedness; to avoid confusion, we will not use the similar notion of columnreducedness in this paper, and thus all further mentions of reducedness refer to row reducedness.For shifted reduced forms, we follow the definitions in [4, 5]. Let A ∈ K [ x ] m × n and s ∈ Z n ,and let t = ( t , . . . , t m ) = rdeg s ( A ). Then, the s -leading matrix of A is the matrix lm s ( A ) ∈ K m × n whose entry ( i , j ) is the coe ffi cient of degree t i − s j of the entry ( i , j ) of A , or 0 if t i = −∞ .Equivalently, lm s ( A ) is the coe ffi cient of degree zero of X − t AX s , whose entries are in K [ x − ].The matrix A is said to be in s -reduced form if its s -leading matrix has full row rank. In particular,a matrix in s -reduced form must have full row rank.For a matrix M ∈ K [ x ] k × m , we have rdeg s ( MA ) ≤ rdeg t ( M ) and this is an equality when nocancellation of leading terms occurs in this left-multiplication. The predictable degree propertystates that A is s -reduced if and only if rdeg s ( MA ) = rdeg t ( M ) holds for any M ∈ K [ x ] k × m . Hereis a useful consequence of this characterization. Lemma 2.2.
Let s ∈ Z n , let A ∈ K [ x ] m × n , and let t = rdeg s ( A ) . If A is s -reduced, then theidentity lm s ( MA ) = lm t ( M )lm s ( A ) holds for any M ∈ K [ x ] k × m .Proof. Let d = rdeg s ( MA ). By definition, lm s ( MA ) is the coe ffi cient of degree 0 of the matrix X − d MAX s = X − d MX t X − t AX s , whose entries are in K [ x − ]. Besides, since rdeg s ( A ) = t andsince the predictable degree property gives rdeg t ( M ) = d , the matrices X − d MX t and X − t AX s are over K [ x − ] and their coe ffi cients of degree 0 are lm t ( M ) and lm s ( A ), respectively.Another characterization of matrices in s -reduced form is that they have minimal s -row de-gree among all matrices which represent the same K [ x ]-module [52, Def. 2.13]; in this paper, wewill use the following consequence of this minimality.10 emma 2.3. Let M be a submodule of K [ x ] × n of rank m, let s ∈ Z n , and let t ∈ Z m be the s -rowdegree of some s -reduced basis of M . Without loss of generality, assume that t is nondecreasing.Let B ∈ K [ x ] m × n be a matrix of rank m whose rows are in M , and let d ∈ Z m be its s -row degreesorted in nondecreasing order. If d ≤ t , then B is an s -reduced basis of M , and d = t .Proof. Up to permuting the rows of B , we assume that rdeg s ( B ) = d without loss of generality.Let A ∈ K [ x ] m × n be an s -reduced basis of M such that rdeg s ( A ) = t . Since the rows of B are in M , there exists a matrix U ∈ K [ x ] m × m such that B = UA ; and U is nonsingular since B and A have rank m . Since A is s -reduced, the predictable degree property applies, ensuring that d = rdeg s ( B ) = rdeg s ( UA ) = rdeg t ( U ) . This means that deg( U i , j ) ≤ d i − t j for all 1 ≤ i , j ≤ m .Now, assume by contradiction that d = t does not hold. Thus, d k < t k for some 1 ≤ k ≤ m .Then, for i ≤ k and j ≥ k we have d i ≤ d k < t k ≤ t j , hence deg( U i , j ) <
0. Thus, the submatrix U { ,..., k } , { k ,..., m } is zero, which implies that U is singular; this is a contradiction, hence d = t .Since t is nondecreasing, the inequality deg( U i , j ) ≤ t i − t j implies that U is a block lowertriangular matrix whose diagonal blocks have degree 0; hence these blocks are invertible matricesover K , and U is unimodular [see 42, Lemma 6 for similar degree considerations, starting fromstronger assumptions on A and B ]. Thus, B is a basis of M .Furthermore, it is easily observed that lm d ( U ) ∈ K m × m is block lower triangular with the sameinvertible diagonal blocks as U ; hence lm d ( U ) is invertible. On the other hand, Lemma 2.2 statesthat lm s ( B ) = lm d ( U )lm s ( A ). Thus lm s ( B ) has rank m = rank(lm s ( A )), and B is s -reduced. For algorithmic purposes, it is often convenient to work with reduced forms that satisfy someadditional requirements, called weak Popov forms . These are intrinsically related to the notionof pivot of a polynomial matrix.For a nonzero vector p ∈ K [ x ] × n and a shift s ∈ Z n , the s -pivot of p is its rightmost entry p j such that deg( p j ) + s j = rdeg s ( p ) [4, 36]; it corresponds to the rightmost nonzero entry of lm s ( p ).The index j = π and the degree deg( p π ) = δ of this entry are called the s -pivot index and s -pivotdegree, respectively. For brevity, in this paper the pair ( π, δ ) is called the s -pivot profile of p . Byconvention, the zero vector in K [ x ] × n has s -pivot index 0 and s -pivot degree −∞ . These notionsare extended to matrices A ∈ K [ x ] m × n by forming row-wise lists. For example, the s -pivot indexof A is π = ( π , . . . , π m ) ∈ Z m > where π i is the s -pivot index of the row A i , ∗ . The s -pivot degree δ and the s -pivot profile ( π i , δ i ) ≤ i ≤ m of A are defined similarly.Then, A is said to be in s -weak Popov form if it has no zero row and π is strictly increasing;if this holds up to row permutation, i.e. the entries of π are pairwise distinct, then A is said tobe in s -unordered weak Popov form. Furthermore, a matrix is in s -Popov form if it is in s -weakPopov form, its s -pivots are monic, and each of these s -pivots has degree strictly larger than theother entries in the same column. For a given K [ x ]-submodule M of K [ x ] n , there is a uniquebasis of M which is in s -Popov form [4].For a given matrix B , the matrix A is said to be an s -reduced (resp. s -weak Popov, s -Popov)form of B if A is a row basis of B and A is in s -reduced (resp. s -weak Popov, s -Popov) form.Like for s -reducedness, the property of a matrix A ∈ K [ x ] m × n to be in s -weak Popov formdepends only on its s -leading matrix lm s ( A ) ∈ K m × n , namely on the fact that it has a staircaseshape. Indeed, A is in s -weak (resp. s -unordered weak) Popov form if and only if lm s ( A ) has no11ero row and J m lm s ( A ) J n is in row echelon form (resp. in row echelon form up to row permu-tation); this was used as a definition by Beckermann et al. [4, 5]. In particular, for any constantmatrix C ∈ K m × n , we have lm ( C ) = C and therefore C is in -weak (resp. -unordered weak)Popov form if and only if it has no zero row and J m CJ n is in row echelon form (resp. in rowechelon form up to row permutation). Taking C = lm s ( A ), the next lemma follows. Lemma 2.4.
Let A ∈ K [ x ] m × n and let s ∈ Z n . Then, A is in s -weak (resp. s -unordered weak)Popov form if and only if lm s ( A ) is in -weak (resp. -unordered weak) Popov form. Furthermore, if A is in s -weak Popov form and ( j , . . . , j m ) is the list of indices of pivotcolumns in the row echelon form J m lm s ( A ) J n (in other words, this list is the column rank profileof that matrix), then the s -pivot index of A is equal to ( n + − j m , . . . , n + − j ). This leadsto the following lemma which states that the s -pivot profile is an invariant of left-unimodularlyequivalent s -weak Popov forms, generalizing the fact that for matrices over K the set of indicesof pivot columns is an invariant of left-equivalent row echelon forms. Lemma 2.5.
Let s ∈ Z n and let A ∈ K [ x ] m × n be in s -unordered weak Popov form with s -pivotprofile ( π i , δ i ) ≤ i ≤ m . Then, the s -pivot profile of the s -Popov form of A is ( π σ ( i ) , δ σ ( i ) ) ≤ i ≤ m , where σ : { , . . . , m } → { , . . . , m } is the permutation such that ( π σ ( i ) ) ≤ i ≤ m is strictly increasing.Proof. Without loss of generality we assume that A is in s -weak Popov form, implying also σ ( i ) = i for 1 ≤ i ≤ m . Let P ∈ K [ x ] m × n be the s -Popov form of A : we want to prove that A and P have the same s -pivot index and the same s -pivot degree. Let U be the unimodular matrix suchthat P = UA ; then Lemma 2.2 yields lm s ( P ) = lm t ( U )lm s ( A ), where t = rdeg s ( A ). Since bothlm s ( P ) and lm s ( A ) have full row rank, lm t ( U ) ∈ K m × m is invertible. Then J m lm s ( P ) J n = J m lm t ( U ) J m J m lm s ( A ) J n holds, and thus the row echelon forms J m lm s ( P ) J n and J m lm s ( A ) J n have the same pivot columnssince J m lm t ( U ) J m ∈ K m × m is invertible. It follows from the discussion preceding this lemma that P has the same s -pivot index as A .As a consequence, P has the same s -pivot degree as A if and only if rdeg s ( P ) = rdeg s ( A ).Suppose by contradiction that there exists an index i such that rdeg s ( P i , ∗ ) < rdeg s ( A i , ∗ ). Then,build the matrix B ∈ K [ x ] m × n which is equal to A except for its i th row which is replaced by P i , ∗ . By construction, B has rank m (since it is in s -weak Popov form) and its rows are in the rowspace of A . Writing d for the tuple rdeg s ( B ) sorted in nondecreasing order, and u for the tuple t sorted in nondecreasing order, we have d ≤ u and d (cid:44) u , which contradicts Lemma 2.3. Hencethere is no such index i , and since this proof by contradiction is symmetric in A and P , there isno index i such that rdeg s ( A i , ∗ ) < rdeg s ( P i , ∗ ) either. Thus rdeg s ( A ) = rdeg s ( P ).We will also use the following folklore fact, which is a corollary of Lemma 2.2, and hasoften been used in algorithms for approximant bases or kernel bases in order to preserve thereducedness of matrices during the computation. Lemma 2.6.
Let A ∈ K [ x ] m × n and B ∈ K [ x ] k × m , and let s ∈ Z n and t = rdeg s ( A ) ∈ Z m . Then, • if A is s -reduced and B is t -reduced, then BA is s -reduced; • if A is in s -weak Popov form and B is t -weak Popov form, then BA is s -weak Popov formProof. Since A is s -reduced, Lemma 2.2 states that lm s ( BA ) = ML where L = lm s ( A ) ∈ K m × n and M = lm t ( B ) ∈ K k × m . The first item then follows from the fact that if M has rank k and L m , then ML has rank k . Similarly, the second item reduces to prove that, assuming M and L are in row echelon form with full row rank, then ML is also in row echelon form. Let( a , . . . , a k ) (resp. ( b , . . . , b m )) be the pivot indices of M (resp. L ). Then the i th row of ML is anonzero multiple of row a i of L combined with multiples of rows of L of index greater than a i .Consequently, the pivot indices of the rows of ML are b a < · · · < b a k , which proves that ML isin row echelon form.Finally, under assumptions that generalize the situation encountered in our determinant algo-rithm below, we show that the pivot entries of a kernel basis [ K K ] are located in its rightmostcolumns, that is, in K . Lemma 2.7.
Let t ∈ Z n , let F ∈ K [ x ] n × n be in t -weak Popov form, and let u = rdeg t ( F ) . Let G ∈ K [ x ] m × n and v ∈ Z m be such that v ≥ rdeg t ( G ) , and let K = [ K K ] ∈ K [ x ] m × ( m + n ) be a ( u , v ) -weak Popov basis of K ([ FG ]) , where K and K have m and n columns, respectively. Then,the ( u , v ) -pivot entries of K are all located in K ; in particular, K is in v -weak Popov form.Proof. Since the ( u , v )-pivot entry of a row is the rightmost entry of that row which reaches the( u , v )-row degree, it is enough to prove that rdeg v ( K ) ≥ rdeg u ( K ). First, from v ≥ rdeg t ( G ), weobtain rdeg v ( K ) ≥ rdeg rdeg t ( G ) ( K ). Now, by definition, rdeg rdeg t ( G ) ( K ) ≥ rdeg t ( K G ). Sincethe rows of K are in K ([ FG ]), we have K G = − K F , hence rdeg t ( K G ) = rdeg t ( K F ). Since F is t -reduced, we can apply the predictable degree property: rdeg t ( K F ) = rdeg u ( K ). Thisproves the sought inequality. For the last point, note that the ( u , v )-pivot entries of [ K K ]located in K correspond to v -pivot entries in K . Thus, since [ K K ] is in ( u , v )-weak Popovform with all ( u , v )-pivot entries in K , it follows that the v -pivot index of K is increasing. To conclude these preliminaries, we recall known fast algorithms for three polynomial matrixsubroutines used in our determinant algorithm: multiplication with unbalanced degrees, minimalapproximant bases, and minimal kernel bases; we give the corresponding complexity estimatesadapted to our context and in particular using our framework stated in Section 1.1.
Unbalanced multiplication.
Polynomial matrix algorithms often involve multiplication with ma-trix operands whose entries have degrees that may be unbalanced but still satisfy properties thatcan be exploited to perform the multiplication e ffi ciently. Here we will encounter products ofreduced matrices with degree properties similar to those discussed in [57, Sec. 3.6], where ane ffi cient approach for computing such products was given. Lemma 2.8.
There is an algorithm U nbalanced M ultiplication which takes as input a matrix B ∈ K [ x ] k × m with k ≤ m, a matrix A ∈ K [ x ] m × n with n ≤ m, and an integer D greater than orequal to both the sum of the positive entries of rdeg ( A ) and that of rdeg rdeg ( A ) ( B ) , and returnsthe product BA using O ( m ω M ( D / m )) operations in K , assuming H sm and H ω .Proof. Zhou et al. [57, Sec. 3.6] gave such an algorithm, yet with a cost analysis which hideslogarithmic factors; because these factors are our main concern here we will rely on the versionin [27, Sec. 4]. In this reference, Algorithm U nbalanced M ultiplication was described for squarematrices. One could adapt it to the case of rectangular A and B as in the statement above.However, for the sake of conciseness and with no impact on the asymptotic cost bound, weconsider the more basic approach of forming the square m × m matrices D = [ B0 ] and C = [ A 0 ],13omputing DC using the above-cited algorithm, and retrieving BA from it. Now, by construction,both the sum of the positive entries of rdeg ( C ) and that of rdeg rdeg ( A ) ( D ) are at most D , hence[27, Prop. 4.1] applies: defining ¯ m and d as the smallest powers of 2 greater than or equal to m and D / m , it states that the computation of DC costs O ( (cid:80) ≤ i ≤ log ( ¯ m ) i (2 − i ¯ m ) ω M (2 i d )) operations in K .Using H sm and H ω , which ensure respectively that M (2 i d ) ≤ M (2 i ) M ( d ) and M (2 i ) ∈ O (2 i ( ω − − ε ) )for some ε >
0, we obtain that this bound is in O ( ¯ m ω M ( d ) (cid:80) ≤ i ≤ log ( ¯ m ) − i ε ) ⊆ O ( ¯ m ω M ( d )). Thisis in O ( m ω M ( D / m )), since ¯ m ∈ Θ ( m ) and d ∈ Θ (1 + D / m ). Minimal approximant basis.
The second basic tool we will use is approximant bases computa-tion; for this, we will use the algorithm PM-B asis , originally described in [21]. Precisely, werely on the slightly modified version presented in [28] which ensures that the computed basis isin shifted weak Popov form.
Lemma 2.9.
There is an algorithm
PM-B asis which takes as input a tuple γ ∈ Z n > , a matrix F ∈ K [ x ] m × n with cdeg( F ) < γ , and a shift s ∈ Z m , and returns a basis of A γ ( F ) in s -weak Popovform using O (( m + n ) m ω − M (cid:48) (max( γ ))) operations in K .Proof. The algorithm is [28, Algo. 2]; to accommodate non-uniform order γ , it is called withinput order Γ = max( γ ) and input matrix FX ( Γ ,..., Γ ) − γ as explained in [28, Rmk. 3.3]. According to[28, Prop. 3.2], this costs O ((1 + nm ) (cid:80) ≤ i ≤(cid:100) log ( Γ ) (cid:101) i m ω M (2 − i Γ )) operations in K , which is preciselythe claimed bound by definition of M (cid:48) ( · ). Minimal kernel basis.
We will make use of the algorithm of Zhou et al. [57], which itself relieson unbalanced products and approximant bases, and returns a kernel basis in shifted reducedform e ffi ciently for input matrices with small average row degree. Lemma 2.10.
There is an algorithm K ernel B asis which takes as input a full column rank matrix F ∈ K [ x ] m × n with m ≥ n and m ∈ O ( n ) , and a shift s ∈ Z m ≥ such that s ≥ rdeg ( F ) , and returns abasis of K ( F ) in s -reduced form using O ( m ω M (cid:48) ( D / m )) operations in K , assuming H sl , H sm , H ω .Here D is the sum of the entries of s , and the sum of the s -row degree of this basis is at most D.Proof. The algorithm of Zhou et al. [57] computes an s -reduced basis K ∈ K [ x ] k × m of K ( F );precisely, this reference is about computing a basis of the right kernel in column reduced form,yet this naturally translates into left kernels and row reduced forms by taking suitable transposes.Furthermore, the last claim in the lemma follows from [57, Thm. 3.4], which states that any suchbasis K is such that | rdeg s ( K ) | ≤ | s | = D . For the complexity, we rely on the analysis in [27,Prop. B.1] which shows that, defining ¯ m and d as the smallest powers of 2 greater than or equalto m and D / m , this computation costs O log ( ¯ m ) (cid:88) j = j log (2 − j ¯ m ) (cid:88) i = i (2 − i − j ¯ m ) ω M (cid:16) i + j d (cid:17) + log (2 j d ) (cid:88) i = i (2 − j ¯ m ) ω M (cid:16) j − i d (cid:17) operations in K . Now the same analysis as in the proof of Lemma 2.8 shows that, assuming H sm and H ω , the first inner sum is in O ((2 − j ¯ m ) ω M (2 j d )), and by definition of M (cid:48) ( · ) the second innersum is in O ((2 − j ¯ m ) ω M (cid:48) (2 j d )). Thus the total cost is in O ( (cid:80) ≤ j ≤ log ( ¯ m ) j (2 − j ¯ m ) ω M (cid:48) (2 j d )), which isin O ( ¯ m ω M (cid:48) ( d ) (cid:80) ≤ j ≤ log ( ¯ m ) j (1 − ω ) M (cid:48) (2 j )) since H sl and H sm ensure that M (cid:48) ( · ) is submultiplicative.Similarly to the proof of Lemma 2.8, this bound is in O ( m ω M (cid:48) ( D / m )) thanks to H ω .14 . Determinant algorithm for reduced matrices In this section, we present the main algorithm in this paper, which computes the determinantof a matrix in reduced form using the subroutines listed in Section 2.5 as well as the algorithmsR educed T o W eak P opov and W eak P opov T o P opov from Theorems 1.3 and 1.4. Taking for grantedthe proof of these theorems in Sections 4 and 5, here we prove the correctness of our determinantalgorithm in Section 3.2 and analyse its complexity in Section 3.3, thus proving Theorem 1.2. ffi cient of the determinant. All bases of a given submodule of K [ x ] n of rank n havethe same determinant up to a constant factor, i.e. up to multiplication by an element of K \ { } .Many algorithms operating on polynomial matrices such as PM-B asis and K ernel B asis computesuch bases, so that their use in a determinant algorithm typically leads to obtaining the soughtdeterminant up to a constant factor; then finding the actual determinant requires to e ffi cientlyrecover this constant [see e.g. 32, Sec. 4]. Since in this paper we seek determinants of matricesin reduced form, this issue is easily handled using the next result. Lemma 3.1.
Let s ∈ Z n and A ∈ K [ x ] n × n . If A is in s -reduced form, the leading coe ffi cientof det( A ) is det(lm s ( A )) . In particular, if A is in s -weak Popov form, the leading coe ffi cient of det( A ) is the product of the leading coe ffi cients of the diagonal entries of A .Proof. The second claim is a direct consequence of the first, since for A in s -weak Popov form,lm s ( A ) is lower triangular with diagonal entries equal to the leading coe ffi cients of the diagonalentries of AX s , which are the leading coe ffi cients of the diagonal entries of A . For the first claimin the case s = , we refer to [29, Sec. 6.3.2], and in particular Eq. (23) therein. Now, for anarbitrary s and A in s -reduced form, we consider the nonnegative shift t = s − (min( s ) , . . . , min( s ))and observe that lm s ( A ) = lm ( AX t ), hence AX t is -reduced, and thus the leading coe ffi cientof det( AX t ) = det( A ) det( X t ) (which is the same as that of det( A )) is equal to det(lm ( AX t )) = det(lm s ( A )). From shifted to non-shifted.
In Section 1.5, we explained that one step of our algorithm consistsin finding det( K ) for a matrix K in v -weak Popov form, and achieves this by computing the − d -Popov form of K T , which is already in − d -weak Popov form. The next lemma substantiates this;note that in Section 1.5 we had left out the reversal matrix J n for the sake of exposition. Lemma 3.2.
Let v ∈ Z n , let K ∈ K [ x ] n × n , and let d = rdeg v ( K ) .(a) if lm v ( K ) has no zero column, then lm − d ( K T ) = lm v ( K ) T and rdeg − d ( K T ) = − v ;(b) if K is in v -reduced form, then K T is in − d -reduced form;(c) if K is in v -weak Popov form, then J n K T J n is in − d J n -weak Popov form;(d) if furthermore P is the − d J n -Popov form of J n K T J n , then P T is in -weak Popov form and det( K ) = det(lm v ( K )) det( P T ) .Proof. By definition, lm v ( K ) T is the coe ffi cient of degree 0 of ( X − d KX v ) T = X v K T X − d , whichis a matrix over K [ x − ]. The assumption on lm v ( K ) implies that this coe ffi cient of degree 0 hasno zero row. It follows that rdeg − d ( K T ) = − v and that this coe ffi cient of degree 0 is lm − d ( K T ).Item (b) follows from Item (a) by definition of shifted reduced forms.From now on, we assume that K is in v -weak Popov form. Then lm v ( K ) is invertible andlower triangular, and in particular lm − d ( K T ) = lm v ( K ) T . Since J n is a permutation matrix, we15btain lm − d J n ( J n K T J n ) = J n lm − d ( K T ) J n = J n lm v ( K ) T J n , which is invertible and lower triangular.Hence J n K T J n is in − d J n -weak Popov form.For Item (d), since P is n × n and in − d J n -Popov form, we have lm ( P T ) = I n , hence P T is in -weak Popov form. Furthermore, since J n K T J n is unimodularly equivalent to P , its determinantis det( K ) = det( J n K T J n ) = λ det( P ) for some λ ∈ K \ { } . Applying Lemma 3.1 to P showsthat det( P ) is monic, hence λ is the leading coe ffi cient of det( K ); applying the same lemma to K yields λ = det(lm v ( K )). Our main determinant algorithm is D eterminant O f W eak P opov (Algorithm 2), which takesas input a matrix in -weak Popov form and computes its determinant using recursive calls onmatrices of smaller dimension. Then, we compute the determinant of a -reduced matrix by firstcalling R educed T o W eak P opov to find a -weak Popov matrix which has the same determinantup to a nonzero constant, and then calling the previous algorithm on that matrix. This is detailedin Algorithm 1. Algorithm 1 D eterminant O f R educed ( A ) Input: a matrix A ∈ K [ x ] m × m in -reduced form. Output: the determinant of A . P ∈ K [ x ] m × m ← R educed T o W eak P opov ( A , ) ∆ ← D eterminant O f W eak P opov ( P ); (cid:96) ∆ ∈ K \ { } ← leading coe ffi cient of ∆ return det(lm ( A )) ∆ /(cid:96) ∆ The correctness of Algorithm 1 is obvious: according to Theorem 1.3 and Proposition 3.3, P is a -weak Popov form of A , and ∆ is the determinant of P up to multiplication by some elementof K \ { } . Thus det( A ) = (cid:96) ∆ for some (cid:96) ∈ K \ { } , and Lemma 3.1 yields (cid:96) = det(lm ( A )) /(cid:96) ∆ .Concerning the cost bound, Theorem 1.3 states that the first step uses O ( m ω (1 + D / m )) oper-ations in K , where D = | rdeg ( A ) | ; since A is -reduced, this is D = deg(det( A )) [29, Sec. 6.3.2].In the last step, the determinant computation costs O ( m ω ) operations, while scaling ∆ by a con-stant costs O ( D ) operations. The second step uses O ( m ω M (cid:48) ( D / m )) operations according to Propo-sition 3.3, under the assumptions H sl , H sm , and H ω .We now describe the main algorithm of this paper (Algorithm 2) and focus on its correctness.We also mention cost bounds for all steps of the algorithm that are not recursive calls, but wedefer the core of the complexity analysis to Section 3.3. Proposition 3.3.
Algorithm 2 is correct, and assuming that H sl , H sm , and H ω hold, it usesO ( m ω M (cid:48) ( D / m )) operations in K .Proof of correctness. The fact that A is in -weak Popov form has two consequences on the tuple s computed at Line 1: first, it is the -pivot degree of A (i.e. its diagonal degrees), and second,the sum D = | s | is equal to the degree of the determinant of A [29, Sec. 6.3.2].The main base case of the recursion is when m = K . We use a second base case at Line 3: if D =
0, then A is an m × m matrix over K . Since it is in -weak Popov, it is invertible and lower triangular, hence det( A ) is the productof its diagonal entries, which is computed in O ( m ) multiplications in K . This base case is notnecessary for obtaining the correctness and the cost bound in Proposition 3.3; still, not using itwould incur a cost of O ( m ω ) operations in the case D = lgorithm 2 D eterminant O f W eak P opov ( A ) Input: a matrix A ∈ K [ x ] m × m in -weak Popov form. Output: the determinant of A , up to multiplication by some element of K \ { } . s = ( s , . . . , s m ) ∈ Z m ≥ ← rdeg ( A ); D ← s + · · · + s m (cid:46) D = degree of det( A ) if m = then return the polynomial f such that A = [ f ] (cid:46) base case: × matrix if D = then return the product of diagonal entries of A (cid:46) base case: matrix over K if D < m then (cid:46) handle constant rows to reduce to dimension ≤ D B ← A lm ( A ) − from which rows and columns with indices in { i | s i = } are removed return D eterminant O f W eak P opov ( B ) if s + · · · + s (cid:98) m / (cid:99) > D / then return D eterminant O f W eak P opov ( J m A lm ( A ) − J m ) Write A = [ A A A A ], with A of size (cid:98) m / (cid:99) × (cid:98) m / (cid:99) and A of size (cid:100) m / (cid:101) × (cid:100) m / (cid:101) K ∈ K [ x ] (cid:100) m / (cid:101)× m ← K ernel B asis ([ A A ] , s ) [ K K ] ∈ K [ x ] (cid:100) m / (cid:101)× m ← R educed T o W eak P opov ( K , s ), where K is (cid:100) m / (cid:101) × (cid:100) m / (cid:101) B ← U nbalanced M ultiplication ([ K K ] , [ A A ] , D ) (cid:46) B = K A + K A ∆ ← D eterminant O f W eak P opov ( B ) (cid:46) first recursive call ∆ ← D eterminant O f W eak P opov ( A ) (cid:46) second recursive call P ← W eak P opov T o P opov ( J (cid:100) m / (cid:101) K T J (cid:100) m / (cid:101) , rdeg ( s (cid:98) m / (cid:99) + ,..., s m ) ( K ) J (cid:100) m / (cid:101) ) ∆ ← D eterminant O f W eak P opov ( P T ) (cid:46) third recursive call return ∆ ∆ / ∆ For the recursion, we proceed inductively: we assume that the algorithm correctly computesthe determinant for all -weak Popov matrices of dimension less than m , and based on this weshow that it is also correct for any -weak Popov matrix A of dimension m . Case 1: D < m . Then A has at least one constant row; using linear algebra we reduce to thecase of a matrix B of dimension at most D with all rows of degree at least 1. Since s = rdeg ( A ),we can write A = X s lm ( A ) + R for a matrix R ∈ K [ x ] m × m such that rdeg ( R ) < s . Since A is -reduced, lm ( A ) is invertible and A lm ( A ) − = X s + R lm ( A ) − with rdeg ( R lm ( A ) − ) < s ,which implies lm ( A lm ( A ) − ) = I m . In particular, A lm ( A ) − is in -weak Popov form, and foreach i such that the row A i , ∗ is constant, i.e. s i =
0, both the i th row and i th column of A lm ( A ) − are the i th row and i th column of the identity matrix. Therefore the matrix B at Line 5 is in -weakPopov form, has the same determinant as A up to a constant, and has dimension { i | s i (cid:44) } ≤ D .Hence the correctness in this case. In terms of complexity, computing B essentially amounts tocomputing the product A lm ( A ) − , which is done by row-wise expanding A into a ( m + D ) × m matrix over K , right-multiplying by lm ( A ) − , and compressing the result back into a polynomialmatrix: this costs O ( m ω (1 + D / m )) ⊆ O ( m ω ) operations. Case 2: s + · · · + s (cid:98) m / (cid:99) > D /
2. Then we modify the input A so as to reduce to Case 3 . Aswe have seen above, lm ( A lm ( A ) − ) = I m . We now reverse the diagonal entries by reversing theorder of rows and columns: let B = J m A lm ( A ) − J m . Then lm ( B ) = J m lm ( A lm ( A ) − ) J m = I m ,hence B is in -weak Popov form: Line 7 calls the algorithm on this matrix to obtain det( B ) upto a constant, and this yields det( A ) since it is equal to det(lm ( A )) det( B ). To conclude the proofof correctness in that case (assuming correctness in Case 3 ), it remains to observe that B has thesame matrix dimension m as A , and that the matrix B has degrees such that calling the algorithmwith input B does not enter Case 2 but
Case 3 . Indeed, we have rdeg ( B ) = s J m , hence the sumof the first (cid:98) m / (cid:99) entries of the tuple rdeg ( B ) is s m + · · · + s (cid:100) m / (cid:101) + = D − ( s + · · · + s (cid:100) m / (cid:101) ),which is at most D / A lm ( A ) − , which costs O ( m ω (1 + D / m )) operations as we have seen above; this is in O ( m ω M (cid:48) ( D / m )). Case 3: s + · · · + s (cid:98) m / (cid:99) ≤ D /
2. Then, Line 7 performs no action, and Line 8 definessubmatrices of A . By construction, [ A A ] has full column rank and s ≥ rdeg ([ A A ]) holds. Thus,according to Lemma 2.10, Line 9 uses O ( m ω M (cid:48) ( D / m )) operations to compute an s -reduced basis K of K ([ A A ]), with | rdeg s ( K ) | ≤ D . Then, Theorem 1.3 states that Line 10 transforms K intoan s -weak Popov basis [ K K ] of this kernel at a cost of O ( m ω (1 + D / m )) operations, since | rdeg s ( K ) | ≤ D and min( s ) ≥
0. Since all s -reduced bases of K ( F ) have the same s -row degree upto permutation, | rdeg s ([ K K ]) | ≤ D holds, hence the assumptions of Lemma 2.8 are satisfiedand Line 11 uses O ( m ω M ( D / m )) operations to compute B = K A + K A .The important observation at this stage is the identity (cid:34) I (cid:100) m / (cid:101) K (cid:35) (cid:34) A A A A (cid:35) = (cid:34) A A (cid:35) (2)which, provided that K is nonsingular, implies det( A ) = det( B ) det( A ) / det( K ). We are goingto show that this is the formula used in Line 16 to compute det( A ).First, A has dimension less than m and, being a principal submatrix of the -weak Popovmatrix A , it is also in -weak Popov form. Hence the recursive call at Line 13 is sound and ∆ isequal to det( A ) up to a constant.Since A is in -weak Popov form and [ I (cid:100) m / (cid:101) K ] is in rdeg ( A )-weak Popov form, their product[ A A ] is in -weak Popov form; indeed lm ([ A A ]) is invertible and lower triangular accordingto Lemma 2.2. It follows that B is in -weak Popov form and has dimension less than m : Line 12recursively computes ∆ , equal to det( B ) up to a constant.It remains to prove that ∆ computed at Lines 14 and 15 is equal to det( K ) up to a constant.Let v = rdeg ( A ) = ( s (cid:98) m / (cid:99) + , . . . , s m ) be the shift used at Line 14, and let d = rdeg v ( K ) = rdeg s ([ K K ]). Applying Lemma 2.7 (with F = A , G = A , t = , and v as above) showsthat [ K K ] has all its s -pivot entries in K , and in particular K is in v -weak Popov form. Let δ ∈ Z n ≥ be the v -pivot degree of K , where n = (cid:100) m / (cid:101) , and note that d = δ + v ≥ δ since v ≥ .Then, Lemma 3.2 states that J n K T J n is in − d J n -weak Popov form; its − d J n -pivot degree is thelist of degrees of its diagonal entries, that is, δ J n . Since d J n ≥ δ J n , we can apply Theorem 1.4,which implies that Line 14 computes the − d J n -Popov form P of J n K T J n using O ( m ω M (cid:48) ( | d | / m ))operations; as we have seen above, | d | = | rdeg s ([ K K ]) | ≤ D . Then, from the last itemof Lemma 3.2, P T is in -weak Popov form and det( K ) = det(lm v ( K )) det( P T ), hence Line 15correctly computes det( K ) up to a constant.To conclude this presentation of our determinant algorithm, we note that it would be benefi-cial, in a practical implementation, to add an early exit. Precisely, just after computing det( B ) atLine 12, one could perform the following action before proceeding to the next steps: if deg( ∆ ) = D then return ∆ (cid:46) early exit Indeed, recall that ∆ is det( C ) up to a constant: for a generic matrix A , we have deg( ∆ ) = D and det( B ) is det( A ) up to a constant, so it is not necessary to perform the remaining recursivecalls of the algorithm (see also [21, Sec. 4.2.2]). Let us prove this claim: if deg( ∆ ) = D , thendet( B ) is det( A ) up to a constant. Since [ K K ] is a kernel basis, it has unimodular column bases1822, Lem. 2.2], and thus it can be completed into a unimodular matrix U = [ U U K K ] ∈ K [ x ] m × m [56, Lem. 2.10]. Therefore UA = (cid:34) U U K K (cid:35) (cid:34) A A A A (cid:35) = (cid:34) B B (cid:35) where [ B B ] = [ U U ] A . Since det( U ) is in K \ { } , det( A ) is equal to det( B ) det( B ) up to aconstant. Now, our assumption deg(det( B )) = D = deg(det( A )) implies that det( B ) is in K \ { } ,hence the correctness of this early exit. We have seen above that all computations in Algorithm 2 other than recursive calls have anarithmetic cost in O ( m ω M (cid:48) ( D / m )) operations in K ; here, we complete the proof of the cost boundin Proposition 5.14. In this section, we use the assumptions H sl , H sm , and H ω as well as theirconsequences stated in Section 1.1.Let C ( m , D ) denote the arithmetic cost of Algorithm 2; recall that D is the degree of thedeterminant of the input, which is also the sum of its row degrees. First consider the case m ≤ D .If s + · · · + s (cid:98) m / (cid:99) > D /
2, the reduction to the case s + · · · + s (cid:98) m / (cid:99) ≤ D / m and D performed at Line 7 costs O ( m ω (1 + D / m )). Once we are in the latter case, there are threerecursive calls with input matrices having the following dimensions and degrees: • At Line 12, the matrix B is (cid:100) m / (cid:101) × (cid:100) m / (cid:101) , and applying the predictable degree propertyon Eq. (2) gives in particular rdeg ( B ) = rdeg s ([ K K ]), hence | rdeg ( B ) | ≤ D . • At Line 13, the matrix A is (cid:98) m / (cid:99)×(cid:98) m / (cid:99) and the sum of its row degrees is s + · · · + s (cid:98) m / (cid:99) ,which is at most D / • At Line 15, the matrix P T is (cid:100) m / (cid:101)×(cid:100) m / (cid:101) and its -pivot degree is rdeg ( P T ) = δ J m . Recallindeed that this is the list of diagonal degrees of P T , which is the same as that of P , andthus the same as that of J n K T J n according to Lemma 2.5. Now, from | δ + v | = | δ | + | v | ≤ D and the assumption | v | = s (cid:98) m / (cid:99) + + . . . + s m > D /
2, we obtain | rdeg ( P T ) | = | δ | ≤ D / m is a power of 2. If it is not, a given input matrixcan be padded with zeros, and ones on the main diagonal, so as to form a square matrix withdimension the next power of two and the same determinant. According to the three items above,the cost bound then satisfies: C ( m , D ) ≤ C ( m / , (cid:98) D / (cid:99) ) + C ( m / , D ) + O (cid:0) m ω M (cid:48) ( D / m ) (cid:1) . Letting the O ( · ) term aside, we illustrate this recurrence relation in Fig. 1.Let µ = log ( m ) and let K be the constant of the O ( · ) term above. Recalling that m ≤ D and (cid:98)(cid:98) D / j (cid:99) / (cid:99) = (cid:98) D / j + (cid:99) , unrolling this recurrence to the i th recursion level for 0 ≤ i ≤ µ yields C ( m , D ) ≤ i (cid:88) j = a i , j C (cid:18) m i , (cid:22) D j (cid:23)(cid:19) + K i − (cid:88) k = k (cid:88) j = a k , j (cid:18) m k (cid:19) ω M (cid:48) (cid:32) D / j m / k (cid:33) ≤ i (cid:88) j = a i , j C (cid:18) m i , (cid:22) D j (cid:23)(cid:19) + Km ω M (cid:48) (cid:18) Dm (cid:19) i − (cid:88) k = k (cid:88) j = a k , j − k ω M (cid:48) (cid:16) k − j (cid:17) , , Dm / , D m / , (cid:98) D / (cid:99) m / , D m / , (cid:98) D / (cid:99) m / , (cid:98) D / (cid:99) m / , D m / , (cid:98) D / (cid:99) m / , (cid:98) D / (cid:99) m / , (cid:98) D / (cid:99) , D . . . , (cid:98) D / j (cid:99) . . . , (cid:98) D / µ (cid:99)
11 21 4 41 6 12 81 2 j (cid:16) µ j (cid:17) µ Figure 1: Directed acyclic graph of recursive calls, of depth µ = log ( m ). Each boxed node shows the matrix dimensionsand the determinantal degree of a recursive call. Beginning with one call in dimension and determinantal degree ( m , D ),for a given node the number above it indicates the number of times a recursive call corresponding to this node is made,and the numbers of recursive sub-calls this node generates are indicated on both arrows starting from this node. where the last inequality comes from the submultiplicativity of M (cid:48) ( · ) and the coe ffi cients a i , j satisfy a i , = , a i , i = i , a i , j = a i − , j + a i − , j − for 0 < j < i . In Fig. 1, one can observe the similarity between Pascal’s triangle and the number of calls withparameters ( m / i , D / j ). This translates as a connection between the a i , j ’s and the binomialcoe ffi cients: one can prove by induction that a i , j = j (cid:16) ij (cid:17) .Now, by assumption M (cid:48) ( d ) ∈ O ( d ω − − ε ) for some ε >
0; for such an ε >
0, let ˜ K be a constantsuch that M (cid:48) ( d ) ≤ ˜ Kd ω − − ε for all d ≥
0. Then i − (cid:88) k = k (cid:88) j = a k , j − k ω M (cid:48) (cid:16) k − j (cid:17) ≤ ˜ K i − (cid:88) k = − (1 + ε ) k k (cid:88) j = a k , j j ( ω − − ε ) ≤ ˜ K i − (cid:88) k = − (1 + ε ) k k (cid:88) j = (cid:32) kj (cid:33) = ˜ K i − (cid:88) k = − ε k and, defining the constant ˆ K = K ˜ K (cid:80) + ∞ k = − ε k , for i = µ we obtain C ( m , D ) ≤ µ (cid:88) j = a µ, j C (cid:18) , (cid:22) D j (cid:23)(cid:19) + ˆ Km ω M (cid:48) (cid:18) Dm (cid:19) . As we have seen above, parameters (1 , d ) for any d ∈ Z ≥ correspond to base cases with a 1 × O (1) operations20ach; then the total cost of these base cases is bounded asymptotically by (cid:80) µ j = a µ, j ≤ m log (3) ).Thus we obtain C ( m , D ) = O ( m ω M (cid:48) ( D / m )), under the assumption m ≤ D .For the case D < m handled at Line 4, Section 3.2 showed that C ( m , D ) = C ( ˆ m , D ) + O ( m ω )where ˆ m is the number of non-constant rows of the input matrix. Since ˆ m ≤ D , our proof aboveshows that C ( ˆ m , D ) is in O ( ˆ m ω M (cid:48) ( D / ˆ m )). Now our assumptions on M ( · ) imply in particularthat M (cid:48) ( · ) is subquadratic, hence this bound is in O ( ˆ m ω ( D / ˆ m ) ) = O ( ˆ m ω − D ) ⊆ O ( m ω ). Thisconcludes the complexity analysis.
4. Shifted forms: from reduced to weak Popov
This section proves Theorem 1.3 by generalizing the approach of [42, Sec. 2 and 3], whichfocuses on the non-shifted case s = . It first uses Gaussian elimination on the -leading matrixof A to find a unimodular matrix U such that UA is in -weak Popov form, and then exploits thespecific form of U to compute UA e ffi ciently. Here we extend this approach to arbitrary shiftsand show how to take into account the possible unbalancedness of the row degree of A .First, we generalize [42, Lem. 8] from s = to an arbitrary s , by describing how U can beobtained by computing a -weak Popov form of the s -leading matrix of A . Lemma 4.1.
Let s ∈ Z n and let A ∈ K [ x ] m × n be s -reduced with rdeg s ( A ) nondecreasing. Thereexists an invertible lower triangular T ∈ K m × m which can be computed in O ( m ω − n ) operationsin K and is such that T lm s ( A ) is in -unordered weak Popov form. For any such matrix T , • U = X t TX − t has polynomial entries and is unimodular, where t = rdeg s ( A ) , • rdeg s ( UA ) = rdeg s ( A ) and lm s ( UA ) = T lm s ( A ) , • UA is in s -unordered weak Popov form.Proof. Consider the matrix lm s ( A ) J n and its generalized Bruhat decomposition lm s ( A ) J n = CPR as defined in [35]: C ∈ K m × m is in column echelon form, R ∈ K m × n is in row echelon form, and P ∈ K m × m is a permutation matrix. Therefore J m RJ n is in -weak Popov form (see the paragraphbefore Lemma 2.4 in Section 2.4), and PRJ n is in -unordered weak Popov form. Since lm s ( A )has full row rank, C is lower triangular and invertible, hence PRJ n = C − lm s ( A ) J n which provesthe existence of T = C − . Computing the decomposition costs O ( m ω − n ) operations [13, Cor. 25]while inverting C costs O ( m ω ) operations. Alternatively, [42, Sec. 3] shows how to compute T within the same cost bound using an LUP decomposition with a modified pivoting strategy.For any such matrix T , write T = ( T i j ) i j and t = ( t i ) i . Then the entry ( i , j ) of U is T i j x t i − t j .Thus, U is lower triangular in K ( x ) m × m with diagonal entries in K \ { } , and since t = rdeg s ( A )is nondecreasing, U is a unimodular matrix in K [ x ] m × m .Now, consider the nonnegative shift u = t − (min( t ) , . . . , min( t )) ∈ Z m ≥ , and note that U = X u TX − u . (The introduction of u is to circumvent the fact that we have not defined the rowdegree of a matrix over the Laurent polynomials, a notion which would be needed if we used X t rather than X u in Eq. (4).) Since A is in s -reduced form, the predictable degree property yieldsrdeg s ( UA ) = rdeg t ( U ) = rdeg u ( U ) + (min( t ) , . . . , min( t )) . (3)On the other hand,rdeg u ( U ) = rdeg ( UX u ) = rdeg ( X u T ) = u = rdeg s ( A ) − (min( t ) , . . . , min( t )) , (4)where the third equality follows from the fact that T is a constant matrix with no zero row. Then,from Eqs. (3) and (4), we obtain rdeg s ( UA ) = rdeg s ( A ) = t .21hen, lm s ( UA ) is formed by the coe ffi cients of nonnegative degree of X − rdeg s ( UA ) UAX s = X − t UAX s = TX − t AX s . Since T is constant and lm s ( A ) is formed by the coe ffi cients of nonnegative degree of X − t AX s ,we obtain that lm s ( UA ) = T lm s ( A ). The third item then directly follows from Lemma 2.4.Knowing T , and therefore U , the remaining di ffi culty is to e ffi ciently compute UA . For this,in the case s = , the approach in [42] has a cost bound which involves the maximum degree d = deg( A ) = max(rdeg ( A )) [42, Thm. 13], and uses the following steps: • first compute x d X − rdeg ( A ) UA = T ( X ( d ,..., d ) − rdeg ( A ) A ); • then scale the rows via the left multiplication by x − d X rdeg ( A ) .While the scaling does not use arithmetic operations, the first step asks to multiply the constantmatrix T by an m × n polynomial matrix of degree d : this costs O ( m ω − nd ) operations in K .Such a cost would not allow us to reach our target bound for the computation of characteristicpolynomials, since d may be too large, namely when A has unbalanced row degrees.Below, we refine the approach to better conform with the row degrees of A . This leads toan improvement of the above cost bound to O ( m ω − n (1 + D / m )) where D = | rdeg ( A ) | , thusinvolving the average row degree of A instead of its maximum degree d . For this, we follow thestrategy of splitting the rows of A into subsets having degrees in prescribed intervals; when theseintervals get further away from the average row degree of A , the corresponding subset containsa smaller number of rows. Earlier works using a similar strategy such as [53], [57, Sec. 2.6], and[22, Sec. 4], are not directly applicable to the problem here. Furthermore, we exploit the fact that A has nondecreasing row degree and that T has a lower triangular shape to avoid logarithmicfactors in the cost bound. This results in Algorithm 3. Algorithm 3 R educed T o W eak P opov ( A , s ) Input: a matrix A ∈ K [ x ] m × n , a shift s ∈ Z n such that A is in s -reduced form. Output: an s -weak Popov form of A . (cid:46) Step 1: Ensure nonnegative shift and nondecreasing s -row degree (cid:47) ˆ s ∈ Z n ≥ ← s − (min( s ) , . . . , min( s )) ( B , t ) ∈ K [ x ] m × n × Z m ≥ ← matrix B obtained from A by row permutation such that the tuple t = rdeg ˆ s ( B ) is nondecreasing (cid:46) Step 2: Compute the factor T in the unimodular transformation U = X t TX − t (cid:47) L ∈ K m × n ← lm ˆ s ( B ), that is, the entry ( i , j ) of L is the coe ffi cient of degree t i − s j of theentry ( i , j ) of B , where ˆ s = ( s , . . . , s n ) and t = ( t , . . . , t m ) T ∈ K m × m ← invertible lower triangular matrix such that TL is in -unordered weak Popovform (cid:46) can be computed via a generalized Bruhat decomposition, see Lemma 4.1 (cid:46) Step 3: Compute the product P = UB (cid:47) D ← t + · · · + t m ; K ← (cid:98) log ( mt m / D ) (cid:99) + (cid:46) K = min { k ∈ Z > | t m < k D / m } i ← i k ← max { i | t i < k D / m } for 1 ≤ k ≤ K (cid:46) = i < i ≤ · · · ≤ i K − < i K = m P ← zero matrix in K [ x ] m × n for k from 1 to K do R ← { i k − + , . . . , m } ; C ← { i k − + , . . . , i k } ; θ ← t i k = max( t C ) P R , ∗ ← P R , ∗ + X t R − ( θ,...,θ ) T R , C X ( θ,...,θ ) − t C B C , ∗ return the row permutation of P which has increasing ˆ s -pivot index22 roof of Theorem 1.3. The first step builds a nonnegative shift ˆ s which only di ff ers from s byan additive constant, and builds a matrix B which is a row permutation of A ; hence any ˆ s -weakPopov form of B is an s -weak Popov form of A . Since the ˆ s -row degree t of B is nondecreasing,the construction of T at Step 2 ensures that the matrix U = X t TX − t ∈ K [ x ] m × m is unimodularand such that UB is in ˆ s -unordered weak Popov form, according to Lemma 4.1. Thus, for thecorrectness, it remains to prove that the matrix P computed at Step 3 is P = UB .Using notation from the algorithm, define C k = { i k − + , . . . , i k } for 1 ≤ k ≤ K . The C k ’sare pairwise disjoint and such that { , . . . , m } = C ∪ · · · ∪ C K . Then, slicing the columns of T according to these sets of column indices, we obtain UB = X t TX − t B = (cid:88) ≤ k ≤ K X t T ∗ , C k X − t C k B C k , ∗ . Furthermore, since T is lower triangular, all rows of T ∗ , C k with index not in R k = { i k − + , . . . , m } are zero. Hence, more precisely, UB = (cid:88) ≤ k ≤ K (cid:34) i k − × n X t R k T R k , C k X − t C k B C k , ∗ (cid:35) . This formula corresponds to the slicing of the product P = X t TX − t B in blocks as follows: ∗∗ ∗ X t R k ∗∗ ∗∗ ∗∗ ∗ T R k , C k ∗ ∗ ∗∗ ∗ ∗ ∗ ∗ ∗ X − t C k ∗ ∗ ∗∗ B C k , ∗ ∗∗ . The for loop at Step 3 computes P by following this formula, hence P = UB . (Note indeedthat the scaling by ( θ, . . . , θ ) in the algorithm can be simplified away and is just there to avoidcomputing with Laurent polynomials.)Concerning the complexity bound, we first note that the quantity D defined in the algorithmis the same as that in Theorem 1.3; indeed, rdeg ˆ s ( B ) = rdeg s ( B ) − (min( s ) , . . . , min( s )), hence D = t + · · · + t m = | rdeg ˆ s ( B ) | = | rdeg s ( B ) | − m · min( s ) = | rdeg s ( A ) | − m · min( s ) . By Lemma 4.1, Step 2 uses O ( m ω − n ) operations. As for Step 3, its cost follows from bounds onthe cardinalities of C k and R k . Precisely, we have C k ⊆ R k , and min( R k ) > i k − implies that eachentry of the subtuple t R k is at least 2 k − D / m . Hence R k · k − D / m ≤ (cid:80) i ∈R k t i ≤ t + · · · + t m = D ,and C k ≤ R k ≤ m / k − .Then, in the product X t R k − ( θ,...,θ ) T R k , C k X ( θ,...,θ ) − t C k B C k , ∗ , the left multiplication by X t R k − ( θ,...,θ ) does not use arithmetic operations; the matrix T R k , C k is over K and has at most m / k − rows andat most m / k − columns; and the matrix X ( θ,...,θ ) − t C k B C k , ∗ is over K [ x ] and has n columns and atmost m / k − rows. Furthermore, the latter matrix has degree at most θ : indeed, ˆ s ≥ impliesthat t = rdeg ˆ s ( B ) ≥ rdeg ( B ), hence in particular t C k ≥ rdeg ( B C k , ∗ ). Recall that θ = t i k ≤ k D / m holds, by definition of θ and i k . From these bounds on the dimensions and the degrees of theinvolved matrices, and using the fact that m / k − ≤ m ≤ n , it follows that computing the product23 t R k − ( θ,...,θ ) T R k , C k X ( θ,...,θ ) − t C k B C k , ∗ uses O (( m / k − ) ω − n ( θ + ⊆ O (( m / k − ) ω − n (2 k D / m + K . Thus, since ω >
2, the cost of Step 3 is O (cid:88) ≤ k ≤ K (cid:18) m k − (cid:19) ω − n (cid:32) k Dm + (cid:33) ⊆ O m ω − nD (cid:88) ≤ k ≤ K k (2 − ω ) + m ω − n (cid:88) ≤ k ≤ K k (1 − ω ) ⊆ O (cid:16) m ω − nD + m ω − n (cid:17) .
5. Shifted forms: from weak Popov to Popov
This section culminates in Section 5.5 with the description of Algorithm W eak P opov T o P opov and a proof of Theorem 1.4. Based on [42, Lem. 14] (which extends to the shifted case), the resultin Theorem 1.4 can easily be used to solve the same problem in the rectangular case with an m × n matrix A ; while this is carried out in Algorithm W eak P opov T o P opov , it is out of the main scopeof this paper and thus for conciseness we only give a cost analysis in the case m = n .Our approach is to obtain the − s -Popov form of A from a shifted reduced kernel basis ofsome matrix F built from A . This fact is substantiated in Section 5.1, which also proves that wehave precise a priori information on the degrees of this sought kernel basis.A folklore method for kernel basis computation is to find an approximant basis at an ordersu ffi ciently large so that it contains a kernel basis as a submatrix. More precisely, assuming weknow a list of bounds s such that there exists a basis of K ( F ) with column degree bounded by s ,the following algorithm computes such a kernel basis which is furthermore − s -reduced: • γ ← cdeg s ( F ) + • M ← basis of A γ ( F ) in − s -reduced form; • return the submatrix of M formed by its rows with nonpositive − s -degree.The idea is that any row p of M is such that pF = mod X γ , and if it further satisfies cdeg( p ) ≤ s then cdeg( pF ) ≤ cdeg s ( F ) < γ , so that pF = holds. Here the complexity mainly depends on | s | and | cdeg s ( F ) | , quantities that are often large in which case the algorithm of Zhou et al. [57]is more e ffi cient. Nevertheless there are cases, such as the one arising in this section, where bothsums are controlled, and elaborating over this approach leads to an e ffi cient algorithm.In order to propose Algorithm K nown D egree K ernel B asis in Section 5.4, e ffi ciently com-puting the kernel basis using essentially a single call to PM-B asis , we transform the input intoone with a balanced shift and a balanced order. Section 5.2 deals with the shifts by describinga transformation of the input inspired from [45, Sec. 3], allowing us to reduce to the case of ashift s whose entries are balanced, meaning that max( s ) is not much larger than the average of s .Section 5.3 deals with balancing the order γ by performing the overlapping partial linearizationof [45, Sec. 2].For the latter transformation, as noted above, we assume that there exists a basis of the consid-ered kernel whose column degree is bounded by s , or equivalently that − s -reduced kernel baseshave nonpositive − s -row degree. On the other hand, for the first transformation we must ensurethat − s -reduced kernel bases have nonnegative − s -row degree. Thus our algorithm works underthe requirement that − s -reduced bases of K ( F ) have − s -row degree exactly , hence its name.This is a restriction compared to [55, Algo. 1] which only assumes that − s -reduced bases of K ( F )have nonpositive − s -row degree and has to perform several approximant basis computations torecover the whole kernel basis, as outlined in the introduction. In short, we have managed to ex-ploit the fact that we have better a priori knowledge of degrees in kernel bases than in the context24f [55], leading to a faster kernel computation which, when used in our determinant algorithm,brings no logarithmic factor in the complexity. We normalize the matrix A into its − s -Popov form P using a kernel basis computation, an ap-proach already used in a context similar to ours in the non-shifted case in [42, Lem. 19]. Roughly,this stems from the fact that the identity UA = P with a unimodular U can be rewritten as (cid:104) U P (cid:105) (cid:34) A − I m (cid:35) = ;for a well-chosen shift, one retrieves [ U P ] as a shifted reduced kernel basis. The next statementgives a choice of shift suited to our situation, and describes the degree profile of such kernelbases. The focus on the shift − δ comes from the fact that any − δ -reduced form R of A is only aconstant transformation away from being the − s -Popov form P of A [see 26, Lem. 4.1]. Lemma 5.1.
Let s ∈ Z m , let A ∈ K [ x ] m × m be in − s -weak Popov form, let δ ∈ Z m ≥ be the − s -pivot degree of A , and assume that s ≥ δ . Let R be a − δ -weak Popov form of A and let U be theunimodular matrix such that UA = R . Let further d = ( s − δ , δ ) ∈ Z m and F = (cid:104) A − I m (cid:105) ∈ K [ x ] m × m .Then, • the − d -pivot profile of [ U R ] is ( m + j , δ j ) ≤ j ≤ m , • [ U R ] is in − d -weak Popov form with rdeg − d ([ U R ]) = and cdeg([ U R ]) ≤ d , • [ U R ] is a basis of K ( F ) .Proof. First, we prove that R has − δ -pivot degree δ ; note that this implies rdeg − δ ( R ) = andcdeg( R ) = δ since R is in − δ -weak Popov form. Lemma 2.5 shows that the − s -Popov form P of A has the same − s -pivot degree as A , that is, δ . Hence, by definition of Popov forms,cdeg( P ) = δ . Then, [26, Lem. 4.1] states that P is also in − δ -Popov form. Since R is a − δ -weakPopov form of P , by Lemma 2.5 it has the same − δ -pivot degree as P , that is, δ .Now, by the predictable degree property and since rdeg − s ( A ) = − s + δ ,rdeg − s + δ ( U ) = rdeg − s ( UA ) = rdeg − s ( R ) ≤ rdeg − δ ( R ) = , where the inequality holds because − s ≤ − δ . Thus, by choice of d , the − d -pivot entries of [ U R ]are the − δ -pivot entries of its submatrix R ; this proves the first item.Then, the second item follows: the matrix [ U R ] is in − d -weak Popov form since its − d -pivot index is increasing; its − d -row degree is equal to the − δ -row degree of R which is ; andrdeg − d ([ U R ]) = implies cdeg([ U R ]) ≤ d by Lemma 2.1.Let [ K K ] ∈ K [ x ] m × m be a basis K ( F ) (it has m = m − m rows since F is 2 m × m and hasrank m ). Since UA = R , the rows of [ U R ] are in this kernel. As a result, there exists a matrix V ∈ K [ x ] m × m such that [ U R ] = V [ K K ]. In particular, U = VK , and since U is unimodular,this implies that V is unimodular as well. Thus, the basis [ K K ] is unimodularly equivalent to[ U R ], and the latter matrix is also a basis of K ( F ).One may note similarities with [27, Lem. 5.1], which is about changing the shift of reducedforms via kernel basis computations. Here, we consider A in − s -reduced form and are interestedin its − δ -reduced forms R : we change − s into − δ = − s + ( s − δ ), with a nonnegative di ff erence s − δ ≥ . Still, the above lemma does not follow from [27] since here the origin shift − s isnonpositive and thus we cannot directly incorporate it in the matrix F by considering AX − s .25he next corollary uses notation from Lemma 5.1 and shows how to obtain the − s -Popovform of A via a − d -reduced basis of the kernel of F . Then Sections 5.2 to 5.4 focus on thee ffi cient computation of such a kernel basis. Corollary 5.2.
Let [ ˆU ˆR ] ∈ K [ x ] m × m be a − d -reduced basis of K ( F ) . Then, ˆR is a − δ -reducedform of A , and P = (lm − δ ( ˆR )) − ˆR is the − s -Popov form of A .Proof. It su ffi ces to prove that ˆR is a − δ -reduced form of A ; then, the conclusion follows from[26, Lem. 4.1]. Both matrices [ ˆU ˆR ] and [ U R ] are − d -reduced and thus have the same − d -rowdegree up to permutation. Lemma 5.1 yields rdeg − d ([ U R ]) = , hence rdeg − d ([ ˆU ˆR ]) = . Inparticular, since − d = ( δ − s , − δ ), we have rdeg − δ ( ˆR ) ≤ .We conclude by applying Lemma 2.3 to the row space of A , which has rank m , and whichhas a basis R in − δ -reduced form with − δ -row degree . From ˆUA = ˆR , we obtain that the rowsof ˆR are in this row space. This identity also implies that ˆR has rank m , otherwise there wouldexist a nonzero vector in the left kernel of ˆR , which would also be in the left kernel of ˆU since A is nonsingular, hence it would be in the left kernel of the full row rank matrix [ ˆU ˆR ]. Theassumptions of the lemma are satisfied, and thus ˆR is a − δ -reduced form of A . Here, we show that the output column partial linearization, used previously in algorithms forapproximant bases and generalizations of them [45, Sec. 3], [54, 26], can be applied to kernelcomputations when the sought kernel basis has nonnegative shifted row degree. The main e ff ectof this transformation is to make the shift and output degrees more balanced, while preservingmost other properties of ( s , F ). The transformation itself, defined below, is essentially a column-wise x δ -adic expansion for a well-chosen integer parameter δ . Definition 5.3.
Let s = ( s , . . . , s m ) ∈ Z m ≥ , and let δ ∈ Z > . For ≤ j ≤ m, write s j = ( α j − δ + β j with α j = (cid:100) s j /δ (cid:101) and ≤ β j ≤ δ if s j (cid:44) , and α j = and β j = if s j = . Definem = α + · · · + α m , and the expansion-compression matrix E ∈ K [ x ] m × m as E = x δ ... x ( α − δ ... x δ ... x ( α m − δ . Define also s = ( δ, . . . , δ, β (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) α , . . . , δ, . . . , δ, β m (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) α m ) ∈ Z m ≥ . In this context, for a matrix K ∈ K [ x ] k × m , we define the column partial linearization of K as theunique matrix K ∈ K [ x ] k × m such that K = KE and all the columns of K whose index is not in { α + · · · + α j , ≤ j ≤ m } have degree less than δ .
26e use the notation in this definition in all of Section 5.2. More explicitly, K = [ K · · · K m ]where K j ∈ K [ x ] k × α j is the unique matrix such that K ∗ , j = K j x δ ... x ( α j − δ and the first α j − K j have degree less than δ .This construction originates from [45, Sec. 3], where it was designed for approximant basiscomputations, in order to make the shift more balanced at the cost of a small increase of thedimension; this is stated in Lemma 5.4. As mentioned above, several slightly di ff erent versionsof this column partial linearization have been given in the literature: each version requires someminor adaptations of the original construction to match the context. Here, in order to benefit fromproperties proved in [28, Sec. 5.1], we follow the construction in [28, Lem. 5.2]: one can checkthat Definition 5.3 is a specialization of the construction in that reference, for the shift − s ∈ Z m ≤ and taking the second parameter to be t = max( − s ) so that t = − s − max( − s ) + t = − s . Lemma 5.4.
The entries of s are in { , , . . . , δ } . Furthermore, if δ ≥ | s | / m, then m ≤ m ≤ m.Proof. The first remark is obvious. For the second one, note that 1 ≤ α j ≤ + s j /δ holds byconstruction, for 1 ≤ j ≤ m . Hence m ≤ m ≤ m + | s | /δ ≤ m .Importantly, this column partial linearization behaves well with respect to shifted row degreesand shifted reduced forms. Lemma 5.5.
Let K ∈ K [ x ] k × m with rdeg − s ( K ) ≥ , and let K ∈ K [ x ] k × m be its column partiallinearization. Then, row degrees are preserved: rdeg − s ( K ) = rdeg − s ( K ) . Furthermore, if K is in − s -weak Popov form, then K is in − s -weak Popov form.Proof. We show that this follows from the second item in [28, Lem. 5.2]; as noted above, theshift t = ( t , . . . , t m ) in that reference corresponds to − s = ( − s , . . . , − s m ) here. Let ( π i , δ i ) ≤ i ≤ k denote the − s -pivot profile of K . By definition of the − s -pivot index and degree, rdeg − s ( K ) = ( δ i − s π i ) ≤ i ≤ k , hence our assumption rdeg − s ( K ) ≥ means that δ i ≥ s π i for 1 ≤ i ≤ k . Thus, wecan apply [28, Lem. 5.2] to each row of K , from which we conclude that K has − s -pivot profile( α + · · · + α π i , δ i − s π i + β π i ) ≤ i ≤ k .First, since the entry of − s at index α + · · · + α π i is − β π i , this implies thatrdeg − s ( K ) = ( δ i − s π i + β π i − β π i ) ≤ i ≤ k = ( δ i − s π i ) ≤ i ≤ k = rdeg − s ( K ) . Furthermore, if K is in − s -weak Popov form, then ( π i ) ≤ i ≤ k is increasing, hence ( α + · · · + α π i ) ≤ i ≤ k is increasing as well and K is in − s -weak Popov form. Remark 5.6.
One may further note the following properties: • Writing lm − s ( K ) = [ L ∗ , · · · L ∗ , m ] ∈ K k × m , we have lm − s ( K ) = [ · · · ∗ , (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) α · · · · · · ∗ , m (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) α m ] ∈ K k × m . If K is in − s -Popov form, then K is in − s -Popov form.These properties are not used here, but for reference we provide a proof in Appendix A. We will also need properties for the converse operation, going from some matrix P ∈ K [ x ] m × m to its compression PE . Lemma 5.7.
Let P ∈ K [ x ] k × m have − s -pivot profile ( π i , δ i ) ≤ i ≤ k and assume π i = α + · · · + α j i for some j i ∈ Z > , for ≤ i ≤ k. Then, rdeg − s ( PE ) = rdeg − s ( P ) , and PE has − s -pivot profile ( j i , δ i + ( α j i − δ ) ≤ i ≤ k = ( j i , δ i + s j i − β j i ) ≤ i ≤ k . If P is in − s -weak Popov form, then PE is in − s -weak Popov form.Proof. The − s -pivot profile of PE is directly obtained by applying the first item in [28, Lem. 5.2]to each row of P . From this − s -pivot profile, we getrdeg − s ( PE ) = ( δ i + s j i − β j i − s j i ) ≤ i ≤ k = ( δ i − β j i ) ≤ i ≤ k . On the other hand, since the entry of − s at index α + · · · + α j i is − β j i , the − s -pivot profile of P yields rdeg − s ( P ) = ( δ i − β j i ) ≤ i ≤ k . Thus, rdeg − s ( PE ) = rdeg − s ( P ). Now, if P is in − s -weak Popovform, then ( α + · · · + α j i ) ≤ i ≤ k is increasing, which implies that ( j i ) ≤ i ≤ k is increasing. As a result, PE is in − s -weak Popov form, since ( j i ) ≤ i ≤ k is the − s -pivot index of PE .Our approach for computing the kernel of F is based on the fact that if K is a basis of K ( F )and K is its column partial linearization, then KF = KEF = . This identity shows that thekernel K ( EF ) contains the rows of K , so we may hope to recover K , and thus K = KE , from abasis of K ( EF ); the main advantage is that the latter basis is computed with the balanced shift − s . Note that a basis of K ( EF ) does not straightforwardly yield K , at least because this kernelalso contains K ( E ). In Lemma 5.8 we exhibit a basis S for K ( E ), and then in Lemma 5.9 weshow that K ( EF ) is generated by K and S . We also give properties which allow us, from a basisof K ( EF ), to easily recover a basis K of K ( F ) which has the sought form (see Lemma 5.11). Lemma 5.8.
The matrix S = diag( S , . . . , S m ) ∈ K [ x ] ( m − m ) × m , where S j = x δ − . . . . . . x δ − ∈ K [ x ] ( α j − × α j for ≤ j ≤ m, is the -Popov basis of the kernel K ( E ) . Furthermore, S is also in − s -Popov form,it has − s -row degree , and its − s -pivot profile is ( α + · · · + α j − + i , δ ) ≤ i <α j , ≤ j ≤ m .Proof. By construction, SE = and S is in -Popov form. Besides, S has rank m − m , which isthe rank of K ( E ) since E has rank m . Now, observe that there is no nonzero vector of degree lessthan δ in the left kernel of the vector [1 x δ · · · x ( α j − δ ] T , and thus there is no nonzero vectorof degree less than δ in K ( E ). It follows that S is a basis of K ( E ). Indeed, if K ∈ K [ x ] ( m − m ) × m isa basis of K ( E ) in -reduced form, then rdeg( K ) ≥ ( δ, . . . , δ ). Since SE = , we have S = UK for some nonsingular U ∈ K [ x ] ( m − m ) × ( m − m ) . By the predictable degree property,( δ, . . . , δ ) = rdeg( S ) = rdeg( UK ) = rdeg rdeg( K ) ( U ) ≥ rdeg( U ) + ( δ, . . . , δ ) , U is constant. Thus U is unimodular, and S is a basis of K ( E ).Now consider j such that α j >
1, and write t = ( − δ, . . . , − δ, − β j ) ∈ Z α j . Then, the definitionof α j and β j , notably the fact that − β j <
0, ensures that S j is in t -Popov form with t -pivot degree( δ, . . . , δ ) and t -pivot index (1 , . . . , α j − S follows. Lemma 5.9.
Let F ∈ K [ x ] m × n . The following properties hold.(i) Column degrees are preserved: cdeg s ( EF ) = cdeg s ( F ) .(ii) The kernels of F and EF are related by K ( EF ) E = K ( F ) .(iii) Let K ∈ K [ x ] k × m be a basis of K ( F ) (hence k = m − rank( F ) ), and let K ∈ K [ x ] k × m be anymatrix such that K = KE . Then, B = (cid:34) KS (cid:35) ∈ K [ x ] ( k + m − m ) × m (5) is a basis of K ( EF ) , where S is the matrix defined in Lemma 5.8.(iv) Assume that rdeg − s ( K ) ≥ and that K is the column partial linearization of K . Then, rdeg − s ( B ) = (rdeg − s ( K ) , ) , and if K is in − s -weak Popov form, then B is in − s -unorderedweak Popov form.Proof. (i) By definition, cdeg s ( EF ) = cdeg( X s EF ) and cdeg s ( F ) = cdeg( X s F ). Since the matrix X s E = X δ ... X ( α − δ X s ... X δ ... X ( α m − δ X sm has column degree s and contains X s as a subset of its rows, we get cdeg( X s EF ) = cdeg( X s F ).(ii) For any vector p ∈ K ( EF ), we have pEF = , which means pE ∈ K ( F ). Conversely,from any p ∈ K ( F ) we can construct q ∈ K [ x ] × m such that qE = p since E has the identity as asubset of its rows; then, qEF = pF = , which means q ∈ K ( EF ).(iii) We recall that, by Lemma 5.8, S is a basis of K ( E ). The rows of B are in K ( EF ), since SEF = and KEF = KF = . Now, we want to prove that any u ∈ K ( EF ) is a K [ x ]-linearcombination of the rows of B . By Item (ii), uE ∈ K ( F ); thus uE = vK = vKE for some v ∈ K [ x ] × k . Therefore u − vK ∈ K ( E ), and it follows that u − vK = wS for some w ∈ K [ x ] × m .This yields u = [ v w ] B .(iv) Since rdeg − s ( K ) ≥ , we can apply Lemma 5.5, which yields rdeg − s ( K ) = rdeg − s ( K ),hence rdeg − s ( B ) = (rdeg − s ( K ) , ). Now, if K is in − s -weak Popov form, then Lemma 5.5 showsthat K is in − s -weak Popov form and that all entries of its − s -pivot index are in { α + · · · + α i , ≤ i ≤ m } . Besides, by Lemma 5.8, S is in − s -Popov form with a − s -pivot index which is disjointfrom { α + · · · + α i , ≤ i ≤ m } . Hence B is in − s -unordered weak Popov form. Remark 5.10.
Similarly to Item (iv), one may observe that if K is in − s -reduced form, then B is in − s -reduced form; and that if K is in − s -Popov form, then B is in − s -Popov form up to rowpermutation. These points will not be used here; for reference a proof is given in Appendix A. K ( F ) bycomputing a − s -weak Popov basis of K ( EF ) and taking a submatrix of it. Lemma 5.11.
Let F ∈ K [ x ] m × n and let r = m − rank( F ) be the rank of K ( EF ) . Assume that − s -reduced bases of K ( F ) have nonnegative − s -row degree. Let Q ∈ K [ x ] r × m be a − s -weakPopov basis of K ( EF ) . Let P ∈ K [ x ] k × m be the submatrix of the rows of Q whose − s -pivot indexis in { α + · · · + α j , ≤ j ≤ m } . Then, PE is a − s -weak Popov basis of K ( F ) .Proof. We first prove that the number of rows of PE is the rank of K ( F ), that is, k = m − rank( F ).Indeed, by Item (iv) of Lemma 5.9, the − s -pivot index of the − s -Popov basis of K ( EF ) containsthe − s -pivot index of S . By Lemma 5.8, the latter is the tuple formed by the integers in the set { , . . . , m } \ { α + · · · + α j , ≤ j ≤ m } sorted in increasing order. Since P is the submatrix of therows of Q whose − s -pivot index is not in this set, P has k = r − ( m − m ) = m − rank( F ) rows.Now, by construction, P is in − s -weak Popov form and its − s -pivot index has entries in { α + · · · + α i , ≤ i ≤ m } . Thus we can apply Lemma 5.7; it ensures that PE is in − s -weak Popovform and that rdeg − s ( PE ) = rdeg − s ( P ).It remains to prove that PE is a basis of K ( F ). Let K ∈ K k × m be the − s -Popov basis of K ( F ), and let K ∈ K k × m be its column partial linearization. Let d = rdeg − s ( K ), which hasnonnegative entries by assumption. Then, according to Lemma 5.5, K is in − s -weak Popov form,and rdeg − s ( K ) = d . Then, we define B ∈ K ( k + m − m ) × m as in Eq. (5); by Item (iv) of Lemma 5.9,the matrix B is a − s -unordered weak Popov basis of K ( EF ).Then, Lemma 2.5 shows that Q has the same − s -pivot profile as the row permutation of B which is in − s -weak Popov form. Since P (resp. K ) is the submatrix of the rows of Q (resp. B )whose − s -pivot index is in { α + · · · + α j , ≤ j ≤ m } , and since both P and K are in − s -weakPopov form, it follows that P and K have the same − s -pivot profile. In particular, we haverdeg − s ( P ) = rdeg − s ( K ) = d , from which we get rdeg − s ( PE ) = d . According to Lemma 2.3, sincethe rows of PE are in K ( F ), this implies that PE is a basis of K ( F ). Now we apply the overlapping partial linearization of [45, Sec. 2], more precisely the versionin [28, Sec. 5.2] which supports arbitrary γ as showed in the definition below that we recall forcompleteness. In the next lemma, we will show that the sought kernel basis can be retrieved as asubmatrix of an approximant basis for the transformed problem. Definition 5.12 ([45, 28]) . Let γ = ( γ , . . . , γ n ) ∈ Z n > , let F ∈ K [ x ] m × n with cdeg( F ) < γ , andlet µ ∈ Z > . Then, for ≤ i ≤ n, let γ i = α i µ + β i with α i = (cid:108) γ i µ − (cid:109) and ≤ β i ≤ µ . Let alson = max( α − , + · · · + max( α n − , , and define L µ ( γ ) = ( γ , . . . , γ n ) ∈ Z n + n > , where γ i = (2 µ, . . . , µ, µ + β i ) ∈ Z α i > if α i > and γ i = γ i otherwise. Considering the ith columnof F , we write its x µ -adic representation as F ∗ , i = F (0) ∗ , i + F (1) ∗ , i x µ + · · · + F ( α i ) ∗ , i x α i µ where cdeg([ F (0) ∗ , i F (1) ∗ , i · · · F ( α i ) ∗ , i ]) < ( µ, . . . , µ, β i ) . If α i > , we define F ∗ , i = (cid:104) F (0) ∗ , i + F (1) ∗ , i x µ F (1) ∗ , i + F (2) ∗ , i x µ · · · F ( α i − ∗ , i + F ( α i ) ∗ , i x µ (cid:105) ∈ K [ x ] m × α i nd J i = [ α i − ] ∈ K [ x ] ( α i − × α i , and otherwise we let F ∗ , i = F ∗ , i and J i ∈ K [ x ] × . Then, L γ ,µ ( F ) = F ∗ , F ∗ , · · · F ∗ , n J J . . . J n ∈ K [ x ] ( m + n ) × ( n + n ) is called the overlapping linearization of F with respect to γ and µ . Lemma 5.13.
Let F ∈ K [ x ] m × n , and let s ∈ Z m ≥ be such that there exists a basis of K ( F ) withnonpositive − s -row degree. Let µ ∈ Z > with µ > max( s ) , t = ( s , µ − , . . . , µ − ∈ Z m + n ≥ , and γ ∈ Z n > with γ ≥ cdeg s ( F ) + . Let M ∈ K [ x ] ( m + n ) × ( m + n ) be a − t -reduced basis of A L µ ( γ ) ( L γ ,µ ( F )) .Then, exactly k = m − rank( F ) rows of M have nonpositive − t -degree, and the first m columns ofthese rows form a matrix K ∈ K [ x ] k × m which is a − s -reduced basis of K ( F ) . Besides, if M is in − t -weak Popov form, then K is in − s -weak Popov form.Proof. Let [
K Q ] be the submatrix of M formed by its rows of nonpositive − t -degree, where K ∈ K [ x ] k × m and Q ∈ K [ x ] k × n ; we have 0 ≤ k ≤ m + n . By choice of t , from rdeg − t ([ K Q ]) ≤ we get deg( Q ) < µ and rdeg − s ( K ) ≤ .In particular, deg( K ) ≤ max( s ) < µ : applying the second item in [28, Lem. 5.6] to eachrow of [ K Q ] shows that rdeg( Q ) < rdeg( K ) and that the rows of K are in A γ ( F ), that is, KF = mod X γ . On the other hand, cdeg( K ) ≤ s implies that cdeg( KF ) ≤ cdeg s ( F ) < γ , hence KF = , i.e. the rows of K are in K ( F ). This implies that the rank of K is at most the rank of themodule K ( F ), i.e. rank( K ) ≤ m − r where r = rank( F ).Furthermore, from rdeg( Q ) < rdeg( K ) and max( s ) < µ we obtainrdeg − s ( K ) ≥ rdeg( K ) − max( s ) > rdeg( Q ) − µ + = rdeg ( − µ + ,..., − µ + ( Q ) , hence by choice of t we have lm − t ([ K Q ]) = [lm − s ( K ) ]. Since this matrix is a subset ofthe rows of the nonsingular matrix lm − t ( M ), it has full row rank, and thus lm − s ( K ) has full rowrank. This shows that K is in − s -reduced form, and that k = rank( K ). If M is furthermorein − t -weak Popov form, then [ K Q ] is in − t -weak Popov form as well; then, the identitylm − t ([ K Q ]) = [lm − s ( K ) ] shows that the − t -pivot entries of [ K Q ] are all located in K ,hence K is in − s -weak Popov form.It remains to prove that k = m − r and that the rows of K generate K ( F ).By assumption, there exists a basis P ∈ K [ x ] ( m − r ) × m of K ( F ) such that cdeg( P ) ≤ s . Inparticular, the rows of P are in A γ ( F ), and applying the first item of [28, Lem. 5.6] to each ofthese rows shows that there exists a matrix R ∈ K [ x ] ( m − r ) × n with rdeg( R ) < rdeg( P ) and suchthat the rows of [ P R ] are in A L µ ( γ ) ( L γ ,µ ( F )). Thus, [ P R ] is a left multiple of M .A key remark now is that rdeg − t ([ P R ]) ≤ . Indeed, by Lemma 2.1 rdeg − s ( P ) ≤ followsfrom cdeg( P ) ≤ s , and we havedeg( R ) < deg( P ) = max(cdeg( P )) ≤ max( s ) < µ. Thus, since M is − t -reduced, the predictable degree property ensures that [ P R ] is a left multipleof M which does not involve the rows of M of positive − t -degree, i.e. a left multiple of [ K Q ].In particular, P is a left multiple of K : there exists a matrix V ∈ K [ x ] ( m − r ) × k such that VK = P .31ince P has rank m − r , we obtain k ≥ m − r , hence k = m − r . On the other hand, since therows of K are in K ( F ), there exists a matrix W ∈ K [ x ] k × k such that WP = K . It follows that P = VWP , and since P has full row rank this implies VW = I k . This means that P and K are leftunimodularly equivalent, hence K is a basis of K ( F ), which concludes the proof. After applying the transformations presented in Sections 5.2 and 5.3, we are left with thecomputation of an approximant basis for a balanced order L µ ( γ ) and a balanced shift − s : this ise ffi ciently by PM-Basis, designed in [21] as an improvement of [3, Algo. SPHPS]. Here, we usethe version in [28, Algo. 2] which ensures that the output basis is in − s -weak Popov form. Algorithm 4 K nown D egree K ernel B asis ( F , s ) Input: a matrix F ∈ K [ x ] m × n , and a nonnegative shift s ∈ Z m ≥ . Requirement: − s -reduced bases of K ( F ) have − s -row degree . Output: a − s -weak Popov basis of K ( F ). (cid:46) Step 1: Output column partial linearization (cid:47) δ ← (cid:100) D / m (cid:101) ∈ Z > , where D = max( | s | , | cdeg s ( F ) | , Apply Definition 5.3 to ( s , δ ) to obtain the parameters and expansion-compression matrix:( α , . . . , α m ) ∈ Z m > , m ∈ Z > , s ∈ Z m ≥ , E ∈ K [ x ] m × m (cid:46) Step 2: Overlapping partial linearization (cid:47) γ ← cdeg s ( F ) + ∈ Z n > (cid:46) order for approximation, equal to cdeg s ( EF ) + Apply Definition 5.12 to ( γ , EF , δ +
1) to obtain the order L δ + ( γ ) ∈ Z n + n > and the matrix L γ ,δ + ( EF ) ∈ K [ x ] ( m + n ) × ( n + n ) (cid:46) Step 3: Compute − t -weak Popov basis of A L δ + ( γ ) ( L γ ,δ + ( EF )) (cid:47) t ← ( s , δ, . . . , δ ) ∈ Z m + n ≥ Γ ← max( L δ + ( γ )); G = L γ ,δ + ( EF ) X ( Γ ,..., Γ ) −L δ + ( γ ) (cid:46) use uniform order ( Γ , . . . , Γ ) M ∈ K [ x ] ( m + n ) × ( m + n ) ← PM-B asis ( Γ , G , − t ) (cid:46) Step 4: Deduce first − s -weak Popov basis of K ( EF ) , then − s -weak Popov basis of K ( F ) (cid:47) Q ∈ K [ x ] k × m ← first m columns of the rows of M which have nonpositive − t -degree P ∈ K [ x ] k × m ← the rows of Q whose − s -pivot index is in { α + · · · + α j , ≤ j ≤ m } return PEProposition 5.14. Algorithm 4 is correct. Let D = max( | s | , | cdeg s ( F ) | , . Then, assuming m ≥ nand using notation from the algorithm, its cost is bounded by the sum of: • the cost of performing PM-B asis at order at most (cid:100) D / m (cid:101) + on an input matrix of rowdimension m + n ≤ m and column dimension n + n ≤ m; • O ( m ) extra operations in K .Thus, Algorithm 4 uses O ( m ω M (cid:48) ( D / m )) operations in K .Proof. Using the assumption that − s -reduced bases of K ( F ) have − s -row degree along withItem (iv) of Lemma 5.9 shows that the − s -reduced bases of K ( EF ) have − s -row degree . Thus,we can apply Lemma 5.13 to ( EF , s , µ, γ ) with µ = δ + > max( s ) and γ = cdeg s ( EF ) + γ = cdeg s ( F ) + M is a − t -weak Popovbasis of A ( Γ ,..., Γ ) ( G ) = A L δ + ( γ ) ( L γ ,δ + ( EF )) (see e.g. [28, Rmk. 3.3] for this approach to make theorder uniform). Then, Lemma 5.13 states that the matrix Q at Line 12 has k = m − rank( EF ) = − rank( F ) rows and is a − s -weak Popov basis of K ( EF ). Then, by Lemma 5.11, PE is a − s -weak Popov basis of K ( F ), hence the correctness.For the cost analysis, we assume m ≥ n , and we start by summarizing bounds on the dimen-sions and degrees at play. Lemma 5.4 yields m ≤ m ≤ m , while Item (i) of Lemma 5.9 ensurescdeg s ( EF ) = cdeg s ( F ). Each entry of L δ + ( γ ) is at most 2( δ + = (cid:100) D / m (cid:101) +
2, by construction.Writing γ = ( γ , . . . , γ n ), by Definition 5.12 we have n = max (cid:18)(cid:24) γ δ + − (cid:25) , (cid:19) + · · · + max (cid:18)(cid:24) γ n δ + − (cid:25) , (cid:19) ≤ γ δ + + · · · + γ n δ + = | γ | δ + . Since | γ | = | cdeg s ( F ) | + n ≤ D + m , and since δ + ≥ ( D + m ) / m , it follows that n ≤ m . Besides,rdeg − s ( P ) ≤ by construction, so that cdeg( P ) ≤ s by Lemma 2.1.The only steps that involve operations in K are the call to PM-B asis at Line 10 and themultiplication PE at Line 14. The construction of E and the inequality cdeg( P ) ≤ s imply thatthe product PE mainly involves concatenating vectors of coe ffi cients; concerning operations in K , there are m − m columns of P for which we may add the constant term of that column tothe term of degree δ of the previous column. Therefore Line 14 has cost bound O ( mk ); since m ≤ m and k ≤ m this is in O ( m ). At Line 10, we call the approximant basis subroutine PM-B asis discussed in Section 2.5; since Γ = max( L δ + ( γ )) is at most 2 (cid:100) D / m (cid:101) + ∈ O (1 + D / m ),Lemma 2.9 states that this call uses O (cid:16) ( m + n + n + n )( m + n ) ω − M (cid:48) ( D / m ) (cid:17) operations in K . Since n + n ≤ m + n ≤ m , this yields the claimed cost bound. Algorithm 5 W eak P opov T o P opov ( A , s ) Input: a matrix A ∈ K [ x ] m × n , a shift s ∈ Z n such that A is in − s -weak Popov form. Requirement: s π ≥ δ , where ( π , δ ) is the − s -pivot profile of A . Output: the − s -Popov form of A . (cid:46) Step 1: Find unimodular transformation and − δ -reduced form of A ∗ , π (cid:47) ( π , δ ) ← the − s -pivot profile of A [ U R ] ∈ K [ x ] m × m ← K nown D egree K ernel B asis (cid:32)(cid:34) A ∗ , π − I m (cid:35) , ( s π − δ , δ ) (cid:33) (cid:46) Step 2: Deduce − s -Popov form of A (cid:47) P ← zero matrix in K [ x ] m × n P ∗ , π ← lm − δ ( R ) − R P ∗ , { ,..., n }\ π ← lm − δ ( R ) − U A ∗ , { ,..., n }\ π return P For proving Theorem 1.4, we describe Algorithm W eak P opov T o P opov (Algorithm 5) and wefocus on the square case, m = n . Then, by definition of the − s -weak Popov form, δ is the tupleof degrees of the diagonal entries of A , and π = (1 , . . . , m ). Furthermore, in this case, s π = s , A ∗ , π = A , and P ∗ , π = P ; in particular, we can discard the step at Line 7 since the submatrices itinvolves are empty.First note that the shift d = ( s π − δ , δ ) = ( s − δ , δ ) used at Line 3 is nonnegative. Besides,Lemma 5.1 shows that − d -reduced bases of K ( F ) have − d -row degree . Thus the requirements33f Algorithm K nown D egree K ernel B asis are met, and Proposition 5.14 shows that the matrix[ U R ] computed at Line 3 is a − d -weak Popov basis of K ( F ). Then, Corollary 5.2 shows thatthe matrix P = P ∗ , π = lm − δ ( R ) − R computed at Line 6 is the − s -Popov form of A . This provesthat Algorithm 5 is correct.Concerning the cost bound, we focus on the case | s | >
0. Indeed, since s ≥ δ ≥ , if | s | = s = δ = . In this case, no computation needs to be done: the − s -Popov form of A is I m , theunique matrix in Popov form whose pivot degree is .The cost of Line 6 is that of multiplying lm − δ ( R ) − ∈ K m × m by R ∈ K [ x ] m × m , which hascolumn degree δ as explained in Section 5.1; this computation corresponds to the second item inTheorem 1.4. This multiplication can be done by first performing a column linearization of R into a m × ( m + | δ | ) matrix R over K , computing lm − δ ( R ) − R , and finally compressing the resultback into a polynomial matrix. This uses O ( m ω (1 + | δ | / m )) operations in K .Concerning Line 3, we rely on Proposition 5.14. Here, the matrix we give as input to Al-gorithm K nown D egree K ernel B asis has dimensions 2 m × m , hence the dimensions in Proposi-tion 5.14 satisfy m ≤ m and n ≤ m (the latter bound comes from the proof of that proposition).Then, Proposition 5.14 states that Line 3 costs: • O ( m ) operations in K , which is the third item in Theorem 1.4, • one call to PM-B asis on a matrix of row dimension m + n ≤ m , column dimension m + n ≤ m , and at order at most 2 (cid:100) D / (2 m ) (cid:101) +
2, where D = max( | d | , | cdeg d ( F ) | ,
1) and F is the input matrix [ A T − I m ] T .Besides, Proposition 5.14 also implies that Line 3 uses O ( m ω M (cid:48) ( D / m )) operations in K .We are going to prove that D = | s | , which concludes the proof. Indeed, the previous paragraphthen directly gives the overall cost bound O ( m ω M (cid:48) ( | s | / m )) in Theorem 1.4, and using2 (cid:24) D m (cid:25) + < (cid:18) D m + (cid:19) + = | s | / m + , the previous paragraph also gives the first item in that theorem.To observe that D = | s | , we first use the definition of d to obtain | d | = | s − δ | + | δ | = | s | .Since | s | ≥
1, this gives D = max( | s | , | cdeg d ( F ) | ). Now, cdeg d ( F ) is the entry-wise maximum ofcdeg s − δ ( A ) and cdeg δ ( − I m ) = δ . Since A is in − s -row reduced form with rdeg − s ( A ) = − s + δ ,Lemma 3.2 ensures that A T is in s − δ -reduced form with cdeg s − δ ( A ) = rdeg s − δ ( A T ) = s . Then,since s ≥ δ we obtain cdeg d ( F ) = s , and thus D = | s | . This concludes the proof of Theorem 1.4.As for the rectangular case m < n , the correctness of Algorithm 5 follows from the aboveproof in the square case, which shows that the algorithm correctly computes the − s π -Popov form P ∗ , π of A ∗ , π , and from [37, Lem. 5.1] concerning the computation of P ∗ , { ,..., n }\ π and the fact thatlm − δ ( R ) − U is the unimodular matrix which transforms A into P . The cost bound can be derivedusing the degree bounds on rectangular shifted Popov forms given in [5, 37] (we do not detailthis here since this would add technical material unrelated to the main results of this article). References [1] Abdeljaoued, J., Malaschonok, G.I., 2001. E ffi cient algorithms for computing the characteristic polynomial in adomain. Journal of Pure and Applied Algebra 156, 127–145. doi: .[2] Baur, W., Strassen, V., 1983. The complexity of partial derivatives. Theoretical computer science 22, 317–330.doi: .[3] Beckermann, B., Labahn, G., 1994. A uniform approach for the fast computation of matrix-type Pad´e approximants.SIAM J. Matrix Anal. Appl. 15, 804–823. doi: .
4] Beckermann, B., Labahn, G., Villard, G., 1999. Shifted normal forms of polynomial matrices, in: ISSAC’99,ACM. pp. 189–196. doi: .[5] Beckermann, B., Labahn, G., Villard, G., 2006. Normal forms for general polynomial matrices. J. SymbolicComput. 41, 708–737. doi: .[6] Berkowitz, S.J., 1984. On computing the determinant in small parallel time using a small number of processors.Information Processing Letters 18, 147–150. doi: .[7] Bunch, J.R., Hopcroft, J.E., 1974. Triangular factorization and inversion by fast matrix multiplication. Mathematicsof Computation 28, 231–236. doi: .[8] B¨urgisser, P., Clausen, M., Shokrollahi, A., 1997. Algebraic Complexity Theory. 1st ed., Springer-Verlag BerlinHeidelberg. doi: .[9] Cantor, D.G., Kaltofen, E., 1991. On fast multiplication of polynomials over arbitrary algebras. Acta Inform. 28,693–701. doi: .[10] Cook, S.A., 1966. On the minimum computation time of functions. Ph.D. thesis.[11] Csanky, L., 1975. Fast parallel matrix inversion algorithms, in: 16th Annual Symposium on Foundations ofComputer Science (sfcs 1975), pp. 11–12. doi: .[12] Danilevskij, A.M., 1937. The numerical solution of the secular equation. Matem. Sbornik 44, 169–171. In Russian.[13] Dumas, J.G., Pernet, C., Sultan, Z., 2017. Fast computation of the rank profile matrix and the generalized Bruhatdecomposition. J. Symbolic Comput. 83, 187–210. doi: .[14] Dumas, J.G., Pernet, C., Wan, Z., 2005. E ffi cient computation of the characteristic polynomial, in: ISSAC’05,ACM. pp. 140–147. doi: .[15] Dummit, D.S., Foote, R.M., 2004. Abstract Algebra. John Wiley & Sons.[16] Faddeev, D., Sominskii, I., 1949. Collected Problems in Higher Algebra, Problem n°979.[17] Forney, Jr., G.D., 1975. Minimal Bases of Rational Vector Spaces, with Applications to Multivariable LinearSystems. SIAM Journal on Control 13, 493–520. doi: .[18] Frame, J., 1949. A simple recurrent formula for inverting a matrix (abstract). Bull. of Amer. Math. Soc. 55, 1045.[19] Gathen, J.v.z., Gerhard, J., 2013. Modern Computer Algebra (third edition). Cambridge University Press. doi: .[20] Giesbrecht, M., 1995. Nearly optimal algorithms for canonical matrix forms. SIAM Journal on Computing 24,948–969. doi: .[21] Giorgi, P., Jeannerod, C.P., Villard, G., 2003. On the complexity of polynomial matrix computations, in: ISSAC’03,ACM. pp. 135–142. doi: .[22] Giorgi, P., Neiger, V., 2018. Certification of minimal approximant bases, in: ISSAC’18, ACM. pp. 167–174.doi: .[23] Gupta, S., Sarkar, S., Storjohann, A., Valeriote, J., 2012. Triangular x -basis decompositions and derandomizationof linear algebra algorithms over K [ x ]. J. Symbolic Comput. 47, 422–453. doi: .[24] Harvey, D., Van Der Hoeven, J., Lecerf, G., 2017. Faster polynomial multiplication over finite fields. J. ACM 63.doi: .[25] Ibarra, O.H., Moran, S., Hui, R., 1982. A generalization of the fast LUP matrix decomposition algorithm andapplications. Journal of Algorithms 3, 45–56. doi: .[26] Jeannerod, C.P., Neiger, V., Schost, E., Villard, G., 2016. Fast computation of minimal interpolation bases in Popovform for arbitrary shifts, in: ISSAC’16, ACM. pp. 295–302. doi: .[27] Jeannerod, C.P., Neiger, V., Schost, E., Villard, G., 2017. Computing minimal interpolation bases. J. SymbolicComput. 83, 272–314. doi: .[28] Jeannerod, C.P., Neiger, V., Villard, G., 2020. Fast computation of approximant bases in canonical form. J.Symbolic Comput. 98, 192–224. doi: .[29] Kailath, T., 1980. Linear Systems. Prentice-Hall.[30] Kaltofen, E., Villard, G., 2005. On the complexity of computing determinants. Computational Complexity 13,91–130. doi: .[31] Keller-Gehrig, W., 1985. Fast algorithms for the characteristic polynomial. Theoretical Computer Science 36,309–317. doi: .[32] Labahn, G., Neiger, V., Zhou, W., 2017. Fast, deterministic computation of the Hermite normal form and determi-nant of a polynomial matrix. J. Complexity 42, 44–71. doi: .[33] Le Gall, F., 2014. Powers of tensors and fast matrix multiplication, in: ISSAC’14, ACM. pp. 296–303. doi: .[34] Le Verrier, U., 1840. Sur les variations s´eculaires des ´el´ements elliptiques des sept plant`etes principales. Journaldes Math´ematiques Pures et Appliqu´ees 5, 220–254.[35] Manthey, W., Helmke, U., 2007. Bruhat canonical form for linear systems. Linear Algebra Appl. 425, 261–282.doi: .[36] Mulders, T., Storjohann, A., 2003. On lattice reduction for polynomial matrices. J. Symbolic Comput. 35, 377–401. oi: .[37] Neiger, V., Rosenkilde, J., Solomatov, G., 2018. Computing Popov and Hermite Forms of Rectangular PolynomialMatrices, in: ISSAC’18, ACM. pp. 295–302. doi: .[38] Neiger, V., Vu, T.X., 2017. Computing canonical bases of modules of univariate relations, in: ISSAC’17, ACM.pp. 357–364. doi: .[39] Pernet, C., Storjohann, A., 2007. Faster Algorithms for the Characteristic Polynomial, in: ISSAC’07, ACM. pp.307–314. doi: .[40] Popov, V.M., 1972. Invariant description of linear, time-invariant controllable systems. SIAM Journal on Control10, 252–264. doi: .[41] Samuelson, P.A., 1942. A method of determining explicitly the coe ffi cients of the characteristic equation. Annalsof Mathematical Statistics 13, 424–429.[42] Sarkar, S., Storjohann, A., 2011. Normalization of row reduced matrices, in: ISSAC’11, ACM. pp. 297–304.doi: .[43] Souriau, J.M., 1948. Une m´ethode pour la d´ecomposition spectrale et l’inversion des matrices. Comptes-Rendusde l’Acad´emie des Sciences 227, 1010–1011.[44] Storjohann, A., 2003. High-order lifting and integrality certification. J. Symbolic Comput. 36, 613–648. doi: .[45] Storjohann, A., 2006. Notes on computing minimal approximant bases, in: Challenges in Symbolic ComputationSoftware. URL: http://drops.dagstuhl.de/opus/volltexte/2006/776 .[46] Strassen, V., 1969. Gaussian elimination is not optimal. Numer. Math. 13, 354–356. doi: .[47] The FFLAS-FFPACK Group, 2019. FFLAS-FFPACK: Finite Field Linear Algebra Subroutines / Package, version2.4.3. http://github.com/linbox-team/fflas-ffpack .[48] The LinBox Group, 2019. Linbox: Linear algebra over black-box matrices, version 1.6.3. https://github.com/linbox-team/linbox/ .[49] Toom, A.L., 1963. The complexity of a scheme of functional elements realizing the multiplication of integers.Soviet Mathematics Doklady 3, 714–716.[50] Van Barel, M., Bultheel, A., 1992. A general module theoretic framework for vector M-Pad´e and matrix rationalinterpolation. Numer. Algorithms 3, 451–462. doi: .[51] Wolovich, W.A., 1974. Linear Multivariable Systems. volume 11 of
Applied Mathematical Sciences . Springer-Verlag New-York. doi: .[52] Zhou, W., 2012. Fast Order Basis and Kernel Basis Computation and Related Problems. Ph.D. thesis. Universityof Waterloo. URL: http://hdl.handle.net/10012/7326 .[53] Zhou, W., Labahn, G., 2009. E ffi cient computation of order bases, in: ISSAC’09, ACM. pp. 375–382. doi: .[54] Zhou, W., Labahn, G., 2012. E ffi cient algorithms for order basis computation. J. Symbolic Comput. 47, 793–819.doi: .[55] Zhou, W., Labahn, G., 2013. Computing column bases of polynomial matrices, in: ISSAC’13, ACM. pp. 379–386.doi: .[56] Zhou, W., Labahn, G., 2014. Unimodular completion of polynomial matrices, in: ISSAC’14, ACM. pp. 413–420.doi: .[57] Zhou, W., Labahn, G., Storjohann, A., 2012. Computing minimal nullspace bases, in: ISSAC’12, ACM. pp.366–373. doi: . Appendix A.
In this appendix, we give proofs for Remarks 5.6 and 5.10. We use notation from Section 5.2:a shift s ∈ Z m ≥ , an integer δ ∈ Z > , and a matrix F ∈ K [ x ] m × n are given; K ∈ K [ x ] k × m is abasis of K ( F ); ( α j , β j ) ≤ j ≤ m and s ∈ Z m ≥ and K ∈ K [ x ] k × m are as described in Definition 5.3.Following the context of Remarks 5.6 and 5.10, we assume rdeg − s ( K ) ≥ , hence in particularrdeg − s ( K ) = rdeg − s ( K ) (see Lemma 5.5). We start with the claims in Remark 5.6. Lemma A.1.
Writing lm − s ( K ) = [ L ∗ , · · · L ∗ , m ] ∈ K k × m , we have lm − s ( K ) = [ · · · ∗ , (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) α · · · · · · ∗ , m (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) α m ] ∈ K k × m . f K is in − s -Popov form, then K is in − s -Popov form.Proof. For given i ∈ { , . . . , k } and j ∈ { , . . . , m } , we rely on the definition of a leading matrix(see Section 2.3): the entry ( i , j ) of lm − s ( K ) is the coe ffi cient of degree d i + s j of K i , j , where d i = rdeg − s ( K i , ∗ ) = rdeg − s ( K i , ∗ ) ≥
0. First consider the case j (cid:60) { α + · · · + α π , ≤ π ≤ m } .Then the j th entry of s is s j = δ , and by definition of the output column partial linearization K i , j has degree less than δ , hence its coe ffi cient of degree d i + δ ≥ δ must be zero. This proves thatall columns of lm − s ( K ) with index not in { α + · · · + α π , ≤ π ≤ m } are indeed zero. It remainsto prove that in the case j = α + · · · + α π for some 1 ≤ π ≤ m , then the entry ( i , j ) of lm − s ( K )is equal to L i ,π . This holds, since in this case we have s j = β π , and by construction of K thecoe ffi cient of degree d i + β π of K i , j is equal to the coe ffi cient of degree d i + ( α π − δ + β π = d i + s π of K i ,π , which itself is equal to L i ,π by definition of a leading matrix.Now assume K is in − s -Popov form. We have showed in Lemma 5.5 that K is in − s -weakPopov form and that the row K i , ∗ has its − s -pivot at index α + · · · + α π i where π i is the − s -pivotindex of K i , ∗ . By construction, the column K ∗ ,α + ··· + α π i is the part of nonnegative degree of thecolumn x − ( α π i − δ K ∗ ,π i = x − s π i + β i K ∗ ,π i . It follows first that the − s -pivot entry K i ,α + ··· + α π i is monicsince it is a high degree part of the (monic) − s -pivot entry K i ,π i ; and second that K i ,α + ··· + α π i hasdegree strictly larger than all other entries in the column K ∗ ,α + ··· + α π i since K i ,π i has degree strictlylarger than all other entries in the column K ∗ ,π i . Hence K is in − s -Popov form.Now we prove the claims in Remark 5.10 concerning variants of Item (iv) of Lemma 5.9,using further notation from Section 5.2: S is the basis of K ( E ) described in Lemma 5.8 and B = [ KS ] ∈ K [ x ] ( k + m − m ) × m is the basis of K ( EF ) given in Eq. (5). As above we assume rdeg − s ( K ) ≥ ,and therefore rdeg − s ( B ) = (rdeg − s ( K ) , ) as stated in Item (iv) of Lemma 5.9. Lemma A.2.
With the above notation and assumptions, • If K is in − s -reduced form, then B is in − s -reduced form. • If K is in − s -Popov form, then B is in − s -Popov form up to row permutation.Proof. Suppose first that K is in − s -reduced form. We apply Lemma A.1:lm − s ( K ) = [ · · · ∗ , · · · · · · ∗ , m ]where lm − s ( K ) = [ L ∗ , · · · L ∗ , m ] ∈ K k × m . The latter matrix has full rank, since K is in − s -reduced form. On the other hand, for matrices such as the diagonal blocks of S we havelm ( − δ,..., − δ, − β ) x δ − . . . . . . x δ − = . . . . . . β ≥
1. As a result,lm − s ( B ) = (cid:34) lm − s ( K )lm − s ( S ) (cid:35) = | | L ∗ , · · · L ∗ , m | | . . . . . . . . . has full rank as well, which means that B is in − s -reduced form.Now, assume further that K is in − s -Popov form. Then Lemma A.1 states that K is in − s -Popov form, and the above description of lm − s ( B ) shows that K and S have disjoint − s -pivotindices, and that all their − s -pivots are monic. Hence B is in − s -unordered weak Popov formwith monic − s -pivots. It remains to show that each of these − s -pivots has degree strictly largerthan the other entries in the same column. Each row of the submatrix S of B has − s -pivot entry x δ and the other entries in the same column of B are either −
1, which has degree 0 < δ , or are ina column of K whose index is not in { α + · · · + α j , ≤ j ≤ m } , which has degree less than δ bydefinition of the column partial linearization. The j th row of the submatrix K of B has − s -pivotindex α + · · · + α j and strictly positive − s -pivot degree since rdeg − s ( K ) ≥ and the entry of − s at index α + · · · + α j is − β j <
0. This concludes the proof since the column S ∗ ,α + ··· + α j has degree0, and since K is in − ss