Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
34,147 characters · 7 sections · 47 citation commands
Implementing Convex Optimization in R: Two Econometric Examples
Key words: big data, convex optimization, high-dimensional model, numerical solver
JEL code: C13, C55, C61, C87
Zhan Gao: [email removed]. Address: Department of Economics, University of Southern California, 3620 South Vermont Ave. Kaprielian Hall, 300 Los Angeles, CA 90089-0253, USA. Zhentao Shi (corresponding author): [email removed]. Address: Department of Economics, 912 Esther Lee Building, the Chinese University of Hong Kong, Sha Tin, New Territories, Hong Kong SAR, China. Tel: (852) 3943-1432. Fax (852) 2603-5805. We thank Roger Koenker for inspiration and hospitality during the second author's visit to University of Illinois.
Equipped with tremendous growth of computing power over the last few decades, econometricians endeavor to tackle high-dimensional real world problems that we could hardly have imagined before. Along with the development of modern asymptotic theory, computation has gradually ascended onto the central stage. Today, discussion of numerical algorithms is essential for new econometric procedures.
Optimization is at the heart of estimation, and convex optimization is the best understood category. Convex problems are ubiquitous in econometric textbooks. The least square problem is convex, and the classical normal regression is also convex after straightforward reparametrization. Given a linear single-index form, the Logit or Probit binary regression, the Poisson regression and the regressions with a censored or truncated normal distributions are all convex. Another prominent example is the quantile regression koenker1978regression, motivated from its robustness to non-Gaussian errors and outlier contamination.
With the advent of big data, practitioners attempt to build general models that involve hundreds or even more parameters in the hope to capture complex heterogeneity in empirical economic studies. Convex optimization techniques lay out the foundation of estimating these high-dimensional models. Recent years witnesses bajari2015machine, gu2017empirical and doudchenko2016balancing, to name a few, exploring new territories by taking advantage of convexity.
To facilitate practical implementation, cvxinr summarize the operation in R by MOSEK via Rmosek to solve linear programming, conic quadratic programming, quadratic programming, etc. R is open-source software, MOSEK is a proprietary convex optimization solver but offers free academic license, and Rmosek is the \texttt{R} interface that communicates with \texttt{MOSEK}. \texttt{MOSEK} specializes in convex problems with reliable performance, and is competitive in high-dimensional problems.
This paper complements cvxinr's work. We replicate by Rmosek two examples of high-dimensional estimators, namely su2016identifying's classifier-Lasso (C-Lasso) and REL's relaxed empirical likelihood (REL). In addition, we replicate by Rmosek an application of C-Lasso that re-examines China's GDP growth rate chinagdp. These exercises highlight two points. Firstly, the R environment is robust in numerical accuracy for high-dimensional convex optimization and Rmosek takes the lead in computational speed. Second, we showcase the ease of creating new econometric estimators\textemdash often no more than a few lines of code\textemdash by the code snippets in the Appendix. Such convenience lowers the cost of turning an idea into a prototype, and enables researchers to glean valuable insights about their archetypes by experimenting new possibilities. All code in this note is hosted at \url{https://github.com/zhentaoshi/convex_prog_in_econometrics}.
It is common practice to assume in linear fixed-effect panel data models that the cross-sectional units are heterogeneous in terms of the time-invariant individual intercept, while they all share the same slope coefficient. This pooling assumption can be tested and is often rejected in real-world applications. In recent years panel data group structure has been attracting attention. bonhomme2015grouped allow group structure in the intercept and use the $k$-means algorithm for classification. When the slope coefficients exhibit group structure, su2016identifying propose Classifier-Lasso (C-Lasso) to identify the latent group pattern.
We illustrate the penalized least square (PLS), a simple special case of C-Lasso.\footnote{The profile log-likelihood function $Q_{1,nT}\left(\beta\right)=\sum_{i=1}^{n}\sum_{t=1}^{T}\psi\left(w_{it},\beta_{i},\hat{\mu}_{i}\left(\beta_{i}\right)\right)$ for nonlinear models can be reformulated into a separable form, while penalized GMM (PGMM) can be handled under the same optimization framework as PLS. They are discussed in Appendix (ref) and Appendix (ref), respectively.} Given a tuning parameter $\lambda$ and the number of groups $K$, PLS is defined as the solution to \[ \min_{\boldsymbol{\beta},\left(\alpha_{k}\right)_{k=1}^{K}}\,\frac{1}{nT}\sum_{i=1}^{n}\sum_{t=1}^{T}\left(y_{it}-x_{it}^{\prime}\beta_{i}\right)^{2}+\frac{\lambda}{n}\sum_{i=1}^{n}\prod_{k=1}^{K}\left\Vert \beta_{i}-\alpha_{k}\right\Vert _{2} \] where $\boldsymbol{\beta}=\left(\beta_{i}\right)_{i=1}^{n}$. The additive-multiplicative penalty pushes the individual slope coefficients $\beta_{i}$ in the same group toward a common coefficient $\alpha_{k}$. This is not a convex problem, but the optimization with the additive-multiplicative penalty can be approximated by an iterative algorithm, as is explained in the Supplement of su2016identifying. Procedures based on such an iteration have been successfully applied to su2017interactive, su2017numberofgroup and su2017sieve. The iterative algorithm initiates at the within-group estimator, which is consistent when $T$ is large. In the $k$-th sub-step of the $r$-th iteration, $(\boldsymbol{\beta},\alpha_{\tilde{k}})$ is chosen to minimize
where $\gamma_{i}=\prod_{k=1}^{\tilde{k}-1}\Vert\hat{\beta}_{i}^{(r,k)}-\hat{\alpha}_{k}^{(r)}\Vert_{2}\cdot\prod_{k=\tilde{k}+1}^{K}\Vert\hat{\beta}_{i}^{(r-1,k)}-\hat{\alpha}_{k}^{(r-1)}\Vert_{2}.$ The iteration proceeds until the $K$-convex problem numerically converges.
Given the multiplier $\gamma_{i}$, the above optimization problem is convex in $(\boldsymbol{\beta},\alpha_{\tilde{k}})$ and the structure is very close to Lasso. Though the R packages lars and glmnet packages can carry out the standard Lasso, however, it is not straightforward how to modify these functions to accommodate ((ref)), where $\alpha_{\tilde{k}}$ is also an unknown parameter to be optimized. A quick review of cvxinr approach to Lasso will be helpful.
Applying the techniques in the Lasso formulation, we can transform the $l_{2}$-norm terms in ((ref)) and formulate the problem into a conic programming:
where $\theta=\left\{ \beta_{i},\nu_{i},\mu_{i},s_{i},r_{i},t_{i},w_{i}\right\} _{i=1}^{n}$. The convexity is manifest when we write the problem in matrix form, as is displayed in Appendix (ref).
We replicate\textcolor{black}{ }the simulation studies in su2016identifying in R via Rmosek and compare the performance of different numerical optimization approaches. su2016identifying conduct their numerical work in MATLAB via CVX cvx. CVX is a MATLAB add-on package for \emph{disciplined convex optimization} grant2006disciplined. It provides an interface to communicate with commercial or open-source solvers. In the \texttt{R }\textcolor{black}{environment, the de facto solver is }\texttt{\textcolor{black}{optimx}}\textcolor{black}{ optimx; another option is the interface }\texttt{\textcolor{black}{nloptr}}\textcolor{black}{ nloptr that hooks optimization solver }\texttt{\textcolor{black}{NLopt}}\textcolor{black}{ nlopt. They are general-purpose optimization solvers not tailored for convexity. Most recently, CVXR are actively developing }\texttt{\textcolor{black}{CVXR}}\textcolor{black}{, }\texttt{\textcolor{black}{CVX}}\textcolor{black}{'s counterpart in }\texttt{\textcolor{black}{R}}\textcolor{black}{. At this stage, it is integrated with the open-source solver }\texttt{\textcolor{black}{ECOS}}\textcolor{black}{ ECOS. We also consider the counterpart of }\texttt{\textcolor{black}{CVX}}\textcolor{black}{ in }\texttt{\textcolor{black}{Python}}\textcolor{black}{ environment, }\texttt{\textcolor{black}{CVXPY}}\textcolor{black}{ cvxpy, to verify the stability of the algorithm across platforms.}\footnote{\textcolor{black}{In the latest version (Version 0.99), }\texttt{\textcolor{black}{CVXR}}\textcolor{black}{ supports }\texttt{\textcolor{black}{MOSEK}}\textcolor{black}{ by sending the problem to }\texttt{\textcolor{black}{MOSEK}}\textcolor{black}{ in the }\texttt{\textcolor{black}{Python}}\textcolor{black}{ environment. In our experiment, large-scale problems like the C-Lasso cause errors in the communication between }\texttt{\textcolor{black}{R}}\textcolor{black}{ and }\texttt{\textcolor{black}{Python}}\textcolor{black}{. In addition, }\texttt{\textcolor{black}{CVXR}}\textcolor{black}{ with }\texttt{\textcolor{black}{MOSEK}}\textcolor{black}{ currently cannot incorporate problems with nonlinear objective functions and hence cannot be used for REL in Section (ref). }}
We follow DGP 1 in su2016identifying as a benchmark. Table (ref) reports under various combinations of the cross sectional units $n$ and the time length $T$, the root-mean-square error (RMSE) of $\widehat{\alpha}_{1}$ and the probability of correct group classification (correct ratio). The DGP, simulation settings and the indicators are relegated to Appendix (ref) to save space.
Within the R environment, the numerical results of estimation error and classification correct ratio by Rmosek are almost identical to CVXR up to rounding errors.\footnote{The de facto solver optimx breaks down when solving such high dimensional problems. nloptr takes more than a few hours to finish one estimation, which makes the full-scale simulation exercise computational infeasible. In addition, nloptr fails to attain an accurate solution in most cases according to our experiments.} We also implement the simulation in \texttt{MATLAB} via \texttt{CVX} and in \texttt{Python} via \texttt{CVXPY}, the results are largely similar, which demonstrates the robustness of the numerical performance of C-Lasso across different computing platforms.
Practitioners may need to try out different specifications for robustness check in real applications. Without fast optimization solvers, computational cost can become a bottleneck. On the same computing platform of Intel(R) Core(TM) i7-8750H CPU @ 2.20GHz, each case is executed in a single thread and we record the running time in the lower panel in Table (ref). Rmosek significantly outperforms all alternatives. CVX in MATLAB is about $3$ to $5.2$ times slower and CVXPY is about $1.4$ to $3.1$ times slower than Rmosek. Although CVX and \texttt{CVXPY} are also powered by \texttt{MOSEK}, the DCP system takes time to check the convexity of the input problem and automate the formulation. For similar reasons, \texttt{CVXR} is $2.7$ to $4$ times slower than \texttt{Rmosek}. According to CVXR, we can skip the DCP formulation steps with \texttt{CVXR} and the advantage of \texttt{Rmosek} becomes around $1.5$ to $1.8$, which illustrates the advantages of \texttt{MOSEK} over the open-source solver \texttt{ECOS}. In summary, DCP is useful when we are uncertain about the convexity and solvability of a problem. However, for problems that are mathematically verified to be convex, directly calling \texttt{MOSEK} saves much computational time.
While China is now the second largest economy in the world in terms of aggregate GDP, the accuracy of its reported national income accounting has been a topic of constant debate over the years. Most recently, chinagdp utilize local economic indicators that are directly associated with economic activities to estimate China's local and aggregate GDP in order to assess the quality of these numbers. Different regions of this continent-size country are growing at varying pace, thereby resulting in tremendous heterogeneity among its provinces. To control the hidden heterogeneity, chinagdp specify a linear fixed effect model with latent group structure \[ y_{it}=x_{it}^{\prime}\beta_{i}+v_{i}+\varepsilon_{it}, \] where $y_{it}$ is the logarithm of GDP for province $i$ at year $t$, $x_{it}$ includes local economic indicators of interests, $v_{i}$ characterizes the fixed effect of province $i$, and $\varepsilon_{it}$ is the idiosyncratic error. The heterogeneous slope coefficients $\beta_{i}$ captures latent group structures across regions to be determined by C-Lasso. The data span from year 2000 to 2007, i.e. $T=8$. Five indicators are employed as regressors to control observable heterogeneity, namely satellite night lights, national tax revenue, exports, imports, and electricity consumption. These indicators are less susceptible to local officials' manipulation and thus more robustly reflect economic activities for real businesses. Two alternative specifications, one with no satellite night lights and the other with neither \emph{satellite night lights} nor \emph{national tax revenue}, are also considered.
In our implementation, the number of groups $K$ and tuning parameter $\lambda$ are determined by the information criterion proposed in su2016identifying. Trials with different specifications and tuning parameters can be time-consuming, particularly for our panel data of a short $T$ as the algorithm has to iterate many times until numerical convergence. chinagdp estimate the model in MATLAB via CVX and the classification results are reported in chinagdp. We replicate the classification results by Rmosek and compare the accuracy and speed to CVX.
The same as in the original paper, in all three specifications the information criterion determines two groups. Displayed in Table (ref), the classification results by Rmosek are identical to those by CVX reported in chinagdp in the second and the third specifications. However, in the first specification where all five indicators are included, we observe deviance across computing platforms. Liaoning, Zhejiang, Shandong and Yunnan are moved into group 2 according to Rmosek whereas they are left in the group 1 in CVX results. Since this is a short $T$ panel data with $T=8$, the numerical stability is more fragile when we include more regressors.
A key observation in chinagdp is that group 1 gathers Beijing, Shanghai and Hainan, the three provinces with the highest GDP shares of the tertiary sector associated with the provision of services. The Rmosek result retains this feature. As shown in Figure (ref), the tertiary industry GDP shares of the four provinces excluded from the Beijing-Shanghai-Hainan group doe not reach the national aggregate share, except for Liaoning in a single year 2004. Given their relatively lower shares, it is sensible that Rmosek removes them out of the high-share group.
In the parameter tuning process, we compute the information criterion for $K=1,2,3,4$ and 10 candidate $\lambda$ values.\footnote{We generate $\lambda=c\mathrm{var}\left(y\right)T^{-\frac{1}{3}}$ where the constant $c$ varies from $0.001$ to $0.01$.} We report the CPU time consumed by the parameter tuning process in Table (ref). Rmosek is about $8$ to $9.6$ times faster than CVX. It demonstrates the speed gain of Rmosek in real applications.
Besides the regression setting in Section (ref), convex programming is also useful in structural econometric estimation. Consider the models with a “true” parameter $\beta_{0}$ satisfying the unconditional moment condition $\mathbb{E}\left[g\left(Z_{i},\beta_{0}\right)\right]=\boldsymbol{0}_{m}$, where $\left\{ Z_{i}\right\} _{i=1}^{n}$ is the observed data, $\beta\in\mathcal{B\subset}\mathbb{R}^{D}$ is a finite dimensional vector in the parameter space $\mathcal{B}$, and $g$ is an $\mathbb{R}^{m}$-valued moment function. GMM hansen1982GMM and empirical likelihood (EL) owen1988empirical,qin1994empirical are two workhorses dealing with moment restriction models. In particular, EL solves \[ \max_{\beta\in\mathcal{B},\pi\in\Delta_{n}}\,\sum_{i=1}^{n}\log\pi_{i}\quad\text{s.t.}\quad\sum_{i=1}^{n}\pi_{i}g\left(Z_{i},\beta\right)=\boldsymbol{0}_{m} \] where $\Delta_{n}=\left\{ \pi\in\left[0,1\right]^{n}:\sum_{i=1}^{n}\pi_{i}=1\right\} $ is the $n$-dimensional probability simplex. However, neither GMM nor EL can be used to estimate a model with more moment equalities than observations, i.e. $m>n$. To make the optimization feasible, REL relaxes the equality restriction $\sum_{i=1}^{n}\pi_{i}g_{i}\left(\beta\right)=\boldsymbol{0}_{m}$ in EL. REL is defined as the solution to \[ \max_{\beta\in\mathcal{B}}\max_{\pi\in\Delta_{n}^{\lambda}\left(\beta\right)}\,\sum_{i=1}^{n}\log\pi_{i} \] where \[ \Delta_{n}^{\lambda}\left(\beta\right)=\left\{ \pi\in\Delta_{n}:\big|\sum_{i=1}^{n}\pi_{i}h_{ij}\left(\beta\right)\big|\leq\lambda,\:j=1,2,\cdots,m\right\} \] is a relaxed simplex, $\lambda\geq0$ is a tuning parameter, $h_{ij}\left(\beta\right)=g_{j}\left(Z_{i},\beta\right)/\hat{\text{\ensuremath{\sigma}}}_{j}\left(\beta\right)$, $g_{j}\left(Z_{i},\beta\right)$ is the $j$-th component of $g\left(Z_{i},\beta\right)$, and $\hat{\text{\ensuremath{\sigma}}}_{j}\left(\beta\right)$ is the sample standard deviation of $\left\{ g_{j}\left(Z_{i},\beta\right)\right\} _{i=1}^{n}$. The formulation of REL is inspired by Dantzig selector dantzig.
Dantzig selector slacks the sup-norm of the first-order condition for optimality. REL borrows the idea to estimate a finite-dimensional parameter in a structural economic model defined by many moment equalities. Comparing to Dantzig selector, REL uses a nonlinear objective function. It is still convex (in minus likelihood) but quantreg that deals with linear programming problems is no longer applicable.
Similar to standard EL, REL's optimization involves an inner loop and an outer loop. The outer loop for $\beta$ is a general low-dimensional nonlinear optimization, which can be solved by Newton-type methods. With the linear constraints and the logarithm objective, the inner loop is convex in $\pi=\left(\pi_{i}\right)_{i=1}^{n}$. For each $\beta$, the inner problem can be formulated as a separable convex optimization problem in the matrix form
and it is readily solvable in Rmosek by translating the mathematical expression into computer code.
We follow the simulation design in REL, which is described in Appendix (ref). Table (ref) reports the bias and RMSE of the estimation of $\hat{\beta}_{1}$, implemented purely in R with the inner loop by Rmosek and the outer loop by nloptr. The results are close to those in REL, where the code is written in MATLAB with the outer loop handled by the function fmincon and the inner loop by CVX solved by \texttt{MOSEK}.
We also experiment with other numerical alternatives. Since the scale of the optimization problems here is much smaller than C-Lasso, the inner loop can be correctly solved by Rmosek, CVXR, CVX in MATLAB, or even nloptr. These four methods produce virtually identical inner loop results up to rounding errors. This finding confirms the robustness of the R environment in high-dimensional optimization. The difference in Table (ref), therefore, is attributed to the outer loop between the function \texttt{nloptr} in \texttt{R} and the function \texttt{fmincon} in \texttt{MATLAB}.
To evaluate the computational cost, we record the time spent in the inner loop. With $100$ sets of data generated by the same DGP for each sample size, we fix $\beta=\left(0.9,\,0.9\right)$ and only numerically solve the inner loop. Since four approaches have virtually identical inner loop results, we only report the running time of each method in Table (ref). Although CVXR and nloptr are able to correctly solve the problem thanks to its small scale, Rmosek remains $4$ to $30$ times faster than these alternatives. We conjecture that bigger speed gain would be observed in a problem of larger scale.
In this note, we demonstrate numerical implementation via Rmosek of two examples of high-dimensional econometric estimators. The convenience and reliability of high-dimensional convex optimization in R will open new possibilities to create estimation procedures. In the era of big data, we are looking forward to witnessing more algorithms blossoming and flourishing along with theoretical research of high-dimensional models.
\singlespacing \onehalfspacing