I don't think your question is a mathematical one, for the question about what do all inequalities eventually reduce to has a simple answer: axioms. I interpret it as a metamathematical question and still I believe the closest answer is the suggestion above about using everything you know.
An inequality is a fairly general mathematical term, which can be attributed to any comparison. One example is complexity hierarchies where you compare which of two problems has the highest complexity, can be solved faster etc. Another one is studying convergence of series, that is comparing a quantity and infinity, here you find Tauberian theory etc. Even though you did not specify in your question which kind of inequalities are you interested in primarily, I am assuming that you are talking about comparing two functions of several real/complex variables. I would be surprised if there is a list of exclusive methods that inequalities of this sort follow from. It is my impression that there is a plethora of theorems/principles/tricks available and the proof of an inequality is usually a combination of some of these. I will list a few things that come to my mind when I'm trying to prove an inequality, I hope it helps a bit.
First I try to see if the inequality will follow from an equality. That is to recognize the terms in your expression as part of some identity you are already familiar with. I disagree with you when you say this shouldn't be counted as a method to prove inequalities. Say you want to prove that $A\geq B$, and you can prove $A=B+C^2$, then, sure, the inequality follows from using "squares are nonnegative", but most of the time it is the identity that proves to be the hardest step. Here's an example, given reals $a_1,a_2,\dots, a_n$, you want to prove that $$\sum_{i,j=1}^n \frac{a_ia_j}{1+|i-j|} \geq 0.$$
After you realize that sum is just equal to $$\frac{1}{2\pi}\cdot\int_{0}^{2\pi}{\int_{0}^{1}{\frac{1-r^{2}}{1-2r\cos(x)+r^{2}}\cdot |\sum_{k=1}^{n}{a_{k}e^{-ikx}}|^{2}dx dr}}$$ then, yes, everything is obvious, but spotting the equality is clearly the nontrivial step in the proof.
In some instances it might be helpful to think about combinatorics, probability, algebra or geometry. Is the quantity $x$ enumerating objects you are familiar with, the probability of an event, the dimension of a vector space, or the area/volume of a region? There is plenty of inequalities that follow this way. Think of Littlewood-Richardson coeficients for example.
Another helpful factor is symmetry. Is your inequality invariant under permuting some of its variables? While I don't remember right now the paper, Polya has an article where he talks about the "principle of nonsufficient reason", which basically boils down to the strategy that if your function is symmetric enough, then so are it's extremal points (there is no sufficient reason to expect assymetry in the maximal/minimal points, is how he puts it). This is similar in vein to using Langrange multipliers. Note however that sometimes it is the oposite of this that comes in handy. Schur's inequality, for example is known to be impossible to prove using "symmetric methods", one must break the symmetry by assuming an arbitrary ordering on the variables. (I think it was sent by Schur to Hardy as an example of a symmetric polynomial inequality that doesn't follow from Muirhead's theorem, see below.)
Majorization theory is yet another powerful tool. The best reference that comes to mind is Marshall and Olkin's book "Inequalities: Theory of Majorization and Its Applications". This is related to what you call convexity and some other notions. Note that there is a lot of literature devoted to inequalities involving "almost convex" functions, where a weaker notion than convexity is usually used. Also note the concepts of Schur-convexity, quasiconvexity, pseudoconvexity etc. One of the simplest applications of majorization theory is Muirhead's inequality which generalizes already a lot of classical inequalities and inequalities such as the ones that appear in competitions.
Sometimes you might want to take advantage of the duality between discrete and continuous. So depending on which tools you have at your disposal you may choose to prove, say the inequality $$\sum_{n=1}^{\infty}\left(\frac{a_1+\cdots+a_n}{n}\right)^p\le \left(\frac{p}{p-1}\right)^p \sum_{n=1}^{\infty}a_n^p$$ or it's continuous/integral version
$$\int_{0}^{\infty}\left(\frac{1}{x}\int_0^x f(t)dt\right)^p dx \le \left(\frac{p}{p-1}\right)^p \int_{0}^{\infty} f(x)^p dx$$ I've found this useful in different occasions (in both directions).
Other things that come to mind but that I'm too lazy to describe are "integration preserves positivity", uncertainity principle, using the mean value theorem to reduce the number of variables etc. What also comes in handy, sometimes, is searching if others have considered your inequality before. This might prevent you from spending too much time on an inequality like $$\sum_{d|n}d \le H_n+e^{H_n}\log H_n$$ where $H_n=\sum_{k=1}^n \frac{1}{k}$.
OK, having spent about 20 hours on the search of a nice proof (which extinguished my passion for beauty for the next several days at least), I'm resorting to the brute force. I will love to see someone else to avenge this pitiful defeat of mine...
As I mentioned in the comment, the key to the solution is the inequality
$$
(1+a)\log(1+b)+(1+b)\log(1+a)\ge 2(1+c)\log(1+c)
$$
where $a,b,c>0, ab=c^2$.
Consider the function
$$
F(t)=(1+a)\log(1+bt)+(1+b)\log(1+at)-2(1+c)\log(1+ct)\,.
$$
We have
$$
\begin{aligned}
F'(t)&=\frac{(1+a)b}{1+bt}+\frac{(1+b)a}{1+at}-2\frac{(1+c)c}{1+ct}
\\
&=(a+b-2c)\frac{(c^2t-1)(ct-1)}{(1+at)(1+bt)(1+ct)}
\end{aligned}
$$
While the derivation of the last equality is tedious, to check it, it suffices to show that $F'(c^{-2})=F'(c^{-1})=0$ and to compute the free term in the numerator to get the right normalization factor, so I'll skip the explicit computations.
It will be convenient to denote $p=c^{-1}$. Then we need to show that
$$
\int_0^1\frac{(t-p^2)(t-p)}{(1+at)(1+bt)(1+ct)}\,dt\ge 0\,.
$$
If $p>1$, there is nothing to do because the integrand is non-negative. So, we'll assume that $p<1$ from now on.
Since the denominator is increasing in $t$, and the numerator is positive on $(0,p^2)$, negative on $(p^2,p)$ and then positive again on $(p,1)$, we would be in good shape if we had
$$
\int_0^p(t-p^2)(t-p)\,dt=-\frac{p^3}6+\frac{p^4}2\ge 0\,,
$$
which is true for $p\ge \frac 13$. Thus, we may assume that $p<\frac 13$.
Now we cannot just neglect the interval $[p,1]$. We shall neglect the interval $[0,p^2]$ instead. Assuming $a\ge c\ge b$, we write
$$
\begin{multline}
\int_{p^2}^p\frac{(t-p^2)(t-p)}{(1+at)(1+bt)(1+ct)}\,dt=
\int_p^1\frac{(t-p)(t-1)p^2}{(p^{-1}+at)(1+bpt)(1+cpt)}\,dt
\\
\ge
\int_p^1\frac{(t-p)(t-1)p^2}{(1+at)(1+bpt)(1+t)}\,dt
\end{multline}
$$
Combining this with $\int_p^1$, we get the lower bound
$$
\int_p^1\frac{(t-p)}{(1+at)}\left[\frac{(t-1)p^2}{(1+bpt)(1+t)}
+\frac{t-p^2}{(1+bt)(1+ct)}\right]\,dt
$$
Note now that $bp\le 1$ and the higher $b$ is, the harder it is for the last square bracket to be non-negative. Thus, it is enough to estimate the bracket for $bp=1$, i.e., to demonstrate that
$$
(1-t)(p+t)^2\le (t-p^2)(1+t)^2\,.
$$
However
$$
t-p^2\ge t(1-t)
$$
and
$$
(p+t)^2=2pt+p^2+t^2<t+2t^2<t(1+t)^2\,.
$$
I would appreciate it if someone checks this ugly monster before I post the remaining (almost trivial) part of the proof :).
Edit: Assuming that those who upvoted took trouble to read and verify the above part, the end is as follows.
Let $X=X_j$, $Y=Y_j$, and $K$ be as before. Then, as I said, $\sqrt{XY}\le K\le\sqrt{XY}+\sqrt{(1-X)(1-Y)}$ (Cauchy). We have $Z=Z_j=\frac{\sqrt{XY}}K$, so it will suffice to show that
$$
2K^2 Z\log\frac 1Z=2K\sqrt{XY}\log\frac{K}{\sqrt{XY}}\le X\log\frac 1X+Y\log\frac 1Y\,.
$$
Note that the left hand side is convex in $K$ for fixed $X,Y$. If $K=\sqrt{XY}$, the inequality is trivial ($0$ is less than or equal to a non-negative number). Thus, we need only consider the case $K=\sqrt{XY}+\sqrt{(1-X)(1-Y)}$. Dividing by $XY$ and putting $X^{-1}=1+x, Y^{-1}=1+y$, so that $\frac{1-X}X=x,\frac{1-Y}Y=y$, we see that this case reduces exactly to the above inequality.
P.S. There is an alternative proof of the main inequality on AoPS. Alas, it also requires some tedious computations...
Best Answer
Take any $r\in(-1,0)$, any vector $(p_i)_{i=1}^n$ with $p_i>0$ for all $i$, and any vector $(h_i)_{i=1}^n\in\mathbb R^n$. For all real $t$ close enough to $0$, let $$g(t):=M(p+th).$$ Then $$g''(0)=(1+r)\Big(\sum_{i=1}^n p_i^{r-1}h_i^2\,\sum_{i=1}^n p_i^{1+r} +\frac{1-r^2}{r^2}\,\Big(\sum_{i=1}^n p_i^r h_i\Big)^2\Big) \Big(\sum_{i=1}^n p_i^{1+r}\Big)^{1/r-2},$$ which is manifestly $\ge0$.
It follows that $M$ is indeed convex (actually, convex on the entire positive orthant of $\mathbb R^n$).