[R] lsqlin in R package pracma

[Berend Hasselman]

Nice and interesting! Something to remember.

Lesson (for me):

Always first look in Task Views on CRAN.
Choose Optimization and look at Mathematical Programming Solvers.

Berend

> On 28 Aug 2015, at 12:56, Hans W Borchers <hwborchers at gmail.com> wrote:
> 
> I got interested in enabling the full funcionality that MATLAB's
> lsqlin() has, that is with equality and bound constraints. To replace
> an equality constraint with two inequality constraints will not work
> with solve.QP() because it requires positive definite matrices. I will
> use kernlab::ipop() instead.
> 
> To handle the full MATLAB example, add the following simple linear
> equality constraint  3x1 + 5x2 + 7x3 + 9x4 = 4 to the example above,
> plus lower and upper bounds -0.1 and 2.0 for all x_i.
> 
>    C <- matrix(c(
>        0.9501,   0.7620,   0.6153,   0.4057,
>        0.2311,   0.4564,   0.7919,   0.9354,
>        0.6068,   0.0185,   0.9218,   0.9169,
>        0.4859,   0.8214,   0.7382,   0.4102,
>        0.8912,   0.4447,   0.1762,   0.8936), 5, 4, byrow=TRUE)
>    d <- c(0.0578, 0.3528, 0.8131, 0.0098, 0.1388)
>    A <- matrix(c(
>        0.2027,   0.2721,   0.7467,   0.4659,
>        0.1987,   0.1988,   0.4450,   0.4186,
>        0.6037,   0.0152,   0.9318,   0.8462), 3, 4, byrow=TRUE)
>    b <- c(0.5251, 0.2026, 0.6721)
> 
>    # Add the equality constraint to matrix A
>    Aeq <- c(3, 5, 7, 9)
>    beq <- 4
>    A1 <- rbind(A ,  c(3, 5, 7, 9))
>    b1 <- c(b, 4)
>    lb <- rep(-0.1, 4)   # lower and upper bounds
>    ub <- rep( 2.0, 4)
>    r1 <- c(1, 1, 1, 0)  # 0 to force an equality constraint
> 
>    # Prepare for a quadratic solver
>    Dmat <- t(C) %*% C
>    dvec <- (t(C) %*% d)
>    Amat <- -1 * A1
>    bvec <- -1 * b1
> 
>    library(kernlab)
>    s <- ipop(-dvec, Dmat, Amat, bvec, lb, ub, r1)
>    s
>    # An object of class "ipop"
>    # Slot "primal":
>    # [1] -0.09999885 -0.09999997  0.15990817  0.40895991
>    # ...
> 
>    x <- s at primal           # [1] -0.1000  -0.1000  0.1599  0.4090
>    A1 %*% x - b1 <= 0      # i.e., A x <= b and 3x[1] + ... + 9x[4] = 4
>    sum((C %*% x - d)^2)    # minimum: 0.1695
> 
> And this is exactly the solution that lsqlin() in MATLAB computes.
> 
> 
> On Thu, Aug 27, 2015 at 6:06 PM, Raubertas, Richard
> <richard_raubertas at merck.com> wrote:
>> Is it really that complicated?  This looks like an ordinary quadratic programming problem, and 'solve.QP' from the 'quadprog' package seems to solve it without user-specified starting values:
>> 
>> library(quadprog)
>> Dmat <- t(C) %*% C
>> dvec <- (t(C) %*% d)
>> Amat <- -1 * t(A)
>> bvec <- -1 * b
>> 
>> rslt <- solve.QP(Dmat, dvec, Amat, bvec)
>> sum((C %*% rslt$solution - d)^2)
>> 
>> [1] 0.01758538
>> 
>> Richard Raubertas
>> Merck & Co.
>> 
> 
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.