Question:
What is the state of the art on analytic approximations of $\min$ and $\max$? My hunch is that numerical analysts probably have a better solution than the one I propose here.
For any $\epsilon$, I'd like analytic approximations $f$ that
Are accurate to within $\epsilon$, in the sense that $$|f(V)-\max(V)|<\epsilon$$ for any vector $V$ in $[-1000,1000]^{100}$, and
Minimize the expected relative error:
$$E\left[\left|\frac{f(V)-\max(V)}{\max(V)}\right|\right]$$ where $V$ is a random vector uniformly distributed in $[-1000,1000]^{100}$.
Motivation:
Within the context of optimisation, differentiable approximations of the min and max operators are very useful. In particular I am looking for $f_N,g_N \in C^\infty$:
\begin{equation} f_N: \mathbb{R}^n \rightarrow \mathbb{R} \end{equation}
\begin{equation} g_N: \mathbb{R}^n \rightarrow \mathbb{R} \end{equation}
where $\forall x \in \mathbb{R}^n \forall \epsilon >0 \exists N \in \mathbb{N}$:
\begin{equation} \forall m > N, \max(\lvert f_m(x)-\max(x) \rvert,\lvert g_m(x)-\min(x) \rvert) < \epsilon \end{equation}
I found a few proposed solutions to a related question on MathOverflow but when I tested these methods I found that none of them were numerically stable with respect to the relative error $\frac{\lvert \delta x \rvert}{\lvert x \rvert}$. Without much loss of generality, I shall focus on approximations of the max operator.
In total, I analysed three different analytic approximations to the max operator $\forall X \in \mathbb{R}^n, \forall N \in \mathbb{N}$:
- The generalised mean:
\begin{equation} GM(X,N)= \big(\frac{1}{n}\sum_{i=1}^n x_i^N\big)^{\frac{1}{N}} \tag{1} \end{equation}
- Exponential generalised mean:
\begin{equation} EM(X,N)= \frac{1}{N} \cdot \log \big(\frac{1}{n}\sum_{i=1}^n e^{N \cdot x_i}\big) \tag{2} \end{equation}
- The smooth max:
\begin{equation} SM(X,N)= \frac{\sum_{i=1}^n x_i \cdot e^{N \cdot x_i}}{\sum_{i=1}^n e^{N \cdot x_i}} \tag{3} \end{equation}
and found that all of these methods were vulnerable to overflow errors. In fact, I created an IJulia notebook where I analysed each method.
My proposed method:
This motivated me to come up with my own solution inspired by the properties of the infinity norm where I first rescale the vectors so they have zero mean and unit variance:
\begin{equation} AM(\hat{X},N) = \sigma_X \cdot \big(\frac{1}{N} \log \big(\sum_{i=1}^n e^{N\cdot \hat{x_i}} \big)\big) + \mu_{X} \tag{4a} \end{equation}
\begin{equation} \hat{X} = \frac{X - 1_n\cdot \mu_X}{\sigma_X} \tag{4b} \end{equation}
whose partial derivative with respect to $\hat{x_i}$ is simply the softmax:
\begin{equation} \frac{\partial}{\partial \hat{x_i}} AM(\hat{X},N) = \frac{e^{N \cdot \hat{x_i}}}{\sum_{i=1}^n e^{N\cdot \hat{x_i}}} \tag{5} \end{equation}
This may be be used to approximate both the min and max operators on $\mathbb{R}^n$ in the Julia programming language:
using Statistics
function analytic_min_max(X::Array{Float64, 1},N::Int64,case::Int64)
"""
An analytic approximation to the min and max operators
Inputs:
X: a vector from R^n where n is unknown
N: an integer such that the approximation of max(X)
improves with increasing N.
case: If case == 1 apply analytic_min(), otherwise
apply analytic_max() if case == 2
Output:
An approximation to min(X) if case == 1, and max(X) if
case == 2
"""
if (case != 1)*(case != 2)
return print("Error: case isn't well defined")
else
## q is the degree of the approximation:
q = N*(-1)^case
mu, sigma = mean(X), std(X)
## standardise the vector so it has zero mean and unit variance:
Z_score = (X.-mu)./sigma
exp_sum = sum(exp.(-Z_score*q))
log_ = log(exp_sum)/q
return (log_*sigma)+mu
end
end
and as expected it passed the numerical stability test that I defined:
function numerical_stability(method,type::Int64)
"""
A simple test for numerical stability with respect to the relative error.
Input:
method: the approximation used
type: 1 for min() and 2 for max()
Output:
Check that the average relative error is less than 10%.
"""
## test will be run 100 times
relative_errors = zeros(100)
for i = 1:100
## a vector sampled uniformly from [-1000,1000]^100
X = (2*rand(100).-1.0)*1000
## the test for min operators
if type == 1
min_ = minimum(X)
relative_errors[i] = abs(min_-method(X,i))/abs(min_)
## the test for max operators
else
max_ = maximum(X)
relative_errors[i] = abs(max_-method(X,i))/abs(max_)
end
end
return mean(relative_errors) < 0.1
end
References:
Wikipedia contributors. Smooth maximum. Wikipedia, The Free Encyclopedia. March 25, 2019, 21:07 UTC. Available at: https://en.wikipedia.org/w/index.php?title=Smooth_maximum&oldid=889462421. Accessed February 12, 2020.
alwayscurious (https://stats.stackexchange.com/users/194748/alwayscurious), What is the reasoning behind standardization (dividing by standard deviation)?, URL (version: 2019-03-18): https://stats.stackexchange.com/q/398116
Sergey Ioffe & Christian Szegedy. Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift. 2015.
J. Cook. Basic properties of the soft maximum. Working Paper Series 70, UT MD Anderson Cancer Center Department of Biostatistics, 2011. https://www.johndcook.com/soft_maximum.pdf
M. Lange, D. Zühlke, O. Holz, and T. Villmann, "Applications of lp-norms and their smooth approximations for gradient based learning vector quantization," in Proc. ESANN, Apr. 2014, pp. 271-276.
Aidan Rocke. analytic_min-max_operators(2020).GitHub repository, https://github.com/AidanRocke/analytic_min-max_operators
