8
$\begingroup$

This is a "common" problem from what I've seen, but with a different spin.

I have a function I use often that finds a fit of an expression to another expression with some free parameters (e.g. for approximation purposes):

ExprSquareError[expr1_, expr2_, vars_, strategy_: Automatic, maxpnts_: Automatic] :=
  Re[NIntegrate @@ Prepend[Join[vars, {Method -> {strategy, "SymbolicProcessing" -> False}, MaxPoints -> maxpnts}], ((expr1) - (expr2))^2]]

(* NMinimize doesn't hold so ExprSquareError will be expanded *)
ExprMinimize[expr1_, expr2_, vars_, params_, constr_: {}, strategy_: Automatic, intstrategy_: Automatic, intmaxpnts_: Automatic] :=
  NMinimize[{ExprSquareError[expr1, expr2, vars, intstrategy, intmaxpnts], constr}, params, Method -> strategy]

Example use:

ExprMinimize[x^2.2, x*x*a + b, {{x, 0, 1}}, {a, b}]

It relies on NIntegrate to find the error between the two expressions with a given assignment of parameters, and NMinimize to find the best choice of parameters to minimize said error.

The problem is that NMinimize doesn't HoldFirst so ExprSquareError is evaluated when you call ExprMinimize, which is fine but generates a warning as NIntegrate tries to expand with some values left unassigned (symbolic) which fails, and in general you get the same issue every time NMinimize wants to evaluate the expression to minimize symbolically for its internal purposes.

This is all fine to me, except that I've found an expression that makes the entire thing crash. Using Monitor[] I see the minimization runs fine but at the end it seems to want to evaluate symbolically again and that last step makes the kernel crash (Tested on Mathematica 8,9 and 10 on Windows)

PointOnPlaneToVector = {x, y, d} \[Function] Normalize[{x, y, d}];
CardKernelFN = {d, x, y, brdf} \[Function] brdf[PointOnPlaneToVector[x, y, d], {0, 0, 1}, {0, 0, 1}];
GGXDh = Function[{HH, NN}, \[Alpha]^2/(Pi*((\[Alpha]^2 - 1)*Max[0, HH . NN]^2 + 1)^2)];
VNormalize = {x} \[Function] x/Norm[x];
GGXD = Function[{LL, VV, NN}, Evaluate[GGXDh[VNormalize[LL + VV], NN]]];
kernelSize = 1;
ggxExprFN = CardKernelFN[d, x, y, GGXD /. {\[Alpha] -> a}];
fitTestParm = {a -> 0.1, d -> 2};

ExprMinimize[ s*ggxExprFN /. {a -> p, d -> 1}, ggxExprFN /. fitTestParm, {{x, -kernelSize, kernelSize}, {y, -kernelSize, kernelSize}}, {p, s}, {(s > 0), (0 <= p <= 1)}, "DifferentialEvolution", "QuasiMonteCarlo", 100]

So now I really have to find a way to restrict the evaluation to be just numerical. If I want to lose the convenience of having ExprMinimize then I can just define the error with pattern matching and make sure everything is NumericQ

minError[p_?NumericQ, s_?NumericQ] = ExprSquareError[s*ggxExprFN /. {a -> p, d -> 1}, ggxExprFN /. fitTestParm, {{x, -kernelSize, kernelSize}, {y, -kernelSize, kernelSize}}, "QuasiMonteCarlo", 100]
NMinimize[{minError[p, s], (s > 0), (0 <= p <= 1)}, {p, s}, Method -> "DifferentialEvolution"]

And this indeed doesn't crash (while still emitting a warning when defining minError as NIntegrate is still being called with symbolic expressions with free parameters).

But the question remains, how could I write a function like ExprMinimize properly? I want the ExprSquareError to be evaluated basically up to NIntegrate, but NIntegrate never be evaluated without fully numeric inputs.

If only it was possible to make pure functions with constrained arguments I could make ExprSquareError return such a function and call it from ExprMinimize:

ExprSquareError[expr1_, expr2_, vars_, params_, strategy_: Automatic, maxpnts_: Automatic] := params \[Function] Evaluate[Re[ NIntegrate @@ Prepend[Join[vars, {Method -> {strategy, "SymbolicProcessing" -> False}, MaxPoints -> maxpnts}], ((expr1) - (expr2))^2]]]
ExprMinimize[expr1_, expr2_, vars_, params_, constr_: {}, strategy_: Automatic, intstrategy_: Automatic, intmaxpnts_: Automatic] := NMinimize[{ExprSquareError[expr1, expr2, vars, params, intstrategy, intmaxpnts] @@ params, constr}, params, Method -> strategy]

But that's not possible, the only way to constrain arguments is with pattern matching it seems. On the other hand, as far as I know, there isn't a way to define a pattern match programmatically, in ExprMinimize I don't know the number or names of the free parameters so I'd need to locally define something like fn[params]=... but that's not possible.

I could maybe try using a CompiledFunction as I think these can be restricted to have only numeric parameters, unlike Function. Any other ideas?

In general Mathematica seems great when you want to use its magic algorithms and prototype, but if you want to execute numerically a specific algorithm on a specific numerical expression it's hard and slow, for example even the fixed method I posted (that doesn't allow for a self-contained ExprMinimize) ends up taking lots of time and memory (I call it many times in a ParallelTable...)

Bonus (malus) point. Even with the "fixed" version using NumeriQ pattern matching the computation ends up running in a second issue, it uses inordinate amounts of memory on the parallel kernels (it always grows, never shrinks) forcing me to generate results in small chunks (quitting kernels every time, save/reload/merge).

E.G.

kernelDistance = 1; kernelSize = 1;
fitDistanceMinMax = {0.1, 5};
fitAlphaMinMax = {0.02, 0.6};

fit2DUnconstrainedError[p_?NumericQ, s_?NumericQ, aa_?NumericQ, dd_?NumericQ] = ExprSquareError[s*ggxExprFN /. {a -> p, d -> 1}, ggxExprFN /. {a -> aa, d -> dd}, {{x, -kernelSize, kernelSize}, {y, -kernelSize, kernelSize}}, "QuasiMonteCarlo", 100]

(* all the code in the comments was tried and didn't help
  ParallelEvaluate[$HistoryLength=1]
  ParallelEvaluate[Share[]]
*)

fitData2DUnconstrained = Flatten[ParallelTable[
  (*Share[]; ClearSystemCache[];*)
  {dd, aa, {p, s} /. Last[NMinimize[{fit2DUnconstrainedError[p, s, aa, dd], (s > 0), (0 <= p <= 1)}, {p, s}, Method -> "DifferentialEvolution"]]},
  {dd, fitDistanceMinMax[[1]], fitDistanceMinMax[[2]], (fitDistanceMinMax[[2]] - fitDistanceMinMax[[1]])/10},
  {aa, fitAlphaMinMax[[1]], fitAlphaMinMax[[2]], (fitAlphaMinMax[[2]] - fitAlphaMinMax[[1]])/10}]
, 1]

P.S. sorry if the code isn't as minimal as possible, I'm sure by bisecting further and experimenting I could have found a more minimal example of the problem, but it's a bit annoying also because everytime Mathematica crashes the kernel the license server takes a while to release the license...

Improve this question
$\endgroup$

1 Answer 1

Reset to default
4
$\begingroup$

I added an extra definition for exprSquareError that will only fire when the parameters are numeric. I redefined exprMinimize to call this version of exprSquareError. Thus, exprMinimize's call of exprSquareError does not expand until the parameters are numeric.

exprSquareError[expr1_, expr2_, vars_, strategy_: Automatic, maxpnts_: Automatic] := Re[NIntegrate @@ Prepend[Join[vars, {Method -> {strategy, "SymbolicProcessing" -> False}, MaxPoints -> maxpnts}], ((expr1) - (expr2))^2]]
exprSquareError[expr1_, expr2_, vars_, {params___?NumericQ}, strategy_: Automatic, maxpnts_: Automatic] := exprSquareError[expr1, expr2, vars, strategy, maxpnts](*note extra definition*)
exprMinimize[expr1_, expr2_, vars_, params_List, constr_: {}, strategy_: Automatic, intstrategy_: Automatic, intmaxpnts_: Automatic] := NMinimize[{exprSquareError[expr1, expr2, vars, params, intstrategy, intmaxpnts], constr}, params, Method -> strategy]

In[191]:= exprMinimize[x^2.2,x*x*a+b,{{x,0,1}},{a,b}]
Out[191]= {0.000129481,{a->0.991587,b->-0.0180289}}

pointOnPlaneToVector = {x, y, d} \[Function] Normalize[{x, y, d}];
cardKernelFN = {d, x, y, brdf} \[Function] brdf[pointOnPlaneToVector[x, y, d], {0, 0, 1}, {0, 0, 1}];
gGXDh = Function[{HH, NN}, \[Alpha]^2/(Pi*((\[Alpha]^2 - 1)*Max[0, HH.NN]^2 + 1)^2)];
vNormalize = {x} \[Function] x/Norm[x];
gGXD = Function[{LL, VV, NN}, Evaluate[gGXDh[vNormalize[LL + VV], NN]]];
kernelSize = 1;
ggxExprFN = cardKernelFN[d, x, y, gGXD /. {\[Alpha] -> a}];
fitTestParm = {a -> 0.1, d -> 2};

In[209]:= exprMinimize[s*ggxExprFN/.{a->p,d->1},ggxExprFN/.fitTestParm,{{x,-kernelSize,kernelSize},{y,-kernelSize,kernelSize}},{p,s},{(s>0),(0<=p<=1)},"DifferentialEvolution","QuasiMonteCarlo",100]
During evaluation of In[209]:= NIntegrate::maxp: The integral failed to converge after 100 integrand evaluations. NIntegrate obtained 170.3335009372765` and 51.646836786706295` for the integral and error estimates. >>
During evaluation of In[209]:= NIntegrate::maxp: The integral failed to converge after 100 integrand evaluations. NIntegrate obtained 104.85984291628401` and 23.239393843835224` for the integral and error estimates. >>
During evaluation of In[209]:= NIntegrate::maxp: The integral failed to converge after 100 integrand evaluations. NIntegrate obtained 174.40018666370233` and 52.26932533157204` for the integral and error estimates. >>
During evaluation of In[209]:= General::stop: Further output of NIntegrate::maxp will be suppressed during this calculation. >>
Out[209]= {0.468343,{p->0.181085,s->3.43887}}
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.