This work is concerned with the structure of bilinear minimization problems arising in recovering subsampled and modulated images in parallel magnetic resonance imaging. By considering a physically reasonable simplified model exhibiting the same fundamental mathematical difficulties, it is shown that such problems suffer from poor gradient scaling and non-convexity, which causes standard optimization methods to perform inadequately. A globalized quasi-Newton method is proposed which is able to reconstruct both image and the unknown modulations without additional a priori information. Thus the present paper serves as a first contribution toward understanding and solving such bilinear optimization problems.