[R] lsqlin in R package pracma

[Hans W Borchers]

On Tue, Sep 1, 2015 at 11:24 PM, Wang, Xue, Ph.D. <Wang.Xue at mayo.edu> wrote:
>
> slsqp in R seems quite slow. Does anyone have some suggestion as how to speed up this?
>

It is no surprise that a general solver like slsqp() takes longer than
specialized quadratic solvers such as solve.QP, ipop(), or
lowRankQP(). You got good advice in this thread, even with code; so
why don't you use what has been suggested in this thread?

If in your timings you include calling lsqlin() in pracma, then this
will increase the time by at least 0.4 ms.

NLopt is an external, compiled program that cannot be accelerated from within R.

I wrote a small wrapper for all the different approaches proposed.
Here is a rough timing estimate (with microbenchmark, in
milliseconds), all derived from the MATLAB example:

                    with constraints    w/o constraints
    ---------------------------------------------------
    qr.solve        --                  0.06 ms
    slsqp           5.5   ms            5.25 ms     [*]
    solve.QP        0.045 ms            --
    ipop            5.0   ms            --
    lowRankQP       0.16  ms            --
    ---------------------------------------------------
    [*] plus 0.4 ms when using lsqlin for start values.

Because solve.QP is fastest, I append a function mimicking the MATLAB
API of lsqlin (without error handling for now) that can be used with
and without equality constraints, utilizing Berend's contribution. If
this is too slow, I guess you will have to look outside R, for example
Julia provides convex solvers that may be fast.

--- cut ----------------------------------------------------------------------
require(quadprog)
constrLsqlin <- function(C, d, A, b, # linear inequalities required
                         Aeq = NULL, beq = NULL, lb = NULL, ub = NULL) {
    m <- nrow(C); n <- ncol(C)
    meq <- nrow(Aeq); neq <- ncol(Aeq)
    if (is.null(meq)) meq = 0

    Dmat <- t(C) %*% C
    dvec <- t(C) %*% d

    if (is.null(Aeq)) {
        Amat <- -A
        bvec <- -b
    } else {
        Amat <- rbind(Aeq, -A)
        bvec <- c(beq, -b)
    }
    if (!is.null(lb)) {
        Amat <- rbind(Amat, diag(n))
        bvec <- c(bvec, lb)
    }
    if (!is.null(ub)) {
        Amat <- rbind(Amat, -diag(n))
        bvec <- c(bvec, -ub)
    }
    rslt <- solve.QP(Dmat, dvec, t(Amat), bvec, meq=meq)
    rslt$solution
}