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.
