1. A method for reconstructing parallel magnetic resonance images (MRI), comprising: providing a set of acquired k-space MR image data y; and finding a target MR image x that minimizes ½∥Fv−y∥ 2 2 +λ∥z∥ 1 wherein v=Sx and z=Wx wherein S is a diagonal matrix containing sensitivity maps of coil elements in an MR receiver array, F is an FFT matrix wherein FF H =I and H denotes a Hermitian adjoint matrix, W is a redundant Haar wavelet matrix satisfying W T W=I, wherein I is an identity 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 ( W ⁢ x k + 1 + b z k , 1 μ 1 ) wherein ⁢ ⁢ 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 ) ] , wherein k is an iteration counter, μ 1 and μ 3 are parameters of an augmented Lagrangian function L ⁡ ( x , z , v , b z , b v ) = 1 2 ⁢  Fv - y  2 2 + λ ⁢  z  1 + μ 1 2 ⁢  Wx - z  2 2 + μ 3 2 ⁢  Sx - v  2 2 + μ 1 ⁢ 〈 b z , Wx - z + y 〉 + μ 3 ⁢ 〈 b v , Sx - v 〉 , and b z and b v are dual variables of the augmented Lagrangian.