1 import math
2 import operator
3 import sets
4
5 import functools
6 from types import *
7
8 """A module providing some basic linear algebra operations.
9
10 The functions cover the standard basic operations on vectors and matrices.
11 If you are looking for an efficient and complete linear algebra module,
12 this isn't it. This module is intended just as a simple (no install,
13 no dependencies) way of double-checking some of the more tedious
14 linear algebra stuff as I worked my way through a linear algebra text.
15
16 The tests embedded in the docstring of each function shows basic usage of the
17 function, and all tests can be run by executing:
18
19 python linearalgebra.py
20
21 No output means success. If you want the output to display what is being
22 tested, run:
23
24 python linearalgebra.py -v
25
26 Some useful math functions and constants are re-exported from the standard
27 Python math module, as a convenience for interactive scripting. These include
28 e, pi, sqrt, cos, sin, tan, acos, asin, and atan.
29
30 Conventions for function parameters:
31 * u and v always refer to a vector;
32 * s and c always refer to a scalar;
33 * A, B, and C always refer to a matrix;
34 * m, n, r always refer to a scalar that is the dimension of a matrix
35 (e.g., m = 3, n = 2 for a 3 x 2 or m x n matrix);
36 * x may be either a vector or a scalar;
37 * any of these may be postfixed by an integer, and is of the same type
38 as the unpostfixed variable name (e.g., m2 is a scalar)."""
39
40
41
42 degrees = math.degrees
43 radians = math.radians
44 e = math.e
45 pi = math.pi
46 sqrt = math.sqrt
47 cos = math.cos
48 sin = math.sin
49 tan = math.tan
50 acos = math.acos
51 asin = math.asin
52 atan = math.atan
68
73 """ Decorates a function such that an exception is raised if the
74 arg at index1 or index 2 is not a vector or if they do not have
75 the same size, and otherwise executes the function as usual. """
76 def decorator(func):
77 @functools.wraps(func)
78 def _decorator(*args):
79 if not is_vector(args[index1]):
80 raise InvalidVectorException("First vector arg to __check_vector_same_size is not a vector.")
81 elif not is_vector(args[index2]):
82 raise InvalidVectorException("Second vector arg to __check_vector_same_size is not a vector.")
83 if len(args[index1]) != len(args[index2]):
84 raise IncompatibleVectorException(
85 "Vectors must be of equal size, not %d and %d."
86 % (len(args[index1]), len(args[index2])))
87 return func(*args)
88 return _decorator
89 return decorator
90
93 """ Decorates a function such that an exception is raised if the
94 arg at index1 or index2 is not a matrix or if they to not have
95 the same dimensions, and otherwise executes the function as usual. """
96 def decorator(func):
97 @functools.wraps(func)
98 def _decorator(*args):
99 (m1, n1) = matrix_dimensions(args[index1])
100 (m2, n2) = matrix_dimensions(args[index2])
101 if m1 != m2 or n1 != n2:
102 msg = "Matrices must have same dimensions. Matrix 1 is %d x %d, but matrix 2 is %d x %d."
103 raise IncompatibleMatrixException(msg % (m1, n1, m2, n2))
104 return func(*args)
105 return _decorator
106 return decorator
107
110 """ Decorates a function such that an exception is raised if the arg at
111 any of the supplied indices is not a scalar, and otherwise executes
112 the function as usual."""
113 def decorator(func):
114 @functools.wraps(func)
115 def _decorator(*args):
116 for arg_index in arg_indices:
117 arg = args[arg_index]
118 if not is_scalar(arg):
119 raise LAValueError("Arg at index %s is not a scalar: %s" % (arg_index, arg))
120 return func(*args)
121 return _decorator
122 return decorator
123
126 """ Decorates a function such that an exception is raised if the vector
127 arg at any of the supplied indices is not actually a vector, and otherwise
128 executes the function as usual. """
129 def decorator(func):
130 @functools.wraps(func)
131 def _decorator(*args):
132 for arg_index in arg_indices:
133 u = args[arg_index]
134 if not is_vector(u):
135 raise InvalidVectorException("Arg should be a vector: %s" % repr(u))
136 return func(*args)
137 return _decorator
138 return decorator
139
142 """ Decorates a function such that an exception is raised if the matrix
143 arg at any of the supplied indices is not a matrix, and otherwise executes
144 the function as usual. """
145 def decorator(func):
146 @functools.wraps(func)
147 def _decorator(*args):
148 for arg_index in arg_indices:
149 A = args[arg_index]
150 if not is_matrix(A):
151 raise InvalidMatrixException("Arg is not a matrix: %s" % repr(A))
152 return func(*args)
153 return _decorator
154 return decorator
155
159 """ Decorates a function such that an exception is raised if
160 the multiplication of matrices (args[index1], args[index2]) is not
161 permitted due to incompatible size. I.e., if the number of the columns
162 in the first matrix is not equal to the number of rows in the second. """
163 def decorator(func):
164 @functools.wraps(func)
165 def _decorator(*args):
166 (m1, n1) = matrix_dimensions(args[index1])
167 (m2, n2) = matrix_dimensions(args[index2])
168 if n1 != m2:
169 msg = "Matrices are not size-compatible. Matrix 1 has %s column(s), but matrix 2 has %s row(s) instead of %s"
170 raise IncompatibleMatrixException(msg % (n1, m2, n1))
171 return func(*args)
172 return _decorator
173 return decorator
174
179 """ Determine whether the given arg is a matrix, which must
180 be a rectangular array of numbers, represented as a list/tuple of
181 list/tuple elements, each of which is a vector of equal size.
182 >>> is_matrix([])
183 True
184 >>> is_matrix([[]])
185 True
186 >>> is_matrix([1,2])
187 False
188 >>> is_matrix([[1,2], (3,4)])
189 True
190 >>> is_matrix([[1,2], [3,4,5]])
191 False
192 >>> is_matrix([[1,2], ['0', 1]])
193 False
194 >>> is_matrix([[.5, -1, 9], [4, 3, 19], [3, 6, 8], [0, 0, 1]])
195 True
196 """
197 if not is_vector_type(A) or not __all_true(A, is_vector_type):
198 return False
199 if len(A) == 0:
200 return True
201 n = len(A[0])
202 return __all_true(A, lambda row: is_vector(row) and len(row) == n)
203
206 """ Calculates the dimensions of the given matrix, returning a pair consisting
207 of the number of rows and the number of columns, raising an exception if any of
208 the rows of the matrix have a different number of elements.
209
210 >>> matrix_dimensions([])
211 (0, 0)
212 >>> matrix_dimensions([[1,2,3]])
213 (1, 3)
214 >>> matrix_dimensions([[1], [2], [3]])
215 (3, 1)
216 >>> matrix_dimensions([[1, 2], [1, 2, 3]])
217 Traceback (most recent call last):
218 ...
219 InvalidMatrixException: Arg is not a matrix: [[1, 2], [1, 2, 3]]
220 """
221 m = len(A)
222 if m == 0:
223 return (0,0)
224 n = len(A[0])
225 for (i, row) in enumerate(A):
226 l = len(row)
227 if n != l:
228 msg = "Matrix has rows of different length: row 0 has %d element(s) and row %d has %d element(s)." % (n, i, l)
229 raise InvalidMatrixException(msg)
230 return (m, n)
231
234 """ Determine whether the given matrix is a square matrix, with the same number of
235 rows and columns.
236
237 >>> matrix_is_square([])
238 True
239 >>> matrix_is_square([[]])
240 False
241 >>> matrix_is_square([[1,1], [2,2], [3,3]])
242 False
243 >>> matrix_is_square([[1,2,3], [4,5,6], [7,8,9]])
244 True
245 >>> matrix_is_square([[1], [1,2]])
246 Traceback (most recent call last):
247 ...
248 InvalidMatrixException: Arg is not a matrix: [[1], [1, 2]]
249 """
250 (m, n) = matrix_dimensions(A)
251 return m == n
252
255 """ Determine whether the given matrix is a diagonal matrix -- i.e., a square matrix with
256 all nondiagonal entries equal to 0.
257
258 >>> matrix_is_diagonal([])
259 True
260 >>> matrix_is_diagonal([[0]])
261 True
262 >>> matrix_is_diagonal([[0, 0], [0, 0], [0, 0]])
263 False
264 >>> matrix_is_diagonal([[1, 0], [0, 2]])
265 True
266 >>> matrix_is_diagonal([[1, 0], [1, 0]])
267 False
268 """
269 if not matrix_is_square(A):
270 return False
271 for (i, row) in enumerate(A):
272 for (j, elem) in enumerate(row):
273 if i != j and elem != 0:
274 return False
275 return True
276
279 """ Determine whether the given matrix is a scalar matrix -- i.e., a diagonal matrix
280 with all diagonal entries equal to each other.
281
282 >>> matrix_is_scalar([])
283 True
284 >>> matrix_is_scalar([[2]])
285 True
286 >>> matrix_is_scalar([[3, 0], [0, 3]])
287 True
288 >>> matrix_is_scalar([[3, 0], [1, 3]])
289 False
290 """
291 if not matrix_is_square(A):
292 return False
293 elif len(A) == 0:
294 return True
295 diag_val = A[0][0]
296 for (i, row) in enumerate(A):
297 for (j, elem) in enumerate(row):
298 if (i != j and elem != 0) or (i == j and elem != diag_val):
299 return False
300 return True
301
304 """ Determine whether the given matrix is an identity matrix -- i.e.,
305 a non-empty scalar matrix with all diagonal entries equal to 1.
306
307 >>> matrix_is_identity([])
308 False
309 >>> matrix_is_identity([[0]])
310 False
311 >>> matrix_is_identity([[1]])
312 True
313 >>> matrix_is_identity([[1, 1], [1, 1]])
314 False
315 >>> matrix_is_identity([[1, 0], [0, 1]])
316 True
317 >>> matrix_is_identity([[1, 0], [0, 1], [0, 0]])
318 False
319 >>> matrix_is_identity([[2, 0], [0, 2]])
320 False
321 """
322 if not matrix_is_square(A):
323 return False
324 elif len(A) == 0:
325 return False
326 for (i, row) in enumerate(A):
327 for (j, elem) in enumerate(row):
328 if (i != j and elem != 0) or (i == j and elem != 1):
329 return False
330 return True
331
334 """ Determine whether the given matrix is symmetric -- i.e., whether it is equal
335 to its own transpose.
336
337 >>> matrix_is_symmetric([[1, 0], [0, 1]])
338 True
339 >>> matrix_is_symmetric([[0, 1], [1, 0]])
340 True
341 >>> matrix_is_symmetric([[1, 2], [2, 0]])
342 True
343 >>> matrix_is_symmetric([[1, 2], [1, 2]])
344 False
345 >>> matrix_is_symmetric([[1, 3, 2], [3, 5, 0], [2, 0, 4]])
346 True
347 """
348 return matrix_equal(A, matrix_transpose(A))
349
351 """ Helper function for matrix_plus and matrix_minus. """
352 (m, n) = matrix_dimensions(A)
353 C = [None] * m
354 for i in range(m):
355 row = [None] * n
356 for j in range(n):
357 row[j] = func(A[i][j], B[i][j])
358 C[i] = row
359 return C
360
361 @__check_matrix_same_size(0, 1)
362 -def matrix_plus(A, B):
363 """ Add the given matrices, which must have the same number of rows and columns.
364
365 >>> matrix_plus([[1, 2], [3, 4]], [[4, 3], [2, 1]])
366 [[5, 5], [5, 5]]
367 >>> matrix_plus([[1,2]], [[3, 4]]), matrix_plus([[1,2]], [[3, 4]])[0] == vector_plus([1,2], [3,4])
368 ([[4, 6]], True)
369 >>> matrix_plus([[1]], [[1, 2]])
370 Traceback (most recent call last):
371 ...
372 IncompatibleMatrixException: Matrices must have same dimensions. Matrix 1 is 1 x 1, but matrix 2 is 1 x 2.
373 """
374 return __matrix_pairwise_op(A, B, operator.add)
375
378 """ Subtract the second matrix from the first, both of which must have the
379 same number of rows and columns.
380
381 >>> matrix_minus([[3, 4], [1, 2]], [[1, 1], [1, 1]])
382 [[2, 3], [0, 1]]
383 >>> matrix_minus([[1,2]], [[3, 4]]), matrix_minus([[1,2]], [[3, 4]])[0] == vector_minus([1,2], [3,4])
384 ([[-2, -2]], True)
385 >>> matrix_minus([[1]], [[1, 2]])
386 Traceback (most recent call last):
387 ...
388 IncompatibleMatrixException: Matrices must have same dimensions. Matrix 1 is 1 x 1, but matrix 2 is 1 x 2.
389 """
390 return __matrix_pairwise_op(A, B, operator.sub)
391
392 @__check_scalar(0)
393 @__check_matrix(1)
394 -def matrix_times(s, A):
395 """ Multiply the matrix A by the scalar s.
396
397 >>> matrix_times(2, [[0,1], [1, 0]])
398 [[0, 2], [2, 0]]
399 >>> matrix_times(-1, [[1,2,3,4]])
400 [[-1, -2, -3, -4]]
401 >>> matrix_times(3, [[]])
402 [[]]
403 >>> matrix_times(1, [])
404 []
405 >>> matrix_times(3, [[1,2,3], [4,5,6], [7,8,9], [10, 11, 12]])
406 [[3, 6, 9], [12, 15, 18], [21, 24, 27], [30, 33, 36]]
407 """
408 return [[s*elem for elem in row] for row in A]
409
412 """ Return the negative of the given matrix, or the matrix multiplied
413 by scalar -1.
414
415 >>> matrix_negative([])
416 []
417 >>> matrix_negative([[]])
418 [[]]
419 >>> matrix_negative(((-1, 2), (2, -3)))
420 [[1, -2], [-2, 3]]
421 """
422 return matrix_times(-1, A)
423
426 """ Generate the transpose of the given m x n matrix, which is the n x m
427 matrix with where each (i,j) element of the original matrix is the (j,i)
428 element of the new matrix.
429
430 >>> matrix_transpose([[1, 2], [3, 4], [5, 6]])
431 [[1, 3, 5], [2, 4, 6]]
432 >>> matrix_transpose([])
433 []
434 >>> matrix_transpose([[]])
435 []
436 >>> matrix_transpose([[1,2,3], [4,5,6], [7,8,9]])
437 [[1, 4, 7], [2, 5, 8], [3, 6, 9]]
438 >>> matrix_transpose(matrix_transpose([[1,2,3], [4,5,6], [7,8,9]]))
439 [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
440 >>> matrix_transpose([[1,2,3]])
441 [[1], [2], [3]]
442 """
443 (m, n) = matrix_dimensions(A)
444 C = [[None] * m for i in range(n)]
445 for i in range(m):
446 for j in range(n):
447 C[j][i] = A[i][j]
448 return C
449
450 @__check_matrix_multipliable(0, 1)
451 @__check_matrix(0, 1)
452 -def matrix_product(A, B):
453 """ Calculate the product of matrices A and B, where A must have the
454 same number of columns and B has rows. The result for A of size m x n
455 and B of size n x r is an m x r matrix.
456
457 >>> matrix_product([[1,1], [1,1]], [[1,1], [1,1]])
458 [[2, 2], [2, 2]]
459 >>> matrix_product([], [])
460 []
461 >>> matrix_product([[1, 3, -1], [-2, -1, 1]], [[-4, 0, 3, -1], [5, -2, -1, 1], [-1, 2, 0, 6]])
462 [[12, -8, 0, -4], [2, 4, -5, 7]]
463 """
464 (m, n) = matrix_dimensions(A)
465 (n, r) = matrix_dimensions(B)
466 C = [[None] * r for i in range(m)]
467 for i in range(m):
468 for j in range(r):
469 C[i][j] = vector_product(A[i], [row[j] for row in B])
470 return C
471
473 """ Determine whether the two scalars are equal to accuracy of epsilon.
474
475 >>> scalar_equal(0.9000000000000002, 0.9000000000000001)
476 True
477 >>> scalar_equal(1, 1)
478 True
479 >>> scalar_equal(1, 1.0)
480 True
481 >>> scalar_equal(1.0, 1.00000001)
482 False
483 """
484 return abs(s1-s2) < eps
485
488 """ Determine whether the two vectors are equal by determining if the
489 pairwise components are equal to accuracy of epsilon.
490
491 >>> vector_equal([2], [2])
492 True
493 >>> vector_equal([1.0, 0.9000000000000002], [1, 0.9000000000000001])
494 True
495 """
496 return __all_true(map(lambda (x,y): scalar_equal(x,y), zip(u, v)))
497
500 """ Determine whether the two matrices are equal by determining if
501 the corresponding components are equal to accuracy of epsilon.
502
503 >>> matrix_equal([], [])
504 True
505 >>> matrix_equal([[1, 2, 3], [3.9000000000000001, 2, 1]], [[1, 2, 3], [3.9, 2, 1]])
506 True
507 """
508 (m1, n1) = matrix_dimensions(A)
509 (m2, n2) = matrix_dimensions(B)
510 if m1 != m2 or n1 != n2:
511 return False
512 for i in range(m1):
513 for j in range(n1):
514 if not scalar_equal(A[i][j], B[i][j]):
515 return False
516 return True
517
519 """ Create a zero vector with same size as the given vector, or
520 having as many elements as the given int.
521
522 >>> vector_zero(3)
523 [0, 0, 0]
524 >>> vector_zero([1,2,3,4])
525 [0, 0, 0, 0]
526 >>> vector_zero([])
527 []
528 """
529 if isinstance(x, int):
530 return [0] * x
531 elif not isinstance(x, list) and not isinstance(x, tuple):
532 raise LAValueError("Arg should be a vector or int, not: %s" % repr(x))
533 else:
534 return [0 for n in x]
535
538 """ Create a unit-length vector in direction of v.
539
540 >>> vector_unit_length([3,3])
541 [0.70710678118654757, 0.70710678118654757]
542 """
543 return vector_times(1.0/vector_norm(v), v)
544
547 """ Multiplies the vector by the given scalar.
548
549 >>> vector_times(3, [1,2])
550 [3, 6]
551 >>> vector_times(5, [4,-2,7])
552 [20, -10, 35]
553 >>> vector_times(0, [1,2,3])
554 [0, 0, 0]
555 """
556 return map(lambda x: s * x, v)
557
559 """ A helper function that applies the given 2-arg func to
560 the pairwise elements of the two vectors. """
561 return map(lambda (x, y): func(x, y), zip(u, v))
562
563 @__check_vector_same_size(0, 1)
564 -def vector_plus(u, v):
565 """ Calculate the vector resulting from the addition of
566 the two given vectors.
567
568 >>> vector_plus([1,2], [3,4])
569 [4, 6]
570 >>> vector_plus([0.5, 1], [2.5, -1])
571 [3.0, 0]
572 """
573 return __vector_pairwise_op(u, v, operator.add)
574
577 """ Subtract the second vector from the first vector.
578
579 >>> vector_minus([2], [-3])
580 [5]
581 >>> vector_minus([1,2,3], [2,3,4])
582 [-1, -1, -1]
583 """
584 return __vector_pairwise_op(u, v, operator.sub)
585
588 """ Calculate the dot product of the given vectors, which is defined
589 as the sum of the pairwise products of the vectors.
590
591 >>> vector_product([0,0], [0,1])
592 0
593 >>> vector_product([1,2,-3], [-3,5,2])
594 1
595 """
596 return sum(__vector_pairwise_op(u, v, operator.mul))
597
600 """ Calculate the length, or norm, of the vector, which
601 is defined to be the square root of the sum of the squares
602 of the components of the vector.
603
604 >>> vector_norm([0,1])
605 1.0
606 >>> vector_norm([1,1]), scalar_equal(vector_norm([1,1]), sqrt(2))
607 (1.4142135623730951, True)
608 >>> vector_norm([2,-1,3]), scalar_equal(vector_norm([2,-1,3]), sqrt(14))
609 (3.7416573867739413, True)
610 >>> vector_norm([sqrt(2), -1, 1])
611 2.0
612 """
613 return math.sqrt(sum(map(lambda x: x * x, v)))
614
617 """ Calculate the distance between the two vectors, which
618 is defined to be the norm of the difference between the
619 vectors.
620 >>> vector_distance([sqrt(2), 1, -1], [0, 2, -2])
621 2.0
622 >>> vector_distance([0,1], [1,0])
623 1.4142135623730951
624 >>> vector_distance([-1,2], [3,1])
625 4.1231056256176606
626 """
627 return vector_norm(vector_minus(u, v))
628
631 """ Determine whether the two vectors are orthogonal, which
632 is defined by whether the dot product is 0.
633
634 >>> vector_is_orthogonal([0,1], [1,0])
635 True
636 >>> vector_is_orthogonal([2,-1,1], [1,-2,-1])
637 False
638 >>> vector_is_orthogonal([1,1,-2], [3,1,2])
639 True
640 """
641 return abs(vector_product(u, v)) < epsilon
642
645 """ Determine the angle between the two vectors, in radians.
646
647 >>> vector_angle([0,1,1], [1,0,1]), scalar_equal(vector_angle([0,1,1], [1,0,1]), pi/3)
648 (1.0471975511965979, True)
649 >>> vector_angle([0.9, 2.1, 1.2], [-4.5, 2.6, -0.8])
650 1.5376292299388497
651 """
652 return math.acos(vector_product(u, v) / (vector_norm(u) * vector_norm(v)))
653
656 """ Calculate the projection of v onto u.
657
658 >>> vector_project([2,1], [-1, 3]), vector_equal(vector_project([2,1], [-1, 3]), [2.0/5, 1.0/5])
659 ([0.40000000000000002, 0.20000000000000001], True)
660 >>> vector_project([0,0,1], [1,2,3]), vector_equal(vector_project( [0,0,1], [1,2,3]), [0,0,3])
661 ([0.0, 0.0, 3.0], True)
662 """
663 return vector_times(float(vector_product(u, v)) / vector_product(u, u), u)
664
669 """ Determine whether arg is a vector type (list or tuple). """
670 return isinstance(u, list) or isinstance(u, tuple)
671
673 """ Determine whether the arg is a vector, or a list or tuple of numbers.
674
675 >>> is_vector('')
676 False
677 >>> is_vector([])
678 True
679 >>> is_vector(())
680 True
681 >>> is_vector((1,))
682 True
683 >>> is_vector([1])
684 True
685 >>> is_vector([1,2,3,4])
686 True
687 >>> is_vector((1,2,3))
688 True
689 >>> is_vector([1, 0.5, -9L])
690 True
691 >>> is_vector([1, 2, "3"])
692 False
693 >>> is_vector([[0]])
694 False
695 >>> is_vector([[1,2]])
696 False
697 """
698 return is_vector_type(u) and __all_true(u, is_scalar)
699
701 """ Determine whether s is a scalar, which is either a float, int,
702 or long. The integral values 0 and 1 may also be given as
703 False and True, respectively.
704
705 >>> is_scalar(0)
706 True
707 >>> is_scalar(True)
708 True
709 >>> is_scalar(-9.2)
710 True
711 >>> is_scalar(100L)
712 True
713 >>> is_scalar((1,2))
714 False
715 """
716 return isinstance(s, int) or isinstance(s, float) or isinstance(s, long)
717
719 """ Determine if all items in list evaluate to True using the optional
720 test function, which defaults to bool.
721
722 >>> __all_true([])
723 True
724 >>> __all_true([1,0])
725 False
726 >>> __all_true([True, 1])
727 True
728 >>> __all_true([True, True, False])
729 False
730 >>> __all_true([1,2,3], lambda s: s % 2 == 0)
731 False
732 >>> __all_true([2,4,6], lambda s: s % 2 == 0)
733 True
734 """
735 for res in lst:
736 if not test(res): return False
737 return True
738
740 import doctest
741 doctest.testmod()
742
743 if __name__ == "__main__":
744 __test()
745