--- title: "Introduction to SOLNP ver 2.0.0" output: rmarkdown::html_vignette: css: custom.css code_folding: hide citation_package: natbib toc: yes bibliography: references.bib bibliography-style: apalike natbiboptions: round vignette: > %\VignetteIndexEntry{Introduction to SOLNP ver 2.0.0} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` ```{r setup, echo=FALSE} library(Rsolnp) ``` ## Introduction The SOLNP algorithm, first proposed by @Ye1989, is a Sequential Quadratic Programming (SQP) approach designed to solve general nonlinear optimization problems with both equality and inequality constraints. In 2010, it was implemented in R as the Rsolnp package, originally to fill the need for a robust and flexible constrained optimizer required by the author’s own GARCH modeling package, [rugarch](https://cran.r-project.org/package=rugarch). Since then, Rsolnp has remained one of the few readily available and user-friendly options for nonlinear constrained optimization in R. Despite the growth of R’s ecosystem in statistics and data science, the development and availability of modern nonlinear constrained optimization solvers in R has lagged behind other languages, with notable exceptions being [nloptr](https://cran.r-project.org/package=nloptr) (which wraps NLopt solvers) and [ipoptr](https://github.com/jyypma/ipoptr) (an R interface to the IPOPT solver)^[the latter is no longer available on CRAN due to licensing issues]. Moreover, in the nonlinear programming (NLP) space, most state-of-the-art solvers remain closed source and commercial, in contrast to other algorithmic advances that have transitioned from academia to open-source availability. Rsolnp does not claim to be competitive with commercial grade nonlinear solvers, but it has generally served its purpose reasonably well for small to medium sized problems. ## Problem Statement The general nonlinear optimization problem addressed by SOLNP can be stated as follows: $$ \begin{align} & \min_{x \in \mathbb{R}^n} \quad f(x) \\\\ & \text{subject to:} \\\\ & \quad l_g \leq g(x) \leq u_g \qquad \text{(general inequality constraints)} \\\\ & \quad h(x) = B \qquad \text{(equality constraints with bound vector \(B\))} \\\\ & \quad x^L \leq x \leq x^U \qquad \text{(variable bounds)} \end{align} $$ where: - $f(x)$: Nonlinear objective function to minimize - $g(x)$: Vector of general nonlinear inequality constraint functions, with lower and upper bounds $l_g$, $u_g$ - $h(x)$: Vector of nonlinear equality constraint functions, required to be equal to vector B - $x^L$, $x^U$: Lower and upper bounds for the variables Note: - Standard equality constraints of the form $h(x) = 0$ can be recovered by setting B = 0. - General inequality constraints allow for two-sided bounds (e.g., $l_g < g(x) < u_g$), not just upper or lower bounds individually. ### Inequalities Each inequality constraint $l_i \leq g_i(x) \leq u_i$ is transformed into two equality constraints via slack variables. Introduce slack variable $s_i \geq 0$, then define: $$ h_i(x, s_i) = g_i(x) - s_i = l_i \quad \text{if } g_i(x) \geq l_i \\ h_i(x, s_i) = g_i(x) + s_i = u_i \quad \text{if } g_i(x) \leq u_i $$ This can be unified into: $$ g_i(x) - s_i = l_i \quad \text{and} \quad g_i(x) + s_i = u_i $$ So every bounded inequality becomes an equality constraint with a slack, and the slack becomes part of the optimization variable vector. In the actual implementation, if a constraint is double-bounded, a "penalty equality constraint" is introduced $h_i(x, s_i) := g_i(x) - s_i = m_i$ and add a penalty term to objective if $s_i < 0$. ### Slack Variable Embedding The augmented optimization vector becomes: $\tilde{x} = \begin{bmatrix} x \\ s \end{bmatrix}$ And the original constraints are reformulated into a system of equalities only: $$ \tilde{g}(x, s) = \begin{bmatrix} g_{\text{eq}}(x) \\ g_{\text{ineq}}(x) - s - l \\ g_{\text{ineq}}(x) + s - u \end{bmatrix} = 0 $$ ### Optimization via Augmented Lagrangian Unlike traditional SQP, SOLNP does not explicitly solve for the Karush–Kuhn–Tucker (KKT) conditions. Instead, it uses a partial augmented Lagrangian: $$ L(x, \lambda, \rho) = f(x) + \lambda^\top g(x) + \frac{\rho}{2} \|g(x)\|^2 $$ where: - $g(x)$ contains only equality constraints (as inequalities are transformed) - $\lambda$ are Lagrange multipliers (estimated iteratively) - $\rho$ is the penalty parameter ### No Stationarity-Based KKT Step Unlike SQP or interior-point methods, SOLNP does not solve the full KKT system at each iteration, instead using: - a quadratic local model - a trust-region style step - penalty-based adjustment (via $\rho$) The stationarity check is only implicit in terms of convergence of $\nabla f + J^\top \lambda \to 0$, not enforced directly. ### Why Higher Tolerances Are Not Preferred Setting the tolerance in SOLNP (or augmented Lagrangian solvers) much tighter than 1e-8—such as 1e-12—can actually degrade solution quality because the algorithm becomes sensitive to numerical round-off and finite precision errors, especially in gradient and Hessian computations or matrix solves. As the tolerance approaches the limits of double-precision arithmetic, the solver may "chase noise" rather than real improvements, leading to erratic or even worse solutions, with unreliable convergence and possible violation of constraints. For most problems, practical and reliable results are achieved with tolerances in the 1e-6 to 1e-8 range. ## C++ Version The original solnp function did not make use of analytic gradients or Jacobians, and was written entirely in R, a direct translation from the Matlab code. Since version 2.0.0, a new function `csolnp`, written in C++^[Using [Rcpp](https://cran.r-project.org/package=Rcpp) and [RcppArmadillo](https://cran.r-project.org/package=RcppArmadillo)] allows the user to provide analytic gradient and/or Jacobians for the inequality and equality constraints. In the absence of analytic functions, finite differences are used using the functions from the [numDeriv](https://cran.r-project.org/package=numDeriv) package. The function signature for the csolnp function is now slightly different from that of solnp in both content and code styling: ```{r} args(csolnp) ``` vs ```{r} args(solnp) ``` Speedup of anywhere up to 10x can be expected for some problems. Additionally, certain enhancements have been made in the C++ code in regards to criteria for early termination to avoid stalling in the line search phase. ## Multi-Start Solver A new function `csolnp_ms` implements a multi-start strategy, similar to the existing `gosolnp` function. However, it differs in the way starting candidate values are calculated, with the original function using rejection sampling to find feasible starting points, whilst the new function uses a combination of buffered box sampling and a fast interior-point feasibility projection. This approach efficiently generates initial values that strictly satisfy parameter bounds and nonlinear constraints by first sampling points slightly away from the boundaries and then rapidly projecting them into the feasible region using a penalized minimization. This results in more reliable and diverse initial candidates for multi-start optimization, especially in problems with complex or tightly constrained feasible regions. ## Test Suite As of version 2.0.0, a new test suite based on @Hock1980 and @Schittkowski2012 has been included, translated from the Fortran codes [here](https://klaus-schittkowski.de/tpnp.htm). Currently, about 60 problems of the 306 have been translated, and it is the intention of the author to eventually translate all the tests. Each test, returns a list with the following information: - fn : the objective function - gr : the analytic gradient function - eq_fn : the equality function - eq_b : the equality bounds - eq_jac : the Jacobian of the equality function - ineq_fn : the inequality function - ineq_lower : the lower bounds for the inequalities - ineq_upper : the upper bounds for the inequalities - ineq_jac : the Jacobian of the inequality function - lower : lower parameter bounds - upper : upper parameter bounds - start : initialization parameters - best_fn : best known optimal objective - best_par : best known optimal parameters - name : test name The suite can be called at once or by reference to a specific test, as illustrated below. ```{r} prob <- solnp_problem_suite(number = 10) sol <- csolnp(pars = prob$start, fn = prob$fn, gr = prob$gr, eq = prob$eq_fn, eq_b = prob$eq_b, eq_jac = prob$eq_jac, ineq_fn = prob$ineq_fn, ineq_lower = prob$ineq_lower, ineq_upper = prob$ineq_upper, ineq_jac = prob$ineq_jac, lower = prob$lower, upper = prob$upper) print(prob$name) print(c("convergence" = sol$convergence)) print(c("csolnp objective" = sol$objective, "best objective" = prob$best_fn)) print(sol$elapsed_time) ``` We also test a more specific problem based on the GARCH(1,1) model with analytic derivatives: ```{r} prob <- solnp_problem_suite(suite = "Other", number = 4) sol <- csolnp(pars = prob$start, fn = prob$fn, gr = prob$gr, eq = prob$eq_fn, eq_b = prob$eq_b, eq_jac = prob$eq_jac, ineq_fn = prob$ineq_fn, ineq_lower = prob$ineq_lower, ineq_upper = prob$ineq_upper, ineq_jac = prob$ineq_jac, lower = prob$lower, upper = prob$upper) print(prob$name) print(c("convergence" = sol$convergence)) print(c("csolnp objective" = sol$objective, "best objective" = prob$best_fn)) print(sol$elapsed_time) ``` Finally, we take a look at a `solver-killer` problem, HS55 (from the Hock-Schittkowski suite), which poses the following challenges for NLP solvers: * Extremely Flat Objective Near Optimum: The objective function in HS55 is nearly constant (very flat) near the solution, so gradients become very small and the solver can struggle to find a meaningful descent direction. This leads to slow progress or stalling. * Tight and Nearly Active Nonlinear Constraints: The feasible region is defined by nonlinear constraints that are nearly active at the optimum, meaning any small move risks violating feasibility. This makes line search and step acceptance delicate. * Ill-conditioning and Sensitivity: The combination of flatness and tight constraints results in an ill-conditioned optimization landscape. Small numerical errors or poorly scaled steps can cause the solver to jump out of the feasible region, or make the algorithm oscillate or stagnate. * Difficult Constraint Jacobians: The constraint gradients in HS55 can be nearly linearly dependent or poorly scaled near the optimum, leading to numerical instability when solving the KKT system or updating multipliers. * Slow Convergence/Failure to Converge: Many solvers take thousands of iterations and may still not converge to high accuracy, or may falsely report convergence with a suboptimal or infeasible solution. Instead of using `csolnp` which fails in this case, we try out the `csolnp_ms` approach. ```{r,warning=FALSE} prob <- solnp_problem_suite(number = 55) sol <- csolnp_ms(prob$fn, gr = prob$gr, eq_fn = prob$eq_fn, eq_b = prob$eq_b, eq_jac = prob$eq_jac, ineq_fn = prob$ineq_fn, ineq_lower = prob$ineq_lower, ineq_upper = prob$ineq_upper, ineq_jac = prob$ineq_jac, lower = prob$lower, upper = prob$upper, n_candidates = 200, penalty = 1e5, control = list(min_iter = 1000, max_iter = 100, tol = 1e-8), seed = 300) print(paste0("Solution : ", round(sol$objective,3), " | Best Objective :", round(prob$best_fn,3))) print(paste0("Equaility Violation : ", round(sol$kkt_diagnostics$eq_violation,3))) ``` The result is is not too bad, though the equality violation is still > 1e-6 which is what we were aiming for. ## Differences with full Active-Set/SQP solvers Below is a short summary of the key differences between the SOLNP algorithm and full Active-set/SQP methods: | Aspect | SOLNP | Active-Set SQP / QP | |-----------------------|------------------------------------------|--------------------------------------| | Inequality handling | Slacks, penalty, all treated equally | Active set, complementarity enforced | | Subproblem | Least squares (QR, no QP) | Quadratic programming | | KKT stationarity | Not always enforced | Enforced at every step | | Complementarity | Not enforced | Enforced at every step | | Scaling sensitivity | Can be significant | Generally better | | Globalization | Penalty parameter tuning crucial | Line search/trust region, more robust| | Numerical behavior | Simple, robust for easy problems | More complex, better for hard ones | ## Conclusion The C++ rewrite of solnp, inclusion of analytic derivative options and the new test suite should hopefully elevate the quality of the Rsolnp package. Future enhancements will include expansion of the test suite and a deeper investigation into some of the underlying approaches used and ways to enhance the solver. ## References