## 3 Materials and Methods
### Implementing the Algorithm  
Let us begin by specifying the value used for all of the parameters involved in the algorithm. Our choices were influenced by some preliminary numerical results and the values proposed in the literature, such as [1] (Chapter 6). In our implementation, the parameters are set to:  
\[\Delta^{0} =1 \text{(initial trust-region radius)},\] \[\Delta_{s}^{0} =0.5 \text{(initial sampling radius)},\] \[\Delta_{\max} =10^{3} \text{(maximal trust-region radius)},\] \[\eta_{1} =0.1 \text{(parameter for accepting the trial point)},\] \[\eta_{2} =0.9 \text{(parameter for the trust-region radius update)},\] \[\gamma =0.5 \text{(parameter for decreasing the trust-region radius)},\] \[\gamma_{inc} =2 \text{(parameter for increasing the trust-region radius)},\] \[\epsilon_{stop} =10^{-5} \text{(parameter for verifying optimality)},\] \[\mu =1 \text{(parameter for verifying the size of the trust-region radius)}.\]  
To build an approximation of the Hessian \(H^{k}\), we compute a simplex Hessian, as defined in Definition 2. Two matrices of directions \(S\) and \(T\) must be chosen; at every iteration \(k\), the matrices \(S^{k}\) and \(T^{k}\) are set to  
\[S^{k} =\frac{\Delta_{s}^{k}}{2}\operatorname{Id}_{n},\] \[\Delta_{s}^{0} =0.5 \text{(initial sampling radius)},\]0  
Multiplying the identity matrix by \(\frac{\Delta_{s}}{2}\) to form \(S^{k}\) and \(T^{k}\) guarantees that the sampling radius for building an approximation of the Hessian is equal to \(\Delta_{s}\). Setting \(S^{k}\) and \(T^{k}\) in this fashion creates \((n+1)(n+2)/2\) distinct sample points poised for quadratic interpolation. This implies that the simplex Hessian is equal to the Hessian of the quadratic interpolation function [59]. Clearly, other matrices \(S^{k}\) and \(T^{k}\) could be used, and \(S^{k}\) does not necessarily need to be equal to \(T^{k}\), nor to always be a multiple of the identity matrix for every iteration \(k\). More details on how to choose \(S\) and \(T\) so that the resulting set of sample points is poised for quadratic interpolation can be found in [59]. To build an approximation of the gradient \(g^{k}\), the gradient of the quadratic interpolation function built using the same \((n+1)(n+2)/2\) sample points used to obtain \(H^{k}\) is simply computed. Therefore, building a model \(m^{k}\) requires \((n+1)(n+2)/2\) function evaluations. Note that the sample points utilized to build \(g^{k}\) and \(H^{k}\) may be outside of the box constraint. This is allowed in our implementation. To solve the trust-region sub-problem, the Matlab command quadprog with the algorithm trust-region-reflective is used. In theory, this method satisfies Assumption 3. Based on the preliminary numerical results, the version of Algorithm 1 that is utilized to conduct the numerical experiments in this section does not attempt to improve the accuracy of the model in Step 4. The results obtained suggest that rebuilding a new model function when \(\rho^{k}<\eta_{1}\) decreases the efficiency of our algorithm. Reducing Numerical Errors  
We next discuss two strategies that decrease the risk of numerical errors. When the sampling radius \(\Delta_{s}\) is sufficiently small, numerical errors occur while computing \(g^{k}\) and \(H^{k}\), and this can cause \(g^{k}\) and \(H^{k}\) to be very bad approximations of the gradient and Hessian at \(x^{k}\). To avoid this situation, a minimal sampling radius \(\Delta_{s\min}\) is defined. Every time thesampling radius \(\Delta_{s}^{k}\) is updated in Algorithm 1 (this may happen in Step 1 or Step 5), the rule implemented is the following:  
\[\Delta_{s}^{0} =0.5 \text{(initial sampling radius)},\]1  
In our implementation, \(\Delta_{s\text{ min}}=10^{-4}\). Numerical errors can occur with relatively large values of the sampling radius \(\Delta_{s}\). The following example illustrates this situation. This motivates our choice of setting \(\Delta_{s\text{ min}}=10^{-4}\). **Example 1**. : _Let_  
\[F(x)=f_{1}(x)\cdot f_{1}(x)=\left(0.5x^{\top} \begin{bmatrix}10&9\\ 9&10\end{bmatrix} x+