Question

1 . Find the row, column and overall means for the following matrix:

m = np.arange(12).reshape((3,4))

# YOUR CODE HERE

m = np.arange(12).reshape((3,4))
print m
print

print "OVerall", m.mean()
print "Row", m.mean(1)
print "Columne", m.mean(0)

[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]

OVerall 5.5
Row [ 1.5  5.5  9.5]
Columne [ 4.  5.  6.  7.]

2 . Find the outer product of the following two vecotrs

u = np.array([1,3,5,7])
v = np.array([2,4,6,8])

Do this in the following ways:

• Using the function outer in numpy
• Using a nested for loop or list comprehension
• Using numpy broadcasting operatoins

# YOUR CODE HERE

u = np.array([1,3,5,7])
v = np.array([2,4,6,8])

print np.outer(u, v)
print
print np.array([[u_ * v_ for v_ in v] for u_ in u])
print
print u[:,None] * v[None,:]

[[ 2  4  6  8]
 [ 6 12 18 24]
 [10 20 30 40]
 [14 28 42 56]]

[[ 2  4  6  8]
 [ 6 12 18 24]
 [10 20 30 40]
 [14 28 42 56]]

[[ 2  4  6  8]
 [ 6 12 18 24]
 [10 20 30 40]
 [14 28 42 56]]

3 . Create a 10 by 6 matrix of random uniform numbers. Set all rows with any entry less than 0.1 to be zero. For example, here is a 4 by 10 version:

array([[ 0.49722235,  0.88833973,  0.07289358,  0.12375223,  0.39659254,
         0.70267114],
       [ 0.3954172 ,  0.889077  ,  0.71286225,  0.06353112,  0.68107965,
         0.17186995],
       [ 0.74821206,  0.92692111,  0.24871227,  0.26904958,  0.80410194,
         0.22304055],
       [ 0.22582605,  0.37671244,  0.96510957,  0.88819053,  0.14654176,
         0.33987323]])

becomes

array([[ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ],
       [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ,
         0.        ],
       [ 0.74821206,  0.92692111,  0.24871227,  0.26904958,  0.80410194,
         0.22304055],
       [ 0.22582605,  0.37671244,  0.96510957,  0.88819053,  0.14654176,
         0.33987323]])

Hint: Use the following numpy functions - np.random.random , np.any as well as Boolean indexing and the axis argument.

# YOUR CODE HERE

xs = np.random.random((10,6))
print xs
print
xs[(xs < 0.1).any(axis=1), :] = 0
print xs

[[ 0.5117  0.9098  0.2184  0.3631  0.855   0.7114]
 [ 0.3929  0.2313  0.3802  0.5492  0.5567  0.0041]
 [ 0.638   0.0576  0.043   0.8751  0.2926  0.7628]
 [ 0.3679  0.8735  0.0294  0.552   0.2402  0.8848]
 [ 0.4602  0.1932  0.2937  0.8179  0.5595  0.6779]
 [ 0.8091  0.8686  0.418   0.0589  0.4785  0.5212]
 [ 0.5806  0.3092  0.9199  0.6553  0.3492  0.5411]
 [ 0.4491  0.2823  0.2959  0.5635  0.7152  0.5176]
 [ 0.352   0.6328  0.8731  0.1679  0.9875  0.3494]
 [ 0.8262  0.0655  0.0054  0.8869  0.9113  0.1994]]

[[ 0.5117  0.9098  0.2184  0.3631  0.855   0.7114]
 [ 0.      0.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.    ]
 [ 0.      0.      0.      0.      0.      0.    ]
 [ 0.4602  0.1932  0.2937  0.8179  0.5595  0.6779]
 [ 0.      0.      0.      0.      0.      0.    ]
 [ 0.5806  0.3092  0.9199  0.6553  0.3492  0.5411]
 [ 0.4491  0.2823  0.2959  0.5635  0.7152  0.5176]
 [ 0.352   0.6328  0.8731  0.1679  0.9875  0.3494]
 [ 0.      0.      0.      0.      0.      0.    ]]

4 . Use np.linspace to create an array of 100 numbers between 0 and \(2\pi\) (includsive).

• Extract every 10th element using slice notation
• Reverse the array using slice notation
• Extract elements where the absolute difference between the sine and cosine functions evaluated at that element is less than 0.1
• Make a plot showing the sin and cos functions and indicate where they are close

# YOUR CODE HERE

xs = np.linspace(0, 2*np.pi, 100)
print xs[::10]
print
print xs[::-1]
print
idx = np.abs(np.sin(xs)-np.cos(xs)) < 0.1
print xs[idx]
print
plt.scatter(xs[idx], np.sin(xs[idx]))
plt.plot(xs, np.sin(xs), xs, np.cos(xs));

[ 0.      0.6347  1.2693  1.904   2.5387  3.1733  3.808   4.4427  5.0773
  5.712 ]

[ 6.2832  6.2197  6.1563  6.0928  6.0293  5.9659  5.9024  5.8389  5.7755
  5.712   5.6485  5.5851  5.5216  5.4581  5.3947  5.3312  5.2677  5.2043
  5.1408  5.0773  5.0139  4.9504  4.8869  4.8235  4.76    4.6965  4.6331
  4.5696  4.5061  4.4427  4.3792  4.3157  4.2523  4.1888  4.1253  4.0619
  3.9984  3.9349  3.8715  3.808   3.7445  3.6811  3.6176  3.5541  3.4907
  3.4272  3.3637  3.3003  3.2368  3.1733  3.1099  3.0464  2.9829  2.9195
  2.856   2.7925  2.7291  2.6656  2.6021  2.5387  2.4752  2.4117  2.3483
  2.2848  2.2213  2.1579  2.0944  2.0309  1.9675  1.904   1.8405  1.7771
  1.7136  1.6501  1.5867  1.5232  1.4597  1.3963  1.3328  1.2693  1.2059
  1.1424  1.0789  1.0155  0.952   0.8885  0.8251  0.7616  0.6981  0.6347
  0.5712  0.5077  0.4443  0.3808  0.3173  0.2539  0.1904  0.1269  0.0635
  0.    ]

[ 0.7616  0.8251  3.8715  3.9349]

5 . Create a matrix that shows the 10 by 10 multiplication table.

• Find the trace of the matrix
• Extract the anto-diagonal (this should be array([10, 18, 24, 28, 30, 30, 28, 24, 18, 10]) )
• Extract the diagnoal offset by 1 upwards (this should be array([ 2, 6, 12, 20, 30, 42, 56, 72, 90]) )

# YOUR CODE HERE

ns = np.arange(1, 11)
m = ns[:, None] * ns[None, :]
print m
print
print m.trace()
print
print np.flipud(m).diagonal()
print
print m.diagonal(offset=1)

[[  1   2   3   4   5   6   7   8   9  10]
 [  2   4   6   8  10  12  14  16  18  20]
 [  3   6   9  12  15  18  21  24  27  30]
 [  4   8  12  16  20  24  28  32  36  40]
 [  5  10  15  20  25  30  35  40  45  50]
 [  6  12  18  24  30  36  42  48  54  60]
 [  7  14  21  28  35  42  49  56  63  70]
 [  8  16  24  32  40  48  56  64  72  80]
 [  9  18  27  36  45  54  63  72  81  90]
 [ 10  20  30  40  50  60  70  80  90 100]]

385

[10 18 24 28 30 30 28 24 18 10]

[ 2  6 12 20 30 42 56 72 90]

6 . Diagonalize the follwoing matrix

A = np.array([
    [1,  2, 1],
    [6, -1, 0],
    [-1,-2,-1]
])

In other words, find the invertible matrix \(P\) and the diagonal matrix \(D\) such that $A = PDP^{-1} $. Confirm by calculating the value of $PDP^{-1} $.

• Do this mnaully
• Then use numpy.linalg functions to do the same

# YOUR CODE HERE

A = np.array([
    [1,  2, 1],
    [6, -1, 0],
    [-1,-2,-1]
])

dotm = lambda *args: reduce(np.dot, args)

u, V = la.eig(A)
P = V
D = np.diag(u)
print P
print
print np.real_if_close(np.round(u))
print

np.real_if_close(np.round(dotm(P, D, la.inv(P)), 6))

[[ 0.4082 -0.4851 -0.0697]
 [-0.8165 -0.7276 -0.418 ]
 [-0.4082  0.4851  0.9058]]

[-4.  3.  0.]

array([[ 1.,  2.,  1.],
       [ 6., -1., -0.],
       [-1., -2., -1.]])

The code below just intorduces some of the symbolic algebra capabilities of Python ...

from sympy import symbols, init_printing, roots, solve, eye
from sympy.matrices import Matrix

init_printing()

x = symbols('x')

M = Matrix([
    [1,  2, 1],
    [6, -1, 0],
    [-1,-2,-1]
])

M

\[\begin{split}\left[\begin{matrix}1 & 2 & 1\\6 & -1 & 0\\-1 & -2 & -1\end{matrix}\right]\end{split}\]

# Find characteristic polynomial
poly = M.charpoly(x)
poly.as_poly()

\[\operatorname{Poly}{\left( x^{3} + x^{2} - 12 x, x, domain=\mathbb{Z} \right)}\]

# eigenvalues are the roots
roots(poly)

\[\begin{split}\begin{Bmatrix}-4 : 1, & 0 : 1, & 3 : 1\end{Bmatrix}\end{split}\]

7 . Use the function provided below to visualize matrix multiplication as a geometric transformation by experiment with differnt values of the matrix \(m\).

• What does a diagonal matrix do to the origianl vectors?
• What does a non-invertible matrix do to the original vectors?
• What property results in matrices that preserves the area of the parallelogram spanned by the two vectors?
• What property results in matrices that also preserve the length and angle of the original vectors?
• What additional property is necessary to preserve the orientation of the original vecotrs?
• What does the transpose of the matrix that preserves the length and angle of the original vectors do?
• Write a function that when given any two non-colinear 2D vectors u, v, finds a transformation that converts u into e1 (1,0) and v into e2 (0,1).

#  Provided function
def plot_matrix_transform(m):
    """Show the geometric effect of m on the standard unit vectors e1 and e2."""

    e1 = np.array([1,0])
    e2 = np.array([0,1])
    v1 = np.dot(m, e1)
    v2 = np.dot(m, e2)

    X = np.zeros((2,2))
    Y = np.zeros((2,2))
    pts = np.array([e1,e2,v1,v2])
    U = pts[:, 0]
    V = pts[:, 1]
    C = [0,1,0,1]

    xmin = min(-1, U.min())
    xmax = max(1, U.max())
    ymin = min(-1, V.min())
    ymax = max(-1, V.max())

    plt.figure(figsize=(6,6))
    plt.quiver(X, Y, U, V, C, angles='xy', scale_units='xy', scale=1)
    plt.axis([xmin, xmax, ymin, ymax]);

### Example usage
m = np.array([[1,2],[3,4]])
plot_matrix_transform(m)

# YOUR CODE HERE

A1 = np.diag([2,3])
A2 = np.array([[2,3],[1,1.5]])
A3 = np.array([[np.cos(1), -np.sin(1)], [np.sin(1), np.cos(1)]])
A4 = np.array([[np.cos(1), np.sin(1)], [np.sin(1), -np.cos(1)]])

print A1, la.det(A1)
print
print A2, la.det(A2)
print
print A3, la.det(A3)
print
print A4, la.det(A4)

[[2 0]
 [0 3]] 6.0

[[ 2.   3. ]
 [ 1.   1.5]] 0.0

[[ 0.5403 -0.8415]
 [ 0.8415  0.5403]] 1.0

[[ 0.5403  0.8415]
 [ 0.8415 -0.5403]] -1.0

print A3.dot(A3.T)

[[ 1.  0.]
 [ 0.  1.]]

print A4.dot(A4.T)

[[ 1.  0.]
 [ 0.  1.]]

# A diagnoal matrix simply scales the vectors
# This gives insight into what the eigendecomposition tells us
plot_matrix_transform(A1)

# A singluar matrix collapses one vector onto another
# The determinant is zero becasue the parallelogram area is zero
plot_matrix_transform(A2)

# An orthogoanl matrix preservees length and angle
# Hence the area is also preserved and the determinant is 1
# In 2D it is etiher a rotation (shown here)
plot_matrix_transform(A3)

# or a refelction
# The reflection does not preserve orietnation
# This is indicated by the determinatn being -1
plot_matrix_transform(A4)

# The tranpose of an orthogonal matrix is its inverse
plot_matrix_transform(A3.T)

def transform(u, v):
    """Retruns a matrix that converts u into e1 (1,0) and v into e2 (0,1)."""
    return la.inv(np.vstack([u, v]).T)

u = np.random.random(2)
v = np.random.random(2)

M = transform(u, v)
print u, M.dot(u)
print v, M.dot(v)

[ 0.0276  0.8173] [ 1.  0.]
[ 0.2418  0.0561] [ 0.  1.]

8 . Find and plot the least squares fit to the given values of \(x\) and \(y\) for the following:

• a constant
• a quadratic equation
• a 5th order polynomial
• a polynomial of order 50

plt.figure(figsize=(12,8))
x = np.load('x.npy')
y = np.load('y.npy')
plt.plot(x, y, 'o')

### YOUR CODE HERE
p0 = np.poly1d(np.polyfit(x, y, 0))
p2 = np.poly1d(np.polyfit(x, y, 2))
p5 = np.poly1d(np.polyfit(x, y, 5))
p50 = np.poly1d(np.polyfit(x, y, 50))

plt.plot(x, p0(x), 'b:', linewidth=2)
plt.plot(x, p2(x), 'b-', linewidth=2)
plt.plot(x, p5(x), 'g--', linewidth=2)
plt.plot(x, p50(x), 'r-.', linewidth=2)
plt.legend(['Data', 'Constant', 'Qaudratic', 'Order 5', 'Order 50'], loc='best', fontsize=20);

/Users/cliburn/anaconda/lib/python2.7/site-packages/numpy/lib/polynomial.py:588: RankWarning: Polyfit may be poorly conditioned
  warnings.warn(msg, RankWarning)

%load_ext version_information

%version_information numpy, scipy

Software	Version
Python	2.7.9 64bit [GCC 4.2.1 (Apple Inc. build 5577)]
IPython	2.3.1
OS	Darwin 13.4.0 x86_64 i386 64bit
numpy	1.9.1
scipy	0.14.0
Thu Jan 22 15:43:33 2015 EST