@rtapia if you are related to Richard Tapia from Rice University, you should ask him about this.
@ctkelley funny, I forget what I was reading but not too long ago I learned of his name.
I have noticed that every nonconvex solver requires upper and lower bounds on all decision variables, which is not inherently required by JuMP/Ipopt. Can this be manually specified using the set_solver_attribute? Wondering if this would make a difference.
If you are referring to Nonconvex.jl, you can use -Inf
and Inf
to signal no lower or upper bound.
Hello @mohamed82008, referring to JuMP in general. I am trying to understand why it cannot solve the problem above, can it solve nonconvex problems given the right solver? Nonconvex.jl worked.
That is inaccurate. Ipopt does not assume that anything is convex. It is simply a local optimizer. As such, it only looks for stationary points.
@rtapia It’s only sheer luck that you obtained the right solution with Ipopt via Nonconvex (if that gave you a different answer than JuMP/Ipopt, it must be because the starting point was different and/or some Ipopt options were different). Most nonconvex solvers are local, which means that they look for stationary points, not minimizers (not even local minimizers). Every normalized eigenvector of P is a stationary point of your problem. The Lagrange multiplier is the associated eigenvalue.
What you’re looking for is a global optimization solver.
KNITRO has a global optimization mode, but it’s commercial software. A cheap way to imitate what it does is to run Ipopt many times from many different starting points. Still, there’s no theoretical guarantee that you’ll find a global minimum, but you’ll increase your chances in probability.
@dpo not at all, the matrices are random and it repeatedly provides the minimum. While there may be multiple stationary points, the gradient is going to continue to point towards a descent direction. This is not a difficult problem to solve at all, I have a small primal-dual algorithm that can do the same, so can Casadi/Ipopt.
Jump/Ipopt states that the local optimal solution is at an infeasible point.
I just realised 2 potential improvements. First, you probably want this to be an equality constraint so add_eq_constraint!
. Second, you can use f1
directly instead of x -> f1(x)
.
@rtapia What I mean is that, even though IPOPT identifies the right eigenvalue in this particular case, the problem is nonconvex (because of the constraints), and there is nothing in its convergence analysis that guarantees that it will. It may be something specific to this problem.
Note that IPOPT does not follow the gradient direction. And it’s not because the negative gradient is a descent direction (regardless of the problem) that you will identify a global (or even local) minimum. Regardless, it’s insufficient to only consider the gradient in the presence of constraints.
I thought I would also show how simple it is to model and solve this problem with the JSO infrastructure:
using LinearAlgebra
using ADNLPModels, NLPModelsIpopt
function build_problem(n = 5)
A = 10 * rand(Float64, (n, n))
P = A' * A
λmin = eigmin(P)
f(x) = dot(x, P * x) / 2
c(x) = [dot(x, x)] # multiplier = ± λmin = 2f value
lcon = [1.]
ucon = [1.]
x0 = 2 * rand(n) .- 1
λmin, ADNLPModel(f, x0, c, lcon, ucon)
end
λmin, model = build_problem()
stats = ipopt(model)
@show stats.objective
@show stats.multipliers
@show stats.solution
Note that, in this particular problem, there are no other local minima (assuming distinct eigenvalues, all the other stationary points are saddle points + one maximum). In principle, local optimizers like Ipopt can converge to stationary points that are not local minima, but in practice IIRC this is pretty unlikely for generic starting points.
That’s why people can use Rayleigh-quotient optimization to find extremal eigenvalues, as in the LOBPCG algorithm — for it to eventually converge to the minimal eigenvalue, with many algorithms it is sufficient for the starting guess to have a nonzero inner product with the corresponding eigenvector (which is true almost certainly for a random starting vector … even if you choose a special starting guess for which it is not true, roundoff errors typically end up pushing you in the right direction).
Would also add that since eigenvectors are directions in the eigenspace, then one would think a Newton step would lead toward the lowest eigenvalue.
I am aware and I noted above that all eigenvectors are stationary points. All I’m saying is that you cannot count on local solvers like IPOPT to provide a local, let alone global, minimum of a nonconvex problem in general. There are initial guesses that will lead to a saddle-point. In this example, they may form a set of measure zero. It’s enough to look at a Newton step on the quadratic -x²/2 + y²/2 starting from (0, 1). Randomization may likely push you away from saddle points, but the convergence properties of solvers like Ipopt do not talk about that. It’s not because a problem has a single local min and other stationary points are saddle points that a local solver is guaranteed to find the local minimum. At any rate, if Ipopt is good enough for you on your problem, great.
That’s my point — I don’t generally worry much about things that happen with probability zero.
The probability is only zero in exact arithmetic, though.
I don’t worry too much about things with probability 10^{-16}, either.
(Especially since once you get into floating-point arithmetic it can help too, via rounding errors that push you off the unstable manifold leading to a saddle point.)
Certainly it seems like hyperbole to say things like “It’s only sheer luck that you obtained the right solution”.