Convex.jl matrix constraints interface : How can i supply this matrix?

Specific Domains Optimization (Mathematical)
1

Hi,

I want to constraint some matrix that is given as a linear functions of the variables to be semi-definite positive. However, I found something strange :

function get_hankels(s,a,b)
    m = length(s)-1
    if iseven(m)
        # m = 2n
        n = Int(m//2)
        H_low = [s[i+j-1] for i in 1:n+1,j in 1:n+1]
        H_up = [(a+b)*s[i+j] - s[i+j+1] - a*b*s[i+j-1] for i in 1:n,j in 1:n]
    else
        # m = 2n+1
        n = Int((m-1)//2)
        H_low = [s[i+j] - a*s[i+j-1] for i in 1:n+1,j in 1:n+1]
        H_up = [b*s[i+j-1] - s[i+j] for i in 1:n+1,j in 1:n+1]
    end
    return H_low,H_up
end

# To see the matrices I want to use : 
using DynamicPolynomials
m=6
@polyvar x[1:m+1]
H_low,H_up = get_hankels(x,0,1)
display(H_low)
display(H_up)

# Now with convex variables, I want to build the constraints : 
using Convex
m=6
s = Variable(m+1)
H_low, H_up = get_hankels(s,0,1) # works

cstr1 = (H_low âȘ° 0) # does not work. 
cstr1 = ([s[1] s[2] s[3] s[4]; s[2] s[3] s[4] s[5]; s[3] s[4] s[5] s[6]; s[4] s[5] s[6] s[7]] âȘ° 0) # works

# Same for H_up 
cstr2 = (H_up âȘ° 0)

I have the follwoing error message :

julia> cstr1 = (H_low âȘ° 0) # does not work.
ERROR: MethodError: no method matching âȘ°(::Matrix{Convex.IndexAtom}, ::Int64)
Closest candidates are:
  âȘ°(::Convex.AbstractExpr, ::Union{Number, AbstractArray}) at C:\Users\u009192\.julia\packages\Convex\ubaUN\src\constraints\sdp_constraints.jl:115
  âȘ°(::Union{Number, AbstractArray}, ::Convex.AbstractExpr) at C:\Users\u009192\.julia\packages\Convex\ubaUN\src\constraints\sdp_constraints.jl:124
Stacktrace:
 [1] top-level scope
   @ REPL[58]:1

julia>

Indeed, the type is not the same :

julia> typeof([s[1] s[2] s[3] s[4]; s[2] s[3] s[4] s[5]; s[3] s[4] s[5] s[6]; s[4] s[5] s[6] s[7]])
Convex.TransposeAtom

julia> typeof(H_low)
Matrix{Convex.IndexAtom} (alias for Array{Convex.IndexAtom, 2})

julia>

Is there a way to build this H_low matrix automatically in the right type ? The size of the problem might change from a run to the other


@ericphanson Do you now the right syntax ?

2

Yeah, this is something Convex struggles with, which is that it wants e.g. a 2x2 matrix to be Variable(2,2), not Matrix{Variable}. In other words, Convex works with vector/matrix-valuned variables but can’t work with vectors/matrices of variables directly. One workaround is to do

function convert_to_variable(arr::AbstractMatrix{<:AbstractVariable})
    x = Variable(size(arr))
    for i = 1:size(arr, 1), j = 1:size(arr, 2)
        add_constraint!(x, x[i,j] == arr[i,j])
    end
    return x
end

(untested), and then do H_low_variable = convert_to_variable(H_low) and use H_low_variable from then on.

This might not be very efficient though, as it creates a new variable with one constraint per variable in your matrix. For Convex.jl, a more efficient approach would be to rewrite get_hankels in a “vectorized” way, i.e. try to write it something like a matrix multiplication (of a non-variable matrix) onto your variable vector s. In other words, try to avoid indexing into s and create H_low by a function of s that doesn’t need to index into it.

Since JuMP works on a scalar level, that might be a better approach for this kind of problem.

2 Likes
3

I only changed the AbstractVariable to AbstractExprOrValue as is :

function convert_to_variable(arr::AbstractMatrix{<:Convex.AbstractExprOrValue})
    x = Variable(size(arr))
    for i = 1:size(arr, 1), j = 1:size(arr, 2)
        add_constraint!(x, x[i,j] == arr[i,j])
    end
    return x
end

I had Matrix{Convex.IndexAtom}, and Convex.IndexAtom <: Convex.AbstractVariable is false but Convex.IndexAtom <: Convex.AbstractExprOrValue is true.

Works like a charm, thanks a lot @ericphanson !

Edit: I’ll try rewriting it as matrix-vector product, as it is clearly possible in this case. It was just easier to reason about with indexing. Jump is out-of-scope since I need arbitrary precision :wink:

1 Like