From 7149b9379b8ec50ee093fa01e91a3f4bed741141 Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Mon, 1 Sep 2025 20:13:41 +0200 Subject: [PATCH 1/3] 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/3] 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/3] 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])