Inexact Newton method for minimization of convex piecewise quadratic functions

TL;DR

This paper introduces an inexact Newton method tailored for minimizing convex piecewise quadratic functions, with proven convergence and practical applications in linear algebra and computational geometry.

Contribution

It presents a novel inexact Newton algorithm for convex piecewise quadratic functions, extending its application to linear systems and polyhedral distance problems.

Findings

Method successfully applied to large sparse matrices.

Convergence analysis confirms theoretical robustness.

Effective in practical computational geometry tasks.

Abstract

An inexact Newton type method for numerical minimization of convex piecewise quadratic functions is considered and its convergence is analyzed. Earlier, a similar method was successfully applied to optimizaton problems arising in numerical grid generation. The method can be applied for computing a minimum norm nonnegative solution of underdetermined system of linear equations or for finding the distance between two convex polyhedra. The performance of the method is tested using sample data from NETLIB family of the University of Florida sparse matrix collection as well as quasirandom data.

Tables

Table 1. Table 1: Matrix properties for 11 NetLib problems

name	m	n	nz(A)	$\Vert x_{\ast }\Vert$	${\min }_i{\left(A{A}^T\right)}_{ii}$	${\max }_i{\left(A{A}^T\right)}_{ii}$
afiro	27	51	102	634.029569	1.18490000	44.9562810
addlittle	56	138	424	430.764399	1.00000000	10654.0000
agg3	516	758	4756	765883.022	1.00000001	179783.783
25fv47	821	1876	10705	3310.45652	0.00000000	88184.0358
pds$\_$02	2953	7716	16571	160697.180	1.00000000	91.0000000
cre$\_$a	3516	7248	18168	1162.32987	0.00000000	27476.8400
80bau3b	2262	12061	23264	4129.96530	1.00000000	321739.679
ken$\_$13	28362	42659	97246	25363.3224	1.00000000	170.000000
maros$\_$r7	3136	9408	144848	141313.207	3.05175947	3.37132546
cre$\_$b	9648	77137	260785	624.270129	0.00000000	27476.8400
osa$\_$14	2337	54797	317097	119582.321	18.0000000	845289.908

Table 2. Table 2: Computational results of different solvers for 11 NetLib problems

name	solver	time(sec)	${\Vert Ax-b\Vert }_{\infty }$	#NewtIter	#MVMult
afiro	GNewtEGK	0.001	8.63E–11	17	398
–”–	ssGNewton	0.06	6.39E–14	–	–
–”–	cqpMOSEK	0.31	1.13E–13	–	–
addlittle	GNewtEGK	0.003	6.45E–10	22	1050
–”–	ssGNewton	0.05	2.27E–13	–	–
–”–	cqpMOSEK	0.35	7.18E–11	–	–
agg3	GNewtEGK	0.14	3.93E–07	116	9234
–”–	ssGNewton	0.27	3.59E–08	–	–
–”–	cqpMOSEK	0.40	2.32E–10	–	–
25fv47	GNewtEGK	0.54	7.15E–10	114	32234
–”–	ssGNewton	1.51	3.43E–09	–	–
–”–	cqpMOSEK	1.36	1.91E–11	–	–
pds$\_$02	GNewtEGK	0.32	1.55E–08	75	8559
–”–	ssGNewton	2.30	1.40E–07	–	–
–”–	cqpMOSEK	0.51	8.20E–06	–	–
cre$\_$a	GNewtEGK	3.36	2.64E–09	219	85737
–”–	ssGNewton	1.25	4.13E–06	–	–
–”–	cqpMOSEK	0.61	2.15E–10	–	–
80bau3b	GNewtEGK	0.27	3.33E–09	79	6035
–”–	ssGNewton	0.95	1.18E–12	–	–
–”–	cqpMOSEK	0.80	2.90E–07	–	–
ken$\_$13	GNewtEGK	1.41	2.70E–08	55	6285
–”–	ssGNewton	9.09	4.39E–09	–	–
–”–	cqpMOSEK	2.09	1.71E–09	–	–
maros$\_$r7	GNewtEGK	0.10	1.18E–09	27	535
–”–	ssGNewton	2.86	2.54E–11	–	–
–”–	cqpMOSEK	55.20	3.27E–11	–	–
cre$\_$b	GNewtEGK	9.25	6.66E–10	75	24590
–”–	ssGNewton	13.20	1.62E–09	–	–
–”–	cqpMOSEK	2.31	1.61E–06	–	–
osa$\_$14	GNewtEGK	42.59	8.25E–08	767	104874
–”–	ssGNewton	60.10	4.10E–08	–	–
–”–	cqpMOSEK	4.40	7.82E–05	–	–

Table 3. Table 3: Comparison of JCG stopping criteria γ ( i ) ≤ ε C ​ G 2 ​ γ ( 0 ) superscript 𝛾 𝑖 superscript subscript 𝜀 𝐶 𝐺 2 superscript 𝛾 0 \gamma^{(i)}\leq\varepsilon_{CG}^{2}\gamma^{(0)} (old) and ( 23 ) (new)

criterion	${\epsilon }_{CG}$	geom.mean time	arithm.mean time	geom.mean res.
old	0.05	1.78	12.37	3.64–09
old	0.03	1.45	12.21	3.47–09
old	0.01	1.35	14.91	2.11–09
old	0.003	1.48	21.14	2.14–09
old	0.001	1.82	29.12	3.64–09
new	0.003	1.66	11.12	3.46–09
new	0.002	1.39	11.36	3.64–09
new	0.001	1.26	10.81	4.65–09
new	0.0003	1.26	11.25	4.23–09
new	0.0001	1.33	12.47	3.60–09

Table 4. Table 4: Performance of the generalized Newton method ( ε = 10 − 4 𝜀 superscript 10 4 \varepsilon=10^{-4} ) for the problem of distance between two quasirandom convex polyhedrons with n / 2 𝑛 2 n/2 faces each

n	${\Vert x_1-x_2\Vert }_2$	${\Vert A^T{x}_{\ast }-c\Vert }_{\infty }$	time(sec)	${\Vert g\left(x_{\ast }\right)\Vert }_{\infty }$	#NewtIter
8	0.001815	9.69–09		7.89–13	15
16	0.481528	8.63–05		1.27–13	3
32	0.795116	8.80–05		1.46–12	28
64	1.102286	1.32–04		5.58–13	13
128	1.446262	1.36–04		7.12–13	17
256	1.449913	9.54–05		4.37–13	11
512	1.460197	1.31–04	$0.001$	8.16–13	15
1024	1.460063	1.46–04	$0.002$	1.09–12	14
2048	1.463320	1.04–04	$0.005$	6.58–13	19
4096	1.463766	1.26–04	$0.009$	3.59–13	20
8192	1.463879	1.03–04	$0.009$	8.32–14	12
16384	1.463976	7.58–05	$0.009$	1.64–12	13
32768	1.464046	3.28–05	$0.018$	1.54–12	13

Equations

$$p_{\ast }=\arg \underset{p\in R^m}{\min }\left(\frac{1}{2}\Vert (\hat{x}+A^{\mathrm{T}}p{)}_{+}{\Vert }^2-b^{\mathrm{T}}p\right),$$

$$x_{\ast }=\arg \underset{{}_{x\ge 0}^{Ax=b}}{\min }\frac{1}{2}\Vert x-\hat{x}{\Vert }^2,$$

$$\phi (p)=\frac{1}{2}\Vert (\hat{x}+A^{\mathrm{T}}p{)}_{+}{\Vert }^2-b^{\mathrm{T}}p,$$

$$g(p)=A(\hat{x}+A^{\mathrm{T}}p{)}_{+}-b,$$

$$H(p)=A\mathrm{Diag}\left(\mathrm{sign}(\hat{x}+A^{\mathrm{T}}p{)}_{+}\right)A^{\mathrm{T}}.$$

$$\frac{1}{2}((\eta +\zeta {)}_{+}{)}^2-\frac{1}{2}({\eta }_{+}{)}^2-\zeta {\eta }_{+}={\zeta }^2{\int }_0^1\left({\int }_0^1\mathrm{sign}(\eta +st\zeta {)}_{+}ds\right)tdt.$$

$$f(\eta +\zeta )=f(\eta )+\zeta f^{\prime }(\eta )+{\zeta }^2{\int }_0^1\left({\int }_0^1{f}^{\prime \prime }(\eta +st\zeta )ds\right)tdt$$

$$\frac{1}{2}\Vert (y+z{)}_{+}{\Vert }^2-\frac{1}{2}\Vert y_{+}{\Vert }^2-z^{\mathrm{T}}{y}_{+}=\frac{1}{2}{z}^{\mathrm{T}}\mathrm{Diag}(d)z,$$

$$d={\int }_0^1\left({\int }_0^1\mathrm{sign}(y+stz{)}_{+}ds\right)2tdt.$$

$$\phi (p+q)-\phi (p)-q^{\mathrm{T}}g(p)=\frac{1}{2}{q}^{\mathrm{T}}A\mathrm{Diag}(d)A^{\mathrm{T}}q,$$

$$d={\int }_0^1\left({\int }_0^1\mathrm{sign}(\hat{x}+A^{\mathrm{T}}p+st{A}^{\mathrm{T}}q{)}_{+}ds\right)2tdt.$$

$$\mathrm{sign}(\hat{x}+A^{\mathrm{T}}p+\vartheta A^{\mathrm{T}}q{)}_{+}=\mathrm{sign}(\hat{x}+A^{\mathrm{T}}p{)}_{+},$$

$$\phi (p+q)-\phi (p)-q^{\mathrm{T}}g(p)=\frac{1}{2}{q}^{\mathrm{T}}H(p)q,$$

$$|(A^{\mathrm{T}}q{)}_j|\le |(\hat{x}+A^{\mathrm{T}}p{)}_j|\mathrm{whenever}(A^{\mathrm{T}}q{)}_j(\hat{x}+A^{\mathrm{T}}p{)}_j<0;$$

$$p_{k+1}=p_k-{\alpha }_k{M}_k^{-1}g(p_k),$$

$$M_k=H(p_k)+\delta \mathrm{Diag}(A{A}^{\mathrm{T}}).$$

$$d_k^{\mathrm{T}}{g}_k=d_k^{\mathrm{T}}{M}_k{d}_k={\vartheta }_k^2{g}_k^{\mathrm{T}}{M}_k^{-1}{g}_k$$

$$\phi (p_k-\alpha d_k)\le \phi (p_k)-\frac{\alpha }{2}{d}_k^{\mathrm{T}}g(p_k),\alpha \in \{1,\text{leavevmode}1/2,\text{leavevmode}1/4,\text{leavevmode}\dots \},$$

$$\phi (p+q)-\phi (p)-q^{\mathrm{T}}g(p)\le \frac{\gamma }{2}{q}^{\mathrm{T}}Mq$$

$$A{\rm Diag}(d)A^{\rm T}\leq AA^{\rm T}\leq\|\widehat{A}\|^{2}{\rm Diag}(AA^{\rm T})\leq\frac{\|\widehat{A}\|^{2}}{\delta}\Bigl{(}\delta{\rm Diag}(AA^{\rm T})+H(p)\Bigr{)}=\frac{\|\widehat{A}\|^{2}}{\delta}M;$$

$$\gamma =\Vert \hat{A}{\Vert }^2/\delta .$$

$$\phi (\hat{p})\le \phi (p)-\frac{{\vartheta }^2}{4\gamma }{g}^{\mathrm{T}}{M}^{-1}g.$$

$$\sum _{j=0}^{k-1}{g}_j^{\mathrm{T}}{g}_j\le \frac{4\gamma (1+\delta )\Vert A{\Vert }^2(\phi (p_0)-\phi (p_{\ast }))}{{\vartheta }_{\min }^2}.$$

$$\underset{k\to \infty }{\lim }\Vert g(p_k)\Vert =0,$$

$$\phi (p-\beta d)$$

$$\phi (\hat{p})\le \phi (p)-\frac{1}{2}{d}^{\mathrm{T}}Md.$$

$$\phi (p)-\alpha d^{\mathrm{T}}Md<\phi (p-2\alpha d)\le \phi (p)-(2\alpha -2\gamma {\alpha }^2)d^{\mathrm{T}}Md,$$

$$\phi (p-\alpha d)\le \phi (p)-\frac{\alpha }{2}{d}^{\mathrm{T}}Md,$$

$$\phi (p-\alpha d)\le \phi (p)-\frac{1}{4\gamma }{d}^{\mathrm{T}}Md.$$

$$d^{(i)}=\sum _{j=0}^{i-1}{s}^{(j)},$$

Full text

1 Introduction

The present paper is devoted to theoretical and experimental study of novel techniques for incorporation of preconditioned conjugate gradient linear solver into inexact Newton method. Earlier, similar method was successfully applied to optimizaton problems arising in numerical grid generation GK99 ; GKK04 ; Ka03 , and here we will consider its application to the numerical solution of piecewise-quadratic unconstrained optimization problems GGEN09 ; KMPN17 ; Ma02 ; Ma04 . The latter include such problems as finding the projection of a given point onto the set of nonnegative solutions of an underdetermined system of linear equations GGE18 or finding a distance between two convex polyhedra Bo89 (and both are tightly related to the standard linear programming problem). The paper is organized as follows. In Section 2, a typical problem of minimization of piecewise quadratic function is formulated. In Section 3, certain technical results are given related to objective functions under consideration. Section 4 describes an inexact Newton method adjusted to the optimization problem. In Section 5, a convergence analysis of the proposed algorithm is given with account of special stopping rule of inner linear conjugate gradient iterations. In Section 6, numerical results are presented for various model problems.

2 Optimization problem setting

Consider the piecewise-quadratic unconstrained optimization problem

[TABLE]

where the standard notation $\xi_{+}=\max(0,\xi)=(\xi+|\xi|)/2$ is used. Problem (1) can be viewed as the dual for finding projection of a vector on the set of nonnegative solutions of underdetermined linear systems of equations GGE18 ; GGEN09 :

[TABLE]

the solution of which is expressed via $p_{*}$ as $x_{*}=(\widehat{x}+A^{\rm T}p_{*})_{+}$. Therefore, we are considering piecewise quadratic function $\varphi:R^{m}\rightarrow R^{1}$ determined as

[TABLE]

which is convex and differentiable. Its gradient $g(p)={\rm grad}\leavevmode\nobreak\ p$ is given by

[TABLE]

and it has generalized Hessian HSN84

[TABLE]

The relation of $H(p)$ to $\varphi(p)$ and $g(p)$ will be explained later in Remark 1.

3 Taylor expansion of $(\cdot)_{+}^{2}$ function

The following result is a special case of Taylor expansion with the residual term in integral form.

LEMMA 1. For any real scalars $\eta$ and $\zeta$ it holds

[TABLE]

PROOF. Consider $f(\xi)=\frac{1}{2}(\xi_{+})^{2}$ and note that $f^{\prime}(\xi)=\xi_{+}$ and $f^{\prime\prime}(\xi)={\rm sign}(\xi_{+})$ (note that $f^{\prime\prime}(0)$ can formally be set equal to any finite real number, and w.l.o.g. we use $f^{\prime\prime}(0)=0$). Inserting this into the Taylor expansion

[TABLE]

readily gives the desired result.

LEMMA 2. For any real $n$-vectors $y$ and $z$ it holds

[TABLE]

where

[TABLE]

PROOF. Setting in (5) $\eta=y_{j}$, $\zeta=z_{j}$, and summing over all $j=1,\ldots,n$ obviously yields the required formula. Note that the use of scalar multiple 2 within the integral provides for the estimate $\|{\rm Diag}(d)\|\leq 1$.

LEMMA 3. Function (2) and its gradient (3) satisfy the identity

[TABLE]

where

[TABLE]

PROOF. Setting in (6) and (7) $y=\widehat{x}+A^{\rm T}p$ and $z=A^{\rm T}q$ readily yields the required result (with account of cancellation of linear terms involving $b$ in the left hand side of (8)).

REMARK 1. As is seen from (9), if the condition

[TABLE]

holds true for any $0\leq\vartheta\leq 1$, then (8) is simplified as

[TABLE]

where the generalized Hessian matrix $H(p)$ is defined in (4). This explains the key role of $H(p)$ in the organization of the Newton-type method considered below. Note that a sufficient condition for (10) to hold is

[TABLE]

that is, if certain components of the increment $q$ are relatively small, then $\varphi$ is exactly quadratic (11) in the corresponding neighborhood of $p$.

4 Inexact Newton method for dual problem

As suggests condition (10) and its consequence (11), one can try to numerically minimize $\varphi$ using Newton type method $p_{k+1}=p_{k}-d_{k}$, where $d_{k}=H(p_{k})^{-1}g(p_{k})$. Note that by (11) this will immediately give the exact minimizer $p_{*}=p_{k+1}$ if the magnitudes of $d_{k}$ components are sufficiently small to satisfy (12) taken with $p=p_{k}$ and $q=-d_{k}$. However, initially $p_{k}$ may be rather far from solution, and only gradual improvements are possible. First, a damping factor $\alpha_{k}$ must be used to guarantee monotone convergence (with respect to the decrease of $\varphi(p_{k})$ as $k$ increases). Second, $H(p_{k})$ must be replaced by some appropriate approximation $M_{k}$ in order to provide its invertibility with a reasonable bound for the inverse. Therefore, we propose the following prototype scheme

[TABLE]

where

[TABLE]

The parameters $0<\alpha_{k}\leq 1$ and $0\leq\delta\ll 1$ must be defined properly for better convergence. Furthermore, at initial stages of iteration, the most efficient strategy is to use approximate Newton directions $d_{k}\approx M_{k}^{-1}g(p_{k})$, which can be obtained using preconditioned conjugate gradient (PCG) method for the solution of Newton equation $M_{k}d_{k}=g(p_{k})$. As will be seen later, it suffices to use any vector $d_{k}$ which satisfies conditions

[TABLE]

with $0<\vartheta_{k}<1$ sufficiently separated from zero. For any preconditioning, the approximations constructed by the PCG method satisfy (14), see Section 5.3 below. With account of the Armijo type criterion

[TABLE]

where the maximum steplength $\alpha$ satisfying (15) is used, the inexact Newton algorithm can be presented as follows:

Algorithm 1.

Input:

$A\in R^{m\times n}$,  $b\in R^{m}$,  $\widehat{x}\in R^{m}$;

Initialization:

$\delta=10^{-6}$,  $\varepsilon=10^{-12}$,  $\tau=10^{-15}$

$k_{\max}=2000$,  $l_{\max}=10$;  $p_{0}=0$,

Iterations:

for $k=0,1,\ldots,k_{\max}-1$:

$x_{k}=(\widehat{x}+A^{\rm T}p_{k})_{+}$

$\varphi_{k}=\frac{1}{2}\|x_{k}\|^{2}-b^{\rm T}p_{k}$

$g_{k}=Ax_{k}-b$

if $(\|g_{k}\|\leq\varepsilon\|b\|)$ return $\{{x_{k},p_{k},g_{k}\}}$

find $d_{k}\in R^{m}$  such that

$d_{k}^{\rm T}g_{k}=d_{k}^{\rm T}M_{k}d_{k}=\vartheta_{k}^{2}g_{k}M_{k}^{-1}g_{k}$,

where $M_{k}=A\leavevmode\nobreak\ {\rm Diag}({\rm sign}(x_{k}))\leavevmode\nobreak\ A^{\rm T}+\delta{\rm Diag}(AA^{\rm T})$

$\alpha^{(0)}=1$

$p_{k}^{(0)}=p_{k}-d_{k}$

for $l=0,1,\ldots,l_{\max}-1$:

$x_{k}^{(l)}=(\widehat{x}+A^{\rm T}p_{k}^{(l)})_{+}$

$\varphi_{k}^{(l)}=\frac{1}{2}\|x_{k}^{(l)}\|^{2}-b^{\rm T}p_{k}^{(l)}$

$\zeta_{k}^{(l)}=\left(\frac{1}{2}\alpha^{(l)}d_{k}^{\rm T}g_{k}+\varphi_{k}^{(l)}\right)-\varphi_{k}$

if ($\zeta_{k}^{(l)}>\tau|\varphi_{k}|$) then

$\alpha^{(l+1)}=\alpha^{(l)}/2$

$p_{k}^{(l+1)}=p_{k}-\alpha^{(l+1)}d_{k}$

else

$p_{k}^{(l+1)}=p_{k}^{(l)}$

go to NEXT

end if

end for

NEXT: $p_{k+1}=p_{k}^{(l+1)}$

end for

Next we explore the convergence properties of this algorithm.

5 Convergence analysis of inexact Newton method

It appears that Algorithm 1 exactly conforms with the convergence analysis presented in Ka03 (see also AK01 ). For the completeness of presentation and compatibility of notations, we reproduce here the main results of Ka03 .

5.1 Estimating convergence of inexact Newton method

The main assumptions we need for the function $\varphi(p)$ under consideration are that it is bounded from below, have gradient $g(p)\in R^{m}$, and satisfies

[TABLE]

for the symmetric positive definite $m\times m$ matrix $M=M(p)$ defined above in (13) and some constant $\gamma\geq 1$. Note that the exact knowledge of $\gamma$ is not necessary for actual calculations. The existence of such $\gamma$ follows from (13) and Lemma 3. Indeed, denoting $D=({\rm Diag}(AA^{\rm T}))^{1/2}$ and $\widehat{A}=D^{-1}A$, for the right hand side of (8) one has, with account of $\|{\rm Diag}(d)\|\leq 1$ and $H(p)\geq 0$,

[TABLE]

therefore, (16) holds with

[TABLE]

The latter formula explains our choice of $M$ which is more appropriate in cases of large variations in norms of rows in $A$ (see the examples and discussion in Section 6).

Next we will estimate the reduction in the value of $\varphi$ attained by the descent along the direction $(-d)$ satisfying (14). One can show the following estimate for the decrease of objective function value at each iteration (here, simplified notations $p=p_{k}$, $\hat{p}=p_{k+1}$ etc. are used) $\hat{p}=p-\alpha d$ with $\alpha=2^{-l}$, where $l=0,1,\ldots$, as evaluated according to (15):

[TABLE]

In particular, if the values of $\vartheta^{2}$ are separated from zero by a positive constant $\vartheta_{\min}^{2}$ (lower estimate for $\vartheta$ follows from Section 5.3 and an upper bound for $\kappa={\rm cond}(CM)$), then, with account for $M\leq(1+\delta)\|A\|^{2}I$ and the boundedness of $\varphi$ from below, it follows

[TABLE]

Noting that the right hand side of the latter estimate does not depend on $k$, it finally follows that

[TABLE]

where $k$ is the number of the outer (nonlinear) iteration.

Estimate (18) can be verified as follows (note that quite similar analysis can be found in Ma02 ). Using $q=-\beta d$, where $0<\beta<2/\gamma$, one can obtain from (16) and (14) the following estimate for the decrease of $\varphi$ along the direction $(-d)$:

[TABLE]

The following two cases are possible.

Case 1. If the condition (15) is satisfied at once for $\alpha=1$, this means that (recall that the left equality of (14) holds)

[TABLE]

Case 2. Otherwise, if at least one bisection of steplength was performed (and the actual steplength is $\alpha$), then, using (19) with $\beta=2\alpha$, it follows

[TABLE]

which readily yields $\alpha>1/(2\gamma)$. Since we also have

[TABLE]

it follows

[TABLE]

Joining these two cases, taking into account that $\gamma\geq 1$, and using the second equality in (14) one obtains the required estimate (18).

It remains to notice that as soon as the norms of $g$ attain sufficiently small values, the resulting directions $d$ will also have small norms. Therefore, the case considered in Remark 1 will take place, and finally the convergence of the Newton method will be much faster than at its initial stage.

5.2 Linear CG approximation of Newton directions

Next we relate the convergence of inner linear Preconditioned Conjugate Gradient (PCG) iterations to the efficiency of Inexact Newton nonlinear solver. Similar issues were considered in AK01 ; GKK04 ; KA94 ; Ka03 .

An approximation $d^{(i)}$ to the solution of the problem $Md=g$ generated on the $i$th PCG iteration by the recurrence $d^{(i+1)}=d^{(i)}+s^{(i)}$ (see Algorithm 2 below) can be written as follows (for our purposes, we always set the initial guess for the solution $d^{(0)}$ to zero):

[TABLE]

where the PCG direction vectors are pairwise $M$-orthogonal: $(s^{(j)})^{{\mbox{\rm\tiny T}}}Ms^{(l)}=0$,  $j\neq l$. Let also denote the $M$-norms of PCG directions as $\eta^{(j)}=(s^{(j)})^{{\mbox{\rm\tiny T}}}Ms^{(j)}$, $j=0,1,\ldots,i-1$. Therefore, from (21), one can determine

[TABLE]

and estimate (20) takes the form

[TABLE]

where $k$ is the Newton iteration number. Summing up the latter inequalities for $0\leq k\leq m-1$, we get

[TABLE]

On the other hand, the cost measure related to the total time needed to perform $m$ inexact Newton iterations with $i_{k}$ PCG iterations at each Newton step, can be estimated as proportional to

[TABLE]

Here $\epsilon_{\rm CG}$ is a small parameter reflecting the ratio of one linear PCG iteration cost to the cost of one Newton iteration (in particular, including construction of preconditioning and several $\varphi$ evaluations needed for backtracking) plus possible efficiency loss due to early PCG termination. Thus, introducing the function $\psi(i)=(\epsilon_{\rm CG}^{-1}+i)/\zeta^{(i)}$, (here, we omit the index $k$) one obtains a reasonable criterion to stop PCG iterations in the form $\psi(i)>\psi(i-1)$. Here, the use of smaller values $\varepsilon_{\rm CG}$ generally corresponds to the increase of the resulting iteration number bound. Rewriting the latter condition, one obtains the final form of the PCG stopping rule:

[TABLE]

Note that by this rule, the PCG iteration number is always no less than 2.

Finally, we explicitly present the resulting formulae for the PCG algorithm incorporating the new stopping rule. Following GGEN09 , we use the Jacobi preconditioning

[TABLE]

Moreover, the reformulation KM11 of the CG algorithm HS52 ; Ax76 is used. This may give a more efficient parallel implementation, see, e.g., GGEN09 .

Following KM11 , recall that at each PCG iteration the $M^{-1}$-norm of the $(i+1)$-th residual $r^{(i+1)}=g-Md^{(i+1)}$ attains its minimum over the corresponding Krylov subspace. Using the standard PCG recurrences (see Section 5.3 below) one can find $d^{(i+1)}=d^{(i)}+Cr^{(i)}\alpha^{(i)}+s^{(i-1)}\alpha^{(i)}\beta^{(i-1)}$. Therefore, the optimum increment $s^{(i)}$ in the recurrence $d^{(i+1)}=d^{(i)}+s^{(i)}$, where $s^{(i)}=V^{(i)}h^{(i)}$ and $V^{(i)}=[Cr^{(i)}\leavevmode\nobreak\ |\leavevmode\nobreak\ s^{(i-1)}]$, can be determined via the solution of the following 2-dimensional linear least squares problem:

[TABLE]

By redefining $r^{(i)}:=-r^{(i)}$ and introducing vectors $t^{(i)}=Ms^{(i)}$, the required PCG reformulation follows:

Algorithm 2.

$r^{(0)}=-g,$ $d^{(0)}=s^{(-1)}=t^{(-1)}=0,$  $\zeta^{(-1)}=0;$

$i=0,1,\ldots,it_{\max}:$

$w^{(i)}=Cr^{(i)},$

$z^{(i)}=Mw^{(i)},$

$\gamma^{(i)}=(r^{(i)})^{{\mbox{\rm\tiny T}}}w^{(i)},$  $\xi^{(i)}=(w^{(i)})^{{\mbox{\rm\tiny T}}}z^{(i)},$  $\eta^{(i-1)}=(s^{(i-1)})^{{\mbox{\rm\tiny T}}}t^{(i-1)},$

$\zeta^{(i)}=\zeta^{(i-1)}+\eta^{(i-1)},$

if  $((\varepsilon_{\rm CG}^{-1}+i)\eta^{(i-1)}\leq\zeta^{(i)})$  or  $(\gamma^{(i)}\leq\varepsilon_{\rm CG}^{2}\gamma^{(0)})$  return $\{d^{(i)}\}$;

if   ($k=0$)   then

$\alpha^{(i)}=-\gamma^{(i)}/\xi^{(i)},\quad\beta^{(i)}=0;$

else

$\delta^{(i)}=\gamma^{(i)}/(\xi^{(i)}\eta^{(i-1)}-(\gamma^{(i)})^{2})$,

$\alpha^{(i)}=-\eta^{(i-1)}\delta^{(i)}$,    $\beta^{(i)}=\gamma^{(i)}\delta^{(i)}$;

end if

$t^{(i)}=z^{(i)}\alpha^{(i)}+t^{(i-1)}\beta^{(i)},$    $r^{(i+1)}=r^{(i)}+t^{(i)},$

$s^{(i)}=w^{(i)}\alpha^{(i)}+s^{(i-1)}\beta^{(i)},$   $d^{(i+1)}=d^{(i)}+s^{(i)}.$

For maximum reliability, the new stopping rule (23) is used along with the standard one; however, in almost all cases the new rule provides for an earlier CG termination.

Despite of somewhat larger workspace and number of vector operations compared to the standard algorithm, the above version of CG algorithm enables more efficient parallel implementation of scalar product operations. At each iteration of the above presented algorithm, it suffices to use one MPI$\_$AllReduce($*$,$*$,3,…) operation instead of two MPI$\_$AllReduce($*$,$*$,1,…) operation in the standard PCG recurrences. This is especially important when many MPI processes are used and the start-up time for MPI$\_$AllReduce operations is relatively large. For another equivalent PCG reformulations allowing to properly reorder the scalar product operations, see DE03 and references cites therein.

5.3 Convergence properties of PCG iterations

Let us recall some basic properties of the PCG algorithm, see, e.g. Ax76 . The standard PCG algorithm (algebraically equivalent to Algorithm 2) for the solution of the problem $Md=g$ can be written as follows (the initial guess for the solution $d_{0}$ is set to zero):

[TABLE]

The scaling property (14) (omitting the upper and lower indices at $d$, it reads $d^{\rm T}g=d^{\rm T}Md$) can be proved as follows. Let $d=d^{(i)}$ be obtained after $i$ iterations of the PCG method applied to $Md=g$ with zero initial guess $d^{(0)}=0$. Therefore, $d\in K_{i}={\rm span}\{Cg,CMCg,\dots,(CM)^{i-1}Cg\}$, and, by the PCG optimality property, it holds

[TABLE]

Since $\alpha d\in K_{i}$ for any scalar $\alpha$, one gets

[TABLE]

Setting here $\alpha=d^{T}g/d^{T}Md$, one can easily transform this inequality as $0\geq(-d^{T}g+d^{T}Md)^{2}$, which readily yields (14). Furthermore, by the well known estimate of the PCG iteration error Ax76 using Chebyshev polynomials, one gets

[TABLE]

where

[TABLE]

By the scaling condition, this gives

[TABLE]

Hence, $0<\theta<1$ and $\theta^{2}\rightarrow 1$ as the PCG iteration number $i$ grows.

6 Numerical test results

Below we consider two families of test problems which can be solved via minimization of piecewise quadratic problems. The first one was described above in Section 2 (see also GGE18 ), while the second coincides with the problem setting for the evaluation of distance between two convex polyhedra used in Bo89 . The latter problem is of key importance e.g., in robotics and computer animation.

6.1 Test results for 11 NETLIB problems

Matrix data from the following 11 linear programming problems (this is the same selection from NETLIB collection as considered in KMPN17 ), were used to form test problems (1). Note that further we only consider the case $\widehat{x}=0$. Recall also the notation $x_{*}=(\widehat{x}+A^{T}p_{*})_{+}$. The problems in Table 1 below are ordered by the number of nonzero elements ${\rm nz}(A)$ in $A\in R^{m\times n}$.

It is readily seen that 3 out of 11 matrices have null rows, and more than half of them have rather large variance of row norms. This explains the proposed Hessian regularization (13) instead of the earlier construction GGE18 ; KMPN17 $M_{k}=H(p_{k})+\delta I_{m}$. The latter is a proper choice only for matrices with rows of nearly equal length, such as  maros$\_$r7  example or various matrices with uniformly distributed quasirandom entries, as used for testing in GGEN09 ; KMPN17 . In particular, estimate (17) with $D=I$ would take the form $\gamma=\|A\|^{2}/\delta$, so the resulting method appears to be rather sensitive to the choice of $\delta$.

In Table 2, the results presented in KMPN17 are reproduced along with similar data obtained with our version of Generalized Newton method. It must be stressed that we used the fixed set of tunung parameters

[TABLE]

for all problems. Note that In KMPN17 the parameter choice for the Armijo procedure was not specified.

In KMPN17 , the calculations were performed on 5GHz AMD 64 Athlon X2 Dual Core. In our experiments, one core of 3.40 GHz x8 Intel (R) Core (TM) i7-3770 CPU was used, which is likely somewhat slower.

Note that the algorithm of KMPN17 is based on direct evaluation of $M_{k}$ and its sparse Cholesky factorization, while our implementation, as was proposed in GGEN09 , uses the Jacobi preconditioned Conjugate Gradient iterations for approximate evaluation of Newton directions. Thus, the efficiency of our implementation critically depends on the CG iteration convergence, which is sometimes slow. On the other hand, since the main computational kernels of the algorithm are presented by matrix-vector multiplications of the type  $x=Ap$  or  $q=A^{T}y$, its parallel implementation can be sufficiently efficient.

In Table 2, the abbreviation cqpMOSEK refers to MOSEK Optimization Software package for convex quadratic problems, see KMPN17 . The abbreviation ssGNewton denotes the method implemented and tested in KMPN17 , while GNewtEGK stands for the method proposed in the present paper.

Despite the use of slower computer, our GNewtEGK demonstrates considerably faster performance in 8 cases of 11. Otherwise, one can observe that smaller computational time of cqpMOSEK goes along with much worse residual norm, see the results for problems  cre$\_$b  and  osa$\_$14 .

Thus, in most cases the presented implementation of Generalized Newton method takes not too large number of Newton iterations using approximate Newton directions generated by CG iterations with diagonal preconditioning (24) and special stopping rule (23).

A direct comparison of efficiency for the standard CG iterations stopping rule $\gamma^{(i)}\leq\varepsilon_{\rm CG}^{2}\gamma^{(0)}$ (see Algorithm 2 for the notations) and the new one (23) is given in Table 3, where the timing (in seconds) and precision results averaged over the same 11 problems are given. One can see that nearly the same average residual norm $\|Ax-b\|_{\infty}$ can be obtained considerably faster and with less critical dependence on $\varepsilon_{\rm CG}$ when using the new PCG iteration stopping rule.

6.2 Evaluating the distance between convex polyhedra

Let the two convex polyhedra ${\cal X}_{1}$ and ${\cal X}_{2}$ be described by the following two systems of linear inequalities:

[TABLE]

where $A_{1}\in R^{s\times n_{1}}$,  $A_{2}\in R^{s\times n_{2}}$, and the vectors $x_{1},x_{2},b_{1},b_{2}$ are of compatible dimensions. The original problem of evaluating the distance between ${\cal X}_{1}$ and ${\cal X}_{2}$ is (cf. Bo89 , where it was solved by the projected gradient method)

[TABLE]

We will use the following regularized/penalized approximate reformulation of the problem in terms of unconstrained convex piecewise quadratic minimization. Introducing the matrices

[TABLE]

and the vector $b=[b_{1}^{T},b_{2}^{T}]^{T}\in R^{n_{1}+n_{2}}$, we consider the problem

[TABLE]

where the regularization/penalty parameter $\varepsilon$ is a sufficiently small positive number (we have used $\varepsilon=10^{-4}$). The latter problem can readily be solved by adjusting the above described Algorithm 1 using $\delta=0$ and the following explicit expressions for the gradient and the generalized Hessian:

[TABLE]

where $D(x)={\rm Diag}({\rm sign}(A^{T}x-b)_{+})$. When solving practical problems of evaluating the distance between two 3D convex polyhedra determined by their faces (so that $s=3$), the inexact Newton iterations are performed in $R^{6}$, and the cost of each iteration is proportional to the total number $n=n_{1}+n_{2}$ of the faces determining the two polyhedra. In this case, the explicit evaluation of $H(x)$ and the use of its Cholesky factorization is more preferable than the use of the CG method.

Test polyhedrons with $n/2$ faces each were centered at the points $e=[1,1,1]^{T}$ or $-e$ and defined as $A_{1}^{T}(x-e)\leq b_{1}$ and $A_{2}^{T}(x+e)\leq b_{2}$ with $b_{1}=b_{2}=[1,\ldots,1]\in R^{n/2}$, respectively. The columns of matrices $A_{1}$ and $A_{2}$ were determined by $n/2$ quasirandom unit $3$-vectors generated with the use of logistic sequence (see, e.g.YBZS10 and references cited therein) $\xi_{0}=0.4$,  $\xi_{k}=1-2\xi_{k-1}^{2}$. We used $A_{1}(i,j)=\xi_{20(i-1+3(j-1))}$ and similar for $A_{2}$; then the columns of these matrices were normalized to the unit length. The corresponding performance results seem quite satisfactory, see Table 4. Note that the lower bound $\|x_{1}-x_{2}\|\geq 2\sqrt{3}-2\approx 1.464101$ always holds for the distance (since the two balls $\|x_{1}-e\|=1$ and $\|x_{2}+e\|=1$ are inscribed in the corresponding polyhedrons).

Acknowledgements.

This work was supported by the Russian Foundation for Basic Research grant No. 17-07-00510 and by Program No. 26 of the Presidium of the Russian Academy of Sciences. The authors are grateful to the anonymous referees and to Prof. V.Garanzha for many useful comments which greatly improved the exposition of the paper.