diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py
index 64445ab2..297270ef 100644
--- a/SimPEG/Inversion.py
+++ b/SimPEG/Inversion.py
@@ -32,12 +32,6 @@ class BaseInversion(object):
             self.opt.printers.insert(2,IterationPrinters.phi_d)
             self.opt.printers.insert(3,IterationPrinters.phi_m)
 
-        if not hasattr(opt, '_bfgsH0') and hasattr(opt, 'bfgsH0'): # Check if it has been set by the user and the default is not being used.
-            #TODO: I don't think that this if statement is working...
-            print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.'
-            opt.bfgsH0 = SimPEG.Solver(objFunc.reg.modelObj2Deriv())
-
-
     #TODO: Move this to the data class?
     @property
     def phi_d_target(self):
diff --git a/SimPEG/ObjFunction.py b/SimPEG/ObjFunction.py
index 47354b4b..2b1912b8 100644
--- a/SimPEG/ObjFunction.py
+++ b/SimPEG/ObjFunction.py
@@ -97,7 +97,7 @@ class BaseObjFunction(object):
         if return_H:
             def H_fun(v):
                 phi_d2Deriv = self.dataObj2Deriv(m, v, u=u)
-                phi_m2Deriv = self.reg.modelObj2Deriv()*v
+                phi_m2Deriv = self.reg.modelObj2Deriv(m, v=v)
 
                 return phi_d2Deriv + self.beta * phi_m2Deriv
 
diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py
index 8e08250b..8bb5c95c 100644
--- a/SimPEG/Optimization.py
+++ b/SimPEG/Optimization.py
@@ -648,10 +648,16 @@ class BFGS(Minimize, Remember):
 
             Must be a SimPEG.Solver
         """
-        _bfgsH0 = getattr(self,'_bfgsH0',None)
-        if _bfgsH0 is None:
-            return Solver(sp.identity(self.xc.size).tocsc(), flag='D')
-        return _bfgsH0
+        if getattr(self,'_bfgsH0',None) is None:
+            # Check if it has been set by the user and the default is not being used.
+            if self.parent is None:
+                self._bfgsH0 = Solver(sp.identity(self.xc.size).tocsc(), flag='D')
+            else:
+                print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.'
+                objFunc = self.parent.objFunc
+                self._bfgsH0 = Solver(objFunc.reg.modelObj2Deriv(objFunc.m_current))
+        return self._bfgsH0
+
     @bfgsH0.setter
     def bfgsH0(self, value):
         assert type(value) is Solver, 'bfgsH0 must be a SimPEG.Solver'
diff --git a/SimPEG/Parameters.py b/SimPEG/Parameters.py
index 9b12d8d1..5a6760e6 100644
--- a/SimPEG/Parameters.py
+++ b/SimPEG/Parameters.py
@@ -141,7 +141,7 @@ class BetaEstimate(Parameter):
 
         x0 = np.random.rand(*m.shape)
         t = x0.dot(objFunc.dataObj2Deriv(m,x0,u=u))
-        b = x0.dot(objFunc.reg.modelObj2Deriv()*x0)
+        b = x0.dot(objFunc.reg.modelObj2Deriv(m, v=x0))
         return self.beta0_ratio*(t/b)
 
 
diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py
index 4fbf5d4b..8dc32b79 100644
--- a/SimPEG/Regularization.py
+++ b/SimPEG/Regularization.py
@@ -79,12 +79,18 @@ class BaseRegularization(object):
             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
 
         """
-        return self.W.T * ( self.W * self.model.transform(m - self.mref) )
+        mTd = self.model.transformDeriv(m - self.mref)
+        return mTd.T * ( self.W.T * ( self.W * self.model.transform(m - self.mref) ) )
 
     @Utils.timeIt
-    def modelObj2Deriv(self):
+    def modelObj2Deriv(self, m, v=None):
         """
 
+            :param numpy.array m: geophysical model
+            :param numpy.array v: vector to multiply
+            :rtype: scipy.sparse.csr_matrix or numpy.ndarray
+            :return: WtW or WtW*v
+
         The regularization is:
 
         .. math::
@@ -98,7 +104,12 @@ class BaseRegularization(object):
             R(m) = \mathbf{W^\\top W}
 
         """
-        return self.W.T * self.W
+        mTd = self.model.transformDeriv(m - self.mref)
+        if v is None:
+            return mTd.T * self.W.T * self.W * mTd
+
+        return mTd.T * ( self.W.T * ( self.W * ( mTd * v) ) )
+