9
$\begingroup$

First I present the specific programming question, and then I provide further background details.

Given an expression in Mathematica, say f=3x+7y, there are two common ways to substitute information for x and y. The first is to use the /. rule, so something like

 t=f/.{x->0,y->0};

would make the value of t equal to 0.

And the other is a more general substitution for x and y everywhere they appear, such as

 x=0; y=0;
 t=f;

This also makes the value of t equal to 0. The difference in memory usage, however, becomes apparent with more complicated polynomials. I have uploaded a text file with a complicated polynomial to MediaFire : http://www.mediafire.com/?0ea3mhtqmc8c3r2 (it is about 4MB, and I'm not aware of a better way of sharing it). If you download it and then load it into Mathematica as I will demonstrate, you can see for yourself the significant memory loss:

Substitution Method

 In[1]:= << "/home/username/Downloads/p";
 In[2]:= MaxMemoryUsed[]
 Out[2]= 79833424
 In[3]:= MemoryInUse[]
 Out[3]= 76373112
 In[4]:= t = p /. {x -> 0, y -> 0, z -> 0};
 In[5]:= MaxMemoryUsed[]
 Out[5]= 152124264
 In[6]:= MemoryInUse[]
 Out[6]= 112876080

General Assignment Method

 In[1]:= << "/home/username/Downloads/p";
 In[2]:= MaxMemoryUsed[]
 Out[2]= 79831632
 In[3]:= MemoryInUse[]
 Out[3]= 75415960
 In[4]:= x = 0; y = 0; z = 0;
 In[5]:= MaxMemoryUsed[]
 Out[5]= 79831632
 In[6]:= MemoryInUse[]
 Out[6]= 75420744
 In[7]:= t = p;
 In[8]:= MaxMemoryUsed[]
 Out[8]= 112888664
 In[9]:= MemoryInUse[]
 Out[9]= 112875336

As you can see, during the substitution operation there was about 40MB of additional memory used to perform the substitution. The substitution method also takes much longer than the general method. The difference in the two methods becomes much higher (several GB) as the complexity of the polynomial increases, and ends up maxing out the available RAM. The problem is that for my program I cannot set x, y, z equal to zero in general.

This seems to be an inefficiency with the way Mathematica performs substitutions. I would appreciate any suggestions for working around this problem.

Further Background

The reason I cannot simply employ the general method is that the next step requires that I take derivatives of the original polynomial with respect to the variables x, y, z. The process would ideally be iterated as follows:

Step 1: Given p, write the result of p/.{x->0,y->0,z->0} to a file.

Step 2: Take derivatives of p with respect to x, y, z.

Step 3: Perform some algebraic operations with the derivatives, creating the more complicated polynomials p1 and p2.

Step 4: Repeat the process on each of p1 and p2.

Since the process runs out of RAM when using the substitution rule, I have been working around this by performing several iterations of the process without writing the output of p/.{x->0,y->0,z->0} to a file, and then applying the "General Assignment Method" (setting x=0; y=0; z=0; and writing the evaluation of each polynomial to a file).

This is an unsatisfactory workaround, however, because it is impossible to know the proper number of iterations which will be required, and so it is possible that after performing all the iterations and writing the results to a file that I will need to run it again, this time performing more iterations. If the results could be written out to a file in the intermediate steps, I would be able to check them and determine if I should let the iteration continue or if I can stop it.

Improve this question
$\endgroup$
1
  • $\begingroup$ It's worth noting that both the Clear and Block solutions do not work when the substitution involves an indexed variable, like "a[8]". $\endgroup$
    – mboratko
    Commented Mar 8, 2012 at 17:56

3 Answers 3

Reset to default
10
$\begingroup$

While I am not able to explain why there is a difference between the two methods, I would like to suggest using

t = Block[{x=0, y=0}, p];

Block will set a value to x and y only temporarily, and avoids the use of Clear.

Some more things that can be useful for you:

  1. When you start working, use $HistoryLength = 0 to avoid remembering past results and wasting memory (note that Mathematica sometimes remembers full results even when the input was ended in a semicolon!)

  2. Use Share[] occasionally. It seems that with your polynomial this will help a lot more than in other cases. After reading in the polynomial, it frees 60 MB for me immediately.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thanks for the tip. Much cleaner than my solution. I had two separate questions about the Share[] command - do you know if you have to run it more than once if you have created new complicated expressions since you originally ran it? Also, do algorithms (say the GroebnerBasis algorithm) take more time when you have used the Share[] command because they have to, in a sense, reconstruct the information? $\endgroup$
    – mboratko
    Commented Mar 8, 2012 at 2:53
  • $\begingroup$ @process91 As I understand, Share[] optimizes the way expressions are stored. Every time a new expression is created, there will be a chance that Share[] can reduce memory usage. According to the documentation, Share[] might improve performance. There is no need to reconstruct anything: the information is all there except duplicate parts of expressions are not stored twice. I do not have much practical experience with the performance effects of Share[] because I don't usually work with expressions where it has such a significant effect, like in your case. $\endgroup$
    – Szabolcs
    Commented Mar 8, 2012 at 5:43
4
$\begingroup$

The basic issue is that when using:

expr /. x->1

Mathematica will need to examine the entire expression tree to determine what parts can be replaced. After replacing all of these parts, it then evaluates the result. Contrast this with using a Block approach:

Block[{x=1}, expr]

In this case, Mathematica sets x=1 and then evaluates expr. When x is encountered during the evaluation of expr, it evaluates to 1. In contrast to using ReplaceAll, Mathematica does not look for parts that match x. So, it is not surprising that the Block approach can be much faster.

Here are three more observations:

Evaluation

With the Block approach, x only gets replaced with 1 when the evaluator encounters x. Here is an example where in the Block approach the evaluator never sees x:

Hold[x] /. x->1
Block[{x=1}, Hold[x]]

Hold[1]

Hold[x]

Rules with patterns

If the LHS of the rule is a pattern and not a literal, then the naive Block approach won't work.

{x[1], x[2.2]} /. x[_Integer] -> 1
Block[{x[_Integer]=1}, {x[1], x[2.2]}]

{1, x[2.2]}

Block::lvset: Local variable specification {x[_Integer]=1} contains x[_Integer]=1, which is an assignment to x[_Integer]; only assignments to symbols are allowed.

Block[{x[_Integer] = 1}, {x[1], x[2.2]}]

The Block approach can be made to work by doing something like:

Block[{x}, x[_Integer]=1; {x[1], x[2.2]}]

{1, x[2.2]}

but one should know what they are doing, as differences in how DownValues and ReplaceAll/Replace work, as well as possible evaluation issues complicate things.

Replace

The Rule approach can be sped up by only doing replacements at the desired level. For example:

poly = Expand[(x+y+z)^100];
r1 = poly /. {x->1, y->2, z->3}; //RepeatedTiming
r2 = Replace[poly, {x->1, y->2, z->3}, {-1}]; //RepeatedTiming
r3 = Block[{x=1, y=2, z=3}, poly]; //RepeatedTiming

{0.024, Null}

{0.015, Null}

{0.0099, Null}

By using Replace with level {-1}, only the leaves are replaced, while using ReplaceAll, all parts of the expression are replaced, if possible.

Improve this answer
$\endgroup$
3
$\begingroup$

I have solved this problem, in a very simple way, as I will demonstrate here:

 In[1]:= f = 3 x + 7 y
 Out[1]= 3 x + 7 y
 In[2]:= x = 0;
 In[3]:= t = f
 Out[3]= 7 y
 In[4]:= Clear[x];
 In[5]:= t
 Out[5]= 7 y
 In[6]:= f
 Out[6]= 3 x + 7 y

Now I only wonder why Mathematica's substitute command does not use this method for substitutions. I'm sure there are situations where more complicated substitution rules may fail with this approach, but in that case some heuristics should determine if this much faster and less RAM intense method would be applicable and, if so, apply it.

Improve this answer
$\endgroup$
3
  • $\begingroup$ The method here is "dangerous" in that you cannot readily tell, without examining previous input, that the variable x has (at least temporarily) been given a value. Using ReplaceAll (/.)with Rule (->) immediately makes the substitution obvious. Moreover, a number of functions such as Solve and DSolve return results as rules rather than as values, hence the common need for ReplaceAll. $\endgroup$
    – murray
    Commented Mar 8, 2012 at 19:51
  • $\begingroup$ Bear in mind that Mathematica is a term-rewriting system and direct substitutions work because the input is evaluated as a program. Thus, you ask here for a heuristic that can efficiently and in general determine whether two different programs (the one evaluated directly, and the one obtained using ReplaceAll to perform substitutions and then evaluated) are equivalent. $\endgroup$
    – Oleksandr R.
    Commented Mar 8, 2012 at 20:15
  • $\begingroup$ @OleksandrR. I would think that simply checking if the substitutions requested were for constants would be adequate. Perhaps I'm overlooking something, but from what I can tell using Block[{x=c1,y=c2,...},f] where the c1,c2,... are constants will always yield the same result as f/.{x->c1,y->c2,...}, and will do so faster and with less RAM required. Thus a simple form of the heuristic I'm looking for would be to see if the substitutions were all constants and, if so, use the Block method. $\endgroup$
    – mboratko
    Commented Mar 17, 2012 at 18:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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