Math 480: Lecture 11
Numpy - The fundamental package for scientific computing with Python
Documentation
- The Guide to NumPy is Travis Oliphant's excellent book about Numpy.
- The Numpy website is filled with useful links to make learning Numpy easy.
Some Quick Teasers
Mandelbrot Example (see http://www.scipy.org/Tentative_NumPy_Tutorial/Mandelbrot_Set_Example). This is how very clever advanced users of Numpy can do amazing things incredibly quickly in just a few lines of code.
{{{id=207|
%python
import numpy
def mandelbrot(h,w, maxiter=20):
y,x = numpy.ogrid[-1.4:1.4:h*1j, -2:0.8:w*1j]
c = x + y*1j
z = c
divtime = maxiter + numpy.zeros(z.shape, dtype=int)
for i in xrange(maxiter):
z = z**2 + c
diverge = z*numpy.conj(z) > 2**2 # who is diverging
div_now = diverge & (divtime==maxiter) # who is diverging now
divtime[div_now] = i # note when
z[diverge] = 2 # avoid diverging too much
return divtime
///
}}}
{{{id=166|
@interact
def _(size=(300,(100..600)), maxiter=(20,(1..40))):
import pylab
time im = mandelbrot(size, size, maxiter)
pylab.imshow(im)
time pylab.savefig('a.png')
///
}}}
{{{id=2|
///
}}}
Image Compression
{{{id=209|
import pylab; import numpy
A_image = numpy.mean(pylab.imread(DATA + 'quad.png'), 2)
u,s,v = numpy.linalg.svd(A_image)
S = numpy.zeros( A_image.shape )
S[:len(s),:len(s)] = numpy.diag(s)
n = A_image.shape[0]
@interact
def svd_image(i = ("Eigenvalues (quality)",(20,(1..A_image.shape[0]//2)))):
A_approx = numpy.dot(numpy.dot(u[:,:i], S[:i,:i]), v[:i,:])
g = graphics_array([matrix_plot(A_approx), matrix_plot(A_image)])
show(g, axes=False, figsize=(10,6))
html("Compressed to %.1f%% of size using %s eigenvalues."%(
100*(2.0*i*n+i)/(n*n), i))
///
}}}
{{{id=208|
///
}}}
{{{id=210|
///
}}}
{{{id=6|
///
}}}
{{{id=12|
///
}}}
Sage Matrices and Numpy Arrays
Sage has a notion of matrices, which model the matrices that come up mainly in pure mathematics (e.g., linear algebra).
{{{id=172|
set_random_seed(2)
A = random_matrix(QQ, 3, 3); A
///
[ -1 -1 -1/2]
[-1/2 2 0]
[ 0 0 0]
}}}
{{{id=171|
A.echelon_form()
///
[ 1 0 2/5]
[ 0 1 1/10]
[ 0 0 0]
}}}
{{{id=173|
A * A.transpose()
///
[ 9/4 -3/2 0]
[-3/2 17/4 0]
[ 0 0 0]
}}}
{{{id=169|
///
}}}
Numpy arrays are vastly different than Sage matrices, to put it mildly. Essentially every rule, constraint, convention, etc., is different. Numpy arrays serve a very different purpose than Sage matrices. They are the ultimate Python general n-dimensional dense array object. In contrast, Sage matrices model very mathematical objects with clear mathematically meaningful semantics and fast algorithms mainly for certain types of algebraic computations. (In some cases Sage matrices are implemented using Numpy under the hood.)
Numpy array have similar functionality to Matlab arrays, given that they have very similar computational domains. See Numpy for Matlab Users for a table that compares and contrasts their functions. That page also has an excellent table showing equivalent commands in Matlab and Numpy.
{{{id=0|
import numpy
B = numpy.array(A); B
///
array([[-1. , -1. , -0.5],
[-0.5, 2. , 0. ],
[ 0. , 0. , 0. ]])
}}}
{{{id=180|
///
}}}
Echelon form?
On Tue, Nov 18, 2008 at 14:21, Robert Young wrote:
> Hi,
>
> Is there a method in NumPy that reduces a matrix to it's reduced row echelon
> form? I'm brand new to both NumPy and linear algebra, and I'm not quite sure
> where to look.
No, we don't have a function to do that. What do you need it for?
While many linear algebra books talk about it, I have never seen a
practical use for it outside of manual matrix computations.
--
Robert Kern
(See http://mail.scipy.org/pipermail/numpy-discussion/2008-November/038704.html.)
In doing linear algebra over the rational numbers (or any field with exact entries) reduced row echelon form is a critically important algorithm. For matrix algebra with floating point numbers (reals, complexes, etc.) reduced row echelon form is basically pointless. Numpy doesn't even have it.
Numpy array multiplication (with a *) is componentwise multiplication of entries:
{{{id=177|
B * B.T
///
array([[ 1. , 0.5, -0. ],
[ 0.5, 4. , 0. ],
[-0. , 0. , 0. ]])
}}}
Use dot for matrix multiplication:
{{{id=178|
numpy.dot(B, B.T)
///
array([[ 2.25, -1.5 , 0. ],
[-1.5 , 4.25, 0. ],
[ 0. , 0. , 0. ]])
}}}
(NOTE: numpy does have numpy.matrix, where * does to matrix multiplication...)
{{{id=185|
C = numpy.matrix(A)
C * C.T
///
matrix([[ 2.25, -1.5 , 0. ],
[-1.5 , 4.25, 0. ],
[ 0. , 0. , 0. ]])
}}}
{{{id=175|
///
}}}
Create a numpy array as follows. For numeric types, one can and should explicitly
give the datatype (dtype), and it can be any numerical type or generic
Python objects. This is already much different than Sage's very mathematical
matrices.
{{{id=25|
Z = numpy.array([[sqrt(2) + 3, random_matrix(QQ, 2)], [1.5, 2/3]])
Z
///
array([[sqrt(2) + 3, [1/2 -1]
[ -1 2]],
[1.50000000000000, 2/3]], dtype=object)
}}}
{{{id=26|
Z + Z
///
array([[2*sqrt(2) + 6, [ 1 -2]
[-2 4]],
[3.00000000000000, 4/3]], dtype=object)
}}}
Symbolics look funny above because Numpy arrays use the "wrong" formatting option when converting them to strings. The following function is a pretty printer for Numpy arrays for the Sage notebook...
{{{id=30|
def shownp(A):
if len(A.shape) > 2:
for i in range(A.shape[0]):
shownp(A[i])
return
if len(A.shape) == 1:
A = A.reshape(A.shape[0], 1)
nrows, ncols = A.shape
s = r'\left(\begin{array}{'
s += 'r'*ncols + '}\n'
for i in range(nrows):
s += ' & '.join([latex(A[i,j]) for j in range(ncols)]) + '\\\\\n'
s += r'\end{array}\right)'
s = s.replace('\\texttt','')
html('$%s$'%s)
///
}}}
{{{id=31|
shownp(Z)
///
\left(\begin{array}{rr} \sqrt{2} + 3 & \left(\begin{array}{rr} \frac{1}{2} & -1 \\ -1 & 2 \end{array}\right)\\ 1.50000000000000 & \frac{2}{3}\\ \end{array}\right)
}}}
{{{id=66|
# A 1-d array
A = numpy.array([1..5])
shownp(A)
///
\left(\begin{array}{r} 1\\ 2\\ 3\\ 4\\ 5\\ \end{array}\right)
}}}
Numpy arrays can have any number of dimensions! Think of these higher dimensional arrays as arrays of arrays (of arrays...). They are useful, e.g., in image processing, where there is one array for each color channel.
{{{id=41|
M = numpy.array([[[1,2], [3,4]],
[[1/3,5], [10,2]], [[8,1], [-5,2]]], dtype=int)
M
///
array([[[ 1, 2],
[ 3, 4]],
[[ 0, 5],
[10, 2]],
[[ 8, 1],
[-5, 2]]])
}}}
{{{id=82|
# The actual array entries are *NOT* Python int's!
type(M[0,0,0])
///
}}}
{{{id=83|
M[0,0,0].item()
///
1
}}}
{{{id=50|
shownp(M)
///
\left(\begin{array}{rr} 1 & 2\\ 3 & 4\\ \end{array}\right)
\left(\begin{array}{rr} 0 & 5\\ 10 & 2\\ \end{array}\right)
\left(\begin{array}{rr} 8 & 1\\ -5 & 2\\ \end{array}\right)
}}}
{{{id=73|
# dtype: A data-type object that fully describes each
# fixed-length item in the array.
# The dtype attribute can be set to anything that can be
# interpreted as a data type (numpy includes 21 of these).
# Setting this attribute allows you to change the interpretation of
# the data in the array. The new data type must be compatible with the array's
# current data type.
///
}}}
{{{id=74|
M.dtype
///
dtype('int64')
}}}
{{{id=75|
M.dtype = float
///
}}}
{{{id=76|
shownp(M)
///
\left(\begin{array}{rr} {4.94065645841e-324} & {9.88131291682e-324}\\ {1.48219693752e-323} & {1.97626258336e-323}\\ \end{array}\right)
\left(\begin{array}{rr} 0.0 & {2.47032822921e-323}\\ {4.94065645841e-323} & {9.88131291682e-324}\\ \end{array}\right)
\left(\begin{array}{rr} {3.95252516673e-323} & {4.94065645841e-324}\\ {nan} & {9.88131291682e-324}\\ \end{array}\right)
}}}
{{{id=77|
M.dtype = numpy.int32
///
}}}
{{{id=78|
shownp(M)
///
\left(\begin{array}{rrrr} 1 & 0 & 2 & 0\\ 3 & 0 & 4 & 0\\ \end{array}\right)
\left(\begin{array}{rrrr} 0 & 0 & 5 & 0\\ 10 & 0 & 2 & 0\\ \end{array}\right)
\left(\begin{array}{rrrr} 8 & 0 & 1 & 0\\ -5 & -1 & 2 & 0\\ \end{array}\right)
}}}
{{{id=84|
shownp(M.astype(float))
///
\left(\begin{array}{rrrr} 1.0 & 0.0 & 2.0 & 0.0\\ 3.0 & 0.0 & 4.0 & 0.0\\ \end{array}\right)
\left(\begin{array}{rrrr} 0.0 & 0.0 & 5.0 & 0.0\\ 10.0 & 0.0 & 2.0 & 0.0\\ \end{array}\right)
\left(\begin{array}{rrrr} 8.0 & 0.0 & 1.0 & 0.0\\ -5.0 & -1.0 & 2.0 & 0.0\\ \end{array}\right)
}}}
{{{id=48|
# flags: (not settable) special array-connected dictionary-like object with
# attributes showing the state of flags in this array;
# only the flags WRITEABLE, ALIGNED, and
# UPDATEIFCOPY can be modified by setting
# attributes of this object
M.flags
///
C_CONTIGUOUS : True
F_CONTIGUOUS : False
OWNDATA : True
WRITEABLE : True
ALIGNED : True
UPDATEIFCOPY : False
}}}
{{{id=43|
# shape: (settable) tuple showing the array shape; setting this
# attribute re-shapes the array
M.shape
///
(3, 2, 4)
}}}
{{{id=67|
M.shape = (3,4,2)
shownp(M)
///
\left(\begin{array}{rr} 1 & 0\\ 2 & 0\\ 3 & 0\\ 4 & 0\\ \end{array}\right)
\left(\begin{array}{rr} 0 & 0\\ 5 & 0\\ 10 & 0\\ 2 & 0\\ \end{array}\right)
\left(\begin{array}{rr} 8 & 0\\ 1 & 0\\ -5 & -1\\ 2 & 0\\ \end{array}\right)
}}}
{{{id=68|
M.shape = (3,2,4); M
///
array([[[ 1, 0, 2, 0],
[ 3, 0, 4, 0]],
[[ 0, 0, 5, 0],
[10, 0, 2, 0]],
[[ 8, 0, 1, 0],
[-5, -1, 2, 0]]], dtype=int32)
}}}
{{{id=71|
# Reshape makes a reshaped copy:
shownp(M.reshape((4,3,2)))
///
\left(\begin{array}{rr} 1 & 0\\ 2 & 0\\ 3 & 0\\ \end{array}\right)
\left(\begin{array}{rr} 4 & 0\\ 0 & 0\\ 5 & 0\\ \end{array}\right)
\left(\begin{array}{rr} 10 & 0\\ 2 & 0\\ 8 & 0\\ \end{array}\right)
\left(\begin{array}{rr} 1 & 0\\ -5 & -1\\ 2 & 0\\ \end{array}\right)
}}}
{{{id=72|
M.shape
///
(3, 2, 4)
}}}
{{{id=44|
# strides: (settable) tuple showing how many bytes must be jumped in
# the data segment to get from one entry to the next
M.strides
///
(32, 16, 4)
}}}
{{{id=51|
# ndim: (not settable) number of dimensions in array
M.ndim
///
3
}}}
{{{id=52|
# data: (settable) buffer object loosely wrapping the array data
# (only works for single-segment arrays)
d = M.data; d
///
}}}
{{{id=54|
d[0:10]
///
'\x01\x00\x00\x00\x00\x00\x00\x00\x02\x00'
}}}
{{{id=55|
type(d)
///
}}}
{{{id=56|
# size: (not settable) total number of elements
M.size
///
24
}}}
{{{id=46|
# itemsize: (settable) one-dimensional, indexable iterator
# object that acts somewhat like a 1-d array
M.flat
///
}}}
{{{id=57|
len(M.flat)
///
24
}}}
{{{id=47|
list(M.flat)
///
[1, 0, 2, 0, 3, 0, 4, 0, 0, 0, 5, 0, 10, 0, 2, 0, 8, 0, 1, 0, -5, -1, 2, 0]
}}}
{{{id=63|
M.flat[0] = 20090283
///
}}}
{{{id=64|
shownp(M)
///
\left(\begin{array}{rrrr} 20090283 & 0 & 2 & 0\\ 3 & 0 & 4 & 0\\ \end{array}\right)
\left(\begin{array}{rrrr} 0 & 0 & 5 & 0\\ 10 & 0 & 2 & 0\\ \end{array}\right)
\left(\begin{array}{rrrr} 8 & 0 & 1 & 0\\ -5 & -1 & 2 & 0\\ \end{array}\right)
}}}
{{{id=61|
# __array_interface__: (not settable) dictionary with keys (data,
# typestr, descr, shape, strides) for compliance with Python side
# of array protocol
M.__array_interface__
///
{'descr': [('', '
}}}
{{{id=65|
# __array_priority__: (not settable) always 0.0 for base type ndarray
M.__array_priority__
///
0.0
}}}
{{{id=42|
///
}}}
A 2d Numpy array:
{{{id=4|
A = numpy.array([[1/3,2,3], [4,5,6], [7,8,9]], dtype=float)
shownp(A)
///
\left(\begin{array}{rrr} 0.333333333333 & 2.0 & 3.0\\ 4.0 & 5.0 & 6.0\\ 7.0 & 8.0 & 9.0\\ \end{array}\right)
}}}
{{{id=8|
type(A[0,0])
///
}}}
{{{id=11|
A[0,0] = 10; shownp(A)
///
\left(\begin{array}{rrr} 10.0 & 2.0 & 3.0\\ 4.0 & 5.0 & 6.0\\ 7.0 & 8.0 & 9.0\\ \end{array}\right)
}}}
{{{id=5|
A.flags
///
C_CONTIGUOUS : True
F_CONTIGUOUS : False
OWNDATA : True
WRITEABLE : True
ALIGNED : True
UPDATEIFCOPY : False
}}}
{{{id=7|
A.flags['WRITEABLE'] = False
///
}}}
{{{id=9|
A.flags
///
C_CONTIGUOUS : True
F_CONTIGUOUS : False
OWNDATA : True
WRITEABLE : False
ALIGNED : True
UPDATEIFCOPY : False
}}}
{{{id=10|
A[0,0] = 10
///
Traceback (most recent call last):
File "", line 1, in
File "_sage_input_324.py", line 10, in
exec compile(u'open("___code___.py","w").write("# -*- coding: utf-8 -*-\\n" + _support_.preparse_worksheet_cell(base64.b64decode("QVswLDBdID0gMTA="),globals())+"\\n"); execfile(os.path.abspath("___code___.py"))' + '\n', '', 'single')
File "", line 1, in
File "/tmp/tmp2vGkEn/___code___.py", line 3, in
exec compile(u'A[_sage_const_0 ,_sage_const_0 ] = _sage_const_10 ' + '\n', '', 'single')
File "", line 1, in
RuntimeError: array is not writeable
}}}
{{{id=13|
# Mutability (writeability) is not strictly enforced once set,
# since you can trivially change it to true, make a write, then
# change the flag back.
before = A.flags['WRITEABLE']
A.flags['WRITEABLE'] = True
A[0,0] = 20
A.flags['WRITEABLE'] = before
shownp(A)
///
\left(\begin{array}{rrr} 20.0 & 2.0 & 3.0\\ 4.0 & 5.0 & 6.0\\ 7.0 & 8.0 & 9.0\\ \end{array}\right)
}}}
{{{id=19|
# This is different than Sage matrix mutability!
B = matrix(RDF, 3, 3, [1..9])
B.set_immutable()
# B *can't* be changed back to being mutable...
B[0,0] = 20
# Sage is much more about enforcing rigor, etc., since a lot of subtle
# mathematical errors in complicated algebraic algorithms can arise
# from the lax (but more user-friendly?) numpy approach.
///
Traceback (most recent call last): # mathematical errors in complicated algebraic algorithms can arise
File "", line 1, in
File "/tmp/tmpk6hZHR/___code___.py", line 6, in
B[_sage_const_0 ,_sage_const_0 ] = _sage_const_20
File "matrix0.pyx", line 1281, in sage.matrix.matrix0.Matrix.__setitem__ (sage/matrix/matrix0.c:5702)
File "matrix0.pyx", line 404, in sage.matrix.matrix0.Matrix.check_mutability (sage/matrix/matrix0.c:3750)
ValueError: matrix is immutable; please change a copy instead (i.e., use copy(M) to change a copy of M).
}}}
{{{id=17|
# numpy arrays can be serialized and saved to disk
save(M, '/tmp/M.sobj')
print M.dumps()
(M == loads(M.dumps())).all()
///
cnumpy.core.multiarray
_reconstruct
qcnumpy
ndarray
qK�UbRq(KKKKcnumpy
dtype
qUi4K�KRq(KU\left(\begin{array}{rrrr} 20090283 & 0 & 2 & 0\\ 3 & 0 & 4 & 0\\ \end{array}\right)