Skip to content

Quadratic Optimization

4 messages · Martin Maechler, Brian Ripley, Spencer Graves

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

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

      S = Lagrange objective = a'x+lam*(x'Ax-c). 
      dS/dx = a + 2*lam*Ax. 
     
      Given lam, x1 = solve(A, a)/(2*lam) 
        
      Then x = c*x1/(x'Ax) 

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

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

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

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

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

      How's this? 
      Spencer Graves
amit soni wrote:
Martin Maechler
# 
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
Brian Ripley
#
On Sat, 2 Dec 2006, Martin Maechler wrote:

            

      
      
      
    

        
But that is a single equality quadratic constraint, and I believe 
'quadratic constraints' (note, plural) conventionally means multiple 
inequality constraints.  That meaning is a hard problem that needs 
specialized software (most likely using interior-point methods).

            

      
      
      
    

        
Not I believe the usual meaning (nor what Googling 'quadratic constraints' 
came up with for me).

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Spencer Graves
    

    

      
      #
    

  


  

      
Hi, Prof. Ripley: 

<snip>

            

      
      
      
    

        
Thanks for the clarification.  Spencer Graves