Handling nonlinear objective involving the inverse of a matrix

Specific Domains Optimization (Mathematical)
1

so the output of:

MOI.get(Ipopt.Optimizer(), MOI.ListOfSupportedNonlinearOperators())

includes inv, so I thought it should handle the inverse of a matrix. however, in my minimal example, it does not:

	n = 5
	
	model = Model(Ipopt.Optimizer)
	
	#=
	decision variables
	=#
	#  the transition matrix
	@variable(model, P[1:n, 1:n] .≥ 0)
	#  the stationary distn
	@variable(model, π★[1:n] .≥ 0.0)

	# the fundamental matrix
	Z = @expression(model, inv(I(n) - (P - ones(n) * π★')))

giving an error MethodError: no method matching JuMP.NonlinearExpr(::JuMP.NonlinearExpr)

I saw your trick here to include Z as a variable with a constraint, but this causes the solver not to converge. is this the only solution?

@variable(model, Z[1:n, 1:n])
@constraint(model, Z * (I(n) - (P - ones(n) * π★')) .== I(n))

I’ve actually had success with inv by registering the function myself, via the legacy way. see below. so, it should be possible in principle, no? note: I’m trying to simplify my code but also add another variable, which is proving difficult to do by modifying my idleness_cost function.

	function idleness_cost(P_flat...)
		# rebuild transition matrix
		P = reshape(collect(P_flat), n, n)
		
		# compute the fundamenatal matrix
		Z = inv(I(n) - (P - ones(n) * π★'))
		
		# first hitting times starting at arbitrary state to all other states
		m₁ⱼ = [((1 == j ? 1.0 : 0.0) - Z[1, j] + Z[j, j]) / π★[j] for j = 1:n]

		# return cost of idleness costs
		return sum(m₁ⱼ[j] * ϕ[j] for j = 1:n)
	end
	
	model = Model(Ipopt.Optimizer)
	
	register(model, :idleness_cost, n * n, idleness_cost, autodiff=true)

thanks!

2

OR, a follow-up, this would help me too. how can I make π★ be an additional variable?!

	function idleness_cost(π★, P_flat...)
		# rebuild transition matrix
		P = reshape(collect(P_flat), n, n)
		
		# compute the fundamenatal matrix
		Z = inv(I(n) - (P - ones(n) * π★'))
		
		# first hitting times starting at arbitrary state to all other states
		m₁ⱼ = [((1 == j ? 1.0 : 0.0) - Z[1, j] + Z[j, j]) / π★[j] for j = 1:n]

		# return cost of idleness costs
		return sum(m₁ⱼ[j] * ϕ[j] for j = 1:n)
	end
	
	model = Model(Ipopt.Optimizer)

	#=
	decision variables
	=#
	#  the transition matrix
	@variable(model, P[1:n, 1:n] .≥ 0)
	#  the stationary distn
	@variable(model, π★[1:n] .≥ 0.0)
         
        register(model, :idleness_cost, n + n * n, idleness_cost, autodiff=true)

        @NLobjective(model, Min, idleness_cost(π★, P...))

I get error Unexpected array JuMP.VariableRef[π★[1], π★[2], π★[3], π★[4], π★[5], π★[6], π★[7], π★[8], π★[9], π★[10], π★[11], π★[12]] in nonlinear expression. Nonlinear expressions may contain only scalar expressions.

3

So you can do this:

julia> using JuMP, LinearAlgebra, Ipopt

julia> n = 5
5

julia> function idleness_cost(args...)
           ϕ = ones(n)
           π★, P = collect(args[1:n]), reshape(collect(args[n+1:end]), n, n)
           Z = inv(LinearAlgebra.I(n) - (P - ones(n) * π★'))
           return sum(
               ((1 == j ? 1 : 0) - Z[1, j] + Z[j, j]) / π★[j] * ϕ[j]
               for j in 1:n
           )
       end
idleness_cost (generic function with 1 method)

julia> begin
           model = Model(Ipopt.Optimizer)
           @variable(model, P[1:n, 1:n] >= 0)
           @variable(model, π★[1:n] >= 0.0)
           @operator(model, op_idleness, n^2+n, idleness_cost)
           @objective(model, Min, op_idleness(π★..., P...))
           optimize!(model)
       end
This is Ipopt version 3.14.16, running with linear solver MUMPS 5.7.3.

Number of nonzeros in equality constraint Jacobian...:        0
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:        0

Total number of variables............................:       30
                     variables with only lower bounds:       30
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        0
        inequality constraints with only lower bounds:        0
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls

It works but it’s not great. It relies on ForwardDiff.jl being able to differentiate through inv(::Matrix).

I’ll make a better error message for inv(::Matrix) in JuMP. We support inv but only for scalars.

The @NL` interface is legacy. You get the last error because you’re passing in a vector as the first argument to a user-defined function.

2 Likes