How to eliminate torsion

I asked on sci.math.research, and got an answer which works and is easy if you have a Smith Normal Form function available. The packages SAGE and gp mentioned below are freeware, and SAGE is in Python, which I wish I had known or I would have mentioned it before.

Victor S. Miller

Suppose that we're working in Z^n, and that we have an m-dimensional
lattice. Form the m x n integer matrix m whose rows are a basis of
your lattice, and call it H. The Smith normal form of H is an n x n
diagonal integer matrix D, whose top m entries are positive (and the
others are zero) and each such entry is an integer multiple of the one
below it and to the right. Furthermore there are matrices U in
SL(m,Z) and V in SL(n,Z) such that

H = U D V

Your lattice will be gcd-reduced if and only if all the non-zero
elements of D are 1. To form the saturation of your lattice (that's
what I think the standard terminology is) calculate

H' = U D' V

where D' is the matrix obtained from D by replacing all of its
non-zero entries by 1.

Most computer algebra packages have a function for Smith normal form
(such as maple, SAGE, Magma, gp, etc.).

Victor

The only caveat is that you need to clean up the result with something like LLL or Hermite.

--- In tuning-math@yahoogroups.com, "genewardsmith" <genewardsmith@...> wrote:

> H' = U D' V
>
> where D' is the matrix obtained from D by replacing all of its
> non-zero entries by 1.

Using Victor's method, I wrote code to pass from a list of vals to a corresponding list of monzos or vice-versa by the P --> I-P interchange method. It all works like a charm. By using Graham's Cangwu-style matrix definitions and the pseudoinverse defined projection map P, it looks like wedge products could be avoided if one wanted. P could be used internally as a way of uniquely characterizing a temperament, though it's not going to work for human consumption.

Great work Gene! -Carl

At 12:50 PM 7/13/2010, you wrote:
>I asked on sci.math.research, and got an answer which works and is
>easy if you have a Smith Normal Form function available. The packages
>SAGE and gp mentioned below are freeware, and SAGE is in Python, which
>I wish I had known or I would have mentioned it before.
>
>Victor S. Miller
>
>Suppose that we're working in Z^n, and that we have an m-dimensional
>lattice. Form the m x n integer matrix m whose rows are a basis of
>your lattice, and call it H. The Smith normal form of H is an n x n
>diagonal integer matrix D, whose top m entries are positive (and the
>others are zero) and each such entry is an integer multiple of the one
>below it and to the right. Furthermore there are matrices U in
>SL(m,Z) and V in SL(n,Z) such that
>
>H = U D V
>
>Your lattice will be gcd-reduced if and only if all the non-zero
>elements of D are 1. To form the saturation of your lattice (that's
>what I think the standard terminology is) calculate
>
>H' = U D' V
>
>where D' is the matrix obtained from D by replacing all of its
>non-zero entries by 1.
>
>Most computer algebra packages have a function for Smith normal form
>(such as maple, SAGE, Magma, gp, etc.).
>
>Victor
>
>The only caveat is that you need to clean up the result with something
>like LLL or Hermite.
>

On 13 July 2010 20:50, genewardsmith <genewardsmith@sbcglobal.net> wrote:
> I asked on sci.math.research, and got an answer which works and is easy if you have a Smith Normal Form function available. The packages SAGE and gp mentioned below are freeware, and SAGE is in Python, which I wish I had known or I would have mentioned it before.

Sage is a big bag of everything knitted together with Python, so
there's no telling if the function we want is in Python and I can move
it to my website. Anyway, I've got it installed and I'm going through
the tutorial now.

Graham

--- In tuning-math@yahoogroups.com, "genewardsmith" <genewardsmith@...> wrote:

Some email:

After thinking about it I realized that the
answer is even slightly simpler. If d is the dimension of your
lattice, then a basis of the saturation is formed by first d rows of
the matrix V (which is n x n). You are right that some reduction
(probably LLL) might be advisable.

Victor

> Victor S. Miller
>
> Suppose that we're working in Z^n, and that we have an m-dimensional
> lattice. Form the m x n integer matrix m whose rows are a basis of
> your lattice, and call it H. The Smith normal form of H is an n x n
> diagonal integer matrix D, whose top m entries are positive (and the
> others are zero) and each such entry is an integer multiple of the one
> below it and to the right. Furthermore there are matrices U in
> SL(m,Z) and V in SL(n,Z) such that
>
> H = U D V
>
> Your lattice will be gcd-reduced if and only if all the non-zero
> elements of D are 1. To form the saturation of your lattice (that's
> what I think the standard terminology is) calculate
>
> H' = U D' V
>
> where D' is the matrix obtained from D by replacing all of its
> non-zero entries by 1.
>
> Most computer algebra packages have a function for Smith normal form
> (such as maple, SAGE, Magma, gp, etc.).
>
> Victor
>
> The only caveat is that you need to clean up the result with something like LLL or Hermite.
>

On 14 July 2010 23:15, genewardsmith > Some email:
>
> After thinking about it I realized that the
> answer is even slightly simpler.  If d is the dimension of your
> lattice, then a basis of the saturation is formed by first d rows of
> the matrix V (which is n x n).  You are right that some reduction
> (probably LLL) might be advisable.
>
> Victor

Currently, with Sage, it's looking even simpler than that. For a
matrix M, with each row a val, some unison vectors come out from

kernel(transpose(M)).basis()

I haven't tested to see if they ever have spurious torsion, but they
certainly shouldn't. They're given in something that looks like my
version of Hermite normal form, which isn't appropriate for unison
vectors, but that can be fixed.

Graham

On 13 July 2010 22:22, genewardsmith <genewardsmith@sbcglobal.net> wrote:

> Using Victor's method, I wrote code to pass from
> a list of vals to a corresponding list of monzos
> or vice-versa by the P --> I-P interchange method.
> It all works like a charm. By using Graham's
> Cangwu-style matrix definitions and the
> pseudoinverse defined projection map P, it
> looks like wedge products could be avoided if
> one wanted. P could be used internally as a way
> of uniquely characterizing a temperament, though
> it's not going to work for human consumption.

I don't see why we need P either. The Hermite normal form of the
mapping works fine as a unique key.

I'm getting further with torsion now. This paper, which is related to
the Sage documentation, has the clue:

http://www.wstein.org/papers/hnf/pernet-stein-fast_computation_of_hnf_of_random_integer_matrices.pdf

"If A is a basis matrix for M , and H is the Hermite form of the
transpose of A with any 0 rows at the bottom deleted (so H is square),
then H^−1 A is a matrix whose rows are a basis for the saturation of
M."

That actually doesn't work. That's what you get for reading
pre-prints. But it does seem to work if A is also in Hermite normal
form, which is easy enough to arrange as you need that algorithm
anyway. I can't say if this is simpler than Smith normal form or not
because I still don't understand the latter. But I seem to already
have Hermite normal form (even if I don't trust it where there are
zero rows) and inverses for rank 2 and 3 cases, even without a
third-party library.

What also seems to follow is that the determinant of this square
matrix H gives the torsion. I'm crunching through hours of CPU time
to verify it against Sage's index_in_saturation. And for rank>2 it
looks like a simpler way of checking for contorsion than wedge
products.

For the rank 2 case, wedgies are pretty simple. You can generate them
easily and once you have a GCD of 1 (as you will most of the time) you
can stop the calculation.

Graham

--- In tuning-math@yahoogroups.com, Graham Breed <gbreed@...> wrote:

> I don't see why we need P either. The Hermite normal form of the
> mapping works fine as a unique key.

Agreed.

> For the rank 2 case, wedgies are pretty simple. You can generate them
> easily and once you have a GCD of 1 (as you will most of the time) you
> can stop the calculation.

It seems to me that they are also pretty simple for rank 3 or 4. It does get harder because of the rth order (rank = r) polynomial growth in terms of the dimension of the vals, but that's not too terrible.

On 18 July 2010 03:25, genewardsmith <genewardsmith@sbcglobal.net> wrote:
>
> --- In tuning-math@yahoogroups.com, Graham Breed <gbreed@...> wrote:
>
>> For the rank 2 case, wedgies are pretty simple.  You can generate them
>> easily and once you have a GCD of 1 (as you will most of the time) you
>> can stop the calculation.
>
> It seems to me that they are also pretty simple
> for rank 3 or 4. It does get harder because of the
< rth order (rank = r) polynomial growth in terms of
> the dimension of the vals, but that's not too terrible.

If they're are simple formulas, let's see them. I don't have any beyond rank 2.

The size of the vals doesn't matter if you don't have to generate the
whole thing, as I said above.

Graham

--- In tuning-math@yahoogroups.com, Graham Breed <gbreed@...> wrote:

> > It seems to me that they are also pretty simple
> > for rank 3 or 4. It does get harder because of the
> < rth order (rank = r) polynomial growth in terms of
> > the dimension of the vals, but that's not too terrible.
>
> If they're are simple formulas, let's see them. I don't have any beyond rank 2.

Herman's method makes it pretty easy. For rank 3, for instance, you have three vals, call them a, b, and c. Take three nested do loops, i from 1 to n-2, j from i+1 to n-1, and k from j+1 to n. Inside the loops you first increment the indexing number telling you how far up the wedge product you've gotten, and then sum over the six products you get times the parity, and that is the trival coefficent for that index. In other words, take a[i]b[j]c[k] + a[j]b[k]c[i] + a[k]b[i]c[j]
- a[j]b[i]c[k] - a[k]b[j]c[i] - a[i]b[k]c[j]. With added effort you can covert this idea into a more general program good for any list of vals of the same size.

On 18 July 2010 15:31, genewardsmith <genewardsmith@sbcglobal.net> wrote:

> Herman's method makes it pretty easy. For rank 3,
> for instance, you have three vals, call them
> a, b, and c. Take three nested do loops,  i from
> 1 to n-2, j from i+1 to n-1, and k from j+1 to n.
> Inside the loops you first increment the indexing
> number telling you how far up the wedge product
> you've gotten, and then sum over the six products
> you get times the parity, and that is the trival
> coefficent for that index. In other words, take
> a[i]b[j]c[k] + a[j]b[k]c[i] + a[k]b[i]c[j]
> - a[j]b[i]c[k] - a[k]b[j]c[i] - a[i]b[k]c[j]. With added
> effort you can covert this idea into a more general
> program good for any list of vals of the same size.

That's good, so rank 3's fairly easy, but not as simple as rank 2. I
still don't know how to make it more general, except by the old method
of implementing all the wedge products with a multivector object.
Whereas the Hermite normal form code is ugly, but the principle behind
it's simple enough, and it's used for other things. So there's not
much point adding these wedgies as well.

Graham

--- In tuning-math@yahoogroups.com, Graham Breed <gbreed@...> wrote:

> That's good, so rank 3's fairly easy, but not as simple as rank 2. I
> still don't know how to make it more general, except by the old method
> of implementing all the wedge products with a multivector object.

You can write a routine incrementing monotone increasing lists of integers 1<= i <= n. Then you do things in the same way, taking a sum over n! things times the parity of the permutation.

On 19 July 2010 17:27, genewardsmith <genewardsmith@sbcglobal.net> wrote:

> You can write a routine incrementing monotone
> increasing lists of integers 1<= i <= n. Then you
> do things in the same way, taking a sum over
> n! things times the parity of the permutation.

Maybe a simple function could be written, but I don't know how to
write it. And I'm not really interested in putting a lot of energy
into this. I know how to find contorsion without wedge products, and
I know of a good library for integer matrices, so I have made
progress.

The code for dealing with multivectors isn't that complicated, but
looks more so than the Hermite normal form code. The advantage is
that you can use it to do a load of different calculations, including
finding mappings from unison vectors and returning the error and
complexity. It compares favorably with writing my own matrix code.
But I don't know how to calculate the TOP-RMS tuning or the Cangwu
badness without matrices, and I do have a pure Python library that
does these well enough, so wedgies are the more complicated option.

In terms of efficiency, this code will always lose. It's written in
Python and I haven't moved it to any other languages because it uses
Python's syntactic support for dictionaries. I could re-write it in
C, but I'd depend on a hash table library, and that would make the
code more complicated. It also wouldn't buy a huge efficiency gain
because it would spend most of its time dealing with the hash tables.
Python's dictionaries are already implemented using efficient hash
tables coded in c. The wedgies end up as quite large objects for
higher prime limits and that slows things down. If there are
efficient libraries for exterior algebra I don't know of them. So
wedgies generally don't look like the pragmatic option.

Wedgies are still useful for finding the contorsion of a rank 2
temperament (or something that would be a temperament if it didn't
have contorsion). All you need to do is generate the wedgie as a lazy
list. The function's so simple that you don't need to understand
exterior algebra to justify it, and other approaches are likely to
give the same result. It's efficient because it's lazy and it should
translate easily to more efficient languages.

The algorithm you gave for rank 3 is something I could implement. If
you're interested in rank 3 but no higher, it's an option. But, given
that I already have to find the Hermite normal form(HNF), using that
looks easier.

The HNF code has to be there because it's used to give a unique key to
each temperament class. It really is better than wedgies for this.
Wedgies are bigger. Some elements are redundant, but I don't know
which ones in the general case, and reconstructing them adds
complexity. Different things like temperaments but with contorsion
have the same wedgie but different HNFs. And the HNF is close to the
generator/period mapping. So HNFs are clearly earning their keep.

I surely don't have the most efficient HNF algorithm. The difficult
part is the Gaussian elimination to get an echelon matrix. But it can
also be used in a lot of different places. The determinants used for
scalar complexity and badness are easy to calculate for an echelon
matrix. Solving the least squares problem to get the optimal tuning
isn't that difficult once you have an echelon matrix. The simplest
algorithm for finding the kernel uses echelon form. The HNF can be
used not only for detecting torsion but reliably removing it.

It happens that I need different functions for Gaussian elimination of
real and integer matrices. Maybe I could fix that but it doesn't
matter because I'm moving to a third party library for the weighted
calculations. (I still don't know how to get PARI working on my
website, and it's probably not worth the hassle given I'll have to
work it all out again if I move hosts. Note, Carl, this is also why I
keep using cgi-bin.) The same multivector code does work for either
integers or floating point, as it happens.

For the rank 2 case, if you have a library to deal with floating point
matrices, the easiest thing to uniquely identify a temperament class
by is the octave-equivalent generator, which is the same as the
octave-equivalent part of the wedgie, and is easiest to calculate.
There are problems with this but it basically works. So for rank 2
temperaments (which are an important special case) you don't need HNF
and the easiest way to detect torsion is to use wedgies.

The simplest way of dealing with contorsion is to ignore it, which is
what my website currently does. Next is to detect but not remove it,
where wedgies are useful. To detect and remove contorsion, beyond
rank 2, I still need HNF or a third party library.

So that's it, a long answer to a simple question. If you want to show
how simple general wedgie code could be, go ahead an implement it.
It's not something I have a need for.

Graham