Table Of Content

Quantifier elimination for approximate BK-factorization E. Kartashova1, S. McCallum2, 1 RISC, J.Kepler University, Linz, Austria 7 0 2 Macquarie University, Sydney, Australia 0 2 e-mail: [email protected], [email protected] n a J 7 1 Introduction 1 v Factorization of linear partial differential operators (LPDOs) is a very well- 9 studied problem and a lot of pure existence theorems are known. The only 1 known constructive factorization algorithm - Beals-Kartashova (BK) factor- 0 1 ization - is presented in [1]). Its comparison with Hensel descent which is 0 sometimesregardedasconstructive, isgivenin[2],wheretheideatouseBK- 7 factorization for approximate factorization is also discussed. It originates in 0 / oneofthemostinterestingfeatures of BK-factorization: atthebeginningall h p the first-order factors are constructed and afterwards the factorization con- - dition(s) should be checked. This leads to the important application area - h t namely, numerical simulations which could be simplified substantially if in- a m stead of computation with one LPDE of order n we will be able to proceed computations with n LPDEs all of order 1. In numerical simulations it is : v not necessary to fulfill factorization conditions exactly but with some given i X accuracy, which we call approximate factorization. r The idea of the present paper is to look into the feasibility of solving a problems of this kind using quantifier elinination by cylindrical algebraic decomposition [3]. In this paper we are going to apply this approach to a hyperbolic LPDO of order 2 with polynomial coefficients. 2 Hyperbolic LPDO of order 2 A bivariate operator of second order has general form A = a ∂2+a ∂ ∂ +a ∂2+a ∂ +a ∂ +a (1) 2 20 x 11 x y 02 y 10 x 01 y 00 1 and is factorizable [1] iff ωa +a −L(2a ω+a ) ωa +a −L(2a ω+a ) 10 01 20 11 10 01 20 11 a = L + × 00 2a ω+a 2a ω+a (cid:26) 20 11 (cid:27) 20 11 a (a −L(a ω+a ))+(a ω+a )(a −La ) 20 01 20 11 20 11 10 20 × . (2) 2a ω+a 20 11 Here coefficients a = a (x,y) are functions on two variables x and y; i,j i,j ω is a distinct root of the following polynomial P (z) := a z2 + a z + 2 20 11 a , P (ω) = 0; and L is a linear differential operator of the form L = 02 2 ∂ −ω∂ . x y Let us introduce a function of two variables x,y ωa +a −L(2a ω+a ) ωa +a −L(2a ω+a ) 10 01 20 11 10 01 20 11 R= L + × 2a ω+a 2a ω+a (cid:26) 20 11 (cid:27) 20 11 a (a −L(a ω+a ))+(a ω+a )(a −La ) 20 01 20 11 20 11 10 20 × , 2a ω+a 20 11 and rewrite factorization condition (2) as a = R. 00 Now suppose that (1) is a hyperbolic operator in canonical form, i.e. H = ∂2−∂2+a ∂ +a ∂ +a (3) 2 x y 10 x 01 y 00 which corresponds to a (x,y) = 1, a (x,y) = 0, a (x,y) = −1. In this 20 11 02 case we have two roots ω = 1 and ω = −1, and function R takes a form 1 2 a ±a (a ±a )2 10 01 10 01 R= L + , 2 4 (cid:26) (cid:27) where ”+” corresponds to ω = 1 and ”−” corresponds to ω = −1. We 1 2 rewrite (4) is a slightly different form which will more convenient for further use: (a +a )/2, ω = 1; R = L{S}+S2 with S = 10 01 (4) ((a10−a01)/2, ω = −1. 3 Polynomial coefficients Let us suppose that operator H has polynomial coefficients a and regard 2 ij cases. 3.1 Polynomials of first degree We have 3polynomials a of firstdegree withtwo variables x,y: a (x,y) = ij 00 b x+b y+b , a (x,y) = c x+c y+c , a (x,y) = d x+d y+d , then 3 2 1 10 3 2 1 01 3 2 1 2 (1) For the first root ω = 1 we have L = ∂ −∂ and 1 x y a +a = (c +d )x+(c +d )y+(c +d ) = f x+f y+f 10 01 3 3 2 2 1 1 3 2 1 withf =(c +d ), f = (c +d ), f =(c +d ).ThenL(a +a )= f −f 1 1 1 2 2 2 3 3 3 10 01 3 2 and f −f (f x+f y+f )2 (1) 3 2 3 2 1 R = + (5) 1 2 4 (2) For the second root ω = −1 we have L = ∂ +∂ and 2 x y a −a = (c −d )x+(c −d )y+(c −d )= h x+h y+h 10 01 3 3 2 2 1 1 3 2 1 with h = (c −d ), h = (c −d ), h = (c −d ). Then L(a −a ) = 1 1 1 2 2 2 3 3 3 10 01 h −h and 3 2 h −h (h x+h y+h )2 (1) 3 2 3 2 1 R = + . (6) 2 2 4 Remark 1. Notice that functions R and R coincide symbolically: 1 2 s −s (s x+s y+s )2 R(1) = 3 2 + 3 2 1 (7) 2 4 but of course, the form of s as functions of coefficients c , b , d will be i j j j different. Remark2. Factorization conditionfortheoperator(3)hasnowverysimple form s −s (s x+s y+s )2 R(1) = a (x,y) ⇒ 3 2 + 3 2 1 = b x+b y+b 00 3 2 1 2 4 which yields s = 0, s s = 0, s = 0, s s = 2b , 3 2 3 2 3 1 3 (s2s1 = 2b2, s21+2(s3−s2) = 4b1 For instance, in case of the second root this system of equations has form (c −d )2 =0, (c −d )(c −d ) = 0, (c −d )2 = 0, 3 3 3 3 2 2 2 2 (c −d )(c −d )= 2b , (c −d )(c −d )= 2b ,  3 3 1 1 3 2 2 1 1 2 (c −d )2+(2(c −d )−(c −d )) = 4b 1 1 3 3 2 2 1  and its solution gives all exactly factorizable operators of this type: a (x,y) = (c −d )2/4 00 1 1 a (x,y) = c x+c y+c  10 3 2 1 a (x,y) = c x+c y+d 01 3 2 1   3 3.2 Polynomials of second degree Now we have 3 polynomials a of second degree with two variables x,y: ij a (x,y) = b x2+b xy+b y2+b x+b y+b , 00 6 5 4 3 2 1 a (x,y) = c x2+c xy+c y2+c x+c y+c , 10 6 5 4 3 2 1 a (x,y) = d x2+d xy+d y2+d x+d y+d , 01 6 5 4 3 2 1 then 1. For the first root ω = 1 we have L = ∂ −∂ and 1 x y a +a = (c +d )x2+(c +d )xy+(c +d )y2+(c +d )x+(c +d )y+(c +d ) 10 01 6 6 5 5 4 4 3 3 2 2 1 1 = f x2+f xy+f y2+f x+f y+f 6 5 4 3 2 1 with f = (c +d ), ∀i= 1,2,...,6. Then i i i L(a +a )= 2(f x−f y)+f (y−x)+f −f 10 01 6 4 5 3 2 and 2(f x−f y)+f (y−x)+f −f (2) 6 4 5 3 2 R = + 1 2 (f x2+f xy+f y2+f x+f y+f )2 6 5 4 3 2 1 + (8) 4 2. For the second root ω = −1 we have L = ∂ +∂ and 2 x y a −a = (c −d )x2+(c −d )xy+(c −d )y2+(c −d )x+(c −d )y+(c −d ) 10 01 6 6 5 5 4 4 3 3 2 2 1 1 = h x2+h xy+h y2+h x+h y+h 6 5 4 3 2 1 with h = (c −d ), ∀i= 1,2,...,6. Then i i i L(a −a )= 2(h x−h y)+h (y−x)+h −h 10 01 6 4 5 3 2 and 2(h x−h y)+h (y−x)+h −h (2) 6 4 5 3 2 R = + 2 2 (h x2+h xy+h y2+h x+h y+h )2 6 5 4 3 2 1 + . (9) 4 As above we have in fact one function 2(s x−s y)+s (y−x)+s −s R(2) = 6 4 5 3 2 + 2 (s x2+s xy+s y2+s x+s y+s )2 6 5 4 3 2 1 + . (10) 4 4 Remark 3. As direct corollaries of linear differentiation of polynomials one can conclude that for a polynomial any finite degree n, function R(n) has form n−1 R(n) = R(k)+(other terms) k=1 X and is a polynomial of degree r with r ≤ 2n, with necessary condition of exact factorization being deg(a )≤ deg(R(n)). 00 4 Problem setting for QE We use standard formal language of elementary real algebra, that is, Tarski algebra [3] and formulate a first simple case of approximate factorization of LPDO as a quantifier elimination (QE) problem. Namely, we shall consider the case of a hyperbolicLPDO of order 2 with polynomial coefficients of the first order. In this case we have (1) three polynomials a (x,y) of first degree in the two variables x,y; ij the coefficients of these polynomials are also variables: a (x,y) = b x+b y+b 00 3 2 1 a (x,y) = c x+c y+c 10 3 2 1 a (x,y) = d x+d y+d ; 01 3 2 1 (2) one function s −s (s x+s y+s )2 R(1)(x,y) = 3 2 + 3 2 1 2 4 withs givenby(5)orby(6),i.e. s = c +d forthefirstrootands = c −d i i i i i i i for the second root; (3) a constant ε; (4) constants m and n, which definea boundedrectangular region in the plane: −m< x < m, −n< y <n. Remark 4. Notice that the special form of the factorization condition allowed us to reduce the number of variables needed for this QE problem. Initially we had 9 variables b ,b ,b ,c ,c ,c ,d ,d ,d , but in fact it 3 2 1 3 2 1 3 2 1 is enough to consider only the 6 variables s ,s ,s ,b ,b ,b . 1 2 3 1 2 3 With all this given, let us consider the quantified formula of elementary ∗ ∗ realalgebraφ = φ (b ,s )whichassertsthat“forallxandyinthebounded i j region −m < x < m, −n < y < n, we have −ε< a (x,y)−R(1)(x,y) < ε.” 00 ∗ We wishto eliminate thequantifiers from φ (b ,s ). More precisely, we wish i j ′ ′ tofindaformulaofelementaryrealalgebraφ = φ(b ,s ),freeofquantifiers, i j ′ ∗ such that if φ(b ,s ) is true then φ (b ,s ) is true. That is, we wish to find i j i j 5 conditions on thecoefficients of the initial polynomialsa (x,y) which imply ij thatthefunctionR(1)(x,y)differsnottoomuchfromonethesepolynomials, namely a (x,y), throughout the bounded region −m < x < m, −n < y < 00 n. 5 Synopsis of QE by CAD Let A be a set of integral polynomials in x ,x ...,x , where r ≥ 1. An 1 2 r A-invariant cylindrical algebraic decomposition (CAD)of Rr, r-dimensional real space, is a decomposition D of Rr into nonempty connected subsets called cells such that 1. the cells of D are cylindrically arranged with respect to the variables x ,x ,...,x ; 1 2 r 2. every cell of D is a semialgebraic set (that is, a set defined by means of boolean combinations of polynomial equations and inequalities); and 3. every polynomial in A is sign-invariant throughout each cell of D. TheCADalgorithm as originally conceived [3,4] has inputsandoutputs as follows. Given such a set A of r-variate polynomials and a nonnegative integer f with f < r, the algorithm produces as its output a description of an A-invariant CAD D of Rr, in which explicit semialgebraic defining formulas are provided only for the cells of the CAD D of Rf induced (that f is, implicitly determined) by D. The description of D comprises lists of indices and sample points for the cells of D. (Every cell is assigned an index which indicates its position within the cylindrical structure of D.) The working of the original CAD algorithm can be summarized as follows. If r = 1, an A-invariant CAD of R1 is constructed directly, using polynomial real root isolation. If r > 1, then the algorithm computes a projection set P of (r −1)-variate polynomials (in x1,...,xr−1) such that any P-invariant CAD D′ of Rr−1 can be extended to a CAD D of Rr. If ′ ′ f = r we set f ← f−1 and otherwise set f ← f. Then the algorithm calls ′ ′ ′ itself recursively on P and f to get such a D . Finally D is extended to D. In order to produce semialgebraic defining formulas for the cells of D the f algorithm must be used in a mode called augmented projection. Thus for r > 1, if we trace the algorithm, we see that it computes a first projection set P, eliminating x , then computes the projection of P, r eliminating xr−1, and so on, until the (r − 1)-st projection set has been obtained, which is a set of polynomials in the variable x only. This is 1 called the projection phase of the algorithm. The construction of a CAD of R1 invariant with respect to the (r−1)-st projection set is called the base phase. The successive extensions of the CAD of R1 to a CAD of R2, the CAD of R2 to a CAD of R3, and so on, until an A-invariant cad of Rr is obtained, constitute the extension phase of the algorithm. Nowweconsiderthequantifier elimination (QE)problemfortheelemen- 6 tary theory of the reals: given a quantified formula (known as a QE problem instance) of elementary real algebra ∗ φ = (Q x )...(Q x )φ(x ,...,x ) f+1 f+1 r r 1 r where φ is a formula involving the variables x ,x ,...,x which is free of 1 2 r ′ ′ quantifiers, find a formula φ(x ,...,x ), free of quantifiers, such that φ is 1 f ∗ equivalent to φ . The QE problem can be solved by constructing a certain CAD of Rr. The method is described as follows. 1. Extract from φ the list A of distinct non-zero r-variate polynomials occurring in φ. 2. Constructlists S and I of sample points and cell indices, respectively, for an A-invariant CAD D of Rr, together with a list F of semialgebraic defining formulas for the cells of the CAD D of Rf induced by D. f ∗ 3. Using S, evaluate the truth value of φ in each cell of D . (By f ∗ construction of D, the truth value of φ is constant throughout each cell c ∗ of D , hence can be determined by evaluating φ at the sample point of c.) f ′ 4. Constructφ(x ,...,x ) as the disjunction of the semialgebraic defin- 1 f ∗ ing formulas of those cells of D for which the value of φ has been deter- f mined to be true. TheabovealgorithmsolvesanygivenparticularinstanceoftheQEprob- leminprinciple. Howeverthecomputingtimeofthealgorithmgrowssteeply as the number r of variables occurring in the input formula φ increases. Collins and Hong[7]introduced themethod ofpartial CAD construction for QE. This method, named with the acronym QEPCAD, is based upon the simple observation that we can often solve a QE problem by means of a partially built CAD. The QEPCAD algorithm was originally implemented by Hong. A recent implementation, denoted by QEPCAD-B, contains im- provements byBrown, Collins,McCallum, andothers–see[5]. QEPCAD-B has solved a range of reasonably interesting problems for which the original QE algorithm takes too much time. Nevertheless the worst case computing time of QEPCAD-B remains large (that is, it depends doubly-exponentially on r). 6 Application of QEPCAD to BK-factorization Weconsideronlythefirstsimplecaseofapproximatefactorization described in Section 4. Using the notation of Section 4, we suppose that ε, m and n havebeengivenspecificconstantvalues,sayε = m = n = 1,andweconsider ∗ the formula φ (b ,s ) which asserts that i j (∀x)(∀y)[(|x| < 1∧|y|< 1) ⇒ |a (x,y)−R(1)(x,y)| < 1]. (11) 00 7 ′ ′ We wish to find a formula φ(b ,s ), free of quantifiers, such that φ(b ,s ) i j i j ∗ implies φ (b ,s ). i j Remark 5. Itwouldbeofgreatestinteresttofindthemostgeneralsuch ′ ′ ∗ φ(b ,s ) – that is, to find quantifier-free φ(b ,s ) equivalent to φ (b ,s ). i j i j i j But as we’ll see it seems that the time and space resources needed to do this areprohibitive. We’ll also seethatitisnotastimeconsuming, yet hopefully ∗ still of interest, to find quantifier-free conditions merely sufficient for φ to be true. We attempted to find a solution to the above QE problem instance by running the program QEPCAD-B with the quantified formula 11 (rewrit- ten so that the variables b ,s appear explicitly, and the denominator 4 is i j cleared from the right hand side of the implication, see 12 below) as its input. The variable ordering used was (s ,s ,s ,b ,b ,b ,x,y). The com- 3 2 1 3 2 1 puter used for this and subsequent experiments was a Sun server having a 292 MHz ultraSPARC risc processor. Forty megabytes of memory were made available for list processing. However the program ran out of memory after approximately one hour and forty minutes. The program was exe- cuting the projection phase of the algorithm when it stopped. The first three projection steps – that is, successive elimination of y, x and b – were 1 complete. Increasing the amount of memory to eighty megabytes did not help – the program still ran out of memory during the fourth projection step (that is, during elimination of b ). 2 6.1 Searching for quantifier-free sufficient conditions Of course a very special, but completely trivial, quantifier-free sufficient condition for our QE problem instance is the formula ′ φ(b ,s ) := [b = 0∧b = 0∧b = 0∧s = 0∧s = 0∧s = 0]. i j 1 2 3 1 2 3 Itcouldbeofsomeinteresttolookforpartialsolutionsto(thatis,quantifier- free sufficient conditions for) our QE problem instance in which some but not all of the variables b ,s are equal to zero. For example, recall that i j the given quantified (11) – after rewriting so that the variables b ,s appear i j explicitly and the denominator 4 is cleared from the right hand side of the implication – is: (∀x)(∀y)[(|x| < 1∧|y| < 1) ⇒ |4b x+4b y+4b −2(s −s )−(s x+s y+s )2| < 4]. 3 2 1 3 2 3 2 1 (12) Suppose that we put b = s = 0 in (12). We obtain: 2 2 (∀x)(∀y)[(|x| < 1∧|y|< 1) ⇒ |4b x+4b −2s −(s x+s )2|< 4] 3 1 3 3 1 8 which is equivalent to: (∀x)[(|x| < 1) ⇒ |4b x+4b −2s −(s x+s )2|< 4], (13) 3 1 3 3 1 ∗ which we shall denote by ψ (b ,s ). i j The following theorem shows that a partial solution to the special QE ∗ problem instance ψ (b ,s ) (that is, a quantifier-free sufficient condition for i j ∗ ∗ ψ ) leads to a partial solution to the QE problem instance φ (that is, a ∗ quantifier-free sufficient condition for φ ). ′ Theorem 1 Suppose that ψ (b ,s ) is a quantifier-free formula, involving i j ∗ only b ,b ,s ,s , which implies ψ (b ,s ). Then the quantifier-free formula 1 3 1 3 i j ′ ∗ ψ (b ,s )∧b = 0∧s = 0 implies φ (b ,s ). i j 2 2 i j ◮ Let b ,s be real numbers. Assume ψ′(b ,s )∧b = 0∧s = 0. Then i j i j 2 2 ∗ ψ (b ,s )∧b = 0∧s = 0 is true, by hypothesis. Take real numbers x and i j 2 2 y, with |x| < 1 and |y|< 1. Then |4b x+4b y+4b −2(s −s )−(s x+s y+s )2|= |4b x+4b −2s −(s x+s )2|< 4, 3 2 1 3 2 3 2 1 3 1 3 3 1 by virtue of (13) (since |x| < 1). Hence (12) is true. (cid:4) The above discussion suggests that it would be worthwhile to try to ∗ find a solution to the simplified, special QE problem instance ψ using the program QEPCAD-B. Putting (13) into a slightly more general form, and hence reducing by 1 the number of variables in the formula, we obtain: (∀x)[(|x| < 1) ⇒ |ax2+bx+c| < 4]. (14) ′ Apartialsolution θ (a,b,c) to(14)couldeasily betransformedintoapartial solution ψ′(b ,b ,s ,s ) to ψ∗ by setting a = −s2, b = 4b − 2s s and 1 3 1 3 3 3 1 3 c= 4b −2s −s2. 1 3 1 We ran program QEPCAD-B with (14) as its input. Eighty megabytes of memory were made available for list processing. After 191 seconds the program produced the following quantifier-free formula equivalent to (14): c - b + a + 4 >= 0 /\ c - b + a - 4 <= 0 /\ c + b + a + 4 >= 0 /\ c + b + a - 4 <= 0 /\ [ 4 a c - b^2 + 16 a > 0 \/ 4 a c - b^2 - 16 a > 0 \/ [ b^2 - 16 a = 0 /\ b^2 + 16 a > 0 ] \/ [ b^2 - 16 a < 0 /\ b - 2 a >= 0 ] \/ [ b^2 - 16 a < 0 /\ b + 2 a <= 0 ] \/ [ b^2 - 16 a > 0 /\ b + 2 a >= 0 ] \/ [ b^2 - 16 a > 0 /\ b - 2 a <= 0 ] \/ [ b^2 - 16 a = 0 /\ c - b + a + 4 > 0 /\ c - b + a - 4 < 0 ] ]. Since a= −s2, we have a≤ 0. We ran QEPCAD-B a second time, this time 3 using the command 9 assume [a <= 0]. After 60 seconds the program produced the following somewhat simpler quantifier-free formula equivalent to (14) under the assumption a ≤ 0: c - b + a + 4 >= 0 /\ c - b + a - 4 <= 0 /\ c + b + a + 4 >= 0 /\ c + b + a - 4 <= 0 /\ [ 4 a c - b^2 - 16 a > 0 \/ [ b > 0 /\ b + 2 a >= 0 ] \/ [ b < 0 /\ b - 2 a <= 0 ] \/ [ b^2 + 16 a = 0 /\ c - b + a + 4 > 0 /\ c - b + a - 4 < 0 ] ]. It is possible to induce the program to produce an arguably even simpler solution formula using less computing time by making two separate runs of QEPCAD-B. The first run uses the command assume [a < 0]. After just 1.9 seconds the program produced the following quantifier-free formula equivalent to (14) under the assumption a< 0: c - b + a + 4 >= 0 /\ c - b + a - 4 <= 0 /\ c + b + a + 4 >= 0 /\ c + b + a - 4 <= 0 /\ [ b - 2 a <= 0 \/ b + 2 a >= 0 \/ 4 a c - b^2 - 16 a > 0 ]. (15) The above formula is perhaps the most elegant and understandable of those obtained by applying QEPCAD-B to Formula 14. For it is a slight improve- mentof(thatis,slightly morecompactthan)aformulaseentobeequivalent to it (under assumption a < 0) which is quite straightforward to derive by hand from (14) using elementary properties of the parabola y = ax2+bx+c on the interval (−1,+1): [ 2 a - b >= 0 /\ a + b + c + 4 >= 0 /\ a - b + c - 4 <= 0] \/ [2 a + b >= 0 /\ a - b + c + 4 >= 0 /\ a + b + c -4 <= 0] \/ [2 a - b < 0 /\ 2 a + b < 0 /\ 4 a c - b^2 - 16 a > 0 /\ a - b + c + 4 >= 0 /\ a + b + c + 4 >= 0]. (16) Remark 6. To derive by hand (16) from (14) under the assumption a < 0, one has to notice that function f(x) = ax2 + bx + c has its max- ′ imum value for f (x) = 2ax + b = 0, that is, for x = −b/(2a), and consider three cases separately: (1) −b/(2a) ≤ −1, (2) −b/(2a) ≥ +1, and (3) −1 < −b/(2a) < +1. For each of the above three cases one can then write down necessary and sufficient conditions for (14) to be true. For example, in Case 1, (14) is clearly equivalent to −4≤ a+b+c∧a−b+c ≤ 4. After treating each of the above cases, we obtain (16) by forming the disjunction of the formulas corresponding to the cases. Of course, (15) for a < 0 is not quite a complete solution to the QE problem instance of (14) under assumption a ≤ 0. To obtain a complete 10