Mapping landforms on the Moon is of great interest and importance for future human settlements and resources exploration. One of the first steps is to map the topography in great detail and resolution. However, data from the Lunar Orbiter Laser Altimeter (LOLA) provide low-resolution elevation maps in comparison to the size of detailed geological features. To improve resolution, we developed a new method to upscale topographic maps to a higher resolution using images from the Lunar Reconnaissance Orbiter Camera (LROC). Our method exploits the relation between topographic gradients and degrees of shading of incoming sunlight. In contrast to earlier published methods, our approach is based on probabilistic, linear inverse theory, and its computational efficiency is very high due to its formulation through the Sylvester Equation. The method operates on multiple images and incorporates albedo variations. A further advantage of the method is that we avoid/reduce the use of arbitrary tuning parameters through a probabilistic formulation where all weighting of data and model parameters is based on prior information about data uncertainties and reasonable bounds on the model. Our results increase the resolution of the topography from -60 m per pixel to 0.9 m per pixel, bringing it to the same pixel resolution as the optical images from LROC, allowing in some cases detection of craters as small as -3 m of diameter. We estimate uncertainties of the topographic model due to noise in the images, and in the low-resolution (LOLA) model.