Hi,
I have a simple problem in which I want to
\max_{x \in \mathbb{R}^{N}} f(x) \text{s.t}\hspace{5pt}x \geq 0
where f(x) is continuous and globally concave. Which is the best solver for nonlinear, globally concave maximization problems? I do not have constraints, besides non-negativity constraints. The only problems are that
(1) N is very large (N should be around 20,000)
(2) The solution might have a large number of zeros.
I have analytical gradients and Hessians!
The reason speed is of first order importance for me is that I have to optimize the function many many many times for different configurations of parameters.
Also, maybe I do not know, what the definition of sparsity is in this context. In the optimum there will be a large amount of 0s.
What is the recommendation?
I am adding a MWE below.
As you can see, I am maximizing the following function
\max_{x_{ij}} \sum_{j=1}^{J} (\sum_{i=1}^{I}x_{ij})^{\frac{p-1}{p}} - \sum_{i=1}^{I}(\sum_{j=1}^{J}t_{ij}x_{ij})^{\frac{1}{k}}
with p>1 and k<1. Hence, 0<(p-1)/p<1 and 1/k>0. The derivative of the typical element x_{ij} is
\frac{\partial f}{\partial x_{ij}} = \frac{(p-1)}{p}(\sum_{i=1}^{I}x_{ij})^{-\frac{1}{p}}-\frac{t_{ij}}{k}(\sum_{j=1}^{J}t_{ij}x_{ij})^{\frac{1-k}{k}}
Moreover, notice I define two other objects:
X_{j} = \sum_{i=1}^{I} x_{ij}
Y_{i} = \sum_{j=1}^{J} t_{ij}x_{ij}
Clearly, it is easy to compute the Hessian, which will be 0 for many entries. Although, it is going to be a pretty large matrix.
\frac{\partial^{2}f}{\partial x_{ij}^{2}} = -\frac{(p-1)}{p^{2}} (\sum_{i=1}^{I}x_{ij})^{-\frac{1-p}{p}} -\frac{(1-k)}{k^{2}}t_{ij}^{2}(\sum_{j=1}^{J}t_{ij}x_{ij})^{\frac{1-2k}{k}}
\frac{\partial^{2}f}{\partial x_{ij}x_{i^{\prime}j^{\prime}}} = 0
\frac{\partial^{2}f}{\partial x_{ij} x_{ij^{\prime}} } = -\frac{(1-k)}{k^{2}}t_{ij}t_{ij^{\prime}}(\sum_{j=1}^{J}t_{ij}x_{ij})^{\frac{1-2k}{k}}
\frac{\partial^{2}f}{\partial x_{ij}x_{i^{\prime}j}} = -\frac{(p-1)}{p^{2}} (\sum_{i=1}^{I}x_{ij})^{-\frac{1-p}{p}}
using JuMP, Ipopt, LinearAlgebra, Random, Gurobi, BenchmarkTools, MosekTools, NLopt
const Ι = 250 #This is Greek Letter Iota
const J = 50
const p = 5
const k = 0.75
function distmat(J,Ι)
Distances = zeros(J,Ι)
Random.seed!(7777)
coordinates_x = hcat([x for x = 1:J],fill(0.0,J))
coordinates_y = hcat([x for x = 1:Ι].*J/Ι,fill(0.0,Ι))
for j = 1:J, l=1:Ι
Distances[j,l] = sqrt((coordinates_x[j,1]-coordinates_y[l,1])^2+(coordinates_x[j,2]-coordinates_y[l,2])^2)
end
return 1 .+ .1*Distances
end
const t = distmat(J,Ι)
function solution_numerical_Ipopt(p,k,t)
(J,Ι) = size(t)
primal = Model(Ipopt.Optimizer)
set_silent(primal)
set_optimizer_attribute(primal, "nlp_scaling_method", "none")
@variable(primal, x[1:J,1:Ι] >= 0)
@NLobjective(primal, Max,
sum( sum(x[j,i] for i = 1:Ι)
^((p-1)/p) for j = 1:J) - sum( sum( t[j,i]*x[j,i] for j =1:J )^(1/k) for i = 1:Ι) )
JuMP.optimize!(primal)
println(objective_value(primal))
sol = value.(x)
return sol
end
sol_Ipopt = @time solution_numerical_Ipopt(p,k,t) ;
sum(sol_Ipopt,dims=2)
function func_grad(x,p,k,t)
(J,Ι) = size(t)
X = sum(sol_Ipopt,dims=2)
Y = sum(t.*sol_Ipopt,dims=1)'
res = sum(X.^((p-1)/p)) - sum(Y.^(1/k))
grad = zeros(J,Ι)
for i= 1:Ι
for j = 1:J
grad[j,i] = (p-1)/p * X[j]^(-1/p) - t[j,i]/k * Y[i]^((1-k)/k)
end
end
return res, grad
end
(f,g) = func_grad(sol_Ipopt,p,k,t)
