---
title: "Quick Tour of Package **optimCheck**"
author: "Martin Lysy"
date: "`r Sys.Date()`"
output:
  html_vignette:
    toc: true
bibliography: references.bib
csl: taylor-and-francis-harvard-x.csl
link-citations: true
vignette: >
  %\VignetteIndexEntry{optimCheck: A Quick Tour}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, eval = FALSE, echo = FALSE}

rmarkdown::render("optimCheck-quicktut.Rmd") # R code to render vignette

```

<!-- latex macros -->
\newcommand{\xx}{\boldsymbol{x}}
\newcommand{\hxx}{\hat{\boldsymbol{x}}}
\newcommand{\txx}{\tilde{\boldsymbol{x}}}
\newcommand{\bb}{\boldsymbol{b}}
\renewcommand{\AA}{\boldsymbol{A}}

## Introduction

The **optimCheck** package provides a set of tools to check that output of an optimization algorithm is indeed at a local mode of the objective function.  The tools include both visual and numerical checks, the latter serving to automate formalized unit tests with e.g., the **R** packages [**testthat**](https://CRAN.R-project.org/package=testthat) or [**RUnit**](https://CRAN.R-project.org/package=RUnit).  

## Example: Quadratic Objective Function

A brief overview of the package functionality is illustrated with the following example.  Let
$$
Q(\xx) = \xx'\AA\xx - 2 \bb'\xx
$$
denote a quadratic objective function in $\xx \in \mathbb R^d$.  If $\AA_{d \times d}$ is a positive-definite matrix, then the unique minimum of $Q(\xx)$ is $\hxx = \AA^{-1}\bb$.  Let us now ignore this information and try to minimize $Q(\xx)$ using **R**'s simplest built-in mode-finding routine, provided by the **R** function `stats::optim()`.

In its simplest configuration, `stats::optim()` requires only the objective function and a starting value $\xx_0$ to initialize the mode-finding procedure.  Let's consider a difficult setting for `stats::optim()`, with a relatively large $d = 12$ and a starting value $\xx_0$ which is far from the optimal value $\hxx$.
```{r, echo = -1}
set.seed(2608) # set seed to fix output
d <- 12 # dimension of optimization problem

# create the objective function: Q(x) = x'Ax - 2b'x
A <- crossprod(matrix(rnorm(d^2), d, d)) # positive definite matrix
b <- rnorm(d)
objfun <- function(x) crossprod(x, A %*% x)[1] - 2 * crossprod(b, x)[1]

xhat <- solve(A, b) # analytic solution

# numerical mode-finding using optim
xfit <- optim(fn = objfun,                    # objective function
              par = xhat * 5,                 # initial value is far from the solution
              control = list(maxit = 1e5))    # very large max. number of iterations

```

### Visual Checks with `optim_proj()`

Like most solvers, `stats::optim()` utilizes various criteria to determine whether its algorithm has converged, which can be assess with the following command:
```{r}
# any value other than 0 means optim failed to converge
xfit$convergence 
```
Here `stats::optim()` reports that its algorithm has converged.  Now let's check this visually with **optimCheck** using *projection plots*.  That is, let $\txx$ denote the potential optimum of $Q(\xx)$.  Then for each $i = 1,\ldots,d$, we plot 
$$
Q_i(x_i) = Q(x_i, \txx_{-i}), \qquad \txx_{-i} = (\tilde x_1, \ldots, \tilde x_{i-1}, \tilde x_{i+1}, \ldots, \tilde x_d).
$$
In other words, projection plot $i$ varies only $x_i$, while holding all other elements of $\xx$ fixed at the value of the potential solution $\txx$.  These plots are produced with the **optimCheck** function `optim_proj()`:
```{r, fig.width = 10, fig.height = 6, out.width = "97%"}
require(optimCheck) # load package

# projection plots
xnames <- parse(text = paste0("x[", 1:d, "]")) # variable names
oproj <- optim_proj(fun = objfun,              # objective function
                    xsol = xfit$par,           # potential solution
                    maximize = FALSE,          # indicates that a local minimum is sought
                    xrng = .5,                 # range of projection plot: x_i +/- .5*|x_i|
                    xnames = xnames)

```
In each of the projection plots, the potential solution $\tilde x_i$ is plotted in red.  The `xrng` argument to `optim_proj()` specifies the plotting range.  Among various ways of doing this, perhaps the simplest is a single scalar value indicating that each plot should span $x_i \pm$ `xrng` $\cdot |x_i|$.  Thus we can see from these plots that `stats::optim()` was sometimes up to 10% away from the local mode of the projection plots.

### Quantification of Projection Plots

Projection plots are a powerful method of assessing the convergence of mode-finding routines to a local mode.  While great for interactive testing, plots are not well-suited to automated unit testing as part of an **R** package development process.  To this end, **optimCheck** provides a means of quantifying the result of a call to `optim_proj()`.  Indeed, a call to `optim_proj()` returns an object of class `optproj` with the following elements:
```{r}
sapply(oproj, function(x) dim(as.matrix(x)))
```
As described in the function documentation, `xproj` and `yproj` are matrices of which each column contains the $x$-axis and $y$-axis coordinates of the points contained in each projection plot.  The `summary()` method for `optproj` objects coverts these to absolute and relative errors in both the potential solution and the objective function.  The `print()` method conveniently displays these results:
```{r}
oproj # same print method as summary(oproj)
```
The documentation for `summary.optproj()` describes the various calculations it provides.  Perhaps the most useful of these are the elementwise absolute and relative differences between the potential solution $\tilde{\xx}$ and $\hxx_\mathrm{proj}$, the vector of optimal 1D solutions in each projection plot. For convenience, these can be extracted with the `diff()` method:
```{r}
diff(oproj) # equivalent to summary(oproj)$xdiff

# here's exactly what these are
xsol <- summary(oproj)$xsol # candidate solution
xopt <- summary(oproj)$xopt # optimal solution in each projection plot
xdiff <- cbind(abs = xopt-xsol, rel = (xopt-xsol)/abs(xsol))
range(xdiff - diff(oproj))
```
Thus it is proposed that a combination of `summary()` and `diff()` methods for projection plots can be useful for constructing automated unit tests.  In this case, plotting itself can be disabled by passing `optim_proj()` the argument `plot = FALSE`.  See the `optimCheck/tests` folder for **testthat** examples featuring: 

- Logistic Regression (`stats::glm()` function).
- Quantile Regression (`quantreg::rq()` function in [**quantreg**](https://CRAN.R-project.org/package=quantreg))
- *Multivariate normal mixtures* (`mclust::emEEE()` in [**mclust**](https://CRAN.R-project.org/package=mclust)).

You can run these tests with the command

```{r eval = FALSE}
testthat::test_package("optimCheck", reporter = "progress")
```

<!-- ```{r, echo = -1, cache = FALSE} -->
<!-- set.seed(151) -->
<!-- # run tests and output timings -->
<!-- tresults <-  testthat::test_package("optimCheck", reporter = "list") -->
<!-- invisible(sapply(tresults, function(tres) { -->
<!--   message("Context: ", tres$context, " [", length(tres$results), " tests -- ", -->
<!--           round(tres$real, 1), " s]") -->
<!-- })) -->
<!-- ``` -->

## `optim_refit()`: A Numerical Alternative to Projection Plots

There are some situations in which numerical quantification of projection plots leaves to be desired:

Generating all projection plots requires `N = 2 * npts * length(xsol)` evaluations of the objective function (where the default value is `npts = 100`), which can belabor the process of automated unit testing.  A different test for mode-finding routines is to recalculate the optimal solution with an "very good" starting point: the current potential solution.  This is the so-called "**refi**ne op**t**izimation" -- or `refit` -- strategy.  

The `optim_refit()` function refines the optimization with a call to **R**'s built-in general-purpose optimizer: the function `stats::optim()`.  In particular, it selects the default Nelder-Mead simplex method with a simplified parameter interface.  As seen in the unit tests above, the `refit` checks are 2-3 times faster than their projection plot counterparts.  Consider now the example of refining the original `stats::optim()` solution to the quadratic objective function:
```{r}
orefit <- optim_refit(fun = objfun,        # objective function
                      xsol = xfit$par,     # potential solution
                      maximize = FALSE)    # indicates that a local minimum is sought
summary(orefit) # same print method as orefit
```
Thus we can see that the first and second run of `stats::optim()` are quite different. 

Of course, this does not mean that the refit solution produced by `stats::optim()` is a local mode:
```{r, fig.width = 10, fig.height = 6, out.width = "97%"}
# projection plots with refined solution
optim_proj(xsol = orefit$xopt, fun = objfun,
           xrng = .5, maximize = FALSE)
```
Indeed, the default `stats::optim()` method is only accurate when initialized close to the optimal solution.  Therefore, one may wish to run the refit test with a different optimizer.  This can be done externally to `optim_refit`, prior to passing the refit solution to the function via its argument `xopt`.  This is illustrated below using `stats::optim()`'s gradient-based quasi-Newton method:
```{r}
# gradient of the objective function
objgrad <- function(x) 2 * drop(A %*% x - b)

# mode-finding using quasi-Newton method
xfit2 <- optim(fn = objfun,                    # objective function
               gr = objgrad,                   # gradient
               par = xfit$par,                 # initial value (first optim fit)
               method = "BFGS")

# external refit test with optimizer of choice
orefit2 <- optim_refit(fun = objfun,
                       xsol = xfit$par,        # initial value (first optim fit)
                       xopt = xfit2$par,       # refit value (2nd fit with quasi-Newton method
                       maximize = FALSE)

# project plot test on refit solution
optim_proj(xsol = orefit2$xopt, fun = objfun,
           xrng = .5, maximize = FALSE, plot = FALSE)
```

<!-- the elementwise optima  -->
<!-- - It is difficult to combine local optima of different projection plots into a global optimal solution -->
<!-- for package development. (as with (http://r-pkgs.had.co.nz/tests.html)[**`testthat`**]),  -->

## Future Work: Constrained Optimization

Many constrained statistical optimization problems, seek a "sparse" solution, i.e., one for which some of the elements of the optimal solution are equal to zero.  In such cases, the relative difference between potential and optimal solution is an unreliable metric.  A working proposal is to flag these "true zeros" in `optim_proj()` and `optim_refit()`, so as to add a 1 to the relative difference denominators.  Other suggestions on this and **optimCheck** in general are [welcome](mailto:mlysy@uwaterloo.ca).
