Minor fixes

CHANGELOG.md

@@ -4,14 +4,24 @@ 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-09-01
 
 ### Added 
 
 ### 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.
+- rank_project working with non-symmetric matrices.
+
 ## [0.0.5] - 2025-05-27
 
 ### Added 

_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()

_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()

cert_tools/__init__.py

@@ -4,4 +4,4 @@
 from .linalg_tools import *
 from .sdp_solvers import *
 
-__version__ = "0.0.5"
+__version__ = "0.0.6"

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:
@@ -87,13 +87,16 @@ 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)
+            p = int(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])),
@@ -104,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:])),
@@ -134,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]
@@ -178,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))
@@ -196,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)

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["error_on_fail"] = False
+    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:
+        S = cas.nlpsol("S", method, nlp, options)
+
     # Run Program
     sol_input = dict(lbg=g_rhs, ubg=g_rhs)
     if x_cand is not None:

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" }]