Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Let me begin by saying that this is my first post on StackExchange. I apologize in advance if I unwittingly break any of its unwritten rules of etiquette.

Recently, I've been trying to understand an algorithm by coding it into Mathematica. One of the main steps in the algorithm requires me to find a integer-valued row echelon form for a matrix with integer entries. There are several constraints that the input matrices satisfy. The matrices represent sparse underdetermined (never inconsistent or trivial) homogeneous linear systems. The number of columns may only take particular values fixed by the algorithm but, for the purposes of this question, it is enough to know that the number of columns is always much, much larger than the number of rows. The number of rows is also typically somewhat large and it may be the case that some number of rows are duplicates or linear combinations of other rows. Another key point is that the integer entries could generically be quite large. Therefore Compile was not an option for me (however if there is some magic way to do multi-precision arithmetic with Compile please let me know!).

My main motivation for asking this question is that this is the first time that it has actually mattered to me how efficient my Mathematica code is and I am having a very difficult time optimizing my solution attempts. First a word about data structures. I represent my sparse arrays using their rule-based form on a row-by-row basis. As a representative example, consider the following 15 X 100 matrix:

{{{1} -> 1, {21} -> -1, {62} -> -1, {92} -> 1}, {{2} -> 
1, {22} -> -1, {64} -> -1, {94} -> 1}, {{3} -> 
1, {23} -> -1, {65} -> -1, {95} -> 1}, {{2} -> 
1, {22} -> -1, {64} -> -1, {94} -> 1}, {{4} -> 
1, {24} -> -1, {67} -> -1, {97} -> 1}, {{5} -> 
1, {25} -> -1, {68} -> -1, {98} -> 1}, {{3} -> 
1, {23} -> -1, {65} -> -1, {95} -> 1}, {{5} -> 
1, {25} -> -1, {68} -> -1, {98} -> 1}, {{6} -> 
1, {26} -> -1, {69} -> -1, {99} -> 1}, {{55} -> -2, {56} -> 
1987, {62} -> -8, {63} -> -2, {65} -> -4, {66} -> 
1987, {72} -> -4, {73} -> -2, {82} -> -4, {95} -> 
2}, {{58} -> -2, {59} -> 
1987, {64} -> -8, {65} -> -2, {68} -> -4, {69} -> 
1987, {74} -> -4, {75} -> -2, {84} -> -4, {98} -> 
2}, {{59} -> -2, {60} -> 
1987, {65} -> -8, {66} -> -2, {69} -> -4, {70} -> 
1987, {75} -> -4, {76} -> -2, {85} -> -4, {99} -> 
2}, {{5} -> -2, {6} -> 1987, {13} -> -2, {16} -> 
1987, {22} -> -4, {23} -> -2, {32} -> 4, {45} -> -2, {55} -> 
4, {56} -> -3974, {62} -> 16, {63} -> 4, {64} -> -8, {65} -> 
4, {66} -> -3974, {72} -> 8, {73} -> 4, {82} -> 8, {94} -> 
8}, {{8} -> -2, {9} -> 1987, {15} -> -2, {19} -> 
1987, {24} -> -4, {25} -> -2, {34} -> 4, {48} -> -2, {58} -> 
4, {59} -> -3974, {64} -> 16, {65} -> 4, {67} -> -8, {68} -> 
4, {69} -> -3974, {74} -> 8, {75} -> 4, {84} -> 8, {97} -> 
8}, {{9} -> -2, {10} -> 1987, {16} -> -2, {20} -> 
1987, {25} -> -4, {26} -> -2, {35} -> 4, {49} -> -2, {59} -> 
4, {60} -> -3974, {65} -> 16, {66} -> 4, {68} -> -8, {69} -> 
4, {70} -> -3974, {75} -> 8, {76} -> 4, {85} -> 8, {98} -> 8}}

My main motivation for using this representation of my data is that it seems to make the implementation of an efficient pivoting strategy incredibly simple. I'm not overly wedded to this choice of data structure and would be happy to switch to something else if it allows for a better solution to the problem.

One obvious question is: Why not just use built-in functions? It seems to me that, at least naively, the out-of-the-box solution based on RowReduce

IntegerRowEchelon1[L_, n_] := (temp = 
DeleteCases[RowReduce[SparseArray[#, n] & /@ L], 0 Range[n]]; 
Drop[ArrayRules[#], -1] & /@ 
Expand[(LCM @@ Denominator[Flatten[temp]]) temp])

is overkill in the extreme. First of all, it seems impossible to force RowReduce to deliver a row echelon as opposed to a reduced row echelon form. The back-solving step of the elimination process is unnecessary for me and therefore should be avoided. Furthermore, since I am working over the integers, I have to waste time by clearing the denominators and turning the reduced row echelon form that RowReduce delivers into an acceptable integer-valued output. So far, however, I haven't been able to improve upon the above solution.

The essence of my question to you, the experts, is whether or not I should just go with the ridiculous work-around IntegerRowEchelon1 or if there is a more natural lower-level solution that would work as fast or faster. There are clearly many computations done by IntegerRowEchelon1 that could in principle be avoided but, on the other hand, RowReduce has undoubtedly been highly optimized by the developers at WRI. So it is not obvious to me whether or not one can do significantly better.

Finally, let me share with you two solution attempts which lose to IntegerRowEchelon1 by factors of roughly 2 and 3 respectively. I tested them both on 1000 15 X 100 matrices with 35% of their entries filled in using RandomInteger. The correctness of the output was verified by applying IntegerRowEchelon1 to both the input and the output. These solution attempts were based on a classical algorithm due to Camille Jordan. Since this algorithm was not widely known until the work of Erwin Bareiss circa 1968 (The original paper, E. H. Bareiss, "Sylvester's Identity and Multistep Integer-Preserving Gaussian Elimination", is googleable (search "Erwin Bareiss Sylvester's Identity", first link) if you want to have a look at it. I also was able to get quite a bit out of the Google Books preview for Chapter 9 of "Algorithms For Computer Algebra" by Geddes et. al. http://books.google.es/books?id=9fOUwkkRxT4C&printsec=frontcover&hl=es&source=gbs_ge_summary_r&cad=0#v=onepage&q&f=false), it is usually called the Bareiss one-step fraction-free Gaussian elimination algorithm. This algorithm is exposited in many textbooks on computer algebra and is probably the simplest one that does what I want.

Although I couldn't find an answer in the documentation, I would guess that this algorithm is what Mathematica is using to forward-solve when you select the method "OneStepRowReduction." Needless to say, I tried forcing RowReduce to use this method but it appears to make no difference whatsoever, in either the timing or in the output. Of course, I've also considered the possibility that the reason I can't beat Mathematica is that it is using an algorithm more efficient than Bareiss's. For now, I will proceed assuming that Bareiss's algorithm is at least no worse than the algorithm(s) upon which RowReduce is ultimately based. For those of you with insider information, any pearls of wisdom you can share about how RowReduce actually works would be more than welcome.

My first attempt at coding Bareiss's algorithm looks pretty ugly but only loses to IntegerRowEchelon1 by a factor of about 1.9:

IntegerRowEchelon2[L_] := (Aold = SortBy[L, #[[1, 1, 1]] &]; 
rnum = Length[Aold]; i = 1; lastpivot = 1; Anew = Aold; flag = 0; 
While[i <= rnum - 1, 
If[Aold[[i, 1, 1]] =!= Aold[[i + 1, 1, 1]], i++; lastpivot = 1; 
Continue[]]; 
For[a = i + 1, a <= rnum, a++, 
If[Aold[[i, All, 1]] =!= Aold[[a, All, 1]], 
 tempRRi = 
  Rule[#, 0] & /@ 
   Complement[Aold[[a, All, 1]], Aold[[i, All, 1]]]; 
 tempRRa = 
  Rule[#, 0] & /@ 
   Complement[Aold[[i, All, 1]], Aold[[a, All, 1]]]; 
 Aold[[a]] = Sort[Join[Aold[[a]], tempRRa]]; 
 Aold[[i]] = Sort[Join[Aold[[i]], tempRRi]]; 
 Anew[[i]] = Aold[[i]]; 
 Anew[[a]] = 
  If[Aold[[a, 1, 1]] === Aold[[i, 1, 1]], Drop[Aold[[a]], 1], 
   Aold[[a]]], 
 Anew[[a]] = 
  If[Aold[[a, 1, 1]] === Aold[[i, 1, 1]], Drop[Aold[[a]], 1], 
   Aold[[a]]]]; 
Which[Aold[[a, 1, 1]] === Aold[[i, 1, 1]], 
 Anew[[a, 1 ;; Length[Aold[[a]]] - 1, 2]] = 
  Drop[1/lastpivot (Aold[[i, 1, 2]] Aold[[a, 
         1 ;; Length[Aold[[a]]], 2]] - 
      Aold[[i, 1 ;; Length[Aold[[a]]], 2]] Aold[[a, 1, 2]]), 1], 
 True, 
 Anew[[a, 1 ;; Length[Aold[[a]]], 2]] = 
  1/lastpivot Aold[[i, 1, 2]] Aold[[a, 1 ;; Length[Aold[[a]]], 
     2]]]; If[
 Anew[[a, 1 ;; Length[Aold[[a]]] - 1, 2]] =!= 
  0 Range[Length[Aold[[a]]] - 1], 
 Anew[[{i, a}]] = DeleteCases[#, {x_} -> 0] & /@ Anew[[{i, a}]]; 
 Aold[[{i, a}]] = DeleteCases[#, {x_} -> 0] & /@ Aold[[{i, a}]], 
 flag = 1; Anew = Delete[Anew, a]; rnum--; Break[]]]; 
If[flag == 1, flag = 0; If[i != Anew[[i + 1, 1, 1, 1]], i++]; 
 Aold = Anew; Continue[]]; lastpivot = Aold[[i, 1, 2]]; i++; 
 Anew = SortBy[Anew, #[[1, 1, 1]] &]; Aold = Anew]; Anew)

I was sure that my second attempt was going to perform much better but it turned out to lose by a factor of about 2.8 to IntegerRowEchelon1

IntegerRowEchelon3[L_] :=
(A = SortBy[L /. Rule[{x_}, y_] :> List[x, y], #[[1, 1]] &]; 
rnum = Length[A]; i = 1; lastpivot = 1; 
While[i < rnum, 
If[A[[i, 1, 1]] != A[[i + 1, 1, 1]], i++; lastpivot = 1; 
Continue[]]; 
pivotrow = A[[i]]; (activerow = A[[#]]; 
  A[[#]] = 
   If[Length[#] == 1, {#[[1, 1]], #[[1, 2]]/
       lastpivot}, {#[[1, 1]], (#[[1, 2]] + #[[2, 2]])/
       lastpivot}] & /@ 
    SplitBy[Sort[
      Join[Map[{#[[1]], pivotrow[[1, 2]] #[[2]]} &, 
        If[activerow[[1, 1]] == pivotrow[[1, 1]], 
         Drop[activerow, 1], activerow]], 
       Map[{#[[1]], -If[activerow[[1, 1]] == pivotrow[[1, 1]], 
             activerow[[1, 2]], 0] #[[2]]} &, 
        Drop[pivotrow, 1]]]], #[[1]] &]) & /@ Range[i + 1, rnum]; 
A = SortBy[
 DeleteCases[DeleteCases[A, {_, 0}, {2}], {}], #[[1, 1]] &]; 
 lastpivot = pivotrow[[1, 2]]; rnum = Length[A]; i++]; 
A /. List[x_, y_] :> Rule[{x}, y])

At least to me,IntegerRowEchelon3 looks much nicer than IntegerRowEchelon2. But it's way slower! This is quite counter-intuitive to me because, in coding IntegerRowEchelon3, I made a concerted effort to get rid of as much fluff as possible and be as functional as possible which I thought was "good Mathematica programming style". It was very disappointing to discover that IntegerRowEchelon3 loses so badly to IntegerRowEchelon1. I must admit that I feel a bit lost at this point about issues of performance tuning and look forward to seeing alternative solutions to my problem. To be clear, at the end of the day speed is what matters and I would be more than happy to accept "high-level" solutions that beat IntegerRowEchelon1 as well.

share|improve this question
1  
Hi! If you don't mind, could you please link to the book or article describing the algorithm you say you're trying to implement? –  J. M. Jun 1 '12 at 16:33
    
There are some unmatched brackets in your definition of IntegerRowEchelon1. –  Heike Jun 1 '12 at 16:43
    
Thanks for the help! Better now? –  Rob2181 Jun 1 '12 at 16:55
add comment

1 Answer 1

up vote 14 down vote accepted

Actually what you want is HermiteDecomposition. It is the integer ring form of RowReduce; that latter, while working largely over the integers for the forward elimination, actually is an echelon form over the rational field.

As for efficiency, you'll have to experiment to see if it meets your needs. If not, feel free to post or send examples so I'll have something concrete to test.

--- edit ---

Here is some testing. I defined your input to be m1, and do not repeat show that definition below. First I turn it into a bona fide mathematica SparseArray object.

m2 = Flatten[MapIndexed[{#2[[1]], #1[[1, 1]]} -> #1[[2]] &, m1, {2}], 1];
mat = SparseArray[m2];

Now we get the HNF.

In[455]:= Timing[{uu, hh} = HermiteDecomposition[mat];]

Out[455]= {0.01, Null}

We have the normal form (as uu) and the conversion matrix (as hh). I'll show the former.

In[457]:= uu

Out[457]= {{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 
  0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 1, 
  0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, -1988, 
  0, -1, -994}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -2, 0, 0, -1}, {0, 
  0, 0, 0, 0, -2, 0, 0, 1987, -2, 0, 0, -1, 0, 0}, {0, 0, 0, 0, 0, 0, 
  0, 0, 0, -1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0, -1, -994, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 
  0}, {0, -1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, -1, 0, 0,
   0, 1, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 
  0, 0, 0, 0}}

Here is a substantially larger example.

randomSparseIntegerMatrix[m_, n_, p_, max_] := 
 RandomChoice[{p, 1 - p} -> {1, 0}, {m, n}]*
  RandomInteger[{-max, max}, {m, n}]

I create a 50x250 matrix with inputs up to 10^8 in magnitude, and sparsity of around 80%.

biggermat = randomSparseIntegerMatrix[50, 250, 1/5, 10^20];

In[465]:= Timing[{uubig, hhbig} = HermiteDecomposition[biggermat];]

Out[465]= {1.68, Null}

Hope that helps to give some idea of possibilities for using HermiteDecomposition.

--- end edit ---

share|improve this answer
    
thanks so much for the prompt response! I didn't know about HermiteDecomposition and will have a look at it. Actually, the example 15 X 100 matrix in my post, while in a somewhat unorthodox form, is quite representative of what I'm working with (apart from the fact that the integers in it are a bit on the smaller side) and represents a good test case. Of course you'd probably need to copy it a 100 times or so to get a timing above the machine's timing resolution but that's easy enough to do. –  Rob2181 Jun 1 '12 at 17:37
    
See edit for examples. –  Daniel Lichtblau Jun 1 '12 at 18:33
    
I tested Daniel's solution in the same way that I tested IntegerRowEchelon2 and IntegerRowEchelon3. I found that it outperformed IntegerRowEchelon1 by a factor of about 1.2. I guess this means that the extra operations performed by IntegerRowEchelon1 are actually quite inexpensive and I shouldn't have lost any sleep over them. –  Rob2181 Jun 2 '12 at 13:55
    
However, I still think that it would be nice if there was a switch in RowReduce that would instruct the function to try and deliver a row echelon form with entries in the domain that the input matrix had entries in. The reason for this is that the Hermite normal form approach, to the best of my knowledge, will only work for principal ideal domains. On the other hand, something like Bareiss's algorithm should work in any integral domain since all one needs is a notation of exact division. To give an example, I would expect Hermite normal form to work for matrices with univariate polynomial –  Rob2181 Jun 2 '12 at 13:59
    
entries but NOT matrices with multivariate polynomial entries. Bareiss's algorithm should work fine for any polynomial ring. –  Rob2181 Jun 2 '12 at 14:00
show 3 more comments

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.