Nonlinear optimization with constraints

Let $w$ be some vector. I want to solve: $$w=\arg\min \sum f_i^2(w) \\s.t. w\geq 0, \sum w_i = 1$$ That can be formulated as follows: $$w=\arg\min ||{\bf f}(w)||^2 \\ s.t. \quad H(w)\geq 0, G(w)= 0$$ where $H(w)=w,G(w)={\bf 1}^T w-1=\sum w_i-1$.

This seems to me a standard problem. The unconstrained problem can be solved using the Gauss-Newton (GN) method. However, I can't seem to find how to solve the constrained problem. Which method should I use? can GN be adapted to these constraints?

Thank you.

EDIT:

a. I will add a formula for $f_i(w)$: $$f_i(w) = w^T A_i w - z_i$$ where $A_i$ is a symmetric matrix and $z_i$ are known scalars.

b. Trying to find my answer, I encountered the Levenberg–Marquardt algorithm, which appears almost suitable here, but I do have a matrix $A_i$ and not a scalar that varies with $i$.

c. I thought of doing the following: since in each iteration of solving the Gauss-Newton method I solve a linear LS, I would solve instead a constrained LS with my constraints. Is that a right way to go?

Thanks again.

• Do you have a formula for $f_i$? – LinAlg Dec 14 '16 at 22:25
• Thanks, I added it + some clarifications. – yoki Dec 15 '16 at 7:39
• Are $A_i$ positive semidefinite? – LinAlg Dec 16 '16 at 9:26
• Are the scalars $z_i$ positive? – littleO Dec 21 '16 at 9:14
• Since your objective is a sum of squares of polynomials, I suggest you take a look at some polynomial optimization methods (see homepages.cwi.nl/~monique/files/moment-ima-update-new.pdf) . Yet, the objective can be written in a tensor form (See also researchgate.net/publication/…). – Alex Silva Dec 21 '16 at 9:42

Your optimization problem has the form \begin{align} \tag{$\spadesuit$}\operatorname{minimize}_w &\quad F(w) \\ \text{subject to} & \quad w \in C \end{align} where $$F(w) = \sum_i f_i(w)^2$$ and $C$ is the probability simplex (which is a closed, convex set).

You can solve problems of the form ($\spadesuit$) using the projected gradient method: $$w^{k+1} = \Pi_C(w^k - t\nabla F(w^k)).$$ Here $\Pi_C(y)$ is the projection of $y$ onto $C$. Projecting onto the probability simplex is a very well studied operation and can be performed very efficiently with existing algorithms. (It takes only a few lines of Matlab. See algorithm 1 here, for example: https://arxiv.org/abs/1309.1541 .)

The projected gradient method has a step size restriction, so be sure to select the step size $t$ small enough so that the step size restriction is satisfied. There is also a line search procedure you can use to select step sizes(which I think is probably a good idea).

You could also try using an accelerated version of the projected gradient method, such as FISTA, although I'm not sure how they will perform for nonconvex problems.

Edit: The projected gradient method is a special case of the proximal gradient method, and a line search version of the proximal gradient method is described in equations 1.26 and 1.27 (p. 18) of the paper "Gradient-Based Algorithms with Applications to Signal Recovery Problems" by Beck and Teboulle. This line search version of the proximal gradient method is also described, with slightly different notation, on slide 6-20 of Vandenberghe's 236c notes: http://www.seas.ucla.edu/~vandenbe/236C/lectures/proxgrad.pdf.

For convenience, I'll write the proximal gradient method with line search in detail, using notation that I think is a little more clear than the notation in the above references. This algorithm minimizes $f(x) + g(x)$, where $f$ is assumed to be differentiable (with a Lipschitz continuous gradient) and $g$ is a convex, lower semicontinuous function. The proximal operator of $g$ with parameter $t > 0$ is denoted by $\operatorname{prox}_{tg}$. If you don't know what the proximal operator is, that's fine, just replace $\operatorname{prox}_{tg}$ with $\Pi_C$, and you have the projected gradient method. Here is the algorithm: \begin{align} &\text{Initialize } x_0 \text{ and select } t_0 > 0, 0 < \alpha < 1, \beta > 1\\ &\textbf{for } k = 1,2,\ldots \textbf{do:}\\ & \quad t := \beta \, t_{k-1} \\ &\quad \textbf{repeat:} \\ &\qquad \quad x := \operatorname{prox}_{tg}(x_{k-1} - t \nabla f(x_{k-1})) \\ &\qquad \quad \textbf{break if } f(x) \leq f(x_{k-1}) + \langle \nabla f(x_{k-1}),x - x_{k-1} \rangle + \frac{1}{2t} \| x - x_{k-1} \|_2^2 \\ &\qquad \quad t := \alpha t \\ & \quad t_k := t \\ & \quad x_k := x\\ &\textbf{end for} \end{align}

I usually select $\alpha = .5$ and $\beta = 1.25$. I think $t_0 = 1$ is usually not a bad choice. You could initialize $x_0$ to be your best available guess of the optimal value of $x$, or perhaps just initialize $x_0$ randomly or to all zeros.

Here is some python code for the proximal gradient method with line search:

def proxGrad(gradf,proxg,evalf,evalg,x0,params):
# minimize f(x) + g(x). f is smooth, g has easy prox-operator

maxIter = params['maxIter']
t = params['stepSize'] # initial step size
showTrigger = params['showTrigger']
reduceFactor = .5
increaseFactor = 1.25
costs = np.zeros((maxIter,1))

xk = x0
for k in np.arange(maxIter):

gxk = evalg(xk)
costs[k] = fxk + gxk
if k % showTrigger == 0:
print "Iteration: " + str(k) + "    cost: " + str(costs[k])

t = t*increaseFactor
acceptFlag = False
while acceptFlag == False:
xkp1 = proxg(xk - t*gradf_xk, t)
fxkp1 = evalf(xkp1)
diff = xkp1 - xk
if fxkp1 <= fxk + np.vdot(gradf_xk,diff) + (.5/t)*np.sum(diff**2):
acceptFlag = True
else:
t = t*reduceFactor
print "Reducing t. t is now: " + str(t)

xk = xkp1

return (xkp1,costs)


Here's a little more info about the relationship between the proximal gradient method and the projected gradient method. The projected gradient method for minimizing $f(x)$ subject to the constraint that $x \in C$, is just the special case of the proximal gradient method where $g$ is the convex indicator function of $C$: $$g(x) = \begin{cases} 0 & \quad \text{if } x \in C \\ \infty & \quad \text{otherwise.} \end{cases}$$ With this choice of $g$, the proximal operator of $g$ is given by the formula $$\operatorname{prox}_{tg}(x) = \Pi_C(x) = \text{projection of x onto C}.$$ (I'm assuming that $C$ is a closed convex set.)

• How do I calculate the cost function for the line search? when there are no constraints it's not a problem, but here should I calculate the cost function $F(w^k_t)$ (i.e. the cost for some $w$ with a specific step size $t$) when $w^k$ is after the projection or before? – yoki Dec 21 '16 at 9:35
• @yoki I updated my answer to give details on the line search procedure. – littleO Dec 21 '16 at 10:42
• Your function $F$ is not convex. I do not see how your algorithm can guarantee optimality. – LinAlg Dec 21 '16 at 11:04
• You're right that this approach does not guarantee global optimality. The Beck and Teboulle paper I referenced does give a convergence result in the case where $f$ is non-convex (Theorem 1.3, p. 21), but convergence to a global minimizer is not guaranteed. – littleO Dec 21 '16 at 11:20
• @littleO Thank you very much for the explanation! – yoki Dec 21 '16 at 12:33

As noted in the comments, this answer does not give an exact solution.

Using an extra variable $t_i$, which is used to model $f_i(w)$, the problem can be represented as a quadratically constrained quadratic optimization problem (QCQP): \begin{align*} \min \quad & t_i^2 \\ \mathrm{s.t.} \quad & t_i \geq w^T A_i w - z_i \\ & \sum w_i = 1 \\ & w \geq 0 \end{align*} You may wonder why the first constraint uses "$\geq$" instead of "$=$". The reason is that it ensures convexity (given two feasible solutions, the line segment between those solutions is also feasible). However, as noted in the comments, the formulation is not equivalent.

The advantage of writing this problem as a QCQP is that there is a plethora of algorithms to solve it. I will focus on interior point methods, which find the (globally) optimal solution in 40-60 iterations independent of the size of your problem. Matlab has quadprog. Commercial software like CPLEX, Mosek, Gurobi has bindings for many languages, including Matlab and Python. QCQP is a special case of convex optimization, so in Python you could use scipy.optimize.

• Is there a problem because the objective function $\sum_i t_i^2$ can never push $t_i$ to a value that is less than $0$? So it is possible that at optimality the constraint $t_i \geq w^T A_i w - z_i$ does not hold with equality? – littleO Dec 21 '16 at 11:14
• You are right, if $w^T A_i w - z_i<0$, my formulation can take $t_i=0$. – LinAlg Dec 21 '16 at 11:18