2
$\begingroup$

I have probability values $p \in \mathbb{R}^n$. Given $A \in \mathbb{R}^{m\times n}$, $b \in \mathbb{R}^m$, I want to minimize the following objective function. $||Ap - b||_2^2 + \sum_{i=1}^{n-2} (\log(p_i) - 2\log(p_{i + 1}) + \log(p_{i+2}))^2$

subject to $\sum p_i = 1$. Unfortunately, the second term is non-convex. I have tried to do change of variables $y_i = \exp(x_i)$ to make the second term convex, but then the first term is non-convex.

I have two questions.

  1. What is a good way to optimize this problem?
  2. I believe that this objective function doesn't have any local minimums besides the global minimum (perhaps with some weak assumptions). Is this true? If so, how could I solve this?
Improve this question
$\endgroup$
9
  • 1
    $\begingroup$ Are you assuming that the p_i values are sorted in increasing order? $\endgroup$
    – Kevin Dalmeijer
    Apr 22 at 4:09
  • 4
    $\begingroup$ How big is n? You could just try throw it in a local nonlinear optimizer and see what happens.. Don't forget the nonnegativity constraints on p; or better yet, $p_i \ge $ small_positive_number. If n is small enough, a (sort of) rigorous global optimization solver such as BARON might be able to solve to global optimality. Why do you think there are no non-globally optimal local minima? $\endgroup$
    – Mark L. Stone
    Apr 22 at 13:56
  • $\begingroup$ Are you interested in heuristic solutions (not guaranteed to find an optimum)? $\endgroup$
    – prubin ♦
    Apr 22 at 19:56
  • $\begingroup$ @KevinDalmeijer no I’m not assuming they are increasing in order $\endgroup$
    – JEK
    Apr 23 at 1:47
  • $\begingroup$ @MarkL.Stone I think using the log implies positive constraint on the probabilities. I’m guessing there is no non global local minimum because both first term and the second terms are convex (after change of variables). $\endgroup$
    – JEK
    Apr 23 at 1:50

2 Answers 2

Reset to default
2
$\begingroup$

You could just throw this at a global NLP solver. For "small" instances they can prove global optimality. I tried the following using GAMS/Baron:

Set
  i /i1*i25/
  j /j1*j25/
;

parameter A(i,j), b(i);
A(i,j) = uniform(0,1);
b(i) = uniform(0,1);

variable p(j),z;
p.lo(j) = 0.001;
p.up(j) = 1;

equations obj,e;

obj.. z =e= sum(i,sqr(sum(j,A(i,j)*p(j))-b(i))) +
            sum(j$(ord(j)<=card(j)-2), sqr(log(p(j))-2*log(p(j+1))+log(p(j+2))));
            
e.. sum(j,p(j)) =e= 1;

model m /all/;
option nlp=baron, threads=0;
solve m minimizing z using nlp;

This gives:

===========================================================================
 BARON version 23.3.11. Built: WIN-64 Sat Mar 11 18:01:20 EST 2023

 BARON is a product of The Optimization Firm.
 For information on BARON, see https://minlp.com/about-baron

 If you use this software, please cite publications from
 https://minlp.com/baron-publications, such as: 

 Khajavirad, A. and N. V. Sahinidis,
 A hybrid LP/NLP paradigm for global optimization relaxations,
 Mathematical Programming Computation, 10, 383-421, 2018.
===========================================================================
 This BARON run may utilize the following subsolver(s)
 For LP/MIP/QP: CLP/CBC, ILOG CPLEX                             
 For NLP: MINOS, SNOPT, External NLP, IPOPT, FILTERSQP
===========================================================================
 Doing local search
 Preprocessing found feasible solution with value 2.10908
 Solving bounding LP
 Starting multi-start local search
 Done with local search
===========================================================================
  Iteration    Open nodes         Time (s)    Lower bound      Upper bound
          1             1             0.59      1.72152          2.10908
       1526            83            29.75      2.00708          2.10908
       4691            73            59.41      2.00708          2.10908
       7903            74            89.16      2.00708          2.10908
      10797            62           118.78      2.00708          2.10908
      13761            51           148.61      2.00820          2.10908
      16724            49           178.17      2.00820          2.10908
      19677            41           207.70      2.00820          2.10908
      22105            36           237.20      2.00820          2.10908
      24746            26           267.09      2.00820          2.10908
      27238            18           296.91      2.00820          2.10908
      29698            17           326.59      2.04487          2.10908
      31844            10           356.42      2.04487          2.10908
      34535            14           386.00      2.05357          2.10908
      37169            10           415.53      2.06959          2.10908
      39835             5           445.19      2.07546          2.10908
      42543             7           474.88      2.08354          2.10908
      43439             0           484.69      2.10886          2.10908

 Calculating duals

                         *** Normal completion ***            

 Wall clock time:                   490.83
 Total CPU time used:               484.69

 Total no. of BaR iterations:   43439
 Best solution found at node:      -1
 Max. no. of nodes in memory:      88
 
 All done
===========================================================================

Solution      = 2.10907512234556  best solution found during preprocessing
Best possible = 2.10886421617
Absolute gap  = 0.00021090617555819  optca = 1E-9
Relative gap  = 9.99993662262895E-5  optcr = 0.0001

Baron finds the global optimum during preprocessing but proving global optimality takes a bit of time. This could be a confirmation that the local optimum is indeed globally optimal. The penalty term seems rather small (for my random data) compared to the LS term.

Note: you can use https://neos-server.org/neos/ to solve this. This can also be useful to verify results found by your own algorithms. NEOS gives you access to different global NLP solvers (both under AMPL and GAMS).

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ Thanks. I suggested this in my first comment, but was too lazy to try it out. $\endgroup$
    – Mark L. Stone
    Apr 24 at 12:16
  • 1
    $\begingroup$ "The penalty term seems rather small (for my random data) compared to the LS term." Given that there is no multiplier on the penalty term, its relative contribution is determined by the scaling of the input data. $\endgroup$
    – Mark L. Stone
    Apr 24 at 12:43
1
$\begingroup$

The following is a heuristic, meaning I cannot prove it converges to an optimum. It generates a sequence of solutions with monotonically decreasing objective values, so it will at least converge to a local minimum.

In addition to the original constraint $\sum_i p_i = 1,$ we also require that $p_i \ge \epsilon\, \forall i$ for some positive constant $\epsilon,$ to avoid have the logs blow up. Let $f(p)$ denote the original objective function. We will use a change of variables. The new variables $x_0,\dots,x_n$ are defined by $x_0=0$ and $x_j = \sum_{i=1}^j p_i$ for $j=1,\dots,n$ (implying $x_n = 1$). The reverse transformation is $p_i = x_i - x_{i-1},$ which we denote by $p = \pi(x).$ To ensure $p_i \ge \epsilon,$ we require $x_{i-1} + \epsilon \le x_i \le x_{i+1} - \epsilon.$

Now let $g(x)=f(\pi(x)).$ We want to minimize $g(x)$ subject to $x_0=1 < x_1 < \dots < x_n = 1$ with the $\epsilon$ difference requirements described above. The heuristic proceeds as follows.

  1. Generate an initial feasible solution and record it as an incumbent. This could be done randomly, or you could start with $p=(1/n,1/n,\dots,1/n)$ (and convert it to the equivalent $x$).

  2. Repeat the following loop until a full pass through the loop produces no change in the incumbent solution.

    a. Shuffle the indices $1,\dots,n-1,$ calling the resulting list $L.$

    b. For each $i\in L,$ search for the value $\xi$ in the interval $[x_{i-1} + \epsilon, x_{i+1} - \epsilon]$ that minimizes $g(x_0,\dots, x_{i-1}, \xi, x_{i+1},\dots, x_n).$ If the new objective value is less than that of the incumbent, make the new $x$ vector the new incumbent.

Step 2b requires an algorithm to optimize a function of one variable over a closed interval. Depending on what language you are using, there may well be library routines for this available. (I did some tests in R, using the optimize method from the standard stats package.) Worst case, you can code something like golden section search very easily.

The virtue of using a random starting solution is that, in the event that the global minimum is not the only local minimum, you can run the heuristic multiple times with different starts and hope to get lucky.

Update: I also tried gradient projection in R, using the BB library. It generated a better solution to my test problem, in considerably less time. The time difference is in part due to BB being coded in a compiled language by better coders than I, but the improved solution is attributable to the method.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Interesting algorithm. What is the intuition behind it? What would be the benefit over simple gradient descent? $\endgroup$
    – JEK
    Apr 23 at 17:01
  • $\begingroup$ Aside from the annoyance of deriving the gradient (or estimating it numerically), there are the constraints to deal with. The requirement that the probabilities sum to 1 can be dealt with either by adding a penalty term to the objective for violations or by replacing $p_n$ with $1 - p_1 - \dots - p_{n-1}.$ The bounds $\epsilon \le p_i $ are a little trickier (involving deflection of the gradient if it points outside the feasible region), although I think some gradient descent algorithms/libraries handle that. There's also possible numerical issues (partial derivative $\approx 1/\epsilon$). $\endgroup$
    – prubin ♦
    Apr 23 at 18:21
  • $\begingroup$ All that said, gradient descent is certainly worth considering. As for intuition, the change in variables is to finesse the constraint issue. Rather than tweaking $p_i$ and worrying how it affects all the other $p_j,$ we can nudge $x_i$ while keeping it in a well defined interval. The downside relative to gradient descent is that we only change one variable at a time, so potentially (likely?) more iterations are needed. $\endgroup$
    – prubin ♦
    Apr 23 at 18:23
  • $\begingroup$ As noted in my update, gradient projection outperformed my heuristic. $\endgroup$
    – prubin ♦
    Apr 23 at 22:24
  • 1
    $\begingroup$ "Aside from the annoyance of deriving the gradient (or estimating it numerically),. Gradients can be calculated by (reverse mode) automatic differentiation. Some modeling systems do this automatically under the hood, such as AMPL and YALMIP. The BARON global optimization solver computes all derivatives (by automatic differentiation) itself, whether BRON is called directly, or though an optimization ,modeling system. There are many good nonlinear optimization solvers which will horrid linear and nonlinear constraints - no need for horrid gradient "descent" (its name is false advertising). $\endgroup$
    – Mark L. Stone
    Apr 24 at 12:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.