We illustrate several robust pooling operations. Denote by $\phi$ the penalty function. Then robust pooling finds the solution to the optimization problem
$$ y \in \text{arg min}_u \sum_{i=1}^{n} \phi(u - x_i; \alpha) $$where $\alpha$ is a parameter of the penalty function.
# define various penalty functions
import numpy as np
phi = {
'Quadratic': lambda z, alpha: 0.5 * np.power(z, 2.0),
'Huber': lambda z, alpha: np.where(np.abs(z) <= alpha, 0.5 * np.power(z, 2.0), alpha * np.abs(z) - 0.5 * alpha * alpha),
'Pseudo-Huber': lambda z, alpha: np.power(alpha, 2.0) * (np.sqrt(1.0 + np.power(z, 2.0) / np.power(alpha, 2.0)) - 1.0),
'Welsch': lambda z, alpha: 1.0 - np.exp(-0.5 * np.power(z / alpha, 2.0)),
'Trunc. Quad.': lambda z, alpha: np.minimum(0.5 * np.power(z, 2.0), 0.5 * alpha * alpha)
}
# plot the penalty functions
%matplotlib inline
import matplotlib.pyplot as plt
alpha = 1.0
z = np.linspace(-5.0, 5.0, num=101)
for name in sorted(phi.keys()):
plt.plot(z, phi[name](z, alpha))
plt.title("$\phi(z, a)$")
plt.legend(sorted(phi.keys()))
plt.show()
By plotting the sum of two penalty functions centered at $x_1 = -2\alpha$ and $x_2 = 2\alpha$, respectively, we can observe some of their properties: such as the convexity of Huber, the strong convexity of pseudo-Huber and the non-convexity of Welsch.
# plot the sum of two penalty functions
alpha = 1.0
delta = 2.0 * alpha
u = np.linspace(-5.0, 5.0, num=101)
for name in sorted(phi.keys()):
plt.plot(u, phi[name](u - delta, alpha) + phi[name](u + delta, alpha))
plt.title("$f(x, u) = \phi(u - x_1, a) + \phi(u - x_2, a)$")
plt.legend(sorted(phi.keys()))
plt.show()
# plot the sum of two penalty functions scaled
for name in sorted(phi.keys()):
f = phi[name](u - delta, alpha) + phi[name](u + delta, alpha)
plt.plot(u, f / (np.max(f) - np.min(f)))
plt.title("$f(x, u) = \kappa \phi(u - x_1, a) + \kappa \phi(u - x_2, a)$")
plt.legend(sorted(phi.keys()), loc='lower right')
plt.show()
We now explore how the solutions of the sum of two penalty functions vary as a function of its inputs. We fix $x_1$ to zero and varying $x_2$. Our objective (with $\alpha = 1$) is $$ f(x, u) = \phi(u) + \phi(u - x) $$ for scalar $x$ and $u$.
As we vary $x$ we see that the non-convex penalty functions can produce multiple isolated solutions, while the strongly convex quadratic and pseudo-Huber always produce a single solution. Note, however, that for a large number of points (sum of $n >> 2$ penalty functions) the effect from changing a single point would not be as pronounced as we demonstrate here for just two points.
# prepare data for animation and plots
def find_local_minima(v):
"""Finds the local minima of a one-dimensional curve and returns their index. In the
case of a plateau (such as for Huber) returns the boundary points of the plateau."""
return [i for i in range(1, len(v)-1) if ((v[i] < v[i-1]) and (v[i] <= v[i+1])) or ((v[i] <= v[i-1]) and (v[i] < v[i+1]))]
x = np.linspace(-4.0, 4.0, num=100)
u = np.linspace(-5.0, 5.0, num=501)
y = {}
for name in sorted(phi.keys()):
y[name] = [[] for xi in x]
for i in range(len(x)):
for name in sorted(phi.keys()):
# find (global) minimia of objective function
f = phi[name](u, 1.0) + phi[name](u - x[i], 1.0)
f_min = np.min(f)
y[name][i] = [j for j in find_local_minima(f) if f[j] <= f_min + 1.0e-3]
%%capture
# animate solution for sum of two penalty functions
from matplotlib import animation
from IPython.display import HTML
def animate(fnum, u, x, y):
plt.clf()
for a in ax:
a.clear()
for name in sorted(phi.keys()):
# plot the (normalized) objective
f = phi[name](u, 1.0) + phi[name](u - x[fnum], 1.0)
f = f / np.max(f)
p = ax[0].plot(u, f)
# indicate minima
for i in y[name][fnum]:
ax[0].plot([u[i], u[i]], [-0.1, f[i]], color=p[0].get_color(), ls='--')
# plot y (global minima) as a function of x
y_min = [min(u[i]) for i in y[name][:fnum+1]]
y_max = [max(u[i]) for i in y[name][:fnum+1]]
ax[1].plot(x[:fnum+1], y_min, color=p[0].get_color())
ax[1].plot(x[:fnum+1], y_max, color=p[0].get_color())
ax[0].set_xlim(u[0], u[-1]); ax[0].set_ylim(-0.1, 1.0)
ax[0].set_xlabel(r"$u$"); ax[0].set_ylabel(r"$f$ (scaled)")
ax[0].set_title("Objective")
ax[1].set_xlim(x[0], x[-1]); ax[1].set_ylim(-4.0, 4.0)
ax[1].set_xlabel(r"$x$"); ax[1].set_ylabel(r"$y$")
ax[1].set_title("Solution")
return ax
fig = plt.figure()
ax = [plt.subplot(1, 2, 1), plt.subplot(1, 2, 2)]
ani = animation.FuncAnimation(fig, animate, fargs=(u, x, y), interval=100, frames=len(x), blit=False, repeat=False)
plt.close(fig)
# display using video or javascript
HTML(ani.to_html5_video())
#HTML(ani.to_jshtml())
<Figure size 432x288 with 0 Axes>
We now plot the solution sets separately for each penalty function.
cm = plt.rcParams["axes.prop_cycle"].by_key()["color"]
for indx, name in enumerate(sorted(phi.keys())):
y_min = [min(u[i]) for i in y[name]]
y_max = [max(u[i]) for i in y[name]]
fig = plt.figure()
plt.plot(x, y_min, color=cm[indx % len(cm)])
plt.plot(x, y_max, color=cm[indx % len(cm)])
for i in range(len(x)):
if len(y[name][i]) > 2:
for uj in [u[jy] for jy in y[name][i] if ((u[jy] != y_min[i]) and (u[jy] != y_max[i]))]:
plt.plot(x[i], uj, color=cm[indx % len(cm)], marker='.')
if name == "Huber":
plt.fill_between(x, y_max, y_min, color=cm[indx % len(cm)])
plt.xlabel(r"$x$"); plt.ylabel(r"$y$"); plt.title(name)
plt.show()