Factorization of Saddle-point Matrices in Dynamical Systems Optimization---Reusing Pivots

TL;DR

This paper introduces a method for efficiently solving sequences of saddle-point systems in dynamical systems optimization by reusing pivots in LDL^T factorizations, leading to faster computations.

Contribution

The paper presents a novel approach to reuse pivots in LDL^T factorizations for saddle-point matrices, improving computational efficiency in dynamical systems optimization.

Findings

Method outperforms Bunch-Parlett in speed

Reuses pivots to reduce computation time

Maintains same accuracy as traditional methods

Abstract

In this paper we consider the application of direct methods for solving a sequence of saddle-point systems. Our goal is to design a method that reuses information from one factorization and applies it to the next one. In more detail, when we compute the pivoted $LDL^T$ factorization we speed up computation by reusing already computed pivots and permutations. We develop our method in the frame of dynamical systems optimization. Experiments show that the method improves efficiency over Bunch-Parlett while delivering the same results.

Figures

Tables

Table 1. Table 1: The results for the linear benchmark problem  [ 12 , Sec. 6.2] . On the left hand side there are the parameters, in the middle the results of Alg.  4 with no monitoring of pivots and on the right hand side the results of Alg.  4 with the monitoring.

$k$	$N$	$\#IT$	$\#LD{L}^T$	$\#B$-$P$/$B$-$K$	$\#IT$	$\#LD{L}^T$	$\#$upd of $P$	$R$
$10$	$5$	$47$	$5$	$42$	$47$	$5$	$1$	$42$
	$10$	$26$	$19$	$7$	$26$	$19$	$1$	$7$
	$15$	$31$	$17$	$14$	$31$	$17$	$1$	$14$
	$20$	$400$	$3$	$418/408$	$111/258$	$3$	$2/5$	-
	$25$	$45$	$17$	$28$	$45$	$17$	$1$	$28$
	$30$	$131/125$	$11$	$120/116$	$89/78$	$11$	$8/4$	-
$20$	$5$	$25$	$27$	$0$	$25$	$27$	$0$	-
	$10$	$27$	$27$	$0$	$27$	$27$	$0$	-
	$15$	$24$	$9$	$15$	$24$	$9$	$1$	$15$
	$20$	$23$	$3$	$20$	$23$	$3$	$1/3$	$20/6$
	$25$	$27$	$5$	$22$	$27$	$5$	$1/2$	$22/11$
	$30$	$34$	$20$	$14$	$34$	$20$	$1$	$14$
$30$	$5$	$36$	$36$	$0$	$36$	$36$	$0$	-
	$10$	$24$	$24$	$0$	$24$	$24$	$0$	-
	$15$	$26$	$10$	$16$	$26$	$10$	$1$	$16$
	$20$	$28$	$23$	$5$	$28$	$23$	$1$	$5$
	$25$	$24$	$10$	$14$	$24$	$10$	$1$	$14$
	$30$	$24$	$6$	$18$	$24$	$6$	$1$	$18$
$40$	$5$	$41$	$41$	$0$	$41$	$41$	$0$	-
	$10$	$24$	$24$	$0$	$24$	$24$	$0$	-
	$15$	$24$	$3$	$21$	$24$	$3$	$1/2$	$21/10$
	$20$	$20$	$5$	$15$	$20$	$5$	$1$	$15$
	$25$	$37$	$11$	$26$	$37$	$11$	$1$	$26$
	$30$	$59$	$9$	$50$	$59$	$9$	$3$	$16$

Table 2. Table 2: The results for the nonlinear benchmark problem  [ 12 , Sec. 6.3] . On the left hand side there are the parameters, in the middle the results of Alg.  4 with no monitoring of pivots and on the right hand side the results of Alg.  4 with the monitoring.

$k$	$N$	$\#IT$	$\#LD{L}^T$	$\#B$-$P$/$B$-$K$	$\#IT$	$\#LD{L}^T$	$\#$upd of $P$	$R$
$10$	$5$	$40$	$40$	$0$	$40$	$40$	$0$	-
	$10$	$246$	$175$	$71$	$246$	$175$	$1$	$71$
	$15$	$115$	$36$	$79$	$115$	$36$	$2/1$	$39/79$
	$20$	$98$	$17$	$81$	$98$	$17$	$1/2$	$81/40$
	$25$	$400$	$14$	$386$	$400$	$14$	$4/2$	$96/193$
	$30$	$168$	$11$	$157$	$168$	$11$	$1$	$157$
$20$	$5$	$67$	$17$	$50$	$67$	$17$	$1$	$50$
	$10$	$50$	$19$	$31$	$50$	$19$	$1$	$31$
	$15$	$97$	$42$	$55$	$97$	$42$	$1/2$	$55/27$
	$20$	$34$	$8$	$26$	$34$	$8$	$1/2$	$26/13$
	$25$	$289$	$10$	$279$	$289$	$10$	$2$	$139$
	$30$	$109$	$22$	$87$	$109$	$22$	$1$	$87$
$30$	$5$	$40$	$30$	$10$	$40$	$30$	$1$	$10$
	$10$	$47$	$30$	$17$	$47$	$30$	$1$	$17$
	$15$	$400$	$12$	$388$	$400$	$12$	$4$	$97$
	$20$	$51$	$3$	$48$	$51$	$3$	$1/3$	$48/16$
	$25$	$109/107$	$3$	$106/104$	$108/113$	$3$	$1/3$	-
	$30$	$400/177$	$20$	$380/157$	$400/374$	$20$	$2/4$	-
$40$	$5$	$113$	$49$	$64$	$113$	$49$	$1$	$64$
	$10$	$400$	$97$	$303$	$400$	$97$	$1$	$303$
	$15$	$210/211$	$10$	$200/201$	$207/213$	$10$	$1/2$	-
	$20$	$58$	$25$	$33$	$58$	$25$	$1$	$33$
	$25$	$134$	$1$	$133$	$134$	$1$	$1/3$	$133/44$
	$30$	$146$	$12$	$134$	$146/155$	$12$	$4/5$	-

Equations

$$K_i{x}_i=y_i,i=1,2,\dots$$

$$PK{P}^T$$

$$\dot{x}(t)=f(x(t)),x(0)=x_0,$$

$$\left[\begin{array}{cc}H& B\\ B^T& -C\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}f\\ g\end{array}\right]\text{or}Ku=b,$$

$$H=\left[\begin{array}{ccc}{H}_1& & \\ & \ddots & \\ & & H_N\end{array}\right],B=\left[\begin{array}{ccccccc}v& -M_1^T& & & & & \\ & -v_1^T& & & & & \\ & I& -M_2^T& & & & \\ & & -v_2^T& & & & \\ & & I& & & & \\ & & & \ddots & & & \\ & & & & & & \\ & & & & & -M_{N-1}^T& \\ & & & & & -v_{N-1}^T& \\ & & & & & I& w\\ & & & & & 0& \beta \end{array}\right].$$

$$\left[\begin{array}{cc}H& B\\ B^T& -C\end{array}\right]$$

$$=\left[\begin{array}{cc}{L}_H& 0\\ B^T{L}_H^{-T}{D}_H^{-1}& L_S\end{array}\right]\left[\begin{array}{cc}{D}_H& 0\\ 0& D_S\end{array}\right]\left[\begin{array}{cc}{L}_H^T& D_H^{-1}{L}_H^{-1}B\\ 0& L_S^T\end{array}\right]$$

$$=LD{L}^T.$$

$$L_H{D}_H{L}_H^T=\left[\begin{array}{ccc}{L}_{H,1}& & \\ & \ddots & \\ & & L_{H,N}\end{array}\right]\left[\begin{array}{ccc}{D}_{H,1}& & \\ & \ddots & \\ & & D_{H,N}\end{array}\right]\left[\begin{array}{ccc}{L}_{H,1}^T& & \\ & \ddots & \\ & & L_{H,N}^T\end{array}\right],$$

$$B^T{L}_H^{-T}{D}_H^{-1}=\left[\begin{array}{cccccc}{s}_1^T& & & & & \\ X_1& Y_2& & & & \\ & X_2& Y_3& & & \\ & & \ddots & & & \\ & & & X_{N-2}& Y_{N-1}& \\ & & & & X_{N-1}& Y_N\\ & & & & & s_2^T\end{array}\right],$$

$$s_1$$

$$X_i$$

$$Y_i$$

$$s_2$$

$$S=-\left[\begin{array}{ccccccc}{\alpha }_1& w_1^T& & & & & \\ w_1& V_1& W_1^T& & & & \\ & W_1& V_2& & & & \\ & & & \ddots & & & \\ & & & & V_{N-2}& W_{N-2}^T& \\ & & & & W_{N-2}& V_{N-1}& w_2\\ & & & & & w_2^T& {\alpha }_2\end{array}\right],$$

$${\alpha }_1$$

$$w_1$$

$$V_i$$

$$W_i$$

$$w_2$$

$${\alpha }_2$$

$$-S=\left[\begin{array}{cccc}{\alpha }_1& w_1^T& & \\ w_1& V_1& W_1^T& \\ & W_1& V_2& w_2\\ & & w_2^T& {\alpha }_2\end{array}\right]⟶\left[\begin{array}{cccc}1& & & \\ w_1/{\alpha }_1& {\hat{V}}_1& W_1^T& \\ & W_1& V_2& w_2\\ & & w_2^T& {\alpha }_2\end{array}\right],$$

$$\left[\begin{array}{cccc}1& & & \\ w_1/{\alpha }_1& {\hat{V}}_1& W_1^T& \\ & W_1& V_2& w_2\\ & & w_2^T& {\alpha }_2\end{array}\right]⟶\left[\begin{array}{cccc}1& & & \\ w_1/{\alpha }_1& {\hat{L}}_1& & \\ & W_1{\hat{L}}_1^{-T}{\hat{D}}_1^{-1}& {\hat{V}}_2& w_2\\ & & w_2^T& {\alpha }_2\end{array}\right],$$

$$\left[\begin{array}{cccc}1& & & \\ w_1/{\alpha }_1& {\hat{L}}_1& & \\ & W_1{\hat{L}}_1^{-T}{\hat{D}}_1^{-1}& {\hat{V}}_2& w_2\\ & & w_2^T& {\alpha }_2\end{array}\right]⟶\left[\begin{array}{cccc}1& & & \\ w_1/{\alpha }_1& {\hat{L}}_1& & \\ & W_1{\hat{L}}_1^{-T}{\hat{D}}_1^{-1}& {\hat{L}}_2& \\ & & w_2^T{\hat{L}}_2^{-T}{\hat{D}}_2^{-1}& {\hat{\alpha }}_2\end{array}\right],$$

$$L_S=\left[\begin{array}{cccccc}1& & & & & \\ l_1& {\hat{L}}_1& & & & \\ & W_1{\hat{L}}_1^{-T}{D}_{S,1}^{-1}& {\hat{L}}_2& & & \\ & & \ddots & & & \\ & & & W_{N-2}{\hat{L}}_{N-2}^{-T}{D}_{S,N-2}^{-1}& {\hat{L}}_{N-1}& \\ & & & & l_2^T& 1\end{array}\right],$$

$$-D_S=\left[\begin{array}{ccccc}{d}_1& & & & \\ & D_{S,1}& & & \\ & & \ddots & & \\ & & & D_{S,N-1}& \\ & & & & d_N\end{array}\right],$$

$$A=\left[\begin{array}{cc}{A}_{1,1}& A_{1,2}\\ A_{1,2}^T& A_{2,2}\end{array}\right],$$

$$A=\left[\begin{array}{cc}{I}_j& 0\\ A_{1,2}^T{A}_{1,1}^{-1}& I_{n-j}\end{array}\right]\left[\begin{array}{cc}{A}_{1,1}& 0\\ 0& A_{2,2}-A_{1,2}^T{A}_{1,1}^{-1}{A}_{1,2}\end{array}\right]\left[\begin{array}{cc}{I}_j& A_{1,1}^{-T}{A}_{1,2}\\ 0& I_{n-j}\end{array}\right],$$

$$m:=\underset{i}{\max }\{|m_i|;m_i\in M\}$$

$${\mu }_0^{(n-1)}:=\underset{i,j}{\max }\{|a_{i,j}|;a_{i,j}\in A^{(n-1)}\}$$

$$m:=\underset{i,j}{\max }\{|m_{i,j}|;m_{i,j}\in M\}$$

$${\mu }_0^{(n-2)}:=\underset{i,j}{\max }\{|a_{i,j}|;a_{i,j}\in A^{(n-2)}\}$$

$${\epsilon }_1<|\det \beta |$$

$$\dot{x}=\left[\begin{array}{ccc}\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]& & \\ & \ddots & \\ & & \left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]\end{array}\right]x,$$

$$\dot{x}=\left[\begin{array}{ccc}\left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]& & \\ & \ddots & \\ & & \left[\begin{array}{cc}0& 1\\ -1& 0\end{array}\right]\end{array}\right]x+\left[\begin{array}{c}\sin (x_k)\\ ⋮\\ \sin (x_1)\end{array}\right],$$

$$R:=\lfloor \frac{\text{\#B-P/B-K}}{\text{\#upd of P}}\rfloor .$$

Full text

Keywords: saddle-point matrix; symmetric indefinite factorization; dynamical systems; sequential quadratic programming

1 Introduction

Consider a sequence of saddle-point systems arising, for example, in Sequential Quadratic Programming (SQP) [16, 18], that is

[TABLE]

where $K_{i}\in\mathbb{R}^{n\times n}$, $x_{i}\in\mathbb{R}^{n}$ and $y_{i}\in\mathbb{R}^{n}$. Here, we assume matrices $K_{i}$ to have the same structure of nonzero entries. Based on this we study how to solve the system $K_{i}x_{i}=y_{i}$, benefiting from information from previous iterations. Our goal is to develop a strategy that reduces the amount of work spent on searching the matrix for pivots in direct methods. In more detail we describe our strategy in Section 6 and compare its stability with Bunch-Parlett [7].

We develop this strategy in the frame of direct methods for solving saddle-point systems that arise in a certain class of optimization problems arising in verification of dynamical systems [12]. Here, one seeks a solution of a dynamical system that originates in a given set of initial states and reaches another set of states that are to be avoided, unsafe states. Saddle-point matrices $K_{i}$ that arise have specific structure which we try to exploit.

Such optimization problem occur for example in control and verification of hybrid systems [5, 8, 17, 27] and in motion planning [14, 19]. In addition, the techniques described in this paper apply to general underdetermined boundary value problems for ordinary differential equations [1]. Moreover, similar saddle-point matrices to ours arise, for example, in the mixed and hybrid finite element discretizations [22], in a class of interior point methods [24, 25, 26], and time-harmonic Maxwell equations [10, 11].

The saddle-point matrices that arise in such optimization problems and applications are sparse, whether one uses the SQP method [16] or the interior-point method [15]. Hence, direct methods for solving the saddle-point system look promising. However, the naive application of straightforward $LDL^{T}$ factorization often results in a failure due to ill-conditioning [12, 24, 25, 26] and the singularity of the $(1,1)$ block [10, 11].

Denoting the saddle-point matrix and its factorization by

[TABLE]

in more detail, the main contributions of this paper are: description of a strategy for selecting and reusing pivots in the symmetric and indefinite $PKP^{T}=LDL^{T}$ factorization; analysis of the growth factor in the reduced matrices; numerical comparison with Bunch-Parlett and Bunch-Kaufman on a series of benchmarks from dynamical systems optimization; description and exploitation of a specific structure of the saddle-point matrix $K$ and the factor $L$ in dynamical systems optimization problem [13]; Alg. 4 that switches from unpivoted to pivoted factorization of the matrix $K$, balancing the speed with the stability of computation.

The outline of the paper is as follows. In Section 2 we briefly review the optimization problem we try to solve [12]. In Section 3 we describe the structure of the saddle-point matrices. In Section 4 we will compute the $LDL^{T}$ factorization and prove that the factor $L$ has a banded structure of nonzero elements. Then the discussion about the implementation of the $LDL^{T}$ factorization follows and a hybrid method for solving the saddle-point system is described in Section 5. Sections 6 contains a detailed description of reusing pivots and its effect on the stability of our method. Furthermore we include numerical results in Section 7. The whole paper is concluded with a summary and a brief discussion of results in Section 8.

2 Motivation

Our motivation originates from the field of computer aided verification [5, 14, 27]. Consider a system of ordinary differential equations such that

[TABLE]

where $x:\mathbb{R}\to\mathbb{R}^{k}$ is a function of variable $t\geq 0$, $x_{0}\in\mathbb{R}^{k}$ and $f:\mathbb{R}^{k}\to\mathbb{R}^{k}$ is continuously differentiable. We denote the flow of the vector field $f$ in (1) by $\Phi:\mathbb{R}\times\mathbb{R}^{k}\to\mathbb{R}^{k}$ and for the fixed $x_{0}$ one has the solution $x(t)$ of (1), where $x(t)=\Phi(t,x_{0})$ for $t\geq 0$.

Denote the set of initial states by $\mathrm{Init}$ and the set of states we try to avoid by $\mathrm{Unsafe}$. Our goal is to find any solution $x(t)$ of (1) such that $x_{0}\in\mathrm{Init}$ and $\Phi(t_{f},x_{0})\in\mathrm{Unsafe}$ for some $t_{f}>0$, if it exists.

In the previous work [12] we solve this boundary value problem by the multiple-shooting method [2]. That is, one computes a solution of (1) from shorter solution segments. Suppose we have $N$ solution segments of (1) such that their initial states are denoted by $x_{0}^{i}$ and their lengths by $t_{i}>0$ for $1\leq i\leq N$. Then the desired solution to our problem satisfies: $x_{0}^{1}\in\mathrm{Init}$, $x_{0}^{i+1}=\Phi(t_{i},x_{0}^{i})$ for $1\leq i\leq N-1$ (these are the matching conditions), and $\Phi(t_{N},x_{0}^{N})\in\mathrm{Unsafe}$.

Boundary conditions $x_{0}^{1}\in\mathrm{Init}$ and $\Phi(t_{N},x_{0}^{N})\in\mathrm{Unsafe}$ can be formulated either as equalities (points belong to the boundaries of the sets), or it can be given as inequalities (points are insides of sets). Either way there are infinitely many solutions [13], therefore, one needs to introduce a regularization. In the paper [13] we formulate an objective function in the form $\sum t_{i}^{2}$, where $t_{i}$ is the length of the $i$-th solution segment, that drives the solution segments to have the same lengths.

In the end one solves a general nonlinear programming problem with $N(k+1)$ parameters, where those parameters are lengths of solution segments $t_{i}>0$ and initial states $x_{0}^{i}\in\mathbb{R}^{k}$, $1\leq i\leq N$. From now onwards we denote the number of parameters by $n=N(k+1)$ and the number of constraints by $m=(N-1)k+2$.

3 Blocks of Saddle-point Matrix

The Line-search SQP method described in [12] requires in each iteration the solution of the saddle-point system (2)

[TABLE]

where $H\in\mathbb{R}^{n\times n}$, $B\in\mathbb{R}^{n\times m}$ with $n\geq m$ and $C\in\mathbb{R}^{m\times m}$.

The elements of the matrix $B$ and $C$ are recomputed and the matrix $H$ is updated block by block by the BFGS scheme. The structure of nonzero elements of the matrix $B$ in (3) remains the same throughout the iterations.

Since we are interested in the $LDL^{T}$ factorization of the saddle-point matrix $K$ that arises in the solution of reachability problems for dynamical systems [12], the blocks $H$ and $B$ have the form

[TABLE]

The matrix $H$ consists of blocks $H_{i}\in\mathbb{R}^{(k+1)\times(k+1)}$, $1\leq i\leq N$. The matrix $C=diag(\gamma_{1},0,\ldots,0,\gamma_{2})$, where $\gamma_{i}\geq 0$ for $i=1,2$. In the matrix $B$, there are blocks $M_{i}\in\mathbb{R}^{k\times k}$ and vectors $v_{i}\in\mathbb{R}^{k}$, $1\leq i\leq N-1$. The matrix $I\in\mathbb{R}^{k\times k}$ is the identity matrix of order $k$, the vectors $v$ and $w$ are nonzero and belong to $\mathbb{R}^{k}$, and $\beta$ is a nonzero scalar. The first and the last columns of the matrix $B$ correspond to the boundary constraints $x_{0}^{1}\in\mathrm{Init}$ and $\Phi(t_{N},x_{0}^{N})\in\mathrm{Unsafe}$.

A similar banded structure of nonzero elements to the one in the matrix $B$ arises, for example, when one solves boundary value problems for ordinary differential equations [1, Sec. 5]. However, we have additional entries $v_{i}$ in (3) because we consider the lengths of time intervals $t_{i}$ (the lengths of solution segments) to be parameters and not fixed values.

The saddle-point matrix $K$ satisfies the following conditions [12]: the matrix $H$ is symmetric positive definite (BFGS approximations of the Hessian), and the matrix $B$ has full column rank. Under these conditions the saddle point matrix $K$ is: nonsingular [4, Th. 3.1], indefinite [4, Th. 3.5], strongly factorizable [23, Th. 2.1].

4 $LDL^{T}$ Factorization

In this section we give formulas for the $LDL^{T}$ factorization of the saddle-point matrix $K$ with blocks (3). In addition, we describe the structure of nonzero elements of the unit lower triangular factor $L$ for which $K=LDL^{T}$.

A standard approach [9, Ch. 4.1] to solving the linear system $Ku=b$, where $K$ is symmetric and nonsingular, is Alg. 1.

Denote by $L_{H}D_{H}L_{H}^{T}$ the $LDL^{T}$ factorization of the matrix $H$ and by $L_{S}D_{S}L_{S}^{T}$ the $LDL^{T}$ factorization of the Schur complement $S=K/H$, where $S=-C-B^{T}H^{-1}B$. When forming the factor $L$ one factorizes $H=L_{H}D_{H}L_{H}^{T}$ and then $S=L_{S}D_{S}L_{S}^{T}$. Then

[TABLE]

Computation of the factor $L$ is summarized in Alg. 2 and follows the framework of the solution of equilibrium systems [9, p. 170].

Matrix $H$ from (3) is block diagonal, then

[TABLE]

where $L_{H,i}\in\mathbb{R}^{(k+1)\times(k+1)}$ is unit lower triangular and $D_{H,i}\in\mathbb{R}^{(k+1)\times(k+1)}$ is diagonal, $1\leq i\leq N$.

Lemma 1**.**

Let $K=LDL^{T}$ be the $LDL^{T}$ factorization of the saddle-point matrix (2) with blocks given in (3). Then the $(1,1)$ block $L_{H}$ of the factor $L$ is block diagonal and its blocks $L_{H,i}\in\mathbb{R}^{(k+1)\times(k+1)}$,$1\leq i\leq N$, are unit lower triangular.

Proof.

The matrix $H$ (3) is symmetric positive definite, therefore, each $H_{i}$, $1\leq i\leq N$, is symmetric positive definite. When one carries out the $LDL^{T}$ factorization of $H_{i}$, then one obtains factors $L_{H,i}$ and $D_{H,i}$ such that $H_{i}=L_{H,i}D_{H,i}L_{H,i}^{T}$ for $1\leq i\leq N$. The $LDL^{T}$ factorization of $H$ can be written in the matrix form as in (5), where $L_{H,i}$ is unit lower triangular and $D_{H,i}$ is diagonal with the $1\times 1$ pivots on the diagonal, $1\leq i\leq N$. ∎

We proceed with the computation of $B^{T}L_{H}^{-T}D_{H}^{-1}$ in the factor $L$ in (4).

Lemma 2**.**

Let $K=LDL^{T}$ be the $LDL^{T}$ factorization of the saddle-point matrix (2) with blocks given in (3). Then

[TABLE]

where

[TABLE]

where $[v^{T}\ 0]^{T}\in\mathbb{R}^{k+1}$, $[w^{T}\ \beta]^{T}\in\mathbb{R}^{k+1}$, $[I\ 0]\in\mathbb{R}^{k\times(k+1)}$ and $[M_{i}\ v_{i}]\in\mathbb{R}^{k\times(k+1)}$ for $1\leq i\leq N-1$.

Proof.

The result follows from the direct computation of the matrix product $BL_{H}^{-T}D_{H}^{-1}$. The matrix $B$ is given in (3) and factors $L_{H}$ and $D_{H}$ in (5). Since $L_{H}$ is block diagonal, then its inverse $L_{H}^{-1}$ is also block diagonal with blocks of the same size. ∎

To finish the description of the factor $L$ in (4) one needs to compute the Schur complement $S=K/H$ and factorize it. Lemma 3 shows the block $3$-diagonal structure of the Schur complement $S$.

Lemma 3**.**

The Schur complement $S=-C-B^{T}H^{-1}B$ has the form

[TABLE]

where

[TABLE]

with $[v^{T}\ 0]^{T}\in\mathbb{R}^{k+1}$, $[w^{T}\ \beta]^{T}\in\mathbb{R}^{k+1}$, $[I\ 0]\in\mathbb{R}^{k\times(k+1)}$, $\gamma_{1}\geq 0$ and $\gamma_{2}\geq 0$ from the matrix $C$, and $[M_{i}^{T}\ v_{i}]\in\mathbb{R}^{k\times(k+1)}$ for $1\leq i\leq N-1$.

Proof.

The result follows from the matrix product $BH^{-1}B^{T}$, where $H^{-1}$ is a block diagonal matrix and $B$ is given in (3). The matrix $C=diag(\gamma_{1},0,\ldots,0,\gamma_{2})$, therefore, it only affects the values $\alpha_{1}$ and $\alpha_{2}$ in (7). ∎

The Schur complement $S$ is tri-block diagonal. We shall illustrate the process of $LDL^{T}$ factorization for $N=3$. Then

[TABLE]

where $\hat{V}_{1}=V_{1}-w_{1}w_{1}^{T}/\alpha_{1}$. Once we form the $LDL^{T}$ factorization of $\hat{V}_{1}=\hat{L_{1}}\hat{D_{1}}\hat{L}_{1}^{T}$, then

[TABLE]

where $\hat{V}_{2}=V_{2}-W_{1}\hat{L}_{1}^{-T}\hat{D_{1}}^{-1}\hat{L}_{1}^{-1}W_{1}^{T}=V_{2}-W_{1}\hat{V}_{1}^{-1}W_{1}^{T}$. Once more we form the $LDL^{T}$ factorization of $\hat{V}_{2}=\hat{L}_{2}\hat{D}_{2}\hat{L}_{2}^{T}$, then

[TABLE]

where $\hat{\alpha}_{2}=\alpha_{2}-w_{2}^{T}\hat{L}_{2}^{-T}\hat{D}_{2}^{-1}\hat{L}_{2}^{-1}w_{2}=\alpha_{2}-w_{2}^{T}\hat{V}_{2}^{-1}w_{2}$. Finally, we put $\hat{\alpha}_{2}$ into the diagonal matrix $D_{S}$ and set the last elemental on the diagonal in $L_{S}$ to one.

Lemma 4**.**

Let $S=L_{S}D_{S}L_{S}^{T}$ be the $LDL^{T}$ factorization of the Schur complement $S=-C-B^{T}H^{-1}B$, then

[TABLE]

where $\hat{L}_{i}\in\mathbb{R}^{k\times k}$, $1\leq i\leq N-1$, are unit lower-triangular, $l_{1}=w_{1}/\alpha_{1}\in\mathbb{R}^{k}$, $l_{2}=\hat{D}_{N-1}^{-T}\hat{L}_{N-1}^{-1}w_{2}\in\mathbb{R}^{k}$ and the diagonal matrix $D_{S}$ is such that

[TABLE]

where $d_{1}\in\mathbb{R}$, $d_{2}\in\mathbb{R}$ and $D_{S,i}\in\mathbb{R}^{k\times k}$ for $1\leq i\leq N-1$. Here $d_{1}=\alpha_{1}$ from (7) and $d_{N}=\alpha_{2}-w_{2}^{T}\hat{L}_{N-1}^{-T}D_{S,N-1}^{-1}\hat{L}_{N-1}^{-1}w_{2}$. The scalar $\alpha_{2}$ and the vector $w_{2}$ come from (7).

Lemmas 1, 2 and 4 describe the block structure of the factor $L$. Note that both matrices $B^{T}L_{H}^{-T}D_{H}^{-1}$ and $L_{S}$ are banded and their width is independent of $N$. In Fig. 1, it is illustrated what the structures of nonzero entries of the saddle-point matrix $K$ and its factor $L$ are.

We conclude this section by the observation that for the block $L_{S}$ the option for parallel computation of blocks $V_{i}$ and $W_{i}$ at the same time in (7) is no longer available as in the case of the block $L_{H}$ in (5). However, we do not need to keep the whole Schur complement $S$ in memory to get $L_{S}$ as it is shown in the discussion preceding Lemma 4 and summarized in Alg. 3.

5 Meeting the Ill-conditioned $H$

The accuracy of the computed solution $u$ of the saddle-point system (2) with blocks from (3) by Alg. 1 depends on the condition number of $H$ [9, p. 171]. For the matrix $K$ is indefinite the condition number of $K$ may be much smaller than the condition number of $H$ [9, p. 171]. Therefore, one may try to find a permutation matrix $P$ such that the factorization $PKP^{T}=LDL^{T}$ gives better numerical results.

One can also find an ill-conditioned $(1,1)$ block in a class of interior point methods, electrical networks modelling and in the finite elements for a heat application [24, 25, 26]. We also encountered the ill-conditioned $(1,1)$ block of $K$ in dynamical system optimization [12]. Note that there are applications where the $(1,1)$ block is singular such as time-harmonic Maxwell equations [10, 11], and linear dynamical systems in the paper [12].

One such approach to computation of $PKP^{T}=LDL^{T}$ is Bunch-Parlett [7]. However, this leads to a dense factor $L$ as illustrated in Fig. 2. In addition, finding elements for a pivoting strategy is very expensive and one needs to search a matrix for its maximal off-diagonal element. For a nonsingular symmetric matrix $A\in\mathbb{R}^{n\times n}$, the pivoting strategy requires between $n^{3}/12$ and $n^{3}/6$ comparisons [7, Sec. 6.3]. On the other hand, the Bunch-Parlett factorization is nearly as stable as the LU factorization with complete pivoting [7].

BFGS approximations of the matrix $H$ (3) approach a singular matrix as the iteration process progresses [12]. Therefore, the accuracy of solutions computed by Alg. 1 deteriorates [9, p. 171]. This happened to us in several cases when we used benchmark problems from [12].

However, we observed that the application of Bunch-Parlett instead of Alg. 1 did not fail and delivered desired results. This leads us to the formulation of a hybrid method Alg. 4 that switches at some point from the straightforward $LDL^{T}$ factorization without pivoting to Bunch-Parlett. Our idea is to use Alg. 1 as long as possible until the condition number of $D_{H}$ gets large. When this behaviour is detected the method switches for Bunch-Parlett to finish. From numerical testing we found a suitable rule for switching to be that condition number of the diagonal matrix $D_{H}$ is greater than $1/\sqrt[3]{\varepsilon}$, where the machine precision $\varepsilon\approx 10^{-16}$.

Our goal in the next section is to reduce the amount of work spent on searching for pivots in Bunch-Parlett.

6 Updating the Matrix $P$

We do not need to compute a new permutation for the Bunch-Parlett method in every iteration and may try to use and update the one from the previous iteration in Alg. 4. Let $P_{i}K_{i}P_{i}^{T}=L_{i}D_{i}L_{i}^{T}$ be the Bunch-Parlett factorization of the saddle-point matrix $K_{i}$ from the $i$-th iteration in Alg. 4. In the very next iteration we try to compute $P_{i}K_{i+1}P_{i}^{T}=L_{i+1}D_{i+1}L_{i+1}^{T}$, where $D_{i+1}$ has the same pattern of $1\times 1$ and $2\times 2$ pivots as the matrix $D_{i}$. Such a factorization of $K_{i+1}$ may not exist and we need to monitor the pivots and update the permutation matrix $P_{i}$ if necessary.

In the paper [21] there is an algorithm for updating a matrix factorization. However, the method works for matrices of the form $K_{i+1}=K_{i}+\sigma vv^{T}$, where $\sigma\in\mathbb{R}$, $v\in\mathbb{R}^{n+m}$ and $\sigma vv^{T}$ is a rank-one matrix. This is not the case in the problem we try to solve.

We employ the following monitoring strategy as we follow the pattern of pivots of $D_{i}$ and factorize the permuted saddle-point matrix $P_{i}K_{i+1}P_{i}^{T}$. If the $1\times 1$ pivot $\beta$ satisfies $|\beta|>\varepsilon_{1}$, where $\varepsilon_{1}>0$, we use it and leave $P_{i}$ unchanged. If not, we apply the Bunch-Parlett method on the reduced matrix and update the permutation. In the case of the $2\times 2$ pivot $\beta=\begin{bmatrix}a&b\\ b&c\end{bmatrix}$ we accept it if $|ac-b^{2}|>\varepsilon_{1}$ and $\|\beta\|<\varepsilon_{2}$. If these conditions do not hold we apply the Bunch-Parlett method on the reduced matrix and update the permutation $P_{i}$.

The $1\times 1$ pivot is useful if and only if $|\beta|$ is large relative to the largest off-diagonal element in absolute value [7, Sec. 4.2]. Therefore, we only bound its modulus from below and the larger it is the better. However, it may happen that $|\beta|$ is small compared to off-diagonal elements and that causes the increase in the magnitudes of elements in the factor $L$ and the reduced matrix. That is the reason for introducing the second condition on the $2\times 2$ pivots for which we require $\|\beta\|<\varepsilon_{2}$. Since we divide by the determinant of the $2\times 2$ pivot $\beta$ we need it to be bounded away from zero [7, Sec. 4.3]. With this monitoring strategy one is not restricted to Bunch-Parlett and we also tried it for the Bunch-Kaufman method.

In the rest of this section, we compare the growth factor of the elements in the reduced matrices for the pivot monitoring strategy with Bunch-Parlett [7]. We use the same notation as in that paper. Let $A\in\mathbb{R}^{n\times n}$ be symmetric and nonsingular such that

[TABLE]

where $A_{1,1}\in\mathbb{R}^{j\times j}$, $c\in\mathbb{R}^{j\times(n-j)}$ and $A_{2,2}\in\mathbb{R}^{(n-j)\times(n-j)}$. If $A_{1,1}^{-1}$ exists, then

[TABLE]

where $I_{j}$ is the identity matrix of order $j$, and $I_{n-j}$ or order $n-j$ respectively. The elements of the matrix $M=A_{1,2}^{T}A_{1,1}^{-1}$ are called multipliers and we consider $j=1$ or $2$. We denote by $A^{(n)}=A$ and let $A^{(k)}$ be the reduced matrix of order $k$. In the end, let $\mu_{0}=\max_{i,j}\{|a_{i,j}|\,;a_{i,j}\in A\}$ and $\mu_{1}=\max_{i}\{|a_{i,i}|\,;a_{i,i}\in A\}$.

Suppose that the pivot $\beta=A_{1,1}$ is of order $1$, that is $j=1$. Under our monitoring strategy we accept $\beta$ for the pivot if $|\beta|>\varepsilon_{1}$. Then the reduced matrix is $A^{(n-1)}=A_{2,2}-A_{1,2}^{T}\beta^{-1}A_{1,2}$.

Lemma 5**.**

Let $\varepsilon_{1}\in(0,1)$. If $|\beta|>\varepsilon_{1}$, then

[TABLE]

Proof.

We follow [7, Lem. 1] and replace $\mu_{1}$ with our lower bound $\varepsilon_{1}$ on the magnitude of the pivot $\beta$. ∎

We observe in Lemma 5 that the bound is more pessimistic than the bound in [7, Lem. 1] for Bunch-Parlett. In more detail the bound $m<1.562$, see [7, Lem. 5], on multipliers may not hold under our monitoring strategy. The reason behind this is that we take the pivot $\beta$ even if $|\beta|<\mu_{0}(1+\sqrt{17})/8$ as long as $|\beta|>\varepsilon_{1}$. Therefore, during the factorization the elements in the reduced matrix may grow in magnitude rapidly.

Suppose the pivot $\beta$ is of order $2$, that is $j=2$. We accept $\beta$ if $|\det\beta|>\varepsilon_{1}$ and $\|\beta\|<\varepsilon_{2}$. Then the matrix $M\in\mathbb{R}^{(n-2)\times 2}$.

Lemma 6**.**

Let $\varepsilon_{1}\in(0,1)$. If $|\det\beta|>\varepsilon_{1}$, then

[TABLE]

Proof.

The first and the second inequality follows from [7, Lem. 2], where we use the lower bound $\varepsilon_{1}$ instead of $\det\beta$ in the denominator. The last chain of inequalities follows partially from our assumption that $|\det\beta|>\varepsilon_{1}$ and [7, Lem. 3]. ∎

Note that we accept the pivot $\beta$ even if $\mu_{1}>\mu_{0}$ as long as $|\det\beta|>\varepsilon_{1}$ and $\|\beta\|<\varepsilon_{2}$ hold. Due to this fact we cannot derive the bound $\mu_{0}^{(k)}<(2.57)^{n-k}\mu_{0}$ as in the paper [7]. In our case the bound $m$ on multipliers depends on $1/\varepsilon_{1}$, therefore, the overall bound on the growth of elements in the reduced matrix $A^{(k)}$ will contain powers of $1/\varepsilon_{1}$. Because of this very reason we place the condition $\|\beta\|<\varepsilon_{2}$ on the $2\times 2$ pivot trying to meet the undesired growth of elements in the reduced matrices.

A less expensive method than Bunch-Parlett is Bunch-Kaufman [6], since it requires only $O(n^{2})$ comparisons when searching the matrix for pivots. It is accepted as the algorithm of choice when solving symmetric indefinite linear systems. Similarly to Bunch-Parlett one can show that the growth of elements in the reduced matrices is bounded [6], however, this time there is no bound on the entries of the factor $L$ [3, Sec. 2.2]. Therefore, it gives lower accuracy and can even be unstable [3]. There are ways around this as described in [3], however, then the modified Bunch-Kaufman requires a higher number of comparisons lying somewhere between the number of comparisons of Bunch-Kaufman and Bunch-Parlett.

Our proposed heuristic also does not bound entries in the factor $L$, however, it tries to skip the search for pivots. In Section 7 we compare the pivot monitoring strategy with Bunch-Parlett and Bunch-Kaufman on a series of benchmarks.

7 Computational Experiments

In this section we apply our method 4 to two benchmark problems from the paper [12]. We test Alg. 4 with and without the monitoring strategy for updating the permutation matrix and compare the results. Both following benchmark problems are described in detail in the paper [12]. In this paper, we consider only equality constraints, hence, the block $C$ in (2) is the zero matrix. For the reader’s convenience we describe the governing differential equations of those dynamical systems.

The first benchmark problem [12, Sec. 6.2] is a linear dynamical system given by

[TABLE]

where the statespace dimension is $k\in\mathbb{N}$. It was shown in [12] that the Hessian matrix is singular and observed that the BFGS approximation of the Hessian approaches a singular matrix. This leads to a saddle-point matrix $K$ that has an ill-conditioned $(1,1)$ block.

The second benchmark problem [12, Sec. 6.3] is a nonlinear dynamical system such that

[TABLE]

where the statespace dimension is $k\in\mathbb{N}$. Similar to the benchmark above the BFGS approximations of the Hessian approach singular matrix [12].

In both benchmark problems the sets $\mathrm{Init}$ and $\mathrm{Unsafe}$ are balls of radius $1/4$. The stopping criteria on the norm of the gradient, the norm of the vector of constraints, the maximum number of iterations and the minimal step-size are the same as in the paper [12]. For our monitoring strategy in Alg. 4 we set $\varepsilon_{1}=10^{-3}$ and $\varepsilon_{2}=10^{6}$. The results for the first benchmark problem are shown in the Tab. 1 and for the second benchmark problem in the Tab. 2 respectively. In all instances we were able to find a desired solution from $N$ solution segments for which $x_{0}^{1}\in\mathrm{Init}$ and $\Phi(\sum t_{i},x_{0}^{1})\in\mathrm{Unsafe}$.

Both tables 1 and 2 have four parts. The first part consists of two columns denoted by $k$ – the statespace dimension of the dynamical system and by $N$ – the number of solution segments. The second part corresponds to Alg. 4 with no monitoring of pivots and there are three columns: $\#IT$ – the number of iterations in Alg. 4, $\#LDL^{T}$ – the number of straightforward $LDL^{T}$ factorizations, and $\#B$-$P$/$B$-$K$ – the number of Bunch-Parlett/Bunch-Kaufman factorizations. The third part shows the results of Alg. 4 with the monitoring of pivots. The meaning of columns $\#IT$ and $\#LDL^{T}$ remains the same. However, the column denoted by “$\#$upd of $P$” shows how many times the matrix $P$ was computed and updated.

In the end, the last column denoted by $R$ gives the ratio of the number of Bunch-Parlett/Bunch-Kaufman factorizations to the number of updates of $P$ using the pivot monitoring, that is

[TABLE]

One can interpret $R$ in the following way. As $R$ approaches $\#B$-$P$/$B$-$K$, then Alg. 4 reuses pivots almost all the time. Especially, when $R=\#B$-$P$/$B$-$K$, then we carry the search for pivots only once. When $R$ approaches $1$, then Alg 4 searches $K$ for pivots more frequently, and ultimately when $R=1$ then it uses standard Bunch-Parlett/Bunch-Kaufman throughout.

It may happen that the sum of numbers from $\#LDL^{T}$ and $\#B$-$P$/$B$-$K$ columns is greater than the number in the column $\#IT$ because of the restarts in the LS-SQP method [12]. Whenever there is a single value in a column, then both Bunch-Parlett and Bunch-Kaufman yield the same results. If the values in columns $\#IT$ in one row differ we do not compute $R$. The same applies when only the $LDL^{T}$ factorization with no pivoting was used.

We read the results in Tab. 1 and 2 in the following way. For example the last row of Tab. 1 is: the statespace dimension $k=40$ and the number of solution segments $N=30$ result in the optimization problem with $30(40+1)=1230$ parameters and $(30-1)40+2=1162$ equality constraints. Then the saddle-point matrix $K$ is of order $2392$; Alg. 4 with no monitoring of pivots took $59$ iterations from which the matrix $K$ was factorized $9$ times by $LDL^{T}$ with no pivoting and $50$ times by Bunch-Parlett/Bunch-Kaufman. When we used Alg. 4 with the monitoring of pivots the matrix $P$ was computed once in the $10$th iteration. From this point onwards it was updated twice, therefore, the matrix $P$ was reused $47$ times. The ratio $R$ is then $\lfloor 50:3\rfloor=16$.

We demonstrated in Tab. 1 and 2 that we can switch between a cheap factorization ($LDL^{T}$ without pivoting) and Bunch-Parlett. In addition we can minimize the cost of finding the pivots in the Bunch-Parlett method. The monitoring strategy that allows us to reuse the permutation matrices in the Bunch-Parlett method is independent of our application and may be used in other problems as well.

We can compare the results in Tab. 1 and 2 with the results in the paper [12], where the preconditioned projected conjugate gradient (PPCG) method [16, 18] was used. For the linear benchmark problem the Alg. 4 required less iterations of the LS-SQP. However, in the nonlinear case, the results are inconclusive.

All the computations were carried out in Scilab 5.5.2 [20] on a computer Intel(R) Xeon(R) CPU X5680 @ 3.33GHz with the operating system Cent OS 6.8. We used the built in ode solver ode in the default settings and the backslash operator for solving systems of linear equations in Alg. 1–3 and 4.

8 Conclusion

We proposed and tested a pivot monitoring strategy that allows us to reuse and update permutation matrices. Therefore, we reduced the cost of finding the pivots in solving a sequence of a saddle-point systems. Numerical experiments show that this successfully speeds up computation in the frame of dynamical systems optimization.

The result is a method that is less stable than the unmodified Bunch-Parlett, as shown in Section 6. However, practice has shown this is very often not a big concern. For example Ashcraft at al. [3, p. 552], observe that for sparse matrices and symmetric factorizations “very often less stable algorithms appear to perform numerically just as well as more reliable algorithms.” Another observation is that the unpivoted factorization requires less time and storage than pivoted factorizations with stability guarantees [3, p. 552]. Our experiments confirm those observations.

Another observation that we find useful is the following. When matrices $K_{i}$ have fixed structure of nonzero entries and the matrix $P$ remains unchanged, then the structure of nonzero entries in the factor $L$ remains the same. This becomes interesting for memory allocation of sparse matrices. Keeping the same matrix $P$ is a sort of data preprocessing. We can arrange matrices so that the pivots are on the diagonal during the factorization.