diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py
index c8dd9caf..68998a3f 100644
--- a/SimPEG/Directives.py
+++ b/SimPEG/Directives.py
@@ -122,8 +122,8 @@ class BetaEstimate_ByEig(InversionDirective):
 
         if self.debug: print 'Calculating the beta0 parameter.'
 
-        m = self.invProb.m_current
-        u = self.invProb.u_current or self.prob.fields(m)
+        m = self.invProb.curModel
+        u = self.prob.fields(m)
 
         x0 = np.random.rand(*m.shape)
         t = x0.dot(self.dmisfit.eval2Deriv(m,x0,u=u))
diff --git a/SimPEG/InvProblem.py b/SimPEG/InvProblem.py
index 259e7e5b..d75362c2 100644
--- a/SimPEG/InvProblem.py
+++ b/SimPEG/InvProblem.py
@@ -18,9 +18,22 @@ class BaseInvProblem(object):
     reg     = None   #: Regularization
     opt     = None   #: Optimization program
 
-    u_current = None #: The most current evaluated field
-    m_current = None #: The most current model
+    deleteTheseOnModelUpdate = [] # List of strings, e.g. ['_MeSigma', '_MeSigmaI']
 
+    @property
+    def curModel(self):
+        """
+            Sets the current model, and removes dependent properties
+        """
+        return getattr(self, '_curModel', None)
+    @curModel.setter
+    def curModel(self, value):
+        if value is self.curModel:
+            return # it is the same!
+        self._curModel = value
+        for prop in self.deleteTheseOnModelUpdate:
+            if hasattr(self, prop):
+                delattr(self, prop)
 
     def __init__(self, dmisfit, reg, opt, **kwargs):
         Utils.setKwargs(self, **kwargs)
@@ -48,23 +61,21 @@ class BaseInvProblem(object):
         self.phi_d = np.nan
         self.phi_m = np.nan
 
-        self.m_current = m0
+        self.curModel = m0
 
         print 'SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv. \n    ***Done using direct methods***'
-        self.opt.bfgsH0 = Solver(self.reg.eval2Deriv(self.m_current))
+        self.opt.bfgsH0 = Solver(self.reg.eval2Deriv(self.curModel))
 
     @Utils.timeIt
     def evalFunction(self, m, return_g=True, return_H=True):
         """evalFunction(m, return_g=True, return_H=True)
         """
 
-        self.u_current = None
-        self.m_current = m
-        forward = self.prob
+        #TODO: check for warmstart
+        self.curModel = m
         gc.collect()
 
         u = self.prob.fields(m)
-        self.u_current = u
 
         phi_d = self.dmisfit.eval(m, u=u)
         phi_m = self.reg.eval(m)
diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py
index 691f143e..3fe9cfd0 100644
--- a/SimPEG/Inversion.py
+++ b/SimPEG/Inversion.py
@@ -59,7 +59,7 @@ class BaseInversion(object):
         """
         self.invProb.startup(m0)
         self.directiveList.call('initialize')
-        self.m = self.opt.minimize(self.invProb.evalFunction, m0)
+        self.m = self.opt.minimize(self.invProb.evalFunction, self.invProb.curModel)
         self.directiveList.call('finish')
 
         return self.m
diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py
index fc7941d8..a2d65e28 100644
--- a/SimPEG/Optimization.py
+++ b/SimPEG/Optimization.py
@@ -687,7 +687,7 @@ class BFGS(Minimize, Remember):
     def bfgsrec(self,k,n,nn,S,Y,d):
         """BFGS recursion"""
         if k < 0:
-            d = self.bfgsH0 * (d)
+            d = self.bfgsH0 * d
         else:
             khat    = 0 if nn is 0 else np.mod(n-nn+k,nn)
             gamma   = np.vdot(S[:,khat],d)/np.vdot(Y[:,khat],S[:,khat])