2

I have a large 4x4 symbolic matrix whose components are scalar functions with arguments (t,r,$\theta$, $\phi$) and various parameters. The LeafCount of the purely symbolic expression is around 500,000 and, after plugging in the values of the parameters, it takes 0.3 seconds to evaluate the matrix on a single co-ord point. The problem is, this function needs to be numerically integrated, so 0.3 seconds for a single point is far too slow. I've simplified all the expression that enter into the derivation of the full symbolic matrix, however executing Simplify or Together on the symbolic matrix itself crashes the kernel (it runs out of memory on my 10Gb-of-ram machine).

Unfortunately, compilation to C or WVM is not feasible since the matrix contains InterpolatingFunctions. I've tried implementing this solution to extract the piecewise polynomial from the InterpolatingFunction, but when trying to compile to C, the kernel crashes from memory usage and when compiling to WVM, the evaluation on a co-ordinate point fails with code

CompiledFunction::cfex: Could not complete external evaluation at instruction 1; proceeding with uncompiled evaluation.

For clarity, I'm compiling the function only after inserting the parameter values so its only a function of (t,r,$\theta$, $\phi$), which are the arguments for the Compile function.

The code that generates this expression is far too long to upload, so I exported it to a .mx file, which can be found here. The parameters are (\[chi], m, n, \[omega], \[nu]) and the arguments are (t,r, \[theta], \[phi]). There are functions S[\[theta]] and R[r[\chi]] and their derivatives within the symbolic expression. Numerical evaluation of the expression can be done using

tmp[t, r, \[Theta], \[Phi]] /. S -> Function[{\[Theta]}, 1] /. 
          R -> Function[{r}, 1] /. \[Omega] -> 1 /. \[Nu] -> 
         1 /. \[Chi] -> 0.5 /. m -> 1 /. n -> 1 /. t -> 0 /. 
   r -> 2 /. \[Theta] -> 1 /. \[Phi] -> 0

for a quick output.

How then can I speed up the numerical evaluation of such a large expression? Is simplification the best route or somehow compiling to WVM or C? Any help is appreciated.

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
shanedrum
  • 577
  • 2
  • 8
  • 5
    I can think of two things, OptimizeExpression and passing the interpolating function values as arguments to the compiled function. – Michael E2 Jun 10 '22 at 15:42
  • 1
    Indeed, I tried to load the matrix into the symbol A. Then LeafCount[Experimental`OptimizeExpression[A]] is only $3655$ instead of $581590$. That tells me that you were quite a bit careless in generating this expression. – Henrik Schumacher Jun 11 '22 at 07:48
  • Thats quite surprising how much it simplifies. The quantity is the result of many tensor contractions in a kerr spacetime, so I would normally expect quite a huge expression, but it seems Mathematica is able to reduce it quite a bit. How would I convert the result from Experimental`OptimizeExpression into a function I can then numerically evaluate? Such as, inserting the parameter values and having a function like Function[{t,r,theta,phi}, result] with which I can evaluate over a co-ordinate grid? – shanedrum Jun 11 '22 at 19:11
  • One idea may be found in hf1 in this answer: myfunc = Experimental`OptimizeExpression[expr] /. Experimental`OptimizedExpression[f_] :> Function @@ Hold[{...args...}, f] – Michael E2 Jun 12 '22 at 22:51

0 Answers0