Mathematica syntax, fishing expedition

Developers
1

I wanted to pull some values from mathematica for testing gradients for 3F2 and related functions. It’s been years since I’ve used it and the best I could come up with was

Hypergeometric3F2[a1_,a2_,a3_,b1_,b2_,z_]:=HypergeometricPFQ[{a1,a2,a3},{b1,b2},z]
NumberForm[Hypergeometric3F2[1.1, 1.2, 1.3, 4.1, 4.2, 1.0] - Hypergeometric3F2[1.1, 1.2, 1.3, 4.1, 4.2, 1.0-10^-10],20]

Anybody know the syntax for the more generic Grad function for evaluating the gradient at a point or something that automatically does decent finite differences in Mathematica (or one of the open source tools)? All I could get was a promise for the function and it never let me substitute in actual values.

2

Not sure I understand, you want to compute

d/dz(z->Hypergeometric3F2[a1,a2,a3,b1,b2,z])

The Mathematica syntax for that is D[ f, z ]

For example:

In[1]:= D[Hypergeometric3F2[a1, a2, a3, b1, b2, z], z]

Out[1]= (a1 a2 a3 HypergeometricPFQ[{1 + a1, 1 + a2, 1 + a3}, {1 + b1, 1 + b2}, z])/(b1 b2)

For gradient, the fastest way is to use the Grad build-in function

Grad[f[x, y, z], {x, y, z}]

will return the list of partial derivatives

{ D[ f[x,y,z],x ], D[ f[x,y,z],y ], D[ f[x,y,z],z ] }

I am using Mathematica and I would be happy to help if I can

-Vincent

3

I wanted Grad[f[x, y, z], {x, y, z}] evaluated at a specific set of points but then I realized there were some tests for grad_F32 that I hadn’t noticed before. I’ll figure that out and check back here. Thanks for offering to help.

4

Ok, so the tests in test/unit/math/mix/scal/fun/grad_F32_test.cpp do test some of the gradients of grad_F32, but the values of grad_F32 are never tested. Our stan::math::grad_F32 is supposed to give us—pretty sure about this—the partial derivative of F32 w.r.t each of the six input parameters at their current values. So here’s what I did, line 1 is just defining 3F2, lines 2 & 3 are the finite difference gradient for the z parameter. Then lines 4 & 5 are just demonstrating my confusion about the Mathematica syntax for Grad. Does this help? I need the gradient at some points (I have a set that exercises all the branches of the logic in grad_F32) for all six parameters.

Hypergeometric3F2[a1_,a2_,a3_,b1_,b2_,z_]:=HypergeometricPFQ[{a1,a2,a3},{b1,b2},z]
NumberForm[(Hypergeometric3F2[1.1, 1.2, 1.3, 4.1, 4.2,1.0] - 
  Hypergeometric3F2[1.1, 1.2, 1.3, 4.1, 4.2, 1.0-10^-10])/10^-10,20]

This one just returns more or less a statement of the input:

Grad[Hypergeometric3F2[a1,a2,a3,b1,b2,z], {1.1,1.2,1.3,4.1,4.2,1.0}]

This one returns a vector of zeros:

Grad[Hypergeometric3F2[ 1.1,1.2,1.3,4.1,4.2,0.1],{a1,a2,a3,b1,b2,z}]
5

Ok, thank you. I should have guessed that from your first post.
I think that what you want is something like:

gradAtPoint[f_, atPoint_, precision_: MachinePrecision] /;
VectorQ[atPoint, NumberQ] :=

Block[{atPointPrec, var, placeHolder, grad, gradValue},
atPointPrec =
SetPrecision[atPoint, Max[precision, Precision[atPoint]]];
var = ReplacePart[
ConstantArray[atPointPrec, Length@atPoint], {i_, i_} →
placeHolder];
grad = Map[D[f @@ #, placeHolder] &, var];
gradValue =
Table[grad[[i]] /. placeHolder → atPointPrec[[i]], {i,
Length@atPoint}];
Return[gradValue];
]

For example:

gradAtPoint[Hypergeometric3F2, {1.1, 1.2, 1.3, 4.1, 4.2, 1.0}]

gives you (computed with MachinePrecision)

{0.136454, 0.125884, 0.116885, -0.0399597, -0.0390586, 0.17696}

For arbitrary precision, here 30 digits (slower):

gradAtPoint[Hypergeometric3F2, {1.1, 1.2, 1.3, 4.1, 4.2, 1.0},30]

gives you

{0.1364541527605578885150865535281,
0.1258842738952212691686044022322, 0.1168846455204336835620604842355,
-0.03995972722956410382059540493864,
-0.03905863436957377012633711359080, 0.17695990197399216889220148097}

Hope this solved your problem
-Vincent

6

I also think that it is possible to write Mathematica scripts to automatically generates c++ unit tests code to check special function gradient values.
Just tell me if you are interested in such kind of approaches.
-Vincent

7

No need for guessing, this is perfect, thank you!

8

It would be really cool to use something like that to improve our test coverage in the math library more systematically but I’m unlikely to have the time to put into it for months into the future (I have a few other Stan projects in mind once my survival model woes are solved). I wonder if this is something a student could do as a discrete project (not sure if the level of guidance required would be worth the returns there).

9

I was thinking of a possible use of Mathematica to generate c++ test codes. The advantage is that Mathematica can also provide the numerical values to check. However code creation is basically string manipulation, perhaps not the most interesting student project. (An alternative would be to call directly the C++ function to check from Mathematica, but this requires Mathematica installed which is not a practical solution here).
I can investigate this for a peculiar special function just to see what happens: can someone suggest a critical/problematic one to check?

10

Yesterday I quickly provided a Mathematica code. However this code is ugly and overcomplicated, I recommend to use this cleaned new version.

gradAtPoint[f_, atPoint_, precision_: MachinePrecision] /;
VectorQ[atPoint, NumberQ] := Block[{var, grad},
var = Table [Unique, {Length@atPoint}];
grad = Map[D[f @@ var, #] &, var] /.
Thread[var →
SetPrecision[atPoint, Max[precision, Precision[atPoint]]]];
Return[grad];
]