The association of seismic‐wave arrivals with causative earthquakes becomes progressively more challenging as arrival detection methods become more sensitive, and particularly when earthquake rates are high. For instance, seismic waves arriving across a monitoring network from several sources may overlap in time, false arrivals may be detected, and some arrivals may be of unknown phase (e.g., P or S waves). We propose an automated method to associate arrivals with earthquake sources and obtain source locations applicable to such situations. To do so, we use a pattern detection metric based on the principle of backprojection to reveal candidate sources followed by graph‐theory‐based clustering and an integer linear optimization routine to associate arrivals with the minimum number of sources necessary to explain the data. This method solves for all sources and phase assignments simultaneously, rather than in a sequential greedy procedure as is common in other association routines. We demonstrate our method on both synthetic and real data from the Integrated Plate Boundary Observatory Chile seismic network of northern Chile. For the synthetic tests, we report results for cases with varying complexity, including rates of 500 earthquakes/day and 500 false arrivals/station/day, for which we measure true positive detection accuracy of >95%. For the real data, we develop a new catalog between 1 January 2010 and 31 December 2017 containing 817,548 earthquakes, with detection rates on average 279 earthquakes/day and a magnitude‐of‐completion of M∼1.8⁠. A subset of detections are identified as sources related to quarry and industrial site activity, and we also detect thousands of foreshocks and aftershocks of the 1 April 2014 M_w 8.2 Iquique earthquake. During the highest rate of aftershock activity, >600 earthquakes/day are detected in the vicinity of the Iquique earthquake rupture zone.