The Kitchin Research Group

Numpy has some gotcha features for linear algebra purists. The first is that a 1d array is neither a row, nor a column vector. That is, \(a\) = \(a^T\) if \(a\) is a 1d array. That means you can take the dot product of \(a\) with itself, without transposing the second argument. This would not be allowed in Matlab.

import numpy as np

a = np.array([0, 1, 2])
print a.shape
print a
print a.T

print
print np.dot(a, a)
print np.dot(a, a.T)

>>> >>> (3L,)
[0 1 2]
[0 1 2]
>>>
5
5

Compare the previous behavior with this 2d array. In this case, you cannot take the dot product of \(b\) with itself, because the dimensions are incompatible. You must transpose the second argument to make it dimensionally consistent. Also, the result of the dot product is not a simple scalar, but a 1 × 1 array.

b = np.array([[0, 1, 2]])
print b.shape
print b
print b.T

print np.dot(b, b)    # this is not ok, the dimensions are wrong.
print np.dot(b, b.T)
print np.dot(b, b.T).shape

(1L, 3L)
[[0 1 2]]
[[0]
 [1]
 [2]]
>>> Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
ValueError: objects are not aligned
[[5]]
(1L, 1L)

Try to figure this one out! x is a column vector, and y is a 1d vector. Just by adding them you get a 2d array.

x = np.array([[2], [4], [6], [8]])
y = np.array([1, 1, 1, 1, 1, 2])
print x + y

>>> [[ 3  3  3  3  3  4]
 [ 5  5  5  5  5  6]
 [ 7  7  7  7  7  8]
 [ 9  9  9  9  9 10]]

Or this crazy alternative way to do the same thing.

x = np.array([2, 4, 6, 8])
y = np.array([1, 1, 1, 1, 1, 1, 2])

print x[:, np.newaxis] + y

>>> >>> [[ 3  3  3  3  3  3  4]
 [ 5  5  5  5  5  5  6]
 [ 7  7  7  7  7  7  8]
 [ 9  9  9  9  9  9 10]]

In the next example, we have a 3 element vector and a 4 element vector. We convert \(b\) to a 2D array with np.newaxis, and compute the outer product of the two arrays. The result is a 4 × 3 array.

a = np.array([1, 2, 3])
b = np.array([10, 20, 30, 40])

print a * b[:, np.newaxis]

>>> >>> [[ 10  40  90]
 [ 20  80 180]
 [ 30 120 270]
 [ 40 160 360]]

These concepts are known in numpy as array broadcasting. See http://www.scipy.org/EricsBroadcastingDoc and http://docs.scipy.org/doc/numpy/user/basics.broadcasting.html for more details.

These are points to keep in mind, as the operations do not strictly follow the conventions of linear algebra, and may be confusing at times.

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source