Showing posts with label math. Show all posts
Showing posts with label math. Show all posts

2015-05-27

Guide to CVXOPT's quadprog() for row-major and/or MATLAB-speaking minds

I am using CVXOPT now. My mind is very row-major, because I mostly program in C and Python. And my mind is also MATLAB-style, because my friend send me a code in MATLAB. Now I need to convert his code to CVXOPT in Python. Here are some thing you wanna pay attention.

1. Matrixes in CVXOPT are column-major


By the code:
G=cvxopt.matrix([[1.,1],[-1,2],[2,1],[-1,0],[0,-1]])
I thought I created a 5-by-2 matrix (5 horizontal rows and 2 vertical columns). However, I got a 2-by-5 matrix.

In [72]: G
Out[72]: <2x5 matrix, tc='d'>

In [73]: print G
[ 1.00e+00 -1.00e+00  2.00e+00 -1.00e+00  0.00e+00]
[ 1.00e+00  2.00e+00  1.00e+00  0.00e+00 -1.00e+00]

2. The matrixes must be of type float.


Make sure your matrixes are of type float, even though your coefficients are integers.

Earlier I created a matrix:
In [53]: q=cvxopt.matrix([[-2],[-6]])

Later I got the error like
TypeError: 'q' must be a 'd' matrix with one column

This is because the matrix is not of type float:
In [54]: q
Out[54]: <1x2 matrix, tc='i'>

To fix, simply make at least one element float, e.g., "1." Then the type code will become d.
In [89]: q
Out[89]: <1x2 matrix, tc='d'>


3. Generating the matrixes needed for your optimization problem.


First of all, the cvxopt formulates quadratic programming so much simpler than MATLAB. Here is how cvxopt formulates the problem:

$$ \begin{array}{ll}
\mbox{minimize} & (1/2) x^TPx + q^T x \\
\mbox{subject to} & G x \preceq h, \\
& Ax = b.
\end{array}
$$

The boundaries for all variables can be defined as part of the matrix G - just use eyes. E.g., the lower boundary constraints

$$x_1 > 1, x_2>2, x_3> 3
$$
is just
$$
\begin{bmatrix}
-1 & 0 & 0 \\
0 & -1 & 0 \\
0 & 0 & -1
\end{bmatrix}
\begin{bmatrix}
x_1\\
x_2 \\
x_ 3
\end{bmatrix}
< \begin{bmatrix} -1\\ -2\\ -3 \end{bmatrix} $$ Note that in either MATLAB's or cvxopt's formulation, the linear terms are always on the side of "smaller than or equal to". Hence, the minus one above. MATLAB formulations have two additional column arrays for upper and lower boundaries: $$ \begin{array}{ll} \mbox{minimize} & (1/2) x^THx + f^T x \\ \mbox{subject to} & A x \preceq b, \\ & A_{eq}x = b_{eq}, \\ & lb \preceq x \preceq ub. \end{array} $$ Apparently, there are lazy people who wants to have a straightforward conversion from MATLAB's formulation to cvxopt's. However, I am not sure whether that translation is out-dated because cvxopt changed their APIs. Here is an easier one.

Let the MATLAB version be
x = quadprog(H,f,A,b,Aeq,beq,lb,ub, ...)
This is how you do it in cvxopt:
import numpy
import cvxopt
import cvxopt.solvers
n = H.shape[1]   # n is the number of variables 

P = H
q = f
G = numpy.vstack([A, -numpy.eye(n), numpy.eye(n)])
h = numpy.vstack([b, -lb, ub])
A = Aeq
b = beq

sol = cvxopt.solvers.qp(cvxopt.matrix(P), cvxopt.matrix(q), cvxopt.matrix(G), cvxopt.matrix(h), cvxopt.matrix(A), cvxopt.matrix(b))
x = sol['x']

If you don't have lb nor/and ub, no need to have it/them in G and h.

Put things together

So, let's see an example. The example in MATLAB Optimization Toolbox can be solved in cvxopt as follows (in iPython shell):

In [81]: G=cvxopt.matrix([[1.,1],[-1,2],[2,1],[-1,0],[0,-1]])

In [82]: P=cvxopt.matrix([[1,-1],[-1,2.]])

In [83]: q=cvxopt.matrix([[-2,-6.]])

In [84]: h=cvxopt.matrix([2,2,3,0,0.])

In [85]: Solv=cvxopt.solvers.qp(P.T, q, G.T, h)
     pcost       dcost       gap    pres   dres
 0: -1.0389e+01 -8.2778e+00  2e+01  9e-01  1e+00
 1: -7.2856e+00 -9.9286e+00  3e+00  1e-16  4e-16
 2: -8.1161e+00 -8.6188e+00  5e-01  8e-17  2e-16
 3: -8.2068e+00 -8.2359e+00  3e-02  6e-17  3e-15
 4: -8.2220e+00 -8.2224e+00  3e-04  5e-17  1e-15
 5: -8.2222e+00 -8.2222e+00  3e-06  8e-17  6e-16
Optimal solution found.

In [86]: print Solv['x']
[ 6.67e-01]
[ 1.33e+00]

2014-06-19

Why SVD can be used to solve EVD (eigenvalue decomposition)

I use SVD (Singular Vector Decomposition) to solve eigenvalue problems (EVD) (note that you don't do so when the matrix is very large) often. But thru years, I gradually forget why we can do so.

I did some revisiting today and found a good explanation. Equation 11 on this webpage (http://mathworld.wolfram.com/EigenDecomposition.html) explains everything clearly. The right-hand side is the form of SVD. Eq. (11) also reveals the "remarkable relationship between a diagonalized matrix, eigenvalues, and eigenvectors follows from the beautiful mathematical identity (the eigen decomposition)" [See also]. And the reason Eq.(11) holds is from Equations 1 to 10.

For more info on the connection between SVD and EVD, you can read this section in Wikipedia:http://en.wikipedia.org/wiki/Singular_value_decomposition#Relation_to_eigenvalue_decomposition

See also: Matrix Diagonalization, Wolfram MathWorld, http://mathworld.wolfram.com/MatrixDiagonalization.html

What a determinant really determines

The determinant is an important concept in linear algebra. Since a determinant is defined for a square matrix, many people would think it as a property of a matrix. But what does it really determine? In other words, why did mathematicians invent it? And why is it defined so? For example, why $\ \det\left ( \left [ \array{ a & b \cr c & d} \right ] \right) = ad - bc $?

Answer: It determines whether a system of linear equations has a unique solution.

The concept of matrix does not come from nowhere. It is strongly related to linear equations. Let's consider a system of linear equations: $$ \begin{align}
ax + by & = C_1\\
cx + dy & = C_2
\label{eq:two_eqs}
\end{align} $$

If $C_1=C_2=0$, eliminating $x$ will result in this: $(ad-bc)y=0$. Pay attention to the coefficient: $ad-bc$. Does it look like $\ \det\left ( \left [ \array{ a & b \cr c & d} \right ] \right) $?

Now if $ad-bc=0$, then $y$ can take whatever value to be a solution of the equations. Namely, the equation has unlimited amount of solutions.

If $C_1 \not = 0$ and $C_2 \not = 0$, the result of eliminating $x$ will have a non-zero constant on the right-hand side: $(ad-bc)y=C_3$. Now if $ad-bc = 0$, there is no way for the equations to have a solution.

Therefore, the determinant actually determines whether a system of linear equations has a unique solution. A system of linear equations can be represented as a matrix. So the determinant of the matrix defines the property of the linear system that the set of equations defines.

References:
1. System of Linear Equations, http://www.math.oregonstate.edu/home/programs/undergrad/CalculusQuestStudyGuides/vcalc/system/system.html
2. Determinant, Wolfram MathWorld, http://mathworld.wolfram.com/Determinant.html
3. Determinant, Wikipedia, http://en.wikipedia.org/wiki/Determinant