| \section{Sets and Relations} |
| |
| \begin{definition}[Polyhedral Set] |
| A {\em polyhedral set}\index{polyhedral set} $S$ is a finite union of basic sets |
| $S = \bigcup_i S_i$, each of which can be represented using affine |
| constraints |
| $$ |
| S_i : \Z^n \to 2^{\Z^d} : \vec s \mapsto |
| S_i(\vec s) = |
| \{\, \vec x \in \Z^d \mid \exists \vec z \in \Z^e : |
| A \vec x + B \vec s + D \vec z + \vec c \geq \vec 0 \,\} |
| , |
| $$ |
| with $A \in \Z^{m \times d}$, |
| $B \in \Z^{m \times n}$, |
| $D \in \Z^{m \times e}$ |
| and $\vec c \in \Z^m$. |
| \end{definition} |
| |
| \begin{definition}[Parameter Domain of a Set] |
| Let $S \in \Z^n \to 2^{\Z^d}$ be a set. |
| The {\em parameter domain} of $S$ is the set |
| $$\pdom S \coloneqq \{\, \vec s \in \Z^n \mid S(\vec s) \ne \emptyset \,\}.$$ |
| \end{definition} |
| |
| \begin{definition}[Polyhedral Relation] |
| A {\em polyhedral relation}\index{polyhedral relation} |
| $R$ is a finite union of basic relations |
| $R = \bigcup_i R_i$ of type |
| $\Z^n \to 2^{\Z^{d_1+d_2}}$, |
| each of which can be represented using affine |
| constraints |
| $$ |
| R_i = \vec s \mapsto |
| R_i(\vec s) = |
| \{\, \vec x_1 \to \vec x_2 \in \Z^{d_1} \times \Z^{d_2} |
| \mid \exists \vec z \in \Z^e : |
| A_1 \vec x_1 + A_2 \vec x_2 + B \vec s + D \vec z + \vec c \geq \vec 0 \,\} |
| , |
| $$ |
| with $A_i \in \Z^{m \times d_i}$, |
| $B \in \Z^{m \times n}$, |
| $D \in \Z^{m \times e}$ |
| and $\vec c \in \Z^m$. |
| \end{definition} |
| |
| \begin{definition}[Parameter Domain of a Relation] |
| Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation. |
| The {\em parameter domain} of $R$ is the set |
| $$\pdom R \coloneqq \{\, \vec s \in \Z^n \mid R(\vec s) \ne \emptyset \,\}.$$ |
| \end{definition} |
| |
| \begin{definition}[Domain of a Relation] |
| Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation. |
| The {\em domain} of $R$ is the polyhedral set |
| $$\domain R \coloneqq \vec s \mapsto |
| \{\, \vec x_1 \in \Z^{d_1} \mid \exists \vec x_2 \in \Z^{d_2} : |
| (\vec x_1, \vec x_2) \in R(\vec s) \,\} |
| . |
| $$ |
| \end{definition} |
| |
| \begin{definition}[Range of a Relation] |
| Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation. |
| The {\em range} of $R$ is the polyhedral set |
| $$ |
| \range R \coloneqq \vec s \mapsto |
| \{\, \vec x_2 \in \Z^{d_2} \mid \exists \vec x_1 \in \Z^{d_1} : |
| (\vec x_1, \vec x_2) \in R(\vec s) \,\} |
| . |
| $$ |
| \end{definition} |
| |
| \begin{definition}[Composition of Relations] |
| Let $R \in \Z^n \to 2^{\Z^{d_1+d_2}}$ and |
| $S \in \Z^n \to 2^{\Z^{d_2+d_3}}$ be two relations, |
| then the composition of |
| $R$ and $S$ is defined as |
| $$ |
| S \circ R \coloneqq |
| \vec s \mapsto |
| \{\, \vec x_1 \to \vec x_3 \in \Z^{d_1} \times \Z^{d_3} |
| \mid \exists \vec x_2 \in \Z^{d_2} : |
| \vec x_1 \to \vec x_2 \in R(\vec s) \wedge |
| \vec x_2 \to \vec x_3 \in S(\vec s) |
| \,\} |
| . |
| $$ |
| \end{definition} |
| |
| \begin{definition}[Difference Set of a Relation] |
| Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation. |
| The difference set ($\Delta \, R$) of $R$ is the set |
| of differences between image elements and the corresponding |
| domain elements, |
| $$ |
| \diff R \coloneqq |
| \vec s \mapsto |
| \{\, \vec \delta \in \Z^{d} \mid \exists \vec x \to \vec y \in R : |
| \vec \delta = \vec y - \vec x |
| \,\} |
| $$ |
| \end{definition} |
| |
| \section{Simple Hull}\label{s:simple hull} |
| |
| It is sometimes useful to have a single |
| basic set or basic relation that contains a given set or relation. |
| For rational sets, the obvious choice would be to compute the |
| (rational) convex hull. For integer sets, the obvious choice |
| would be the integer hull. |
| However, {\tt isl} currently does not support an integer hull operation |
| and even if it did, it would be fairly expensive to compute. |
| The convex hull operation is supported, but it is also fairly |
| expensive to compute given only an implicit representation. |
| |
| Usually, it is not required to compute the exact integer hull, |
| and an overapproximation of this hull is sufficient. |
| The ``simple hull'' of a set is such an overapproximation |
| and it is defined as the (inclusion-wise) smallest basic set |
| that is described by constraints that are translates of |
| the constraints in the input set. |
| This means that the simple hull is relatively cheap to compute |
| and that the number of constraints in the simple hull is no |
| larger than the number of constraints in the input. |
| \begin{definition}[Simple Hull of a Set] |
| The {\em simple hull} of a set |
| $S = \bigcup_{1 \le i \le v} S_i$, with |
| $$ |
| S : \Z^n \to 2^{\Z^d} : \vec s \mapsto |
| S(\vec s) = |
| \left\{\, \vec x \in \Z^d \mid \exists \vec z \in \Z^e : |
| \bigvee_{1 \le i \le v} |
| A_i \vec x + B_i \vec s + D_i \vec z + \vec c_i \geq \vec 0 \,\right\} |
| $$ |
| is the set |
| $$ |
| H : \Z^n \to 2^{\Z^d} : \vec s \mapsto |
| S(\vec s) = |
| \left\{\, \vec x \in \Z^d \mid \exists \vec z \in \Z^e : |
| \bigwedge_{1 \le i \le v} |
| A_i \vec x + B_i \vec s + D_i \vec z + \vec c_i + \vec K_i \geq \vec 0 |
| \,\right\} |
| , |
| $$ |
| with $\vec K_i$ the (component-wise) smallest non-negative integer vectors |
| such that $S \subseteq H$. |
| \end{definition} |
| The $\vec K_i$ can be obtained by solving a number of |
| LP problems, one for each element of each $\vec K_i$. |
| If any LP problem is unbounded, then the corresponding constraint |
| is dropped. |
| |
| \section{Parametric Integer Programming} |
| |
| \subsection{Introduction}\label{s:intro} |
| |
| Parametric integer programming \parencite{Feautrier88parametric} |
| is used to solve many problems within the context of the polyhedral model. |
| Here, we are mainly interested in dependence analysis \parencite{Fea91} |
| and in computing a unique representation for existentially quantified |
| variables. The latter operation has been used for counting elements |
| in sets involving such variables |
| \parencite{BouletRe98,Verdoolaege2005experiences} and lies at the core |
| of the internal representation of {\tt isl}. |
| |
| Parametric integer programming was first implemented in \texttt{PipLib}. |
| An alternative method for parametric integer programming |
| was later implemented in {\tt barvinok} \cite{barvinok-0.22}. |
| This method is not based on Feautrier's algorithm, but on rational |
| generating functions \cite{Woods2003short} and was inspired by the |
| ``digging'' technique of \textcite{DeLoera2004Three} for solving |
| non-parametric integer programming problems. |
| |
| In the following sections, we briefly recall the dual simplex |
| method combined with Gomory cuts and describe some extensions |
| and optimizations. The main algorithm is applied to a matrix |
| data structure known as a tableau. In case of parametric problems, |
| there are two tableaus, one for the main problem and one for |
| the constraints on the parameters, known as the context tableau. |
| The handling of the context tableau is described in \autoref{s:context}. |
| |
| \subsection{The Dual Simplex Method} |
| |
| Tableaus can be represented in several slightly different ways. |
| In {\tt isl}, the dual simplex method uses the same representation |
| as that used by its incremental LP solver based on the \emph{primal} |
| simplex method. The implementation of this LP solver is based |
| on that of {\tt Simplify} \parencite{Detlefs2005simplify}, which, in turn, |
| was derived from the work of \textcite{Nelson1980phd}. |
| In the original \parencite{Nelson1980phd}, the tableau was implemented |
| as a sparse matrix, but neither {\tt Simplify} nor the current |
| implementation of {\tt isl} does so. |
| |
| Given some affine constraints on the variables, |
| $A \vec x + \vec b \ge \vec 0$, the tableau represents the relationship |
| between the variables $\vec x$ and non-negative variables |
| $\vec y = A \vec x + \vec b$ corresponding to the constraints. |
| The initial tableau contains $\begin{pmatrix} |
| \vec b & A |
| \end{pmatrix}$ and expresses the constraints $\vec y$ in the rows in terms |
| of the variables $\vec x$ in the columns. The main operation defined |
| on a tableau exchanges a column and a row variable and is called a pivot. |
| During this process, some coefficients may become rational. |
| As in the \texttt{PipLib} implementation, |
| {\tt isl} maintains a shared denominator per row. |
| The sample value of a tableau is one where each column variable is assigned |
| zero and each row variable is assigned the constant term of the row. |
| This sample value represents a valid solution if each constraint variable |
| is assigned a non-negative value, i.e., if the constant terms of |
| rows corresponding to constraints are all non-negative. |
| |
| The dual simplex method starts from an initial sample value that |
| may be invalid, but that is known to be (lexicographically) no |
| greater than any solution, and gradually increments this sample value |
| through pivoting until a valid solution is obtained. |
| In particular, each pivot exchanges a row variable |
| $r = -n + \sum_i a_i \, c_i$ with negative |
| sample value $-n$ with a column variable $c_j$ |
| such that $a_j > 0$. Since $c_j = (n + r - \sum_{i\ne j} a_i \, c_i)/a_j$, |
| the new row variable will have a positive sample value $n$. |
| If no such column can be found, then the problem is infeasible. |
| By always choosing the column that leads to the (lexicographically) |
| smallest increment in the variables $\vec x$, |
| the first solution found is guaranteed to be the (lexicographically) |
| minimal solution \cite{Feautrier88parametric}. |
| In order to be able to determine the smallest increment, the tableau |
| is (implicitly) extended with extra rows defining the original |
| variables in terms of the column variables. |
| If we assume that all variables are non-negative, then we know |
| that the zero vector is no greater than the minimal solution and |
| then the initial extended tableau looks as follows. |
| $$ |
| \begin{tikzpicture} |
| \matrix (m) [matrix of math nodes] |
| { |
| & {} & 1 & \vec c \\ |
| \vec x && |(top)| \vec 0 & I \\ |
| \vec r && \vec b & |(bottom)|A \\ |
| }; |
| \begin{pgfonlayer}{background} |
| \node (core) [inner sep=0pt,fill=black!20,right delimiter=),left delimiter=(,fit=(top)(bottom)] {}; |
| \end{pgfonlayer} |
| \end{tikzpicture} |
| $$ |
| Each column in this extended tableau is lexicographically positive |
| and will remain so because of the column choice explained above. |
| It is then clear that the value of $\vec x$ will increase in each step. |
| Note that there is no need to store the extra rows explicitly. |
| If a given $x_i$ is a column variable, then the corresponding row |
| is the unit vector $e_i$. If, on the other hand, it is a row variable, |
| then the row already appears somewhere else in the tableau. |
| |
| In case of parametric problems, the sign of the constant term |
| may depend on the parameters. Each time the constant term of a constraint row |
| changes, we therefore need to check whether the new term can attain |
| negative and/or positive values over the current set of possible |
| parameter values, i.e., the context. |
| If all these terms can only attain non-negative values, the current |
| state of the tableau represents a solution. If one of the terms |
| can only attain non-positive values and is not identically zero, |
| the corresponding row can be pivoted. |
| Otherwise, we pick one of the terms that can attain both positive |
| and negative values and split the context into a part where |
| it only attains non-negative values and a part where it only attains |
| negative values. |
| |
| \subsection{Gomory Cuts} |
| |
| The solution found by the dual simplex method may have |
| non-integral coordinates. If so, some rational solutions |
| (including the current sample value), can be cut off by |
| applying a (parametric) Gomory cut. |
| Let $r = b(\vec p) + \sp {\vec a} {\vec c}$ be the row |
| corresponding to the first non-integral coordinate of $\vec x$, |
| with $b(\vec p)$ the constant term, an affine expression in the |
| parameters $\vec p$, i.e., $b(\vec p) = \sp {\vec f} {\vec p} + g$. |
| Note that only row variables can attain |
| non-integral values as the sample value of the column variables is zero. |
| Consider the expression |
| $b(\vec p) - \ceil{b(\vec p)} + \sp {\fract{\vec a}} {\vec c}$, |
| with $\ceil\cdot$ the ceiling function and $\fract\cdot$ the |
| fractional part. This expression is negative at the sample value |
| since $\vec c = \vec 0$ and $r = b(\vec p)$ is fractional, i.e., |
| $\ceil{b(\vec p)} > b(\vec p)$. On the other hand, for each integral |
| value of $r$ and $\vec c \ge 0$, the expression is non-negative |
| because $b(\vec p) - \ceil{b(\vec p)} > -1$. |
| Imposing this expression to be non-negative therefore does not |
| invalidate any integral solutions, while it does cut away the current |
| fractional sample value. To be able to formulate this constraint, |
| a new variable $q = \floor{-b(\vec p)} = - \ceil{b(\vec p)}$ is added |
| to the context. This integral variable is uniquely defined by the constraints |
| $0 \le -d \, b(\vec p) - d \, q \le d - 1$, with $d$ the common |
| denominator of $\vec f$ and $g$. In practice, the variable |
| $q' = \floor{\sp {\fract{-f}} {\vec p} + \fract{-g}}$ is used instead |
| and the coefficients of the new constraint are adjusted accordingly. |
| The sign of the constant term of this new constraint need not be determined |
| as it is non-positive by construction. |
| When several of these extra context variables are added, it is important |
| to avoid adding duplicates. |
| Recent versions of {\tt PipLib} also check for such duplicates. |
| |
| \subsection{Negative Unknowns and Maximization} |
| |
| There are two places in the above algorithm where the unknowns $\vec x$ |
| are assumed to be non-negative: the initial tableau starts from |
| sample value $\vec x = \vec 0$ and $\vec c$ is assumed to be non-negative |
| during the construction of Gomory cuts. |
| To deal with negative unknowns, \textcite[Appendix A.2]{Fea91} |
| proposed to use a ``big parameter'', say $M$, that is taken to be |
| an arbitrarily large positive number. Instead of looking for the |
| lexicographically minimal value of $\vec x$, we search instead |
| for the lexicographically minimal value of $\vec x' = \vec M + \vec x$. |
| The sample value $\vec x' = \vec 0$ of the initial tableau then |
| corresponds to $\vec x = -\vec M$, which is clearly not greater than |
| any potential solution. The sign of the constant term of a row |
| is determined lexicographically, with the coefficient of $M$ considered |
| first. That is, if the coefficient of $M$ is not zero, then its sign |
| is the sign of the entire term. Otherwise, the sign is determined |
| by the remaining affine expression in the parameters. |
| If the original problem has a bounded optimum, then the final sample |
| value will be of the form $\vec M + \vec v$ and the optimal value |
| of the original problem is then $\vec v$. |
| Maximization problems can be handled in a similar way by computing |
| the minimum of $\vec M - \vec x$. |
| |
| When the optimum is unbounded, the optimal value computed for |
| the original problem will involve the big parameter. |
| In the original implementation of {\tt PipLib}, the big parameter could |
| even appear in some of the extra variables $\vec q$ created during |
| the application of a Gomory cut. The final result could then contain |
| implicit conditions on the big parameter through conditions on such |
| $\vec q$ variables. This problem was resolved in later versions |
| of {\tt PipLib} by taking $M$ to be divisible by any positive number. |
| The big parameter can then never appear in any $\vec q$ because |
| $\fract {\alpha M } = 0$. It should be noted, though, that an unbounded |
| problem usually (but not always) |
| indicates an incorrect formulation of the problem. |
| |
| The original version of {\tt PipLib} required the user to ``manually'' |
| add a big parameter, perform the reformulation and interpret the result |
| \parencite{Feautrier02}. Recent versions allow the user to simply |
| specify that the unknowns may be negative or that the maximum should |
| be computed and then these transformations are performed internally. |
| Although there are some application, e.g., |
| that of \textcite{Feautrier92multi}, |
| where it is useful to have explicit control over the big parameter, |
| negative unknowns and maximization are by far the most common applications |
| of the big parameter and we believe that the user should not be bothered |
| with such implementation issues. |
| The current version of {\tt isl} therefore does not |
| provide any interface for specifying big parameters. Instead, the user |
| can specify whether a maximum needs to be computed and no assumptions |
| are made on the sign of the unknowns. Instead, the sign of the unknowns |
| is checked internally and a big parameter is automatically introduced when |
| needed. For compatibility with {\tt PipLib}, the {\tt isl\_pip} tool |
| does explicitly add non-negativity constraints on the unknowns unless |
| the \verb+Urs_unknowns+ option is specified. |
| Currently, there is also no way in {\tt isl} of expressing a big |
| parameter in the output. Even though |
| {\tt isl} makes the same divisibility assumption on the big parameter |
| as recent versions of {\tt PipLib}, it will therefore eventually |
| produce an error if the problem turns out to be unbounded. |
| |
| \subsection{Preprocessing} |
| |
| In this section, we describe some transformations that are |
| or can be applied in advance to reduce the running time |
| of the actual dual simplex method with Gomory cuts. |
| |
| \subsubsection{Feasibility Check and Detection of Equalities} |
| |
| Experience with the original {\tt PipLib} has shown that Gomory cuts |
| do not perform very well on problems that are (non-obviously) empty, |
| i.e., problems with rational solutions, but no integer solutions. |
| In {\tt isl}, we therefore first perform a feasibility check on |
| the original problem considered as a non-parametric problem |
| over the combined space of unknowns and parameters. |
| In fact, we do not simply check the feasibility, but we also |
| check for implicit equalities among the integer points by computing |
| the integer affine hull. The algorithm used is the same as that |
| described in \autoref{s:GBR} below. |
| Computing the affine hull is fairly expensive, but it can |
| bring huge benefits if any equalities can be found or if the problem |
| turns out to be empty. |
| |
| \subsubsection{Constraint Simplification} |
| |
| If the coefficients of the unknown and parameters in a constraint |
| have a common factor, then this factor should be removed, possibly |
| rounding down the constant term. For example, the constraint |
| $2 x - 5 \ge 0$ should be simplified to $x - 3 \ge 0$. |
| {\tt isl} performs such simplifications on all sets and relations. |
| Recent versions of {\tt PipLib} also perform this simplification |
| on the input. |
| |
| \subsubsection{Exploiting Equalities}\label{s:equalities} |
| |
| If there are any (explicit) equalities in the input description, |
| {\tt PipLib} converts each into a pair of inequalities. |
| It is also possible to write $r$ equalities as $r+1$ inequalities |
| \parencite{Feautrier02}, but it is even better to \emph{exploit} the |
| equalities to reduce the dimensionality of the problem. |
| Given an equality involving at least one unknown, we pivot |
| the row corresponding to the equality with the column corresponding |
| to the last unknown with non-zero coefficient. The new column variable |
| can then be removed completely because it is identically zero, |
| thereby reducing the dimensionality of the problem by one. |
| The last unknown is chosen to ensure that the columns of the initial |
| tableau remain lexicographically positive. In particular, if |
| the equality is of the form $b + \sum_{i \le j} a_i \, x_i = 0$ with |
| $a_j \ne 0$, then the (implicit) top rows of the initial tableau |
| are changed as follows |
| $$ |
| \begin{tikzpicture} |
| \matrix [matrix of math nodes] |
| { |
| & {} & |(top)| 0 & I_1 & |(j)| & \\ |
| j && 0 & & 1 & \\ |
| && 0 & & & |(bottom)|I_2 \\ |
| }; |
| \node[overlay,above=2mm of j,anchor=south]{j}; |
| \begin{pgfonlayer}{background} |
| \node (m) [inner sep=0pt,fill=black!20,right delimiter=),left delimiter=(,fit=(top)(bottom)] {}; |
| \end{pgfonlayer} |
| \begin{scope}[xshift=4cm] |
| \matrix [matrix of math nodes] |
| { |
| & {} & |(top)| 0 & I_1 & \\ |
| j && |(left)| -b/a_j & -a_i/a_j & \\ |
| && 0 & & |(bottom)|I_2 \\ |
| }; |
| \begin{pgfonlayer}{background} |
| \node (m2) [inner sep=0pt,fill=black!20,right delimiter=),left delimiter=(,fit=(top)(bottom)(left)] {}; |
| \end{pgfonlayer} |
| \end{scope} |
| \draw [shorten >=7mm,-to,thick,decorate, |
| decoration={snake,amplitude=.4mm,segment length=2mm, |
| pre=moveto,pre length=5mm,post length=8mm}] |
| (m) -- (m2); |
| \end{tikzpicture} |
| $$ |
| Currently, {\tt isl} also eliminates equalities involving only parameters |
| in a similar way, provided at least one of the coefficients is equal to one. |
| The application of parameter compression (see below) |
| would obviate the need for removing parametric equalities. |
| |
| \subsubsection{Offline Symmetry Detection}\label{s:offline} |
| |
| Some problems, notably those of \textcite{Bygde2010licentiate}, |
| have a collection of constraints, say |
| $b_i(\vec p) + \sp {\vec a} {\vec x} \ge 0$, |
| that only differ in their (parametric) constant terms. |
| These constant terms will be non-negative on different parts |
| of the context and this context may have to be split for each |
| of the constraints. In the worst case, the basic algorithm may |
| have to consider all possible orderings of the constant terms. |
| Instead, {\tt isl} introduces a new parameter, say $u$, and |
| replaces the collection of constraints by the single |
| constraint $u + \sp {\vec a} {\vec x} \ge 0$ along with |
| context constraints $u \le b_i(\vec p)$. |
| Any solution to the new system is also a solution |
| to the original system since |
| $\sp {\vec a} {\vec x} \ge -u \ge -b_i(\vec p)$. |
| Conversely, $m = \min_i b_i(\vec p)$ satisfies the constraints |
| on $u$ and therefore extends a solution to the new system. |
| It can also be plugged into a new solution. |
| See \autoref{s:post} for how this substitution is currently performed |
| in {\tt isl}. |
| The method described in this section can only detect symmetries |
| that are explicitly available in the input. |
| See \autoref{s:online} for the detection |
| and exploitation of symmetries that appear during the course of |
| the dual simplex method. |
| |
| Note that the replacement of the $b_i(\vec p)$ by $u$ may lose |
| information if the parameters that occur in $b_i(\vec p)$ also |
| occur in other constraints. The replacement is therefore currently |
| only applied when all the parameters in all of the $b_i(\vec p)$ |
| only occur in a single constraint, i.e., the one in which |
| the parameter is removed. |
| This is the case for the examples from \textcite{Bygde2010licentiate} |
| in \autoref{t:comparison}. |
| The version of {\tt isl} that was used during the experiments |
| of \autoref{s:pip:experiments} did not take into account |
| this single-occurrence constraint. |
| |
| \subsubsection{Parameter Compression}\label{s:compression} |
| |
| It may in some cases be apparent from the equalities in the problem |
| description that there can only be a solution for a sublattice |
| of the parameters. In such cases ``parameter compression'' |
| \parencite{Meister2004PhD,Meister2008} can be used to replace |
| the parameters by alternative ``dense'' parameters. |
| For example, if there is a constraint $2x = n$, then the system |
| will only have solutions for even values of $n$ and $n$ can be replaced |
| by $2n'$. Similarly, the parameters $n$ and $m$ in a system with |
| the constraint $2n = 3m$ can be replaced by a single parameter $n'$ |
| with $n=3n'$ and $m=2n'$. |
| It is also possible to perform a similar compression on the unknowns, |
| but it would be more complicated as the compression would have to |
| preserve the lexicographical order. Moreover, due to our handling |
| of equalities described above there should be |
| no need for such variable compression. |
| Although parameter compression has been implemented in {\tt isl}, |
| it is currently not yet used during parametric integer programming. |
| |
| \subsection{Postprocessing}\label{s:post} |
| |
| The output of {\tt PipLib} is a quast (quasi-affine selection tree). |
| Each internal node in this tree corresponds to a split of the context |
| based on a parametric constant term in the main tableau with indeterminate |
| sign. Each of these nodes may introduce extra variables in the context |
| corresponding to integer divisions. Each leaf of the tree prescribes |
| the solution in that part of the context that satisfies all the conditions |
| on the path leading to the leaf. |
| Such a quast is a very economical way of representing the solution, but |
| it would not be suitable as the (only) internal representation of |
| sets and relations in {\tt isl}. Instead, {\tt isl} represents |
| the constraints of a set or relation in disjunctive normal form. |
| The result of a parametric integer programming problem is then also |
| converted to this internal representation. Unfortunately, the conversion |
| to disjunctive normal form can lead to an explosion of the size |
| of the representation. |
| In some cases, this overhead would have to be paid anyway in subsequent |
| operations, but in other cases, especially for outside users that just |
| want to solve parametric integer programming problems, we would like |
| to avoid this overhead in future. That is, we are planning on introducing |
| quasts or a related representation as one of several possible internal |
| representations and on allowing the output of {\tt isl\_pip} to optionally |
| be printed as a quast. |
| |
| Currently, {\tt isl} also does not have an internal representation |
| for expressions such as $\min_i b_i(\vec p)$ from the offline |
| symmetry detection of \autoref{s:offline}. |
| Assume that one of these expressions has $n$ bounds $b_i(\vec p)$. |
| If the expression |
| does not appear in the affine expression describing the solution, |
| but only in the constraints, and if moreover, the expression |
| only appears with a positive coefficient, i.e., |
| $\min_i b_i(\vec p) \ge f_j(\vec p)$, then each of these constraints |
| can simply be reduplicated $n$ times, once for each of the bounds. |
| Otherwise, a conversion to disjunctive normal form |
| leads to $n$ cases, each described as $u = b_i(\vec p)$ with constraints |
| $b_i(\vec p) \le b_j(\vec p)$ for $j > i$ |
| and |
| $b_i(\vec p) < b_j(\vec p)$ for $j < i$. |
| Note that even though this conversion leads to a size increase |
| by a factor of $n$, not detecting the symmetry could lead to |
| an increase by a factor of $n!$ if all possible orderings end up being |
| considered. |
| |
| \subsection{Context Tableau}\label{s:context} |
| |
| The main operation that a context tableau needs to provide is a test |
| on the sign of an affine expression over the elements of the context. |
| This sign can be determined by solving two integer linear feasibility |
| problems, one with a constraint added to the context that enforces |
| the expression to be non-negative and one where the expression is |
| negative. As already mentioned by \textcite{Feautrier88parametric}, |
| any integer linear feasibility solver could be used, but the {\tt PipLib} |
| implementation uses a recursive call to the dual simplex with Gomory |
| cuts algorithm to determine the feasibility of a context. |
| In {\tt isl}, two ways of handling the context have been implemented, |
| one that performs the recursive call and one, used by default, that |
| uses generalized basis reduction. |
| We start with some optimizations that are shared between the two |
| implementations and then discuss additional details of each of them. |
| |
| \subsubsection{Maintaining Witnesses}\label{s:witness} |
| |
| A common feature of both integer linear feasibility solvers is that |
| they will not only say whether a set is empty or not, but if the set |
| is non-empty, they will also provide a \emph{witness} for this result, |
| i.e., a point that belongs to the set. By maintaining a list of such |
| witnesses, we can avoid many feasibility tests during the determination |
| of the signs of affine expressions. In particular, if the expression |
| evaluates to a positive number on some of these points and to a negative |
| number on some others, then no feasibility test needs to be performed. |
| If all the evaluations are non-negative, we only need to check for the |
| possibility of a negative value and similarly in case of all |
| non-positive evaluations. Finally, in the rare case that all points |
| evaluate to zero or at the start, when no points have been collected yet, |
| one or two feasibility tests need to be performed depending on the result |
| of the first test. |
| |
| When a new constraint is added to the context, the points that |
| violate the constraint are temporarily removed. They are reconsidered |
| when we backtrack over the addition of the constraint, as they will |
| satisfy the negation of the constraint. It is only when we backtrack |
| over the addition of the points that they are finally removed completely. |
| When an extra integer division is added to the context, |
| the new coordinates of the |
| witnesses can easily be computed by evaluating the integer division. |
| The idea of keeping track of witnesses was first used in {\tt barvinok}. |
| |
| \subsubsection{Choice of Constant Term on which to Split} |
| |
| Recall that if there are no rows with a non-positive constant term, |
| but there are rows with an indeterminate sign, then the context |
| needs to be split along the constant term of one of these rows. |
| If there is more than one such row, then we need to choose which row |
| to split on first. {\tt PipLib} uses a heuristic based on the (absolute) |
| sizes of the coefficients. In particular, it takes the largest coefficient |
| of each row and then selects the row where this largest coefficient is smaller |
| than those of the other rows. |
| |
| In {\tt isl}, we take that row for which non-negativity of its constant |
| term implies non-negativity of as many of the constant terms of the other |
| rows as possible. The intuition behind this heuristic is that on the |
| positive side, we will have fewer negative and indeterminate signs, |
| while on the negative side, we need to perform a pivot, which may |
| affect any number of rows meaning that the effect on the signs |
| is difficult to predict. This heuristic is of course much more |
| expensive to evaluate than the heuristic used by {\tt PipLib}. |
| More extensive tests are needed to evaluate whether the heuristic is worthwhile. |
| |
| \subsubsection{Dual Simplex + Gomory Cuts} |
| |
| When a new constraint is added to the context, the first steps |
| of the dual simplex method applied to this new context will be the same |
| or at least very similar to those taken on the original context, i.e., |
| before the constraint was added. In {\tt isl}, we therefore apply |
| the dual simplex method incrementally on the context and backtrack |
| to a previous state when a constraint is removed again. |
| An initial implementation that was never made public would also |
| keep the Gomory cuts, but the current implementation backtracks |
| to before the point where Gomory cuts are added before adding |
| an extra constraint to the context. |
| Keeping the Gomory cuts has the advantage that the sample value |
| is always an integer point and that this point may also satisfy |
| the new constraint. However, due to the technique of maintaining |
| witnesses explained above, |
| we would not perform a feasibility test in such cases and then |
| the previously added cuts may be redundant, possibly resulting |
| in an accumulation of a large number of cuts. |
| |
| If the parameters may be negative, then the same big parameter trick |
| used in the main tableau is applied to the context. This big parameter |
| is of course unrelated to the big parameter from the main tableau. |
| Note that it is not a requirement for this parameter to be ``big'', |
| but it does allow for some code reuse in {\tt isl}. |
| In {\tt PipLib}, the extra parameter is not ``big'', but this may be because |
| the big parameter of the main tableau also appears |
| in the context tableau. |
| |
| Finally, it was reported by \textcite{Galea2009personal}, who |
| worked on a parametric integer programming implementation |
| in {\tt PPL} \parencite{PPL}, |
| that it is beneficial to add cuts for \emph{all} rational coordinates |
| in the context tableau. Based on this report, |
| the initial {\tt isl} implementation was adapted accordingly. |
| |
| \subsubsection{Generalized Basis Reduction}\label{s:GBR} |
| |
| The default algorithm used in {\tt isl} for feasibility checking |
| is generalized basis reduction \parencite{Cook1991implementation}. |
| This algorithm is also used in the {\tt barvinok} implementation. |
| The algorithm is fairly robust, but it has some overhead. |
| We therefore try to avoid calling the algorithm in easy cases. |
| In particular, we incrementally keep track of points for which |
| the entire unit hypercube positioned at that point lies in the context. |
| This set is described by translates of the constraints of the context |
| and if (rationally) non-empty, any rational point |
| in the set can be rounded up to yield an integer point in the context. |
| |
| A restriction of the algorithm is that it only works on bounded sets. |
| The affine hull of the recession cone therefore needs to be projected |
| out first. As soon as the algorithm is invoked, we then also |
| incrementally keep track of this recession cone. The reduced basis |
| found by one call of the algorithm is also reused as initial basis |
| for the next call. |
| |
| Some problems lead to the |
| introduction of many integer divisions. Within a given context, |
| some of these integer divisions may be equal to each other, even |
| if the expressions are not identical, or they may be equal to some |
| affine combination of other variables. |
| To detect such cases, we compute the affine hull of the context |
| each time a new integer division is added. The algorithm used |
| for computing this affine hull is that of \textcite{Karr1976affine}, |
| while the points used in this algorithm are obtained by performing |
| integer feasibility checks on that part of the context outside |
| the current approximation of the affine hull. |
| The list of witnesses is used to construct an initial approximation |
| of the hull, while any extra points found during the construction |
| of the hull is added to this list. |
| Any equality found in this way that expresses an integer division |
| as an \emph{integer} affine combination of other variables is |
| propagated to the main tableau, where it is used to eliminate that |
| integer division. |
| |
| \subsection{Experiments}\label{s:pip:experiments} |
| |
| \autoref{t:comparison} compares the execution times of {\tt isl} |
| (with both types of context tableau) |
| on some more difficult instances to those of other tools, |
| run on an Intel Xeon W3520 @ 2.66GHz. |
| These instances are available in the \lstinline{testsets/pip} directory |
| of the {\tt isl} distribution. |
| Easier problems such as the |
| test cases distributed with {\tt Pip\-Lib} can be solved so quickly |
| that we would only be measuring overhead such as input/output and conversions |
| and not the running time of the actual algorithm. |
| We compare the following versions: |
| {\tt piplib-1.4.0-5-g0132fd9}, |
| {\tt barvinok-0.32.1-73-gc5d7751}, |
| {\tt isl-0.05.1-82-g3a37260} |
| and {\tt PPL} version 0.11.2. |
| |
| The first test case is the following dependence analysis problem |
| originating from the Phideo project \parencite{Verhaegh1995PhD} |
| that was communicated to us by Bart Kienhuis: |
| \begin{lstlisting}[flexiblecolumns=true,breaklines=true]{} |
| lexmax { [j1,j2] -> [i1,i2,i3,i4,i5,i6,i7,i8,i9,i10] : 1 <= i1,j1 <= 8 and 1 <= i2,i3,i4,i5,i6,i7,i8,i9,i10 <= 2 and 1 <= j2 <= 128 and i1-1 = j1-1 and i2-1+2*i3-2+4*i4-4+8*i5-8+16*i6-16+32*i7-32+64*i8-64+128*i9-128+256*i10-256=3*j2-3+66 }; |
| \end{lstlisting} |
| This problem was the main inspiration |
| for some of the optimizations in \autoref{s:GBR}. |
| The second group of test cases are projections used during counting. |
| The first nine of these come from \textcite{Seghir2006minimizing}. |
| The remaining two come from \textcite{Verdoolaege2005experiences} and |
| were used to drive the first, Gomory cuts based, implementation |
| in {\tt isl}. |
| The third and final group of test cases are borrowed from |
| \textcite{Bygde2010licentiate} and inspired the offline symmetry detection |
| of \autoref{s:offline}. Without symmetry detection, the running times |
| are 11s and 5.9s. |
| All running times of {\tt barvinok} and {\tt isl} include a conversion |
| to disjunctive normal form. Without this conversion, the final two |
| cases can be solved in 0.07s and 0.21s. |
| The {\tt PipLib} implementation has some fixed limits and will |
| sometimes report the problem to be too complex (TC), while on some other |
| problems it will run out of memory (OOM). |
| The {\tt barvinok} implementation does not support problems |
| with a non-trivial lineality space (line) nor maximization problems (max). |
| The Gomory cuts based {\tt isl} implementation was terminated after 1000 |
| minutes on the first problem. The gbr version introduces some |
| overhead on some of the easier problems, but is overall the clear winner. |
| |
| \begin{table} |
| \begin{center} |
| \begin{tabular}{lrrrrr} |
| & {\tt PipLib} & {\tt barvinok} & {\tt isl} cut & {\tt isl} gbr & {\tt PPL} \\ |
| \hline |
| \hline |
| % bart.pip |
| Phideo & TC & 793m & $>$999m & 2.7s & 372m \\ |
| \hline |
| e1 & 0.33s & 3.5s & 0.08s & 0.11s & 0.18s \\ |
| e3 & 0.14s & 0.13s & 0.10s & 0.10s & 0.17s \\ |
| e4 & 0.24s & 9.1s & 0.09s & 0.11s & 0.70s \\ |
| e5 & 0.12s & 6.0s & 0.06s & 0.14s & 0.17s \\ |
| e6 & 0.10s & 6.8s & 0.17s & 0.08s & 0.21s \\ |
| e7 & 0.03s & 0.27s & 0.04s & 0.04s & 0.03s \\ |
| e8 & 0.03s & 0.18s & 0.03s & 0.04s & 0.01s \\ |
| e9 & OOM & 70m & 2.6s & 0.94s & 22s \\ |
| vd & 0.04s & 0.10s & 0.03s & 0.03s & 0.03s \\ |
| bouleti & 0.25s & line & 0.06s & 0.06s & 0.15s \\ |
| difficult & OOM & 1.3s & 1.7s & 0.33s & 1.4s \\ |
| \hline |
| cnt/sum & TC & max & 2.2s & 2.2s & OOM \\ |
| jcomplex & TC & max & 3.7s & 3.9s & OOM \\ |
| \end{tabular} |
| \caption{Comparison of Execution Times} |
| \label{t:comparison} |
| \end{center} |
| \end{table} |
| |
| \subsection{Online Symmetry Detection}\label{s:online} |
| |
| Manual experiments on small instances of the problems of |
| \textcite{Bygde2010licentiate} and an analysis of the results |
| by the approximate MPA method developed by \textcite{Bygde2010licentiate} |
| have revealed that these problems contain many more symmetries |
| than can be detected using the offline method of \autoref{s:offline}. |
| In this section, we present an online detection mechanism that has |
| not been implemented yet, but that has shown promising results |
| in manual applications. |
| |
| Let us first consider what happens when we do not perform offline |
| symmetry detection. At some point, one of the |
| $b_i(\vec p) + \sp {\vec a} {\vec x} \ge 0$ constraints, |
| say the $j$th constraint, appears as a column |
| variable, say $c_1$, while the other constraints are represented |
| as rows of the form $b_i(\vec p) - b_j(\vec p) + c$. |
| The context is then split according to the relative order of |
| $b_j(\vec p)$ and one of the remaining $b_i(\vec p)$. |
| The offline method avoids this split by replacing all $b_i(\vec p)$ |
| by a single newly introduced parameter that represents the minimum |
| of these $b_i(\vec p)$. |
| In the online method the split is similarly avoided by the introduction |
| of a new parameter. In particular, a new parameter is introduced |
| that represents |
| $\left| b_j(\vec p) - b_i(\vec p) \right|_+ = |
| \max(b_j(\vec p) - b_i(\vec p), 0)$. |
| |
| In general, let $r = b(\vec p) + \sp {\vec a} {\vec c}$ be a row |
| of the tableau such that the sign of $b(\vec p)$ is indeterminate |
| and such that exactly one of the elements of $\vec a$ is a $1$, |
| while all remaining elements are non-positive. |
| That is, $r = b(\vec p) + c_j - f$ with $f = -\sum_{i\ne j} a_i c_i \ge 0$. |
| We introduce a new parameter $t$ with |
| context constraints $t \ge -b(\vec p)$ and $t \ge 0$ and replace |
| the column variable $c_j$ by $c' + t$. The row $r$ is now equal |
| to $b(\vec p) + t + c' - f$. The constant term of this row is always |
| non-negative because any negative value of $b(\vec p)$ is compensated |
| by $t \ge -b(\vec p)$ while and non-negative value remains non-negative |
| because $t \ge 0$. |
| |
| We need to show that this transformation does not eliminate any valid |
| solutions and that it does not introduce any spurious solutions. |
| Given a valid solution for the original problem, we need to find |
| a non-negative value of $c'$ satisfying the constraints. |
| If $b(\vec p) \ge 0$, we can take $t = 0$ so that |
| $c' = c_j - t = c_j \ge 0$. |
| If $b(\vec p) < 0$, we can take $t = -b(\vec p)$. |
| Since $r = b(\vec p) + c_j - f \ge 0$ and $f \ge 0$, we have |
| $c' = c_j + b(\vec p) \ge 0$. |
| Note that these choices amount to plugging in |
| $t = \left|-b(\vec p)\right|_+ = \max(-b(\vec p), 0)$. |
| Conversely, given a solution to the new problem, we need to find |
| a non-negative value of $c_j$, but this is easy since $c_j = c' + t$ |
| and both of these are non-negative. |
| |
| Plugging in $t = \max(-b(\vec p), 0)$ can be performed as in |
| \autoref{s:post}, but, as in the case of offline symmetry detection, |
| it may be better to provide a direct representation for such |
| expressions in the internal representation of sets and relations |
| or at least in a quast-like output format. |
| |
| \section{Coalescing}\label{s:coalescing} |
| |
| See \textcite{Verdoolaege2015impact} for details on integer set coalescing. |
| |
| \section{Transitive Closure} |
| |
| \subsection{Introduction} |
| |
| \begin{definition}[Power of a Relation] |
| Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation and |
| $k \in \Z_{\ge 1}$ |
| a positive number, then power $k$ of relation $R$ is defined as |
| \begin{equation} |
| \label{eq:transitive:power} |
| R^k \coloneqq |
| \begin{cases} |
| R & \text{if $k = 1$} |
| \\ |
| R \circ R^{k-1} & \text{if $k \ge 2$} |
| . |
| \end{cases} |
| \end{equation} |
| \end{definition} |
| |
| \begin{definition}[Transitive Closure of a Relation] |
| Let $R \in \Z^n \to 2^{\Z^{d+d}}$ be a relation, |
| then the transitive closure $R^+$ of $R$ is the union |
| of all positive powers of $R$, |
| $$ |
| R^+ \coloneqq \bigcup_{k \ge 1} R^k |
| . |
| $$ |
| \end{definition} |
| Alternatively, the transitive closure may be defined |
| inductively as |
| \begin{equation} |
| \label{eq:transitive:inductive} |
| R^+ \coloneqq R \cup \left(R \circ R^+\right) |
| . |
| \end{equation} |
| |
| Since the transitive closure of a polyhedral relation |
| may no longer be a polyhedral relation \parencite{Kelly1996closure}, |
| we can, in the general case, only compute an approximation |
| of the transitive closure. |
| Whereas \textcite{Kelly1996closure} compute underapproximations, |
| we, like \textcite{Beletska2009}, compute overapproximations. |
| That is, given a relation $R$, we will compute a relation $T$ |
| such that $R^+ \subseteq T$. Of course, we want this approximation |
| to be as close as possible to the actual transitive closure |
| $R^+$ and we want to detect the cases where the approximation is |
| exact, i.e., where $T = R^+$. |
| |
| For computing an approximation of the transitive closure of $R$, |
| we follow the same general strategy as \textcite{Beletska2009} |
| and first compute an approximation of $R^k$ for $k \ge 1$ and then project |
| out the parameter $k$ from the resulting relation. |
| |
| \begin{example} |
| As a trivial example, consider the relation |
| $R = \{\, x \to x + 1 \,\}$. The $k$th power of this map |
| for arbitrary $k$ is |
| $$ |
| R^k = k \mapsto \{\, x \to x + k \mid k \ge 1 \,\} |
| . |
| $$ |
| The transitive closure is then |
| $$ |
| \begin{aligned} |
| R^+ & = \{\, x \to y \mid \exists k \in \Z_{\ge 1} : y = x + k \,\} |
| \\ |
| & = \{\, x \to y \mid y \ge x + 1 \,\} |
| . |
| \end{aligned} |
| $$ |
| \end{example} |
| |
| \subsection{Computing an Approximation of $R^k$} |
| \label{s:power} |
| |
| There are some special cases where the computation of $R^k$ is very easy. |
| One such case is that where $R$ does not compose with itself, |
| i.e., $R \circ R = \emptyset$ or $\domain R \cap \range R = \emptyset$. |
| In this case, $R^k$ is only non-empty for $k=1$ where it is equal |
| to $R$ itself. |
| |
| In general, it is impossible to construct a closed form |
| of $R^k$ as a polyhedral relation. |
| We will therefore need to make some approximations. |
| As a first approximations, we will consider each of the basic |
| relations in $R$ as simply adding one or more offsets to a domain element |
| to arrive at an image element and ignore the fact that some of these |
| offsets may only be applied to some of the domain elements. |
| That is, we will only consider the difference set $\Delta\,R$ of the relation. |
| In particular, we will first construct a collection $P$ of paths |
| that move through |
| a total of $k$ offsets and then intersect domain and range of this |
| collection with those of $R$. |
| That is, |
| \begin{equation} |
| \label{eq:transitive:approx} |
| K = P \cap \left(\domain R \to \range R\right) |
| , |
| \end{equation} |
| with |
| \begin{equation} |
| \label{eq:transitive:path} |
| P = \vec s \mapsto \{\, \vec x \to \vec y \mid |
| \exists k_i \in \Z_{\ge 0}, \vec\delta_i \in k_i \, \Delta_i(\vec s) : |
| \vec y = \vec x + \sum_i \vec\delta_i |
| \wedge |
| \sum_i k_i = k > 0 |
| \,\} |
| \end{equation} |
| and with $\Delta_i$ the basic sets that compose |
| the difference set $\Delta\,R$. |
| Note that the number of basic sets $\Delta_i$ need not be |
| the same as the number of basic relations in $R$. |
| Also note that since addition is commutative, it does not |
| matter in which order we add the offsets and so we are allowed |
| to group them as we did in \eqref{eq:transitive:path}. |
| |
| If all the $\Delta_i$s are singleton sets |
| $\Delta_i = \{\, \vec \delta_i \,\}$ with $\vec \delta_i \in \Z^d$, |
| then \eqref{eq:transitive:path} simplifies to |
| \begin{equation} |
| \label{eq:transitive:singleton} |
| P = \{\, \vec x \to \vec y \mid |
| \exists k_i \in \Z_{\ge 0} : |
| \vec y = \vec x + \sum_i k_i \, \vec \delta_i |
| \wedge |
| \sum_i k_i = k > 0 |
| \,\} |
| \end{equation} |
| and then the approximation computed in \eqref{eq:transitive:approx} |
| is essentially the same as that of \textcite{Beletska2009}. |
| If some of the $\Delta_i$s are not singleton sets or if |
| some of $\vec \delta_i$s are parametric, then we need |
| to resort to further approximations. |
| |
| To ease both the exposition and the implementation, we will for |
| the remainder of this section work with extended offsets |
| $\Delta_i' = \Delta_i \times \{\, 1 \,\}$. |
| That is, each offset is extended with an extra coordinate that is |
| set equal to one. The paths constructed by summing such extended |
| offsets have the length encoded as the difference of their |
| final coordinates. The path $P'$ can then be decomposed into |
| paths $P_i'$, one for each $\Delta_i$, |
| \begin{equation} |
| \label{eq:transitive:decompose} |
| P' = \left( |
| (P_m' \cup \identity) \circ \cdots \circ |
| (P_2' \cup \identity) \circ |
| (P_1' \cup \identity) |
| \right) \cap |
| \{\, |
| \vec x' \to \vec y' \mid y_{d+1} - x_{d+1} = k > 0 |
| \,\} |
| , |
| \end{equation} |
| with |
| $$ |
| P_i' = \vec s \mapsto \{\, \vec x' \to \vec y' \mid |
| \exists k \in \Z_{\ge 1}, \vec \delta \in k \, \Delta_i'(\vec s) : |
| \vec y' = \vec x' + \vec \delta |
| \,\} |
| . |
| $$ |
| Note that each $P_i'$ contains paths of length at least one. |
| We therefore need to take the union with the identity relation |
| when composing the $P_i'$s to allow for paths that do not contain |
| any offsets from one or more $\Delta_i'$. |
| The path that consists of only identity relations is removed |
| by imposing the constraint $y_{d+1} - x_{d+1} > 0$. |
| Taking the union with the identity relation means that |
| that the relations we compose in \eqref{eq:transitive:decompose} |
| each consist of two basic relations. If there are $m$ |
| disjuncts in the input relation, then a direct application |
| of the composition operation may therefore result in a relation |
| with $2^m$ disjuncts, which is prohibitively expensive. |
| It is therefore crucial to apply coalescing (\autoref{s:coalescing}) |
| after each composition. |
| |
| Let us now consider how to compute an overapproximation of $P_i'$. |
| Those that correspond to singleton $\Delta_i$s are grouped together |
| and handled as in \eqref{eq:transitive:singleton}. |
| Note that this is just an optimization. The procedure described |
| below would produce results that are at least as accurate. |
| For simplicity, we first assume that no constraint in $\Delta_i'$ |
| involves any existentially quantified variables. |
| We will return to existentially quantified variables at the end |
| of this section. |
| Without existentially quantified variables, we can classify |
| the constraints of $\Delta_i'$ as follows |
| \begin{enumerate} |
| \item non-parametric constraints |
| \begin{equation} |
| \label{eq:transitive:non-parametric} |
| A_1 \vec x + \vec c_1 \geq \vec 0 |
| \end{equation} |
| \item purely parametric constraints |
| \begin{equation} |
| \label{eq:transitive:parametric} |
| B_2 \vec s + \vec c_2 \geq \vec 0 |
| \end{equation} |
| \item negative mixed constraints |
| \begin{equation} |
| \label{eq:transitive:mixed} |
| A_3 \vec x + B_3 \vec s + \vec c_3 \geq \vec 0 |
| \end{equation} |
| such that for each row $j$ and for all $\vec s$, |
| $$ |
| \Delta_i'(\vec s) \cap |
| \{\, \vec \delta' \mid B_{3,j} \vec s + c_{3,j} > 0 \,\} |
| = \emptyset |
| $$ |
| \item positive mixed constraints |
| $$ |
| A_4 \vec x + B_4 \vec s + \vec c_4 \geq \vec 0 |
| $$ |
| such that for each row $j$, there is at least one $\vec s$ such that |
| $$ |
| \Delta_i'(\vec s) \cap |
| \{\, \vec \delta' \mid B_{4,j} \vec s + c_{4,j} > 0 \,\} |
| \ne \emptyset |
| $$ |
| \end{enumerate} |
| We will use the following approximation $Q_i$ for $P_i'$: |
| \begin{equation} |
| \label{eq:transitive:Q} |
| \begin{aligned} |
| Q_i = \vec s \mapsto |
| \{\, |
| \vec x' \to \vec y' |
| \mid {} & \exists k \in \Z_{\ge 1}, \vec f \in \Z^d : |
| \vec y' = \vec x' + (\vec f, k) |
| \wedge {} |
| \\ |
| & |
| A_1 \vec f + k \vec c_1 \geq \vec 0 |
| \wedge |
| B_2 \vec s + \vec c_2 \geq \vec 0 |
| \wedge |
| A_3 \vec f + B_3 \vec s + \vec c_3 \geq \vec 0 |
| \,\} |
| . |
| \end{aligned} |
| \end{equation} |
| To prove that $Q_i$ is indeed an overapproximation of $P_i'$, |
| we need to show that for every $\vec s \in \Z^n$, for every |
| $k \in \Z_{\ge 1}$ and for every $\vec f \in k \, \Delta_i(\vec s)$ |
| we have that |
| $(\vec f, k)$ satisfies the constraints in \eqref{eq:transitive:Q}. |
| If $\Delta_i(\vec s)$ is non-empty, then $\vec s$ must satisfy |
| the constraints in \eqref{eq:transitive:parametric}. |
| Each element $(\vec f, k) \in k \, \Delta_i'(\vec s)$ is a sum |
| of $k$ elements $(\vec f_j, 1)$ in $\Delta_i'(\vec s)$. |
| Each of these elements satisfies the constraints in |
| \eqref{eq:transitive:non-parametric}, i.e., |
| $$ |
| \left[ |
| \begin{matrix} |
| A_1 & \vec c_1 |
| \end{matrix} |
| \right] |
| \left[ |
| \begin{matrix} |
| \vec f_j \\ 1 |
| \end{matrix} |
| \right] |
| \ge \vec 0 |
| . |
| $$ |
| The sum of these elements therefore satisfies the same set of inequalities, |
| i.e., $A_1 \vec f + k \vec c_1 \geq \vec 0$. |
| Finally, the constraints in \eqref{eq:transitive:mixed} are such |
| that for any $\vec s$ in the parameter domain of $\Delta$, |
| we have $-\vec r(\vec s) \coloneqq B_3 \vec s + \vec c_3 \le \vec 0$, |
| i.e., $A_3 \vec f_j \ge \vec r(\vec s) \ge \vec 0$ |
| and therefore also $A_3 \vec f \ge \vec r(\vec s)$. |
| Note that if there are no mixed constraints and if the |
| rational relaxation of $\Delta_i(\vec s)$, i.e., |
| $\{\, \vec x \in \Q^d \mid A_1 \vec x + \vec c_1 \ge \vec 0\,\}$, |
| has integer vertices, then the approximation is exact, i.e., |
| $Q_i = P_i'$. In this case, the vertices of $\Delta'_i(\vec s)$ |
| generate the rational cone |
| $\{\, \vec x' \in \Q^{d+1} \mid \left[ |
| \begin{matrix} |
| A_1 & \vec c_1 |
| \end{matrix} |
| \right] \vec x' \,\}$ and therefore $\Delta'_i(\vec s)$ is |
| a Hilbert basis of this cone \parencite[Theorem~16.4]{Schrijver1986}. |
| |
| Note however that, as pointed out by \textcite{DeSmet2010personal}, |
| if there \emph{are} any mixed constraints, then the above procedure may |
| not compute the most accurate affine approximation of |
| $k \, \Delta_i(\vec s)$ with $k \ge 1$. |
| In particular, we only consider the negative mixed constraints that |
| happen to appear in the description of $\Delta_i(\vec s)$, while we |
| should instead consider \emph{all} valid such constraints. |
| It is also sufficient to consider those constraints because any |
| constraint that is valid for $k \, \Delta_i(\vec s)$ is also |
| valid for $1 \, \Delta_i(\vec s) = \Delta_i(\vec s)$. |
| Take therefore any constraint |
| $\spv a x + \spv b s + c \ge 0$ valid for $\Delta_i(\vec s)$. |
| This constraint is also valid for $k \, \Delta_i(\vec s)$ iff |
| $k \, \spv a x + \spv b s + c \ge 0$. |
| If $\spv b s + c$ can attain any positive value, then $\spv a x$ |
| may be negative for some elements of $\Delta_i(\vec s)$. |
| We then have $k \, \spv a x < \spv a x$ for $k > 1$ and so the constraint |
| is not valid for $k \, \Delta_i(\vec s)$. |
| We therefore need to impose $\spv b s + c \le 0$ for all values |
| of $\vec s$ such that $\Delta_i(\vec s)$ is non-empty, i.e., |
| $\vec b$ and $c$ need to be such that $- \spv b s - c \ge 0$ is a valid |
| constraint of $\Delta_i(\vec s)$. That is, $(\vec b, c)$ are the opposites |
| of the coefficients of a valid constraint of $\Delta_i(\vec s)$. |
| The approximation of $k \, \Delta_i(\vec s)$ can therefore be obtained |
| using three applications of Farkas' lemma. The first obtains the coefficients |
| of constraints valid for $\Delta_i(\vec s)$. The second obtains |
| the coefficients of constraints valid for the projection of $\Delta_i(\vec s)$ |
| onto the parameters. The opposite of the second set is then computed |
| and intersected with the first set. The result is the set of coefficients |
| of constraints valid for $k \, \Delta_i(\vec s)$. A final application |
| of Farkas' lemma is needed to obtain the approximation of |
| $k \, \Delta_i(\vec s)$ itself. |
| |
| \begin{example} |
| Consider the relation |
| $$ |
| n \to \{\, (x, y) \to (1 + x, 1 - n + y) \mid n \ge 2 \,\} |
| . |
| $$ |
| The difference set of this relation is |
| $$ |
| \Delta = n \to \{\, (1, 1 - n) \mid n \ge 2 \,\} |
| . |
| $$ |
| Using our approach, we would only consider the mixed constraint |
| $y - 1 + n \ge 0$, leading to the following approximation of the |
| transitive closure: |
| $$ |
| n \to \{\, (x, y) \to (o_0, o_1) \mid n \ge 2 \wedge o_1 \le 1 - n + y \wedge o_0 \ge 1 + x \,\} |
| . |
| $$ |
| If, instead, we apply Farkas's lemma to $\Delta$, i.e., |
| \begin{verbatim} |
| D := [n] -> { [1, 1 - n] : n >= 2 }; |
| CD := coefficients D; |
| CD; |
| \end{verbatim} |
| we obtain |
| \begin{verbatim} |
| { rat: coefficients[[c_cst, c_n] -> [i2, i3]] : i3 <= c_n and |
| i3 <= c_cst + 2c_n + i2 } |
| \end{verbatim} |
| The pure-parametric constraints valid for $\Delta$, |
| \begin{verbatim} |
| P := { [a,b] -> [] }(D); |
| CP := coefficients P; |
| CP; |
| \end{verbatim} |
| are |
| \begin{verbatim} |
| { rat: coefficients[[c_cst, c_n] -> []] : c_n >= 0 and 2c_n >= -c_cst } |
| \end{verbatim} |
| Negating these coefficients and intersecting with \verb+CD+, |
| \begin{verbatim} |
| NCP := { rat: coefficients[[a,b] -> []] |
| -> coefficients[[-a,-b] -> []] }(CP); |
| CK := wrap((unwrap CD) * (dom (unwrap NCP))); |
| CK; |
| \end{verbatim} |
| we obtain |
| \begin{verbatim} |
| { rat: [[c_cst, c_n] -> [i2, i3]] : i3 <= c_n and |
| i3 <= c_cst + 2c_n + i2 and c_n <= 0 and 2c_n <= -c_cst } |
| \end{verbatim} |
| The approximation for $k\,\Delta$, |
| \begin{verbatim} |
| K := solutions CK; |
| K; |
| \end{verbatim} |
| is then |
| \begin{verbatim} |
| [n] -> { rat: [i0, i1] : i1 <= -i0 and i0 >= 1 and i1 <= 2 - n - i0 } |
| \end{verbatim} |
| Finally, the computed approximation for $R^+$, |
| \begin{verbatim} |
| T := unwrap({ [dx,dy] -> [[x,y] -> [x+dx,y+dy]] }(K)); |
| R := [n] -> { [x,y] -> [x+1,y+1-n] : n >= 2 }; |
| T := T * ((dom R) -> (ran R)); |
| T; |
| \end{verbatim} |
| is |
| \begin{verbatim} |
| [n] -> { [x, y] -> [o0, o1] : o1 <= x + y - o0 and |
| o0 >= 1 + x and o1 <= 2 - n + x + y - o0 and n >= 2 } |
| \end{verbatim} |
| \end{example} |
| |
| Existentially quantified variables can be handled by |
| classifying them into variables that are uniquely |
| determined by the parameters, variables that are independent |
| of the parameters and others. The first set can be treated |
| as parameters and the second as variables. Constraints involving |
| the other existentially quantified variables are removed. |
| |
| \begin{example} |
| Consider the relation |
| $$ |
| R = |
| n \to \{\, x \to y \mid \exists \, \alpha_0, \alpha_1: 7\alpha_0 = -2 + n \wedge 5\alpha_1 = -1 - x + y \wedge y \ge 6 + x \,\} |
| . |
| $$ |
| The difference set of this relation is |
| $$ |
| \Delta = \Delta \, R = |
| n \to \{\, x \mid \exists \, \alpha_0, \alpha_1: 7\alpha_0 = -2 + n \wedge 5\alpha_1 = -1 + x \wedge x \ge 6 \,\} |
| . |
| $$ |
| The existentially quantified variables can be defined in terms |
| of the parameters and variables as |
| $$ |
| \alpha_0 = \floor{\frac{-2 + n}7} |
| \qquad |
| \text{and} |
| \qquad |
| \alpha_1 = \floor{\frac{-1 + x}5} |
| . |
| $$ |
| $\alpha_0$ can therefore be treated as a parameter, |
| while $\alpha_1$ can be treated as a variable. |
| This in turn means that $7\alpha_0 = -2 + n$ can be treated as |
| a purely parametric constraint, while the other two constraints are |
| non-parametric. |
| The corresponding $Q$~\eqref{eq:transitive:Q} is therefore |
| $$ |
| \begin{aligned} |
| n \to \{\, (x,z) \to (y,w) \mid |
| \exists\, \alpha_0, \alpha_1, k, f : {} & |
| k \ge 1 \wedge |
| y = x + f \wedge |
| w = z + k \wedge {} \\ |
| & |
| 7\alpha_0 = -2 + n \wedge |
| 5\alpha_1 = -k + x \wedge |
| x \ge 6 k |
| \,\} |
| . |
| \end{aligned} |
| $$ |
| Projecting out the final coordinates encoding the length of the paths, |
| results in the exact transitive closure |
| $$ |
| R^+ = |
| n \to \{\, x \to y \mid \exists \, \alpha_0, \alpha_1: 7\alpha_1 = -2 + n \wedge 6\alpha_0 \ge -x + y \wedge 5\alpha_0 \le -1 - x + y \,\} |
| . |
| $$ |
| \end{example} |
| |
| The fact that we ignore some impure constraints clearly leads |
| to a loss of accuracy. In some cases, some of this loss can be recovered |
| by not considering the parameters in a special way. |
| That is, instead of considering the set |
| $$ |
| \Delta = \diff R = |
| \vec s \mapsto |
| \{\, \vec \delta \in \Z^{d} \mid \exists \vec x \to \vec y \in R : |
| \vec \delta = \vec y - \vec x |
| \,\} |
| $$ |
| we consider the set |
| $$ |
| \Delta' = \diff R' = |
| \{\, \vec \delta \in \Z^{n+d} \mid \exists |
| (\vec s, \vec x) \to (\vec s, \vec y) \in R' : |
| \vec \delta = (\vec s - \vec s, \vec y - \vec x) |
| \,\} |
| . |
| $$ |
| The first $n$ coordinates of every element in $\Delta'$ are zero. |
| Projecting out these zero coordinates from $\Delta'$ is equivalent |
| to projecting out the parameters in $\Delta$. |
| The result is obviously a superset of $\Delta$, but all its constraints |
| are of type \eqref{eq:transitive:non-parametric} and they can therefore |
| all be used in the construction of $Q_i$. |
| |
| \begin{example} |
| Consider the relation |
| $$ |
| % [n] -> { [x, y] -> [1 + x, 1 - n + y] | n >= 2 } |
| R = n \to \{\, (x, y) \to (1 + x, 1 - n + y) \mid n \ge 2 \,\} |
| . |
| $$ |
| We have |
| $$ |
| \diff R = n \to \{\, (1, 1 - n) \mid n \ge 2 \,\} |
| $$ |
| and so, by treating the parameters in a special way, we obtain |
| the following approximation for $R^+$: |
| $$ |
| n \to \{\, (x, y) \to (x', y') \mid n \ge 2 \wedge y' \le 1 - n + y \wedge x' \ge 1 + x \,\} |
| . |
| $$ |
| If we consider instead |
| $$ |
| R' = \{\, (n, x, y) \to (n, 1 + x, 1 - n + y) \mid n \ge 2 \,\} |
| $$ |
| then |
| $$ |
| \diff R' = \{\, (0, 1, y) \mid y \le -1 \,\} |
| $$ |
| and we obtain the approximation |
| $$ |
| n \to \{\, (x, y) \to (x', y') \mid n \ge 2 \wedge x' \ge 1 + x \wedge y' \le x + y - x' \,\} |
| . |
| $$ |
| If we consider both $\diff R$ and $\diff R'$, then we obtain |
| $$ |
| n \to \{\, (x, y) \to (x', y') \mid n \ge 2 \wedge y' \le 1 - n + y \wedge x' \ge 1 + x \wedge y' \le x + y - x' \,\} |
| . |
| $$ |
| Note, however, that this is not the most accurate affine approximation that |
| can be obtained. That would be |
| $$ |
| n \to \{\, (x, y) \to (x', y') \mid y' \le 2 - n + x + y - x' \wedge n \ge 2 \wedge x' \ge 1 + x \,\} |
| . |
| $$ |
| \end{example} |
| |
| \subsection{Checking Exactness} |
| |
| The approximation $T$ for the transitive closure $R^+$ can be obtained |
| by projecting out the parameter $k$ from the approximation $K$ |
| \eqref{eq:transitive:approx} of the power $R^k$. |
| Since $K$ is an overapproximation of $R^k$, $T$ will also be an |
| overapproximation of $R^+$. |
| To check whether the results are exact, we need to consider two |
| cases depending on whether $R$ is {\em cyclic}, where $R$ is defined |
| to be cyclic if $R^+$ maps any element to itself, i.e., |
| $R^+ \cap \identity \ne \emptyset$. |
| If $R$ is acyclic, then the inductive definition of |
| \eqref{eq:transitive:inductive} is equivalent to its completion, |
| i.e., |
| $$ |
| R^+ = R \cup \left(R \circ R^+\right) |
| $$ |
| is a defining property. |
| Since $T$ is known to be an overapproximation, we only need to check |
| whether |
| $$ |
| T \subseteq R \cup \left(R \circ T\right) |
| . |
| $$ |
| This is essentially Theorem~5 of \textcite{Kelly1996closure}. |
| The only difference is that they only consider lexicographically |
| forward relations, a special case of acyclic relations. |
| |
| If, on the other hand, $R$ is cyclic, then we have to resort |
| to checking whether the approximation $K$ of the power is exact. |
| Note that $T$ may be exact even if $K$ is not exact, so the check |
| is sound, but incomplete. |
| To check exactness of the power, we simply need to check |
| \eqref{eq:transitive:power}. Since again $K$ is known |
| to be an overapproximation, we only need to check whether |
| $$ |
| \begin{aligned} |
| K'|_{y_{d+1} - x_{d+1} = 1} & \subseteq R' |
| \\ |
| K'|_{y_{d+1} - x_{d+1} \ge 2} & \subseteq R' \circ K'|_{y_{d+1} - x_{d+1} \ge 1} |
| , |
| \end{aligned} |
| $$ |
| where $R' = \{\, \vec x' \to \vec y' \mid \vec x \to \vec y \in R |
| \wedge y_{d+1} - x_{d+1} = 1\,\}$, i.e., $R$ extended with path |
| lengths equal to 1. |
| |
| All that remains is to explain how to check the cyclicity of $R$. |
| Note that the exactness on the power is always sound, even |
| in the acyclic case, so we only need to be careful that we find |
| all cyclic cases. Now, if $R$ is cyclic, i.e., |
| $R^+ \cap \identity \ne \emptyset$, then, since $T$ is |
| an overapproximation of $R^+$, also |
| $T \cap \identity \ne \emptyset$. This in turn means |
| that $\Delta \, K'$ contains a point whose first $d$ coordinates |
| are zero and whose final coordinate is positive. |
| In the implementation we currently perform this test on $P'$ instead of $K'$. |
| Note that if $R^+$ is acyclic and $T$ is not, then the approximation |
| is clearly not exact and the approximation of the power $K$ |
| will not be exact either. |
| |
| \subsection{Decomposing $R$ into strongly connected components} |
| |
| If the input relation $R$ is a union of several basic relations |
| that can be partially ordered |
| then the accuracy of the approximation may be improved by computing |
| an approximation of each strongly connected components separately. |
| For example, if $R = R_1 \cup R_2$ and $R_1 \circ R_2 = \emptyset$, |
| then we know that any path that passes through $R_2$ cannot later |
| pass through $R_1$, i.e., |
| \begin{equation} |
| \label{eq:transitive:components} |
| R^+ = R_1^+ \cup R_2^+ \cup \left(R_2^+ \circ R_1^+\right) |
| . |
| \end{equation} |
| We can therefore compute (approximations of) transitive closures |
| of $R_1$ and $R_2$ separately. |
| Note, however, that the condition $R_1 \circ R_2 = \emptyset$ |
| is actually too strong. |
| If $R_1 \circ R_2$ is a subset of $R_2 \circ R_1$ |
| then we can reorder the segments |
| in any path that moves through both $R_1$ and $R_2$ to |
| first move through $R_1$ and then through $R_2$. |
| |
| This idea can be generalized to relations that are unions |
| of more than two basic relations by constructing the |
| strongly connected components in the graph with as vertices |
| the basic relations and an edge between two basic relations |
| $R_i$ and $R_j$ if $R_i$ needs to follow $R_j$ in some paths. |
| That is, there is an edge from $R_i$ to $R_j$ iff |
| \begin{equation} |
| \label{eq:transitive:edge} |
| R_i \circ R_j |
| \not\subseteq |
| R_j \circ R_i |
| . |
| \end{equation} |
| The components can be obtained from the graph by applying |
| Tarjan's algorithm \parencite{Tarjan1972}. |
| |
| In practice, we compute the (extended) powers $K_i'$ of each component |
| separately and then compose them as in \eqref{eq:transitive:decompose}. |
| Note, however, that in this case the order in which we apply them is |
| important and should correspond to a topological ordering of the |
| strongly connected components. Simply applying Tarjan's |
| algorithm will produce topologically sorted strongly connected components. |
| The graph on which Tarjan's algorithm is applied is constructed on-the-fly. |
| That is, whenever the algorithm checks if there is an edge between |
| two vertices, we evaluate \eqref{eq:transitive:edge}. |
| The exactness check is performed on each component separately. |
| If the approximation turns out to be inexact for any of the components, |
| then the entire result is marked inexact and the exactness check |
| is skipped on the components that still need to be handled. |
| |
| It should be noted that \eqref{eq:transitive:components} |
| is only valid for exact transitive closures. |
| If overapproximations are computed in the right hand side, then the result will |
| still be an overapproximation of the left hand side, but this result |
| may not be transitively closed. If we only separate components based |
| on the condition $R_i \circ R_j = \emptyset$, then there is no problem, |
| as this condition will still hold on the computed approximations |
| of the transitive closures. If, however, we have exploited |
| \eqref{eq:transitive:edge} during the decomposition and if the |
| result turns out not to be exact, then we check whether |
| the result is transitively closed. If not, we recompute |
| the transitive closure, skipping the decomposition. |
| Note that testing for transitive closedness on the result may |
| be fairly expensive, so we may want to make this check |
| configurable. |
| |
| \begin{figure} |
| \begin{center} |
| \begin{tikzpicture}[x=0.5cm,y=0.5cm,>=stealth,shorten >=1pt] |
| \foreach \x in {1,...,10}{ |
| \foreach \y in {1,...,10}{ |
| \draw[->] (\x,\y) -- (\x,\y+1); |
| } |
| } |
| \foreach \x in {1,...,20}{ |
| \foreach \y in {5,...,15}{ |
| \draw[->] (\x,\y) -- (\x+1,\y); |
| } |
| } |
| \end{tikzpicture} |
| \end{center} |
| \caption{The relation from \autoref{ex:closure4}} |
| \label{f:closure4} |
| \end{figure} |
| \begin{example} |
| \label{ex:closure4} |
| Consider the relation in example {\tt closure4} that comes with |
| the Omega calculator~\parencite{Omega_calc}, $R = R_1 \cup R_2$, |
| with |
| $$ |
| \begin{aligned} |
| R_1 & = \{\, (x,y) \to (x,y+1) \mid 1 \le x,y \le 10 \,\} |
| \\ |
| R_2 & = \{\, (x,y) \to (x+1,y) \mid 1 \le x \le 20 \wedge 5 \le y \le 15 \,\} |
| . |
| \end{aligned} |
| $$ |
| This relation is shown graphically in \autoref{f:closure4}. |
| We have |
| $$ |
| \begin{aligned} |
| R_1 \circ R_2 &= |
| \{\, (x,y) \to (x+1,y+1) \mid 1 \le x \le 9 \wedge 5 \le y \le 10 \,\} |
| \\ |
| R_2 \circ R_1 &= |
| \{\, (x,y) \to (x+1,y+1) \mid 1 \le x \le 10 \wedge 4 \le y \le 10 \,\} |
| . |
| \end{aligned} |
| $$ |
| Clearly, $R_1 \circ R_2 \subseteq R_2 \circ R_1$ and so |
| $$ |
| \left( |
| R_1 \cup R_2 |
| \right)^+ |
| = |
| \left(R_2^+ \circ R_1^+\right) |
| \cup R_1^+ |
| \cup R_2^+ |
| . |
| $$ |
| \end{example} |
| |
| \begin{figure} |
| \newcounter{n} |
| \newcounter{t1} |
| \newcounter{t2} |
| \newcounter{t3} |
| \newcounter{t4} |
| \begin{center} |
| \begin{tikzpicture}[>=stealth,shorten >=1pt] |
| \setcounter{n}{7} |
| \foreach \i in {1,...,\value{n}}{ |
| \foreach \j in {1,...,\value{n}}{ |
| \setcounter{t1}{2 * \j - 4 - \i + 1} |
| \setcounter{t2}{\value{n} - 3 - \i + 1} |
| \setcounter{t3}{2 * \i - 1 - \j + 1} |
| \setcounter{t4}{\value{n} - \j + 1} |
| \ifnum\value{t1}>0\ifnum\value{t2}>0 |
| \ifnum\value{t3}>0\ifnum\value{t4}>0 |
| \draw[thick,->] (\i,\j) to[out=20] (\i+3,\j); |
| \fi\fi\fi\fi |
| \setcounter{t1}{2 * \j - 1 - \i + 1} |
| \setcounter{t2}{\value{n} - \i + 1} |
| \setcounter{t3}{2 * \i - 4 - \j + 1} |
| \setcounter{t4}{\value{n} - 3 - \j + 1} |
| \ifnum\value{t1}>0\ifnum\value{t2}>0 |
| \ifnum\value{t3}>0\ifnum\value{t4}>0 |
| \draw[thick,->] (\i,\j) to[in=-20,out=20] (\i,\j+3); |
| \fi\fi\fi\fi |
| \setcounter{t1}{2 * \j - 1 - \i + 1} |
| \setcounter{t2}{\value{n} - 1 - \i + 1} |
| \setcounter{t3}{2 * \i - 1 - \j + 1} |
| \setcounter{t4}{\value{n} - 1 - \j + 1} |
| \ifnum\value{t1}>0\ifnum\value{t2}>0 |
| \ifnum\value{t3}>0\ifnum\value{t4}>0 |
| \draw[thick,->] (\i,\j) to (\i+1,\j+1); |
| \fi\fi\fi\fi |
| } |
| } |
| \end{tikzpicture} |
| \end{center} |
| \caption{The relation from \autoref{ex:decomposition}} |
| \label{f:decomposition} |
| \end{figure} |
| \begin{example} |
| \label{ex:decomposition} |
| Consider the relation on the right of \textcite[Figure~2]{Beletska2009}, |
| reproduced in \autoref{f:decomposition}. |
| The relation can be described as $R = R_1 \cup R_2 \cup R_3$, |
| with |
| $$ |
| \begin{aligned} |
| R_1 &= n \mapsto \{\, (i,j) \to (i+3,j) \mid |
| i \le 2 j - 4 \wedge |
| i \le n - 3 \wedge |
| j \le 2 i - 1 \wedge |
| j \le n \,\} |
| \\ |
| R_2 &= n \mapsto \{\, (i,j) \to (i,j+3) \mid |
| i \le 2 j - 1 \wedge |
| i \le n \wedge |
| j \le 2 i - 4 \wedge |
| j \le n - 3 \,\} |
| \\ |
| R_3 &= n \mapsto \{\, (i,j) \to (i+1,j+1) \mid |
| i \le 2 j - 1 \wedge |
| i \le n - 1 \wedge |
| j \le 2 i - 1 \wedge |
| j \le n - 1\,\} |
| . |
| \end{aligned} |
| $$ |
| The figure shows this relation for $n = 7$. |
| Both |
| $R_3 \circ R_1 \subseteq R_1 \circ R_3$ |
| and |
| $R_3 \circ R_2 \subseteq R_2 \circ R_3$, |
| which the reader can verify using the {\tt iscc} calculator: |
| \begin{verbatim} |
| R1 := [n] -> { [i,j] -> [i+3,j] : i <= 2 j - 4 and i <= n - 3 and |
| j <= 2 i - 1 and j <= n }; |
| R2 := [n] -> { [i,j] -> [i,j+3] : i <= 2 j - 1 and i <= n and |
| j <= 2 i - 4 and j <= n - 3 }; |
| R3 := [n] -> { [i,j] -> [i+1,j+1] : i <= 2 j - 1 and i <= n - 1 and |
| j <= 2 i - 1 and j <= n - 1 }; |
| (R1 . R3) - (R3 . R1); |
| (R2 . R3) - (R3 . R2); |
| \end{verbatim} |
| $R_3$ can therefore be moved forward in any path. |
| For the other two basic relations, we have both |
| $R_2 \circ R_1 \not\subseteq R_1 \circ R_2$ |
| and |
| $R_1 \circ R_2 \not\subseteq R_2 \circ R_1$ |
| and so $R_1$ and $R_2$ form a strongly connected component. |
| By computing the power of $R_3$ and $R_1 \cup R_2$ separately |
| and composing the results, the power of $R$ can be computed exactly |
| using \eqref{eq:transitive:singleton}. |
| As explained by \textcite{Beletska2009}, applying the same formula |
| to $R$ directly, without a decomposition, would result in |
| an overapproximation of the power. |
| \end{example} |
| |
| \subsection{Partitioning the domains and ranges of $R$} |
| |
| The algorithm of \autoref{s:power} assumes that the input relation $R$ |
| can be treated as a union of translations. |
| This is a reasonable assumption if $R$ maps elements of a given |
| abstract domain to the same domain. |
| However, if $R$ is a union of relations that map between different |
| domains, then this assumption no longer holds. |
| In particular, when an entire dependence graph is encoded |
| in a single relation, as is done by, e.g., |
| \textcite[Section~6.1]{Barthou2000MSE}, then it does not make |
| sense to look at differences between iterations of different domains. |
| Now, arguably, a modified Floyd-Warshall algorithm should |
| be applied to the dependence graph, as advocated by |
| \textcite{Kelly1996closure}, with the transitive closure operation |
| only being applied to relations from a given domain to itself. |
| However, it is also possible to detect disjoint domains and ranges |
| and to apply Floyd-Warshall internally. |
| |
| \LinesNumbered |
| \begin{algorithm} |
| \caption{The modified Floyd-Warshall algorithm of |
| \protect\textcite{Kelly1996closure}} |
| \label{a:Floyd} |
| \SetKwInput{Input}{Input} |
| \SetKwInput{Output}{Output} |
| \Input{Relations $R_{pq}$, $0 \le p, q < n$} |
| \Output{Updated relations $R_{pq}$ such that each relation |
| $R_{pq}$ contains all indirect paths from $p$ to $q$ in the input graph} |
| % |
| \BlankLine |
| \SetAlgoVlined |
| \DontPrintSemicolon |
| % |
| \For{$r \in [0, n-1]$}{ |
| $R_{rr} \coloneqq R_{rr}^+$ \nllabel{l:Floyd:closure}\; |
| \For{$p \in [0, n-1]$}{ |
| \For{$q \in [0, n-1]$}{ |
| \If{$p \ne r$ or $q \ne r$}{ |
| $R_{pq} \coloneqq R_{pq} \cup \left(R_{rq} \circ R_{pr}\right) |
| \cup \left(R_{rq} \circ R_{rr} \circ R_{pr}\right)$ |
| \nllabel{l:Floyd:update} |
| } |
| } |
| } |
| } |
| \end{algorithm} |
| |
| Let the input relation $R$ be a union of $m$ basic relations $R_i$. |
| Let $D_{2i}$ be the domains of $R_i$ and $D_{2i+1}$ the ranges of $R_i$. |
| The first step is to group overlapping $D_j$ until a partition is |
| obtained. If the resulting partition consists of a single part, |
| then we continue with the algorithm of \autoref{s:power}. |
| Otherwise, we apply Floyd-Warshall on the graph with as vertices |
| the parts of the partition and as edges the $R_i$ attached to |
| the appropriate pairs of vertices. |
| In particular, let there be $n$ parts $P_k$ in the partition. |
| We construct $n^2$ relations |
| $$ |
| R_{pq} \coloneqq \bigcup_{i \text{ s.t. } \domain R_i \subseteq P_p \wedge |
| \range R_i \subseteq P_q} R_i |
| , |
| $$ |
| apply \autoref{a:Floyd} and return the union of all resulting |
| $R_{pq}$ as the transitive closure of $R$. |
| Each iteration of the $r$-loop in \autoref{a:Floyd} updates |
| all relations $R_{pq}$ to include paths that go from $p$ to $r$, |
| possibly stay there for a while, and then go from $r$ to $q$. |
| Note that paths that ``stay in $r$'' include all paths that |
| pass through earlier vertices since $R_{rr}$ itself has been updated |
| accordingly in previous iterations of the outer loop. |
| In principle, it would be sufficient to use the $R_{pr}$ |
| and $R_{rq}$ computed in the previous iteration of the |
| $r$-loop in Line~\ref{l:Floyd:update}. |
| However, from an implementation perspective, it is easier |
| to allow either or both of these to have been updated |
| in the same iteration of the $r$-loop. |
| This may result in duplicate paths, but these can usually |
| be removed by coalescing (\autoref{s:coalescing}) the result of the union |
| in Line~\ref{l:Floyd:update}, which should be done in any case. |
| The transitive closure in Line~\ref{l:Floyd:closure} |
| is performed using a recursive call. This recursive call |
| includes the partitioning step, but the resulting partition will |
| usually be a singleton. |
| The result of the recursive call will either be exact or an |
| overapproximation. The final result of Floyd-Warshall is therefore |
| also exact or an overapproximation. |
| |
| \begin{figure} |
| \begin{center} |
| \begin{tikzpicture}[x=1cm,y=1cm,>=stealth,shorten >=3pt] |
| \foreach \x/\y in {0/0,1/1,3/2} { |
| \fill (\x,\y) circle (2pt); |
| } |
| \foreach \x/\y in {0/1,2/2,3/3} { |
| \draw (\x,\y) circle (2pt); |
| } |
| \draw[->] (0,0) -- (0,1); |
| \draw[->] (0,1) -- (1,1); |
| \draw[->] (2,2) -- (3,2); |
| \draw[->] (3,2) -- (3,3); |
| \draw[->,dashed] (2,2) -- (3,3); |
| \draw[->,dotted] (0,0) -- (1,1); |
| \end{tikzpicture} |
| \end{center} |
| \caption{The relation (solid arrows) on the right of Figure~1 of |
| \protect\textcite{Beletska2009} and its transitive closure} |
| \label{f:COCOA:1} |
| \end{figure} |
| \begin{example} |
| Consider the relation on the right of Figure~1 of |
| \textcite{Beletska2009}, |
| reproduced in \autoref{f:COCOA:1}. |
| This relation can be described as |
| $$ |
| \begin{aligned} |
| \{\, (x, y) \to (x_2, y_2) \mid {} & (3y = 2x \wedge x_2 = x \wedge 3y_2 = 3 + 2x \wedge x \ge 0 \wedge x \le 3) \vee {} \\ |
| & (x_2 = 1 + x \wedge y_2 = y \wedge x \ge 0 \wedge 3y \ge 2 + 2x \wedge x \le 2 \wedge 3y \le 3 + 2x) \,\} |
| . |
| \end{aligned} |
| $$ |
| Note that the domain of the upward relation overlaps with the range |
| of the rightward relation and vice versa, but that the domain |
| of neither relation overlaps with its own range or the domain of |
| the other relation. |
| The domains and ranges can therefore be partitioned into two parts, |
| $P_0$ and $P_1$, shown as the white and black dots in \autoref{f:COCOA:1}, |
| respectively. |
| Initially, we have |
| $$ |
| \begin{aligned} |
| R_{00} & = \emptyset |
| \\ |
| R_{01} & = |
| \{\, (x, y) \to (x+1, y) \mid |
| (x \ge 0 \wedge 3y \ge 2 + 2x \wedge x \le 2 \wedge 3y \le 3 + 2x) \,\} |
| \\ |
| R_{10} & = |
| \{\, (x, y) \to (x_2, y_2) \mid (3y = 2x \wedge x_2 = x \wedge 3y_2 = 3 + 2x \wedge x \ge 0 \wedge x \le 3) \,\} |
| \\ |
| R_{11} & = \emptyset |
| . |
| \end{aligned} |
| $$ |
| In the first iteration, $R_{00}$ remains the same ($\emptyset^+ = \emptyset$). |
| $R_{01}$ and $R_{10}$ are therefore also unaffected, but |
| $R_{11}$ is updated to include $R_{01} \circ R_{10}$, i.e., |
| the dashed arrow in the figure. |
| This new $R_{11}$ is obviously transitively closed, so it is not |
| changed in the second iteration and it does not have an effect |
| on $R_{01}$ and $R_{10}$. However, $R_{00}$ is updated to |
| include $R_{10} \circ R_{01}$, i.e., the dotted arrow in the figure. |
| The transitive closure of the original relation is then equal to |
| $R_{00} \cup R_{01} \cup R_{10} \cup R_{11}$. |
| \end{example} |
| |
| \subsection{Incremental Computation} |
| \label{s:incremental} |
| |
| In some cases it is possible and useful to compute the transitive closure |
| of union of basic relations incrementally. In particular, |
| if $R$ is a union of $m$ basic maps, |
| $$ |
| R = \bigcup_j R_j |
| , |
| $$ |
| then we can pick some $R_i$ and compute the transitive closure of $R$ as |
| \begin{equation} |
| \label{eq:transitive:incremental} |
| R^+ = R_i^+ \cup |
| \left( |
| \bigcup_{j \ne i} |
| R_i^* \circ R_j \circ R_i^* |
| \right)^+ |
| . |
| \end{equation} |
| For this approach to be successful, it is crucial that each |
| of the disjuncts in the argument of the second transitive |
| closure in \eqref{eq:transitive:incremental} be representable |
| as a single basic relation, i.e., without a union. |
| If this condition holds, then by using \eqref{eq:transitive:incremental}, |
| the number of disjuncts in the argument of the transitive closure |
| can be reduced by one. |
| Now, $R_i^* = R_i^+ \cup \identity$, but in some cases it is possible |
| to relax the constraints of $R_i^+$ to include part of the identity relation, |
| say on domain $D$. We will use the notation |
| ${\cal C}(R_i,D) = R_i^+ \cup \identity_D$ to represent |
| this relaxed version of $R^+$. |
| \textcite{Kelly1996closure} use the notation $R_i^?$. |
| ${\cal C}(R_i,D)$ can be computed by allowing $k$ to attain |
| the value $0$ in \eqref{eq:transitive:Q} and by using |
| $$ |
| P \cap \left(D \to D\right) |
| $$ |
| instead of \eqref{eq:transitive:approx}. |
| Typically, $D$ will be a strict superset of both $\domain R_i$ |
| and $\range R_i$. We therefore need to check that domain |
| and range of the transitive closure are part of ${\cal C}(R_i,D)$, |
| i.e., the part that results from the paths of positive length ($k \ge 1$), |
| are equal to the domain and range of $R_i$. |
| If not, then the incremental approach cannot be applied for |
| the given choice of $R_i$ and $D$. |
| |
| In order to be able to replace $R^*$ by ${\cal C}(R_i,D)$ |
| in \eqref{eq:transitive:incremental}, $D$ should be chosen |
| to include both $\domain R$ and $\range R$, i.e., such |
| that $\identity_D \circ R_j \circ \identity_D = R_j$ for all $j\ne i$. |
| \textcite{Kelly1996closure} say that they use |
| $D = \domain R_i \cup \range R_i$, but presumably they mean that |
| they use $D = \domain R \cup \range R$. |
| Now, this expression of $D$ contains a union, so it not directly usable. |
| \textcite{Kelly1996closure} do not explain how they avoid this union. |
| Apparently, in their implementation, |
| they are using the convex hull of $\domain R \cup \range R$ |
| or at least an approximation of this convex hull. |
| We use the simple hull (\autoref{s:simple hull}) of $\domain R \cup \range R$. |
| |
| It is also possible to use a domain $D$ that does {\em not\/} |
| include $\domain R \cup \range R$, but then we have to |
| compose with ${\cal C}(R_i,D)$ more selectively. |
| In particular, if we have |
| \begin{equation} |
| \label{eq:transitive:right} |
| \text{for each $j \ne i$ either } |
| \domain R_j \subseteq D \text{ or } \domain R_j \cap \range R_i = \emptyset |
| \end{equation} |
| and, similarly, |
| \begin{equation} |
| \label{eq:transitive:left} |
| \text{for each $j \ne i$ either } |
| \range R_j \subseteq D \text{ or } \range R_j \cap \domain R_i = \emptyset |
| \end{equation} |
| then we can refine \eqref{eq:transitive:incremental} to |
| $$ |
| R_i^+ \cup |
| \left( |
| \left( |
| \bigcup_{\shortstack{$\scriptstyle\domain R_j \subseteq D $\\ |
| $\scriptstyle\range R_j \subseteq D$}} |
| {\cal C} \circ R_j \circ {\cal C} |
| \right) |
| \cup |
| \left( |
| \bigcup_{\shortstack{$\scriptstyle\domain R_j \cap \range R_i = \emptyset$\\ |
| $\scriptstyle\range R_j \subseteq D$}} |
| \!\!\!\!\! |
| {\cal C} \circ R_j |
| \right) |
| \cup |
| \left( |
| \bigcup_{\shortstack{$\scriptstyle\domain R_j \subseteq D $\\ |
| $\scriptstyle\range R_j \cap \domain R_i = \emptyset$}} |
| \!\!\!\!\! |
| R_j \circ {\cal C} |
| \right) |
| \cup |
| \left( |
| \bigcup_{\shortstack{$\scriptstyle\domain R_j \cap \range R_i = \emptyset$\\ |
| $\scriptstyle\range R_j \cap \domain R_i = \emptyset$}} |
| \!\!\!\!\! |
| R_j |
| \right) |
| \right)^+ |
| . |
| $$ |
| If only property~\eqref{eq:transitive:right} holds, |
| we can use |
| $$ |
| R_i^+ \cup |
| \left( |
| \left( |
| R_i^+ \cup \identity |
| \right) |
| \circ |
| \left( |
| \left( |
| \bigcup_{\shortstack{$\scriptstyle\domain R_j \subseteq D $}} |
| R_j \circ {\cal C} |
| \right) |
| \cup |
| \left( |
| \bigcup_{\shortstack{$\scriptstyle\domain R_j \cap \range R_i = \emptyset$}} |
| \!\!\!\!\! |
| R_j |
| \right) |
| \right)^+ |
| \right) |
| , |
| $$ |
| while if only property~\eqref{eq:transitive:left} holds, |
| we can use |
| $$ |
| R_i^+ \cup |
| \left( |
| \left( |
| \left( |
| \bigcup_{\shortstack{$\scriptstyle\range R_j \subseteq D $}} |
| {\cal C} \circ R_j |
| \right) |
| \cup |
| \left( |
| \bigcup_{\shortstack{$\scriptstyle\range R_j \cap \domain R_i = \emptyset$}} |
| \!\!\!\!\! |
| R_j |
| \right) |
| \right)^+ |
| \circ |
| \left( |
| R_i^+ \cup \identity |
| \right) |
| \right) |
| . |
| $$ |
| |
| It should be noted that if we want the result of the incremental |
| approach to be transitively closed, then we can only apply it |
| if all of the transitive closure operations involved are exact. |
| If, say, the second transitive closure in \eqref{eq:transitive:incremental} |
| contains extra elements, then the result does not necessarily contain |
| the composition of these extra elements with powers of $R_i$. |
| |
| \subsection{An {\tt Omega}-like implementation} |
| |
| While the main algorithm of \textcite{Kelly1996closure} is |
| designed to compute and underapproximation of the transitive closure, |
| the authors mention that they could also compute overapproximations. |
| In this section, we describe our implementation of an algorithm |
| that is based on their ideas. |
| Note that the {\tt Omega} library computes underapproximations |
| \parencite[Section 6.4]{Omega_lib}. |
| |
| The main tool is Equation~(2) of \textcite{Kelly1996closure}. |
| The input relation $R$ is first overapproximated by a ``d-form'' relation |
| $$ |
| \{\, \vec i \to \vec j \mid \exists \vec \alpha : |
| \vec L \le \vec j - \vec i \le \vec U |
| \wedge |
| (\forall p : j_p - i_p = M_p \alpha_p) |
| \,\} |
| , |
| $$ |
| where $p$ ranges over the dimensions and $\vec L$, $\vec U$ and |
| $\vec M$ are constant integer vectors. The elements of $\vec U$ |
| may be $\infty$, meaning that there is no upper bound corresponding |
| to that element, and similarly for $\vec L$. |
| Such an overapproximation can be obtained by computing strides, |
| lower and upper bounds on the difference set $\Delta \, R$. |
| The transitive closure of such a ``d-form'' relation is |
| \begin{equation} |
| \label{eq:omega} |
| \{\, \vec i \to \vec j \mid \exists \vec \alpha, k : |
| k \ge 1 \wedge |
| k \, \vec L \le \vec j - \vec i \le k \, \vec U |
| \wedge |
| (\forall p : j_p - i_p = M_p \alpha_p) |
| \,\} |
| . |
| \end{equation} |
| The domain and range of this transitive closure are then |
| intersected with those of the input relation. |
| This is a special case of the algorithm in \autoref{s:power}. |
| |
| In their algorithm for computing lower bounds, the authors |
| use the above algorithm as a substep on the disjuncts in the relation. |
| At the end, they say |
| \begin{quote} |
| If an upper bound is required, it can be calculated in a manner |
| similar to that of a single conjunct [sic] relation. |
| \end{quote} |
| Presumably, the authors mean that a ``d-form'' approximation |
| of the whole input relation should be used. |
| However, the accuracy can be improved by also trying to |
| apply the incremental technique from the same paper, |
| which is explained in more detail in \autoref{s:incremental}. |
| In this case, ${\cal C}(R_i,D)$ can be obtained by |
| allowing the value zero for $k$ in \eqref{eq:omega}, |
| i.e., by computing |
| $$ |
| \{\, \vec i \to \vec j \mid \exists \vec \alpha, k : |
| k \ge 0 \wedge |
| k \, \vec L \le \vec j - \vec i \le k \, \vec U |
| \wedge |
| (\forall p : j_p - i_p = M_p \alpha_p) |
| \,\} |
| . |
| $$ |
| In our implementation we take as $D$ the simple hull |
| (\autoref{s:simple hull}) of $\domain R \cup \range R$. |
| To determine whether it is safe to use ${\cal C}(R_i,D)$, |
| we check the following conditions, as proposed by |
| \textcite{Kelly1996closure}: |
| ${\cal C}(R_i,D) - R_i^+$ is not a union and for each $j \ne i$ |
| the condition |
| $$ |
| \left({\cal C}(R_i,D) - R_i^+\right) |
| \circ |
| R_j |
| \circ |
| \left({\cal C}(R_i,D) - R_i^+\right) |
| = |
| R_j |
| $$ |
| holds. |