Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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...

share|improve this question

This question has an open bounty worth +50 reputation from Angelo Pesce ending in 4 days.

This question has not received enough attention.

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.