Newtonian machine learning (NML) inversion has been shown to accurately recover the low-to-intermediate wavenumber information of subsurface velocity models. This method uses the wave-equation inversion kernel to invert the skeletonized data that is automatically learned by an autoencoder. The skeletonised data is a one-dimensional latent-space representation of the seismic trace. However, for a complicated dataset, the decoded waveform could lose some details if the latent space dimension is set to one, which leads to a low-resolution NML tomogram. To mitigate this problem, an autoencoder with a higher dimensional latent space is needed to encode and decode the seismic data. In this paper, we present a wave equation inversion that inverts the multi-dimensional latent variables of an autoencoder for the subsurface velocity model. The multi-variable implicit function theorem is used to determine the perturbation of the multi-dimensional skeletonised data with respect to the velocity perturbations. In this case, each dimension of the latent variable is characterized one gradient and the velocity model is updated by the weighted sum of all these gradients. Numerical results suggest that the multidimensional NML inverted result can achieve a higher resolution in the tomogram compared to the conventional single dimensional NML inversion.