How to: Use Python to Solve Optimization Problems

Some Python programmers may be interested in learning how to use Python, and various supporting packages, to solve mathematics problems frequently encountered by social scientists. Depending on the level of interest, I will make more of these tutorials to tackle other problems.
Optimization problems are a class frequently encountered by social scientists, and therefore this is where I will begin. Optimization techniques can answer questions such as, “if we spend x on schools, y on hospitals and z on infrastructure in a attempt gain confidence from some population P, what is the optimal amount to spend in each to gain maximum confidence?” In a step-by-step fashion, this tutorial will explain how to use Python and SymPy to solve these types of problems.

1. If you have not already, you will need to download and install the Python IDLE and SymPy compatible with your operating system. Each of these software packages has their own installation instructions, review them carefully before installing. Also, if you are totally new to Python, I recommend the Dive Into Python tutorial. Once you have the software up and running we are ready to begin!
2. Problem: Find the critical points for the function below, then classify these critical points as either local maximum, local minimum, saddle point, or neither.
   
   
   First, we must find the first-order conditions of this function by setting up the gradient. To do this using Python, we must first import SymPy into the IDLE and then add the appropriate variables to memory.
3. Open the Python IDLE and type the following commands (all Python commands will be in block-quotes):
   

   >>> import sympy as S
   >>> x,y=S.symbols(‘xy’)

   The first line loads the entire SymPy package into the variable S, while the second line adds the variables x and y into memory using SymPy’s symbolic format so they can be manipulated.
4. The gradient of this function will be the partial derivatives of both variables, x and y.
   
   To solve this with Python, we will first store the function as a variable, then use SymPy’s diff function to perform partial differentiation on f.
5. We store the function in a variable to prevent ourselves from having to repeatedly enter the function into IDLE.
   

   >>> f=x**4+x**2-6*x*y+3*y**2
   >>> f
   -6*x*y + x**2 + 3*y**2 + x**4

   You will notice two things; first, the ** is Python’s is exponent operator, that is, if you are calculating 4² you would type 4**2 into the interpreter; second, SymPy has rearranged the function. This rearranging occurs because SymPy automatically reorders variables from most negative to most positive rank. Now we perform the partial differentiation.
   

   >>> a=S.diff(f,x)
   >>> S.pprint(a)
   -6*y + 2*x + 4*x³
   >>> b=S.diff(f,y)
   >>> S.pprint(b)
   -6*x + 6*y

   SymPy’s S.diff command will differentiate a function with respect to all or some of its variables. In our example, S.diff(f,x) is differentiating f with respect to x and likewise S.diff(f,y) is differentiating f with respect to y. To find the full derivative; however, simply execute S.diff(f,[x,y]). Another very useful command in SymPy is S.pprint or “pretty print”. This command will take SymPy’s symbolic representation of a function and output it in a more readable manner. To see the difference, inspect our new variables a or b without using S.pprint. With the gradient now stored in memory, we can use these variables to solve the system of equations and find the critical points.
6. To do this, we can use SymPy’s built-in equations solver.
   

   >>> S.solve_system([a,b],[x,y])
   [(-1, -1), (0, 0), (1, 1)]

   Our critical points are (-1,-1), (0,0) and (1,1). We now use the second-order conditions to find if they are local maximum/minimum, saddle points or neither. To do so, we must find the Hessian matrix of f and find out the “definiteness” of this matrix for all of our critical points.
7. SymPy has commands that will easily allow us to calculate the Hessian, as well as the principal minors of those matrices, which will tell us the definiteness.
   

   >>> H=S.hessian(f,[x,y])

   The translation from Python matrix form to HTML is ugly, so for clarity, we find the Hessian of f is equal to the equation below.
   
   
   Given that the Hessian has only one element containing a variable (12x²+2), we will only have to plug in the x-values of our critical points to be able to classify them; however, in order to find the definiteness of these matrices we will have to store three matrices in memory for SymPy to operate on.
8. To do so, we use SymPy’s internal matrix package, S.Matrix, to store these matrices.
   

   >>> M1=S.Matrix([14,-6],[-6,6])
   >>> M2=S.Matrix([2,-6],[-6,6])

   Note, since 2+12*(+/-1)² are equivalent, M1 represents both critical points (-1,-1) and (1,1) and M2 represents (0,0). Our final step is to find the principal minors of these matrices to determine their definiteness.
9. SymPy has a built in function that will return all N principal minors of an NxN matrix, so we need only to check the sign of the results of this command to get our answers.
   

   >>> S.Matrix.berkowitz_minors(M1)
   (14, 48)
   >>> S.Matrix.berkowitz_minors(M2)
   (2, -24)

   M1 is positive definite because all principal minors are &ge 0, therefore (-1,-1) and (1,1) are both local minimum. Conversely, because the principal minors of M2 flip sign from positive to negative it is indefinite, therefore (0,0) is a saddle-point. 

To subscribe to the "Guy WhoSteals" feed, click here.
Stolen from: http://www.drewconway.com/zia/?p=274
You can add yourself to the GuyWhoSteals fanpage on Facebook or follow GuyWhoSteals on Twitter.