[Python] The Semidefinite Programming solver code with arbitrary precision

I'm trying to find a SDP solver with arbitrary precision support, to solve the standard semidefinite program as shown in

https://irene.readthedocs.io/en/latest/sdp.html
minmize sum_{i=1}^m b_i x_i
subject to
sum_{i=1}^m A_{ij} x_i - C_j >= 0 j = 1,...,k


where A_{ij} and C_j were some square matrices, and j the index for the set of constraints. The x_i were some arbitrary variables(scalar) to be optimized.

However, it's been challenging since most SDP solvers did not support arbitrary precision, which is necessary. The SDPA for example, had stopped maintenance around 2015, and even the website stopped working.

I tried to use the famous SDPB solver

https://github.com/davidsd/sdpb/tree/master
#the document used x for A_{ij} and C_j dependence, and y_n for the scalar variables.


However, when the matrices A_{ij} and C_j are constant matrices instead of the polynomial matrices, the SDPB were not very stable.

I then found a SDPJ solver, which the findFeasible solve the question

https://github.com/FishboneChiang/SDPJSolver.jl


But it's native in Julia, and was very hard to use with the python environment.

I therefore converted it using LLM. One of those models worked

import mpmath
import numpy as np

def matrix_trace(matrix):
"""Calculate the trace of a matrix (sum of diagonal elements)."""
n = matrix.rows
trace = mpmath.mpf(0)
for i in range(n):
trace += matrix[i, i]
return trace

def solve_sdp_feasibility(A_list, b, gamma=1.0, max_iter=1000, tolerance=1e-10, dps=50, X0=None):
"""
Solves the SDP feasibility problem using the Douglas-Rachford splitting method.

Parameters:
A_list (list): List of constraint matrices (each as a list of lists or mpmath matrix).
b (list): List of constraint values.
gamma (float): Step size for the Douglas-Rachford update.
max_iter (int): Maximum number of iterations.
tolerance (float): Convergence tolerance.
dps (int): Decimal places precision for mpmath.
X0 (mpmath.matrix): Initial guess for the solution matrix.

Returns:
mpmath.matrix: Solution matrix X that satisfies the constraints and is PSD.
"""
mpmath.mp.dps = dps
n = len(A_list[0])
m = len(b)

# Convert A_list and b to mpmath matrices and numbers
A_mpmath = []
for A in A_list:
A_mpmath.append(mpmath.matrix(A))
b_mpmath = [mpmath.mpf(b_i) for b_i in b]

# Precompute M matrix: M[j] = trace(A_i * A_j)
M = mpmath.zeros(m, m)
for i in range(m):
for j in range(m):
# Calculate trace(A_i * A_j)
product = A_mpmath * A_mpmath[j]
M[i, j] = matrix_trace(product)

# Precompute the inverse of M
try:
M_inv = mpmath.inverse(M)
except mpmath.MatrixError:
raise ValueError("Matrix M is singular. Cannot compute inverse.")

# Initialize X
if X0 is None:
X = mpmath.zeros(n, n)
else:
X = mpmath.matrix(X0)

# Define projection functions
def proj_affine(X):
r = mpmath.zeros(m, 1)
for i in range(m):
r = matrix_trace(A_mpmath * X) - b_mpmath
lambda_vec = M_inv * r
update = mpmath.zeros(n, n)
for i in range(m):
update += lambda_vec * A_mpmath
return X - update

def proj_psd(X):
# Ensure symmetry
X = (X + X.T) / 2

# Convert to numpy for eigenvalue decomposition (more reliable)
X_np = np.array(X.tolist(), dtype=float)
D_np, V_np = np.linalg.eigh(X_np)

# Set negative eigenvalues to zero
D_plus_np = np.maximum(D_np, 0)

# Reconstruct the matrix
X_psd_np = V_np @ np.diag(D_plus_np) @ V_np.T

# Convert back to mpmath matrix
return mpmath.matrix(X_psd_np.tolist())

# Main iteration loop
for it in range(max_iter):
Y = proj_affine(X)
U = 2 * Y - X
Z = proj_psd(U)
X_next = X + gamma * (Z - Y)

# Calculate Frobenius norm of the difference
diff = X_next - X
frob_norm = mpmath.sqrt(matrix_trace(diff.T * diff))

if frob_norm < tolerance:
print(f"Converged after {it+1} iterations")
break

X = X_next
else:
print(f"Reached maximum iterations ({max_iter}) without convergence")

return X

# Example usage with more complex constraints:
if __name__ == "__main__":
# Example: Find a 3x3 PSD matrix X such that:
# trace(A1 * X) = b1
# trace(A2 * X) = b2
# trace(A3 * X) = b3

# Define constraint matrices A_i (these are fixed)
A1 = [[1, 2, 3],
[2, 4, 5],
[3, 5, 6]]

A2 = [[0, 1, 0],
[1, 0, 1],
[0, 1, 0]]

A3 = [[2, 0, 1],
[0, 3, 0],
[1, 0, 2]]

# Right-hand side values
b = [10, 5, 8]

# Solve the SDP
X_sol = solve_sdp_feasibility([A1, A2, A3], b, dps=50)

print("Solution matrix X:")
for i in range(3):
print(f"[{X_sol[i, 0]}, {X_sol[i, 1]}, {X_sol[i, 2]}]")

# Verify all constraints
A_matrices = [A1, A2, A3]
for i in range(3):
A_mp = mpmath.matrix(A_matrices)
constraint_value = matrix_trace(A_mp * X_sol)
print(f"Constraint {i+1}: {constraint_value} (expected: {b})")

# Check if X is PSD (all eigenvalues non-negative)
X_np = np.array(X_sol.tolist(), dtype=float)
eigenvalues = np.linalg.eigvalsh(X_np)
print(f"Eigenvalues of X: {eigenvalues}")
print(f"Is X PSD? {np.all(eigenvalues >= 0)}")

# Additional example with more constraints
print("\nAdditional example with 4 constraints:")

A4 = [[1, 0, 0],
[0, 0, 0],
[0, 0, 0]]

A5 = [[0, 0, 0],
[0, 1, 0],
[0, 0, 0]]

A6 = [[0, 0, 0],
[0, 0, 0],
[0, 0, 1]]

A7 = [[0, 1, 0],
[1, 0, 0],
[0, 0, 0]]

b2 = [1, 2, 3, 4]

X_sol2 = solve_sdp_feasibility([A4, A5, A6, A7], b2, dps=50)

print("Solution matrix X:")
for i in range(3):
print(f"[{X_sol2[i, 0]}, {X_sol2[i, 1]}, {X_sol2[i, 2]}]")

# Verify all constraints
A_matrices2 = [A4, A5, A6, A7]
for i in range(4):
A_mp = mpmath.matrix(A_matrices2)
constraint_value = matrix_trace(A_mp * X_sol2)
print(f"Constraint {i+1}: {constraint_value} (expected: {b2})")

# Check if X is PSD
X_np2 = np.array(X_sol2.tolist(), dtype=float)
eigenvalues2 = np.linalg.eigvalsh(X_np2)
print(f"Eigenvalues of X: {eigenvalues2}")
print(f"Is X PSD? {np.all(eigenvalues2 >= 0)}")


The code outputted the results to be

Reached maximum iterations (1000) without convergence
Solution matrix X:
[235.48002885886930412540835581970375114319575650077, -308.36933515559355018981090516603668191002666830055, 199.79293247284757166979960523751779258037136206704]
[-308.36933515559355018981090516603668191002666830055, 433.68754529824929370872769794228559241615115634857, -209.69987778500855636981793146901967731772005678184]
[199.79293247284757166979960523751779258037136206704, -209.69987778500855636981793146901967731772005678184, 400.93037961268435538239059000611706001255961675798]
Constraint 1: 2244.0939640925982768160373836963117355480790138246 (expected: 10)
Constraint 2: -1036.1384258812042131192576732701127184554934501648 (expected: 5)
Constraint 3: 2973.4693177835503434813801959535339847207069396973 (expected: 8)
Eigenvalues of X: [ 5.34290205 219.83433164 844.92072007]
Is X PSD? True

Additional example with 4 constraints:
Reached maximum iterations (1000) without convergence
Solution matrix X:
[350.2235623147925538400926370741217397153377532959, -270.6655180158920781519071852017077617347240447998, 0.0]
[-270.66551801589207792986258027667645364999771118164, 214.89080330684650510519873023440595716238021850586, 0.0]
[0.0, 0.0, 3.0]
Constraint 1: 350.2235623147925538400926370741217397153377532959 (expected: 1)
Constraint 2: 214.89080330684650510519873023440595716238021850586 (expected: 2)
Constraint 3: 3.0 (expected: 3)
Constraint 4: -541.33103603178415608176976547838421538472175598145 (expected: 4)
Eigenvalues of X: [ 3. 3.56155281 561.55281281]
Is X PSD? True


It's using Douglas-Rachford splitting method for solving semidefinite programming (SDP) , but I don't quite understand the code. Especially, the X_i was supposed to be some unknown variables, not the unknown matrices, i.e. each X_i was supposed to be a scalar.

Could you explain what the code was doing? Why all the eigenvalues of X_sol are non negative led to positive semidefinite condition? What if some of the eigenvalues are non negative, some are zeros, and some are positive?

How to convert the SDPJ or SDPB to a standard SDP solver where A_i are constant matrices and X_i are scalars?(Somehow I felt the Douglas-Rachford splitting method seemed to be more stable than the Newtonian gradient based algorithms, but I'm not sure.) How to write a SDP solver with arbitrary precision support in python?

Continue reading...