Solving nonlinear constrained optimization problem

Specific Domains Optimization (Mathematical)
,
1

Hello,

I am trying to solve the following nonconvex problem in Julia using Optim.jl:
min x’Px
s.t: 1 -x’*x <=0

where P is a positive definite matrix. I have defined the following

using JuMP, Optim

n = 1500;
A = 10*rand(Float64,(n,n));
P = A'*A;

m = Model(Optim.Optimizer)
set_optimizer_attribute(m, "method", BFGS())

@variable(m, x[1:n]);
@objective(m,Min, x'*P*x);
f = 1 - x'*x;
@constraint(m, f1, f <= 0)

optimize!(m)

But I get a a stackoverflowerror, so I am not sure what I am missing. I know there are multiple ways of defining the problem in Optim. But the only constrained example I see is the boxconstraint: Interior point Newton - Optim.jl (julianlsolvers.github.io)

Thanks

2

Isn’t this analytically solvable? According to the min–max theorem, your minimum will be the smallest eigenvalue of P, achieved when x is a corresponding eigenvector (normalized to unit length).

If you are just using this as a test case for JuMP, I don’t think the BFGS algorithm in Optim.jl supports arbitrary nonlinear constraints? You may need some other engine, like Nonconvex.jl or NLopt.jl or Ipopt.jl or …

1 Like
3

Yes this is just a test case, according to Optim.jl optim can handle generic nonlinear constraints using IPNewton.

I have since corrected the code to:

set_optimizer_attribute(m, "method", IPNewton());
x0=ones(n,1);

@variable(m, x[1:n]);
#@variable(m, x[i=1:n], start = x0[i]);
set_start_value.(x,1);

I will try Nonconvex.jl, I have tried NLopt and Ipopt. Just getting familiar with the available packages in Julia.

4

Hello,

I have formed the following problem using Optim.jl with a nonconvex constraint:

min x’Px
1-x’*x <=0

According to: Optim.jl
Optim should be able to solve it using IPNewton?

the code is as follows:

using JuMP, Optim, LineSearches

n = 10;
A = 10*rand(Float64,(n,n));
P = A'*A;

# cost function
f0(x) = x'*P*x;

# gradient of cost function
function f0_grad!(gf0,x)
	for j=1:n
		gf0[j] = 2*P[j,:]'*x;
	end
end

# Hessian is a constant
function f0_hess!(h,x)
	for j=1:n
		for i=1:n
			h[i,j] = 2*P[i,j];
		end
	end
end
# nonlinear constraint
con_c!(c,x) = (c[1] = 1 - x'*x);

# Jacobian of constraint
function con_jacobian!(J,x)
	for j=1:n
		J[j] = -2*x[j];
	end
end

# Hessian of Constraint
function con_h!(h2,x,λ)
	#h = zeros(n,n);
        for j=1:n
			h2[j,j] = -2*λ[1];
		end
end

#upper and lower bounds
lx = -10*ones(n); ux = 10*ones(n);
lc = [-Inf]; uc = [-1];

#initial condition
x0 = 1.1*ones(n);

df = TwiceDifferentiable(f0, f0_grad!, f0_hess!, x0)	
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!, lx, ux, lc, uc)

res = optimize(df, dfc, x0, IPNewton())

The solution states “success” but the solution is incorrect. The solution of x should correspond to the eigenvector (normalized) associated with the lowest eigenvalue.

Thanks,

5

How far off is the minimum from the minimal eigenvalue? I notice that you don’t set a tolerance or any other convergence criteria — maybe you just need to run for more iterations / lower the tolerance?

6

Its far off, I increased the number iterations and lowered the tolerance, but it seems to make it worse. The value should be small, but I get a solution ~e3.

7

You could do:

using JuMP, Ipopt
function main(A)
    n = size(A, 1)
    model = Model(Ipopt.Optimizer)
    x0 = zeros(n)
    _, i = findmin(vec(sum(A; dims = 2)))
    x0[i] = 1.0
    y0 = A * x0
    @variable(model, x[i in 1:n], start = x0[i]);
    @variable(model, y[i in 1:n], start = y0[i]);
    @objective(model, Min, sum(y.^2))
    @constraint(model, A * x .== y)
    @constraint(model, 1 <= sum(x.^2))
    optimize!(model)
    @assert is_solved_and_feasible(model)
    return value.(x)
end

n = 1_500
A = 10 * rand(Float64, (n, n))
main(A)

But note that Ipopt assumes the problem is convex, so you are not guaranteed to find a global minima.

8

Well that is working much better, I modified it to use @NLconstraint, that should technically work. I am trying to see if increasing tolerance could have an impact on the solution, but I am not sure what is the correct way to implement.

using JuMP, Ipopt, LinearAlgebra

x1 = zeros(n);

function main(A)
    n = size(A, 1)
    model = Model(Ipopt.Optimizer)
    set_attribute(model,"max_iter", 10000)
    model.xtol_rel = 1e-6;
    set_optimizer_attribute(model,"acceptable_constr_viol_tol", 1e-10)
    set_optimizer_attribute(model,"constr_viol_tol", 1e-6)
    x0 = zeros(n)
    x0 = ones(n)
    y0 = A * x0
    @variable(model, x[i in 1:n], start = x0[i]);
    @variable(model, y[i in 1:n], start = y0[i]);
    @objective(model, Min, sum(y.^2))
    f = x'*x;
    @constraint(model, A * x .== y)
    @NLconstraint(model, 1 <= f)
    JuMP.optimize!(model)
    @assert is_solved_and_feasible(model)
    global x1 =value.(x);
    return JuMP.objective_value(model)
end

n = 10;
A = 10 * rand(Float64, (n, n));
A = A'*A;
main(A)

V = eigvals(A);
print("solution should be = ",V[1]);
9

Please don’t :laughing: the @NL interface is now legacy (see ANN: JuMP v1.15 is released) and it should not be encouraged for new uses. Also, it means the constraint is added as a generic nonlinear constraint. If you use @constraint JuMP and Ipopt can detect it as quadratic.

global x1 =value.(x);

Also, don’t use global. Return the value from the function with

    return value.(x), JuMP.objective_value(model)

and then call it as

x1, obj = main(A)

model.xtol_rel = 1e-6;

This I not valid JuMP syntax.

You can see the full list of Ipopt options here: Ipopt: Ipopt Options

The way to set an option is set_attribute(model, "tol", 1e-6).

10

Just a guess, but you are mainly tackling the Rayleigh quotient here (up to scaling, and your constraints).

Your constraints basically say you want t norm (squared) larger than 1. So the minimum will lie on the boundary (you could also take norme equal) to one.

If you try to solve the (unconstraint) Rayleigh problem with any Newton based iteration, one can compute that the update will boil down to multiplying your (nearly) eigenvector with a (positive) scalar. Might that be a cause here?

11

Clearly a noob here :rofl:, thanks for the corrections. I learned based off older versions.

@kellertuer I have solved this with a Primal-Dual interior point solver

2 Likes
12

for anyone interested I got the exact solution (or near exact) using Nonconvex.jl

using Nonconvex, NonconvexIpopt, LinearAlgebra, Printf

n = 10;
A = 10*rand(Float64,(n,n));
P = A'*A;


f(x) = x'*P*x;
f1(x) = 1- x'*x;
lx = -500*ones(n); ux = 500*ones(n);
model = Model(f)
addvar!(model, lx,ux)
add_ineq_constraint!(model, x -> f1(x))

x0 = 1.5*ones(n);
#alg = NLoptAlg(:LD_MMA)
alg = IpoptAlg()
#alg = MMA02()
options = IpoptOptions(tol=1e-10, first_order = false)
#options = MMAOptions()
res = optimize(model, alg, x0, options =options, convcriteria = KKTCriteria())

V = eigvals(P);
@printf("Answer should be = %0.4f\n",V[1])
@printf("Answer is = %0.4f\n", res.minimum)
1 Like
13

If you want something scalable, try using MadNLP.

using MadNLP, LinearAlgebra

struct MyProblem{T, VT, MT} <: MadNLP.AbstractNLPModel{T, VT}
    Q::MT
    counters::MadNLP.NLPModels.Counters
    meta::MadNLP.NLPModelMeta{T,VT}
end

MyProblem(Q) = MyProblem(
    Q,
    MadNLP.NLPModels.Counters(),
    MadNLP.NLPModelMeta(
        size(Q,1);
        x0 = fill!(similar(Q,size(Q,1)), 10),
        ncon = 1,
        ucon = fill!(similar(Q,1), 0)
    )
)

MadNLP.obj(m::MyProblem, x::AbstractVector) = dot(x,m.Q,x)
MadNLP.grad!(m::MyProblem, x::AbstractVector, y::AbstractVector) = (mul!(y,m.Q,x); y.*=2)
MadNLP.cons!(m::MyProblem, x::AbstractVector, y::AbstractVector) = fill!(y, 1 - dot(x,x))
MadNLP.jac_dense!(m::MyProblem, x::AbstractVector, j::AbstractMatrix) = (j .= x'; j.*=-2)
function MadNLP.hess_dense!(m::MyProblem, x::AbstractVector, y::AbstractVector, h::AbstractMatrix; obj_weight= one(eltype(x)))
    copyto!(h, m.Q)
    h .*= obj_weight
    view(h,1: n+1: n^2) .-= y[1]
end

n = 100;
A = 10*rand(Float64,(n,n));
P = A'*A;

m = MyProblem(P)
madnlp(m; lapack_algorithm=MadNLP.CHOLESKY, max_iter= 10000)


# got an Nvidia gpu?
using MadNLPGPU, CUDA

MadNLP.obj(m::MyProblem, x::CuVector) = dot(x,m.Q * x)
function MadNLP.hess_dense!(m::MyProblem, x::CuVector, y::CuVector, h::CuMatrix; obj_weight= one(eltype(x)))
    copyto!(h, m.Q)
    h .*= obj_weight
    CUDA.@allowscalar view(h,1: n+1: n^2) .-= y[1]
end

m = MyProblem(CuArray(P))
madnlp(m; lapack_algorithm=MadNLP.CHOLESKY, max_iter= 10000)
14

I’ve had issues with MadNLP, get the following error. When I search madnlp function it does not have the option to specify the algorithm.

image
image1265×808 39.5 KB

15

Could you try with a more recent version of MadNLP? This seems to be using v0.6.0 or older. The most recent version is v0.8.0.

16

I get the following error when trying to update:

image
image1223×1020 100 KB

17

Could you try to install MadNLP in a separate project environment? Most likely some other package is blocking the upgrade.

cd to/some/tmp/directory
julia --project=.
]add MadNLP
18

That worked, but the output is incorrect, answer should be the smallest eigenvalue of P.

I am a little confused on the line:
MadNLP.grad!(m::MyProblem, x::AbstractVector, y::AbstractVector) = (mul!(y,m.Q,x); y.*=2)
is this stating the gradient is [ Px; 2Px] ? or just 2Px?

19

It first put m.Q*x into y, and then multiply the entries of y by 2, which has the same effect of
y .= 2 .* m.Q *x, but just reduces memory allocation.

MadNLP is a local solver, like Ipopt and many other NLP solvers. So, it doesn’t aim to find the global minimum.

Just curious, what’s the motivation for solving this as an optimization problem rather than using more specialized methods (e.g., power iteration)?

20

Hello,

The objective value was smaller than the possible minimum. I am just testing, then I will scale to a larger optimization problem that is nonconvex. The two problems are not related but would like to seek a “reliable” solver.

Thanks