Recently I tried to do what I thought was a fairly small (relative to the 6 GB of RAM that I have on my machine) row reduction calculation on a matrix representing an undetermined linear system and was very surprised by how much trouble I had with it. The 2590 X 12155 matrix, A (available for download via Mediafire), takes a little over 17 MB to store, has rational entries, and is very sparse (~3% of the entries non-zero). The data structure I used to store A is Mathematica's rule-based representation for sparse arrays on a row-by-row basis (if this is unclear please see the example matrix from my last post).
Curiously, I find that using RowReduce to find a row echelon form for A leaves a huge memory footprint well beyond what Mathematica reports with MemoryInUse. I was able to observe this effect without leaving Mathematica in an (to me miraculous) allegedly platform-independent way by using the answer that @WReach gave here in response to a question about memory monitoring. I began with
In[1] := $HistoryLength = 0;
to turn off caching and then set up a pair of commands which give respectively my machine's and Mathematica's estimates of how much memory is currently in use using the ideas of @WReach:
In[2] := Needs["JLink`"]
InstallJava[];
LoadJavaClass["java.lang.management.ManagementFactory"];
In[5] := temp = JavaBlock[{#,java`lang`management`ManagementFactory
`getOperatingSystemMXBean[]
@#[]} & /@ {getFreePhysicalMemorySize,
getTotalPhysicalMemorySize}]; memstart1 =
temp[[2, 2]] - temp[[1, 2]]
Out[5] := 3096420352
In[6] := memstart2 = MemoryInUse[]
Out[6] := 18286664
In[7] := temp = JavaBlock[{#,
java`lang`management`ManagementFactory`getOperatingSystemMXBean[]
@#[]} & /@ {getFreePhysicalMemorySize,
getTotalPhysicalMemorySize}]; memnow1 =
temp[[2, 2]] - temp[[1, 2]] - memstart1
Out[7] := 94208
In[8] := memnow2 = MemoryInUse[] - memstart2
Out[8] := 1808
I've found that copying and running the preceding two lines of code wherever I want a memory sanity-check seems to give a fairly consistent way to compare how much memory Mathematica has actually used since the start of the session and how much memory Mathematica thinks it has used since the start of the session (the numbers 94208 and 1808 from the initialization above should be thought of as being effectively zero). Although it is probably obvious, I have of course assumed that the amount of RAM used by all other programs is constant throughout the session (which should be approximately correct if they are just running idly in the background).
Now I read in the file containing the matrix and assign it to A:
In[9] := f = OpenRead["rrdemo"];
A = Read[f];
Close[f];
I also go ahead and define an array that will be used momentarily (recall that 12155 is the number of columns in A)
In[12] := W = Array[w, 12155];
Actually, although it's a bit off-topic, it's worth pointing out that there already seems to be some missing memory associated with the process of reading in the file "rrdemo" (in a separate session I checked that there is no missing memory associated with the definition of W):
In[13] := temp = JavaBlock[{#,
java`lang`management`ManagementFactory`getOperatingSystemMXBean[]
@#[]} & /@ {getFreePhysicalMemorySize,
getTotalPhysicalMemorySize}]; memnow1 =
temp[[2, 2]] - temp[[1, 2]] - memstart1
Out[13] := 425725952
In[14] := memnow2 = MemoryInUse[] - memstart2
Out[14] := 268374384
This is actually nothing at all compared to what happens with RowReduce. To obtain a row echelon form for my matrix I use RowReduce in a straightforward way on a sparse array. To show that the problem is not with my generation of the sparse array, I first check the memory footprint of that step of the calculation:
In[15] := AbsoluteTiming[C1 = CoefficientArrays[Plus @@@ (#[[2]] w @@
#[[1]] & /@ # & /@ A), W][[2]];]
Out[15] := {5.447995, Null}
In[16] := temp = JavaBlock[{#,java`lang`management`ManagementFactory
`getOperatingSystemMXBean[]@#[]} & /@ {getFreePhysicalMemorySize,
getTotalPhysicalMemorySize}]; memnow1 =
temp[[2, 2]] - temp[[1, 2]] - memstart1
Out[16] := 539504640
In[17] := memnow2 = MemoryInUse[] - memstart2
Out[17] := 385560664
In[18] := ByteCount[C1]
Out[18] := 100604400
The performance of Mathematica here seems reasonable. The actual and reported memory usage go up from what they were before in proportion to the ByteCount of C1.
Now look what happens if we RowReduce C1:
In[19] := AbsoluteTiming[C2 = Drop[ArrayRules[#], -1] & /@
RowReduce[C1, Method -> "OneStepRowReduction"];]
Out[19] := {402.172945, Null}
In[20] := temp = JavaBlock[{#,
java`lang`management`ManagementFactory`getOperatingSystemMXBean[]
@#[]} & /@ {getFreePhysicalMemorySize,
getTotalPhysicalMemorySize}]; memnow1 =
temp[[2, 2]] - temp[[1, 2]] - memstart1
Out[20] := 2461675520
In[21] := memnow2 = MemoryInUse[] - memstart2
Out[21] := 731499200
In[22] := ByteCount[C2]
Out[22] := 280479264
Something seems to be seriously wrong here. The ByteCount of C2 is only ~ 280 Mb and this is close to the increase in the total memory usage reported by Mathematica. But the actual memory usage is ~3.3 times what Mathematica thinks it is! I'm forced to conclude that Mathematica fails to garbage collect more than 1.5 Gb of memory at the end of the calculation.
In summary, I have a couple of related questions. Does anybody have any insights into the anomalous behavior of RowReduce that I just illustrated? Am I doing something wrong? More generally, does anybody have any idea as to whether or not Mathematica has special algorithms for symbolic sparse linear algebra? Memory issues aside, it makes me suspicious that RowReduce takes so long to process C1 (I have a 2.5 GHz processor) and that it returns a result that is no longer a sparse array. At this point I'm wondering whether Mathematica is even the right tool for the job (symbolic linear algebra on large, very sparse undetermined linear systems). Any advice from more experienced users would be greatly appreciated. If it makes any difference, I'm using 64 bit Linux and am running Mathematica v.8.0.4.
RowReduceonC1twice in a row. $\endgroup$