# Affect of approximating a non-differentiable function on optimisation of minimisation

I am looking at a problem of constrained minimization, where the function to be minimized contains the Heaviside function, and as such is not twice continuously differentiable.

My question is what effect would the use of a twice continuously differentiable approximation to the Heaviside function have on the accuracy and efficiency of the optimization?

the problem takes the form

$min \sum_{i}^{N} k_{i}\left [ H\left ( x_{i} \right ) -T\right ],$

$s.t. \sum_{i}^{N} k_{i}\left [ x_{i} -R\right ] \le 0.01,$

where $H\left(x \right )$ is the Heaviside function

and $x,k,T,R\in\mathbb{R}$

Would the use of the approximate Heaviside

$H\left ( x \right )\approx \frac{1}{2}+\frac{1}{2}tanh(kx)=\frac{1}{1+e^{-2kx}}$

help with the minimization and give a sufficiently accurate result

or would a Legendre polynomial expansion (similar to http://www.phys.ufl.edu/~fry/6346/legendrestep.pdf) be more successful?

This approach will need to be extended into multiple dimension before implementation, for example in 3 dimensions the minimization becomes

$\sum_{i}\sum_{j}\sum_{l}k_{1,i}k_{2,j}k_{3,l}\left [ H\left ( x_{ijk}\right )-T \right ],$

and the constraint

$\sum_{i}\sum_{j}\sum_{l} k_{1,i }k_{2,j}k_{3,l}\left [x_{ijk} - R \right] \le 0.01,$

The minimization is covered by the question constrained minimization in N dimensions, this question is regarding the effect of using an approximation to the step function in the algorithm (and ultimately the choice of algorithm to use)

• Can you say which are your optimization variables and which are coefficients whose values you know? – Wolfgang Bangerth Jan 10 '14 at 13:39
• I am favoring the approximation to the Heaviside mentioned above, as this has a fairly simple form, can be twice differentiated easily, and given that I only need H=1 for $x\ge0.01$ it is sufficently accurate. – MikeW Jan 10 '14 at 13:39
• @WolfgangBangerth, apologies, I am looking to optimize wrt k. The coefficients x, T and R are all known. I now spot the error of my notation and I should have had x and k the other way round to follow convention. – MikeW Jan 10 '14 at 13:40
• Intuitively it seems to me that the tanh approximation will be much better behaved for optimization than the polynomial expansion which is going to introduce Gibbs phenomena (oscillations) creating very misleading gradients that could screw up an optimization algorithm. – Aurelius Jan 10 '14 at 18:34
• Do you know the locations of all discontinuities? If so, would it be practical to subdivide the feasible region of your problem into subregions, and optimize over each of those subregions? – Geoff Oxberry Jan 10 '14 at 19:43

If you're optimizing with respect to the $k_{i,j}$, $k_{i,j}$ isn't an argument to the Heaviside function (for each $i, j$), and the $x_{i}$ and $T$ are known parameters, your problem is continuous and (twice continuously) differentiable with respect to $k$. Your problem still looks to be nonconvex at first glance, so deterministic global optimization would require a branch-and-bound method with convex relaxations.

If the $x_{i}$ are decision variables (with or without the $k_{i,j}$), then your problem isn't just non-differentiable: it's discontinuous. Consequently, again, it's nonconvex and requires a branch-and-bound method. Generally speaking, as Wolfgang points out, if you know locations of discontinuities, you can subdivide your problem into a number of smaller subproblems on which you may be able to use more efficient methods (read: not branch-and-bound, preferable some second-order optimization algorithm). The number of subproblems may be unacceptably large (as Wolfgang points out, in this case, $2^{n}$, where $n$ is the number of decision variables that are arguments to Heaviside functions).

If you replace the Heaviside function with the hyperbolic tangent approximation, you gain twice-differentiability, but given that the hyperbolic tangent is nonconvex over any subinterval containing 0, you may also introduce nonconvexity into your problem. As a result, I don't think you gain much by introducing that approximation, and I would only use it in the event that it is impractical to divide the problem into many subproblems over which all functions are convex, continuous, and twice-differentiable.

If the problem were just about nondifferentiability and not about discontinuities, I'd say you could try looking at subgradient methods or bundle methods. These are slower than methods that require some degree of differentability. It's worth noting that even derivative-free methods require some degree of differentiability -- they merely assume that derivative information is unavailable.

• You rightly point out that all x are known, so the result of the Heaviside function is known for all values of i. As such this is a simple pair of nonlinear functions of the same varibles k. Really I am just looking for the most efficient way of optimising the minimisation subject to the constraint and generalise this to an arbitrary number of dimensions. – MikeW Jan 11 '14 at 9:51
• If that's the case, then at this point, you're probably best off trying to use a library to implement at least a local optimization of your problem. Broadly speaking, using a second-order gradient-based optimization method is probably the best that you can do. Use a library like IPOPT, or a modeling language like GAMS or AMPL. From there, you can get some sense of what obstacles might be, and we could provide further assistance. – Geoff Oxberry Jan 11 '14 at 10:12

If the number of variables you have is modest, then you can just subdivide the domain into the half-lines/quadrants/octants/etc of the space. In each of them, the Heaviside function is constant. For example, for a single optimization variable, your original problem was $$\min_x k[H(x)-T] \\ \text{s.t.} \ k[x-R] \le c.$$ You can split this into two optimization problems: first, solve $$x_1^* = \arg \min_x -kT \\ \qquad \text{s.t.} \ k[x-R] \le c \quad \text{and} \quad x\le 0.$$ then solve $$x_2^* = \arg \min_x k[1-T] \\ \qquad \text{s.t.} \ k[x-R] \le c \quad \text{and} \quad x\ge 0.$$ Note that both of these are linear optimization problems and are, consequently rather simple to solve. It is possible that one of these problems has no solution because the feasible set is empty.

If only one of the problems has a solution, then you're done. If both have a solution, you simply have to compare the value of the objective function at these points to find which one is better.

If you have more than one variable $x_i$, then the problem becomes slightly more complicated. For example, if you have two, then you need to consider four cases of signs of the variables $x_i$. For general, $N$-dimensional problems, the number of cases to compare is $2^N$. This will quickly become impractical, but for modest numbers you should be able to deal with it. In particular, you can use the usual tricks of integer programming, e.g. branch and bound: if you have already found a solution $x^\ast$ from one of the sub-problems, then this provides you with a value of the objective function that you need to get below to find a better solution. I believe you can use this in the dual formulation of linear programs to rule out that a given one of the other subproblems has a better solution, without even solving it. In effect, you may get away with solving many fewer than $2^N$ linear programs.

• The number of variables in 1 dimenson is relatively modest (~10) but I am looking at problems where this is expanded to around 6 dimensions with up to 100 variables in each, while this is still plausable I am looking to explore other approaches, mainly for these nonlinear cases – MikeW Jan 11 '14 at 9:07
• I still believe that that may be the best you can do. Yes, the smooth approximation is something you can stick into general optimizer, but you end up with ${\cal O}(2^N)$ local minima that you may have to search through. Geoff explains this in more detail. – Wolfgang Bangerth Jan 11 '14 at 13:00