Optimizing vector values for maximum correlation
This is one approach you can follow:


Setup a linear regression system, with each match as a row, and each feature corresponding to a team.
The unknown feature coefficient for each team is the team 'strength' that we try to determine.
The feature values will be one of 0, 1, or -1, depending on whether the team did not play, was the column team, or the row team respectively in that match.
The regression target will be the score differential (column team score - row team score) in that match.


Eg: For the result matrix above, the system would be:

$$
(-1) * x_1 + (1) * x_2 + (0) * x_3 = 3\\
(-1) * x_1 + (0) * x_2 + (1) * x_3 = 6\\
(0) * x_1 + (0-1) * x_2 + (1) * x_3 = 3
$$


One solution to the above system is: $x_1=-6; x_2=-3; x_3=0$
Multiple solutions are possible, which can be viewed as translations of each other (adding a constant to all team strengths).
If there are 'n' teams, then there are only 'n-1' linearly independent columns in the regression. (In R, one of the coefficients comes out as NA. This can be treated as 0, or dropped from the regression, effectively making it 0).