== name == Sheafhom == category == Mathematics Packages == author == Mark McConnell mailto(mmcconnell17704@yahoo.com) == author-image == mmcconnell17704 == short-description == Sheafhom is a package for large-scale mathematical computations. Its front end is a language for problems in algebraic topology and number theory. These problems come down to large sparse systems of linear equations over integral domains, especially the integers. Sheafhom's back end solves the sparse systems. == long-description == Sheafhom is a package for large-scale mathematical computations. Its front end is a language for problems in algebraic topology and number theory. These problems come down to large sparse systems of linear equations over integral domains, especially the integers. Sheafhom's back end solves the sparse systems. Sheafhom 1.italic(x) was developed in 1993-99 in Common Lisp. Sheafhom 2.0 was written in 2001-2004 in Java. Sheafhom 2.1 is again in Common Lisp, which is ideal for several reasons. Arbitrary-precision integers are built into the language and can be very fast. It is easier to write the code in layers by using macros, as well as CLOS. bold(The back end.) Fill-in is a concern in sparse linear algebra using either exact or inexact (floating-point) numbers. Imagine row-reducing the following matrix, where the letters are arbitrary non-zero numbers (all italic(x)'s distinct). a x 0 0 x x x 0 x x 0 0 x 0 x x 0 x b 0 0 x 0 x We want to choose one row and add multiples of it to the other rows until the first column has only one non-zero entry. If we choose the top row (italic(a) is the italic(pivot)), the result will look like this, where * is a new non-zero entry. a x 0 0 x x 0 * x x * * 0 * x x * x 0 * 0 x * x If italic(b) is the pivot, the result will be 0 x 0 * x x 0 0 x x 0 * 0 0 x x 0 x b 0 0 x 0 x The second result has less fill-in, two *'s rather than seven. Fill-in compounds on itself as we reduce row after row. A sparse solver must reduce fill-in by choosing the pivots cleverly. Sheafhom currently uses the well-established Markowitz pivoting algorithm (see the Book section). It is also a platform for experiments with pivoting strategies that may be better for topology. Floating-point work must balance fill-in with numerical stability. In the example, if italic(b) is much smaller than italic(a), using italic(b) as pivot may introduce too much round-off error. We may be forced to use italic(a), with its larger fill-in. Over the integers, numerical stability is replaced with the opposite problem: integer explosion. In integer multiplication, the length (number of digits) of the product is roughly the sum of the lengths of the factors. By the end of a row-reduction, matrix entries often have hundreds or thousands of digits. Sheafhom has several strategies for reducing integer explosion. Sheafhom 2.italic(x) offers two kinds of graphical tools to study fill-in and integer explosion. * Line graphs show the growth of the sparsity percentage and number of row/column operations. * Sheafhom can generate movies showing the sparsity pattern changing in real time. This slows down the computation, but is a strong asset for testing pivoting algorithms. bold(The front end.) The front end of Sheafhom 1.italic(x) was a language for sheaf theory. A topological space of dimension italic(n) can be glued together from cells of dimensions 0 through italic(n), represented in code as a data structure of face relations. A sheaf is a choice of sparse matrices, one for each face relation, with various commuting properties. A linear map of sheaves is a choice of sparse matrix for each cell commuting with the matrices from the face relations. A complex of sheaves is a sequence of linear maps of sheaves, and so on. CLOS is a good tool for handling these layers of abstractions piled on abstractions. The full front end of Sheafhom 1.italic(x) has not been ported over to 2.1. Version 1.italic(x) did not handle integers modulo italic(n), and the front end will be rewritten in 2.2 to fix this gap. Nevertheless, 2.1 supports good code for topology, as the Tutorial section shows. bold(Language features.) Sheafhom 2.1's lowest layer is a language of macros that expresses matrices in the bare bones of Lisp, integers for entries and lists for sparse vectors. Macros optimize the low-level code, with prolix declarations and with macroexpand and disassemble for checking the work. Above the lower levels, Sheafhom treats speed as a second--though high--priority. Avoiding fill-in and integer explosion is the first priority. In the Tutorial, we will reduce a matrix of 26 rows and over 173,000 columns. The macro with-splicer embodies a mini-language for surgery on lists. It iterates down a list and offers insert, delete, and larger splicing commands, all without disturbing the iteration. The core sparse vector routines rely on with-splicer. Sheafhom 2.1 is strictly in ANSI CL. The only exception is the file gui.lisp with the graphical tools. This uses Allegro's jlinker to call out to Java. It could be adapted to other platforms. The back end of Sheafhom 2.1 uses integers. However, the low-level macro language has been carefully written so it can be extended to finite fields, rings of integers in algebraic number fields, and other types of "numbers". == examples == Please see the Tutorial section. == instructions == * If you're not using Allegro CL with jlinker, please open the file gui.lisp and edit it as indicated there. * Open an ANSI Common Lisp. Change to the directory where the lisp files are located--in ACL, type the equivalent of :cd \mark\math\shh2 Then type (load "make") == tutorial == bold(Definition of Smith normal form.) The Smith normal form (SNF) of an integer matrix is a way of breaking it down into basic pieces while respecting the integer structure. Let's begin with an example. International Allegro CL Professional Edition 7.0 [Windows] (Oct 20, 2004 22:05) Copyright (C) 1985-2004, Franz Inc., Oakland, CA, USA. All Rights Reserved. CL-USER(2): :cd \mark\math\shh2 c:\mark\math\shh2\ CL-USER(3): (load "make") [list of files being compiled and loaded] Enter a 3 by 3 matrix, listing the entries row by row. CL-USER(4): (defvar mm (input-csparse 3 3 '(5 -1 0 0 2 10 7 5 32))) MM CL-USER(5): mm 5 -1 0 0 2 10 7 5 32 This matrix is an instance of the class bold(csparse). The name stands for "column sparse". A csparse stores a matrix by columns, and each column is a sparse vector. Now compute the SNF of the matrix and take a look at its internals. CL-USER(6): (defvar snf-mm (make-snf mm)) SNF-MM CL-USER(7): (describe snf-mm) [SNF 3 by 3, [one unit] 2 [one 0]] is an instance of #: The following slots have :INSTANCE allocation: MAT 5 -1 0 0 2 10 7 5 32 D -1 0 0 0 2 0 0 0 0 USE-P T USE-Q T P-LIST ((T 1 0 -2) (T 2 0 -5) (T 2 1 3) (SHEAFHOM2::P 1 2) (T 2 1 5)) Q-LIST ((T 1 2 1) (T 0 1 -5) (SHEAFHOM2::P 0 1)) This say mm is a product of matrices * mm = italic(P0) italic(P1) italic(P2) italic(P3) italic(P4) italic(D) italic(Q0) italic(Q1) italic(Q2) The italic(P)'s and italic(Q)'s are italic(elementary matrices). For instance, italic(P0), the 0th element of p-list, is a italic(transvection) matrix, built by taking the identity matrix and inserting -2 in the 1,0 position: 1 0 0 -2 1 0 0 0 1 For a given matrix italic(A), the product italic(P0 A) alters the middle row of italic(A) by subtracting off two times the top row. italic(P3) is the 2 by 2 italic(permutation) matrix 1 0 0 0 0 1 0 1 0 The product italic(P3 A) interchanges the middle and bottom rows of italic(A). Multiplying by elementary matrices on the right performs the analogous column operations. The italic(P)'s and italic(Q)'s are isomorphisms representing changes of basis. The essential structure of mm is in the central matrix italic(D), the diagonal matrix -1 0 0 0 2 0 0 0 0 Because there are only two non-zero diagonal entries, mm has rank 2. Viewing mm for the moment over the real numbers, rank 2 means the image of mm has only two dimensions instead of three. The kernel of mm is one-dimensional, generated by the following vector. CL-USER(8): (shh2::ker snf-mm) -1 -5 1 (We need the prefix "shh2::" because the current implementations of kernel, cokernel, and so on are provisional. We expect them to be made public in the next release.) Over the real numbers, the quotient of three-dimensional space by the image of mm is one-dimensional. Over the integers, there is more to the story, because the middle entry of italic(D) is 2 rather than +/- 1. The quotient is bold(Z) + (bold(Z) / 2 bold(Z)) including the italic(torsion) subgroup bold(Z) / 2 bold(Z). Printing an SNF summarizes what's on the diagonal of italic(D). CL-USER(9): snf-mm [SNF 3 by 3, [one unit] 2 [one 0]] A "unit" is an invertible element in a ring. Over the integers, the units are +/- 1. After the units, we list the torsion coefficients and the number of trailing zeroes. bold(Operating on csparses.) We'll need the next function later in the tutorial. CL-USER(10): (defun csparse-incf (matrix i j) "Adds 1 to the entry in the i,j position." (let* ((entry (csparse-get matrix i j)) (val (if (null entry) 0 (sp-v entry))) (new-val (1+ val)) (new-entry (make-sp i new-val))) (csparse-set matrix new-entry j))) CSPARSE-INCF Let's walk through this function line by line, applying its steps to the csparse mm. CL-USER(11): mm 5 -1 0 0 2 10 7 5 32 We'll add a 1 to the 1,2 position, changing the 10 to 11. A csparse is stored by columns, and each column is a sparse vector. The "10" in column 2 must maintain its index in the vector (1) in addition to its value (10). Sheafhom wraps this data into an instance of the type bold(sparse-elt). The function csparse-get extracts a sparse-elt. CL-USER(12): (csparse-get mm 1 2) (1 . 10) Over the integers, a sparse-elt is just a cons holding the index and value. For other rings, the representation could be more complicated. We extract the index of a sparse-elt with the macro sp-i, and the value with sp-v. CL-USER(13): (sp-v *) 10 To replace the 10 with 11, we need to create a new sparse-elt with index 1 and value 11. The macro make-sp creates a sparse-elt. CL-USER(14): (make-sp 1 11) (1 . 11) The function csparse-set puts the new sparse-elt into column 2 of the csparse. CL-USER(15): (csparse-set mm * 2) 5 -1 0 0 2 11 7 5 32 This change raises the rank of mm to 3. There's also a lot more torsion: the quotient by the image is cyclic of order 32. CL-USER(16): (make-snf mm) [SNF 3 by 3, [two units] 32 [zero 0s]] Sheafhom has many more operators. Check the documentation or the source files for more information. Operators on sparse-elts have "sp" in their names. Operators on sparse vectors (type bold(sparse-v)) have "v", csparses have "csparse", and general matrix operations have "m". An operator named "make-" is a basic constructor, and "input-" is a user-friendly constructor for use at the command prompt. "Alter" and "swap" provide elementary row/column operations. bold(Modulo the dictionary.) The funny article href(Homophonic quotients of free groups, http://www.expmath.org/restricted/2/2.3/mestre.ps) looks at the free group on the letters A through Z modulo equivalence of words that have the same pronunciation in French. For instance, italic(soie) = italic(soi) implies italic(e) = 1. The authors prove the quotient is trivial. Let's solve a related problem: find the quotient of the free abelian group on A through Z by the relation italic(WWWW) = 0, where italic(WWWW) runs through a complete English dictionary. If the quotient is too boring, we can also try subsets of the dictionary. We'll use href(this dictionary,http://www.itasoftware.com/careers/WORD.LST), which comes from a site that's of interest to Lispers for several reasons. :-) In the dictionary, no hyphens or other punctuation marks appear, just the letters a through z. Words of length one are excluded because they make combinatorial problems like this one too easy. We load the words into a list for flexibility. CL-USER(17): (defconstant +words+ (let ((result '())) (with-open-file (stream "WORD.LST" :direction :input) (loop (let ((symbol (read stream nil 'end-of-file))) (if (eq symbol 'end-of-file) (return (nreverse result)) (push (symbol-name symbol) result))))))) +WORDS+ CL-USER(18): +words+ ("AA" "AAH" "AAHED" "AAHING" "AAHS" "AAL" "AALII" "AALIIS" "AALS" "AARDVARK" ...) CL-USER(19): (length +words+) 173528 CL-USER(20): (nth 100000 +words+) "NONJURY" What's the range of word lengths? CL-USER(21): (remove-duplicates (sort (mapcar #'length +words+) #'<)) (2 3 4 5 6 7 8 9 10 11 ...) While we're at it, what is the longest word? CL-USER(22): (last *) (28) CL-USER(23): (remove-if-not #'(lambda (w) (= (length w) 28)) +words+) ("ETHYLENEDIAMINETETRAACETATES") ANTIDISESTABLISHMENTARIANISM is also 28 letters long, though it isn't in this dictionary for some reason. Back to the problem on free abelian groups. We'll create a 26 by 173528 csparse with one column for each word. The column for DAD is 1 0 0 2 0 ... 0 standing for one A and two D's. Here is a utility for converting A through Z to the row indices 0 through 25. CL-USER(24): (defconstant +code-A+ (char-code #\A)) +CODE-A+ CL-USER(25): (defun az-to-int (ch) "Converts ASCII A to 0, B to 1, ..., Z to 25." (let ((result (- (char-code ch) +code-A+))) (assert (<= 0 result 25) (ch) "Not between A and Z: ~:C" ch) result)) AZ-TO-INT The next function takes a list of words and returns the csparse. CL-USER(26): (defun az-list-to-csparse (words) (let ((result (make-csparse-zero 26 (length words))) (j 0)) (dolist (word words result) (map nil #'(lambda (ch) (csparse-incf result (az-to-int ch) j)) word) (incf j)))) AZ-LIST-TO-CSPARSE Let's apply this to the whole dictionary. CL-USER(27): (compile 'az-to-int) ;; to speed things up AZ-TO-INT NIL NIL CL-USER(28): (compile 'az-list-to-csparse) AZ-LIST-TO-CSPARSE NIL NIL CL-USER(29): (defvar words-all (az-list-to-csparse +words+)) WORDS-ALL ;; after a little while CL-USER(30): words-all [CSPARSE 26 by 173528] Finally, we find the SNF. CL-USER(31): (let ((*show-stats* t)) (make-snf words-all nil nil t)) [SNF 26 by 173528, [twenty-six units] [zero 0s]] Since the SNF has nothing but units down the diagonal, we've proved a theorem: the free abelian group on A through Z modulo the dictonary is trivial! We need to explain several items in the last form. When *show-stats* is true, make-snf pops up a line graph showing how the pattern of sparsity changed during the computation. The graph shows that the matrix started with 27.8% of the entries non-zero. The sparsity decreased slightly while the first 20 of the 26 rows were being reduced, then increased till the end. When *show-prog* is t, make-snf prints timing data to standard output. Another flag *show-csw* will be explained below. The last argument t to make-snf means the SNF was computed destructively. This is crucial for saving memory with large matrices. The entries of words-all were modified in place, with sparse-elts thrown to the garbage collector whenever possible. If the last argument had been nil (the default), words-all would have been preserved, and the snf would have been computed using a copy. As another way to save space, we gave the arguments "nil nil" to make-snf. These stop the program from storing the italic(P)'s and italic(Q)'s. Only the main italic(D) portion was recorded. As we've suggested, the trivial group is boring. What is the quotient of A through Z modulo the two-letter words? The quotient will be at least bold(Z) / 2 bold(Z) for parity reasons: no word with an odd number of letters can be equivalent to 0. But are there other relations? CL-USER(32): (remove-if-not #'(lambda (w) (= (length w) 2)) +words+) ("AA" "AB" "AD" "AE" "AG" "AH" "AI" "AL" "AM" "AN" ...) CL-USER(33): (length *) 96 CL-USER(34): (az-list-to-csparse **) [CSPARSE 26 by 96] With *show-csw* true, we see a movie of how the sparsity pattern changes as the SNF is computed. Csw stands for "csparse window". CL-USER(35): (let ((*show-csw* t)) (make-snf * nil nil t)) [SNF 26 by 96, [twenty-one units] -2 [four 0s]] The torsion coefficient 2 (up to sign) is the one we predicted. But the matrix has rank 22. Roughly speaking, there are four linearly independent "words" that are not combinations of two-letter words. The quotient is bold(Z)^4 + (bold(Z) / 2 bold(Z)). What about four-letter words? CL-USER(36): (remove-if-not #'(lambda (w) (= (length w) 4)) +words+) ("AAHS" "AALS" "ABAS" "ABBA" "ABBE" "ABED" "ABET" "ABLE" "ABLY" "ABOS" ...) CL-USER(37): (az-list-to-csparse **) [CSPARSE 26 by 3919] CL-USER(38): (make-snf * nil nil t) [SNF 26 by 3919, [twenty-five units] 4 [zero 0s]] We have torsion of order 4, again for parity reasons. Apart from that, the quotient by four-letter words is a small as possible. What about lengths greater than or equal to some cut-off value? The same methods show that the quotient by words of length 20 or more is trivial. But for length 21 or more, the quotient is bold(Z). CL-USER(39): (remove-if-not #'(lambda (w) (>= (length w) 21)) +words+) ("ACETYLCHOLINESTERASES" "ADRENOCORTICOSTEROIDS" "ADRENOCORTICOTROPHINS" "ANTHROPOMORPHIZATIONS" "ANTIAUTHORITARIANISMS" "ANTIFERROMAGNETICALLY" "BUCKMINSTERFULLERENES" "CARBOXYMETHYLCELLULOSE" "CARBOXYMETHYLCELLULOSES" "CLINICOPATHOLOGICALLY" ...) CL-USER(40): (make-snf (az-list-to-csparse *) t nil t) [SNF 26 by 118, [twenty-five units] [one 0]] Note the arguments "t nil t", meaning we remember the italic(P)'s but not the italic(Q)'s, and work destructively. Because we remembered the italic(P)'s, we can use shh::coker-section, which picks out a generator for the quotient. CL-USER(41): (shh2::coker-section *) The result is a column vector with all entries 0 except for a 1 standing for the letter Q. We've proved that CARBOXYMETHYLQCELLULOSE is not a word. bold(Topology in dimension two.) Sheafhom was written to solve problems in algebraic topology and derived categories. The file two-complex.lisp sets up a simplicial complex X. It's restricted to dimension two to make things easier. The method make-d turns the boundary maps of X into csparses. The method print-homology produces the homology of X with coefficients in bold(Z), including the torsion. Here is a brief example. Consider a 3 by 3 grid of squares with vertices numbered as shown. 0 2 1 0 6 7 8 6 3 4 5 3 0 1 2 0 We list the nine squares. CL-USER(42): *klein-bottle-cells* ((6 7 2 0) (7 8 1 2) (8 6 0 1) (3 4 7 6) (4 5 8 7) (5 3 6 8) (0 1 4 3) (1 2 5 4) (2 0 3 5)) Now attach the left side of the grid to the right side by gluing 0 to 0, 3 to 3, and 6 to 6. The result is a cylindrical tube. Next, glue the bottom circle of the cylinder to the top circle, 0 to 0, 1 to 1, and 2 to 2. The trouble is that the order is reversed, 0-1-2-0 on the bottom and 0-2-1-0 on the top. After squinting, you'll see you can't do the gluing in three dimensions. You have to bend the tube into the fourth dimension and turn one end "upside down" to match up the circles in the proper order. The resulting space is a manifold of dimension two, the Klein bottle. The file two-complex.lisp lets you construct the Klein bottle this way. CL-USER(43): (defvar *klein-bottle* (make-two-complex *)) *KLEIN-BOTTLE* CL-USER(44): (print-homology *klein-bottle*) H_0 is free of rank 1 H_1 has free part of rank 1 and torsion coeffs (2) H_2 is free of rank 0 NIL == home-url == http://www.geocities.com/mmcconnell17704/math.html == doc-url == http://www.geocities.com/mmcconnell17704/shh21-manual.txt == license == Public domain. == book == Direct Methods for Sparse Matrices, I. S. Duff, A. M. Erisman and J. K. Reid, Clarendon Press, 1987. The authors are experts in sparse matrices over the real numbers. The book covers the Markowitz algorithm and other pivoting algorithms. Cohomology of Sheaves, Birger Iversen, Universitext, Springer-Verlag, 1986. The inspiration for the abstract front end of Sheafhom: exact categories, abelian categories, complexes, and the derived category. Methods of Homological Algebra, 2nd ed., S. I. Gelfand and Yu. I. Manin, Springer Monographs in Mathematics, Springer, 2003. A modern approach to derived categories. Basic Algebra I, 2nd ed., Nathan Jacobson, W. H. Freeman and Co., 1985. A serious abstract algebra text. Where I learned Smith normal form. A Course in Computational Algebraic Number Theory, Henri Cohen, Graduate Text in Mathematics 138, Springer, 1993. The bible in its field. Discusses Smith normal form. == references == The following papers are related to Sheafhom or its ancestors. The rational homology of toric varieties is not a combinatorial invariant, Mark McConnell, Proc. AMS 105 (1989), 986-991. Cohomology at infinity and the well-rounded retract for general linear groups, Avner Ash and Mark McConnell, Duke Math. Jour. 90 (1997), 549-576. Cohomology of congruence subgroups of SL(4,Z), Avner Ash, Paul Gunnells and Mark McConnell, J. Number Theory 94 (2002), no. 1, 181-212. == source-fooball == href(Source code and documentation,http://www.geocities.com/mmcconnell17704/math.html) for Sheafhom 1.italic(x), 2.0 and 2.1. You can also download the source code href(here,/source/shh21.zip). == release-date == March 17, 2005 == release-version == 2.1 == status == stable == history == Sheafhom 1.italic(x), Common Lisp, 1993-99. Sheafhom 2.0, Java, 2001-04. Sheafhom 2.1, Common Lisp, 2004-05. == acl-dependencies == None. == other-dependencies == Sheafhom 2.1 is strictly in ANSI CL. The only exception is the file gui.lisp, which uses ACL's jlinker to call out to Java. This file can be edited for other Common Lisps--see the instructions in the file. == platform == No restrictions. == ad ==