Internal versus external balancing in the evaluation of graph-based number types

TL;DR

This paper compares internal and external balancing techniques for graph-based number types, demonstrating how restructuring and balancing can improve evaluation efficiency through theoretical analysis and experiments.

Contribution

It provides a comparative analysis of internal versus external balancing methods for graph-based number types, highlighting their respective advantages and disadvantages.

Findings

External balancing can significantly improve evaluation performance.

Internal balancing offers a different approach to error management.

Theoretical and experimental results illustrate the trade-offs between the methods.

Abstract

Number types for exact computation are usually based on directed acyclic graphs. A poor graph structure can impair the efficency of their evaluation. In such cases the performance of a number type can be drastically improved by restructuring the graph or by internally balancing error bounds with respect to the graph's structure. We compare advantages and disadvantages of these two concepts both theoretically and experimentally.

Figures

Tables

Table 1. Table 1 : Operation-dependent constants c ⁡ ( e l ) c subscript 𝑒 𝑙 \operatorname{c}({e_{l}}) and c ⁡ ( e r ) c subscript 𝑒 𝑟 \operatorname{c}({e_{r}}) for an accuracy-driven evaluation in Real_algebraic , with x h ​ i ​ g ​ h , y h ​ i ​ g ​ h subscript 𝑥 ℎ 𝑖 𝑔 ℎ subscript 𝑦 ℎ 𝑖 𝑔 ℎ {x_{high}},{y_{high}} upper bounds and x l ​ o ​ w , y l ​ o ​ w subscript 𝑥 𝑙 𝑜 𝑤 subscript 𝑦 𝑙 𝑜 𝑤 {x_{low}},{y_{low}} lower bounds on the child values.

	negation	add./sub.	multipl.	division	$d$-th root
$c\left(e_l\right)$	$1$	$1$	$y_{high}$	$\frac{1}{y_{low}}$	$\frac{1}{d}{\left(x_{low}\right)}^{\frac{1-d}{d}}$
$c\left(e_r\right)$	$0$	$1$	$x_{high}$	$\frac{1}{y_{low}^{}{}_{}{}^2}$	$0$

Equations

$$2^{q+i_v}+c_l{2}^{q+i_l}+c_r{2}^{q+i_r}\le 2^q\text{ or, equivalently, }{2}^{i_v}+c_l{2}^{i_l}+c_r{2}^{i_r}\le 1$$

$$\cos {\mathrm{t}}_{\mathrm{i}}(v)=-(q+i_v+w_l{i}_l+w_r{i}_r)$$

$$(1+w_l)z_l+w_l{z}_r-w_l$$

$$(1+w_r)z_r+w_r{z}_l-w_r$$

$$i_l$$

$$i_r$$

$$i_v$$

$$\cos {\mathrm{t}}_{\mathrm{v}}(P)=-{\sum }_{e\in P}(i(e)+\log (c(e)))-i(v)\text{ and }\cos {\mathrm{t}}_{\mathrm{f}}(P)={\sum }_{e\in P}\log (c(e))$$

$$\mathrm{cost}(E)=n\log (n)+{\sum }_{v\in E}\log \left({\sum }_{P\in P(v)}{2}^{\cos {\mathrm{t}}_{\mathrm{f}}(P)}\right)-nq$$

$$w(e)=\frac{{\mathrm{c}}_{\mathrm{f}}(v,e)}{{\mathrm{c}}_{\mathrm{f}}(v)}w(v)=\frac{{\sum }_{P_e\in P(v),e\in P_e}{2}^{\cos {\mathrm{t}}_{\mathrm{f}}(P_e)}}{{\sum }_{P\in P(v)}{2}^{\cos {\mathrm{t}}_{\mathrm{f}}(P)}}w(v)$$

$$\cos {\mathrm{t}}_{\mathrm{r}}(P)=\cos {\mathrm{t}}_{\mathrm{v}}(P)+\cos {\mathrm{t}}_{\mathrm{f}}(P)$$

$$=-{\sum }_{e\in P}(i(e)+\log (c(e)))-i(v)+\cos {\mathrm{t}}_{\mathrm{f}}(P)$$

$$=-{\sum }_{j=0}^k(\log (w(e_j))-\log (w(v_j))+\log (w(v))+\cos {\mathrm{t}}_{\mathrm{f}}(P)$$

$$=\log (w(v_0))-{\sum }_{j=0}^k\left(\log ({\mathrm{c}}_{\mathrm{f}}(v_j))+\log (c(e_j))-\log ({\mathrm{c}}_{\mathrm{f}}(v_{j+1})\right)+\cos {\mathrm{t}}_{\mathrm{f}}(P)$$

$$=\log (w(v_0))-\log ({\mathrm{c}}_{\mathrm{f}}(v))$$

$$\delta (v)\le -\delta (e_l)w(e_l)-\delta (e_r)w(e_r)=-{\sum }_{v^{\prime }\in E}\left(\delta (e_l)\frac{{\mathrm{c}}_{\mathrm{f}}(v^{\prime },e_l)}{{\mathrm{c}}_{\mathrm{f}}(v^{\prime })}+\delta (e_r)\frac{{\mathrm{c}}_{\mathrm{f}}(v^{\prime },e_r)}{{\mathrm{c}}_{\mathrm{f}}(v^{\prime })}\right)$$

$$\Delta \mathrm{cost}(E)$$

$$=-\sum _{v\in E}\underset{P\in P(v)}{\min }\sum _{e\in P}\delta (e)-\sum _{v\in E}\delta (v)$$

$$\ge -\sum _{v\in E}\underset{P\in P(v)}{\min }\sum _{e\in P}\delta (e)+\sum _{v\in E}\sum _{e\in \text{E}(E)}\delta (e)\frac{{\mathrm{c}}_{\mathrm{f}}(v,e)}{{\mathrm{c}}_{\mathrm{f}}(v)}$$

$$=-\sum _{v\in E}\underset{P\in P(v)}{\min }\sum _{e\in P}\delta (e)+\sum _{v\in E}\sum _{P\in P(v)}\sum _{e\in P}\delta (e)\frac{2^{\cos {\mathrm{t}}_{\mathrm{f}}(P)}}{{\mathrm{c}}_{\mathrm{f}}(v)}$$

$$\ge -\sum _{v\in E}\underset{P\in P(v)}{\min }\sum _{e\in P}\delta (e)+\sum _{v\in E}\underset{P\in P(v)}{\min }\sum _{e\in P}\delta (e)=0$$

$$\mathrm{cost}(E)=\sum _{v\in E}\underset{P\in P(v)}{\max }\cos {\mathrm{t}}_{\mathrm{r}}(P)-nq=\sum _{v\in E}(\log (n)-\log ({\mathrm{c}}_{\mathrm{f}}(v)))-nq$$

$$\mathrm{cost}(T)=n\log n+\sum _{v\in V}\cos {\mathrm{t}}_{\mathrm{f}}(\mathrm{path}(v))-nq$$

$$\log (1+r)\le \log (1)+r\frac{d}{dx}\log (x){|}_1=\frac{r}{\ln (2)}$$

$$\mathrm{cost}(P)=k\log k+\cos {\mathrm{t}}_{\mathrm{C}}(P)$$

$$\mathrm{cost}(P)$$

$$=k\log k+\frac{k(k-1)}{2}+\cos {\mathrm{t}}_{\mathrm{C}}(P)$$

$$\mathrm{cost}(\mathrm{cp}(L))-d_l{i}_l=\mathrm{cost}(\mathrm{cp}(R))-d_r{i}_r$$

$$\mathrm{cost}(\mathrm{cp}(X))=\mathrm{cost}(\mathrm{cp}(L))-d_l{i}_l-i_v$$

$$\frac{-z-c{d}_f{z}^{d_f}}{1-z-c{z}^{d_f}}-d_l=0⟺c\frac{d_r-1}{d_r}{z}^{d_f}+\frac{d_l-1}{d_l}z-1=0$$

$$i_v$$

Full text

1 Introduction

Inexact computation causes many problems when algorithms are implemented, ranging from slightly wrong results to crashes or invalid program states. This is especially prevalent in the field of computational geometry, where real number computations and combinatorical properties intertwine [10]. In consequence, various exact number types have been developed [6, 8, 16]. It is an ongoing challenge to make these number types sufficiently efficient to be an acceptable alternative to floating-point primitives in practical applications. Number types based on the Exact Computation Paradigm recompute the value of complex expressions if the currently stored error bound is not sufficient for an exact decision [15]. Hence, they store the computation history of a value in a directed acyclic graph, which we call an expression dag. The structure of the stored graph is then determined by the order in which the program executes the operations. It lies in the nature of iterative programming that values are often computed step by step, resulting in list-like graph structures.

Re-evaluating expressions in an unbalanced graph is more expensive than in a balanced one [4, 11]. We discuss two general approaches on reducing the impact of graph structure on the evaluation time. Prior to the evaluation, the expression dag can be restructured. Originally proposed by Yap [15], restructuring methods with varying degrees of invasiveness were developed [11, 14]. Root-free expression trees can be restructured to reach optimal depth as shown by Brent [3]. In Section 2.1 we introduce a weighted version of Brent’s algorithm applied on maximal subtrees inside an expression dag. Besides restructuring, which can be considered ‘external’ with respect to the evaluation process, we can make ‘internal’ adjustments during the evaluation to compensate for bad structure. Error bounds occuring during an evaluation can be balanced to better reflect the structure of the graph [4]. Doing so requires a switch from an integer to a floating-point error bound representation, leading to numerical issues that need to be taken into consideration [7, 13]. In Section 2.2 we show how error bounds can be balanced optimally in both the serial and the parallel case and compare several heuristics. Finally, in Section 3 we experimentally highlight strengths and weaknesses of each approach.

2 Concepts

An expression dag is a rooted ordered directed acyclic graph in which each node is either a floating-point number, a unary operation ($\sqrt[d]{\phantom{x}}$,$-$) with one child, or a binary operation ($+$,$-$,$*$,$/$) with two children. We call an expression dag $E^{\prime}$, whose root is part of another expression dag $E$ a subexpression of $E$. We write $v\in E$ to indicate that $v$ is a node in $E$ and we write $|E|$ to represent the number of operator nodes in $E$. In an accuracy-driven evaluation the goal is to evaluate the root node of an expression dag with absolute accuracy $q$, i.e., to compute an approximation ${\tilde{x}}$ for the value $x$ of the represented expression, such that $|{\tilde{x}}-x|\leq 2^{q}$ (cf. [15]). To reach this goal, sufficiently small error bounds for the (up to two) child nodes and for the operation error are set and matching approximations are computed recursively for the children. Let $v\in E$ be a node with outgoing edges ${e_{l}}$ to the left and ${e_{r}}$ to the right child. Let $\operatorname{i}({e_{l}}),\operatorname{i}({e_{r}})$ be the increase in accuracy for the left and the right child of $v$ and $\operatorname{i}(v)$ be the increase in accuracy for the operation (i.e. the increase in precision) at $v$. Depending on the operation in $v$ we assign constants ${\operatorname{c}({e_{l}})},{\operatorname{c}({e_{r}})}$ to its outgoing edges as depicted in Table 1. If the node $v$ is known from the context, we shortly write ${i_{v}},{i_{l}},{i_{r}}$ for the accuracy increases at $v,{e_{l}},{e_{r}}$ and ${c_{l}},{c_{r}}$ for the respective constants. To guarantee an accuracy of $q$ at $v$, the choice of ${i_{v}},{i_{l}},{i_{r}}$ must satisfy the inequality

[TABLE]

Aside from this condition, the choice of ${i_{v}},{i_{l}},{i_{r}}$ is arbitrary and usually done by a symmetric distribution of the error. In the exact number type Real_algebraic they are chosen such that ${c_{l}}2^{i_{l}}\leq 0.25$, ${c_{r}}2^{i_{r}}\leq 0.25$ and $2^{i_{v}}=0.5$ (and adjusted accordingly for one or zero children). Let the depth of a node $v$ in an expression dag be the length of the longest path from the root to $v$. In general, the precision ${p_{v}}$ needed to evaluate a node $v$ increases linearly with the depth of the node due to the steady increase through ${i_{l}},{i_{r}}$. The approximated value of each node is stored in a multiple-precision floating-point type (bigfloat). The cost of evaluating a node is largely dominated by the cost of the bigfloat operation, which is linear in $|{p_{v}}|$ in case of addition and subtraction and linear up to a logarithmic factor in case of multiplication, division and roots. So the precision ${p_{v}}$ is a good indicator for the total evaluation cost of a node (except for negations).

Let $E$ be an expression dag. We define the cost of a node $v\in E$ to be $|{p_{v}}|$ and the cost of $E$, denoted by $\operatorname{cost}(E)$, as the sum of the cost of all nodes in $E$. We set the depth of $E$ to the maximum depth of all nodes in $E$. Let ${E_{list}}$ be a list-like expression dag, i.e., an expression dag with depth $\Theta(n)$, where $n$ is the number of its nodes and let ${E_{bal}}$ be a balanced expression dag, i.e., an expression dag with depth $\Theta(\log(n))$. Since the precision increases linearly with the depth, we have $\operatorname{cost}({E_{list}})=\Theta(n^{2})$ and $\operatorname{cost}({E_{bal}})=\Theta(n\log(n))$, assuming that the operation constants can be bounded (cf. [11]). In a parallel environment the cost of the evaluation is driven by dependencies between the nodes. For an expression dag $E$ with $n$ nodes let the cost of a path in $E$ be the sum of the cost of the nodes along the path. Let $\operatorname{cp}(E)$ be a path in $E$ with the highest cost. We call $\operatorname{cp}(E)$ a critical path in $E$. Then the cost of evaluating $E$ in parallel is $\Theta(\operatorname{cost}(\operatorname{cp}(E)))$ with $O(n)$ processors. Let ${E_{list}},{E_{bal}}$ be defined as before. Then obviously $\operatorname{cost}(\operatorname{cp}({E_{list}}))=\Theta(n)$ and $\operatorname{cost}(\operatorname{cp}({E_{bal}}))=\Theta(\log n)$ (cf. [14]). So in both the serial and the parallel case, balanced graph structures are superior.

2.1 Graph Restructuring

By definition, exact number types that use accuracy-driven evaluation act lazy, i.e., expressions are not evaluated until a decision needs to be made. Before their first evaluation, underlying graph structures are lightweight and can be changed at low cost. Therefore graph restructuring algorithms ideally take place when the first decision is demanded. While it is not impossible to restructure graphs that have already been evaluated, it comes with several downsides. Since subexpressions will change during restructuring, all approximations and error bounds associated with these subexpressions are lost, although they could be reused in later evaluations. Since stored data may depend on data in subexpressions, the internal state of the whole expression dag may be invalidated. Those effects can make restructuring expensive if many decisions are requested without significant changes to the graph in between. Let $E$ be an expression dag. We call a connected, rooted subgraph of $E$ an operator tree if it consists solely of operator nodes, does not contain root operations and does not contain nodes with two or more parents (not necessarily in $E$), except for its root. We restructure each maximal operator tree in $E$ according to a weighted version of Brent’s algorithm. Let $T$ be an operator tree in $E$. We call the children of the leaves of $T$ the operands of $T$ and associate a positive weight with each of those operands. We define a weight function, such that for each node $v\in T$ the weight of $v$ is greater or equal than the weight of its children. The main difference between the original algorithm and the weighted variation lies in the choice of the split node. We give a brief outline of the algorithm.

The algorithm builds upon two operations, compress and raise. The operation compress takes an expression tree $E$ and returns an expression tree of the form $F/G$ and raise takes an expression tree $E$ and a subtree $X$ and returns an expression tree of the form $(AX+B)/(CX+D)$, where $A,B,C,D,F,G$ are division-free expression trees with logarithmic depth.

Let ${v_{r}}$ be the root node of $T$. We choose ${v_{s}}$ as a node with maximal weight in $T$ such that both children have either weight $<\frac{1}{2}\operatorname{weight}({v_{r}})$ or are operands. Note that this implies $\operatorname{weight}({v_{s}})\geq\frac{1}{2}\operatorname{weight}({v_{r}})$. We then recursively call compress on ${v_{s}}$ and raise ${v_{s}}$ to the root by repeating the following steps:

Search for a new split node ${v_{s}^{\prime}}$ on the path from ${v_{r}}$ to ${v_{s}}$ that splits at a weight of $\frac{1}{2}(\operatorname{weight}({v_{r}})+\operatorname{weight}({v_{s}}))$. 2. 2.

Recursively raise ${v_{s}^{\prime}}$ to ${v_{r}}$ and ${v_{s}}$ to the respective child node in ${v_{s}^{\prime}}$. 3. 3.

Substitute ${v_{s}^{\prime}}$ and its children into ${v_{r}}$ by incorporating the operation at ${v_{s}^{\prime}}$.

Let $R$ be the expression at ${v_{r}}$, let $Y$ be the expression at ${v_{s}^{\prime}}$ and let $X$ be the expression at ${v_{s}}$. After the second step, $R=\frac{A^{\prime}Y+B^{\prime}}{C^{\prime}Y+D^{\prime}}$ and we have $Y=Y_{L}\circ Y_{R}$ with $Y_{L}=\frac{A^{\prime\prime}X+B^{\prime\prime}}{C^{\prime\prime}X+D^{\prime\prime}}$ and $Y_{R}=F^{\prime\prime}/G^{\prime\prime}$ or vice versa. Substituting $Y$ (with respect to the operation $\circ$ at $Y$) then gives the desired $R=\frac{AX+B}{CX+D}$. Substituting $X=F^{\prime}/G^{\prime}$ finally leads to a balanced expression of the form $R=F/G$.

The new split operation is shown in Algorithm 2. If unit weight is chosen, there will never be an operand that does not satisfy the split condition. If furthermore the weight function is chosen as the number of operands in a subtree, satisfying the split condition implies having a bigger weight than the sibling. Therefore the algorithm is identical to Brent’s original algorithm applied to subtrees of the expression dag and guarantees logarithmic depth for the new operator tree. Regarding the overall expression dag, nodes which contain root operations, have more than one parent or have been evaluated before are treated equally to the other operands in this case and therefore act as ‘blocking nodes’ for the balancing process. Let $k$ be the number of these blocking nodes in $E$. If the number of incoming edges for each blocking node is bounded by a constant, the depth of $E$ after applying the algorithm to each operator tree is in $O(k\log(\frac{n}{k}))$. This depth can be reduced by applying appropriate weights to the blocking nodes. From a conceptual perspective, a sensible choice for the weight of an operand (as well as for the weights of the inner nodes) would be the number of operator nodes in the subexpression rooted at the operand. Note that we are actually interested in the number of bigfloat operations. However, it is very expensive to compute the number of descendants for a node in a DAG, since one has to deal with duplicates [2]. Ignoring duplicates, we could choose the number of operators we would get by expanding the DAG to a tree. While computable in linear time, the number of operators can get exponential (cf. [4, 11]) and therefore we cannot store the exact weight in an integer data type anymore. There are ways of managing such weights, as we discuss in Section 2.2, but they are imprecise and less efficient than relying on primitives. Both weight functions behave identical to the unit weight case when there are no blocking nodes present. If there are blocking nodes on the other hand, these nodes get weighted accordingly and expensive nodes are risen to the top of the operator tree. The depth after restructuring for $k$ blocking nodes therefore becomes $O(k+\log n)$.

The weight functions described above are optimal, but hard to compute. Let the weight of both operators and operands be the depth of the subexpression rooted at the operand or operator in the underlying expression dag. Then the algorithm subsequently reduces the length of the longest paths in the expression dag. Note that this strategy does not necessarily lead to an optimal result. Nevertheless, computing the depth of a subexpression in an expression dag can be done fast and the depth can be represented efficiently. Therefore this strategy might prove to be a good heuristic to combine advantages of the unit weight algorithm and the weighted approach.

2.2 Error Bound Balancing

As described at the start of this section, the additional cost of unbalanced graph structures originates in the increase in accuracy associated with each node. A more careful choice of ${i_{v}},{i_{l}},{i_{r}}$ in (1) may compensate for an unfavorable structure. If set correctly, linear depth still only leads to a logarithmic increase in accuracy aside from operation constants [4].

An increase in accuracy at an operator node only affects the operation itself. An increase in accuracy for a child node affects all operations in the subexpression of the child. We associate a non-negative weight $\operatorname{w}(e)$ with each edge $e$ in an expression dag $E$, representing the impact a change in $\operatorname{i}(e)$ has on the total cost of $E$. For a node $v\in E$ with outgoing edges ${e_{l}},{e_{r}}$ let ${w_{l}}=\operatorname{w}({e_{l}})$ and ${w_{r}}=\operatorname{w}({e_{r}})$. We then say that, for an evaluation to accuracy $q$, the cost induced on $E$ by the choice of parameters in $v$ is given by

[TABLE]

whereas $\operatorname{cost}(E)=\sum_{v\in E}\operatorname{cost_{i}}(v)$. To minimize the total cost we want to minimize the cost induced by each node while maintaining the condition in (1). Let ${z_{l}}={c_{l}}2^{i_{l}}$, ${z_{r}}={c_{r}}2^{i_{r}}$ and let ${w_{all}}=1+{w_{l}}+{w_{r}}$. With an optimal choice of the parameters, (1) is an equality and we have ${i_{v}}=\log(1-{z_{l}}-{z_{r}})$. Substituting ${i_{v}}$ into (2) and setting $\frac{\partial}{\partial{i_{l}}}\operatorname{cost_{i}}(v)=\frac{\partial}{\partial{i_{r}}}\operatorname{cost_{i}}(v)=0$ we get

[TABLE]

leading to ${z_{l}}=\frac{{w_{l}}}{{w_{all}}}$ and ${z_{r}}=\frac{{w_{r}}}{{w_{all}}}$. Resubstituting ${z_{l}}$ and ${z_{r}}$, the optimal choice of the parameters for error bound distribution inside a node is

[TABLE]

We can show that this parameter choice makes the cost of the evaluation to some degree independent of the structure of the graph. For a node $v\in E$ we denote the set of paths between the root node of $E$ and $v$ by $\operatorname{\textbf{P}}(v)$. For a path $P\in\operatorname{\textbf{P}}(v)$ we write $e\in P$ to indicate that $e$ is an edge along $P$. The precision requested at $v$ along $P$ can be expressed as $\operatorname{cost_{r}}(P)=\operatorname{cost_{v}}(P)+\operatorname{cost_{f}}(P)$ where

[TABLE]

denote the variable cost induced by the choice of ${i_{l}},{i_{r}},{i_{v}}$ and the fixed cost induced by the operation constants along the path.

Theorem 2.1

Let $E$ be an expression dag consisting of $n$ unevaluated operator nodes. Then the cost of evaluating $E$ with accuracy $q\leq 0$ and with an optimal choice of parameters is

[TABLE]

Proof

We define weights for each node $v$ and each edge $e$ in $E$ with respect to (2). Let $\textstyle\operatorname{c_{f}}(v)=\sum_{P\in\operatorname{\textbf{P}}(v)}2^{\operatorname{cost_{f}}(P)}\textrm{ and }\operatorname{c_{f}}(v,e)=\sum_{P_{e}\in\operatorname{\textbf{P}}(v),e\in P_{e}}2^{\operatorname{cost_{f}}(P_{e})}$. Then we set $\operatorname{w}(v)={w_{all}}=1+{w_{l}}+{w_{r}}$ if $v$ is an operator node and $\operatorname{w}(v)=0$ otherwise. For an edge $e$ leading to $v$ we set

[TABLE]

We show that choosing the parameters as in (2.2) with this weight function is optimal and that it leads to the desired total evaluation cost. For a node $v\in E$ let $P\in\operatorname{\textbf{P}}(v)$ be any path to $v$ of the form $P=(v_{0},e_{0},...,v_{k},e_{k},v_{k+1}=v)$, then

[TABLE]

In particular, the precision requested at $v$ along each path is the same. Assume that the parameter choice is not optimal. For an edge $e$ let $\delta(e)$ be the difference in $\operatorname{i}(e)$ between the optimal value and the value resulting from (2.2) with weights as defined in (6) and let $\delta(v)$ be the respective difference in $\operatorname{i}(v)$ for a node $v$. Due to the optimization that led to (2.2), the slope of $\operatorname{i}(v)$ is $-\operatorname{w}({e_{l}})$ in direction of $\operatorname{i}({e_{l}})$ and $-\operatorname{w}({e_{r}})$ in direction of $\operatorname{i}({e_{r}})$ when keeping (1) equal. So the difference in $\operatorname{i}(v)$ can be bounded through

[TABLE]

Denote the difference in cost by preceeding it with $\Delta$ and let ${\textbf{E}(E)}$ be the set of edges in $E$. For our parameter choice, the precision requested at a node $v$ is the same along each path as shown in (7), so $\Delta\max_{P\in\operatorname{\textbf{P}}(v)}\operatorname{cost_{r}}(P)=\max_{P\in\operatorname{\textbf{P}}(v)}\Delta\operatorname{cost_{r}}(P)$. We then get

[TABLE]

and therefore our parameter choice is optimal. It remains to calculate the total cost for evaluating $E$. Since $\operatorname{w}(v_{0})=n$ and each path $P\in\operatorname{\textbf{P}}(v)$ leads to the same requested precision, the desired equation follows directly from (7) with

[TABLE]

Choosing the parameters as in (2.2) leads to an optimal distribution of error bounds under the assumption that the weights ${w_{l}},{w_{r}}$ accurately reflect the impact of an increase in ${i_{l}},{i_{r}}$ on the total cost. Computing the exact weight shown in (6) is hard since we have to know and to maintain the cost along all paths leading to a node. We discuss several heuristic approaches. From Theorem 2.1 we can immediately conclude:

Corollary 1

Let $T$, $|T|=n$, be an expression tree, i.e., an expression dag where each node has at most one parent. Then the optimal weight choice for an edge leading to a node $v$ is the number of operator nodes in the subexpression rooted at $v$ and the cost of an evaluation of $T$ to accuracy $q\leq 0$ is

[TABLE]

where $\operatorname{path}(v)$ denotes the unique path $P\in\operatorname{\textbf{P}}(v)$. ∎

So a natural choice for the weight of an edge is the number of operator nodes in the respective subexpression of the target node. Then the optimality condition holds for tree-like expression dags but fails when common subexpressions exist. Figure 1 shows a graph for which the optimal distribution (1(a)) differs from the distribution achieved through counting the operators (1(b)). In the example the weights for the middle node are ${w_{l}}={w_{r}}=1$. Since the lower addition is a common child of the left and the right path, it gets evaluated only once. The optimal weights would therefore be ${w_{l}}={w_{r}}=0.5$. When constants are present it may even occur that a common subexpression already needs to be evaluated at a much higher accuracy and therefore the weight can be set close to zero.

Computing the actual number of operators without duplicates in an expression dag is already a difficult task. As in Section 2.1, we can set the weight of an edge to the number of operators in the subexpression, counting duplicates, in which case we need to deal with a possible exponential increase in weight size. This leads to an additional loss in optimality (cf. Figure 1(c)), but makes it algorithmically feasible to compute the weights. This approach is largely identical to the one of van der Hoeven, who defined the weights as the number of leaves in the left and right subexpression [4]. Regarding the exponential weight increase, van der Hoeven suggested the use of a floating-point representation. Effectively managing correct floating-point bounds can get expensive. We use a different approach. In the definition of ${i_{v}},{i_{l}},{i_{r}}$ the actual value of the weights is never needed. This enables us to store the weight in a logarithmic representation from the start. The downside of this approach is that an exact computation of the weight is not possible even for small values. Note that an overestimation of the weights will never lead us to violate the condition in (1) and therefore maintains exact computation. When computing the weights, we need to compute terms of the form $\log(2^{a}+2^{b})$. Let $a\geq b$, then we have $\log(2^{a}+2^{b})=a+\log(1+2^{b-a})$ with $2^{b-a}\leq 1$. An upper bound on the logarithm can be obtained through repeated squaring [5]. For $a\gg b$ squaring $1+2^{b-a}$ is numerically unstable. In this case we can approximate the logarithm by linearization near $1$. Then

[TABLE]

and therefore $\log(2^{a}+2^{b})\leq a+\frac{1}{\ln(2)}2^{b-a}$. This approximation works well for a large difference between $a$ and $b$. For small values of $a-b$ we can use repeated squaring. Otherwise we simply set the result to $1$ for $a-b\leq\log(\ln(2))$. One way to efficiently compute an upper bound to the power term is to compute the product $2^{2^{d_{1}}}\cdots 2^{2^{d_{k}}}$ with $d_{min}\leq d_{i}\leq 0$ for $1\leq i\leq k$ where $d_{1},...,d_{k}\in\mathbb{Z}$ are the digits set to one in the binary representation of $b-a$. Since the number of possible factors is finite, we can store upper bounds for them in a lookup table.

Error bound balancing does not alter the structure of the expression dag and therefore does not change its parallelizability. The maximum cost of a critical path is reduced from $\Theta(n^{2})$ to $\Theta(n\log n)$, but multiple threads cannot be utilized effectively. If an arbitrary number of processors is available, the total cost of the evaluation reduces to the cost of evaluating a critical path. We can therefore choose the error bounds in such a way that the highest cost of a path from the root to a leaf is minimized. A lower bound on the cost of a critical path $P=(v_{0},e_{0},...,e_{k-1},v_{k})$ with $k$ operators can be obtained by isolating it, i.e., by assuming that each other edge in the expression dag leads to an operand. Let $\operatorname{cost_{C}}(P)=\sum_{i=0}^{k}\operatorname{cost_{f}}(\operatorname{path}(v_{i}))-kq$ be the cost induced by the constants and the initial accuracy along $P$. Then Corollary 1 gives

[TABLE]

If $k=n$ the weight choice is already optimal. Let ${E_{bal}}$ be an expression dag that resembles a perfectly balanced tree with depth $k$ and $2^{k}-1$ operator nodes. Since we do not have common subexpressions, $|\operatorname{\textbf{P}}(v)|=1$ for each $v\in{E_{bal}}$ and with (7) the total cost of any path $P$ in ${E_{bal}}$ is $\operatorname{cost}(P)=k\log(2^{k}-1)+\operatorname{cost_{C}}(P)=\Theta(k^{2})$. When minimizing the total cost of ${E_{bal}}$, the precision increase ${i_{v}}$ at a node $v\in{E_{bal}}$ is weighted against the cost induced in all operators in its subexpression and therefore logarithmic in their number. The cost induced on the critical path, however, depends on the depth of the subexpression. Building upon this observation, the cost of the critical path in ${E_{bal}}$ can be reduced. For a node $v\in{E_{bal}}$ with subexpression depth $j$ and outgoing edges ${e_{l}},{e_{r}}$ we set ${i_{v}}=-\log(j)$ and ${i_{l}}-\log({c_{l}})={i_{r}}-\log({c_{r}})=\log(j-1)-\log(j)-1$ (cf. (2.2)). Then the cost of the critical path $P$ in ${E_{bal}}$ is

[TABLE]

It can be shown that this parameter choice is optimal, aside from taking the operation constants into account. Although not an asymptotic improvement, the cost of the critical path was cut nearly in half. In the derivation of the chosen parameters, we made use of the symmetry of the expression. In general it is hard to compute the optimal parameters for minimizing the critical path. Let $v$ be the root node of an expression dag $X$ with outgoing edges ${e_{l}},{e_{r}}$ where the left subexpression $L$ has depth ${d_{l}}\geq 1$ and the right subexpression $R$ has depth ${d_{r}}\geq 1$. In an optimal parameter choice we have

[TABLE]

Otherwise, ${i_{l}}$ or ${i_{r}}$ could be decreased without increasing the cost of the critical path of $X$ and ${i_{v}}$ could be increased, reducing its cost. Let ${d_{f}}=\frac{{d_{l}}}{{d_{r}}}$, let ${c_{f}}=\frac{\operatorname{cost}(\operatorname{cp}(L))-\operatorname{cost}(\operatorname{cp}(R))}{{d_{r}}}$ and let ${c}=2^{c_{f}}$. Then ${i_{r}}={d_{f}}{i_{l}}+{c_{f}}$ and with (1) and $z=2^{i_{l}}$ we get ${i_{v}}=\log(1-z-{c}z^{d_{f}})$. Due to (10) there is a critical path through ${e_{l}}$ and therefore

[TABLE]

Substituting ${i_{v}}$ and forming the derivative with respect to ${i_{l}}$ we get

[TABLE]

Solving this equation yields an optimal choice for ${i_{l}}$ (and hence with (10) and (1) for ${i_{r}}$ and ${i_{v}}$). Note that for ${d_{f}}=1$, ${c_{f}}=0$ and $d_{l}=d_{r}$ we get the parameters used for ${E_{bal}}$. Unfortunately, there is no closed form for the solution of (11) for arbitrary ${d_{f}}$. Thus, for an implementation a numerical or a heuristic approach is needed.

The cost induced by operation constants and the initial accuracy usually increases with a higher depth. So it is plausible to assume for a node $v$ that the child with the higher subexpression depth will contain a more expensive path in the evaluation, if the difference in accuracy increase at $v$ is relatively small. We can use this observation in the following heuristic. We set

[TABLE]

Figure 2 shows an example for the differences between the critical path optimization, total cost optimization and the depth heuristic. The heuristic reduces the weight of the critical path compared to the previous strategies.

3 Experiments

We present experiments to underline differences between restructuring (Section 2.1) and error bound balancing (Section 2.2). For the comparison, the policy-based exact-decisions number type Real_algebraic with multithreading is used [8, 12]. We compare several different strategies. In our default configuration for Real_algebraic we use boost::interval as floating-point filter and mpfr_t as bigfloat data type. Furthermore we always enable topological evaluation, bottom-up separation bound representation and error representation by exponents [9, 13]. We call the default strategy without balancing def. For restructuring we use the weighted version of Brent’s algorithm with unit weights (bru) and with setting the weights to the expression depth (brd). For error bound balancing we use the weight function counting all operators without removing duplicates (ebc) and the depth-based approach for reducing the length of critical paths (ebd). We furthermore test combinations of internal and external balancing as described in the respective sections. For every strategy we use a variant with and without multithreading (m). The experiments are performed on an Intel i7-4700MQ with 16GB RAM under Ubuntu 18.04, using g++ 7.3.0, Boost 1.62.0 and MPFR 4.0.1. All data points are averaged over twenty runs if not specified otherwise. All expressions are evaluated to an accuracy of $q=-10000$.

3.1 List-like expression dags

List-like expression dags with linear depth have quadratic cost (cf. Section 2). Both restructuring and error bound balancing should reduce the cost significantly in this case. We build an expression dag ${E_{list}}$ by computing $res:=res\circ a_{i}$ in a simple loop starting with $res=a_{0}$, where $\circ\in\{+,-,*,/\}$ is chosen randomly and uniformly and $a_{i}$ are operands ($0\leq i\leq n$). For the operands we choose random rationals, i.e., expressions of the form $a_{i}=d_{i,1}/d_{i,2}$ where $d_{i,j}\neq 0$ are random double numbers exponentially distributed around $1$. By using exact divisions we assure that the operands have sufficient complexity for our experiments. To prevent them from being affected by restructuring, we assign an additional (external) reference to each operand. Figure 3 shows the results for evaluating ${E_{list}}$.

Both balancing methods lead to a significant reduction in running time compared to the default configuration (note the logarithmic scale). For large numbers of operators, restructuring is superior to error bound balancing. While error bound balancing optimizes the variable precision increase, it does not reduce the cost associated with the operation constants. The precision increase due to operation constants affects more nodes in an unbalanced structure than in a balanced one, which gives restructuring an advantage. For small numbers of operators, error bound balancing leads to better results than restructuring, since the cost of evaluating additional operators created through restructuring becomes more relevant. The structure of ${E_{list}}$ is highly detrimental to efficient parallelization. Consequently, neither the default evaluation nor the error bound balanced evaluation show significant cost reduction when run on multiple processors. With Brent’s algorithm a speedup of about $1.7$, i.e., a runtime reduction of about $40\text{\,}\mathrm{\char 37\relax}$, can be observed. Since ${E_{list}}$ does not contain any common subexpressions or other barriers, the results for other restructuring or error bound balancing strategies are indistinguishable from their counterparts. Interestingly, the evaluation does not benefit from a combination of both balancing strategies. Instead the results closely resemble the results obtained by using only restructuring and even get a bit worse in the multithreaded case. Since through restructuring a perfectly balanced dag is created, the default error bounds are already close to optimal (cf. Section 3.3).

3.2 Blocking nodes

Restructuring gets difficult as soon as ‘blocking nodes’, such as nodes with multiple parents, occur in the expression dag (cf. Section 2.1). We repeat the experiment from Section 3.1, but randomly let about $30\text{\,}\mathrm{\char 37\relax}$ of the operator nodes be blocking nodes by adding an additional parent (which is not part of our evaluation). Nodes with such a parent cannot be part of a restructuring process, since the subexpressions associated with them might be used somewhere else and therefore cannot be destroyed.

Both the default and the internal balancing method are not affected by the change and thus show the same results as before. Restructuring on the other hand performs worse and falls back behind error bound balancing (cf. Figure 4). The depth heuristic leads to fewer losses for both total and parallel running time. It reduces the running time by about $10\text{\,}\mathrm{\char 37\relax}$ in single-threaded and about $20\text{\,}\mathrm{\char 37\relax}$ in multithreaded execution compared to using unit weights. Combining internal and external balancing combines the advantages of both strategies in this case. Exemplarly, a combination of brd and ebc, named cmb, is shown in Figure 4. For serial evaluation the running time of the combined approach mostly resembles the running time of error bound balancing, getting slightly faster for a large number of operators (about $9\text{\,}\mathrm{\char 37\relax}$ for $N=50000$). In parallel, however, it strongly increases parallelizability leading to a speedup of $1.6$ and a total runtime reduction of up to $85\text{\,}\mathrm{\char 37\relax}$ compared to the default strategy.

3.3 Balanced expression dags

When an expression dag is already balanced, there is not much to gain by either balancing method. In a perfectly balanced expression dag ${E_{bal}}$, restructuring cannot reduce the depth and therefore does not reduce its cost, neither in serial nor in parallel. Brent’s algorithm still creates a normal form, which adds additional operations and might even increase the maximum depth. Error bound balancing on the other hand can potentially make a difference. For a balanced expression dag the total cost is strongly influenced by the operation constants, which is reflected in a high variance when choosing the operators at random. In the experiment shown in Figure 5, we increase the number of test sets for each data point from $20$ to $50$ and use the same test data for each number type. The single data points lie in a range of about $\pm$20\text{,}\mathrm{\char 37\relax}$$ of the respective average.

As expected, restructuring performs worse than the default number type, doubling the depth and replacing each division by, on average, two multiplications. Error bound balancing performs worse than not balancing as well. Neither the total operator count, nor the depth-based strategy have a significant impact on the running time of the bigfloat operations, since the cost decrease per operation is at most logarithmic in the number of operators. For the same reason and due to the limited number of processors, the expected cost reduction between ebcm and ebdm in the multithreaded case (cf. Section 2.2) can not be observed in the experimental data.

3.4 Common Subexpressions

Random expression dags, created by randomly applying operations on a forest of operands until it is reduced to a single DAG, tend to be balanced and therefore behave similarly to a perfectly balanced tree. This changes if common subexpressions are involved. With error bound balancing, common subexpressions can be recognized and the error bounds at the parent nodes can be adjusted, such that both request the same accuracies (cf. Theorem 2.1). The two implemented heuristics to some degree take common subexpressions into account, since they contribute the same weight to all of the subexpression’s parents. We test the behavior of error bound balancing strategies by randomly reusing a certain percentage of subtrees during randomized bottom-up construction of the graph. To avoid zeros, ones, or an exponential explosion of the expression’s value we only use additions if two subtrees are identical during construction. While error bound balancing still cannot outperform the default strategy due to the balanced nature, it moves on par with it. If $5\text{\,}\mathrm{\char 37\relax}$ of the operations have more than one parent, the error bound balancing strategies improve the single-threaded running time by about $1\text{\,}\mathrm{\char 37\relax}5\text{\,}\mathrm{\char 37\relax}$. In a parallel environment, it still performs worse with ebdm being slightly superior to ebcm.

If common subexpressions lead to a large difference between the actual number of operators and the number of operators in a tree expansion, ebc significantly overestimates the optimal weight of its edges (cf. Figure 1). Figure 6 shows results for evaluating a sequence of additions where the left and the right summand is the result of the previous addition. The full operator count heuristic does not reduce the running time and even performs worse than the default strategy for large numbers of operators, whereas the depth-based heuristic clearly outperforms the other strategies. Note that in this case, the depth-based heuristic leads to the optimal error distribution for both total and critical path cost.

3.5 A note on floating-point primitives

Error bound balancing requires the use of floating-point error bounds. While IEEE 754 requires that floating point computations must be exactly rounded, it is surprisingly difficult to find an adequate upper or lower bound to the result of such an operation. IEEE 754 specifies four rounding modes: Round to nearest, Round to positive/negative infinty and Round to zero [1]. For the last three modes, which are commonly referred to as directed rounding, it is easy to obtain a lower or upper bound by negating the operands adequately. Unfortunately, most systems implement Round to nearest. Switching the rounding mode is expensive. While double operations with appropriate negations for directed rounding are about two times slower, switching to an appropriate rounding mode can increase the running time of a single operation by a factor of $100$. The same factor applies if we manually jump to the next (or previous) representable double value.

Handling floating-point primitives correctly can, depending on the architecture, be very expensive. In most cases, however, the computed error bounds massively overestimate the actual error. Moreover, for the actual bigfloats computations the error bounds are rounded up to the next integer. It is therefore almost impossible that floating-point rounding errors make an actual difference in any computation. For our experiments we refrained from handling those bounds correctly to make the results more meaningful and less architecture-dependent.

4 Conclusion

We have shown, theoretically and experimentally, that both external and internal balancing methods are useful tools to mitigate the impact of badly balanced expression dags. Restructuring has a higher potential on reducing the cost, but can become useless or even detrimental if the graph has many common subexpressions or is already balanced. In a parallel environment, restructuring is necessary to make use of multiple processors in an unbalanced graph. Error bound balancing is more widely applicable, but is limited in its effectivity. If the graph is small or already sufficiently balanced, neither of the methods has a significant positive impact on the evaluation cost. A general purpose number type should therefore always check whether the structure generally requires balancing before applying either of the algorithms. For both strategies we have described optimal weight functions. In both cases implementations require heuristics to be practicable. Our experiments show that carefully chosen heuristics are in most cases sufficient to increase the performance of exact number types.