Module linearalgebra
[frames] | no frames]

Source Code for Module linearalgebra

  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  # We re-export these from the math module for convenience. 
 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 
52 53 # The value we use for determining whether a float is equal to an 54 # int: if the difference is less than epsilon, we deem them equal. 55 epsilon = 0.000000000000001 56 57 -class LAValueError(ValueError):
58 pass
59 -class InvalidVectorException(LAValueError):
60 pass
61 -class IncompatibleVectorException(LAValueError):
62 pass
63 -class InvalidMatrixException(LAValueError):
64 pass
65 -class IncompatibleMatrixException(LAValueError):
66 pass
67
68 ## First, some decorators that are used for verifying arguments, such as 69 ## that a vector is a vector or that two matrices have dimensions that 70 ## allow them to be multiplied. 71 72 # This decorator verifies args to the other decorators, all of which 73 # take integer arguments indicating the index of the arg to be checked 74 # before the function that is decorated is called. 75 -def __intArgs(func):
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
88 # Verify that two vectors have the same size. 89 @__intArgs 90 -def __samesizeV(index1, index2):
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
116 # Verify that two matrices have the same number of rows and columns. 117 @__intArgs 118 -def __samesizeM(index1, index2):
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 # Verify that one or more args of a function is a scalar. 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 # Verify that one or more args is a vector. 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 # Verify that one or more args is a matrix. 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
191 # Verify that the number of columns of the first matrix equals number of rows 192 # of the second matrix. 193 @__intArgs 194 -def __columnRowCompatibleM(index1, index2):
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
215 ## Actual functions defined below here. 216 217 -def isM(A):
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
254 @__verifyM(0) 255 -def dimensionsM(A):
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
282 @__verifyM(0) 283 -def isSquareM(A):
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
303 @__verifyM(0) 304 -def isDiagonalM(A):
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
327 @__verifyM(0) 328 -def isScalarM(A):
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
352 @__verifyM(0) 353 -def isIdentityM(A):
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
382 @__verifyM(0) 383 -def isSymmetricM(A):
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
400 -def __pairwiseM(A, B, func):
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
462 @__verifyM(0) 463 -def negativeM(A):
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
476 @__verifyM(0) 477 -def transposeM(A):
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
524 -def areEqualS(s1, s2):
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
534 @__samesizeV(0, 1) 535 -def areEqualV(u, v):
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
546 @__verifyM(0, 1) 547 -def areEqualM(A, B):
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)