mcpsolve keeps getting stuck

[djsegal]

For example,

Iter     f(x) inf-norm    Step 2-norm
------   --------------   --------------
     0     3.324599e+00              NaN
     1     1.074657e+00     1.023475e+01
     2     1.074657e+00     8.011869e-31
     3     1.074657e+00     0.000000e+00
     4     1.074657e+00     1.232595e-32
     5     1.074657e+00     8.011869e-31
     6     1.074657e+00     0.000000e+00
     7     7.190746e-01     4.632532e+00
     8     7.190746e-01     5.320852e-01

---

Notice that 7 & 8 have the same inf-norm.

The code then stays there forever

[sglyon]

I've noticed similar things, but without code representing your objective function (and derivatives if you have them analytically) it will be hard for us to do any debugging.

Any chance you can supply that without too much effort and without it being too much code for us to look over?

[djsegal]

Is there a way I could decrease the scope of the problem to make it more manageable?

[djsegal]

@sglyon one question I had was, why do I need:

      if abs(rho) > 1
        cur_vars = [ 0 , 0 , 0 ]
        cur_F = [ 0 , 0 , 0 ]
        return
      end

in my equation_set! if I set the range of rho to [0,1] in:

  cur_solution = mcpsolve(
    @generate_wave_equation_set(cur_solved_R_0, cur_solved_B_0),
    [1.0, 0, 0.5], [3.0, 1.0, 3.0], [2.0, 0.5, 1.5], factor=3.0,
    ftol = 1e-3, xtol = 5e-4, show_trace = true
  ).zero

// i.e. why does rho sometimes escape its given bounds?

[sglyon]

Is there a way I could decrease the scope of the problem to make it more manageable?

@djsegal I'm not sure -- I was going to ask the same thing. while I'd love to help, I personally don't have the time to go through all the code there.

// i.e. why does rho sometimes escape its given bounds?

This is a result of how the mixed complementarity problem (mcp) is obtained from the underlying nonlinear system. Check out this file (https://github.com/JuliaNLSolvers/NLsolve.jl/blob/master/src/mcp.jl) for more details. The basic idea its that the original objective function (call it f) is transformed into a secondary function g that knows something about the bounds. In the minmax formulation the function g looks something like

 g(x) = min(max(f(x), x-upper), x-lower)

where upper and lower are the upper and lower bounds you supply. If you think carefully about this transformation you notice that x can be arbitrary, so your function f can be evaluated anywhere -- even outside the domain [lower, upper]. The min and max operations then make sure that the residual for g (that the actual solver is trying to find a zero for) takes into account the bounds.

Does that make sense?

[djsegal]

Yes,

I think one of the problems is I have this line of code:

cur_n_profile = ( 1 - rho ^ 2 ) ^ 1/2

which prevents real solutions when | rho | > 1

// i.e. i get the following error trace:

 in nan_dom_err at ./math.jl:196 [inlined]
 in ^(::Float64, ::Float64) at ./math.jl:355

---

edit: is there a way to push the solver off from bad corners once it finds one?

[djsegal]

update:

it's nothing to write home about, but this was my solution to the problem:

  cur_solution = nothing

  for cur_attempt in 1:20
    did_work = true

    try
      cur_initial_values = [
        rand(linspace(1.0, 3.0)),
        rand(linspace(0.0, 1.0)),
        rand(linspace(0.5, 3.0))
      ]

      cur_solution = mcpsolve(
        @generate_wave_equation_set(cur_solved_R_0, cur_solved_B_0),
        [0.0, 0.0, 0.0], [20.0, 1.0, 15.0], cur_initial_values, factor=5.0,
        show_trace = false, xtol = 1e-6, ftol = 5e-4, iterations=40
      ).zero

    catch
      did_work = false
    end

    if !did_work ; continue ; end

    print("✓")
    break
  end

  if cur_solution == nothing
    print("x")
    cur_solution = [ 0 , 0 , 0 ]
  end

[sglyon]

I wish there was a more satisfactory answer for you, but I think that as long as you've found something that works -- it might be the best you can do for now.

The main issue is that a mixed complementarity solver is not a constrained root finder. As long as the function f satisfies the sign conditions outside the bounds (see the readme), the final solution for x will be in the constraints you give the mcp solver. But, there is no guarantee that the solver will never evaluate the original objective outside the bounds.

[djsegal]

Ok. Closing

[sglyon]

I think the original issue about the mcpsolver getting stuck is still relevant.

I'll reopen until we have figured that one out. I think it might be related to the issues with trust region reported in #89