partial numerical evaluation

Question

Context

The following is an issue I have come across many times and still do not have a satisfactory solution for.

The context of it is this. Let's say we have an equation, a linear ODE for the function f of the variable x:

equation = c2[x] f''[x] + c1[x] f'[x] + c0[x] f[x] + d0[x] == 0

each of the coefficients depends on x, but they can also depend on another function(s), let's say g[x], or on other parameter(s), let's say p.

We want to solve this equation numerically, therefore we discretize the variable x onto some grid, let's say

grid = {1.,2.,3.,4.};

The derivatives will of course become matrices, but this question is only about the coefficients, which become vectors.

At this point an example is probably useful. The equation above is not important, we can focus on one of the coefficients, which could for example take the form:

c = (Sinh[1/2 (g x^4 + x^2)] + p^3 x/ (g - x^3)) g^2 x^2;

(Note: I hope this is a sufficiently complicated example, but any simplifying structure here is unintended, the dependence on the variable x, the function(s) g and the parameters p can be arbitrary)

Of course x is known ahead of time (grid above), but p and g are not. Imagine a sequence of computations running many times in a loop, that contains the pieces:

(* ... *)
p = 3.; (* actually some computation *)
g = {0.714, 0.429,0.734, 0.921}; (* some other computation *)
cComputed = Block[{x = grid}, c]
(* use this vector c to construct and solve the discretized equation *)

{-47.6352, 78.8935, 1.77501*10^15, 3.19207*10^55}

Done this way, at every pass through the loop everything in c is computed from scratch. However since we know x ahead of time, we can already do a big part of the computations:

cPartlyComputed = g^2({1.`, 4.`, 9.`, 16.`} Sinh[{0.5`, 8.`, 40.5`, 128.`} g + {0.5`, 2.`, 4.5`, 8.`}] 
+ p^3  {1.`, 8.`, 27.`, 64.`}/(g - {1.`, 8.`, 27.`, 64.`}));

Note that I have not merely inserted {1.,2.,3.,4.} for x, but also computed x^2, 1/2 x^2 and 1/2 x^4.

speed comparison

A quick test to show that this gives a significant speedup:

without = Function[{g, p, x}, g^2 x^2 ((p^3 x)/(g - x^3) + Sinh[1/2 (x^2 + g x^4)])];
with = Function[{g, p}, g^2 ({1.`, 4.`, 9.`, 16.`} Sinh[g {0.5`, 8.`, 40.5`, 128.`} + {0.5`, 2.`, 4.5`, 8.`}] + (p^3 {1.`, 8.`, 27.`, 64.`})/(g + {-1.`, -8.`, -27.`, -64.`}))];

gval = RandomReal[{}, {4}];
pval = 3.;

Do[without[gval, pval, grid], {10^5}] // AbsoluteTiming//Print;
Do[with[gval, pval], {10^5}] // AbsoluteTiming//Print;

1.35754

0.719002

About a factor two speedup.

The question

This is the question: how can we create the expression cPartlyComputed above automatically?

naive solution

The naive solution illustrates the problem:

cTry = c /. x -> grid;

{1. g^2 ((1. p^3)/(-1. + g) + Sinh[1/2 (1. + 1. g)]), 4. g^2 ((2. p^3)/(-8. + g) + Sinh[1/2 (4. + 16. g)]), 9. g^2 ((3. p^3)/(-27. + g) + Sinh[1/2 (9. + 81. g)]), 16. g^2 ((4. p^3)/(-64. + g) + Sinh[1/2 (16. + 256. g)])}

The unknown vector g has entered into the elements of a list. When we have computed g and insert it into this, we will get a matrix, not the vector that we want.

Best attempt so far

Since the essence of the problem seems to be the Listable property of Plus and Times when applied to symbolic arguments, the idea was to temporarily remove this property, only for symbolic arguments, or create our own Plus and Times identical to the built-in ones except for this property.

I tried the latter:

ClearAll[CircleTimes, CirclePlus]
Attributes[CirclePlus] = Attributes[CircleTimes] = Complement[Attributes[Plus], {Listable, Protected}];
Verbatim[CircleTimes][a_] := a;
Verbatim[CirclePlus][a_] := a;
numericOrList[a_] := NumericQ[a] || Head[a] === List;
a_\[CircleTimes]b_ := a*b /; numericOrList[a] && numericOrList[b];
a_\[CirclePlus]b_ := a + b /; numericOrList[a] && numericOrList[b];
a_\[CircleTimes](b_\[CirclePlus]c_) := 
a\[CircleTimes]b\[CirclePlus]a\[CircleTimes]c /; numericOrList[a]
a_\[CircleTimes]x1_\[CirclePlus]a_\[CircleTimes]x2_ := 
        a\[CircleTimes](x1 + x2) /; numericOrList[x1] && numericOrList[x2];


evaluate[expression_, var_, grid_] := Hold[Evaluate[
expression /. {Plus -> CirclePlus, Times -> CircleTimes} /. var -> grid
]] /. {CirclePlus -> Plus, CircleTimes -> Times}

some explanation:

CircleTimes and CirclePlus are undefined built-in functions that display nicely. As mentioned we give them all the properties of Times and Plus except Listable (and Protected of course).

We check their arguments with numericOrList which is True for pure numbers or for lists, which in our context will always be lists of pure numbers. On such arguments we do want the Listable property, so all the definitions are to ensure that they reduce to Times and Plus in these cases.

The Verbatim is necessary to prevent an infinite recursion, see Carl Woll's answer to this question.

Finally the function evaluate uses these CircleTimes and CirclePlus together with the naive answer above, holding the result ready to be used in a later computation.

performance

This seems to work:

evaluate[c,x,grid]

Hold[g^2 ({1., 4., 9., 16.} Sinh[ g {0.5, 8., 40.5, 128.} + {0.5, 2., 4.5, 8.}] + ( p^3 {1., 8., 27., 64.})/(g + {-1., -8., -27., -64.}))]

identical to the cPartlyComputed above that I did by hand.

However it is prohibitively slow, presumably because of all the pattern matching. On this other still very simple example:

c2=-5832. - 432. x^2 - 1728. f x^3 - 12. (32. + 72. d1f + 144. f^2 - 27. p^4) x^4
    - 24. f (32. + 72. d1f - 27. p^4) x^5

it creates something that seems optimal:

evaluate[c2,x,grid]

Hold[d1f {-864., -13824., -69984., -221184.} + f (d1f {-1728., -55296., -419904., -1.76947*10^6} + p^4 {648., 20736., 157464., 663552.} + {-768., -24576., -186624., -786432.}) + f {-1728., -13824., -46656., -110592.} + f^2 {-1728., -27648., -139968., -442368.} + p^4 {324., 5184., 26244., 82944.} + {-6648., -13704., -40824., -111048.}]

but takes 34 seconds to do so on my laptop (I don't understand why, the pattern matching is probably slow but that slow?) This is so slow that I cannot even test if it would at least work in a more general case.

Requirements of the solution

So what I am looking for is this:

1. A function evaluate[expression_,var_,grid_] that replaces the symbolic variable var with the numerical list grid in expression, without messing it up by taking any other symbols in expression into the lists.
2. expression can be either a single expression or a (nested) list of expressions, each one being an arbitrary nonlinear function of one or several symbols. A more realistic example on which it should work is given below.
3. The result should contain no, or as few as possible, operations that can be carried out with just the knowledge of var.
4. While for a given problem this has to be done only once and the result will be used many times, it is still important that it be reasonably fast.
5. If the equation, and the functions and parameters, are at machine precision, then inserting packed arrays for the functions should ideally result in a packed array as output as well.
6. Ideally the code should be short and easy to understand.

It doesn't have to be related to my attempt above, perhaps there is a way to "hack" the built in Times and Plus to behave like this in a much faster way, or something else entirely.

An actual example to test on

I've uploaded a more realistic case on which it should work. It depends on the functions f,g,h and their derivatives, which have been rewritten as f'[x] -> d1f etc( if it's easier this can of course be undone) and on the scalar parameter p.

Update

There have been two answers so far, by Carl Woll and xzczd, which are both close to what I'm looking for but not quite. A quick comparison on the example linked above: with a grid size of 4, Carl Woll's solution takes about 0.013 seconds to produce something with a LeafCount of 4573, while xzczd's answer takes about 0.25 seconds to produce something with a LeafCount of 3409*. So the former is a lot faster but also computes less. Computations with the result of the latter are then about 25% faster than with the former.

However both leave computations undone. With Carl Woll's answer I think it's because the rules aren't repeated, resulting in expressions like f {-1728., -13824., -46656., -110592.} + f {0., -51104., -61502., -3.270*^6} that are not evaluated. (it is also becoming harder to read the code) With xzczd's answer the problem is that additions of numbers or a number and a list aren't computed, because for some reason including that it becomes incredibly slow.

Therefore I will set a bounty for a solution that does do really all the computations in a reasonable time (I don't think it should take more than 0.5 seconds for the example linked above, possibly much shorter). And if it's still easily readable a year from now that would be great!

[* Actually to obtain this I've added to his last two patterns also the ones with patt2 instead of patt]

Answer 1

You can set x to be numeric, and then search for numeric subexpressions to be replaced. For example:

evaluate[params_, expr_, v_, grid_] := Module[{x},
    NumericQ[x]=True;
    ReplaceAll[
        Activate @ Inactive[Function][params, expr /. v->x],
        h_?NumericQ :> RuleCondition[h /. x->grid]
    ]
]

For your first example:

c = (Sinh[1/2 (g x^4 + x^2)] + p^3 x/ (g - x^3)) g^2 x^2;
evaluate[{g, p}, c, x, {1, 2, 3, 4}]

Function[{g, p}, g^2 {1, 4, 9, 16} ((p^3 {1, 2, 3, 4})/(g + {-1, -8, -27, -64}) + Sinh[1/2 ({1, 4, 9, 16} + g {1, 16, 81, 256})])]

For your second example:

c2 = -5832.-432. x^2-1728. f x^3-12. (32.+72. d1f+144. f^2-27. p^4) x^4 -24. f (32.+72. d1f-27. p^4) x^5;
evaluate[{f, d1f, p}, c2, x, {1, 2, 3, 4}]

Function[{f, d1f, p}, -5832. + {-432., -1728., -3888., -6912.} - 1728. f {1, 8, 27, 64} - 12. (32. + 72. d1f + 144. f^2 - 27. p^4) {1, 16, 81, 256} - 24. f (32. + 72. d1f - 27. p^4) {1, 32, 243, 1024}]

You'll notice that multiplicands and addends that are both numeric don't necessarily get combined in this approach, but multiplying and adding packed arrays is really fast, and is unlikely to be where most of the computational time is spent.

Still, it is not too difficult to combine these multiplicands and addends, and the above code could be modified to do so:

evaluate[params_, expr_, v_, grid_] := Module[{x, foo, tmp, rule1, rule2},
    SetAttributes[foo, NumericFunction];
    NumericQ[x]=True;
    rule1 = RuleDelayed[
        HoldPattern[Plus[h__]],
        Total @ KeyValueMap[
            If[#1, foo[Plus@@#2], Plus@@(#2/.rule1)]&,
            GroupBy[{h},NumericQ]
        ]
    ];
    rule2 = h_?NumericQ z_:1 :> If[z===1, foo[h], foo[h] ReplaceAll[rule2] @ z];
    Activate[
        Inactive[Function][params, expr] /.
            v->x /.
            rule1 /.
            rule2
    ] /. foo[z_] :> RuleCondition[z /. {foo -> Identity, x -> grid}]
]

Your second example:

c2 = -5832.-432. x^2-1728. f x^3-12. (32.+72. d1f+144. f^2-27. p^4) x^4-24. f (32.+72. d1f-27. p^4) x^5;
evaluate[{f,d1f,p}, c2, x, {1,2,3,4}]

Function[{f, d1f, p}, f {-1728., -13824., -46656., -110592.} + (p^4 (-27.) + 32. + d1f 72. + f^2 144.) {-12., -192., -972., -3072.} + f (p^4 (-27.) + 32. + d1f 72.) {-24., -768., -5832., -24576.} + {-6264., -7560., -9720., -12744.}]

Update

My first couple answers traverse the expression in a top-down approach, isolating the x-dependent parts along the way. Your CirclePlus/CircleTimes approach (as well as @xzxcd) traverse the expression in a bottom-up approach.

Now, there are two things that my approach didn't do that you wanted:

1. Convert $f(x) (a + b + \ldots)$ to $a f(x) + b f(x) + \ldots$

2. Convert $f(x) a + g(x) a + h(x) b + \dots$ to $(f(x) + g(x)) a + h(x) b + \ldots$

This is difficult to do in a single-pass top-down approach, but straightforward in single-pass bottom-up approach. So, rather than processing the expression with two or more passes, I will switch to using a bottom-up approach as well.

Now, the reason for the slowness of the addition rules in the other rules is because they rely on the flat and orderless attributes of CirclePlus. Consider the following definition:

ClearAll[CirclePlus]
SetAttributes[CirclePlus, {Flat,Orderless}]

a_?NumericQ ⊕ b_?NumericQ := a+b

Let's count how many times NumericQ is called when evaluating 1 ⊕ x ⊕ 2

i = 0;
TracePrint[1 ⊕ x ⊕ 2,_NumericQ, TraceAction :> (i++&)]
i

3⊕x

70

We see that NumericQ is called 70 times! This is why including such binary rules causes the other answers to be so slow. Rather than using binary definitions and relying on the Flat and Orderless attributes, it is much faster to use definitions that don't rely on these attributes. Here is an example:

Clear[CirclePlus]

Verbatim[CirclePlus][a__] /; Count[{a}, _?NumericQ]>1 := CirclePlus @@ GroupBy[{a}, NumericQ, Apply[Plus]]

Then:

i = 0;
TracePrint[CirclePlus[1, x, 2], _NumericQ, TraceAction -> (i++&)]
i

3⊕x

7

Much better. Note that the excessive NumericQ calls are exacerbated by the Orderless attribute. I will first present a broad overview of my updated answer.

1. Replace Times and Plus with CircleTimes and CirclePlus.

2. Have CircleTimes separate numeric parts (assuming x is numeric) from non-numeric parts, wrapping the numeric parts in a n wrapper, and the non-numeric parts in a s wrapper.

3. Define distributive and product rules for the n wrapper.

4. Have CirclePlus factor non-numeric parts.

Here is the code:

ClearAll[CircleTimes, CirclePlus, nWrapper]

evaluate[params_, expr_, v_, grid_] := Module[{x, res},
    NumericQ[x]=True;
    res = Inactive[Function][params, expr] /. v->x /. {Times->CircleTimes, Plus->CirclePlus};
    ReplaceAll[
        Activate[res /. sWrapper -> Identity],
        {
        nWrapper[a_] :> RuleCondition @ ReplaceAll[a, {x -> grid, nWrapper -> Identity}],
        n_?NumericQ :> RuleCondition @ ReplaceAll[n, x->grid]
        }
    ]
]

CircleTimes[a__] := Times @@ KeyValueMap[
    If[#, nWrapper[#2], sWrapper[#2]]&,
    GroupBy[{a}, NumericQ, Apply[Times]]
]

CirclePlus[a__] := Plus @@ KeyValueMap[
    Times[#1, nWrapper[#2]]&,
    GroupBy[
        Flatten @ Replace[{a}, Verbatim[Plus][b__] :> {b}, {1}],
        sPart -> nPart,
        Apply[Plus]
    ]
]

sPart[_?NumericQ] := 1
sPart[u_sWrapper] := u
sPart[t_Times] := Replace[t, _?NumericQ->1, {1}]
sPart[u_] := sWrapper[u]

nPart[t_?NumericQ] := nWrapper[t]
nPart[u_sWrapper] := 1
nPart[t_Times] := Replace[t, Except[_?NumericQ]->1, {1}]
nPart[u_] := 1

SetAttributes[nWrapper, {Flat, NumericFunction}]
nWrapper[1] = 1;
nWrapper /: f_nWrapper a_Plus := Distribute @ Unevaluated[f a]
nWrapper /: f_nWrapper sWrapper[a_Plus b_.] := Distribute @ Unevaluated[f a] If[b===1, 1, sWrapper[b]]
nWrapper /: Verbatim[nWrapper][f_] Verbatim[nWrapper][g_] := nWrapper[f g]

In your examples, numbers and lists now distribute, producing your expected outputs:

c = (Sinh[1/2 (g x^4+x^2)]+p^3 x/(g-x^3)) g^2 x^2;
evaluate[{g,p}, c, x, {1, 2, 3, 4}]

Function[{g, p}, g^2 ((p^3 {1, 8, 27, 64})/( g + {-1, -8, -27, -64}) + {1, 4, 9, 16} Sinh[{1/2, 2, 9/2, 8} + g {1/2, 8, 81/2, 128}])]

c2 = -5832.-432. x^2-1728. f x^3-12. (32.+72. d1f+144. f^2-27. p^4) x^4 -24. f (32.+72. d1f-27. p^4) x^5;
evaluate[{f, d1f, p}, c2, x, {1, 2, 3, 4}]

Function[{f, d1f, p}, f {-1728., -13824., -46656., -110592.} + f^2 {-1728., -27648., -139968., -442368.} + d1f {-864., -13824., -69984., -221184.} + p^4 {324., 5184., 26244., 82944.} + f (d1f {-1728., -55296., -419904., -1.76947*10^6} + {-768., -24576., \ -186624., -786432.} + p^4 {648., 20736., 157464., 663552.}) + {-6648., -13704., -40824., -111048.}]

Here's a view of the "collected" expression before x is replaced with a list:

evaluate[{f, d1f, p}, c2, x, x]

Function[{f, d1f, p}, f (-1728. x^3) + f^2 (-1728. x^4) + d1f (-864. x^4) + p^4 (324. x^4) + f (d1f (-1728. x^5) - 768. x^5 + p^4 (648. x^5)) + (-5832. - 432. x^2 - 384. x^4)]

Notice how all of the x-dependent coefficients are grouped together. Basically, the function evaluate acts like Collect, but without restructuring. For instance, compare the above to:

evaluate[{f, d1f, p}, Expand @ c2, x, x]

Function[{f, d1f, p}, f^2 (-1728. x^4) + d1f (-864. x^4) + p^4 (324. x^4) + d1f f (-1728. x^5) + f p^4 (648. x^5) + (-5832. - 432. x^2 - 384. x^4) + f (-1728. x^3 - 768. x^5)]

The function evaluate does not do any factoring or expanding of the parameters. On the other hand, if you used Collect, you would get a restructuring dependent on the order of parameters used:

Activate[
    Inactive[Function][{f, d1f, p}, Collect[c2, {f, d1f, p}, w]]
] /. w[s_] -> s

Function[{f, d1f, p}, f^2 (-1728. x^4) + d1f (-864. x^4) + p^4 (324. x^4) + (-5832. - 432. x^2 - 384. x^4) + f (d1f (-1728. x^5) + p^4 (648. x^5) + (-1728. x^3 - 768. x^5))]

vs:

Activate[
    Inactive[Function][{f, d1f, p}, Collect[c2, {p, d1f, f}, w]]
] /. w[s_] -> s

Function[{f, d1f, p}, f^2 (-1728. x^4) + d1f (-864. x^4 + f (-1728. x^5)) + p^4 (324. x^4 + f (648. x^5)) + (-5832. - 432. x^2 - 384. x^4) + f (-1728. x^3 - 768. x^5)]

However, for your big example, the kind of restructuring done by Collect would be bad, as it will expand all powers, creating many more terms to evaluate.

At any rate, let's see how evaluate does on your big example:

evaluate[
    {d1f,d1g,d1h,d2f,d2g,d2h,f,g,h,p},
    exampleCoefficients,
    x,
    {1,2,3,4}
]; //AbsoluteTiming

{0.048297, Null}

Answer 2

The matching for certain pattern e.g. 1⊕{1, 2} seems to be extremely slow, so I exclude them from the matcher:

ClearAll[CircleTimes, CirclePlus]
Attributes[CirclePlus] = 
  Attributes[CircleTimes] = Complement[Attributes[Plus], {Listable, Protected}];
Verbatim[CircleTimes][a_] := a;
Verbatim[CirclePlus][a_] := a;

With[{patt =(*__?NumericQ|*){__?NumericQ}, patt2 = __?NumericQ}, 
 (a : patt)⊗(b : patt) := a b;
 (a : patt)⊗(b : patt2) := a b;
 (a : patt2)⊗(b : patt2) := a b;
 (a : patt)⊕(b : patt) := a + b;
 (*(a:patt)⊕(b:patt2):=a+b;*)
 (a : patt)⊗(b_⊕c_) := a⊗b⊕a⊗c;
 a_⊗(x1 : patt)⊕a_⊗(x2 : patt) := a⊗(x1 + x2);]

evaluate[expression_, var_, grid_] := 
 Hold[Evaluate[
    expression /. {Plus -> CirclePlus, Times -> CircleTimes} /. 
     var -> grid]] /. {CirclePlus -> Plus, CircleTimes -> Times}

test = evaluate[c2, x, grid]; // AbsoluteTiming
(* {0., Null} *)

Still, the result is very close to the optimized one.

Answer 3

You may use the Trott-Strzebonski in-place evaluation technique which I first learned from here: 19542.

With

c = (Sinh[1/2 (g x^4 + x^2)] + p^3 x/(g - x^3)) g^2 x^2;
grid = {1., 2., 3., 4.};

Then

cPartlyComputed =
 Function[{g, p}, Evaluate[c]] /.
  {
   f_[x, n_] :> With[{v = f[grid, n]}, v /; True],
   x :> With[{v = grid}, v /; True]
   }

Function as attribute HoldAll so it can be treated as a held expression. f_[x, n_] is getting the Power[x, _] entries and x picks up any left over x'es.

Hope this helps.

Answer 4

Here is a different attempt to work around the Listable attribute of Plus and Times.

doPE[c_, x_, grid_] :=
  Internal`InheritedBlock[{Plus, Times},
    ClearAttributes[{Plus, Times}, {Flat, Listable}];
    Hold @@ {
      c /. x -> grid //.
       {
         a___ * (n_?NumericQ | n_List) * L_List :> a * Thread[n * L],
         a___ + (n_?NumericQ | n_List) + L_List :> a + Thread[n + L],
         a___ + f_*L_List + f_*M_List :> a + f Thread[L + M]
       }
     }
  ]

Test:

c = (Sinh[1/2 (g x^4 + x^2)] + p^3 x/(g - x^3)) g^2 x^2;
c2 = -5832. - 432. x^2 - 1728. f x^3 - 
   12. (32. + 72. d1f + 144. f^2 - 27. p^4) x^4 - 
   24. f (32. + 72. d1f - 27. p^4) x^5;
grid = {1., 2., 3., 4.};

doPE[c, x, grid]

doPE[c2, x, grid]

Hold[g^2 {1., 4., 9., 16.}
  ((p^3 {1., 2., 3., 4.})/(g + {-1., -8., -27., -64.}) + 
    Sinh[1/2 ({1., 4., 9., 16.} + g {1., 16., 81., 256.})])]

Hold[{-6264., -7560., -9720., -12744.} +
  f {-1728., -13824., -46656., -110592.} +
  f (32. + 72. d1f - 27. p^4) {-24., -768., -5832., -24576.} +
  (32. + 72. d1f + 144. f^2 - 27. p^4) {-12., -192., -972., -3072.}]