About the Code

This is some logic for this code whose use is described here. It also largely applies to this older code.

Consider an n dimensional abstract vector space V. Consider two sets, A and B, of n vectors each. Suppose that each set spans V. Each member of A can be represented as a linear combination of the vectors in B and vice-versa. The coefficients in these two representations each constitute an n by n matrix and these matrices are inverses of each other.

In these programs a matrix is represented as a list of lists of scalars—a list of rows, each row is a list of scalars. (A scalar is an element of the field.) For an n by n matrix the inversion routine may be understood as taking n column vectors and appending n more columns, each describing a basis vector. For n = 3 the picture is:

a₁₁  a₁₂  a₁₃  1  0  0
a₂₁  a₂₂  a₂₃  0  1  0
a₃₁  a₃₂  a₃₃  0  0  1

The 6 columns correspond to the 6 vectors of A ∪ B. The 3 scalars of a column are the coefficients of the linear combination that forms the corresponding vector in terms of the vectors of set B.

Inside the program this is a list of 3 lists of lists each of 6 scalars.

The named let called pad appends these extra columns. It omits trailing zeros which are implicitly supplied as needed.

We apply a sequence of elementary operations to modify the representation of each of these 2n vectors. We justify the operations of the program in light of a sequence of changes to the basis in which our constant 2n vectors are expressed. These changes are thus to be considered as alias transformations—a sequence of new designations for a constant set of 2n vectors in an abstract space V. Each such operation expresses the 2n old vectors in terms of a new gradually changing vector basis. This changing basis is initially the set B and after all our elementary operations the changing basis is the set A. At the end stage the left n² scalars will form the identity matrix and the right n² scalars will be the matrix inverse of the initial left n² scalars. The new picture is thus:

1  0  0  b₁₁  b₁₂  b₁₃
0  1  0  b₂₁  b₂₂  b₂₃
0  0  1  b₃₁  b₃₂  b₃₃

The b’s are the inverse matrix, for the right n columns have continued to depict the old basis vectors thru a sequence of evolving coordinate systems, while the right n rows described and still describe the original input vectors in that same sequence of coordinate systems. The right n rows thus describe the old basis in terms of the input vectors, which is literally the inverse of describing the input vectors in terms of the old basis.

Alas the code is shorter, more efficient and more obscure than these figures suggest. Not only are zeros at the ends of rows omitted but predestined zeros and ones at the beginning of rows as well. I don’t have a clear explanation of why this is more than a space saving hack but the code is both shorter and faster this way. I hesitate to make it longer, slower and consume more space, in order to make it correspond to the simple process above. The definition of matrix inverse has not changed in 170 years and I thus think the code need not evolve. (Note 2013: I now study this to extract the magic of finding x such that Mx=0 to supply that magic in a language that does not abstract as well as Scheme.)

The elementary operations are:

swapping of rows which is performed by the routine tf. tf is called to insure that the first scalar of the first row is not zero for we shall need to divide by it. (This corresponds to swapping two of the basis vectors.) tf is called for sub matrices as we shall see.
multiplying a row by a constant which is used to produce unit diagonal elements for the new identity matrix. (This corresponds to shrinking one of the basis vectors.)
subtracting a scalar multiple of one row from another, which is performed by routine fms. There are two episodes of this. First it is used to force zeros below the emerging diagonal on the left. These operations alternate with calls to tf which may need to deal with zeros not in the original data. The second episode is to force zeros above the diagonal.

Named let inx takes the padded matrix and recursively:

ensures that top left element isn’t 0 by calling tf,
ensures that it is 1 by multiplying the top row by its reciprocal. The top row is clipped of its first element as it is redundantly 1. (The divided clipped top row is the value of b.)
The padded matrix is now shorn of its top row and left column ((map cdr (cdr A))) and a multiple of the top row, b, is subtracted from each clipped row of the new matrix. The missing left row would have been all zeros, for that is how the multiplier was selected.
What remains is an n-1 by 2n-1 matrix to which inx is recursively applied.

Upon completion inx returns an upper triangular matrix with the leading zeros and ones of each row suppressed. This value is named ut. Of course this value is still padded with the transformed names of the old basis vectors.

Now named let ol, to which this was all prelude, takes ut and eliminates the values above the diagonal by yet more elementary transformations.

This may be relevant.

For the logic of a report on singular matrices, see this.

We speculate now on a generalization to the Penrose inverse. To compute the inverse of A we form the doubled matrix [A, I] and derive a sequence of matrices [X, Y] preserving the invariant X = AY where [A, I] is the first in the sequence and [I', A⁺] is the last. I' may have some zeros on the diagonal. Perhaps X = YA is a more convenient invariant. I must understand this before I proceed but that scheme which depends on some notion of orthogonality may fail for finite fields.