A semi-inclusive measurement of charged hadron multiplicities in deep inelastic muon scattering off an isoscalar target was performed using data collected by the COMPASS Collaboration at CERN. The following kinematic domain is covered by the data: photon virtuality $Q^{2}>1$ (GeV/$c$)$^2$, invariant mass of the hadronic system $W > 5$ GeV/$c^2$, Bjorken scaling variable in the range $0.003 < x < 0.4$, fraction of the virtual photon energy carried by the hadron in the range $0.2 < z < 0.8$, square of the hadron transverse momentum with respect to the virtual photon direction in the range 0.02 (GeV/$c)^2 < P_{\rm{hT}}^{2} < 3$ (GeV/$c$)$^2$. The multiplicities are presented as a function of $P_{\rm{hT}}^{2}$ in three-dimensional bins of $x$, $Q^2$, $z$ and compared to previous semi-inclusive measurements. We explore the small-$P_{\rm{hT}}^{2}$ region, i.e. $P_{\rm{hT}}^{2} < 1$ (GeV/$c$)$^2$, where hadron transverse momenta are expected to arise from non-perturbative effects, and also the domain of larger $P_{\rm{hT}}^{2}$, where contributions from higher-order perturbative QCD are expected to dominate. The multiplicities are fitted using a single-exponential function at small $P_{\rm{hT}}^{2}$ to study the dependence of the average transverse momentum $\langle P_{\rm{hT}}^{2}\rangle$ on $x$, $Q^2$ and $z$. The power-law behaviour of the multiplicities at large $P_{\rm{hT}}^{2}$ is investigated using various functional forms. The fits describe the data reasonably well over the full measured range.