Creative Telescoping on Multiple Sums
CCreative Telescoping on Multiple Sums
Christoph Koutschan and Elaine Wong
Abstract.
We discuss the strategies and difficulties of determining arecurrence which a certain polynomial (in the form of a symbolic mul-tiple sum) satisfies. The polynomial comes from an analysis of integralestimators derived via quasi-Monte Carlo methods.
Keywords.
Symbolic Summation, Creative Telescoping, Holonomic Func-tion, Hypergeometric Series.
1. The Problem
Recently, Wiart and Wong [6] derived a formula for the covariance of anintegral estimator for functions satisfying a certain decay condition, based ona quasi-Monte Carlo framework developed by Wiart, Lemieux, and Dong [5].This formula is written as the following polynomial in x , G s ( x ) ∶= m + s − ∑ k = ⎛⎝ s ∑ r = ( sr )( k − r − ) b − (− b ) r r − − c m ( k ) ∑ i = (− b ) i ( r − i )⎞⎠( bx ) k , (1.1)where c m ( k ) = max ( k − m, ) . The goal in [6] is to show that (1.1) is notpositive for all b, m, s ∈ N , b ≥ x ∈ [ , ) . We approach this problemby employing symbolic computation rather than analysis: using the avail-able computer algebra software for holonomic functions [1, 2, 4], we carryout a guess-and-prove strategy that ultimately leads us to deduce a suitableclosed form for (1.1). The result is an expression in terms of regularized betafunctions which allows us to show the desired nonpositivity statement. Thisarticle serves to highlight the “proving” aspect of the strategy, i.e., the deriva-tion and proof of a third-order recurrence that (1.1) satisfies. At the sametime, we discuss some general solutions to address difficulties in summationproblems of a similar nature. Both authors were supported by the Austrian Science Fund (FWF): F5011-N15. a r X i v : . [ c s . S C ] O c t C. Koutschan and E. Wong
2. Motivation and Background
On a first glance, we note that all constituents of (1.1) have the propertyof being “holonomic”. For the purposes of this paper, we stay slightly in-formal (but still rigorously practical) and use the definition that holonomicfunctions are sequences which satisfy “sufficiently many” linear recurrenceswith polynomial coefficients. One convenience of dealing with such functionscomes from the fact that holonomicity is preserved through basic operations:we will refer to these as “closure properties”. For example, the product of twoholonomic functions is again holonomic [7, Proposition 3.2]. From a compu-tational point of view, if we know the recurrences that the two holonomicfunctions satisfy, we can construct a recurrence that their product satisfiesand even bound its order and the degrees of its polynomial coefficients.Our expression (1.1) consists of summation quantifiers, (products of)binomial coefficients, and polynomial/exponential functions in the parame-ters. The binomial coefficient, for example, can be described completely viatwo recurrences and finitely many initial conditions. From our usual notionof binomial coefficients, we can immediately write down the two recurrences(valid on Z × Z ): ( n − k + )( n + k ) − ( n + )( nk ) = , ( k + )( nk + ) − ( n − k )( nk ) = , (2.1)and specify the initial conditions ( ) = , (− ) = , (− − ) = . By doing so, we interpret the binomial coefficient over the integers in thetraditional (combinatorial) way, namely, it is nonzero only for 0 ≤ k ≤ n . Weexplicitly highlight this because the computer algebra system Mathematicauses a more general notion of the binomial coefficient, which extends itsdefinition to the negative integers. It is important to note that this generalizedbinomial coefficient is defined by the very same recurrences, but just usesdifferent initial conditions: ( ) = , (− ) = , (− − ) = . Unfortunately, the more general view of the binomial coefficient affects thenatural boundaries of our summation limits: the summands containing thesecoefficients behave inappropriately outside of our prescribed bounds, whichis of course irrelevant regarding the definition of G s ( x ) , but which may causeproblems when evaluating boundary terms originating from the telescopingsums. We will have to address this issue later.In the following, we will use this simple bivariate sequence ( nk ) to il-lustrate the main features of the holonomic systems approach. Using thenotation S n resp. S k to represent the forward shift operator in the givenreative Telescoping on Multiple Sums 3variable, we can rewrite (2.1) so that each of the corresponding operators ( n − k + ) S n − ( n + ) , ( k + ) S k − ( n − k ) , (2.2)maps ( nk ) to the zero sequence. We say that these operators annihilate thegiven function. As one can see, the translation between recurrence and oper-ator can be read off immediately. Viewing recurrences as operators enablesthe use of algebraic methods to manipulate them more efficiently. However,we then have to deal with objects that do not always commute. The appro-priate algebraic framework to represent such operators is an Ore algebra. Inthe following technical definition, ∂ serves as a placeholder for any of ouroperator symbols S n or S k . Definition 2.1.
Let R be a ring.1. If σ ∶ R → R is a ring endomorphism and δ ∶ R → R is such that the “skew”Leibniz law is satisfied, that is, δ ( f g ) = δ ( f ) g + σ ( f ) δ ( g ) for all f, g ∈ R , then δ is called a σ -derivation.2. Suppose now that there is an endomorphism σ ∶ R → R and a σ -derivation δ ∶ R → R . Suppose further that a ring structure is defined on the set R [ ∂ ] of univariate polynomials in ∂ with coefficients in R , equipped with theusual addition, and multiplication is such that ∂ i ∂ j = ∂ i + j and ∂f = σ ( f ) ∂ + δ ( f ) for all i, j ∈ N and f ∈ R. Then R [ ∂ ] is an Ore algebra over R . We typically use the symbol O todenote such algebras.3. Suppose that f is in a left- R [ ∂ ] -module with action ⋅ ∶ O × R → R suchthat 1 ⋅ f = f and L ⋅ ( L ⋅ f ) = ( L L ) ⋅ f for all L , L ∈ O . Then wecall ann ( f ) = { L ∈ R [ ∂ ] ∣ L ⋅ f = } the annihilator of f in R [ ∂ ] .For the binomial coefficient ( nk ) , the Ore algebra that we use is O = R [ S k , S n ] with R = Q ( n, k ) . A quick check shows that shift operators satisfythe required commutation properties. In the definition of this Ore algebra,each σ denotes a forward shift operation (clearly a ring endomorphism) and δ ≡ σ -derivation).The reader may have wondered why we did not consider Pascal’s ruleas a potential defining recurrence for the binomial coefficient. The reason isthat (2.1) are quite canonical generators for the set of all recurrences satisfiedby ( nk ) . In algebraic terms, we can formulate this statement precisely: the an-nihilator of ( nk ) , which is a left ideal in O , is generated by the operators (2.2),let’s call them P , P ∈ O :ann ( nk ) = { C ⋅ P + C ⋅ P ∣ C , C ∈ O } . C. Koutschan and E. WongMoreover, the two operators in (2.2) even form a (left) Gr¨obner basis ofann ( nk ) .The key tool used in rigorously deriving a “grand” recurrence for G s ( x ) lies in the highly touted creative telescoping algorithm [8] for symbolic sumsand integrals, as implemented in the HolonomicFunctions package [2]. Inorder to construct a recurrence for a symbolic parametric sum of the form ∑ k summand , (2.3)the algorithm takes as input a list of generators, like the ones in (2.2), for anannihilating ideal of the summand. If the summand is given as a closed-formexpression, then such a list is automatically computed, provided that it isrecognized to be holonomic. The algorithm then identifies lists of operators P and Q (in the form of Ore polynomials in the algebra as described inDefinition 2.1) such that for each P ∈ P and its corresponding Q ∈ Q , theoperator P − ( S k − ) ⋅ Q is an element of the given annihilating ideal. Theset P contains the so-called “telescopers” (all of which are free of k and S k but may contain the other parameters), and the set Q the corresponding“certificates”.How do these objects help us? Summing with respect to k gives relationsof the form ∑ k P ⋅ summand − ∑ k ( S k − ) ⋅ Q ⋅ summand = . (2.4)In a best-case scenario, each of the P ’s commutes with the first summationin (2.4) (allowing us to pull it out of the sum so that we can view the elementsof P being applied to the whole sum and not just the summand) and thesecond summation collapses to zero by telescoping (leaving no trace of thecertificate). From there, we would conclude that P generates a left ideal ofannihilating operators for (2.3), that is, it represents a set of recurrenceswhich are satisfied by the sum. Coming back to our triple sum (1.1), we canrepeatedly apply this process until a recurrence for the outermost (and hencethe whole) sum is deduced.However, life is not always that easy: during the application of thisstrategy to the particular summation problem (1.1), we encountered the fol-lowing difficulties that are somewhat prototypical for the holonomic systemsapproach. This explains why, despite being automatable in principle, it stilllacks a press-the-button implementation that would provide a computer proofof a claimed identity in a completely automatic way and without any humaninteraction.1. The summand, i.e., the expression inside a summation quantifier, maytake on nonzero values outside of the respective summation bounds.Thus, there is no reason to expect a priori that the second summationin (2.4) will evaluate to zero. And indeed, we found that it did not,and such terms constitute some of the “inhomogeneous parts” of theequation. An additional annihilator for them is required in order tohomogenize the recurrence.reative Telescoping on Multiple Sums 52. The upper boundaries contain the variable s , and the operators in P contain shifts in s , causing difficulties with moving P ∈ P to be outsideof the sum.3. Some of the operators in Q contain singularities at the boundary valuesso we were forced to exclude these values (which required compensationelsewhere). This is because the sum in (2.4) containing the certificatesis designed to collapse to only boundary value evaluations and we en-counter problems if the summands are undefined at such values. Furtherissues could surface if those summands were also undefined at some in-termediary value. Luckily, this was not the case here.4. Mathematica, in its symbolic zeal, rewrites the innermost sum in (1.1)as a hypergeometric F series and the second innermost sum as a Dif-ferenceRoot. While the values of the F function match with our sumwithin the domain in question, there are still an infinite number of valuesfor which it doesn’t. The DifferenceRoot is Mathematica’s version of arecurrence together with initial values, but unfortunately not helpful forour purposes because it is incompatible with HolonomicFunctions anddoes not support multivariate recurrences that are needed for creativetelescoping.We illustrate the first three difficulties with a toy example. For thor-oughness, we apply our strategy fully to this example to give a sufficient ideaof the broader behavior. From this point on, we will periodically perform theservice of demonstrating how computers and humans interact, by highlighting(in brackets) when paper-and-pencil reasoning is used and when automationis applied.Suppose we want to rigorously determine an annihilating operator for n ∑ k = ( nk ) (2.5)for n ≥
5. In other words, we would like to identify a recurrence that it satisfies.The creative telescoping algorithm (computer) outputs the telescoper P = S n − Q = kk − n − . Then (2.4) implies n ∑ k = ( S n − )( nk ) − ⎛⎝ kk − n − ( nk )RRRRRRRRRRR k = n + − kk − n − ( nk )RRRRRRRRRRR k = ⎞⎠·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ collapsed sum with singularity at k = n + = , with the summation containing the certificate collapsing to only evaluationsat the boundary values. We note that ( nk ) is nonzero for k = , . . . ,
4, i.e.,outside of the summation bounds. After substituting k = ≤ k = n + n − k = n rather than k = n + n − ∑ k = ( S n − )( nk ) − ( n − n + n − n )·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ inhomogeneous part = . (2.6)This fixes the issue of the singularity (alternatively, we could have rewritten kk − n − ( nk ) = − kn + ( n + k ) to get rid of the pole). Next, we note that the uppersummation bound contains the parameter n while our telescoper P containsthe shift operator S n : applying the operator to the whole sum affects boththe upper bound and ( nk ) . This is fixed with the (human) observation that n − ∑ k = ( S n − )( nk ) = ( S n − ) n ∑ k = ( nk ) −( n + n ) − ( n + n + ) + ( nn )·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ compensated terms = − n . In this situation, we say that the operator and the summation does notcommute. Then (2.6) simplifies to the inhomogeneous recurrence ( S n − ) n ∑ k = ( nk ) = n − n + n − n . If one prefers a homogeneous recurrence, an annihilating operator for theright-hand side can be determined to be ( n − ) S n − ( n + ) . We can thereforeconclude that (( n − ) S n − ( n + )) ⋅ ( S n − ) = ( n − ) S n + ( − n ) S n + ( n + ) is the desired annihilating operator for (2.5). As is expected in such simpleexamples, the recurrence can be solved (by the computer) to obtain a closedform for (2.5). Of course, it agrees with the one that one directly gets frominvoking the binomial theorem.In the above example, we can see that there is a “dance” between humanand the computer and only upon a careful collaboration does it bear fruit.We now proceed to use a similar strategy to attack the big sum G s ( x ) andfurthermore present some alternatives to improve performance. The totalcomputation time largely depends on how complicated the summands andinhomogeneous parts turn out to be after (human) simplification. The nextsection outlines some of these strategies and in particular highlights how wewere able to successfully derive (and prove) a recurrence for G s ( x ) .
3. A Playbook for the Holonomic Approach
This section illustrates how to generally overcome the difficulties listed in theprevious section and how to effectively perform the human-computer danceto prove our main result. We envision that the discussion leads to a deeperunderstanding of the practical issues when applying the holonomic systemsapproach and makes it accessible for future applications. The Mathematicanotebook containing implementations of these strategies can be found in theonline supplementary material [3].reative Telescoping on Multiple Sums 7
Theorem 3.1.
For b, m, s ∈ N , b ≥ , the polynomial given in (1.1) , G s ( x ) ∶= m + s − ∑ k = ⎛⎝ s ∑ r = ( sr )( k − r − ) b − (− b ) r r − − c m ( k ) ∑ i = (− b ) i ( r − i )⎞⎠( bx ) k , satisfies the recurrence ( s + )( bx − ) ⋅ G s + + ( m ( bx − )( x − ) + bsx ( x − ) + bx ( x − ) − s ( x − ) − x + ) ⋅ G s + − ( x − )( bmx + bsx + bx + mx − m + sx − s + x − ) ⋅ G s + + ( x − ) ( m + s + ) ⋅ G s = . We note that this result is already contained in [6, Lemma 15] withcomputational details found in [3]. The following discussion serves to outlinealternate (and in some cases, faster) proof strategies, to provide some expo-sition for technical details that were not mentioned in [6], and to explicitlyresolve some of the issues mentioned in the previous section in as much gen-erality as possible under the context of using our problem as a case study.We hope that this will be useful for future practitioners.
Before we dive in, we make a few remarks about how to view G s ( x ) to makeour life easier. On the one hand, the summation (1.1) can be separated intotwo parts G s = G ( ) s + G ( ) s , in order to remove the max function in the upperlimit of the innermost sum. After a mild simplification (human), these twoparts look as follows: G ( ) s ∶= − m + s − ∑ k = s ∑ r = ( sr )( k − r − ) ( b − b ) r ( bx ) k ,G ( ) s ∶= m + s − ∑ k = m + s ∑ r = ( sr )( k − r − ) − b (− b ) r r − ∑ i = r −( k − m ) (− b ) i ( r − i )( bx ) k . Observe that − G ( ) s is the collection of terms that is added to G s to enablethe sum to collapse to G ( ) s . We write this out to show that initially applyingthe full strategy to “only” the double sum gives a hypothetical lower boundfor the time and effort required to treat the whole thing. We note that if k > m + s −
1, then there is no reason to expect that either of the inner sumswould be zero, which may cause the inhomogeneous parts in (2.4) to survive.The split sums also serve as an example of how to apply closure prop-erties: the sum of holonomic functions is still holonomic [7, Proposition 3.1],so ann ( G ( ) s + G ( ) s ) can be deduced by executing (computer) the corresponding “closure propertyof addition” algorithm after separately computing a respective annihilatingideal for each of the two terms. This closure property can also be applied C. Koutschan and E. Wongin intermediate computations (for example, during the treatment of the in-homogeneous parts). However, the user should be aware that there is a riskthat the recurrence order (more precisely: the holonomic rank) may increaseduring each such application (but not more than the sum of their orders).On the other hand, if we choose to deal with the full triple sum rightfrom the start, then we can do better. By observing that when k − m < , we have the situation where r − < i ≤ r − − ( k − m ) forces the innermostbinomial coefficient ( r − i ) to be 0. Thus, the max function in G s ( x ) can besafely removed. Addition is commutative, so it does not matter if we move allsummations to the front and consider only one summand with three indexedparameters. In other words, G s ( x ) can be rewritten as m + s − ∑ k = s ∑ r = r − −( k − m ) ∑ i = ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ quantifiers grouped ( sr )( k − r − )( r − i ) b − (− b ) r − i ( bx ) k ·„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„‚„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„„¶ one summand . (3.1) Another “preparation” step involves employing the Guess package [1] to pre-dict the recurrence that our polynomial satisfies, by using sufficiently genericevaluations of (3.1). This step gives us confidence that a nice enough recur-rence exists, and it serves as an additional sanity check. For our problem,the (computer) guessing procedure already produced the claimed minimalthird-order recurrence (from Theorem 3.1) in the parameter s . To prove thatthe guess is correct, it is enough to compute the same recurrence (or a higherorder one) via creative telescoping. In the latter case, one also has to verifythat the guessed recurrence (operator) is a right factor of the bigger one, andconsider a sufficient number of initial values.In general, the creative telescoping algorithm works very well and canbe applied directly without adjustments if all conditions are “ideal”. In suchcases, the outputted telescoper corresponds exactly to the recurrence that wewant. If certain issues arise such that (human) adjusting for them producesextraneous terms that are not a part of the original sum and/or come fromcompensations due to the rebuilding of the original sum, then we need tocollect all such terms and find a collective annihilator for them for the purposeof “homogenizing” the recurrence given by the telescoper. We found that some certificates Q contain singularities on the boundaryvalues of the inner sum. This implies that the limits of the sum must beadjusted so that we avoid evaluating at those points. In fact, the summationrange must be adjusted so that there are no singularities at all intermediatevalues.To illustrate this a little more generally, suppose that, upon applyingthe creative telescoping algorithm on the summation with respect to r , thecomputer outputs a certificate Q ∈ Q ( s, r ) containing poles at r i ∈ [ , s + ] ⊂ N for a finite number of i . All parameters besides r are treated symbolically.reative Telescoping on Multiple Sums 9Then we can see that the sum ∑ sr = ( S r − ) ⋅ Q ⋅ F ( s, r ) cannot be determinedbecause evaluations at those poles for the summation range [ , s ] are notpossible and therefore the sum is undefined.Such singularities can be removed from the offending sum so that theevaluation(s) can happen. We also remove the exact same values from thesummations containing the telescopers to match so that (2.4) makes sense.The summations with the telescopers can then be subsequently “filled in”with the removed summands (balanced of course by subtracting those sameterms from the inhomogeneous part). This strategy can be quite effective ifall of the poles are collected contiguously at either of the summation bounds,or if there are only one or two. Note that if the number of compensated termsrequired is too large, it may be be better to consider another strategy, suchas rewriting terms in some alternative (but equivalent) form to avoid the poleentirely (as was suggested in the analysis for (2.5)). We also found that the telescoper P does not commute with our summation.To illustrate this a little more generally, suppose the operator P is in the Orealgebra generated by the shift operator S s . In other words, suppose that P can be written as some polynomial in S s , for example, P = p + p S s + ⋯ + p j S js , where the p i ( i = , . . . , j ) may be rational functions containing the parame-ter s . If we apply such a P to a summation of the form ∑ m + s − k = H ( s, k ) , thenwe face the issue that the application not only affects the parameter s in thesummand H ( s, k ) , but also the upper limit m + s −
1. Then if we apply P tothe whole sum, we get p m + s − ∑ k = H ( s, k ) + p m + s ∑ k = H ( s + , k ) + ⋯ + p jm + s + j − ∑ k = H ( s + j, k ) . (3.2)It is quite obvious that this is not the same as applying P to only the sum-mand H ( s, k ) : m + s − ∑ k = ( p H ( s, k ) + p H ( s + , k ) + ⋯ + p j H ( s + j, k )) . (3.3)However, we can simulate the “factoring out” of the P in (3.3) if we peeloff a sufficient (and finite) number of terms from each sum in (3.2) such thatits upper limits are all m + s −
1. Then (3.3) can be replaced by the peeledversion of (3.2) with P on the outside and the removed terms can then bemerged with the inhomogeneous part. In all of our examples, the adjustment of the summation limits to avoidsingularities in the certificates was completed first. After that, it is a priorinot clear if one should proceed to adjust for the operator commuting withthe summand or to fill in the terms that were removed for the singularities.0 C. Koutschan and E. WongThe (human) decision may depend on the number of singularities, where thesingularities are located, how complicated the telescoper expression is, andwhether or not the lower/upper boundaries are influenced by the telescopers.This makes it difficult to automate adjustments effectively.Once we have collected all of the inhomogeneous parts, we face thequestion of how to process them. In principle, we could just write them downand try to use Mathematica’s symbolic power to simplify them as much aspossible. Unfortunately, this does not work very well on our problem, withthe only progress being that some of the inhomogeneous parts convenientlycollapse to zero (so we remove them). Instead, we take advantage of the factthat we can write all of the inhomogeneous parts as different shifts and sub-stitutions of the given summand. More precisely, the total of some of theseparts can be expressed as an operator applied to the summand, followed bya substitution. Then, an annihilator for this total can be derived by apply-ing the closure properties “application of an operator” and “integer-linearsubstitution”. In this way, we completely avoid dealing with expressions likeMathematica’s DifferenceRoots.We illustrate this strategy with an example that comes from the inho-mogeneous parts for G ( ) . It involves observing patterns in a complicatedexpression to construct an operator that would give the same result when itis applied to some simpler version of the expression. First, we identified allterms that contain a certain hypergeometric F series: − ( b − )( m + s + )( bx ) m + s + b x ⋅ F ( − s, − m − s ; 2; b − b )+ ( b − )( m + bs )( bx ) m + s b ⋅ F ( − s, − m − s ; 2; b − b )( b − )( s + )( bx − )( bx ) m + s b ⋅ F ( − m − s, − s ; 2; b − b ) . (3.4)From there, we can see that selecting the operator ( b − )( m + s + ) b x ( bx ) S m − ( b − )( m + bs ) b ( bx ) S m + ( b − )( s + )( bx − ) b ( bx ) S s (3.5)and applying it to ( bx ) m + s ⋅ F ( − s, − m − s ; 2; b − b ) gives (3.4). The anni-hilator for (3.4) can therefore be obtained by “applying” (3.5) (in the senseof closure properties) to the annihilator of this single expression. This pro-cess is much faster than trying to directly compute an annihilating ideal ofthe sum of hypergeometric series, with the added benefit that the order ofthe recurrence will usually be smaller compared to applying the Annihilatorcommand directly to expressions such as (3.4). In our particular example, thelatter method even caused the program to crash.We can therefore see that directly constructing an operator by closelyinspecting patterns in the inhomogeneous parts can be computationally effec-tive. In a way, the fact that we use F ’s in the previous argument is inconse-quential to the construction of the operator to be applied (it could have easilybe replaced with a symbolic expression that exhibits similar shift behaviors,reative Telescoping on Multiple Sums 11for example). Thus, an annihilating operator for the inhomogeneous termscan be deduced in this way before administering any substitutions. As men-tioned before, this could also involve removing all terms that would collapseto zero anyway and building from scratch the new operator by only using theshifts needed to produce compensation terms that may have resulted fromthe treatment of singularities and commutation. We similarly applied thisstrategy to G ( ) s and it cost us 30 hours. The collection of all inhomogeneous parts and its subsequent removal via itsannihilator can eat up a lot up computation time depending on how compli-cated these parts are. In particular, we experience this in the computation ofthe annihilators for the inhomogeneous part of G ( ) s when applying the clo-sure property of integer-linear substitution. Thus, as an alternative to blindlyapplying the corresponding computer command and not knowing what is go-ing on behind the scenes while waiting patiently for the code to finish, we cantake better control of the process by making a few additional optimizations,resulting in a significant speedup in computation time.With the “application of an operator” closure property we produced anannihilating ideal, with its Gr¨obner basis denoted by U ( ) , for the combinedinhomogeneous parts of G ( ) s , denoted by H ( s, k ) , but without the necessarysubstitution k → m + s according to the upper summation bound. This impliesthat it is necessary to apply the closure property “integer-linear substitution”to U ( ) .The above Gr¨obner basis U ( ) has the set of irreducible monomials { S s , S k , } and hence is of holonomic rank 3. The theory tells us that for H ( s, m + s ) we should expect a recurrence of order 3 in s , however, we know(from the longer computation in Section 3.5) that our final recurrence is oforder 2 in s .For constructing such a recurrence, we want to find an operator T in theleft ideal generated by U ( ) with the support { S s S k , S s S k , } , correspondingto a bivariate recurrence involving the terms H ( s + , k + ) , H ( s + , k + ) , H ( s, k ) , which, after the substitution k → m + s , turns into the desired second-orderrecurrence for H ( s, m + s ) . We therefore make the following ansatz for T : T = c ( k, s ) S k S s + c ( k, s ) S k S s + c ( k, s ) , where the coefficients c ( k, s ) , c ( k, s ) , c ( k, s ) are to be determined. Gr¨obnerbasis theory tells us that this T is an element of the annihilating ideal (inother words: represents a valid recurrence for H ( s, k ) ) if and only if it reducesto 0 by the Gr¨obner basis U ( ) . Reducing the ansatz T by U ( ) results in alinear combination of the basis monomials { S s , S k , } , i.e., an Ore polynomialof the form E ( k, s, c , c , c ) S s + E ( k, s, c , c , c ) S k + E ( k, s, c , c , c ) E , E , E . This polynomial is zero if and only if E = E = E =
0, so we proceed to solve this system for c , c , c .Ultimately, this procedure gives us an Ore polynomial with the supportwe want that annihilates the inhomogeneous parts after substituting k → m + s in the coefficients of T and omitting the shift operator S k . This is essentiallythe same as saying: substitute k → m + s into the recurrence T ⋅ H ( s, k ) .However, since this substitution tends to decrease the size of expressions (itreduces the number of variables), it is desirable to perform it as early aspossible, and not only at the very end.Indeed, we were able to speed up our computation significantly by per-forming the substitution already during the reduction of the monomials of T ,but care has to be taken: to match the leading monomials, we may have tomultiply by (a power of) S k , and for this (noncommutative) multiplicationone needs to keep the variable k . However, it can be substituted immediatelyafterwards. This leads to a less dramatic swell of expressions. Other sourcesof speedup include a manual selection strategy of Gr¨obner basis elements tobe used for the reduction, and the order in which these reductions are made.It might be worthwhile to note that there are places in this processwhere we got lucky: as fate would have it, the coefficient E of S s is zero (butonly for k = m + s ), and this gives us the luck of finding a recurrence of ordertwo instead of three! This procedure has now reduced our total computationtime to 1.4 hours. This next strategy takes a more holistic approach in that we treat the triplesum (3.1) all at once and make adjustments to the single summand in orderto deal with the “unnatural boundary” problems simultaneously. For ourproblem, this is accomplished by making the observation that the issue of“unnatural boundaries” occurs whenever all three binomial coefficients inour summand are nonzero beyond the limits of our summation. Assumingthat m and s are fixed positive integers, we let B denote the collection ofpoints ( k, i, r ) ∈ Z for which the summand in (3.1) is nonzero, and let B denote all points ( k, i, r ) ∈ Z that are inside the summation ranges. While B corresponds to all integer points of a bounded polytope in R , the set B isunbounded (see Fig. 1). We essentially want to sum over all points ( k, i, r ) inthe intersection of B and B (depicted in blue), while we want to avoid thosepoints in B ∖ B (depicted in red). Hence, the set of “bad points” also formsan infinite polytope and we remove these points using gamma functions.We first recall that Γ ( k ) has poles exactly at the non-positive integers,and therefore ( k ) has zeros at k = , − , − , . . . . Then the summand canbe modified by the following gamma functions in order to enforce naturalboundaries: C ( ε, i, k, r ) ∶= ( sr ) ⋅ ( k − r − ) ⋅ ( r − i ) ⋅ Γ ( k + ε ) Γ ( k ) ⋅ Γ ( r − i − ( k − m ) + ε ) Γ ( r − i − ( k − m )) reative Telescoping on Multiple Sums 13 Figure 1. “Bad” points B ∖ B (red) and “good” points B ∩ B (blue) for fixed m =
15 and s = ε . Upon sending ε → C ( ε, i, k, r ) to be zero whenever k ≤ r − i − ( k − m ) ≤
0. In other words,the introduction of the reciprocal gammas is balanced out by the perturbedgammas in the numerators, which conveniently avoids division by zero andgives an equivalent (final) result for the original problem after setting ε = m + s − ∑ k = s ∑ r = r − −( k − m ) ∑ i = C ( ε, i, k, r ) ⋅ b − (− b ) r − i ( bx ) k and afterwards take the limit ε → G s ( x ) . Unfortunately, we do not get the minimal-order recurrence, but a4 C. Koutschan and E. Wongfourth-order one. This has allowed us to reduce our computation time toabout 11 minutes.However, this is not yet the end of the story. In Fig. 1 we were trappedby Mathematica’s definition of the binomial coefficient (cf. the discussion inSection 2): actually, the summand is zero for k ≤ ( nk ) is zero unless 0 ≤ k ≤ n . Theconditions for the summand to be nonzero (implied by the three binomialcoefficients) somehow correspond to the summation bounds (given by thethree summation quantifiers), which is illustrated in the following table (it isactually a curiosity of this problem that we can find such correspondence):Factor inSummand Nonzero Range Summation Bounds( B ) ( B ) ( sr ) ≤ r ≤ s ≤ r ≤ s ( k − r − ) ≤ r − ≤ k − ≤ k ≤ m + s − ( r − i ) ≤ i ≤ r − ≤ i ≤ r − − ( k − m ) After close inspection of this table, it becomes evident that only onegamma correction is actually needed (and hence Fig. 1 does not show thetrue situation). We can therefore redefine C ( ε, i, k, r ) ∶= ( sr ) ⋅ ( k − r − ) ⋅ ( r − i ) ⋅ Γ ( r − i − ( k − m ) + ε ) Γ ( r − i − ( k − m )) . This observation speeds up the computations significantly, and the winningtime is 30 seconds. Moreover, we obtain the minimal recurrence of orderthree. The reader may now wonder how we can tell the HolonomicFunctionspackage that this computation should be executed with a different definitionof the binomial coefficient (that differs from Mathematica’s)? The answer is:we do not have to, since it is completely irrelevant (from the viewpoint of thepackage), because both versions of the binomial coefficient satisfy the verysame recurrence equations, as we have seen in Section 2! The difference onlybecomes relevant when we evaluate the summand at particular values (whichis done outside of the package), e.g., when checking initial conditions.
4. Conclusions and Future Work
In this expository article we demonstrated the usage of the HolonomicFunc-tions package to deal with an intricate triple sum coming from an applicationin quasi-Monte Carlo integration. We had two goals in mind for this paper:first, we felt the need to deliver some technical details for a key lemma in [6],where only the main ideas of the computer algebra proof were mentioned;second, we wanted to provide a somewhat easy-to-digest tutorial for provingspecial function and combinatorial identities with the help of the computerby expounding on the difficulties that may arise in similar applications andreative Telescoping on Multiple Sums 15Approximate Strategies ResultComp. Time Implemented30 hours – split sums fifth-order recurrence inthe ideal generated by theguessed recurrence– sing./comm. corrections– closure properties1.4 hours – split sums– sing./comm. corrections– closure properties– substitution speedup11 minutes – one triple sum fourth-order recurrence inthe ideal generated by theguessed recurrence– two gamma insertions– sing. corrections30 seconds – one triple sum same third-orderrecurrence as the guessedrecurrence– one gamma insertion– sing. corrections
Figure 2.
Results and Comparisonsa highlighting a few creative ways to cure them. We hope that we haveconvinced the reader that it is not always so cut and dry to prove a givenidentity with the holonomic systems approach. While in principle it allowsone to prove holonomic identities in an automated way, we have seen that inpractice, even with using current state-of-the-art software tools, many stepsin the proof require human interaction. At many positions in the provingprocess, we had to make a choice on how to proceed, and the decision mayboth influence the optimality of the final result and the time that is requiredto obtain it. Fig. 2 gives an impressive overview how much difference in run-time such choices can make. In the future we plan to automate some of theproof steps that had to be done “by hand” in this case study, for example,in the analysis of singularities in the certificate(s) and dealing with the issueof commutation.
Acknowledgements
We would like to thank the organizers of CASC 2020 for providing an occasionand opportunity to give a talk about the work on which this article is based.We were encouraged by the positive feedback from the audience, which mo-tivated this post-proceedings contribution. We also express our appreciationto Hao Du and Ali Uncu for their support and helpful commentary.