Algorithms for Computing Sparse Shifts for Multivariate Polynomials

Dima Yu. Grigoriev

IRMAR Université Rennes-1 Campus Beaulieu 35042 Rennes France

Lakshman Y. N.^{1}

Computing Science Research Bell Labs., 600 Mountain Hill Murray Hill, NJ 07974, USA

          e-mail: dima@maths.univ-rennes1.fr,  ynl@research.bell-labs.com

Abstract

In this paper, we investigate the problem of finding t-sparse shifts for multivariate polynomials. Given a polynomial f ∈ F[x_{1}, x_{2}, …, x_{n}] of degree d, and a positive integer t, we consider the problem of representing f(x) as a K-linear combination of the power products of u_{i} where u_{i} = x_{i} − b_{i} for some b_{i} ∈ K, an extension of F, for i = 1, …, n, i.e., f = ∑_{j} f_{j}u^{α_{j}}, in which at most t of the f_{j} are non-zero. We provide sufficient conditions for uniqueness of sparse shifts for multivariate polynomials, prove tight bounds on the degree of the polynomial being interpolated in terms of the sparsity bound t and a bound on the size of the coefficients of the polynomial in the standard representation, and describe two new efficient algorithms for computing sparse shifts for a multivariate polynomial.

1Work supported by NSF grant CCR–9203062

Introduction

In this paper, we consider the problem of computing t-sparse shifts for multivariate polynomials. Given a polynomial f ∈ F[x_{1}, x_{2}, …, x_{n}] of degree d (where F is a field of characteristic 0), consider the representation of f(x) as a K-linear combination of the power products of u_{i} where u_{i} = x_{i} − b_{i} for some b_{i} ∈ K, an extension of F, for i = 1, …, n, i.e., f = ∑_{j} f_{j}u^{α_{j}} where α_{j} denotes the multi-index (α_{j1}, α_{j2}, …, α_{jn}), and u^{α_{j}} indicates the power product u_{1}^{α_{j1}}u_{2}^{α_{j2}} … u_{n}^{α_{jn}}. Let t be a positive integer ≤ \binom{d+n}{n}. We say that b⃗ = (b_{1}, b_{2}, …, b_{n}) is a t-sparse shift for f (or, f is t-sparse in the shifted basis consisting of the power products of the u_{i}) if at most t of the f_{j} in the above representation are non-zero (the term “basis” refers to the fact that the power products of the u_{i} form a F-basis for the polynomial ring F[x_{1}, x_{2}, …, x_{n}]).

The main problem that we address is: given an f and t as above, can we efficiently compute a t-sparse shift for f if one exists? We are particularly interested of in the case polynomials that have rational shifts (each b_{i} ∈ F) and the case of polynomials that have finitely many t-sparse shifts.

Recently, there has been much interest in the design of efficient algorithms for computing sparse representations for various classes of functions such as polynomials, rational functions, and algebraic functions (Grigoriev-Karpinski 1987, Clausen et al 1988, Ben-Or&Tiwari 1988, Kaltofen-Lakshman 1988, Borodin-Tiwari 1990, Grigoriev-Karpinski-Singer 1990,1991,1992, 1993, 1994, Mansour 1992, Lakshman-Saunders 1993, 1994). The problem of finding sparsifying invertible linear tranformations for polynomials in F[x_{1}, x_{2}, …, x_{n}] was first addressed in a recent paper by Grigoriev and Karpinski (Grigoriev-Karpinski 1992) where they provide an algebraic criterion to be satisfied by any sparsifying linear tranformation for a polynomial in F[x_{1}, x_{2}, …, x_{n}] and an algorithm based on the algebraic criterion for computing sparsifying linear transformations. However, the algorithm requires solving systems of polynomial equations and inequalities involving the parameters a_{i,j}, c_{i} of the sparsifying transformations over the algebraic closure of F. While this is possible in principle, it is known to be very hard. For the important special case considered in this paper, we show that one can compute a system of polynomial equations involving the parameters of the unknown transformation which is already “solved” in a sense (the parameters are separated) in time that is polynomially bounded by t. The dependence of the algorithm on d, n is sensitive to how the polynomial f is presented. We will state the precise complexity result a little later. In this paper, we build on the results of two earlier papers, (Grigoriev-Karpinski 1993, Lakshman-Saunders 1994), and we make use of techniques used to deal with zero-dimensional Gröbner bases.

We assume that we are given a straight line program for computing f. Consequently, we can generate straight line programs for some low order partial derivatives of f efficiently (see Baur and Strassen, 1983). Instead of a straight-line program for f, it is enough if we have black boxes for f and certain low order partial derivatives of f. It is indeed possible to use this approach even if we have just a black box for computing f. We make some remarks later on as to how to modify our approach to work in this situation.

The main contributions of this paper are:

• sufficient conditions for uniqueness of sparse shifts for multivariate polynomials; • tight bounds on the degree of the polynomial being interpolated in terms of the sparsity

bound t and a bound M on the size of the coefficients of the polynomial in the standard representation.

• two new efficient algorithms for computing sparse shifts for a multivariate polynomial.

In section 1 we discuss conditions under which a polynomial can have a unique sparse shift. In section 2, we describe our first algorithm for computing sparse shifts. In section 3, we describe our degree bounds and the second algorithm for computing sparse shifts. The first algorithm computes t-sparse shifts for a multivariate polynomial f with finitely many t-sparse shifts in all cases except one – it can fail when deg_{x_{i}}(f) < t for two or more x_{i} and f still has finitely many sparse shifts. It performs (dt)^{O(n)} Q-operations if randomization is not allowed and t^{O(n)} Qoperations if randomization is allowed. When there is a unique shift, the algorithm performs (td^{n})^{O(1)} Q-operations if randomization is not allowed and (nt)^{O(1)} Q-operations if randomization is allowed. The second algorithm computes t-sparse shifts for a multivariate polynomial f without any finiteness restrictions on the number of t-sparse shifts. It has running time bounded by (nt)^{O(n^{2})}. We conclude with a discussion of open problems and possible applications in section 4.

Actually, the algorithms could run when F-operations are admitted, but for the complexity analysis to make the algorithms more realistic we allow just Q-operations.

1. Observations on the Uniqueness of Sparse Shifts

In (Lakshman-Saunders 1994), the following were shown to be true.

Theorem 1 Let f(x) ∈ F[x] be of degree d and let t ≤ (d + 1)/2. If there exists an α in some algebraic extension K of F such that f(x) is t-sparse in the shifted power basis 1, (x−α), (x−α)^{2}, …, then the shift α is unique. □

Corollary 1 Let f(x) ∈ F[x] be of degree d and let t ≤ (d + 1)/2. If α (in any extension K of F) is a t-sparse shift (hence the unique t-sparse shift) for f(x), then α ∈ F. □

The situation is more complicated than this for the multivariate case. The following example illustrates one difficulty. Consider the polynomial h(x, y) = x^{d}y + x^{d} with d > 4. h is 2-sparse with respect to the shift (0, c) for any c ∈ F (h(x, y) = x^{d}(y − c) + (c + 1)x^{d}). Obvious generalizations of this example to polynomials in F[x_{1}, x_{2}, …, x_{n}] lead us to the following conclusion. Let h ∈ F[x_{1}, x_{2}, …, x_{n}] and

B = {b⃗ = (b_{1}, b_{2}, …, b_{n}) ∈ F̄^{n} such that b⃗ is a t-sparse shift for h}.

Note that if b⃗ is a t-sparse shift for h, then b⃗ is a common zero of at least some \binom{d+n}{n} − t partial derivatives of h. It follows that B is an algebraic set in F̄^{n} and can be of any dimension from −1 to n − 1 (as usually we agree that the empty set has the dimension −1).

The next two weaker uniqueness results follow from theorem 1. Let h ∈ F[x_{1}, x_{2}, …, x_{n}] with deg_{x_{i}}(h) = d_{i} and let t < (d_{i} + 1)/2.

Lemma 1 If b⃗ = (b_{1}, b_{2}, …, b_{n}) and c⃗ = (c_{1}, c_{2}, …, c_{n}) are two t-sparse shifts for h, then b_{i} = c_{i} for any 1 ≤ i ≤ n.

Proof: For a fixed i consider a mapping Φ : F[x_{1}, x_{2}, …, x_{n}] ↦ F[x_{i}] with Φ(x_{j}) = a_{j} for i ≠ j where a_{j} ∈ F are chosen to preserve the degree in x_{i}, i.e., deg_{x_{i}}(h) = d_{i} = deg(Φ(h)). Both b_{i} and c_{i} are t-sparse shifts for Φ(h) and using theorem 1, we conclude that b_{i} = c_{i}. □

Lemma 2 Let d = min_{1≤i≤n}{deg_{x_{i}}(h)}. If t < (d + 1)/2, then h has at most one t-sparse shift.

Proof: Apply lemma 1 to each x_{i}. □

Stronger uniqueness results hold for the more general case of sparsifying linear transformations under side conditions.

1.1 Sufficient conditions for Uniqueness

In this subsection, we prove two different sufficient conditions for the uniqueness of sparsifying transformations. Let

f = ∑_{i=1}^{t} f_{i}x^{α_{i}} = ∑_{i=1}^{T} ϕ_{i}u^{β_{i}}, u⃗ = Ax⃗ + d⃗

where

u⃗ = ( \begin{matrix}u_{1} \\ u_{2} \\ … \\ u_{n} \end{matrix} ) ,  A = (a_{i,j})_{1≤i,j≤n},  x⃗ = ( \begin{matrix}x_{1} \\ x_{2} \\ … \\ x_{n} \end{matrix} ) ,  d⃗ = ( \begin{matrix}d_{1} \\ d_{2} \\ … \\ d_{n} \end{matrix} ) ,

a_{i,j}, d_{i} ∈ K and A non-singular. Let a_{1}, a_{2}, …, a_{n} denote the columns of A with A denoting the set {a_{1}, a_{2}, …, a_{n}}. Let d be the degree of f.

Theorem 2 If d⃗ is not in the span of B, for every proper subset B of A, then, t + T > d.

Proof: Let us substitute linear forms in a new variable ζ for the x_{i} as

x_{i} ⟶ k_{i}(ζ + N)

with k_{i}, N ∈ K to be chosen in a certain way. Clearly, the linear forms u_{i} become linear forms in ζ under this substitution. Let us try to choose k_{i}, N is such a way that each of the u_{i}-s becomes a scalar multiple of one and the same linear form in ζ, different from ζ + N. The u_{i} are tranformed as follows:

( \begin{matrix}a_{1,1} & a_{1,2} & … & a_{1,n} \\ a_{2,1} & a_{2,2} & … & a_{2,n} \\ ⋮ & ⋮ & ... & ⋮ \\ a_{n,1} & a_{n,2} & … & a_{n,n} \end{matrix} ) ( k_{1}(ζ + N) ; k_{2}(ζ + N) ; ⋮ ; k_{n}(ζ + N) ) + ( d_{1} ; d_{2} ; ⋮ ; d_{n} ) = ( r_{1} ; r_{2} ; ⋮ ; r_{n} ) ζ + ( Nr_{1} + d_{1} ; Nr_{2} + d_{2} ; ⋮ ; Nr_{n} + d_{n} )

where r_{i} = ∑_{j=1}^{n} a_{i,j}k_{j}. Let

r⃗ = ( \begin{matrix}r_{1} \\ r_{2} \\ … \\ r_{n} \end{matrix} ) , k⃗ = ( \begin{matrix}k_{1} \\ k_{2} \\ … \\ k_{n} \end{matrix} ) .

We want each r_{i}ζ + (r_{i}N + d_{i}) to be a scalar multiple of one and the same linear form in ζ which means we want λr⃗ = Nr⃗ + d⃗ for some non-zero λ ∈ K. So, we have (λ − N)r⃗ = d⃗. Recalling that r⃗ = Ak⃗, we have k⃗ = 1/(λ−N)A^{−1}d⃗. Let us choose λ, N such that λ, N ≠ 0 and λ ≠ N and find the corresponding k_{i} by solving the above system of equations. It follows from the hypothesis (that d⃗ is not in the span of any proper subset of A) that no k_{i} is zero. This implies that under the substitution considered, the degree does not decrease, i.e., deg(f) = deg(Φ(ζ)) where Φ(ζ) = f(x_{i} = k_{i}(ζ + N)). Now, Φ(ζ) has degree d, is t-sparse in the powers of ζ + N and T -sparse in the powers of ζ + λ. Invoking the univariate theorem, we have t + T > d. □

For the special case of sparse shifts, we get the following:

Corollary 2 Let f = ∑_{j=1}^{t} f_{j}x^{α_{j}} = ∑_{j=1}^{T} ϕ_{j}u^{β_{j}} where u_{i} = x_{i} + d_{i}. If no d_{i} = 0, then t + T > d.

Proof: Let A be the identity matrix in the above theorem. □

The second theorem on the number of sparse shifts imposes a different criterion. Let B = A^{−1} = (b_{i,j})_{1≤i,j≤n}. We have x⃗ = Bu⃗ + d⃗′ where d⃗′ = −Bd⃗. Note that

x_{i}∂/∂x_{i}(u_{1}^{s_{1}}u_{2}^{s_{2}} … u_{n}^{s_{n}}) = u_{1}^{s_{1}}u_{2}^{s_{2}} … u_{n}^{s_{n}}(a_{1,i}s_{1}/u_{1} + a_{2,i}s_{2}/u_{2} + … + a_{n,i}s_{n}/u_{n})(∑_{j=1}^{n} b_{i,j}u_{j}) + ϕ

where deg_{u}(ϕ) < s_{1} + s_{2} + … + s_{n}. Let d be the degree of f.

Theorem 3 If for some i, ∏_{j=1}^{n} a_{j,i} ∏_{j=1}^{n} b_{i,j} ≠ 0, then tT > d/n for n ≥ 0.

Assume the contrary. We set up some notation and prove supporting claims and then present the proof of the theorem. Let D = x_{i}∂/∂x_{i}. There exist ϵ_{0}, ϵ_{1}, …, ϵ_{t} ∈ F such that ∑_{0≤j≤t} ϵ_{j}D^{j}f = 0 (see lemma 4). Let t_{0} = max{j : ϵ_{j} ≠ 0} ≤ t.

Let S_{n}(d) ⊂ Z^{n} denote the set of integer points of the (n−1)-dimensional simplex {(z_{1}, z_{2}, …, z_{n}) : z_{i} ≥ 0, 1 ≤ i ≤ n, ∑_{1≤i≤n} z_{i} = d}. Call the set {β_{1}, …, β_{T} } the support of f. For any two vectors s, w ∈ S_{n}(d), let var(s, w) = ∑_{1≤i≤n} max(s_{i} − w_{i}, 0) = − ∑_{1≤i≤n} min(s_{i} − w_{i}, 0). For any vector v = (v_{1}, v_{2}, …, v_{n}) such that ∑_{1≤i≤n} v_{i} = 0, represent v uniquely as v = v^{(+)} + v^{(−)} where v_{i}^{(+)} = max(v_{i}, 0) and v_{i}^{(−)} = min(v_{i}, 0). We have var(v) = ∑_{1≤i≤n} v_{i}^{(+)} = − ∑_{1≤i≤n} v_{i}^{(−)}.

Let T_{0} be the number of points in the intersection of the support of f with S_{n}(d) and denote the points in this intersection by a_{1}, a_{2}, …, a_{T_{0}}. Observe that for any b ∈ Z^{n} satisfying the conditions var(b) = t_{0}, a_{j} + b ∈ S_{n}(d) (see above) the point a_{j} + b belongs to the support of D^{t_{0}}(u^{a^{j}}) (due to the assumption in the statement of the theorem) and does not belong to the support of D^{l}(u^{a^{j}}) for any l < t_{0}.

For each m, 1 ≤ m ≤ n, re-order the a_{i} in non-decreasing order on their m-th coordinates as 0 ≤ a_{1,m} ≤ a_{2,m} ≤ … ≤ a_{T_{0},m} ≤ d.

Lemma 3 There is an m, 1 ≤ m ≤ n, for which either a_{1,m} ≥ t_{0}, or, for a certain l, 1 ≤ l ≤ T_{0}−1, a_{l+1,m} − a_{l,m} ≥ 2t_{0} + 1.

Proof: Suppose not; then for every m we have a_{l,m} ≤ t_{0} − 1 + (l − 1)2t_{0} for each l, 1 ≤ l ≤ T_{0}. Hence,

∑_{1≤l≤T_{0}} a_{l,m} ≤ (t_{0} − 1)T_{0} + 2t_{0}T_{0}(T_{0} − 1)/2 = T_{0}(t_{0}T_{0} − 1) ≤ T_{0}(d/n − 1),

the latter inequality follows from the assumption that the theorem is wrong. On the other hand, ∑_{1≤m≤n} ∑_{1≤l≤T_{0}} a_{l,m} = dT_{0} since a_{l} ∈ S_{n}(d) for 1 ≤ l ≤ T_{0}, we get a contradiction. Hence, the lemma is true. □

Proof: (of theorem 3): Assume that tT ≤ d/n. Fix an m, 1 ≤ m ≤ n satisfying lemma 3.

First, consider the case a_{1,m} ≥ t_{0} (see the beginning of the proof of the theorem). Consider those points among a_{1}, …, a_{T_{0}} that belong to the (n − 2)-dimensional simplex S_{n}(d) ∩ {a_{j,m} = a_{1,m}}. Without loss of generality, assume that they are a_{1}, …, a_{T_{1}}, T_{1} ≤ T_{0}. Among these points, choose the one that is largest in the lexicographic order, without loss of generality, let it be a_{1}. Since a_{1,m} ≥ t_{0}, we have ∑_{k≠m} a_{1,k} ≤ d − t_{0}. Consider the point ā_{1} = (a_{1,1} + t_{0}, a_{1,2}, …, a_{1,m−1}, a_{1,m} − t_{0}, …, a_{1,n}) ∈ S_{n}(d). Clearly, var(a_{1}, ā_{1}) = t_{0}. We will show that for any point a_{j}, 2 ≤ j ≤ T_{0}, var(a_{j}, ā_{1}) > t_{0}, thereby showing that ā_{1} does not belong to the support of D^{l}(u^{a^{j}}) for all 0 ≤ l ≤ t_{0} and 2 ≤ j ≤ T_{0} and proving the theorem for the case a_{1,m} ≥ t_{0}. Suppose that var(a_{j}, ā_{1}) ≤ t_{0}; this means that v = a_{j} − ā_{1} = v^{(+)} + v^{(−)} with v^{(+)} = (0, …, 0, t_{0}, 0, …, 0) (the t_{0} appearing in the m-th position) and a_{j} ∈ {a_{1}, …, a_{T_{1}}}. Let a_{j} = (a_{j,1}, …, a_{j,n}). Since a_{j,1} ≤ a_{1,1}, (as a_{1} is lexicographically the largest among the a_{j}), v^{(−)} = (−t_{0}, 0, …, 0), but then a_{j} = a_{1}, a contradiction! This proves that var(a_{j}, ā_{1}) > t_{0}.

Now, consider the case a_{l+1,m} − a_{l,m} ≥ 2t_{0} + 1 for a certain l, 1 ≤ l T_{0} − 1 (according to lemma 3 we need just a trial of this and of the previous case). Among the points in {a_{1}, …, a_{T_{0}}} ∩ {a_{j,m} = a_{l+1,m}}, choose the largest in lexicographic order, call it a_{l+1}. As before, consider the point ā_{l+1} = a_{l+1} + (t_{0}, 0, …, 0, −t_{0}, 0, …, 0) ∈ S_{n}(d) (the −t_{0} is in position m) and prove that var(a_{j}, ā_{l+1}) > t_{0} for every j, 1 ≤ j ≤ T_{0} and j ≠ l + 1 and complete the proof as in the previous case. □

2. Computing Sparse Multivariate Shifts

We assume in our discussion that f has a finite (possibly zero) number of sparse shifts. If there are infinitely many t-sparse shifts with respect to one or more x_{i} ( when f is seen as an element of F[x_{1}, …, x_{i−1}, x_{i+1}, …, x_{n}][x_{i}] with t > deg_{x_{i}}(f)) but only finitely many of them can be combined to t-sparse shifts for the polynomial f, then the algorithm of this section can fail. This situation is taken care of by the algorithm of next section as the “low degree case”. In this section, we assume that t ≤ deg_{x_{i}}(f) for 1 ≤ i ≤ n. As before, let f = ∑_{j=1}^{t} f_{j}u^{α_{j}}. Recall that α_{j} denotes the multi-index (α_{j1}, α_{j2}, …, α_{jn}), and u^{α_{j}} indicates the power product u_{1}^{α_{j1}}u_{2}^{α_{j2}} … u_{n}^{α_{jn}}. Consider the ideal I ⊂ Q[y_{1}, y_{2}, …, y_{n}] which is the ideal of the points {α_{1}, α_{2}, …, α_{t}}. We shall construct a reduced Gröbner basis for I under any admissible term ordering. For convenience, we choose the lexicographic term ordering with y_{1} ≺ y_{2} ≺ … ≺ y_{n}. In fact, we can construct a triangular set decomposition of the ideal I instead of a Gröbner basis for I. For a good description of zero-dimensional Gröbner bases and triangular sets, we refer the reader to (Lazard 1992, Kapur–Lakshman 1992).

The Gröbner basis G for I under the chosen term ordering looks as follows:

G_{1}^{(1)}(y_{1}) = y_{1}^{δ_{1}} + g_{1,δ_{1}−1}^{(1)}y_{1}^{δ_{1}−1} + … + g_{1,1}^{(1)}y_{1} + g_{1,0}^{(1)},
G_{1}^{(2)}(y_{1}, y_{2}), G_{2}^{(2)}(y_{1}, y_{2}), …, G_{k_{2}}^{(2)}(y_{1}, y_{2}),
G_{1}^{(3)}(y_{1}, y_{2}, y_{3}), G_{2}^{(3)}(y_{1}, y_{2}, y_{3}), …, G_{k_{3}}^{(3)}(y_{1}, y_{2}, y_{3}),
⋮
G_{1}^{(n)}(y_{1}, y_{2}, …, y_{n}), G_{2}^{(n)}(y_{1}, y_{2}, …, y_{n}), …, G_{k_{n}}^{(n)}(y_{1}, y_{2}, …, y_{n}).

We recall a few standard terms and facts from the theory of Gröbner bases. For further details, see (Becker et al 1993).

• the head term of a polynomial h, HeadTerm(h) is the largest term (under the term ordering ≺) appearing in h with a non-zero coefficient.

• a term s is reduced with respect to the Gröbner basis G if it is not divisible by HeadTerm(g) for any g ∈ G.

• The number of reduced terms with respect to the Gröbner basis G is equal to the dimension of the residue class ring Q[y_{1}, …, y_{n}]/I as a Q-vector space and, in our case, is ≤ t.

Suppose h ∈ I and

h = y^{ϵ_{1}} + h_{2}y^{ϵ_{2}} + h_{3}y^{ϵ_{3}} + … + h_{L}y^{ϵ_{L}}

where y^{ϵ_{i}} = y_{1}^{ϵ_{i,1}}y_{2}^{ϵ_{i,2}} … y_{n}^{ϵ_{i,n}} and y^{ϵ_{i+1}} ≺ y^{ϵ_{i}} for 0 ≤ i < L. Let us define f_{ϵ_{i}} for a multi-index ϵ_{i} = (ϵ_{i,1}, ϵ_{i,2}, …, ϵ_{i,n}) with ϵ_{i,1}, ϵ_{i,2}, …, ϵ_{i,n} ≥ 0 as follows:

f_{(0,0,...,0)} = f,
f_{(ϵ_{i,1},ϵ_{i,2},...,ϵ_{i,n})} = (x_{1} − z_{1})∂f_{(ϵ_{i,1}−1,ϵ_{i,2},...,ϵ_{i,n})}/∂x_{1}  if ϵ_{i,1} > 0, else
= (x_{2} − z_{2})∂f_{(ϵ_{i,1},ϵ_{i,2}−1,...,ϵ_{i,n})}/∂x_{2}  if ϵ_{i,2} > 0, else
⋮
= (x_{n} − z_{n})∂f_{(ϵ_{i,1},ϵ_{i,2},...,ϵ_{i,n}−1)}/∂x_{n}  otherwise.

Furthermore, we extend the ordering ≺ to the set of polynomials f_{ϵ_{i}} as

f_{ϵ_{i}} ≺ f_{ϵ_{j}} iff y^{ϵ_{i}} ≺ y^{ϵ_{j}}.

Let

f̂_{ϵ_{i}} = f_{ϵ_{i}}(z_{1} = b_{1}, …, z_{n} = b_{n}, x_{1}, x_{2}, …, x_{n}) ∈ K[x_{1}, …, x_{n}]

where (b_{1}, b_{2}, …, b_{n}) is a t-sparse shift for f.

Lemma 4 Let h ∈ I and h = y^{ϵ_{1}} +h_{2}y^{ϵ_{2}} +h_{3}y^{ϵ_{3}} +…+h_{L}y^{ϵ_{L}}. Then the polynomials f̂_{ϵ_{1}}, f̂_{ϵ_{2}}, …, f̂_{ϵ_{L}} satisfy the Q-linear relation

h_{1}f̂_{ϵ_{1}} + h_{2}f̂_{ϵ_{2}} + h_{3}f̂_{ϵ_{3}} + … + h_{L}f̂_{ϵ_{L}} = 0, with h_{1} = 1.

Proof: We have

f̂_{ϵ_{i}} = ∑_{j=1}^{t} α_{j}^{ϵ_{i}}f_{j}(x − b)^{α_{j}}

where α_{j}^{ϵ_{i}} = α_{j,1}^{ϵ_{i,1}}α_{j,2}^{ϵ_{i,2}} … α_{j,n}^{ϵ_{i,n}},and, (x − b)^{α_{j}} = (x_{1} − b_{1})^{α_{j,1}}(x_{2} − b_{2})^{α_{j,2}} … (x_{n} − b_{n})^{α_{j,n}}. Therefore,

∑_{i=1}^{L} h_{i}f̂_{ϵ_{i}} = ∑_{i=1}^{L} h_{i}(∑_{j=1}^{t} α_{j}^{ϵ_{i}}f_{j}(x − b)^{α_{j}})

= ∑_{j=1}^{t} ∑_{i=1}^{L} h_{i}α_{j}^{ϵ_{i}}f_{j}(x − b)^{α_{j}}

= ∑_{j=1}^{t} h(α_{j})f_{j}(x − b)^{α_{j}}
= 0,  since h ∈ I.

□

For polynomials f_{1}, f_{2}, …, f_{m}, f let us define a generalized Wronskian matrix and the ω-vector as

W_{m}(f_{1}, f_{2}, …, f_{m}) = ( \begin{matrix}f_{1} & f_{2} & … & f_{m} \\ D(f_{1}) & D(f_{2}) & … & D(f_{m}) \\ ⋮ & ⋮ & ... & ⋮ \\ D^{m−1}(f_{1}) & D^{m−1}(f_{2}) & … & D^{m−1}(f_{m}) \end{matrix} ) , ω_{m}(f) = ( f ; D(f) ; ⋮ ; D^{m−1}(f) ) ,

where D is a generic linear combination of ∂/∂x_{i}, i.e., D = ∑_{i=1}^{n} l_{i}∂/∂x_{i}, l_{i} ∈ F. As usual, D^{i} denotes the operator D applied i times. Let

W_{m}(f_{1}, f_{2}, …, f_{m}) = det(W_{m}(f_{1}, f_{2}, …, f_{m})).

Clearly,

W_{L}(f_{ϵ_{1}}, …, f_{ϵ_{L}}) ∈ F[z_{1}, …, z_{n}, x_{1}, …, x_{n}]

and

W_{L}(f̂_{ϵ_{1}}, …, f̂_{ϵ_{L}}) ∈ K[x_{1}, …, x_{n}].

Lemma 5 If a set of polynomials F = {f̂_{ϵ_{1}}, f̂_{ϵ_{2}}, …, f̂_{ϵ_{L}}} satisfies the K-linear relation

f̂_{ϵ_{1}} + h_{2}f̂_{ϵ_{2}} + h_{3}f̂_{ϵ_{3}} + … + h_{L}f̂_{ϵ_{L}} = 0,

and no proper subset of F satifies a K-linear relation, then h⃗ = (h_{2} h_{3} … h_{L})^{T r} is the unique solution to the system of equations

W_{L−1}(f̂_{ϵ_{2}}, …, f̂_{ϵ_{L}})ϕ⃗ = −ω_{L−1}(f̂_{ϵ_{1}}).

Proof: The proof is classical and we give only a brief sketch. By rewriting the linear dependency as

−f̂_{ϵ_{1}} = h_{2}f̂_{ϵ_{2}} + h_{3}f̂_{ϵ_{3}} + … + h_{L}f̂_{ϵ_{L}}

and applying the operator D successively L − 1 times, it follows that

W_{L−1}(f̂_{ϵ_{2}}, …, f̂_{ϵ_{L}})h⃗ = −ω_{L−1}(f̂_{ϵ_{1}}).

If h is not the only solution, then W_{L−1}(f̂_{ϵ_{2}}, …, f̂_{ϵ_{L}}) is singular, i.e., the Wronskian W_{L−1}(f̂_{ϵ_{2}}, …, f̂_{ϵ_{L}}) vanishes identically. This implies that the polynomials f̂_{ϵ_{2}}, …, f̂_{ϵ_{L}} are K-linearly dependent (Kaplanski 1957). But no proper subset of F is supposed to satify a K-linear relation, hence, h is the only solution. □

From lemmas 4,5, it follows that if b⃗ = (b_{1}, b_{2}, …, b_{n}) is a t-sparse shift for f, and if h = y^{ϵ_{1}}+h_{2}y^{ϵ_{2}}+h_{3}y^{ϵ_{3}}+…+h_{L}y^{ϵ_{L}} ∈ I, then W_{L}(f_{ϵ_{1}}, …, f_{ϵ_{L}}) vanishes identically under the substitution z_{1} = b_{1}, z_{2} = b_{2}, …, z_{n} = b_{n}. A partial converse of lemma 4 is also true and it gives us a way to determine sparse shifts of f.

For a polynomial h = y^{ϵ_{1}}+h_{2}y^{ϵ_{2}}+h_{3}y^{ϵ_{3}}+…+h_{L}y^{ϵ_{L}}, let support(h) denote the list of polynomials f_{ϵ_{1}}, f_{ϵ_{2}}, …, f_{ϵ_{L}} and W(h) denote W_{L}(support(h)) where L is the cardinality of support(h).

Lemma 6 Let G = {g_{1}, g_{2}, …, g_{r}} ⊂ Q[y_{1}, y_{2}, …, y_{i}] be a reduced Gröbner basis for a zerodimensional ideal J with dim(Q[y_{1}, y_{2}, …, y_{i}]/J) ≤ t, and b⃗_{i} = (b_{1}, b_{2}, …, b_{i}) ∈ K^{i} such that

W(g_{1})_{z_{j}=b_{j}, j=1,...,i} = W(g_{2})_{z_{j}=b_{j}, j=1,...,i} = … = W(g_{r})_{z_{j}=b_{j}, j=1,...,i} = 0.

Then b⃗_{i} is a partial t-sparse shift for f, i.e., f = ∑_{j=1}^{t} ϕ_{j}(x − b)^{β_{j}}) where (x − b)^{β_{j}}) denotes (x_{1} − b_{1})^{β_{j,1}}(x_{2} − b_{2})^{β_{j,2}} … (x_{i} − b_{i})^{β_{j,i}} and J is the ideal of the points {β_{1}, β_{2}, …, β_{t}}.

Proof: Let f = ∑_{j=1}^{k} ϕ_{j}(x − b)^{β_{j}} with ϕ_{j} ≠ 0. We will show that k ≤ t. Suppose g_{1} = g_{1,1}y^{ϵ_{1}} + g_{1,2}y^{ϵ_{2}} +g_{1,3}y^{ϵ_{3}} +…+g_{1,L}y^{ϵ_{L}}. Since W(g_{1})_{z_{j}=b_{j}, j=1,...,i} = 0, by lemma 5, g_{1,1}f̂_{ϵ_{1}} +g_{1,2}f̂_{ϵ_{2}} +g_{1,3}f̂_{ϵ_{3}} + … + g_{1,L}f̂_{ϵ_{L}} = 0, i.e.,

= ∑_{p=1}^{L} g_{1,p}(∑_{j=1}^{k} β_{j}^{ϵ_{p}}ϕ_{j}(x − b)^{β_{j}})

= ∑_{j=1}^{k} ∑_{p=1}^{L} g_{1,p}β_{j}^{ϵ_{p}}ϕ_{j}(x − b)^{β_{j}}

= ∑_{j=1}^{k} g_{1}(β_{j})ϕ_{j}(x − b)^{β_{j}}
= 0.

Therefore, g_{1}(β_{j}) = 0 for j = 1, …, k. In a similar fashion, we can show that g_{p}(β_{j}) = 0 for each g_{p} ∈ G. In other words, each β_{j} is a zero of the ideal J. Since dim(Q[y_{1}, y_{2}, …, y_{i}]/J) ≤ t, J can have at most t distinct zeros. Therefore, k ≤ t. □

2.1 The Sparse Shift Algorithm – Case of Finitely Many Shifts

Our goal is to construct the Gröbner basis G. Given any term s (in y_{1}, …, y_{n}), it is either

• reduced with respect to G, or

• the head term of some polynomial in G, or

• a multiple of some head term in G.

A term s is called a simple multiple of the term s′ if s = y_{j}s′ for some y_{j}. Every head term in G is a simple multiple of some reduced term. The idea is to enumerate f_{ϵ_{i}} in the increasing order according to ≺ and decide, with the help of a Wronskian test, which of the above properties the term y^{ϵ_{i}} satisfies. The idea of systematic enumeration is borrowed from the well-known FGLM basis conversion algorithm (Faugère et al 1993).

• The function Complete basis constructs a reduced Gröbner basis for the ideal I under the pure lexicographic term order with y_{1} ≺ y_{2} ≺ … ≺ y_{n}. In the process, it also computes the corresponding t-sparse shifts for f. If f has several t-sparse shifts satisfying the assumptions stated earlier, Complete basis constructs reduced Gröbner bases for the ideal I corresponding to each of the sparse shifts. It returns a set of ordered pairs (b⃗, G) where b⃗ is a t-sparse shift for f and G is a reduced Gröbner basis for the corresponding ideal I. It uses the functions Next term and Wronskian test.

• The function Next term takes four parameters:

• current basis, reduced terms and last shift index are passed in unmodified by Complete basis, and,

• new var is a flag that is set by Next term.

Next term returns the smallest term y^{ϵ}( =: s) (according to the ordering ≺) that is neither in reduced terms nor is a multiple of some known head term in current basis; it returns null if no such y^{ϵ} exists. When such a term exists, if it has a new variable, i.e., the number of variables in the term is > last shift index, the flag new var is set to true, else new var is set to false.

• The function Wronskian test takes five parameters: s,reduced terms, list of bis, last shift index and new var which are all passed in unmodified by Complete basis and it operates in two distinct modes.

• If new var = true, then Wronskian test tries to extend the partial t-sparse shift list of bis = (b_{1}, …, b_{k}) to the next variable, i.e., variable whose index is last shift index + 1 = k + 1. If it finds a possible shift, it returns the pair of values [ true, bis or poly ] where bis or poly contains all the zeros of the content of W_{L+1}(f_{ϵ}, f_{ϵ_{1}}, …, f_{ϵ_{L}})_{z_{i}=b_{i}, i=1,...,k} as an element of K[z_{k+1}][x_{1}, x_{2}, …, x_{n}] where y^{ϵ_{1}}, …, y^{ϵ_{L}} ∈ reduced terms, and L = Cardinality(reduced terms). If the content is 1, then Wronskian test returns the pair of values [ false, [ ] ] .

• If new var = false, then Wronskian test attempts to solve the system of linear equations

W_{L}(f̂_{ϵ_{1}}, …, f̂_{ϵ_{L}})ϕ⃗ = −ω_{L}(f̂_{ϵ})

where y^{ϵ_{1}}, …, y^{ϵ_{L}} ∈ reduced terms, and L = Cardinality(reduced terms). If there is a solution ϕ = (g_{1} g_{2} … g_{L})^{T r}, then Wronskian test returns the pair of values [ true, y^{ϵ} + g_{1}y^{ϵ_{1}} + … + g_{L}y^{ϵ_{L}} ] where y^{ϵ} = s. If there is no solution, then Wronskian test returns the pair of values [ false, [ ] ] .

Complete basis( f, current basis, /* set containing a partially constructed Gröbner basis; */ reduced terms, /* set containing terms known to be reduced with respect to current basis; */

list of bis,  /* list containing partial shifts; if list of bis = (b_{1}, …, b_{k}), then the b_{i} are possible t-sparse shifts for x_{i}, i = 1, …, k; */

last shift index /* index of the last variable for which a shift has been com- puted; last shift index = Cardinality(list of bis) always; */

term limit  /* bound on the number of terms (t) in the shifted sparse representation of f; */

)

begin local

set of GBs, /* (b⃗, G) are accumulated in this; */ new var,Wr flag,

bis or poly;  /* flags and place-holders for return values and return status of Next term, Wronskian test; */

set of GBs := { }; if Cardinality(reduced terms) ≤ term limit then

if s := Next term( reduced terms, current basis,

last shift index, new var) then [ Wr flag, bis or poly ] := Wronskian test( s, reduced terms, list of bis, last shift index, new var);

if new var then if Wr flag then for each b ∈ bis or poly do set of GBs := set of GBs ∪ Complete basis(reduced terms,

current basis, [list of bis, b], last shift index+1, term limit);

od ; set of GBs := set of GBs ∪ Complete basis(reduced terms ∪ s,

current basis, list of bis, last shift index, term limit);

return set of GBs  /* new var is true, and the flag Wr flag is set to trueby Wronskian test ; therefore, one or more possible shifts were found for the last shift index+1-th variable; branch out to complete each of the shifts; continue the original branch also in search of other shifts. */

else set of GBs := Complete basis(reduced terms ∪ s,

current basis, list of bis, last shift index, term limit)

return set of GBs; fi /* No possible shift was found for the last shift index+1- th variable; no Q-linear combination of the term s and all the lower terms known to be reduced with respect to cur- rent basis belongs to I; classify s as reduced with respect to current basisand continue.*/

else if Wr flag then set of GBs := Complete basis(reduced terms, current basis ∪

bis or poly, list of bis, last shift index, term limit)

return set of GBs  /* new var is false, and the flag Wr flagis set to trueby Wronskian test , therefore, a Q-linear combination of the term s and all the lower terms known to be reduced with respect to current basis was found to belong to I; update the current basis and continue. */

else set of GBs := Complete basis(reduced terms ∪ s,

current basis, list of bis, last shift index, term limit)

return set of GBs; fi ; /* no Q-linear combination of the term s and all the lower terms known to be reduced with respect to current basis belongs to I; classify s as reduced with respect to cur- rent basisand continue.*/

fi else

return { (list of bis, current basis) } /* Next term failed to return a new term; so, every term is either known to be reduced with respect to current basis, or is a multiple of some head term in current basis. This means that current basis is a zero-dimensional Gröbner basis; return the basis and the corresponding shift. */

fi else

return { } ;  /* The basis being built has more than t reduced terms which means that the shift being computed can not be completed to a t-sparse shift. */

end.

As described earlier, the function Wronskian test returns a pair of values consisting of a flag (to denote what was computed) and either a list of possible shifts or a new element of the Gröbner basis being constructed. It uses two functions, Roots of Content and Lin Sys Solve.

• The function Roots of Content takes 2 parameters, an index k, 0 ≤ k ≤ n − 1, and a polynomial w ∈ K[z_{k+1}][x_{1}, x_{2}, …, x_{n}]. Roots of Content returns a list of all the zeros of the content of the polynomial w.

• The function Lin Sys Solve takes 2 parameters, W (∈ K^{m×m}), ω(∈ K^{m}), and attempts to solve the m × m system of linear equations W ϕ⃗ = ω. If the system has a solution, it returns a g⃗ ∈ K^{m} such that W g⃗ = ω, else it returns a null list. In our case, the system of equations either has a unique solution or no solution.

Wronskian test( s, reduced terms, list of bis, last shift index, new var )

begin local

shift candidates, /* used to accumulate the list of zeroes of the content of W_{L+1}(f_{ϵ}, f_{ϵ_{1}}, …, f_{ϵ_{L}})_{z_{i}=b_{i}, i=1,...,k} as an element of K[z_{k+1}][x_{1}, x_{2}, …, x_{n}] where y^{ϵ_{1}}, …, y^{ϵ_{L}} ∈ reduced terms, and L = Cardinality(reduced terms). */

list of coeffs, /* used to store the list of coefficients returned by Lin Sys Solve. */ basis element, /* used to store a polynomial y^{ϵ} + g_{1}y^{ϵ_{1}} + … + g_{L}y^{ϵ_{L}} that will become part of the Gröbner basis being constructed. */

k;

k := last shift index;

if new var then

shift candidates := Roots of Content(last shift index + 1,

W_{L+1}(f_{ϵ}, f_{ϵ_{1}}, …, f_{ϵ_{L}})_{z_{i}=b_{i}, i=1,...,k});

if shift candidates ≠ [ ] then return [ true, shift candidates ] else

return [ false, [ ] ] fi else

list of coeffs := Lin Sys Solve(W_{L}(f̂_{ϵ_{1}}, …, f̂_{ϵ_{L}}), −ω_{L}(f̂_{ϵ}));

if list of coeffs ≠ [ ] then /* list of coeffs is a list of coefficients [g_{1}, g_{2}, …, g_{l}] */ basis element := y^{ϵ} + g_{1}y^{ϵ_{1}} + g_{2}y^{ϵ_{2}} + … + g_{L}y^{ϵ_{L}}; return [ true, basis element ]

else

return [ false, [ ] ] fi fi

end. Initially, Complete basis is invoked with the following parameter values: the polynomial f, current basis set to { } , reduced terms set to {1}, list of bis set to ( ), last shift index set to 0, and term limit set to t.

Correctness of the Algorithm: The algorithm is in one of two states always, new varbeing true and new varbeing false.

State 1, new var true: In this state, the algorithm is attempting to extend current basis to include a polynomial in F[y_{1}, …, y_{i+1}] (where i = last shift index). Since the Gröbner basis being constructed is a lexicographically ordered basis, and the partial basis current basis is being built from the smallest head term up, current basis is actually a reduced Gröbner basis for a zero-dimensional ideal J in F[y_{1}, …, y_{i}] with dim(Q[y_{1}, y_{2}, …, y_{i}]/J) ≤ t, and the elements of current basis satisfy the Wronskian tests with respect to the candidate shift list of bis. By theorem 6, list of bis is indeed a partial shift for f and a candidate for a complete shift. There are at most finitely many ways in which this partial shift can be extended to variable x_{i+1} (that is our assumption). Each possible way to extend the shift appears as a root of the content of an appropriate Wronskian, the algorithm tries to find the contents of all such Wronskians by exhaustive enumeration. If the content of a Wronskian has more than one root, each root is a possible way to continue the shift list of bis and the algorithm branches into as many branches as the roots and also continues the computations along the parent branch. When the algorithm starts a new branch, it enters state 2.

State 2, new var false: In this state, the algorithm already has a candidate shift for x_{1}, …, x_{i} and is attempting to extend current basis to include relevant elements of the Gröbner basis that are in F[y_{1}, …, y_{i}] (where i = last shift index). If a polynomial h ∈ F[y_{1}, …, y_{i}] belongs to the ideal under construction, then by lemmas 4 and 5, h satisfies the Wronskian test. If h is indeed an element of the target Gröbner basis, then it will be generated by the algorithm because of the particular order in which the Wronskian tests are performed (this is one of the key ideas in the FGLM algorithm; see (Faugere et al) for details). The algorithm enters state 1 when after it generates a basis element whose head term is a pure power of x_{i} (because Next term now sets new var to true).

In every branch, in either state, after a Wronskian test, either a new reduced term is deduced or a new basis polynomial is generated. The number of basis elements is bounded by term limit and the number of reduced terms is bounded by term limit. If the cardinality of reduced terms exceeds term limit, then the ideal J has more than t zeros, i.e., the partial shift list of bis cannot be extended to x_{i+1} and the algorithm terminates that branch. Therefore, each branch terminates after at most (n + 1)t Wronskian tests. When a branch terminates, if the number of reduced terms is ≤ term limit, then the branch must have terminated because Next term returned the empty list; this means that every term is known to be either reduced or a head term or a multiple of a head term with respect to current basis which means that current basis is a reduced Gröbner basis in Q[x_{1}, …, x_{n}].

If there are infinitely many t-sparse shifts with respect to one or more x_{i} ( t > deg_{x_{i}}(f)) but only finitely many of them can be combined to t-sparse shifts for the polynomial f, then the algorithm can encounter an identically vanishing Wronskian in state 1 and will fail. An example of a polynomial for which such a phenomenon happens is (x − 1)^{2}(y − 1)^{2}(z − 1)^{2} + (x − 1)(y − 1)(z − 1) + 2. (1, 1, 1) is a 3-sparse shift for this polynomial and it is the only 3-sparse shift for the polynomial. However, with respect to any single variable, there are infinitely many 3-sparse shifts. This situation is taken care of by the algorithm of the next section as the “low degree case”.

The function Complete basis returns a set of ordered pairs (b⃗, G) where b⃗ is a t-sparse shift for f and G is a reduced Gröbner basis for the corresponding ideal I. From these, it is quite straightforward to compute the shifted sparse representation f = ∑_{i=1}^{t} f_{i}u^{α_{i}} corresponding to the pair (b⃗, G).

• Find the zeros of the ideal (G). We know that the zeros are the multi-indices α_{i} in the above representation. Since we already have a reduced, lexicographically ordered Gröbner basis and know that all the roots are n-tuples of integers, this can be done fast. See (Lazard 1992).

• Once we know the shift b⃗ and the multi-indices α_{i}, we can find the coefficients of f from its values at t selected points by solving a t × t system of linear equations.

Complexity Analysis: In the following analysis, our main goal is to get an upper bound on the number of Q-operations performed by the sparse shift algorithm. Our emphasis is not so much on getting the sharpest possible bounds (as that would depend on the intricate details of how each step is implemented) as on finding the coarse dependence (polynomial or exponential) of the running time (number of Q-operations) on n, t, d. We choose the primitive element methos for computing with algebraic numbers for convenience and the other models that one finds in the literature are polynomially related to this (polynomial in the degree of the extension under consideration).

The polynomial f is assumed to have rational number coefficients. In fact, the algorithm could run over any field of characteristic zero, but in the complexity analysis we assume that the coefficients of f are rationals. The t-sparse shifts may be algebraic over Q as we have seen. We assume that the algebraic numbers that arise in a particular branch of the algorithm are expressed as Q-linear combinations of 1, ζ, ζ^{2}, … where ζ is a primitive element of the smallest degree algebraic extension over Q that contains all the algebraic numbers that arise in that branch.

The main operation in the algorithm is the Wronskian test in Complete basis. We know that in each branch generated by Complete basis, there are at most n(t + 1) Wronskian tests. How many distinct branches can there be? Notice that branching can take place only when Complete basis is in state 1. Branching corresponds to finding more than one root for the content of a Wronskian. Consider a Wronskian test in state 1 with last shift index = i. The Wronskians are determinants of matrices of size no more than (t + 1) × (t + 1) and the degree of each entry in z_{i+1} is no more than t, therefore, the content of the Wronskian, which is a polynomial in z_{i+1}, has degree O(t^{2}) in z_{i+1}. Hence a worst case bound on the number of branches at a time is O(t^{2}). The main branch in state 1 can branch at most t times and the main branches are the only ones that can branch. The branching stops when last shift index becomes equal to n. Therefore, the number of branches is bounded by O(t^{3n}). Conclude that there are O(nt^{3n+1}) Wronskian tests in all.

Consider any branch generated by Complete basis at a time when last shift index = i. The entries of the Wronskian are polynomials in z_{i+1} and x_{1}, …, x_{n} with coefficients from an algebraic extension of the rationals of degree O(t^{2i}) (the contents of the previous Wronskians whose roots form the partial shift along the chosen branch are of degree O(t^{2}) in z_{j}, j < i + 1, and if each of the contents is irreducible over the earlier extensions, the current extension will have degree O(t^{2i}) over Q). Each arithmetic operation in such an extension costs O(t^{4i}) Q-operations.

To compute the content of a Wronskian, we have to compute the Wronskian (a (t + 1) × (t + 1) determinant at most). The Wronskian is a polynomial of degree O(dt) in the x_{i} and O(t^{2}) in z_{i+1}. Computing the Wronskian and then its content costs O(t^{3}(dt)^{2n}t^{4n}) = O(d^{2n}t^{6n+3}) Q-operations if we are to do it deterministically. If we are allowed to use randomization, we can substitute two different sets of random rational numbers for the x_{i} in the matrix corresponding to the Wronskian and compute the univariate gcd of the determinants of the matrix under the two specializations. With high probability, the gcd will be the content of the Wronskian. The cost of doing this is O(t^{3}.t^{4}.t^{4n}) = O(t^{4n+7}) Q-operations.

The last two steps in the sparse shift algorithm are much simpler. Finding the integer roots of a zero-dimensional lexicographic Gröbner basis with at most t zeros can be done by finding all the integer roots of a univariate polynomial of degree at most t and evaluating n − 1 other univariate polynomials of degree < t at the roots this polynomial and solving an n × n linear system (see (Lakshman 1990)). The total cost is bounded by O(t^{3} + nt^{2} + n^{3}) Q-operations. This has to be done at most O(nt^{3n+1}) times (once for each branch of Complete basis).

Setting up and solving a t × t linear system to compute the coefficients in the shifted sparse representation costs O(t^{3}) Q-operations and this too has to be done at most O(nt^{3n+1}) times, once for each branch of Complete basis. Adding up the costs of all the steps, we have:

Theorem 4 The algorithm of this section computes all shifted t-sparse representations for f provided deg_{x_{i}}(f) ≥ t for each x_{i}. If randomization is not allowed, the algorithm performs (n(dt)^{n})^{O(1)} Q-operations. If randomization is allowed, the algorithm performs (nt^{n})^{O(1)} Q-operations. □

For the special case of t ≤ (deg_{x_{i}}(f) + 1)/2 for each x_{i}, the polynomial f has at most one t-sparse shift by lemma 2, the algorithm runs much faster. In this case, there is essentially no branching and all individual shifts are rational. For this case, we have:

Theorem 5 If t ≤ (deg_{x_{i}}(f) + 1)/2 for each x_{i}, the algorithm of this section computes a shifted t-sparse representations for f (if it has one) in time polynomial in t. More specifically, if randomization is not allowed, the algorithm performs (ntd^{n})^{O(1)} Q-operations. If randomization is allowed, the algorithm performs (nt)^{O(1)} Q-operations. □

Remark: If we do not have the derivarives of f available, but only a black box that evaluates f, we can still use ideas close to the ones describes in this section to construct a Gröbner basis for the ideal of points p^{α_{1}}, p^{α_{2}}, …, p^{α_{t}} with p_{j}^{α_{j}} denoting the power product p_{1}^{α_{j,1}}p_{2}^{α_{j,2}} … p_{t}^{α_{j,t}} where the p_{i} are distinct prime integers. In place of f_{ϵ}, one can use f(p_{1}^{ϵ_{1}} + z_{1}, p_{2}^{ϵ_{2}} + z_{2}, …, p_{n}^{ϵ_{n}} + z_{n}).

3. Sparse Shift Algorithm – The General Case

In this section, we deal with the case of finding t-sparse shifts for a polynomial f ∈ F[x_{1}, …, x_{n}] (we suppose that F ⊂ R) for which deg_{x_{i}}(f) < 2t for one or more x_{i} and consequently, there might be several (possibly infinitely many) t-sparse shifts with respect to x_{i} alone. This covers the case that causes the algorithm in the previous section to fail. We assume that the polynomial f is given to us as a black box or a straight line program. Assuming that we are given a bound M on the sizes of the coefficients of f in the standard representation (and without knowing anything about the degree of f), we identify three kinds of variables in f:

• variables that have very high degree (Ω((M + t^{2})^{1+ϵ})), these appear with a unique t-sparse shift that is either 0, or 1 or -1. We make repeated use of the Ben-Or and Tiwari algorithm (Ben-Or, Tiwari 1988) in conjunction with some bounds proved in this section to identify the high degree variables and to find the corresponding shifts.

• variables that have moderate degree (2t < degree < c(M +t^{2})^{1+ϵ} for some constant c), these appear with a unique, rational shift in any t-sparse representation of f. We make repeated use of the sparse-shift algorithm in (Lakshman-Saunders 1994) to identify the moderate degree variables.

• variables that have low degree (degree < 2t), these may appear with different, algebraic shifts in different t-sparse representation of f. We make use of the algorithm in (Grigoriev-Karpinski 1993) to identify the low degree variables and find their sparse shifts.

The algorithm used to find sparse shifts for the low degree variables (Grigoriev-Karpinski 1993) is very general and in fact solves the problem of computing sparse shifts completely without making such distinctions as listed here. However, we apply it selectively, to polynomials of low degree (d ≤ 2tn), and as a result, the complexity of the algorithm comes down from O(d^{O(n^{4})}) to O((nt)^{O(n^{2})}) operations.

3.1 Degree Bounds on Shifted Sparse Polynomials

We establish some bounds on the degrees of f in each x_{i} in terms of t, n, M. These bounds are used in the main algorithm which unfolds in the rest of the section. The main steps of the algorithm are collected together at the end of the section.

Let f = ∑_{i=0}^{d} f_{i}x^{i} = ∑_{j=1}^{t} f_{j}(x − b)^{δ_{j}}. For a rational number p/q ∈ F, define size(p/q) = ⌈log_{2}(p+1)⌉+⌈log_{2}(q+1)⌉ and for the polynomial f ∈ F[x] as above, define size(f) = max{size(f_{i})}. Let size(f) ≤ M for an integer M.

Lemma 7 For any ϵ > 0 there exists c such that if f has two sparse representations as above where b ∉ {0, 1, −1}, then d ≤ c(M + t^{2})^{1+ϵ}.

Proof: If d ≤ 3t, the lemma is obvious. If d > 3t, then b is unique and rational (Lakshman-Saunders 1994) and we have f = ∑_{i=0}^{d} f_{i}x^{i} = ∑_{j=1}^{t} f_{j}(x − b)^{δ_{j}} with d = δ_{t} > δ_{t−1} > … > δ_{1}. Rewrite this as

f⃗B_{t}CD_{d+1} = f⃗

where f⃗ = (f_{δ_{1}} f_{δ_{2}} … f_{δ_{t}}), B_{t} is the t × t diagonal matrix with B_{t}(j, j) = b^{δ_{j}}, C is a t × (d + 1) matrix with C(i, j) = \binom{δ_{i}}{d+1−j}, D_{d+1} the (d + 1) × (d + 1) diagonal matrix with B_{t}(j, j) = b^{−d+j−1}, and, f⃗ = (f_{d} f_{d−1} … f_{0}). Let u⃗ = (u_{1} u_{2} … u_{t}) = f⃗B_{t}. Let C_{t} be the t × t submatrix of C consisting of the last t columns of C.

We know that C_{t} is non-singular (see (Lakshman-Saunders 1994)) and hence,

u⃗ = (f_{t−1} f_{t−2} … f_{0}) ( \begin{matrix}b^{t−1} & 0 & … & 0 \\ 0 & b^{t−2} & … & 0 \\ ⋮ & ⋮ & ... & ⋮ \\ 0 & 0 & … & b^{0} \end{matrix} ) C_{t}^{−1}.

Therefore, size(u_{1}) < M + tsize(b) + ct^{2} log d for some constant c, the ct^{2} log d term coming from C_{t}^{−1}. Since f_{d} = u_{1}b^{−d}, we get

size(u_{1}) > size(b^{d}) − M > d(size(b) − 2) − M

Since b ∉ {0, 1, −1}, we have size(b) ≥ 3 and therefore, d−3t ≤ (d−3t)(size(b)−2) ≤ c′(M+t^{2} log d) for some constant c′. It follows that d/logd ≤ c′(M + t^{2}) and d = O((M + t^{2})^{1+ϵ}) for any ϵ > 0. □

Lemma 8 Let g ∈ F[x_{1}, x_{2}, …, x_{n}] be a non-zero polynomial that is shifted t-sparse and S = {a_{1} < a_{2} < … < a_{2t}} ⊂ R. Then g cannot vanish everywhere on S^{n}.

Proof: The lemma follows from successive applications of Descartes’ rule to each variable x_{1}, x_{2}, …, x_{n}. □

Choose a constant ĉ and an arbitrarily small ϵ > 0 such that deg_{x_{i}}(f) < ĉ(M + t^{2})^{1+ϵ} for 1 ≤ i ≤ n in lemma 7. Denote this degree bound ĉ(M + t^{2})^{1+ϵ} by D.

Let 1 ≤ i ≤ n. Fix b to be one of 0, −1, −1 and let X_{i} = x_{i} −b. For each W = (w_{1}, w_{2}, …, w_{n}) ∈ S^{n−1}, consider the univariate polynomial f_{W} (X_{i}) = f(w_{1}, …, w_{i−1}, x_{i}, w_{i+1}, …, w_{n}). By lemma 8, it follows that there is a W such that deg(f_{W} (X_{i})) = deg_{x_{i}}(f). Fix such a W.

Let A_{1} = {−2, (−2)^{2}, (−2)^{3}, …, (−2)^{2t}}, A_{2} = {−3, (−3)^{2}, (−3)^{3}, …, (−3)^{2t}}, A_{3} = {−5, (−5)^{2}, (−5)^{3}, …, (−5)^{2t}}, and, A_{4} = {−7, (−7)^{2}, (−7)^{3}, …, (−7)^{2t}}. Apply the Ben-Or and Tiwari interpolation algorithm (Ben-Or, Tiwari 1988) to f_{W} (X_{i}) at the four sets of points A_{1}, A_{2}, A_{3}, A_{4} respectively. If the algorithm succeeds, it returns a t-sparse polynomial in F[X_{i}] for each set of evaluation points A_{1}, A_{2}, A_{3}, A_{4}. Let us denote the polynomials returned by the Ben-Or and Tiwari algorithm by g_{1}, g_{2}, g_{3}, g_{4} ∈ F[X_{i}] respectively.

Lemma 9 deg_{x_{i}}(f) > D if and only if g_{1} = g_{2} = g_{3} = g_{4} and deg(g_{1}) > D.

Proof: (⟹) : Since deg_{x_{i}}(f) > D, then by lemma 7, we have b_{i} ∈ {−1, 0, 1} in any t-sparse shift (b_{1}, b_{2}, …, b_{n}) for f. Therefore, f_{W} (X_{i}) is t-sparse for any W ∈ S^{n−1} and since the interpolating polynomial produced by the Ben-Or and Tiwari algorithm (Ben-Or, Tiwari 1988) is unique, we have g_{1} = g_{2} = g_{3} = g_{4} = f_{W} (X_{i}). Since we choose a W such that deg(f_{W} (X_{i})) = deg_{x_{i}}(f), it follows that deg(g_{1}) > D.

(⟸) : Assume that for a certain W_{0}, the Ben-Or and Tiwari algorithm returns four t-sparse polynomials in X_{i} such that g_{1} = g_{2} = g_{3} = g_{4} and deg(g_{1}) > D. The polynomial f_{W_{0}}(X_{i}) is shifted t-sparse and coincides with the t-sparse polynomial g_{1} at 4t positive points and 4t negative points. Therefore, by theorem 5 in (Lakshman-Saunders 1994), we have f_{W_{0}}(X_{i}) = g_{1}(X_{i}). □

If the Ben-Or, Tiwari algorithm fails to return, for each b = 0, 1, −1, and for each W four t-sparse polynomials in X_{i} g_{1} = g_{2} = g_{3} = g_{4} with deg_{X_{i}}(g_{1}) > D, then f_{W} (x_{i}) has a non-zero t-sparse shift and deg_{x_{i}}(f) < D. Denote the set of all indices i such that deg_{x_{i}}(f) > D by I_{B}. The set I_{B} can be determined by performing ((log d)t^{n})^{O(1)} arithmetic operations (for each i, for each b, and for each W ∈ S^{n−1}, we have to perform the Ben-Or and Tiwari algorithm 4 times).

Let j ∈ {1, 2, …, n} \ I_{B}. For each W ∈ S^{n−1}, interpolate f_{W} (x_{j}) as a dense univariate polynomial. If for some W_{0}, deg(f_{W_{0}}(x_{j}) > 2t, then we can use the sparse-shift algorithm from (Lakshman-Saunders 1994) to find the unique b_{j} ∈ F such that f_{W_{0}}(X_{j}) is t-sparse relative to x_{j} − b_{j}.

Lemma 10 Let j ∈ {1, 2, …, n}\I_{B}. Then deg_{x_{j}}(f) > 2t iff for a certain W_{0} ∈ S^{n−1}, deg(f_{W_{0}}(x_{j})) > 2t. If the latter is true, then there is unique b_{j} ∈ F such that f_{W_{0}}(x_{j}) is t-sparse relative to x_{j} − b_{j}. Moreover,in any t-sparse shift (b′_{1}, b′_{2}, …, b′_{n}) of f, we have b′_{j} = b_{j}.

Proof: The existence of W_{0} follows from lemma 8. Since f_{W_{0}}(x_{j}) is shifted t-sparse, and its degree > 2t, the shift is unique (theorem 1, Lakshman-Saunders, 1994). If (b′_{1}, b′_{2}, …, b′_{n}) is a t-sparse shift for f, then f_{W_{0}}(x_{j}) is t-sparse with respect to x_{j} − b′_{j}. From the uniqueness of the shift, it follows that b′_{j} = b_{j}. □

Denote the set of all indices j ∈ {1, 2, …, n} \ I_{B} such that deg_{x_{j}}(f) > 2t by I_{C}. The set I_{C} can be determined by performing (Mt^{n})^{O(1)} arithmetic operations (for each j and for each W ∈ S^{n−1}, we have to perform a dense univariate interpolation and the sparse-shift algorithm of Lakshman-Saunders once; since the degree with respect to x_{j} is bounded by D, we have the above bound.)

From lemmas 9, 10, it follows that for every i ∈ I_{B} only the powers of x_{i} and for every j ∈ I_{C}, only the powers of x_{j} − b_{j} appear in any t-sparse representation of f. Moreover, any two shifted t-sparse representations of f contain the same power products of x_{i} and x_{j} − b_{j}. The last statement follows from lemma 8 by considering suitable f_{W} ’s. Consider the representation of f as

f = ∑_{1≤i≤t} f_{i}X_{I_{B}}^{α_{i}}(x − b)_{I_{C}}^{β_{i}}

where f_{i} ∈ F[{x_{k}}], k ∈ {1, 2, …, n}\(I_{B}∪I_{C}), X_{I_{B}}^{α_{i}} denotes a power product of X_{i}, i ∈ I_{B}, and (x− b)_{I_{C}}^{β_{i}}. Note that the degrees of the f_{i} in any x_{k} are less than 2t. Therefore, the dense representations of f_{i} as F-linear combinations of power products of x_{k} for k ∈ {1, 2, …, n} \ (I_{B} ∪ I_{C}) have O(t^{n}) terms. We explicitly compute the f_{i} by dense interpolation from their values at O(t^{n}) points. The values of the f_{i} for a particular specialization S taking x_{k} to v_{k} ∈ F for k ∈ {1, 2, …, n} \ (I_{B} ∪ I_{C}) are obtained by constructing S(f) ∈ F[I_{B} ∪I_{C}] (S(f) denotes the image of f under the substitution S). For any S, S(f) is obtained by sparse interpolation using the Ben-Or and Tiwari algorithm (since we know the non-zero shifts b_{i}, this can be achieved by direct application of the Ben-Or and Tiwari algorithm). The cost of constructing the f_{i} this way is ((log d)t^{n})^{O(1)} F-operations.

Once we have the f_{i}, the problem is to find shifts b_{k} for the x_{k} for k ∈ {1, 2, …, n} \ (I_{B} ∪ I_{C}) such that the total number of terms in all the f_{i} represented in the power products of x_{k} − b_{k} is at most t. This is done by finding a shifted t-sparse representation for the polynomial Ψ(x_{k}, z) = f_{1}z^{2t} + f_{2}z^{2t−1} + … + f_{t}z^{t+1} where z is a new unknown. In any t-sparse shift of the above polynomial, the shift with respect to z has to be 0 since its degree in z is greater than 2t − 1. To find a t-sparse shift for Ψ, we apply the algorithm from (Grigoriev-Karpinski 1993) which finds the variety V of all t-sparse shifts of f, i.e., set of r-tuples (b_{k_{1}}, b_{k_{2}}, …, b_{k_{r}}) (r = n − cardinality(I_{B} ∪ I_{C}) and {k_{1}, …, k_{r}} = {1, …, n} \ (I_{B} ∪ I_{C})) such that the total number of terms in all the f_{i} represented in the power products of x_{k} − b_{k} is at most t. The algorithm from (Grigoriev-Karpinski 1993) returns the variety V as a union of its irreducible components ∪_{l}V_{l} and for each V_{l} and i, the algorithm returns a set of exponent vectors e_{1}, …, e_{s} and a set of rational functions g_{1}(b_{k}|k ∈ {1, …, n} \ (I_{B} ∪ I_{C})), …, g_{s}(b_{k}|k ∈ {1, …, n} \ (I_{B} ∪ I_{C})) such that f_{i} = ∑_{l=1}^{s} g_{l}x^{e_{l}} for any (b_{k_{1}}, b_{k_{2}}, …, b_{k_{r}}) ∈ V_{l}. For this input, the number of operations performed by the algorithm of Grigoriev and Karpinski is bounded by O((nt)^{O(n^{2})}) since deg(f_{i}) ≤ 2tn. We now collect the main steps of the algorithm together and summarize its asymptotic time complexity in theorem 6:

Algorithm to find all multivariate sparse shifts

• Compute the index set I_{B} of variables of high degree in f. If the Ben-Or and Tiwari algorithm fails to produce t-sparse g_{1}, g_{2}, g_{3}, g_{4} satisfying lemma 9 for a variable X_{i} = x_{i} − b, for every b ∈ {0, 1, −1}, then that variable has degree ≤ D.

• For each variable x_{j} whose index j is not in I_{B}, construct f_{W} (x_{j}) for all possible W by dense interpolation assuming that deg(f_{W} (x_{j})) ≤ D. Denote the interpolant by g_{W} (x_{j}) and perform the sparse shift algorithm of Lakshman-Saunders on g_{W} (x_{j}) whenever deg(g_{W} (x_{j})) > 2t. If each time, we discover the same shift, note that j ∈ I_{C}.

• Consider the representation of f as

f = ∑_{1≤i≤t} f_{i}X_{I_{B}}^{α_{i}}(x − b)_{I_{C}}^{β_{i}}.  (1)

Determine the f_{i} by dense interpolation as outlined in the text and compute sparse shifts for Ψ(x_{k}, z) = f_{1}z^{2t}+f_{2}z^{2t−1}+…+f_{t}z^{t+1} by the algorithm of Grigoriev and Karpinski. For each such shift of the f_{i}, return the corresponding t-sparse representation f obtained by substituting the shifted sparse representations of the f_{i} into the representation (1) above.

Theorem 6 The algorithm of this section computes all shifted t-sparse representations for f. The algorithm performs O(M^{O(1)}(nt)^{O(n^{2})}) operations. □

4. Discussion

In this paper, we have investigated the problem of finding t-sparse shifts for multivariate polynomials. The first algorithm has the advantage that the unknown shifts are obtained as zeros of univariate polynomials over algebraic extensions of Q (as opposed to having to solve general systems of polynomial equations). The algorithm uses Gröbner basis techniques. However, this algorithm cannot handle all polynomials with t-sparse shifts. The second algorithm is a complete algorithm slower than the first, but improves significantly over previously known algorithms. Several interesting issues remain unresolved at this time:

• find a necessary and sufficient condition for the uniqueness of t-sparse shifts and sparsifying linear transformations for multivariate polynomials.

• find an algorithm that handles the low degree case more efficiently than our algorithm. It is intriguing to see what might happen if we try to construct Gröbner bases as in the first algorithm, but with respect to other term orderings to handle the low degree case.

• find efficient algorithms for finding sparsifying linear tranformations (see Grigoriev and Karpinski, 1993). The Gröbner basis techniques can be extended to find sparsifying linear transformations for bivariate polynomials efficiently (polynomial in t, d).

It is also interesting to consider more general transformations such as ones leading to sparse decompositions of polynomials.

References

Baur, W., and Strassen, V., (1983), ”The complexity of partial derivatives,” Theoretical Computer Science, Vol. 22, pp. 317–330. Becker, T., and Weispfenning, V., (1993), “Gröbner bases – a computational approach to commutative algebra,” Graduate texts in Mathematics, Vol. 141, Springer-Verlag, New York. Ben-Or, M. and Tiwari, P. (1988), “A deterministic algorithm for sparse multivariate polynomial interpolation,” Proc. 20th Symp. Theory of Computing, ACM Press, pp. 301–309. Borodin, A. and Tiwari, P. (1990), “On the decidability of sparse univariate polynomial interpolation,” Proc. 22nd Symp. Theory of Computing, ACM Press, pp. 535–545. Clausen, M., Dress, A., Grabmeier, J., Karpinski, M. (1988), “On zero testing and interpolation of k-sparse multivariate polynomials over finite fields,” TR 88.06.006, IBM Germany, Heidelberg Scientific Center, June 1988. Faugère J.C., Gianni P., Lazard D., Mora T., (1994) “ Efficient Computation of Zero-dimensional Gröbner Bases by Change of Ordering ”, Jour. Symb. Computation, 1994.

Grigoriev, D. Yu. and Karpinski, M. (1987), “The matching problem for bipartite graphs with polynomially bounded permanents is in NC,” Proc. 28th IEEE Symp. Foundations Comp. Sci., pp. 166–172. Grigoriev, D., Karpinski, M., and Singer, M. (1990), “Fast parallel algorithms for sparse multivariate polynomial interpolation over finite fields,” SIAM J.Comp, Vol. 19, pp. 1059–1063. Grigoriev, D., Karpinski, M., and Singer, M. (1991), “The interpolation problem for k-sparse sums of eigenfunctions of operators,” Advances in Appl Math, Vol 12, pp. 76-81. Grigoriev, D., Karpinski, M., and Odlyzko, A.M.(1992), Existence of short proofs of non-divisibility of sparse polynomials under the extended Riemann hypothesis,” Proc. ISSAC 92, ACM Press, pp. 117-122. Grigoriev, D., Karpinski, M., and Singer, M. (1993a), “Computational complexity of sparse real algebraic function interpolation,” Proc. MEGA ’92, Progress in Mathematics, Birkhauser, Vol. 109, pp. 91–104. Grigoriev, D., and Karpinski, M., (1993), “A zero-test and an interpolation algorithm for the shifted sparse polynomials,” Proc. AAECC-93, Lect. Notes in Comp. Sci., Vol. 673, pp. 162–169. Grigoriev, D., Karpinski, M., and Singer, M. (1994), “Computational complexity of sparse rational interpolation,” SIAM J. Comp., Vol. 23, pp. 1-11. Kaltofen, E. (1987), “Single-factor Hensel lifting and its application to the straight-line complexity of certain polynomials,” Proc. 19th Symp. Theory of Computing, ACM Press, pp. 443–452. Kaltofen, E. and Lakshman, Y.N. (1988), “Improved sparse multivariate polynomial interpolation algorithms,” Proc. ISSAC 1988, Rome, Italy, Springer-Verlag LNCS 358, pp. 467–474. Kaltofen, E., and Trager, B. (1990), “Computing with polynomials given by black boxes for their evaluations: Greatest common divisors, factorization, separation of numerators and denominators,” J. Symb. Comp. 9 (1990), pp. 301-320. Kaplanski, I. (1957), “An introduction to differential algebra,” Hermann, Paris. Kapur, D., and Lakshman Y.N., (1992), “ Elimination methods: An introduction,” in Symbolic and Numerical Computation for Artificial Intelligence, (Ed. Bruce Donald, Deepak Kapur, Joe Mundy), Academic Press, 1992, pp. 45–89. Lakshman Y.N., (1990) “ On the Complexity of Computing Gröbner Bases for Zero Dimensional Polynomial Ideals ”, Ph.D. Thesis, Rensselaer Polytechnic Institute, Troy, New York, December 1990. Lakshman, Y.N., and Saunders, B.D., (1993), “Sparse polynomial interpolation in non-standard bases,” to appear in SIAM J. Comp. Lakshman, Y.N., and Saunders, B.D., (1994), “On computing sparse shifts for univariate polynomials,” Proc. Int. Symp. Symb. Alg. Comp., 1994 (ISSAC 94), Oxford, UK, ACM Press. Lazard D., (1989), “ Solving Zero–Dimensional Algebraic Systems ”, Tech. Report.no.89–48, LITP, Universite Paris VI, June 1989. Mansour, Y. (1992), “Randomized interpolation and approximation of sparse polynomials ,” Proc. ICALP ’92. Zippel, R. (1990), “Interpolating polynomials from their values,” Jour. Symb. Comp., Vol. 9, No. 3, pp. 375–403.