1math::linearalgebra(n) Tcl Math Library math::linearalgebra(n)
2______________________________________________________________________________
6
8 math::linearalgebra - Linear Algebra
9
11 package require Tcl ?8.4?
12 package require math::linearalgebra ?1.0.1?
14 ::math::linearalgebra::mkVector ndim value
16 ::math::linearalgebra::mkUnitVector ndim ndir
18 ::math::linearalgebra::mkMatrix nrows ncols value
20 ::math::linearalgebra::getrow matrix row
22 ::math::linearalgebra::setrow matrix row newvalues
24 ::math::linearalgebra::getcol matrix col
26 ::math::linearalgebra::setcol matrix col newvalues
28 ::math::linearalgebra::getelem matrix row col
30 ::math::linearalgebra::setelem matrix row ?col? newvalue
32 ::math::linearalgebra::show obj ?format? ?rowsep? ?colsep?
34 ::math::linearalgebra::dim obj
36 ::math::linearalgebra::shape obj
38 ::math::linearalgebra::conforming type obj1 obj2
40 ::math::linearalgebra::symmetric matrix ?eps?
42 ::math::linearalgebra::norm vector type
44 ::math::linearalgebra::norm_one vector
46 ::math::linearalgebra::norm_two vector
48 ::math::linearalgebra::norm_max vector
50 ::math::linearalgebra::normMatrix matrix type
52 ::math::linearalgebra::dotproduct vect1 vect2
54 ::math::linearalgebra::unitLengthVector vector
56 ::math::linearalgebra::normalizeStat mv
58 ::math::linearalgebra::axpy scale mv1 mv2
60 ::math::linearalgebra::add mv1 mv2
62 ::math::linearalgebra::sub mv1 mv2
64 ::math::linearalgebra::scale scale mv
66 ::math::linearalgebra::rotate c s vect1 vect2
68 ::math::linearalgebra::transpose matrix
70 ::math::linearalgebra::matmul mv1 mv2
72 ::math::linearalgebra::angle vect1 vect2
74 ::math::linearalgebra::crossproduct vect1 vect2
76 ::math::linearalgebra::matmul mv1 mv2
78 ::math::linearalgebra::mkIdentity size
80 ::math::linearalgebra::mkDiagonal diag
82 ::math::linearalgebra::mkHilbert size
84 ::math::linearalgebra::mkDingdong size
86 ::math::linearalgebra::mkOnes size
88 ::math::linearalgebra::mkMoler size
90 ::math::linearalgebra::mkFrank size
92 ::math::linearalgebra::mkBorder size
94 ::math::linearalgebra::mkWilkinsonW+ size
96 ::math::linearalgebra::mkWilkinsonW- size
98 ::math::linearalgebra::solveGauss matrix bvect
100 ::math::linearalgebra::solveTriangular matrix bvect
102 ::math::linearalgebra::solveGaussBand matrix bvect
104 ::math::linearalgebra::solveTriangularBand matrix bvect
106 ::math::linearalgebra::determineSVD A eps
108 ::math::linearalgebra::eigenvectorsSVD A eps
110 ::math::linearalgebra::leastSquaresSVD A y qmin eps
112 ::math::linearalgebra::choleski matrix
114 ::math::linearalgebra::orthonormalizeColumns matrix
116 ::math::linearalgebra::orthonormalizeRows matrix
118 ::math::linearalgebra::to_LA mv
120 ::math::linearalgebra::from_LA mv
122_________________________________________________________________
124
126 This package offers both low-level procedures and high-level algorithms
127 to deal with linear algebra problems:
128 · robust solution of linear equations or least squares problems
130 · determining eigenvectors and eigenvalues of symmetric matrices
132 · various decompositions of general matrices or matrices of a spe‐
134 cific form
135 · (limited) support for matrices in band storage, a common type of
137 sparse matrices It arose as a re-implementation of Hume's LA
138 package and the desire to offer low-level procedures as found in
139 the well-known BLAS library. Matrices are implemented as lists
140 of lists rather linear lists with reserved elements, as in the
141 original LA package, as it was found that such an implementation
142 is actually faster.
143 It is advisable, however, to use the procedures that are offered, such
145 as setrow and getrow, rather than rely on this representation explic‐
146 itly: that way it is to switch to a possibly even faster compiled
147 implementation that supports the same API.
148 Note: When using this package in combination with Tk, there may be a
150 naming conflict, as both this package and Tk define a command scale.
151 See the NAMING CONFLICT section below.
152
154 The package defines the following public procedures (several exist as
155 specialised procedures, see below):
156 Constructing matrices and vectors
158 ::math::linearalgebra::mkVector ndim value
160 Create a vector with ndim elements, each with the value value.
161 ndim integer Dimension of the vector (number of components)
163 value double Uniform value to be used (default: 0.0)
165 ::math::linearalgebra::mkUnitVector ndim ndir
168 Create a unit vector in ndim-dimensional space, along the ndir-
169 th direction.
170 ndim integer Dimension of the vector (number of components)
172 ndir integer Direction (0, ..., ndim-1)
174 ::math::linearalgebra::mkMatrix nrows ncols value
177 Create a matrix with nrows rows and ncols columns. All elements
178 have the value value.
179 nrows integer Number of rows
181 ncols integer Number of columns
183 value double Uniform value to be used (default: 0.0)
185 ::math::linearalgebra::getrow matrix row
188 Returns a single row of a matrix as a list
189 matrix list Matrix in question
191 row int Index of the row to return
193 ::math::linearalgebra::setrow matrix row newvalues
196 Set a single row of a matrix to new values (this list must have
197 the same number of elements as the number of columns in the
198 matrix)
199 matrix list name of the matrix in question
201 row int Index of the row to update
203 newvalues list List of new values for the row
205 ::math::linearalgebra::getcol matrix col
208 Returns a single column of a matrix as a list
209 matrix list Matrix in question
211 col int Index of the column to return
213 ::math::linearalgebra::setcol matrix col newvalues
216 Set a single column of a matrix to new values (this list must
217 have the same number of elements as the number of rows in the
218 matrix)
219 matrix list name of the matrix in question
221 col int Index of the column to update
223 newvalues list List of new values for the column
225 ::math::linearalgebra::getelem matrix row col
228 Returns a single element of a matrix/vector
229 matrix list Matrix or vector in question
231 row int Row of the element
233 col int Column of the element (not present for vectors)
235 ::math::linearalgebra::setelem matrix row ?col? newvalue
238 Set a single element of a matrix (or vector) to a new value
239 matrix list name of the matrix in question
241 row int Row of the element
243 col int Column of the element (not present for vectors)
245 Querying matrices and vectors
247 ::math::linearalgebra::show obj ?format? ?rowsep? ?colsep?
249 Return a string representing the vector or matrix, for easy
250 printing. (There is currently no way to print fixed sets of
251 columns)
252 obj list Matrix or vector in question
254 format string Format for printing the numbers (default: %6.4f)
256 rowsep string String to use for separating rows (default: new‐
258 line)
259 colsep string String to use for separating columns (default:
261 space)
262 ::math::linearalgebra::dim obj
265 Returns the number of dimensions for the object (either 0 for a
266 scalar, 1 for a vector and 2 for a matrix)
267 obj any Scalar, vector, or matrix
269 ::math::linearalgebra::shape obj
272 Returns the number of elements in each dimension for the object
273 (either an empty list for a scalar, a single number for a vector
274 and a list of the number of rows and columns for a matrix)
275 obj any Scalar, vector, or matrix
277 ::math::linearalgebra::conforming type obj1 obj2
280 Checks if two objects (vector or matrix) have conforming shapes,
281 that is if they can be applied in an operation like addition or
282 matrix multiplication.
283 type string Type of check:
285 · "shape" - the two objects have the same shape (for
287 all element-wise operations)
288 · "rows" - the two objects have the same number of
290 rows (for use as A and b in a system of linear
291 equations Ax = b
292 · "matmul" - the first object has the same number of
294 columns as the number of rows of the second
295 object. Useful for matrix-matrix or matrix-vector
296 multiplication.
297 obj1 list First vector or matrix (left operand)
299 obj2 list Second vector or matrix (right operand)
301 ::math::linearalgebra::symmetric matrix ?eps?
304 Checks if the given (square) matrix is symmetric. The argument
305 eps is the tolerance.
306 matrix list Matrix to be inspected
308 eps float Tolerance for determining approximate equality
310 (defaults to 1.0e-8)
311 Basic operations
313 ::math::linearalgebra::norm vector type
315 Returns the norm of the given vector. The type argument can be:
316 1, 2, inf or max, respectively the sum of absolute values, the
317 ordinary Euclidean norm or the max norm.
318 vector list Vector, list of coefficients
320 type string Type of norm (default: 2, the Euclidean norm)
322 ::math::linearalgebra::norm_one vector
324 Returns the L1 norm of the given vector, the sum of absolute
325 values
326 vector list Vector, list of coefficients
328 ::math::linearalgebra::norm_two vector
330 Returns the L2 norm of the given vector, the ordinary Euclidean
331 norm
332 vector list Vector, list of coefficients
334 ::math::linearalgebra::norm_max vector
336 Returns the Linf norm of the given vector, the maximum absolute
337 coefficient
338 vector list Vector, list of coefficients
340 ::math::linearalgebra::normMatrix matrix type
343 Returns the norm of the given matrix. The type argument can be:
344 1, 2, inf or max, respectively the sum of absolute values, the
345 ordinary Euclidean norm or the max norm.
346 matrix list Matrix, list of row vectors
348 type string Type of norm (default: 2, the Euclidean norm)
350 ::math::linearalgebra::dotproduct vect1 vect2
353 Determine the inproduct or dot product of two vectors. These
354 must have the same shape (number of dimensions)
355 vect1 list First vector, list of coefficients
357 vect2 list Second vector, list of coefficients
359 ::math::linearalgebra::unitLengthVector vector
362 Return a vector in the same direction with length 1.
363 vector list Vector to be normalized
365 ::math::linearalgebra::normalizeStat mv
368 Normalize the matrix or vector in a statistical sense: the mean
369 of the elements of the columns of the result is zero and the
370 standard deviation is 1.
371 mv list Vector or matrix to be normalized in the above sense
373 ::math::linearalgebra::axpy scale mv1 mv2
376 Return a vector or matrix that results from a "daxpy" operation,
377 that is: compute a*x+y (a a scalar and x and y both vectors or
378 matrices of the same shape) and return the result.
379 Specialised variants are: axpy_vect and axpy_mat (slightly
381 faster, but no check on the arguments)
382 scale double The scale factor for the first vector/matrix (a)
384 mv1 list First vector or matrix (x)
386 mv2 list Second vector or matrix (y)
388 ::math::linearalgebra::add mv1 mv2
391 Return a vector or matrix that is the sum of the two arguments
392 (x+y)
393 Specialised variants are: add_vect and add_mat (slightly faster,
395 but no check on the arguments)
396 mv1 list First vector or matrix (x)
398 mv2 list Second vector or matrix (y)
400 ::math::linearalgebra::sub mv1 mv2
403 Return a vector or matrix that is the difference of the two
404 arguments (x-y)
405 Specialised variants are: sub_vect and sub_mat (slightly faster,
407 but no check on the arguments)
408 mv1 list First vector or matrix (x)
410 mv2 list Second vector or matrix (y)
412 ::math::linearalgebra::scale scale mv
415 Scale a vector or matrix and return the result, that is: compute
416 a*x.
417 Specialised variants are: scale_vect and scale_mat (slightly
419 faster, but no check on the arguments)
420 scale double The scale factor for the vector/matrix (a)
422 mv list Vector or matrix (x)
424 ::math::linearalgebra::rotate c s vect1 vect2
427 Apply a planar rotation to two vectors and return the result as
428 a list of two vectors: c*x-s*y and s*x+c*y. In algorithms you
429 can often easily determine the cosine and sine of the angle, so
430 it is more efficient to pass that information directly.
431 c double The cosine of the angle
433 s double The sine of the angle
435 vect1 list First vector (x)
437 vect2 list Seocnd vector (x)
439 ::math::linearalgebra::transpose matrix
442 Transpose a matrix
443 matrix list Matrix to be transposed
445 ::math::linearalgebra::matmul mv1 mv2
448 Multiply a vector/matrix with another vector/matrix. The result
449 is a matrix, if both x and y are matrices or both are vectors,
450 in which case the "outer product" is computed. If one is a vec‐
451 tor and the other is a matrix, then the result is a vector.
452 mv1 list First vector/matrix (x)
454 mv2 list Second vector/matrix (y)
456 ::math::linearalgebra::angle vect1 vect2
459 Compute the angle between two vectors (in radians)
460 vect1 list First vector
462 vect2 list Second vector
464 ::math::linearalgebra::crossproduct vect1 vect2
467 Compute the cross product of two (three-dimensional) vectors
468 vect1 list First vector
470 vect2 list Second vector
472 ::math::linearalgebra::matmul mv1 mv2
475 Multiply a vector/matrix with another vector/matrix. The result
476 is a matrix, if both x and y are matrices or both are vectors,
477 in which case the "outer product" is computed. If one is a vec‐
478 tor and the other is a matrix, then the result is a vector.
479 mv1 list First vector/matrix (x)
481 mv2 list Second vector/matrix (y)
483 Common matrices and test matrices
485 ::math::linearalgebra::mkIdentity size
487 Create an identity matrix of dimension size.
488 size integer Dimension of the matrix
490 ::math::linearalgebra::mkDiagonal diag
493 Create a diagonal matrix whose diagonal elements are the ele‐
494 ments of the vector diag.
495 diag list Vector whose elements are used for the diagonal
497 ::math::linearalgebra::mkHilbert size
500 Create a Hilbert matrix of dimension size. Hilbert matrices are
501 very ill-conditioned with respect to eigenvalue/eigenvector
502 problems. Therefore they are good candidates for testing the
503 accuracy of algorithms and implementations.
504 size integer Dimension of the matrix
506 ::math::linearalgebra::mkDingdong size
509 Create a "dingdong" matrix of dimension size. Dingdong matrices
510 are imprecisely represented, but have the property of being very
511 stable in such algorithms as Gauss elimination.
512 size integer Dimension of the matrix
514 ::math::linearalgebra::mkOnes size
517 Create a square matrix of dimension size whose entries are all
518 1.
519 size integer Dimension of the matrix
521 ::math::linearalgebra::mkMoler size
524 Create a Moler matrix of size size. (Moler matrices have a very
525 simple Choleski decomposition. It has one small eigenvalue and
526 it can easily upset elimination methods for systems of linear
527 equations.)
528 size integer Dimension of the matrix
530 ::math::linearalgebra::mkFrank size
533 Create a Frank matrix of size size. (Frank matrices are fairly
534 well-behaved matrices)
535 size integer Dimension of the matrix
537 ::math::linearalgebra::mkBorder size
540 Create a bordered matrix of size size. (Bordered matrices have a
541 very low rank and can upset certain specialised algorithms.)
542 size integer Dimension of the matrix
544 ::math::linearalgebra::mkWilkinsonW+ size
547 Create a Wilkinson W+ of size size. This kind of matrix has
548 pairs of eigenvalues that are very close together. Usually the
549 order (size) is odd.
550 size integer Dimension of the matrix
552 ::math::linearalgebra::mkWilkinsonW- size
555 Create a Wilkinson W- of size size. This kind of matrix has
556 pairs of eigenvalues with opposite signs, when the order (size)
557 is odd.
558 size integer Dimension of the matrix
560 Common algorithms
562 ::math::linearalgebra::solveGauss matrix bvect
564 Solve a system of linear equations (Ax=b) using Gauss elimina‐
565 tion. Returns the solution (x) as a vector or matrix of the
566 same shape as bvect.
567 matrix list Square matrix (matrix A)
569 bvect list Vector or matrix whose columns are the individual
571 b-vectors
572 ::math::linearalgebra::solveTriangular matrix bvect
575 Solve a system of linear equations (Ax=b) by backward substitu‐
576 tion. The matrix is supposed to be upper-triangular.
577 matrix list Upper-triangular matrix (matrix A)
579 bvect list Vector or matrix whose columns are the individual
581 b-vectors
582 ::math::linearalgebra::solveGaussBand matrix bvect
584 Solve a system of linear equations (Ax=b) using Gauss elimina‐
585 tion, where the matrix is stored as a band matrix (cf. STORAGE).
586 Returns the solution (x) as a vector or matrix of the same shape
587 as bvect.
588 matrix list Square matrix (matrix A; in band form)
590 bvect list Vector or matrix whose columns are the individual
592 b-vectors
593 ::math::linearalgebra::solveTriangularBand matrix bvect
596 Solve a system of linear equations (Ax=b) by backward substitu‐
597 tion. The matrix is supposed to be upper-triangular and stored
598 in band form.
599 matrix list Upper-triangular matrix (matrix A)
601 bvect list Vector or matrix whose columns are the individual
603 b-vectors
604 ::math::linearalgebra::determineSVD A eps
607 Determines the Singular Value Decomposition of a matrix: A = U S
608 Vtrans. Returns a list with the matrix U, the vector of singu‐
609 lar values S and the matrix V.
610 A list Matrix to be decomposed
612 eps float Tolerance (defaults to 2.3e-16)
614 ::math::linearalgebra::eigenvectorsSVD A eps
617 Determines the eigenvectors and eigenvalues of a real symmetric
618 matrix, using SVD. Returns a list with the matrix of normalized
619 eigenvectors and their eigenvalues.
620 A list Matrix whose eigenvalues must be determined
622 eps float Tolerance (defaults to 2.3e-16)
624 ::math::linearalgebra::leastSquaresSVD A y qmin eps
627 Determines the solution to a least-sqaures problem Ax ~ y via
628 singular value decomposition. The result is the vector x.
629 Note that if you add a column of 1s to the matrix, then this
631 column will represent a constant like in: y = a*x1 + b*x2 + c.
632 To force the intercept to be zero, simply leave it out.
633 A list Matrix of independent variables
635 y list List of observed values
637 qmin float Minimum singular value to be considered (defaults to
639 0.0)
640 eps float Tolerance (defaults to 2.3e-16)
642 ::math::linearalgebra::choleski matrix
645 Determine the Choleski decomposition of a symmetric positive
646 semidefinite matrix (this condition is not checked!). The result
647 is the lower-triangular matrix L such that L Lt = matrix.
648 matrix list Matrix to be decomposed
650 ::math::linearalgebra::orthonormalizeColumns matrix
653 Use the modified Gram-Schmidt method to orthogonalize and nor‐
654 malize the columns of the given matrix and return the result.
655 matrix list Matrix whose columns must be orthonormalized
657 ::math::linearalgebra::orthonormalizeRows matrix
660 Use the modified Gram-Schmidt method to orthogonalize and nor‐
661 malize the rows of the given matrix and return the result.
662 matrix list Matrix whose rows must be orthonormalized
664 Compability with the LA package
666 ::math::linearalgebra::to_LA mv
668 Transforms a vector or matrix into the format used by the origi‐
669 nal LA package.
670 mv list Matrix or vector
672 ::math::linearalgebra::from_LA mv
674 Transforms a vector or matrix from the format used by the origi‐
675 nal LA package into the format used by the present implementa‐
676 tion.
677 mv list Matrix or vector as used by the LA package
679
681 While most procedures assume that the matrices are given in full form,
682 the procedures solveGaussBand and solveTriangularBand assume that the
683 matrices are stored as band matrices. This common type of "sparse"
684 matrices is related to ordinary matrices as follows:
685 · "A" is a full-size matrix with N rows and M columns.
687 · "B" is a band matrix, with m upper and lower diagonals and n
689 rows.
690 · "B" can be stored in an ordinary matrix of (2m+1) columns (one
692 for each off-diagonal and the main diagonal) and n rows.
693 · Element i,j (i = -m,...,m; j =1,...,n) of "B" corresponds to
695 element k,j of "A" where k = M+i-1 and M is at least (!) n, the
696 number of rows in "B".
697 · To set element (i,j) of matrix "B" use:
699 setelem B $j [expr {$N+$i-1}] $value
701 (There is no convenience procedure for this yet)
703
705 There is a difference between the original LA package by Hume and the
706 current implementation. Whereas the LA package uses a linear list, the
707 current package uses lists of lists to represent matrices. It turns out
708 that with this representation, the algorithms are faster and easier to
709 implement.
710 The LA package was used as a model and in fact the implementation of,
712 for instance, the SVD algorithm was taken from that package. The set of
713 procedures was expanded using ideas from the well-known BLAS library
714 and some algorithms were updated from the second edition of J.C. Nash's
715 book, Compact Numerical Methods for Computers, (Adam Hilger, 1990) that
716 inspired the LA package.
717 Two procedures are provided to make the transition between the two
719 implementations easier: to_LA and from_LA. They are described above.
720
722 Odds and ends: the following algorithms have not been implemented yet:
723 · determineQR
725 · certainlyPositive, diagonallyDominant
727
729 If you load this package in a Tk-enabled shell like wish, then the com‐
730 mand
731 namespace import ::math::linearalgebra
732 results in an error message about "scale". This is due to the fact that
733 Tk defines all its commands in the global namespace. The solution is to
734 import the linear algebra commands in a namespace that is not the
735 global one:
736 package require math::linearalgebra
738 namespace eval compute {
739 namespace import ::math::linearalgebra::*
740 ... use the linear algebra version of scale ...
741 }
742 To use Tk's scale command in that same namespace you can rename it:
744 namespace eval compute {
746 rename ::scale scaleTk
747 scaleTk .scale ...
748 }
749
752 least squares, linear algebra, linear equations, math, matrices, vec‐
753 tors
754
756 Copyright (c) 2004-2006 Arjen Markus <arjenmarkus@users.sourceforge.net>
757 Copyright (c) 2004 Ed Hume <http://www.hume.com/contact.us.htm>
758math 1.0.1 math::linearalgebra(n)