Matrices in Python

BMI 214
Something weird to start us off, which explains why youre reading this additional document:
Lets define a 2 dimensional matrix M M will contain 2 lists, with 3 zeros in each. >>> M=[[0]*3]*2 >>> print M [[0, 0, 0], [0, 0, 0]] >>> M[0][0]=5 >>> M [[5, 0, 0], [5, 0, 0]]
Huh? When I update the value of position 0 in the first list in the matrix, it also updates the value of position 0 in the second list. WHY?? An important point you must be aware of: Python variables represent pointers to objects in memory. In fact, multiple variables can point to the same object in memory, and so if you manipulate that object via one of the variables, those changes will be reflected even when you look at the object (which is the same object in memory) through a different variable. This is what is happening here. The list defined by [0]*3 is actually multiplied when creating M, but both lists in M actually point to the same object in memory. That is, M[0] and M[1] are actually the same object. How do we get around this, in order to create a matrix? You have a few options: 1. You can build up the matrices you need in project1 as you go along. If you build them by appending the values to the list once youve calculated them, youre actually adding doubles to the list, and so wont run into this problem. So, if were building the second row of the score matrix, we can do it by calculating each columns value, and using append or extend to add that value to M[1]. For example: M = [[0]*6] # initialize M to have one row with 6 zeros, as first row in matirx M.append([]) # initialize second row of M to be empty M[1].append(3.0) # append the value of first element in second row M[1].append(5) # append the value of second element in second row print M prints out: [[0, 0, 0, 0, 0, 0], [3.0, 5]] Equivalent to matrix:
2. You can choose to learn and use NumPy (http://www.scipy.org/). It's much faster than using lists-of-lists-of-lists for multi-dimensional arrays, and you can do things like adding them or multiplication directly. It does have a learning curve, but it is considered superior to anything else out there. 3. If you want to do it in pure Python, with initialization of a matrix of defined size, to zero values, you can do the following:
def zeros(*shape): if len(shape) == 0: return 0 car = shape[0] cdr = shape[1:] return [zeros(*cdr) for i in range(car)] >>> zeros(3) [0, 0, 0] >>> zeros(2,3) [[0, 0, 0], [0, 0, 0]] >>> zeros(1,2,3) [[[0, 0, 0], [0, 0, 0]]] >>> zeros(1,3,4) [[[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]]
Something we didnt go over in Section: You can create a list like this, using the + (concatenation): a = [0] + [0] + [0] This result in a is equal to [0, 0, 0]
With thanks to Zach Pincus

COMPARISON OF MATLAB, OCTAVE, AND PYTHON/SCIPY/PYLAB EXAMPLES, FOR MATH 694

ED BUELER

Trefethen & Baus Numerical Linear Algebra uses Matlab (http://www.mathworks.com) for very good reason. Matlab was designed by Cleve Moler for teaching numerical linear algebra. Matlab has, since its invention, become a powerful programming language and engineering tool. The best single tool, in my opinion. There are open source alternatives to Matlab, and theyll work for this course. I like free, open source software. Among other open source alternatives, GNU Octave is intended to be a Matlab clone, and the examples below work in an identical way in Matlab and in Octave. Incompatibilities between Octave and Matlab are now treated by the Octave developers as a reportable bug. I will mostly use Octave myself for teaching, but Ill make an eort to test most examples in both Octave and Matlab. To download Octave, go to http://www.gnu.org/software/octave/. The general purpose language Python has developed in the direction of Matlab functionality with the scipy (http://www.scipy.org/) and pylab (http://matplotlib.sourceforge.net/) projects. The ipython interactive shell gives, with pylab, the most Matlab-like experience. I like Python. The examples below hint at the computer language dierences, and the dierent modes of thought, between Matlab/Octave and Python. I suspect that only students who already use Python will nd Python/scipy/pylab eective for this course. Some brief how-to comments might help compare examples below. The Matlab/Octave examples legendre.m and hello.m are scripts. These are run by starting Matlab/Octave, making sure that the path includes the directory containing the examples. Then type the name of the script at the prompt, without the.m: >> legendre or >> hello. For the rst two Python examples type run legendre.py or run hello.py at the ipython prompt or python legendre.py, python hello.py at a ordinary shell prompt. The last example mgs.m, mgs.py is a function which needs an input, and Ill address that later. legendre.m
(Matlab & Octave version)

legendre.py

(Python version)
% Trefethen & Bau, page 64 x = (-128:128)/128; A = [x.^0 x.^1 x.^2 x.^3]; [Q,R] = qr(A,0); scale = Q(257,:); Q = Q * diag(1./ scale); plot(Q)
# Trefethen & Bau, page 64 from pylab import * x = linspace(-1.0,1.0,257).reshape((257,1)) A = concatenate((x**0, x**1, x**2, x**3),axis=1) [Q,R] = qr(A) scale = Q[256,:] Q = dot(Q,diag(1.0 / scale)) plot(Q), show()
Date: January 22, 2009. For these examples, follow links at http://www.dms.uaf.edu/~bueler/Math694S09.htm.

hello.m

% assembles HELLO matrix % see Trefethen&Bau lecture 9, exer bl = ones(8,6); H = bl; H(1:3,3:4) = zeros(3,2); H(6:8,3:4) = zeros(3,2); E = bl; E(3,3:6) = zeros(1,4); E(6,3:6) = zeros(1,4); L = bl; L(1:6,3:6) = zeros(6,4); O = bl; O(3:6,3:4) = zeros(4,2); HELLO = zeros(15,40); HELLO(2:9,2:7) = H; HELLO(3:10,10:15) = E; HELLO(4:11,18:23) = L; HELLO(5:12,26:31) = L; HELLO(6:13,34:39) = O; spy(HELLO)

hello.py

# assembles HELLO matrix # see Trefethen&Bau lecture 9, exer from pylab import ones, zeros, spy, show bl = ones((8,6)) H = bl.copy() H[0:3,2:4] = zeros((3,2)) H[5:8,2:4] = zeros((3,2)) E = bl.copy() E[2,2:6] = zeros((1,4)) E[5,2:6] = zeros((1,4)) L = bl.copy() L[0:6,2:6] = zeros((6,4)) O = bl.copy() O[2:6,2:4] = zeros((4,2)) HELLO = zeros((15,40)) HELLO[1:9,1:7] = H HELLO[2:10,9:15] = E HELLO[3:11,17:23] = L HELLO[4:12,25:31] = L HELLO[5:13,33:39] = O spy(HELLO,marker=.); show()
function [Q,R] = mgs(A); % MGS computes QR factors twoeps = 2*eps; [m n] = size(A); R = zeros(n,n); scal = max(max(abs(A))); if scal == 0 Q = eye(m,n); return, end Q = A; for i = 1:n r = norm(Q(:,i),2); R(i,i) = r; if abs(r)/scal < twoeps error(nearly zero col), end w = Q(:,i)/r; Q(:,i) = w; for j = i+1:n r = w*Q(:,j); R(i,j) = r; Q(:,j) = Q(:,j)-r*w; end end