From 2eb3d933a0db7f5ca522d84cfcf4364e91b72a57 Mon Sep 17 00:00:00 2001 From: Enzo Busseti Date: Sat, 14 Sep 2024 14:56:33 +0800 Subject: [PATCH] minor edits; planning to try a new formulation with QR nullspace projection now --- project_euromir/equilibrate.py | 2 +- project_euromir/solver_new.py | 73 +++++++++++++++++++++++++++------- 2 files changed, 60 insertions(+), 15 deletions(-) diff --git a/project_euromir/equilibrate.py b/project_euromir/equilibrate.py index 5e3ab65..1ac702d 100644 --- a/project_euromir/equilibrate.py +++ b/project_euromir/equilibrate.py @@ -34,7 +34,7 @@ def _cones_separation_matrix(zero, nonneg, second_order): def hsde_ruiz_equilibration( # pylint: disable=too-many-arguments matrix, b, c, dimensions, d=None, e=None, rho=1., sigma=1., - eps_rows=1E-4, eps_cols=1E-4, max_iters=25): + eps_rows=1E-1, eps_cols=1E-1, max_iters=25): """Ruiz equilibration of problem matrices for the HSDE system. :param matrix: Problem matrix. diff --git a/project_euromir/solver_new.py b/project_euromir/solver_new.py index 1610819..994dc17 100644 --- a/project_euromir/solver_new.py +++ b/project_euromir/solver_new.py @@ -26,7 +26,7 @@ from project_euromir.direction_calculator import ( CGNewton, DenseNewton, DiagPreconditionedCGNewton, ExactDiagPreconditionedCGNewton, LSMRLevenbergMarquardt, - LSQRLevenbergMarquardt, MinResQLPTest, WarmStartedCGNewton, + LSQRLevenbergMarquardt, MinResQLPTest, WarmStartedCGNewton, _densify, nocedal_wright_termination) from project_euromir.line_searcher import (BacktrackingLineSearcher, LogSpaceLineSearcher, @@ -42,7 +42,7 @@ logger = logging.getLogger(__name__) -QR_PRESOLVE = False +QR_PRESOLVE = True def solve(matrix, b, c, zero, nonneg, soc=(), # xy = None, # need to import logic for equilibration @@ -67,15 +67,19 @@ def solve(matrix, b, c, zero, nonneg, soc=(), equilibrate.hsde_ruiz_equilibration( matrix, b, c, dimensions={ 'zero': zero, 'nonneg': nonneg, 'second_order': soc}, - max_iters=25) + max_iters=1000)#, eps_rows=1E-5, eps_cols=1E-5,) if QR_PRESOLVE: - q, r = np.linalg.qr( - np.vstack([matrix_transf.todense(), c_transf.reshape((1, n))])) - matrix_transf = q[:-1].A - c_transf = q[-1].A1 + # q, r = np.linalg.qr( + # np.vstack([matrix_transf.todense(), c_transf.reshape((1, n))])) + # matrix_transf = q[:-1].A + # c_transf = q[-1].A1 + q, r = np.linalg.qr(matrix_transf.todense()) + matrix_transf = q.A + c_transf = np.linalg.solve(r, c_transf) + sigma_qr = np.linalg.norm( - b_transf) / np.mean(np.linalg.norm(matrix_transf, axis=1)) + b_transf) #/ np.mean(np.linalg.norm(matrix_transf, axis=1)) b_transf = b_transf/sigma_qr workspace = create_workspace(m, n, zero) @@ -100,9 +104,9 @@ def _local_hessian_x_nogap(x): return hessian_x_nogap( x, m=m, n=n, zero=zero, matrix=matrix_transf, b=b_transf) - def _local_hessian_y_nogap(x): + def _local_hessian_y_nogap(y): return hessian_y_nogap( - x, m=m, n=n, zero=zero, matrix=matrix_transf) + y, m=m, n=n, zero=zero, matrix=matrix_transf) def _local_residual(xy): return residual( @@ -137,7 +141,7 @@ def _local_derivative_residual(xy): # direction_calculator = WarmStartedCGNewton( # # warm start causes issues if null space changes b/w iterations # hessian_function=_local_hessian, - # rtol_termination=lambda x, g: min(0.5, np.linalg.norm(g)**0.5), + # rtol_termination=lambFalseda x, g: min(0.5, np.linalg.norm(g)**0.5), # max_cg_iters=None, # minres=False, # regularizer=1e-8, # it seems 1e-10 is best, but it's too sensitive to it :( @@ -186,7 +190,7 @@ def _local_derivative_residual(xy): # direction_calculator = LSMRLevenbergMarquardt( # residual_function=_local_residual, # derivative_residual_function=_local_derivative_residual, - # warm_start=True, # also doesn't work with warm start + # warm_start=True, # also doesn't work with warm start # ) # direction_calculator = DenseNewton( #WarmStartedCGNewton( @@ -201,6 +205,8 @@ def _local_derivative_residual(xy): _start = time.time() # extra_iters=5 # all_losses = [] + # all_dirnorms = [] + # all_dirnorms_times_steplen = [] for newton_iterations in range(1000): @@ -222,16 +228,51 @@ def _local_derivative_residual(xy): # logger.info('Converged in %d iterations.', newton_iterations) # break + # dense_hess = _densify(_local_hessian(xy)) + # dense_hessx_nogap = _densify(_local_hessian_x_nogap(xy[:n])) + # dense_hessy_nogap = _densify(_local_hessian_y_nogap(xy[n:])) + # eivals = np.linalg.eigh(dense_hess)[0] + + # #diag_precond = np.diag(1./(np.diag(dense_hess))) + # #dense_hess_diag_precond = dense_hess @ diag_precond + # #eivals_diag_precond = np.linalg.eigh(dense_hess_diag_precond)[0] + + # import matplotlib.pyplot as plt + # plt.plot(eivals, label='hess eivals') + # # plt.plot(eivals_diag_precond) + # plt.plot(np.diag(dense_hess), label='hess diag') + # diag_nogap = np.concatenate([np.diag(dense_hessx_nogap),np.diag(dense_hessy_nogap)]) + # plt.plot(diag_nogap, label='hess diag nogap') + + # gap = np.concatenate([c_transf, b_transf]) + # approx_hess = np.diag(diag_nogap) + np.outer(gap, gap) + np.eye(n+m) * 1e-2 + # eivals_approx = np.linalg.eigh(approx_hess)[0] + # plt.plot(eivals_approx, label='eivals approx hess') + + # pinv_approx_hess = np.linalg.pinv(approx_hess) + # precond_hess = dense_hess @ pinv_approx_hess + # eivals_precond = np.linalg.eigh(precond_hess)[0] + # plt.plot(eivals_precond, label='eivals pinv precond hess') + + # plt.legend() + # plt.show() + # breakpoint() + direction = direction_calculator.get_direction( current_point=xy, current_gradient=grad_xy) + logger.info('direction norm %.2e', np.linalg.norm(direction)) + # all_dirnorms.append(np.linalg.norm(direction)) + # oldxy = np.copy(xy) + # all_losses.append(loss_xy) + xy, loss_xy, grad_xy = \ line_searcher.get_next(current_point=xy, current_loss=loss_xy, current_gradient=grad_xy, direction=direction) - # all_losses.append(loss_xy) + # all_dirnorms_times_steplen.append(np.linalg.norm(xy-oldxy)) # import matplotlib.pyplot as plt # iter_x = xy[:n] @@ -268,8 +309,12 @@ def _local_derivative_residual(xy): f'Solver did not converge in {newton_iterations} iterations.') # import matplotlib.pyplot as plt - # plt.semilogy(all_losses) + # plt.semilogy(all_dirnorms) + # plt.semilogy(all_dirnorms_times_steplen) + # plt.semilogy(np.sqrt(all_losses)) + # plt.show() + # breakpoint() if loss_xy > np.finfo(float).eps: raise NotImplementedError(