I'm looking to find a computationally efficient solution to a large system of nonlinear equations. I'm trying to maximize the following function:
$$ f(\vec{x}) = \sum_i^N C_i (x_i-A_i)x_i^{\epsilon_{ii}} \prod_{j\not = i}^N x_j^{\epsilon_{ij}}$$
Subject to the constraints: $$ C_i > 0 \quad\forall\: i$$ $$ x_i > A_i > 0 \quad\forall\: i$$ $$ {\epsilon}_{ii} < -1 \quad\forall\: i$$ $$ {\epsilon}_{ij} \ge 0 \quad\forall\: i,j \quad\text{s.t.}\; i \not= j$$
I've been able to write a few solvers that can solve this, but I'm struggling to find/write a solver that can scale to $N=1000000$ (one million). Unfortunately, the $\epsilon_{ij}$ matrix is dense.
Things I've tried:
- scipy.optimize BFGS algorithm: takes 8 hours when N = 2000
- a fixed point solver: takes 3 hours for N = 10000, but seems to get stuck in local maxima
I would really appreciate advice on what to investigate next. I was gonna start going down the path of either (a) hand deriving the Jacobian and Hessian and taking another stab at using a scipy.optimize solver or (b) trying to coerce this to look more like a geometric programming problem and using gpkit. I'm pretty new to optimization so I could very well be missing something obvious - please don't hesitate to point out anything I've overlooked!
If it helps, I do have distributed computing resources available (spark, dask, kubernetes, etc).