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. -----Original Message----- From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of Hans W Borchers Sent: Wednesday, August 26, 2015 6:22 AM To: r-help at stat.math.ethz.ch Subject: Re: [R] lsqlin in R package pracma On Mon Aug 24 Wang, Xue, Ph.D. Wang.Xue at mayo.edu wrote> I am looking for a R version of Matlab function lsqlin. I came across > R pracma package which has a lsqlin function. Compared with Matlab lsqlin, > the R version does not allow inequality constraints. > I am wondering if this functionality will be available in future. And also > like to get your opinion on which R package/function is the best forsolving> least square minimization problem with linear inequality constraints. > Thanks very much for your time and attention!Solving (linear) least-squares problems with linear inequality constraints is more difficult then one would expect. Inspecting the MATLAB code reveals that it employs advanced methods such as active-set (linear inequality constraints) and interior-point (for bounds constraints). Function nlsLM() in package *minpack.lm* supports bound constraints if that is sufficient for you. The same is true for *nlmrt*. Convex optimization might be a promising approach for linear inequality constraints, but there is no easy-to-handle convex solver in R at this moment. So the most straightforward way would be to use constrOptim(), that is optim with linear constraints. It requires a reasonable starting point, and keeping your fingers crossed that you are able to find such a point in the interior of the feasible region. I someone wants to try: Here is the example from the MATLAB "lsqlin" page: 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) The least-square function to be minimized is ||C x - d||_2 , and the constraints are A x <= b : f <- function(x) sum((C %*% x - d)^2) The solution x0 returned by MATLAB has a minimum of f(x0) = 0.01759204 . This point does not lie in the interior and cannot be used for a start. [[alternative HTML version deleted]] ______________________________________________ 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. Notice: This e-mail message, together with any attachme...{{dropped:11}}
Hi Richard, It is good to know that solve.QP could solve quadratic programming problem. The difficulty here is that the objective function might not be in quadratic form. It is not in the form of t(X)QX, where Q is an n by n symmetric matrix. Thanks, Xue -----Original Message----- From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of Raubertas, Richard Sent: Thursday, August 27, 2015 12:07 PM To: Hans W Borchers; r-help at stat.math.ethz.ch Subject: Re: [R] lsqlin in R package pracma 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. -----Original Message----- From: R-help [mailto:r-help-bounces at r-project.org] On Behalf Of Hans W Borchers Sent: Wednesday, August 26, 2015 6:22 AM To: r-help at stat.math.ethz.ch Subject: Re: [R] lsqlin in R package pracma On Mon Aug 24 Wang, Xue, Ph.D. Wang.Xue at mayo.edu wrote> I am looking for a R version of Matlab function lsqlin. I came across > R pracma package which has a lsqlin function. Compared with Matlab > lsqlin, the R version does not allow inequality constraints. > I am wondering if this functionality will be available in future. And > also like to get your opinion on which R package/function is the best > forsolving> least square minimization problem with linear inequality constraints. > Thanks very much for your time and attention!Solving (linear) least-squares problems with linear inequality constraints is more difficult then one would expect. Inspecting the MATLAB code reveals that it employs advanced methods such as active-set (linear inequality constraints) and interior-point (for bounds constraints). Function nlsLM() in package *minpack.lm* supports bound constraints if that is sufficient for you. The same is true for *nlmrt*. Convex optimization might be a promising approach for linear inequality constraints, but there is no easy-to-handle convex solver in R at this moment. So the most straightforward way would be to use constrOptim(), that is optim with linear constraints. It requires a reasonable starting point, and keeping your fingers crossed that you are able to find such a point in the interior of the feasible region. I someone wants to try: Here is the example from the MATLAB "lsqlin" page: 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) The least-square function to be minimized is ||C x - d||_2 , and the constraints are A x <= b : f <- function(x) sum((C %*% x - d)^2) The solution x0 returned by MATLAB has a minimum of f(x0) = 0.01759204 . This point does not lie in the interior and cannot be used for a start. [[alternative HTML version deleted]] ______________________________________________ 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. Notice: This e-mail message, together with any attachme...{{dropped:11}} ______________________________________________ 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.
> On 27 Aug 2015, at 23:12, Wang, Xue, Ph.D. <Wang.Xue at mayo.edu> wrote: > > Hi Richard, > > It is good to know that solve.QP could solve quadratic programming problem. The difficulty here is that the objective function might not be in quadratic form. It is not in the form of t(X)QX, where Q is an n by n symmetric matrix.Unless I?m very mistaken the objective function is in the form you mention. The quadratic part is t(x) %*% t(C) %*% C %*% x so your Q is simply equivalent to t(C) %*% C. Berend
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.
>
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.