Scalable Preconditioning of Block-Structured Linear Algebra Systems using ADMM

TL;DR

This paper introduces a scalable iterative method combining ADMM and GMRES to efficiently solve large, block-structured linear systems common in optimization, outperforming traditional Schur complement methods in speed and robustness.

Contribution

The paper presents a novel ADMM-preconditioned GMRES approach that overcomes scalability issues in solving large, coupled linear systems in optimization problems.

Findings

ADMM-GMRES is nearly ten times faster than Schur decomposition.

The method handles systems with up to 2 million variables efficiently.

It is robust to the choice of the augmented Lagrangian penalty parameter.

Abstract

We study the solution of block-structured linear algebra systems arising in optimization by using iterative solution techniques. These systems are the core computational bottleneck of many problems of interest such as parameter estimation, optimal control, network optimization, and stochastic programming. Our approach uses a Krylov solver (GMRES) that is preconditioned with an alternating method of multipliers (ADMM). We show that this ADMM-GMRES approach overcomes well-known scalability issues of Schur complement decomposition in problems that exhibit a high degree of coupling. The effectiveness of the approach is demonstrated using linear systems that arise in stochastic optimal power flow problems and that contain up to 2 million total variables and 4,000 coupling variables. We find that ADMM-GMRES is nearly an order of magnitude faster than Schur complement decomposition. Moreover, we demonstrate that the approach is robust to the selection of the augmented Lagrangian penalty parameter, which is a key advantage over the direct use of ADMM.

Figures

Tables

Table 1. Table 1: Performance for ADMM-GMRES (ADGM), ADMM, and Schur decomposition (Schur) on problem ( 22 ) with n z ≥ 100 subscript 𝑛 𝑧 100 n_{z}\geq 100 .

Pglib-matpower Case	$n_x$	$n_z$	Solution Time (s)			Iteration Count		$\Vert Hu-r\Vert$
			Schur	ADMM	ADGM	ADMM	ADGM	Schur	ADMM	ADGM
P1. case13659_pegase	2115450	4091	1670.692	723.158	183.298	2000	325	2.798E-09	3.466E-03	1.602E-08
P2. case9241_pegase	1408950	1444	369.058	507.632	81.278	2000	165	3.841E-09	3.869E-05	2.981E-09
P3. case6470_rte	849800	760	124.416	336.881	51.106	2000	139	2.981E-10	3.932E-03	1.984E-10
P4. case6515_rte	845950	683	112.975	333.421	40.883	2000	109	3.206E-10	7.631E-04	2.627E-10
P5. case6495_rte	843650	679	111.035	320.943	37.302	2000	108	2.427E-10	7.128E-04	4.291E-10
P6. case2868_rte	389850	560	49.152	174.801	23.113	2000	105	3.228E-10	4.195E-04	1.404E-10
P7. case2848_rte	382250	510	43.584	168.000	22.508	2000	100	1.631E-10	2.329E-04	4.032E-11
P8. case2869_pegase	423500	509	48.424	189.774	38.293	2000	154	4.435E-10	1.135E-03	6.830E-10
P9. case6468_rte	813250	398	62.343	326.375	24.449	2000	65	2.774E-10	5.404E-04	4.416E-10
P10. case3012wp_k	367050	378	32.732	174.999	18.422	2000	74	7.032E-11	1.356E-05	1.338E-11
P11. case1951_rte	263900	365	24.864	132.780	19.566	2000	80	2.764E-10	2.836E-05	3.227E-10
P12. case2383wp_k	295900	319	24.458	152.102	17.499	2000	58	6.209E-11	7.444E-07	1.632E-10
P13. case1888_rte	249900	289	19.114	129.940	16.811	2000	66	1.534E-10	2.708E-06	2.935E-10
P14. case3120sp_k	367900	272	23.729	175.397	12.998	2000	56	3.402E-11	1.214E-08	2.220E-10
P15. case1354_pegase	193200	259	15.185	112.178	6.520	2000	47	1.453E-10	3.649E-07	2.728E-10
P16. case240_pserc	48650	142	5.262	64.665	4.737	2000	59	3.942E-10	1.902E-08	7.815E-10
P17. case2746wp_k	311600	103	8.455	100.004	11.080	1287	35	4.171E-11	9.925E-10	4.428E-10

Table 2. Table 2: Performance for ADMM-GMRES (ADGM), ADMM and Schur decomposition (Schur) on problem ( 22 ) with n z ≤ 100 subscript 𝑛 𝑧 100 n_{z}\leq 100 .

Pglib-matpower Case	$n_x$	$n_z$	Solution Time (s)			Iteration Count		$\Vert Hu-r\Vert$
			Schur	ADMM	ADGM	ADMM	ADGM	Schur	ADMM	ADGM
P18. case73_ieee_rts	19200	95	3.161	21.751	2.914	772	38	1.001E-11	9.795E-10	6.237E-11
P19. case2746wop_k	311100	84	7.033	78.842	10.125	1022	34	2.763E-11	9.911E-10	1.306E-09
P20. case2736sp_k	308400	81	6.758	81.199	9.796	1054	35	2.538E-11	9.978E-10	4.499E-10
P21. case300_ieee	41200	56	2.212	41.976	6.157	1305	38	6.341E-11	9.891E-10	6.844E-10
P22. case2737sop_k	305650	53	4.668	54.118	9.791	696	34	1.808E-11	9.905E-10	3.925E-10
P23. case200_pserc	26000	37	1.440	13.064	3.150	428	28	3.524E-12	9.645E-10	2.619E-11
P24. case24_ieee_rts	6250	31	1.096	6.219	2.059	231	30	3.218E-12	9.918E-10	5.142E-11
P25. case118_ieee	17050	18	0.799	10.767	2.212	375	32	1.329E-11	9.662E-10	2.874E-10
P26. case162_ieee_dtc	23450	11	0.618	4.629	1.829	148	20	7.489E-12	8.752E-10	1.347E-10
P27. case89_pegase	16100	11	0.588	3.927	1.936	132	23	4.015E-12	9.344E-10	1.228E-10
P28. case39_epri	5200	9	0.476	3.471	1.516	124	18	8.340E-12	8.749E-10	1.543E-10
P29. case30_as	4100	5	0.360	2.386	0.573	81	11	1.632E-13	9.865E-09	5.271E-11
P30. case30_fsr	4100	5	0.362	2.047	0.577	73	11	1.997E-13	8.138E-09	3.175E-11
P31. case5_pjm	1000	4	0.319	2.707	0.474	103	8	2.434E-13	9.123E-09	3.137E-10
P32. case57_ieee	7200	3	0.314	2.382	0.747	82	7	7.034E-13	9.208E-09	2.171E-09
P33. case14_ieee	1850	1	0.237	1.343	0.346	48	3	1.355E-13	9.339E-09	3.045E-11
P34. case30_ieee	3700	1	0.243	1.357	0.354	47	3	5.790E-14	9.653E-09	6.383E-12
P35. case3_lmbd	450	1	0.235	1.233	0.342	44	3	1.515E-13	7.080E-09	1.520E-12

Equations

$$\underset{x_i,z}{\min }$$

$$J_i{x}_i=b_i,({\lambda }_i)i\in \mathcal{P}$$

$$A_i{x}_i+B_{i}z=0,(y_i)i\in \mathcal{P}.$$

$$\mathcal{L}(x,z,\lambda ,y)=\sum _{i\in \mathcal{P}}\frac{1}{2}{x}_i^T{D}_i{x}_i+c_i^T{x}_i+y_i^T(A_i{x}_i+B_{i}z)+{\lambda }_i^T(J_i{x}_i-b_i),$$

$${\nabla }_{x_i}\mathcal{L}$$

$${\nabla }_{{\lambda }_i}\mathcal{L}$$

$${\nabla }_{y_i}\mathcal{L}$$

$${\nabla }_z\mathcal{L}$$

$$\underset{H}{\underbrace{\left[\begin{array}{ccc}K& & A^T\\ & & B^T\\ A& B& \end{array}\right]}}\underset{u}{\underbrace{\left[\begin{array}{c}\upsilon \\ z\\ y\end{array}\right]}}=\underset{r}{\underbrace{\left[\begin{array}{c}\gamma \\ 0\\ 0\end{array}\right]}}.$$

$$K_i=\left[\begin{array}{cc}{D}_i& J_i^T\\ J_i& \end{array}\right]{\tilde{A}}_i=\left[\begin{array}{cc}{A}_i& 0\end{array}\right],i\in \mathcal{P}.$$

$$\left[\begin{array}{ccccccccc}{D}_1& J_1^T& & & & & A_1^T& & \\ J_1& & & & & & & & \\ & & \ddots & & & & & \ddots & \\ & & & D_P& J_P^T& & & & A_P^T\\ & & & J_P& & & & & \\ & & & & & & B_1^T& \cdots & B_P^T\\ A_1& & & & & B_1& & & \\ & & \ddots & & & ⋮& & & \\ & & & A_P& & B_P& & & \end{array}\right]\left[\begin{array}{c}{x}_1\\ {\lambda }_1\\ ⋮\\ x_P\\ {\lambda }_P\\ z\\ y_1\\ ⋮\\ y_P\end{array}\right]=\left[\begin{array}{c}-c_1\\ b_1\\ ⋮\\ -c_P\\ b_P\\ 0\\ 0\\ ⋮\\ 0\end{array}\right]$$

$$\left[\begin{array}{ccccccccc}{D}_1& J_1^T& A_1^T& & & & & & \\ J_1& & & & & & & & \\ A_1& & & & & & & & B_1\\ & & & \ddots & & & & & ⋮\\ & & & & D_P& J_P^T& A_P^T& & \\ & & & & J_P& & & & \\ & & & & A_P& & & & B_P\\ & & B_1^T& \cdots & & & B_P^T& & \end{array}\right]\left[\begin{array}{c}{x}_1\\ {\lambda }_1\\ y_1\\ ⋮\\ x_P\\ {\lambda }_P\\ y_P\\ z\end{array}\right]=\left[\begin{array}{c}-c_1\\ b_1\\ 0\\ ⋮\\ -c_P\\ b_P\\ 0\\ 0\end{array}\right]$$

$$\left[\begin{array}{cc}{K}_s& B_s\\ B_s^T& \end{array}\right]\left[\begin{array}{c}{\upsilon }_s\\ z\end{array}\right]=\left[\begin{array}{c}{\gamma }_s\\ 0\end{array}\right]$$

$$K_{s_i}=\left[\begin{array}{ccc}{D}_i& J_i^T& A^T\\ J_i& & \\ A& & \end{array}\right]B_{s_i}={\left[\begin{array}{ccc}0& 0& B_i^T\end{array}\right]}^T,i\in \mathcal{P}.$$

$$(B_s^T{K}_s^{-1}{B}_s)z=B_s^T{K}_s^{-1}{\gamma }_s.$$

$$K_s{\upsilon }_s={\gamma }_s-B_{s}z.$$

$${\mathcal{L}}_{\rho }(x,z,\lambda ,y)=\sum _{i\in \mathcal{P}}{x}_i^T{D}_i{x}_i+c_i^T{x}_i+{\left(A_i{x}_i+B_{i}z\right)}^T{y}_i+{\lambda }_i^T(J_i{x}_i-b_i)+\frac{\rho }{2}\Vert A_i{x}_i+B_{i}z{\Vert }^2$$

$$K_{{\rho }_i}=\left[\begin{array}{cc}{D}_i+\rho A_i^T{A}_i& J_i^T\\ J_i& 0\end{array}\right].$$

$$\underset{x_i,z}{\min }$$

$$J_i{x}_i=b_i,({\lambda }_i)i\in \mathcal{P}$$

$$A_i{x}_i+B_{i}z=0,(y_i)i\in \mathcal{P}.$$

$$\underset{H_{\rho }}{\underbrace{\left[\begin{array}{ccc}{K}_{\rho }& \rho A^{T}B& A^T\\ \rho B^{T}A& \rho B^{T}B& B^T\\ A& B& \end{array}\right]}}\underset{u}{\underbrace{\left[\begin{array}{c}\upsilon \\ z\\ y\end{array}\right]}}=\underset{r}{\underbrace{\left[\begin{array}{c}\gamma \\ 0\\ 0\end{array}\right]}}.$$

$$\underset{H_{\rho }}{\underbrace{\left[\begin{array}{ccc}{K}_{\rho }& \rho A^{T}B& A^T\\ \rho B^{T}A& \rho B^{T}B& B^T\\ A& B& \end{array}\right]}}=\underset{M_{\rho }}{\underbrace{\left[\begin{array}{ccc}{K}_{\rho }& & \\ \rho B^{T}A& \rho B^{T}B& \\ A& B& -\frac{1}{\rho }I\end{array}\right]}}-\underset{N_{\rho }}{\underbrace{\left[\begin{array}{ccc}& -\rho A^{T}B& -A^T\\ & & -B^T\\ & & -\frac{1}{\rho }I\end{array}\right]}}$$

$$T_{\rho }(u):=G_{\rho }u+f_{\rho }$$

$$M_{\rho }^{-1}{H}_{\rho }u=M_{\rho }^{-1}r.$$

$$M_{\rho }^{-1}{H}_{\rho }h$$

$$=h-\left[G_{\rho }h+f_{\rho }-f_{\rho }\right]$$

$$=h-G_{\rho }h$$

$$=(I-G_{\rho })h.$$

$$\underset{x_i,z}{\min }$$

$$J{x}_i=b_i,({\lambda }_i)i\in \mathcal{P}$$

$$A_i{x}_i-z=0,(y_i)i\in \mathcal{P}$$

$$\underset{x_i,z}{\min }$$

$$\sum _{j\in {\Omega }_{G_i}}{P}_{G_{j,s}}-P_{L_{i,s}}=\sum _{j\in {\Omega }_i}{P}_{F_{i,j,s}},i\in {\Omega }_{\mathcal{B}},s\in \mathcal{P}$$

$$P_{F_{i,j,s}}=\frac{{\theta }_{i,s}-{\theta }_{j,s}}{X_{i,j}},i,j\in {\Omega }_L,s\in \mathcal{P}$$

$$P_{G_{j,s}}-z_j=0,j\in {\Omega }_G,s\in \mathcal{P}$$

$${\theta }_{0,s}=0,s\in \mathcal{P}.$$

Full text

Keywords: Schur complement decomposition; ADMM; iterative; linear algebra; large-scale

1 Introduction

The scalability of optimization solvers relies quite heavily on the solution of the underlying linear algebra systems. Advances in direct sparse linear algebra solvers have been instrumental in the widespread use of quadratic programming and nonlinear programming solvers such as Ipopt, OOQP, and Knitro [27, 1, 2]. Specialized direct solution techniques have also been developed to tackle large-scale and block-structured systems (using variants of Schur complement decomposition techniques) [38, 37, 10, 20, 21]. Block structures appear in many important applications such as parameter estimation, stochastic programming, network optimization, and optimal control. Schur decomposition techniques can also leverage parallel computing architectures and have enabled the solution of problems with millions to billions of variables and constraints. Unfortunately, many applications of interest still remain inaccessible due to fundamental scalability limitations of Schur complement techniques. Specifically, Schur complement decomposition does not scale well in problems that exhibit high degrees of block coupling. This is because high degrees of coupling require assembling and factorizing large Schur complement matrices (which are often highly dense).

Iterative solution techniques [33, 31, 13, 12, 6, 5] and associated preconditioning strategies [8, 32, 19, 35, 17, 43] have been proposed to address fundamental scalability issues of direct linear algebra strategies. In the context of block-structured problems, attempts have been made to solve the Schur complement system by using iterative solution techniques (to avoid assembling and factorizing the Schur complement). Preconditioners for Schur complements arising in special problem classes such as multi-commodity networks and stochastic programs have been developed [8]. Unfortunately, preconditioning strategies for general problem classes are still lacking. Another important issue that arises in this context is that the implementation of advanced linear algebra strategies is non-trivial (e.g., it requires intrusive modifications of optimization solvers).

Along a separate line of research, significant advances have been made in the development of problem-level decomposition techniques such as the alternating direction method of multipliers (ADMM) and Lagrangian dual decomposition [23, 22, 26, 23, 18, 25, 34]. Such approaches are flexible and rather easy to implement but suffer from slow convergence. Recently, it has been proposed to use ADMM as a preconditioner for Krylov-based iterative solvers such as GMRES [40, 41]. In this work, we provide a detailed derivation of this ADMM-GMRES approach and test its performance in the context of block-structured linear algebra systems. We demonstrate that this approach overcomes the scalability issues of Schur complement decomposition. We also demonstrate that this approach is significantly more effective than using ADMM directly. Our tests are facilitated by the use of PyNumero, a recently-developed Python framework that enables the implementation and benchmarking of optimization algorithms. We use the proposed framework to tackle problems with hundreds of thousands to millions of variables generated from standard benchmark sets and power grid applications.

The paper is structured as follows. In Section 2 we define the problem of interest and provide preliminary information on the use of Schur complement decomposition and ADMM approaches. In Section 3 we provide a detailed derivation of the ADMM-GMRES approach and in Section 4.1 we provide benchmark results.

2 Preliminaries

We study the solution of block-structured quadratic programs (QP) of the form:

[TABLE]

Here, $\mathcal{P}:=\{1,\dots,P\}$ is a set of block variable partitions. Each partition contains a vector of primal variables $x_{i}\in\mathbb{R}^{n_{x_{i}}}$ and the vector $z\in\mathbb{R}^{n_{z}}$ contains the primal variables that the couple partitions. The total number of primal variables is $n:=n_{z}+\sum_{i\in\mathcal{P}}n_{x_{i}}$. Equation (1b) are the partition constraints with their respective dual variables $\lambda_{i}\in\mathbb{R}^{m_{i}}$. Equation (1c) are the constraints that link partitions across set $\mathcal{P}$ and have associated dual variables $y_{i}\in\mathbb{R}^{l_{i}}$. We assume that the partition matrices $J_{i}\in\mathbb{R}^{m_{i}\times n_{i}}$ have full rank and that the right-hand-side coefficients $b_{i}\in\mathbb{R}^{m_{i}}$ are in the column space of $J_{i}$. The total number of partition constraints is $m:=\sum_{i\in\mathcal{P}}m_{i}$. We refer to $A_{i}\in\mathbb{R}^{n_{z}\times n_{i}}$ and $B_{i}\in\mathbb{R}^{n_{z}\times n_{z}}$ as linking matrices and we assume them to have full rank. The total number of linking constraints is $l:=\sum_{i\in\mathcal{P}}l_{i}$. The QP under study is the main computational kernel behind nonlinear programming strategies because it is used for computing primal-dual search steps.

We make the blanket assumption that the block-structured QP is strongly convex and that the combined Jacobian matrix (obtained by assembling partition and coupling constraints) has full rank. Strong convexity can be obtained by ensuring that all block Hessian matrices $D_{i}$ are positive definite. Strong convexity and full-rank conditions guarantee that the primal-dual solution of the QP exists and is unique. Moreover, these assumptions guarantee that the QP solution is a unique minimizer and that this can be found by solving the first-order stationarity conditions. Additional assumptions will also be needed on the nature of the building blocks of the QP (associated with each partition). Such assumptions are needed to ensure that proposed decomposition schemes are well-defined and will be stated as we proceed (in order to maintain clarity in the presentation).

The Lagrange function of (1) can be expressed as:

[TABLE]

where $x:=(x_{1},\cdots,x_{P})$, $\lambda:=(\lambda_{1},\cdots,\lambda_{P})$ and $y:=(y_{1},\cdots,y_{P})$. The first-order optimality conditions of (1) are given by:

[TABLE]

These conditions form a block-structured linear system of the form shown in (6). For the sake of compactness and ease of notation, we rewrite (6) as:

[TABLE]

where $\upsilon=(\upsilon_{1},\cdots,\upsilon_{P})$, $\upsilon_{i}=(x_{i},\lambda_{i})$, $\gamma=(\gamma_{1},\cdots,\gamma_{P})$, $\gamma_{i}=(-c_{i},b_{i})$, $u=(v,z,y)$, $r=(\gamma,0,0)$, and $K=\text{blkdiag}\{K_{1},K_{2},\cdots,K_{P}\}$. We also have $A=\text{blkdiag}\{\tilde{A}_{1},\tilde{A}_{2},\cdots,\tilde{A}_{P}\}$, $B=\text{rowstack}\{B_{1},B_{2},\cdots,B_{P}\}$ with:

[TABLE]

[TABLE]

2.1 Solution using Schur Decomposition

One can solve large instances of the block-structured QP by using a Schur-complement decomposition method (we refer to this approach simply as Schur decomposition) [39]. This approach decomposes (1) by using block Gaussian elimination on a permuted version of the linear system (6). The permuted system has the structure:

[TABLE]

This system can be expressed in compact form as:

[TABLE]

where $\upsilon_{s}=(\upsilon_{s_{1}},\cdots,\upsilon_{s_{P}})$, $\upsilon_{s_{i}}=(x_{i},\lambda_{i},y_{i})$, $\gamma_{s}=(\gamma_{s_{1}},\cdots,\gamma_{s_{P}})$, and $\gamma_{i}=(-c_{i},b_{i},0)$. We also have $K_{s}=\text{blkdiag}\{K_{s_{1}},K_{s_{2}},\cdots,K_{s_{P}}\}$ and $B_{s}=\text{rowstack}\{B_{s_{1}},B_{s_{2}},\cdots,B_{s_{P}}\}$ with:

[TABLE]

Because we have assumed that $D_{i}$ is positive definite and the combined constraint Jacobian $(J_{i},A)$ has full rank, we have that $K_{s_{i}}$ is nonsingular. This also implies $K_{s}$ is nonsingular and, as a result, we can form the Schur complement system:

[TABLE]

We refer to the coefficient matrix of (10) as the Schur complement. Since $K_{s}$ is block-diagonal, one can assemble the Schur complement by factorizing the blocks $K_{s_{i}}$ independently. By using this assembled Schur complement, one can then factorize the Schur complement matrix and solve system (10) to find a solution for the coupling variables $z$. Having $z$, one then proceeds to find solutions for the partition variables $\upsilon_{s}$ by solving the following system:

[TABLE]

Here, again, one can solve for each element $\upsilon_{s_{i}}$ independently because $K_{s}$ is block-diagonal. The Schur decomposition method is summarized in Algorithm 1. We refer the reader to [29] for details on the implementation of Schur decomposition approaches.

Schur decomposition is a flexible approach that enables the solution of problems with many block partitions. A fundamental limitation of this approach, however, is that one needs to assemble and factorize the Schur complement matrix $S$ (which is often highly dense). As a result, Schur decomposition does not scale well with the number of coupling variables $z$. Iterative approaches can in principle be used to solve the Schur system but effective preconditioning strategies do not currenlty exist for general block-structured systems.

2.2 Solution using ADMM

The block-structured QP can also be decomposed and solved by using ADMM. This approach seeks to minimize the augmented Lagrangian function:

[TABLE]

by performing alternating minimization with respect to the block variables $(x_{i},y_{i})$ and the coupling variables $z$. A standard implementation of ADMM for solving the QP of interest is presented in Algorithm 2.

In the above algorithm, $\mathcal{X}_{i}=\{x\,|\,J_{i}x-b_{i}=0\}$ is used to denote the feasible set of each partition (the inner block constraints are satisfied exactly). The ADMM algorithm can be implemented by solving the first-order conditions of each subproblem directly. This is because each block subproblem is strongly convex and the block Jacobian has full rank. This approach is sketched in Algorithm 3.

In the above algorithm we have that $K_{\rho}=\text{blkdiag}\{K_{\rho_{1}},K_{\rho_{2}},\cdots,K_{\rho_{P}}\}$ with

[TABLE]

We note that the update of the coupling variables still requires forming and factorizing the matrix $B^{T}B$ (but this can be done in blocks by forming and factorizing $B_{i}^{T}B_{i}$ individually). Moreover, this operation only needs to be performed once. Note also that $B^{T}B$ is invertible since $B$ has full rank. As a result, the update step for the coupling variables in ADMM is cheaper than that of Schur decomposition and can thus overcome the main computational bottleneck of the latter. Unfortunately, it is well-known that ADMM exhibits slow convergence and thus the ability to perform fast operations might be shadowed by the need to perform many iterations.

3 Solution using ADMM-GMRES

The key observation that motivates our work is that ADMM can be used as a preconditioner for iterative linear algebra techniques such as GMRES [40, 41]. To derive the ADMM preconditioning strategy, we consider the regularized QP (1):

[TABLE]

The solution of this problem is also a solution of (1) (since the penalization term vanishes at the solution). The optimality conditions of the regularized QP are given by:

[TABLE]

We refer to (15) as the KKT system and to $H_{\rho}$ as the KKT matrix. ADMM can be interpreted as a Gauss-Seidel (alternating) minimization of the block and coupling variables and the dual variables [7, 11, 34]. This induces a splitting operator $H_{\rho}=M_{\rho}-N_{\rho}$ satisfying:

[TABLE]

Applying splitting (16) to (15) gives the operator:

[TABLE]

where $G_{\rho}{=}M_{\rho}^{-1}N_{\rho}$ and $f_{\rho}{=}M_{\rho}^{-1}r$. Note that any $u$ satisfying the fixed-point $T_{\rho}(u){=}u$ also satisfies $(I-G_{\rho})u{=}f_{\rho}$ and is a solution of the preconditioned KKT system:

[TABLE]

This follows from $M_{\rho}^{-1}H_{\rho}=M_{\rho}^{-1}\left(M_{\rho}-N_{\rho}\right)=I-G_{\rho}$. This motivates the development of a Richardson recursion of the form $u^{k+1}{=}G_{\rho}u^{k}+f_{\rho}$, which converges to a $u$ satisfying $(I-G_{\rho})u{=}f_{\rho}$ and $M_{\rho}^{-1}H_{\rho}u{=}M_{\rho}^{-1}r$ (provided that the eigenvalues of $G_{\rho}$ are inside the unit circle). In Appendix A we show that the operator $T_{\rho}(u)$ can be computed by performing one ADMM iteration (using $u$ as starting point). In other words, we have that $T_{\rho}(u){=}$ADMM($u$, $N{=}1$, $\rho$). This also implies that the Richardson recursion can be written as $u^{k+1}{=}$ADMM($u^{k}$, $N{=}1$, $\rho$). Consequently, the Richardson recursion (and thus ADMM) are consistent preconditioner choices.

The key idea behind ADMM-GMRES is to solve the system $M_{\rho}^{-1}H_{\rho}u{=}M_{\rho}^{-1}r$ by using the Krylov solver GMRES. This is equivalent to solving $(I-G_{\rho})u{=}f_{\rho}$. The right-hand side of this system can be computed as $f_{\rho}=T_{\rho}(0)$. GMRES requires the computation of matrix-vector products with the preconditioned coefficient matrix of the form $M_{\rho}^{-1}H_{\rho}h{=}(I-G_{\rho})h$. This can be done by using the operator (17) as:

[TABLE]

This follows from the observation that:

[TABLE]

From (3) we note that asking the ADMM oracle ADMM($h$, $N$, $\rho$) to iterate until reaching convergence will deliver $T_{\rho}(h){=}h$ satisfying $M_{\rho}^{-1}H_{\rho}h{=}f_{\rho}$. In such a case, the ADMM preconditioner is perfect (since it solves the actual preconditioned KKT system). Consequently, the quality of the ADMM preconditioner will improve as one increases $N$. For details on the properties of the preconditioner, we refer the reader to [41, 40]. The ADMM-GMRES strategy is summarized in Algorithm 4.

4 Numerical Results

In this section we discuss the implementation of ADMM-GMRES and present results for different benchmark problems. All numerical experiments were performed using PyNumero, which is an open-source framework written in Python and C++ that combines modeling capabilities of the algebraic modeling language Pyomo [24] with efficient libraries like the AMPL solver library [14], the Harwell Subroutine Library (HSL), and NumPy/SciPy [28]. It uses object-oriented principles that facilitate the implementation of algorithms and problem formulations that exploit block-structures via polymorphism and inheritance. All these features facilitate the implementation of ADMM, Schur decomposition, and ADMM-GMRES. The optimization models were implemented in Pyomo/PyNumero and all linear algebra operations were performed in compiled code. Within PyNumero, we used an MA27 interface to perform all direct linear algebra operations. We use the GMRES implementation available in Scipy to perform all iterative linear algebra operations. Iterative linear algebra routines available in KRYPY [15] were also used to validate results. To implement the power grid models we used EGRET 222https://github.com/grid-parity-exchange/Egret, a Pyomo-based package that facilitates the formulation if optimization problems that arise in power systems. The convergence criterion for GMRES and ADMM requires that the norm of the KKT system residual $Hu-r$ is smaller than $\epsilon=1\cdot 10^{-8}$. If the convergence criterion is not satisfied after 2,000 iterations, the algorithm was aborted and we report the final residual achieved. The linear solver MA27 was used with a pivoting tolerance of $1\cdot 10^{-8}$

4.1 Standard Benchmark Problems

We first conducted tests with randomly generated instances to study qualitatively the performance of ADMM-GMRES on block-structured optimization problems. This section focuses on two-stage stochastic programs of the form:

[TABLE]

where $\mathcal{P}$ is the scenario set, $x_{i}$ are the second-stage (recourse) variables, and $z$ are first-stage (coupling) variables. We defined a nominal vector $b$ and create scenarios with right-hand-side vector $b_{i}$ using the nominal vector $b$ as the mean and a standard deviation $\sigma{=}0.5b$. We first demonstrate the scalability of Algorithm 4 when solving instances of problem (21) with high dimensionality in the coupling variables $z$. The stochastic problem was constructed in the following manner: $D$ was set to a $4{,}800\times 4{,}800$ block diagonal matrix with $16$ dense symmetric blocks. Each block was generated following Algorithm 14 from [41] with log-standard-deviation $s{=}0.5$ (see [41] for details). The random matrix $J$ has dimensions of $100\times 4{,}800$. The number of scenarios was set to $50$, giving an initial problem with $240{,}000$ variables and $5{,}000$ constraints. To investigate the scalability of Algorithm 4, the number of complicating variables was varied from $100$ to $4{,}000$. Note that, as $n_{z}$ increases, the number of constraints of (21) increases as $50n_{z}$. The largest problem solved contained $244{,}000$ total variables and $205{,}000$ total constraints.

We solved (21) using four different strategies. Figure 1 summarizes the results obtained with Schur decomposition, GMRES (without preconditioner), ADMM, and ADMM-GMRES. These results confirm the observations of Section 2. Specifically, Schur decomposition does not scale well as $n_{z}$ increases. The main reason for this behavior is that, as $n_{z}$ increases, the number of operations required to form the Schur-complement increase. In addition, because the Schur-complement matrix is a dense $n_{z}\times n_{z}$ matrix, the factorization time increases cubically as $n_{z}$ increases. We observe that GMRES takes the longest time to solve the problem. For ADMM and ADMM-GMRES we see more favorable scalability as $n_{z}$ increases. Specifically, we see that ADMM-GMRES converges faster and that the savings increase as the number of coupling variables increases. We also note that ADMM-GMRES mimics the trend in performance of ADMM but is significantly faster.

Figure 2 shows the residuals for the iterative approaches for a problem with $1{,}000$ complicating variables. We can see that all methods exhibit linear convergence but that ADMM-GMRES outperforms both ADMM and GMRES (unpreconditioned). Notably, ADMM-GMRES converges in just $35$ iterations while ADMM and GMRES require over $300$ iterations and $1{,}000$ iterations, respectively. Moreover, we note that ADMM-GMRES can reach high accuracy levels (of $1\cdot 10^{-8}$), which is a desirable feature of iterative solution strategies.

An important drawback of ADMM is the need to tune the penalty parameter $\rho$. The work in [16] shows that an optimal value for $\rho$ can be chosen based on the smallest and largest eigenvalues of the matrix $A^{T}D^{-1}A$. In principle, this selection of $\rho$ is optimal for ADMM-GMRES as well. However, for large-scale structured problems such as the ones considered here, computing the eigenvalues of $A^{T}D^{-1}A$ is expensive. Heuristic approaches have also been proposed to select $\rho$ at every ADMM iteration with the objective of accelerating convergence [36]. Unfortunately these heuristics do not provide guarantees and might incur additional overhead. In particular, for the QP problems considered here, varying $\rho$ at every ADMM iteration will require forming and factorizing the $K_{\rho_{i}}$ repetitively. Interestingly, we now proceed to demonstrate that ADMM-GMRES is fairly insensitive to the choice of $\rho$.

Figure 3 compares the performance of ADMM against that of ADMM-GMRES for a stochastic program with $n_{z}{=}1{,}000$ and different values of $\rho$ ranging from $10^{-3}$ to $10^{3}$. Here we measure performance in terms of the number of iterations. Our results on the block-structured problem are in agreement with those in [40], where a single block problem is solved. We see that ADMM-GMRES is remarkably robust to the choice of $\rho$ (the number of iterations remain below 30). ADMM, on the other hand, fails to converge within 2,000 iterations for small and large values of $\rho$. The superior performance of ADMM-GMRES is attributed to the fact that the selection of $\rho$ only has an effect on the preconditioner and not on the convergence properties of GMRES. Nevertheless, we do observe that an optimal selection of $\rho$ improves performance of both ADMM and ADMM-GMRES. The robustness of ADMM-GMRES is a desirable feature when using the solver within more advanced SQP-based solvers. In particular, recent developments of augmented lagrangian interior-point approaches [9, 30] provide promising frameworks for ADMM-GMRES because the selection of $\rho$ is made based on information from the outer-iteration of the interior-point [3, 4]. Finally, recent developments of inertia-free methods for nonconvex nonlinear optimization enable the use of iterative linear solvers for the Newton subproblem.

4.2 Optimal Power Flow Problems

We now demonstrate the computational benefits of using ADMM-GMRES by solving stochastic optimal power flow problems. The optimal power flow problem is frequently used in power networks to determine an efficient dispatch of power generators that satisfy demand and maintains feasible operation conditions. This approach assumes that demand forecasts are accurate and determines a nominal operation point for power generation, power flow in transmission lines and voltage angle at each bus in the power grid. We solve the DC-power flow problem for 35 benchmark problems available in the PGLIB_OPF library [42] and determine nominal operation points for each of them. The solution of each benchmark problem was also obtained using Ipopt.

To asses the computational performance of ADMM-GMRES, we formulated a set-point problem that uses the nominal solution of the DC-power flow problem but seeks to minimize the effect of potential uncertainty in electricity demand values. The optimization solved is the quadratic problem:

[TABLE]

Here, $\Omega_{\mathcal{B}}$ and $\Omega_{L}$ denote the set of buses and transmission lines in the network, $\Omega_{G}$ the set of generators, $\Omega_{G_{i}}$ the set of generators at bus $i$, and $\Omega_{i}$ the set of buses connected to bus $i$. The variables in the model are the generator outputs $P_{G}$ , the flow in the transmission lines $P_{F}$, and the voltage angles $\theta_{j}$. As parameters we have the reactance of the lines $X_{i,j}$, the loads $P_{L}$, and the set-point values $P_{G}^{{\dagger}}$, $P_{F}^{{\dagger}}$, and $\theta^{{\dagger}}$ obtained from the DC-power flow solution. We denote the reference bus as $\theta_{0}$ and define objective weight values $w$. The goal of formulation (22) is to find the closest feasible operation to the optimal DC-power flow solution while accounting for potential uncertainty in the demands. In this problem the first-stage variables are the output of the generators for which we use Equation (22d) to enforce the same power generation across the set of scenarios. The dimensionality of the first-stage of this problem is given by the number of generators in the power network. Hence, this number varies from three to $4{,}092$ in the $35$ different benchmark problems considered in our study. For each benchmark we generated $50$ random scenarios with normal random distributed noise on the load $P_{L}$

Tables 1 and 2 summarize the results for the $35$ benchmarks. The problems were sorted according to their number of coupling variables. Table 1 presents the results for benchmarks with a first-stage dimension greater than $100$. Results for the smaller benchmarks are shown in Table 2. We see that, for the smaller problems, Schur decomposition is the best alternative as it can give exact solutions in about the same time as ADMM. However, for all benchmarks shown in Table 1, ADMM-GMRES finds an $\epsilon$-accurate solution in less time than ADMM and Schur decomposition. In particular, for case13659_pegase (with $4{,}091$ first-stage variables), ADMM-GMRES solved the problem almost an order of magnitude faster than Schur decomposition. By observing the trend for the rest of the problems, we conclude that these favorable scalability results can be expected to hold as the number of coupling variables increases. We highlight that, for many of the problems shown in Table 1, ADMM does not converge after $2{,}000$ iterations; ADMM-GMRES, on the other hand, consistently achieves $\epsilon$-accurate solutions in few iterations and regardless of the choice of $\rho$. In summary, our results demonstrate that ADMM-GMRES provides a plausible approach to overcome the limitations of Schur complement decomposition.

[FIGURE:]

[FIGURE:]

5 Conclusions and Future Work

We have demonstrated that ADMM provides an effective mechanism to precondition iterative linear solvers and with this overcome scalability limitations of Schur complement decomposition. Our results also demonstrate that the approach is robust to the choice of the penalty parameter. As part of future work, we will investigate the performance of ADMM-GMRES within a nonlinear interior-point framework. Here, it will be necessary to relax our assumptions on strong convexity and on the full rank of the Jacobian. Preliminary results reported in the literature indicate that different types of primal-dual regularized KKT systems can be used to compute search steps within interior-point methods under such relaxed conditions [9]. For instance, the primal-dual regularized system correspond to the optimality conditions of the QP problem:

[TABLE]

We will investigate ADMM variants to precondition such systems. The effectiveness of using ADMM as a preconditioner makes us wonder whether other approaches can be used for preconditioning as well. For instance, inexact dual Newton strategies can potentially be used to precondition structured KKT systems. This is an interesting direction of future work. We will also investigate advanced ADMM strategies that use second-order multiplier updates to accelerate the preconditioner.

Disclaimer

Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia, LLC., a wholly owned subsidiary of Honeywell International, Inc., for the U.S. Department of Energy’s National Nuclear Security Administration under contract DE-NA-0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the U.S. Department of Energy or the United States Government.This work was conducted as part of the Institute for the Design of Advanced Energy Systems (IDAES) with funding from the Office of Fossil Energy, Cross-Cutting Research, U.S. Department of Energy. V. M. Zavala acknowledges funding from the National Science Foundation under award NSF- EECS-1609183.

Appendix A Computing Operator $T_{\rho}(u)$ using ADMM

Here we prove that the operator $T_{\rho}(u)$ can be computed by applying one ADMM iteration. Consider, without loss of generality (and in order to simplify the presentation), the case of a single block problem of the form:

[TABLE]

The results that we derive next can be extended to multiple blocks using induction. The KKT system for problem (24) is:

[TABLE]

where $K_{\rho}=D+\rho A^{T}A$. Applying a Gauss-Seidel splitting to this system at a point $u=(x,z,y)$ leads to the update $u^{+}=T_{\rho}(u)=G_{\rho}u+f_{\rho}$, where $G_{\rho}=M_{\rho}^{-1}N_{\rho}$ and $f_{\rho}=M_{\rho}^{-1}r$. The explicit form of $M_{\rho}^{-1}$ is given by:

[TABLE]

where $\Sigma:=B^{T}B$. Having $M_{\rho}^{-1}$ we construct:

[TABLE]

and the right-hand-side-vector

[TABLE]

By defining $Q:=AK_{\rho}^{-1}A^{T}$ we can write the explicit form of the update $u^{+}$ as:

[TABLE]

Upon expansion we obtain the update $u^{+}=(x^{+},z^{+},y^{+})$:

[TABLE]

We now show that ADMM delivers the same updates after one iteration. We use the augmented Lagrange function:

[TABLE]

Initializing at $u=(x,z,y)$, the update $x^{+}$ is given by:

[TABLE]

For which the optimality conditions are

[TABLE]

and thus,

[TABLE]

We note that (34) and (30a) are equivalent. The update for the coupling variables $z^{+}$ is given by:

[TABLE]

The optimality conditions are given by:

[TABLE]

and thus,

[TABLE]

We thus have that (37) and (30b) are equivalent. Finally, the dual variables are updated as $y^{+}=y+\rho(Ax^{+}+Bz^{+})$. Substituting (34) and (37) in this expression leads to (30c) .