Linear mixture model (LMM) is widely used in remote sensing image analysis. In the original LMM, it is assumed that each pixel is a linear combination of all the endmember signatures. For an image scene covering a large geospatial area, the number of endmembers is quite large but only a small subset of these endmembers may actually participate in the composition of a specific pixel. If all the endmembers are used in the spectral unmixing process, the result may contain additional estimation error. The Multiple Endmember Spectral Mixture Analysis (MESMA) approach has been proposed which allows the endmember subset to be varied from pixel to pixel. We have developed approaches that can automatically search for an endmember subset for each pixel from a large number of endmembers for the entire image scene, without the limitation on the number of endmembers to be included. Due to high computational complexity, in this paper, we will develop a fast searching algorithm by considering the neighborhood relationship among pixels. We will show that the resulting endmember sets can reduce pixel reconstruction error and improve the quality of estimated abundances (i.e., all the pixels have nonnegative abundances and most of pixels have sum-to-one abundances).