# Handling nonlinear objective involving the inverse of a matrix

**URL:** <https://discourse.julialang.org/t/handling-nonlinear-objective-involving-the-inverse-of-a-matrix/123726>\
**Category:** Optimization (Mathematical)\
**Created:** [December 11, 2024, 8:29pm UTC](https://discourse.julialang.org/t/handling-nonlinear-objective-involving-the-inverse-of-a-matrix/123726 "2024-12-11T20:29:50Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Cory\_Simon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cory_simon/32/3950_2.png) [@Cory\_Simon](https://discourse.julialang.org/u/Cory_Simon)\
**Post date:** [December 11, 2024, 8:29pm UTC](https://discourse.julialang.org/t/handling-nonlinear-objective-involving-the-inverse-of-a-matrix/123726/1 "2024-12-11T20:29:50Z")

</div>

so the output of:

```julia
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:

```julia
	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](https://discourse.julialang.org/t/matrix-inverse-in-constraints-in-nlconstraint/88462/4) to include `Z` as a variable with a constraint, but this causes the solver not to converge. is this the only solution?

```julia
@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 `register`ing 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.

```julia
	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!

---

<div class="post-metadata">

**Author:** ![Cory\_Simon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cory_simon/32/3950_2.png) [@Cory\_Simon](https://discourse.julialang.org/u/Cory_Simon)\
**Post date:** [December 11, 2024, 8:35pm UTC](https://discourse.julialang.org/t/handling-nonlinear-objective-involving-the-inverse-of-a-matrix/123726/2 "2024-12-11T20:35:42Z")

</div>

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

```julia
	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.`

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [December 11, 2024, 8:55pm UTC](https://discourse.julialang.org/t/handling-nonlinear-objective-involving-the-inverse-of-a-matrix/123726/3 "2024-12-11T20:55:40Z")

</div>

So you can do this:

```Julia
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.