NAME

math::linearalgebra - Linear Algebra

SYNOPSIS

package require Tcl ?8.4?

package require math::linearalgebra ?1.1?

::math::linearalgebra::mkVector ndim value

::math::linearalgebra::mkUnitVector ndim ndir

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

::math::linearalgebra::getrow matrix row ?imin? ?imax?

::math::linearalgebra::setrow matrix row newvalues ?imin? ?imax?

::math::linearalgebra::getcol matrix col ?imin? ?imax?

::math::linearalgebra::setcol matrix col newvalues ?imin? ?imax?

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

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

::math::linearalgebra::swaprows matrix irow1 irow2 ?imin? ?imax?

::math::linearalgebra::swapcols matrix icol1 icol2 ?imin? ?imax?

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

::math::linearalgebra::dim obj

::math::linearalgebra::shape obj

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

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

::math::linearalgebra::norm vector type

::math::linearalgebra::norm_one vector

::math::linearalgebra::norm_two vector

::math::linearalgebra::norm_max vector ?index?

::math::linearalgebra::normMatrix matrix type

::math::linearalgebra::dotproduct vect1 vect2

::math::linearalgebra::unitLengthVector vector

::math::linearalgebra::normalizeStat mv

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

::math::linearalgebra::add mv1 mv2

::math::linearalgebra::sub mv1 mv2

::math::linearalgebra::scale scale mv

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

::math::linearalgebra::transpose matrix

::math::linearalgebra::matmul mv1 mv2

::math::linearalgebra::angle vect1 vect2

::math::linearalgebra::crossproduct vect1 vect2

::math::linearalgebra::matmul mv1 mv2

::math::linearalgebra::mkIdentity size

::math::linearalgebra::mkDiagonal diag

::math::linearalgebra::mkRandom size

::math::linearalgebra::mkTriangular size ?uplo? ?value?

::math::linearalgebra::mkHilbert size

::math::linearalgebra::mkDingdong size

::math::linearalgebra::mkOnes size

::math::linearalgebra::mkMoler size

::math::linearalgebra::mkFrank size

::math::linearalgebra::mkBorder size

::math::linearalgebra::mkWilkinsonW+ size

::math::linearalgebra::mkWilkinsonW- size

::math::linearalgebra::solveGauss matrix bvect

::math::linearalgebra::solvePGauss matrix bvect

::math::linearalgebra::solveTriangular matrix bvect ?uplo?

::math::linearalgebra::solveGaussBand matrix bvect

::math::linearalgebra::solveTriangularBand matrix bvect

::math::linearalgebra::determineSVD A eps

::math::linearalgebra::eigenvectorsSVD A eps

::math::linearalgebra::leastSquaresSVD A y qmin eps

::math::linearalgebra::choleski matrix

::math::linearalgebra::orthonormalizeColumns matrix

::math::linearalgebra::orthonormalizeRows matrix

::math::linearalgebra::dger matrix alpha x y ?scope?

::math::linearalgebra::dgetrf matrix

::math::linearalgebra::det matrix

::math::linearalgebra::largesteigen matrix tolerance maxiter

::math::linearalgebra::to_LA mv

::math::linearalgebra::from_LA mv

DESCRIPTION

This package offers both low-level procedures and high-level algorithms to deal with linear algebra problems:

• 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.

Note: When using this package in combination with Tk, there may be a naming conflict, as both this package and Tk define a command scale. See the NAMING CONFLICT section below.

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.

::math::linearalgebra::mkUnitVector ndim ndir

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

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

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

::math::linearalgebra::getrow matrix row ?imin? ?imax?

Returns a single row of a matrix as a list

::math::linearalgebra::setrow matrix row newvalues ?imin? ?imax?

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)

::math::linearalgebra::getcol matrix col ?imin? ?imax?

Returns a single column of a matrix as a list

::math::linearalgebra::setcol matrix col newvalues ?imin? ?imax?

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)

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

Returns a single element of a matrix/vector

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

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

::math::linearalgebra::swaprows matrix irow1 irow2 ?imin? ?imax?

Swap two rows in a matrix completely or only a selected part

::math::linearalgebra::swapcols matrix icol1 icol2 ?imin? ?imax?

Swap two columns in a matrix completely or only a selected part

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)

::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)

::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)

::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.

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

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

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.

::math::linearalgebra::norm_one vector

Returns the L1 norm of the given vector, the sum of absolute values

::math::linearalgebra::norm_two vector

Returns the L2 norm of the given vector, the ordinary Euclidean norm

::math::linearalgebra::norm_max vector ?index?

Returns the Linf norm of the given vector, the maximum absolute coefficient

::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.

::math::linearalgebra::dotproduct vect1 vect2

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

::math::linearalgebra::unitLengthVector vector

Return a vector in the same direction with length 1.

::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.

::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.

Specialised variants are: axpy_vect and axpy_mat (slightly faster, but no check on the arguments)

::math::linearalgebra::add mv1 mv2

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

Specialised variants are: add_vect and add_mat (slightly faster, but no check on the arguments)

::math::linearalgebra::sub mv1 mv2

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

Specialised variants are: sub_vect and sub_mat (slightly faster, but no check on the arguments)

::math::linearalgebra::scale scale mv

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

Specialised variants are: scale_vect and scale_mat (slightly faster, but no check on the arguments)

::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.

::math::linearalgebra::transpose matrix

Transpose a matrix

::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.

::math::linearalgebra::angle vect1 vect2

Compute the angle between two vectors (in radians)

::math::linearalgebra::crossproduct vect1 vect2

Compute the cross product of two (three-dimensional) vectors

::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.

Common matrices and test matrices

::math::linearalgebra::mkIdentity size

Create an identity matrix of dimension size.

::math::linearalgebra::mkDiagonal diag

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

::math::linearalgebra::mkRandom size

Create a square matrix whose elements are uniformly distributed random numbers between 0 and 1 of dimension size.

::math::linearalgebra::mkTriangular size ?uplo? ?value?

Create a triangular matrix with non-zero elements in the upper or lower part, depending on argument uplo.

::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.

::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.

::math::linearalgebra::mkOnes size

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

::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.)

::math::linearalgebra::mkFrank size

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

::math::linearalgebra::mkBorder size

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

::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.

::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.

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.

::math::linearalgebra::solvePGauss matrix bvect

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

::math::linearalgebra::solveTriangular matrix bvect ?uplo?

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

::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.

::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.

::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.

::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.

::math::linearalgebra::leastSquaresSVD A y qmin eps

Determines the solution to a least-sqaures problem Ax ~ y via singular value decomposition. The result is the vector x.

Note that if you add a column of 1s to the matrix, then this column will represent a constant like in: y = a*x1 + b*x2 + c. To force the intercept to be zero, simply leave it out.

::math::linearalgebra::choleski matrix

Determine the Choleski decomposition of a symmetric positive semidefinite matrix (this condition is not checked!). The result is the lower-triangular matrix L such that L Lt = matrix.

::math::linearalgebra::orthonormalizeColumns matrix

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

::math::linearalgebra::orthonormalizeRows matrix

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

::math::linearalgebra::dger matrix alpha x y ?scope?

Perform the rank 1 operation A + alpha*x*y' inline (that is: the matrix A is adjusted). For convenience the new matrix is also returned as the result.

::math::linearalgebra::dgetrf matrix

Computes an LU factorization of a general matrix, using partial, pivoting with row interchanges. Returns the permutation vector.

The factorization has the form

   P * A = L * U
    

where P is a permutation matrix, L is lower triangular with unit diagonal elements, and U is upper triangular. Returns the permutation vector, as a list of length n-1. The last entry of the permutation is not stored, since it is implicitely known, with value n (the last row is not swapped with any other row). At index #i of the permutation is stored the index of the row #j which is swapped with row #i at step #i. That means that each index of the permutation gives the permutation at each step, not the cumulated permutation matrix, which is the product of permutations.

::math::linearalgebra::det matrix

Returns the determinant of the given matrix, based on PA=LU decomposition, i.e. Gauss partial pivotal.

::math::linearalgebra::largesteigen matrix tolerance maxiter

Returns a list made of the largest eigenvalue (in magnitude) and associated eigenvector. Uses iterative Power Method as provided as algorithm #7.3.3 of Golub & Van Loan. This algorithm is used here for a dense matrix (but is usually used for sparse matrices).

Compability with the LA package Two procedures are provided for compatibility with Hume's LA package:

::math::linearalgebra::to_LA mv

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

::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.

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:

• determineQR
• certainlyPositive, diagonallyDominant

NAMING CONFLICT

If you load this package in a Tk-enabled shell like wish, then the command

namespace import ::math::linearalgebra

results in an error message about "scale". This is due to the fact that Tk defines all its commands in the global namespace. The solution is to import the linear algebra commands in a namespace that is not the global one:

package require math::linearalgebra
namespace eval compute {


    namespace import ::math::linearalgebra::*


    ... use the linear algebra version of scale ...
}

To use Tk's scale command in that same namespace you can rename it:

namespace eval compute {


    rename ::scale scaleTk


    scaleTk .scale ...
}

BUGS, IDEAS, FEEDBACK

This document, and the package it describes, will undoubtedly contain bugs and other problems. Please report such in the category math :: linearalgebra of the Tcllib SF Trackers [http://sourceforge.net/tracker/?group_id=12883]. Please also report any ideas for enhancements you may have for either package and/or documentation.

KEYWORDS

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

CATEGORY

Mathematics

COPYRIGHT

Copyright (c) 2004-2008 Arjen Markus <arjenmarkus@users.sourceforge.net>
Copyright (c) 2004 Ed Hume <http://www.hume.com/contact.us.htm>
Copyright (c) 2008 Michael Buadin <relaxkmike@users.sourceforge.net>