mixsqp

https://stephenslab.github.io/mixsqp/
Admin
mixsqp provides algorithms based on sequential quadratic programming for fast maximum-likelihood...
Tokens:12,840
Snippets:99
Trust Score:9.7
Update:1 week ago
Show doc for...
Tokens:
Raw
### Install and get started with mixsqp in R

Source: https://stephenslab.github.io/mixsqp/index.html

Steps to install and explore the mixsqp package. First install from CRAN, then load the library contact, run a built-in example, and open the help page.

```R
install.packages("mixsqp")
```

```R
library(mixsqp)
```

```R
example("mixsqp")
```

```R
help("mixsqp")
```

--------------------------------

### Build and test mixsqp package with R CMD commands

Source: https://stephenslab.github.io/mixsqp/index.html

Build, install, and check the mixsqp package from the command line in a shell. Requires dependencies to be installed.

```shell
R CMD build mixsqp
R CMD INSTALL mixsqp_0.3-18.tar.gz
R CMD check --as-cran mixsqp_0.3-18.tar.gz
```

--------------------------------

### Install mixsqp from GitHub with devtools

Source: https://stephenslab.github.io/mixsqp/index.html

Install the latest mixsqp version from GitHub using devtools. This command automatically installs required packages.

```R
install.packages("devtools")
library(devtools)
install_github("stephenslab/mixsqp",build_vignettes = TRUE)
```

--------------------------------

### Install local mixsqp repository

Source: https://stephenslab.github.io/mixsqp/index.html

Install mixsqp from a local clone of the repository using the install_local function from devtools

```R
library(devtools)
install_local("mixsqp")
```

--------------------------------

### Build and view mixsqp vignette

Source: https://stephenslab.github.io/mixsqp/index.html

Build the vignette and view it locally after installing REBayes and Rmosek for the optional 'mixkwdual' solver.

```R
devtools::install_github("stephenslab/mixsqp",build_vignettes = TRUE)
library(mixsqp)
vignette("mixsqp-intro")
```

--------------------------------

### Fit IP algorithm with mixkwdual

Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html

Computes the solution using the IP algorithm from the REBayes package. Requires REBayes to be installed; otherwise this step can be skipped.

```R
fit.ip <- mixkwdual(L)
```

--------------------------------

### tacks

Source: https://stephenslab.github.io/mixsqp/reference/index.html

Beckett & Diaconis tack rolling example.

```APIDOC
## tacks

### Description
Beckett & Diaconis tack rolling example.

```

--------------------------------

### Generate likelihood matrix with simulatemixdata

Source: https://stephenslab.github.io/mixsqp/reference/simulatemixdata.html

Example usage of `simulatemixdata` to generate a likelihood matrix for 1,000 samples and 20 mixture components. The result is stored in `dat`.

```R
# Generate the likelihood matrix for a data set with 1,000 samples
# and a nonparametric Empirical Bayes model with 20 mixture
# components.
dat <- simulatemixdata(1000,20)
```

--------------------------------

### Fit mixture proportions with mixsqp() and evaluate with mixobjective()

Source: https://stephenslab.github.io/mixsqp/reference/mixsqp.html

Basic call to mixsqp on a simulated 100000 x 10 matrix, then evaluates the objective with mixobjective. The run prints algorithm settings and a per-iteration table, ending when convergence criteria are met in about 0.5 seconds. The optional comparison with mixkwdual requires MOSEK and REBayes, so it is wrapped in if (FALSE).

```R
__
set.seed(1)
n <- 1e5
m <- 10
w <- rep(1,n)/n
L <- simulatemixdata(n,m)$L
out.mixsqp <- mixsqp(L,w)
#> Running mix-SQP algorithm 0.3-17 on 100000 x 10 matrix
#> convergence tol. (SQP):     1.0e-08
#> conv. tol. (active-set):    1.0e-10
#> zero threshold (solution):  1.0e-08
#> zero thresh. (search dir.): 1.0e-10
#> l.s. sufficient decrease:   1.0e-02
#> step size reduction factor: 7.5e-01
#> minimum step size:          1.0e-08
#> max. iter (SQP):            1000
#> max. iter (active-set):     11
#> number of EM iterations:    10
#> Computing SVD of 100000 x 10 matrix.
#> Matrix is not low-rank; falling back to full matrix.
#> iter        objective max(rdual) nnz stepsize max.diff nqp nls
#>    1 +6.649406794e-01  -- EM --   10 1.00e+00 7.79e-02  --  --
#>    2 +6.336404241e-01  -- EM --   10 1.00e+00 2.61e-02  --  --
#>    3 +6.201431305e-01  -- EM --   10 1.00e+00 2.41e-02  --  --
#>    4 +6.131770842e-01  -- EM --   10 1.00e+00 2.06e-02  --  --
#>    5 +6.091788165e-01  -- EM --   10 1.00e+00 1.71e-02  --  --
#>    6 +6.066948308e-01  -- EM --   10 1.00e+00 1.42e-02  --  --
#>    7 +6.050509636e-01  -- EM --   10 1.00e+00 1.18e-02  --  --
#>    8 +6.039062025e-01  -- EM --   10 1.00e+00 9.95e-03  --  --
#>    9 +6.030755992e-01  -- EM --   10 1.00e+00 8.45e-03  --  --
#>   10 +6.024525768e-01  -- EM --   10 1.00e+00 7.25e-03  --  --
#>    1 +6.024525768e-01 +2.365e-02  10  ------   ------   --  --
#>    2 +6.058450302e-01 +4.073e+00   4 1.00e+00 3.91e-01  11   1
#>    3 +5.999690521e-01 +1.894e+00   5 1.00e+00 3.21e-02  10   1
#>    4 +5.997818561e-01 +7.929e-01   5 1.00e+00 2.95e-04   2   1
#>    5 +5.997500584e-01 +2.749e-01   5 1.00e+00 1.98e-03   2   1
#>    6 +5.996811464e-01 +1.482e-01   5 1.00e+00 5.64e-03   2   1
#>    7 +5.994718234e-01 +6.121e-02   5 1.00e+00 9.54e-03   2   1
#>    8 +5.991459366e-01 +1.505e-02   5 1.00e+00 8.23e-03   2   1
#>    9 +5.989768704e-01 +1.193e-03   5 1.00e+00 2.88e-03   2   1
#>   10 +5.989591648e-01 +8.215e-06   5 1.00e+00 2.42e-04   2   1
#>   11 +5.989590364e-01 +3.904e-10   5 1.00e+00 1.63e-06   2   1
#> Optimization took 0.54 seconds.
#> Convergence criteria met---optimal solution found.
f <- mixobjective(L,out.mixsqp$x,w)
print(f,digits = 16)
#> [1] 0.5989590363728419

# We can also compare this result with solution found from an 
# interior-point approach called via the "KWDual" function from the 
# REBayes package. (This requires installation of the MOSEK 
# optimization library as well as the REBayes package, so we have
# made this step optional.)
if (FALSE) {
out.kwdual <- mixkwdual(L,w)
print(mixobjective(L,out.kwdual$x,w),digits = 16)
}
```

--------------------------------

### Build pkgdown site for mixsqp

Source: https://stephenslab.github.io/mixsqp/index.html

Build the pkgdown website for mixsqp while connected to the Internet. Version 1.4.1 of pkgdown was used.

```R
pkgdown::build_site()
```

--------------------------------

### Load mixsqp package

Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html

Loads the mixsqp package so its functions can be used.

```R
library(mixsqp)
```

--------------------------------

### Compare SQP and IP objective values

Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html

Prints the objective values from the SQP and IP solutions and the difference between them, showing they are nearly identical.

```R
cat(sprintf("Objective at SQP solution: %0.16f\n",fit.sqp$value,digits = 16))
cat(sprintf("Objective at IP solution:  %0.16f\n",fit.ip$value,digits = 16))
cat(sprintf("Difference in objectives:  %0.4e\n",fit.sqp$value - fit.ip$value))
```

--------------------------------

### Run IP algorithm with mixkwdual

Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html

Runs the IP method using the mixkwdual function from the REBayes package, timing the computation and printing the elapsed time.

```R
timing <- system.time(fit.ip  <- mixkwdual(L))
cat(sprintf("Computation took %0.2f seconds\n",timing["elapsed"]))
```

--------------------------------

### mixsqp-package

Source: https://stephenslab.github.io/mixsqp/reference/index.html

mixsqp: Sequential Quadratic Programming for Fast Maximum-Likelihood Estimation of Mixture Proportions.

```APIDOC
## mixsqp-package

### Description
mixsqp: Sequential Quadratic Programming for Fast Maximum-Likelihood Estimation of Mixture Proportions.

```

--------------------------------

### Fit mixture model with mixsqp

Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html

Runs the SQP algorithm on the likelihood matrix. The output shows default settings and progress, including objective, dual residual, and number of nonzeros. Convergence is achieved when the dual residual is below conv.tol (default 1e-8).

```R
fit.sqp <- mixsqp(L)
# Running mix-SQP algorithm 0.3-17 on 1000 x 20 matrix
# convergence tol. (SQP):     1.0e-08
# conv. tol. (active-set):    1.0e-10
# zero threshold (solution):  1.0e-08
# zero thresh. (search dir.): 1.0e-10
# l.s. sufficient decrease:   1.0e-02
# step size reduction factor: 7.5e-01
# minimum step size:          1.0e-08
# max. iter (SQP):            1000
# max. iter (active-set):     20
# number of EM iterations:    10
# Computing SVD of 1000 x 20 matrix.
# Matrix is not low-rank; falling back to full matrix.
# iter        objective max(rdual) nnz stepsize max.diff nqp nls
#    1 +6.825854577e-01  -- EM --   20 1.00e+00 3.43e-02  --  --
#    2 +6.608901451e-01  -- EM --   20 1.00e+00 1.12e-02  --  --
#    3 +6.501638059e-01  -- EM --   20 1.00e+00 8.83e-03  --  --
#    4 +6.441429919e-01  -- EM --   20 1.00e+00 7.64e-03  --  --
#    5 +6.405380236e-01  -- EM --   20 1.00e+00 6.44e-03  --  --
#    6 +6.382624101e-01  -- EM --   20 1.00e+00 5.36e-03  --  --
#    7 +6.367521107e-01  -- EM --   20 1.00e+00 4.46e-03  --  --
#    8 +6.357010188e-01  -- EM --   20 1.00e+00 3.75e-03  --  --
#    9 +6.349367202e-01  -- EM --   20 1.00e+00 3.18e-03  --  --
#   10 +6.343585101e-01  -- EM --   20 1.00e+00 2.73e-03  --  --
#    1 +6.343585101e-01 +2.161e-02  20  ------   ------   --  --
#    2 +6.331717425e-01 +8.779e-03   4 1.00e+00 4.36e-01  20   1
#    3 +6.282364710e-01 +8.715e-05   4 1.00e+00 3.39e-03   2   1
#    4 +6.281978492e-01 +9.817e-09   4 1.00e+00 4.13e-05   2   1
# Optimization took 0.01 seconds.
# Convergence criteria met---optimal solution found.
```

--------------------------------

### mixsqp_control_default()

Source: https://stephenslab.github.io/mixsqp/reference/index.html

Default control parameters for mixsqp.

```APIDOC
## mixsqp_control_default()

### Description
Default control parameters for mixsqp.

### Usage
mixsqp_control_default()

```

--------------------------------

### Plot optimal solution for tacks data

Source: https://stephenslab.github.io/mixsqp/reference/tacks.html

Loads the tacks dataset and plots the optimal solution, which is noted for being extremely sparse.

```R
data(tacks)
plot(tacks$x,type = "l",col = "royalblue")
```

--------------------------------

### mixsqp

Source: https://stephenslab.github.io/mixsqp/reference/mixsqp.html

Solves the convex optimization problem of minimizing a weighted sum of negative log-likelihoods subject to a simplex constraint, using a sequential quadratic programming (SQP) algorithm.

```APIDOC
## mixsqp

### Description
Solves the convex optimization problem of minimizing a weighted sum of negative log-likelihoods subject to a simplex constraint, using a sequential quadratic programming (SQP) algorithm.

### Method
Function call

### Endpoint
mixsqp(L, w, control = list())

### Parameters
#### Arguments
- **L** (matrix) - Required - A matrix with n rows and m columns containing only non-negative entries.
- **w** (numeric vector) - Required - A vector of non-negative weights of length n.
- **control** (list) - Optional - A list of control parameters that override the default optimization algorithm settings. See Details for available components.

### Return Value
A list with the following elements:
- **x** (numeric vector) - The solution to the convex optimization problem if the SQP algorithm converges; otherwise, the best estimate of the solution achieved by the algorithm. Note that if the algorithm fails to converge, the entries of `x` may not sum to 1.
- **value** (numeric) - The value of the objective function, f(x), at `x`.
- **status** (character) - A character string describing the status of the algorithm upon termination.
- **progress** (data.frame) - A data frame containing detailed information about the algorithm's progress, with one row per SQP iteration. See the description of the `verbose` control parameter in Details for an explanation of the columns.

### Details
The `mixsqp` function solves the following optimization problem. Let L be a matrix with n rows and m columns containing only non-negative entries, and let w=(w1,…,wn) be a vector of non-negative "weights". `mixsqp` computes the value of vector x=(x1,…,xm) minimizing the objective function f(x)=−∑j=1n wj log(∑k=1m Ljk xk), subject to the constraint that x lie within the simplex; that is, the entries of x are non-negative and sum to 1. Implicitly, there is an additional constraint L∗x>0 to ensure the objective has a finite value.

If all weights are equal, solving this optimization problem corresponds to finding the maximum-likelihood estimate of the mixture proportions x given n independent data points drawn from a mixture model with m components. In this case, Ljk is the likelihood for mixture component k and data point j.

The Expectation Maximization (EM) algorithm can be used to solve this optimization problem, but it is intolerably slow in many interesting cases, and mixsqp is much faster.

A special feature of this optimization problem is that the gradient of the objective does not change with re-scaling; for example, if all the entries of matrix `L` are multiplied by 100, the gradient does not change. A practical benefit of this property is that the optimization algorithm will behave similarly irrespective of the scale of `L`; for example, the same value for the convergence tolerance `convtol.sqp` will have the same effect at different scales.

A related feature is that the solution to the optimization problem is invariant to rescaling the rows of `L`; for example, the solution will remain the same after all the entries in a row of `L` are multiplied by 10. A simple normalization scheme divides each row by the largest entry in the row so that all entries of `L` are at most 1: `L <- L / apply(L,1,max)` Occasionally, it can be helpful to normalize the rows when some of the entries are unusually large or unusually small. This can help to avoid numerical overflow or underflow errors.

The SQP algorithm is implemented using the Armadillo C++ linear algebra library, which can automatically take advantage of multithreaded matrix computations to speed up `mixsqp` for large `L` matrices, but only when R has been configured with a multithreaded BLAS/LAPACK library (e.g., OpenBLAS).

The `control` argument is a list in which any of the following named components will override the default optimization algorithm settings (as they are defined by `mixsqp_control_default`):

- **normalize.rows** (logical) - When `normalize.rows = TRUE`, the rows of the data matrix `L` are automatically scaled so that the largest entry in each row is 1. This is the recommended setting for better stability of the optimization. When `log = TRUE`, this setting is ignored because the rows are already normalized. Note that the objective is computed on the original (unnormalized) matrix to make the results easier to interpret.
- **force.sparse.init** (logical) - The computational complexity of the SQP updates grows rapidly by the number of nonzeros in the current solution estimate. When `force.sparse.init = TRUE`, and there are more than 20 nonzeros in the solution before beginning the SQP updates, only the 20 largest entries are kept, the rest are forced to zero.
- **tol.svd** (numeric) - Setting used to determine rank of truncated SVD approximation for L. The rank of the truncated singular value decomposition is determined by the number of singular values surpassing `tol.svd`. When `tol.svd = 0` or when `L` has 4 or fewer columns, all computations are performed using full L matrix.
```

--------------------------------

### mixobjective()

Source: https://stephenslab.github.io/mixsqp/reference/index.html

Compute objective optimized by mixsqp.

```APIDOC
## mixobjective()

### Description
Compute objective optimized by mixsqp.

### Usage
mixobjective()

```

--------------------------------

### Generate small data set

Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html

Generates a small data set using simulatemixdata and prints its dimensions.

```R
L <- simulatemixdata(1000,20)$L
dim(L)
```

--------------------------------

### mixsqp()

Source: https://stephenslab.github.io/mixsqp/reference/index.html

Maximum-likelihood estimation of mixture proportions using SQP.

```APIDOC
## mixsqp()

### Description
Maximum-likelihood estimation of mixture proportions using SQP.

### Usage
mixsqp()

```

--------------------------------

### Show session information

Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html

Displays the R session information, including the version of R, platform, and packages used to generate the vignette results.

```R
sessionInfo()
# R version 3.6.2 (2019-12-12)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS Catalina 10.15.2
# 
# Matrix products: default
# BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
# 
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# 
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# 
# other attached packages:
# [1] mixsqp_0.3-18
# 
# loaded via a namespace (and not attached):
#  [1] Rcpp_1.0.3       knitr_1.26       magrittr_1.5     REBayes_1.8     
#  [5] MASS_7.3-51.4    lattice_0.20-38  R6_2.4.1         rlang_0.4.2     
#  [9] stringr_1.4.0    tools_3.6.2      grid_3.6.2       xfun_0.11       
# [13] irlba_2.3.3      htmltools_0.4.0  yaml_2.2.0       assertthat_0.2.1
# [17] digest_0.6.23    rprojroot_1.3-2  pkgdown_1.4.1    crayon_1.3.4    
# [21] Matrix_1.2-18    fs_1.3.1         memoise_1.1.0    evaluate_0.14   
# [25] rmarkdown_2.0    stringi_1.4.3    compiler_3.6.2   desc_1.2.0      
# [29] backports_1.1.5  Rmosek_9.0.96
```

--------------------------------

### simulatemixdata()

Source: https://stephenslab.github.io/mixsqp/reference/index.html

Create likelihood matrix from simulated data set.

```APIDOC
## simulatemixdata()

### Description
Create likelihood matrix from simulated data set.

### Usage
simulatemixdata()

```

--------------------------------

### Set random seed

Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html

Initializes the sequence of pseudorandom numbers for reproducible results.

```R
set.seed(1)
```

--------------------------------

### mixobjective(L, x, w = rep(1, nrow(L)))

Source: https://stephenslab.github.io/mixsqp/reference/mixobjective.html

Computes the value of the objective function optimized by the mix-SQP algorithm at the point x. The objective is defined by the matrix L and optional weights w. Returns Inf if any entry of L %*% x is non-positive.

```APIDOC
## mixobjective(L, x, w = rep(1, nrow(L)))

### Description
Computes the value of the objective function optimized by the mix-SQP algorithm at the point x. The objective is defined by the matrix L and optional weights w. Returns Inf if any entry of L %*% x is non-positive.

### Usage
mixobjective(L, x, w = rep(1, nrow(L)))

### Arguments
- **L** (matrix) - Required - Matrix specifying the optimization problem. In the context of mixture-model fitting, L[j,k] should be the value of the kth mixture component density at the jth data point. L should be a numeric matrix with at least two columns, with all entries being non-negative and finite (and not missing). No column should be entirely zeros. For large matrices, it is preferable that the matrix is stored in double-precision; see storage.mode.
- **x** (numeric) - Required - The point at which the objective is evaluated; see argument x0 in mixsqp for details.
- **w** (numeric vector) - Optional - An optional numeric vector, with one entry for each row of L, specifying the "weights" associated with the rows of L. All weights must be finite, non-negative and not missing. Internally, the weights are normalized to sum to 1, which does not change the problem, but does change the value of the objective function reported. By default, all weights are equal.

### Value
The value of the objective at x. If any entry of L %*% x is less than or equal to zero, Inf is returned.

### See Also
mixsqp
```

--------------------------------

### Plot objective progress from mixsqp fit

Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html

Uses the progress data stored in fit.sqp$progress to plot the objective value at each SQP iteration.

```R
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))
```

--------------------------------

### mixkwdual()

Source: https://stephenslab.github.io/mixsqp/reference/index.html

Maximum-likelihood estimation of mixture proportions using MOSEK.

```APIDOC
## mixkwdual()

### Description
Maximum-likelihood estimation of mixture proportions using MOSEK.

### Usage
mixkwdual()

```

--------------------------------

### Simulate mixture data with simulatemixdata

Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html

Creates an n×m conditional likelihood matrix for a mixture of zero-centered normals. By default, rows are normalized so the maximum entry in each row is 1.

```R
L <- simulatemixdata(1000,20)$L
dim(L)
# [1] 1000   20
```

--------------------------------

### Simulate larger mixture data set

Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html

Creates a larger likelihood matrix with 100,000 samples and 100 mixture components for performance comparison.

```R
L <- simulatemixdata(1e5,100)$L
dim(L)
# [1] 100000    100
```

--------------------------------

### Time SQP fit on large data set

Source: https://stephenslab.github.io/mixsqp/articles/mixsqp-intro.html

Measures the computation time for fitting the SQP model on the large data set and prints the elapsed time.

```R
timing <- system.time(fit.sqp <- mixsqp(L))
cat(sprintf("Computation took %0.2f seconds\n",timing["elapsed"]))
```

--------------------------------

### simulatemixdata function signature

Source: https://stephenslab.github.io/mixsqp/reference/simulatemixdata.html

The function signature for `simulatemixdata`, which simulates data and computes a likelihood matrix. It takes parameters for sample size, number of mixture components, simulation type, log transformation, and row normalization.

```R
simulatemixdata(
  n,
  m,
  simtype = c("n", "nt"),
  log = FALSE,
  normalize.rows = !log
)
```
Tokens:
Raw