14
$\begingroup$

Does anybody know what is the algorithm for generating random unimodular matrices (integer matrices with determinant $\pm 1$) whose elements do not exceed a given bound? Such an algorithm is mentioned here, and the following reference is provided:

Jürgen Hausen, Generating Problems in Linear Algebra, MapleTech, Vol. 1, No. 2, 1994.

However, this paper seems to be no longer accessible online. If the algorithm is based on the Hermite normal form, how do they ensure that the elements of the generated matrix are bounded by a given positive integer? Many thanks in advance for any insights :-)

$\endgroup$
9
  • $\begingroup$ I don't know how but I remember my old linear algebra lecturer saying this was one of the most handy things he knew. $\endgroup$
    – muzzlator
    Mar 21, 2013 at 13:26
  • $\begingroup$ Damn, I should have bountied this question at first sight; it was languishing for quite a while... $\endgroup$ Jun 1, 2013 at 18:20
  • $\begingroup$ This should be an answer, but it involves lattice reduction, which is a rich subject and I wanted to comment in case I do not get around to the work required for a good explanation. Generate a random matrix (so a random determinant). Choose a row to change. The corresponding column in the inverse gives the values with which to reduce using lattice reduction techniques. This will give you the small values for the new row. $\endgroup$
    – adam W
    Jun 4, 2013 at 17:04
  • $\begingroup$ Using the Hermite Normal form also works, it just does not give the most reduced values (though not near as large as Kirill's answer-I think maybe because I use upper not lower triangular, either that or in his application he did not reduce fully to the unique form...) $\endgroup$
    – adam W
    Jun 4, 2013 at 17:05
  • $\begingroup$ adam W: can you please provide some reference to that magic lattice thing? thankx $\endgroup$ Jun 4, 2013 at 18:20

3 Answers 3

6
+100
$\begingroup$

I tried the following approach

  1. generate a unimodular matrix $A$ (possibly with large entries)
  2. fix it to reduce the abs size of the entries

In 1, I generate $n\times n$ matrices $L$ and $U$ with 1 on the diagonal; multiply them and randomly permute rows and columns.

In 2, repeat until "convergence" the following step. Pick the largest element in absolute value, assume it's $a_{ij}$. Pick a random non zero entry $a_{kj}$ on the $j$-th column, $k\ne i$, to reduce the size of $a_{ij}$ with $a'_{i,.} = a_{i,.} + a_{k,.}$.

genunimod <- function(n,k=5){
r <- trunc(sqrt(k))+1
A <- matrix(sample(-r:r,n**2,rep=T),n,n)
L <- matrix(0,n,n)
U <- matrix(0,n,n)
l <- outer(1:n,1:n,">")
u <- outer(1:n,1:n,"<")
L[l] <- A[l]
U[u] <- A[u]
diag(U) <- 1
diag(L) <- 1

p1 <- sample(1:n) #permutation shuffle rows
p2 <- sample(1:n) #permutation shuffle cols
P <- diag(rep(1,n)) # identity

A <- P[p1,] %*% L %*% U %*% P[p2,]
# this A is unimodular but may have large entries

m <- k+1
useless <- 0
steps <- 0
while(m>k & useless<=n**2){
    steps <- steps+1
    mm <- max(abs(A))
    v <- which(abs(A)==mm,arr.ind=T)[1,]        
    # use a random element to reduce element v by row
    k <- sample(which(A[,v[2]]!=0)[-v[1]],1)
    s <- sign(A[v[1],v[2]]*A[k,v[2]])
    A[v[1],] <- A[v[1],]-s*A[k,]
    m <- max(abs(A)) # new max term
    if(m>=mm) useless <- useless+1
    }
list(A=A,LU=P[p1,] %*% L %*% U %*% P[p2,],steps=steps,useless=useless)
}

The code can fail in several ways: the largest element can, in principle, be located in an otherwise null column; a single step can be unsuccessful in reducing the largest entry and I stop computations if more than $n^2$ such useless steps occur.

Efficiency can be vastly improved using less primitive methods (say, using "pivoting" ideas instead of randomly picking $k$, using columns and rows and so on).

It takes 1.2 seconds (on a MacBook Air) to generate 1000 $5\times5$ unimodular matrices with entries $\leq 10$. I never had an "otherwise zero" column (good!) and experienced 69 "failures" (7\%, less good!) due to too many useless steps. Hence, 1.2 seconds were indeed needed to create 931 matrices. The digits have the following bell-shaped distribution. It should be symmetric and more weight is given to small numbers.

enter image description here

$\endgroup$
1
  • $\begingroup$ I suppose a nice alternative would be to grab the unimodular factor from the Hermite decomposition of a random integer matrix, and then apply your "fix" to the unimodular matrix thus generated. Thanks for the elaborate answer! $\endgroup$ Jun 7, 2013 at 1:35
5
$\begingroup$

This isn't a complete answer, but it's a start. The citation is wrong as given, in that the author's name is misspelled. See here (go to vol.1 no.2), which unfortunately has the desired link broken, but here is a fixed version of that link. Alas, the paper is absent, but its abstract is:

This article presents a method for creating exercises in linear algebra by making use of Maple's random number generator. The main tool for doing that is a Maple procedure which generates unimodular integral matrices. The presentation also considers standard exercises for basic topics such as linear independence, bases and the rank of a matrix.

Most importantly, the author's name is spelled Hausen (not Hansen). He is currently a professor at the University of Tübingen. His website is here, and his email address can be found there. Although his research interests appear to have departed from Maple programming, I would suggest contacting him directly with your inquiry.

$\endgroup$
1
2
$\begingroup$

First, the unimodular matrix $U$ that you get from doing a Hermite decomposition $A=U R$ ($U$ is unimodular, $R$ is lower-triangular) might have some very large entries. Numerically, I found that for a matrix $A$ of dimension $n\times n$, with its entries chosen uniformly at random from $\{0,\pm1\}$, the maximum entry of $U$ behaves roughly as $$ \log \max(U) \approx 2(n-10) \log 2, \qquad n\geq 20, $$ so for large dimensions, the straightforward approach doesn't work so well.

Second, if you don't care all that much about getting a specific posterior distribution over the unimodular matrices, there is a very simple algorithm. Let $U$ be an upper-triangular matrix with ones on the diagonal and whose above-diagonal entries are chosen uniformly at random from an integer interval $I$, with maximum magnitude $a = \max\limits_{x\in I}|x|$.

Clearly, $\det U=1$. So, if you form a matrix $A = U_1 U_2^\top$, it will have determinant $1$, and its entries will be at most $(n-1)^2 a^2$. However, its entries will not be uniformly distributed; for instance, the top-left entry will always be $1$. It is also possible to take matrices such as $U_1 U_2^\top U_3 U_4^\top$, which will still have an uneven distribution of entries, but a little less so.

To satisfy an upper bound on matrix elements, it should be enough to pick a sufficiently small value of $a$, and possibly generate several matrices until one of them has sufficiently small elements.

Third, another way is to let $J_k$ be the identity matrix with the nonzero entries in its $k$-th column replaced with random integers. Then $\det J_k=1$ and the matrix $A=J_1J_2\cdots J_n$ will have $\det A=1$, although different columns will have different distributions.

If you wish to have determinant $-1$, choose randomly an odd permutation with a matrix $P$, $\det P=-1$, and use $AP$.

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .