The goal

Implement some Python code to find the extrema points of a function that is strongly oscillating.

The background

Let $G$ be a connected graph with $n$ points with Laplacian matrix $L(G)$. We may represent each vertex of $G$ as a canonical basis vector $|krangle$ for $kin {1, …, n}$, spanning a position Hilbert space $mathcal{H}$. A notable element of such a space is the flat superposition
$$|srangle=frac{1}{sqrt n}(|1rangle+|2rangle+…|nrangle).$$

Within such a model, assume you are given the Hamiltonian
$$ H(gamma, w, G)=gamma L(G)-|wranglelangle w|$$
where $|wrangle$ is one of the position basis vectors and $gamma$ a real parameter, and consider the probability function $$P_{G,n}(t,gamma)=|langle w| U_H(t,gamma)| srangle|^2$$ for a system prepared in the flat state to evolve under $H$ into $|wrangle$.

---

For a general $G$, $P$ is a strongly oscillating function. Calculating the right maximum for a pair $(gamma, t)$ given $G, n$ and $|wrangle$ as fixed parameters is challenging, because most algorithms will rather find something local and then halt.

My implementation

The flat state $|srangle$ and the target state $|wrangle$ can be implemented in Python as

def flat(n):
    s = []
    for j in range(n):
        s.append(1/m.sqrt(n))
    return np.array(s)

def target(n,p):
    w=[]
    for j in range(n):
        if j==p:
            w.append(1)
        else: w.append(0)
    return np.array(w)

The probability function is then implemented as

def probability(x,G,n,p):
    t,gamma=x 
    w=target(n,p)
    U=expm(-1j*t*(gamma*nx.laplacian_matrix(G).todense()-np.outer(w,w)))
    psi=np.dot(U,flat(n))
    
    return -abs(np.dot(w,psi))*abs(np.dot(w,psi))

where I have used an intermediate state $|psi(t)rangle=U_H(t)|srangle$.

The Scipy library has a function minimize in the scipy.optimize module (note: the reason why the probability function returns a negative quantity is that there is no direct maximizing function, so we find the maxima of $f$ by finding the minima of $-f$). The documentation can be found here. I have used it to implement the minimizer function:

def minimizer(G,n,p,x0):
    res=minimize(probability,x0,args=(G,n,p))
    return [-res.fun, res.x[0]]

returning an array with $-P_{text{max}}$ and the corresponding $(t_{text{max}}, gamma_{text{max}})$. The argument x0 is an array containing the initial guess on $t$ and $gamma$ respectively.

The issue

scipy.optimize.minimize works only after setting an initial guess x0, and the result of the optimization strongly depends on such a choice. As a consequence, obtaining a graph of the parameters that maximize probability as a function of the number of vertices $n$ becomes impossible, as any x0 is only good in a range $(n_1, n_2)$ with endpoints depending on the chosen graph.

So here I ask for help: is there any way to make sure minimizer finds the global minima of probability (and in particular the first of them, if there are many) consistently in a dimension range $n in (0, Nge 100)$, regardless of the underlying graph?

I hope everything is as clear as possible, and thanks in advance to whoever might answer.

An example: the complete graph

For $n=20$, the probability function for the complete graph looks like this:

The complete graph is a good benchmark for testing the code as we know that $P_{text{max}}=1$ for $t_{text{max}}=pisqrt{n}/2$ and $gamma=1/n$ regardless of the chosen vertex $|wrangle$.

I have written a snippet of code that plots either of these quantities as functions of $n$, taking $|wrangle=|0rangle$:

#i=0 plots probability_max,  #i=1 plots t_max, #i=2 plots gamma_max
#plots in the interval [n=2, n=maxi) 
def plotComplete(x0,maxi,p,i):
    zPoints = []
    yPoints = []
    refX=[]
    refY=[]

    for z in range(2,maxi):
        y=round(minimizer(nx.complete_graph(z),z,p,x0)[i],5)
        zPoints.append(z)
        yPoints.append(y)

    pyplot.plot(zPoints, yPoints, 'b')
    
    if i==1:
        for z in range(2,maxi): 
            y=m.pi*m.sqrt(z)/2
            refX.append(z)
            refY.append(y)

        pyplot.plot(refX, refY, 'r')
    
    pyplot.show() 

Calling plotComplete([0.1,1],100,0,0) (notice that the initial guess on $t$ is $0.1$; changing it to, for example, $0$ or $1$ gives wildly different results) gives as expected the flat line

for the first $100$ complete graphs. Calling plotComplete([0.1,1],100,0,1) plots $t_{text{max}}$ with the same initial guess:

Here, I have made two overlapping plots: the red one is the theoretical value $pisqrt{n}/2$, the blue one is the implemented minimizer. For the most part up to around $n=60$ the two plots are in almost perfect agreement, then after that the initial guess $0.1$ is no more good enough to produce the correct result. The blue spikes appear any time the wrong minimum is found.