knitr::opts_chunk$set(collapse = TRUE,results = "hold",comment = "#", fig.align = "center",eval = FALSE)
In this vignette, we illustrate the use of the sequential quadratic
programming (SQP) algorithm implemented in mixsqp
.
Load the mixsqp
package.
library(mixsqp)
Next, initialize the sequence of pseudorandom numbers.
set.seed(1)
We begin with a small example to show how mixsqp
works.
L <- simulatemixdata(1000,20)$L dim(L)
This call to simulatemixdata
created an $n \times m$ conditional
likelihood matrix for a mixture of zero-centered normals, with $n =
1000$ and $m = 20$. By default, simulatemixdata
normalizes the rows
of the likelihood matrix so that the maximum entry in each row is 1.
Now we fit the mixture model using the SQP algorithm:
fit.sqp <- mixsqp(L)
In this example, the SQP algorithm converged to a solution in a small number of iterations.
By default, mixsqp
outputs information on its progress. It begins by
summarizing the optimization problem and the algorithm settings used.
(Since we did not change these settings in the mixsqp
call, all the
settings shown here are the default settings.)
After that, it outputs, at each iteration, information about the current solution, such as the value of the objective ("objective") and the number of nonzeros ("nnz").
The "max(rdual)" column shows the quantity used to assess convergence.
It reports the maximum value of the "dual residual"; the SQP solver
terminates when the maximum dual residual is less than conv.tol
,
which by default is $10^{-8}$. In this example, we see that the dual
residual shrinks rapidly toward zero.
Another useful indicator of convergence is the "max.diff" column---it reports the maximum difference between the solution estimates at two successive iterations. We normally expect these differences to shrink as we approach the solution, which is precisely what we see in this example.
This information is also provided in the return value, which we can use, for example, to create a plot of the objective value at each iteration of the SQP algorithm:
numiter <- nrow(fit.sqp$progress) plot(1:numiter,fit.sqp$progress$objective,type = "b", pch = 20,lwd = 2,xlab = "SQP iteration", ylab = "objective",xaxp = c(1,numiter,numiter - 1))
This next code chunk gives information about the computing environment used to generate the results contained in this vignette, including the version of R and the packages used.
sessionInfo()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.