** Sums products and linear algebra notation - avoiding loops where possible
:PROPERTIES:
:categories: Linear algebra
:date: 2013/02/26 09:00:00
:updated: 2013/02/27 13:12:15
:END:
[[http://matlab.cheme.cmu.edu/2012/01/03/sums-products-and-linear-algebra-notation-avoiding-loops-where-possible/][Matlab comparison]]
Today we examine some methods of linear algebra that allow us to
avoid writing explicit loops in Matlab for some kinds of
mathematical operations.
Consider the operation on two vectors $\bf{a}$
and $\bf{b}$.
$$y=\sum\limits_{i=1}^n a_ib_i$$
a = [1 2 3 4 5]
b = [3 6 8 9 10]
*** Old-fashioned way with a loop
We can compute this with a loop, where you initialize y, and then
add the product of the ith elements of a and b to y in each
iteration of the loop. This is known to be slow for large vectors
#+BEGIN_SRC python
a = [1, 2, 3, 4, 5]
b = [3, 6, 8, 9, 10]
sum = 0
for i in range(len(a)):
sum = sum + a[i] * b[i]
print sum
#+END_SRC
#+RESULTS:
: 125
This is an old fashioned style of coding. A more modern, pythonic approach is:
#+BEGIN_SRC python
a = [1, 2, 3, 4, 5]
b = [3, 6, 8, 9, 10]
sum = 0
for x,y in zip(a,b):
sum += x * y
print sum
#+END_SRC
#+RESULTS:
: 125
*** The numpy approach
The most compact method is to use the methods in numpy.
#+BEGIN_SRC python
import numpy as np
a = np.array([1, 2, 3, 4, 5])
b = np.array([3, 6, 8, 9, 10])
print np.sum(a * b)
#+END_SRC
#+RESULTS:
: 125
*** Matrix algebra approach.
The operation defined above is actually a dot product. We an directly compute the dot product in numpy. Note that with 1d arrays, python knows what to do and does not require any transpose operations.
#+BEGIN_SRC python
import numpy as np
a = np.array([1, 2, 3, 4, 5])
b = np.array([3, 6, 8, 9, 10])
print np.dot(a, b)
#+END_SRC
#+RESULTS:
: 125
*** Another example
Consider $y = \sum\limits_{i=1}^n w_i x_i^2$. This operation is like a weighted sum of squares.
The old-fashioned way to do this is with a loop.
#+BEGIN_SRC python
w = [0.1, 0.25, 0.12, 0.45, 0.98];
x = [9, 7, 11, 12, 8];
y = 0
for wi, xi in zip(w,x):
y += wi * xi**2
print y
#+END_SRC
#+RESULTS:
: 162.39
Compare this to the more modern numpy approach.
#+BEGIN_SRC python
import numpy as np
w = np.array([0.1, 0.25, 0.12, 0.45, 0.98])
x = np.array([9, 7, 11, 12, 8])
y = np.sum(w * x**2)
print y
#+END_SRC
#+RESULTS:
: 162.39
We can also express this in matrix algebra form. The operation is equivalent to $y = \vec{x} \cdot D_w \cdot \vec{x}^T$ where $D_w$ is a diagonal matrix with the weights on the diagonal.
#+BEGIN_SRC python
import numpy as np
w = np.array([0.1, 0.25, 0.12, 0.45, 0.98])
x = np.array([9, 7, 11, 12, 8])
y = np.dot(x, np.dot(np.diag(w), x))
print y
#+END_SRC
#+RESULTS:
: 162.39
This last form avoids explicit loops and sums, and relies on fast linear algebra routines.
*** Last example
Consider the sum of the product of three vectors. Let $y = \sum\limits_{i=1}^n w_i x_i y_i$. This is like a weighted sum of products.
#+BEGIN_SRC python
import numpy as np
w = np.array([0.1, 0.25, 0.12, 0.45, 0.98])
x = np.array([9, 7, 11, 12, 8])
y = np.array([2, 5, 3, 8, 0])
print np.sum(w * x * y)
print np.dot(w, np.dot(np.diag(x), y))
#+END_SRC
#+RESULTS:
: 57.71
: 57.71
*** Summary
We showed examples of the following equalities between traditional
sum notations and linear algebra
$$\bf{a}\bf{b}=\sum\limits_{i=1}^n a_ib_i$$
$$\bf{x}\bf{D_w}\bf{x^T}=\sum\limits_{i=1}^n w_ix_i^2$$
$$\bf{x}\bf{D_w}\bf{y^T}=\sum\limits_{i=1}^n w_i x_i y_i$$
These relationships enable one to write the sums as a single line of
python code, which utilizes fast linear algebra subroutines, avoids
the construction of slow loops, and reduces the opportunity for
errors in the code. Admittedly, it introduces the opportunity for
new types of errors, like using the wrong relationship, or linear
algebra errors due to matrix size mismatches.