1 import math
2 import operator
3 import sets
4 from types import *
5
6 '''A module providing some basic linear algebra operations.
7
8 The functions cover the standard basic operations on vectors and matrices.
9 If you are looking for an efficient and complete linear algebra module,
10 this isn't it. This module is intended just as a simple (no install,
11 no dependencies) way of double-checking some of the more tedious
12 linear algebra stuff as I worked my way through a linear algebra text.
13
14 The tests embedded in the docstring of each function shows basic usage of the
15 function, and all tests can be run by executing:
16
17 python linearalgebra.py
18
19 No output means success. If you want the output to display what is being
20 tested, run:
21
22 python linearalgebra.py -v
23
24 Some useful math functions and constants are re-exported from the standard
25 Python math module, as a convenience for interactive scripting. These include
26 e, pi, sqrt, cos, sin, tan, acos, asin, and atan.
27
28 Conventions for function parameters:
29 * u and v always refer to a vector;
30 * s and c always refers to a scalar;
31 * A, B, and C always refer to a matrix;
32 * m, n, r always refer to a scalar that is the dimension of a matrix
33 (e.g., m = 3, n = 2 for a 3 x 2 or m x n matrix);
34 * x may be either a vector or a scalar;
35 * any of these may be postfixed by an integer, and is of the same type
36 as the unpostfixed variable name.'''
37
38
39
40
41 degrees = math.degrees
42 radians = math.radians
43 e = math.e
44 pi = math.pi
45 sqrt = math.sqrt
46 cos = math.cos
47 sin = math.sin
48 tan = math.tan
49 acos = math.acos
50 asin = math.asin
51 atan = math.atan
67
76 ''' Verifies that the indices that were passed to the given decorator for
77 the given function are all ints. '''
78 def wrapper(*args):
79 for arg in args:
80 if not isinstance(arg, int):
81 msg = "Arg for '%s' decorator should be an int: %s"
82 raise LAValueError(msg % (func.__name__, repr(arg)))
83 return func(*args)
84 wrapper.__doc__ = func.__doc__
85 wrapper.__name__ = func.__name__
86 return wrapper
87
91 ''' Decorates a function such that an exception is raised if
92 the vector arg at position index1 does not have the same size as
93 the vector arg at position index2, starting at 0 for the first arg.
94 If the vectors have the same size, then the function is executed
95 as it would have been.'''
96 def decorator(func):
97 def _decorator(*args):
98 vec_u = args[index1]
99 vec_v = args[index2]
100 if not isV(vec_u):
101 raise InvalidVectorException("First vector arg is not a vector.")
102 elif not isV(vec_v):
103 raise InvalidVectorException("Second vector arg is not a vector.")
104 len_u = len(vec_u)
105 len_v = len(vec_v)
106 if len_u != len_v:
107 raise IncompatibleVectorException(
108 "Vectors must be of equal size, not %d and %d."
109 % (len_u, len_v))
110 return func(*args)
111 _decorator.__doc__ = func.__doc__
112 _decorator.__name__ = func.__name__
113 return _decorator
114 return decorator
115
119 ''' Decorates a function such that an exception is raised if
120 the matrix arg at position index1 does not have the same dimensions
121 as the matrix arg at position index2, starting at 0 for the first arg.
122 If the matrices have the same number of rows and columns, then the
123 function is executed as it would have been.'''
124 def decorator(func):
125 def _decorator(*args):
126 (m1, n1) = dimensionsM(args[index1])
127 (m2, n2) = dimensionsM(args[index2])
128 if m1 != m2 or n1 != n2:
129 msg = "Matrices must have same dimensions. Matrix 1 is %d x %d, but matrix 2 is %d x %d."
130 raise IncompatibleMatrixException(msg % (m1, n1, m2, n2))
131 return func(*args)
132 _decorator.__doc__ = func.__doc__
133 _decorator.__name__ = func.__name__
134 return _decorator
135 return decorator
136
137
138 @__intArgs
139 -def __verifyS(*argIndices):
140 ''' Decorates a function such that an exception is raised if the scalar arg at
141 any of the supplied indices is not actually a vector, and otherwise executes
142 the function as usual.'''
143 def decorator(func):
144 def _decorator(*args):
145 for argIndex in argIndices:
146 arg = args[argIndex]
147 if not isS(arg):
148 raise LAValueError("Arg at index %s is not a scalar: %s" % (argIndex, arg))
149 return func(*args)
150 _decorator.__doc__ = func.__doc__
151 _decorator.__name__ = func.__name__
152 return _decorator
153 return decorator
154
155
156 @__intArgs
157 -def __verifyV(*argIndices):
158 ''' Decorates a function such that an exception is raised if the vector
159 arg at any of the supplied indices is not actually a vector, and otherwise
160 executes the function as usual. '''
161 def decorator(func):
162 def _decorator(*args):
163 for argIndex in argIndices:
164 u = args[argIndex]
165 if not isV(u):
166 raise InvalidVectorException("Arg should be a vector: %s" % repr(u))
167 return func(*args)
168 _decorator.__doc__ = func.__doc__
169 _decorator.__name__ = func.__name__
170 return _decorator
171 return decorator
172
173
174 @__intArgs
175 -def __verifyM(*argIndices):
176 ''' Decorates a function such that an exception is raised if the matrix
177 arg at any of the supplied indices is not a matrix, and otherwise executes
178 the function as usual. '''
179 def decorator(func):
180 def _decorator(*args):
181 for argIndex in argIndices:
182 A = args[argIndex]
183 if not isM(A):
184 raise InvalidMatrixException("Arg is not a matrix: %s" % repr(A))
185 return func(*args)
186 _decorator.__doc__ = func.__doc__
187 _decorator.__name__ = func.__name__
188 return _decorator
189 return decorator
190
195 ''' Decorates a function such that an exception is raised if
196 the matrix arg at position index1 is not column-row-compatible
197 with the matrix arg at position index2, starting at 0 for the first arg.
198
199 Matrices A (of size a x b) and B (of size c x d) are column-row-compatible
200 if and only if b == c; i.e., if under variable name reassignment,
201 A has size m x n and B has size n x r.'''
202 def decorator(func):
203 def _decorator(*args):
204 (m1, n1) = dimensionsM(args[index1])
205 (m2, n2) = dimensionsM(args[index2])
206 if n1 != m2:
207 msg = "Matrices are not size-compatible. Matrix 1 has %s column(s), but matrix 2 has %s row(s) instead of %s"
208 raise IncompatibleMatrixException(msg % (n1, m2, n1))
209 return func(*args)
210 _decorator.__doc__ = func.__doc__
211 _decorator.__name__ = func.__name__
212 return _decorator
213 return decorator
214
218 ''' Determine whether the given arg is a matrix, which must
219 be a rectangular array of numbers, represented as a list/tuple of
220 list/tuple elements, each of which is a vector of equal size.
221
222 >>> isM([])
223 True
224 >>> isM([[]])
225 True
226 >>> isM([1,2])
227 False
228 >>> isM([[1,2], (3,4)])
229 True
230 >>> isM([[1,2], [3,4,5]])
231 False
232 >>> isM([[1,2], ['0', 1]])
233 False
234 >>> isM([[.5, -1, 9], [4, 3, 19], [3, 6, 8], [0, 0, 1]])
235 True
236 '''
237 if not isinstance(A, list) and not isinstance(A, tuple):
238 return False
239 m = len(A)
240 if m == 0:
241 return True
242 n = None
243 for row in A:
244 if not isinstance(row, list) and not isinstance(row, tuple):
245 return False
246 if not isV(row):
247 return False
248 if n is None:
249 n = len(row)
250 elif n != len(row):
251 return False
252 return True
253
256 ''' Calculates the dimensions of the given matrix, returning a pair consisting
257 of the number of rows and the number of columns, raising an exception if any of
258 the rows of the matrix have a different number of elements.
259
260 >>> dimensionsM([])
261 (0, 0)
262 >>> dimensionsM([[1,2,3]])
263 (1, 3)
264 >>> dimensionsM([[1], [2], [3]])
265 (3, 1)
266 >>> dimensionsM([[1, 2], [1, 2, 3]])
267 Traceback (most recent call last):
268 ...
269 InvalidMatrixException: Arg is not a matrix: [[1, 2], [1, 2, 3]]
270 '''
271 m = len(A)
272 if m == 0:
273 return (0,0)
274 n = len(A[0])
275 for (i, row) in enumerate(A):
276 l = len(row)
277 if n != l:
278 msg = "Matrix has rows of different length: row 0 has %d element(s) and row %d has %d element(s)." % (n, i, l)
279 raise InvalidMatrixException(msg)
280 return (m, n)
281
284 ''' Determine whether the given matrix is a square matrix, with the same number of
285 rows and columns.
286
287 >>> isSquareM([])
288 True
289 >>> isSquareM([[]])
290 False
291 >>> isSquareM([[1,1], [2,2], [3,3]])
292 False
293 >>> isSquareM([[1,2,3], [4,5,6], [7,8,9]])
294 True
295 >>> isSquareM([[1], [1,2]])
296 Traceback (most recent call last):
297 ...
298 InvalidMatrixException: Arg is not a matrix: [[1], [1, 2]]
299 '''
300 (m, n) = dimensionsM(A)
301 return m == n
302
305 ''' Determine whether the given matrix is a diagonal matrix -- i.e., a square matrix with
306 all nondiagonal entries equal to 0.
307
308 >>> isDiagonalM([])
309 True
310 >>> isDiagonalM([[0]])
311 True
312 >>> isDiagonalM([[0, 0], [0, 0], [0, 0]])
313 False
314 >>> isDiagonalM([[1, 0], [0, 2]])
315 True
316 >>> isDiagonalM([[1, 0], [1, 0]])
317 False
318 '''
319 if not isSquareM(A):
320 return False
321 for (i, row) in enumerate(A):
322 for (j, elem) in enumerate(row):
323 if i != j and elem != 0:
324 return False
325 return True
326
329 ''' Determine whether the given matrix is a scalar matrix -- i.e., a diagonal matrix
330 with all diagonal entries equal to each other.
331
332 >>> isScalarM([])
333 True
334 >>> isScalarM([[2]])
335 True
336 >>> isScalarM([[3, 0], [0, 3]])
337 True
338 >>> isScalarM([[3, 0], [1, 3]])
339 False
340 '''
341 if not isSquareM(A):
342 return False
343 elif len(A) == 0:
344 return True
345 diagVal = A[0][0]
346 for (i, row) in enumerate(A):
347 for (j, elem) in enumerate(row):
348 if (i != j and elem != 0) or (i == j and elem != diagVal):
349 return False
350 return True
351
354 ''' Determine whether the given matrix is an identity matrix -- i.e., a scalar matrix
355 with all diagonal entries equal to 1.
356
357 >>> isIdentityM([])
358 False
359 >>> isIdentityM([[0]])
360 False
361 >>> isIdentityM([[1]])
362 True
363 >>> isIdentityM([[1, 1], [1, 1]])
364 False
365 >>> isIdentityM([[1, 0], [0, 1]])
366 True
367 >>> isIdentityM([[1, 0], [0, 1], [0, 0]])
368 False
369 >>> isIdentityM([[2, 0], [0, 2]])
370 False
371 '''
372 if not isSquareM(A):
373 return False
374 elif len(A) == 0:
375 return False
376 for (i, row) in enumerate(A):
377 for (j, elem) in enumerate(row):
378 if (i != j and elem != 0) or (i == j and elem != 1):
379 return False
380 return True
381
384 ''' Determine whether the given matrix is symmetric -- i.e., whether it is equal
385 to its own transpose.
386
387 >>> isSymmetricM([[1, 0], [0, 1]])
388 True
389 >>> isSymmetricM([[0, 1], [1, 0]])
390 True
391 >>> isSymmetricM([[1, 2], [2, 0]])
392 True
393 >>> isSymmetricM([[1, 2], [1, 2]])
394 False
395 >>> isSymmetricM([[1, 3, 2], [3, 5, 0], [2, 0, 4]])
396 True
397 '''
398 return areEqualM(A, transposeM(A))
399
401 ''' Helper function for plusM and minusM. '''
402 (m, n) = dimensionsM(A)
403 C = [None] * m
404 for i in range(m):
405 row = [None] * n
406 for j in range(n):
407 row[j] = func(A[i][j], B[i][j])
408 C[i] = row
409 return C
410
411 @__verifyM(0, 1)
412 @__samesizeM(0, 1)
413 -def plusM(A, B):
414 ''' Add the given matrices, which must have the same number of rows and columns.
415
416 >>> plusM([[1, 2], [3, 4]], [[4, 3], [2, 1]])
417 [[5, 5], [5, 5]]
418 >>> plusM([[1,2]], [[3, 4]]), plusM([[1,2]], [[3, 4]])[0] == plusV([1,2], [3,4])
419 ([[4, 6]], True)
420 >>> plusM([[1]], [[1, 2]])
421 Traceback (most recent call last):
422 ...
423 IncompatibleMatrixException: Matrices must have same dimensions. Matrix 1 is 1 x 1, but matrix 2 is 1 x 2.
424 '''
425 return __pairwiseM(A, B, operator.add)
426
427 @__verifyM(0, 1)
428 @__samesizeM(0, 1)
429 -def minusM(A, B):
430 ''' Subtract the second matrix from the first, both of which must have the
431 same number of rows and columns.
432
433 >>> minusM([[3, 4], [1, 2]], [[1, 1], [1, 1]])
434 [[2, 3], [0, 1]]
435 >>> minusM([[1,2]], [[3, 4]]), minusM([[1,2]], [[3, 4]])[0] == minusV([1,2], [3,4])
436 ([[-2, -2]], True)
437 >>> minusM([[1]], [[1, 2]])
438 Traceback (most recent call last):
439 ...
440 IncompatibleMatrixException: Matrices must have same dimensions. Matrix 1 is 1 x 1, but matrix 2 is 1 x 2.
441 '''
442 return __pairwiseM(A, B, operator.sub)
443
444 @__verifyS(0)
445 @__verifyM(1)
446 -def timesM(s, A):
447 ''' Multiply the matrix A by the scalar s.
448
449 >>> timesM(2, [[0,1], [1, 0]])
450 [[0, 2], [2, 0]]
451 >>> timesM(-1, [[1,2,3,4]])
452 [[-1, -2, -3, -4]]
453 >>> timesM(3, [[]])
454 [[]]
455 >>> timesM(1, [])
456 []
457 >>> timesM(3, [[1,2,3], [4,5,6], [7,8,9], [10, 11, 12]])
458 [[3, 6, 9], [12, 15, 18], [21, 24, 27], [30, 33, 36]]
459 '''
460 return [[s*elem for elem in row] for row in A]
461
464 ''' Return the negative of the given matrix, or the matrix multiplied
465 by scalar -1.
466
467 >>> negativeM([])
468 []
469 >>> negativeM([[]])
470 [[]]
471 >>> negativeM(((-1, 2), (2, -3)))
472 [[1, -2], [-2, 3]]
473 '''
474 return timesM(-1, A)
475
478 ''' Generate the transpose of the given m x n matrix, which is the n x m
479 matrix with where each (i,j) element of the original matrix is the (j,i)
480 element of the new matrix.
481
482 >>> transposeM([[1, 2], [3, 4], [5, 6]])
483 [[1, 3, 5], [2, 4, 6]]
484 >>> transposeM([])
485 []
486 >>> transposeM([[]])
487 []
488 >>> transposeM([[1,2,3], [4,5,6], [7,8,9]])
489 [[1, 4, 7], [2, 5, 8], [3, 6, 9]]
490 >>> transposeM(transposeM([[1,2,3], [4,5,6], [7,8,9]]))
491 [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
492 >>> transposeM([[1,2,3]])
493 [[1], [2], [3]]
494 '''
495 (m, n) = dimensionsM(A)
496 C = [[None] * m for i in range(n)]
497 for i in range(m):
498 for j in range(n):
499 C[j][i] = A[i][j]
500 return C
501
502 @__columnRowCompatibleM(0, 1)
503 @__verifyM(0, 1)
504 -def productM(A, B):
505 ''' Calculate the product of matrices A and B, where A must have the
506 same number of columns and B has rows. The result for A of size m x n
507 and B of size n x r is an m x r matrix.
508
509 >>> productM([[1,1], [1,1]], [[1,1], [1,1]])
510 [[2, 2], [2, 2]]
511 >>> productM([], [])
512 []
513 >>> productM([[1, 3, -1], [-2, -1, 1]], [[-4, 0, 3, -1], [5, -2, -1, 1], [-1, 2, 0, 6]])
514 [[12, -8, 0, -4], [2, 4, -5, 7]]
515 '''
516 (m, n) = dimensionsM(A)
517 (n, r) = dimensionsM(B)
518 C = [[None] * r for i in range(m)]
519 for i in range(m):
520 for j in range(r):
521 C[i][j] = productV(A[i], [row[j] for row in B])
522 return C
523
525 ''' Determine whether the two scalars are equal to accuracy of epsilon.
526
527 >>> areEqualS(0.9000000000000002, 0.9000000000000001)
528 True
529 >>> areEqualS(1, 1)
530 True
531 '''
532 return abs(s1-s2) < epsilon
533
536 ''' Determine whether the two vectors are equal by determining if the
537 pairwise components are equal to accuracy of epsilon.
538
539 >>> areEqualV([2], [2])
540 True
541 >>> areEqualV([1.0, 0.9000000000000002], [1, 0.9000000000000001])
542 True
543 '''
544 return __allTrue(map(lambda (x,y): areEqualS(x,y), zip(u, v)))
545
548 ''' Determine whether the two matrices are equal by determining if
549 the corresponding components are equal to accuracy of epsilon.
550
551 >>> areEqualM([], [])
552 True
553 >>> areEqualM([[1, 2, 3], [3.9000000000000001, 2, 1]], [[1, 2, 3], [3.9, 2, 1]])
554 True
555 '''
556 (m1, n1) = dimensionsM(A)