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