The latest Martian gravity and topographic models derived from the Mars Orbiter Laser Altimeter and the Mars Global Surveyor spacecraft tracking data are used to compute the sub-crustal stress field on Mars. For this purpose, we apply the method for a simultaneous determination of the horizontal sub-crustal stress component and the crustal thickness based on solving the Navier–Stokes problem and incorporating the Vening Meinesz–Moritz inverse problem of isostasy. Results reveal that most of the Martian sub-crustal stress is concentrated in the Tharsis region, with the most prominent signatures attributed to a formation of Tharsis major volcanoes followed by crustal loading. The stress distribution across the Valles Marineris rift valleys indicates extensional tectonism. This finding agrees with more recent theories of a tectonic origin of Valles Marineris caused, for instance, by a crustal loading of the Tharsis bulge that resulted in a regional trusting and folding. Aside from these features, the Martian stress field is relatively smooth with only a slightly enhanced pattern of major impact basins. The signatures of active global tectonics and polar ice load are absent. Whereas the signature of the hemispheric dichotomy is also missing, the long-wavelength spectrum of the stress field comprises the signature of additional dichotomy attributed to the isostatically uncompensated crustal load of Tharsis volcanic accumulations. These results suggest a different origin of the Earth's and Martian sub-crustal stress. Whereas the former is mainly related to active global tectonics, the latter is generated by a crustal loading and regional tectonism associated with a volcanic evolution on Mars. The additional sub-crustal stress around major impact basins is likely explained by a crustal extrusion after impact followed by a Moho uplift.