Lesion to symptom mapping (LSM) is a crucial tool for understanding the causality of brain-behavior relationships. The analyses are typically performed by applying statistical methods on individual brain voxels (VLSM), a method called the mass-univariate approach. Several authors have shown that VLSM suffers from limitations that may decrease the accuracy and reliability of the findings, and have proposed the use of multivariate methods to overcome these limitations. In this study, we propose a multivariate optimization technique known as sparse canonical correlation analysis for neuroimaging (SCCAN) for lesion to symptom mapping. To validate the method and compare it with mass-univariate results, we used data from 131 patients with chronic stroke lesions in the territory of the middle cerebral artery, and created synthetic behavioral scores based on the lesion load of 93 brain regions (putative functional units). LSM analyses were performed with univariate VLSM or SCCAN, and the accuracy of the two methods was compared in terms of both overlap and displacement from the simulated functional areas. Overall, SCCAN produced more accurate results - higher dice overlap and smaller average displacement - compared to VLSM. This advantage persisted at different sample sizes (N = 20–131) and different multiple comparison corrections (false discovery rate, FDR; Bonferroni; permutation-based family wise error rate, FWER). These findings were replicated with a fully automated SCCAN routine that relied on cross-validated predictive accuracy to find the optimal sparseness value. Simulations of one, two, and three brain regions showed a systematic advantage of SCCAN over VLSM; under no circumstance could VLSM exceed the accuracy obtained with SCCAN. When considering functional units composed of multiple brain areas VLSM identified fewer areas than SCCAN. The investigation of real scores of aphasia severity (aphasia quotient and picture naming) showed that SCCAN could accurately identify known language-critical areas, while VLSM either produced diffuse maps (FDR correction) or few scattered voxels (FWER correction). Overall, this study shows that a multivariate method, such as, SCCAN, outperforms VLSM in a number of scenarios, including functional dependency on single or multiple areas, different sample sizes, different multi-area combinations, and different thresholding mechanisms (FWER, Bonferroni, FDR). These results support previous claims that multivariate methods are in general more accurate than mass-univariate approaches, and should be preferred over traditional VLSM approaches. All the methods described in this study are available in the newly developed LESYMAP package for R.