If you ignore the relationships between the $l_i$, you can still use the union bound to get an upper bound on the probability that the maximum is at least $x$ that is off by at most a factor of $n$, and then you can use this to show the expected value of the maximum is $H_n(1+o(1))$.

Since each $l_i$ is a sum of Bernoulli random variables, the indicators of positions of potential minima, we can use [Hoeffding's inequality][1] for sums of bounded random variables. The average value of $l_i$ is $H_{n-i+1}$. It is a sum of $n-i+1$ Bernoulli random variables, so for $x \ge 0$, $P[l_i - H_{n-i+1} \gt t] \le \exp(-2(n-i+1)t^2)$, or for $i'=n-i+1$, $P[l_{i'} - h_i \gt t] \le \exp(-2 i t^2)$. Choose any $\epsilon \gt 0$.

$$\begin{eqnarray} E[\max(l_1,...,l_n)] & \le & (1+\epsilon)H_n + \sum_{x \gt (1+ \epsilon)H_n} P(\max(l_1,...,l_n) \gt x) \newline & \le & (1+\epsilon)H_n + \sum_{x \gt (1+\epsilon)H_n} \sum _{i=1}^n P(l_i \gt x) \newline & \le & (1+\epsilon)H_n + \sum_{x \gt (1+\epsilon)H_n} \sum _{i=1}^n \exp(-2 (n-i+1)(x-H_{n-i+1})^2) \newline & = & (1+\epsilon)H_n + \sum_{x \gt (1+\epsilon)H_n} \sum _{i=1}^n \exp(-2 i(x-H_{i})^2) \end{eqnarray}$$

We can change the order of summation and estimate very coarsely

$$\begin{eqnarray}\sum_{x \gt (1+\epsilon)H_n} \exp (-2i(x-H_i)^2) & \le & \sum_{x \ge 0} \exp(-2i(x+(1+\epsilon)H_n-H_i)^2) \newline & \le & \sum_{x \ge 0} \exp(-2(x+\epsilon H_n)^2) \newline & \le & \exp(-2\epsilon^2 H_n^2) \frac{1}{1-\exp(-2(2\epsilon H_n+1))}\end{eqnarray}$$

by comparison with a geometric series with the same first two terms. Since $H_n \gt \log n$, 

$$ \begin{eqnarray} & & &\le &\exp(-2 \epsilon^2 (\log n)^2) \frac{1}{1- \exp (-4 \epsilon H_n - 2)} \newline & & &\le & \exp(-2 \epsilon^2 (\log n)^2) \frac{1}{1- \exp (- 2)} \newline & & &\le& 2 \exp(-2 \epsilon^2 (\log n)^2) \newline & & & \le & 2 n^{-2 \epsilon^2 \log n} \end{eqnarray}$$

This estimate doesn't depend on $i$. 

$$\begin{eqnarray} E[\max(l_1,...,l_n)] & \le & (1+\epsilon) H_n + \sum_{i=1}^n 2 n^{-2 \epsilon^2 \log n} \newline &=& (1+\epsilon)H_n + 2n^{-2 \epsilon^2\log n + 1} \newline & = & (1+\epsilon)H_n + o(1).\end{eqnarray} $$

So, $E[\max(l_1,...l_n)] \sim H_n$. The estimates here were coarse, and while they could be improved to estimate the difference $ E[\max(l_1,...l_n)] - H_n$, using some of the dependencies between the $\lbrace l_i \rbrace$ should be necessary to get the true order.


  [1]: http://en.wikipedia.org/wiki/Hoeffding%27s_inequality