The p‐median problem (PMP) consists of locating p facilities (medians) in order to minimize the sum of distances from each client to the nearest facility. The interest in the large‐scale PMP arises from applications in cluster analysis, where a set of patterns has to be partitioned into subsets (clusters) on the base of similarity. In this paper we introduce a new heuristic for large‐scale PMP instances, based on Lagrangean relaxation. It consists of three main components: subgradient column generation, combining subgradient optimization with column generation; a ‘core’ heuristic, which computes an upper bound by solving a reduced problem defined by a subset of the original variables chosen on a base of Lagrangean reduced costs; and an aggregation procedure that defines reduced size instances by aggregating together clients with the facilities. Computational results show that the proposed heuristic is able to compute good quality lower and upper bounds for instances up to 90,000 clients and potential facilities.