An Interior-Point Lagrangian Decomposition Method for Separable Convex Optimization

Introduction

Can self-concordance and interior-point methods be incorporated into a Lagrangian dual decomposition framework? This paper presents a decomposition algorithm that incorporates the interior-point method into augmented Lagrangian decomposition technique for solving large-scale separable convex problems. Separable convex problems, i.e. optimization problems with a separable convex objective function but with coupling constraints, arise in many fields: networks (communication networks, multicommodity network flows) , process system engineering (e.g. distributed model predictive control) , stochastic programming , etc. There has been considerable interest in parallel and distributed computation methods for solving this type of structured optimization problems and many methods have been proposed: dual subgradient methods , alternating direction methods , proximal method of multipliers , proximal center method , interior-point based methods , etc.

The methods mentioned above belong to the class of augmented Lagrangian or multiplier methods , i.e. they can be viewed as techniques for maximizing an augmented dual function. For example in the alternating direction method a quadratic penalty term is added to the standard Lagrangian to obtain a smooth dual function and then using a steepest ascent update for the multipliers. However, the quadratic term destroys the separability of the given problem. Moreover, the performance of these methods is very sensitive to the variations of their parameters and some rules for choosing these parameters were given e.g. in . In the proximal center method we use smoothing techniques in order to obtain a well-behaved Lagrangian, i.e. we add a separable strongly convex term to the ordinary Lagrangian. This technique leads to a smooth dual function, i.e. with Lipschitz continuous gradient, which preserves separability of the problem, the corresponding parameter is selected optimally and moreover the multipliers are updated using an optimal gradient based scheme. In interior-point methods are proposed for solving special classes of separable convex problems with a particular structure of the coupling/local constraints. In those papers the Newton direction is used to update the primal variables and/or multipliers obtaining polynomial-time complexity for the proposed algorithms. In the present paper we use a similar smoothing technique as in in order to obtain a well-behaved augmented dual function. Although we relax the coupling constraints using the Lagrangian dual framework as in , the main difference here is that the smoothing term is a self-concordant barrier, while in the main property of the smoothing term was strong convexity. Therefore, using the properties of self-concordant functions we show that the augmented dual function becomes under mild assumptions also self-concordant. Hence the Newton direction can be used instead of gradient based directions as it is done in most of the augmented Lagrangian methods. Furthermore, we develop a specialized interior-point method to maximize the augmented dual function which takes into account the special structure of our problem. We present a parallel algorithm for computing the Newton direction of the dual function and we also prove global convergence of the proposed method.

The main contributions of the paper are the following: (i) We consider a more general model for separable convex problems that includes local equality and inequality constraints, and linear coupling constraints which generalizes the models in . (ii) We derive sufficient conditions for self-concordance of augmented Lagrangian and we prove self-concordance for the corresponding family of augmented dual functions. (iii) We provide an interior-point based algorithm for solving the dual problem with proofs of global convergence and polynomial-time complexity. (iv) We propose a practical implementation of the algorithm based on solving approximately the subproblems and on parallel computations of the Newton directions.

Note that item (ii) generalizes the results of . However, the consideration of general convex problems with local equality constraints requires new proofs with more involved arguments in order to prove self-concordance for the family of dual functions.

This paper is organized as follows. In Section 2 we formulate the separable convex problem followed by a brief description of some of the existing decomposition methods for this problem. The main results are given in Sections 3 and 4. In Section 3 we show that the augmented Lagrangian obtained by adding self-concordant barrier terms to the ordinary Lagrangian forms a self-concordant family of dual functions. Then an interior-point Lagrangian decomposition algorithm with polynomial complexity is proposed in Section 4. The new algorithm makes use of the special structure of our problem so that it is highly parallelizable and it can be effectively implemented on parallel processors. We conclude the paper with some possible applications.

Throughout the paper we use the following notations. For a function $\psi$ with two arguments, scalar parameter $t$ and decision variable $x$, i.e. $\psi(t,x)$, we use “ $\prime$” to denote the partial derivative of $\psi(t,x)$ with respect to $t$ and “$\nabla$” with respect to $x$: e.g. $\nabla\psi^{\prime}(t,x)=\frac{\partial^{2}}{\partial t\partial x}\psi(t,x)$. For a function $\phi$, three times differentiable, i.e. $\phi\in{\cal C}^{3}(\text{dom}\!\ \phi)$, $\nabla^{3}\phi(x)[h_{1},h_{2},h_{3}]$ denotes the third differential of $\phi$ at $x$ along directions $h_{1},h_{2}$ and $h_{3}$. We use the notation $A\preceq B$ if $B-A$ is positive semidefinite. We use $D_{A}$ to denote the block diagonal matrix having on the main diagonal the matrices $A_{1},\cdots,A_{N}$. We use $\text{int}(X)$ to denote the interior of a set $X$.

Problem Formulation

We consider the following general separable convex optimization problem:

In order to obtain a smooth dual function we need to use smoothing techniques applied to the ordinary Lagrangian $L_{0}$. One approach is the augmented Lagrangian obtained e.g. by adding a quadratic penalty term to the Lagrangian $L_{0}$: $t\|Bx-a\|^{2}$. In the alternating direction method the minimization of the augmented Lagrangian is performed by alternating minimization in a Gauss-Seidel fashion followed by a steepest ascent update for the multipliers.

In we proposed the proximal center method in which we added to the standard Lagrangian a smoothing term $t\sum_{i=1}^{N}g_{X_{i}}(x_{i})$, where each function $g_{X_{i}}$ is strongly convex and depends on the set $X_{i}$ so that the augmented Lagrangian takes the following form:

Therefore, the augmented Lagrangian $L_{t}^{\text{prox}}$ is strongly convex, preserves separability of the problem like $L_{0}$ and the associated augmented dual function

is differentiable and has also a Lipschitz continuous gradient. In an accelerated gradient based method is used to maximize the augmented dual function $d_{t}^{\text{prox}}$, while the corresponding minimization problems are solved in parallel. Moreover, the smoothing parameter $t$ is selected optimally.

Note that the methods discussed above use only the gradient directions of the augmented dual function in order to update the multipliers. Therefore, in the absence of more conservative assumptions like strong convexity, the global convergence rate of these methods is slow, in general sub-linear. In this paper we propose to smoothen the Lagrangian by adding instead of strongly convex terms $g_{X_{i}}$, self-concordant barrier terms $\phi_{X_{i}}$ associated with the sets $X_{i}$, in order to obtain the self-concordant Lagrangian:

In the next section we show, using the theory of self-concordant barrier functions , that for a relatively large class of convex functions $f_{i}$ (see also Section 5), we can obtain a self-concordant augmented dual function:

This opens the possibility of deriving an interior-point method using Newton directions for updating the multipliers to speed up the convergence rate of the proposed algorithm.

Sufficient Conditions for Self-Concordance of the Augmented Dual Function

In this section we derive sufficient conditions under which the family of augmented dual functions is self-concordant. A key property that allows to prove polynomial convergence for barrier type methods is the property of self-concordance (see Definition 2.1.1 in ):

Note that (5) is equivalent to (see , pp. 14):

Moreover, if Hessian $\nabla^{2}\phi(x)$ is positive definite, then the inequality (6) is equivalent to

Next lemma provides some basic properties of self-concordant functions:

Note that a self-concordant function which is also a barrier for its domain is called strongly self-concordant. The next lemma gives some helpful composition rules for self-concordant functions.

then $\bar{\psi}_{t}(x)=\psi(x)-t\sum_{i=1}^{n}\log(u_{i}-x_{i})(x_{i}-l_{i})$ is $2(1+\beta)/\sqrt{t}$-self concordant.

(i) and (ii) can be found in , pp. 13. (iii) Denote $\phi_{\text{box}}(x)=-\sum_{i=1}^{n}\log(u_{i}-x_{i})(x_{i}-l_{i})$. Note that

and using Cauchy-Schwarz inequality it follows that $\phi_{\text{box}}$ is $2$-self-concordant function on $\text{int}(X_{\text{box}})$. Let us denote

Using hypothesis (9) and 2-self-concordance of $\phi_{\text{box}}$ we have the following inequalities:

With some computations we can observe that

Note that condition (9) is similar to $\psi$ is $\beta$-compatible with $\phi_{\text{box}}$ on $X_{\text{box}}$, defined in . The following assumptions will be valid throughout this section:

We analyze the following prototype minimization problem:

In the following four lemmas we derive the main properties of the family of augmented dual functions $\{d(t,\cdot)\}_{t>0}$. We start with a linear algebra result:

and thus $x=0$ which is a contradiction. ∎

If Assumption 3.1 holds, then for any $t>0$ the function $d(t,\cdot)$ is $M_{t}$-self-concordant, where $M_{t}$ is either $2/\sqrt{t}$ or $\max\{M_{f},2/\sqrt{t}\}$ or $2(1+\beta)/\sqrt{t}$.

Since $f$ is assumed to be either linear or convex quadratic or $M_{f}$-self-concordant or $X$ is a box and $f$ satisfies condition (9) it follows from Lemma 3.1 that $f+t\phi$ is also $M_{t}$-self concordant (where $M_{t}$ is either $2/\sqrt{t}$ or $\max\{M_{f},2/\sqrt{t}\}$ or $2(1+\beta)/\sqrt{t}$, respectively) and with positive definite Hessian (according to our assumptions and Proposition 3.1). Moreover, $f+t\phi$ is strongly self-concordant since $\phi$ is a barrier function for $X$. Since $A$ has full row rank and $p<n$, then there exists some vectors $u_{i},i=1\cdots n-p$, that form a basis of the null space of this matrix. Let $U$ be the matrix having as columns the vectors $u_{i}$ and $x_{0}$ a particular solution of $Ax=a$. Then, for a fixed $t$, the feasible set of (10) can be described as

which is an open convex set. Using that self-concordance is affine invariant it follows that the functions $\bar{f}(y)=f(x_{0}+Uy)$, $\bar{\phi}(y)=\phi(x_{0}+Uy)$ have the same properties as the functions $f$, $\phi$, respectively, that $\bar{f}+t\bar{\phi}$ is also $M_{t}$-self concordant and that

From our assumptions and Proposition 3.1 it follows that the Hessian of $\phi$ and $\bar{\phi}$ are positive definite. Since $f$ is convex it follows that the Hessian of $\bar{f}+t\bar{\phi}$ is also positive definite and thus invertible. Let

Since \left[\begin{array}[]{c}A\\ B\end{array}\right] has full row rank, then from Lemma 3.2 $BU$ has full row rank. Moreover, since $\nabla^{2}\bar{F}(t,\cdot)$ is positive definite and

it follows that $\nabla^{2}d(t,\cdot)$ is positive definite on its domain

Moreover, since self-concordance is affine invariant it follows that $d(t,\cdot)$ is also $M_{t}$-self-concordant on the domain $X_{d(t,\cdot)}$. ∎

First we determine the formula for the Hessian. It follows immediately from (11) that

For simplicity, we drop the dependence of all the functions on $x(t,\lambda)$ and $(t,\lambda)$. Differentiating (11) with respect to $\lambda$ we arrive at the following system in $\nabla x$ and $\nabla\nu$:

Since $H$ is positive definite and according to our assumption $A$ is full row rank, it follows that the system matrix is invertible. Using the well-known formula for inversion of partitioned matrices we find that:

Differentiating the first part of (11) with respect to $t$ and using the same procedure as before we arrive at a similar system as above in the unknowns $x^{\prime}$ and $\nu^{\prime}$. We find that

We also introduce the following notation: $F:=H^{-1}A^{T}(AH^{-1}A^{T})^{-1}AH^{-1}$ and $G:=H^{-1}-F$, which are positive semidefinite. Using a similar reasoning as in and Cauchy-Schwarz inequality we obtain:

We recall that $H(t,\lambda)=\nabla^{2}f(x(t,\lambda))+t\nabla^{2}\phi(x(t,\lambda))$. Therefore

We again drop the dependence on $(t,\lambda)$ and after some straightforward algebra computations we arrive at the following expression:

Taking into account the expression of $H^{\prime}$ derived above we obtain:

Using the self-concordance property (7) for $f+t\phi$ we obtain that:

Moreover, since $f$ is convex, $\nabla^{2}f$ is positive semidefinite and thus:

Combining the last two inequalities we obtain:

With some algebra we can check that the following identity holds: $FH(H^{-1}~-~F)~=~0$. Based on this identity we can compute $u^{T}Hu$ and $(x^{\prime})^{T}Hx^{\prime}$. Indeed,

The inequality from lemma follows then by replacing the last two relations in (13). ∎

The main result of this section is summarized in the next theorem.

Under the Assumption 3.1, $\{d(t,\lambda)\}_{t>0}$ is a strongly self-concordant family in the sense of Definition Note that according to Definition 3.1.1 in $\gamma_{t}=1$ and $\mu_{t}=1$ for our case. 3.1.1 in with parameters $\alpha_{t}=M_{t},\xi_{t}=(M_{t}/2)\sqrt{N_{\phi}/t}$ and $\eta_{t}=(M_{t}/2)\sqrt{N_{\phi}/t}+(1/2t)$, where $M_{t}$ is defined in Lemma 3.3.

Basically, from Definition 3.1.1 in we must check three properties: self-concordance of $d(t,\lambda)$ (Lemma 3.3) and that the first and second order derivative of $d(t,\cdot)$ vary with $t$ at a rate proportional to the derivative itself (Lemmas 3.4 and 3.5). In conclusion, the Lemmas 3.3–3.5 prove our theorem. ∎

It is known that self-concordant families of functions can be minimized by path-following methods in polynomial time. Therefore, this type of family of augmented dual functions $\{d(t,\cdot)\}_{t>0}$ plays an important role in the algorithm of the next section.

Parallel Implementation of an Interior-Point Based Decomposition Method

In this section we develop an interior-point Lagrangian decomposition method for the separable convex problem given by (1)–(2). Our previous Theorem 3.1 is the major contribution of our paper since it allows us to effectively utilize the Newton method for tracing the trajectory of optimizers of the self-concordant family of augmented dual functions (4).

The following assumptions for optimization problem (1)–(2) will be valid in this section:

Note that boundedness of the set $X_{i}$ can be relaxed to $X_{i}$ does not contain straight lines and the set of optimal solutions to problem (1)–(2) is bounded. Note also that the rank assumption (iii) is not restrictive since we can eliminate the redundant equalities (see also Lemma 3.2 for other less restrictive conditions). The constraint qualification condition from Assumption 4.1 (iii) guarantees that strong duality holds for problem (1)–(2) and thus there exists a primal-dual optimal solution $(x^{*},\lambda^{*})$.

Note that the function $d(t,\cdot)$ can be computed in parallel by decomposing the original large optimization problem (1)–(2) into $N$ independent small convex subproblems.

(i) The family $\{d_{i}(t,\cdot)\}_{t>0}$ is strongly self-concordant with the parameters $\alpha_{i}(t)=M_{i}(t),\xi_{i}(t)=(M_{i}(t)/2)\sqrt{N_{i}/t}$ and $\eta_{i}(t)=(M_{i}(t)/2)\sqrt{N_{i}/t}+(1/2t)$, where $M_{i}(t)$ is either $2/\sqrt{t}$ or $\max\{M_{f_{i}},2/\sqrt{t}\}$ or $2(1+\beta)/\sqrt{t}$ for all $i=1\cdots N$. (ii) The family $\{d(t,\cdot)\}_{t>0}$ is strongly self-concordant with parameters $\alpha(t)=\alpha/\sqrt{t}$, $\xi(t)=\xi/t$ and $\eta(t)=\eta/t$, for some fixed positive constants $\alpha,\xi$ and $\eta$ depending on $(N_{i},M_{f_{i}},\beta)$.

(i) is a straightforward consequence of Assumption 4.1 and Theorem 3.1.

From Assumption 4.1 and the discussion from previous section, it follows that the optimizer of each maximization is unique and denoted by

The central path $\{(x(t,\lambda(t)),\lambda(t)):t>0\}$ converges to the optimal solution $(x^{*},\lambda^{*})$ as $t\to 0$ and $\{x(t,\lambda(t)):t>0\}$ is feasible for the problem (1)–(2).

Let $x(t):=\arg\min_{x}\{f(x)+t\phi_{X}(x):Bx=b,x_{i}\in\text{int}(X_{i}),A_{i}x_{i}=a_{i}\;\forall i\}$, then it is known that $x(t)\to x^{*}$ as $t\to 0$. It is easy to see that the Hessian of $f+t\phi_{X}$ is positive definite and thus $f+t\phi_{X}$ is strictly convex and $x(t)$ is unique. From Assumption 4.1 it also follows that strong duality holds for this barrier function problem and therefore

In conclusion, $x(t)=x(t,\lambda(t))$ and thus $x(t,\lambda(t))\to x^{*}$ as $t\to 0$. As a consequence it follows that $x(t,\lambda(t))$ is feasible for the original problem, i.e. $Bx(t,\lambda(t))~=~b$, $A_{i}x_{i}(t,\lambda(t))=a_{i}$ and $x_{i}(t,\lambda(t))\in\text{int}(X_{i})$. It is also clear that $\lambda(t)\to\lambda^{*}$ as $t\to 0$.∎

The next theorem describes the behavior of the central path:

For $x(t)=x(t,\lambda(t))$ the following bound holds for the central path: given any $0<\tau<t$ then,

It follows immediately that $\langle\nabla f(x(s)),x^{\prime}(s)\rangle=-s\langle\nabla\phi_{X}(x(s)),x^{\prime}(s)\rangle$. Since $0<\tau<t$, then there exists $\tau\leq s\leq t$ such that

Using a similar reasoning as in Lemma 3.4 we have:

where we denote with $H(s)=\nabla^{2}f(x(s))+s\nabla^{2}\phi_{X}(x(s))$. Using (8), the expression for $x^{\prime}(s)$ and since $0\prec\nabla^{2}\phi_{X}(x(s))\preceq 1/sH(s)$ and $H^{-1}(s)\preceq 1/s\big(\nabla^{2}\phi_{X}(x(s))\big)^{-1}$ we obtain:

It follows immediately that $f(x(t))-f(x(\tau))\leq N_{\phi}(t-\tau)$. ∎

A simple consequence of Theorem 4.1 is that the following bounds on the approximation of the optimal value function $f^{*}$ hold:

It is easy to see that the gradient of the self-concordant function $d(t,\cdot)$ is given by

For every $(t,\lambda)$ let us define the positive definite matrix

The Hessian of function $d_{i}(t,\cdot)$ is positive definite and from (12) it has the form

In conclusion, the Hessian of dual function $d(t,\cdot)$ is also positive definite and given by:

Denote the Newton direction associated to self-concordant function $d(t,\cdot)$ at $\lambda$ with

For every $t>0$, we define the Newton decrement of the function $d(t,\cdot)$ at $\lambda$ as:

replace $r$ by $r+1$ and go to Step 1 Step 3. output $(t^{0},\lambda^{0})=(t_{0},\lambda_{r_{f}})$.

Note that Algorithm 4.1 approximates the optimal Lagrange multiplier $\lambda(t_{0})$ of the dual function $d(t_{0},\cdot)$, i.e. the sequence $(t_{0},\lambda_{r})$ moves into the neighborhood $V(t,\epsilon_{V})=\{(t,\lambda):\delta(t,\lambda)\leq\epsilon_{V}\}$ of the trajectory $\{(t,\lambda(t)):t>0\}$.

(Path-Following Algorithm) Step 0. input: $(t^{0},\lambda^{0})$ satisfying $\delta(t^{0},\lambda^{0})\leq\epsilon_{V}$ , $k=0$, $0<\tau<1$ and $\epsilon>0$ Step 1. if $t^{k}N_{\phi}\leq\epsilon$, then $k_{f}=k$ and go to Step 5 Step 2. (outer iteration) let $t^{k+1}=\tau t^{k}$ and go to inner iteration (Step 3) Step 3. (inner iteration) initialize $\lambda=\lambda^{k}$, $t=t^{k+1}$ and $\delta=\delta(t^{k+1},\lambda^{k})$ while $\delta>\epsilon_{V}$ do

Step 3.1 compute $x_{i}=x_{i}(t,\lambda)\;\forall i$, determine a step size $\sigma$ and compute

$\lambda^{+}=\lambda+\sigma\Delta\lambda(t,\lambda)$

Step 3.2 compute $\delta^{+}=\delta(t,\lambda^{+})$ and update $\lambda=\lambda^{+}$ and $\delta=\delta^{+}$ Step 4. $\lambda^{k+1}=\lambda$ and $x_{i}^{k+1}=x_{i}$; replace $k$ by $k+1$ and go to Step 1 Step 5. output: $(x_{1}^{k_{f}},\cdots,x_{N}^{k_{f}},\lambda^{k_{f}})$.

In Algorithm 4.2 we trace numerically the trajectory $\{(t,\lambda(t)):t>0\}$ from a given initial point $(t^{0},\lambda^{0})$ close to this trajectory. The sequence $\{(x_{1}^{k},\cdots,x_{N}^{k},\lambda^{k})\}_{k>0}$ lies in a neighborhood of the central path and each limit point of this sequence is primal-dual optimal. Indeed, since $t^{k+1}=\tau t^{k}$ with $\tau<1$, it follows that $\lim_{k\to\infty}t^{k}=0$ and using Theorem 4.1 the convergence of the sequence $x^{k}=[(x_{1}^{k})^{T}\cdots(x_{N}^{k})^{T}]^{T}$ to $x^{*}$ is obvious.

The step size $\sigma$ in the previous algorithms is defined by some line search rule. There are many strategies for choosing $\tau$. Usually, $\tau$ can be chosen independent of the problem (long step methods), e.g. $\tau=0.5$, or depends on the problem (short step methods). The choice for $\tau$ is crucial for the performance of the algorithm. An example is that in practice long step interior-point algorithms are more efficient than short step interior-point algorithms. However, short step type algorithms have better worst-case complexity iteration bounds than long step algorithms. In the sequel we derive a theoretical strategy to update the barrier parameter $\tau$ which follows from the theory described in and consequently we obtain complexity bounds for short step updates. Complexity iteration bounds for long step updates can also be derived using the same theory (see Section 3.2.6 in ). The next lemma estimates the reduction of the dual function at each iteration.

(ii) If $\delta\leq\delta_{*}$, then defining the Newton iterate $\lambda^{+}=\lambda+\Delta\lambda$ we have

(iii) If $\delta\leq\delta_{*}/2$, then defining $t^{+}=\frac{2c}{2c+1}t$, where $c=1/4+2\xi/\delta_{*}+\eta$, we have

(i) and (ii) follow from Theorem 2.2.3 in and Lemma 4.1 from above.

(iii) is based on the result of Theorem 3.1.1 in . In order to apply this theorem, we first write the metric defined by (3.1.4) in for our problem: given $0<t^{+}<t$ and using Lemma 4.1 we obtain

Since $\delta\leq\delta_{*}/2<\delta_{*}$ and since for $t^{+}=\frac{2c}{2c+1}t$, where $c$ is defined above, one can verify that $\rho_{\delta_{*}/2}(t,t^{+})=c\log(1+1/2c)\leq 1/2\leq 1-\delta/\delta_{*}$, i.e. $t^{+}$ satisfies the condition (3.1.5) of Theorem 3.1.1 in , it follows that $\delta(t^{+},\lambda)\leq\delta_{*}$. ∎

Define the following step size: $\sigma(\delta)=1/(1+\delta)$ if $\delta>\delta_{*}$ and $\sigma(\delta)=1$ if $\delta\leq\delta_{*}$. With Algorithm 4.1 for a given $t^{0}$ and $\epsilon_{V}=\delta_{*}/2$, we can find $(t^{0},\lambda^{0})$ satisfying $\delta(t^{0},\lambda^{0})\leq\delta_{*}/2$ using the step size $\sigma(\delta)$ (see previous lemma). Based on the analysis given in Lemma 4.3 it follows that taking in Algorithm 4.2 $\epsilon_{V}=\delta_{*}/2$ and $\tau=2c/(2c+1)$, then the inner iteration stage (step 3) reduces to only one iteration:

Step 3. compute $\lambda^{k+1}=\lambda^{k}+\Delta\lambda(t^{k+1},\lambda^{k})$.

However, the number of outer iterations is larger than in the case of long step algorithms.

2 Practical Implementation

for some $\epsilon_{x}>0$. Note however that even when such approximations are considered, the vector $\Delta\lambda$ still defines a search direction in the $\lambda$-space. Moreover, the cost of computing an extremely accurate maximizer of (14) as compared to the cost of computing a good maximizer of (14) is only marginally more, i.e. a few Newton steps at most (due to quadratic convergence of the Newton method close to the solution). Therefore, it is not unreasonable to assume even exact computations in the proposed algorithms.

In the rest of this section we discuss the complexity of our method and parallel implementations for solving efficiently the Newton direction $\Delta\lambda$. At each iteration of the algorithms we need to solve basically a linear system of the following form:

where $G_{i}=B_{i}[H_{i}^{-1}-H_{i}^{-1}A_{i}^{T}\big(A_{i}H_{i}^{-1}A_{i}^{T}\big)^{-1}A_{i}H_{i}^{-1}]B_{i}^{T}$, the positive definite matrix $H_{i}$ denotes the Hessian of $f_{i}+t\phi_{X_{i}}$ and some appropriate vector $g$. In order to obtain the matrices $H_{i}$ we can solve in parallel $N$ small convex optimization problems of the form (14) by Newton method, each one of dimension $n_{i}$ and with self-concordant objective function. The cost to solve each subproblem (14) by Newton method is ${\cal O}(n_{i}^{3}(n_{\lambda}+\log\log 1/t\epsilon_{x}))$, where $n_{\lambda}$ denotes the number of Newton iterations before the iterates $x_{i}$ reaches the quadratic convergence region (it depends on the update $\lambda$) and $t\epsilon_{x}$ is the required accuracy for the approximation of (14). Note that using the Newton method for solving (14) we automatically obtain also the expression for $H_{i}^{-1}$ and $A_{i}H_{i}^{-1}A_{i}^{T}$. Assuming that a Cholesky factorization for $A_{i}H_{i}^{-1}A_{i}^{T}$ is used to solve the Newton system corresponding to the optimization subproblem (14), then this factorization can also be used to compute in parallel the matrix of the linear system (15). Finally, we can use a Cholesky factorization of this matrix and then forward and backward substitution to obtain the Newton direction $\Delta\lambda$. In conclusion, we can compute the Newton direction $\Delta\lambda$ in ${\cal O}(\sum_{i=1}^{N}n_{i}^{3})$ arithmetic operations.

Note however that in many applications the matrices $H_{i}$, $A_{i}$ and $B_{i}$ are very sparse and have special structures. For example in network optimization (see Section 5.2 below for more details) the $H_{i}$’s are diagonal matrices, $B_{i}$’s are the identity matrices and the matrices $A_{i}$’s are the same for all $i$ (see (17)), i.e. $A_{i}=A$. In this case the Cholesky factorization of $AH_{i}^{-1}A^{T}$ can be done very efficiently since the sparsity pattern of those matrices is the same in all iterations and coincides with the sparsity pattern of $AA^{T}$, so the analyse phase has to be done only once, before optimization.

For large problem instances we can also solve the linear system (15) approximately using a preconditioned conjugate gradient algorithm. There are different techniques to construct a good preconditioner and they are spread across optimization literature. Detailed simulations for the method proposed in this paper and comparison of different techniques to solve the Newton system (15) will be given elsewhere.

Let us also note that the number of Newton iterations performed in Algorithm 4.1 can be determined via Lemma 4.3 (i). Moreover, if in Algorithm 4.2 we choose $\epsilon_{V}=\delta_{*}/2$ and $\tau=2c/(2c+1)$ we need only one Newton iteration at the inner stage. It follows that for this particular choice for $\epsilon_{V}$ and $\tau$ the total number of Newton iterations of the algorithm is given by the number of outer iterations, i.e. the algorithm terminates in polynomial-time, within ${\cal O}\big(\frac{1}{\log(\tau^{-1})}\log(N_{\phi}t^{0}/\epsilon)\big)$ iterations. This choice is made only for a worst-case complexity analysis. In a practical implementation one may choose larger values using heuristic considerations.

Applications with Separable Structure

In this section we briefly discuss some of the applications to which our method can be applied: distributed model predictive control and network optimization. Note that for these applications our Assumption 4.1 holds.

A first application that we will discuss here is the control of large-scale systems with interacting subsystem dynamics. A distributed model predictive control (MPC) framework is appealing in this context since this framework allows us to design local subsystem-base controllers that take care of the interactions between different subsystems and physical constraints. We assume that the overall system model can be decomposed into $N$ appropriate subsystem models:

A similar formulation of distributed MPC for coupled linear subsystems with decoupled costs was given in , but without state constraints. In , the authors proposed to solve the optimization problem (16) in a decentralized fashion, using the Jacobi algorithm . But, there is no theoretical guarantee of the Jacobi algorithm about how good the approximation to the optimum is after a number of iterations and moreover one needs strictly convex functions $f_{i}$ to prove asymptotic convergence to the optimum.

2 Network Optimization

Network optimization furnishes another area in which our algorithm leads to a new method of solution. The optimization problem for routing in telecommunication data networks has the following form :

Each function $f_{j}\in{\cal C}^{3}\big([0,\ d_{j})\big)$ is convex and $f_{j}$ is $3$-compatible with the self-concordant barrier $\phi_{j}(y_{j})=-\log(y_{j}(d_{j}-y_{j}))$ on the interval $(0,\ d_{j})$.

Therefore, we can solve this network optimization problem with our method. Note that the standard dual function $d_{0}$ is not differentiable since it is the sum of a differentiable function (corresponding to the variable $y$) and a polyhedral function (corresponding to the variable $x$). In a bundle-type algorithm is developed for maximizing the non-smooth function $d_{0}$, in the dual subgradient method is applied for maximizing $d_{0}$, while in alternating direction methods were proposed.

3 Preliminary Numerical Results

In the table we display the CPU time (seconds) and the number of calls of the dual function (i.e. the total number of outer and inner iterations) for our dual interior-point algorithm (DIP) and an algorithm based on alternating direction method (ADI) for different values of $m_{1},n_{1},N$ and fixed accuracy $\epsilon=10^{-4}$. For two problems the ADI algorithm did not produce the result after running one day. All codes are implemented in Matlab version 7.1 on a Linux operating system for both methods. The computational time can be considerably reduced, e.g. by treating sparsity using more efficient techniques as explained in Section 4.2 and programming the algorithm in C. There are primal-dual interior-point methods that treat sparsity very efficiently but most of them specialized to block-angular linear programs . For different data but with the same dimension and structure we observed that the number of iterations does not vary much.

Conclusions

A new decomposition method in convex programming is developed in this paper using dual decomposition and interior-point framework. Our method combines the fast local convergence rates of the Newton method with the efficiency of structural optimization for solving separable convex programs. Although our algorithm resembles augmented Lagrangian methods, it differs both in the computational steps and in the choice of the parameters. Contrary to most augmented Lagrangian methods that use gradient based directions to update the Lagrange multipliers, our method uses Newton directions and thus the convergence rate of the proposed method is faster. The reason for this lies in the fact that by adding self-concordant barrier terms to the standard Lagrangian we proved that under appropriate conditions the corresponding family of augmented dual functions is also self-concordant. Another appealing theoretical advantage of our interior-point Lagrangian decomposition method is that it is fully automatic, i.e. the parameters of the scheme are chosen as in the path-following methods, which are crucial for justifying its global convergence and polynomial-time complexity.

References