cvxpy

Write quadratic fractional program in CVXPY


I am trying the write the following quasiconvex function in CVXPY, which I eventually want to minimize.

# Import packages
import cvxpy as cp
import numpy as np

# parameters
A1 = np.matrix([[2,1],[1,3]])
A2 = np.matrix([[1,1],[1,4]])

# create model
x = cp.Variable(2)
f = cp.quad_form(x,A1,assume_PSD=True)
h = -cp.quad_form(x,A2,assume_PSD=True) + 10
quadFrac = f / h

# check the function type
print("f is convex: ",f.is_convex())
print("h is concave: ",h.is_concave())
print("f/h is quasiconvex: ",quadFrac.is_quasiconvex())
print("f/h is ",quadFrac.curvature)

I know that this function f/h is quasiconvex for all x where h(x)>0 based on Boyd's textbook exercise 4.7. However, when I check the function curvature it says it is unknown. I think this is due to not following explicit DQCP composition rules.

So, how can I write this function in CVXPY for it to be recognized as quasiconvex?

Based on the DQCP atoms, the Ratio atom seems to be able to accommodate my function where it says

...the ratio of a nonnegative convex function and a positive concave function is quasiconvex.


Solution

  • Please note that this is OUTSIDE the standard use of CVXPY. So be careful about any decisions based on the results

    I discussed this with some people on Discord. Taking into account DQCP, wraps, signs, and curvature, this an appropriate (hacky) way to solve the problem

    # import
    import cvxpy as cp
    import numpy as np
    
    # parameters
    A1 = np.matrix([[2,1],[1,3]])
    A2 = np.matrix([[1,1],[1,4]])
    
    # create variables and functions
    x = cp.Variable(2)
    f = cp.quad_form(x,A1,assume_PSD=True)
    h = -cp.quad_form(x,A2,assume_PSD=True) + 10
    
    # assertions
    h = cp.atoms.affine.wraps.nonneg_wrap(h)
    
    # function
    F = f/h
    
    # check assertions
    print("f is ",f.sign,f.curvature)
    print("h is ",h.sign,h.curvature)
    print("t is ",t.sign,t.curvature)
    print("F is ",F.sign,F.curvature)
    
    # constraints
    constraints = [h >= 1e-4]
    
    # problem
    problem = cp.Problem(cp.Minimize(F),constraints)
    
    # solve
    assert problem.is_dqcp()
    problem.solve(qcp=True)
    
    # print results
    print(x.value)