diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py
index a35f97c0..69e85e84 100644
--- a/SimPEG/Directives.py
+++ b/SimPEG/Directives.py
@@ -144,6 +144,7 @@ class BetaSchedule(InversionDirective):
             if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
             self.invProb.beta /= self.coolingFactor
 
+
 class TargetMisfit(InversionDirective):
 
     chifact = 1.
@@ -242,12 +243,6 @@ class SaveOutputDictEveryIteration(SaveEveryIteration):
         # Save the file as a npz
         np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
 
-
-# class UpdateReferenceModel(Parameter):
-
-#     mref0 = None
-
-#     def nextIter(self):
 #         mref = getattr(self, 'm_prev', None)
 #         if mref is None:
 #             if self.debug: print 'UpdateReferenceModel is using mref0'
@@ -258,56 +253,138 @@ class SaveOutputDictEveryIteration(SaveEveryIteration):
 class Update_IRLS(InversionDirective):
 
     eps_min = None
+    eps_p = None
+    eps_q = None
+    norms = [2.,2.,2.,2.]
     factor = None
     gamma = None
     phi_m_last = None
     phi_d_last = None
+    f_old = None
+    f_min_change = 1e-2
+    beta_tol = 5e-2
 
+    # Solving parameter for IRLS (mode:2)
+    IRLSiter   = 0
+    minGNiter = 5
+    maxIRLSiter = 10
+    iterStart = 0
+
+    # Beta schedule
+    coolingFactor = 2.
+    coolingRate = 1
+
+    mode = 1
+
+    @property
+    def target(self):
+        if getattr(self, '_target', None) is None:
+            self._target = self.survey.nD*0.5
+        return self._target
+    @target.setter
+    def target(self, val):
+        self._target = val
 
     def initialize(self):
 
-        # Scale the regularization for changes in norm
-        if getattr(self, 'phi_m_last', None) is not None:
-
-            self.reg.curModel = self.invProb.curModel
-            self.reg.gamma = 1.
-            phim_new = self.reg.eval(self.invProb.curModel)
-            self.gamma = self.phi_m_last / phim_new
-
-        self.reg.curModel = self.invProb.curModel
-        self.reg.gamma = self.gamma
-
-        if getattr(self, 'phi_d_last', None) is None:
-            self.phi_d_last = self.invProb.phi_d
+        if self.mode == 1:
+            self.reg.norms = [2., 2., 2., 2.]
 
     def endIter(self):
-        # Cool the threshold parameter if required
-        if getattr(self, 'factor', None) is not None:
-            eps = self.reg.eps / self.factor
 
-            if getattr(self, 'eps_min', None) is not None:
-                self.reg.eps = np.max([self.eps_min,eps])
+        # After reaching target misfit with l2-norm, switch to IRLS (mode:2)
+        if self.invProb.phi_d < self.target and self.mode == 1:
+            print "Convergence with smooth l2-norm regularization: Start IRLS steps..."
+
+            self.mode = 2
+            print self.eps_p, self.eps_q, self.norms
+            self.reg.eps_p = self.eps_p
+            self.reg.eps_q = self.eps_q
+            self.reg.norms = self.norms
+            self.coolingFactor = 1.
+            self.coolingRate = 1
+            self.iterStart = self.opt.iter
+            self.phi_d_last = self.invProb.phi_d
+            self.phi_m_last = self.invProb.phi_m_last
+
+            self.reg.l2model = self.invProb.curModel
+            self.reg.curModel = self.invProb.curModel
+
+            if getattr(self, 'f_old', None) is None:
+                self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False)
+
+        # Beta Schedule
+        if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
+            if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
+            self.invProb.beta /= self.coolingFactor
+
+
+        # Only update after GN iterations
+        if (self.opt.iter-self.iterStart) % self.minGNiter == 0 and self.mode==2:
+
+            self.IRLSiter += 1
+
+            phim_new = self.reg.eval(self.invProb.curModel)
+            self.f_change = np.abs(self.f_old - phim_new) / self.f_old
+
+            print "Regularization decrease: %6.3e" % (self.f_change)
+
+            # Check for maximum number of IRLS cycles
+            if self.IRLSiter == self.maxIRLSiter:
+                print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter
+                self.opt.stopNextIteration = True
+                return
+
+            # Check if the function has changed enough
+            if self.f_change < self.f_min_change and self.IRLSiter > 1:
+                print "Minimum decrease in regularization. End of IRLS"
+                self.opt.stopNextIteration = True
+                return
             else:
-                self.reg.eps = eps
+                self.f_old = phim_new
 
-        # Get phi_m at the end of current iteration
-        self.phi_m_last = self.invProb.phi_m_last
+            # Cool the threshold parameter if required
+            if getattr(self, 'factor', None) is not None:
+                eps = self.reg.eps / self.factor
 
-         # Update the model used for the IRLS weights
-        self.reg.curModel = self.invProb.curModel
+                if getattr(self, 'eps_min', None) is not None:
+                    self.reg.eps = np.max([self.eps_min,eps])
+                else:
+                    self.reg.eps = eps
 
-        # Temporarely set gamma to 1. to get raw phi_m
-        self.reg.gamma = 1.
+            # Get phi_m at the end of current iteration
+            self.phi_m_last = self.invProb.phi_m_last
 
-        # Compute new model objective function value
-        phim_new = self.reg.eval(self.invProb.curModel)
+            # Reset the regularization matrices so that it is
+            # recalculated for current model
+            self.reg._Wsmall = None
+            self.reg._Wx = None
+            self.reg._Wy = None
+            self.reg._Wz = None
 
-        # Update gamma to scale the regularization between IRLS iterations
-        self.reg.gamma = self.phi_m_last / phim_new
+             # Update the model used for the IRLS weights
+            self.reg.curModel = self.invProb.curModel
 
-        # Set the weighting matrix to None so that it is recomputed next time
-        # it is called in the inversion
-        self.reg._W = None
+            # Temporarely set gamma to 1. to get raw phi_m
+            self.reg.gamma = 1.
+
+            # Compute new model objective function value
+            phim_new = self.reg.eval(self.invProb.curModel)
+
+            # Update gamma to scale the regularization between IRLS iterations
+            self.reg.gamma = self.phi_m_last / phim_new
+
+            # Reset the regularization matrices again for new gamma
+            self.reg._Wsmall = None
+            self.reg._Wx = None
+            self.reg._Wy = None
+            self.reg._Wz = None
+
+            # Check if misfit is within the tolerance, otherwise scale beta
+            val = self.invProb.phi_d / (self.survey.nD*0.5)
+
+            if np.abs(1.-val) > self.beta_tol:
+                self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
 
 class Update_lin_PreCond(InversionDirective):
     """
@@ -360,19 +437,3 @@ class Update_Wj(InversionDirective):
             JtJdiag = JtJdiag / max(JtJdiag)
 
             self.reg.wght = JtJdiag
-
-class Scale_Beta(InversionDirective):
-    """
-        Instead of a linear cooling schedule, beta is allowed to change based
-        on the ratio between the target misfit and the current data misfit. The
-        update is done only if the misfit is outside some threshold bounds.
-    """
-    tol = 0.05
-
-    def endIter(self):
-
-        # Check if misfit is within the tolerance, otherwise adjust beta
-        val = self.invProb.phi_d / (self.survey.nD*0.5)
-
-        if np.abs(1.-val) > self.tol:
-            self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py
index 06ef4be4..afd90525 100644
--- a/SimPEG/Examples/Inversion_IRLS.py
+++ b/SimPEG/Examples/Inversion_IRLS.py
@@ -1,7 +1,7 @@
 from SimPEG import *
 
 
-def run(N=200, plotIt=True):
+def run(N=100, plotIt=True):
     """
         Inversion: Linear Problem
         =========================
@@ -18,6 +18,8 @@ def run(N=200, plotIt=True):
     mesh = Mesh.TensorMesh([N])
 
     m0 = np.ones(mesh.nC) * 1e-4
+    mref = np.zeros(mesh.nC)
+
     nk = 10
     jk = np.linspace(1.,nk,nk)
     p = -2.
@@ -50,57 +52,47 @@ def run(N=200, plotIt=True):
     wr = np.sum(prob.G**2.,axis=0)**0.5
     wr = ( wr/np.max(wr) )
 
-    reg = Regularization.Simple(mesh)
-    reg.wght = wr
-
+#    reg = Regularization.Simple(mesh)
+#    reg.mref = mref
+#    reg.cell_weights = wr
+#
     dmis = DataMisfit.l2_DataMisfit(survey)
     dmis.Wd = 1./wd
-
-    opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4)
-    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
-    invProb.curModel = m0
-
-    beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
-    target = Directives.TargetMisfit()
-
+#
+#    opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4)
+#    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
+#    invProb.curModel = m0
+#
+#    beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
+#    target = Directives.TargetMisfit()
+#
     betaest = Directives.BetaEstimate_ByEig()
-    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
-
-
-    mrec = inv.run(m0)
-    ml2 = mrec
-    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
-
-    # Switch regularization to sparse
-    phim = invProb.phi_m_last
-    phid =  invProb.phi_d
+#    inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
+#
+#
+#    mrec = inv.run(m0)
+#    ml2 = mrec
+#    print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
+#
+#    # Switch regularization to sparse
+#    phim = invProb.phi_m_last
+#    phid =  invProb.phi_d
 
     reg = Regularization.Sparse(mesh)
+    reg.mref = mref
+    reg.cell_weights = wr
 
-#==============================================================================
-#     fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
-#     dmdx = reg.mesh.cellDiffxStencil * mrec
-#     plt.plot(np.sort(dmdx))
-#==============================================================================
-
-    #reg.recModel = mrec
-    reg.wght = np.ones(mesh.nC)
     reg.mref = np.zeros(mesh.nC)
-    reg.eps_p = 5e-2
-    reg.eps_q = 1e-2
-    reg.norms   = [0., 0., 2., 2.]
-    reg.wght = wr
+    eps_p = 5e-2
+    eps_q = 5e-2
+    norms   = [0., 0., 2., 2.]
 
-    opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3)
-    invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.)
-    beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
-    #betaest = Directives.BetaEstimate_ByEig()
-    target = Directives.TargetMisfit()
-    IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
+    opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
+    invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
+    update_Jacobi = Directives.Update_lin_PreCond()
+    IRLS = Directives.Update_IRLS( norms=norms,  eps_p=eps_p, eps_q=eps_q)
 
-    inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
-
-    m0 = mrec
+    inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
 
     # Run inversion
     mrec = inv.run(m0)
@@ -117,7 +109,7 @@ def run(N=200, plotIt=True):
         axes[0].set_title('Columns of matrix G')
 
         axes[1].plot(mesh.vectorCCx, mtrue, 'b-')
-        axes[1].plot(mesh.vectorCCx, ml2, 'r-')
+        axes[1].plot(mesh.vectorCCx, reg.l2model, 'r-')
         #axes[1].legend(('True Model', 'Recovered Model'))
         axes[1].set_ylim(-1.0,1.25)
 
diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py
index 706124c8..c66242aa 100644
--- a/SimPEG/Optimization.py
+++ b/SimPEG/Optimization.py
@@ -1008,4 +1008,4 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
             indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
             delx[indx] = 0.
 
-            return delx
+            return delx
\ No newline at end of file
diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py
index 96beabcb..2ce5614c 100644
--- a/SimPEG/Regularization.py
+++ b/SimPEG/Regularization.py
@@ -1,4 +1,6 @@
-import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp
+import Utils, Maps, Mesh
+import numpy as np
+import scipy.sparse as sp
 
 class RegularizationMesh(object):
     """
@@ -403,7 +405,238 @@ class BaseRegularization(object):
 
         return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
 
-class Tikhonov(BaseRegularization):
+class Simple(BaseRegularization):
+    """
+    Simple regularization that does not include length scales in the derivatives.
+    """
+
+    mrefInSmooth = False  #: include mref in the smoothness?
+    alpha_s      = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
+    alpha_x      = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'],     "Weight for the first derivative in the x direction")
+    alpha_y      = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'],     "Weight for the first derivative in the y direction")
+    alpha_z      = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'],     "Weight for the first derivative in the z direction")
+    cell_weights = 1.
+
+    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
+        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+
+        if isinstance(self.cell_weights,float):
+            self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights
+
+    @property
+    def Wsmall(self):
+        """Regularization matrix Wsmall"""
+        if getattr(self,'_Wsmall', None) is None:
+            self._Wsmall = Utils.sdiag((self.alpha_s*self.cell_weights)**0.5)
+        return self._Wsmall
+
+    @property
+    def Wx(self):
+        """Regularization matrix Wx"""
+        if getattr(self, '_Wx', None) is None:
+            self._Wx = Utils.sdiag((self.alpha_x * (self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.regmesh.cellDiffxStencil
+        return self._Wx
+
+    @property
+    def Wy(self):
+        """Regularization matrix Wy"""
+        if getattr(self, '_Wy', None) is None:
+            self._Wy = Utils.sdiag((self.alpha_y * (self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.regmesh.cellDiffyStencil
+        return self._Wy
+
+    @property
+    def Wz(self):
+        """Regularization matrix Wz"""
+        if getattr(self, '_Wz', None) is None:
+            self._Wz = Utils.sdiag((self.alpha_z * (self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.regmesh.cellDiffzStencil
+        return self._Wz
+
+#    @property
+#    def Wsmooth(self):
+#        """Full smoothness regularization matrix W"""
+#        print 'wtf why are we using Wsmooth'
+#        raise NotImplementedError
+#        if getattr(self, '_Wsmooth', None) is None:
+#            wlist = (self.Wx,)
+#            if self.regmesh.dim > 1:
+#                wlist += (self.Wy,)
+#            if self.regmesh.dim > 2:
+#                wlist += (self.Wz,)
+#            self._Wsmooth = sp.vstack(wlist)
+#        return self._Wsmooth
+#
+#    @property
+#    def W(self):
+#        """Full regularization matrix W"""
+#        print 'wtf why are we using W'
+#        if getattr(self, '_W', None) is None:
+#            wlist = (self.Wsmall, self.Wx)
+#            if self.regmesh.dim > 1:
+#                wlist += (self.Wy,)
+#            if self.regmesh.dim > 2:
+#                wlist += (self.Wz,)
+#            self._W = sp.vstack(wlist)
+#        return self._W
+
+
+    @Utils.timeIt
+    def _evalSmall(self, m):
+        r = self.Wsmall * ( self.mapping * (m - self.mref) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmallDeriv(self, m):
+        r = self.Wsmall * ( self.mapping * (m - self.mref) )
+        return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
+
+    @Utils.timeIt
+    def _evalSmall2Deriv(self, m, v = None):
+        rDeriv = self.Wsmall * ( self.mapping.deriv(m - self.mref) )
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothx(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wx * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wx * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmoothy(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wy * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wy * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmoothz(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wz * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wz * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmooth(self, m):
+        phiSmooth = self._evalSmoothx(m)
+        if self.regmesh.dim > 1:
+            phiSmooth += self._evalSmoothy(m)
+        if self.regmesh.dim > 2:
+            phiSmooth += self._evalSmoothz(m)
+        return phiSmooth
+
+    @Utils.timeIt
+    def _evalSmoothxDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wx * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wx * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wx * ( self.mapping * m )
+            return r.T * ( self.Wx * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothx2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wx * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wx * ( self.mapping.deriv(m) )
+
+        if v is not None:
+            return rDeriv.T * ( rDeriv * v )
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothyDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wy * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wy * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wy * ( self.mapping * m )
+            return r.T * ( self.Wy * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothy2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wy * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wy * ( self.mapping.deriv(m) )
+
+        if v is not None:
+            return rDeriv.T * ( rDeriv * v )
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothzDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wz * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wz * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wz * ( self.mapping * m )
+            return r.T * ( self.Wz * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothz2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wz * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wz * ( self.mapping.deriv(m) )
+
+        if v is not None:
+            return rDeriv.T * ( rDeriv * v )
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothDeriv(self, m):
+        deriv = self._evalSmoothxDeriv(m)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothyDeriv(m)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothzDeriv(m)
+        return deriv
+
+    @Utils.timeIt
+    def _evalSmooth2Deriv(self, m, v=None):
+        deriv = self._evalSmoothx2Deriv(m, v)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothy2Deriv(m, v)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothz2Deriv(m, v)
+        return deriv
+
+
+    @Utils.timeIt
+    def eval(self, m):
+        return self._evalSmall(m) + self._evalSmooth(m)
+
+    @Utils.timeIt
+    def evalDeriv(self, m):
+        """
+        The regularization is:
+
+        .. math::
+
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+        So the derivative is straight forward:
+
+        .. math::
+
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+        """
+        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
+
+    @Utils.timeIt
+    def eval2Deriv(self, m, v=None):
+        return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v)
+
+
+
+class Tikhonov(Simple):
     """
     L2 Tikhonov regularization with both smallness and smoothness (first order
     derivative) contributions.
@@ -493,56 +726,131 @@ class Tikhonov(BaseRegularization):
             self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz
         return self._Wzz
 
+
     @property
-    def Wsmooth(self):
+    def Wsmooth2(self):
         """Full smoothness regularization matrix W"""
         if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx, self.Wxx)
+            wlist = (self.Wxx)
             if self.regmesh.dim > 1:
-                wlist += (self.Wy, self.Wyy)
+                wlist += (self.Wyy)
             if self.regmesh.dim > 2:
-                wlist += (self.Wz, self.Wzz)
+                wlist += (self.Wzz)
             self._Wsmooth = sp.vstack(wlist)
         return self._Wsmooth
 
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-
     @Utils.timeIt
-    def _evalSmall(self, m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmooth(self, m):
+    def _evalSmoothxx(self, m):
         if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
+            r = self.Wxx * ( self.mapping * (m - self.mref) )
         elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * (m) )
+            r = self.Wxx * ( self.mapping * (m) )
         return 0.5 * r.dot(r)
 
+    @Utils.timeIt
+    def _evalSmoothyy(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wyy * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wyy * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmoothzz(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wzz * ( self.mapping * (m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wzz * ( self.mapping * (m) )
+        return 0.5 * r.dot(r)
+
+    @Utils.timeIt
+    def _evalSmooth2(self, m):
+        phiSmooth2 = self._evalSmoothxx(m)
+        if self.regmesh.dim > 1:
+            phiSmooth2 += self._evalSmoothyy(m)
+        if self.regmesh.dim > 2:
+            phiSmooth2 += self._evalSmoothzz(m)
+        return phiSmooth2
+
+    @Utils.timeIt
+    def _evalSmoothxxDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wxx * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wxx * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wxx * ( self.mapping * m )
+            return r.T * ( self.Wxx * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothyyDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wyy * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wyy * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wyy * ( self.mapping * m )
+            return r.T * ( self.Wyy * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothzzDeriv(self, m):
+        if self.mrefInSmooth == True:
+            r = self.Wzz * ( self.mapping * ( m - self.mref ) )
+            return r.T * ( self.Wzz * self.mapping.deriv(m - self.mref) )
+        elif self.mrefInSmooth == False:
+            r = self.Wzz * ( self.mapping * m )
+            return r.T * ( self.Wzz * self.mapping.deriv(m) )
+
+    @Utils.timeIt
+    def _evalSmoothxx2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wxx * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wxx * self.mapping.deriv(m)
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothyy2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wyy * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wyy * self.mapping.deriv(m)
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothzz2Deriv(self, m, v=None):
+        if self.mrefInSmooth == True:
+            rDeriv = self.Wzz * ( self.mapping.deriv( m - self.mref ) )
+        elif self.mrefInSmooth == False:
+            rDeriv = self.Wzz * self.mapping.deriv(m)
+        if v is not None:
+            return rDeriv.T * (rDeriv * v)
+        return rDeriv.T * rDeriv
+
+    @Utils.timeIt
+    def _evalSmoothDeriv2(self, m):
+        deriv = self._evalSmoothxxDeriv(m)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothyyDeriv(m)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothzzDeriv(m)
+        return deriv
+
+    @Utils.timeIt
+    def _evalSmooth2Deriv2(self, m, v=None):
+        deriv = self._evalSmoothxx2Deriv(m, v)
+        if self.regmesh.dim > 1:
+            deriv += self._evalSmoothyy2Deriv(m, v)
+        if self.regmesh.dim > 2:
+            deriv += self._evalSmoothzz2Deriv(m, v)
+        return deriv
+
+
     @Utils.timeIt
     def eval(self, m):
-        return self._evalSmall(m) + self._evalSmooth(m)
-
-    @Utils.timeIt
-    def _evalSmallDeriv(self,m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
-
-    @Utils.timeIt
-    def _evalSmoothDeriv(self,m):
-        if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * ( m - self.mref ) )
-            return r.T * ( self.Wsmooth * self.mapping.deriv(m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * m )
-            return r.T * ( self.Wsmooth * self.mapping.deriv(m) )
+        return self._evalSmall(m) + self._evalSmooth(m) + self._evalSmooth2(m)
 
     @Utils.timeIt
     def evalDeriv(self, m):
@@ -560,184 +868,134 @@ class Tikhonov(BaseRegularization):
             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
 
         """
-        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
+        return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmoothDeriv2(m)
+
+    def eval2Deriv(self, m, v=None):
+        """
+        The regularization is:
+
+        .. math::
+
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+        So the derivative is straight forward:
+
+        .. math::
+
+            R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+        """
+        return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v) + self._evalSmooth2Deriv2(m, v)
 
 
-class Simple(Tikhonov):
+
+class Sparse(Simple):
     """
-    Simple regularization that does not include length scales in the derivatives.
+        The regularization is:
+    
+        .. math::
+    
+            R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top R^\\top R W(m-m_\\text{ref})}
+    
+        where the IRLS weight
+    
+        .. math::
+    
+            R = \eta TO FINISH LATER!!!
+    
+        So the derivative is straight forward:
+    
+        .. math::
+    
+            R(m) = \mathbf{W^\\top R^\\top R W (m-m_\\text{ref})}
+    
+        The IRLS weights are recomputed after each beta solves.
+        It is strongly recommended to do a few Gauss-Newton iterations
+        before updating.
     """
-
-    mrefInSmooth = False  #: SMOOTH and SMOOTH_MOD_DIF options
-    alpha_s     = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight")
-    alpha_x     = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
-    alpha_y     = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
-    alpha_z     = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
-    wght        = 1.
+        
+    # set default values
+    eps_p = 1e-1        # Threshold value for the model norm
+    eps_q = 1e-1        # Threshold value for the model gradient norm
+    curModel = None     # Requires model to compute the weights
+    l2model = None
+    gamma = 1.          # Model norm scaling to smooth out convergence
+    norms = [0., 2., 2., 2.] # Values for norm on (m, dmdx, dmdy, dmdz)
+    cell_weights = 1.        # Consider overwriting with sensitivity weights
 
     def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
+        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
 
-        if isinstance(self.wght,float):
-            self.wght = np.ones(self.regmesh.nC) * self.wght
+        if isinstance(self.cell_weights,float):
+            self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights
 
     @property
     def Wsmall(self):
         """Regularization matrix Wsmall"""
         if getattr(self,'_Wsmall', None) is None:
-            self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
+            if getattr(self, 'curModel', None) is None:
+                self.Rs = Utils.speye(self.regmesh.nC)
+
+            else:
+                f_m = self.mapping * (self.curModel - self.reg.mref)
+                self.rs = self.R(f_m , self.eps_p, self.norms[0])
+                self.Rs = Utils.sdiag( self.rs )
+
+            self._Wsmall = Utils.sdiag((self.alpha_s*self.gamma*self.cell_weights)**0.5)*self.Rs
+
         return self._Wsmall
 
     @property
     def Wx(self):
         """Regularization matrix Wx"""
-        if getattr(self, '_Wx', None) is None:
-            self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil
+        if getattr(self,'_Wx', None) is None:
+            if getattr(self, 'curModel', None) is None:
+                self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
+
+            else:
+                f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel)
+                self.rx = self.R( f_m , self.eps_q, self.norms[1])
+                self.Rx = Utils.sdiag( self.rx )
+
+            self._Wx = Utils.sdiag(( self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
+
         return self._Wx
 
     @property
     def Wy(self):
         """Regularization matrix Wy"""
-        if getattr(self, '_Wy', None) is None:
-            self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil
+        if getattr(self,'_Wy', None) is None:
+            if getattr(self, 'curModel', None) is None:
+                self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
+
+            else:
+                f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel)
+                self.ry = self.R( f_m , self.eps_q, self.norms[2])
+                self.Ry = Utils.sdiag( self.ry )
+
+            self._Wy = Utils.sdiag((self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
+
         return self._Wy
 
     @property
     def Wz(self):
         """Regularization matrix Wz"""
-        if getattr(self, '_Wz', None) is None:
-            self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil
+        if getattr(self,'_Wz', None) is None:
+            if getattr(self, 'curModel', None) is None:
+                self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
+
+            else:
+                f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel)
+                self.rz = self.R( f_m , self.eps_q, self.norms[3])
+                self.Rz = Utils.sdiag( self.rz )
+
+            self._Wz = Utils.sdiag((self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
+
         return self._Wz
 
-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        if getattr(self, '_Wsmooth', None) is None:
-            wlist = (self.Wx,)
-            if self.regmesh.dim > 1:
-                wlist += (self.Wy,)
-            if self.regmesh.dim > 2:
-                wlist += (self.Wz,)
-            self._Wsmooth = sp.vstack(wlist)
-        return self._Wsmooth
-
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-
-    @Utils.timeIt
-    def _evalSmall(self, m):
-        r = self.Wsmall * ( self.mapping * (m - self.mref) )
-        return 0.5 * r.dot(r)
-
-    @Utils.timeIt
-    def _evalSmooth(self, m):
-        if self.mrefInSmooth == True:
-            r = self.Wsmooth * ( self.mapping * (m - self.mref) )
-        elif self.mrefInSmooth == False:
-            r = self.Wsmooth * ( self.mapping * m)
-        return 0.5 * r.dot(r)
-
-
-class Sparse(Simple):
-
-    # set default values
-    eps_p      = 1e-1
-    eps_q      = 1e-1
-    curModel = None # use a model to compute the weights
-    gamma    = 1.
-    norms    = [0., 2., 2., 2.]
-    wght     = 1.
-
-    def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
-        Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
-
-        if isinstance(self.wght,float):
-            self.wght = np.ones(self.regmesh.nC) * self.wght
-
-    @property
-    def Wsmall(self):
-        """Regularization matrix Wsmall"""
-        if getattr(self, 'curModel', None) is None:
-            self.Rs = Utils.speye(self.regmesh.nC)
-
-        else:
-            f_m = self.curModel - self.reg.mref
-            self.rs = self.R(f_m , self.eps_p, self.norms[0])
-            #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
-            self.Rs = Utils.sdiag( self.rs )
-
-        return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
-
-
-    @property
-    def Wx(self):
-        """Regularization matrix Wx"""
-
-        if getattr(self, 'curModel', None) is None:
-            self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
-
-        else:
-            f_m = self.regmesh.cellDiffxStencil * self.curModel
-            self.rx = self.R( f_m , self.eps_q, self.norms[1])
-            self.Rx = Utils.sdiag( self.rx )
-
-        return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
-
-    @property
-    def Wy(self):
-        """Regularization matrix Wy"""
-
-        if getattr(self, 'curModel', None) is None:
-            self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0])
-
-        else:
-            f_m = self.regmesh.cellDiffyStencil * self.curModel
-            self.ry = self.R( f_m , self.eps_q, self.norms[2])
-            self.Ry = Utils.sdiag( self.ry )
-
-        return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
-
-    @property
-    def Wz(self):
-        """Regularization matrix Wz"""
-
-        if getattr(self, 'curModel', None) is None:
-            self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0])
-
-        else:
-            f_m = self.regmesh.cellDiffzStencil * self.curModel
-            self.rz = self.R( f_m , self.eps_q, self.norms[3])
-            self.Rz = Utils.sdiag( self.rz )
-
-        return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
-
-    @property
-    def Wsmooth(self):
-        """Full smoothness regularization matrix W"""
-        #if getattr(self, '_Wsmooth', None) is None:
-        wlist = (self.Wx,)
-        if self.regmesh.dim > 1:
-            wlist += (self.Wy,)
-        if self.regmesh.dim > 2:
-            wlist += (self.Wz,)
-        #self._Wsmooth = sp.vstack(wlist)
-        return sp.vstack(wlist)
-
-    @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Wsmall, self.Wsmooth)
-            self._W = sp.vstack(wlist)
-        return self._W
-
     def R(self, f_m , eps, exponent):
 
+        # Eta scaling is important for mix-norms...do not mess with it
         eta = (eps**(1.-exponent/2.))**0.5
         r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.)