next up previous Next: Optimization Up: ALIAS-Maple Previous: The solving procedures

  • La page de J-P. Merlet
  • J-P. Merlet home page
  • La page Présentation de HEPHAISTOS
  • HEPHAISTOS home page
  • La page "Présentation" de l'INRIA
  • INRIA home page

    Subsections


    Simplification procedures

    Introduction

    For a specific problem the end-user may have some knowledge of the problem that implies that he may determine that for some intervals for the variables the problem has no solution or that the intervals for some variables may be reduced.

    ALIAS-Maple allows one to introduce such knowledge in the generated C++ code. The mechanism that is used is to allow one to specify a C++ routine, that will be called a simplification procedure. We will present in this chapter simplification procedures devoted to system solving. A simplification procedure takes as input a set of intervals for the variables and will return an integer:

    Procedures returning -1,0 or 1 will be called filtering procedures while procedures returning 11 are called root procedures.

    This routine will be called by all the solving procedure of ALIAS-Maple whenever a new set of intervals is considered in the algorithm i.e. right after the bisection process (see the ALIAS-C++ manual). System solving procedures in ALIAS C++ are designed to discard box for which a simplification procedure returns -1 and to update it if the simplification procedure returns 1.

    The syntax of such a routine Simp is: 

     
int Simp(INTERVAL_VECTOR &v_IS)
{
....
return 0;
}

    Alternatively ALIAS-Maple allows to automatically produce simplification procedures according to classical methods used in interval analysis.

    Note that numerical round-off errors in Maple may produce solving errors as the filtering procedures uses frequently an expanded version of the expressions. To reduce the impact of this problem it is possible to adjust the number of digits that will be used when doing numerical computation in the filtering procedure by using the ALIAS variable `ALIAS/digits` which is set by default to 40.


    Filtering simplification procedures

    An important point: some of these procedures have as argument a list of variable names and this list must be the same than the one that is used for the solving procedure.


    The HullConsistency, Simp2B and HullIConsistency procedures

    The idea underlying this simplification procedure (known in the community of constraint programming as the hull-consistency or 2B-consistency) is to rewrite the equation and check if it is consistent at each time. For example imagine that one of the equation is $x^2-2x+1=0$. The procedure will introduce a new variable $X$ such that $X=2x-1$ and compute its interval evaluation. As $X$ should be equal to $x^2$ if $X$ has a negative upper bound the simplification procedure will return -1, which indicate to the solving algorithm that it may discard the current input interval. If the upper bound $U$ of $X$ is positive the procedure will check if $x$ lie in $[-\sqrt{U},\sqrt{U}]$ and update the interval for $x$ if this is not the case.

    In a second step the procedure will introduce another new variable $Y$ such that $Y=(x^2+1)/2$ and compute its interval evaluation. As $Y$ should be equal to $x$ the procedure will check if the interval value for $Y$ is consistent with the interval value for $x$ and if not update the interval for $x$ accordingly.

    The HullConsistency and HullIConsistency procedures of ALIAS-Maple partly implements this principle respectively for equalities and inequalities. They may be used if any of the equation has linear or square terms in the unknown (e.g. $x*y, x^2*y$). The procedure will isolate the linear and square terms of each unknown and produce the C++ code that will check is their interval evaluation is consistent.

    The syntax of the HullConsistency procedure is: 

     
HullConsistency([x^2+y^2-1,x*y-1],[x,y],"Simp");
    The first argument is a list of equations, the second a list of variable names and the third argument is the name of a file. In the above example the procedure will create the C++ program Simp.C that may be later be given to a solving procedure.

    Until version 2.4 was a HullConsistencyStrong procedure available, that has been contributed by G. Chabert. The procedure was intended to deal with terms that were different from linear and square terms that are managed by HullConsistency. Its drawback was that it was based on a semantic tree of the expression and was not able to decrease the dependency problem. For example if the expression is

    \begin{displaymath} x^2+x*y+x \end{displaymath}

    and the final node is $x^2$, then the 2B will be implemented using

    \begin{displaymath} x^2 = -x*y-x \end{displaymath}

    while HullConsistency will use

    \begin{displaymath} x^2= -x*(y-1) \end{displaymath}

    which is more efficient. For example in the equation

    \begin{displaymath} \sin(x_1+x_2)\cos(x_1^2)*x_2+x_1=0 \end{displaymath}

    we will use

    \begin{displaymath} \sin(x_1+x_2)=\frac{-x_1}{\cos(x_1^2)*x_2}~~~x_2=\frac{-x_1}{\sin(x_1+x_2)*\cos(x_1^2)} ~~~x_1=-\sin(x_1+x_2)\cos(x_1^2)*x_2 \end{displaymath}

    to update successively the range for $(x_1, x_2),(x_2),(x_1)$. Note that there is not an extensive check of the possibility of interval evaluating the right hand side term of these new equations and thus the simplification procedure may lead to an error during this evaluation. The syntax of Simp2B is similar to the one of Hullconsistency.

    The new procedure Simp2B manages in a better way the dependency problem but is less complete than HullConsistencyStrong, which is no more available for Maple 9.5, as it manages only the following expressions

    \begin{displaymath} \sin, \cos, \tan, \log, exp, x^n (n~being~an~integer) \end{displaymath}

    It is also usually more computer intensive than the HullConsistencyStrong procedure. Note that you have some control over the simplification rule that will be used. Consider for example the equation

    \begin{displaymath} x^2+ e^z y+1=0 \end{displaymath}

    By default Simp2B will use the following equalities: $x^2 = -1 -e^z y$, $e^z = (-1-x^2)/y$, $y= (-1-x^2)/e^z$ to reduce the range respectively for $x, z,y$. But if you specify in the list ALIAS_FORBIDDEN_TERMS: 
     
ALIAS_FORBIDDEN_TERM:=[e^z]:
    then only the first and third rule will be used.

    The syntax of the HullIConsistency procedure is: 

     
HullIConsistency([x^2+y^2-1<=0,x*y-1<=0],[x,y],0,"Simp");
    Here the third argument is an integer which should be either 0 or 1. It 1 it indicates that the derivatives of the left side of the inequalities are used for their interval evaluation.

    There is an optional last argument to these procedures which is a floating point number $f$. If during the simplification there is change in the interval for a variable such that the width of the new range is lower by more than $f$ from the width of the initial range, then the simplification process will be repeated until either the procedure detect that there is no solution to the system or that no change occur in the variable or that the amplitude of the change is lower than $f$.

    If you have numerous equations that have such terms the C++ procedure may be quite large. To reduce this size you may choose to look only at the linear term or only at the square term by setting the flag `ALIAS/hullC` to 1 or 2 (its default value is 0). Note that this option is not valid for the HullConsistencyStrong and HullIConsistency procedures.

    A typical use of this procedure is given below, assuming that you have a list of equations EQ and a list of variable VAR: 

     
HullConsistency(EQ,VAR,"SS"):
HessianSolve(EQ,VAR,Init,"SS");
    Here the procedure generate the C++ procedure SS (which will be written in the file SS.C) and the solving algorithm HessianSolve will use this procedure during the computation.

    The efficiency of the hull consistency approach is usually quite good but decreases quite quickly after the first improvements. Hence it is in general not a good idea to repeat the procedure when only small improvements on the range are obtained.

    A test to determine if an expression involved in the equations and inequalities can be interval-valuated may be performed if the Verify_Problem_Expression procedure has been called before the call to HullConsistency or HullIConsistency (see section 2.1.5). If different calls to the consistency procedures are made with different set of expressions it is necessary to specify a different string `ALIAS/ID` string before each call to Verify_Problem_Expression. For example 

     
 `ALIAS/ID`:="Simp1":
 Verify_Problem_Expression(EQ1,VAR):
 HullIConsistency(EQ1,VAR,"Simp1"):
 `ALIAS/ID`:="Simp2":
 Verify_Problem_Expression(EQ2,VAR):
 HullConsistency(EQ2,VAR,"Simp2"):
 `ALIAS/ID`:="Simp3":
 Verify_Problem_Expression(EQ2,VAR):
 HullConsistencyStrong(EQ3,VAR,"Simp3"):
 `ALIAS/ID`:="other":
    is valid.

    A special key-word may be used for the procedures HullConsistency and HullIConsistency in the case of an optimization problem. In that case the ALIAS-C++ variable ALIAS_Optimize is set to -1 (minimum problem) or 1 (maximum problem), the flag ALIAS_Has_Optimum is set to 1 as soon as an estimation of the optimum is found while the interval ALIAS_Optimum is set to the estimation of the optimum. If an equation or an inequality has the key-word "Optimum" in it and if ALIAS_Has_Optimum is 1, then "Optimum" will be substituted by Sup(ALIAS_Optimum) (minimum problem) or Inf(ALIAS_Optimum) (maximum problem).

    For example assume that you deal with a problem involving the minimum of $x\sin(x)+y\cos(y)$: you may then use the following call 

     
HullIConsistency([x*sin(x)+y*cos(y)-Optimum<=0],[x,y],0,"Simp");
    to create a simplification procedure that may reduce the current box. See another example in section 5.2

    The HullConsistencyTaylor procedure

    This procedure uses a third order Taylor extension of an equation and an interval evaluation of the remainder. This leads to a second order equation in each of the variables, whose coefficients is computed. The analytic solution of this second order equation is used to improve the variable range. This filter is very efficient for non algebraic equations and if the width of the range for the variable is not too large.

    The syntax is: 

     
HullConsistencyTaylor(func,vars,procname)

    An optional fourth argument may be used. It consists in a list of one or two number. If this argument is [0.1,0.01], then the filter will be used only if the maximal width for the ranges is less than 0.1 and the filter will be used again if the gain in a range if larger than 0.01. If only one number is given in the list, then it will indicate the maximal width for the ranges.

    The BiCenteredForm procedure

    This procedure is based on the interval evaluation of equations using centered form with as centers two possible centers which may be optimal (see the ALIAS C++ manual for the calculation of these centers). The centered form uses the fact the interval evaluation of a function $f$ in the $n$ unknown $X=\{x_1,\ldots,x_n\}$ is included in $f(c)+J(X)(X-c)$ where $c$ is point called the center. As mentioned earlier two possible centers are calculated. Furthermore the calculation of the product $J(X)(X-c)$ may involve several occurrences of the same variables, that are not detected when computing the product numerically. This procedure allows, as an option, to compute symbolically the product and to re-arrange it in order to try to reduce multiple occurrences.

    The syntax of this procedure is: 

     
BiCenteredForm(func,funcproc,Jfuncproc,vars,custom,procname)
    where In addition the procedure admits a 7th optional argument "prod". In that case the procedure will compute symbolically the product $J(X)(X-c)$ and re-arrange its terms. Both the numerical form and the one resulting from the symbolic calculation will be used in the simplification procedure.

    The GlobalConsistencyTaylor procedure

    This procedure uses a second order Taylor extension of an equation and an interval evaluation of the remainder. The equation system is transformed into a linear part and an interval remainder and an interval linear system solver is used to determine the bounds for the solutions. Hence this filter is global (i.e. it takes into account all equations of a system) while most of the other filter are local (i.e. they deal only with one equation at a time). This filter is usually very efficient for non algebraic equations and if the width of the range for the variable is not too large.

    The syntax is: 

     
GlobalConsistencyTaylor(func,vars,procname)

    An optional fourth argument may be used. It consists in a list of one or two number. If this argument is [0.1,0.01], then the filter will be used only if the maximal width for the ranges is less than 0.1 and the filter will be used again if the gain in a range if larger than 0.01. If only one number is given in the list, then it will indicate the maximal width for the ranges.

    The SimpAngle procedure

    This procedure may be used if one of the equation indicates that the scalar product of two vectors must have some given value. Hence such equation may be written as:

    \begin{displaymath} \sum (x_i-x_j)(x_k-x_m) =0 \end{displaymath}

    where the $x_i, x_k$ are unknowns and the $x_j, x_m$ unknowns or constants.

    A typical use of this procedure is given below, assuming that you have a list of equations EQ and a list of variable VAR: 

     
eq1:=(x1-3)*(x2-4)+(y1-1)*(y2-4)-10:
SimpAngle([eq1],VAR,"SS"):
GeneralSolve(EQ,VAR,Init,"SS"):
    Here the procedure generate the C++ procedure SS (which will be written in the file SS.C) and the solving algorithm GeneralSolve will use this procedure during the computation.

    There is an optional last argument to this procedure which is a floating point number $f$. If during the simplification there is change in the interval for a variable such that the width of the new range is lower by more than $f$ from the width of the initial range, then the simplification process will be repeated until either the procedure detect that there is no solution to the system or that no change occur in the variable or that the amplitude of the change is lower than $f$.

    The SimplexConsistency procedure

    This simplification procedure may be used to solve systems of equations as soon as the system has at least two equations with polynomial terms in at least one variable. The principle used for this simplification procedure is explained in the ALIAS-C++ manual section of Solve_Simplex.

    The syntax of this procedure is: 

     
\SimplexConsistency(EQ,VAR,"Simp"):
    where EQ is a list of equation, VAR a list of variable name and Simp is the name of the simplification procedure that will be written in the file Simp.C.

    The behavior of this simplification procedure may be modified by using the following variables:

    Two remarks about this procedure:

    1. it is usually very effective
    2. it may be computer intensive for large systems: hence its use should in general be avoided within the 3B method (see section 4.5)
    3. it is not numerically safe with the current implementation of the simplex i.e. in some cases it may lead to miss some solutions

    The IntervalNewton procedure

    This procedure allows one to take into account the specific structure of the system to solve when using the interval Newton method. In this method the classical algorithm, for example used by default in GradientSolve or HessianSolve, computes numerically for a box $X$ the interval evaluation of $M=J^{-1}(X_m)J(X)$ where $X_m$ is the mid-point of $X$. This numerical calculation does not allow to to take into account the dependency between the elements of $J$ when computing $M$. This procedure compute symbolically the matrix $M$ and re-arrange its elements (with MinimalCout) in order to reduce the dependency problem.

    The syntax of this procedure is 

    IntervalNewton(func,nfunc,funcproc,njfunc,Jfuncproc,typemid,grad,grad3B,incl,vars,procname)
    where

    Roots simplification procedures

    A special type of simplification procedure is the one that are used not for filtering boxes but for finding roots of a system of equations. The general solving algorithm of ALIAS accept such type of procedure that returns 11 if a root is provided by the simplification procedure.

    TryNewton

    This procedure will create a C++ program that is intended to be used within an equation solving procedure as a simplification procedure. For each call to this program the program will run at most newton_iteration iterations of the classical Newton scheme with as initial guess for the solution the mid-point of the interval vector that is the argument of a simplification procedure. If Newton seems to converge (i.e. the residues have a width lower than fepsilon) the program first checks that the approximation of the root lies within a given search domain and then uses the Kantorovitch theorem and the epsilon-inflation scheme to certify that there is a single solution within a given ball. In that case the program will return 11, otherwise it returns 0. The general use of this procedure is for finding quickly a root of the system although in some cases it may speed up the process of finding all solutions of the system.

    The syntax of this procedure is 

     
TryNewton(func,vars,Init,procname)

    This procedure may be used with GeneralSolve and GradientSolve. It should not be used with HessianSolve as the corresponding C++ program is already embedded in the C++ code of this procedure.

    Rouche

    The Rouche procedure has to be used for a system of equations and will create a simplification procedure that may return 11, i.e. it provide a ball that include a single solution of the system, that may be found using a Newton iteration. A key-point for using this procedure is to consider the matrix constituted of the $k$-th derivatives of the equations with respect to the variable, that will be denoted $F^(k)$ and to be able to find $k_0$ such that there is a $k_j$ in $[2,k_0]$

    \begin{displaymath} \frac{\vert\vert(F^{(1)}(X_0))^{-1}F^{(k)}(X_0)\vert\vert^{1... ...F^{(1)}(X_0))^{-1}F^{(k_j)}(X_0)\vert\vert^{1/(k_0-1)}}{k_j!} \end{displaymath}

    for all $X_0$ in the search space and all possible $k\ge 2$. The value $k_0$ will be called the order of the Rouche method. Two different cases will be considered here: For algebraic equations there will be clearly a $k_1$ such that $F^{(k)}=0$ for all $k\ge k_1$. Hence we may use $k1-1$ as order for the Rouche method.

    As soon as there are non algebraic terms in at least one equation, then we cannot determine in advance the order of the Rouche method without an in-depth analysis of the system.

    The syntax of this procedure is: 

     
Rouche(Func,nfunc,Order,Vars,ProcName)
    where A good example is 
     
EQ:=[x^2+y^2+z^2-9,(x-1)^2+(y-1)^2+z^2-9,(x-1)^2+(y+1)^2+z^2-16]]:
VAR:=[x,y,z]:
Rouche(EQ,"p",VAR,"rouche");

    This procedure will be called only if the width of the box is lower than ALIAS/maxnewton and uses at most ALIAS/newton_iteration of the Newton scheme to compute an approximation of the root such that the residues of the system are lower than ALIAS/fepsilon. The generated C++ program will also use the epsilon inflation method to enlarge as much as possible the ball that includes the root.

    The Rouche procedure may be quite powerful to find a ball that includes a single root of the system (even more powerful than the Kantorovitch scheme). For example by using Rouche in combination with DeflationUP (see section 8.1.2.1) we have been able to solve the Wilkinson polynomial of order 19 (see section 12.3), while the general procedure failed certifying roots starting at order 13.


    Concatenation of simplification procedures

    In some cases you may have the opportunity to write more than one simplification procedure (for example by using the HullConsistency procedure and writing you own simplification procedure). Currently the procedures in the C++ library and in the Maple library accept only one simplification procedure. But you may concatenate an unlimited number of simplification procedures by using the procedure CatSimp. The first argument of this procedure is the name of the procedure that will be used as unique simplification procedures by ALIAS-C++. The following arguments are names of simplification procedures. Note that the name of the file containing the simplification procedure must always be the name of the procedure with the extension .C. Hence 

    CatSimp("WholeSimp","Simp1","Simp2")
    will produce the file WholeSimp.C defining the procedure WholeSimp which in turn calls the procedures Simp1,Simp2 defined in the file Simp1.C, Simp2.C. Hence the arguments of the CatSimp procedure are usually strings.

    This procedure accepts a numerical value as optional last argument: in that case the whole chain of simplification procedures will be repeated if the width of at least one range has decreased by more than the specified numerical value during the previous run of the simplification procedures.

    Some of the simplification procedures may be computer intensive and therefore their use in the 3B method (see section 4.5) may not be recommended. For that purpose it is possible to specify within the CatSimp procedure that a given simplification procedure shall not be used within the 3B method: for that purpose just add the string NO_3B right after the procedure name. Hence: 

     
CatSimp("Simp","Simp1","Simp2","NO_3B"):
    will create the simplification procedure Simp that will usually calls Simp1 and Simp2 except for the 3B method where only Simp1 will be used.


    Using the 3B method

    Principle

    The 3B method is an usually efficient method to reduce the width of the variable intervals. This method can be used for any algorithm of ALIAS and is embedded in the code of the solving procedures of the C++ library.

    The principle of the method is quite simple. Let assume that we deal with a problem with $n$ variables $x_1,\ldots,x_n$ and that the range for the variable $x_i$ is [$a_1,b_1$] while $\epsilon$ is a small number. The procedure will examine if the problem may have a solution if the range for $x_i$ is reduced to [ $a_1,a_1+\epsilon$]. If it is not the case the range for $x_i$ may be reduced to [$a_2,b_2$]=[ $a_1+\epsilon,b_1$]. The procedure will then be repeated by using $2\epsilon$ as a new value for $\epsilon$ until the algorithm either determine that there is no solution to the problem or fails to show that that is no solution for the range [ $a_n,a_n+n\epsilon$]. The algorithm will then try to reduce the range for $x_i$ on the "right" side i.e. by examining if the problem has no solution for the range [ $b_1-\epsilon, b_1$]. This process is repeated in turn for each variable.

    Within Maple the value of $\epsilon$ is given by the variable `ALIAS/Delta3B` which set an identical value of $\epsilon$ for all variables. Alternatively you may specify a specific value of $\epsilon$ for each variable by assigning a list to the variable `ALIAS/Delta3B_ARRAY` with as many elements as unknowns: each element of this list is the value of $\epsilon$ for each variable.

    The 3B method will not be used for variables whose range have a width greater than a given threshold. This threshold is defined by the variable `ALIAS/Max3B` (default value: 5). You may assign a new value for this threshold by: 

     
`ALIAS/Max3B`:=10:
    Alternatively you may specify a specific threshold for each variable by using the array `ALIAS/Max3B_ARRAY`. For example: 
     
`ALIAS/Max3B_ARRAY:=array([0.1,0.1,0.05]):
    indicates that the 3B method will be used only if the width or the input intervals is lower than 0.1 for the 2 first unknowns and 0.005 for the third one.

    As soon as the Max3B, Delta3B values have been set you activate the 3B method by setting the `ALIAS/3B` variable to 1 or 2. At iteration $n$ of the 3B method the algorithm will try to reduce the range for the variable by $n\epsilon$. Assume that $n$ is larger than 1 and the algorithm fails:

    Repeating the 3B method

    An improvement may be obtained for the variable $x_i$ with $i>1$ while no improvement was obtained for the variable $x_1,\ldots,x_{i-1}$. But during the calculation for these variables the range for the variable $x_i$ was the initial one and not the improved value. Hence it may be interesting to repeat the procedure as soon as one the variable is improved.

    The 3B procedure will be repeated if:

    The 3B method usually exhibit an asymptotically slow convergence speed and it is hence bad policy of repeating the process for very small improvements in the range.

    Other parameters for the 3B method

    At one step of the algorithm the 3B method will determine that there is no solution either by considering the simplification procedures provided by the end-user or by a direct evaluation of the expressions involved in the problem.

    As the simplification procedures may be computer intensive you may avoid using them within the 3B method by setting the flag `ALIAS/Use_Simp_3B` to 0. Alternatively you may specify that only some of the simplification procedures are to be used within the 3B method by using the mechanism described in the CatSimp procedure (section 4.4). You may also set `ALIAS/Use_Simp_3B` to 2: in that case if some variable ranges have been improved the algorithm will check if the whole improved set of ranges may not be a solution of the problem by using all available simplification procedures (including the one having a "NO_3B" in CatSimp).

    As for the interval evaluation of the expressions involved in the problem some of the procedures proposed in ALIAS-Maple may use the derivatives of the expression to improve the interval evaluation. This may be computer intensive and you may avoid using these derivatives within the 3B method by setting the variable `ALIAS/Grad_3B` to a list with as many elements as number of variables. The values in this list are thresholds for the width of the variable ranges such that the derivatives are not used if at least one range width exceed the value.


    Simplification procedures for matrix

    For some procedure involving matrices there are specific simplification procedures with a different syntax: 
     
SIMP(int dimA,INTERVAL_MATRIX & A)
    This procedure shall return -1 if the set of matrices A satisfies a given property, 0 otherwise (dimA is the dimension of the matrices).

    next up previous Next: Optimization Up: ALIAS-Maple Previous: The solving procedures
  • J-P. Merlet home page
  • La page Présentation de HEPHAISTOS
  • HEPHAISTOS home page
  • La page "Présentation" de l'INRIA
  • INRIA home page

    jean-pierre merlet
    2018-07-25