[R] Quadratic Optimization

[Martin Maechler]

>>>>> "SpG" == Spencer Graves <spencer.graves at pdf.com>
>>>>>     on Fri, 01 Dec 2006 17:29:56 -0800 writes:

    SpG>       Unless I'm missing something, optimizing a linear
    SpG> function with quadratic constraints is almost trivial
    SpG> with Langrange multipliers.

yes. Good point, let's hope we're not solving someone's homework
here :-)

    SpG>       Maximize a'x subject to x'Ax=c.

    SpG>       S = Lagrange objective = a'x+lam*(x'Ax-c).
    SpG>         dS/dx = a + 2*lam*Ax.

    SpG>       Given lam, x1 = solve(A, a)/(2*lam)

  [you forgot a "-" , but it's irrelevant, since 'lam' arbitrarily varies;
   and we can use  x1 = solve(A,a) in the following ]

    SpG>       Then x = c*x1/(x'Ax)

but the last equation above is not useful ("x" on the RHS).
Rather, we know that  x = s * x1 (for some scalar 's') and from
the constraint x'Ax = c (or equivalently, from dS/d{lam} = 0),
we find that

   s = +/- sqrt(c / (x1' A x1))

and can even see that  x1' A x1 = x1' A A^{-1} a =  x1' a

Note the "+/-" : we get the minimizing and the maximizing values

    SpG>       In R, you need to know that "t" = transpose of a
    SpG> matrix.

In other words, as an "algorithm" in R



solveLinQuad <- function(a,A,c)
{
    ## Purpose: maximize a'x  under constraint  x'Ax = c
    ##       we return x = arg.max{...}, but note that  arg.min{..} = -x
    ## -------------------------------------------------------------------
    ## Author: Martin Maechler, Date:  2 Dec 2006

    ## argument checking ...
    stopifnot(is.matrix(A), nrow(A) == ncol(A), length(a) == nrow(A),
              is.numeric(A), is.numeric(a), is.numeric(c))
    x1 <- solve(A,a)
    ## note that  x1' A x1 = x1' A A^{-1} a = x1' a { == sum(x1 * a) }
    x <- sqrt(c / sum(x1 * a)) * x1
    ## we return the arg.max :
    if(sum(x * a) >= 0) x else -x
}

## Test:
set.seed(1)
A <- crossprod(matrix(round(rnorm(24),2), 6,4))
a <- c(-1, 1, -2, 6)
c <- 3
(x <- solveLinQuad(a,A,c))
##    -------------
sum(a * x) # 4.677
## constraint
sum(x * (A %*% x)) ## -> 3

    SpG>       I thought I had seen mention of a contributed
    SpG> package for optimization with nonlinear constraints.
    SpG> However, you don't need that here.

    SpG>       In case this does not solve your problem, my
    SpG> crude engineer's approach to constrained optimization
    SpG> includes the following:

    SpG>       (1) Find transformation to send the constraints
    SpG> to +/-Inf.

    SpG>       (2) If that fails, add the constraints as a
    SpG> penalty.  Start with a low penalty and gradually
    SpG> increase it if necessary until you solve the problem.
    SpG> Of course, to do this, you have to make sure your
    SpG> objective function returns valid numbers outside the
    SpG> constrained region.

    SpG>       How's this?  Spencer Graves

    SpG> amit soni wrote:
    >> Hi,
    >>
    >> I need to solve an optimization problem in R having
    >> linear objective function and quadratic
    >> constraints(number of variables is around 80).  What are
    >> the possible choices to do this in R.  optim() function
    >> only allows box constrained problems. Is it possible in
    >> nlm()? Or please tell me if there is any other routine.
    >>
    >> Thanks
    >>
    >> Amit