** Potential gotchas in linear algebra in numpy
:PROPERTIES:
:categories: linear algebra, gotcha
:date: 2013/03/12 22:19:53
:updated: 2013/03/12 22:19:53
:END:
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.
#+BEGIN_SRC python :session
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)
#+END_SRC
#+RESULTS:
:
: >>> >>> (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 \times 1 array.
#+BEGIN_SRC python :session
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
#+END_SRC
#+RESULTS:
#+begin_example
(1L, 3L)
[[0 1 2]]
[[0]
[1]
[2]]
>>> Traceback (most recent call last):
File "", line 1, in
ValueError: objects are not aligned
[[5]]
(1L, 1L)
#+end_example
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.
#+BEGIN_SRC python :session
x = np.array([[2], [4], [6], [8]])
y = np.array([1, 1, 1, 1, 1, 2])
print x + y
#+END_SRC
#+RESULTS:
:
: >>> [[ 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.
#+BEGIN_SRC python :session
x = np.array([2, 4, 6, 8])
y = np.array([1, 1, 1, 1, 1, 1, 2])
print x[:, np.newaxis] + y
#+END_SRC
#+RESULTS:
:
: >>> >>> [[ 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 \times 3 array.
#+BEGIN_SRC python :session
a = np.array([1, 2, 3])
b = np.array([10, 20, 30, 40])
print a * b[:, np.newaxis]
#+END_SRC
#+RESULTS:
:
: >>> >>> [[ 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.