[prev in list] [next in list] [prev in thread] [next in thread] 

List:       r-help
Subject:    Re: [R] lsqlin in R package pracma
From:       Berend Hasselman <bhh () xs4all ! nl>
Date:       2015-08-28 14:48:52
Message-ID: 7238F988-D7F5-4C42-9931-F0B1C67242D3 () xs4all ! nl
[Download RAW message or body]



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@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@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@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@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.

______________________________________________
R-help@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.


[prev in list] [next in list] [prev in thread] [next in thread]