In elastohydrodynamic lubrication (EHL), the lubricant experiences pressures in excess of 500 MPa and strain rates larger than $$10^5$$
10
5
$$\text{s}^{-1}$$
s
-
1
. The high pressures lead to a dramatic rise in the Newtonian viscosity and the high rates lead to large shear stresses and pronounced shear thinning. The extraction of accurate rheological properties using non-equilibrium molecular dynamics simulations (NEMD) has played a key role in improving our understanding of lubricant flow in EHL conditions. However, the high dimensionality of the output data generated by NEMD simulations often makes a deeper interrogation of the link between molecular-scale features and rheological properties challenging. In this paper, we explore the use of machine learning to analyze and visualize the high-dimensional output data generated in typical NEMD simulations. We show that dimension reduction of NEMD simulation data describing the shear flow of squalane enables a clear visualization of the transition from Newtonian to non-Newtonian shear thinning with increasing shear rate and provides a reliable assessment of the correlation between shear thinning and the evolution in molecular alignment. The end-to-end atom pairs dominate the largest variations in pair orientation tensor components for low-pressure systems (0.1, 100 MPa) and provide the clearest separation of the orientation tensors with rate. On the other hand, the side atom pairs dominate the largest variation in the tensor components for the high-pressure systems ($$P\ge 400$$
P
≥
400
MPa) which exhibit an overall limited evolution in orientation tensors as a function of rate. Dimension reduction using all the six components of the orientation tensors of all 435 pairs associated with a squalane molecule shows that the decrease in viscosity with rate for low pressures is strongly correlated with changes in molecular alignment. However, for high pressures, shear thinning occurs at saturated orientational order.