+ + + + +
+A mathematical (quadratic) optimization problem, or just optimization problem, has the form +$$ +\begin{align*} + &\mathrm{min}_{\lambda}\hspace{0.2cm} \frac{1}{2}\boldsymbol{\lambda}^T\boldsymbol{P}\boldsymbol{\lambda}+\boldsymbol{q}^T\boldsymbol{\lambda},\\ \nonumber + &\mathrm{subject\hspace{0.1cm}to} \hspace{0.2cm} \boldsymbol{G}\boldsymbol{\lambda} \preceq \boldsymbol{h} \wedge \boldsymbol{A}\boldsymbol{\lambda}=f. +\end{align*} +$$ + +subject to some constraints for say a selected set \( i=1,2,\dots, n \). +In our case we are optimizing with respect to the Lagrangian multipliers \( \lambda_i \), and the +vector \( \boldsymbol{\lambda}=[\lambda_1, \lambda_2,\dots, \lambda_n] \) is the optimization variable we are dealing with. + +
+In our case we are particularly interested in a class of optimization problems called convex optmization problems. +In our discussion on gradient descent methods we discussed at length the definition of a convex function. + +
+Convex optimization problems play a central role in applied mathematics and we recommend strongly Boyd and Vandenberghe's text on the topics. + +
+
+ +
+ + +
+ + + + +
+If we use Python as programming language and wish to venture beyond +scikit-learn, tensorflow and similar software which makes our +lives so much easier, we need to dive into the wonderful world of +quadratic programming. We can, if we wish, solve the minimization +problem using say standard gradient methods or conjugate gradient +methods. However, these methods tend to exhibit a rather slow +converge. So, welcome to the promised land of quadratic programming. + +
+The functions we need are contained in the quadratic programming package CVXOPT and we need to import it together with numpy as + +
+ + +
import numpy
+import cvxopt
++This will make our life much easier. You don't need t write your own optimizer. + +
+
+ +
+ + +
+ + + + +
+We remind ourselves about the general problem we want to solve +$$ +\begin{align*} + &\mathrm{min}_{x}\hspace{0.2cm} \frac{1}{2}\boldsymbol{x}^T\boldsymbol{P}\boldsymbol{x}+\boldsymbol{q}^T\boldsymbol{x},\\ \nonumber + &\mathrm{subject\hspace{0.1cm} to} \hspace{0.2cm} \boldsymbol{G}\boldsymbol{x} \preceq \boldsymbol{h} \wedge \boldsymbol{A}\boldsymbol{x}=f. +\end{align*} +$$ + +
+Let us show how to perform the optmization using a simple case. Assume we want to optimize the following problem +$$ +\begin{align*} + &\mathrm{min}_{x}\hspace{0.2cm} \frac{1}{2}x^2+5x+3y \\ \nonumber + &\mathrm{subject to} \\ \nonumber + &x, y \geq 0 \\ \nonumber + &x+3y \geq 15 \\ \nonumber + &2x+5y \leq 100 \\ \nonumber + &3x+4y \leq 80. \\ \nonumber +\end{align*} +$$ + +The minimization problem can be rewritten in terms of vectors and matrices as (with \( x \) and \( y \) being the unknowns) +$$ +\frac{1}{2}\begin{bmatrix} x\\ y \end{bmatrix}^T \begin{bmatrix} 1 & 0\\ 0 & 0 \end{bmatrix} \begin{bmatrix} x \\ y \end{bmatrix} + \begin{bmatrix}3\\ 4 \end{bmatrix}^T \begin{bmatrix}x \\ y \end{bmatrix}. +$$ + +Similarly, we can now set up the inequalities (we need to change \( \geq \) to \( \leq \) by multiplying with \( -1 \) on bot sides) as the following matrix-vector equation +$$ +\begin{bmatrix} -1 & 0 \\ 0 & -1 \\ -1 & -3 \\ 2 & 5 \\ 3 & 4\end{bmatrix}\begin{bmatrix} x \\ y\end{bmatrix} \preceq \begin{bmatrix}0 \\ 0\\ -15 \\ 100 \\ 80\end{bmatrix}. +$$ + +We have collapsed all the inequalities into a single matrix \( \boldsymbol{G} \). We see also that our matrix +$$ +\boldsymbol{P} =\begin{bmatrix} 1 & 0\\ 0 & 0 \end{bmatrix} +$$ + +is clearly positive semi-definite (all eigenvalues larger or equal zero). +Finally, the vector \( \boldsymbol{h} \) is defined as +$$ +\boldsymbol{h} = \begin{bmatrix}0 \\ 0\\ -15 \\ 100 \\ 80\end{bmatrix}. +$$ + +
+Since we don't have any equalities the matrix \( \boldsymbol{A} \) is set to zero +The following code solves the equations for us +
+ + +
# Import the necessary packages
+import numpy
+from cvxopt import matrix
+from cvxopt import solvers
+P = matrix(numpy.diag([1,0]), tc=’d’)
+q = matrix(numpy.array([3,4]), tc=’d’)
+G = matrix(numpy.array([[-1,0],[0,-1],[-1,-3],[2,5],[3,4]]), tc=’d’)
+h = matrix(numpy.array([0,0,-15,100,80]), tc=’d’)
+# Construct the QP, invoke solver
+sol = solvers.qp(P,q,G,h)
+# Extract optimal value and solution
+sol[’x’]
+sol[’primal objective’]
++
+ +
+ + +
+ + + + +
+We are now ready to return to our setup of the optmization problem for a more realistic case. Introducing the slack parameter \( C \) we have +$$ +\frac{1}{2} \boldsymbol{\lambda}^T\begin{bmatrix} y_1y_1K(\boldsymbol{x}_1,\boldsymbol{x}_1) & y_1y_2K(\boldsymbol{x}_1,\boldsymbol{x}_2) & \dots & \dots & y_1y_nK(\boldsymbol{x}_1,\boldsymbol{x}_n) \\ +y_2y_1K(\boldsymbol{x}_2,\boldsymbol{x}_1) & y_2y_2K(\boldsymbol{x}_2,\boldsymbol{x}_2) & \dots & \dots & y_1y_nK(\boldsymbol{x}_2,\boldsymbol{x}_n) \\ +\dots & \dots & \dots & \dots & \dots \\ +\dots & \dots & \dots & \dots & \dots \\ +y_ny_1K(\boldsymbol{x}_n,\boldsymbol{x}_1) & y_ny_2K(\boldsymbol{x}_n\boldsymbol{x}_2) & \dots & \dots & y_ny_nK(\boldsymbol{x}_n,\boldsymbol{x}_n) \\ +\end{bmatrix}\boldsymbol{\lambda}-\mathbb{I}\boldsymbol{\lambda}, +$$ + +subject to \( \boldsymbol{y}^T\boldsymbol{\lambda}=0 \). Here we defined the vectors \( \boldsymbol{\lambda} =[\lambda_1,\lambda_2,\dots,\lambda_n] \) and +\( \boldsymbol{y}=[y_1,y_2,\dots,y_n] \). +With the slack constants this leads to the additional constraint \( 0\leq \lambda_i \leq C \). + +
+ +
+ +
+ + +