From 7149b9379b8ec50ee093fa01e91a3f4bed741141 Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Mon, 1 Sep 2025 20:13:41 +0200 Subject: [PATCH 1/6] Add success flag to ipopt output --- cert_tools/sdp_solvers.py | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/cert_tools/sdp_solvers.py b/cert_tools/sdp_solvers.py index c04a352..707af64 100644 --- a/cert_tools/sdp_solvers.py +++ b/cert_tools/sdp_solvers.py @@ -164,19 +164,27 @@ def solve_low_rank_sdp( S = cas.nlpsol("S", "ipopt", nlp, options) else: S = cas.nlpsol("S", method, nlp, options) - + # Run Program sol_input = dict(lbg=g_rhs, ubg=g_rhs) if x_cand is not None: sol_input["x0"] = x_cand.reshape((-1, 1)) r = S(**sol_input) - Y_opt = r["x"] + + # Get solver info + stats = S.stats() + success = stats.get("success", False) + info = {"solver_stats": stats, "success": success} + # Reshape and generate SDP solution + Y_opt = r["x"] Y_opt = np.array(Y_opt).reshape((n, rank), order="F") X_opt = Y_opt @ Y_opt.T + # Get cost scale, offset = adjust cost = np.trace(Q @ X_opt) * scale + offset + # Construct certificate mults = np.array(r["lam_g"]) H = Q @@ -186,10 +194,11 @@ def solve_low_rank_sdp( H = H + A * mults[i, 0] # Return - info = {"X": X_opt, "H": H, "cost": cost} + info = {"X": X_opt, "H": H, "cost": cost, "success": success} return Y_opt, info + def solve_sdp_mosek( Q, Constraints, From bd5b4e206d2fbed8270a3653116ce02018c027d7 Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Mon, 8 Sep 2025 13:02:03 +0200 Subject: [PATCH 2/6] Implement low-rank projection algorithm --- cert_tools/linalg_tools.py | 41 ++++++++++++++++++++++++++++++++++++-- cert_tools/sdp_solvers.py | 3 +-- 2 files changed, 40 insertions(+), 4 deletions(-) diff --git a/cert_tools/linalg_tools.py b/cert_tools/linalg_tools.py index 6d93bb1..af65ed1 100644 --- a/cert_tools/linalg_tools.py +++ b/cert_tools/linalg_tools.py @@ -1,3 +1,4 @@ +import warnings from copy import deepcopy import numpy as np @@ -84,7 +85,7 @@ def project_so3(X): return U @ Vh -def rank_project(X, p=1, tolerance=1e-10): +def rank_project(X, p=1, tolerance=1e-10) -> tuple[np.ndarray, dict]: """Project matrix X to matrix of rank p.""" try: assert la.issymmetric(X, atol=tolerance) @@ -113,7 +114,43 @@ def rank_project(X, p=1, tolerance=1e-10): "error eigs": np.sum(np.abs(E[p:])), "EVR": abs(E[p - 1] / E[p]), # largest over second-largest } - return x, info + return np.asarray(x), info + + +def extract_lower_rank_solution(X, Constraints, tol=1e-8, max_trials=10): + """Given a psd matrix X of rank r, returns a matrix X' of rank r-1""" + r_start = 0 + for k in range(max_trials): + E, U = np.linalg.eigh(X) + r = int(np.sum(E > tol)) + if r == 1: + return X, {"success": True, "rank": 1, "iter": k} + elif r == r_start - 1: + return X, {"success": True, "rank": r, "iter": k} + if r_start == 0: + r_start = r + V = U[:, -r:] @ np.diag(np.sqrt(E[-r:])) + np.testing.assert_allclose(X, V @ V.T, atol=tol) + + AV = np.vstack([svec(V.T @ Ai @ V) for Ai, _ in Constraints]) # type: ignore + + # below, the rows of bi have the nullspace vectors. + # AV @ bi = 0 + b, info_null = get_nullspace(AV, tolerance=tol) + + # chose one of the vectors + idx = np.random.choice(b.shape[0]) + bi = b[idx] + + assert np.all(AV @ bi < tol) + Delta = smat(bi) + + lambdas = np.abs(np.linalg.eigvalsh(Delta)) + # maximum magnitude eigenvalue + alpha = -1 / np.max(lambdas) + X = V @ (np.eye(r) + alpha * Delta) @ V.T + warnings.warn(f"Did not find a lower-rank solution in {max_trials} trials.") + return X, {"success": True, "rank": r, "iter": k} def find_dependent_columns(A_sparse, tolerance=1e-10, verbose=False, debug=False): diff --git a/cert_tools/sdp_solvers.py b/cert_tools/sdp_solvers.py index 707af64..8fecfef 100644 --- a/cert_tools/sdp_solvers.py +++ b/cert_tools/sdp_solvers.py @@ -164,7 +164,7 @@ def solve_low_rank_sdp( S = cas.nlpsol("S", "ipopt", nlp, options) else: S = cas.nlpsol("S", method, nlp, options) - + # Run Program sol_input = dict(lbg=g_rhs, ubg=g_rhs) if x_cand is not None: @@ -198,7 +198,6 @@ def solve_low_rank_sdp( return Y_opt, info - def solve_sdp_mosek( Q, Constraints, From e438415d48b18360c4b854646e3cb18151b3e922 Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Tue, 9 Sep 2025 09:03:47 +0200 Subject: [PATCH 3/6] Change Constraints argument to A_list (b not used) --- cert_tools/linalg_tools.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/cert_tools/linalg_tools.py b/cert_tools/linalg_tools.py index af65ed1..396b1de 100644 --- a/cert_tools/linalg_tools.py +++ b/cert_tools/linalg_tools.py @@ -117,7 +117,7 @@ def rank_project(X, p=1, tolerance=1e-10) -> tuple[np.ndarray, dict]: return np.asarray(x), info -def extract_lower_rank_solution(X, Constraints, tol=1e-8, max_trials=10): +def extract_lower_rank_solution(X, A_list, tol=1e-8, max_trials=10): """Given a psd matrix X of rank r, returns a matrix X' of rank r-1""" r_start = 0 for k in range(max_trials): @@ -132,11 +132,14 @@ def extract_lower_rank_solution(X, Constraints, tol=1e-8, max_trials=10): V = U[:, -r:] @ np.diag(np.sqrt(E[-r:])) np.testing.assert_allclose(X, V @ V.T, atol=tol) - AV = np.vstack([svec(V.T @ Ai @ V) for Ai, _ in Constraints]) # type: ignore + AV = np.vstack([svec(V.T @ Ai @ V) for Ai in A_list]) # type: ignore # below, the rows of bi have the nullspace vectors. # AV @ bi = 0 b, info_null = get_nullspace(AV, tolerance=tol) + if b.shape[0] == 0: + return X, {"success": True, "rank": r, "iter": k} + break # chose one of the vectors idx = np.random.choice(b.shape[0]) From 326339cda0bf86f4504f11cc682dc25719281c10 Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Tue, 9 Sep 2025 09:11:02 +0200 Subject: [PATCH 4/6] Expand doc, add points to CHANGELOG --- CHANGELOG.md | 3 +++ cert_tools/linalg_tools.py | 8 +++++++- 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index d982e56..5d5fc57 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,8 +7,11 @@ and this project adheres to [Semantic Versioning](http://semver.org/). ## [Unreleased] - 2025-09-01 ### Added +- linalg_tools.extract_lower_rank_solution: new algorithm to extract + lower-rank solution from a solution to an SDP. ### Changed +- sdp_solvers.solve_low_rank_sdp: now returns the actuall success status of IPOPT. ### Fixed diff --git a/cert_tools/linalg_tools.py b/cert_tools/linalg_tools.py index 396b1de..2bee1dd 100644 --- a/cert_tools/linalg_tools.py +++ b/cert_tools/linalg_tools.py @@ -118,7 +118,13 @@ def rank_project(X, p=1, tolerance=1e-10) -> tuple[np.ndarray, dict]: def extract_lower_rank_solution(X, A_list, tol=1e-8, max_trials=10): - """Given a psd matrix X of rank r, returns a matrix X' of rank r-1""" + """Given a psd matrix X of rank r, returns a matrix X' of rank r-1. + + This is an implementation of Algorithm 2.1 by [1]. + + [1] Lemon et al. -- Low-Rank Semidefinite Programming Theory and Applications + Foundations and Trends in Optimization, 2015. + """ r_start = 0 for k in range(max_trials): E, U = np.linalg.eigh(X) From 424ea344caa8a14091dc56e23426796d06af4a6a Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Tue, 9 Sep 2025 09:12:12 +0200 Subject: [PATCH 5/6] Bump version --- CHANGELOG.md | 10 +++++++++- cert_tools/__init__.py | 2 +- setup.cfg | 2 +- 3 files changed, 11 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 5d5fc57..47d30a9 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,7 +4,15 @@ All notable changes to this project will be documented in this file. The format is based on [Keep a Changelog](http://keepachangelog.com/) and this project adheres to [Semantic Versioning](http://semver.org/). -## [Unreleased] - 2025-09-01 +## [Unreleased] + +### Added + +### Changed + +### Fixed + +## [0.0.7] - 2025-09-09 ### Added - linalg_tools.extract_lower_rank_solution: new algorithm to extract diff --git a/cert_tools/__init__.py b/cert_tools/__init__.py index 7e0abbd..fe2a767 100644 --- a/cert_tools/__init__.py +++ b/cert_tools/__init__.py @@ -4,4 +4,4 @@ from .linalg_tools import * from .sdp_solvers import * -__version__ = "0.0.6" +__version__ = "0.0.7" diff --git a/setup.cfg b/setup.cfg index 81583ce..e83fb68 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = cert_tools -version = 0.0.6 +version = 0.0.7 authors = [ {name = "Frederike Dümbgen", email = "frederike.dumbgen@utoronto.ca" }, {name = "Connor Holmes", email = "connor.holmes@mail.utoronto.ca" }] From 3f72c8ba399afa5154dfda082bec1a6f74d487db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Frederike=20D=C3=BCmbgen?= Date: Tue, 9 Sep 2025 09:40:38 +0200 Subject: [PATCH 6/6] Improve interface, add low-rank extraction (#18) * Add success flag to ipopt output * Implement low-rank projection algorithm * Change Constraints argument to A_list (b not used) --- cert_tools/linalg_tools.py | 44 ++++++++++++++++++++++++++++++++++++-- cert_tools/sdp_solvers.py | 12 +++++++++-- 2 files changed, 52 insertions(+), 4 deletions(-) diff --git a/cert_tools/linalg_tools.py b/cert_tools/linalg_tools.py index 6d93bb1..396b1de 100644 --- a/cert_tools/linalg_tools.py +++ b/cert_tools/linalg_tools.py @@ -1,3 +1,4 @@ +import warnings from copy import deepcopy import numpy as np @@ -84,7 +85,7 @@ def project_so3(X): return U @ Vh -def rank_project(X, p=1, tolerance=1e-10): +def rank_project(X, p=1, tolerance=1e-10) -> tuple[np.ndarray, dict]: """Project matrix X to matrix of rank p.""" try: assert la.issymmetric(X, atol=tolerance) @@ -113,7 +114,46 @@ def rank_project(X, p=1, tolerance=1e-10): "error eigs": np.sum(np.abs(E[p:])), "EVR": abs(E[p - 1] / E[p]), # largest over second-largest } - return x, info + return np.asarray(x), info + + +def extract_lower_rank_solution(X, A_list, tol=1e-8, max_trials=10): + """Given a psd matrix X of rank r, returns a matrix X' of rank r-1""" + r_start = 0 + for k in range(max_trials): + E, U = np.linalg.eigh(X) + r = int(np.sum(E > tol)) + if r == 1: + return X, {"success": True, "rank": 1, "iter": k} + elif r == r_start - 1: + return X, {"success": True, "rank": r, "iter": k} + if r_start == 0: + r_start = r + V = U[:, -r:] @ np.diag(np.sqrt(E[-r:])) + np.testing.assert_allclose(X, V @ V.T, atol=tol) + + AV = np.vstack([svec(V.T @ Ai @ V) for Ai in A_list]) # type: ignore + + # below, the rows of bi have the nullspace vectors. + # AV @ bi = 0 + b, info_null = get_nullspace(AV, tolerance=tol) + if b.shape[0] == 0: + return X, {"success": True, "rank": r, "iter": k} + break + + # chose one of the vectors + idx = np.random.choice(b.shape[0]) + bi = b[idx] + + assert np.all(AV @ bi < tol) + Delta = smat(bi) + + lambdas = np.abs(np.linalg.eigvalsh(Delta)) + # maximum magnitude eigenvalue + alpha = -1 / np.max(lambdas) + X = V @ (np.eye(r) + alpha * Delta) @ V.T + warnings.warn(f"Did not find a lower-rank solution in {max_trials} trials.") + return X, {"success": True, "rank": r, "iter": k} def find_dependent_columns(A_sparse, tolerance=1e-10, verbose=False, debug=False): diff --git a/cert_tools/sdp_solvers.py b/cert_tools/sdp_solvers.py index c04a352..8fecfef 100644 --- a/cert_tools/sdp_solvers.py +++ b/cert_tools/sdp_solvers.py @@ -170,13 +170,21 @@ def solve_low_rank_sdp( if x_cand is not None: sol_input["x0"] = x_cand.reshape((-1, 1)) r = S(**sol_input) - Y_opt = r["x"] + + # Get solver info + stats = S.stats() + success = stats.get("success", False) + info = {"solver_stats": stats, "success": success} + # Reshape and generate SDP solution + Y_opt = r["x"] Y_opt = np.array(Y_opt).reshape((n, rank), order="F") X_opt = Y_opt @ Y_opt.T + # Get cost scale, offset = adjust cost = np.trace(Q @ X_opt) * scale + offset + # Construct certificate mults = np.array(r["lam_g"]) H = Q @@ -186,7 +194,7 @@ def solve_low_rank_sdp( H = H + A * mults[i, 0] # Return - info = {"X": X_opt, "H": H, "cost": cost} + info = {"X": X_opt, "H": H, "cost": cost, "success": success} return Y_opt, info