From c547ab2c65011fb06420bbd3c5b802834c709899 Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Mon, 16 Jun 2025 16:26:02 +0200 Subject: [PATCH 1/7] Fix rank_project for p>1 --- CHANGELOG.md | 3 ++- cert_tools/linalg_tools.py | 5 ++++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a59ac4a..ebd4685 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,13 +4,14 @@ 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] - YYYY-MM-DD +## [Unreleased] - 2025-06-16 ### Added ### Changed ### Fixed +- rank_project working with p>1. ## [0.0.5] - 2025-05-27 diff --git a/cert_tools/linalg_tools.py b/cert_tools/linalg_tools.py index 1323bb8..7e99079 100644 --- a/cert_tools/linalg_tools.py +++ b/cert_tools/linalg_tools.py @@ -93,7 +93,10 @@ def rank_project(X, p=1, tolerance=1e-10): p = np.sum(np.abs(E) > tolerance) x = V[:, -p:] * np.sqrt(E[-p:]) - X_hat = np.outer(x, x) + if p == 1: + X_hat = np.outer(x, x) + else: + X_hat = x @ x.T info = { "error X": np.linalg.norm(X_hat - X), "error eigs": np.sum(np.abs(E[:-p])), From 18dbfa8a37b3e6209fbb5a3c44eabe39e9c60e2c Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Mon, 1 Sep 2025 18:16:20 +0200 Subject: [PATCH 2/7] Add more options to IPOPT --- cert_tools/sdp_solvers.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/cert_tools/sdp_solvers.py b/cert_tools/sdp_solvers.py index 8236d9b..ab08785 100644 --- a/cert_tools/sdp_solvers.py +++ b/cert_tools/sdp_solvers.py @@ -125,6 +125,7 @@ def solve_low_rank_sdp( options={}, limit_constraints=False, verbose=False, + method="ipopt", ): """Use the factorization proposed by Burer and Monteiro to solve a fixed rank SDP. @@ -150,9 +151,20 @@ def solve_low_rank_sdp( g_rhs = cas.vertcat(*g_rhs) # Define Low Rank NLP nlp = {"x": Y.reshape((-1, 1)), "f": f, "g": g_lhs} - options["ipopt.print_level"] = int(verbose) options["print_time"] = int(verbose) - S = cas.nlpsol("S", "ipopt", nlp, options) + options["print_level"] = 5 + options["error_on_fail"] = False + options["warm_start_init_point"] = "yes" + options["mu_init"] = 1e2 + options["mu_strategy"] = "adaptive" + options["warm_start_bound_push"] = 1e-6 + options["warm_start_mult_bound_push"] = 1e-6 + if method == "ipopt": + options["ipopt.print_level"] = int(verbose) + 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: From 2ec1209d67b2fb1f344e1f6af087f7d22d7b3bfb Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Mon, 1 Sep 2025 18:17:24 +0200 Subject: [PATCH 3/7] Add tolerance to issymmetric check --- cert_tools/linalg_tools.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cert_tools/linalg_tools.py b/cert_tools/linalg_tools.py index 7e99079..775e26e 100644 --- a/cert_tools/linalg_tools.py +++ b/cert_tools/linalg_tools.py @@ -87,7 +87,7 @@ def project_so3(X): def rank_project(X, p=1, tolerance=1e-10): """Project matrix X to matrix of rank p.""" try: - assert la.issymmetric(X) + assert la.issymmetric(X, atol=tolerance) E, V = np.linalg.eigh(X) if p is None: p = np.sum(np.abs(E) > tolerance) From 3906000fe759985b37e6cd0e1f217f9459d667f7 Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Mon, 1 Sep 2025 19:16:43 +0200 Subject: [PATCH 4/7] Fix ipopt options --- _test/test_sdp_solvers.py | 8 ++++---- cert_tools/sdp_solvers.py | 12 ++++++------ 2 files changed, 10 insertions(+), 10 deletions(-) diff --git a/_test/test_sdp_solvers.py b/_test/test_sdp_solvers.py index acddb44..dec49c2 100644 --- a/_test/test_sdp_solvers.py +++ b/_test/test_sdp_solvers.py @@ -6,10 +6,10 @@ import numpy as np from cert_tools import ( - solve_sdp_mosek, solve_low_rank_sdp, - solve_sdp_fusion, solve_sdp_cvxpy, + solve_sdp_fusion, + solve_sdp_mosek, ) root_dir = os.path.abspath(os.path.dirname(__file__) + "/../") @@ -178,5 +178,5 @@ def test_sdp_solvers( continue # test_sdp_solvers(f"test_prob_{i}.pkl", primal_list=[False], verbose=True) test_sdp_solvers(f"test_prob_{i}.pkl") - # test_p1_low_rank() - # test_p3_low_rank() + test_p1_low_rank() + test_p3_low_rank() diff --git a/cert_tools/sdp_solvers.py b/cert_tools/sdp_solvers.py index ab08785..c04a352 100644 --- a/cert_tools/sdp_solvers.py +++ b/cert_tools/sdp_solvers.py @@ -152,14 +152,14 @@ def solve_low_rank_sdp( # Define Low Rank NLP nlp = {"x": Y.reshape((-1, 1)), "f": f, "g": g_lhs} options["print_time"] = int(verbose) - options["print_level"] = 5 options["error_on_fail"] = False - options["warm_start_init_point"] = "yes" - options["mu_init"] = 1e2 - options["mu_strategy"] = "adaptive" - options["warm_start_bound_push"] = 1e-6 - options["warm_start_mult_bound_push"] = 1e-6 if method == "ipopt": + options["ipopt.warm_start_init_point"] = "yes" + options["ipopt.warm_start_bound_push"] = 1e-6 + options["ipopt.warm_start_mult_bound_push"] = 1e-6 + options["ipopt.print_level"] = 5 + options["ipopt.mu_init"] = 1e2 + options["ipopt.mu_strategy"] = "adaptive" options["ipopt.print_level"] = int(verbose) S = cas.nlpsol("S", "ipopt", nlp, options) else: From 0331cde16acf4e7bfa1d87cddcef8df025e7cc33 Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Mon, 1 Sep 2025 19:17:02 +0200 Subject: [PATCH 5/7] Fix and test rank_project --- _test/test_linalg_tools.py | 68 ++++++++++++++++++++++++++++++++++---- cert_tools/linalg_tools.py | 13 ++++---- 2 files changed, 69 insertions(+), 12 deletions(-) diff --git a/_test/test_linalg_tools.py b/_test/test_linalg_tools.py index 8d8acc4..25b713a 100644 --- a/_test/test_linalg_tools.py +++ b/_test/test_linalg_tools.py @@ -1,10 +1,27 @@ import unittest -import matplotlib.pyplot as plt import numpy as np import scipy.sparse as sp -from cert_tools.linalg_tools import smat, svec +from cert_tools.linalg_tools import rank_project, smat, svec + + +def best_rotation_nd(X, Y): + """ + X, Y: (N, d) arrays with corresponding row points. + Returns: R (dxd orthogonal matrix), + such that Y ≈ (X @ R.T) + """ + cx = X.mean(axis=0) + cy = Y.mean(axis=0) + Xc = X - cx + Yc = Y - cy + + H = Xc.T @ Yc # d x d + U, S, Vt = np.linalg.svd(H) + V = Vt.T + R = V @ U.T + return R class TestLinAlg(unittest.TestCase): @@ -30,7 +47,7 @@ def test_svec(self): np.testing.assert_almost_equal(smat(s2), S2) # products should be equal prod_mat = np.trace(S1 @ S2) - prod_vec = np.dot(s1, s2) + prod_vec = np.dot(s1, s2) # type: ignore assert abs(prod_mat - prod_vec) < 1e-10, "PSD Inner product not equal" def test_svec_sparse(self): @@ -51,18 +68,57 @@ def test_svec_sparse(self): s2 = svec(S2) s1_dense = svec(S1_dense) s2_dense = svec(S2_dense) - np.testing.assert_almost_equal(s1_dense, s1.toarray().squeeze(0)) - np.testing.assert_almost_equal(s2_dense, s2.toarray().squeeze(0)) + np.testing.assert_almost_equal(s1_dense, s1.toarray().squeeze(0)) # type: ignore + np.testing.assert_almost_equal(s2_dense, s2.toarray().squeeze(0)) # type: ignore # test mapping np.testing.assert_almost_equal(smat(s1), S1.toarray()) np.testing.assert_almost_equal(smat(s2), S2.toarray()) # products should be equal prod_mat = np.trace(S1.toarray() @ S2.toarray()) - prod_vec = (s1 @ s2.T).toarray() + prod_vec = (s1 @ s2.T).toarray() # type: ignore assert abs(prod_mat - prod_vec) < 1e-10, "PSD Inner product not equal" + def test_rank_project(self): + # test symmetric rank-one matrix + x_gt = np.random.randn(5, 1) + X = x_gt @ x_gt.T + x_test, info_rank = rank_project(X, p=1) + if x_gt[0] * x_test[0] < 0: + x_test = -x_test + np.testing.assert_allclose(x_test, x_gt) + + # test symmatric rank-two matrix + # below we need to compensate for a rotation because + # X = U E U' = U sqrt(E) sqrt(E) U' = U sqrt(E) R' R sqrt(E) U' + # for any orthogonal R. + x_gt = np.random.randn(5, 2) + X = x_gt @ x_gt.T + x_test, info_rank = rank_project(X, p=2) + R = best_rotation_nd(x_test, x_gt) + x_test = x_test @ R.T # align for comparison + np.testing.assert_allclose(x_test, x_gt) + + # test non-symmetric rank-one matrix + x_gt = np.random.randn(5, 1) + X = x_gt @ x_gt.T + X[0, 2] += 1e-8 + x_test, info_rank = rank_project(X, p=1) + if x_gt[0] * x_test[0] < 0: + x_test = -x_test + np.testing.assert_allclose(x_test, x_gt, atol=1e-5) + + # test non-symmetric rank-two matrix + x_gt = np.random.randn(5, 2) + X = x_gt @ x_gt.T + X[0, 2] += 1e-8 + x_test, info_rank = rank_project(X, p=2) + R = best_rotation_nd(x_test, x_gt) + x_test = x_test @ R.T # align for comparison + np.testing.assert_allclose(x_test, x_gt, atol=1e-5) + if __name__ == "__main__": test = TestLinAlg() test.test_svec() test.test_svec_sparse() + test.test_rank_project() diff --git a/cert_tools/linalg_tools.py b/cert_tools/linalg_tools.py index 775e26e..6d93bb1 100644 --- a/cert_tools/linalg_tools.py +++ b/cert_tools/linalg_tools.py @@ -8,7 +8,7 @@ NULL_THRESH = 1e-5 -def svec(S, order="C"): +def svec(S, order="C") -> np.ndarray | sp.csc_array: """Convert symmetric matrix to vectorized form Args: @@ -90,7 +90,7 @@ def rank_project(X, p=1, tolerance=1e-10): assert la.issymmetric(X, atol=tolerance) E, V = np.linalg.eigh(X) if p is None: - p = np.sum(np.abs(E) > tolerance) + p = int(np.sum(np.abs(E) > tolerance)) x = V[:, -p:] * np.sqrt(E[-p:]) if p == 1: @@ -107,7 +107,7 @@ def rank_project(X, p=1, tolerance=1e-10): if p is None: p = np.sum(np.abs(E) > tolerance) X_hat = U[:, :p] @ np.diag(E[:p]) @ Vh[:p, :] - x = U[:, :p] @ np.diag(E[:p]) + x = U[:, :p] @ np.diag(np.sqrt(E[:p])) info = { "error X": np.linalg.norm(X_hat - X), "error eigs": np.sum(np.abs(E[p:])), @@ -137,6 +137,7 @@ def find_dependent_columns(A_sparse, tolerance=1e-10, verbose=False, debug=False print(f"clean_constraints: keeping {rank}/{A_sparse.shape[1]} independent") bad_idx = list(range(A_sparse.shape[1])) + assert E is not None good_idx_list = sorted(E[sort_inds[:rank]])[::-1] for good_idx in good_idx_list: del bad_idx[good_idx] @@ -181,8 +182,8 @@ def get_nullspace(A_dense, method=METHOD, tolerance=NULL_THRESH): # Based on Section 5.5.5 "Basic Solutions via QR with Column Pivoting" from Golub and Van Loan. # assert A_dense.shape[0] >= A_dense.shape[1], "only tall matrices supported" - Q, R, P = la.qr(A_dense, pivoting=True, mode="economic") - if Q.shape[0] < 1e4: + Q, R, P = la.qr(A_dense, pivoting=True, mode="economic") # type: ignore + if Q.shape[0] < 1e4: # type: ignore np.testing.assert_allclose(Q @ R, A_dense[:, P], atol=1e-5) S = np.abs(np.diag(R)) @@ -199,7 +200,7 @@ def get_nullspace(A_dense, method=METHOD, tolerance=NULL_THRESH): Pinv[p] = k LHS = R1[:, Pinv] - info["Q1"] = Q[:, :rank] + info["Q1"] = Q[:, :rank] # type: ignore info["LHS"] = LHS basis = np.zeros(N.T.shape) From 85f25e0556e0c05305d5e3386ab01cb017a7e910 Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Mon, 1 Sep 2025 19:17:12 +0200 Subject: [PATCH 6/7] Update CHANGELOG --- CHANGELOG.md | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ebd4685..fff181d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,14 +4,17 @@ 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-06-16 +## [Unreleased] - 2025-09-01 ### Added +- sdp_solvers.solve_low_rank_sdp: expose options of ipopt. ### Changed ### Fixed - rank_project working with p>1. +- rank_project working with almost-symmetric matrices. +- rank_project working with non-symmetric matrices. ## [0.0.5] - 2025-05-27 From 5cf3d1a85b32ad7f4017f42c7fe332d62d0bf389 Mon Sep 17 00:00:00 2001 From: Frederike Duembgen Date: Mon, 1 Sep 2025 19:35:47 +0200 Subject: [PATCH 7/7] Bump version to 0.0.6 --- CHANGELOG.md | 8 +++++++- cert_tools/__init__.py | 2 +- setup.cfg | 2 +- 3 files changed, 9 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fff181d..d982e56 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,10 +7,16 @@ and this project adheres to [Semantic Versioning](http://semver.org/). ## [Unreleased] - 2025-09-01 ### Added -- sdp_solvers.solve_low_rank_sdp: expose options of ipopt. ### Changed +### Fixed + +## [0.0.6] - 2025-09-01 + +### Added +- sdp_solvers.solve_low_rank_sdp: expose options of ipopt. + ### Fixed - rank_project working with p>1. - rank_project working with almost-symmetric matrices. diff --git a/cert_tools/__init__.py b/cert_tools/__init__.py index 59fea62..7e0abbd 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.5" +__version__ = "0.0.6" diff --git a/setup.cfg b/setup.cfg index 92320b7..81583ce 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = cert_tools -version = 0.0.5 +version = 0.0.6 authors = [ {name = "Frederike Dümbgen", email = "frederike.dumbgen@utoronto.ca" }, {name = "Connor Holmes", email = "connor.holmes@mail.utoronto.ca" }]