Home > Analytics, Optimization > Generating random unimodular matrices with a column of ones

Generating random unimodular matrices with a column of ones

In an effort to broaden my blog readership, today’s topic is generating random unimodular matrices with a column of ones. Here’s the tweet from a fellow Hawkeye that got me thinking about this:

A unimodular matrix is a matrix whose determinant is +/-1, denoted det(M) for a matrix M. Two important facts are that det(A)det(B)=det(AB), and the determinant of a triangular matrix is the product of the diagonal entries.

So randomly generating a unimodular matrix is easy:

  • Generate a lower triangular matrix L and upper triangular matrix U with 1s on the diagonal and random entries elsewhere.
    The problem is that the product M=LU does not have a column of ones. But we can fix this. Let’s assume that we want the last column of M should contain all 1s. Then the inner product of row I of L and the last column of U should be 1. We can fiddle with the last element of U in this inner product so that this works out: changing this element will not change the determinant. There is a special case for the lower right element because we do not want to change the diagonals, so in that case we should adjust L instead.

Here is the code in python, assuming the NumPy package has been loaded. I do not use “dot” to compute val so I don’t have to special case i=0. Too long for a tweet, and matrices that are produced are a little funny looking, and if we get unlucky we may divide but zero, but good enough for a blog post.

def rand_unimod(n):
    l = tril(random.randint(-10, 10, size=(n,n))).astype('float')
    u = triu(random.randint(-10, 10, size=(n,n))).astype('float')
    for i in range(0, n):
        l[i, i] = u[i, i] = 1.0
        if i < n - 1:
            val = sum([l[i, j] * u[j, n-1] for j in range(0, i)])
            u[i, n-1] = (1 - val) / l[i, i]
        else:
            val = sum([l[i, j] * u[j, n-1] for j in range(1, i+1)])
            l[n-1, 0] = (1 - val) / u[0, n-1]
    return dot(l, u), l, u
About these ads
  1. Jon Lee
    14 September 2013 at 11:58 pm

    *totally* unimodular

    • 15 September 2013 at 12:01 am

      Yes! Post to be updated over the weekend…

  2. Rongqin Sheng
    24 September 2013 at 2:42 am

    I have a formula to share. Suppose that the 1st column contains all 1’s. Then the matrix can be expressed as

    [ 1 x^T ]
    [ -------------------- ]
    [ e ex^T + B ]

    where e is a vector of all 1’s, x is a random vector of integers, and B is a random unimodular matrix.

  3. Rongqin Sheng
    24 September 2013 at 3:11 am

    Oops! it’s not TUM for the formula. It’s interesting.

  1. No trackbacks yet.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 40 other followers

%d bloggers like this: