Looking for different proofs of "Discrete Liouville's Theorem".

Good day.

There is a question I have already encountered twice, in very different contexts, that is relatively simple looking, but both solutions I know involve some pretty advanced theorems from the respective fields. I would like to ask if someone could think of any different proofs.

So, here is the question: Let $f: \mathbb Z^2 \rightarrow \mathbb R$ satisfy $f(m,n) = \frac14(f(m+1,n) + f(m-1,n) + f(m,n+1) + f(m,n-1))$. i.e. $f$ is defined on a grid, at a value at a point is the average of its four neighbors. Assume that $f$ is bounded (in some versions, bounded from below only). Prove that $f$ is constant.

One proof of the question uses martingales, and relies, I think on the martingale convergence theorem. In the proof, you define a two-dimensional random walk $S_n$ and look at the process $f(S_n)$.

The other proof comes from functional analysis, and uses the Krein-Milman theorem on the set $K_1 = \{f \in L_{\infty}(\mathbb Z^2) | f(m,n) = \frac14(...), ||f||\leq1 \} $, after finding its extremal set (which is the constant functions $\{+1,-1\}$).

So, any other proofs of this seemingly simple question? I should suspect there is something related to complex analysis, since for analytic functions this would follow from the Liouville theorem. Also, perhaps a combinatorial solution?

This is the discrete equivalent of Liouville's principle for harmonic function: every bounded function with a null laplacian is constant.

Without loss of generality you can assume that $f(0,0)=0$ (if not so, replace $f$ with $f-f(0,0)$). Now consider the square $B_N=\{(m,n):|m|+|n|\leq N\}$. Over $B_N$, $f$ attains its maximum and minimum on $\partial B_N$, and $$\Gamma_N = \max_{(m,n)\in B_N}f(m,n)=\max_{(m,n)\in \partial B_N}f(m,n)$$ is non decreasing as a function of $N$. $\Gamma_N$ is also bounded by $\Gamma$, and satisfies the almost-midpoint-convexity condition: $$\Gamma_N\leq\frac{3\Gamma_{N+1}+\Gamma_{N-1}}{4}$$ since the values of $f$ over $\partial B_N$ are completely determined by the values of $f$ over $\partial B_{N+1}$ and $\partial B_{N-1}$. Since any non-decreasing, bounded, midpoint-convex function is constant, we would expect that the same applies here. We have: $$\Gamma_{N+1}-\Gamma_{N}\geq\frac{1}{3}(\Gamma_{N}-\Gamma_{N-1})\geq\frac{\Gamma_1}{3^N}$$ so for any $\epsilon>0$ we can find a point $(a_1,b_1)$ such that $f(a_1,b_1)\geq\frac{3}{2}\Gamma_1-\varepsilon>\frac{10}{7}\Gamma_1.$ This gives that the difference of the values of $f$ in two neighbour points cannot be bounded, hence the function itself cannot be bounded.

We have just used a Borel-Caratheodory-type argument to provide lower bounds for $\Gamma_N$.