How to maximize just with respect two parameters out of four?

Question

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 
((0<Eigenvalues[SparseArray[Automatic,<<3>>],Method->Banded]<1.)+1.4427
Eigenvalues[SparseArray[<<1>>],<<1>>]
Log[Eigenvalues[<<1>>,Method->Banded]])-0.5
((0<Eigenvalues[SparseArray[Automatic,{11,11},0,
{1,{{0,2,5,8,11,14,17,20,23,26,29,31},{{2},{1},{1},{3},{2},{2},{4},{3},{3},{5},{4},{4},{6},{5},{5},{7},{6},{6},{8},{7},{7},{9},{8},{8},{10},{9},{9},{11},{10},{10},{11}}},
{Times[<<4>>],Times[<<2>>],Times[<<4>>],Times[<<4>>],Times[<<2>>],Times[<<4>>],Times[<<4>>],Times[<<2>>],Times[<<4>>],Times[<<4>>],Times[<<2>>],Times[<<4>>],Times[<<4>>],Times[<<2>>],Times[<<4>>],Times[<<4>>],Times[<<2>>],Times[<<4>>],Times[<<4>>],Times[<<2>>],Times[<<4>>],Times[<<4>>],Times[<<2>>],Times[<<4>>],Times[<<4>>],Times[<<2>>],Times[<<4>>],Times[<<4>>],Times[<<2>>],Times[<<4>>],Times[<<2>>]}}],Method->Banded]<1.)+1.4427
Eigenvalues[SparseArray[<<1>>],Method->Banded]
Log[Eigenvalues[SparseArray[Automatic,{11,11},0,{1,{<<2>>},
{<<31>>}}],Method->Banded]]) is not a number at {phi,th} =
{6.20181,2.61304}.

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

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]

-8

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.