NLsolve cannot solve this linear complementarity problem

[chkwon]

I've tried to solve the following simple linear complementarity problem using NLsolve.

using NLsolve

M = [0  0 -1 -1 ;
     0  0  1 -2 ;
     1 -1  2 -2 ;
     1  2 -2  4 ]

q = [2; 2; -2; -6]

function f!(x, fvec)
    fvec = M * x + q
end

r = mcpsolve(f!, [0., 0., 0., 0.], [Inf, Inf, Inf, Inf],
             [1.25, 0., 0., 0.5], reformulation = :smooth, autodiff = true)

x = r.zero  # [1.25, 0.0, 0.0, 0.5]
@show dot( M*x + q, x )  # 0.5

sol = [2.8, 0.0, 0.8, 1.2]
@show dot( M*sol + q, sol )  # 0.0

While the solution is know to [2.8, 0.0, 0.8, 1.2], I obtained [1.25, 0.0, 0.0, 0.5], just same as the initial solution provided.

julia> @show r
r = Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [1.25,0.0,0.0,0.5]
 * Zero: [1.25,0.0,0.0,0.5]
 * Inf-norm of residuals: 0.000000
 * Iterations: 0
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 1
 * Jacobian Calls (df/dx): 1

It seems that it was not run properly with 0 iterations.

[sglyon]

Hmm this is odd. I don't have time to take a closer look right now, but thank you for reporting.

@KristofferC if you have any time to look into this that'd be great -- if not that's ok too!

[KristofferC]

Since I know nothing of this thing you call "complementary problems" I am not the right person to ask :)

All I can see is that it says |f(x)| < 1.0e-08: true so it considers the equation solved.

[pkofod]

@chkwon was this ever solved?

[chkwon]

@pkofod It seems not.

[pkofod]

I'm not at a computer but it seems as if f! Has a 🐜. You need to assign into fvec, copy! or similar. Fvec isn't updated the way you're doing it

[sglyon]

Good catch @pkofod.

Defining

function f!(x, fvec)
    copy!(fvec, M * x + q)
end

results in

julia> r = mcpsolve(f!, [0., 0., 0., 0.], [Inf, Inf, Inf, Inf],
                    [1.25, 0., 0., 0.5], reformulation = :smooth, autodiff = true)
Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [1.25,0.0,0.0,0.5]
 * Zero: [2.8,-4.72171e-16,0.8,1.2]
 * Inf-norm of residuals: 0.000000
 * Iterations: 6
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 7
 * Jacobian Calls (df/dx): 7

julia> x = r.zero  # [1.25, 0.0, 0.0, 0.5]
4-element Array{Float64,1}:
  2.8
 -4.72171e-16
  0.8
  1.2

julia> @show dot( M*x + q, x )  # 0.5
dot(M * x + q,x) = -4.468660707835697e-14
-4.468660707835697e-14

julia> sol = [2.8, 0.0, 0.8, 1.2]
4-element Array{Float64,1}:
 2.8
 0.0
 0.8
 1.2

julia> @show dot( M*sol + q, sol )  # 0.0
dot(M * sol + q,sol) = 3.552713678800501e-16
3.552713678800501e-16

[chkwon]

Brilliant! Thanks a lot!