0
$\begingroup$

I have a bottleneck problem involving a function call with symbolic derivatives.

I don't have much experience with Mathematica (or any dynamic language) and I'm certain that the way I've set this up is about as basic as it gets, so anyone Mathematica -fluent will probably spot several "obvious" changes that should speed things up. Please tell me what these are! There are probably a host of RTFM suggestions for speeding up calls to this function (which I welcome) but specific suggestions would be great too.

The quantity mainTerm[l, m, r, Theta, Phi, FactorA, FactorB] gets called many times in a nested loop of (r,Theta,Phi) coordinates and (l,m) iterations; so this involves calling spherical harmonics/Hankel functions (and calculating their derivatives and double derivatives) ad nauseum. When I expanded out the derivatives/double derivatives down to the level of hypergeometric basis functions (to avoid doing calculus on the fly), this seemed to slow things down, which has me baffled.

The block shown below is mainly to illustrate the complexity of the function I need to compute (i.e., I'm hoping for advice on how to compile or otherwise speed up functions when analytic derivatives might need to be calculated on the fly, not quantitative // AbsoluteTiming comparisons):

k = somevalue; preAAA := FactorA*(SphericalHankelH1[l, k*r]* D[D[SphericalHarmonicY[l,m,Theta,Phi],Theta]*Sin[Theta],Theta] + SphericalHankelH1[l, k*r]*D[D[SphericalHarmonicY[l,m,Theta,Phi],Phi]]); AAA := Expand[preRadial]; BBB := (FactorB*SphericalHankelH1[l, k*r])*D[SphericalHarmonicY[l,m,Theta,Phi],Phi] - FactorA*D[SphericalHarmonicY[l,m,Theta,Phi],Theta]* D[r*SphericalHankelH1[l, k*r],r]; CCC := (FactorB*SphericalHankelH1[l, k*r])*D[SphericalHarmonicY[l,m,Theta,Phi],Theta] - FactorA*D[SphericalHarmonicY[l,m,Theta,Phi],Phi]* D[r*SphericalHankelH1[l, k*r],r]; mainTerm[l_,m_,r_,Theta_,Phi_,FactorA_,FactorB_] = AAA + BBB + CCC;

Subequently, inside the (l, m) loop I set storeValue[[l, m]] = mainTerm[l, m, ...]. When I compile AAA, BBB and CCC individually and set storeValue this way, the code runs slower (which I am confused about).

I'm aware that compiling my component terms could (should?) speed things up, and I've tried various approaches to that end (with/without compiling to C, with/without explicitly defining Reals and Integers, etc.) but none of them had much effect and some slowed things down.

Any suggestions on how better to set this up for a large number of nested function calls?

Improve this question
$\endgroup$

1 Answer 1

Reset to default
3
$\begingroup$

Probably the main reason why your calculation is slow is that you're using SetDelayed (:=) which means that the derivatives are calculated again and again every time a quantity such as PreAAA is called. Instead, just replace all those assignments with Set (=) to evaluate the derivatives symbolically once and for all:

Clear[l, r, m, Theta, Phi, k]; preAAA = FactorA*(SphericalHankelH1[l, k*r]* D[D[SphericalHarmonicY[l, m, Theta, Phi], Theta]*Sin[Theta], Theta] + SphericalHankelH1[l, k*r]* D[D[SphericalHarmonicY[l, m, Theta, Phi], Phi]])

and so on for the other quantities as well. As to your compilation problems, it's not certain without more information, but you probably have some expressions in there that simply cannot be compiled, so that the compilation doesn't actually work fully. See, e.g., here.

Edit

Another speed-up tip would be to thread the special functions over lists of arguments as much as possible. I.e., generate lists of the possible values of l, m, and maybe a grid of the continuous parameters k, r, and the angles for which you'd like to evaluate the expressions. Then use the fact that you can do one-time calls like this to get many function values:

SphericalHankelH1[Range[10], 1.] (* ==> {0.301169 - 1.38177 I, 0.0620351 - 3.60502 I, 0.00900658 - 16.6433 I, 0.00101102 - 112.898 I, 0.0000925612 - 999.44 I, 7.15694*10^-6 - 10880.9 I, 4.79013*10^-7 - 140453. I, 2.8265*10^-8 - 2.09591*10^6 I, 1.49138*10^-9 - 3.549*10^7 I, 7.11655*10^-11 - 6.72215*10^8 I} *)

This is just the basic idea. The list over which you thread can actually have sub-lists too, so one can make pretty complex structures with one single call to a special function.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thanks for pointing out my inappropriate use of SetDelayed (I thought this was necessary to force evaluation of AAA, BBB & CCC with the new arguments in mainTerm). I just made this change, but unfortunately it does not appear to speed things up appreciably (I'm running a very reduced number of iterations to check speed). Is there an obvious way to pre-compile the spherical harmonics/Hankel function derivatives, and could this improve run time? Thanks! $\endgroup$ Commented Mar 21, 2013 at 4:10
  • $\begingroup$ Oh... silly me... seeing as SphericalHarmonicY and SphericalHankelH1 do not appear anywhere on that list, I might just be bum out of luck. There has to be a way to speed this up... The loop iterates at least 40 million times and it is taking aeons. Any ideas would be hugely appreciated. $\endgroup$ Commented Mar 21, 2013 at 4:27
  • $\begingroup$ So you are suggesting pre-populating a 2D lookup matrix (say matrixH1) for the SphericalHankel case threaded over my range of l values and a mesh of r values, then calling matrixH1[[l,krIndex]] in each of my AAA/BBB/CCC contributions to mainTerm. And then obviously repeating the process for preallocating the D[] and D[D[]] terms in each of AAA/BBB/CCC. I tried this last week, but it slowed my algorithm down - I assumed this was something to do with the cost of looking up elements in a 4D matrix? I'll try again, more carefully, and see if it was just my implementation. $\endgroup$ Commented Mar 21, 2013 at 5:35
  • $\begingroup$ Yes, that's pretty much the extent of it. And note that you can feed the list of kr points to these functions as a list too. Which of the arguments to thread over depends on the application, I guess. $\endgroup$ Commented Mar 21, 2013 at 5:48
  • $\begingroup$ Thank you for your suggestions, Jens. The speed of my implementation is still woefully inadequate, but you were very helpful in pointing out the limited scope of Compile[] and confirming the wisdom of pre-populating matrices for reiterated elements. Unrelated: I love your MathematicaGraphics page, and now for an Honours talk I'm trying to reproduce the appearance of your spherical coordinate system plot using your 3D text labels and vector styles. I don't think mine will look any where near as nice, but thank you for this resource! $\endgroup$ Commented Mar 22, 2013 at 2:32

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.