Comparison of CAD-based Methods for Computation of Rational Function Limits

Adam Strzeboński

Wolfram Research Inc.
100 Trade Centre Dr.
Champaign, IL 61820
U.S.A.

adams@wolfram.com
http://members.wolfram.com/adams

Abstract

I will present five methods for computation of limits of real multivariate rational functions. The methods do not require any assumptions about the rational function and compute the lower limit and the upper limit. All methods are based on the cylindrical algebraic decomposition (CAD) algorithm, but use different formulations of the problem. I will show results of an empirical comparison of the methods on a large set of examples.

Problem statement

Let  g∈ℚ[(Subscript[x, 1]),…,(Subscript[x, n])],  h∈ℚ[(Subscript[x, 1]),…,(Subscript[x, n])]\{0},  D:={u∈(ℝ^n):h(u)≠0},  f:D∋u⟶((g(u))/(h(u)))∈ℝ, and let c∈(ℚ^n).

Problem 1: Find l∈(OverBar[ℝ])  such that  l=(Subscript[lim, x->c])f(x)  or prove that such l does not exist.

Problem 2: Find (Subscript[l, 1]),(Subscript[l, 2])∈(OverBar[ℝ])  such that  (Subscript[l, 1])=(Subscript[min lim, x->c])f(x)  and  (Subscript[l, 2])=(Subscript[max lim, x->c])f(x) .

Example

f=((4 x y)/((x^2)+2 (y^2)));

p1=Plot3D[f,{x,-1,1},{y,-1,1},PlotPoints->100,PlotStyle->{Opacity[0.7]}];
p2=Plot3D[{-Sqrt[2],Sqrt[2]},{x,-1,1},{y,-1,1},PlotStyle->{{Opacity[0.5],Red},{Opacity[0.5],Green}}];
Show[{p2,p1}]

(Underscript[\[Limit], {x,y}->0])f

(Underscript[\[MinLimit], {x,y}->0])f

(Underscript[\[MaxLimit], {x,y}->0])f

Recent related work

• P. Alvandi, M. Kazemi, and M. Moreno Maza. 2016. Computing Limits of Real Multivariate Rational Functions. In Proceedings of the International Symposium on Symbolic and Algebraic Computation, ISSAC 2016. ACM, 39–46.
• C. Cadavid, S. Molina, and J. D. Velez. 2013. Limits of quotients of bivariate real analytic functions. J. Symbolic Comp. 50 (2013), 197–207.
• J. D. Velez, J. P. Hernandez, and C. A. Cadavid. 2015. Limits of quotients of real polynomial functions of three variables. arXiv 1505.04121.

Solves Problem 1 under the additional assumption that c is an isolated zero of h . The authors use the theory of Lagrange multipliers to reduce the problem to computing the limit along a real algebraic set, and solve the reduced problem using regular chains methods.

• S. J. Xiao and G. X. Zeng. 2014. Determination of the limits for multivariate rational functions. Science China Mathematics 57, 2 (2014), 397–416.

Solves a modified version of Problem 1, namely decides whether the limit exists and is finite. The negative answer includes both the case when the limit does not exist and the case when the limit exists and is infinite. The algorithm uses Wu’s elimination method, rational univariate representations, and requires adjoining two infinitesimal elements to the field.

Simple solution using quantifier elimination

Eliminate quantifiers from the definition.

Problem 1

The limit of f at c is equal to l∈ℝ iff

(Subscript[∀, ϵ>0])(Subscript[∃, δ>0])(Subscript[∀, x∈(ℝ^n)])(0<∥x-c∥<δ∧h(x)≠0)⇒l-ϵ<f(x)<l+ϵ

The limit of f at c is equal to -∞ iff

(Subscript[∀, M>0])(Subscript[∃, δ>0])(Subscript[∀, x∈(ℝ^n)])(0<∥x-c∥<δ∧h(x)≠0)⇒f(x)<-M

The limit of f at c is equal to ∞ iff

(Subscript[∀, M>0])(Subscript[∃, δ>0])(Subscript[∀, x∈(ℝ^n)])(0<∥x-c∥<δ∧h(x)≠0)⇒f(x)>M

Problem 2

The lower limit of f at c is equal to  l∈ℝ iff

[(Subscript[∀, ϵ>0])(Subscript[∀, δ>0])(Subscript[∃, x∈(ℝ^n)])0<∥x-c∥<δ∧h(x)≠0 ∧f(x)<l+ϵ]∧
[(Subscript[∀, ϵ>0])(Subscript[∃, δ>0])(Subscript[∀, x∈(ℝ^n)])(0<∥x-c∥<δ∧h(x)≠0)⇒f(x)>l-ϵ]

The lower limit of f at c is equal to -∞.

(Subscript[∀, M>0])(Subscript[∀, δ>0])(Subscript[∃, x∈(ℝ^n)])0<∥x-c∥<δ∧h(x)≠0∧f(x)<-M

The lower limit of f at c is equal to ∞.

(Subscript[∀, M>0])(Subscript[∃, δ>0])(Subscript[∀, x∈(ℝ^n)])(0<∥x-c∥<δ∧h(x)≠0)⇒f(x)>M

How this can be improved

Quantifier elimination algorithms are highly sensitive to the number of variables -- the worst-case complexity of CAD algorithm is doubly exponential in the number of variables.  Hence it is very important in practice to find a formulation of the problem which introduces as few variables as possible. Direct application of the definitions results in CAD problems with n+3 variables. We propose five formulations that work with n+1 or n+2 variables and use simpler versions of the CAD algorithm customized for limit computation. We compare their efficiency on examples form literature and a large set of examples extracted from inputs submitted to the Wolfram Alpha website.

For simplicity I will talk about algorithms computing the lower limit. This is sufficient to solve both Problem 1 and Problem 2.  When both the lower limit and the upper limit are needed, certain steps of the computation can be shared. If only a solution of Problem 1 is required, it can often be obtained in a simpler way. The proceedings paper shows how the algorithms can be adjusted to efficiently obtain the required goal.

Without loss of generality let us assume that c=0, h(0)=0, and g and h are relatively prime.

Methods based on optimization

These methods use the definition of the lower limit in the following formulation

(Subscript[min lim, u->0])f(u)=(Subscript[lim, ϵ->(SuperPlus[0])])(inf_(∥u∥<ϵ∧u∈D))f(u)

where ∥u∥ := (max_(1≤i≤n))|(Subscript[u , i])|.

Algorithm OLIM1

This algorithm makes the most direct use of the definition.

Algorithm OLIM1:
Input:  g∈Q[x],  h∈Q[x]\{0}.
Output: (Subscript[min lim, u->0])((g(u))/(h(u))).
(1) Set  S(ϵ,y,x) := ϵ>0∧∥x∥<ϵ∧h≠0∧g=y h
[y is constrained to be equal to f(x)]
(2) Compute the CAD projection for S , with variables ordered ϵ<y<(Subscript[x, 1])<...<(Subscript[x, n]). 
Let (Subscript[F, v]) be the set of projection polynomials with the “largest” variable v.
(3) Choose a positive rational number e smaller than all positive roots of elements of (Subscript[F, ϵ]) .
(4) Construct cells extending the open cell with sample point (e), starting with the smallest values of y, until the first cell C on which S is true is found. Let r(ϵ) be the lower bound on y in C.
(5) If r(ϵ)=−∞, return −∞.
(6) Otherwise r(ϵ) is an algebraic function of ϵ. Return (Subscript[lim, ϵ->(SuperPlus[0])])r(ϵ).

The algorithm introduces two variables ϵ and y in addition to the variables x , however it constructs only one cell in ϵ and treats the variables x as existentially quantified. The formula S for which CAD is constructed always contains an equational constraint and hence the variant of CAD making use of equational constraints is used. The variable y is constrained to be equal to f(x) with ∥x∥<ϵ,  hence (Subscript[min lim, u->0])f(u)=(Subscript[lim, ϵ->(SuperPlus[0])])r(ϵ).

Limits of univariate real algebraic functions can be computed using a specialized version of the CAD cell adjacency computation algorithm described in:
A. Strzeboński.  CAD Adjacency Computation Using Validated Numerics. In Proceedings of the International Symposium on Symbolic and Algebraic Computation, ISSAC 2017. ACM, 413–420.

Algorithm OLIM2

This algorithm re-formulates the problem so that it only contains strict inequalities. This allows to use the algorithm GCAD -- a specialized version of CAD which uses a simpler projection operator and constructs only open cells.

Algorithm OLIM2:
Input:  g∈Q[x],  h∈Q[x]\{0}.
Output: (Subscript[min lim, u->0])((g(u))/(h(u))).
(1) Set  S(ϵ,y,x) := ϵ>0∧ ∥x∥<ϵ∧((h>0∧g<y h)∨(h<0∧g>y h))
[y is constrained to be greater than f(x)]
(2) Compute the GCAD projection for S , with variables ordered ϵ<y<(Subscript[x, 1])<...<(Subscript[x, n]). 
Let (Subscript[F, v]) be the set of projection polynomials with the “largest” variable v.
(3) Choose a positive rational number e smaller than all positive roots of elements of (Subscript[F, ϵ]) .
(4) Find rational numbers (Subscript[a, 0]),...,(Subscript[a, m]) such that (Subscript[a, 0])<(Subscript[r, 1])<(Subscript[a, 1])<...<(Subscript[r, m])<(Subscript[a, m]), where (Subscript[r, 1]),...,(Subscript[r, m]) are all the real roots of p(e,y) for p∈(Subscript[F, y]).
(5) Perform the binary search to find (Subscript[i, 0]): the minimum i for which (Subscript[∃, x])S(e,(Subscript[a, i]),x) is true.
[Note that if i<j then (Subscript[∃, x])S(e,(Subscript[a, i]),x)⇒(Subscript[∃, x])S(e,(Subscript[a, j]),x)]
(5) If (Subscript[i, 0])=0, return −∞.
(6) Otherwise let r(ϵ):=(Subscript[Root, y,λ])p(ϵ) be the algebraic function such that  p∈(Subscript[F, y]) and r(e)=(Subscript[r, (Subscript[i, 0])]).  Return (Subscript[lim, ϵ->(SuperPlus[0])])r(ϵ).

The algorithm introduces two variables ϵ and y in addition to the variables x , however it constructs only one cell in ϵ, performs a binary search on cells in y, and treats the variables x as existentially quantified. The variable y is constrained to be greater than f(u) with ∥u∥<ϵ,  hence (Subscript[min lim, u->0])f(u)=(Subscript[lim, ϵ->(SuperPlus[0])])r(ϵ).

The algorithm GCAD is described in:
A. Strzeboński. Solving Systems of Strict Polynomial Inequalities. J. Symbolic Comp. 29 (2000), 471–480.

Methods based on topological properties

These methods are based on the following characterization of the lower limit.

Lemma:
Let 
l=(Subscript[min lim, u->0])f(u),
S(y,x)=(h>0∧g<y h)∨(h<0∧g>y h),
Z={(a,u)∈(ℝ^(n+1)):S(a,u)},
Z(a)={u∈(ℝ^n):S(a,u)}.
Then 
(1) if  a>l  then 0∈(OverBar[Z(a)]),
(2) if  a<l  then (a,0)∉(OverBar[Z]).

[Z(a) is the set of u for which f(u)<a]

Regular projection

Let p∈ℚ[y,(Subscript[x, 1]),...,(Subscript[x, n])]. We say that p is regular in (Subscript[x, k]) if either p does not depend on (Subscript[x, k]) or the leading coefficient of p treated as a polynomial in (Subscript[x, k]) belongs to ℚ[y]. We say that a finite set P⊂ℚ[y,(Subscript[x, 1]),...,(Subscript[x, n])] is regular in (Subscript[x, k]) if all elements of P are regular in (Subscript[x, k]).

For a finite set of polynomials F⊂ℚ[y,(Subscript[x, 1]),...,(Subscript[x, k+1])] let GP(F, (Subscript[x, k+1])) ⊂ℚ[y,(Subscript[x, 1]),...,(Subscript[x, k])] denote the set of irreducible factors of the leading coefficients, discriminants, and pairwise resultants of the elements of F treated as polynomials in (Subscript[x, k+1]).

For a matrix M∈(ℤ^(n×n)) let (Subscript[M, i]) denote the i-th row of M, so that (Subscript[M, i])x is the i-th component of the vector M x.

Algorithm RPROJ:
Input:  A finite set F⊂ℚ[y,(Subscript[x, 1]),...,(Subscript[x, n])].
Output:  ((Subscript[P, 0]),...,(Subscript[P, n])) and an invertible matrix M∈(ℤ^(n×n)) such that
• (Subscript[P, n]) is the set of irreducible factors of F(y,(Subscript[M, 1])x,...,(Subscript[M, n])x),
• for 1≤k≤n, (Subscript[P, k]) is regular in (Subscript[x, k]) and (Subscript[P, k-1])=GP((Subscript[P, k]),(Subscript[x, k])).

Algorithm TLIM

Algorithm TLIM:
Input:  g∈Q[x],  h∈Q[x]\{0}.
Output: (Subscript[min lim, u->0])((g(u))/(h(u))).
(1)  Set F={g−y h,h,(Subscript[x, 1]),...,(Subscript[x, n])} and compute (((Subscript[P, 0]),...,(Subscript[P, n])),M)=RPROJ(F;y,(Subscript[x, 1]),...,(Subscript[x, n])).
(2) Compute the real roots  (Subscript[r, 1])<...<(Subscript[r, m]) of elements of (Subscript[P, 0]).
(3) Find rational numbers (Subscript[a, 0]),...,(Subscript[a, m]) such that -∞=(Subscript[r, 0])<(Subscript[a, 0])<(Subscript[r, 1])<(Subscript[a, 1])<...<(Subscript[r, m])<(Subscript[a, m])<(Subscript[r, m+1])=∞.
(3) Perform the binary search to find (Subscript[i, 0]): the minimal 0≤i≤m such that 0∈(OverBar[Z((Subscript[a, i]))])  or m + 1 if 0∉(OverBar[Z((Subscript[a, i]))]) for all 0≤i≤m.
(5) Return (Subscript[r, (Subscript[i, 0])]).

The algorithm introduces one variable y in addition to the variables x. The variable y is used only in one CAD projection computation needed to determine a finite number of possible values for the lower limit. The lower limit is then determined by a binary search among the possible values. The algorithm relies on the fact that (Subscript[min lim, u->0])f(u) is equal to one of the (Subscript[r, i]).

Depending on the method used to decide whether 0∈(OverBar[Z((Subscript[a, i]))])  we obtain three versions of the algorithm TLIM.

Outline of proof that (Subscript[min lim, u->0])f(u) is equal to one of the (Subscript[r, i]).

We show that for each  0≤i≤m either for all (Subscript[r, i])<b<(Subscript[r, i+1]) 0∈(OverBar[Z(b)])  or for all (Subscript[r, i])<b<(Subscript[r, i+1]) 0∉(OverBar[Z(b)]). To show this we perform a rational variable replacement y->σ(t) with σ:ℝ⟶((Subscript[r, i]),α] with an arbitrary fixed (Subscript[r, i])<α<(Subscript[r, i+1]). Let ((Subscript[Q, 0]),...,(Subscript[Q, n])) be the polynomial sets obtained from ((Subscript[P, 0]),...,(Subscript[P, n])) after the replacement and clearing of denominators. The leading coefficient of an element of (Subscript[Q, k]) is a univariate polynomial in t with no real roots. Hence ((Subscript[Q, 0]),...,(Subscript[Q, n])) is a valid CAD projection and a CAD Δ constructed using this projection is well-based in the sense of (Schwartz & Sharir, 1983). The set (Subscript[C, 0])=ℝ×{(0,...,0)} is a cell of Δ and the solution set W of S(σ(t),M x) is a union of cells of Δ. Since Δ is well-based, either (Subscript[C, 0])⊂(OverBar[W]) or (Subscript[C, 0])⋂(OverBar[W])=∅. In the first case 0∈(OverBar[Z(σ(a))])   for all a, in the second case  0∉(OverBar[Z(σ(a))])   for all a. This concludes the proof, since (Subscript[r, i])<α<(Subscript[r, i+1])  was arbitrary.

Deciding whether zero belongs to the closure of a semialgebraic set

Algorithm ZCQ1

Algorithm ZCQ1:
Input: A system S(x) of strict polynomial inequalities.
Output: true or false depending on whether 0 belongs to the closure of the solution set of S(x).
(1) If S(0) is true return true.
(2) Let T(ϵ,x):=∥x∥<ϵ∧S(x), and let (Subscript[F, (Subscript[x, n])]) be the set of all irreducible factors of polynomials that appear in T.
(3) For i=n-1 to 1 set (Subscript[F, (Subscript[x, i])]):=GP((Subscript[F, (Subscript[x, i+1])]),(Subscript[x, i+1])).
(4) Set (Subscript[F, ϵ]):=GP((Subscript[F, (Subscript[x, 1])]),(Subscript[x, 1])).
(5) Choose a positive rational number e smaller than all positive roots of elements of (Subscript[F, ϵ]).
(6) Run the lifting phase of the algorithm GCAD with the projection set ((Subscript[F, ϵ]),(Subscript[F, (Subscript[x, 1])]),...,(Subscript[F, (Subscript[x, n])])) for the open cell with the sample point (e).
[GCAD constructs only open cells]
(7) Return true if an open cell on which T(ϵ,x) is true is found, and return false otherwise.

The algorithm introduces a variable ϵ  in addition to the variables x, however it constructs only one cell in ϵ , constructs only open cells, and treats the variables x as existentially quantified. Moreover, the algorithm TLIM can call it with the original S((Subscript[a, i]),x), without transforming the variables x.

The algorithm GCAD is described in:
A. Strzeboński. Solving Systems of Strict Polynomial Inequalities. J. Symbolic Comp. 29 (2000), 471–480.

Algorithm ZCQ2

Algorithm ZCQ2:
Input: A system S(x) of polynomial inequalities.
Output: true or false depending on whether 0 belongs to the closure of the solution set of S(x).
(1) If S(0) is true return true.
(2) Compute a well-based CAD Δ of (ℝ^n) such that the solution set of S(x) is a union of cells of Δ.
(3) Use the algorithm CADAdjacency to compute adjacencies between cells of Δ.
(4) Return true if the cell containing 0 is adjacent to a cell on which S(x) is true, and return false otherwise.

The algorithm does not introduce variables, however it  requires the ability to compute a well-based CAD and hence the algorithm TLIM needs to call it with S((Subscript[a, i]),M x) as the input.

The algorithm CADAdjacency is described in:
A. Strzeboński.  CAD Adjacency Computation Using Validated Numerics. In Proceedings of the International Symposium on Symbolic and Algebraic Computation, ISSAC 2017. ACM, 413–420.

Comparison of the algorithms

We compare five CAD-based methods for computation of limits of real multivariate rational functions: OLIM1, OLIM2, TLIM1 (TLIM using ZCQ1), TLIM2 (TLIM using ZCQ2), and TLIM3 (TLIM using ZCQ2 if the matrix M is an identity matrix or a permutation matrix, and using ZCQ1 otherwise).

The motivation for the criterion used in TLIM3 is that ZCQ2 needs to use a transformed input, so ZCQ1 may be preferable if the linear transformation represented by M may destroy the sparsity of the input.

The algorithms have been implemented in Mathematica, as a part of the Limit command, using the built-in C implementation of the CAD algorithm.

Wolfram Alpha problems

In this example we used inputs selected from 120000 multivariate limit problems submitted to Wolfram Alpha. We have run the problems through the Mathematica built-in Limit and MinLimit commands, and selected the multivariate rational function limit inputs for which Mathematica was using a CAD-based method. Mathematica first tries to determine the limit by applying structural criteria and by computing univariate limits along curves. CAD-based methods are used only if the simpler methods fail. We have obtained 4818 distinct rational function inputs for which solving of Problem 2 required a CAD-based method. In 4793 of these at least one of the methods finished in under 600 seconds. For Problem 1 we obtained 612 inputs that required a CAD-based method. The large difference in the number of examples comes from the fact that Mathematica uses criteria that can prove non-existence of limit, but do not provide the value of the lower limit.

Timings for Problem 1:

Timings for Problem 2:

TLIM1 has the best average timing for 6 sets of Problem 1 examples and 3 sets of Problem 2 examples. TLIM3 has the best average timing for 2 sets of Problem 1 examples and 5 sets of Problem 2 examples. OLIM1 has the best average time in one set of Problem 2 examples. TLIM3 is at least as fast as TLIM2 for all sets of examples.

Timeout percentages for Problem 1:

Timeout percentages for Problem 2:

OLIM1 and OLIM2 timed out more frequently than TLIM1, TLIM2, and TLIM3, likely due to the fact that they require projection computation with one variable more.

Win counts for Problem 1:

Win counts for Problem 2:

TLIM1 and OLIM2 were the fastest methods for the largest number of Problem 1 examples. TLIM3 and OLIM1 were the fastest methods for the largest number of Problem 2 examples. However, each of the methods was the fastest for a significant number of examples of each problem.

Examples from literature

The functions in this example were selected from examples in

• P. Alvandi, M. Kazemi, and M. Moreno Maza. 2016. Computing Limits of Real Multivariate Rational Functions. In Proceedings of the International Symposium on Symbolic and Algebraic Computation, ISSAC 2016. ACM, 39–46.

• S. J. Xiao and G. X. Zeng. 2014. Determination of the limits for multivariate rational functions. Science China Mathematics 57, 2 (2014), 397–416.

using the same procedure as for Wolfram Alpha examples, that is functions for which Mathematica solves both Problem 1 and Problem 2 without using CAD-based methods were not selected. We obtained one instance of Problem 1 and six instances of Problem 2.

Here methods TLIM2 and TLIM3 were identical (because the matrix M was always an identity matrix or a permutation matrix) and were faster than all other methods.

Conclusions

Experimental results suggest that for problems with at least three variables algorithm TLIM3 gives the best average timings. However, since each of the methods (except TLIM2 when different than TLIM3) was the fastest for a significant number of examples, finding a good automatic method selection heuristic is a nontrivial problem for future research.