DM554/DM545 – Linear and Integer Programming

Using IPython for interactive work

Preamble for working with Arrays and Fractions

See this tutorial Tutorial for an introduction to IPython. Here we show how to use IPython for doing, for example, matrix operations with fractions. Copy this preamble in your script to load the needed modules and implement a function printm and verbatim for pretty printing.

In [1]: 
from numpy import *
from fractions import Fraction as f
set_printoptions(precision=3,suppress=True)

def printm(a):
    """Prints the array as strings 
    :a: numpy array
    :returns: prints the array
    """
    def p(x):
        return str(x)
    p = vectorize(p,otypes=[str])
    print p(a)
    
def tableau(a,W=7):
    """Returns a string for verbatim printing
    :a: numpy array
    :returns: a string
    """
    if len(a.shape) != 2:
        raise ValueError('verbatim displays two dimensions')
    rv = []
    rv+=[r'|'+'+'.join('{:-^{width}}'.format('',width=W) for i in range(a.shape[1]))+"+"]
    rv+=[r'|'+'|'.join(map(lambda i: '{0:>{width}}'.format("x"+str(i+1)+" ",width=W), range(a.shape[1]-2)) )+"|"+
         '{0:>{width}}'.format("-z ",width=W)+"|"
         '{0:>{width}}'.format("b ",width=W)+"|"]
    rv+=[r'|'+'+'.join('{:-^{width}}'.format('',width=W) for i in range(a.shape[1]))+"+"]
    for i in range(a.shape[0]-1):
        rv += [r'| '+' | '.join(['{0:>{width}}'.format(str(a[i,j]),width=W-2) for j in range(a.shape[1])])+" |"]
    rv+=[r'|'+'+'.join('{:-^{width}}'.format('',width=W) for i in range(a.shape[1]))+"+"]
    i = a.shape[0]-1
    rv += [r'| '+' | '.join(['{0:>{width}}'.format(str(a[i,j]),width=W-2) for j in range(a.shape[1])])+" |"]
    rv+=[r'|'+'+'.join('{:-^{width}}'.format('',width=W) for i in range(a.shape[1]))+"+"]
    print '\n'.join(rv)

Elementary row operations

To do elementary row operations on a matrix one can do as follows:

In [2]: 
A = array( [[1, f(9,4), 2, f(4,5)],
             [f(3,4), 1, f(5,4), 0]])
A[0,:] = f(4,9) * A[0,:]
printm(A)

[['4/9' '1' '8/9' '16/45']
 ['3/4' '1' '5/4' '0']]
In [3]: 
tableau(A)

|-------+-------+-------+-------+
|    x1 |    x2 |    -z |     b |
|-------+-------+-------+-------+
|   4/9 |     1 |   8/9 | 16/45 |
|-------+-------+-------+-------+
|   3/4 |     1 |   5/4 |     0 |
|-------+-------+-------+-------+

Note that indices start from 0, hence a vector of size \(n\) has indices from \(0\) to \(n-1\).

Matrix Operations

Construct an array of zeros

In [4]: 
zeros((3,1))
Out[4]:
array([[ 0.],
       [ 0.],
       [ 0.]])

Concatenate an identity matrix to a matrix

In [5]: 
A_1 = concatenate([A,identity(2)],axis=1)
printm(A_1)

[['4/9' '1' '8/9' '16/45' '1.0' '0.0']
 ['3/4' '1' '5/4' '0' '0.0' '1.0']]

Insert a row in between other two:

In [6]: 
A_2 = concatenate([ A[0:1,:], array( [[0,2,1,0]] ), A[1:2,:] ],axis=0)
printm(A_2)

[['4/9' '1' '8/9' '16/45']
 ['0' '2' '1' '0']
 ['3/4' '1' '5/4' '0']]

Ask for the dimensions of an array:

In [7]: 
A.shape
Out[7]:
(2, 4)

Multiply matrices

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 [8]: 
A = ones((4,2)) # a matrix of 1s of size 4x2
b = array([1,2])
A*b
Out[8]:
array([[ 1.,  2.],
       [ 1.,  2.],
       [ 1.,  2.],
       [ 1.,  2.]])

which is element-wise and hence not matrix-vector multiplication. The correct way to do matrix-vector and matrix-matrix multiplication is:

In [9]: 
dot(A,b)
Out[9]:
array([ 3.,  3.,  3.,  3.])

Inverse

In [10]: 
A = random.randint(1,10,size=(4,4))
A_i = linalg.inv(A)
around( dot(A,A_i), decimals=2) # let's verify
Out[10]:
array([[ 1., -0.,  0.,  0.],
       [ 0.,  1.,  0., -0.],
       [ 0., -0.,  1., -0.],
       [ 0., -0.,  0.,  1.]])

Transpose

In [11]: 
A.T
Out[11]:
array([[8, 5, 4, 5],
       [8, 6, 5, 5],
       [7, 3, 3, 3],
       [6, 4, 7, 4]])

Rank

In [12]: 
linalg.matrix_rank(A)
Out[12]:
4

Reduced row echelon form

To find the reduced row echelon form we need to import the module for symbolic computation sympy.

In [13]: 
A=array([[3,1,1,1],
           [2, 4, 1, 1],
           [2,1,2,2]])
b=array([ 2, 2, 1])
import sympy as sy
# np.linalg.solve(A,b)
AA=column_stack([A,b])
sy.Matrix(AA).rref()
Out[13]:
(Matrix([
 [1, 0, 0, 0,  9/13],
 [0, 1, 0, 0,  3/13],
 [0, 0, 1, 1, -4/13]]), [0, 1, 2])