In [1]:
# Preamble for rendering into this notebook.
from util import style
style()
# Get `numpy` loaded up so we can begin.
import numpy as np
In [2]:
# "Normal" array creation
print np.array([0,1,2,3,4])
In [3]:
# `python` auto ranged list creation
print range(5) # end only
print range(5, 10) # start and end
print range(5, 20, 3) # start, end, and step
In [4]:
# `numpy` equivalent for creating `ndarray`s.
# (From here on out we'll refer to `numpy`'s `ndarray`s simply as "array"s.)
print np.arange(5) # end only
print np.arange(5, 10) # start and end
print np.arange(5, 20, 3) # start, end, and step
In [5]:
# Basic arrays have a "flat" shape.
x = np.arange(100)
print x
print
print 'x shape:', x.shape
In [6]:
# We can reshape them to be 2d (i.e. a 10 x 10 matrix)
y = x.reshape(10, 10)
print y
print
print 'y shape:', y.shape
In [7]:
# ... or we can even make them tensors.
# (Note the shape must match the data size.)
z = x.reshape(5, 4, 5)
print z
print
print 'z shape:', z.shape
numpy also has a special "matrix" type that is a thin wrapper around an array. It's guaranteed to be 2d. It changes just a couple methods so matrix multiplication has a nicer syntax. For more info on the differences, see: this question on Stack Overflow. I'll use them throughout just because we'll be multiplying a lot of matrices.
In [8]:
# I'll use capital variable names for matrices. Here we'll turn a 1D `numpy`
# array into a 2D `numpy` matrix.
# Reminder: this is our array `x`.
print x
print
print 'x shape:', x.shape
print
print 'x type:', type(x)
# We can turn it into a matrix `M`.
M = np.matrix(x)
print 'M type:', type(M)
print
print M
print
print 'M shape:', M.shape
In [9]:
# Note that M by default is just a row vector form of x (that is, all
# of x is put in a single row).
#
# We can make it a different shape.
R = M.reshape((10, 10))
print R
print 'R shape:', R.shape
# We can also create a `numpy` matrix from a 2D `numpy` array directly and
# it will maintain its shape.
print
print 'y shape:', y.shape
N = np.matrix(y)
print N
print 'N shape:', N.shape
In [10]:
# Just to be sure this doesn't work, let's try making a matrix from a
# tensor. (Uncomment the following line to see the error).
# T = np.matrix(z)
In [11]:
# This function displays a matrix as a heatmap of relative values.
# Check out util.py for its source.
from util import display_matrix
In [12]:
# Let's play with some small matrices and visualize them!
A = np.matrix('2 3 ; 2 2')
print 'A'
print A
display_matrix(A, base_size=200)
In [13]:
B = A.I
print 'B = A^-1'
print B
display_matrix(B, base_size=200)
In [14]:
print 'BA = AB = A^-1A = AA^-1 = I'
R = B*A
print R
display_matrix(R, base_size=200)
In [15]:
# Here's a zero matrix
Z = np.asmatrix(np.zeros((4,4)))
print 'Z (0 matrix)'
print Z
display_matrix(Z, base_size=300)
In [16]:
# Let's make a slightly larger matrix to visualize.
M = np.matrix(np.arange(15).reshape(5, 3))
print M
display_matrix(M)
In [17]:
# ...and here's its transpose
print M.T
display_matrix(M.T)
In [18]:
# OK, what about our bigger (10 x 10) matrix, N?
display_matrix(N)
In [19]:
# ...and its transpose
display_matrix(N.T)
In [ ]:
Here—during bootcamp—we got to talking about this matrix (N) and why it doesn't have an invserse. It is singular. You can verify that you can create columns as linear combinations of others (which we'll do right here).
In [20]:
# We can create, e.g., column 4 as a linear cobination of
# columns 2, 3, and 5 (note that we're using 0-based column
# indexing here).
c4 = N[:][2] - N[:][3] + N[:][5]
print c4
print np.allclose(c4, N[:][4])
In [21]:
# N does, however, have a pseudoinverse, though. Let's compute and
# print it.
N_plus = np.linalg.pinv(N)
print N_plus
In [22]:
# That's nice... but what does it look like?
display_matrix(N_plus)
As a side note, the whole reason I decided to try to visualize matrices as heatmaps was to check out what an inverse or pseduoinverse looks like. If you're inspired, Wikipedia has lots more about these subjects; e.g., here's the pseudoinverse page.
In [23]:
# Now that we've got a pseudoinverse, let's see what it looks like
# when multiplied with the original.
print 'Left-multiply'
display_matrix(N_plus * N)
In [24]:
# Neat, it's kind of like it's rotated 90 degrees and scaled.
# (Exercise for the reader: figure out what's going on here.)
# What about right-multiplying?
print 'Right-multiply:'
display_matrix(N * N_plus)
In [25]:
# Those looked the same. Are the left- and right-multiplications of
# the inverses with the original matrix the same?
print np.allclose(N*N_plus, N_plus*N)
In [26]:
# Pop quiz: what is this?
display_matrix(N * N_plus * N)
In [27]:
# That's right---it's the original!
print np.allclose(N * N_plus * N, N)
In [28]:
# A noise matrix!
#
# This `numpy` function takes the
# - mean (mu)
# - standard deviation (sigma)
# - and number of points to generate (n)
#
# ... and gives back a `numpy` array of length n with values
# drawn from Normal(mu, sigma).
noise = np.random.normal(0.0, 10, 30) # noise
print noise
In [29]:
# What do you expect it to look like?
display_matrix(np.matrix(noise))
In [30]:
# All right, let's return to our good old matrix M.
display_matrix(M)
In [31]:
# What happens if we scale it? Guess how you think the visualization
# will change.
display_matrix(M * 5)
In [32]:
# Very good---it doens't change at all!
In [33]:
# What about eigenvalue decompositions? Here's a simple matrix.
A = np.diag((1,2,3))
display_matrix(A)
In [34]:
# This is the eigenvalue decomposition.
# Pop quiz: Can you guess what the eigenvalues and eigenvectors will be?
w, V = np.linalg.eig(A)
In [35]:
# The eigenvalues are scalars; w is a vector of all of them.
display_matrix(np.asmatrix(w))
In [36]:
# The eigenvectors are vectors (well, duh); V is a matrix of
# them (one in each row)
display_matrix(V)
In [37]:
# How about an SVD (singular value decomposition)? Let's create a random
# matrix first. We'll use 9 rows and 6 columns. (Note that this example
# is taken right from `numpy`'s docs on SVD.)
B = np.matrix(np.random.randn(9, 6))
display_matrix(B)
In [38]:
# The SVD gives us a matrix, vector, and matrix.
# We'll perform it and verify their shapes here.
U, s, V = np.linalg.svd(B, full_matrices=True)
print U.shape, s.shape, V.shape
In [39]:
# Here's the first matrix (9 x 9).
display_matrix(U)
In [40]:
# Here's the length-6 vector. Note that the values are of decreasing
# magnitude--this is done by `numpy` by design.
display_matrix(np.asmatrix(s))
In [41]:
# ... and here's the final matrix, a 6 x 6.
display_matrix(V)
In [42]:
# Cool. Let's try to reconstruct B (our original matrix) from the SVD.
# Remember that B is 9 x 6. We'll make a new matrix S with the correct
# dimensions and fill it with zeros to start.
S = np.matrix(np.zeros((9,6)))
display_matrix(S)
In [43]:
# Now, we fill the main diagonal with s.
S[:6, :6] = np.diag(s)
display_matrix(S)
In [44]:
# To reconstruct B, we do:
#
# U x S x V
# 9 x 9 9 x 6 6 x 6
#
# .. which you can verify will result in a 9 x 6
# matrix.
#
# We'll reconstruct into a new matrix 'R' (for
# "Reconstruction").
R = U*S*V
display_matrix(R)
In [45]:
# Does that look like B?
display_matrix(B)
In [46]:
# Yep! We can check if each of the elements are close:
print np.isclose(B, R)
In [47]:
# ... or use a shortcut to verify that all of the elements are close:
print np.allclose(B, R)
In [48]:
# As a little exercise, remember our matrix N from above, which is
# a 10 x 10 with values 0..99?
display_matrix(N)
In [49]:
# We decided it (N) is singular. We might wonder, since its columns
# aren't linearly independent, what in fact *is* its rank?
#
# We can find out with the SVD.
U, s, V = np.linalg.svd(N, full_matrices=True)
display_matrix(np.asmatrix(s), base_size=900)
In [50]:
# It looks like there are only two nonzero values. We can verify this with
# `numpy`'s builtin `matrix_rank` function---which also uses the SVD!
print np.linalg.matrix_rank(N)