3
$\begingroup$

As a very minimal example, consider a sum of fractions such as:

$A=\frac{a}{s_{123}bc}+\frac{d}{e}$

In practice, I have many hundreds of thousands of these fractions generated from a recursion relation. My desired result comes around from the prescription in my theory, where I am required to multiply through by $s_{123}$ and then take $s_{123}\to 0$. Upon replicating in MMA, I use a code similar to:

A*s[1,2,3]//Expand
%/.{s[1,2,3]->0}

This Expand is necessary since by simply multiplying A by the relevant factor, MMA will first interpret it as $s_{123} \big (\frac{a}{s_{123}bc}+\frac{d}{e}\big)$. I think this is what is slowing my code down substantially, and crashes MMA due to a lack of memory.

To get around this, I have had a few ideas;

  1. Use something like DeleteCase as a simple method (or something perhaps similar to Filter out all terms not involving a given variable which might be a bit overkill for my needs) to delete all terms not involving $s_{123}$, and then set the remaining factors of it in the leftover terms to 1.
  2. Instead of using Expand, define a new function which brings common factors inside to each fraction.

Which of these, alongside any other methods, would you suggest? How could I go about implementing them? I think 2 would be easier to do, but at this point, I am needing efficiency in my code on such a large scale.

$\endgroup$
  • $\begingroup$ You might find that working with a List of terms is better (e.g. effectively List@@A) than working with a sum of terms. Presumably Expand is comparing a large number of terms, looking for common subexpressions. $\endgroup$ – mikado Mar 2 '19 at 12:15
  • $\begingroup$ I would have thought Cancel[] would work better here than Expand[]. $\endgroup$ – J. M.'s technical difficulties Mar 2 '19 at 20:16
  • 1
    $\begingroup$ Cancel actually crashes my code in this case. It simply runs out of memory on my laptop, but nonetheless; Expand was my first thought, and did the trick at the time when efficiency wasn't as important. $\endgroup$ – Brad Mar 2 '19 at 20:26
4
$\begingroup$

Pick out only those terms in A that contain exactly one s123 in the denominator:

A = a/(s123 b c) + d/e;

Select[A, Exponent[#, s123] == -1 &] /. s123 -> 1

a/(b c)

or with pattern matching using a Default pattern (*):

Total[Cases[A, _./s123, {1}]] /. s123 -> 1

a/(b c)

or with pattern-based case deletion:

DeleteCases[A, Except[_./s123]] /. s123 -> 1

a/(b c)

You'll have to check what's fastest in your case.

Timing comparison

Faking a random input:

SeedRandom[1234];
A = Total[Times @@@ ({a, b, c, s}^# & /@ RandomInteger[{-10, 10}, {10^4, 4}])];

Method 1: Select is quite fast:

First[T1 = RepeatedTiming[Select[A, Exponent[#, s] == -1 &] /. s -> 1]]
(* 0.042 *)

Method 2: Cases is very fast:

First[T2 = RepeatedTiming[Total[Cases[A, _./s, {1}]] /. s -> 1]]
(* 0.006 *)

Method 3: DeleteCases is about equally fast, depending on the exact input:

First[T3 = RepeatedTiming[DeleteCases[A, Except[_./s]] /. s -> 1]]
(* 0.0081 *)

@CarlWoll's method is very slow:

First[T4 = RepeatedTiming[SeriesCoefficient[A, {s, 0, -1}]]]
(* 5.12 *)

check that the results agree:

T1[[2]] == T2[[2]] == T3[[2]] == T4[[2]] // Expand
(* True *)

Footnote

(*) We need a special Default pattern for 1/s123 because it does not match the pattern _/s123. Why not? The full form of x/s123 (with anything for x) has head Times:

x/s123 //FullForm
(* Times[Power[s123, -1], x] *)

The full form of 1/s123, on the other hand, has head Power:

1/s123 //FullForm
(* Power[s123, -1] *)
| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for your very helpful and informative comment. Only off the top of your head, can you think of any other alternative methods that may be faster which I have not considered before? $\endgroup$ – Brad Mar 2 '19 at 15:22
  • 1
    $\begingroup$ The pattern (_ : 1)/s instead of _/s | 1/s works. $\endgroup$ – Somos Mar 2 '19 at 19:44
  • $\begingroup$ I've added DeleteCases as another fast method. $\endgroup$ – Roman Mar 3 '19 at 2:43
  • $\begingroup$ @Somos the default pattern _./s works too and looks even simpler. $\endgroup$ – Roman Mar 15 '19 at 2:44
2
$\begingroup$

Why not just use SeriesCoefficient?

A = a/(s123 b c) + d/e;
SeriesCoefficient[A, {s123, 0, -1}]

a/(b c)

| improve this answer | |
$\endgroup$
  • $\begingroup$ Because it's terribly slow, see my answer. $\endgroup$ – Roman Mar 2 '19 at 16:38
  • 1
    $\begingroup$ In my extended case, it took around 41 seconds, compared to my original method which took around 5. The improvements from Roman both took sub-1 second when applied to the practical case. Thank you for your answer; although it's not quite what I needed in this situation, it serves as a pleasant reminder of the in-built capabilities of MMA. $\endgroup$ – Brad Mar 2 '19 at 16:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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