Expressing positive semidefiniteness of block matrices in CVXPY

652 Views Asked by At

I am trying to solve a robust optimization problem using CVXPY. I am new to this package and also a beginner of Python. I don't know how to implement the constraint that makes a block decision variable matrix to be positive semidefinite.


enter image description here


Can anyone help with it?

1

There are 1 best solutions below

2
On BEST ANSWER

First, you are lucky because a similar code took me a lot of time! The following code misses the second term in the objective function (so you should add that one) but it definitely has the overall structure of the code. Feel free to edit my response if there are issues (or comment so I can edit it), I just did not want you to suffer looking for syntax as I did. At the same time, unfortunately, I do not have time to debug.

def solving_SDP(mu,delta,lambd,S_bar,S_under):
    N = np.shape(S_bar)[1] #number of columns of Sigma_bar
    # construct the variables
    G_bar = cp.Variable(S_bar.shape, symmetric = True)  
    G_under = cp.Variable(S_under.shape, symmetric = True)  
    w = cp.Variable(shape=(mu.shape[0],1)) 
    V = cp.Variable((N+1,N+1), symmetric=True)
    V = cp.bmat([[G_bar-G_under, w], [w.T, np.eye(1)]])
    objective =  cp.trace(mu.T@w)-lambd*(cp.trace(G_bar@S_bar)-cp.trace(G_under@S_under))
    constraints = [V >> 0]
    constraints += [G_bar >> 0]
    constraints += [G_under >> 0]
    constraints += [cp.norm([email protected](shape=(N,1)), 1) == 1]
    prob = cp.Problem(cp.Minimize(objective), constraints)
    prob.solve()
    #print("Status of problem:", prob.status)
    #print("Optimal value of problem is", prob.value)
    return w,G_bar,G_under
    
#Test it:
N = 10
S_bar = make_spd_matrix(N, random_state=None)
S_under = make_spd_matrix(N, random_state=None)
mu = np.random.randn(N,1)
delta = np.random.randn(N,1)
lambd = 1
w,G_bar,G_under = solving_SDP(mu,delta,lambd,S_bar,S_under)
print("w = \n ", w)
eigVals, eigVecs = LA.eig(G_bar)
eigVals = -np.sort(-eigVals) #sort eigenvalues
print("Minimum eigenvalue of G_bar = ", eigVals[N-1])