There are some posts explaining how do we maximize a function of two parameters with respect to one of those and how to plot the resulting function: Plot a function after taking the supremum with respect to one variable.

Now I have a quite complicated function which has four parameters: $(n,t,\theta,\phi)$ where $n\in \mathbb{N}$ is natural, $t\in (-1,1)$ and $(\theta,\phi)\in S^2$ are angles on the sphere.

The function is highly non trivial. It is the output of computations involving matrices, eigenvalues, and a bunch of things like that. The $n$ is the dimension of the matrix, which acts as a cuttoff to get an approximation to the real computation with operators on a Hilbert space (so one is actually truncating to dimension $n$).

Now I need to maximize with respect to $(\theta,\phi)$ and get a function of $(n,t)$ that for each $n$ fixed can be plotted.

I believe that numerically there sould be no issue in doing that. The method of the question didn't work however. Calling preClassicalCorrelationsAB[nmax_Integer, t_, th_, phi_] the four-argument function, I tried:

f[nmax_Integer, t_?NumericQ] := NMaximize[{preClassicalCorrelationsAB[nmax, t, th, phi], 0 <= th <= Pi, 0 <= phi <= 2 Pi}, t]

And tried plotting f for fixed nmax = 100. It didn't work. In particular I got the error:

0.05102712269051777` is not a valid variable.

The same happend with MaxValue.

I mean, I know the function is highly complicated and so on, but numerically it should be possible to handle it.

So, how can I proceed to maximize this function with respect to $(\theta,\phi)$ and get a function of $(n,t)$ which I can plot for any value of $n$?

Edit: The comment showed I was using NMaximize in the wrong way. Still, I have a problem due to the kind of function I'm considering. The code defining is as follows. First two helper functions are defined:

cs = Compile[{{x, _Real}}, If[0. < x < 1., x Log[2, x], 0.], CompilationTarget -> "C", RuntimeAttributes -> {Listable}, Parallelization -> True]
myEigenvalues[mat_SparseArray?SquareMatrixQ] := If[mat["NonzeroPositions"] === {} && mat["Background"] == 0., ConstantArray[0., Length[mat]], Eigenvalues[N[mat], Method -> "Banded"]]

Next one defines the matrices rhoBPostAPlus and rhoBPostAMinus. Finally one defines:

probPlus[nmax_Integer, t_, th_, phi_] :=Tr[rhoBPostAPlus[nmax, t, th, phi]]
probMinus[nmax_Integer, t_, th_, phi_] := Tr[rhoBPostAMinus[nmax, t, th, phi]]
rhoBPostAPlusNormalized[nmax_Integer, t_, th_, phi_] := rhoBPostAPlus[nmax, t, th, phi]/probPlus[nmax, t, th, phi]
rhoBPostAMinusNormalized[nmax_Integer, t_, th_, phi_] := rhoBPostAMinus[nmax, t, th, phi]/probMinus[nmax, t, th, phi]

entropyBPostAPlus[nmax_Integer, t_, th_, phi_] :=Total[cs@myEigenvalues[rhoBPostAPlusNormalized[nmax, t, th, phi]]]
entropyBPostAMinus[nmax_Integer, t_, th_, phi_] := Total[cs@myEigenvalues[rhoBPostAMinusNormalized[nmax, t, th, phi]]]

preClassicalCorrelationsAB[nmax_Integer, t_, th_, phi_] := entropyB[nmax, t] + probPlus[nmax, t, th, phi]*entropyBPostAPlus[nmax, t, th, phi] + probMinus[nmax, t, th, phi]*entropyBPostAMinus[nmax, t, th, phi]

This preClassicalCorrelationsAB works perfectly if I try to evaluate it. We can even plot it without error.

Now when we evaluate:

f[nmax_Integer, t_?NumericQ] := NMaximize[{preClassicalCorrelationsAB[nmax, t, th, phi],0 <= th <= Pi, 0 <= phi <= 2 Pi}, {th, phi}]

we get the error:

NMaxValue::nnum: The function value -1.90526-0.5 
{<<31>>}}],Method->Banded]]) is not a number at {phi,th} =

Why is that? NMaximize can't handle a function like that?

  • 3
    $\begingroup$ I would think that your call to NMaximize within the definition of f should read: NMaximize[{func, cons}, {th, phi}], i.e. you should NOT be trying to maximize with respect to t, but with respect to $\theta$ and $\phi$. t is a numerical value by the time the NMaximize call is evaluated. $\endgroup$
    – MarcoB
    Mar 12, 2019 at 17:55
  • $\begingroup$ What is rhoBPostAPlus [ ] in your code? The function is undefined. Same for entropyB and rhoBPostAMinus. $\endgroup$
    – MikeY
    Mar 13, 2019 at 16:56
  • $\begingroup$ rhoBPostAPlus is a matrix which depends on the four parameters. I didn't add it to avoid making the post even lengthier, since its definition is rather complicated. The same for rhoBPostAMinus. Intuitively I thought the problem couldn't be with them, since the preClassicalCorrelationsAB function works on its own, being even possible to plot it. $\endgroup$ Mar 13, 2019 at 17:03
  • 1
    $\begingroup$ Try forcing the preClassicalCorrelationsAB function to not run until it has numeric values, so use preClassicalCorrelationsAB[nmax_Integer, t_?NumericQ, th_?NumericQ, phi_?NumericQ] $\endgroup$
    – MikeY
    Mar 13, 2019 at 17:12
  • $\begingroup$ If you can simulate the behavior with a simplified rhoBPostAPlus we can probably solve this in minutes. $\endgroup$
    – MikeY
    Mar 13, 2019 at 17:47

1 Answer 1


Define a function to play with...

preClassicalCorrelationsAB[nmax_, t_, th_, phi_] := -(nmax^2 + t^2 + th^2 + phi^2)

Then define

f[nmax_Integer, t_?NumericQ] := NMaximize[preClassicalCorrelationsAB[nmax, t, th, phi], {th, phi}] // First

We are optimizing over th and phi with n and t fixed, then picking off the value of the function with the First call.

Try it.

f[2, 2]


EDIT: just saw @MarcoB's comment, which is the answer too.

Another EDIT: bringing the answer down from the comments.

When you get these sorts of issues, they can often be addressed by forcing the kernel function of the optimization to only operate with numerical values.

preClassicalCorrelationsAB[nmax_Integer, t_?NumericQ, th_?NumericQ, phi_?NumericQ]

If you are having trouble using a Plot command on the function, and know that the function is reasonably smooth, you can try ListPlot3D@Table[etc.,{ }] as a workaround that lets you control the sampling interval.


Your Answer

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

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