[R] lsqlin in R package pracma

[Hans W Borchers]

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