new model, using newton-CG;

about 3x faster than before already

project_euromir/loss_nullspace_proj.py

@@ -1,3 +1,19 @@
+# Copyright 2024 Enzo Busseti
+#
+# This file is part of Project Euromir.
+#
+# Project Euromir is free software: you can redistribute it and/or modify it
+# under the terms of the GNU General Public License as published by the Free
+# Software Foundation, either version 3 of the License, or (at your option) any
+# later version.
+#
+# Project Euromir is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+# details.
+#
+# You should have received a copy of the GNU General Public License along with
+# Project Euromir. If not, see <https://www.gnu.org/licenses/>.
 """Loss and related functions for nullspace projection model."""
 
 import numpy as np

project_euromir/solver_new.py

@@ -33,6 +33,7 @@
                                            ScipyLineSearcher)
 from project_euromir.loss_no_hsde import (Dresidual, create_workspace, hessian,
                                           loss_gradient, residual)
+from project_euromir.loss_nullspace_proj import NullSpaceModel
 from project_euromir.minamide import (MinamideTest, hessian_x_nogap,
                                       hessian_y_nogap)
 from project_euromir.refinement import refine as hsde_refine
@@ -43,6 +44,7 @@
 logger = logging.getLogger(__name__)
 
 QR_PRESOLVE = True
+NULLSPACE = True
 
 def solve(matrix, b, c, zero, nonneg, soc=(),
         # xy = None, # need to import logic for equilibration
@@ -88,10 +90,10 @@ def solve(matrix, b, c, zero, nonneg, soc=(),
     # nullspace projection idea
     y0 = -c_transf @ matrix_transf.T
     # y is y0 plus vector in span of nullspace_projector
-    assert np.allclose(matrix_transf.T @ (y0 + nullspace_projector @ np.random.randn(m-n)) + c_transf, 0.)
+    assert np.allclose(
+        matrix_transf.T @ (
+            y0 + nullspace_projector @ np.random.randn(m-n)) + c_transf, 0.)
 
-    variable = np.zeros(m)
-    #
     def new_nullspace_loss(variable):
         x = variable[:n]; y_null = variable[n:]
         y = y0 + nullspace_projector @ y_null
@@ -103,11 +105,16 @@ def new_nullspace_loss(variable):
         return (y_loss + s_loss_zero + s_loss_nonneg + gap_loss) / 2.
 
     # import scipy.optimize as opt
+    # variable = np.zeros(m)
     # result = opt.minimize(new_nullspace_loss, variable, tol=1e-32)
     # print(result)
 
     # breakpoint()
 
+    ns_model = NullSpaceModel(
+        m=m, n=n, zero=zero, nonneg=nonneg, matrix_transfqr=matrix_transf,
+        b=b_transf, c=c_transf, nullspace_projector=nullspace_projector)
+
     workspace = create_workspace(m, n, zero)
 
     def _local_loss(xy):
@@ -144,9 +151,19 @@ def _local_derivative_residual(xy):
             xy, m=m, n=n, zero=zero, matrix=matrix_transf, b=b_transf,
             c=c_transf, nonneg=nonneg,  soc=soc)
 
-    xy = np.zeros(n+m)
-    loss_xy = _local_loss(xy)
-    grad_xy = _local_grad(xy)
+    if NULLSPACE:
+        xy = np.zeros(m)
+        loss_func = ns_model.loss
+        grad_func = ns_model.gradient
+        hess_func = ns_model.hessian
+    else:
+        xy = np.zeros(n+m)
+        loss_func = _local_loss
+        grad_func = _local_grad
+        hess_func = _local_hessian
+
+    loss_xy = loss_func(xy)
+    grad_xy = grad_func(xy)
 
     # line_searcher = LogSpaceLineSearcher(
     #     loss_function=_local_loss,
@@ -161,7 +178,7 @@ def _local_derivative_residual(xy):
     line_searcher = BacktrackingLineSearcher(
         # Scipy searcher is not stable enough, breaks on numerical errors
         # with small steps
-        loss_function=_local_loss,
+        loss_function=loss_func,
         max_iters=1000)
 
     # direction_calculator = WarmStartedCGNewton(
@@ -190,7 +207,7 @@ def _local_derivative_residual(xy):
     # direction_calculator = ExactDiagPreconditionedCGNewton(
     direction_calculator  = CGNewton(
         # warm start causes issues if null space changes b/w iterations
-        hessian_function=_local_hessian,
+        hessian_function=hess_func,
         rtol_termination=lambda x, g: min(0.5, np.linalg.norm(g)**0.5), #,**2),
         max_cg_iters=None,
         # minres=True, # less stable, needs more stringent termination,
@@ -236,7 +253,7 @@ def _local_derivative_residual(xy):
 
     for newton_iterations in range(1000):
 
-        grad_xy = _local_grad(xy)
+        grad_xy = grad_func(xy)
 
         logger.info(
             'Iteration %d, current loss %.2e, current inf norm grad %.2e',
@@ -354,6 +371,12 @@ def _local_derivative_residual(xy):
     print('DirectionCalculator statistics', direction_calculator.statistics)
     print('LineSearcher statistics', line_searcher.statistics)
 
+    if NULLSPACE:
+        x, y_null = xy[:n], xy[n:]
+        xy = np.empty(n+m)
+        xy[:n] = x
+        xy[n:] = y0 + nullspace_projector @ y_null
+
     if not HSDE_REFINEMENT:
         # switch to refinement
         x = xy[:n]