Objective: We propose and develop a novel biclustering (N-BiC) approach for performing N-way biclustering of neuroimaging data. Our approach is applicable to an arbitrary number of features from both imaging and behavioral data (e.g., symptoms). We applied it to structural MRI data from patients with schizophrenia. Methods: It uses a source-based morphometry approach [i.e., independent component analysis of gray matter segmentation maps] to decompose the data into a set of spatial maps, each of which includes regions that covary among individuals. Then, the loading parameters for components of interest are entered to an exhaustive search, which incorporates a modified depth-first search technique to carry out the biclustering, with the goal of obtaining submatrices where the selected rows (individuals) show homogeneity in their expressions of selected columns (components) and vice versa. Results: Findings demonstrate that multiple biclusters have an evident association with distinct brain networks for the different types of symptoms in schizophrenia. The study identifies two components: inferior temporal gyrus (16) and brainstem (7), which are related to positive (distortion/excess of normal function) and negative (diminution/loss of normal function) symptoms in schizophrenia, respectively. Conclusion: N-BiC is a data-driven method of biclustering MRI data that can exhaustively explore relationships/substructures from a dataset without any prior information with a higher degree of robustness than earlier biclustering applications. Significance: The use of such approaches is important to investigate the underlying biological substrates of mental illness by grouping patients into homogeneous subjects, as the schizophrenia diagnosis is known to be relatively nonspecific and heterogeneous.