A method for reconstructing parallel magnetic resonance images includes providing a set of acquired k-space MR image data y, and finding a target MR image x that minimizes ½∥Fv−y∥22+λ∥z∥1 where v=Sx and z=Wx where S is a diagonal matrix containing sensitivity maps of coil elements in an MR receiver array, F is an FFT matrix, W is a redundant Haar wavelet matrix, and λ≧0 is a regularization parameter, by updating x k + 1 = ( μ 1  I + μ 3  S H  S ) - 1  [ μ 1  W H  ( z k - b z k ) + μ 3  S H  ( v k - b v k ) ] ,  z k + 1 = soft  ( Wx k + 1  b z k , 1 μ 1 )   where soft  ( x , T ) = { x + T if   x ≤ - T , 0 if    x  ≤ T , x - T if   x ≥ T ,   and   v k + 1 = ( F H  F + μ 3  I ) - 1  [ F H  y + μ 3  ( Sx k + 1 + b v k ) ] , where k is an iteration counter, μ1 and μ3 are parameters of an augmented Lagrangian function, and bz and bv are dual variables of the augmented Lagrangian.