Total variation (TV) and wavelet L1 norms have often been used as regularization terms in image restoration and reconstruction problems. However, TV regularization can introduce staircase effects and wavelet regularization some ringing artifacts, but hybrid TV and wavelet regularization can reduce or remove these drawbacks in the reconstructed images. We need to compute the proximal operator of hybrid regularizations, which is generally a sub-problem in the optimization procedure. Both TV and wavelet L1 regularisers are nonlinear and non-smooth, causing numerical difficulty. We propose a dual iterative approach to solve the minimization problem for hybrid regularizations and also prove the convergence of our proposed method, which we then apply to the problem of MR image reconstruction from highly random under-sampled k-space data. Numerical results show the efficiency and effectiveness of this proposed method.