expressions that involve operations on small vectors

[schrimpf]

Is there a recommended way to write constraints and objectives that involve operations on vectors?

For a simple example, suppose we want the constraint Ax >= 0 , where A is a given n by d matrix, and x is a d vector of variables. Some obvious ways to express this constraint do not work, which is somewhat related to #102

import ExaModels

c = ExaModels.ExaCore()
d = 10
n = 10_000
x = ExaModels.variable(c, d, start = 0.0)
A = randn(n,d)
Arows = eachrow(A)
#  A*x ≥ 0
ExaModels.constraint(c, dot(a,x) for a in Arows; ucon=Inf)  # MethodError: no method matching iterate
ExaModels.constraint(c, a'*x for a in Arows; ucon=Inf)  # MethodError: no method matching adjoint
ExaModels.constraint(c, sum(x.*a) for a in Arows; ucon=Inf)  # MethodError: no method matching length

A work around is to unroll the loop to avoid iteration or length.

@generated function unrolldot(x,y, ::Val{N}) where N
  ex = :(x[1]*y[1])
  for i in 2:N
    ex = :( $ex + x[$i]*y[$i])
  end
  return ex
end
ExaModels.constraint(c, unrolldot(x,a, Val(d))  for a in Arows; ucon=Inf)

Of course, this doesn't scale great with d, but it seems to be quite performant for moderate d up to around 100.

To complicate things further, suppose the constraint or objective is a function of a collection multiple low dimensional vectors. As an example, consider a least squares problem.

# min ||A*x - b||^2
b = randn(n)
data = [(a=a, b=bi) for (a,bi) in zip(eachrow(A), b)]
ExaModels.objective(c, (unrolldot(x,di.a, Val(d)) - di.b)^2 for di in data)  # MethodError: no method matching getindex

This fails due to the two levels of indexing in di.a[1] (one to get a inside di and then another to get [1] inside di.a.

A workaround is to organize data so that everything relevant is available with one level of indexing. My approach has been to convert the vector a to a named tuple (a1=a[1], a2=a[2], ... ) and then use another generated function to create the dot product with x.

function vectotuple(x::AbstractVector{T}, prefix, ::Val{N}) where {T, N}
  NamedTuple{ntuple(i -> Symbol(prefix,i), N)}(ntuple(i -> x[i], N))
end

@generated function unrolldot(x,y,prefix::Symbol, ::Val{N}) where N
  ex = :(x[1]*y[Symbol(prefix,1)])
  for i in 2:N
    ex = :( $ex + x[$i]*y[Symbol(prefix,$i)])
  end
  return ex
end

data= [(vectotuple(a, :a, Val(d))..., b=bi) for (a,bi) in zip(eachrow(A), b)]
ExaModels.objective(c, (unrolldot(x,di,:a, Val(d)) - di.b)^2 for di in data)

This works and is performant, but the code starts to feel quite arcane. I am wondering whether there is a better way or maybe what I'm doing is just a bad idea.