Lipschitz Optimisation for Lipschitz Interpolation

TL;DR

This paper introduces a novel method for Lipschitz-based nonparametric learning that optimizes metric parameters to improve prediction accuracy, addressing robustness and local minima issues in existing techniques.

Contribution

It proposes a parameter optimization approach for Lipschitz metrics using validation errors and leverages Lipschitz properties to ensure global optimization in nonparametric learning.

Findings

Demonstrates improved prediction accuracy on artificial data

Shows robustness to noise in real-world datasets

Outperforms existing Lipschitz-based methods

Abstract

Techniques known as Nonlinear Set Membership prediction, Kinky Inference or Lipschitz Interpolation are fast and numerically robust approaches to nonparametric machine learning that have been proposed to be utilised in the context of system identification and learning-based control. They utilise presupposed Lipschitz properties in order to compute inferences over unobserved function values. Unfortunately, most of these approaches rely on exact knowledge about the input space metric as well as about the Lipschitz constant. Furthermore, existing techniques to estimate the Lipschitz constants from the data are not robust to noise or seem to be ad-hoc and typically are decoupled from the ultimate learning and prediction task. To overcome these limitations, we propose an approach for optimising parameters of the presupposed metrics by minimising validation set prediction errors. To avoid poor performance due to local minima, we propose to utilise Lipschitz properties of the optimisation objective to ensure global optimisation success. The resulting approach is a new flexible method for nonparametric black-box learning. We provide experimental evidence of the competitiveness of our approach on artificial as well as on real data.

Tables

Table 1. TABLE I: Sample statistics.

Data set	Sample mean	Sample std	Size	Input dim.
UCI-CCPP	454.37	17.06	9568	4
PumaDyn-8hn	1.16	5.62	4915	8

Table 2. TABLE II: Measured abs. test error statistics.

Data set:	UCI-CCPP			PumaDyn-8hn
Statistic:	Mean	Std		Mean	Std	Median
Lin. Regr.	$3.64$	$2.81$		3.68	2.49	3.39
GP1	$206.18$	$144.98$		2.63	2.15	2.14
LACKI	$3.40$	$3.44$		3.28	2.55	2.69
POKI-LC	$3.25$	$2.98$		3.25	2.43	2.76
POKI-ARD	2.60	$2.79$		2.81	2.07	2.44

Equations

$$\widehat{{\mathcal{E}}_1}({\hat{\mathfrak{f}}}_{\mathfrak{n}};f):=\frac{1}{\left|{\mathcal{X}}_{sample}\right|}\sum _{x\in {\mathcal{X}}_{sample}}\left|{\hat{\mathfrak{f}}}_{\mathfrak{n}}(x)-f(x)\right|$$

$$\widetilde{{\mathcal{E}}_1}({\hat{\mathfrak{f}}}_{\mathfrak{n}};\tilde{f}):=\frac{1}{\left|{\mathcal{X}}_{sample}\right|}\sum _{x\in {\mathcal{X}}^{\text{test}}}\left|{\hat{\mathfrak{f}}}_{\mathfrak{n}}(x)-\tilde{f}(x)\right|.$$

$$\displaystyle\,\mathfrak{\hat{f}_{n}}\bigl{(}x;\theta(n),\mathcal{D}_{n}\bigr{)}:=\frac{1}{2}\mathfrak{u}_{n}(x;\theta(n)\bigr{)}+\frac{1}{2}\mathfrak{l}_{n}(x;\theta(n)\bigr{)}.$$

$$\hat{\ell }(n):=\underset{i\ne j}{\max }\frac{\left|{\tilde{f}}_i-{\tilde{f}}_j\right|}{\Vert s_i-s_j\Vert }.$$

$$\hat{}\theta :={\text{argmin}}_{\theta \in \Theta }\mathfrak{l}(\theta ;{\mathcal{D}}^{\text{cond}},{\mathcal{D}}^{\text{eval}})$$

$$\mathfrak{l}(\theta ;{\mathcal{D}}^{\text{cond}},{\mathcal{D}}^{\text{eval}})=\frac{1}{\left|{\mathcal{X}}^{\text{eval}}\right|}\sum _{x\in {\mathcal{X}}^{\text{eval}}}\left|\tilde{f}(x)-{\hat{\mathfrak{f}}}_{\mathfrak{n}}(x;\theta ,{\mathcal{D}}^{\text{cond}})\right|.$$

$$\Omega :\{\begin{array}{l}\Theta \to \mathbb{R}\\ \theta \mapsto \mathfrak{l}(\theta ;{\mathcal{D}}_n^{\text{cond}},{\mathcal{D}}_n^{\text{eval}})\end{array}$$

$$L(\Omega )\le \max \{L({\phi }_i)|i=1,\dots ,Q_n\}$$

$${\hat{\mathfrak{f}}}_{\mathfrak{n}}(q;\theta ;{\mathcal{D}}_n^{\text{cond}})$$

$$L({\phi }_i)$$

$$L(\Omega )$$

$$D_j:=\max \{a_1,\dots ,a_d,b_1,\dots ,b_d\}$$

$$\underline{M}:=\underset{j=1...,N}{\min }{f}_j-L{D}_j^p\le \underset{x\in I}{\min }f(x)$$

$$\overline{M}:=\underset{j=1,...,N}{\max }{f}_j+L{D}_j^p\ge \underset{x\in I}{\max }f(x).$$

Full text

I Introduction

Supervised machine learning methods are algorithms for inductive inference. On the basis of a sample, they construct (learn) a computable model of a data generating process that facilitates inference over the underlying ground truth function and aims to predict its function values at unobserved inputs.

Among supervised learning methods, nonparametric algorithms tend to offer greater flexibility to learn rich function classes. Unfortunately, many classical techniques for nonparametric regression, such as the Nadaraya-Watson estimator [21, 14] or the LOESS method, [6] suffer from a practical limitation: their regression performance depends on the choice of hyperparameters. While in principle, it would be possible to tune these to the data (in manner similar in spirit to the one we propose in this work), to the best of our knowledge, currently there is little understanding on how to do so with a global optimiser that offers theoretical performance guarantees on the optimisation solution. This means that in practice, one is left to engineer these hyperparameters (or the settings of an optimiser) by manual tuning in order to ensure good performance on a particular learning problem. Of course, this stands in opposition to the motivation for utilising nonparametric learning, especially in system identification: which is to facilitate flexible and fully automated black-box learning that does not require manual intervention.

Perhaps one of the most popular nonparametric machine learning method is Bayesian inference with Gaussian processes (GPs) [16]. GPs offer a flexible and principled probabilistic method for nonparametric regression and have evolved into one of the chief work-horses for learning dynamic systems [9, 10, 7, 15, 8, 11] in the research communities related to artificial intelligence and, more recently, also in control. However, they suffer from several limitations, including scalability to large data sets, a lack of understanding of how to bound the closed-loop dynamics resulting from controlling on the basis of a GP state-space model and, quite similar to the other classical regression methods, the question of how to choose a good prior and its hyper-parameters in a principled yet computable manner. To alleviate the last problem, it is common practice to tune hyper-parameters of a chosen (typically universal) kernel to explain the data via the marginal log-likelihood [16]. While often successful on many data sets, the result can be highly sensitive to the choice of optimiser, initialisations, data sets and computational budget. Unfortunately, little theoretical understanding of the important interplay between these components in the resulting inference mechanism seems to exist.

In contrast to such Bayesian methods this work builds on nonparametric regression techniques that harnesses Lipschitz regularity of the target function to provide bounds on the predictions of the target function at unobserved inputs. Applied to machine learning the basic idea is that Lipschitz continuity constrains the set of possible function values of a target function at a query input, dependent on the distance between the query and the previously observed training examples. A prediction is then made by choosing a function value in the middle of the set of possible function values. This idea, at least going back to the era of “Russian mathematics” [19], has been redeveloped and advanced under different headlines including Lipschitz Interpolation [23, 1], Nonlinear Set Membership (NSM) interpolation methods [12] and Kinky Inference (KI) [4]. The presupposed Lipschitz constant as well as the assumed input space metric are crucial hyper-parameter choices of these methods that can drastically affect the predictive behaviour of trained inference rule. While a variety of Lipschitz constant estimators are known, they are designed independently from the prediction task and tend to be sensitive to noise. As an alternative, [12] proposed a method where the Lipschitz constant is estimated from a parametric model that is fitted to the data (e.g. a linear model or a neural network). However, it remains unclear in how far this Lipschitz constant will aide the predictive performance of the Lipschitz interpolation model.

In this work, we propose a different approach. Firstly, we consider the more general problem of optimising for parameters of a chosen pseudo-metric. (As we will see the Lipschitz constant determination problem can be cast as a special case of this). We then determine these free parameters by minimising an empirical estimate of the prediction error directly of the KI rule. When purposefully stating the error as an $\ell_{1}$-loss, we can derive a Lipschitz constant for the optimisation objective. This allows us to employ Lipschitz optimisation. In contrast, to other hyperparameter tuning approaches utilised in nonparametric machine learning, this global optimisation approach offers bounds on the optimisation success and can avoid falling into suboptimal solutions.

The result of this merger of Lipschitz optimisation for parameter optimisation and KI (based on the newly found parameters) is a reliable and fast nonparametric machine learning approach which we will refer to as Parameter Optimised Kinky Inference (POKI).

Apart from the application to the automated determination of the Lipschitz constant in nonlinear set membership methods, we discuss other settings where our POKI approach uncovers periodicities or inherent low-dimensionality (automated relevance determination) to learn good predictive models from data.

II Kinky Inference and Lipschitz Interpolation

In this section, we will rehearse the class of learning rules sometimes referred to as Kinky Inference. They encompass a host of other methods such as Lipschitz Interpolation and Nonlinear Set Interpolation.

Setting. Let $f:\mathcal{X}\to\mathcal{Y}$ be a target or ground-truth function we desire to learn and let $\mathcal{X}$, $\mathcal{Y}$ be two spaces endowed with (pseudo-) metrics $\,\mathfrak{d}:\mathcal{X}^{2}\to\mathbb{R}_{\geq 0},\,\mathfrak{d}_{\mathcal{Y}}:\mathcal{Y}^{2}\to\mathbb{R}_{\geq 0}$, respectively.

For simplicity, we restrict our exposition to real-valued targets with $\,\mathfrak{d}_{\mathcal{Y}}(y,y^{\prime})=\left|y-y^{\prime}\right|$.

Assume that, at time step $n$, we have access to a sample or data set \mathcal{D}_{n}:=\{\bigl{(}s_{i},\tilde{f}_{i}\bigr{)}\,|\,i=1,\ldots,N_{n}\} containing $N_{n}\in\mathbb{N}$ (possibly corrupted) sample values $\tilde{f}_{i}\in\mathcal{Y}$ of target function $f$ at sample input $s_{i}\in\mathcal{X}$. The sampled function values are allowed to have observational error given by a (potentially stochastic) error function $\mathfrak{e}:\mathcal{X}\to\mathbb{R}^{m}_{\geq 0}$. That is, all we know is that $\tilde{f}_{i}=f(s_{i})+\mathfrak{e}(s_{i})$. For convenience, we may also write $\mathcal{D}_{n}=(\mathcal{X}_{n},\mathcal{Y}_{n})$ where $\mathcal{X}_{n}=\{s_{i}|i=1,...,N_{n}\}\subset\mathcal{X}$ is the collection of sample inputs and $\mathcal{Y}_{n}=\{\tilde{f}_{i}|i=1,...,N_{n}\}\subset\mathcal{Y}$ is the sequence of observed function values. It is our aim to learn target function $f$ in the sense that utilise the available data $\mathcal{D}_{n}$ to infer predictions $\,\mathfrak{\hat{f}_{n}}(x)$ of $f(x)$ at unobserved query inputs $x\notin\mathcal{X}_{n}$. In our context, the evaluation of $\,\mathfrak{\hat{f}}_{n}$ is what we refer to as (inductive) inference or prediction.

The entire function $\,\mathfrak{\hat{f}_{n}}$ that is learned to facilitate predictions is referred to as the predictor.

To ground the inference, a priori assumptions are necessary. In our context, we will generally assume that the target can be arbitrarily well approximated by Lipschitz continuous function with some Lipschitz constant. Remember, relative to our chosen pseudo-metrics, a real-valued function $\phi$ is Lipschitz continuous (with Lipschitz constant $\ell\geq 0$) on domain $I\subset\mathcal{X}$ if $\,\mathfrak{d}_{\mathcal{Y}}(\phi(x),\phi(x^{\prime}))\leq\ell\,\mathfrak{d}(x,x^{\prime}),\forall x,x^{\prime}\in I$. Note the metrics, may depend on a parameter $\theta$. In this case, we highlight this dependency explicitly by writing $\,\mathfrak{d}(\cdot,\cdot;\theta)$ instead of $\,\mathfrak{d}(\cdot,\cdot)$. Of course, the Lipschitz parameter can be absorbed into the parameter of a chosen pseudo-metric.

The approximation success of $\,\mathfrak{\hat{f}_{n}}$ relative to a target can be measure by various metrics. In this work, we will be most interested in the $\mathcal{L}_{1}$-prediction error $\mathcal{E}_{1}(\,\mathfrak{\hat{f}_{n}};f):=\left\|\,\mathfrak{\hat{f}_{n}}-f\right\|_{1}=\int_{\mathcal{X}}\left|\,\mathfrak{\hat{f}_{n}}(x)-f(x)\right|\text{ d}x$. In practice, this error can be estimated by the * empirical prediction error* estimate

[TABLE]

where the sample $\mathcal{X}_{sample}$ is a finite set of sample inputs.

The average test set prediction error serves as a surrogate measure for the true prediction error; if the test set is chosen sufficiently dense (and the predictor and target are continuous) then $\hat{\mathcal{E}_{1}}\approx\mathcal{E}_{1}$. In the case of i.i.d. samples, we can also construe Eq. 1 as a Monte-Carlo estimate of the $\mathcal{L}_{1}$- error.

In case we do not have access to a the target (or a noise-free sample), we might have to base our assessment of the method on the *empirical sample prediction error *

[TABLE]

Learning rule. In this work we will expand on the basis of a special case of the class of kinky inference predictors [4] to perform learning as inference over unobserved function values. The prediction rule can be stated as follows:

Definition II.1** (Kinky inference (KI) rule (simplified) ).**

Given access to a sample set $\mathcal{D}_{n}$ and an input space pseudo-metric $\,\mathfrak{d}(\cdot,\cdot;\theta(n)):\mathcal{X}^{2}\to\mathbb{R}$ parameterised by $\theta(n)$, we define the KI predictor by \,\mathfrak{\hat{f}_{n}}\bigl{(}\cdot;\theta(n),\mathcal{D}_{n}\bigr{)}:\mathcal{X}\to\mathcal{Y} to perform inference over function values as per:

[TABLE]

Here, \mathfrak{u}_{n}\bigl{(}\cdot;\theta(n)\bigr{)},\mathfrak{l}_{n}\bigl{(}\cdot;\theta(n)\bigr{)}:\mathcal{X}\to\mathbb{R}^{m} are called ceiling and floor functions, respectively. They are given by \mathfrak{u}_{n}\bigl{(}x;\theta(n)\bigr{)}:=\min_{i=1,\ldots,N_{n}}\tilde{f}_{i}+\,\mathfrak{d}(x,s_{i};\theta(n)) and \mathfrak{l}_{n}\bigl{(}x;\theta(n)\bigr{)}:=\max_{i=1,\ldots,N_{n}}\tilde{f}_{i}-\,\mathfrak{d}(x,s_{i};\theta(n)), respectively.

In the literature, various generalisations and special cases exist. For instance, Calliess [4] proposes a generalised framework called Kinky Inference, that also allows for the specification of additional parameters. These include functions that allow the incorporation of a priori knowledge about upper and lower bounds on the target as well as upper and lower bounds on observational noise.

A special case arises for the choice of $\,\mathfrak{d}(x,y;\theta(n))=\ell(n)\left\|x-y\right\|$ which is referred to as Lipschitz Interpolation [1] or as Nonlinear Set Interpolation [12]. Here the parameter $\theta(n)=\ell(n)$ is the supposed Lipschitz constant of the target. Typically this constant is assumed to be either known a priori or estimated lazily from the data, e.g.  [18] as follows:

[TABLE]

Unfortunately, the latter estimate has the problem of being unbounded in the presence of observational noise. In the case of bounded noise, [4] proposed to utilise the alternative estimator $\ell(n):=\max_{i\neq j}\frac{\left|\tilde{f}_{i}-\tilde{f}_{j}\right|-2\bar{\mathfrak{e}}}{\left\|s_{i}-s_{j}\right\|}$ where $\bar{\mathfrak{e}}:=\sup_{x}\left|\mathfrak{e}(x)\right|$ is an upper bound on the (zero-mean) noise. While this prevents the estimates to blow up in the bounded noise case, generally it is not clear how to obtain $\bar{\mathfrak{e}}$. And, if $\bar{\mathfrak{e}}$ is chosen too conservatively large then the prediction quality is questionable. As an alternative approach, we will employ Lipschitz optimisation to find a parameter $\theta(n)$ that minimises empirical prediction error on a validation data set.

III Parameter Optimised Kinky Inference (POKI)

Given some separate data sets $\mathcal{D}^{\text{cond}}$ and $\mathcal{D}^{\text{eval}}$, we will aim to choose a parameter of the pseudo-metric that minimises the empiricial sample prediction error. That is, we choose our parameter to be

[TABLE]

where $\Theta$ denotes a predefined parameter space and the loss quantifies the empirical sample prediction error quantified by the loss function $\,{\mathfrak{l}}(\cdot;\mathcal{D}^{\text{cond}},\mathcal{D}^{\text{eval}}):\mathcal{X}\to\mathbb{R}$ with

[TABLE]

Here, we refer to \mathcal{D}^{\text{eval}}=\bigl{(}\,\mathcal{X}^{\text{eval}},\,\mathcal{Y}^{\text{eval}}\bigr{)} as the available parameter evaluation data and $\mathcal{D}^{\text{cond}}$ as the conditioning data.

To avoid overfitting and exact interpolation of the noise, we will generally ensure that the conditioning and evaluation data are different, i.e. at least $\mathcal{D}^{\text{eval}}\neq\mathcal{D}^{\text{cond}}$. Ideally, they should be drawn independently from each other. We will give an illustration of the utility of this approach below. For now, we will confine ourselves to describe the approach we propose:

Def. (Parameter Optimised Kinky Inference (POKI)): Assume at trial $n$, we are given access to a set of training examples $\mathcal{D}_{n}$ of size $N_{n}$.

We will undergo the following steps:

We randomly partition the data into two sub sets: a conditioning data set $\mathcal{D}^{\text{cond}}_{n}$ and a parameter evaluation data set $\mathcal{D}^{\text{eval}}_{n}$, $\mathcal{D}^{\text{tex}}_{n}=\mathcal{D}^{\text{cond}}_{n}\,{\dot{\cup}}\,\mathcal{D}^{\text{eval}}_{n}$. Unless explicitly stated otherwise, both sets will be made close to equal in size. 2. 2.

Utilising the conditioning and evaluation data sets, we compute the minimum-loss parameter estimate $\hat{}\theta_{n}:=\text{argmin}_{\theta\in\Theta}\,{\mathfrak{l}}(\theta;\mathcal{D}^{\text{cond}}_{n},\mathcal{D}^{\text{eval}}_{n})$ (cf. Eq. 4). 3. 3.

For prediction of future query inputs $q\in\mathcal{X}$, we use the KI prediction $\,\mathfrak{\hat{f}_{n}}(q;\hat{}\theta_{n},\mathcal{D}_{n})$ (cf. Def. II.1), utilising the learned newly identified parameter estimate $\hat{}\theta_{n}$ and conditioning on the full set of available data.

We refer to the resulting predictor $\,\mathfrak{\hat{f}_{n}}(\cdot;\hat{}\theta_{n},\mathcal{D}^{\text{tex}}_{n})$ as a Parameter Optimised Kinky Inference (POKI) rule.

III-A Lipschitz optimisation and constant determination

In order to guarantee the quality of the parameter estimate $\hat{}\theta_{n}$ even in the presence of non- differentiabilities and local optima of the loss function, we propose to utilise Lipschitz optimisation methods. In one dimensional settings we could for instance utilise Shubert’s method [17]. Various extensions to the multi-dimensional case exist including stochastic methods that offer probabilistic bounds that allow for guarantees with computational effort scaling linearly with the number of dimensions [24] as well as recent local approaches based on regret analysis [13]. In this work we use the simple approach sketched in Sec. 2.5 of [4]. While the effort to reduce the given worst-case bound scales exponentially in the number of dimensions, the approach has the advantage of being simple and giving deterministic bounds. However, we intend to try out the more advanced, scalable approaches to Lipschitz optimisation in the course of future work.

In our chosen Lipschitz optimisation method, if we know a Lipschitz constant $L(\Omega)$ of the objective function

[TABLE]

then we can find a guaranteed global minimum up to a worst-case error bound that can be pre-specified in advance, thereby avoiding the pitfall of local minima.

To this end, we need to determine the Lipschitz constant of $\Omega$, i.e. of the loss as a function of the parameter $\theta$. We denote the Lipschitz constant of a function $f$ by $L(f)$. For simplicity, we choose the metric induced by the maximum-norm to define Lipschitz continuity of $\Omega$. Hence, for a given loss we desire to derive a nonnegative number $L(\Omega)\in\mathbb{R}_{+}$ such that: $\forall\theta,\theta^{\prime}\in\Theta:\left|\Omega(\theta)-\Omega(\theta^{\prime})\right|\leq L(\Omega)\,\left\|\theta-\theta^{\prime}\right\|_{\infty}.$

Let the conditioning data be given by $\mathcal{D}^{\text{cond}}_{n}=\{(\mathfrak{x}_{i},\mathfrak{y}_{i})|i=1,\ldots,M_{n}\}$ and the evaluation data be given by $\mathcal{D}^{\text{eval}}_{n}=\{(x_{i},y_{i})|i=1,\ldots,Q_{n}\}$ such that the full data is $\mathcal{D}_{n}=\{(s_{i},\tilde{f}_{i})|i=1,...,N_{n}\}=\mathcal{D}^{\text{cond}}_{n}\cup\mathcal{D}^{\text{eval}}_{n}$.

Appealing to the definition of the loss and the rules stated in Lem. .1 (refer to the appendix), we see that

[TABLE]

where $\varphi_{i}:\theta\mapsto\,\mathfrak{\hat{f}_{n}}(x_{i};\theta,\mathcal{D}^{\text{cond}}_{n})$ is the prediction of the $i$th evaluation example input as a function of the parameter $\theta$. We will now bound the Lipschitz constant of each $\varphi_{i}$:

By definition of the predictor we have

[TABLE]

for some indices $k,j\in\{1,\ldots,M_{n}\}$ and all query inputs $q\in\mathcal{X}$.

Appealing to Lem. .1.8), we thus see that that for all $i$ we have

[TABLE]

where each $\,\mathfrak{d}(x_{i},\mathfrak{x}_{k};\cdot)$ is a function of the parameter of the pseudo-metric $\,\mathfrak{d}$. Since the Lipschitz constant of the loss is bounded from above by the maximum of the Lipschitz constants of all the $\varphi_{i}$ (cf. Eq. 6) and since $\mathcal{D}^{\text{cond}}_{n},\mathcal{D}^{\text{eval}}_{n}\subset\mathcal{D}_{n}$, we obtain the bound:

[TABLE]

where, for every pair $s_{j},s_{k}$ of inputs in the available data, \,\mathfrak{d}(s_{j},s_{k};\cdot)\bigr{)} is a function of hyper-parameter $\theta$. So, the determination of a bound on the desired Lipschitz constant depends on the conditioning data and the arithmetic definition the parameter $\theta$ within the chosen pseudo-metric $\,\mathfrak{d}$. For several cases of general interest, we will next provide some Lipschitz constant derivations. The best Lipschitz constant of a function $\phi$ will be denoted by $L(\phi)$.

• •

Pseudo-metric for periodic functions: If we know the target to be periodic with some frequency we may choose the pseudo-metric $\,\mathfrak{d}(x,x^{\prime};\theta)=\left|\sin(\pi\theta\left|x-x^{\prime}\right|)\right|$. Knowing that taking the absolute value does not change a Lipschitz constant and that the Lipschitz constant of a differentiable function is the supremum of the absolute value of its derivative, we see that L\bigl{(}\,\mathfrak{d}(x,x^{\prime};\cdot)\bigr{)}\leq\pi\left|x-x^{\prime}\right|. With reference to Eq. 9 it follows that in our case, a Lipschitz bound on the loss is given by $\pi$ times the diameter of the sample input set: $L(\Omega)\leq\pi\max\{\left|s_{j}-s_{k}\right|\,|\,j,k=1,\dots,N_{n}\}.$

Note, the parameter $\theta$ determines the frequency with which the pseudo-metric cannot distinguish two inputs in a periodic input space. Choosing a parameter $\theta$ that matches the frequency of the target would be important to achieve adequate predictive performance. Conversely, if one presupposes a false frequency $\theta\neq 2$ then the prediction performance can suffer dramatically. Therefore, when optimising the parameter finding avoiding local minima is critical. Here, our global Lipschitz optimisation approach shines. For an illustration, refer to Fig. 1.

• •

Automated Lipschitz constant determination: As mentioned above, we might seek to perform Lipschitz interpolation with automated Lipschitz constant determination. To this end, we might choose the metric $\,\mathfrak{d}(x,y;\theta_{n})=\theta_{n}\left\|x-y\right\|$ where parameter $\theta$ acts as the assumed Lipschitz constant of the predictor. Clearly, we have $L(\,\mathfrak{d}(s_{i},s_{k};\cdot))\leq\left|s_{j}-s_{k}\right|$. Thus, the Lipschitz constant of the optimisation objective $\Omega$ can be bounded by the diameter of the sample input set $\mathcal{X}_{n}$: $L(\Omega)\leq\max\{\left|s_{j}-s_{k}\right|\,|\,j,k=1,\dots,N_{n}\}$.

We refer to the resulting POKI inference rule as POKI-LC (where LC stands for Lipschitz constant). Illustrations of the performance of this approach in contrast other methods are provided in the experimental section.

• •

Automated Relevance Determination (ARD). A more general case of Lipschitz parameter determination is Automated Relevance Determination (ARD). Here, the goal is to find weights encoding the extent to which input space parameters are relevant for the prediction. This is interesting in high-dimensional input space where the input vectors might be features whose predictive role might be unclear a priori. To facilitate ARD, we might choose the metric $\,\mathfrak{d}(x,x^{\prime};\theta)=\max_{i=1,\ldots,d}\theta_{i}\left|x_{i}-x_{i}^{\prime}\right|$. The weight $\theta_{i}$ quantifies the degree to which the $i$th input component should contribute to the prediction. Large weights suggest a strong degree of importance: i.e. small deviations of a query $x_{i}$ from the closest example $s_{j,i}$ in the $i$th component will result in large uncertainty of the prediction. Conversely, if $\theta_{i}=0$ effectively disables any influence of the $i$th input dimension in the prediction process.

Analogously to the kernel literature, inferring $\theta$ from the data will be referred to as automated relevance determination (ARD) and the resulting POKI rule will be referred to as POKI-ARD.

To facillitate Lipschitz optimisation of the parameter vector, we need to derive a Lipschitz constant bound. Again appealing to the rules of Lipschitz arithmetic, we see that $L(\,\mathfrak{d}(s_{i},s_{k};\cdot))\leq\left\|s_{i}-s_{k}\right\|_{\infty}$ and thus, L(\Omega)\leq\max\bigl{\{}\left\|s_{j}-s_{k}\right\|_{\infty}\,|\,j,k=1,\dots,N_{n}\bigr{\}}. In other words, the Lipschitz constant of the loss is bounded from above by the diameter of the data relative to the standard maximum-norm.

As a first illustration, consider learning the dynamics of a friction-less pendulum based on a sample. Denoting the angle of a pendulum by $q$, its dynamics are given by the ordinary differential equation $\ddot{q}=f(x)$ where [m;b;l;g] $x=[q;\dot{q}]$ is called the state and $f(x)=-9.81\sin(x_{1})$. Here, we learn the dynamics by learning the function $f$ based on a sample of triplets of angles, velocities and accelerations. We can see that this target function only depends on the first component of its input, namely the angle. So, if we incorporate this in our KI method by choosing the metric $\,\mathfrak{d}(x,x^{\prime};\theta)=\max\{\theta_{1}\left|x_{1}-x_{1}^{\prime}\right|,\theta_{2}\left|x_{2}-x_{2}^{\prime}\right|\}$ with $\theta=[1;0]$ we should expect superior predictive power over using the standard maximum norm (where $\theta=[1;1]$). Therefore, we expect ARD method to automatically uncover this advantageous parameter. As an example affirming this intuition, refer to Fig. 2. The example of Fig. 2(a) - Fig. 2(c) illustrate how the ARD method can uncover the fact that the target function only depends on one of the input components and how this can affect the predictive performance.

IV Experiments

In this section, we compare our approach to a number of well-established machine learning methods both on artificial and on real data.

IV-A Artifical data.

In order to have access to the ground-truth, we first tested the method on a sequence of artificial benchmark regression tasks. As the ground truth, we chose the target function $f:\mathcal{X}:=[0,1]^{d}\to\mathbb{R},x\mapsto\left|-\cos(2\pi\,x_{1})\right|+x_{1}$. Its values have a linear trend as well as a strongly nonlinear one, varying only with the first input component.

The data was obtained by artificially sampling uniformly $\mathcal{X}$ and superimposing the function values with i.i.d. zero-mean Gaussian noise with standard deviation $\frac{1}{4}$.

That is the sample was obtained from randomly selected inputs of the noisy data-generating “test function” $\tilde{f}(x)=f(x)+\nu(x)$ where $\nu(x)\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\mathcal{N}(0,\frac{1}{16})$.

As mentioned in the introduction most general nonparametric regression techniques’ predictive performance hinges on the choice of a set of hyper-parameters. For instance, in the Loess method [6] there is a smoothness parameter to be set a priori. For an illustration of the strong impact of this parameter to the predictor refer to Fig. 3.

We can see the behaviour ranges between over-smoothing and overfitting to the noise (e.g. the yellow curve). Unfortunately, for the Loess method we are not aware of a systematic method for finding the parameter in question. Of course, in principle, our approach we use for the KI regression approach might be applicable as well. That said, currently we are not sure how to determine a Lipschitz constant (or even gradient) of the test error as a function of the parameter that would yield an analogous approach to the POKI method we propose in this paper. Therefore, for the time being, we will restrict our comparisons to other methods where automatic (hyper-) parameter determination techniques are available.

In particular, we trained predictors based on the following learning methods:

Linear regression model (Lin. Mod.) trained on the data with the least-squares method. 2. 2.

A Gaussian process (GP) provided with the correct observational noise level. The prior was based on a zero-mean function $\mu(\cdot)\equiv 0$ and a squared-exponential (SE) ARD-kernel k(x,x^{\prime})=\omega\exp\bigl{(}-\sum_{i=1}^{d}\frac{(x_{i}-x_{i})^{2}}{\ell_{i}}\bigr{)} with relevance length scale parameters $\ell_{i}$ chosen uniformly to be $\ell_{i}=0.5,\forall i$ and output-scale parameter chosen to be $\omega=1$. 3. 3.

A kinky inference predictor (LACKI) with $\,\mathfrak{d}(x,x^{\prime};\theta)=\theta\left\|x-x^{\prime}\right\|_{\infty}$ where the parameter was the Lipschitz constant $\theta=L(n)$ lazily estimated as per Eq. 3. This is equivalent to the LACKI method described in [3] with hyper-parameter choices $\underline{L}=0$ and $\lambda=0$. 4. 4.

Our parameter optimised kinky inference approach (POKI-LC) where again $\,\mathfrak{d}(x,x^{\prime};\theta)=\theta\left\|x-x^{\prime}\right\|_{\infty}$. 5. 5.

Our parameter optimised kinky inference approach (POKI-LC2) trained as POKI-LC but where the parameter was optimised with Brent’s method instead of Lipschitz optimisation. 6. 6.

Our parameter optimised kinky inference approach (POKI-ARD) where the metric was an ARD metric with relevance parameters optimised as described above. 7. 7.

A Gaussian process (GP-ARD) with a SE-ARD kernel as above, but whose relevance hyper-parameters were optimised with the conjugate-gradient approach to maximise the marginal log-likelihood [16].

Exp.1 : As a first illustration, we considered regression on one-dimensional input space ($d=1$) based on a sample of 84 noise-corrupted data points. The results for a selection of the methods are depicted in Fig. 4(a). We can see that the linear regression method accurately fitted the linear trend in the data but does not explain the nonlinear variation. For the given hyper-parameter settings, the GP was clearly over-smoothing the data. This might well have to do with the choice of assumed observational noise and length scale. To obtain better results we have also tried GP hyper-parameter optimisation on the data (GP-ARD). The conjugate gradient optimiser managed to determine sensible hyper-parameters, thereby reducing the over-smoothing markedly. In contrast, the LACKI approach uncovered the nonlinearity well but overfitted to the noise in the data. Finally, our POKI-LC approach was able to smooth out the noise sufficiently well to offer a good fit to both the linear as well as to the non-linear trend in the target function. Being of comparable prediction accuracy as the solution found by GP-ARD, POKI’s predictor was markedly less smooth than the GP’s predictor with the chosen kernel and thereby, being able to fit the non-smooth parts of the target function better.111Of course, we could have handpicked a less smooth kernel class for the GP to potentially match the non-smooth parts of target better. However, we have opted for the SE class as this is the universal kernel class most commonly utilised in “black-box” learning.

Exp.2 : We desired to gain an impression of the performance of our POKI methods with (i) increasing input space dimensionality and (ii) in the limit of increasing sample size. In a third set of experiments (iii), we also the repeated the experiments in (ii) but for training data corrupted by noise $\nu$ drawn i.i.d. from a uniform over the intervall $[-0.5,0.5]$. The results are depicted in Fig. 5, Fig. 6 and Fig. 7, respectively. As performance metrics of interest, we considered the quality of the prediction and the computational effort required for training as well as for predicting with the trained models. The predictive accuracy was estimated by the empirical absolute prediction error means (cf. Eq. 1) estimated on the basis of a test set sample of 4000 inputs drawn i.i.d. from a uniform over the input space. The test sample was drawn independently from the training examples. Runtime measurements were based on code written in Julia 4.7 running on a single core of a 2.5 GHz i7 processor on a laptop with 16GB RAM. Our evaluation of the GPs was based on the implementation provided by the GaussianProcesses.jl library.

Examining the plots we note the following observations: Firstly, for all POKI methods, the prediction error seemed to vanish with increasing sample size. Secondly, the ARD methods were able to yield better prediction accuracy. We attribute this to their ability to uncover the inherently one-dimensional structure of the functional relationship. In comparison to the tested alternatives, POKI-ARD seemed to be more reliable in finding good parameters (and hence, to yield superior prediction accuracy) than its competitors. We attribute this to the use of the global Lipschitz optimiser. This benefit comes at the cost of increased training duration in the current implementation of the optimiser in the ARD case, when the dimensionality of the data increases. Addressing this issue is deferred to future work. In Exp.2 (ii) and (iii), it appeared that the Lipschitz optimisation-based hyper-parameter tuning seemed to be particularly beneficial for a small- to mid- range training data density (e.g. compare the prediction error plots of POKI-LC vs POKI-LC2 as well as POKI-ARD vs GP-ARD). Furthermore, in Exp.2 (i), it emerged that both GP-based predictors performed comparably poorly as the dimensionality of the input space grew (cf. Fig. 5).

IV-B Real data.

We also tested the predictive power of our methods on two real data sets. With no access to the true target the performance estimates had to rely on the empirical sample prediction error (cf. Eq. 2). The first data set was the UCI-CCPP data set [20] where the task was to predict the output of a power plant based on sensor measurements and control inputs. The available data was partitioned randomly into a training and a test data set containing 10% and 90% of all available data, respectively. The second data set was the Pumadyn data set [2]. Simulating the dynamics of the robot arm of a 560 Puma robot, the data has been used to showcase the potency of Gaussian processes for learning dynamical systems [16]. The available data was partitioned randomly into a training and a test data set containing 60% and 40% of all available data, respectively.

The results are documented in Tab. II. While the GP with hyperparameter optimisation (GP1) outperformed all other methods on the Pumadyn data set, POKI-ARD exhibits similarly low test set error mean. Furthermore, on the CCPP data the GP model performed very poorly, while our POKI methods were estimated to have the lowest test set error.

IV-C Conclusions

We developed POKI, a method for nonparametric regression that optimises parameters of a chosen pseudo-metric with the aim to maximise predicitve performance. Our approach can be applied to automated determination of Lipschitz constants in nonlinear set membership methods thereby addressing an open problem in data-driven control in a more principled manner than previously proposed methods [12, 5]. In addition, we have also given examples where the automatically learned (hyper-) parameters were utilised for automated relevance determination and detected and leveraged periodicity in the data.

We have found that our POKI methods work. They can compensate for noise in the data without the necessity to know the correct noise distributions. They are flexible to learn nonlinear, non-smooth functions. While the performance will not always be better than competing methods (such as GPs), an advantage of our approach is that the inference only involves basic computational steps that we would expect to be numerically robust and executable even on reduced instruction set micro-controllers. Furthermore, the particular mathematical properties of the Lipschitz interpolation rule allow us to easily determine Lipschitz constants of the empirical prediction errors. This in turn allows us to determine hyper-parameters by global Lipschitz optimisation, avoiding the pitfalls of local optima. Therefore, our approach can be more flexible and reliable to perform well in a large variety of problem domains and data sets without any manual tweaking. This is in contrast to many other methods (such as GPs or deep learning methods), where the predictive performance can be sensitive to a priori choices, e.g. of priors, hyper-parameters, optimisers and initialisations.

With the present choice of Lipschitz optimiser, problems involving high-dimensional hyper-parameters can render the determination of the approximately optimal hyper-parameter computationally intractable. Future work intends to explore the utilisation or development of alternative Lipschitz optimisation approaches that promise improved scalability (for a fixed error bound) in the dimensionality of the optimisation problem (e.g. [24]).

We observed that our approach can smooth out noise and hence, can yield accurate prediction of the ground truth even in the presence of additive observational noise. A theoretical challenge that remains to be investigated is the asymptotics of the true $\mathcal{L}_{1}$-prediction error in the limit of increasing number of noisy data points. We would hope that these results would allow us to derive probabilistic guarantees for a data-driven controller that combines our learning method with stochastic MPC approaches. This would provide an important extension to existing work of NSM-based MPC that had to rely on the knowledge of the Lipschitz constant [5, 4].

Finally, we would like to point out that Lipschitz interpolation rules (i.e. KI or NSM predictors) have a mathematical structure that made it particularly easy to calculate a Lipschitz bound on the empirical $\ell_{1}$ prediction loss, facilitating global Lipschitz optimisation of the hyper-parameters. However, it might be worthwhile exploring in how far the Lipschitz optimisation approach can also be successfully employed in automated (hyper-) parameter tuning of other machine learning methods.

V Acknowledgements

I would like to thank Carl Rasmussen and Jan Maciejowski for helpful discussions and encouraging feedback. Also, funds via EPSRC NMZR/031 RG64733 are gratefully acknowledged.

-A *Hölder arithmetic and optimisation *

This work utilised global Lipschitz optimisation and, to facilitate this, we derived Lipschitz constants for various functions pointing to the algorithm and results stated in [4]. For convenience, we restate these results at this point followed by an outline of a simple method for bounded Hölder optimisation we utilised in this work.

-A1 Hölder arithmetic

As a slight generalisation of Lipschitz continuity, we say that a function $f:\mathcal{X}\to\mathcal{Y}$ is Hölder continuous on $I\subseteq\mathcal{X}$ with (Hölder) constant $L$ and (Hölder) exponent $p$ if $\,\mathfrak{d}_{\mathcal{Y}}(f(x),f(x^{\prime}))\leq L\,\mathfrak{d}_{\mathcal{X}}(x,x^{\prime})\forall x,x^{\prime}\in I$.

Lemma .1** (Hölder arithmetic).**

Let, $I,J\subset\mathcal{X}$ where $\mathcal{X}$ is a metric space endowed with metric $\,\mathfrak{d}$. Let $f:\mathcal{X}\to\mathbb{R}$ be Hölder on $I$ with constant $L_{I}(f)\in\mathbb{R}_{+}$ and $g:\mathcal{X}\to\mathbb{R}$ be Hölder on $J$ with constant $L_{J}(g)\in\mathbb{R}_{+}$. Assume both functions have the same Hölder exponent $p\in(0,1]$. That is, $\forall x,x^{\prime}\in\mathcal{X}:\left|f(x)-f(x^{\prime})\right|\leq L(f)\,\mathfrak{d}(x,x^{\prime})^{p}$ and $\forall x,x^{\prime}\in\mathcal{X}:\left|g(x)-g(x^{\prime})\right|\leq L(g)\,\mathfrak{d}(x,x^{\prime})^{p}$. We have:

Mapping $x\mapsto|f(x)|$ is Hölder on $I$ with constant $L_{I}(f)$ and exponent $p_{f}$. 2. 2.

If $g$ is Hölder on all of $J=f(I)$ the concatenation $g\circ f:t\mapsto g(f(t))$ is Hölder on $I$ with constant $L_{I}(g\circ f)\leq$ $L_{J}(g)\,L_{I}^{p}(f)$ and exponent $p^{2}$. 3. 3.

Let $r\in\mathbb{R}$. $r\,f:x\mapsto r\,f(x)$ is Hölder on $I$ having a constant $L_{I}(r\,f)=|r|\,L_{I}(f)$. 4. 4.

$f+g:x\mapsto f(x)+g(x)$* has Hölder constant at most $L_{I}(f)+L_{J}(g)$.* 5. 5.

Let $m_{f}=\sup_{x\in\mathcal{X}}f(x)$ and $m_{g}=\sup_{x\in\mathcal{X}}g(x)$. Product function $f\cdot g:x\mapsto f(x)\,g(x)$ has Hölder exponent $p$ and a Hölder constant on $I$ which is at most $(m_{f}\,L_{J}(g)+m_{g}\,L_{I}(f))$. 6. 6.

For some countable index set ${\mathcal{I}_{t}}$, let the sequence of functions $f_{i}$ be Hölder with exponent $p$ and constant $L(f_{i})$ each. Furthermore, let $H(x)=\sup_{i\in{\mathcal{I}_{t}}}f_{i}(x)$ and $h(x):=\inf_{i\in{\mathcal{I}_{t}}}f_{i}(x)$ be finite for all $x$. Then $H,h$ are also Hölder with exponent $p$ and have a Hölder constant which is at most $\sup_{i\in{\mathcal{I}_{t}}}L(f_{i})$. 7. 7.

Let $b:=\inf_{x\in\mathcal{X}}|f(x)|>0$ and let $\phi(x)=\frac{1}{f(x)},\forall x\in\mathcal{X}$ be well-defined. Then $L_{I}(\phi)\leq b^{-2}\,L_{I}(f)$. 8. 8.

Let $p=1$ (that is we consider the Lipschitz case), let $I$ be convex and $\,\mathfrak{d}(x,x^{\prime})=\left\|x-x^{\prime}\right\|$ where $\left\|\cdot\right\|$ is a norm that induces a sub-multiplicative matrix norm (e.g. all $p-$ norms are valid). $f$ cont. differentiable on I $\Rightarrow$ $L_{I}(f)\leq\sup_{x\in I}\left\|\nabla f(x)\right\|.$ For one-dimensional input space, $\mathcal{X}=\mathbb{R}$, $L_{I}(f)=\sup_{x\in I}\left|\nabla f(x)\right|$ is the smallest Lipschitz number. 9. 9.

Let $c\in\mathbb{R}$, $f(t)=c,\forall x\in I$. Then $f$ is Hölder continuous with constant $L_{I}(f)=0$ and for any coefficient $p_{f}\in\mathbb{R}$. 10. 10.

$L_{I}(f^{2})\leq 2\,L_{I}(f)\,\sup_{t\in I}f\,$. 11. 11.

With conditions as in 8), and input space dimension one, we have $\forall q\in\mathbb{Q}:L_{I}(f^{q})=|q|\,\sup_{\tau\in\mathcal{I}}|f^{q-1}(\tau)\,\dot{f}(\tau)|$.

Proof.

Refer to [4]. ∎

-A2 Hölder continuity for optimisation with bounds

In the presence of a known Lipschitz constant, optimisation of a Lipschitz continuous function can be done with Shubert’s algorithm [17]. Of importance to us at various parts of the thesis is, that it also provides non-asymptotic error bounds. Unfortunately, it is limited to minimising functions on one-dimensional input domains. For multi-dimensional domains DIRECT [direct:93] is a global optimiser that is popular, not least because no Lipschitz constant is required. On the flip side, it does not provide the error bounds afforded by Shubert’s method. Furthermore, we will sometimes be interested in bounds on the maxima and minima for the more general case, where the functions to be optimised are Hölder continuous with respect to normed spaces. We are currently working on such a method, as well as on a corresponding cubature algorithm. As it is work in progress and its exposition is beyond the scope of this work, we will merely sketch a naive version that provides some desired bounded optimisation for Hölder continuous functions:

For simplicity, assume the input and output space metrics are induced by the maximum norms. For example, $\,\mathfrak{d}_{\mathcal{X}}(x,x^{\prime})=\left\|x-x^{\prime}\right\|_{\infty}$.222Due to norm-equivalence relationships we can bound each norm by each other norm by multiplication with a well-known constant (see e.g. [koenigsberger:2000]). Therefore, knowing a Hölder-constant with respect to one norm immediately yields a constant with respect to another, e.g. $\left\|\cdot\right\|_{\infty}\leq\left\|\cdot\right\|_{2}\leq\sqrt{d}\left\|\cdot\right\|_{\infty}$ where $d$ is the dimension of the input space. Let $f:I\subset\mathbb{R}^{d}\to\mathbb{R}$ be a Hölder continuous function with Hölder constant $L$ and exponent $p$. Let $S=\{(s_{i},f_{i})|i=1,...,N\}$ be a sample of target function $f$ with $f_{i}=f(s_{i})$, $s_{i}\in I,\forall i$.

Assume $I$ is a hyper-rectangle and let $J_{1},...,J_{N}\subset I$ be a partition of the domain $I$ such that each sub-domain $J_{j}$ is a hyperrectangle containing $s_{j}$. (More generally, to avoid $I$ having to be a hyperrectangle itself, one could alternatively assume that the $J_{1},\ldots,J_{N}$ are a covering of $I$. That is, $I\subset\cup_{i}J_{i}$.)

By Hölder continuity $\forall x\in J_{i}:\left|f_{i}-f(x)\right|\leq L\left\|x-s_{i}\right\|_{\infty}^{p}$for some $p\in(0,1],L\geq 0$. Let $a_{1},\ldots,a_{d}>0$, $b_{1},\ldots,b_{d}>0$ be defined such that $J_{i}=s_{i}+([-a_{1},b_{1}]\times...\times[-a_{d},b_{d}])$ where the addition is defined point-wise. Furthermore, let

[TABLE]

be the maximal distance of $s_{i}$ to the boundary of $J_{j}$.

Then, $\forall x\in J_{j}:f(x)\leq f_{j}+L\max_{x\in J_{j}}\left\|x-s_{j}\right\|_{\infty}^{p}=f_{j}+LD_{j}^{p}=:q_{i}$.

Since, the hyperrectangles formed a covering of domain $i$, the desired upper bound on the maximum is $\max_{j=1}^{N}q_{j}\geq\max_{x\in I}f(x)$. By the same argument, a lower bound on the maximum based on the finite sample can be obtained as: $\min_{j=1}^{N}f_{j}-LD_{j}^{p}\leq\min_{x\in I}f(x)$.

In conclusion, we have

[TABLE]

A batch algorithm would get a sample, construct a partition $J_{1},...,J_{N}$ for that sample and then compute the minimum and maximum bounds $\underline{M},\overline{M}$ given above.

For our tests in this work, we have implemented an adaptive algorithm for bounded optimisation which incrementally expands the given sample. That is, it incrementally partitions domain $I$ into increasing small sub-hyperrectangles until the computational budget is expended or the bounds on the maxima and minima have shrunk satisfactorily.

The exploration criterion of where to obtain the next function sample is tailored to finding $\overline{M}$ or $\underline{M}$, depending on the task at hand.

For instance, in a minimisation problem we chose to obtain a new sample point in the hyperrectangle $J_{j_{*}}$ where $j_{*}=\arg\min_{j=1\,...,N}f_{j}-L\,D_{j}^{p}$ is the current rectangle where our present bounds allow for the lowest value to be.