DM554/DM545 – Linear and Integer Programming

Ipython

To install ipython under Macosx: install via pip

pip install "ipython[notebook]"

To install under Linux:

sudo apt-get install ipython ipython-qtconsole ipython-notebook

To install under Windows follow these instructions.

Ipython can be used from the shell by calling:

ipython

or with a console that allows the inlining of plots and text highlighting by typing:

ipython qtconsole

or, finally, as notebook from a borwser started by typing

ipython notebook

You can read documentation on ipython at these links:

Plotting

Plotting in python works with the Matplotlib plotting library. Matplotlib integrates with IPython's display system and event loop handling (there is also a feature to load Matplotlib examples %load URL_to_py)

To make plots using Matplotlib, you must first enable IPython's matplotlib mode. To do this, run the %matplotlib command to enable plotting in the current Notebook. This takes an optional argument that specifies which Matplotlib backend should be used. Most of the time, in the Notebook, you will want to use the inline backend, which will embed plots inside the Notebook:

In [1]: 
%matplotlib inline

Functions with % are extra functions of Ipython that add functionalities to the environment. They are called magic functions. Hence the command above does not load the library in python. As usual we need:

In [2]: 
import matplotlib.pyplot as plt
import numpy as np

Let's plot the following polynomial of degree 3:

\[P_3(x)=x^3-7x+6=(x-1)(x-2)(x+3)\]

We use the numpy.poly1d function to represent the polynomial. It takes an array of coefficients a of length \(n+1\) such that the polynomial is described by:

    a[0] * x**n + a[1] * x**(n-1) + ... + a[n-1]*x + a[n]

(Try numpy.polyfit? to know more).

In [3]: 
a=[1,0,-7,6]
P=np.poly1d(a)
print P

   3
1 x - 7 x + 6
In [4]: 
x = np.linspace(-3.5, 3.5, 500)
plt.plot(x, P(x), '-')
plt.axhline(y=0)
plt.title('A polynomial of order 3');

Roots of a polynomial of degree n

Let's \(P_n(x)\) be a general polynomial of degree \(n\):

\[P_n(x)=p_n+p_{n-1}x+p_{n-2}x^2+\cdots+p_0x^n\]

numpy.roots(p) Return the roots of a polynomial with coefficients given in P. For the polynomial of 3rd degree seen above:

In [5]: 
np.roots(P)
Out[5]:
array([-3.,  2.,  1.])

Arrays and Matrices

In python we shall use np.arrays to represent matrices and vectors. Althoguh very similar to lists they differ in that they store only the same numeric type and there is no append method.

In [6]: 
b=np.array([1,2,3])
print 2*b
print b/2
#np.norm(b)

[2 4 6]
[0 1 1]

Indexing and slices

Indices start at 0.

In [7]: 
A=np.array([[1.,2,3],
            [4,5,6],
            [7,8,9]])
print "A = ", A
print A[0,0]
print A[1:] # all rows starting from the second
print A[1] # the second row of the matrix
print A[:,2] # a vector filled by the 3nd column of the matrix 
print "Type: ", A.dtype
print "Dimension:", A.shape, A.ndim

A =  [[ 1.  2.  3.]
 [ 4.  5.  6.]
 [ 7.  8.  9.]]
1.0
[[ 4.  5.  6.]
 [ 7.  8.  9.]]
[ 4.  5.  6.]
[ 3.  6.  9.]
Type:  float64
Dimension: (3, 3) 2
In [8]: 
A.reshape(1,9)
Out[8]:
array([[ 1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9.]])

Caution: reshape is just a view on the array, it does not copy the data. Hence:

In [9]: 
M=A.reshape(1,9)
M[0,7]=10 # now A[2,1] is 10
A
Out[9]:
array([[  1.,   2.,   3.],
       [  4.,   5.,   6.],
       [  7.,  10.,   9.]])
In [10]: 
print A.base
print M.base

None
[[  1.   2.   3.]
 [  4.   5.   6.]
 [  7.  10.   9.]]

Caution with casting:

In [11]: 
B=A.astype('int')
A[0,0]=1.5
B[0,1]=2.5
A, B
Out[11]:
(array([[  1.5,   2. ,   3. ],
       [  4. ,   5. ,   6. ],
       [  7. ,  10. ,   9. ]]),
 array([[ 1,  2,  3],
       [ 4,  5,  6],
       [ 7, 10,  9]]))

Functions to construct arrays

The following is a list of functions to construct matrices and vectors.

In [12]: 
np.zeros((3,4))
Out[12]:
array([[ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.]])
In [13]: 
np.ones((3,4))
np.identity(4)
np.diag([1,2,3,4])
np.random.rand(3,4)
np.random.randint(1,10,size=(3,4))
Out[13]:
array([[5, 4, 6, 7],
       [1, 1, 5, 7],
       [4, 4, 5, 7]])
In [14]: 
np.arange(5)
Out[14]:
array([0, 1, 2, 3, 4])
In [15]: 
np.linspace(0,5,11)
Out[15]:
array([ 0. ,  0.5,  1. ,  1.5,  2. ,  2.5,  3. ,  3.5,  4. ,  4.5,  5. ])

To concatenate arrays you have to use the concatenate function

In [16]: 
np.concatenate([A,A],axis=1) # same as np.column_stack()
Out[16]:
array([[  1.5,   2. ,   3. ,   1.5,   2. ,   3. ],
       [  4. ,   5. ,   6. ,   4. ,   5. ,   6. ],
       [  7. ,  10. ,   9. ,   7. ,  10. ,   9. ]])
In [17]: 
np.concatenate([A,A],axis=0) # same as np.hstack()
Out[17]:
array([[  1.5,   2. ,   3. ],
       [  4. ,   5. ,   6. ],
       [  7. ,  10. ,   9. ],
       [  1.5,   2. ,   3. ],
       [  4. ,   5. ,   6. ],
       [  7. ,  10. ,   9. ]])

Operations with arrays

The operations +,*,/,- are elementwise. For the operations of scalar multiplication and matrix multiplication the dot method must be used. Compare the following two cells:

In [18]: 
b=np.array([1,2,3])
v=np.array([2,0,1])
b*v
Out[18]:
array([2, 0, 3])
In [19]: 
np.dot(b,v)
Out[19]:
5
In [20]: 
A*b
Out[20]:
array([[  1.5,   4. ,   9. ],
       [  4. ,  10. ,  18. ],
       [  7. ,  20. ,  27. ]])
In [21]: 
np.dot(A,b)
Out[21]:
array([ 14.5,  32. ,  54. ])
In [22]: 
A.T # Matrix transpose
Out[22]:
array([[  1.5,   4. ,   7. ],
       [  2. ,   5. ,  10. ],
       [  3. ,   6. ,   9. ]])

Beware that like reshape, transpose does not copy the data. Vectors have one dimension hence transposing does not do anything particular. In this case it is more appropriate to use reshape:

In [23]: 
b.reshape(-1,1) # yelds a column vector, -1 tells to guess the size
Out[23]:
array([[1],
       [2],
       [3]])

Most mathematical functions act element-wise. They are also called universal functions.

In [24]: 
np.cos(np.array([1,2,3]))
Out[24]:
array([ 0.54030231, -0.41614684, -0.9899925 ])
In [25]: 
2**np.arange(4) # arange includes the 0
Out[25]:
array([1, 2, 4, 8])
In [26]: 
np.arange(4)**np.arange(4)
Out[26]:
array([ 1,  1,  4, 27])

Some functions do not act component-wise, eg, min, max, sum. These functions may operate on the whole matrix or column-wise or row-wise.

In [27]: 
print A
np.sum(A) # on the whole matrix
np.sum(A,axis=1) # row-wise
np.sum(A,axis=0) # column-wise (this is also the result of `sum` from base python)

[[  1.5   2.    3. ]
 [  4.    5.    6. ]
 [  7.   10.    9. ]]
Out[27]:
array([ 12.5,  17. ,  18. ])

If it is not, it is possible to make a function element-wise via vectorize. This yields often an improvement in performance with respect to for loops. Compare for example:

In [28]: 
def f(x):
    return 0 if x<=5 else 1
# f(A) # error
np.vectorize(f)(A)
Out[28]:
array([[0, 0, 0],
       [0, 0, 1],
       [1, 1, 1]])

Matrix inverse

The inverse of a matrix can be calculated with a function from the set of linear algebra functions of numpy.

In [29]: 
np.linalg.det(A) # returns the determinant of A, $det(A) \neq$ if matrix invertible 
np.linalg.inv(A)
Out[29]:
array([[-3.33333333,  2.66666667, -0.66666667],
       [ 1.33333333, -1.66666667,  0.66666667],
       [ 1.11111111, -0.22222222, -0.11111111]])

Plotting vectors

To plot vectors we can use matplotlib.pyplot.quiver. In the basic form quiver(X,Y,U,V), X and Y are the coordinates of the tail of the vector and U and V the components of the vector. The plot needs some tweak for the scales.

In [30]: 
plt.quiver(0,0,3,4,scale=1,scale_units='xy',angles='xy')
plt.gca().set_xlim([0,5])
plt.gca().set_ylim([0,5]);
In [31]: 
soa =np.array( [ [3,2], [1,1], [9,9]]) 
U,V = zip(*soa)
plt.quiver([0,0,0],[0,0,0],U,V,angles='xy',scale_units='xy',scale=1)
ax=plt.gca()
ax.set_xlim([-1,10])
ax.set_ylim([-1,10]);

Solving a Linear System

If \(A\) is a mtrix and \(b\) is a vector, the system of linear equations \[Ax=b\] is solved using the solve method from numpy.linalg or scipy.linalg (the second being more efficient):

In [32]: 
import scipy.linalg as sl
x=sl.solve(A, b)
x
Out[32]:
array([ -1.11022302e-16,   7.77156117e-17,   3.33333333e-01])
In [33]: 
np.allclose(np.dot(A, x), b)
Out[33]:
True

Easy printing for LaTex

In [34]: 
def bmatrix(a):
    """Returns a LaTeX bmatrix
    :a: numpy array
    :returns: LaTeX bmatrix as a string
    """
    if len(a.shape) > 2:
        raise ValueError('bmatrix can at most display two dimensions')
    lines = str(a).replace('[', '').replace(']', '').splitlines()
    rv = [r'\begin{bmatrix}']
    rv += [r'  ' + ' & '.join(l.split()) + '\\\\' for l in lines]
    rv +=  [r'\end{bmatrix}']
    return '\n'.join(rv)
In [35]: 
A=np.random.randint(4,size=(2,2))
print bmatrix(A), bmatrix(A.T)

\begin{bmatrix}
  3 & 1\\
  0 & 2\\
\end{bmatrix} \begin{bmatrix}
  3 & 0\\
  1 & 2\\
\end{bmatrix}