robust solution of linear equations or least squares problems
determining eigenvectors and eigenvalues of symmetric matrices
various decompositions of general matrices or matrices of a specific form
(limited) support for matrices in band storage, a common type of sparse matrices

It arose as a re-implementation of Hume's LA package and the desire to offer low-level procedures as found in the well-known BLAS library. Matrices are implemented as lists of lists rather linear lists with reserved elements, as in the original LA package, as it was found that such an implementation is actually faster.

It is advisable, however, to use the procedures that are offered, such as setrow and getrow, rather than rely on this representation explicitly: that way it is to switch to a possibly even faster compiled implementation that supports the same API.

PROCEDURES

The package defines the following public procedures (several exist as specialised procedures, see below):

Constructing matrices and vectors

::math::linearalgebra::mkVector ndim value

Create a vector with ndim elements, each with the value value.

integer ndim: Dimension of the vector (number of components)
double value: Uniform value to be used (default: 0.0)

::math::linearalgebra::mkUnitVector ndim ndir

Create a unit vector in ndim-dimensional space, along the ndir-th direction.

integer ndim: Dimension of the vector (number of components)
integer ndir: Direction (0, ..., ndim-1)

::math::linearalgebra::mkMatrix nrows ncols value

Create a matrix with nrows rows and ncols columns. All elements have the value value.

integer nrows: Number of rows
integer ncols: Number of columns
double value: Uniform value to be used (default: 0.0)

::math::linearalgebra::getrow matrix row

Returns a single row of a matrix as a list

list matrix: Matrix in question
int row: Index of the row to return

::math::linearalgebra::setrow matrix row newvalues

Set a single row of a matrix to new values (this list must have the same number of elements as the number of columns in the matrix)

list matrix: name of the matrix in question
int row: Index of the row to update
list newvalues: List of new values for the row

::math::linearalgebra::getcol matrix col

Returns a single column of a matrix as a list

list matrix: Matrix in question
int col: Index of the column to return

::math::linearalgebra::setcol matrix col newvalues

Set a single column of a matrix to new values (this list must have the same number of elements as the number of rows in the matrix)

list matrix: name of the matrix in question
int col: Index of the column to update
list newvalues: List of new values for the column

::math::linearalgebra::getelem matrix row col

Returns a single element of a matrix/vector

list matrix: Matrix or vector in question
int row: Row of the element
int col: Column of the element (not present for vectors)

::math::linearalgebra::setelem matrix row ?col? newvalue

Set a single element of a matrix (or vector) to a new value

list matrix: name of the matrix in question
int row: Row of the element
int col: Column of the element (not present for vectors)

Querying matrices and vectors

::math::linearalgebra::show obj ?format? ?rowsep? ?colsep?

Return a string representing the vector or matrix, for easy printing. (There is currently no way to print fixed sets of columns)

list obj: Matrix or vector in question
string format: Format for printing the numbers (default: %6.4f)
string rowsep: String to use for separating rows (default: newline)
string colsep: String to use for separating columns (default: space)

::math::linearalgebra::dim obj

Returns the number of dimensions for the object (either 0 for a scalar, 1 for a vector and 2 for a matrix)

any obj: Scalar, vector, or matrix

::math::linearalgebra::shape obj

Returns the number of elements in each dimension for the object (either an empty list for a scalar, a single number for a vector and a list of the number of rows and columns for a matrix)

any obj: Scalar, vector, or matrix

::math::linearalgebra::conforming type obj1 obj2

Checks if two objects (vector or matrix) have conforming shapes, that is if they can be applied in an operation like addition or matrix multiplication.

string type

Type of check:

"shape" - the two objects have the same shape (for all element-wise operations)
"rows" - the two objects have the same number of rows (for use as A and b in a system of linear equations Ax = b
"matmul" - the first object has the same number of columns as the number of rows of the second object. Useful for matrix-matrix or matrix-vector multiplication.

list obj1

First vector or matrix (left operand)

list obj2

Second vector or matrix (right operand)

::math::linearalgebra::symmetric matrix ?eps?

Checks if the given (square) matrix is symmetric. The argument eps is the tolerance.

list matrix: Matrix to be inspected
float eps: Tolerance for determining approximate equality (defaults to 1.0e-8)

Basic operations

::math::linearalgebra::norm vector type

Returns the norm of the given vector. The type argument can be: 1, 2, inf or max, respectively the sum of absolute values, the ordinary Euclidean norm or the max norm.

list vector: Vector, list of coefficients
string type: Type of norm (default: 2, the Euclidean norm)

::math::linearalgebra::normMatrix matrix type

Returns the norm of the given matrix. The type argument can be: 1, 2, inf or max, respectively the sum of absolute values, the ordinary Euclidean norm or the max norm.

list matrix: Matrix, list of row vectors
string type: Type of norm (default: 2, the Euclidean norm)

::math::linearalgebra::dotproduct vect1 vect2

Determine the inproduct or dot product of two vectors. These must have the same shape (number of dimensions)

list vect1: First vector, list of coefficients
list vect2: Second vector, list of coefficients

::math::linearalgebra::unitLengthVector vector

Return a vector in the same direction with length 1.

list vector: Vector to be normalized

::math::linearalgebra::normalizeStat mv

Normalize the matrix or vector in a statistical sense: the mean of the elements of the columns of the result is zero and the standard deviation is 1.

list mv: Vector or matrix to be normalized in the above sense

::math::linearalgebra::axpy scale mv1 mv2

Return a vector or matrix that results from a "daxpy" operation, that is: compute a*x+y (a a scalar and x and y both vectors or matrices of the same shape) and return the result.

double scale: The scale factor for the first vector/matrix (a)
list mv1: First vector or matrix (x)
list mv2: Second vector or matrix (y)

::math::linearalgebra::add mv1 mv2

Return a vector or matrix that is the sum of the two arguments (x+y)

list mv1: First vector or matrix (x)
list mv2: Second vector or matrix (y)

::math::linearalgebra::sub mv1 mv2

Return a vector or matrix that is the difference of the two arguments (x-y)

list mv1: First vector or matrix (x)
list mv2: Second vector or matrix (y)

::math::linearalgebra::scale scale mv

Scale a vector or matrix and return the result, that is: compute a*x.

double scale: The scale factor for the vector/matrix (a)
list mv: Vector or matrix (x)

::math::linearalgebra::rotate c s vect1 vect2

Apply a planar rotation to two vectors and return the result as a list of two vectors: c*x-s*y and s*x+c*y. In algorithms you can often easily determine the cosine and sine of the angle, so it is more efficient to pass that information directly.

double c: The cosine of the angle
double s: The sine of the angle
list vect1: First vector (x)
list vect2: Seocnd vector (x)

::math::linearalgebra::transpose matrix

Transpose a matrix

list matrix: Matrix to be transposed

::math::linearalgebra::matmul mv1 mv2

Multiply a vector/matrix with another vector/matrix. The result is a matrix, if both x and y are matrices or both are vectors, in which case the "outer product" is computed. If one is a vector and the other is a matrix, then the result is a vector.

list mv1: First vector/matrix (x)
list mv2: Second vector/matrix (y)

::math::linearalgebra::angle vect1 vect2

Compute the angle between two vectors (in radians)

list vect1: First vector
list vect2: Second vector

::math::linearalgebra::matmul mv1 mv2

list mv1: First vector/matrix (x)
list mv2: Second vector/matrix (y)

Common matrices and test matrices

::math::linearalgebra::mkIdentity size

Create an identity matrix of dimension size.

integer size: Dimension of the matrix

::math::linearalgebra::mkDiagonal diag

Create a diagonal matrix whose diagonal elements are the elements of the vector diag.

list diag: Vector whose elements are used for the diagonal

::math::linearalgebra::mkHilbert size

Create a Hilbert matrix of dimension size. Hilbert matrices are very ill-conditioned with respect to eigenvalue/eigenvector problems. Therefore they are good candidates for testing the accuracy of algorithms and implementations.

integer size: Dimension of the matrix

::math::linearalgebra::mkDingdong size

Create a "dingdong" matrix of dimension size. Dingdong matrices are imprecisely represented, but have the property of being very stable in such algorithms as Gauss elimination.

integer size: Dimension of the matrix

::math::linearalgebra::mkOnes size

Create a square matrix of dimension size whose entries are all 1.

integer size: Dimension of the matrix

::math::linearalgebra::mkMoler size

Create a Moler matrix of size size. (Moler matrices have a very simple Choleski decomposition. It has one small eigenvalue and it can easily upset elimination methods for systems of linear equations.)

integer size: Dimension of the matrix

::math::linearalgebra::mkFrank size

Create a Frank matrix of size size. (Frank matrices are fairly well-behaved matrices)

integer size: Dimension of the matrix

::math::linearalgebra::mkBorder size

Create a bordered matrix of size size. (Bordered matrices have a very low rank and can upset certain specialised algorithms.)

integer size: Dimension of the matrix

::math::linearalgebra::mkWilkinsonW+ size

Create a Wilkinson W+ of size size. This kind of matrix has pairs of eigenvalues that are very close together. Usually the order (size) is odd.

integer size: Dimension of the matrix

::math::linearalgebra::mkWilkinsonW- size

Create a Wilkinson W- of size size. This kind of matrix has pairs of eigenvalues with opposite signs, when the order (size) is odd.

integer size: Dimension of the matrix

Common algorithms

::math::linearalgebra::solveGauss matrix bvect

Solve a system of linear equations (Ax=b) using Gauss elimination. Returns the solution (x) as a vector or matrix of the same shape as bvect.

list matrix: Square matrix (matrix A)
list bvect: Vector or matrix whose columns are the individual b-vectors

::math::linearalgebra::solveTriangular matrix bvect

Solve a system of linear equations (Ax=b) by backward substitution. The matrix is supposed to be upper-triangular.

list matrix: Upper-triangular matrix (matrix A)
list bvect: Vector or matrix whose columns are the individual b-vectors

::math::linearalgebra::solveGaussBand matrix bvect

Solve a system of linear equations (Ax=b) using Gauss elimination, where the matrix is stored as a band matrix (cf. STORAGE). Returns the solution (x) as a vector or matrix of the same shape as bvect.

list matrix: Square matrix (matrix A; in band form)
list bvect: Vector or matrix whose columns are the individual b-vectors

::math::linearalgebra::solveTriangularBand matrix bvect

Solve a system of linear equations (Ax=b) by backward substitution. The matrix is supposed to be upper-triangular and stored in band form.

list matrix: Upper-triangular matrix (matrix A)
list bvect: Vector or matrix whose columns are the individual b-vectors

::math::linearalgebra::determineSVD A eps

Determines the Singular Value Decomposition of a matrix: A = U S Vtrans. Returns a list with the matrix U, the vector of singular values S and the matrix V.

list A: Matrix to be decomposed
float eps: Tolerance (defaults to 2.3e-16)

::math::linearalgebra::eigenvectorsSVD A eps

Determines the eigenvectors and eigenvalues of a real symmetric matrix, using SVD. Returns a list with the matrix of normalized eigenvectors and their eigenvalues.

list A: Matrix whose eigenvalues must be determined
float eps: Tolerance (defaults to 2.3e-16)

::math::linearalgebra::orthonormalizeColumns matrix

Use the modified Gram-Schmidt method to orthogonalize and normalize the columns of the given matrix and return the result.

list matrix: Matrix whose columns must be orthonormalized

::math::linearalgebra::orthonormalizeRows matrix

Use the modified Gram-Schmidt method to orthogonalize and normalize the rows of the given matrix and return the result.

list matrix: Matrix whose rows must be orthonormalized

Compability with the LA package

::math::linearalgebra::to_LA mv

Transforms a vector or matrix into the format used by the original LA package.

list mv: Matrix or vector

::math::linearalgebra::from_LA mv

Transforms a vector or matrix from the format used by the original LA package into the format used by the present implementation.

list mv: Matrix or vector as used by the LA package

STORAGE

While most procedures assume that the matrices are given in full form, the procedures solveGaussBand and solveTriangularBand assume that the matrices are stored as band matrices. This common type of "sparse" matrices is related to ordinary matrices as follows:

"A" is a full-size matrix with N rows and M columns.
"B" is a band matrix, with m upper and lower diagonals and n rows.
"B" can be stored in an ordinary matrix of (2m+1) columns (one for each off-diagonal and the main diagonal) and n rows.
Element i,j (i = -m,...,m; j =1,...,n) of "B" corresponds to element k,j of "A" where k = M+i-1 and M is at least (!) n, the number of rows in "B".

To set element (i,j) of matrix "B" use:

    setelem B $j [expr {$N+$i-1}] $value

(There is no convenience procedure for this yet)

REMARKS ON THE IMPLEMENTATION

There is a difference between the original LA package by Hume and the current implementation. Whereas the LA package uses a linear list, the current package uses lists of lists to represent matrices. It turns out that with this representation, the algorithms are faster and easier to implement.

The LA package was used as a model and in fact the implementation of, for instance, the SVD algorithm was taken from that package. The set of procedures was expanded using ideas from the well-known BLAS library and some algorithms were updated from the second edition of J.C. Nash's book, Compact Numerical Methods for Computers, (Adam Hilger, 1990) that inspired the LA package.

Two procedures are provided to make the transition between the two implementations easier: to_LA and from_LA. They are described above.

TODO

Odds and ends: the following algorithms have not been implemented yet:

normalizeStat (!)
determineQR, cholesky
leastSquaresSVD, leastSquares
certainlyPositive, diagonallyDominant

KEYWORDS

least squares , linear algebra , linear equations , math , matrices , vectors

math::linearalgebra(n) 1.0 "Tcl Math Library"

Copyright © 2005 for compilation: ActiveState