[JuMP] Creating a Nonlinear Least Squares model

Specific Domains Optimization (Mathematical)
1

Related to this PR and this post . We are trying to make a Nonlinear Least Squares model using JuMP, so that the matrices are sparse. From @miles.lubin suggestion, I’m building a model passing the residuals as NLexpressions, and creating another model where those expressions are constraints. I can also have constraints, and those are kept in another model.

function problem()
  m = Model()
  @variable(m, x[1:2])
  @NLexpression(m, F1, 10 * (x[2] - x[1]^2)) # Residual
  @NLexpresison(m, F2, x[1] - 1.0) # Residual
  @NLconstraint(m, x[1]^3 - x[2] <= 1.0) # Constraint
  
  mF, mc = NLSModel(m, [F1, F2])
end

function NLSModel(mc, F :: Vector)
  @NLobjective(mc, Min, sum(Fi^2 for Fi in F))
  ev = JuMP.NLPEvaluator(mc)
  MathProgBase.initialize(ev, [:ExprGraph])
  mF = JuMP.Model()
  @objective(mF, Min, 0.0)
  @variable(mF, x[1:2]) # I Want to avoid this also
  for Fi in F
    expr = ev.subexpressions_as_julia_expressions[Fi.index]
    expr = :($expr == 0)
    JuMP.addNLconstraint(mF, expr) # ERROR
  end

  return mF, mc
end

I get the following error:

ERROR: LoadError: Unrecognized expression x[2]. JuMP variable
objects and input coefficients should be spliced directly into expressions.

Thanks in advance.

2

You need to walk the expression and substitute in the actual variables. Something like this seems to work for me.

Note that JuMP will rename all of the variables in the expression to x[i] in the order they were added regardless of their name.

function problem()
  m = Model()
  @variable(m, a)
  @variable(m, b)
  @NLexpression(m, F1, 10 * (b - a^2)) # Residual
  @NLexpression(m, F2, a - 1.0) # Residual
  @NLconstraint(m, a^3 - b <= 1.0) # Constraint

  mF, mc = NLSModel(m, [F1, F2])
end

function NLSModel(mc, F :: Vector)
  @NLobjective(mc, Min, sum(Fi^2 for Fi in F))
  ev = JuMP.NLPEvaluator(mc)
  MathProgBase.initialize(ev, [:ExprGraph])
  mF = JuMP.Model()
  @objective(mF, Min, 0.0)
  @variable(mF, x[1:MathProgBase.numvar(mc)])
  for Fi in F
    expr = ev.subexpressions_as_julia_expressions[Fi.index]
    replace!(expr, x)
    expr = :($expr == 0)
    JuMP.addNLconstraint(mF, expr)
  end

  return mF, mc
end

function replace!(ex, x)
    if isa(ex, Expr)
        for (i,arg) in enumerate(ex.args)
            if isa(arg, Expr)
                if arg.head == :ref && arg.args[1] == :x
                    ex.args[i] = x[arg.args[2]]
                else
                    replace!(arg, x)
                end
            end
        end
    end
end

mF, mc = problem()

@show mF[:x]
@show mc[:a]
1 Like
3

Thanks @odow, that solved the problem.