Overfitting Reduction in Convex Regression

Overfitting Reduction in Convex Regression

\def\spacingset#1{ {#1}} \spacingset{1}

\if11 \fi

\if01 {

} \fi

{\it Keywords:} Machine learning, Convex regression, Overfitting, Lipschitz regularization, Weight restriction

\spacingset{1.8}

Introduction

Convex regression (CR), a classical nonparametric regression method with shape constraints dating back to hildreth1954point, has attracted growing interest in operations research lee2013a, balazs2015near, econometrics kuosmanen2008representation, yagi2020shape, statistical learning blanchet2019multivariate, bertsimas2021sparse, and many other areas. One of the great advantages of CR over other nonparametric techniques is that it does not require any tuning parameters (e.g., the bandwidth parameter as in kernel estimation). Consequently, many applications of CR can be found in various fields such as decision analysis lim2012consistency, portfolio selection hannah2013multivariate, productivity analysis kuosmanen2012stochastic2, energy regulation kuosmanen2020conditional, and environmental policy kuosmanen2021shadow.

Even though the CR estimator is a natural candidate for an estimator of convex functions, it produces inconsistent estimates of the underlying true convex function and extremely large subgradients near the boundary ghosal2017univariate. This undesirable overfitting problem has been observed by, e.g., lim2012consistency, mazumder2019a, and liao2024convex. Recently, an empirical application of CR to the estimation of the cost frontier function for Finnish energy distribution firms yielded extremely large subgradients near the boundary kuosmanen2022, which deteriorated the performance of CR in evaluating their incentive regulatory model for electricity firms. Although these studies have shown that the overfitting problem tends to occur in CR, they do not provide any theoretical evidence that explains the existence of overfitting. This motivated us to theoretically investigate the overfitting problem in the framework of CR.

Several attempts have been made to avoid overfitting in the context of CR. For instance, bertsimas2021sparse proposed penalized convex regression (PCR) by adding a penalty term to the objective function of the optimization problem. This formulation allows us to limit the overall behavior of CR, but we cannot directly control the magnitude of each subgradient. mazumder2019a proposed using Lipschitz convex regression (LCR) by adding additional constraints to the subgradient at each data point and, thus, restricting overfitting in a straightforward sense. However, little is known about how LCR behaves on the boundary of the domain and whether it does not exhibit overfitting behavior near the boundary. Indeed, whether placing a hard bound on the subgradients eliminates the overfitting problem remains undetermined. We answer this question affirmatively by proving that the hard bound on the norm of the subgradients generates an estimator that converges to the underlying function uniformly over the domain, including the boundary, with probability one as the sample size increases to infinity. Furthermore, this estimator has subgradients that converge to the gradient of the underlying function uniformly, including the boundary, with probability one as the sample size increases to infinity; see Theorems (ref) and (ref) of this paper.

We further propose two practical guides on how to reduce overfitting by placing a bound on the subgradients. We first proposed an augmented extension of LCR, which originated from the work of mazumder2019a. In the first method, we noticed that we could find a good “reference” for the estimator. For example, one can imagine that the linear regression estimator could be a good starting point without exhibiting any overfitting behavior when estimating a convex function. Alternatively, the percentiles of the subgradients of the convex regression estimator can be another good starting point for the search. The second remedy for the overfitting problem is to place the lower and upper bounds on the subgradients, which are referred to as weight restrictions in the literature (see, e.g., podinovski2016optimal). The easiest way to obtain the hard bounds is to get input from decision-makers. The decision-makers are often able to specify the upper and lower bounds of the subgradients. For example, in the application context where a decision-maker is trying to estimate a cost function, the subgradients are interpreted as the marginal costs, and the decision-maker may have a good understanding of how high the marginal cost can be kuosmanen2021shadow. Thus, the prespecified values of the upper and lower bounds on the subgradients can be readily obtained from the decision-maker. If such information is unavailable, one can rely on the original subgradients obtained from the CR estimator and use their percentiles as the bounds. Adding an explicit bound helps decision-makers understand how the model works, so they trust the model when they use it.

Our main contributions can be summarized as follows.

Theoretical evidence of the overfitting behavior: We study the theoretical properties of the CR estimator near the boundary. Although the consistency of the CR estimator in the interior of the domain has been proven by seijo2011nonparametric and lim2012consistency, its convergence near the boundary of the domain is not guaranteed. Furthermore, the subgradient of the CR estimator is observed to be very large at the boundary mazumder2019a. ghosal2017univariate have shown that the subgradient of the CR estimator is unbounded in probability at the boundary in the univariate case, but no such result has been obtained in the multivariate case. In this paper, we establish this result in the multivariate setting and prove the unboundedness of the subgradient of the CR estimator at the boundary of its domain. To the best of our knowledge, this is the first attempt to systematically study the overfitting behavior of CR under the multivariate setting.

Our proposed estimators: We propose two new estimators, the augmented Lipschitz convex regression (ALCR) and weight-restricted convex regression (WRCR) estimators. Each of them provides a simple and practical way to bound the subgradients of the fitted convex function and, therefore, can alleviate overfitting. We further prove that these estimators and their subgradients are strongly consistent uniformly over the domain. Numerical results in Sections (ref) and (ref) display superior performance of the ALCR and WRCR estimators compared to some existing estimators. While other techniques for avoiding overfitting exist, our proposed estimators provide both theoretical and practical evidence supporting the performance of the proposed hard-bounded approaches.

Application: In Finland, the Energy Authority (EV) has systematically applied CR to implement economic incentives for Finnish electricity distribution firms since 2012 kuosmanen2012stochastic,kuosmanen2020conditional. However, the conventional CR models resulted in high estimated marginal costs (or subgradients) for several big firms; see, e.g., the descriptive statistics in Table 2 of kuosmanen2012stochastic. Moreover, the incentives for these firms in future periods always rely on the predictive performance of the CR model based on historical data. That is, the effects of the incentive program are subject to the predictive power of the CR model. The empirical study shows that we can address the overfitting problem using the proposed estimators and create the right incentives for regulated firms.

The rest of this paper is organized as follows. In Section (ref), We introduce the problems of convex regression and overfitting reduction, along with relevant notation and definitions. In Section (ref), we provide some theoretical evidence on the overfitting behavior of the CR estimator. In Section (ref), we introduce the proposed ALCR estimator and investigate its statistical properties. The WRCR estimator and its properties are developed and analyzed in Section (ref). In Section (ref), we perform Monte Carlo studies to compare the performance of the proposed estimators to that of some existing estimators. An empirical application of our proposed methods to Finnish electricity distribution firms is presented in Section (ref). In Section (ref), we conclude this paper with suggestions on future research avenues.

Preliminaries

Convex regression

CR can be formally described as follows. Given data $(\boldsymbol{x}_1,y_1), \cdots, (\boldsymbol{x}_n, y_n)$, we assume

for $i = 1, \cdots, n$, where $\boldsymbol{x}_i \in \Omega \subset \mathbb{R}^d$, $f_0:\Omega \rightarrow \mathbb{R}$ is the underlying function to be estimated, and the $\varepsilon_i$'s are the error terms satisfying $\mathbb{E}[\varepsilon_i]=0$ and $\text{Var}[\varepsilon_i]=\sigma^2 < \infty$ for $i = 1, \cdots, n$. Even though $f_0$ cannot be observed exactly, it is known to be convex. In the context of CR, our goal is to estimate $f_0$ by fitting a convex function to the given $n$ observations. In particular, we estimate $f_0$ by minimizing the sum of squared errors:

over $f \in \mathcal{F}$, where $\mathcal{F} \triangleq \{f:\Omega\rightarrow \mathbb{R} \mbox{ such that } f \mbox{ is convex}\}$ is the class of all real-valued convex functions over $\Omega$.

In ((ref)), we try to fit a convex function $f$ to the data points $(\boldsymbol{x}_1, y_1), \cdots, (\boldsymbol{x}_n, y_n)$ and find the one with the least squares, so ((ref)) is an optimization problem over functions and it appears to be infinite-dimensional. However, the infinite-dimensional problem ((ref)) can be reduced to a finite-dimensional quadratic program where the decision variables are the function values and subgradients of the fitted convex function $f$ at the $\boldsymbol{x}_i$'s. The following quadratic program, in the decision variable $f_1, \cdots, f_n \in \mathbb{R}$ and $\boldsymbol{\beta}_1, \cdots, \boldsymbol{\beta}_n \in \mathbb{R}^d$, is one such formulation:

see, e.g., seijo2011nonparametric and lim2012consistency. In ((ref)), $f_i$ and $\boldsymbol{\beta}_i$ represent the value of the fitted function $f$ and a subgradient of $f$, respectively, for $i= 1, \cdots, n$. The set of constraints in ((ref)), namely shape constraints, enforces the fitted function $f$ to be convex.

In the cases of production/cost function estimation, an alternative formulation is commonly used in the literature as the returns to scale need to be specified; see, e.g., afriat1972efficiency, banker1993maximum, and kuosmanen2010data. Formally, the following quadratic program, in the decision variables $\alpha_1, \cdots, \alpha_n \in \mathbb{R}$, and $\boldsymbol{\beta}_1, \cdots, \boldsymbol{\beta}_n$ $\in \mathbb{R}^d$, can be used:

see kuosmanen2008representation for details. In ((ref)), $\alpha_i$ is a constant representing the intercept and $\boldsymbol{\beta}_i$ denotes the slope coefficients (or subgradient of $f$), at $\boldsymbol{x}_i$ for $i= 1, \cdots, n$. By solving ((ref)) or ((ref)), we can numerically compute a solution to ((ref)), which is referred to as the CR estimator. For simplicity of notation, we will use the formulation in (ref) in the following text, and we will consider problem (ref) in Section (ref) for cost function estimation.

Overfitting reduction

Overfitting is a common issue in machine learning where the fitted models are overly complex and fail to capture the underlying true patterns lever2016points. This leads to a model that performs perfectly on training data but poorly on unseen data, resulting in poor out-of-sample performance. For example, high-order polynomial regression is more flexible than its low-order counterpart and, hence, is more prone to suffer from the overfitting problem.

To avoid overfitting in CR, we often need to restrict the model flexibility so that the estimator does not lead to an extreme model with large subgradients, especially near the boundary. A widely used way to do so is to add the penalty (regularization) term to the optimization problem (ref). For example, bertsimas2021sparse proposed PCR by adding $L_2$-norm regularization (see also in lim2021consistency). That is, we minimize the sum of squared errors and squared norms of the subgradients as follows:

over $f \in \mathcal{F}$ for a given sequence of nonnegative real numbers $(\lambda_n: n \geq 1)$. The infinite-dimensional problem ((ref)) can be reformulated as the following finite-dimensional quadratic program in the decision variables $f_1, \cdots, f_n \in \mathbb{R}$ and $\boldsymbol{\beta}_1, \cdots, \boldsymbol{\beta}_n \in \mathbb{R}^d$:

The solution $\hat{f}_{1, p}, \cdots, \hat{f}_{n, p}, \hat{\boldsymbol{\beta}}_{1, p}, \cdots, \hat{\boldsymbol{\beta}}_{n, p}$ to ((ref)) exists uniquely lim2021consistency. We define $\hat{f}_{n, p}:\Omega \rightarrow \mathbb{R}$ by

for $\boldsymbol{x} \in \Omega$ and refer to $\hat{f}_{n, p}$ as the PCR estimator.

Instead of imposing $L_2$-norm on all subgradients $\sum_{i = 1}^n \|\boldsymbol{\beta}_i\|^2$, mazumder2019a apply the $L_2$-norm Lipschitz regularization on each subgradient to manage the overfitting problem. For a given constant $L>0$, LCR is then formulated as

Similarly, we can define a representation function like (ref) by replacing the solution $\hat{f}_{1}, \cdots, \hat{f}_{n}, \hat{\boldsymbol{\beta}}_{1}, \cdots, \hat{\boldsymbol{\beta}}_{n}$ with those obtained from problem (ref).

Leveraging the idea of placing hard bounds on the subgradients in mazumder2019a, we first consider an extension of the $L_2$-norm Lipschitz regularization in the current paper. That is, we add a reference to the $L_2$-norm such that we aim to regularize the subgradients to the benchmark reference. One can use these “reference” values and add the inequality

for $i = 1, \cdots, n$ to the constraints of ((ref)), where $\boldsymbol{b}_0$ represents the reference vector mentioned so far. A key advantage of the new constraints is that they provide a practical way of integrating the prior knowledge of decision-makers (see an empirical study for illustration in Section (ref)).

Another practical way to impose a hard bound on the subgradients is adding upper and lower bounds to each subgradient directly as follows:

for $n = 1, \cdots, n$, where $\boldsymbol{l}_0 \in \mathbb{R}^d$ and $\boldsymbol{u}_0 \in \mathbb{R}^d$ indicate lower and upper bounds on the subgradient, respectively, and $\leq$ is used componentwise. We will revisit these two hard bounds and investigate their statistical properties when applying to CR estimators in Sections (ref) and (ref).

Notations

For $\boldsymbol{x} \in \mathbb{R}^d$, $x(k)$ represents its $k$-th element for $k = 1, \cdots, d$, so $\boldsymbol{x} = (x(1), \cdots, x(k))$. We write $\|\boldsymbol{x}\| = \left\{x(1)^2 + \cdots + x(d)^2\right\}^{1/2}$. The transpose of $\boldsymbol{x} \in \mathbb{R}^d$ is denoted by $\boldsymbol{x}^T$. $\mathbf{0} \in \mathbb{R}^d$ denotes the vector with zeros in all entries. For $\boldsymbol{a}, \boldsymbol{b} \in \mathbb{R}^d$, we write $\boldsymbol{a} \leq \boldsymbol{b}$ if and only if $a(i) \leq b(i)$ for $i = 1, \cdots, d$. We also write $\boldsymbol{a} < \boldsymbol{b}$ if and only if $a(i) < b(i)$ for $i = 1, \cdots, d$.

Let $\Lambda\subset \mathbb{R}^d$ be a convex set. For a convex function $f:\Lambda \rightarrow \mathbb{R}$, we call $\boldsymbol{\beta} \in \mathbb{R}^d$ a subgradient of $f$ at $\boldsymbol{x} \in \Lambda$ if $\boldsymbol{\beta}^T(\boldsymbol{y} - \boldsymbol{x}) \leq f(\boldsymbol{y}) - f(\boldsymbol{x})$ for any $\boldsymbol{y} \in \Lambda$. We denote a subgradient of $f$ at $\boldsymbol{x} \in \Lambda$ by $\mbox{subgrad }f(\boldsymbol{x})$. We call the set of all subgradients at $\boldsymbol{x}$ the subdifferential at $\boldsymbol{x}$, and denote it by $\partial f(\boldsymbol{x})$. For any differentiable function $f:\Lambda \rightarrow \mathbb{R}$, $\nabla f(\boldsymbol{x})$ denotes the derivative of $f$ at $\boldsymbol{x} \in \Lambda$.

The overfitting problem

In this section, we investigate the behavior of the CR estimator near the boundary of its domain and theoretically explain why it shows the overfitting behavior.

We start by formally defining the CR estimator. In CR, we minimize the sum of squared errors in ((ref)) over $f \in \mathcal{F}$. Problem ((ref)) determines the values of the fitted convex function $f$ at the $\boldsymbol{x}_i$'s (i.e., $f_i$'s) uniquely, but it does not determine the subgradients of the fitted function at the $\boldsymbol{x}_i$'s (i.e., $\boldsymbol{\beta}_i$'s) uniquely; see Lemma 2.5 of seijo2011nonparametric. To define the CR estimator uniquely over $\Omega$, we take, among all the $\boldsymbol{\beta}_i$'s solving problem ((ref)), the one with the minimum norm. More precisely, to determine the subgradients of the fitted function at the $\boldsymbol{x}_i$'s uniquely, we solve the following finite-dimensional quadratic program in the decision variables $\boldsymbol{\beta}_1, \cdots, \boldsymbol{\beta}_n \in \mathbb{R}^d$:

where $\hat{f}_1, \cdots, \hat{f}_n$ are the minimizing values of ((ref)). The solution $\tilde{\boldsymbol{\beta}}_1, \cdots, \tilde{\boldsymbol{\beta}}_n$ to ((ref)) exists uniquely since ((ref)) is a minimization problem with continuous, strictly convex, and coercive objective function over a non-empty and closed set (Proposition 7.3.1 and Theorem 7.3.7 in kurdila2006convex). We now define $\hat{f}_n: \Omega \rightarrow \mathbb{R}$ by

for $\boldsymbol{x} \in \Omega$ and refer to $\hat{f}_n$ as the CR estimator.

To analyze the properties of the CR estimator near the boundary of $\Omega$, we need the following assumptions:

Now, we state the main results of this section. Theorem (ref) shows that assuming the convexity of $f_0$, the CR estimator is inconsistent in estimating $f_0$ near the boundary of $\Omega$ as $n \rightarrow \infty$. Theorem (ref) states that assuming the differentiability of $f_0$ over $\Omega$, the subgradients of the CR estimator are unbounded in probability near the boundary of $\Omega$ as $n \rightarrow \infty$. The detailed proofs for Theorems (ref) and (ref) are available in the supplementary material to this paper.

The unboundedness property in Theorem (ref) suggests that the estimated subgradients of CR could take very large values near zero, which is one of the boundary points. In the following sections, we propose restricting the estimated subgradients with a known Lipschitz bound (Section (ref)) and with known upper and lower bounds (Section (ref)).

Augmented Lipschitz convex regression

In this section, we propose the ALCR estimator, which bounds the norm of the subgradients around a reference vector $\boldsymbol{b}_0$. This method becomes useful when some prior knowledge on the subgradients of $f_0$ is available. We will then show that the ALCR estimator does not show the overfitting behavior by proving its uniform consistency over the entire domain $\Omega$.

Given $\boldsymbol{b}_0 \ge \mathbf{0}$ and $L_0 > 0$, we consider a class of convex functions whose subgradients are uniformly bounded by $L_0$ around $\boldsymbol{b}_0$:

In our ALCR, we minimize the sum of squared errors:

over $f \in \mathcal{F}_A$. The infinite-dimensional problem ((ref)) can be reduced to the following finite-dimensional convex program in the decision variables $f_1, \cdots, f_n \in \mathbb{R}$ and $\boldsymbol{\beta}_1, \cdots, \boldsymbol{\beta}_n \in \mathbb{R}^d$:

The solution $\hat{f}_{1, A}, \cdots, \hat{f}_{n, A}, \hat{\boldsymbol{\beta}}_{1, A}, \cdots, \hat{\boldsymbol{\beta}}_{n, A}$ to ((ref)) exists because ((ref)) is an optimization problem with continuous and coercive objective function over a nonempty convex subset (Proposition 7.3.1 and Theorem 7.3.7 in kurdila2006convex). Furthermore, the minimizing values $\hat{f}_{1, A}, \cdots, \hat{f}_{n, A}$ are unique because the objective function is strictly convex, but the $\hat{\boldsymbol{\beta}}_{i, A}$'s are not unique. We now define $\hat{f}_{n, A}: \Omega \rightarrow \mathbb{R}$ by

for $\boldsymbol{x} \in \Omega$ and refer to $\hat{f}_{n, A}$ as the ALCR estimator.

The following theorem, Theorem (ref), establishes the strong uniform consistency of the ALCR estimator and its subgradients over the entire domain $\Omega$ as $n \rightarrow \infty$. The proof of Theorem (ref) is provided in the supplementary material to this paper.

Theorem (ref) ensures that the ALCR estimator and its subgradients converge uniformly on $\Omega$ to $f_0$ and $\nabla f_0$, respectively, with probability one as $n \rightarrow \infty$. This shows that the overfitting behavior is successfully eliminated in the ALCR estimator.

Weight-restricted convex regression

In this section, we propose the WRCR estimator by imposing explicit upper and lower bounds to the subgradients. This will help avoid very large values of the subgradients estimated from CR.

Given $\boldsymbol{l}_0 \in \mathbb{R}^d$ and $\boldsymbol{u}_0 \in \mathbb{R}^d$, we consider the class of convex functions whose subgradients are bounded between $\boldsymbol{l}_0$ and $\boldsymbol{u}_0$:

In our proposed WRCR, we minimize the sum of squared errors:

over $f \in \mathcal{F}_B$. The infinite-dimensional problem ((ref)) can be reduced to the following finite-dimensional convex program in the decision variables $f_1, \cdots, f_n \in \mathbb{R}$ and $\boldsymbol{\beta}_1, \cdots, \boldsymbol{\beta}_n \in \mathbb{R}^d$:

The solution $\hat{f}_{1, B}, \cdots, \hat{f}_{n, B}, \hat{\boldsymbol{\beta}}_{1, B}, \cdots, \hat{\boldsymbol{\beta}}_{n, B}$ to ((ref)) exists because ((ref)) is an optimization problem with continuous and coercive objective function over a nonempty convex subset (Proposition 7.3.1 and Theorem 7.3.7 in kurdila2006convex). Furthermore, the minimizing values $\hat{f}_{1, B}, \cdots, \hat{f}_{n, B}$ are unique because the objective function is strictly convex, but the $\hat{\boldsymbol{\beta}}_{i, B}$'s are not unique. We now define $\hat{f}_{n, B}: \Omega \rightarrow \mathbb{R}$ by

for $\boldsymbol{x} \in \Omega$ and refer to $\hat{f}_{n, B}$ as the WRCR estimator.

The following theorem, Theorem (ref), establishes the strong uniform consistency of the WRCR estimator over the entire domain $\Omega$. The proof of Theorem (ref) is provided in the supplementary material to this paper.

Theorem (ref) ensures that the WRCR estimator and its subgradients converge uniformly on $\Omega$ to $f_0$ and $\nabla f_0$, respectively, with probability one as $n \rightarrow \infty$. This shows that the overfitting behavior is successfully eliminated in the WRCR estimator.

Monte Carlo study

The goal of this section is to observe the numerical behavior of our proposed ALCR and WRCR estimators compared to some existing estimators. We start with Section (ref), which illustrates the CR estimator's overfitting behavior numerically.

Illustration of the overfitting behavior

In this section, we observe the overfitting behavior of the CR estimator numerically. We consider $f_0:[0.2, 1.0] \rightarrow \mathbb{R}$ defined by $f_0(\boldsymbol{x}) = 1/\boldsymbol{x}$ for $\boldsymbol{x} \in [0.2, 1.0]$. The $\boldsymbol{x}_i$'s are evenly distributed, so we let $\boldsymbol{x}_i = 0.2 + 0.8i/n -0.8/(2n)$ for $i = 1, \cdots, n$. We generate the $y_i$'s from $y_i = f_0(\boldsymbol{x}_i) + \varepsilon_i$ for $i = 1, \cdots, n$, where the $\varepsilon_i$'s are i.i.d. normal random variables with mean 0 and variance 1. Once the $(\boldsymbol{x}_i, y_i)$'s are obtained, we compute the CR estimator $\hat{f}_n$ by solving ((ref)) and ((ref)) with Mosek.

To observe how the CR estimator behaves at the boundary of $\Omega = [0.2, 1.0]$, we compute

We repeat this 100 times independently, generating 100 replications of ((ref)), and use these 100 values to compute the 95% confidence interval of $\mathbb{E}[\max_{1 \leq i \leq n} |\hat{f}_n(\boldsymbol{x}_i) - f_0(\boldsymbol{x}_i)|]$. Table (ref) reports these 95% confidence intervals for a wide range of $n$.

Another way to observe the overfitting behavior is by looking at the accuracy of the estimated subgradients. We thus compute

repeat this 100 times independently, generating 100 replications of ((ref)), and use these 100 values to compute the 95% confidence interval of $\mathbb{E} [\max_{1 \leq i \leq n} |\mbox{subgrad }\hat{f}_n(\boldsymbol{x}_i) - \nabla f_0(\boldsymbol{x}_i)|]$. Table (ref) reports these 95% confidence intervals for a wide range of $n$.

The results in Table (ref) show that the CR estimator produces inconsistent estimators of $f_0$ near the boundary of $\Omega$ as $n$ increases. It also generates subgradients whose magnitudes are increasing to infinity as $n \rightarrow \infty$.

Comparisons of statistical performance

This section is concerned with the numerical performance of our proposed ALCR and WRCR estimators compared with three existing estimators, CR, PCR, and LCR estimators. The CR estimator is defined in ((ref)). The LCR estimator is a special case of the ALCR estimator ((ref)) when $\boldsymbol{b}_0 = \boldsymbol{0} \in \mathbb{R}^d$.

Setup

Consider the following two test functions:

where $y_i = f_0(\boldsymbol{x}_i)+\varepsilon_i$ for $i=1,\cdots,n$. We generate the $\boldsymbol{x}_i$'s, independently of one another, from the uniform distribution over $[1, 10]^d$ and draw the $\varepsilon_i$'s from the normal distribution with mean 0 and variance $\sigma^2$, where $\sigma^2$ is determined by the signal-to-noise ratio (SNR), $\text{SNR}= \text{Var}[f_0 (\boldsymbol{x}_1)]/\sigma^2$.

In all simulations, we generate $2n$ independent observations and equally split them into the training and validation sets. We first use the training set to compute the CR estimator and its subgradients $\hat{\boldsymbol{\beta}}_i$'s by solving ((ref)). Next, we use the validation set to compute $\boldsymbol{b}_0$ and all the tuning parameters. We find $\boldsymbol{b}_0$ by computing the slope of the linear regression estimator with the validation set. The tuning parameters we need to compute from the validation set are $L$ for the LCR estimator, $\lambda_n$ for the PCR estimator, $L_0$ for the ALCR estimator, and $q$ (defined in Remark (ref)) for the WRCR estimator. We use the 5-fold cross-validation method to compute these tuning parameters. For example, to find $\lambda_n$ for the PCR estimator using the 5-fold cross-validation method, we split the validation set into five equally sized sets, say $S_1, S_2, \cdots, S_5$. For each $\lambda \in \mathbb{R}$, we define the following cross-validation function: \[\mbox{CV}(\lambda) \triangleq \frac{1}{5} \sum_{k = 1}^5 \sum_{(\boldsymbol{x}_j, y_j) \in S_k} (y_j - f_{\lambda, n}^{[k]}(\boldsymbol{x}_j))^2,\] where $f_{\lambda, n}^{[k]}$ is the PCR estimator in ((ref)) and ((ref)) computed with $\lambda_n = \lambda$ and the $(\boldsymbol{x}_i, y_i)$'s replaced by $\{(\boldsymbol{x}_j, y_j) \;| \; (\boldsymbol{x}_j, y_j) \mbox{ is in the validation set, but } (\boldsymbol{x}_j, y_j) \notin S_k\}$. We then evalulate $\mbox{CV}(\lambda)$ for each value $\lambda$ in the set of candidate values, say $\{\lambda_1, \cdots, \lambda_m\}$, and select the value with the least value of $\mbox{CV}(\lambda)$. We find $L$, $L_0$, and $q$ using the 5-fold cross-validation method in a similar fashion. For the set of candidate values, we use $\{0.01, 0.02, \cdots, 0.49\}$ to search for $q$. To search for $L_0$, we use 50 equally spaced values from $[0,\overline{L}]$, where $\overline{L} = \max_{1 \leq i \leq n}\|\hat{\boldsymbol{\beta}}_i - \boldsymbol{b}_0\|$. For $L$, we use 50 equally spaced values from $[0.1, 5.0]$.

In the following experiments, we use the standard solver Mosek (9.2.44) within the Julia/JuMP package to compute optimization problems. All computations are performed on a computing cluster with Xeon @2.8 GHz processors, 10 CPUs, and 8 GB RAM per CPU.

Statistical performance

We first consider the Type A function with $n=\{50, 100, 200\}$ and $\text{SNR}=3$ to compare our proposed estimators to the CR, PCR, and LCR estimators.

As described in Section (ref), we generate the training and validation sets and use them to compute all the tuning parameters. Next, we generate 1000 independent replications of $(\boldsymbol{x}_1, y_1)$, say $((\boldsymbol{x}_i, y_i): 1 \leq i \leq 1000)$ and use the data along with the tuning parameters to compute prediction errors of the LCR, PCR, ALCR, and WRCR estimators.

To measure the accuracy of the ALCR estimator $\hat{f}_{n, A}$, we compute the empirical mean square error (MSE) value as follows:

We repeat this 50 times independently, generating 50 replications of ((ref)). Table (ref) reports the 95% confidence intervals of the expected MSE that are computed using these 50 values for various values of $n$. We repeat this procedure for each of the PCR, LCR, and WRCR estimators and report the confidence intervals of their expected MSEs in Table (ref).

Although all of the PCR, LCR, ALCR, and WRCR estimators seem to benefit from the additional restrictions on subgradients, our ALCR and WRCR estimators yield more significant improvements in the predictive power over the PCR and LCR estimators.

We next repeat this procedure with the Type B function, generating 50 replications of the MSE value for each of the LCR, PCR, ALCR, and WRCR estimators, and create box plots each with 50 MSE values. In Figure (ref), each box plot shows the median (in the middle of each box), the first quantile (lower end of each box), the third quartile (higher end of each box), and outliers as small circles. Figure (ref) shows the box plots for a wide range of $d$ when $n = 100$ and $n = 500$. Our proposed estimators, the ALCR and WRCR estimators, perform better than other estimators in most cases, especially when the dimension is higher, e.g., $d = 10$.

We next use the Type B function to explore how the noise level affects the predictive performance of the proposed estimators. We repeat the previous experiment with the Type B function. Figure (ref) reports the box plots each with 50 MSE values for a wide range of SNR$\in\{1.0, 1.5, 2.0, 2.5, 3.0\}$ when $n = 100, d = 3$ and when $n=500, d = 10$.

LCR performs similarly to ALCR when the sample size is small ($n=100, d=3$). However, the new method shows better predictive power in large sample settings ($n=500, d=10$). In these settings, WRCR also performs reasonably well compared to other benchmarks but does not show such an advantage at high noise levels. The supplementary material to this paper provides additional numerical results with different underlying functions $f_0$.

Regulation on electricity distribution firms in Finland

In Finland, electricity distribution firms typically enjoy a natural local monopoly due to prohibitively expensive construction fees. This forces governments to establish regulatory agencies to monitor the electricity distribution firms, reduce their local monopoly power, and provide incentives to those adopting the best technology. The EV is one of the pioneers in practically implementing energy regulation programs, which involve cost-efficiency analysis for distribution firms. The traditional frontier estimation techniques, such as data envelopment analysis and stochastic frontier analysis, were applied by EV from 2005 to 2011, and a more reliable CR model was adopted from 2012 onward kuosmanen2012stochastic. Recently, EV considered developing new techniques to reduce the overfitting behavior in CR and increase the accuracy of its incentive regulatory model for Finnish electricity distribution firms kuosmanen2022. In our analysis, we examine the performance of our proposed estimators on the cost frontier estimation problem for Finnish electricity distribution firms.

Frontier estimation model

We consider a general cost frontier model in production economics kuosmanen2020conditional:

for $i = 1, \cdots, n$, where $C$ is a nondecreasing and convex cost function, and the $\varepsilon_i$'s are the error terms representing random effects. $\boldsymbol{x}_i$ is the vector of performance variables for firm $i$, and $y_i$ is the variable cost of firm $i$. The variables are specified as follows:

We collect data from the EV over three regulation periods with 86 Finnish electricity distribution firms, covering the years from 2008 to 2020. The total number of $(\boldsymbol{x}_i, y_i)$ pairs in our data set is $1042$, with some companies going bankrupt after 2010. The original dataset is applied to build the incentive regulation model for Finnish electricity distribution firms kuosmanen2022, and its earlier version has been widely used; see kuosmanen2012stochastic, kuosmanen2013what, and kuosmanen2020conditional for details.

To estimate the cost function $C$ in (ref), we consider the following variant of CR that incorporates the monotonicity of $C$:

The last constraint in (ref) guarantees that the estimated cost function is non-decreasing. The monotonicity assumption on production/cost function has been identified in various areas such as operations research and economics; see, e.g., kuosmanen2010data and kuosmanen2012stochastic2. However, the additional monotonicity constraint in (ref) does not eliminate the overfitting behavior.

Using our data set $((\boldsymbol{x}_i, y_i): 1 \leq i \leq 1042)$, we compute the CR estimator and its subgradients $\hat{\boldsymbol{\beta}}_i = (\hat{\beta}_i(1), \hat{\beta}_i(2), \hat{\beta}_i(3), \hat{\beta}_i(4))$ at $\boldsymbol{x}_i$. Table (ref) reports some summary statistics, such as the mean, standard deviation, minimum, 10th percentile, median, 90th percentile, and maximum, for $(\hat{\beta}_i(k): 1 \leq i \leq 1042)$ for each $k \in \{1, 2, 3, 4\}$.

For each variable, the range between the minimum and maximum estimated subgradients is quite large. Those large subgradients usually occur near the boundary of the domain and, hence, lead to the overfitting behavior in the CR estimator. The overfitting problem may yield invalid results in frontier estimation and produce misleading benchmarks for the acceptable level of cost in future regulation periods.

The 90th percentiles indicate that the estimated subgradients within the 90th percentile should not be very large. We notice that the changes in the estimated subgradients are much smaller within the range between the 10th and 90th percentiles compared to the changes in the range between the minimum and the maximum. This difference is caused by a few big firms with extreme data values that lead to overfitting behavior in the estimation of the cost function using CR.

The descriptive statistics in Table (ref) provide some information that we can use in our proposed estimators. For example, the median values in Table (ref) are natural candidates for $\boldsymbol{b}_0$ in the ALCR estimator.

Evaluating the predictive performance

For the next regulation period between 2024 and 2031, improving the predictive accuracy of the cost frontier model is one of the main concerns of EV regulators. In this section, we compare the two proposed estimators with several key benchmark methods that have been used on the Finnish electricity distribution firm data. For this purpose, we need to slightly modify the CR, PCR, ARCR, and WRCR estimators. The CR estimator (ref) is a variant of the CR estimator where the monotonicity of the underlying function is incorporated as $\boldsymbol{\beta}_i \geq 0$ for $i = 1, \cdots, n$ in the constraint. We also add these conditions to the constraints of PCR's, LCR's, ARCR's, and WRCR's formulations.

In all computations, we use the observations in the first 9 years $((\boldsymbol{a}_i, b_i): i = 1, \cdots, 734)$ as the training data set and the observations in the last 4 years $((\boldsymbol{c}_i, d_i): i = 1, \cdots, 308)$ as the test data set. We set $\boldsymbol{b}_0 = (1.8, 0.5, 34.3, 9.1)$ as the median values reported in Table (ref). We use the training data set to compute all other tuning parameters, $\lambda$ in PCR, $L$ in LCR, $L_0$ in ARCR, and $q$ (defined in Remark (ref)) in WRCR. We use the 5-fold cross-validation method with the following set of candidate values: 50 equally spaced values from $[1, 500]$ for each of $\lambda_n$, $L$ and $L_0$, and $\{0.01,0.02, \cdots, 0.49\}$ for $q$. Next, we use the test data set to measure the accuracy of each estimator. For example, to measure the accuracy of the CR estimator, we define in-sample root mean square error (in-sample RMSE) as follows: \[\left\{\frac{1}{734}\sum_{i = 1}^{734} (b_i - \mbox{ CR at }\boldsymbol{a}_i)^2\right\}^{1/2}\] and out-of-sample root mean square error (out-of-sample RMSE) as follows: \[\left\{\frac{1}{308}\sum_{i = 1}^{308} (d_i - \mbox{CR at }\boldsymbol{c}_i)^2\right\}^{1/2}.\]

Table (ref) reports the in-sample and out-of-sample RMSEs for various estimators.

The CR estimator results in a severe overfitting problem, where the in-sample RMSE is very small while the out-of-sample RMSE appears to be very large. The PCR, LCR, and our proposed estimators (ALCR and WRCR) yield results that are significantly better than those of CR, demonstrating that an additional bound on the subgradients helps improve the predictive performance of CR. We observe that the proposed estimators work best in terms of the prediction accuracy. The WRCR estimator achieves the smallest out-of-sample RMSEs, indicating that it outperforms other CR estimators in terms of predictive power.

From the regulator's perspective, it can be hard to interpret the tuning parameter $\lambda_n$ in PCR and $L$ in LCR. This makes the model hard to be trusted by practitioners. By contrast, the tuning parameter $L_0$ in ALCR can be interpreted as the distance between the estimated marginal costs and the reference vector $\boldsymbol{b}_0$ (the median values of the subgradients of the CR estimator). The $q$th and $(1-q)$th percentiles used in the WRCR can be seen as the proportion of extreme values obtained from CR that the regulators would like to exclude from the subsequent estimation.

Conclusions

In this paper, we provide theoretical evidence of the overfitting behavior of the CR estimator near the boundary of the domain. To eliminate the overfitting behavior, we consider placing bounds on the subgradients of the fitted convex function. In particular, we propose two practical ways to place the bounds, leading our proposed estimators, the ALCR and WRCR estimators. Monte Carlo studies have shown the superior performance of our proposed estimators.

This work was initially motivated by the observed overfitting problem in the regulatory models previously used by the Finnish EV and the need to improve the prediction accuracy of their models in the 5th and 6th regulation periods from 2024 to 2031. The two newly proposed estimators allow EV to have a model that does not show the overfitting behavior and that outperforms the existing CR estimators.

There are several open research avenues. In this paper, we conducted numerical experiments in the cases where $d\ll n$. Understanding and investigating the behavior of the CR estimator when the number of variables is large could be a direction of interest. Furthermore, the predictive performance of our approaches relies on the choice of tuning parameters. It would be interesting to develop a tuning-free CR model that does not show overfitting behavior.

Acknowledgements

The authors acknowledge the computational resources provided by the Aalto Science-IT project. Zhiqiang Liao gratefully acknowledges financial support from the HSE Support Foundation [grant no. 18-3419] and the Jenny and Antti Wihuri Foundation [grant no. 00220201]. Sheng Dai gratefully acknowledges financial support from the OP Group Research Foundation [grant no. 20230008].

\baselineskip 20pt