Hadamard's three circle theorem
Let $\lambda=\frac{\log(b/r)}{\log(b/a)}$. Then $1-\lambda=\frac{\log(r/a)}{\log(b/a)}$. Dividing both sides of your equation by $\log(b/a)$ gives:
$$\log(M(r))\leq \lambda \log(M(a))+(1-\lambda)\log(M(b)).$$
Notice that $a^{\lambda}b^{1-\lambda}=\exp(\lambda\log(a)+(1-\lambda)\log(b))=r$ (verify this!). Thus,
$$\log(M(a^\lambda b^{1-\lambda})\leq \lambda\log(M(a))+(1-\lambda)\log(M(b))$$
which is saying that $M$ is log-convex, in that $\log(M(\exp(z))$ is convex in $z$. It is not hard to show that the only times you get strict equality for convex functions is iff the function is affine (linear): $\log(M(\exp(z))=Az+B$. Thus, $M(\exp(z))=Ce^{Az}$, or $M(z)=Cz^A$
Addendum: to show that all strictly equal convex functions are affine, write:
$f''(z)=\lim_{h\rightarrow 0}\frac{f(z+h)+f(z-h)-2f(z)}{h^2}=\lim_{h\rightarrow 0}\frac{2f(z)-2f(z)}{h^2}=0$
where we used $z=\frac{1}{2}(z+h)+\frac{1}{2}(z-h)$ and equal convexity.
If $f=cz^\lambda$ for some $c\in\mathbb C$ and $\lambda\in\mathbb Z$, then $f$ is holomorphic on $\mathbb C-\{0\}$, and $M(r)=|c|r^\lambda$. Thus $\log M(r)=\log c+\lambda \log r$, and
$$\log\frac br\log M(a)+\log\frac ra\log M(b)=\log\frac br\left(\log c+\lambda \log a\right)+\log\frac ra\left(\log c+\lambda \log b\right)$$ $$=(\log b\log c-\log r\log c+\lambda\log b\log a-\lambda\log r\log a)$$ $$\qquad+\ (\log r\log c-\log a\log c+\lambda\log r\log b-\lambda\log a\log b)$$ $$=\log b\log c-\lambda\log r\log a-\log a\log c+\lambda\log r\log b$$ $$=\log\frac ba\left(\log c+\lambda \log r\right)=\log\frac ba\log M(r).$$
Thus the inequality is saturated. Conversely, if $\log\frac br\log M(a)+\log\frac ra\log M(b)=\log\frac ba\log M(r)$ for all $a<r<b$, then
$$\log r\log\frac{M(b)}{M(a)}+(\log b\log M(a)-\log a\log M(b))=\log\frac ba\log M(r)$$
so if we let $\lambda\log\frac ba=\log\frac{M(b)}{M(a)}$ and $c'\log\frac ba=\log b\log M(a)-\log a\log M(b)$, then
$$\lambda\log r+c'=\log M(r)\Rightarrow M(r)=cr^\lambda$$
where $c=e^{c'}$. Now since $M(r)=\max_{|z|=r}|f(z)|$, for each point $r_0$ there is a corresponding point $z_0$ such that $|f(z_0)|=M(r_0)$. Since $f$ is continuous on the circle of radius $r_0$, $\frac d{d\theta}|f(z)|=0$ at $z_0$. Also, $\frac d{dr}|f(z)|=M'(r)$. Thus, $\frac d{dz}f(z)=\psi z/|z|M'(|z|)=c\psi\lambda z^{\lambda-1}$ at $z_0$ (where $|\psi|=1$) by the Cauchy-Riemann relations (which is to say, the "full" derivative of $f$ at $z_0$ is consistent with the derivative in the radial direction being $M'(|z_0|)$.
The $|f(z_0)|\leq M(|z_0|)$ constraint is strong enough to constrain the higher derivatives in the same manner, so that in fact $f(z)$ looks like $g(z)=\alpha z^\lambda$ (where $\alpha=c\psi$) to all orders. To show this, note that we have here a tight inequality on "one side" of $f$ (the outer radial direction) against another analytic function $g$ in some neighborhood of $z_0$ - we are going to show that any perturbation of $f$ will push up against this wall somewhere.
If the derivatives of $f$ and $g$ do not coincide, we can zoom in close enough that the lowest different derivative is the dominant factor in the difference $f-g$. Then $f(z+z_0)\approx g(z+z_0)+dz^n+O(z^{n+1})$ for some $n\in\Bbb N$ and $d\in\Bbb C$. We choose a direction $\phi$ such that $d\phi^n$ has the same phase as $g(z_0)$, and then taking $z=\epsilon\phi$ for sufficiently small $\epsilon$ yields $|f(z+z_0)|\approx|g(z+z_0)|+|d||z|^n+O(|z|^{n+1})$, which violates the bound $|f(z+z_0)|\le |g(z+z_0)|$.
Thus, $f(z)=\alpha z^\lambda$, but in analytically continuing around the pole at zero, one may find that $f$ does not match up, i.e. a branch cut is necessary. This contradicts the assumption that $f$ is holomorphic on $A$, so this constrains the set of possible $f$'s satisfying the criteria from $\{cz^\lambda:c,\lambda\in\mathbb C\}$ to $\{cz^\lambda:c\in\mathbb C,\lambda\in\mathbb Z\}$.