Why the sequence of Bernstein polynomials of $\sqrt x$ is increasing? Bernstein polynomials preserves nicely several global properties of the function to be approximated: if e.g. $f:[0,1]\to\mathbb R$ is non-negative, or  monotone, or convex; or if it has, say, non-negative $17$-th derivative, on $[0,1]$, it is easy to see that the same holds true for the polynomials $B_nf$. In particular, since  all $B_n$ fix all affine functions, if $f\le ax+b$, also $B_nf(x)\le ax+b$, whence it follows immediately $B_nf\le f$ for concave $f$. On the other hand, comparing $B_nf$ and $B_{n+1}f$ turns out to be harder, due to the different choice of nodes where $f$ is evaluated. Consider for instance the Bernstein polynomials of the function $\sqrt{x}$,
$$p_n(x):=\sum_{k=0}^n{n\choose k}\sqrt{\frac kn}\,x^k(1-x)^{n-k}.$$

Question: Is this sequence of polynomial increasing? More generally, when is $B_nf$ increasing?

Some tentative approaches and remarks.
1. To compare $p_{n+1}$ with $p_n$ we may write the binomial coefficients in the expression for $p_{n+1}(x)$ as ${n+1\choose k}={n\choose k}+{n\choose k-1}$; splitting correspondingly the sum into two sums, and shifting the index in the latter,
we finally get
$$p_{n+1}(x)-p_n(x)=\sum_{k=0}^n{n\choose k}\bigg[x\sqrt{\frac{k+1}{n+1}}+(1-x)\sqrt{\frac k{n+1}}-\sqrt{\frac kn}\,\bigg]x^k(1-x)^{n-k},$$
which at least has non-negative terms approximatively for $\frac kn<x$, which is still not decisive.
2. Monotonicity of the sequence $B_nf$ is somehow reminiscent of that of the real exponential  sequence $\big(1+\frac xn\big)^n$. Precisely, let $\delta_n:f\mapsto \frac{f(\cdot+\frac1n)-f(\cdot)}{\frac1n}$ denote the discrete difference operator, and $e_0:f\mapsto f(0)$ the evaluation at $0$. Then the Bernstein operator $f\mapsto B_nf$ can be written as $B_n=e_0\displaystyle \Big({\bf 1} + \frac{x\delta_n}n\Big)^n$ (which, at least for analytic functions, converges to the Taylor series $e^{xD}$ at $0$). Unfortunately, the analogy seems to stop here.
3. The picture below shows the graphs of  $\sqrt x$ and of the first ten $p_n(x)$.
(The convergence is somehow slow; indeed it is $O(n^{-1/4})$, as per Kac' general estimate $O(n^{-\frac \alpha2})$ for $\alpha$-Hölder functions). The picture leaves some doubts about the endpoints; yet there should be no surprise, since $p_n(0)=0$, $p_n'(0)=\sqrt{n}\uparrow+\infty$, $p_n(1)=1$, $p_n'(1)=\frac1{1+\sqrt{1-\frac1n}}\downarrow\frac12$.

 A: As noted by Paata Ivanishvili,  if $f$ is concave on $[0,1]$, then the Bernstein polynomials $B_n(f,p)$ are increasing in $n$.
Here is a probabilistic proof:
Let $I_j$ for $j \ge 1$ be independent variables taking value 1 with probability $p$ and $0$ with probability $1-p$. Then  $X_n:=\sum_{j=1}^n I_j$ has a binomial Bin$(n,p)$ distribution and the Bernstein polynomial can be written as $B_n(f,p)=E[f(X_n/n)]$. Now for every $j \in [1,n+1]$, the random variable $Y_j=Y_j(n)=X_{n+1}-I_j$ also has a Bin$(n,p)$  distribution
and $$ {X_{n+1}} = {\sum_{j=1}^{n+1} (Y_j/n)} \, .$$
For concave $f$, Jensen's inequality gives
$$  f \left(\frac{\sum_{j=1}^{n+1} (Y_j/n)}{n+1} \right) \ge  \left(\frac{\sum_{j=1}^{n+1} f(Y_j/n)}{n+1} \right)  $$
whence
$$B_{n+1}(f,p)=E f \left(\frac{X_{n+1}}{n+1}\right)=E f \left(\frac{\sum_{j=1}^{n+1} (Y_j/n)}{n+1} \right) \ge  E \left(\frac{\sum_{j=1}^{n+1} f(Y_j/n)}{n+1} \right) =B_n(f,p) $$
A: We have
$$B_n(f,\frac{u}{u+v}) = \frac{1}{(u+v)^n} \sum_{k=0}^n \binom{n}{k} f(\frac{k}{n})\,u^k v^{n-k}$$
so
$$B_n(f,\frac{u}{u+v})- B_{n+1}(f,\frac{u}{u+v})=\\=
\frac{1}{(u+v)^{n+1}}\left((u+v)\sum_{k=0}^n \binom{n}{k} f(\frac{k}{n})\,u^k v^{n-k}- \sum_{k=0}^{n+1} \binom{n+1}{k} f(\frac{k}{n+1})\,u^k v^{n+1-k}\right)$$
Now it's enough to see that the coefficient of each $u^k v^{n+1-k}$ is positive.  Indeed, it equals
$$  \binom{n}{k-1} f( \frac{k-1}{n}) + \binom{n}{k} f(\frac{k}{n}) - \binom{n+1}{k} f(\frac{k}{n+1})$$
Now note that $\binom{n}{k-1} + \binom{n}{k} = \binom{n+1}{k}$ and
$$\frac{k-1}{n} \binom{n}{k-1} + \frac{k}{n} \binom{n}{k} = \binom{n-1}{k-2} + \binom{n-1}{k-1} = \binom{n}{k-1} = \frac{k}{n+1}\binom{n+1}{k}$$
so everything checks out: if $f$ is convex then the sequence of functions  $(B_n(f,\cdot))$ on $[0,1]$ is decreasing in $n$. We also note that each $B_n(f, x)$ is convex in $x$. Indeed, we can show the formula
$$B^{(2)}_{n+2}(f,\frac{u}{u+v}) = \frac{(n+2)(n+1)}{(u+v)^n}\cdot \sum_{k=0}^{n}  \Delta^2 c(k)\binom{n}{k} u^k v^{n-k}$$
where $c(k) = f(\frac{k}{n+2})$ and $\Delta$ is the forward discrete derivative.
By the way, this formula generalizes for any derivative, with the consequence that if $f^{(l)}\ge 0$ ( even in the weaker sense) on $[0,1]$, then the same holds for all the associated Bernstein polynomials.
A: While nothing will beat the brilliant probabilistic proof given in Yuval Peres's answer, a more conventional argument goes as follows. Write $$a_{n,k} = \tbinom nk x^k (1-x)^{n-k} $$ and $$ p_{n,k} = \tfrac{k}{n} $$ Observe that
$$ \tfrac{k}{n} = (1 - p_{n,k}) \tfrac{k}{n-1} + p_{n,k} \tfrac{k-1}{n-1} . $$
Thus, if $f$ is concave, then
$$ f(\tfrac{k}{n}) \geqslant (1 - p_{n,k}) f(\tfrac{k}{n-1}) + p_{n,k} f(\tfrac{k-1}{n-1}) . $$
Multiply this by $a_{n,k}$ and add up to get
$$ \begin{aligned}
p_n(x) & = \sum_{k=0}^n a_{n,k} f(\tfrac{k}{n}) \\
& \geqslant \sum_{k=0}^n a_{n,k} \Bigl( (1 - p_{n,k}) f(\tfrac{k}{n-1}) + p_{n,k} f(\tfrac{k-1}{n-1}) \Bigr) \\
& = \sum_{k=0}^{n-1} a_{n,k} (1 - p_{n,k}) f(\tfrac{k}{n-1}) + \sum_{k=1}^n a_{n,k} p_{n,k} f(\tfrac{k-1}{n-1}) \\
& = \sum_{k=0}^{n-1} \bigl(a_{n,k} (1 - p_{n,k}) + a_{n,k+1} p_{n,k+1} \bigr) f(\tfrac{k}{n-1}) \\
& = \sum_{k=0}^{n-1} a_{n-1,k}   f(\tfrac{k}{n-1})=p_{n-1}(x).
\end{aligned} $$

(I make this answer CW, since I expect this is exactly the solution Paata had in mind.)
A: I would also be happy to know the converse implication, that you quoted in the last remark of your notes
Let $f \in C([0,1])$. Then $p_{n+1}(f,x) \geq p_{n}(f,x)$ for all $n \in \mathbb{N}$ and all $x \in [0,1]$ if and only if $f$ is concave. Yuval gave a nice proof of one implication. To show the converse, assume contrary that $f$ is not concave. By adding a linear function if necessary, and scaling the result by a positive constant we can assume that there exists $x \in [a',b'] \subseteq [0,1]$ such that $f(a')=f(b')>1$, and
$$0=\min_{[a',b']}f = f(x) < \frac{x-a'}{b'-a'}f(b')+\frac{b'-x}{b'-a'}f(a').$$
By continuity we can further find $a<b$ such that $0\leq a' <a <x<b <b' \leq 1$ so that $\min\limits _{[a'b']\setminus [a,b]} f \geq 1$. Then, I claim that for $n$ large enough we must have $p_{n}(f,x)>0$. Indeed,
\begin{eqnarray*}
p_{n}(f,x) & = & \sum_{0\leq k \leq n} f\left(\frac{k}{n}\right)\binom{n}{k}x^{k}(1-x)^{n-k} \\
& \geq & -\|f\|_{C}\,  n \max\limits_{0\leq k \leq [a'n]}\binom{n}{k}x^{k}(1-x)^{n-k}+\binom{n}{[an]}x^{[an]}(1-x)^{n-[an]} \\
& & + \binom{n}{[bn]}x^{[bn]}(1-x)^{n-[bn]} -\|f\|_{C}\,  n\max\limits_{[b'n]\leq k \leq n}\binom{n}{k}x^{k}(1-x)^{n-k}.
\end{eqnarray*}
Let us show that the sum of the  first two terms in the right hand side of the inequality is positive (for $n$ large enough), similarly the sum of the last two terms will be positive. Indeed, to verify
$$
\binom{n}{[an]}x^{[an]}(1-x)^{n-[an]} >\|f\|_{C}\,  n \max\limits_{0\leq k \leq [a'n]}\binom{n}{k}x^{k}(1-x)^{n-k}
$$
Take $1/n$ power and let $n \to \infty$. Since $\binom{n}{sn}^{1/n} \to s^{-s}(1-s)^{s-1}$ as $n \to \infty$ it suffices to verify that
$$
\max_{0\leq s \leq a'}s^{-s}(1-s)^{s-1}x^{s} (1-x)^{1-s} < a^{-a}(1-a)^{a-1}x^{a}(1-x)^{1-a}
$$
Notice that $\varphi(s) := s^{-s}(1-s)^{s-1}x^{s} (1-x)^{1-s}$ is increasing on $[0,a]$. Indeed,
$$
(\ln \varphi (s) )' =  \ln \left(\frac{1-s}{s} \cdot \frac{x}{1-x}\right) >0,
$$
where the last inequality follows because $x>a\geq s$.
Thus for $n$ large enough $p_{n}(f,x)>0$. On the other hand, as Yuval pointed out,  $p_{n}(f,x) = \mathbb{E}f(\frac{\xi_{1}+\cdots+\xi_{n}}{n}) \to f(\mathbb{E} \xi_{1})=f(x)=0$ by the strong law of large numbers. Since $n \mapsto p_{n}(f,x)$ increasing we have $p_{n}(f,x) \leq 0$ which is in contradiction with $p_{n}(f,x)>0$ for $n$ large enough.
