refinement without hsde plugged in;

makes very little difference, just simplifies design,
for infeas/unbound will have to see

project_euromir/refinement_no_hsde.py

@@ -16,6 +16,8 @@
 # Project Euromir. If not, see <https://www.gnu.org/licenses/>.
 """Define residual and Dresidual for use by refinement loop."""
 
+import time
+
 import numpy as np
 import scipy as sp
 
@@ -128,6 +130,46 @@ def _rmatvec(dr):
         matvec = _matvec,
         rmatvec = _rmatvec)
 
+def refine(xz, m, n, zero, nonneg, matrix, b, c, soc=(), max_iters=None):
+    """Refinement without using HSDE."""
+
+    res = residual(
+        xz, m=m, n=n, zero=zero, nonneg=nonneg, matrix=matrix, b=b, c=c,
+        soc=soc)
+
+    oldnormsq = np.linalg.norm(res)**2
+    print('residual norm sq before LSQR', oldnormsq)
+
+    DR = Dresidual(
+        xz, m=m, n=n, zero=zero, nonneg=nonneg, matrix=matrix, b=b, c=c,
+        soc=soc)
+
+    # call LSQR
+    start = time.time()
+    result = sp.sparse.linalg.lsqr(
+        DR, res, atol=0., btol=0.,
+        iter_lim=(n+m)*2 if max_iters is None else max_iters)
+    print('LSQR result[1:-1]', result[1:-1])
+    print('LSQR took', time.time() - start)
+    dxz = result[0]
+
+    # recompute problem variables
+    for i in range(20):
+        xz1 = xz - dxz * (0.5) ** i
+        res = residual(
+            xz1, m=m, n=n, zero=zero, nonneg=nonneg, matrix=matrix, b=b, c=c,
+            soc=soc)
+        newnormsq = np.linalg.norm(res)**2
+
+        if newnormsq < oldnormsq:
+            break
+    else:
+        print('REFINEMENT FAILED!')
+    print("residual norm sq after LSQR", newnormsq)
+
+    return xz1
+
+
 if __name__ == '__main__': # pragma: no cover
 
     from scipy.optimize import check_grad
@@ -170,3 +212,10 @@ def my_Dresidual_dense(xz):
         DRT = _densify_also_nonsquare(my_Dresidual(xy).T)
         assert np.allclose(DR.T, DRT)
     print('\tOK!')
+
+    print('\nCHECKING REFINE')
+    for i in range(1):
+        xz = np.random.randn(n+m)
+        for i in range(2):
+            xz = refine(
+                xz, m, n, zero, nonneg, matrix, b, c, soc=(), max_iters=None)

project_euromir/solver_new.py

@@ -37,7 +37,10 @@
                                           loss_gradient, residual)
 from project_euromir.minamide import (MinamideTest, hessian_x_nogap,
                                       hessian_y_nogap)
-from project_euromir.refinement import refine
+from project_euromir.refinement import refine as hsde_refine
+from project_euromir.refinement_no_hsde import refine
+
+HSDE_REFINEMENT = False
 
 logger = logging.getLogger(__name__)
 
@@ -281,33 +284,52 @@ def _local_derivative_residual(xy):
     print('DirectionCalculator statistics', direction_calculator.statistics)
     print('LineSearcher statistics', line_searcher.statistics)
 
-    # create HSDE variables for refinement
-    u = np.empty(n+m+1, dtype=float)
-    u[:-1] = xy
-    u[-1] = 1.
-    v = np.zeros_like(u)
-    v[n:-1] = -matrix_transf @ u[:n] + b_transf
-
-    for _ in range(3):
-        u, v = refine(
-            z=u-v, matrix=matrix_transf, b=b_transf, c=c_transf, zero=zero,
-            nonneg=nonneg, soc=soc)
-
-    if u[-1] < 1e-8:
-        raise FloatingPointError(
-            "Refinement failed, Newton-CG solution wasn't good enough.")
-
-    # Transform back into problem format
-    u1, u2, u3 = u[:n], u[n:n+m], u[-1]
-    v2, v3 = v[n:n+m], v[-1]
+    if not HSDE_REFINEMENT:
+        # switch to refinement
+        x = xy[:n]
+        y = xy[n:]
+        y[zero:] = np.maximum(y[zero:], 0.)
+        s = -matrix_transf @ x + b_transf
+        s[:zero] = 0.
+        s[zero:] = np.maximum(s[zero:], 0.)
+        z = y-s
+        xz = np.concatenate([x, z])
+        for _ in range(3):
+            xz = refine(xz, matrix=matrix_transf, b=b_transf, c=c_transf, zero=zero,
+                nonneg=nonneg, soc=soc, m=m, n=n)
+        x = xz[:n]
+        y = xz[n:]
+        y[zero:] = np.maximum(y[zero:], 0.)
 
-    if v3 > u3:
-        raise NotImplementedError('Certificates not yet implemented.')
+    else:
 
-    # Apply HSDE scaling
-    x = u1 / u3
-    y = u2 / u3
-    s = v2 / u3
+        # create HSDE variables for refinement
+        u = np.empty(n+m+1, dtype=float)
+        u[:-1] = xy
+        u[-1] = 1.
+        v = np.zeros_like(u)
+        v[n:-1] = -matrix_transf @ u[:n] + b_transf
+
+        for _ in range(3):
+            u, v = hsde_refine(
+                z=u-v, matrix=matrix_transf, b=b_transf, c=c_transf, zero=zero,
+                nonneg=nonneg, soc=soc)
+
+        if u[-1] < 1e-8:
+            raise FloatingPointError(
+                "Refinement failed, Newton-CG solution wasn't good enough.")
+
+        # Transform back into problem format
+        u1, u2, u3 = u[:n], u[n:n+m], u[-1]
+        v2, v3 = v[n:n+m], v[-1]
+
+        if v3 > u3:
+            raise NotImplementedError('Certificates not yet implemented.')
+
+        # Apply HSDE scaling
+        x = u1 / u3
+        y = u2 / u3
+        s = v2 / u3
 
     if QR_PRESOLVE:
         x = np.linalg.solve(r, x) * sigma_qr