Mean field approximation convergence 
*

*The last sentence of Christopher M. Bishop, Pattern Recognition and Machine Learning Section 10.1.1 Factorized distributions on p.466, states, referring to Equation $(10.9)$, that "Convergence is guaranteed because bound is convex with respect to each of the factors $q_i(Z_i)$ (Stephen Boyd and Lleven Vandenberghe (2004), Convex Optimization. Cambridge University Press)."

*The same claim is made in the Wikipedia entry, Variational Bayesian Methods Section 2. Mean field Approximation with the same reference to the Boyd & Vandenberghe book above: "An algorithm of this sort is guaranteed to converge."

I see that $q^*_j(Z_j)$ minimizes $-\mathcal L(q)$ and thus $-\mathcal L(q^*)$ is nonincreasing with respect to the iteration on $q^*$. However, I do not see how the convexity of $-\mathcal L$ with respect to $q$ leads to the convergence of $q^*$. On the contrary, I can imagine counterexamples satisfying the conditions of decreasing functional value and convexity of the function. Neither do I know where in the referenced Boyd book to find the proof. Is there a fixed-point theorem to use?

I haven’t read too much into section 10.1.1 in Bishop’s book, but the full paragraph you are referring to is (emphasis in bold mine)

The set of equations given by (10.9) for $j = 1,...,M$ represent a set of consistency conditions for the maximum of the lower bound subject to the factorization constraint. However, they do not represent an explicit solution because the expression on the right-hand side of (10.9) for the optimum $q^*_j(Z_j)$ depends on expectations computed with respect to the other factors $q_i(Z_i)$ for $i \neq j$. We will therefore seek a consistent solution by first initializing all of the factors $q_i(Z_i)$ appropriately and then cycling through the factors and replacing each in turn with a revised estimate given by the right-hand side of (10.9) evaluated using the current estimates for all of the other factors. Convergence is guaranteed because bound is convex with respect to each of the factors $q_i(Z_i)$ (Boyd and Vandenberghe, 2004).

The emphasis in bold suggests that the author means that the negative KL divergence $\mathcal L(q)$ is maximized using coordinate ascent, where one component of $q$ is adjusted while all others are fixed. The proof of convergence of coordinate ascent in the convexity case is not as straightforward as the author makes it out to be. See for example the references here.
