Claim 1: If $A$ is positive semidefinite then $\langle Ax, x \rangle \langle Ay,y\rangle \geq \langle Ax,y\rangle^{2}$.
This Cauchy--Schwartz inequality can be proved by considering quadratic function $\lambda \mapsto \langle x+\lambda y, A(x+\lambda y) \rangle \geq 0$ and writing down that its discriminant $D = 4 (\langle x, y\rangle ^{2} - \langle x, Ax\rangle \langle y, Ay\rangle$) must be non-positive.
Claim 2: if $A$ is positive definite then $\langle A^{-1}c,c \rangle =\max_{x}\left\{ 2\langle x,c \rangle - \langle Ax,x\rangle \right\}$ for all vectors $c$.
Indeed, the direction $\langle A^{-1}c,c \rangle \geq 2\langle x,c \rangle - \langle Ax,x\rangle$ can be obtained as follows:
$\langle A^{-1}c,c \rangle +\langle Ax,x\rangle \geq 2 \sqrt{\langle A^{-1}c,c \rangle \langle Ax,x\rangle} \geq 2|\langle x,c \rangle|$, where the last inequality with $c=Ay$ coincides with Claim 1. Since $A$ is invertible it covers all c's. On the other hans the choice $x=A^{-1}c$ in $\max_{x}\left\{ 2\langle x,c \rangle - \langle Ax,x\rangle \right\}$ gives $\langle A^{-1}c,c \rangle$.
Claim 3: $(\frac{A+B}{2})^{-1} \leq \frac{A^{-1}+B^{-1}}{2}$ holds for all positive definite matrices.
We have $\langle c, (\frac{A+B}{2})^{-1}c \rangle = \max_{x} \left\{ 2\langle x,c \rangle - \langle (\frac{A+B}{2})x,x\rangle \right\} \leq \frac{\max_{x} \left\{ 2\langle x,c \rangle - \langle (A)x,x\rangle \right\}}{2} +\frac{\max_{x} \left\{ 2\langle x,c \rangle - \langle Bx,x\rangle \right\}}{2} = \frac{\langle A^{-1}c,c\rangle + \langle B^{-1}c,c\rangle}{2}$ $\square$.
- the idea of representing convex function as a sup of linear functions (or sometimes even sup of concave functions with comparable quantity) is a standard technique in stochastic optimal control theory.
- the link that I provided in the comments, I think the second answer using derivatives looks correct to me.