# Solve a Quadratic Programming Problem via Dykstra's Algorithm This function uses Dykstra's cyclic projection algorithm to solve quadratic programming problems of the form [REMOVE_ME]$$ - d^T x + (1/2) x^T D x [REMOVE_ME_2]$$ subject to $A^T x >= b$ where $D$ is a positive definite (or positive semidefinite) matrix. ## Description This function uses Dykstra's cyclic projection algorithm to solve quadratic programming problems of the form $$ - d^T x + (1/2) x^T D x $$ subject to $A^T x >= b$ where $D$ is a positive definite (or positive semidefinite) matrix. ```r dykstra(Dmat, dvec, Amat, bvec, meq = 0, factorized = FALSE, maxit = NULL, eps = NULL) ``` ## Arguments - `Dmat`: Quadratic program matrix $D$ of order $n \times n$. - `dvec`: Quadratic program vector $d$ of length $n$. - `Amat`: Constraint matrix $A$ of order $n \times r$. - `bvec`: Constraint vector $b$ of length $r$. Defaults to vector of zeros. - `meq`: First `meq` constraints are equality constraints (remaining are inequality constraints). Defaults to zero. - `factorized`: If `TRUE`, argument `Dmat` is $R^{-1}$ where $R^T R = D$. - `maxit`: Maximum number of iterations (cycles). Defaults to $30n$. - `eps`: Numeric tolerance. Defaults to `n * .Machine$double.eps`. ## Details Arguments 1-6 of the `dykstra` function are inspired by (and identical to) the corresponding arguments of the `solve.QP` function in the quadprog package. ## Returns - **solution**: Vector $x$ that minimizes quadratic function subject to constraints. - **value**: Value of quadratic function at `solution`. Will be `NA` if `factorized = TRUE`. - **unconstrained**: Vector $x_0 = D^{-1} d$ that minimizes quadratic function ignoring constraints. - **iterations**: Number of iterations (cycles) of the algorithm. - **converged**: `TRUE` if algorithm converged. `FALSE` if iteration limit exceeded. ## References Dykstra, Richard L. (1983). An algorithm for restricted least squares regression. **Journal of the American Statistical Association**, Volume 78, Issue 384, 837-842. doi: 10.1080/01621459.1983.10477029 ## Author(s) Nathaniel E. Helwig ## Note For positive semidefinite $D$, a small constant is added to each eigenvalue of $D$ before solving the quadratic programming problem. ## Examples ```r ### EXAMPLE 1: Generic Quadratic Programming Problem ### # constraint 1 (equality): coefficients sum to 1 # constraints 2-4 (inequality): coefficients non-negative # define QP problem Dmat <- diag(3) dvec <- c(1, 1.5, 1) Amat <- cbind(rep(1, 3), diag(3)) bvec <- c(1, 0, 0, 0) # solve QP problem dykstra(Dmat, dvec, Amat, bvec, meq = 1) # solve QP problem (factorized = TRUE) dykstra(Dmat, dvec, Amat, bvec, meq = 1, factorized = TRUE) ### EXAMPLE 2: Regression with Non-Negative Coefficients ### # generate regression data set.seed(1) nobs <- 100 nvar <- 5 X <- matrix(rnorm(nobs*nvar), nobs, nvar) beta <- c(0, 1, 0.3, 0.7, 0.1) y <- X %*% beta + rnorm(nobs) # define QP problem Dmat <- crossprod(X) dvec <- crossprod(X, y) Amat <- diag(nvar) # solve QP problem dykstra(Dmat, dvec, Amat) # solve QP problem (factorized = TRUE) Rmat <- chol(Dmat) Rinv <- solve(Rmat) dykstra(Rinv, dvec, Amat, factorized = TRUE) ### EXAMPLE 3: Isotonic Regression ### # generate regression data set.seed(1) n <- 50 x <- 1:n y <- log(x) + rnorm(n) # define QP problem Dmat <- diag(n) Amat <- Dmat[, 2:n] - Dmat[, 1:(n-1)] # solve QP problem dyk <- dykstra(Dmat, y, Amat) dyk # plot results plot(x, y) lines(x, dyk$solution) ### EX 4: Large Non-Negative Quadratic Program ### # define QP problem set.seed(1) n <- 1000 Dmat <- Amat <- diag(n) dvec <- runif(n, min = -2) # solve QP problem with dykstra dyk <- dykstra(Dmat, dvec, Amat) dyk ```

dykstra function