In photogrammetry, this is a pretty standard problem, so we can consider it solved. ;-)
But as I already said: It is not simple.
You can probably get away with a reasonably good approximation by doing the following:
1) Solve the equation system for each pair of projected and unprojected points. This yields a transformation matrix for each pair.
2) Choose a good one from these by projecting all the points with each matrix and pick the one that produces the least error.
Do you mean I should compute such a transformation matrix which is one part of the modelview matrix? If yes, then the projection matrix and factors should be considered fixed for this problem, it is right?