\[ \newcommand{\argmin}{\mathop{\mathrm{argmin}}} \newcommand{\argmax}{\mathop{\mathrm{argmax}}} \newcommand{\minimize}{\mathop{\mathrm{minimize}}} \newcommand{\st}{\mathop{\mathrm{subject\,\,to}}} \]

Problem setup

Consider the problem \[\begin{alignat*}{2} &\minimize_\alpha \quad && \sum_{k=1}^r \sum_{i=1}^n w_i \psi_{\tau_k} \bigg(y_i - \sum_{j=1}^p \alpha_j q_{ijk} \bigg) \\ &\st && \sum_{j=1}^p \alpha_j = 1, \; \alpha_j \geq 0. \end{alignat*}\] Here \(\tau_k\), \(k=1,\ldots,r\) is a set of quantile levels, assumed to be in increasing order and each \(q_{ijk}\) is an estimate of the quantile of \(y_i\) at the level \(\tau_k\), from ensemble component member \(j\) Also, \[ \psi_\tau(v) = \max\{\tau v, (\tau-1) v)\}, \] often called the “pinball” or “tilted \(\ell_1\)” loss, for a quantile level \(\tau \in (0,1)\), and \(w_i\), \(i=1,\ldots,n\) are observation weights. A more flexible approach would be to estimate a separate ensemble weight \(\alpha_{jk}\) per component method \(j\) and quantile level \(k\): \[\begin{alignat*}{2} &\minimize_\alpha \quad && \sum_{k=1}^r \sum_{i=1}^n w_i \psi_{\tau_k} \bigg(y_i - \sum_{j=1}^p \alpha_{jk} q_{ijk} \bigg) \\ &\st && \sum_{j=1}^p \alpha_{jk} = 1, \; \alpha_{jk} \geq 0. \end{alignat*}\] As a form of regularization, we can additionally incorporate noncrossing constraints into the above optimization, which take the form: \[ \alpha_{\bullet,k}^T q \leq \alpha_{\bullet,k+1}^T q, \; q \in \mathcal{Q}. \] where \(\mathcal{Q}\) is some collection of points over which to enforce the constraints (for example, the training points, or the training points along with some unlabeled test points).

LP reformulation

Here are the LP formulations of the two quantile stacking approaches. The standard one: \[\begin{alignat*}{2} &\minimize_{\alpha,u} \quad && \sum_{i=1}^n w_i \sum_{k=1}^r u_{ik} \\ &\st \quad && u_{ik} \geq \tau_k \bigg(y_i - \sum_{j=1}^p \alpha_j q_{ijk}\bigg), \\ &&& u_{ik} \geq (\tau_k-1)\bigg(y_i - \sum_{j=1}^p \alpha_j q_{ijk}\bigg), \\ &&& \sum_{j=1}^p \alpha_j = 1 \; \alpha_j \geq 0. \end{alignat*}\] The flexible one: \[\begin{alignat*}{2} &\minimize_{\alpha,u} \quad && \sum_{i=1}^n w_i \sum_{k=1}^r u_{ik} \\ &\st \quad && u_{ik} \geq \tau_k \bigg(y_i - \sum_{j=1}^p \alpha_{jk} q_{ijk}\bigg), \\ &&& u_{ik} \geq (\tau_k-1)\bigg(y_i - \sum_{j=1}^p \alpha_{jk} q_{ijk}\bigg), \\ &&& \sum_{j=1}^p \alpha_{jk} = 1 \; \alpha_{jk} \geq 0, \\ &&& \alpha_{\bullet,k}^T q \leq \alpha_{\bullet,k+1}^T q, \; q \in \mathcal{Q}. \end{alignat*}\]

Heavy-tailed example

We give a simple example of regression data with a skewed error distribution: Gaussian for the left tail, and t-distributed (with 3 degrees of freedom) for the right tail. We use three methods to estimate the conditional quantile function, at 23 quantile levels:

  1. lasso, tuned by cross-validation (CV), plus Gaussian tails;
  2. quantile lasso at 5 quantile levels, also tuned by CV, then extrapolated out to the full set of 23;
  3. quantile lasso “refit” (starting from the previous CV-tuned model) to the full 23 quantile levels.

As we can see from the plots below, the first method often generally does poorly, both in the tails and also in the middle of the distribution (because it models the conditional mean, not the conditional median, and these are quite different due to skewness); the second method often does better in the left tail but underestimates quantiles in the extreme right tail (because it extrapolates using a Gaussian quantile function); and the third method often does better in the right tail and underestimates quantiles the extreme left tail. We can see that flexible stacking roughly learns how to account for these complementary strengths/weaknesses, and sets ensembles weights accordingly.

library(glmnet)
## Loading required package: Matrix
## Loaded glmnet 3.0-2
library(quantgen)

set.seed(33)
n = 300
p = 50
x = matrix(rnorm(n*p), n, p)
mu = function(x) x[1] + x[2]
e = ifelse(runif(n) < 0.5, -abs(rnorm(n)), abs(rt(n, df=3)))
y = apply(x, 1, mu) + e

# Histogram of error distribution: skewed, and heavy-tailed on one side
hist(e, breaks=40, col="lightblue", main="Error distribution", prob=TRUE)
lines(density(e), lwd=2)

# Run CV for usual lasso, and quantile lasso
tau = c(0.1, 0.3, 0.5, 0.7, 0.9)
glmnet_obj = cv.glmnet(x, y, nfolds=5)
quant_obj = cv_quantile_lasso(x, y, tau=tau, nlambda=30, nfolds=5, lp_solver="gurobi", verbose=TRUE, sort=TRUE)
## CV fold 1 ...
## Problems solved (of 150): 5 ... 10 ... 15 ... 20 ... 25 ... 30 ... 35 ... 40 ... 45 ... 50 ... 55 ... 60 ... 65 ... 70 ... 75 ... 80 ... 85 ... 90 ... 95 ... 100 ... 105 ... 110 ... 115 ... 120 ... 125 ... 130 ... 135 ... 140 ... 145 ... 150 ... 
## CV fold 2 ...
## Problems solved (of 150): 5 ... 10 ... 15 ... 20 ... 25 ... 30 ... 35 ... 40 ... 45 ... 50 ... 55 ... 60 ... 65 ... 70 ... 75 ... 80 ... 85 ... 90 ... 95 ... 100 ... 105 ... 110 ... 115 ... 120 ... 125 ... 130 ... 135 ... 140 ... 145 ... 150 ... 
## CV fold 3 ...
## Problems solved (of 150): 5 ... 10 ... 15 ... 20 ... 25 ... 30 ... 35 ... 40 ... 45 ... 50 ... 55 ... 60 ... 65 ... 70 ... 75 ... 80 ... 85 ... 90 ... 95 ... 100 ... 105 ... 110 ... 115 ... 120 ... 125 ... 130 ... 135 ... 140 ... 145 ... 150 ... 
## CV fold 4 ...
## Problems solved (of 150): 5 ... 10 ... 15 ... 20 ... 25 ... 30 ... 35 ... 40 ... 45 ... 50 ... 55 ... 60 ... 65 ... 70 ... 75 ... 80 ... 85 ... 90 ... 95 ... 100 ... 105 ... 110 ... 115 ... 120 ... 125 ... 130 ... 135 ... 140 ... 145 ... 150 ... 
## CV fold 5 ...
## Problems solved (of 150): 5 ... 10 ... 15 ... 20 ... 25 ... 30 ... 35 ... 40 ... 45 ... 50 ... 55 ... 60 ... 65 ... 70 ... 75 ... 80 ... 85 ... 90 ... 95 ... 100 ... 105 ... 110 ... 115 ... 120 ... 125 ... 130 ... 135 ... 140 ... 145 ... 150 ... 
## Computing CV errors and optimum lambdas ...
## Refitting on full training set with optimum lambdas ...
## Problems solved (of 5): 1 ... 2 ... 3 ... 4 ... 5 ...
plot(glmnet_obj)

plot(quant_obj)

# Refit quantile lasso at more quantile levels
tau_new = c(0.01, 0.025, seq(0.05, 0.95, by=0.05), 0.975, 0.99) 
refit_obj = refit_quantile_lasso(quant_obj, x, y, tau_new, lp_solver="gurobi", verbose=TRUE)
## Problems solved (of 23): 5 ... 10 ... 15 ... 20 ...
# Generate test data 
n0 = 300
x0 = matrix(rnorm(n0*p), n0, p)
e0 = ifelse(runif(n0) < 0.5, -abs(rnorm(n0)), abs(rt(n0, df=3)))
y0 = apply(x0, 1, mu) + e0

# Predicted quantiles at test points 
qtrue = outer(apply(x0, 1, mu), ifelse(tau_new < 0.5, qnorm(tau_new), qt(tau_new, df=3)), "+")
qpred1 = outer(as.numeric(predict(glmnet_obj, x0)), qnorm(tau_new), "+")
qpred2_init = predict(quant_obj, x0, sort=TRUE)
qpred2 = quantile_extrapolate(tau, qpred2_init, tau_new, qfun_left=qnorm, qfun_right=qnorm)
qpred3 = predict(refit_obj, x0, sort=TRUE)

par(mfrow=c(1,3))
for (i in 1:9) {
  plot(tau_new, qtrue[i,], type="o", ylim=range(qtrue, qpred1, qpred2, qpred3), ylab="Quantile")
  lines(tau_new, qpred1[i,], col=2, pch=20, type="o")
  lines(tau_new, qpred2[i,], col=3, pch=20, type="o")
  lines(tau_new, qpred3[i,], col=4, pch=20, type="o")
  legend("topleft", legend=c("True", "Lasso", "QLasso (extrap)", "QLasso (refit)"), 
         col=1:4, pch=c(21,20,20,20))
}

# Construct array of predicted quantiles
qarr = combine_into_array(qpred1, qpred2, qpred3)

# Standard stacking: one weight per ensemble member
st_obj1 = quantile_ensemble(qarr, y0, tau_new, lp_solver="gurobi", verbose=TRUE)
## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 13801 rows, 6903 columns and 55203 nonzeros
## Model fingerprint: 0x171bd348
## Coefficient statistics:
##   Matrix range     [4e-05, 6e+00]
##   Objective range  [1e+00, 1e+00]
##   Bounds range     [0e+00, 0e+00]
##   RHS range        [5e-06, 7e+00]
## 
## Concurrent LP optimizer: dual simplex and barrier
## Showing barrier log only...
## 
## Presolve removed 5130 rows and 2100 columns
## Presolve time: 0.03s
## Presolved: 8671 rows, 4803 columns, 34683 nonzeros
## 
## Ordering time: 0.00s
## 
## Barrier statistics:
##  Dense cols : 3
##  AA' NZ     : 2.988e+04
##  Factor NZ  : 3.856e+04 (roughly 6 MBytes of memory)
##  Factor Ops : 1.736e+05 (less than 1 second per iteration)
##  Threads    : 1
## 
##                   Objective                Residual
## Iter       Primal          Dual         Primal    Dual     Compl     Time
##    0   3.47349691e+04  3.74105079e+02  2.54e+01 0.00e+00  5.16e+01     0s
##    1   2.96935223e+04 -4.43163854e+02  3.72e+00 4.61e-01  5.09e+00     0s
##    2   9.34048009e+03  3.47161975e+02  1.71e-01 2.18e-07  6.77e-01     0s
##    3   3.81137137e+03  1.19711986e+03  3.78e-02 4.95e-09  1.95e-01     0s
##    4   2.77592183e+03  1.51240881e+03  1.08e-02 6.18e-09  9.39e-02     0s
##    5   2.62076636e+03  1.70784716e+03  6.16e-03 1.79e-09  6.78e-02     0s
##    6   2.46225904e+03  2.01920800e+03  2.22e-03 2.22e-16  3.29e-02     0s
##    7   2.40852633e+03  2.18752147e+03  1.08e-03 2.22e-16  1.64e-02     0s
##    8   2.35631879e+03  2.28383124e+03  1.35e-04 2.49e-10  5.38e-03     0s
##    9   2.35125050e+03  2.31274058e+03  2.18e-05 1.93e-09  2.86e-03     0s
##   10   2.34966461e+03  2.32378304e+03  4.14e-07 1.59e-09  1.92e-03     0s
##   11   2.34924837e+03  2.33115862e+03  1.54e-06 1.10e-09  1.34e-03     0s
##   12   2.34911543e+03  2.33385382e+03  2.05e-06 9.50e-10  1.13e-03     0s
##   13   2.34921962e+03  2.33631888e+03  1.77e-06 7.88e-10  9.57e-04     0s
##   14   2.34904406e+03  2.34610743e+03  1.03e-08 1.59e-10  2.18e-04     0s
##   15   2.34882823e+03  2.34730545e+03  2.52e-09 8.98e-11  1.13e-04     0s
##   16   2.34877983e+03  2.34785555e+03  2.40e-10 4.99e-11  6.86e-05     0s
##   17   2.34875560e+03  2.34850791e+03  9.56e-10 1.82e-11  1.84e-05     0s
##   18   2.34874967e+03  2.34854536e+03  3.31e-09 3.33e-12  1.52e-05     0s
##   19   2.34874920e+03  2.34874771e+03  7.81e-09 2.22e-16  1.11e-07     0s
##   20   2.34874839e+03  2.34874838e+03  1.70e-11 4.68e-10  2.09e-10     0s
## 
## Barrier solved model in 20 iterations and 0.13 seconds
## Optimal objective 2.34874839e+03
## 
## Crossover log...
## 
##        0 DPushes remaining with DInf 0.0000000e+00                 0s
## 
##        0 PPushes remaining with PInf 0.0000000e+00                 0s
## 
##   Push phase complete: Pinf 0.0000000e+00, Dinf 1.1368684e-12      0s
## 
## Iteration    Objective       Primal Inf.    Dual Inf.      Time
##       16    2.3487484e+03   0.000000e+00   0.000000e+00      0s
## 
## Solved with barrier
## Solved in 16 iterations and 0.15 seconds
## Optimal objective  2.348748384e+03
coef(st_obj1)
## [1] 0.2380858 0.4798837 0.2820304
# Flexible stacking: one weight per ensemble member, per quantile level
st_obj2 = quantile_ensemble(qarr, y0, tau_new, tau_groups=1:length(tau_new), lp_solver="gurobi", verbose=TRUE)
## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 20423 rows, 6969 columns and 94869 nonzeros
## Model fingerprint: 0x72c45be4
## Coefficient statistics:
##   Matrix range     [4e-05, 6e+00]
##   Objective range  [1e+00, 1e+00]
##   Bounds range     [0e+00, 0e+00]
##   RHS range        [5e-06, 7e+00]
## 
## Concurrent LP optimizer: dual simplex and barrier
## Showing barrier log only...
## 
## Presolve removed 5184 rows and 2100 columns
## Presolve time: 0.04s
## Presolved: 15239 rows, 4869 columns, 74025 nonzeros
## 
## Ordering time: 0.00s
## 
## Barrier statistics:
##  Dense cols : 69
##  AA' NZ     : 6.922e+04
##  Factor NZ  : 8.517e+04 (roughly 9 MBytes of memory)
##  Factor Ops : 5.033e+05 (less than 1 second per iteration)
##  Threads    : 1
## 
##                   Objective                Residual
## Iter       Primal          Dual         Primal    Dual     Compl     Time
##    0   3.40773208e+04  3.74105079e+02  2.49e+01 0.00e+00  1.06e+01     0s
##    1   3.13973628e+04 -3.48603373e+02  1.83e+01 7.84e-01  6.52e+00     0s
##    2   2.17050534e+04 -1.17168100e+03  6.81e+00 3.78e-08  2.35e+00     0s
##    3   6.98115532e+03 -3.16895819e+02  1.62e+00 2.74e-08  6.15e-01     0s
##    4   3.52164632e+03  7.32154371e+02  3.76e-01 3.57e-09  1.81e-01     0s
##    5   2.94449502e+03  1.28846703e+03  1.96e-01 2.71e-09  9.81e-02     0s
##    6   2.55967311e+03  1.77738835e+03  7.20e-02 3.12e-09  4.23e-02     0s
##    7   2.43039517e+03  2.08039046e+03  2.96e-02 3.92e-09  1.80e-02     0s
##    8   2.38780785e+03  2.22895125e+03  1.81e-02 5.53e-09  8.05e-03     0s
##    9   2.35917126e+03  2.26728480e+03  9.90e-03 4.24e-09  4.61e-03     0s
##   10   2.35135546e+03  2.28585743e+03  7.56e-03 3.03e-09  3.27e-03     0s
##   11   2.34077763e+03  2.29520036e+03  4.38e-03 2.33e-09  2.27e-03     0s
##   12   2.33528134e+03  2.31071075e+03  2.81e-03 1.11e-09  1.22e-03     0s
##   13   2.33264558e+03  2.31523504e+03  1.98e-03 8.10e-10  8.65e-04     0s
##   14   2.32850824e+03  2.32016272e+03  8.05e-04 1.03e-09  4.14e-04     0s
##   15   2.32712112e+03  2.32330456e+03  4.44e-04 5.77e-10  1.89e-04     0s
##   16   2.32652150e+03  2.32500868e+03  3.57e-04 3.55e-09  7.48e-05     0s
##   17   2.32590813e+03  2.32545646e+03  8.43e-05 7.16e-10  2.23e-05     0s
##   18   2.32565188e+03  2.32560705e+03  9.70e-06 2.61e-10  2.21e-06     0s
##   19   2.32561891e+03  2.32561460e+03  5.43e-07 8.34e-11  2.14e-07     0s
##   20   2.32561691e+03  2.32561682e+03  2.60e-09 2.71e-10  4.46e-09     0s
##   21   2.32561689e+03  2.32561689e+03  2.60e-11 9.25e-11  9.28e-12     0s
## 
## Barrier solved model in 21 iterations and 0.23 seconds
## Optimal objective 2.32561689e+03
## 
## Crossover log...
## 
##        0 DPushes remaining with DInf 0.0000000e+00                 0s
## 
##        1 PPushes remaining with PInf 0.0000000e+00                 0s
##        0 PPushes remaining with PInf 0.0000000e+00                 0s
## 
##   Push phase complete: Pinf 0.0000000e+00, Dinf 8.7383399e-13      0s
## 
## Iteration    Objective       Primal Inf.    Dual Inf.      Time
##       13    2.3256169e+03   0.000000e+00   0.000000e+00      0s
## 
## Solved with barrier
## Solved in 13 iterations and 0.26 seconds
## Optimal objective  2.325616894e+03
coef(st_obj2)
##      0.01     0.025       0.05 0.1       0.15 0.2      0.25       0.3 0.35
## [1,]    0 0.0000000 0.00000000   0 0.00000000   0 0.0000000 0.0000000    0
## [2,]    1 0.8401616 0.95802079   0 0.98590488   1 0.2989672 0.7017906    1
## [3,]    0 0.1598384 0.04197921   1 0.01409512   0 0.7010328 0.2982094    0
##             0.4      0.45       0.5      0.55       0.6      0.65
## [1,] 0.01489847 0.1647208 0.3397269 0.3881592 0.4335043 0.4744474
## [2,] 0.86078398 0.7480508 0.6602731 0.6118408 0.5664957 0.3119742
## [3,] 0.12431755 0.0872284 0.0000000 0.0000000 0.0000000 0.2135784
##            0.7      0.75        0.8      0.85      0.9      0.95     0.975
## [1,] 0.5221814 0.6143174 0.72695885 0.8359713 0.618521 0.4923194 0.0000000
## [2,] 0.0000000 0.0000000 0.17813977 0.1640287 0.381479 0.1629737 0.4065387
## [3,] 0.4778186 0.3856826 0.09490138 0.0000000 0.000000 0.3447069 0.5934613
##           0.99
## [1,] 0.0000000
## [2,] 0.3841211
## [3,] 0.6158789
# Somewhere in the middle: group the extreme 3 quantiles together on either tail, and the middle 
st_obj3 = quantile_ensemble(qarr, y0, tau_new, tau_groups=c(rep(1,3),rep(2,17),rep(3,3)), lp_solver="gurobi", verbose=TRUE)
## Gurobi Optimizer version 9.0.2 build v9.0.2rc0 (mac64)
## Optimize a model with 20483 rows, 6969 columns and 94989 nonzeros
## Model fingerprint: 0xefb5dc64
## Coefficient statistics:
##   Matrix range     [4e-05, 6e+00]
##   Objective range  [1e+00, 1e+00]
##   Bounds range     [0e+00, 0e+00]
##   RHS range        [5e-06, 7e+00]
## 
## Concurrent LP optimizer: dual simplex and barrier
## Showing barrier log only...
## 
## Presolve removed 11221 rows and 2160 columns
## Presolve time: 0.06s
## Presolved: 9262 rows, 4809 columns, 38223 nonzeros
## 
## Ordering time: 0.00s
## 
## Barrier statistics:
##  Dense cols : 9
##  AA' NZ     : 3.342e+04
##  Factor NZ  : 4.273e+04 (roughly 6 MBytes of memory)
##  Factor Ops : 2.027e+05 (less than 1 second per iteration)
##  Threads    : 1
## 
##                   Objective                Residual
## Iter       Primal          Dual         Primal    Dual     Compl     Time
##    0   3.53208152e+04  3.74105079e+02  2.58e+01 0.00e+00  4.73e+01     0s
##    1   3.29628104e+04 -3.90148222e+02  1.33e+01 1.69e+00  1.31e+01     0s
##    2   2.29939910e+04 -3.87300517e+02  2.08e+00 3.43e-07  1.94e+00     0s
##    3   5.83250300e+03  2.33760917e+02  3.20e-01 1.27e-08  4.28e-01     0s
##    4   2.78148151e+03  1.57345778e+03  4.62e-02 9.33e-09  8.70e-02     0s
##    5   2.49503771e+03  2.02501364e+03  1.95e-02 1.71e-08  3.36e-02     0s
##    6   2.42133882e+03  2.14770822e+03  1.07e-02 1.68e-08  1.95e-02     0s
##    7   2.37780706e+03  2.20905439e+03  4.28e-03 1.10e-08  1.20e-02     0s
##    8   2.35426849e+03  2.29862613e+03  1.34e-03 8.13e-09  3.95e-03     0s
##    9   2.34894736e+03  2.31598300e+03  7.10e-04 4.87e-09  2.34e-03     0s
##   10   2.34573096e+03  2.32124917e+03  3.27e-04 3.87e-09  1.74e-03     0s
##   11   2.34505803e+03  2.32650955e+03  2.77e-04 3.05e-09  1.32e-03     0s
##   12   2.34446317e+03  2.32861731e+03  2.25e-04 2.87e-09  1.13e-03     0s
##   13   2.34323113e+03  2.33225109e+03  8.82e-05 2.17e-09  7.80e-04     0s
##   14   2.34298397e+03  2.33388813e+03  6.40e-05 1.82e-09  6.46e-04     0s
##   15   2.34245686e+03  2.33580199e+03  1.39e-05 1.35e-09  4.73e-04     0s
##   16   2.34242055e+03  2.33641437e+03  1.28e-05 1.21e-09  4.27e-04     0s
##   17   2.34236339e+03  2.33756407e+03  1.05e-05 9.53e-10  3.41e-04     0s
##   18   2.34228120e+03  2.33868364e+03  1.80e-06 7.18e-10  2.56e-04     0s
##   19   2.34218806e+03  2.33913108e+03  1.21e-06 6.19e-10  2.17e-04     0s
##   20   2.34219714e+03  2.34026609e+03  5.31e-07 3.83e-10  1.37e-04     0s
##   21   2.34213539e+03  2.34167612e+03  1.47e-07 1.10e-09  3.26e-05     0s
##   22   2.34212385e+03  2.34184461e+03  9.72e-08 9.54e-10  1.98e-05     0s
##   23   2.34212079e+03  2.34187293e+03  9.17e-08 8.06e-10  1.76e-05     0s
##   24   2.34211249e+03  2.34194745e+03  1.78e-08 5.87e-10  1.17e-05     0s
##   25   2.34210854e+03  2.34197707e+03  1.79e-08 4.91e-10  9.35e-06     0s
##   26   2.34210502e+03  2.34210430e+03  1.14e-09 3.72e-12  5.06e-08     0s
##   27   2.34210488e+03  2.34210488e+03  3.92e-12 3.34e-11  5.06e-11     0s
## 
## Barrier solved model in 27 iterations and 0.20 seconds
## Optimal objective 2.34210488e+03
## 
## Crossover log...
## 
##        0 DPushes remaining with DInf 0.0000000e+00                 0s
## 
##        0 PPushes remaining with PInf 0.0000000e+00                 0s
## 
##   Push phase complete: Pinf 0.0000000e+00, Dinf 2.6956215e-12      0s
## 
## Iteration    Objective       Primal Inf.    Dual Inf.      Time
##       13    2.3421049e+03   0.000000e+00   0.000000e+00      0s
## 
## Solved with barrier
## Solved in 13 iterations and 0.22 seconds
## Optimal objective  2.342104878e+03
coef(st_obj3)
##      0.01 0.025 0.05        0.1       0.15        0.2       0.25
## [1,]    0     0    0 0.28966961 0.28966961 0.28966961 0.28966961
## [2,]    1     1    1 0.66645785 0.66645785 0.66645785 0.66645785
## [3,]    0     0    0 0.04387254 0.04387254 0.04387254 0.04387254
##             0.3       0.35        0.4       0.45        0.5       0.55
## [1,] 0.28966961 0.28966961 0.28966961 0.28966961 0.28966961 0.28966961
## [2,] 0.66645785 0.66645785 0.66645785 0.66645785 0.66645785 0.66645785
## [3,] 0.04387254 0.04387254 0.04387254 0.04387254 0.04387254 0.04387254
##             0.6       0.65        0.7       0.75        0.8       0.85
## [1,] 0.28966961 0.28966961 0.28966961 0.28966961 0.28966961 0.28966961
## [2,] 0.66645785 0.66645785 0.66645785 0.66645785 0.66645785 0.66645785
## [3,] 0.04387254 0.04387254 0.04387254 0.04387254 0.04387254 0.04387254
##             0.9      0.95     0.975      0.99
## [1,] 0.28966961 0.3032585 0.3032585 0.3032585
## [2,] 0.66645785 0.1296539 0.1296539 0.1296539
## [3,] 0.04387254 0.5670876 0.5670876 0.5670876
# Make predictions back on qarr
head(predict(st_obj1, newq=qarr))
##           [,1]      [,2]       [,3]        [,4]       [,5]       [,6]
## [1,] -4.428630 -3.999992 -3.7174945 -3.44806317 -3.2705028 -3.0983393
## [2,] -3.121930 -2.303798 -2.0212531 -1.57652884 -1.3994706 -1.1693358
## [3,] -2.068990 -0.881339 -0.4667785 -0.08046935  0.1652078  0.3637540
## [4,] -3.856904 -3.303942 -2.9742296 -2.62165679 -2.4713571 -2.2548354
## [5,] -2.461607 -1.464681 -1.0103869 -0.62866592 -0.3641431 -0.2138155
## [6,] -3.196885 -2.456294 -2.1064217 -1.67639632 -1.5234953 -1.3107861
##             [,7]       [,8]       [,9]      [,10]      [,11]      [,12]
## [1,] -2.95803801 -2.8362678 -2.7188885 -2.6076810 -2.5054108 -2.4280199
## [2,] -1.00693814 -0.8701292 -0.7205397 -0.5529001 -0.2863353 -0.1457398
## [3,]  0.50846671  0.6523020  0.7653205  0.9414211  1.0849621  1.2254508
## [4,] -2.11173599 -1.9989858 -1.8560463 -1.7146942 -1.5922160 -1.4050354
## [5,] -0.04198922  0.1019476  0.2181188  0.3393690  0.4680913  0.6375706
## [6,] -1.18715309 -1.0743326 -1.0209065 -0.9544769 -0.8762442 -0.8242802
##            [,13]      [,14]      [,15]      [,16]       [,17]      [,18]
## [1,] -2.26436745 -2.1017738 -1.9256445 -1.7952051 -1.46290380 -1.0387754
## [2,] -0.02701688  0.1136685  0.2234250  0.3656670  0.66686101  0.9687887
## [3,]  1.27593633  1.3228259  1.3742339  1.4160784  1.51853594  1.6699907
## [4,] -1.19538940 -1.0271787 -0.8212653 -0.6932740 -0.45626770 -0.1220763
## [5,]  0.83518429  1.0066977  1.1895090  1.3176126  1.54944534  1.6579067
## [6,] -0.71554901 -0.5428933 -0.3473636 -0.2106073  0.08545666  0.5415989
##           [,19]     [,20]    [,21]    [,22]    [,23]
## [1,] -0.3300554 0.2081290 1.109665 1.652999 2.050579
## [2,]  1.4024727 1.8407155 2.402407 2.810060 3.207640
## [3,]  1.8569536 2.0126108 2.450446 3.005996 3.403575
## [4,]  0.1681488 0.5211423 1.341834 1.901240 2.298820
## [5,]  1.7995146 1.9344683 2.447252 2.927486 3.325066
## [6,]  0.9137146 1.2875046 1.901907 2.449387 2.846967
head(predict(st_obj2, newq=qarr))
##           [,1]       [,2]       [,3]        [,4]       [,5]       [,6]
## [1,] -4.667375 -4.2298567 -3.9720529 -3.62385053 -3.4816913 -3.3406635
## [2,] -2.668166 -2.3294556 -1.9987868 -1.62333347 -1.4513397 -1.2780631
## [3,] -1.098134 -0.8009529 -0.4200120 -0.05287728  0.1308695  0.3096993
## [4,] -3.763818 -3.4002633 -3.0808933 -2.72049039 -2.5762609 -2.4301663
## [5,] -1.651198 -1.3737083 -0.9723338 -0.60541560 -0.4214081 -0.2442551
## [6,] -2.755741 -2.4469787 -2.0842062 -1.71118901 -1.5731093 -1.4339217
##              [,7]       [,8]       [,9]      [,10]      [,11]      [,12]
## [1,] -3.144768101 -3.0574765 -2.9644884 -2.8465365 -2.5985087 -2.3180888
## [2,] -1.029718785 -0.9327928 -0.7424100 -0.5662331 -0.3300829 -0.1624060
## [3,]  0.565619403  0.6722759  0.8276161  0.9757706  1.1080312  1.2321958
## [4,] -2.201523436 -2.1398422 -1.9859213 -1.8302542 -1.6033763 -1.3541806
## [5,]  0.009363702  0.1169054  0.2516879  0.3712316  0.4985774  0.6296879
## [6,] -1.226113129 -1.1574061 -1.1158069 -1.0703074 -0.9283627 -0.7521391
##            [,13]      [,14]      [,15]       [,16]      [,17]      [,18]
## [1,] -2.11495232 -1.9179965 -1.7043152 -1.52836276 -1.3367500 -0.8563548
## [2,] -0.03088534  0.1007511  0.2280887  0.37071524  0.5487756  0.7966504
## [3,]  1.27981823  1.3391456  1.4128974  1.50147445  1.6288645  1.8526950
## [4,] -1.16393198 -0.9803320 -0.7609197 -0.61787153 -0.4510492 -0.1602102
## [5,]  0.78616700  0.9386384  1.1224471  1.23754558  1.4799457  1.5042040
## [6,] -0.57085265 -0.3940294 -0.2067650 -0.04225404  0.1767576  0.5518117
##            [,19]       [,20]    [,21]    [,22]    [,23]
## [1,] -0.51112692 -0.05695729 1.052448 2.665426 3.160885
## [2,]  1.01339341  1.48902372 2.211461 3.351414 3.809045
## [3,]  2.07868893  2.18477632 2.614851 3.263672 3.726142
## [4,]  0.05382439  0.38940933 1.355816 2.758713 3.249028
## [5,]  1.64614353  1.89301307 2.460396 3.304929 3.765124
## [6,]  0.78454370  1.14294073 1.861163 3.073762 3.546704
head(predict(st_obj3, newq=qarr))
##           [,1]       [,2]       [,3]        [,4]       [,5]       [,6]
## [1,] -4.667375 -4.3009914 -3.9858810 -3.40997700 -3.2366347 -3.0768410
## [2,] -2.668166 -2.3017818 -1.9866714 -1.56638813 -1.3747533 -1.1875709
## [3,] -1.098134 -0.7317501 -0.4166397 -0.08644747  0.1209827  0.3083722
## [4,] -3.763818 -3.3974341 -3.0823237 -2.60024346 -2.4290300 -2.2602240
## [5,] -1.651198 -1.2848144 -0.9697041 -0.63370334 -0.4237606 -0.2442910
## [6,] -2.755741 -2.3893574 -2.0742470 -1.66885812 -1.5011035 -1.3367541
##            [,7]        [,8]       [,9]      [,10]      [,11]      [,12]
## [1,] -2.9289967 -2.78834064 -2.6799139 -2.5742427 -2.4710450 -2.3722289
## [2,] -1.0179186 -0.85655242 -0.6854762 -0.5133875 -0.3269932 -0.1607061
## [3,]  0.4803943  0.64797449  0.7926233  0.9452897  1.0918080  1.2373401
## [4,] -2.1098326 -1.96846773 -1.8200044 -1.6735834 -1.5311816 -1.3792261
## [5,] -0.0684701  0.09870689  0.2316799  0.3636476  0.4956945  0.6335700
## [6,] -1.1932545 -1.05574223 -0.9881213 -0.9202728 -0.8516715 -0.7876679
##            [,13]      [,14]      [,15]      [,16]      [,17]      [,18]
## [1,] -2.21638656 -2.0601975 -1.9008197 -1.7467540 -1.3030142 -0.8406843
## [2,] -0.03167778  0.1012785  0.2305066  0.3665836  0.7027778  1.0433915
## [3,]  1.28429576  1.3312035  1.3798973  1.4288988  1.5569062  1.6968411
## [4,] -1.20146670 -1.0296415 -0.8508681 -0.6824209 -0.4061144 -0.1103844
## [5,]  0.80386613  0.9706135  1.1402016  1.3030747  1.4600256  1.6020904
## [6,] -0.64408623 -0.4900490 -0.3313703 -0.1800390  0.1598617  0.5289696
##           [,19]     [,20]    [,21]    [,22]    [,23]
## [1,] -0.3270905 0.1719314 1.608686 2.382694 2.811804
## [2,]  1.4114942 1.7930363 2.542664 3.043854 3.472964
## [3,]  1.8492926 2.0098556 2.598795 3.397364 3.826474
## [4,]  0.1854992 0.5029868 1.782054 2.588378 3.017488
## [5,]  1.7563045 1.9227762 2.601360 3.248490 3.677600
## [6,]  0.8919994 1.2678542 2.131018 2.913362 3.342472