Let $V$ and $H$ denote the sets of visible and hidden indices. We use $s_V := \{ s_i\}_{i \in V}$ and $s_H := \{s_i \}_{i \in H}$ to denote the entire collections of visible and hidden variables.
Suppose $s_V^\star := \{s_i^\star\}_{i\in V}$ are the values taken by the visible variables $s_V$ in our data. Our objective is to tune the $w_{ij}$ and $b_i$ parameters so as to optimise the log-likelihood of observing this data:
\begin{align*} \log P(s^\star_{V}) &= \log \left(\sum_{s_{H}} P(s^\star_{V}, s_{H}) \right) \\ &= \log \left( \frac{\sum_{s_{H}} e^{-E(s^\star_{ V}, s_{H})} }{\sum_{s} e^{-E(s)} }\right).\end{align*}
We solve this optimisation problem by gradient descent.
Taking derivatives respect to $w_{ij}$, we have
\begin{align*} \frac{\partial \log P(s^\star_V)}{\partial w_{ij}} &= - \frac{\sum_{s_{H}} \frac{\partial E(s^\star_{V}, s_{ H})}{\partial w_{ij}}e^{-E( s^\star_{V}, s_{ H})} }{\sum_{ s_{H}} e^{-E(s^\star_{ V}, s_{ H})} } + \frac{\sum_{s} \frac{\partial E(s)}{\partial w_{ij}}e^{-E(s)} }{\sum_{s} e^{-E(s)} } \\ &= - \sum_{s_{ H}} \left(\frac{\partial E(s^\star_{ V}, s_{ H})}{\partial w_{ij}} P\left( s_{ H} \mid s^\star_{V}\right)\right)+ \sum_{s} \left(\frac{\partial E(s)}{\partial w_{ij}} P\left( s\right)\right) \\
&= \mathbb E_{s_{H} \sim P\left(s_{H} \mid s^\star_{V}\right)} \left[s_i s_j|_{s_V = s^\star_V}\right]- \mathbb E_{s \sim P(s)}\left[s_i s_j \right]
\end{align*}
In your paper's notation,
$$ \langle s_i s_j \rangle_{\rm data} := \mathbb E_{s_{H} \sim P\left(s_{H} \mid s^\star_{V}\right)} \left[s_i s_j|_{s_V = s^\star_V}\right],$$
$$\langle s_i s_j \rangle_{\rm model } := \mathbb E_{s \sim P(s)}\left[s_i s_j \right],
$$
which completes the derivation.
[Note: $$P\left( s_H \mid s_V^\star \right) = \frac{e^{-E(s_V^\star, s_H)}}{\sum_{s_H}e^{-E(s_V^\star, s_H)}}$$ is the conditional distribution for the hidden variables $s_H$, given the values $s_V^\star$ for the visible variables observed in the data. To sample from this distribution, you clamp the values of $s_i$ for $i\in V$ at the values $s_i^\star$ observed in the data, and only perform the Gibbs sampling for the variables $s_i$ with $i\in H$.]
It depends on what you mean by "optimize" or "solve".
Essentially, there are (at least) two distinct notions of a computer accessing a function --- white box, and black box.
Black box access is somewhat more common (within theoretical computer science at least). You view the computer as having some "oracle" $\mathcal{O}(\cdot)$ which it can query. It has no idea of the inner workings of the oracle, but it can send the oracle inputs, and receive outputs back from the oracle. You can view this as having some API for some function, but having absolutely no clue how it is implemented.
White box access means that you give the computer "the code of the function" in some way. Maybe (theoretically) you give it $M$, the description of a Turing Machine (or some other model of computation) which computes the function. Maybe (practically) you give the computer access to the file in which the function is defined.
We'll discretize things to make things simpler (from a computation perspective). Let $f : D \to R$ be the function you want to optimize. Then, for an arbitrary discontinuous function $f$, black-box optimization takes at least $|D|$ queries to the oracle, i.e. you must check all points.
The proof is as follows. Say you have some optimization function which queries the oracle for the function $f$ at points $D' \subsetneq D$. Let $x\in D\setminus D'$ be a point which is not queried. Then the optimization function is incorrect on the function $f'$, which is the same as $f$ on $D'$, but has $f(x) = 1 + \max_{x\in D} f(x)$.
Note that this (simplistic) proof stops working if your optimization function is randomized, but I believe that your general intuition still holds true that one cannot optimize arbitrary discontinuous functions.
In the white box model there are things you can do, but it is highly dependent on how the function is presented. Your function (for example) can be represented via a tree, where each internal node you compute some predicate of the input then branch left or right as a result (this is formalized via the notion of a decision tree). While I do not know the state of art for optimizing functions represented as decision trees, there's the following simple algorithm that shows you get "more power" than the black box query model --- each leaf of the tree corresponds to an output. Sort them. Then loop through the (sorted) leaves, and check if there exists an input which satisfies all predicates on the internal nodes on the path to that leaf. The first leaf for which there does exist such an input is the maximum value.
As for why this is "more power", technically the amount of power depends on the complexity of the underlying decision tree. But for your function things reduce to analyzing a binary tree of depth two (which is quite simple). How long the black-box model takes depends on how exactly you discretize things, but if you wanted to "query all floating point numbers" (for example) it suffices to say it would take a while.
Note that for more general notions of computation (say Turing machines) little is possible for the "standard reasons" (the Turing machines may not even halt).
That being said, if you know that the function is "nice enough", you could imagine doing something similar to what I described with decision trees for the control flow graph of your function. Decision trees control flow graph is especially simple --- it's a tree (instead of a general graph), and there is no mutation of variables.
There may be thorny details when one steps away from this discretized world that I just discussed, but I (unfortunately) cannot comment on those.
Best Answer
In Novikoff's theorem, we assume that
i) The data is linearly separable: $$\forall(\vec{x}, y) \in \mathcal{X} \text{ } \exists \vec{w}_* \exists \gamma > 0: \langle\vec{w}_*, \vec{x}\rangle y \ge \gamma .$$ ii) The weights are updated following Hebb's rule: $$\text{if } \langle\vec{w}_{t-1},\vec{x}\rangle y < 0, \text{ then } \vec{w}_t \leftarrow \vec{w}_{t-1} + y\vec{x} .$$
Thus,
$$\langle\vec{w}_t , \vec{w}_*\rangle^2 = \langle\vec{w}_{t-1}+y\vec{x} , \vec{w}_*\rangle^2 = (\langle\vec{w}_{t-1}, \vec{w}_*\rangle + \langle\vec{w}_*, y\vec{x}\rangle)^2 = (\langle\vec{w}_{t-1}, \vec{w}_*\rangle + \langle\vec{w}_*, \vec{x}\rangle y)^2 \ge (\langle\vec{w}_{t-1} , \vec{w}_*\rangle+\gamma)^2 .$$
$$(\langle\vec{w}_{t-1}, \vec{w}_*\rangle + \langle\vec{w}_*, \vec{x}\rangle y)^2 = (\langle\vec{w}_{t-2}, \vec{w}_*\rangle + 2\langle\vec{w}_*, \vec{x}\rangle y)^2 = \ldots =$$
$$= (\langle\vec{w}_{0}, \vec{w}_*\rangle + t\langle\vec{w}_*, \vec{x}\rangle y)^2 = (\langle0, \vec{w}_*\rangle + t\langle\vec{w}_*, \vec{x}\rangle y)^2 \ge t^2\gamma^2.$$
$$\le ||\vec{w}_{t-1}||^2 + ||\vec{x}||^2 \le ||\vec{w}_{t-1}||^2 + R^2 \le \ldots \le ||\vec{w}_0||^2 + t^2R^2 = t^2R^2.$$