FIX compute p_boj_acc only if acceleration is performed (#41)

Author: Badr-MOUFAD

skglm/solvers/group_bcd_solver.py

@@ -91,13 +91,15 @@ def bcd_solver(X, y, datafit, penalty, w_init=None, p0=10,
             # inplace update of w and Xw
             _bcd_epoch(X, y, w, Xw, datafit, penalty, ws)
 
-            w_acc, Xw_acc = accelerator.extrapolate(w, Xw)
-            p_obj = datafit.value(y, w, Xw) + penalty.value(w)
-            p_obj_acc = datafit.value(y, w_acc, Xw_acc) + penalty.value(w_acc)
+            w_acc, Xw_acc, is_extrapolated = accelerator.extrapolate(w, Xw)
+
+            if is_extrapolated:  # avoid computing p_obj for un-extrapolated w, Xw
+                p_obj = datafit.value(y, w, Xw) + penalty.value(w)
+                p_obj_acc = datafit.value(y, w_acc, Xw_acc) + penalty.value(w_acc)
 
-            if p_obj_acc < p_obj:
-                w, Xw = w_acc, Xw_acc
-                p_obj = p_obj_acc
+                if p_obj_acc < p_obj:
+                    w[:], Xw[:] = w_acc, Xw_acc
+                    p_obj = p_obj_acc
 
             # check sub-optimality every 10 epochs
             if epoch % 10 == 0:
@@ -106,13 +108,15 @@ def bcd_solver(X, y, datafit, penalty, w_init=None, p0=10,
                 stop_crit_in = np.max(opt_in)
 
                 if max(verbose - 1, 0):
+                    p_obj = datafit.value(y, w, Xw) + penalty.value(w)
                     print(
                         f"Epoch {epoch + 1}, objective {p_obj:.10f}, "
                         f"stopping crit {stop_crit_in:.2e}"
                     )
 
                 if stop_crit_in <= 0.3 * stop_crit:
                     break
+        p_obj = datafit.value(y, w, Xw) + penalty.value(w)
         p_objs_out[t] = p_obj
 
     return w, p_objs_out, stop_crit

skglm/tests/test_group.py

@@ -137,7 +137,7 @@ def test_anderson_acceleration():
     w = np.ones(n_features)
     Xw = X @ w
     for i in range(max_iter):
-        w, Xw = acc.extrapolate(w, Xw)
+        w, Xw, _ = acc.extrapolate(w, Xw)
         w = rho * w + 1
         Xw = X @ w
 

skglm/utils.py

@@ -273,7 +273,7 @@ def __init__(self, K):
         self.arr_w_, self.arr_Xw_ = None, None
 
     def extrapolate(self, w, Xw):
-        """Return ``w`` and ``Xw`` extrapolated."""
+        """Return w, Xw, and a bool indicating whether they were extrapolated."""
         if self.arr_w_ is None or self.arr_Xw_ is None:
             self.arr_w_ = np.zeros((w.shape[0], self.K+1))
             self.arr_Xw_ = np.zeros((Xw.shape[0], self.K+1))
@@ -282,19 +282,19 @@ def extrapolate(self, w, Xw):
             self.arr_w_[:, self.current_iter] = w
             self.arr_Xw_[:, self.current_iter] = Xw
             self.current_iter += 1
-            return w, Xw
+            return w, Xw, False
 
         U = np.diff(self.arr_w_, axis=1)  # compute residuals
 
         # compute extrapolation coefs
         try:
             inv_UTU_ones = np.linalg.solve(U.T @ U, np.ones(self.K))
         except np.linalg.LinAlgError:
-            return w, Xw
+            return w, Xw, False
         finally:
             self.current_iter = 0
 
         # extrapolate
         C = inv_UTU_ones / np.sum(inv_UTU_ones)
         # floating point errors may cause w and Xw to disagree
-        return self.arr_w_[:, 1:] @ C, self.arr_Xw_[:, 1:] @ C
+        return self.arr_w_[:, 1:] @ C, self.arr_Xw_[:, 1:] @ C, True