To accelerate the acquisition of J-resolved proton magnetic resonance spectroscopic imaging (1 H-MRSI) data for high-resolution mapping of brain metabolites and neurotransmitters. Methods: The proposed method used a subspace model to represent multidimensional spatiospectral functions, which significantly reduced the number of parameters to be determined from J-resolved 1 H-MRSI data. A semi-LASER-based (Localization by Adiabatic SElective Refocusing) echo-planar spectroscopic imaging (EPSI) sequence was used for data acquisition. The proposed data acquisition scheme sampled k, t 1 , t 2-space in variable density, where t 1 and t 2 specify the J-coupling and chemical-shift encoding times, respectively. Selection of the J-coupling encoding times (or, echo time values) was based on a Cramer-Rao lower bound analysis, which were optimized for gamma-aminobutyric acid (GABA) detection. In image reconstruction, parameters of the subspace-based spatiospectral model were determined by solving a constrained optimization problem. Results: Feasibility of the proposed method was evaluated using both simulated and experimental data from a spectroscopic phantom. The phantom experimental results showed that the proposed method, with a factor of 12 acceleration in data acquisition, could determine the distribution of J-coupled molecules with expected accuracy. In vivo study with healthy human subjects also showed that 3D maps of brain metabolites and neurotransmitters can be obtained with a nominal spatial resolution of 3.0 × 3.0 × 4.8 mm 3 from J-resolved 1 H-MRSI data acquired in 19.4 min. Conclusions: This work demonstrated the feasibility of highly accelerated J-resolved 1 H-MRSI using limited and sparse sampling of k, t 1 , t 2-space and subspace modeling. With further development, the proposed method may enable high-resolution mapping of brain metabolites and neurotransmitters in clinical applications.