10
$\begingroup$

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):

    Mathematica graphics

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

    Mathematica graphics

    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}$.

Improve this question
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – Mark L. Stone
    Commented Jan 2, 2017 at 15:08
  • $\begingroup$ @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. $\endgroup$
    – Szabolcs
    Commented Jan 2, 2017 at 15:21
  • $\begingroup$ I look forward to seeing it, $\endgroup$
    – Mark L. Stone
    Commented Jan 2, 2017 at 16:05

1 Answer 1

Reset to default
1
$\begingroup$

Since there is a single root and there are no constraints, you might have luck posing it as an optimization problem: minimize the sum (along each dimension) of the squares of your original function.

Classical Optimization methods will likely fail, but heuristic methods such as genetic algorithms or CME-ES (covariant etc matrix adaptation - evolutionary strategy) might work.

Improve this answer
$\endgroup$
4
  • $\begingroup$ That is indeed the approach to go. I would particularly look at the SPSA algorithm that was developed specifically for your purpose and is quite robust. $\endgroup$
    – Wolfgang Bangerth
    Commented Jan 3, 2017 at 14:25
  • 3
    $\begingroup$ The OP mentions that the function is very expensive to evaluate (applying Monte Carlo simulation for a function evaluation). Doesn't that pose a very big problem for genetic algorithms and other evolutionary algorithms? They are "trivially parallel" (and is MC usually too) so massive parallel computing could be possible but are they the best way to go here? $\endgroup$
    – GertVdE
    Commented Jan 3, 2017 at 22:10
  • $\begingroup$ @WolfgangBangerth Thank you, as you say it sounds like the right solution. I will look at SPSA. $\endgroup$
    – Szabolcs
    Commented Jan 4, 2017 at 12:01
  • 1
    $\begingroup$ Regarding expensive function evaluations: It is true that genetic algorithms and related heuristic methods tend to require a larger number of function evaluations than traditional methods. The benefit is that heuristic methods can often solve problems that 1) would otherwise require a problem-specific method or 2) would fail due to numerical problems. For this example, it is likely that traditional methods would have trouble because of the stochastic nature of the objective function and the small gradients along some dimensions. SPSA looks like a great candidate method for this problem. $\endgroup$
    – MattKelly
    Commented Jan 4, 2017 at 13:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.