LinearAlgebra Module Documentation

(typeset 22 October 1996)

Just like the FFT module is a wrapper around some low-level FFTPACK routines, the LinearAlgebra module is a wrapper around LAPACK routines (or lapack_lite, if LAPACK was not installed on your system).

Solving Linear Equations

The function solve_linear_equations(a,b) ... [XXX]

Inverting Matrices

inverse(a) returns the inverse of a if possible (if the array a is singular, inverse will (rightly) complain. If a is not singular, then matrixmultiply(a, inverse(a)) will be very close to identity(len(a)).

Finding the Eigenvalues and Eigenvectors of a Matrix

The eigenvalues and eigenvectors of a matrix A are such that for each "matched" pair of eigenvalue u and eigenvector v,

	A.matrixmultiply(v) = u*v

If you only need to find the eigenvalues, then you should use the call eigenvalues(a), which returns an array of the eigenvalues. If you need both the eigenvalues and the eigenvectors, call eigenvectors(a), which returns a 2-tuple (U,V), where the first element of the tuple is the array of eigenvalues and the second is the matrix of eigenvectors.

	>>> a
	 2  3  4  5  6
	 7  8  9 10 11
	12 13 14 15 16
	17 18 19 20 21
	22 23 24 25 26
	>>> evals, evecs = eigenvectors(a)
	>>> for i in range(len(evals)):
	...    eval = evals[i]
	...    evec = evecs[i]
	...    print cumsum(matrixmultiply(a, evec) - eval*evec)
	...
	  7.99360578e-15   1.93178806e-14   2.32591724e-14   2.20379270e-14   2.20379270e-14
	 -3.55271368e-15  -3.55271368e-14  -4.26325641e-14  -3.55271368e-14  -4.26325641e-14
	  4.62284525e-16   4.31902095e-15   8.39325454e-15   1.01401235e-14   1.44328993e-14
	  4.22433695e-16   7.21922331e-16   1.80510920e-15   1.78632987e-15   2.67147415e-15
	  1.12863693e-16   7.36033120e-16   1.51145339e-15   3.39893173e-15   5.13478149e-15

Notice that these are all very small numbers -- proving that all of the solutions were good

Singular Value Decomposition

The function singular_value_decomposition(a) returns, as its name implies, the singular value decomposition of the matrix it gets as an argument. [XXX]

Generalized Inverse

The function generalized_inverse(a) returns, as its name implies, the generalized inverse of the matrix it gets as an argument. [XXX]

Determinant

The determinant of a matrix is given by the function determinant(a), which returns the cumulative product of all the eigenvalues of the matrix:

	>>> a = array([[1,2],[3,4]])
	>>> a                    # for a matrix this simple, it's easy to check:
	1 2                      # det(a) = 1*4 - 3*2 = -2
	3 4
	>>> determinant(a)
	-2.0                     # voila!

Linear Least Squares

The function linear_least_squares(a, b, rcond=1.e-10 takes two arrays a and band a "small number". It returns a 4-tuple (x, residuals, rank, s) such that:

x minimizes 2-norm(|b - Ax|), residuals is the sum square residuals, rank is the rank of the array and s is the rank of the singular values of A in descending order.

If b is a matrix then x is also a matrix with corresponding columns. If the rank of A is less than the number of columns of A or greater than the number of rows, then residuals will be returned as an empty array otherwise resids = sum((b-dot(A,x)**2). Singular values less than s[0]*rcond are treated as zero.