8
$\begingroup$

Bug introduced sometime between 10.0.2 and 11.2 and persisting through 12.3.1 or later

Having just discovered Sum's Method -> "Procedural", I've quickly run into what seems to be a bug in calculating a Hessian matrix of such a sum:

D[Sum[(xx - x[j])^2, {j, nsp}, Method -> "Procedural"], {{xx}, 2}]
(* {{0}} *)

The actual answer should be

(* {{Sum[2, {j, nsp}, Method -> "Procedural"]}} *)

This is strange, because all of the following work:

Non-vectorized second derivative:

D[Sum[(xx - x[j])^2, {j, nsp}, Method -> "Procedural"], {xx, 2}]
(* Sum[2, {j, nsp}, Method -> "Procedural"] *)

Vectorized first derivative:

D[Sum[(xx - x[j])^2, {j, nsp}, Method -> "Procedural"], {{xx}, 1}]
(* {Sum[2 (xx - x[j]), {j, nsp}, Method -> "Procedural"]} *)

Vectorized second derivative (no Method -> "Procedural"):

D[Sum[(xx - x[j])^2, {j, nsp}], {{xx}, 2}]
(* {{2 nsp}} *)

Bug? Any ideas why this happens, and how to work around it?

I should say that the real problem happens within another function, where the function to be differentiated is passed in, and the list of independent variables can be longer than one. Therefore I'm looking for a general solution, not just for this minimal example.

Addition:

Probably the best workaround for me is to just make the Hessian manually with Table:

Table[D[Sum[(xx - x[j])^2, {j, nsp}, Method -> "Procedural"], var1, var2], {var1, {xx}}, {var2, {xx}}]
(* {{Sum[2, {j, nsp}, Method -> "Procedural"]}} *)
$\endgroup$
3
  • $\begingroup$ It seems from a trace that Method -> "X" is parsed as an iterator in the differentiate-the-sum code, which leads to a $Failed being thrown. From there it's a mystery to me why {{0}} is returned. Definitely seems a bug. $\endgroup$
    – Michael E2
    Jul 4 at 0:29
  • $\begingroup$ How frustrating. Does reversing Sum and D help in your more complicated expressions? That allows things to work properly (but maybe more slowly for the more complicated expressions). $\endgroup$
    – JimB
    Jul 4 at 3:00
  • $\begingroup$ @MichaelE2 Thanks, I'll report to WRI. $\endgroup$
    – Chris K
    Jul 4 at 13:29
4
$\begingroup$

The code for D[Sum[..],..] assumes no options to Sum, so Method -> "Procedural" is treated as an iterator. This is a bug. After assuming it's an iterator, the code fails internally because it's a bad iterator. This mysteriously leads to a derivative of {{0}}, which seems an unimportant bug. This last bug happens with D[Sum[(xx - x[j])^2, j], {{xx}, 2}], too, and maybe that makes it an important bug.

A possible workaround would be for the user to remove the options from Sum, differentiate, and then add the options back if differentiating did not remove the sum.

Here's a way to hack the internals to get the right behavior. It basically implements the above procedure, but on the internal sum differentiator. We have to add a down value to catch the case of a Sum with options. (A full fix would include a number of other heads like Inactive[Sum].) We have to add that down value in the right place, before other sums. Maybe it could go first, but I decided it was safer to put it right before the first Sum rule (there are two).

(* save old DownValues *)
olddvs = DownValues@Sum`IterationDDump`iterationD;

newdef = 
  HoldPattern[
     Sum`IterationDDump`iterationD[
      Sum[expr_, its__, sumopts : OptionsPattern[Sum]], vars_, opts_, 
      indspat_]] /; 
    Length[{sumopts}] >= 1 :>
     (Sum`IterationDDump`iterationD[Sum[expr, its], vars, opts, 
       indspat] /. HoldPattern[Sum][args__] :> Sum[args, sumopts]);
newpos = 
 Position[DownValues@Sum`IterationDDump`iterationD, Sum][[1, 1]]

(*  6  *)

DownValues@Sum`IterationDDump`iterationD = 
  Insert[DownValues@Sum`IterationDDump`iterationD,
   newdef,
   newpos
   ];

D[Sum[(xx - x[j])^2, {j, nsp}, Method -> "Procedural"], {{xx}, 2}]

(*  {{2 nsp}}  *)

If you need to fix or change newdef:

(* update definition *)
DownValues@Sum`IterationDDump`iterationD = ReplacePart[
   DownValues@Sum`IterationDDump`iterationD,
   newpos -> newdef
   ];

To go back to the original code, restart the kernel or restore:

(* restore old DownValues *)
DownValues@Sum`IterationDDump`iterationD = olddvs

The above internal hack could be done with WithCleanup in V12.2+ or Internal`InheritedBlock.

$\endgroup$
1
  • $\begingroup$ Wow, that's quite a hack! $\endgroup$
    – Chris K
    Jul 4 at 13:27

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.