We study active array imaging of small but strong scatterers in homogeneous media when multiple scattering between them is important. We use the Foldy-Lax equations to model wave propagation with multiple scattering when the scatterers are small relative to the wavelength. In active array imaging we seek to locate the positions and reflectivities of the scatterers, that is, to determine the support of the reflectivity vector and the values of its nonzero elements from echoes recorded on the array. This is a nonlinear inverse problem because of the multiple scattering. In this paper, we show how to avoid the nonlinearity and form images noniteratively through a two-step process which involves l norm minimization. However, under certain illuminations imaging may be affected by screening, where some scatterers are obscured by multiple scattering. This problem can be mitigated by using multiple and diverse illuminations. In this case, we determine solution vectors that have a common support. The uniqueness and stability of the support of the reflectivity vector obtained with single or multiple illuminations are analyzed, showing that the errors are proportional to the amount of noise in the data with a proportionality factor dependent on the sparsity of the solution and the mutual coherence of the sensing matrix, which is determined by the geometry of the imaging array. Finally, to filter out noise and improve the resolution of the images, we propose an approach that combines optimal illuminations using the singular value decomposition of the response matrix together with sparsity promoting optimization jointly for all illuminations. In this work, an extension of our previous paper [A. Chai, M. Moscoso, and G. Papanicolaou, Inverse Problems, 29 (2013), 025016] on imaging using optimization techniques, we now account for multiple scattering effects.