[image] [image] Ressource protégée - Inria - Institut national de recherche en sciences et technologies du numérique Accéder directement au contenu

[image] [image]

* Documentation

* FR
  Français (FR)

* Anglais (EN)
  Se connecter

* Portail HAL Inria

* Recherche
  Loading...

## Recherche avancée

Information de documents

Titres

○ Titres

○ Sous-titre

○ Titre de l'ouvrage

○ Titre du volume (Série)

○ Champ de recherche par défaut (multicritères) + texte intégral des PDF

○ Résumé

○ Texte intégral indexé des documents PDF

○ Mots-clés

○ Type de document

○ Sous-type de document

○ Tous les identifiants du document

○ Identifiant HAL du dépôt

○ Langue du document (texte)

○ Pays (Texte)

○ Ville

○ À paraître (true ou false)

Ajouter
# Mixing neural networks, continuation and symbolic computation to solve parametric systems of non linear equations

Jean-Pierre Merlet

► To cite this version:

Jean-Pierre Merlet. Mixing neural networks, continuation and symbolic computation to solve parametric systems of non linear equations. Neural Networks, 2024, 176, pp.106316. hal-04609234

HAL Id: hal-04609234

https://inria.hal.science/hal-04609234v1

Submitted on 12 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
## Auteur

### Auteur (multicritères)

○ Auteur (multicritères)

○ Auteur : Nom complet

○ Auteur : Nom de famille

○ Auteur : Prénom

○ Auteur : Complément de nom, deuxième prénom

○ Auteur : Organisme payeur

○ Auteur : IdHal (chaîne de caractères)

○ Auteur : Fonction

○ Auteur : personID (entier)

○ Auteur : Identifiant de l'organisme payeur

○ Auteur : Identifiant de la structure

○ Directeur de thèse

○ Éditeur

○ Éditeur scientifique

○ Éditeur de la série

## Ajouter

### Structure

### Structure (multicritères)

○ Structure (multicritères)

○ Structure : Acronyme

○ Structure : Nom

○ Structure : Code

○ Structure : Pays

○ Structure : Type

○ Structure : État dans le référentiel

○ Structure : Identifiant HAL de la structure

○ Structure/équipe de recherche : Pays

○ Structure/regroupement d'équipes : Pays

○ Structure/laboratoire : Pays

○ Structure/regroupement de laboratoires : Pays

○ Structure/institution : Pays

○ Structure/regroupement d'institutions : Pays

## Ajouter
# Mixing neural networks, continuation and symbolic computation to solve parametric systems of non linear equations

J-P. Merlet

INRIA Sophia-Antipolis, 2004 routes des Lucioles, 06902 Sophia-Antipolis, France,

## Abstract

We consider a square non linear parametric equations system $\mathbf{F}(\mathbf{P}, \mathbf{X}) = \mathbf{0}$ which is constituted of $n$ equations in the $n$ unknowns $\mathbf{X}$ while $\mathbf{P} = \{p_1, p_2, ..., p_m\}$ is a set of parameters that play a role in the definition of the equations $\mathbf{F}$. We are especially interested in systems that have multiple real solutions with the objective of developing a solver for obtaining **all** real solutions **exactly** (a notion that is defined in the paper) for any parameter values within a bounded region. We assume that some numerical methods has allowed us to determine the real solutions for a very limited number of fixed $\mathbf{P}$ (but not necessarily all of them) called the *initial solution set*. We then present a generic approach that allows one to design a fast solver based on the automatic construction of a limited set of multi-layer perceptrons (MLP) whose predictions will be used to obtain exact real solutions by using the Newton scheme. This approach is illustrated on two engineering problems in robotics and chemistry and it will be shown that this solver provides between 99.9 and 100 % of the solutions on verification sets that are totally independent of the sets used for learning the MLPs. In both cases (and especially for the robotics example) the solver very largely outperform others alternate methods especially if the solutions have to be obtained for multiple $\mathbf{P}$. This solver has an architecture that allows for a distributed implementation both for its design and for its use: we will therefore discuss the possibility of such an implementation on a CPU/GPU system<fnref n="1" />.

**Keywords:** parametric equations, system solving, multi-layer perceptron, Newton, robotics, MPFR

## 1. Introduction

## 2. Objectives, issues and assumption

We consider a square system of non linear parametric equations $\mathbf{F}(\mathbf{P}, \mathbf{X}) = \mathbf{0}$ which is constituted of $n$ equations in the $n$ unknowns of $\mathbf{X}$ while $\mathbf{P} = \{p_1, p_2, ..., p_{n_p}\}$ is a set of $n_p$ parameters that play a role in the definition of the system $\mathbf{F}$. We assume that for any parameter vector $\mathbf{P}$ the system $\mathbf{F}$ may have multiple real solutions (but in most cases the number of solutions cannot be determined in advance by inspecting the system), while there is no assumption on the nature of the equations in $\mathbf{F}$ except that at least one of them is non linear and all of them are at least C1.

Using neural networks to solve partial differential equations [1, 2] or other specific system [3] has attracted a lot of attention but works related to the solving of non differential equations are much more sparse. Neural networks have been used to enlarge the domain for the initial guess of the Netwon scheme that will lead to a solution [4, 5]. Goulianas and al. have proposed a specific neural network architecture that may find a solution over a finite domain for the unknowns values [6, 7] and create several such networks to look for all solutions over a possibly large search domain for the value of the unknowns.

<footnote n="1">¹This work has been partly supported by grant ANR-19-P3IA-0002
URL: Jean-Pierre.Merlet@inria.fr (J-P. Merlet)</footnote>
Autres

Champ de recherche par défaut (multicritères)

○ <u>Champ de recherche par défaut (multicritères)</u>

○

○ <u>Revues (multicritères)</u>

○ <u>Revue : Éditeur</u>

○ <u>Revue : Titre abrégé</u>

○ <u>Revue : Titre</u>

○ <u>Revue : Identifiant interne</u>

○ <u>Revue : Date de début de publication</u>

○ <u>Revue : ISSN électronique</u>

○ <u>Revue : ISSN</u>

○ <u>Revue : Éditeur</u>

○ <u>Revue : Couleur dans SHERPA/RoMEO</u>

○ <u>Revue : État dans le référentiel</u>

○

○ <u>Colloque (multicritères)</u>

○ <u>Colloque : Titre</u>

○ <u>Colloque : Organisateur</u>

○ <u>Colloque : date de début (Année)</u>

○ <u>Colloque : date de fin (Année)</u>

○ <u>Conférence invitée (oui/non)</u>

○

○ <u>Projets ANR (multicritères)</u>

○ <u>Projet ANR : Acronyme</u>

○ <u>Projet ANR : Acronyme du programme</u>

○ <u>Projet ANR : Code décision (référence)</u>

○ <u>Projet ANR : Nom</u>

○ <u>Projet ANR : Identifiant interne</u>

○ <u>Projet ANR : État dans le référentiel</u>

○

○ <u>Projets européens (multicritères)</u>

○ <u>Projet européen : Acronyme</u>

○ <u>Projet européen : Identifiant de l'appel à projet</u>

○ <u>Projet européen : Référence</u>

○ <u>Projet européen : Nom</u>

○ <u>Projet européen : Date de fin</u>

○ <u>Projet européen : Financement</u>

○ <u>Projet européen : Date de début</u>
## 2.1. Objectives

Compared to the state of the art the problem we are addressing is somewhat special:

* we are not considering a specific system but an infinite set of equations. We impose however a constraint on the parameters **P**: they should lie in a bounded region **R** defined by a set of intervals for each elements $p_i$ of **P** i.e $p_i \in [\underline{p_i}, \bar{p_i}]$

* there is no limit a priori for the value of the unknowns although we may impose such a limit if needed

* we aim at finding **all** real solutions of **F** for any **P** in **R**

* the calculated solutions will be **exact** in the sense that they can be refined at will to provide an approximation $\mathcal{A}$ of the real solution $\mathcal{S}$ such that for all unknowns we have $|\mathcal{A}_i - \mathcal{S}_i| \le 10^{-j}, j$ being an arbitrary positive integer

The contribution of this paper is a generic methodology to design a fast solver which will provide ideally all the solutions of the system for any **P** located in **R**. A key point of this methodology is the automatic determination of a limited set of multi-layers perceptrons (MLP) whose prediction will be used as a guess for the Newton method so that there is a high probability that for a given **P** the set of solutions determined over all MLPs fully cover the set of all solutions.

They are a few alternate numerical methods for solving this problem for an instance of **P** but they may be extremely computer intensive (as exemplified by the first example in the next section). A practical motivation for developing such a solver is, for example, to determine approximate minimum and maximum values for the unknowns over the bounded domain **R** for **P**. Finding the exact extrema amounts to solve a constrained optimization problem that may be intractable. In that case we may sample the domain **R** for **P**, solve the equations for the samples and retain as optimum the lowest and largest unknowns values over the whole sampling. The quality of the determination of the extrema will depend upon the number of samples and therefore motivate developing a fast solver.

There is no specific motivation for using MLP for this solver beside the facility of implementing and training them (but we will see that we take into account the structure of the problem during the training process). However using MLP raises issues that are mentioned in the next section.

## 2.2. Issues

By design an MLP takes an input (in our case a set of **P**) and provide a fixed size output prediction. As the maximal number of solutions is not known we cannot design a single MLP providing several sets of $n$ dimensional vectors each of them being a prediction for one of the solutions. Hence our approach is based on the use of several MLPs that will provide as output a vector of dimension $n$ that should be a prediction for one of the solution for a given **P**. This raises one issue: how many MLPs will be required to get all solutions as the number of MLPs is an upper bound for the number of solutions that can be obtained by the proposed solver ? We will see however that the proposed approach will automatically provide a significant number of MLPs although it is not sufficient to guarantee to obtain all solutions for all instances of **P** over the **R** domain. Furthermore this choice imposes a constraint on the learning sets: in a learning set we cannot have samples with the same input but different outputs.

## 2.3. Assumption

We assume that a numerical method such as interval analysis or continuation (see section 3.1) has allowed us to calculate the solutions (not necessarily all of them, see section 5.2) for a limited set of $l$ **P** vectors {$\mathbf{P}_1, \mathbf{P}_2, \dots, \mathbf{P}_l$} i.e for each $\mathbf{P}_j$ we will obtain $n_j$ solutions {$S_{j_1}, S_{j_2}, \dots, S_{j_{n_j}}$}. This set will be called the **initial solution set** and it is of good practice to distribute the $\mathbf{P}_j$ over all the bounded domain for **P**. Our approach will be illustrated on two examples that are presented in the next section.
○ Projet européen : État dans le référentiel

○ Projet européen : Identifiant interne

○

○ Date de production : année

○ Date de mise en ligne : année

○ Date de publication : année

○ Date d'écriture : année

○ Date de modification du dépôt : année

○ Date de dépôt : année

○ Date de publication électronique : année

○

○ Collection HAL (multicritères)

○ Collection HAL : catégorie

○ Collection HAL : Code

○ Collection HAL : Nom

○ Collection HAL : Identifiant interne

○

○ Identifiant interne du contributeur/déposant

○ Nom complet du contributeur/déposant

○

○ Domaines

○ Domaine primaire

○ Domaine racine

○ Sous-domaine niveau 1

○ Sous-domaine niveau 2

○ Sous-domaine niveau 3

○

○ Statut du document

○ Version du document

○ Type de dépôt

○ Type de document

○ ISBN

○ Numéro - référence

○ Identifiant DOI

○

○ Classification

○ Audience

○ Vulgarisation

○ Comité de lecture - texte (oui ou non)

○ Actes de colloque
## 2.4. Examples

### 2.4.1. Example 1: cable-driven parallel robot

A cable-driven parallel robot (CDPR) has a moving platform connected to the ground by $n$ cables (figure 1). At the points $A_i$ we have winches that may change the length of the cables that are attached to the platform at points $B_i$. A reference frame $(O, \mathbf{x}, \mathbf{y}, \mathbf{z})$ is defined and the coordinates of the $A_i$ in this frame are supposed to be known. Similarly a mobile frame $(G, \mathbf{x}_r, \mathbf{y}_r, \mathbf{z}_r)$ is attached to the platform and the coordinates of the $B_i$ in this frame are known. The pose $\mathbf{X}$ of the platform may be defined by the coordinates in the reference frame of a specific point $G$ attached to the platform (e.g. its center of mass) and three angles (e.g. roll, pitch, yaw) that define the orientation of the mobile frame with respect to the reference frame. For controlling a CDPR we need to solve the inverse kinematics (IK) problem (i.e

<figure><img src="image_2.png" /><figcaption>Figure 1. A cable-driven parallel robot (CDPR)</figcaption></figure>

finding the lengths $L_0$ of the cable to have the platform in a given pose) and the direct kinematics (DK) (i.e finding the pose of the platform being given the $L_0$). The kinematics equations are dependent upon the assumed behavior of the cables. If they have no mass and no elasticity, then their shape under traction is the straight line between $A_i$ and $B_i$ (such a cable is called ideal). However the assumptions of no mass, no elasticity do not hold for real cables and we have adopted a model presented in the textbook of Irvine [8] that takes into account both the cable mass and its elasticity. Irvine model is planar: the cable is attached at one of its extremities $A$ and is submitted at its other extremity $B$ to a force with $F_x, F_z$ horizontal and vertical components. In the planar frame with origin at $A$, an horizontal axis $x_i$ and vertical axis $z_i$ the coordinates of $B$ are $x_b, z_b$. The Irvine model relates $x_b, z_b$ to $F_x, F_z$ and to the cable length at rest $L_0$. We have:

$$x_b = F_x \left( \frac{L_0}{EA_0} + \frac{\sinh^{-1}\left(\frac{F_z}{F_x}\right) - \sinh^{-1}\left(\frac{-\mu g L_0 + F_z}{F_x}\right)}{\mu g} \right) \quad (1)$$

$$z_b = \frac{F_z L_0}{EA_0} - \frac{\mu g L_0^2}{2EA_0} + \frac{\sqrt{F_x^2 + F_z^2} - \sqrt{F_x^2 + (-\mu g L_0 + F_z)^2}}{\mu g} \quad (2)$$

where $E$ is the Young modulus of the cable matrix (which characterizes its elasticity), $\mu$ the linear density of the material and $A_0$ the area of a cable cross-section. It can be shown that the $x_b, z_b$ of each cable may be expressed as functions of the 6 pose parameters. Hence for a CDPR with $n$ cables the above equations provide $2n$ equations between the pose parameters and the $3n F_{x_i}, F_{z_i}, L_0$. Then we have to consider that the system is mechanically in equilibrium and the cancellation of the forces and torques acting on the platform provide 6 equations relating the pose parameters and the $2n F_{x_i}, F_{z_i}$. If we consider the DK problem the $L_0$ are known and we have to determine the 6 pose parameters and the $2n F_{x_i}, F_{z_i}$: in that case we have $2n + 6$ equations for $2n + 6$ unknowns and the system is always square whatever $n$ is. As for the IK the system is square only if $N = 6$. Our second example deal with a specific CDPR (with known coordinates for the $A_i, B_i$) having 8 cables. Hence for this parametric system the parameters are
○ Référence interne

○ Financement

○ Collaborations

Ajouter

Lancer la recherche

Recherche experte (SolR)

Recherche experte (SolR)

Lancer la recherche
Vers la recherche avancée

• Déposer

Portail INRIA.HAL.SCIENCE
Les publications, logiciels... des scientifiques Inria

• Accueil

• Publications Inria

○ Consulter les derniers dépôts

○ Par équipe de recherche

○ Par auteur

○ Rapports de recherche et techniques

○ Livres blancs

○ Rapports d'activité des équipes

○ Statistiques

• Boite à Outils

○ X2Hal : importer par lot

○ Logiciels : créer le Codemeta.json

○ Créer mon IDHal

○ Créer mon CVHal

○ Haltools : créer ma page de publications

○ Haltools : extraire des données (réservé Inria)

○ Les bonnes pratiques

■ Comment bien s'affilier ?
the 8 L₀ while the unknowns are the 6 pose parameters and the 16 Fₓᵢ, Fᵧᵢ, the system having 22 equations. Solving this system for a given instance of the L₀ (i.e finding all its solutions) is very difficult. To the best of the author knowledge only 2 methods have been proposed: continuation [9] (see section 3.1) and interval analysis [10]. The use of these methods has shown that the DK has usually several solutions (with an example with 39 solutions) but also that they have several drawbacks: both methods may miss solution(s) and they are computer intensive (several hours or even a day of computation time for large CDPR). The IK for n = 6 cables has also usually several solutions but has apparently less solution than the DK (up to now only examples with 5 solutions have been presented). Neural networks have been proposed for the DK but with the purpose of finding a single, approximate, solution [11], or for finding a solution of the IK [12]. For this example, called example 1, the initial solution set has 8 elements with a number of solution between 8 and 24 for a total of 156 solutions. They have been obtained by using a mix of interval analysis and continuation.

### 2.4.2. Example 2: countercurrent reactor

Example 2 is derived from the equations of a biochemical reactor [13]. It has n = 9 equations, 3 parameters (a, b, c) and 9 unknowns (from x₁ to x₉).

$$
\begin{align*}
j = 1 &\quad x_1 - (1 - x_1)x_{j+2} - a(4x_{j+1}) = 0 \\
j = 2 &\quad -(1 - x_1)x_{j+2} - b(1 + 4x_j) = 0 \\
j = 3 &\quad cx_1 - (1 - x_1)x_{j+2} - x_j * (1 + 4x_{j-1}) = 0 \\
3 < j < n-1 &\quad x_1x_{j-2} + (1 - x_1)x_{j+2} - x_j(1 + 4x_{j-1}) = 0 \\
j = n-1 &\quad x_1x_{j-2} + x_j(1 + 4x_{j-1}) = 0 \\
j = n &\quad x_1x_{j-2} + (1 - x_1)x_{j+2} - x_j(1 + 4x_{j-1}) = 0
\end{align*}
$$

and we are interested in finding all solutions of this parametric system when a, b, c are constrained to lie in the range [-0.5,0.5]. The initial solution set has 8 elements that have from 4 to 15 solutions that have been obtained with interval analysis, the computation time being between 3 and 5 seconds. The total number of solutions in the set is 62.

## 3. Designing a solver based only on MLPs

A MLP requires a learning set and we show in the next section a first approach to design such a set.

### 3.1. Generating learning sets

At this stage our only knowledge is the initial solution set obtained for *l* particular **P** with **P**<sub>j</sub> leading to *n*<sub>j</sub> solutions. We will generate learning sets by using as starting point each of the solution *S*<sub>jk</sub> of each **P**<sub>j</sub> (we have therefore Σ<sub>j=1</sub><sup>j=l</sup> *n*<sub>j</sub> starting points). We then select randomly an unit vector **v** in the **P** space. An iterative scheme is used where we have **P** defined at iteration *m* as **P**<sub>j</sub><sup>m</sup> = **P**<sub>j</sub><sup>m-1</sup> + ε(*m*)**v** with ε(0) = 0. Based on the implicit function theorem it is known that each of the solutions of **F**(**P**<sub>j</sub>, **X**) = **0** when **P** follows a continuous path also follows a continuous path that we will call a **branch**. The idea of the **continuation** scheme [14] is to note that at iteration 0 we know the solution of **F** so that if ε(1) is chosen small enough the solution *S*<sub>jm</sub><sup>1</sup> of **F**(**P**<sub>j</sub><sup>m</sup>, **X**) = **0** may be found by using *S*<sub>jk</sub> as initial guess for the Newton scheme. More generally at iteration *m* + 1 we set ε(*m* + 1) = ε(*m*) + α(*m* + 1) where α has a “small” value and the solution *S*<sub>jk</sub><sup>m</sup> is used as a guess for the determining *S*<sub>jk</sub><sup>m+1</sup>. However α(*m* + 1) must be chosen carefully so that the Newton scheme does not converge to a solution on another branch. Fortunately there are tools that allow one to determine α(*m* + 1) to avoid this problem (see section 4) but α(*m* + 1) may become very small when we are close to a singularity of **F** or one of the *p*<sub>i</sub> in **P**<sub>j</sub><sup>m</sup> reaches a domain limit: in that case we stop the process. At the start of the process we have stored the pair (P<sub>j</sub><sup>0</sup>, S<sub>jk</sub> = S<sub>jk</sub><sup>0</sup>) in a list *B*<sub>jk</sub> and during the process we get successive pairs (P<sub>j</sub><sup>m</sup>, S<sub>jk</sub><sup>m</sup>). If *B*<sub>l</sub> denotes the last pair in *B*<sub>jk</sub> a pair at iteration *u* will be stored in *B*<sub>jk</sub> as soon as one of the *p*<sub>i</sub><sup>u</sup> is such that |*p*<sub>i</sub><sup>u</sup> - *p*<sub>iB<sub>l</sub></sub>| is greater than a fixed threshold (a similar strategy may be used for the difference between the solutions). For example 1 we get 146 learning sets and 62 for example 2.

**Note:** In some cases a Newton scheme implemented in C with double floating point accuracy whose stopping criteria is that the residuals have an absolute value lower than a threshold may fail or may provide only a very approximate
* Quelle version puis-je déposer ?

* Choisir le bon « type de document »

* Mettre en valeur ses publications

○ Les tutos IES Inria

* Lier les notices dans HAL

* Umbrhal Inria pour l'accessibilité

* Personnaliser son CVHAL

* Créer et Personnaliser sa collection

• Documentation

○ Aide en ligne HAL

○ Aide en ligne Haltools

○ Aide en ligne X2hal

○ Tutoriels video HAL

○ Documentation API HAL

○ Ajouter des vignettes

○ Identifiant auteur IdHAL et CV

○ Compléter son dossier ORCID avec ses dépôts HAL

○ Inria et la Science Ouverte

X

X

X

Loading...

X

Le document que vous demandez est en access restreint
jusqu'au 2024-12-12 !

Revenir a la notice

Retourner à l'accueil

Contact

•

Ressources
solution. As we will be seen in section 4.1 we will also use a Maple Newton procedure which allows one to obtain the solution with an arbitrary accuracy but which is much slower then the C Newton. A general policy when using continuation is to refine the solution obtained by the C Newton by using it as initial guess for the Maple Newton, the result being the one that is stored in the sets. However in some cases there may be no convergence for the C Newton while the Maple Newton will converge. For example for example 2 a check of 522698 samples has shown that 595 samples have a solution that was not obtained by the C Newton. Hence when the continuation process stops because the C Newton has not converged we also check that the Maple Newton also do not converge otherwise the continuation goes on.

## 3.2. MLPs training and prediction quality

A MLP geometry is defined by the number of hidden layers, the number of neurons in each layer, the activation function used for each neuron, the loss function and the optimizer used to decrease the loss. We have no mean to determine in advance what should be these parameters for obtaining the best prediction. However we have determined that the training time was low (typically around one minute). We have therefore used a systematics approach with the symplifying assumptions of having the same number of neurons for each layer and to use the same activation function for all the neurons. We have tested each MLP having from 1 to 20 hidden layers with a step size of 1, from 10 to 200 neurons with a step size of 10 and any activation function among ReLU, CELU, LeakyReLU, Softplus and the sigmoid while MSE was chosen as the loss function and Adam as the optimizer. As for the learning rate we have an adaptive learning rate that starts at 1e-3 and is divided by 1.1 if the optimizer has not succeeded in decreasing the loss function after 15 iterations. The training process is stopped as soon as the learning rate reaches a value lower than 1e-6 or the loss function is lower than 1e-6. We have tested all the above combinations for each learning set and have retained as best MLPs the 10 MLPs that were exhibiting the lowest loss.

On both examples the final loss for these MLPs was extremely variable ranging from 0.3 to over 1000 and in almost all cases the quality of the prediction was very poor on the learning sets with large difference between the prediction and the exact solution. A detailed analysis of the prediction error will be presented in section 9 but in summary this black box approach does not allow to reach any of our objectives.

# 4. Hybridation of the MLPs

## 4.1. The Newton scheme

As we are looking for exact solutions a natural approach will be to rely on the Newton scheme to improve the result of the MLPs. The Newton scheme is provided with an initial guess $X_0$ for the solution and is an iterative scheme. If $J$ is the jacobian of the system $F$ the Newton scheme is defined as

$$X_{k+1} = X_k - J^{-1}(X_k)F(X_k) \quad (3)$$

where $J^{-1}$ is the inverse of $J$. The scheme is stopped if the absolute value of each component of $F(X_k)$ is lower than a fixed threshold or if $k$ (the number of iteration) is greater than a fixed value. The Newton scheme has interesting properties:

* when converging it is quadratically convergent so that it is a fast method

* on the other hand Newton may converge to a solution that is not the closest solution to $X_0$ and $X_k$ should be not too close to a singularity because of the numerical instability of $J^{-1}$

The later property raises an issue for the continuation process at it may cause a branch jump: if $\alpha$ is chosen too large the continuation process may leave the branch that was followed and jump on another branch. A fundamental theorem to better understand the behavior of the Newton scheme is the Kantorovitch theorem [15].

### Kantorovitch theorem

1. the Jacobian matrix of $F$ has an inverse $\Gamma_0$ at $X_0$ such that $||\Gamma_0|| \le A_0$

2. let $B_0$ be such that $||\Gamma_0 F(X_0)|| \le 2B_0$

3. let $U$ be a ball centered at $X_0$ with diameter $2B_0$
* Documentation

* FAQ

* API

* OAI-PMH

* AuréHAL

## Informations

* À propos

* Données personnelles

* Mentions légales

* Accessibilité

* Conformité RGAA

## Questions juridiques

* Je publie, quels sont mes droits ?

* Loi pour une République numérique

* Stratégie de non-cession des droits

* SHERPA/RoMEO

## Portails

* Portails institutionnels

* HAL

* HAL SHS

* HAL Thèses

* MédiHAL

## CCSD

* CCSD

* Episciences

* Sciencesconf
4. let $C$ be such that the Hessian matrix $\mathbf{H}$ of $\mathbf{F}$ satisfies $\|\mathbf{H}(\mathbf{X})\| \le C$ for all $\mathbf{X}$ in $U$

If $2nA_0B_0C \le 1$, then there is an unique solution of $\mathbf{F} = \mathbf{0}$ in $U$ and the Newton method used with $\mathbf{X}_0$ as estimate of the solution will converge toward this solution. We will denote the quantity $2nA_0B_0C$ the Kantorovitch criteria.

This theorem is useful for us for two purposes:

* determining a safe value of $\alpha$ in the continuation process that is used for building the learning set. As seen previously a solution $S = \mathbf{X}_0$ for $\mathbf{P}_j^m$ has been obtained and Newton is used to obtain a solution on the same branch for $\mathbf{P}_j^m + \alpha\mathbf{v}$. We select an arbitrary value for $\alpha$ and compute the value of $A_0$ and of $B_0 = \|\Gamma_0\mathbf{F}(\mathbf{X}_0)\|$. The calculation of $C$ is more intricate in general (although for example 2 matrix $\mathbf{H}$ is constant) and require the use of interval analysis. If $2nA_0B_0C \le 1$ then the chosen value of $\alpha$ is valid otherwise it is divided by a fixed number larger than 1 and the process is repeated until $\alpha$ becomes valid. Note that the same process may be used to increase the value of $\alpha$ to speed up the continuation. Other simpler method may be used to adjust $\alpha$ [14]

* for computing an exact solution. Here symbolic calculation tool like Maple that allows for extended arithmetics will be used for finding an approximation of the exact solution that differs from it by less than $10^{-a}$ where $a$ is set by the end-user. For example Maple performs its calculation with an accuracy determined by the variables Digits (default value = 10). The Newton calculation is performed iteratively. At iteration 0 we obtain a solution $\mathbf{S}_d$ with Digits set to $d$, $S_d$ is used as initial guess for a Newton where Digits is set to $d+1$, leading to a new solution $\mathbf{S}_{d+1}$. During this process the value of $B_0 = \|\Gamma_0(\mathbf{S}_{d+1})\mathbf{F}(\mathbf{S}_{d+1})\|/2$ will decrease and we stop the process if $B_0 < 10^{-a}$ and $\|\mathbf{S}_{d+1} - \mathbf{S}_d\| < 10^{-a}$. A drawback of this approach is that it is relatively slow. Initially we have mitigated this problem by using a Newton scheme implemented in $C$ whose result is used as a guess for a Maple Newton scheme with high accuracy (this will be called the classical Newton scheme) but unfortunately the $C$ procedure may miss solutions and the data have to be corrected. An interesting alternate that we have used lately is to use the extended arithmetic package MPFR [16]. We have implemented a Newton for example 2 using this package, the stopping criteria being that the residus should have an absolute value lower than $10^{-50}$ and we will see in the next section that it is extremely efficient.

## 4.2. Using Newton for hybridation

It may seem to be natural to look for a MLP whose final loss is minimal and then use the Newton scheme with as initial guess the prediction of this MLP. However by using this strategy if we check all the samples of the learning sets of example 1 we get only about 2% of the total number of solutions as exact solution and we never get more than one solution for any $\mathbf{P}$ in the learning sets. To explain this disappointing result we have first to look at the influence of the unknowns on the Kantorovitch criteria. To exemplify this point we have considered a value of $\mathbf{P}$ for which example 2 has 6 solutions. We have then added a value of 0.01 to each element of $\mathbf{P}$ and we have numerically estimated the derivatives of the Kantorovitch criteria with respect to each of the unknowns for the new system at each of the solution (table 1).

<table><caption>Table 1. Absolute value of the derivatives of the Kantorovitch criteria with respect to each unknowns and for the 6 possible solutions, ordered by increasing influence</caption><thead><tr><th>solution</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td></td><td>x<sub>1</sub> 13.751</td><td>x<sub>1</sub> 3.351</td><td>x<sub>5</sub> 9.772</td><td>x<sub>1</sub> 9.579</td><td>x<sub>7</sub> 2.454</td><td>x<sub>7</sub> 407.668</td></tr><tr><td></td><td>x<sub>5</sub> 13.751</td><td>x<sub>5</sub> 3.392</td><td>x<sub>4</sub> 9.772</td><td>x<sub>6</sub> 9.579</td><td>x<sub>6</sub> 2.454</td><td>x<sub>6</sub> 407.668</td></tr><tr><td></td><td>x<sub>3</sub> 13.751</td><td>x<sub>9</sub> 3.435</td><td>x<sub>1</sub> 9.772</td><td>x<sub>5</sub> 9.579</td><td>x<sub>9</sub> 2.468</td><td>x<sub>9</sub> 407.668</td></tr><tr><td></td><td>x<sub>2</sub> 13.751</td><td>x<sub>7</sub> 3.436</td><td>x<sub>8</sub> 9.772</td><td>x<sub>4</sub> 9.579</td><td>x<sub>8</sub> 2.504</td><td>x<sub>8</sub> 407.668</td></tr><tr><td></td><td>x<sub>7</sub> 13.751</td><td>x<sub>3</sub> 3.460</td><td>x<sub>2</sub> 9.772</td><td>x<sub>2</sub> 9.579</td><td>x<sub>5</sub> 2.506</td><td>x<sub>5</sub> 407.668</td></tr><tr><td></td><td>x<sub>4</sub> 13.751</td><td>x<sub>2</sub> 3.547</td><td>x<sub>3</sub> 9.772</td><td>x<sub>9</sub> 9.579</td><td>x<sub>2</sub> 2.525</td><td>x<sub>2</sub> 407.668</td></tr><tr><td></td><td>x<sub>6</sub> 13.751</td><td>x<sub>8</sub> 3.554</td><td>x<sub>7</sub> 9.850</td><td>x<sub>7</sub> 9.579</td><td>x<sub>1</sub> 2.530</td><td>x<sub>1</sub> 407.668</td></tr><tr><td></td><td>x<sub>8</sub> 13.896</td><td>x<sub>4</sub> 3.568</td><td>x<sub>6</sub> 9.854</td><td>x<sub>3</sub> 9.579</td><td>x<sub>4</sub> 2.598</td><td>x<sub>4</sub> 407.668</td></tr><tr><td></td><td>x<sub>9</sub> 14.953</td><td>x<sub>6</sub> 3.665</td><td>x<sub>9</sub> 10.656</td><td>x<sub>8</sub> 10.405</td><td>x<sub>3</sub> 2.681</td><td>x<sub>3</sub> 418.291</td></tr></tbody></table>
As may be seen from this table the most influential unknown and the influence of the unknowns change with the solution. This shows that the influence on the Kantorovitch criteria (and therefore on the convergence of the Newton method) is not the same for all unknowns and furthermore that this influence changes with the initial guess of the Newton scheme. On the other hand the optimizer basically try to decrease the error just by looking at the gradient of the loss function. It may thus happen that the optimizer at some iteration may obtain a large loss decrease by decreasing largely the errors on the unknowns with minor role on the value of the Kantorovitch criteria while the errors on the Kantorovitch criteria influent unknowns may increase by a small amount. This may lead to a divergence of Newton or a convergence toward a different solution. In other words a decrease of the loss does not necessary imply an increase in the number of occurences of Newton convergence toward the expected solution over the learning set.

Before proposing another strategy we have to define the *success rate* of a MLP over a learning set. If the learning set has *l* samples and the Newton scheme using as initial guess the MLP prediction converges to the expected solution in *f* cases, the the success rate (SR) is defined as 100*f*/l (consequently if the SR is 100 the MLP provides the expected solution for all the learning set samples). The new strategy is based on the calculation of the SR during the MLP training either because the loss has decreased by more than a fixed threshold (or a percentage of the loss) and/or because there has been a number *h* of decrease of the loss since the last check of the SR that is greater or equal to a threshold (typically we use 2). The final MLP will be the one that has provided the largest SR during the training. If the SR of one MLP during the training reaches 100 we store it and if this occurs three times during the training, then the training is stopped.

To show the evolution of the SR during the training we have recorded its every change over a full training session that was stopped when the loss was minimal. Figure 2 shows a typical record for example 1 (and similar data was obtained for example 2).

<figure><img src="image_3.png" /><figcaption>Figure 2. Evolution of the success rate as function of the loss during a training session</figcaption></figure>

This figure shows interesting properties of the SR. We notice first the very oscillatory behavior of the SR. Then we see on the first part of the graph that the SR decreases although the loss is decreasing. the SR is also decreasing until the loss reaches a value around 1450 at which is is minimal (33.33%). Afterward the SR exhibits an increasing trend and reaches 100 for 7 different loss around 600 but never again reaches this value while the loss decreases. At the minimal loss (value =1) the SR is about 56%. This confirms once more that the SR may not be maximal when the loss is minimal.

Regarding the solving problem when using the MLPs designed with the new strategy we get between 20 and 30%
of the exact solutions over all the learning sets for example 1 and between 30 and 100% for example 2. Furthermore for example 1 we have been able to find at most 2 different solutions for the same $\mathbf{P}$ (and only in rare cases) while we obtain up to 4 solutions for example 2. On the other hand for example 1 these MLPs has allowed to determine 2 solutions for the initial solutions set that were not known: this is the first element that shows that it is not necessary to have all solutions for this set. Overall this results show a large progress but we are still far away from our objectives.

The MLP training time will increase because of the computation of the SR during the process but may also be decreased because a SR of 100 is reached well before we come close to the minimal loss. For example 2 we have trained 226 MLPs and for 149 of them (66%) we have obtained a SR of 100 while for example 1 only about 10% of the MLPs have reached a SR of 100. Therefore a key point for the computation time of the program that computes the SR over all the samples of a learning set. This time may be decomposed into 2 components: the time for obtaining the prediction from the MLP, which is approximately constant, and the time for the phase where we use the Newton scheme to possibly obtain an exact solution and check if it corresponds to the solution provided in the learning set. If we use the classical Newton scheme on average the prediction time is only 0.24% of the total computation time. If we use the Newton based on MPFR the prediction time represents 4.36% of the total computation time as the computation time for finding and checking the solution is about 30 time faster than with classical Newton scheme.

A first interpretation of the limitation of the solver is that we do not sufficiently take into the structure of the problem, namely the interconnection between the solutions of the different learning sets. Furthermore the change of the number of solution between two close $\mathbf{P}$ is not managed. For instance for example 2 the initial solutions set has between 4 and 15 solutions but there is no member with 9, 10, 11, 12, 13 and 14 solutions in the set. The next section will an approach to manage these issues by creating samples in a structured way.

## 5. Structuring the learning sets

Let us consider again the notion of *branch* that has been introduced in section 3.1. A branch is a continuous curve $C$ in the solution space that is obtained when $\mathbf{P}$ follows a curve $C_P$ in the parameters space and to each point on this curve corresponds a usually unique point on $C$. We have seen in section 3.1 how the Newton scheme may be used to obtain points on $C$ for a given $C_P$ provided that no *singularity* of $\mathbf{F}$ is encountered on the branch as Newton cannot be used at a singularity. Without going into the details [17] a singularity may be *isolated* and *removable*. For example such a singularity will occur at the intersection of two branches while after the singularity the branches will continue and there are methods to recover the branches [14]. . A singularity may also be non-isolated for example when two branches join at a singularity or when a branch reaches a natural boundary where $\mathbf{F}$ cannot be analytically continued: in both cases if this situation occurs for any path $C_P$, then the number of solutions after the singularity will be lower than the number of solutions before the singularity. Encountering a singularity during the continuation is not an issue. But we may have an opposite situation when working with real solution. If for a given $\mathbf{P} = \mathbf{P}_1$ we have $n$ branches and compute all the solutions with an alternate method such as interval analysis for $\mathbf{P}_1 + \epsilon$ it may happen that we have more than $n$ solutions. Indeed between $\mathbf{P}_1$ and $\mathbf{P}_1 + \epsilon$ the system may have exited from a singularity such as a natural boundary so that real solution(s) may pop-up and we have no mean to detect such a situation during a continuation process. Fortunately we will see that the MLPs we will have created may solve this issue.

Singularities play also a role in the interconnection between different solutions obtained for the same $\mathbf{P}$. Assume that we have defined a curve $C_P$ such that this curve contains several occurences of the same $\mathbf{P} = \mathbf{P}_s$. For each of this occurence $\mathbf{P}_{s_j}$ we have a corresponding point $S_{s_j}$ on $C$ and it may happen that we have $S_{s_k} \neq S_{s_j}$ for some of them. In that case we have found a singularity-free connecting path between $S_{s_k}$ and $S_{s_j}$ (we will show examples of this behavior for example 1).

Interconnection and singularity should play a role when building the learning sets: this may be managed by using the notion of *aspect*.

### 5.1. Aspect

We consider the Cartesian product $\mathcal{U}$ of the parameters space and of the solution space. A point in $\mathcal{U}$ is therefore defined as a pair $(\mathbf{P}_1, \mathbf{S}_1)$ where $\mathbf{S}_1$ is a solution for $\mathbf{F}(\mathbf{P}_1, \mathbf{X}) = \mathbf{0}$. An *aspect* $\mathcal{A}$ is a region of $\mathcal{U}$ such that any two points $(\mathbf{P}_1, \mathbf{S}_1)$, $(\mathbf{P}_2, \mathbf{S}_2)$ of $\mathcal{A}$ are such that there is a continuous path in the parameters space that joins $\mathbf{P}_1$ to $\mathbf{P}_2$ leading to a continuous, singularity-free path in the solution space that connect $\mathbf{S}_1$ to $\mathbf{S}_2$.
An aspect is therefore a variety such that all the solutions are interconnected and calculating aspect is usually a complex task. We will therefore use a weaker concept, the *linear aspect*. Two points ($\mathbf{P}_1, \mathbf{S}_1$), ($\mathbf{P}_2, \mathbf{S}_2$) of $\mathcal{U}$ belong to a linear aspect if the linear path in the parameters space defined by $\mathbf{P} = \mathbf{P}_1 + \lambda(\mathbf{P}_2 - \mathbf{P}_1)$ with $\lambda \in [0, 1]$ leads to a path in the solution space that connect $\mathbf{S}_1$ to $\mathbf{S}_2$. A linear aspect is therefore a subset of a more general aspect.

## 5.2. The solution graph

Aspects involving the points of the initial solutions set may be graphically represented with a planar graph. For that purpose we define a *level* as a particular value for $\mathbf{P}$ and a level will be represented by an horizontal line. For a given level we have various solutions that will be represented by nodes on the level line. For the initial solutions set we have a limited number $l$ of level that may ordered arbitrarily as $\mathcal{L}_1, \mathcal{L}_2, \dots, \mathcal{L}_l$ that are associated to the parameters $\mathbf{P}_1, \dots, \mathbf{P}_l$ (this ordering will be called the *initial ordering*). We will consider in turn each of the solutions $\mathbf{S}_{1k}$ obtained for level $\mathcal{L}_1$ and we will consider a linear path in the parameters space that leads from level 1 to level 2. Using continuation we follow a path in the solution space that starts from $\mathbf{S}_{1k}$ and may possibly lead to solution $\mathbf{S}_{2m}$ at level 2. If this is the case we draw an edge between the corresponding nodes on the level 1 and 2. From level 1 we repeat the process for all levels larger than 1 up to level $l$. More generally the edge construction will be applied from level $m$, examining the interconnection of the solution at level $m$ to all level $L_k$ with $m < k \le l$. Using this process we obtain a graph, the *initial graph*, that connect nodes at different levels. Note that this graph is associated to the ordering that has been chosen, all the edges connecting a level to an upper level.

We will repeat the process for all possible ordering of the initial solutions set. However the new graph will be transposed to be represented in the same ordering than the initial graph meaning that now we have edges connecting one level to a lower level, meaning that the edges are *directed*. We thus obtain an *ordered graph* that will be called the **solutions graph**.

Note that in the process we may discover solutions for some level(s) that were not present in the initial solution set. Discovering such a solution imposes to update the graph but at some point no new solutions will be discovered and the graph is stable. The edges of such a graph may be represented by a set of integer quadruplets $(l_1, k_1, l_2, k_2)$ that means that solution $k_1$ at level $l_1$ is connected to solution $k_2$ at level $l_2$. For example 1 we have 307 quadruplets in the graph and 402 for example 2.

By construction if $(l_1, k_1, l_2, k_2)$ is present in the graph, then $(l_2, k_2, l_1, k_1)$ is also part of the graph. Any set of nodes that are connected by edges are part of the same linear aspect. Figure 3 shows the graph obtained for the initial solutions set of example 1.

<figure><img src="image_4.png" /><figcaption>Figure 3. The solution graph for example 1</figcaption></figure>

We note some interesting properties of this graph:

* some nodes are *isolated* meaning that a solution at a given level has no connection to any other solution at any level. There are 27 such nodes for example 1 and none for example 2
* we have path in the node that connects a node at a given level to another node at the same level. For example the path shown in red on the figure includes the quadruplets (6,1,7,1), (7,1,8,1) and (6,3,8,1) meaning that solution 1 and 3 of level 6 are part of the same linear aspect. The path (6,3,8,1), (8,1,7,1), (7,1,6,1) therefore allows to move from solution 3 of level 6 to solution 1 of the same level without crossing any singularity of F

The solution graph reflects the structure of the system. It contains the information on all the members of the initial solution set that belong to the same linear aspect while the isolated nodes indicates that the initial solution set has linear aspects with a single point within the set. The solution graph is the basic element that will allow us to construct an initial solver, called the *primary solver*, that will be refined incrementally if needed.

# 6. The primary solver

The first step for establishing the primary solver is to create learning sets based on the solution graph.

## 6.1. Solution graph and learning sets

### 6.1.1. Isolated nodes and blocs

The first learning sets we will establish are based on the isolated nodes of the solution graph if any. These nodes are the start points of the learning set that are established using the method described in section 3.1 using random unit vectors $v$ in the parameter space and moving from the $P$ of the starting point along these vectors until a singularity is encountered and storing samples on the way.

We then consider the nodes that are connected by edges. As mentioned previously on a path connecting nodes we may have different nodes at the same level i.e. different solutions for the same set of $P$. Such a situation shall not occur in a learning set for a MLP as the training process assumes that there is a single output for a given input. The solution graph is constituted of $m$ quadruplets $Q_j = (l_1^j, k_1^j, l_2^j, k_2^j)$ and we will deduce from the graph several ordered lists of quadruplets, called a *bloc*. A bloc if constituted of a sequence of quadruplet $\{Q_j, Q_{j+1}, \dots, Q_{j+m}\}$ such that for any $k > 0$ we have $l_1^{j+k} = l_2^{j+k-1}$ (in other words $Q_m, Q_{m+1}$ are in a bloc because in the graph we have an edge that joins solution $k_2^m$ at level $l_2^m$ to solution $k_1^{m+1}$ at level $l_1^{m+1}$). Elements of the bloc must satisfy the *inclusion condition*: a quadruplet may be added to a bloc only if its $l_2$ component is different from any of the $l_1, l_2$ components of the quadruplet already stored in the bloc, the purpose being not having in a bloc two solutions at the same level.

To construct the bloc we will consider in sequence each of the quadruplets in the graph and assume that it is the first element of a bloc. If this quadruplet is $Q_j = (l_1^j, k_1^j, l_2^j, k_2^j)$ we look in the graph for a quadruplet $Q_m$ such that $l_1^m = l_2^j, k_1^m = k_2^j$ satisfying the inclusion condition. If such a quadruplet is found it is added to the bloc and we repeat the process for $Q_m$ until no more quadruplet can be added to the bloc. After completing this process for all quadruplets in the graph we consider the bloc that has the largest number of quadruplets and retain it as *final bloc*. We then repeat the process of finding the largest bloc except that no quadruplet in a final bloc will be considered as first element of a bloc. The process stops when all quadruplets of the graph belong to a final bloc. For instance for example 2 we end up with 75 blocs, each bloc having from 2 to 8 quadruplets.

### 6.1.2. Creating learning sets for the blocs

A bloc establish a connection between nodes $\{n_1, n_2, \dots, n_u\}$ at different level $\{P_1, P_2, \dots, P_u\}$. To obtain a learning set for a bloc we use in sequence a continuation process between two successive nodes with $P = P_k + \lambda(P_{k+1} - P_k)$ and we store sample in the learning set according to the following rules:

* $P_1, S_{n_1}$ is stored in the learning set

* a sample $(P_k, S_{n_k})$ is stored in the learning set if the absolute value of the difference between one element of $P_k$ and the last $P$ stored in the learning set is larger than a given threshold

* a similar rule is applied for $S_{n_k}$

* a sample is not stored if $P_k$ is equal to any $P$ already stored in the learning set
In figure 4 the elements of the bloc ($\mathbf{P}_k, \mathbf{S}_{nk}$) are represented by black square while the samples obtained from the continuation process are represented by black circles.

We aim at having a fixed minimal size $H$ for the learning set (typically $H = 1000$). If after having completed the first part of the process the number of samples in the learning set is less than $H$ we add samples by considering as starting points the samples $(\mathbf{P}a_1, \mathbf{S}a_1), (\mathbf{P}a_2, \mathbf{S}a_2), \dots (\mathbf{P}a_v, \mathbf{S}a_v)$ stored in the learning set for the continuation process of section 3.1 using random vector $v$. These added samples are represented by white circle in figure 4. At the opposite if a block has a size that is much larger than $H$ the learning set is split in different learning sets that have a size that is roughly equal to $H$. By construction the samples of a bloc belong to the same linear aspect.

<figure><img src="image_5.png" /><figcaption>Figure 4. The samples that are stored in the learning set for a given bloc</figcaption></figure>

## 6.2. Training and checking the learning sets

The MLP training is performed using the success rate approach described in section 4.2 and after completing the process for each learning set we obtain a list $L_{mlp}$ of MLPs. The use of the predictions of the MLPs in $L_{mlp}$ combined with a Newton scheme applied on the samples of a specific learning shall ideally provides the exact solution for all samples in the learning set. In that case we will say that we have a *full coverage* of the learning set by the set of MLPs in $L_{mlp}$. Our aim is evidently to obtain a full coverage of all the learning sets.

Some learning sets will have full coverage with a single MLP as during the training this MLP has exhibited a success rate of 100. For some other learning sets the dedicated MLP may have reached a success rate lower than 100, meaning that for some samples the exact solution has not been found. Still full coverage of these learning sets may be obtained as the solution(s) missed by the dedicated MLP may be found with other MLPs in $L_{mlp}$. However it may happen that for some sample(s) in the learning sets the use of all the MLPs does not lead to the exact solution: these samples are stored in a list $\mathcal{L}_{nf}$. We consider the first element $\mathcal{L}_{nf}^1$ of $\mathcal{L}_{nf}$ to create a specific learning set for this sample, using $\mathcal{L}_{nf}^1$ as starting point for the procedure described in section 3.1. This learning set is used to create a MLP $M_1$ that is added to $L_{mlp}$, provided that this MLP has a success rate of 100 (which ensure that $\mathcal{L}_{nf}^1$ is found). Otherwise we create a new learning set using lower storing thresholds, this ensuring that we will find a set leading to a $M_1$ that has a success rate of 100. We also use $M_1$ on the other elements of $\mathcal{L}_{nf}$ as it may also lead to other solution(s). The solutions that have been found are removed from $\mathcal{L}_{nf}$ and we repeat the process until $\mathcal{L}_{nf}$ is empty. Using this strategy the MLPs stored in $L_{mlp}$ ensure a full coverage of all the learning sets. The MLP(s) that have been added during this step are called *specific MLPs*.

For our test examples we obtain full coverage of the learning sets for example 1 with 1154 MLPs (with 16 specific one) and with 160 MLPs for example 2 (with 7 specific one). Figure 5 shows the repartition of a, b in the learning sets for example 2 around some specific value for c while figure 6 presents a 3D view of the a, b, c.

## 6.3. Trimming the number of MLPs

Among all the MLPs in $L_{mlp}$ we may have redundant MLP(s). For example the solutions found by using $MLP_A$ may be included in the union of the solution found by using $MLP_B$ and $MLP_C$ and hence $MLP_A$ is redundant and may be eliminated from $L_{mlp}$. We propose in this section 2 approaches that may lead to a lower number of MLPs while keeping full coverage of all the learning sets.
<figure><img src="image_6.png" /><figcaption>Figure 5. Repartition of a, b in the various learning sets around $c = -0.2, 0, 0.2$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6. Perspective and side view of the repartition of a, b, c in the various learning sets</figcaption></figure>

### 6.3.1. Approach 1

Regarding redundancy the aim should be to find the lowest number of MLPs that are necessary to obtain the full coverage of all the learning sets but this is a difficult problem. We aim at the more modest objective: to obtain a list $L_{mlp}^u$ with full coverage and a number of elements lower than $L_{mlp}$.

Approach 1 is iterative and prior to the first iteration we set $L_{mlp}^u$ to be empty, we store all the samples of all the learning sets in a list $S$ and we determine the set of solutions $S_{mlp}^i$ that are found by each of the MLP (denoted $MLP_i$) in the list $L_{mlp}$. Then we determine if they are $S_{mlp}^i$ that are found only by a single MLP. If this is the case we add the MLPs to $L_{mlp}^u$, remove them from $L_{mlp}$ and remove all their $S_{mlp}^i$ from the $S_{mlp}^i$ of the MLPs still in $L_{mlp}$ and from $S$. We will use an iterative process and at each iteration:

* we determine the MLP $m$ that has the largest number of elements in $S_{mlp}^m$. This MLP is stored in $L_{mlp}^u$ and is removed from $L_{mlp}$

* all elements of $S_{mlp}^m$ are removed from $S$ and from all the $S_{mlp}^i$ resulting from the MLPs still in $L_{mlp}$

The iterations stop when $S$ is empty and we have in $L_{mlp}^u$ a set of MLPs that provide full coverage for all elements in the learning sets. This set has possibly less elements than the initial $L_{mlp}$. Note that the only computer intensive step of this approach is the determination of the initial set of solutions $S_{mlp}^i$ that are found by each of the MLP, the remaining steps being very fast.

For our test examples the result of this trimming approach is very good: for example 1 the number of necessary MLPs decreases from 1154 to 138 and from 160 to 87 for example 2.

### 6.3.2. Approach 2

The drawback of approach 1 is that computing all solutions found by each MLP may be extremely intensive if we have a very large number of samples.
To propose an alternate iterative method we first introduce some notations: $N$ is the total number of samples in $\mathcal{S}$, $M$ is the total number of MLPs, $n_i$ is a flag that is set to 1 when the solution for sample $i$ has been found, $\hat{S}_j$ will be the set of solutions found for MLP $j$ while $\hat{S}$ will be the union of the $\hat{S}_j$, $m_j$ is a flag that is set to 1 when MLP $j$ has been used. As for approach 1 $L_{mlp}^u$ will be the possibly trimmed set of MLPs that will be the result of the algorithm and if the empty set when starting the algorithm. At the first iteration we have all $n_i$, $m_j$ set to 0 while $\hat{S}_j = \emptyset$ for all MLP and consequently $\hat{S} = \emptyset$.

At the difference of approach 1 we will not determine the $\hat{S}_j$ over all the samples in the learning sets but we will work iteratively in the main iteration process and select only a subset $\mathcal{V}$ of $\mathcal{S}$ with around 10 000 samples. These samples will be obtained by picking the elements $\{S_1, S_{1+u}, S_{1+2u}, \dots\}$ of $\mathcal{S}$ where $S_1$ is the first element of $\mathcal{S}$ and $u = \lfloor N/10000 \rfloor + 1$. At each main iteration we first compute the solution set $\hat{S}_j$ for all MLPs in $L_{mlp}$. We then consider the union of the $\hat{S}_j$ for the MLPs that have $m_j = 1$ (and thus belong to $L_{mlp}^u$ as we will see later on) and we decrement $N$ by the number of elements in this union. If this union fully cover the solution of $\mathcal{V}$, then the set $L_{mlp}^u$ fully cover $\mathcal{V}$ and we may proceed to the next iteration in the main loop. Otherwise we remove the elements of this union of all the $\hat{S}_j$ of the MLP that have $m_j = 0$. We will then use a second iterative process:

* we determine the MLP $m$ with $m_m = 0$ that has the largest number of elements in its solution set. This MLP is stored in $L_{mlp}^u$ and $m_m$ is set to 1.

* we set $n_i$ to 1 for all elements of $\hat{S}_m$ and we remove these elements from the $\hat{S}_j$ of the remaining MLPs that have $m_j = 0$. We then decrement $N$ by the number of elements in $\hat{S}_m$

* we exit if the $\hat{S}_j$ of all MLPs in $L_{mlp}$ with $m_j = 0$ are empty

At the end of this second iteration process we have a full coverage of $\mathcal{V}$ by the MLPs in $L_{mlp}^u$. If $L_{mlp}^u = L_{mlp}$ (condition G), then no trimming is possible. Otherwise we perform an iteration of the main iteration process. If condition G is not satisfied after having examined all the samples of $\mathcal{S}$, then $L_{mlp}^u$ has less elements than $L_{mlp}$. Compared to approach 1 this algorithm may determine that no trimming is possible in a faster way as it may not be necessary to check all samples before reaching conition G.

For example we have considered 598069 samples for example 2 that have been obtained by considering 106 MLPs. Approach 2 stops after having considered 523106 samples as condition G is reached. The three MLPs that have provided the largest number of solutions have respectively found 55284 (10.568413 %), 44463 (8.499807 %) and 42713 (8.165267 %) solutions.

# 7. Checking the primary solver

In the previous section we have established a primary solver using MLPs that ensure full coverage of the learning sets that has been used for creating the MLPs. These learning sets were providing a single solution for a given set of parameters. We will now establish verification sets that will provide all the solutions for a given set of parameters.

## 7.1. Generating verification sets

For establishing a valid verification set we must ensure than the **P** of the sample in the verification set are different from the **P** in the learning sets. The most basic way to create a verifications is to use the same continuation process than for the learning sets using as starting points the initial solutions set and with a random choice for the vectors **v** with the following differences:

* we follow simultanously all the solutions from the starting points and stop the process if Newton does not converge for one solution.

* we use a different seed than for the learning sets for initializing the random process of the **v** vector

* the **P** of the stored sample must be at least at a threshold distance from the **P** of the learning sets (for example the starting points are not stored in the verification set)
If the size of the verification set is not sufficient we may repeat the process still using as starting point a member of the initial solution set but with a different **v** or by selecting randomly one of the sample of the verification that will be used as starting point.

There is a drawback in this approach as the number of solutions in all the samples deduced from a specific member *M* of the initial solutions set will have the same number of solution than *M*. For instance for example 2 we have in the initial solutions set starting points with 4, 5, 6, 7, 9 and 15 solutions and therefore none with less than 4 solutions, with 8,10,11,12,13,14 solutions or for more than 15 solutions. But we will see that this drawback may be overcome.

## 7.2. Verification set checking with the primary solver

The process of checking a verification set is to use all the MLPs in $L_{mlp}$ to check if they allow to determine all solution for the samples in the verification set (i.e we have a full coverage of the verification set). Two situation may occur during this verification:

* a new solution for a given sample is discovered by the solver. The verification set is updated by adding the solution to the sample of the verification set

* one (or more) solution for a sample is not found by the solver

For instance we have considered a first verification set for example 2 with 1030 samples for a total of 7980 solutions and have used the primary solver based on 87 MLPs. This solver finds correctly all solutions for 1011 samples (98.15%) and a total of 7961 solutions (99.76%) while 313 new solutions are found. Adding these new solutions to the verification set we have 8293 solutions and 8274 that are found by the solver (99.77%), meaning that the solver has missed 19 solutions. After adding the new solutions we have samples in the verification set with any number of solutions between 4 and 15 except for 14. For example 1 the resultats are similar with 99.25% of the samples in the verification set for which all the solutions are found and overall 99.95% of the solutions are found.

Hence the result is already quite good but we have neither a full coverage of the verification set nor samples in the verification set that cover all number of solutions between the two extrema. The next section will present an approach to improve the result.

# 8. Creating an improved solver

As seen previously we have not a full coverage of the verification set as solutions are not found. The next section will explain how we may obtain a full coverage of the verification set.

## 8.1. Improving the verification set coverage

As some solution(s) of the sample(s) in the verification set are not found we will use a similar strategy than for the one used for obtaining a full coverage for the learning sets (section 6.2). We select one of the solution $S_1$ that has not been found by the solver as starting point to create a specific learning set and create a new MLP from this learning set in such a way that the new MLP will find the missing solution.

Applying this strategy to example 2 we use now 88 MLPs and the new solver determines 8277/8293 solutions (99.807%) meaning that the new MLP has provided not only the solution $S_1$ but also two other solutions as the number of missed solution is now 16. No new solution have been found by the solver and we have still only 1011 samples for which all solutions have been found.

The process may be repeated for each of the solutions that the solver does not find. Table 2 shows for example 2 at each iteration the number of MLPs used for the verification, the percentage of found solutions, the number of samples for which all solutions are found (full solution). It may be seen that the iteration stops for 92 MLPs as we have a full coverage of the verification set and evidently of the learning sets. However we are not fully satisfied by the new solver as the verification set has no sample with 14 solutions and the next section will propose a strategy to correct that.
<table><caption>Table 2. Improvement of the verification set coverage as a function of the number of MLPs. One MLP is added at each iteration of the improvement algorithm.</caption><tr><td>number of MLPs</td><td>88</td><td>89</td><td>90</td><td>91</td><td>92</td></tr><tr><td>% of found solution</td><td>99.807066</td><td>99.819125</td><td>99.855300</td><td>99.867358</td><td>100.000000</td></tr><tr><td>full solution</td><td>1011</td><td>1014</td><td>1015</td><td>1018</td><td>1030</td></tr></table>

## 8.2. Improving the new solver with new verification sets

As seen in section 7.1 the creation of a verification set is basically a random operation that is modified by changing the initial seed of the random generator that is used to calculate the v vectors. As generating a verification set and calculating its coverage are relatively fast operations (see section 10.2) it may be of interest to check the new solver with different verifications sets. We have created a program to automatically generate an arbitrary number of random verification sets and, if needed, create the necessary MLPs that are needed to obtain a full coverage of each verification sets. Typically figure 7 shows for example 2 how checking between 2 and 78 new verification sets influences the number of required MLPs to have full coverage.

<figure><img src="image_8.png" /><figcaption>Figure 7. Number of MLPs that are required to have full coverage of all verification sets as a function of the number of checked verification sets</figcaption></figure>

When we have checked 18 new verification sets (using 102 MLPs) the latest set has samples with between 4 and 16 solutions, all numbers in between being represented in the samples. coverage. After having considered 78 new verification sets we have a total of 64978 samples in the verification sets and 519925 solutions with the following repartition in term of number of samples for a given number of solutions in the samples: 4 (7049), 5 (8744), 6 (7515), 7 (16120), 8 (557), 9 (9707), 10 (3229), 11 (3294), 12 (442), 13 (188), 14 (8), 15 (8124), 16 (1). Let $n_{mlp}^k$ denotes the number of MLPs that are required to have full coverage of all verification set when having k new verification sets. It must be noted that when considering a new verification set while we have already checked k verification sets we have either a stability in the number of MLPs (i.e. $n_{mlp}^{k+1} = n_{mlp}^k$ or an increase only by one (i.e. $n_{mlp}^{k+1} = n_{mlp}^k + 1$). Furthermore in all cases when considering a new verification set k + 1 we have to create a new MLP only because we have a single solution that has not been found with the MLPs used to manage the k previous verification sets. In practice this means that if we have used only the 95 MLPs required to manage 2 new sets to check the 78 new verification sets, then we will have missed 11 solutions which represent 0.00211 % of the total number of solutions. However it may be seen that the number of MLPs required to have full coverage all the verification sets is not stable even after having checked 78 new verification sets. This implies that even having considered a large number of verification sets we still cannot guarantee to be able to find all the solutions for any particular parameters values.

As we cannot consider a too large number of verification sets to improve the solver we need a rule to stop considering new verification sets. Let us assume that after considering k new verification sets we have obtained samples whose number of solutions is in the range [$n_{min}$, $n_{max}$] but with no sample for some values in this range. A rule of the thumb derived from example 2 is to stop considering new verification sets when the number of solutions for all samples lies in the range [$n'_{min}$, $n'_{max}$] with samples that have a number solution equal to any value in this range. For instance for example 2 all number of solutions in the range [4,16] appears when the number of MLPs is 102, obtained
after checking 18 verification sets. This means that if we have stopped the process at this point we will have missed 4 solutions that appears when using 78 new verification sets i.e. 0.000769 % of the total number of solutions, which is quite acceptable.

## 9. MLP prediction quality

We have motivated the hybridation mechanism used for the training of the MLPs and for exploiting their predictions because the Newton method allows for getting exact solutions and in a significant number of cases to reduce the training time. But it is of interest to better quantify the quality of the MLPs predictions. For that purpose we have retained the 98 MLPs that were used for example 2 by one of the solvers. We have considered the latest verification set used to build the solver, which has a total of 8133 solutions. We have then recorded all the MLPs predictions for each sample of the verification set. For each sample we have determined the MLPs that has provided the closest prediction to the exact solution with the additional constraint that the distance between the solution and the prediction shall not be larger than 5. Among the 8133 solutions 7215 predictions satisfy the distance criteria (88.71%). The bounds for the corresponding variables $x_i$ at the solution are given in table 3.

<table><caption>Table 3. Bounds for the unknowns in the verification set</caption><thead><tr><th>x<sub>1</sub></th><th>x<sub>2</sub></th><th>x<sub>3</sub></th><th>x<sub>4</sub></th><th>x<sub>5</sub></th><th>x<sub>6</sub></th><th>x<sub>7</sub></th><th>x<sub>8</sub></th><th>x<sub>9</sub></th></tr></thead><tbody><tr><td>[-7.3, 7.3]</td><td>[-3.3, 2.5]</td><td>[-4.2, 1.3]</td><td>[-4.7, 3.7]</td><td>[-3.8, 6.2]</td><td>[-10.4,11.0]</td><td>[-10.1,38.0]</td><td>[-21.3,459.6]</td><td>[-17.5,52.2]</td></tr></tbody></table>

We have first considered the value of the 9 equations for each of this prediction (meaning that we have checked 64935 equations). Table 4 presents statistics about the equation values. Overall the mean value of the absolute value of

<table><caption>Table 4. Statistical indicators on the absolute value of the equations for the 7215 predictions.</caption><thead><tr><th>Equation</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th></tr></thead><tbody><tr><th>Mean</th><td>0.3398</td><td>0.3224</td><td>0.2598</td><td>0.2543</td><td>0.4259</td><td>1.2638</td><td>0.5090</td><td>1.0739</td><td>5.6274</td></tr><tr><th>Variance</th><td>0.3762</td><td>0.3855</td><td>0.1935</td><td>0.2578</td><td>1.8634</td><td>34.3310</td><td>4.2177</td><td>16.8348</td><td>2952.64</td></tr><tr><th>Min</th><td>1.4e-05</td><td>6.6e-06</td><td>1.4e-05</td><td>5.6e-06</td><td>1.3e-06</td><td>2.2e-05</td><td>1.6e-05</td><td>1.4e-05</td><td>2.9e-05</td></tr><tr><th>Max</th><td>5.1499</td><td>4.2956</td><td>4.1759</td><td>9.4985</td><td>18.5974</td><td>280.0568</td><td>31.1026</td><td>125.4207</td><td>1369.18</td></tr></tbody></table>

the equations is 1.119 (variance 337.18, min 0.1371e-5, max 1369.18). We have 7186 equations that have an absolute value less than 0.01 (11.064 %) and 33293 equations with an absolute value less than 0.1 (51.27%).

If we look at each of the 7125 solutions of the samples and consider the prediction that is the closest to the solution we may determine how many predictions we have such that all the absolute values of the 9 equations are lower than a given threshold. Table 5 provides this number (and the % of the total number of solutions it represent) as a function of various threshold. It must be noted that we have no prediction such that all equation values are lower than 0.0025. Figure 8 shows an histogram of the equation value distribution for the 6276 cases where the equation value is $\le 10$.

<table><caption>Table 5. Number of solution predictions such that all 9 absolute values of the equations are lower than a given threshold</caption><thead><tr><th>Threshold</th><th>0.01</th><th>0.1</th><th>0.2</th><th>0.3</th><th>0.4</th><th>0.5</th><th>0.6</th><th>0.7</th><th>0.88</th><th>0.9</th></tr></thead><tbody><tr><th>solution</th><td>95</td><td>1595</td><td>2280</td><td>2736</td><td>3398</td><td>3742</td><td>3916</td><td>4036</td><td>4198</td><td>4452</td></tr><tr><th>in %</th><td>1.3%</td><td>22.1%</td><td>31.6%</td><td>37.92%</td><td>47.09%</td><td>51.8%</td><td>54.2%</td><td>55.9%</td><td>58.1%</td><td>61.7%</td></tr></tbody></table>

With 98.7% of the predictions having at least one of the 9 equations having an absolute value larger than 0.01 it can be seen from the equations viewpoint the predictions perform rather poorly.

We may now look at the distances between the best prediction and the solution. Their mean value is 0.9736 (variance: 1.2228, min: 0.19 e-2, max: 4.998). Figure 8 shows the distribution of the number of predictions as function of the distance. Among the 7215 data set we have 72 cases where the distance is less than 0.01 (0.997%).
and 1039 for which the distance is less than 0.1 (14.4%). But as seen in table 3 there are large variations in the value for the unknowns at a solution. We have therefore investigated the absolute value $d_i$ of the difference between each unknown value $s_{x_i}$ at a solution and its prediction $p_{x_i}$ over the whole 7215 cases. The mean value of the $d_i$ is 0.9736

<table><thead><tr><th>variable</th><th>d<sub>1</sub></th><th>d<sub>2</sub></th><th>d<sub>3</sub></th><th>d<sub>4</sub></th><th>d<sub>5</sub></th><th>d<sub>6</sub></th><th>d<sub>7</sub></th><th>d<sub>8</sub></th><th>d<sub>9</sub></th></tr></thead><tbody><tr><td>Mean</td><td>0.192</td><td>0.215</td><td>0.14</td><td>0.176</td><td>0.133</td><td>0.242</td><td>0.164</td><td>0.4</td><td>0.279</td></tr><tr><td>Variance</td><td>0.096</td><td>0.1937</td><td>0.051</td><td>0.1188</td><td>0.0368</td><td>0.2789</td><td>0.1237</td><td>0.4497</td><td>0.3423</td></tr><tr><td>Min (e-5)</td><td>0.49</td><td>1.5</td><td>0.6</td><td>0.16</td><td>0.98</td><td>0.0052</td><td>1.7</td><td>0.44</td><td>1.4</td></tr><tr><td>Max</td><td>3.36</td><td>3.11</td><td>2.636</td><td>3.55</td><td>1.995</td><td>4.1</td><td>4.233</td><td>4.96</td><td>4.363</td></tr></tbody></table>

and if we consider $s_{x_i} - p_{x_i}$ instead of the absolute value we get a mean value of -0.00149. We have also calculated $|100(s_{x_i} - p_{x_i})/s_{x_i}|$ when $|s_{x_i}| \ge 0.1$ to get the error in term of percentage: we get a mean value of 36.46% with a minimum of 0.309e-3 % and a maximum of 2430.3 %. As seen from these results the predictions of the MLPs do not satisfy the minimal requirement for the solver (less than 1% of error on the unknowns) in spite of the structuration of the data (without this structuration the mean error for example 1 is around 200%). However it appears that they are in general sufficiently close to the solution to be used as an initial guess for the Newton method.

<figure><img src="image_9.png" /><figcaption>Figure 8. Left: histogram of the number of sample as a function of the maximal absolute value among the equations of F. Right: histogram of the distances between prediction and solutions.</figcaption></figure>

# 10. Software and computation time

## 10.1. Solver implementation

We have implemented a generic framework based on Linux scripts, Maple, C and Python programs for developing a solver whatever F may be. The starting points are

* an initial solution set provided by the end-user

* a *data file* that provides parameters for constructing the solver (such as the number of layers, neurons and activation functions for the MLPs training, the parameters and unknowns names, the box for the parameters, the storing thresholds, ...)

* a *Maple program* that provides the equations
The C and Python programs that are needed are automatically created based on the data file. The parameters for the MLPs have to be determined experimentally during preliminary trials as they may depend upon the problem at hand. After this step the learning sets and MLPs are created automatically to provide the primary solver (section 6). Then verification sets are created automatically (section 7.1) and improving the solver by using additional verification sets (section 8) can be performed incrementally until the end-user estimates that the solver results are satisfactory.

## 10.2. Computation time for establishing the solver

The calculations for establishing the solver are performed on a recent desktop computer. The Newton scheme that has been used for all the phases has been the *classical Newton scheme* (section 4.1). We have timed the various steps for the creation of the solvers for example 2. Creating the solution graph with 402 nodes takes 44 minutes. Establishing the 75 branches from the solution graph is neglectible and it takes a maximum of 2 minutes to establish the learning set for each branch (total time 3h14mn). For the training we have considered 2 possible values for the number of layers and two possible activation functions. The training time for each learning set has a mean value of 18 minutes (min: 1mn, max: 230 mn). The total time for establishing the initial MLPs was 32 hours and 33 minutes (total time: 35h47mn). To have full coverage of the learning sets we have to introduce 28 additional learning sets that are calculated in 3 hours and 20 minutes (total time: 39h07mn). Trimming the MLPs requires 1 hour and 11 minutes (total time: 40h18mn) so that obtaining the primary solver takes about 40 hours.

Now we may move to the verification sets: the mean time value for creating and checking a verification set is 17 minutes and we have considered a total of 106 verification sets. This has required about 30 hours so that the total time to calculate the solver for example 2 is about 70 hours. For example 1 this time is about 80 hours as the equations are more intricate and require more MLPs. A major part of this time is devoted to the MLPs training (about 50% of the total time) partly because of the calculation of the SR with the classical Newton scheme (this scheme having also an important role on the computation time of other steps e.g. establishing the learning and verification sets). Hence a major mean for reducing the computation time will be to use the much faster MPFR Newton scheme. For instance we have implemented this scheme for example 2 and tested it for establishing the learning sets: it appears that using the Newton MPFR is 60% faster than the classical Newton scheme while providing the same accuracy. Even a better result is obtained for creating the solution graph: compared to the 44 minutes it takes with the Maple Newton we have the same result in 3 minutes with the MPFR Newton. Unfortunately for keeping the framework automated we must have a translation program that transforms any equations into an efficient MPFR code. Such a program is not available currently but developing it is a major a priority for this research.

## 10.3. Implementation with GPU

The current framework is sequential while numerous steps of the calculation may benefit from a distributed implementation and we are currently working to implement the solver calculation on an IA dedicated processor with multiple GPUs (we use a JETSON Nano whose low energy consumption is interesting for autonomous robots). This is not an easy task as we will exemplify on a specific step in the construction of the solver namely the calculation of the SR of a MLP during the training. As mentioned in section 4.2 the SR calculation if required to obtain an efficient MLP. It consists in getting predictions for the whole learning set and running a Newton scheme with each of them to check the number of time the Newton scheme converge to the exact solution provided in the learning set. Currently on a CPU the training is performed iteratively by a Python program and every 3 iteration it checks the value of the SR for the current MLP. For calculating the SR the predictions are obtained through a Python program that uses a Pytorch model of the MLP, then a classical Newton scheme is used to compute the SR. A first step was to avoid using Python for obtaining the predictions. The Pytorch model was transformed in a C program so that the SR calculation is fully performed by a C program except the Newton scheme that mixes C and Maple (but we have see that the use MPFR allows to obtain a high accuracy C Newton running in a faster way). Managing a learning set with 50 000 samples requires about 150 seconds, the samples being processed by batch of 100 samples. If we use a pure C Newton this time is reduced to 13s (but in many cases this Newton is not accurate enough to converge toward the exact solution). We have tried to use 2 GPUs, one managing the predictions while the second one is managing a pure C Newton but the computation time is 46 s because it appears that loading the model on the GPU is a lengthy process. Consequently we have run a multi-threaded program on the CPU for calculating the predictions on batches. As soon as a batch has been processed the predictions are copied to the memory of a GPU that will run the Newton scheme while the CPU
process another batch. In this process several GPUs are involved and a synchronisation mechanism on the CPU checks if some of them have completed their job, retrieve their results (in this case only the number of exact solutions that has been found by the Newton scheme) and free them for another job. With this strategy managing 50 000 samples requires only 5 seconds. But this strategy has a drawback: as MPFR has not been ported to GPU we cannot use it to obtain the solutions with high accuracy.

In conclusion although a distributed implementation is clearly a mean to decrease drastically the computation time for establishing the solver the use of GPU will provide a limited gain compared to a classical grid.

### 10.4. Computation time for using the solver

The process that uses the solver for determining all solutions for a given **P** (for short the solver) is written in C. It has several variants but all of them share the first task to launch a Python program that uses in sequence the Pytorch model of each MLP to produce a solution prediction for the current **P** and write them in a file. The C program then read the prediction file and from this point start the variants:

* variant 1: uses a C Newton scheme in sequence on each prediction and sort out a list of possible solutions as the same solution may be found by different MLPs. These solutions are made exact by using a Maple program that uses them as initial guess and runs the high accuracy Maple Newton

* variant 2: uses the prediction as initial guess in a Maple program that runs the high accuracy Maple Newton

* variant 3: same as variant 2 but run the MPFR solver to get exact solutions

* variant 4: same as variant 1 but uses the MPFR Newton to make the solutions exact

Variant 1 and 4 share a potential problem: the C-Newton stops its calculation when the residual are all lower than a threshold (typically 1e-13) but this may not be sufficient to distinguish two very close solutions. Variant 3 may experiment the same problem but it has a much lower threshold (1e-40) while variant 2 is the only one that is safe from this viewpoint.

It may be noted that if the solutions have to be found for multiple **P** the described process is not efficient as we will have the overhead to launch the Python prediction program for each of the **P**. Consequently the process has been implemented to deal with multiple **P** so that the Python program is launched only once as a background task. It may also be expected that variant 3/4 should be faster than variant 2/1 as these variants use the MPFR-Newton instead of the Maple-Newton. As the MPFR-Newton has been implemented only for example 2 table 6 presents the computation time for different number of samples.

<table><caption>Table 6. Mean computation time in seconds for getting all solutions for a given **P** according to the variants and the number of samples that are considered (example 2)</caption><thead><tr><th>sample↓ variant →</th><th>1</th><th>2</th><th>3</th><th>4</th></tr></thead><tbody><tr><td>1</td><td>0.9</td><td>12.41</td><td>0.963</td><td>0.975</td></tr><tr><td>100</td><td>0.15</td><td>5.038945</td><td>0.155</td><td>0.03084</td></tr><tr><td>200</td><td>0.1379</td><td>5.28</td><td>0.158097</td><td>0.023792</td></tr><tr><td>300</td><td>0.1319</td><td>5.613</td><td>0.1651</td><td>0.021977</td></tr><tr><td>400</td><td>0.135986</td><td>6.227732</td><td>0.209614</td><td>0.022202</td></tr></tbody></table>

As seen in the table variants 1,3,4 exhibit approximately the same computation time when dealing with only one **P** while variant 2 exhibits a large computation time. For multiple samples the general trend is that the computation time per sample decreases up to a certain number of samples and somewhat increases for a larger number of samples: according to the variants this number is 300 (variant 1), 100 (variant 2), 100 (variant 3), 300 (variant 4). The best variants are 4, 1, 3, 2 while variant 4 is about 5 times faster than variant 1, 7 times than variant 3 and 160 times than variant 2. But these data are valid only for example 2 and may change for another case. This analysis has not been performed for example 1 as we have not yet a MPFR-Newton for this case but for variant 1, 100 samples, the mean computation for a sample is 1.5 seconds: this time is much larger than for example 2 because we have a
larger number of MLPs, a larger number of equations (22 instead of 9) and because the number of solutions is also larger. Nevertheless the current implementation of variant 4 allows one to obtain multiple solutions for example 2 in an average time which is roughly 100 times faster than the method we have used to determine the initial solution set. For example 1 although we have not used a MPFR Newton the gain is already very large compared to the currently known alternate methods.

Regarding a distributed CPU/GPU implementation our experience described in section 10.3 lead us to believe that we may somewhat increase the gain for a large number of samples by using the Pytorch MLPs translated in C and running in a multithread program on the CPU while the GPU will perform the C-Newton calculation whose result will be sent back to the CPU for refinement by a MPFR-Newton.

## 11. Conclusion

We have presented in this paper a generic approach based on MLPs to provide in a fast way exact solutions of non-linear parametric system with multiple solutions. They are very limited constraints on the system, namely that the **F** must be at least C1 and that it is necessary to be able to find the system solutions, called initial solution sets, for a very limited number of specific values for the parameters (and we have shown that is was not required to find all solutions for these cases). Building learning set (constituted of one solution for a given **P**) or a verification set (all solutions for a given **P**) is easy although it may be necessary to consider extended arithmetics with MPFR as good choice. Although various network parameters have been tested the prediction quality of the MLP is relatively poor. Therefore the first step to improve the solving is to complement the classical loss function by creating MLPs based on the number of time the Newton method allows to determine exact solutions for a given learning set. The second step was to create learning sets based on the structure of the system with the concept of aspect. This has allowed to improve the quality of the MLP predictions and to reduce the number of necessary MLPs. Still some solutions in the learning set were not found and we have proposed an incremental strategy for adding MLPs in order to reach full coverage of the learning sets. We have then proposed a trimming algorithm that allows one to keep a full coverage of the learning sets with a lower number of MLPs. We then moved to verification sets and have shown that considering several such sets allows one to fill the gap in the coverage to end up with a fast solver that will miss very few solutions. This process has been illustrated on two examples: example 1 is especially illustrative as it is a very complex problem and the developed solver largely outperforms existing methods.

As for possible improvements the use of MLP may be discussed. We have looked at Physics-informed neural networks (PINNs) by developing networks whose loss function is computed from the weighted values of the equations thereby allowing to get rid of the learning sets. The results were disappointing as the training time is very large while the quality of the predictions remains poor. We have used weights for the equations in order to produce several MLPs that may provide different solutions for the same **P** but only in a few cases have we got 2 different solutions.

Another possibility may be to change the loss function of the MLP which is currently being built based on the difference between the predictions and the sample exact solution: we may, for example, use the value of the Kantorovitch criteria at the prediction with the purpose of minimizing the average Kantorovitch criteria over the learning sets. Such a loss function will favor the convergence of the Newton scheme toward the expected solution.

The concept of aspect has helped structuring the learning sets but is not completely satisfactory as an aspect may include samples that have the same **P** but different solutions, a situation that must not occur in the learning sets. Furthermore the aspects do not convey a change in the number of solutions for different **P**. For that purpose we may restrict aspect to *discriminant varieties*, i.e. the regions in the parameter space such that for all **P** in the same region we have the same number of solutions [18] and therefore computing aspect within the same discriminant variety may reflect more the structure of the problem. For polynomial systems of reasonable size it is possible to calculate explicitly the discriminant varieties but this not so in the general case. However we are confronted to a weaker problem which is to determine learning sets with samples in the same variety. We are currently investigating this approach.

Eventually we have to consider a distributed implementation of the software first for the design of the solver and then for decreasing the computation time of its exploitation. Some elements for such an implementation has been proposed but it remains a work in progress.
## References

[1] J. Blechschmidt, O. G. Ernst, Three ways to solve partial differential equations with neural networks — a review, GAMM-Mitteilungen 44 (2) (2021) e202100006.

[2] A. Sacchetti, et al., Neural networks to solve partial differential equations: A comparison with finite elements, IEEE Access 10 (2022) 32271–32279.

[3] Y. Zhang, D. Jiang, J. Wang, A recurrent neural network for solving Sylvester equation with time-varying coefficients, IEEE Transactions on Neural Networks 13 (5) (2002) 1053–1063.

[4] G. Li, Z. Zeng, A neural network algorithm for solving nonlinear equation systems, in: Int. Conf. on Computational intelligence and Security, Suzhou, December, 13-17, 2008.

[5] D. Mishra, P. Kaira, Modified hopfield neural network approach for solving nonlinear algebraic equations, Engineering Letters 14 (1) (2007) 135–142.

[6] K. Goulianas, et al., A back propagation-type neural network architecture for solving the nxn nonlinear algebraic algebraic system of equations, Advances in Pure Mathematics 6 (2016) 455–480.

[7] K. Goulianas, et al., Solving polynomial systems using a fast adaptive back propagation-type neural network algorithm, European Journal of Applied Mathematics 29 (2) (2018) 301–337.

[8] H. M. Irvine, Cable Structures, MIT Press, 1981.

[9] J.-P. Merlet, A generic numerical continuation scheme for solving the direct kinematics of cable-driven parallel robot with deformable cables, in: IEEE Int. Conf. on Intelligent Robots and Systems (IROS), Daejeon, October, 9-14, 2016.

[10] J.-P. Merlet, The forward kinematics of cable-driven parallel robots with sagging cables, in: 2nd Int. Conf. on cable-driven parallel robots (CableCon), Duisburg, August, 24-27, 2014, pp. 3–16.

[11] A. Ghasemimi, M. Eghtesad, M. Farid, Neural network solution for forward kinematics problem of cable robots, J. of Intelligent and Robotic Systems 60 (2010) 201–215.

[12] I. Chawla, et al., Neural network-based inverse kineto-static analysis of cable-driven parallel robot considering cable mass and elasticity, in: 5th Int. Conf. on cable-driven parallel robots (CableCon), virtual, July, 7-9, 2021.

[13] L. Luksan, J. Vleck, Sparse and partially separable test problems for unconstrained and equality constrained optimization, Tech. Rep. 767, Institute of Computer Science, Academy of Sciences of the Czech Republic (January 1999).

[14] E. Allgower, Numerical continuation methods, Springer-Verlag, 1990.

[15] R. Tapia, The Kantorovitch theorem for Newton's method, American Mathematic Monthly 78 (1.ea) (1971) 389–392.

[16] L. Fousse, G. Hanrot, V. Lefèvre, P. Pélissier, P. Zimmermann, MPFR: A multiple-precision binary floating-point library with correct rounding, ACM Trans. Math. Softw. 33 (2).

[17] V. Arnold, et al., Singularity theory, Springer-Verlag, 1998.

[18] D. Lazard, F. Rouillier, Solving parametric polynomial systems, Journal of Symbolic Computation 42 (6) (2007) 636–667.