# Optimize least squares penalized by curvature of log pdf

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?
• Are you assuming that the p_i values are sorted in increasing order? Apr 22, 2023 at 4:09
• 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? Apr 22, 2023 at 13:56
• Are you interested in heuristic solutions (not guaranteed to find an optimum)?
– prubin
Apr 22, 2023 at 19:56
• @KevinDalmeijer no I’m not assuming they are increasing in order
– JEK
Apr 23, 2023 at 1:47
• @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).
– JEK
Apr 23, 2023 at 1:50

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). • Thanks. I suggested this in my first comment, but was too lazy to try it out. Apr 24, 2023 at 12:16 • "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. Apr 24, 2023 at 12:43 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. • Interesting algorithm. What is the intuition behind it? What would be the benefit over simple gradient descent? – JEK Apr 23, 2023 at 17:01 • 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$). – prubin Apr 23, 2023 at 18:21 • 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.
– prubin
Apr 23, 2023 at 18:23
• As noted in my update, gradient projection outperformed my heuristic.
– prubin
Apr 23, 2023 at 22:24
• "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). Apr 24, 2023 at 12:22