# Solving MCP with PathSolver with indexing

**URL:** <https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488>\
**Category:** Optimization (Mathematical)\
**Tags:** complementarity\
**Created:** [January 7, 2025, 4:47am UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488 "2025-01-07T04:47:45Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![jmc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmc/32/52649_2.png) [@jmc](https://discourse.julialang.org/u/jmc)\
**Post date:** [January 7, 2025, 4:47am UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/1 "2025-01-07T04:47:45Z")

</div>

Hello,

I am trying to solve an MCP that requires performing operations across some indices.

Here is a MWE:

```julia
using JuMP
using Random
using LinearAlgebra
using BenchmarkTools

using JuMP
import PATHSolver
PATHSolver.c_api_License_SetString("2830898829&Courtesy&&&USR&45321&5_1_2021&1000&PATH&GEN&31_12_2025&0_0_0&6000&0_0")

struct MyArrays
    denom_i::Vector{Float64}
    v::Vector{Float64}
    u::Vector{Float64}
end

function make_arrays(S, J)
    MyArrays(zeros(S), zeros(S), zeros(J))
end

function example_one_go(J::Int, S::Int, A::Matrix, Guess_y::Vector{Float64}, z::Vector{Float64}, C::Matrix, Guess::Matrix, myarrays::MyArrays)
    model = Model(PATHSolver.Optimizer)
    set_silent(model)

    @variable(model, x[i = 1:J, j = 1:J] >= 0, start = Guess[i, j])
    y_fixed_1 = 1.0
    @variable(model, y[i = 2:J] >= 0, start = Guess_y[i])

    last_x = Vector{Float64}()

    function f(x_i, x_j, x...)
        i = round(Int, x_i)
        j = round(Int, x_j)
        new_x_vec = collect(x)

        if last_x === nothing || new_x_vec != last_x
            X = reshape(new_x_vec, J, J)
            denom_i = myarrays.denom_i
            @views mul!(denom_i, A, X[i, :])
            last_x = x
        end

        return sum(A[s, j] / denom_i[s] for s in 1:S) * (1 / S)
    end

    function ∇f(g, x_i, x_j, x...)
        local i = round(Int, x_i)
        local j = round(Int, x_j)
        new_x_vec = collect(x)

        if last_x === nothing || new_x_vec != last_x
            X = reshape(new_x_vec, J, J)
            fill!(g, 0.0)
            denom_i = myarrays.denom_i
            @views mul!(denom_i, A, X[i, :])
            v = myarrays.v
            u = myarrays.u

            @inbounds @simd for s in 1:S
                v[s] = A[s, j] / (denom_i[s]^2)
            end

            @views mul!(u, transpose(A), v)
            for q in 1:J
                local col = 2 + (q - 1) * J + i
                g[col] = -(1 / S) * u[q]
            end

            last_x = x
        end
        return
    end

    get_y(i) = (i == 1) ? y_fixed_1 : y[i]

    @operator(model, op_f, 2 + J * J, f, ∇f)
    @constraint(model, comp1[i in 1:J, j in 1:J], -get_y(i) * op_f(i, j, x...) + get_y(j) * C[i, j] / z[j] ⟂ x[i, j])
    @constraint(model, comp2[i in 2:J], -sum(x[j, i] * C[i, j] for j in 1:J) + z[i] ⟂ y[i])

    JuMP.optimize!(model)

    Y = vcat(1.0, collect(value.(y)))
    return value.(x), Y, ((Y ./ z) .* C) .* value.(x)
end

# Example inputs to run the function
J = 10
S = 100
Guess_y = ones(J)
z = normalize(rand(J), 1) .+ 0.1 # Normalize z to have unit sum
C = ones(J, J)
for i in 1:J
    for j in 1:J
        if i != j
            C[i, j] = 1.25
        end
    end
end
Guess = ones(J, J)
A = rand(S,J) .+ 0.1 # Normalize and shift values to avoid zeros

myarrays = make_arrays(S, J)

# Running the function
result_x, result_Y, transformed_matrix = example_one_go(J, S, A, Guess_y, z, C, Guess, myarrays)

```

I have a couple of questions:

1. Are there any performance gains to be achieved in this code? In my true application, both J and S are in the thousands.

2. Is this code amenable to using other MCP solvers such as GAMS or KNitro? Or the syntax/JuMP structure does not translate well?

3. Is there any way that is useful for performance in which instead of looping over i,j, I just do everything once filling up a large (possible SParse) matrix, and then call that matrix to evaluate the constraints?

---

<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:** [January 7, 2025, 7:04am UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/2 "2025-01-07T07:04:10Z")

</div>

I’d write your model as:

```Julia
using JuMP
import LinearAlgebra
import PATHSolver

function example_one_go(
    J::Int,
    S::Int,
    A::Matrix,
    Guess_y::Vector{Float64},
    z::Vector{Float64},
    C::Matrix,
    Guess::Matrix,
)
    model = Model(PATHSolver.Optimizer)
    @variables(model, begin
        x[i = 1:J, j = 1:J] >= 0, (start = Guess[i, j])
        y[i = 1:J] >= 0, (start = Guess_y[i])
        denom[i in 1:J, s in 1:S] >= 1e-4, (start = A[s, :]' * Guess[i, :])
    end)
    fix(y[1], 1.0; force = true)
    @constraints(model, begin
        [i in 1:J, s in 1:S], A[s, :]' * x[i, :] - denom[i, s] ⟂ denom[i, s]
        [i in 1:J, j in 1:J],
            -y[i] * sum(A[s, j] / denom[i, s] for s in 1:S) / S + y[j] * C[i, j] / z[j] ⟂ x[i, j]
        [i in 2:J], -C[i, :]' * x[:, i] + z[i] ⟂ y[i]
    end)
    optimize!(model)
    @assert is_solved_and_feasible(model)
    X, Y = value.(x), value.(y)
    return X, Y, Y ./ z .* C .* X
end

# Example inputs to run the function
J, S = 10, 100
A = rand(S, J) .+ 0.1
Guess_y = ones(J)
z = LinearAlgebra.normalize(rand(J), 1) .+ 0.1
C = [i == j ? 1.0 : 1.25 for i in 1:J, j in 1:J]
Guess = ones(J, J)
@time ret = example_one_go(J, S, A, Guess_y, z, C, Guess)

```

1. Unless you have other things going on, using a user-defined operator is more hassle than it’s worth. I assume you went that way because JuMP’s default AD system doesn’t exploit common subexpressions in the denominator, so I’ve used a trick of adding a variable to represent the expression. This also lets you add bounds on the denominator which fixes a bunch of iterations like this in PATH:

```julia
    2 1 15 3 nan 1.0e+00 SD nan 5.6e-16 (comp1[1,2)

```

1. I don’t have a GAMS license so I can’t comment. In theory KNITRO supports this, but it really supports only `0 <= x ⟂ y >= 0` where `x` and `y` are variables, so JuMP will add extra variables and nonlinear equality constraints. This might defeat the benefit of using KNITRO.

2. No. Looping over indices is fine.

---

<div class="post-metadata">

**Author:** ![jmc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmc/32/52649_2.png) [@jmc](https://discourse.julialang.org/u/jmc)\
**Post date:** [January 7, 2025, 1:37pm UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/3 "2025-01-07T13:37:46Z")

</div>

Hi @odow , thanks for this suggestion.

I am implementing your suggestion with the function `example_one_go_v2`, in the following MWE:

```julia
using JuMP
using Random
using LinearAlgebra
using BenchmarkTools

using JuMP
import PATHSolver
PATHSolver.c_api_License_SetString("2830898829&Courtesy&&&USR&45321&5_1_2021&1000&PATH&GEN&31_12_2025&0_0_0&6000&0_0")

struct MyArrays
    denom_i::Vector{Float64}
    v::Vector{Float64}
    u::Vector{Float64}
end

function make_arrays(S, J)
    MyArrays(zeros(S), zeros(S), zeros(J))
end

function example_one_go(J::Int, S::Int, A::Matrix, Guess_y::Vector{Float64}, z::Vector{Float64}, C::Matrix, Guess::Matrix, myarrays::MyArrays)
    model = Model(PATHSolver.Optimizer)
    set_silent(model)

    @variable(model, x[i = 1:J, j = 1:J] >= 0, start = Guess[i, j])
    y_fixed_1 = 1.0
    @variable(model, y[i = 2:J] >= 0, start = Guess_y[i])

    last_x = Vector{Float64}()

    function f(x_i, x_j, x...)
        i = round(Int, x_i)
        j = round(Int, x_j)
        new_x_vec = collect(x)

        if last_x === nothing || new_x_vec != last_x
            X = reshape(new_x_vec, J, J)
            denom_i = myarrays.denom_i
            @views mul!(denom_i, A, X[i, :])
            last_x = x
        end

        return sum(A[s, j] / denom_i[s] for s in 1:S) * (1 / S)
    end

    function ∇f(g, x_i, x_j, x...)
        local i = round(Int, x_i)
        local j = round(Int, x_j)
        new_x_vec = collect(x)

        if last_x === nothing || new_x_vec != last_x
            X = reshape(new_x_vec, J, J)
            fill!(g, 0.0)
            denom_i = myarrays.denom_i
            @views mul!(denom_i, A, X[i, :])
            v = myarrays.v
            u = myarrays.u

            @inbounds @simd for s in 1:S
                v[s] = A[s, j] / (denom_i[s]^2)
            end

            @views mul!(u, transpose(A), v)
            for q in 1:J
                local col = 2 + (q - 1) * J + i
                g[col] = -(1 / S) * u[q]
            end

            last_x = x
        end
        return
    end

    get_y(i) = (i == 1) ? y_fixed_1 : y[i]

    @operator(model, op_f, 2 + J * J, f, ∇f)
    @constraint(model, comp1[i in 1:J, j in 1:J], -get_y(i) * op_f(i, j, x...) + get_y(j) * C[i, j] / z[j] ⟂ x[i, j])
    @constraint(model, comp2[i in 2:J], -sum(x[j, i] * C[i, j] for j in 1:J) + z[i] ⟂ y[i])

    JuMP.optimize!(model)

    Y = vcat(1.0, collect(value.(y)))
    return value.(x), Y, ((Y ./ z) .* C) .* value.(x)
end

function example_one_go_v2(
    J::Int,
    S::Int,
    A::Matrix,
    Guess_y::Vector{Float64},
    z::Vector{Float64},
    C::Matrix,
    Guess::Matrix,)
    model = Model(PATHSolver.Optimizer)
    set_silent(model)
    @variables(model, begin
        x[i = 1:J, j = 1:J] >= 0, (start = Guess[i, j])
        y[i = 1:J] >= 0, (start = Guess_y[i])
        denom[i in 1:J, s in 1:S] >= 1e-8, (start = A[s, :]' * Guess[i, :])
    end)
    fix(y[1], 1.0; force = true)
    @constraints(model, begin
        [i in 1:J, s in 1:S], A[s, :]' * x[i, :] - denom[i, s] ⟂ denom[i, s]
        [i in 1:J, j in 1:J],
            -y[i] * sum(A[s, j] / denom[i, s] for s in 1:S) / S + y[j] * C[i, j] / z[j] ⟂ x[i, j]
        [i in 2:J], -C[i, :]' * x[:, i] + z[i] ⟂ y[i]
    end)
    optimize!(model)
    #@assert is_solved_and_feasible(model)
    X, Y = value.(x), value.(y)
    return X, Y, Y ./ z .* C .* X
end

# Example inputs to run the function
J = 20
S = 1000
Guess_y = ones(J)
z = normalize(rand(J), 1) .+ 0.1 # Normalize z to have unit sum

x_coord = [x for x = 1 :J]
y_coord = ones(J)

function distance1(elas,x,y)
    dimJ = size(x)[1]
    DistMat = zeros(size(x)[1],size(x)[1])
        for i=1:dimJ, j=1:dimJ
            DistMat[i,j] = (1+sqrt((x[i]-x[j])^2+(y[i]-y[j])^2))^elas
        end    
    return DistMat    
end

C = distance1(0.1,x_coord,y_coord)

Guess = ones(J, J)
A = rand(S,J) .+ 0.1 # Normalize and shift values to avoid zeros
myarrays = make_arrays(S, J)

# Running the function
ret_v1 = @btime example_one_go(J, S, A, Guess_y, z, C, Guess, myarrays)

ret_v2 = @btime example_one_go_v2(J, S, A, Guess_y, z, C, Guess)

isapprox(ret_v1[2],ret_v2[2])

isapprox(ret_v1[1],ret_v2[1])

```

and at least for these values, the original function is much faster, and also more “robust” which is something that really surprised me, because I agree that iterations in which the denominator is 0 are troublesome.

In particular, if I run with J = 10,

```julia
# Running the function
ret_v1 = @btime example_one_go($J, $S, $A, $Guess_y, $z, $C, $Guess, $myarrays)

ret_v2 = @btime example_one_go_odow($J, $S, $A, $Guess_y, $z, $C, $Guess)

isapprox(ret_v1[2],ret_v2[2])

isapprox(ret_v1[1],ret_v2[1])

```

I get:

```julia
764.341 ms (45807012 allocations: 883.98 MiB)
3.741 s (5831544 allocations: 414.33 MiB)
true
true

```

with J = 15,

```julia
 1.769 s (98376055 allocations: 2.01 GiB)
14.397 s (12009241 allocations: 896.22 MiB)
true
false

```

but this `false` is noise (difference of 1e-8).

But when I do J = 20, for example, then the second code crashes.

```julia
 2.891 s (155912056 allocations: 3.63 GiB)
 55.848 s (20321527 allocations: 1.43 GiB)

```

and when I uncomment the assert the function fails.

I do like that the function you propose makes much less allocations and uses much less data. This might be very important for scaling up purposes. But it seems this is coming at the expense of performance and reliability. This to me is very odd, because your proposed code is very transparent.

Am I missing something, or I messed up the function you proposed in a way I cannot find.

---

<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:** [January 7, 2025, 10:40pm UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/4 "2025-01-07T22:40:00Z")

</div>

Not sure of the reason. Have you looked at the log from PATH? Is the time being spent in JuMP? Or is PATH taking many more iterations to converge? You might consider making the lower bound of the denominator larger.

If you’re going to manually code the gradient, you should almost just use `PATHSolver.solve_mcp` directly:

> <https://github.com/chkwon/PATHSolver.jl/blob/db02fc9f4285f1dec3f3db4ba4faed3acfd53552/src/C_API.jl#L576-L707>

---

<div class="post-metadata">

**Author:** ![jmc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmc/32/52649_2.png) [@jmc](https://discourse.julialang.org/u/jmc)\
**Post date:** [January 7, 2025, 11:51pm UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/5 "2025-01-07T23:51:36Z")

</div>

Thanks, the logs are, for J = 10, S=1000

For `example_one_go`

```julia
Path 5.0.03 (Fri Jun 26 10:05:33 2020)
Written by Todd Munson, Steven Dirkse, Youngdae Kim, and Michael Ferris

Crash Log
major func diff size residual step prox (label)
    0 0 1.2615e+01 0.0e+00 (comp2[2)
    1 1 17 109 6.2502e+00 1.0e+00 0.0e+00 (comp2[2)
pn_search terminated: no progress.

Major Iteration Log
major minor func grad residual step type prox inorm (label)
    0 0 14 2 6.2502e+00 I 0.0e+00 5.3e+00 (comp2[2)
    1 13 15 3 4.7591e-02 1.0e+00 SM 0.0e+00 2.7e-02 (comp1[1,6)
    2 119 16 4 -nan(ind) 1.0e+00 SD 0.0e+00 -nan(ind) (comp1[1,1)
    3 46 17 5 -nan(ind) 1.0e+00 SD -nan(ind) -nan(ind) (comp1[1,1)
    4 11 18 6 -nan(ind) 1.0e+00 SD -nan(ind) -nan(ind) (comp1[1,1)
    5 3 19 7 -nan(ind) 1.0e+00 SD -nan(ind) -nan(ind) (comp1[1,1)
    6 1 20 8 -nan(ind) 1.0e+00 SD -nan(ind) -nan(ind) (comp1[1,1)
    7 1 21 9 -nan(ind) 1.0e+00 SD -nan(ind) -nan(ind) (comp1[1,1)
    8 1 22 10 -nan(ind) 1.0e+00 SD -nan(ind) -nan(ind) (comp1[1,1)
    9 1 23 11 -nan(ind) 1.0e+00 SD -nan(ind) -nan(ind) (comp1[1,1)
   10 1 24 12 -nan(ind) 1.0e+00 SD -nan(ind) -nan(ind) (comp1[1,1)
   11 1 26 13 8.1305e-01 8.0e-01 SW 0.0e+00 2.9e-01 (comp1[6,2)
   12 66 27 14 2.1668e-01 1.0e+00 SO 0.0e+00 7.2e-02 (comp1[6,2)
   13 26 28 15 9.5901e-02 1.0e+00 SO 0.0e+00 4.0e-02 (comp1[9,7)
   14 6 29 16 2.0485e-02 1.0e+00 SO 0.0e+00 8.0e-03 (comp1[9,7)
   15 7 30 17 1.4439e-02 1.0e+00 SO 0.0e+00 1.0e-02 (comp1[1,4)
   16 2 31 18 8.4694e-04 1.0e+00 SO 0.0e+00 4.8e-04 (comp1[1,4)
   17 1 32 19 5.1479e-06 1.0e+00 SO 0.0e+00 3.6e-06 (comp1[1,3)
   18 1 33 20 3.0525e-10 1.0e+00 SO 0.0e+00 2.2e-10 (comp1[1,3)

Major Iterations. . . . 18
Minor Iterations. . . . 307
Restarts. . . . . . . . 0
Crash Iterations. . . . 1
Gradient Steps. . . . . 0
Function Evaluations. . 33
Gradient Evaluations. . 20
Basis Time. . . . . . . 0.000000
Total Time. . . . . . . 0.375000
Residual. . . . . . . . 3.052528e-10
  0.389121 seconds (24.48 M allocations: 475.030 MiB, 8.65% gc time)

```

and for `example_one_go_v2`

```julia
Path 5.0.03 (Fri Jun 26 10:05:33 2020)
Written by Todd Munson, Steven Dirkse, Youngdae Kim, and Michael Ferris

Preprocessed size : 10109

Crash Log
major func diff size residual step prox (label)
    0 0 1.2615e+01 0.0e+00 ()
    1 12 0 10109 1.2614e+01 5.4e-05 0.0e+00 ()
pn_search terminated: no basis change.

Major Iteration Log
major minor func grad residual step type prox inorm (label)
    0 0 13 2 1.2614e+01 I 0.0e+00 1.2e+01 ()
    1 114 14 3 1.5193e-01 1.0e+00 SM 0.0e+00 5.2e-02 ()
    2 17 15 4 7.5712e-02 1.0e+00 SO 0.0e+00 2.6e-02 ()
    3 121 17 5 6.4499e+00 8.0e-01 CB 0.0e+00 6.4e+00 ()
    4 335 19 6 1.3900e+00 8.0e-01 CB 0.0e+00 1.4e+00 ()
    5 161 20 7 1.4274e+00 1.0e+00 CO 0.0e+00 1.4e+00 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 1.8984e-03 10019 50 0 39 9 w[52] z[101]
    6 16 21 8 1.0775e+02 1.0e+00 SD 1.4e-01 9.8e+00 ()
    7 205 22 9 1.0768e+02 1.0e+00 CD 5.7e-02 9.8e+00 ()
    8 264 23 10 1.0768e+02 1.0e+00 CD 2.3e-02 9.8e+00 ()
    9 272 24 11 1.0768e+02 1.0e+00 CD 9.1e-03 9.8e+00 ()
   10 114 25 12 1.0768e+02 1.0e+00 CD 3.7e-03 9.8e+00 ()
   11 308 26 13 1.0768e+02 1.0e+00 CD 1.5e-03 9.8e+00 ()
   12 54 28 14 7.5303e-02 1.0e+00 RG 0.0e+00 2.6e-02 ()
   13 14 29 15 3.7618e-02 1.0e+00 SO 0.0e+00 1.3e-02 ()
   14 403 31 16 7.4297e+00 8.0e-01 CB 0.0e+00 7.3e+00 ()
   15 171 33 17 4.8264e+00 8.0e-01 CB 0.0e+00 4.6e+00 ()
   16 318 34 18 5.6495e+01 1.0e+00 CD 0.0e+00 6.1e+00 ()
   17 85 35 19 2.3467e+02 1.0e+00 CD 0.0e+00 1.1e+02 ()
   18 114 36 20 2.3405e+02 1.0e+00 CD 0.0e+00 1.1e+02 ()
   19 23 38 21 3.7416e-02 1.0e+00 RG 0.0e+00 1.3e-02 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -1.7129e-04 10055 53 0 0 12 w[46] z[56]

Major Iteration Log
major minor func grad residual step type prox inorm (label)
   20 503 45 22 6.2490e+00 1.3e-01 CB 0.0e+00 6.2e+00 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -1.0592e-03 10053 55 0 0 13 w[36] z[26]

Major Iteration Log
major minor func grad residual step type prox inorm (label)
   21 944 46 23 4.8430e+00 1.0e+00 CO 0.0e+00 4.8e+00 ()
   22 23 47 24 4.8430e+00 1.0e+00 CO 0.0e+00 4.8e+00 ()
   23 33 48 25 4.8430e+00 1.0e+00 CO 0.0e+00 4.8e+00 ()
   24 38 63 26 4.8430e+00 3.5e-08 RB 0.0e+00 4.8e+00 ()
   25 24 75 27 3.7404e-02 3.2e-04 CG 0.0e+00 1.3e-02 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -2.0783e-03 10046 62 0 0 12 w[62] z[82]

Major Iteration Log
major minor func grad residual step type prox inorm (label)
   26 610 76 28 2.9577e-02 1.0e+00 CO 0.0e+00 1.0e-02 ()
   27 275 84 29 2.8827e+00 5.5e-02 CB 0.0e+00 2.9e+00 ()
   28 94 107 30 2.9575e-02 2.8e-02 CG 0.0e+00 1.0e-02 ()
   29 38 133 31 2.9573e-02 3.2e-04 RG 0.0e+00 1.0e-02 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -2.9705e-02 10050 58 0 0 35 z[61] w[81]
 1000 -3.6487e-02 10038 70 0 0 71 w[41] w[12]

Restart Log
proximal_perturbation 0
crash_method none
crash_perturb yes
nms_initial_reference_factor 2
lemke_start_type slack
proximal_perturbation 1.0000e-01

Major Iteration Log
major minor func grad residual step type prox inorm (label)
   29 0 160 32 1.2615e+01 R 1.0e-01 1.2e+01 ()
   30 13 161 33 9.4936e+01 1.0e+00 SD 4.0e-02 1.7e+00 ()
   31 13 162 34 3.6983e+00 1.0e+00 SO 1.6e-02 1.5e-01 ()
   32 11 163 35 2.4337e+00 1.0e+00 SO 6.4e-03 6.2e-02 ()
   33 7 164 36 4.0659e-01 1.0e+00 SO 2.6e-03 1.4e-02 ()
   34 6 165 37 1.9530e-01 1.0e+00 SO 1.0e-03 8.2e-03 ()
   35 30 166 38 1.6286e-01 1.0e+00 SO 4.1e-04 2.5e-02 ()
   36 7 167 39 2.9250e-02 1.0e+00 SO 1.6e-04 2.9e-03 ()
   37 2 168 40 3.4643e-03 1.0e+00 SO 6.6e-05 2.5e-04 ()
   38 1 169 41 2.3248e-04 1.0e+00 SO 2.6e-05 1.2e-05 ()
   39 1 170 42 6.7830e-06 1.0e+00 SO 1.0e-05 3.5e-07 ()
   40 1 171 43 1.3597e-08 1.0e+00 SO 6.8e-07 7.0e-10 ()

Major Iterations. . . . 40
Minor Iterations. . . . 7289
Restarts. . . . . . . . 1
Crash Iterations. . . . 1
Gradient Steps. . . . . 6
Function Evaluations. . 171
Gradient Evaluations. . 43
Basis Time. . . . . . . 18.594000
Total Time. . . . . . . 20.593000
Residual. . . . . . . . 1.359654e-08
Postsolved residual: 1.3597e-08
 21.027442 seconds (5.83 M allocations: 414.477 MiB, 0.80% gc time)

```

---

<div class="post-metadata">

**Author:** ![jmc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmc/32/52649_2.png) [@jmc](https://discourse.julialang.org/u/jmc)\
**Post date:** [January 7, 2025, 11:55pm UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/6 "2025-01-07T23:55:24Z")

</div>

as of using the solve\_mcp directly I am little bit lost, do you know where can I find a MWE, or how would you adjust this example?

---

<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:** [January 8, 2025, 12:58am UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/7 "2025-01-08T00:58:16Z")

</div>

See [PATHSolver.jl/README.md at master · chkwon/PATHSolver.jl · GitHub](https://github.com/chkwon/PATHSolver.jl/blob/master/README.md#use-with-the-c-api)

But now that I’ve looked at your example a little more deeply, I’m not sure that it is correct, because it doesn’t account for the different `i` and `j` at the same `x`.

How about something like:

```Julia
function example_one_go(
    J::Int,
    S::Int,
    A::Matrix,
    Guess_y::Vector{Float64},
    z::Vector{Float64},
    C::Matrix,
    Guess::Matrix,
)
    denom_i, v, last = zeros(S), zeros(S), nothing
    function update_cache(i, j, x)
        if (i, j, x) != last
            denom_i .= A * collect(x[i:J:end])
            for s in 1:S
                v[s] = A[s, j] / denom_i[s]^2
            end
            last = (i, j, x)
        end
        return
    end
    function f(x_i, x_j, x...)
        i, j = round(Int, x_i), round(Int, x_j)
        update_cache(i, j, x)
        return sum(A[s, j] / denom_i[s] for s in 1:S) / S
    end
    function ∇f(g, x_i, x_j, x...)
        i, j = round(Int, x_i), round(Int, x_j)
        update_cache(i, j, x)
        fill!(g, 0.0)
        for q in 1:J
            g[2+J*(q-1)+i] = -sum(A[s, q] * v[s] for s in 1:S) / S
        end
        return
    end
    model = Model(PATHSolver.Optimizer)
    @operator(model, op_f, 2 + J * J, f, ∇f)
    @variable(model, x[i = 1:J, j = 1:J] >= 0, start = Guess[i, j])
    @variable(model, y[i = 1:J] >= 0, start = Guess_y[i])
    fix(y[1], 1.0; force = true)
    @constraints(model, begin
        [i in 1:J, j in 1:J], -y[i] * op_f(i, j, x...) + y[j] * C[i, j] / z[j] ⟂ x[i, j]
        # 0.0 ⟂ y[1] # Added by default
        [i in 2:J], -sum(x[j, i] * C[i, j] for j in 1:J) + z[i] ⟂ y[i]
    end)
    optimize!(model)
    @assert is_solved_and_feasible(model)
    X, Y = value.(x), value.(y)
    return X, Y, Y ./ z .* C .* X
end

```

---

<div class="post-metadata">

**Author:** ![jmc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmc/32/52649_2.png) [@jmc](https://discourse.julialang.org/u/jmc)\
**Post date:** [January 8, 2025, 1:19am UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/8 "2025-01-08T01:19:17Z")

</div>

Can you clarify what do you mean by “it doesn’t account for the different `i` and `j` at the same `x”?`

---

<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:** [January 8, 2025, 1:53am UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/9 "2025-01-08T01:53:19Z")

</div>

Your code checks only if `new_x_vec != last_x`, but it might be called with the same value of `x` but a different `i`. In which case, `denom_i` isn’t going to be correct.

You don’t hit the bug because `last_x === nothing || new_x_vec != last_x` is always `true` (check it) because you have `last_x = x` instead of `last_x = new_x_vec`, so your clever caching trick isn’t actually doing anything.

---

<div class="post-metadata">

**Author:** ![jmc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmc/32/52649_2.png) [@jmc](https://discourse.julialang.org/u/jmc)\
**Post date:** [January 8, 2025, 3:21am UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/10 "2025-01-08T03:21:25Z")

</div>

Ah! Great catch. Thanks for this! If I may abuse your time, can I ask four questions?

1. `if (i, j, x) != last` in this line why are you including the `j`. I mean, in your code I understand why you are doing this, because of the `v`. But could it be faster to just update `denom_i` in the `update_cache` which only would be triggered by different `i` and different `x`.
2. Does it make sense to clamp `denom_i`? Or that usually hurts performance?
3. Do you think that using `solve_mcp` will provide large performance improvements?
4. Are there any limits to how many `J` or `S` I can use with `PATH`? Or the limits are given by my machine’s specs?

Again, thanks!

---

<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:** [January 8, 2025, 3:58am UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/11 "2025-01-08T03:58:40Z")

</div>

1. Sure, you could be cleverer with the caching
2. I don’t understand. Clamp to what values?
3. Not really
4. It’s problem dependent. You won’t know until you try

---

<div class="post-metadata">

**Author:** ![jmc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmc/32/52649_2.png) [@jmc](https://discourse.julialang.org/u/jmc)\
**Post date:** [January 13, 2025, 9:11pm UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/12 "2025-01-13T21:11:46Z")

</div>

I am trying to implement the a restriction, such that for `x[i,:]>0`, which is a sufficient condition for `denom[i,s]>0`, but I am struggling to implement it. This is the code I am using:

```julia
function example_one_go(
    J::Int,
    S::Int,
    A::Matrix,
    Guess_y::Vector{Float64},
    z::Vector{Float64},
    C::Matrix,
    Guess::Matrix,
)
    denom_i, v, last = zeros(S), zeros(S), nothing
    function update_cache(i, j, x)
        if (i, j, x) != last
            if sum(collect(x[i:J:end])) == 0.0
            println(sum(collect(x[i:J:end])))
            end
            denom_i .= A * collect(x[i:J:end])
            for s in 1:S
                v[s] = A[s, j] / denom_i[s]^2
            end
            last = (i, j, x)
        end
        return
    end
    function f(x_i, x_j, x...)
        i, j = round(Int, x_i), round(Int, x_j)
        update_cache(i, j, x)
        return sum(A[s, j] / denom_i[s] for s in 1:S) / S
    end
    function ∇f(g, x_i, x_j, x...)
        i, j = round(Int, x_i), round(Int, x_j)
        update_cache(i, j, x)
        fill!(g, 0.0)
        for q in 1:J
            g[2+J*(q-1)+i] = -sum(A[s, q] * v[s] for s in 1:S) / S
        end
        return
    end
    model = Model(PATHSolver.Optimizer)
    set_optimizer_attribute(model, "time_limit", 36000)
    set_optimizer_attribute(model, "major_iteration_limit", 1000)
    set_optimizer_attribute(model, "minor_iteration_limit", 20000)

    @operator(model, op_f, 2 + J * J, f, ∇f)
    @variable(model, x[i = 1:J, j = 1:J] >= 0, start = Guess[i, j])
    @variable(model, y[i = 1:J] >= 0, start = Guess_y[i])
    @variable(model, λ[i=1:J]>= 0, start = zeros(J)[i]) 
    fix(y[1], 1.0; force = true)
    @constraints(model, begin
        [i in 1:J, j in 1:J], -y[i] * op_f(i, j, x...) + y[j] * C[i, j] / z[j] ⟂ x[i, j]
        # 0.0 ⟂ y[1] # Added by default
        [i in 2:J], -sum(x[j, i] * C[i, j] for j in 1:J) + z[i] ⟂ y[i]
        [i in 1:J], 1e-6 - sum(x[i,j] for j in 1:J) ⟂ λ[i]
    end)
    optimize!(model)
    @assert is_solved_and_feasible(model)
    X, Y = value.(x), value.(y)
    return X, Y, Y ./ z .* C .* X
end

```

But I am still getting instances in which `x[i,:]==0`, and therefore I have a bunch of iterations that become `NaN`, i.e, the `if` that prints the output is triggered. Do you know what might be happening?

It also happens if I set a much larger threshold for the `sum(x[i,j] for j in 1:J)`.

Relatedly, there are a number of cases in which I am getting

`Termination status: ITERATION_LIMIT`

with

```julia
Major Iterations. . . . 83
Minor Iterations. . . . 10365

```

As if the option was not being changed.

---

<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:** [January 13, 2025, 11:47pm UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/13 "2025-01-13T23:47:47Z")

</div>

> [@jmc](#):
>
> `[i in 1:J], 1e-6 - sum(x[i,j] for j in 1:J) ⟂ λ[i]`

Do you mean instead: `[i in 1:J], sum(x[i,j] for j in 1:J) - 1e-6 ⟂ λ[i]`?

From the [definition of mixed-complementarity](https://jump.dev/JuMP.jl/stable/manual/constraints/#Complementarity-constraints): if `λ[i]` is at it’s lower bound (here 0), then `sum(x[i,j] for j in 1:J) - 1e-6 >= 0`, otherwise, `sum(x[i,j] for j in 1:J) - 1e-6 = 0`.

---

<div class="post-metadata">

**Author:** ![jmc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmc/32/52649_2.png) [@jmc](https://discourse.julialang.org/u/jmc)\
**Post date:** [January 14, 2025, 3:17am UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/14 "2025-01-14T03:17:00Z")

</div>

Yes, sorry my bad.

```julia
function example_one_go(
    J::Int,
    S::Int,
    A::Matrix,
    Guess_y::Vector{Float64},
    z::Vector{Float64},
    C::Matrix,
    Guess::Matrix,)
    denom_i, v, last = zeros(S), zeros(S), nothing
    function update_cache(i, j, x)
        if (i, j, x) != last
            if sum(collect(x[i:J:end])) == 0.0
                println(sum(collect(x[i:J:end])))
            end
            denom_i .= A * collect(x[i:J:end])
            for s in 1:S
                v[s] = A[s, j] / denom_i[s]^2
            end
            last = (i, j, x)
        end
        return
    end
    function f(x_i, x_j, x...)
        i, j = round(Int, x_i), round(Int, x_j)
        update_cache(i, j, x)
        return sum(A[s, j] / denom_i[s] for s in 1:S) / S
    end
    function ∇f(g, x_i, x_j, x...)
        i, j = round(Int, x_i), round(Int, x_j)
        update_cache(i, j, x)
        fill!(g, 0.0)
        for q in 1:J
            g[2+J*(q-1)+i] = -sum(A[s, q] * v[s] for s in 1:S) / S
        end
        return
    end
    model = Model(PATHSolver.Optimizer)
    set_silent(model)
    set_optimizer_attribute(model, "time_limit", 36000)
    set_optimizer_attribute(model, "major_iteration_limit", 1000)
    set_optimizer_attribute(model, "minor_iteration_limit", 20000)

    @operator(model, op_f, 2 + J * J, f, ∇f)
    @variable(model, x[i = 1:J, j = 1:J] >= 0, start = Guess[i, j])
    @variable(model, y[i = 1:J] >= 0, start = Guess_y[i])
    @variable(model, λ[i=1:J]>= 0, start = zeros(J)[i]) 
    fix(y[1], 1.0; force = true)
    @constraints(model, begin
        [i in 1:J, j in 1:J], -y[i] * op_f(i, j, x...) + y[j] * C[i, j] / z[j] ⟂ x[i, j]
        # 0.0 ⟂ y[1] # Added by default
        [i in 2:J], -sum(x[j, i] * C[i, j] for j in 1:J) + z[i] ⟂ y[i]
        [i in 1:J], sum(x[i,j] for j in 1:J) - 1e-6 ⟂ λ[i]
    end)
    optimize!(model)
    println(termination_status(model))
    #@assert is_solved_and_feasible(model)
    X, Y = value.(x), value.(y)
    return X, Y, Y ./ z .* C .* X
end

```

This code still prints 0s for

```julia
if sum(collect(x[i:J:end])) == 0.0
                println(sum(collect(x[i:J:end])))
            end

```

For example, for J = 10,

```julia
> time_limit 36000
 > minor_iteration_limit 20000
 > major_iteration_limit 1000
Read of options file complete.

Path 5.0.03 (Fri Jun 26 10:05:33 2020)
Written by Todd Munson, Steven Dirkse, Youngdae Kim, and Michael Ferris

Preprocessed size : 119

Crash Log
major func diff size residual step prox (label)
    0 0 3.1214e+00 0.0e+00 ()
    1 4 0 109 1.7606e+00 4.1e-01 0.0e+00 ()
pn_search terminated: no basis change.

Major Iteration Log
major minor func grad residual step type prox inorm (label)
    0 0 5 2 1.7606e+00 I 0.0e+00 5.9e-01 ()
    1 97 6 3 2.1569e-02 1.0e+00 SM 0.0e+00 7.9e-03 ()
    2 45 11 4 1.1930e-02 3.3e-01 RB 0.0e+00 4.4e-03 ()
    3 58 12 5 1.1227e-02 1.0e+00 RO 0.0e+00 4.1e-03 ()
    4 27 13 6 2.3131e+01 1.0e+00 RD 0.0e+00 2.3e+01 ()
    5 17 14 7 2.2672e+01 1.0e+00 RD 0.0e+00 2.3e+01 ()
    6 20 15 8 2.2672e+01 1.0e+00 RD 0.0e+00 2.3e+01 ()
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
    7 44 17 9 1.1226e-02 1.0e+00 RG 0.0e+00 4.1e-03 ()
    8 114 18 10 1.1015e-02 1.0e+00 RO 0.0e+00 4.0e-03 ()
    9 66 19 11 8.3437e-03 1.0e+00 RO 0.0e+00 3.1e-03 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 1.4033e-02 67 51 0 0 9 w[19] z[49]
 1000 4.1136e-05 58 60 0 0 19 z[58] w[48]
 1500 1.6093e-03 55 63 0 0 29 w[78] z[58]

Major Iteration Log
major minor func grad residual step type prox inorm (label)
   10 1596 20 12 8.6068e-03 1.0e+00 CO 8.3e-04 3.0e-03 ()
   11 311 21 13 8.6068e-03 1.0e+00 CO 3.3e-04 3.0e-03 ()
   12 151 22 14 6.5863e-03 1.0e+00 RO 1.3e-04 2.3e-03 ()
   13 75 27 15 6.5527e-03 4.1e-01 RG 5.3e-05 2.3e-03 ()
   14 68 28 16 6.5730e-03 1.0e+00 CO 2.1e-05 2.3e-03 ()
   15 52 29 17 6.5675e-03 1.0e+00 CO 6.6e-04 2.3e-03 ()
   16 178 30 18 6.5731e-03 1.0e+00 CO 2.6e-04 2.3e-03 ()
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
   17 126 31 19 1.1182e-02 1.0e+00 RO 1.1e-04 2.7e-03 ()
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
   18 35 32 20 1.1182e-02 1.0e+00 RD 1.1e-03 2.7e-03 ()
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
   19 165 40 21 6.4795e-03 6.9e-02 CG 5.3e-05 2.3e-03 ()
   20 307 41 22 6.5354e-03 1.0e+00 CO 2.1e-05 2.3e-03 ()
   21 67 42 23 6.5353e-03 1.0e+00 CO 6.5e-04 2.3e-03 ()
   22 164 43 24 8.1009e-03 1.0e+00 RO 2.6e-04 2.1e-03 ()
   23 53 44 25 8.1009e-03 1.0e+00 RO 1.0e-04 2.1e-03 ()
   24 171 45 26 7.1599e-03 1.0e+00 RO 4.2e-05 2.1e-03 ()
   25 91 49 27 6.4647e-03 6.4e-01 RG 5.3e-05 2.3e-03 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -5.6081e-01 58 60 0 0 9 w[23] z[53]

Major Iteration Log
major minor func grad residual step type prox inorm (label)
   26 700 50 28 6.4794e-03 1.0e+00 CO 2.1e-05 2.3e-03 ()
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
   27 67 51 29 -nan(ind) 1.0e+00 RL 8.5e-06 1.9e-03 ()
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
   28 3 52 30 -nan(ind) 1.0e+00 ND -nan(ind) 1.7e-03 ()
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
   29 2 53 31 -nan(ind) 1.0e+00 ND -nan(ind) 1.7e-03 ()
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
   30 2 54 32 -nan(ind) 1.0e+00 ND -nan(ind) 1.7e-03 ()
   31 2 62 33 6.4180e-03 6.9e-02 RG 5.3e-05 2.3e-03 ()
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
   32 74 63 34 -nan(ind) 1.0e+00 RL 2.1e-05 1.9e-03 ()
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
   33 3 64 35 -nan(ind) 1.0e+00 ND -nan(ind) 1.8e-03 ()
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
   34 2 65 36 -nan(ind) 1.0e+00 ND -nan(ind) 1.7e-03 ()
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
   35 2 66 37 -nan(ind) 1.0e+00 ND -nan(ind) 1.7e-03 ()
   36 2 79 38 6.4717e-03 4.4e-05 RW 5.3e-05 2.3e-03 ()
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0
0.0

Restart Log
proximal_perturbation 0
crash_method none
crash_perturb yes
nms_initial_reference_factor 2
lemke_start_type slack
proximal_perturbation 3.1214e-02

Major Iteration Log
major minor func grad residual step type prox inorm (label)
   36 0 87 39 3.1214e+00 R 3.1e-02 9.9e-01 ()
   37 10 88 40 5.8700e-01 1.0e+00 SO 1.2e-02 1.5e-01 ()
   38 10 89 41 2.8095e-01 1.0e+00 SO 5.0e-03 8.2e-02 ()
   39 1 90 42 1.1860e-01 1.0e+00 SO 2.0e-03 4.1e-02 ()
   40 2 91 43 5.8734e-02 1.0e+00 SO 8.0e-04 2.0e-02 ()
   41 6 92 44 1.4371e-02 1.0e+00 SO 3.2e-04 5.5e-03 ()
   42 22 93 45 7.2756e-03 1.0e+00 SO 1.3e-04 2.6e-03 ()
   43 32 94 46 3.1172e-03 1.0e+00 SO 5.1e-05 1.1e-03 ()
   44 12 95 47 1.4272e-03 1.0e+00 SO 2.0e-05 5.8e-04 ()
   45 6 96 48 6.0114e-04 1.0e+00 SO 8.2e-06 3.4e-04 ()
   46 3 97 49 6.8275e-05 1.0e+00 SO 3.3e-06 3.3e-05 ()
   47 2 98 50 3.2237e-06 1.0e+00 SO 1.3e-06 1.8e-06 ()
   48 3 99 51 1.6754e-08 1.0e+00 SO 3.2e-07 9.8e-09 ()

Major Iterations. . . . 48
Minor Iterations. . . . 5487
Restarts. . . . . . . . 1
Crash Iterations. . . . 1
Gradient Steps. . . . . 6
Function Evaluations. . 99
Gradient Evaluations. . 51
Basis Time. . . . . . . 0.078000
Total Time. . . . . . . 2.250000
Residual. . . . . . . . 1.675384e-08
Postsolved residual: 1.6754e-08
LOCALLY_SOLVED

```

and when the problem is large (e.g. for J==33), it hits an iteration limit, that does not correspond with the one I set.

I get the following optimizer output:

```julia
 > time_limit 36000
 > minor_iteration_limit 20000
 > major_iteration_limit 1000
Read of options file complete.

Path 5.0.03 (Fri Jun 26 10:05:33 2020)
Written by Todd Munson, Steven Dirkse, Youngdae Kim, and Michael Ferris

Preprocessed size : 1154

Crash Log
major func diff size residual step prox (label)
    0 0 5.6310e+00 0.0e+00 ()
    1 8 0 1121 5.5478e+00 2.8e-02 0.0e+00 ()
pn_search terminated: no basis change.

Major Iteration Log
major minor func grad residual step type prox inorm (label)
    0 0 9 2 5.5478e+00 I 0.0e+00 9.8e-01 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 4.7949e-02 621 532 0 0 9 w[212] z[247]
 1000 1.2867e-01 384 769 0 0 19 z[639] w[417]
 1500 9.4313e-01 818 335 0 0 29 z[531] w[714]

Major Iteration Log
major minor func grad residual step type prox inorm (label)
    1 1740 10 3 7.6207e-02 1.0e+00 SM 0.0e+00 1.5e-02 ()
    2 1 11 4 3.7881e-02 1.0e+00 SO 0.0e+00 7.4e-03 ()
    3 1 12 5 1.8930e-02 1.0e+00 SO 0.0e+00 3.7e-03 ()
    4 6 13 6 8.1340e+01 1.0e+00 RD 0.0e+00 7.7e+01 ()
    5 6 14 7 7.5953e+01 1.0e+00 SD 8.1e+00 7.2e+01 ()
    6 1 15 8 6.1566e+01 1.0e+00 SD 3.3e+00 5.8e+01 ()
    7 1 16 9 2.8252e+01 1.0e+00 SD 1.3e+00 2.6e+01 ()
    8 18 17 10 5.0035e+00 1.0e+00 SO 5.2e-01 4.2e+00 ()
    9 3 18 11 1.9058e-02 1.0e+00 SO 2.1e-01 1.3e-02 ()
   10 3 19 12 7.4784e-03 1.0e+00 SO 1.9e-03 1.5e-03 ()
   11 6 20 13 6.9905e-03 1.0e+00 SO 7.5e-04 1.4e-03 ()
   12 14 21 14 6.1666e-03 1.0e+00 SO 3.0e-04 1.2e-03 ()
   13 15 22 15 6.0588e-03 1.0e+00 SO 6.2e-04 1.1e-03 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -8.0448e+03 635 518 0 0 9 w[137] z[175]

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -3.2248e+04 639 514 0 0 9 w[236] z[203]

Major Iteration Log
major minor func grad residual step type prox inorm (label)
   14 507 23 16 6.1056e-03 1.0e+00 NO 6.1e-04 1.0e-03 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -1.8464e+03 635 518 0 0 9 w[261] z[129]

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -1.0094e+04 635 518 0 0 9 w[261] z[129]

Major Iteration Log
major minor func grad residual step type prox inorm (label)
   15 683 24 17 6.1060e-03 1.0e+00 NO 6.1e-04 1.0e-03 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -1.8701e+03 634 519 0 0 9 w[129] z[558]

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -1.0099e+04 634 519 0 0 9 w[129] z[558]

Major Iteration Log
major minor func grad residual step type prox inorm (label)
   16 696 25 18 6.1060e-03 1.0e+00 NO 6.1e-04 1.0e-03 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -1.8925e+03 633 520 0 0 9 w[558] z[492]

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -3.1079e+01 578 575 0 0 9 z[165] w[1093]

Major Iteration Log
major minor func grad residual step type prox inorm (label)
   17 817 26 19 6.1060e-03 1.0e+00 RO 6.1e-04 1.0e-03 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -1.8925e+03 633 520 0 0 9 w[558] z[492]

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -3.1079e+01 578 575 0 0 9 z[165] w[1093]

Major Iteration Log
major minor func grad residual step type prox inorm (label)
   18 817 27 20 6.1060e-03 1.0e+00 NO 6.1e-04 1.0e-03 ()

Minor Iteration Log
minor t z w v art ckpts enter leave
  500 -3.3776e-04 499 654 0 0 9 z[714] w[814]
 1000 -1.0582e-04 999 154 0 0 19 z[96] w[784]
 1500 -3.9670e-03 650 503 0 0 29 w[317] z[917]
 2000 -1.9719e-02 194 959 0 0 39 w[478] z[412]

Major Iteration Log
major minor func grad residual step type prox inorm (label)
   19 2137 28 21 5.1411e-03 1.0e+00 RO 2.4e-04 9.9e-04 ()

Major Iterations. . . . 19
Minor Iterations. . . . 11178
Restarts. . . . . . . . 0
Crash Iterations. . . . 1
Gradient Steps. . . . . 0
Function Evaluations. . 28
Gradient Evaluations. . 21
Basis Time. . . . . . . 1.813000
Total Time. . . . . . . 22.235000
Residual. . . . . . . . 5.141060e-03
Postsolved residual: 5.1411e-03
ITERATION_LIMIT
 33.628792 seconds (213.57 M allocations: 11.120 GiB, 2.70% gc time, 31.52% compilation time)

```

---

<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:** [January 14, 2025, 3:24am UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/15 "2025-01-14T03:24:56Z")

</div>

Re the zeros: It’s probable that PATH visits solutions that slightly violate complementarity during the solution process. It might be trying those solutions frequently because that’s the lower bound of `x`, so it starts there and tries searching in a positive direction. I don’t know the details of the algorithm implemented in PATH.

Re the iteration limit: PATH is terminating because of `cumulative_iteration_limit`. Here are the options for PATH:

> **[options.pdf](https://pages.cs.wisc.edu/~ferris/path/options.pdf)**
>
> 16.79 KB

---

<div class="post-metadata">

**Author:** ![jmc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmc/32/52649_2.png) [@jmc](https://discourse.julialang.org/u/jmc)\
**Post date:** [January 14, 2025, 3:32am UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/16 "2025-01-14T03:32:14Z")

</div>

Thanks for both answers. Super helpful. I clearly missed that option.

Yes, that could be the case. I find surprising that if I “ask” `sum(collect(x[i:J:end]))` to be above 1, it still visits that area

---

<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:** [January 14, 2025, 4:06am UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/17 "2025-01-14T04:06:03Z")

</div>

Many solution algorithms visit primal infeasible points on their way to an optimal solution. For example, the dual simplex algorithm starts with an infeasible point and iterates through infeasible points, and as soon as it finds a primal feasible point, that is the optimal solution. You should not expect solvers to visit only feasible primal points.

---

<div class="post-metadata">

**Author:** ![jmc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmc/32/52649_2.png) [@jmc](https://discourse.julialang.org/u/jmc)\
**Post date:** [January 14, 2025, 1:21pm UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/18 "2025-01-14T13:21:09Z")

</div>

Makes sense! Thanks.

From your experience, does setting a floor to `denom_i[s]`, i.e.,

```julia
denom_i .= A * collect(x[i:J:end])
            for s in 1:S
                v[s] = A[s, j] / maximum(denom_i[s],epsilon)^2
            end

```

Is a good idea? Or will this mostly confuse the solver?

---

<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:** [January 14, 2025, 7:30pm UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/19 "2025-01-14T19:30:52Z")

</div>

I don’t have much experience with the “tricks” of MCP modeling. I would suggest that it is a bad idea because it probably means that your Jacobian does not match the function.

---

<div class="post-metadata">

**Author:** ![jmc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmc/32/52649_2.png) [@jmc](https://discourse.julialang.org/u/jmc)\
**Post date:** [January 16, 2025, 11:41am UTC](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488/20 "2025-01-16T11:41:03Z")

</div>

I am getting this error message for large `J`:

```julia
Internal error: encountered unexpected error during compilation of NonlinearOperator:
StackOverflowError()
unknown function (ip: 00007ffc42ef4006)

```

Is this something I can solve using HPC rather than a laptop? Is this related to memory? I understand that the dimension of the problem increases to the order of `J^2`, and having large `J` can be too demanding, but can I address this issue with a better computer?

[Next page](https://discourse.julialang.org/t/solving-mcp-with-pathsolver-with-indexing/124488.md?page=2)