An Average of Probabilities This is an open ended question in that the desired end result is not well posed.  Still, it may be of some interest.
Suppose you have a number of teams which play against each other in two team competitions (all pairings occur with equal likehood, no ties allowed).  Let $P(A,B)$ be the probability that A triumphs in a match between A and B.  Clearly, if you know all the $P(A,B)$ you can compute the team average $P(A)$, the probability that A will win against an unknown opponent.  You can't go the other way, though.  To see this, look at 3 teams $A,B,C$ and suppose $P(A) = P(B) = P(C) = \frac{1}{2}$.  One way we might have gotten this outcome is if $P(X,Y) = \frac{1}{2}$ for any pair of teams X,Y.  In some vague sense, I suppose this is the "most natural" guess.  But it is not the only possibility.  We might, for instance,  be in "Rock, Paper, Scissors" world, where A always beats B, B always beats C, and C always beats A.  
The Question:  Given all the P(A), is there a reasonably natural way to produce a   list of probabilities P(A, *)?
Of course the probabilities must be consistent (so that the list of P(A,B)'s does imply the given list of team averages).
Sample attempt:  attach a "greatness factor" $\lambda_i$ to each team and then define $$P(A_i,A_j) = .5 + \lambda_i - \lambda_j$$
You can even make the $\lambda$'s sum to 1.  This is easily computable and gives plausible results so long as all the team averages are reasonably near .5.  In more extreme cases, though, you get unphysical results (probabilities greater than 1 and so on). For baseball statistics, where this problem first arose, this crude method works fairly well.
Absent any actual ideas, my impulse is to minimize variance.  That is: "choose the solution which is as nearly constant as possible."  That is obviously one way to go, but perhaps someone has a better idea?
 A: Your solution $P(A_i,A_j) = .5 + \lambda_i - \lambda_j$ is the stationary point of $\sum_{ij}(P(A_i,A_j)-.5)^2$ under the given equality constraints. (I guess that's how you derived it, so it's no coincidence that your "greatness factors" look like Lagrange parameters :-)
However, this is not the solution to the complete problem, since we also have the inequality constraints $0\le P(A_i,A_j)\le1$. You mentioned that these are not automatically satisfied, but to solve the problem of minimizing the objective function under the given constraints, we have to make them satisfied.
This would generally be a quadratic programming problem, but I suspect that in the present case you can get away with simply deactivating the variables one by one. I haven't tried this, but the prescription would be: Calculate $P(A_i,A_j) = .5 + \lambda_i - \lambda_j$, determining the $\lambda_i$ up to an additive constant from the equality contraints: $\lambda_i=P(A_i)-.5+\left<\lambda_j\right>_{j\ne i}$. Set all negative probabilities to $0$ (and thus, by symmetry, all that are above $1$ to $1$). Then recalculate the $\lambda_i$; but now the $P(A_i,A_j)$ that hit the inequality constraints are no longer given by $.5+\lambda_i-\lambda_j$, but by their fixed values $0$/$1$. So you get a different system of equations for $\lambda_i$. If you're lucky, solving it and recalculating the probabilities will lead to the same probabilities exceeding the constraints, and then you're done. If not, you can either iterate until you get consistency, or if that doesn't work, use a proper quadratic programming algorithm. In any case, in the end you'll find the solution to your minimization problem under the given constraints.
My personal favorite, by the way, for a "natural" determination of the pairwise probabilities, would be to model a team's performance in any of their games by a Gaussian, with identical variances for all teams but different means. Then you have $n$ parameters (the means) to fit to the $n$ given data (the $P(A_i)$), with one "zero mode" in each (the probabilities have to average to $.5$ and shifting all the means by an additive constant makes no difference). The fitting isn't as nicely linear as in your approach, and I'm not even sure there would always be a solution, but if there is, then it gives you a nice model that to me feels closer to how a game is actually decided than the "greatness factors".
A: The Gaussian performance model of joriki seems nice, but a little expensive computationally. In contrast the additional constraints of the greatness model make it a little unwieldy.
I've previously used a model of the form,
$$
P(A,B) = \frac{\lambda_A}{\lambda_A+\lambda_B}
$$
with the additional constraints that $\lambda_X > 0$ and $\sum\lambda_X =1$
(needed since scaling all $\lambda_X$ by a constant gives the same probabilities).
Note: I was not doing exactly the same as you, I was trying to estimate the $P(A,B)$ from a limited number of observations (and used a maximum likelihood method), so I'm not certain this form will lead to anything solvable analytically, or even produce nice results for your data.
