## nonlinear optimization functions

Librairie C++ de calcul formel/ C++ symbolic computation library

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### nonlinear optimization functions

Hello Bernard,

I have some useful code for Giac that I want to share. I have attached a source file containing few rather basic but useful functions for nonlinear optimization, inspired by the Maple functions of the same name.
1. minimize - minimizes a (multivariate) continuous function on a bounded area. In fact, it solves constrained minimization problems with equalities and inequalities as constraints.
2. maximize - reversed 'minimize'
3. extrema - finds extrema of a (multivariate) differentiable function subject to equality constraints. It uses standard methods of classifying critical points but also examines higher derivatives in both univariate and multivariate cases, thus being able to detect all kinds of saddle points, for example. It is an improvement over the existing function 'extrema' in Giac.
4. minimax - finds minimax polynomial approximating the given continuous univariate function on a certain segment. It is a simple but working implementation of Remez method.
5. tpsolve - transportation problem solver. Despite the fact that transportation problem is in fact integer programming problem which can be solved with 'lpsolve', I thought that it would be useful to have a separate function using specialized method (in this implementation, MODI method is combined with northwest corner method for initial feasible solution).

What do you think about adding this code to Giac source? Together with function 'lpsolve', it makes a decent set of optimization tools.
Pièces jointes
optimization.zip
(18.68 Kio) Téléchargé 122 fois

parisse
Messages : 5139
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: nonlinear optimization functions

Hi,

Great work, thanks!
It compiles nicely but your second example takes forever
extrema(x^3-2x^2+y^2+z^2-2x*y+x*z-y*z+3z,[x,y,z])
I traced it with gdb, it fails while calling minimize line 906 (I did not look why right now).
But here the hessian is [[-4,-2,1],[-2,2,-1],[1,-1,2]], numeric eigenvalues can be easily computed and one can see it's a saddle point. Another option (a little bit more complicated to implement, but exact) would be to call gauss(a2q(H,vlam),vlam) and look at the sign of the coefficients of the squares. Computing minor determinants is probably slower than computing numeric eigenvalues of a symmetric matrix and does not work every time. ker could be called to make sure the hessian is not degenerate.

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

You're right, using eigenvalues with hessian is simpler and faster. But I'm not sure if I'm doing it right. Bordered hessian is computed with respect to variables [x1,x2,...,xn,L1,L2,...,Lm] (in that order) where L1,..,Lm are Lagrange multipliers. Now i look only for first n eigenvalues: if they're all positive -> minimum, if they're all negative -> maximum, if they have mixed signs -> saddle point, if one of them is zero -> undecided. But that gives (1/3,1/3,1/3) as a saddle point of function x*y*z under condition x+y+z=1, but it is in fact local maximum, as one can check by plotting f(x,y)=x*y*(1-x-y). What I'm doing wrong?

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

See the comment below! It contains updated code.
Pièces jointes
optimization-fix1.zip
(18.37 Kio) Téléchargé 129 fois
Dernière modification par lukamar le sam. sept. 02, 2017 12:14 am, modifié 1 fois.

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

I made some changes to 'extreme' function. For problems with no constraints, standard eigenvalue test you proposed is used. For constrained problems, I had to stick with computing determinants of the bordered hessian (i added kernel check as you suggested), as I am not sure how to apply eigenvalues here (see my previous post for an example). There are only n-m determinants to compute, where n is number of variables and m number of constraints. Do you know how to detect saddle points with bordered hessian (for constrained problems)? I haven't figured out that yet.

Also, I deactivated higher derivative test for constrained problems, as I think I haven't applied it properly, since it was given only for unconstrained problems. Of course, Lagrangian holds the objective function as well as constraints, so it seems Lagrange multipliers method turns the problem into an unconstrained one at first, but in fact it can't be treated as such. If it was, critical points would be tested with respect to variables and multipliers together, which may shadow their real nature.

I made one improvement to 'minimax' function. Now an option 'limit=<posint>' can be added as argument, which limits number of iterations of Remez algorithm. That is useful for slow-converging or divergent cases when n is high (due to rounding errors, I think).
Pièces jointes
optimization-fix1.zip
(18.51 Kio) Téléchargé 136 fois

parisse
Messages : 5139
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: nonlinear optimization functions

Sorry, I should have made clear that my post was valid for unconstrained optimization problems only. I don't have experience with constrained problems.

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

One more change: I put test for degeneracy (ker(H)==[]) before computing eigenvalues too, so that they don't have to be computed if one of them will surely be equal to zero (if ker(H)!=[]).
Pièces jointes
optimization-fix2.zip
(18.47 Kio) Téléchargé 129 fois

parisse
Messages : 5139
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: nonlinear optimization functions

Actually one can conclude if there are strictly positive and strictly negative eigenvalues even if dim(ker(H))>0. But perhaps it's a good idea to keep the test as is, because if you are computing numeric eigenvalues, then you must define what strictly means...
By the way, if H does not contain symbolic expressions, it's probably a good idea to evalf it before computing eigenvalues, otherwise the system might return rootofs...

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

Here's probably the final version of the code. Now everything works. For the function 'extrema', I've dropped the kernel test and applied evalf to H before computing its eigenvalues and minors. I've succeeded to detect saddle points during the bordered hessian test. I also fixed few minor bugs and wrote the documentation entries for cascmd_en.tex. The first one (tpsolve.tex) contains a subsection about 'tpsolve' and the other (nonlinear_optimization.tex) contains a section 'Nonlinear optimization' with entries for 'minimize', 'maximize', 'extrema' and 'minimax' functions. The entry about 'tpsolve' should be put in Linear Programmation section, as it is a specialized LP solver.
Pièces jointes
optimization-fix3.zip
(22.63 Kio) Téléchargé 142 fois

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

I have one more update to the code. In the line 876 evalf is used, which may prevent confusion when testing ==0 in some cases. Also, I added one more function 'thiele' for rational interpolation which uses Thiele interpolated continued fraction based on inverse differences. It was quite simple to implement, so i thought it may come handy... Maple has a similar function 'ThieleInterpolation' and I made this one compatible to it.

Here is the documentation entry about 'thiele' function to be put somewhere in cascmd_en.tex. Of course, if you decide so!

Code : Tout sélectionner

\subsection{Rational interpolation : {\tt thiele}}\index{thiele}
\noindent {\tt thiele} takes as the first argument a matrix {\tt data} of type $n\times 2$ where that $i$-th row holds coordinates $x$ and $y$ of $i$-th point, respectively. The second argument is {\tt v}, which may be an identifier, number or any symbolic expression. Function returns $R(v)$ where $R$ is the rational interpolant. Instead of a single matrix {\tt data}, two vectors $\mathbf{x}=(x_1,x_2,\dots,x_n)$ and $\mathbf{y}=(y_1,y_2,\dots,y_n)$ may be given (in this case, {\tt v} is given as the third argument).

This method computes Thiele interpolated continued fraction based on the concept of reciprocal differences.

It is not guaranteed that $R$ is continuous, i.e.~it may have singularities in the shortest segment which contains all components of $\mathbf{x}$.

\subsubsection{Examples}
\noindent Input :
\begin{center}{\tt thiele([[1,3],[2,4],[4,5],[5,8]],x)}\end{center}
Output :
\begin{center}{\tt (19*x\verb|^|2-45*x-154)/(18*x-78)}\end{center}
Input :
\begin{center}{\tt thiele([1,2,a],[3,4,5],3)}\end{center}
Output :
\begin{center}{\tt (13*a-29)/(3*a-7)}\end{center}

In the following example, data is obtained by sampling the function $f(x)=(1-x^4)\,\mathrm{e}^{1-x^3}$.

\noindent Input :
\begin{center}
{\tt data\_x:=[-1,-0.75,-0.5,-0.25,0,}\\
{\tt 0.25,0.5,0.75,1,1.25,1.5,1.75,2];}
{\tt data\_y:=[0.0,2.83341735599,2.88770329586,}\\
{\tt 2.75030303645,2.71828182846,2.66568510781,}\\
{\tt 2.24894558809,1.21863761951,0.0,-0.555711613283,}
{\tt -0.377871362418,-0.107135851128,-0.0136782294833];}
{\tt thiele(data\_x,data\_y,x)}
\end{center}
Output :
\begin{center}
{\tt (-1.55286115659*x\verb|^|6+5.87298387514*x\verb|^|5-5.4439152812*x\verb|^|4}\\
{\tt +1.68655817708*x\verb|^|3-2.40784868317*x\verb|^|2-7.55954205222*x}\\
{\tt +9.40462512097)/(x\verb|^|6-1.24295718965*x\verb|^|5-1.33526268624*x\verb|^|4}\\
{\tt +4.03629272425*x\verb|^|3-0.885419321*x\verb|^|2-2.77913222418*x}\\
{\tt +3.45976823393)}
\end{center}
Pièces jointes
optimization-fix4.zip
(19.95 Kio) Téléchargé 136 fois
Dernière modification par lukamar le dim. sept. 03, 2017 11:49 pm, modifié 5 fois.

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

Forgot to add lines 1790-1792 in optimization.cc. Just replace the file in fix4 with this one.
Pièces jointes
optimization-fix5.zip
(19.55 Kio) Téléchargé 135 fois

parisse
Messages : 5139
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: nonlinear optimization functions

Thanks! I'm adding thiele documentation after Lagrange and spline.

parisse
Messages : 5139
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: nonlinear optimization functions

https://www-fourier.ujf-grenoble.fr/~pa ... casfr.html
The new commands are not highlighted and you may need to clear your browser cache.

lukamar
Messages : 297
Inscription : ven. juin 30, 2017 9:55 am
Localisation : Zagreb, Croatia

### Re: nonlinear optimization functions

I fixed couple of bugs.

1. The function '_minimize' was returning some meaningless results in cases when 'location' was not given because of applying '_feuille' at the end.

2. Functions 'minimize', 'maximize' and 'extrema' must make some assumptions in certain cases, namely when the variables are entered with bounds, such as 'x=a..b'. Up to now, these assumptions affected the given symbols, remaining active after function call and causing (seemingly) wrong output from other, subsequently called functions. I fixed this by introducing a temporary symbol (purged before use) for each variable to make the required assumptions. The result from each of the above functions should not depend on any of these symbols anyway.

Also, I have implemented checking for singularities of the interpolant in 'thiele' by using 'sturmab'. If singularities are detected, a warning message containing their approximate locations is printed.
Pièces jointes
optimization-fix6.zip
(19.73 Kio) Téléchargé 126 fois

parisse
Messages : 5139
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

### Re: nonlinear optimization functions

Your data for sturmab will be approx most of the time? I probably need to make sure it is exactified before Euclide algorithm is called, otherwise roundoff errors might happen.