Optimal Talmudic Zigzag
This is a special case of finding the longest path in a directed acyclic graph. Namely, the vertices of our graph are $1$, $2$, ..., $n$, there is an edge $i \to j$ for each $i<j$ and the length of that edge is $\log \tfrac{p_i + p_j}{2 \sqrt{p_i p_j}}$. We want the longest path from $1$ to $n$.
There are a number of standard efficient algorithms for this; see for example Section 4.7 in Dasgupta, Papadimitriou and Vazirani's Algorithms.
An observation: If we chose to include $z_i$, then we must have $$\frac{p_{z_{i-1}} + p_{z_i}}{2 \sqrt{p_{z_{i-1}} p_{z_i}}}\ \frac{p_{z_{i}} + p_{z_{i+1}}}{2 \sqrt{p_{z_{i}} p_{z_{i+1}}}} \geq \frac{p_{z_{i-1}} + p_{z_{i+1}}}{2 \sqrt{p_{z_{i-1}} p_{z_{i+1}}}}$$ and multiplying this out and rearranging gives $$\left( \frac{p_{z_{i+1}}}{p_{z_i}} - \frac{p_{z_{i}}}{p_{z_{i+1}}}\right) \left( \frac{p_{z_{i}}}{p_{z_{i-1}}} - \frac{p_{z_{i-1}}}{p_{z_{i}}}\right) < 0$$ so either $p_{z_{i-1}} < p_{z_i} > p_{z_{i+1}}$ or $p_{z_{i-1}} > p_{z_i} < p_{z_{i+1}}$. In other words, the sequence $p_{z_1}$, $p_{z_2}$, ..., $p_{z_k}$ is alternating. Finding the longest alternating subsequence of a sequence is a well studied problem, and good algorithms for that might be good heurisitics for this.
Another thought (and I'm going to drop out of this conversation soon): Running times for longest path algorithms depends on the number of edges in the graph. If the optimal solution uses the edge $j \to \ell$, then we must have $p_k$ between $p_j$ and $p_{\ell}$ for all $j < k < \ell$. (Otherwise, inserting $p_k$ would be an improvement.) Let's call $(j, \ell)$ a good pair if this condition holds.
Let's insert a preprocessing step which finds all good pairs $(j, \ell)$ and restrict the graph algorithm to them. Here is how that preprocessing step works. For each $j$, do the following. Suppose $p_{j+1} > p_j$ (if $p_{j+1} < p_j$, reverse all inequalities that follow.) Let $k_j$ be the smallest index $k_j>j$ for which $p_{k_j}<p_j$. If the $p_i$ are independently identically distributed, then $k_j-j=r$ with probability $1/2^{r-1}$ for $r \geq 2$. The good pairs are $(j,\ell)$ with $j < \ell < k_j$ where $p_{\ell}$ is the maximum of $(p_{j+1}, p_{j+2}, \dots, p_{\ell})$. So the expected search time to find $k_j$ is $\sum_{r \geq 2} r/2^{r-1} =3$ and the expected number of such $\ell$ is $\sum_{r \geq 2} 2^{-(r-1)} (1+1/2+\cdots+1/(r-1)) = 2 \log 2 \approx 1.39$.
Thus, if your $p_i$ are independently identically distributed, then in $O(n)$ time you can find the good pairs, and there are only about $1.39 n$ of them. Djisktra's algorithm on a graph with $O(n)$ edges, implemented with a binary heap, runs in time $O(n \log n)$.