Generalization-Aware Structured Regression towards Balancing Bias and Variance

1 Introduction

One of the fundamental challenges in machine learning is to develop models that can learn from empirical evidence and make accurate predictions. Such a challenge requires a design of learning algorithms capable of producing hypotheses (models) that 1) assimilate the empirical evidence, and 2) generalize well to unobserved data. The former is controlled by the notion of empirical risk Remp (training error), while the latter depends on the generalization risk Rgen = |Remp - Rtrue| which determines whether Remp is a valid estimate of the true unknown risk Rtrue (test error). A situation of high empirical risk Remp and low generalization risk Rgen causes underfitting, characterized by the presence of high bias. On the other hand, the converse leads to overfitting indicating high variance. Therefore, in mathematical terms, the problem comes down to minimizing Remp, while maintaining low Rgen. Although, minimizing the empirical risk can be easily attained since it is "measurable" from the observed data, we are not aware of the conditions under which a learning algorithm generalizes. The generalization risk is often impossible to determine since the true risk is unknown. However, there is a broad range of established upper bounds on Rgen for both deterministic and randomized learning algorithms, and in both regression and classification cases. Initially introduced by Vapnik [Vapnik, 1999], generalization risk bounds have been derived on the basis of uniform convergence [Vapnik, 1999], algorithmic stability [Elisseeff et al., 2003; Elisseeff et al., 2005], generic chaining [Audibert and Bousquet, 2007; Talagrand, 1996], the PAC-Bayesian framework [Audibert and Bousquet, 2007; McAllester, 2003], Rademacher and Gaussian complexities [Bartlett and Mendelson, 2002], and robustness-based analysis [Xu and Mannor, 2012]. Some of these bounds are derived on the hypotheses selected by learning algorithms, while others bound the risks of the learning algorithms. For instance, the Vapnik-Chervonenkis (VC) theory [Blumer et al., 1989] provides generalization bounds on hypotheses' risks, while stability bounds [Bousquet and Elisseeff, 2002; Kutin and Niyogi, 2002; Poggio et al., 2004] are derived on learning algorithms' risks. For the theoretical justification in this study, we will focus on the latter.

Inspired by this insight, we designed a bias-variance balancing objective function that aims at tightening an algorithmic stability-based upper bound on the expected true risk in order to yield improved predictive performance. It utilizes distance correlation [Sz´ekely et al., 2007; Sz´ekely et al., 2009], a measure for statistical dependence, to control mutual stability between a model's loss w.r.t. given data and the data itself. This enables learning from empirical evidence, while at the same time enhancing the overall model's generalization performance. Since ensemble construction is a natural way of achieving greater generalization performance, we present the GeneraLization-Aware Collaborative Ensemble Regressor (GLACER), a collaborative ensemble-based model for structured regression that optimizes a generalization-aware objective function. GLACER bags multiple graphical models, namely Gaussian CRFs, and allows them to interact by exchanging examples in a way that minimizes the proposed objective function. Here, we utilize the idea of exchanging examples between ensemble components in a structured regression setting. This decreases GLACER's empirical risk and at the same time accounts for its generalization risk by tightening an upper bound of its expected true risk.

GLACER has been assessed on both synthetic and realworld data, and compared against both traditional and structured regression models. In the conducted experiments, GLACER's predictions have shown to be stable, yielding statistically significant improvements in average MSE.

GLACER's underlying framework can be viewed as a collaborative ensemble that aims to minimize a generalization error bound. Thus, it can be extended beyond regression by replacing its base components by any kind of learning models (e.g. binary classification models [Arsov et al., 2017]).

The main contributions of this work are summarized as follows:

1. Theory and calculation of a stability-based generalization error bound are bridged by leveraging the distance correlation measure;
2. A bias-variance balancing objective function that takes advantage of distance correlation is designed to indirectly control the mutual stability between a regression model's loss on a dataset and the dataset itself, and thus to effectively address the trade-off between underfitting and overfitting;
3. A collaborative ensemble regressor was developed whose components are allowed to interact through a novel exampleexchange-driven optimization in order to optimize the proposed objective.

2 ProblemFormulation

Let Z = X ×Y be an observation space of input-output pairs, where X ∈ Rd and Y ∈ R denote the input and output space, respectively. Further, let D = {z1 = (x1, y1), . . . , zN = (xN, yN)} ∈ Z N be a training set of N input-output pairs (xi, yi) referred to as examples, where Z N represents the space of all training sets of size N. The task is the prediction of an unobserved vector y′ = [y′1 , . . . , y′N ′ ]⊤ of real-valued outputs (targets), given their corresponding inputs X′ = [x′1 , . . . , x′N ′ ]⊤.

3 Main Theoretical Insight

The proposed objective function optimized by GLACER (introduced later in Section 4.2) relies on the following insights:

• An upper bound, derived on the expected true risk ˆRtrue, i.e. the true risk of any learning algorithmL:

where Remp(h, D) denotes the empirical risk of a hypothesis h (inferred by L), while S(ℓ(·, h), ztrn) represents the mutual stability between the loss of h and a random training example ztrn.

• A proposed bias-variance balancing objective function designed to minimize the upper bound in Eq. (1):

where Remp(h,D) is the empirical risk of h and dCorr(ℓ(·, h), ztrn) is the distance correlation, a measure of statistical dependence between ℓ(·, h) and ztrn.

Next, we discuss the derivations of Eqs. (1) and (2) in detail.

3.1 Theoretical Background

The main goal is to address the fundamental learning problem of balancing between the notions of underfitting and overfitting. The former is controlled by the empirical risk Remp, while maintaining a low true risk Rtrue assists in avoiding the latter. A convenient way to take both into account is by analyzing their absolute difference, i.e. the generalization risk Rgen = |Remp - Rtrue|. However, a low Rgen guarantees similar risks, leaving open the possibility of similar but still high risks. Therefore, we consider maintaining a low empirical risk Remp while minimizing the risk difference Rgen.

The empirical risk Remp can be further minimized since it is "measurable", i.e. its value can be calculated with respect to a particular loss function using the observed data. However, the true risk Rtrue is unknown. Therefore, the risk difference Rgen cannot be directly calculated. Instead, one can try to minimize the value of the upper bound of Rgen. There is a broad range of upper bounds on Rgen, but their values are often not directly controllable. In this work, we focus on the bounds of the expected generalization risk ˆRgen which were recently established and are based on the notion of stability.

Here, we define the generalization risk of a learning algorithm L : ∪∞N=1Z N → H, or simply the expected generalization risk ˆRgen. The algorithm L selects a hypothesis h : ∪∞N=1Z N → Y from a hypothesis space H using a subset of the whole observation space. Given a bounded parametric loss function ℓ(·, h) : Z → [0, 1], the generalization risk of L w.r.t. ℓ(·, h) is defined as the absolute difference between the empirical and true risks of L, i.e.

The empirical and true risks of L are defined as the expected empirical and expected true risks of h, i.e.

where Remp(h,D) = 1N ∑Ni=1 ℓ(zi, h); Rtrue(h) = Ez∼P(z)[ℓ(z, h)]. We next present a theorem that outlines a stability-based bound of ˆRgen which the proposed objective function is based on.

Theorem 1 ([Alabdulmohsin, 2017; Alabdulmohsin, 2015]) For any learning algorithm L : ∪∞N=1Z N → H, algorithmic stability is both a necessary and sufficient condition for uniform generalization. Moreover,

where S(ℓ(·, h), ztrn) denotes the mutual stability between ℓ(·, h) and ztrn which essentially represents the overlap between their probability distributions. For the complete proof of the theorem, just follow the reference to its original source.

Next, we justify our choice of the upper bound in Eq. (4) by providing the main two reasons for it:

• In [Alabdulmohsin, 2015] it has been proven that uniform generalization is essentially equivalent to algorithmic stability. Since algorithmic stability is defined by the mutual stability S(ℓ(·, h), ztrn) (on which the bound Eq. (4) is based),

it follows that mutual stability is also connected to uniform generalization. In other words, one needs to be able to control mutual stability in order to improve generalization performance.

• The value of the chosen bound still cannot be directly minimized, but at least it can be indirectly controlled, unlike most of the bounds on Rgen.

We focus solely on 1 - S(ℓ(·, h), ztrn) as it is the tightest bound presented in [Alabdulmohsin, 2015, Trm. 1].

3.2 Bias-Variance Balancing Objective Function

The ultimate goal is to minimize both ˆRemp and the upper bound of ˆRgen at the same time, which can be achieved by tightening the upper bound of ˆRtrue. Expressing the chosen bound (Eq. (4)) in terms of ˆRtrue yields

According to Eq. (3), the upper bound in Eq. (5) can be rewritten as

To control the above upper bound, we define the following objective function:

It is clear that the first term in Eq. (7), which minimizes the empirical error Remp(h,D), will induce a decrease of ˆRemp(L). The second term dCorr(ℓ(·, h), ztrn), defined later in Eq. (8), measures the statistical dependence between the loss ℓ(·, h) and a random training example ztrn. Minimizing dCorr tightens the upper bound of ˆRgen.

In the following section, we provide details on the connection between the notions of mutual stability and statistical dependence, and further explain how distance correlation can indirectly control the stability term S(ℓ(·, h), ztrn) in Eq. (6).

3.3 How Statistical Dependence Controls Stability?

From an information-theoretic perspective, as the variational information I(ℓ(·, h), ztrn) = 1 - S(ℓ(·, h), ztrn) decreases, the empirical loss ℓ(zi, h) becomes more representative to the true loss ℓ(z, h), z ∼ P(z). In terms of stability theory, this means that the learning algorithm that selected h is becoming more stable, since a high value of S(ℓ(·, h), ztrn) indicates that the probability distribution of ℓ(·, h) is not perturbed by a random training example ztrn. Perfect mutual stability S(ℓ(·, h), ztrn) = 1 is achieved when I(ℓ(·, h), ztrn) = 0, i.e. when ℓ(·, h) and ztrn are statistically independent. Therefore, we account for the dependence between them by introducing the second term in our objective function, which relies on the distance correlation measure between ℓ(·, h) and ztrn. Thinking of ℓ(·, h) and ztrn as random vectors L and Z whose observations are ℓ(z1, h), . . . , ℓ(zN, h) and z1, . . . , zN, respectively, the (empirical) distance correlation between them is defined as the square root of

when dCov2N (L,L)dCov2N (Z,Z) > 0. Otherwise, dCorr2N (L,Z) = 0. In Eq. (8), the squared sample distance covariance is calculated as

Here, di· = 1N ∑Nj=1 dij, d·j = 1N ∑Ni=1 dij, and d·· = 1 N 2 ∑Ni=1 ∑Nj=1 dij, where dij is the euclidean distance between the i-th and j-th observations of the random vector for which it is being calculated.

Distance correlation generalizes the idea of correlation between L and Z, while preserving the following properties:

• (i) dCorr(L,Z) is defined for two random vectors of arbitrary, not necessarily equal dimensions;
• (ii) dCorr(L, Z) = 0 iff L and Z are independent;
• (iii) 0 ≤ dCorr(L,Z) ≤ 1.

From (ii), it is clear that minimizing dCorr(L,Z) in the objective's second term may cause a positive change in the mutual stability S(ℓ(·, h), ztrn). The assumption here is that, in the idealistic case, when the value of dCorr(L,Z) reaches zero, L and Z are independent and 1 - S(ℓ(·, h), ztrn) vanishes in the upper bound (Eq. (6)). A direct relation between distance correlation and information-theoretic mutual stability is not claimed, but rather, in this case it is inspired by (indirectly related to) the concept of mutual stability.

4 Methodology

4.1 Preliminaries: Gaussian CRF

A Continuous Conditional Random Field (CCRF) [Qin et al., 2009] models the conditional distribution of the outputs y, given all inputs X, as P(y|X) = 1Z exp {-α ∑Ni=1 (yi - φ(xi))2 - β ∑i∼j Sij (yi - yj )2} , where [α, β]⊤ is a parameter vector, while Z is a normalization constant, calculated as an integral over y of the term in the exponent. The first term in the exponent models the relevance of an "unstructured predictor" φ by assigning a weight α. In the second term, Sij is extracted by a user-defined similarity measure s(xi, xj) assuming that if the two inputs xi and xj are similar, their corresponding outputs yi and yj are also similar. The similarity matrix S = [Sij]N×N can be seen as an adjacency matrix of an undirected weighted graph whose relevance is determined by the β parameter.

Structured Learning. Since the constituents of the exponent in P(y|X) are defined as quadratic functions in terms of y, the conditional probability can be transposed directly onto a multivariate Gaussian distribution, hence the CCRF becomes a Gaussian CRF (GCRF) [Radosavljevic et al., 2010]. Learning a GCRF boils down to determining the precision matrix Q = αI + βL, where L is the Laplacian of S and I is an identity matrix. Note that Q is used to obtain an explicit expression for the inverse covariance matrix Σ-1 = 2Q. The learning is governed by a convex optimization that strives to determine [ˆα, ˆ β]⊤ = arg maxα,β log (P (y|X; α, β)).

Inference. Since the model relies on a multivariate Gaussian distribution, the estimate of a target vector y′ is the distribution's expected value, for which P(y′ |X′ ) is maximized. Hence, the inference task comes down to calculating µ = Σb = αQ-1 [φ(x′1 ), . . . , φ(x′N ′ )]⊤.

4.2 GeneraLization-Aware Collaborative Ensemble Regressor (GLACER)

The proposed model, GLACER, represents a collaborative ensemble model for structured regression that strives to avoid both underfitting and overfitting by optimizing the proposed objective function (Eq. (7)). Since bagging [Andonova et al., 2002] combined with boosting [Friedman, 2001] can be very effective for achieving both variance and bias reduction [B¨uchlmann and Yu, 2002], GLACER sub-bags multiple GCRFs, each utilizing an LSBoost model as an unstructured predictor. Moreover, it employs collaboration between its constituents by allowing them to interactively exchange examples during training. This example-exchange-guided optimization aims to minimize Eq. (7), thus enhancing the overall model's predictive performance. GLACER's training procedure is outlined in Algorithm 1.

Ensemble Construction (Lines 1-5 in Algorithm 1). In order to employ multiple "local" GCRF models to discover different data substructures, a training set D and its corresponding similarity matrix S are sampled uniformly M times without replacement, thus generating M data subsets D1 , . . . ,DM of size ηN, where η ∈ (0, 1), and corresponding M similarity submatrices S1 , . . . , SM . Then, each Dm is used to train a single GCRF component FDm characterized by αm and βm, ∀m = 1, . . . ,M. Since it is obvious that Sm is passed as an input to the m-th GCRF component together with Dm, Sm is left out from the subscript of FDm for simplicity of notation. The final predictions for the output values are made by combining all GCRF components' outcomes in a subbagging fashion: ΦD(X, S) = 1M ∑Mm=1 FDm(X, S).

Next, we proceed by presenting a general interactive optimization procedure that determines which examples should fall in which subsets, that is, modify Dm and Sm for the components FDm so as to decrease the objective function Robj .

Step I (Lines 7-17 in Algorithm 1). In the first step of the optimization, candidate (potential) swaps between GCRF components are made; and the decrease in the objective function is calculated after each of them occurs. Examples are tentatively assorted for exchange, and selection is determined by evaluating the losses of each GCRF component with respect to all examples in its subset. When evaluating the loss of a regression model, the coefficient of determination R2 is often used. The range of R2 values is (-∞, 1]. Since there is a constraint related to the loss function which should fall between 0 and 1, we define a loss function that measures how well a GCRF component FDm fits an example z as follows,

Thereafter, the worst-fit example zm ∗ of the mth GCRF component is determined by zm ∗ =

argmaxz∈Dm ℓ(FDm, z), ∀m = 1, . . . ,M. In order to improve the predictive performance of a single GCRF component, the worst-fit example from its training subset is exchanged with another worst-fit example from another GCRF component's subset. To avoid duplicating examples and compromising the subset size equality among the components, a pair of GCRF components is allowed to exchange examples only if each of them receives an example from the other one that is not present in its subset, and vice versa. Therefore, after an example exchange occurs, ΦD will still represent a valid subbagging ensemble. As per this constraint, the max-loss examples are exchanged between a pair of GCRF components as (Djp , Dk p ) = (Dj \ {zj∗} ∪ {zk ∗ }, Dk \ {zk ∗ } ∪ {zj∗}) , if zk ∗ /∈ Dj ∧ zj∗ /∈ Dk , where Djp and Dk p denote the modified versions of Dj and Dk , respectively. Otherwise, the original

subsets remain the same, i.e. (Djp ,Dk p ) = (Dj , Dk ).

The exchange of zj∗ and zk ∗ is referred to as one "candidate swap". It should be noted that, along with the temporary update of the j-th and k-th subsets, their corresponding similarity matrices Sj and Sk are also updated accordingly. After performing this candidate swap, the outcome of the overall model for arbitrary X and S is

To quantify the contribution of each swap, the difference between the values of the objective function Robj (Eq. (7)), before and after a swap occurs, is calculated:

The value of ∆jk essentially represents the improvement in the objective function value brought by the collaboration between the j-th and k-th GCRF components. This measure is calculated for all pairs of GCRF components.

Step II (Lines 18-20 in Algorithm 1). Once all tentative candidate swaps are examined, examples are exchanged between those two GCRF components that foster the highest decrease in the objective function.

Convergence (Line 6 in Algorithm 1). The aforedescribed steps constitute one iteration of GLACER's training procedure and all of them are repeated until no additional exchange between any pair of GCRF components can further decrease the objective function.

Monotonicity of Robj. Let us assume that τ is the current optimization iteration. If convergence is not the case at this iteration, then there must exist at least one pair (j, k) for which ∆jk > 0. Henceforth, an exchange between a pair of GCRF components must occur, and the fact that ∆jk > 0 only if the exchange fosters an improvement in the value of Robj guarantees that Robj will monotonically decrease as τ increments by one.

Algorithm 1 GLACER

Input: Training set D, #components M, sub-sampling fraction η, similarity matrix S

Procedure: 1: for m= 1, . . . ,M do 2: (Dm, Sm) ← Sub-sample(D, S, η) 3: φ← Train LSBoost(Dm) 4: FDm ← Train GCRF(Dm, Sm, [φ(x)]⊤x∈Dm) 5: Construct ΦD = 1/M ∑Mm=1 FDm(X, S) 6: while at least one ∆jk > 0 do 7: for m = 1, . . . ,M do 8: zm ∗ = argmaxz∈Dm ℓ(FDm, z) 9: for all unique pairs (j, k) ∈ [1,M]2 do 10: if zk ∗ /∈ Dj ∧ zj∗ /∈ Dk then 11: Calculate the objective Robj (ΦD, D) 12: Swap examples between the Dj and Dk 13: Retrain FDj and FDk 14: Construct a modified ensemble ΦD(j,k) 15: Calculate Robj (ΦD(j,k) , D) 16: ∆jk = Robj (ΦD, D) - Robj (ΦD(j,k) , D) 17: Put zj∗ and zk ∗ in their original subsets 18: Find the optimal pair (j ∗ , k∗ ) = argmax(j,k) ∆jk 19: Permanently modify subsets Dj ∗ and Dk∗ 20: In ΦD: replace FDj by FDj∗ , and FDk by FDk∗ Output: Return ΦD

5 Experiments

5.1 Baselines

The baselines considered in this study are listed as follows:

• MLR: Multiple Linear Regression model.
• NN: Feed-Forward Neural Network with 3 layers.
• SVM: SVM regression model with an RBF kernel.
• SUBBAG: Variation of bagging that considers sampling at random, but without replacement to generate training subsets.
• RF: Random Forest ensemble that performs random sampling of examples to generate subsets, followed by a random feature selection within each subset.
• LSB: Gradient-boosting ensemble for additive expansions based on the least-squares fitting criterion.
• (C/NC)NL: Network Lasso, a structured regression model capable of simultaneous clustering and optimization on graphs. Both its convex (C)NL and non-convex (NC)NL variants were trained with λ = 5, same as in [Hallac et al., 2015].

5.2 Setup

In all conducted experiments, the ensemble-based baselines (SUBBAG, RF, LSB) and GLACER were run with M components, while the same sub-sampling fraction η was used to construct training subsets for both SUBBAG and GLACER. Moreover, each baseline other than NL (which is already a structured approach) was run in both traditional (unstructured) and structured mode. A baseline's structured variant is obtained by passing it to a GCRF as its unstructured predictor. Structured variants are distinguished from unstructured ones by the "S-" at the beginning of their names (e.g. S-MLR

denotes a structured MLR). Mean squared error (MSE) was calculated for all models. The average testing MSEs are reported, along with their corresponding two-sided confidence intervals at 90% confidence level.

5.3 Experiments on Synthetic Data

Data Generation. The synthetic data was generated such that it holds a certain structure ([Pavlovski et al., 2017] used a similar generation process). First, N = 3000 examples xi ∈ Rd (d = 5) were generated such that each attribute is normally distributed according to a standard normal distribution. Thereafter, outputs were created as parameterized polynomials with uniformly distributed parameters. Normally distributed noise was applied to these outputs, yielding ∼y. The structure among the examples was created on the basis of an Erd˝os-R´enyi random graph G, each node corresponding to one example. Accordingly, the between-example similarities were calculated as Sij = e-|∼yi-∼yj | in case an edge between i and j exists in G. Otherwise, Sij was set to zero. Finally, GCRF was utilized in a generative manner in order to infer the final outputs as y = α(αI + βL)-1∼y, where α = 1 and β = 5 were chosen, thus making the structure more significant than the input-output relationship within the data.

Parameter Analysis. The predictive performance of GLACER was analyzed under different sets of parameters. Using the aforedescribed data generation procedure, 10 different training and independent test sets were generated. GLACER was then run on each train/test pair with M = 5, 10, 30 components, while for each value of M the subsampling fraction η varied within {0.3, 0.5, 0.7}. The testing MSEs regarding all these different parameter sets are left out due to lack of space. As expected, GLACER stabilizes and achieves greater performance as M increases. When M = 10, GLACER already shows good generalization performance, and η does not seem to play a crucial role. Therefore, in the following experiments we chose to run GLACER with M = 10 and set η to 0.3 for efficiency.

Generalization Capability. Different fractions of data (10%, 50% and the entire training set) were used to further investigate GLACER's generalization performance. For each training data size, the average MSE over 10 repetitions was calculated for GLACER and all baselines. As shown in Table 1, GLACER yields the lowest MSEs under all three training data sizes. Although all models' MSEs decrease with the increased size of training data, GLACER sustains stable predictions when only 50% training data is available, which clearly indicates the generalization performance of GLACER.

Influence of Distance Correlation. To evaluate the influence of distance correlation in the objective function, GLACER's generalization performance was evaluated on different fractions of the same synthetic data in two cases: (1) when the dCorr term is excluded from Robj ; (2) when dCorr is incorporated in Robj. According to Table 2, GLACER manifests lower average MSEs when dCorr is used in Robj . This is consistent as the training data increases. Besides, without dCorr, the average MSE worsens once the training set fraction increases from 50% to 100% which might be an indication of overfitting. On the contrary, incorporating dCorr in Robj prevents from large increases in MSE. This indicates that dCorr plays an important role in the model's generalization performance.

5.4 Real-World Datasets

Sacramento Real-Estate. A collection of 985 real estate transactions were observed in the Greater Sacramento area, California, made over a period of one week in May 2008. Each example refers to a house sale record that contains information about the number of bedrooms and bathrooms, the house area in square feet, and its location in terms of latitude and longitude. The regression task is to predict houses' prices based on their characteristics. Since some attribute values are missing, we used a dataset version that was pre-processed by [Hallac et al., 2015]. An undirected similarity graph was constructed for each training and test set by coupling each house with its 5 nearest houses and the other way around. A weight e-distij was assigned to each existing edge (i, j) based on the geospatial distance between houses i and j.

Medicare Readmissions 1 . This data consists of 1000 hospital records referring to hospitals that have more than ∼150 readmissions. Each record contains information about the number of discharges, the excess readmission ratio, as well as the estimated and expected readmission rates. Given a hospital record, the goal is to predict the number of readmissions at the hospital. The structure among the hospital records was constructed by calculating distij as the Manhattan distance between the attribute values of records i and j, and assigning edge weights in the same way as for the Sacramento dataset.

1 https://data.medicare.gov/data/hospital-compare

Results and Discussion. The Sacramento Real-Estate dataset was split into a training set of 785 house transactions and a test set of 200 transactions (same as in [Hallac et al., 2015]). As for Medicare Readmissions, half of the data was randomly sampled and used for training, while the other half was used for evaluation. The average testing MSEs obtained on both datasets are summarized in Table 3. From these results, it is evident that GLACER outperforms its alternatives by significant margins (the corresponding p-values are smaller than 0.01 for Sacramento, and 0.021 for Medicare). For instance, GLACER is 38.86% (p-value = 5.6 × 10-9 ) and 40.79% (p-value = 2.3 × 10-8 ) more accurate than the convex and non-convex NL on the Sacramento Real-Estate dataset, respectively. GLACER also manifests a lower average MSE than its building blocks, namely SUBBAG and S-LSB on both datasets. As for all other baselines, it obtains substantial MSE decreases on the Sacramento dataset ranging from 9.64% to 56.4%, and from 49.32% up to 96.41% on the Medicare dataset, thus demonstrating a considerable generalization capability.

In addition, GLACER is more stable compared to alternatives, as it has the tightest confidence interval for its average MSE. A possible explanation behind this would be the nature of GLACER's stability-encouraging objective function. Lastly, there is no overlap between the confidence interval for GLACER's average MSE and any other model's interval. Henceforth, GLACER's improvements are statistically significant.

6 Conclusion

In this work, we bridged theory and calculation of a stabilitybased generalization error bound by leveraging the distance correlation measure and proposed a bias-variance balancing objective function which utilizes the properties of this measure to address the trade-off between underfitting and overfitting. In addition, we introduced GLACER, a model for structured regression which optimizes the proposed objective function through an example-exchange-driven optimization. GLACER was assessed on multiple datasets, on which it manifested stable predictions and significantly outperformed a broad range of traditional and structured regression models.

Acknowledgments

This research was supported in part by the Defense Advanced Research Projects Agency (DARPA) and the Army Research Office (ARO) under Contract no. DARPA-THOR-W911NF16-C-0050, an IQVIA Big Data Analytics contract, NSF BIGDATA grant 14476570, and ONR/ONR Global grant N62909-16-1-2222.