# Solving a difficult system of equations numerically

I have a system of $n$ non-linear equations that I want to solve numerically:

$$\mathbf{f}(\mathbf{x})=\mathbf{a}$$ $$\mathbf{f}=(f_1,\dots,f_n)\quad\mathbf{x}=(x_1,\dots,x_n)$$

This system has a number of characteristics that makes it particularly difficult to handle. I am looking for ideas on how to deal with the system more effectively.

Why is the system difficult?

• The functions are similar to this one (but of course in multiple dimensions):

They have flat plateaus separated by a region of smooth change. In 2D, you can imagine something like this for one $f_i$:

Generally, each $f_i$ has two plateaus separated by smooth change around an $n-1$ dimensional hyperplane.

Functions like this are difficult to handle with Newton-like methods because the derivative is effectively zero on the plateaus. In multiple dimensions I cannot easily find a region where none of the $f_i$ have a plateau—if I could that would solve the problem. The bisection method works well for $n=1$, but it does not generalize well to multiple dimensions.

• The functions are very slow to compute. I am looking for a method that can get a reasonable approximation of the root in as few iterations as possible.

• The functions are calculated with a Monte Carlo method. This means that each time they are calculated, I get a slightly different random value. Derivatives are difficult to estimate. Once we are close enough to the root, the noise will start to dominate, and it is necessary to use averaging to increase precision. Ideally it should be possible to generalize the method to an equivalent stochastic approximation version (e.g., Newton → Robbins-Monro).

• The system is high-dimensional. $n$ can be as large as 10-20. When $n=2$, an effective method would probably be the following: try to follow the contours defined by $f_1(x_1,x_2) = 0$ and $f_2(x_1,x_2)=0$ and see where they intersect. It is not clear how this would generalize to high dimensions.

What else do I know about the system?

• There is precisely one root (from theoretical results).

• I know the value of $f_i$ on the plateaus (let's say it's 0 and 1 for any $i$).

• $f_i$ has a special relationship to $x_i$:  $f_i(\dots, x_i, \dots)$ changes monotonically from 1 to 0 as $x_i$ goes from $-\infty$ to $\infty$. This is true for any fixed value of the other $x_{j\ne i}$.

• Do you know lower and upper bounds on all the variables, within which the solution must lie? The tighter those bounds, the better. Can you give a deterministic example, in as high a dimension as you want, which illustrates your plateaus and difficulties, but does not require Monte Carlo simulation and does not have random errors in the functions (bonus points if derivatives can be computed)? The purpose of such a deterministic example is to understand the difficulties of the problem, not to say that the Monte Carlo evaluation will not be used in the ultimate solution of your actual problem. Jan 2 '17 at 15:08
• @MarkL.Stone Bounds: I don't know them. But I can guess. The guesses would have to be pretty broad to be confident that they are correct. Example: I will come up with an example and will edit the question tomorrow. I don't have a much clearer picture about the true form of $f$ than what I described here, so my first example might not be truly representative of the real problem. But I'll put together something made of Fermi functions (sigmoids), and I'll try to make it have as many of the difficulties of the real problem as I can. Jan 2 '17 at 15:21
• I look forward to seeing it, Jan 2 '17 at 16:05