Piecewise linear secant approximation via Algorithmic Piecewise Differentiation

TL;DR

This paper introduces a method for locally approximating piecewise differentiable functions with a bilinear order error using only two sample points, and applies it to develop a generalized Newton's method with proven convergence properties.

Contribution

It presents a novel approach for piecewise linear approximation of functions using minimal samples and derives a stable implementation for the method.

Findings

Discrepancy between function and model is bilinear order O(‖x−x̂‖‖x−x̌‖)

Generalized Newton's method based on piecewise linearization converges under certain conditions

Formulas for numerically stable implementation of the approximation method

Abstract

It is shown how piecewise differentiable functions $F: \mathbb R^n \mapsto \mathbb R^m $ that are defined by evaluation programs can be approximated locally by a piecewise linear model based on a pair of sample points $\check x$ and $\hat x$. We show that the discrepancy between function and model at any point $x$ is of the bilinear order $O(\|x-\check x\| \|x-\hat x\|)$. This is a little surprising since $x \in \mathbb R^n$ may vary over the whole Euclidean space, and we utilize only two function samples $\check F=F(\check x)$ and $\hat F=F(\hat x)$, as well as the intermediates computed during their evaluation. As an application of the piecewise linearization procedure we devise a generalized Newton's method based on successive piecewise linearization and prove for it sufficient conditions for convergence and convergence rates equaling those of semismooth Newton. We conclude with the derivation of formulas for the numerically stable implementation of the aforedeveloped piecewise linearization methods.

Equations

$$F(x)=\mathring{F}+\Delta F(\mathring{x};x-\mathring{x})+\mathcal{O}(\Vert x-\mathring{x}{\Vert }^2).$$

$$F(x)=\mathring{F}+\Delta F(\check{x},\hat{x};x-\mathring{x})+\mathcal{O}(\Vert x-\check{x}\Vert \Vert x-\hat{x}\Vert ).$$

$$\Phi :=\tilde{\Phi }\cup \{\mathbf{abs}()\}.$$

$$\max (u,w)=(u+w+\mathbf{abs}(u-w))/2,\min (u,w)=(u+w-\mathbf{abs}(u-w))/2.$$

$$\Delta F(\mathring{x};\alpha \Delta x)=\alpha \Delta F(\mathring{x};\Delta x)\text{mbox}for0<\alpha \Vert \Delta x\Vert <\rho (\mathring{x}).$$

$$v_i=v_j\circ v_k\text{mbox}or{v}_i={\phi }_i(v_j)\text{mbox}fori=1\dots l.$$

$${\phi }_i\in \Phi \equiv \{\text{mbox}rec,sqrt,\sin ,\cos ,\exp ,\log ,\dots ,\mathbf{abs},\dots \}$$

$$\phi (u)\equiv |u|>0?u^p\cdot \sin (1/u):0\text{mbox}forp\ge 3.$$

$$(x,z,y)=\left(v_{1-n},\dots ,v_{-1},v_0,\dots ,v_{l-m},v_{l-m+1},\dots ,v_l\right)\in {\mathbb{R}}^{n+l}$$

$$v_i(\mathring{x}+\Delta x)-{\mathring{v}}_i\approx \Delta v_i\equiv \Delta v_i(\Delta x).$$

$$\Delta v_i=$$

$$\Delta y=\Delta F(\mathring{x};\Delta x)\equiv F^{\prime }(\mathring{x})\Delta x,$$

$$\Delta v_i=\mathbf{abs}({\mathring{v}}_j+\Delta v_j)-{\mathring{v}}_i\text{mbox}when{v}_i=\mathbf{abs}(v_j).$$

$$\Delta y=\Delta y(\Delta x)=\Delta F(\mathring{x};\Delta x):{\mathbb{R}}^n\to {\mathbb{R}}^m.$$

$${\mathring{v}}_i\equiv ({\check{v}}_i+{\hat{v}}_i)/2\text{mbox}of{\check{v}}_i\equiv v_i(\check{x})\text{mbox}and{\hat{v}}_i\equiv v_i(\hat{x}).$$

$$c_{ij}\equiv\left\{\begin{array}[]{ll}(\check{v}_{i}-\hat{v}_{i})/(\check{v}_{j}-\hat{v}_{j})&\;\mbox{if}\;\check{v}_{j}\neq\hat{v}_{j}\\ \varphi_{i}^{\prime}(\mathring{v}_{j})&\;\mbox{otherwise}\end{array}\,.\right.$$

$${\mathring{v}}_i\equiv \frac{1}{2}({\check{v}}_i+{\hat{v}}_i)=\frac{1}{2}[\mathbf{abs}({\check{v}}_j)+\mathbf{abs}({\hat{v}}_j)].$$

$$\Delta y=\Delta y(\Delta x)=\Delta F(\check{x},\hat{x};\Delta x):{\mathbb{R}}^n\to {\mathbb{R}}^m.$$

$$u\cdot w=\frac{1}{4}\left[(u+w{)}^2-(u-w{)}^2\right].$$

$${◊}_{\mathring{x}}F(x)\equiv F(\mathring{x})+\Delta F(\mathring{x};x-\mathring{x})$$

$${◊}_{\check{x}}^{\hat{x}}F(x)\equiv \frac{1}{2}(F(\check{x})+F(\hat{x}))+\Delta F(\check{x},\hat{x};x-\mathring{x}).$$

$${◊}_{\check{x}}^{\hat{x}}v={◊}_{\check{x}}^{\hat{x}}{x}^2$$

$$=\frac{1}{2}[{\check{x}}^2+{\hat{x}}^2]+2\mathring{x}(x-\mathring{x})=\mathring{v}+2\mathring{x}(x-\mathring{x}).$$

$${◊}_{(\check{u},\check{w})}^{(\hat{u},\hat{w})}(u\cdot w)$$

$$=\frac{1}{2}(\hat{u}\hat{w}+\check{u}\check{w})+\mathring{w}(u-\mathring{u})+\mathring{u}(w-\mathring{w}).$$

$$4\cdot {◊}_{(\check{u},\check{w})}^{(\hat{u},\hat{w})}(u\cdot w)$$

$$=\frac{1}{2}\left([\hat{u}+\hat{w}{]}^2+[\check{u}+\check{w}{]}^2\right)+2(\mathring{u}+\mathring{w})\left(u+w-[\mathring{u}+\mathring{w}]\right)$$

$$-\frac{1}{2}\left([\hat{u}-\hat{w}{]}^2+[\check{u}-\check{w}{]}^2\right)+2(\mathring{u}-\mathring{w})\left(u-w-[\mathring{u}-\mathring{w}]\right)$$

$$=\frac{1}{2}\left([\hat{u}+\hat{w}{]}^2+[\check{u}+\check{w}{]}^2\right)+2(\mathring{u}+\mathring{w})\left([u-\mathring{u}]+[w-\mathring{w}]\right)$$

$$-\frac{1}{2}\left([\hat{u}-\hat{w}{]}^2+[\check{u}-\check{w}{]}^2\right)+2(\mathring{u}-\mathring{w})\left([u-\mathring{u}]-[w-\mathring{w}]\right)$$

$$=4\cdot \left[\frac{1}{2}\left(\check{u}\check{w}+\hat{u}\hat{w}\right)+\mathring{w}[u-\mathring{u}]+\mathring{u}[w-\mathring{w}]\right].$$

$$\Vert F(x)-F(\tilde{x})\Vert$$

$$\max \left(\Vert {◊}_{\check{x}}^{\hat{x}}F(x)-{◊}_{\check{x}}^{\hat{x}}F(\tilde{x})\Vert ,\Vert {◊}_{\mathring{x}}F(x)-{◊}_{\mathring{x}}F(\tilde{x})\Vert \right)$$

Full text

Keywords Automatic differentiation, Stable piecewise linearization, Generalized Newton’s method, Lipschitz continuity, Generalized Hermite interpolation, ADOL-C

MSC2010 65D25, 65K10, 49J52

1 Introduction

In this paper we refine and extend the theory of piecewise linearizations111Our notion of linearity includes nonhomogeneous functions, where the adjective affine or perhaps polyhedral would be more precise. However, such mathematical terminology might be less appealing to computational practicioners and to the best of our knowledge there are no good nouns corresponding to linearity and linearization for the adjectives affine and polyhedral. of piecewise differentiable functions $F:\mathbb{R}^{n}\mapsto\mathbb{R}^{m}$ that was introduced in [Gri13]. Throughout, we assume that any such function is defined by an evaluation procedure consisting of a sequence of elemental functions $v_{i}=\varphi_{i}(v_{j})_{j\prec i}$, where $\varphi_{i}$ is either smooth or the absolute value function. The data dependence relation $\prec$ generates a partial ordering, which yields a directed acyclic graph. In other words, we assume that we have a straight line program without any loops or jumps in the control flow.

In the above reference a piecewise linearization $\Delta F(\mathring{x};\Delta x)$ of $F$ about the point $\mathring{x}$ is constructed such that, for arbitrary $x\in\mathbb{R}^{n}$ with $\mathring{F}=F(\mathring{x})$, it holds

[TABLE]

Since $\Delta F(\mathring{x};\Delta x)$ is constructed by replacing smooth elementals by their tangent lines at $\mathring{x}$, we will refer to it as piecewise tangent linearization. An important property of these piecewise linearizations is that they vary continuously with respect to the sample point or points at which they are developed.

As a generalization of $\Delta F(\mathring{x};\Delta x)$ we construct a piecewise secant linearization $\Delta F(\check{x},\hat{x};\Delta x)$ of $F$ at the pair $\check{x}$ and $\hat{x}$ such that for the midpoints $\mathring{x}=(\check{x}+\hat{x})/2$ and $\mathring{F}=(\check{F}+\hat{F})/2$ with $\check{F}=F(\check{x})$ and $\hat{F}=F(\hat{x})$ it holds

[TABLE]

The new $\Delta F(\check{x},\hat{x};\Delta x)$ reduces to $\Delta F(\mathring{x};\Delta x)$ when the two sample points $\check{x}$ and $\hat{x}$ coalesce at $\mathring{x}$.

Our results in this article are twofold:

• •

Firstly, we present approximation properties and Lipschitz continuity estimates. Here the two major points are that the piecewise linearization is a second order approximation to the underlying function. Moreover, we do not only prove the existence of Lipschitz constants, but provide explicit estimates. These results immediately yield additional statement on conditions for perturbation-stable surjectivity of a piecewise linear model.

• •

Secondly, we develop an application for the piecewise linearization. Namely, we introduce, for each linearization mode, a Newton’s method based on successive piecewise linearization and give sufficient conditions for convergence, as well as statements on the convergence rates.

In an appended section we moreover provide formulas for singularity free and thus numerically stable implementations of the piecewise secant linearizations.

One significant advantage of the piecewise linearization-approach is that, once the piecewise linearization is generated, we also know where its kinks are located. This means that we are liberated of the complications of event handling that the nondifferentiabilities of piecewise smooth functions usually bring about, which was one of the motivations for the development of the techniques presented in [Gri13]. In said reference it is proved that the variations of the piecewise linear model are Lipschitz continuous with respect to perturbations of the development point. The results mentioned in the first bullet point represent a significant improvement over those given in [Gri13] in so far as they are not only sharper, but also provide explicit bounds for the Lipschitz constants for the variations of the model with respect to perturbations of the base point.

Content and Structure

To provide some more background, we will, in the second section, elaborate on the setting outlined above both from the mathematical and the implementation perspective. In Section 2 we will introduce the piecewise linearization framework that is investigated thereafter. In Section 3 we derive the approximation and stability properties of said framework. In the subsequent section the generalized Newton’s methods are derived for both linearization modes. In Section 5 we show that the secant linearization can be computed in a division free, centered form such that it reduces continuously to the tangent mode when the reference points $\check{x}$ and $\hat{x}$ coalesce. The latter is followed by a numerical example in Section 6. The article is concluded by some final remarks.

Preliminary Remarks

Let $\tilde{\Phi}$ be a library of elemental functions that conforms to the condition of elemental differentiability as described in [GW08]. In the following we will consider piecewise differentiable functions that are defined by an evaluation procedure consisting of a sequence of such elemental functions $v_{i}=\varphi_{i}(v_{j})_{j\prec i}$, where the $\varphi_{i}$ are contained in a library

[TABLE]

We remark that the inclusion of the absolute value function into the library means that we can also evaluate $\min()$ and $\max()$ via the identities

[TABLE]

There is a slight implicit restriction, namely, we assume that whenever $\min$ or $\max$ are evaluated both their arguments have well defined finite values so that the same is true for their sum and difference. On the other hand, the expression $\min(1,1/\mathbf{abs}(u))$ makes perfect sense in IEEE arithmetic [iee85], but rewriting it as above leads to a $NaN$ at $u=0$. While this restriction may appear quite technical, it imposes the requirement that all relevant quantities are well defined at least in some open neighborhood, which is exactly in the nature of piecewise differentiability. For an in-depth investigation of piecewise differentiable functions, see, e.g., the books by Kummer [Kum88] and Scholtes [Sch12].

The function $\Delta F(\mathring{x};\Delta x)$ is incremental in that $\Delta F(\mathring{x};0)=0$, but like general piecewise linear continuous functions, it is only locally and positively homogeneous so that

[TABLE]

Here the bound $\rho(\mathring{x})$ is positive everywhere, but generally not continuous with respect to $\mathring{x}$.

In [Gri13] it has been shown that the Jacobians of the linear pieces of $\Delta F(\mathring{x};\alpha\Delta x)$ in the ball of radius $\rho(\mathring{x})$ about $\mathring{x}$ are conically active generalized Jacobians of the underlying nonlinear function $F(x)$. We will not elaborate on this connection here, because even in the smooth case secant approximations need not correspond to exact Jacobians. In contrast, the generalized derivative sets in the sense of Clarke [Cla83] are for piecewise differentiable functions almost everywhere just singletons, containing the classical Jacobian matrix. In floating point arithmetic the user or client algorithm will then quite likely never ’see’ the nonsmoothness or gain any useful information about it.

2 Piecewise Linearization by Tangents and Secants

Suppose the vector function $F:{\cal D}\subset\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ in question is evaluated by a sequence of assignments

[TABLE]

Here $\circ\in\{+,-,*,/\}$ is a polynomial arithmetic operation and

[TABLE]

a univariate function. The user or reader may extend the library by other locally Lipschitz-continuously differentiable functions like the analysis favorites

[TABLE]

To fit into the the framework they then also have to supply an evaluation procedure for both the elemental function $\varphi$ and its derivative $\varphi^{\prime}$, which cannot be based mechanically on the chain rule.

Following the notation from [GW08] we partition the sequence of scalar variables $v_{i}$ into the vector triple

[TABLE]

such that $x\in\mathbb{R}^{n}$ is the vector of independent variables, $y\in\mathbb{R}^{m}$ the vector of dependent variables and $z\in\mathbb{R}^{l-m}$ the (internal) vector of intermediates.

Some of the elemental functions like the reciprocal $rec(u)\equiv 1/u$, the square root and the logarithm are not globally defined. As mentioned above, we will assume that the input variables $x$ are restricted to an open domain ${\cal D}\subset\mathbb{R}^{n}$ such that all resulting intermediate values $v_{i}=v_{i}(x)$ are well defined.

Throughout we will assume that the evaluation procedure for $F$ involves exactly $s\geq 0$ calls to $\mathbf{abs}()$, including $\min$ and $\max$ rewritten or at least reinterpreted as discussed above. Starting from $\mathring{x}$ and an increment $\Delta x=x-\mathring{x}$, we will now construct for each intermediate $v_{i}$ an approximation

[TABLE]

Here the incremental function $\Delta v_{i}(\Delta x)$ is continuous and piecewise linear, with $\mathring{x}$ or $\check{x}$ and $\hat{x}$ considered constant in the tangent and secant case, respectively. Hence, we will often list $\Delta x$, but only rarely use $\mathring{x}$, $\check{x}$ and $\hat{x}$ as arguments of the $\Delta v_{i}$ in proofs.

Defining Relations for Tangent Approximation

We use the reference values $\mathring{v}_{i}=v_{i}(\mathring{x})$ and, assuming that all $\varphi_{i}$ other than the absolute value function are differentiable within the domain of interest, we may use the tangent linearizations

[TABLE]

Here $\mathring{c}_{ij}\equiv\varphi_{i}^{\prime}(\mathring{v}_{j})$, which will be different for the secant linearization.

If no absolute value or other nonsmooth elemental occurs, the function $y=F(x)$ is, at the current point, differentiable and by the chain rule we have the relation

[TABLE]

where $F^{\prime}(x)\in\mathbb{R}^{m\times n}$ is the Jacobian matrix. Thus we observe the obvious fact that smooth differentiation is equivalent to linearizing all elemental functions.

Now, let us move to the piecewise differentiable scenario, where the absolute value function does occur $s>0$ times. We may then obtain a piecewise linearization of the vector function $F(\mathring{x}+\Delta x)-F(\mathring{x})$ by incrementing

[TABLE]

Here $\mathring{v}_{i}=\mathbf{abs}(\mathring{v}_{j})$, which will be slightly different for the secant linearization. In other words, we keep the piecewise linear function $\mathbf{abs}()$ unchanged so that the resulting $\Delta y$ represents, for each fixed $x\in{\cal D}$, the piecewise linear and continuous increment function

[TABLE]

Defining Relations for Secant Approximation

In the tangent approximation the reference point was always the evaluation point $\mathring{x}$ and the resulting values $\mathring{v}_{i}=v_{i}(\mathring{x})$. Now we will make reference to the midpoints

[TABLE]

Consequently, we have the functional dependence $\mathring{v}_{i}\,=\,\mathring{v}_{i}(\check{x},\hat{x})$, which is at least Lipschitz continuous under our assumptions. Now an intriguing observation is that the recurrences (2) and (3) for arithmetic operations can stay just the same, and the recurrence (4) for nonlinear univariates is still formally valid, except that the tangent slope $\varphi^{\prime}(\mathring{v}_{j})$ must be replaced by the secant slope

[TABLE]

Theoretically, some $\check{v}_{i}$ and $\hat{v}_{i}$ may coincide, even if the underlying sample points $\check{x}$ and $\hat{x}$ are not selected identically, in which case the secant based model would reduce to the tangent based model. While exact coincidence of any pair $\check{v}_{i}$ and $\hat{v}_{i}$ is rather unlikely, taking the divided difference over small increments is likely to generate numerical cancellation. Therefore we will develop a division free centered form in Section 5. Finally, the nonsmooth rule (5) can stay unchanged except that we now set

[TABLE]

Hence, it is immediately clear that the new secant approximation reduces to the old tangent approximation when $\check{x}=\hat{x}$. In general, we will denote the mapping between the input increments $\Delta x\in\mathbb{R}^{n}$ and the resulting values $\Delta y\in\mathbb{R}^{m}$ by

[TABLE]

Its piecewise linear structure is very much the same as that of the tangent based model, which is described in detail in [Gri13]. Here we emphasize its quality in approximating the underlying nonlinear and nonsmooth $F$.

In contrast to the tangent model, the secant model is not a priori unique in that it depends quite strongly on the procedural representation of the vector function $F$ and not just its values, i.e. its properties as a mathematical map. For example, one can easily check that applying the above secant modeling rules to $f(x)=\log(\exp(x))$ does not yield the same approximation $\Delta f(\check{x},\hat{x},\Delta x)$ as the one for $f(x)=x$. On the other hand the natural secant linearization rule for the product $v=u\cdot w$ is equivalent to that obtained by applying the Appolonius identity

[TABLE]

Of course, the same is true for the tangent linearization and we may assume without loss of generality that we only have three kinds of elemental functions, the addition, the modulus and smooth univariate functions. That reduction greatly simplifies the thoretical analysis but might not be numerically optimal for actual implementations.

3 Approximation and Stability Properties

In contrast to the presentation in our previous papers we will now also use the nonincremental forms

[TABLE]

and

[TABLE]

For the square as a univariate nonlinear function $v(x)=x^{2}$ we find:

[TABLE]

Lemma 3.1**.**

Plugging the secant approximation for the square into the Appolonius identity, we obtain for the general multiplication $v=u\cdot w$:

[TABLE]

Proof.

[TABLE]

∎

Hereafter we will denote by $\lVert\cdot\rVert\equiv\lVert\cdot\rVert_{\infty}$ the infinity norm. Due to the norm equivalence in finite dimensional spaces all inequalities to be derived take the same form in other norms, provided the constants are adjusted accordingly. The infinity norm is particularly convenient, since we can then prove the following result for the general vector case $m>1$ by considering the absolute values of the individual components $f=F_{i}$ for $i=1\ldots m$. Moreover, we will make use of the Appolonius identity in the following proposition:

Proposition 3.2**.**

Suppose $\tilde{x},\check{x},\hat{x},\check{y},\hat{y},\check{z},\hat{z}\in\mathbb{R}^{n}$ are restricted to a sufficently small closed convex neighboorhod $K\subset\mathbb{R}^{n}$ where the evalution procedure for $F:\mathbb{R}^{n}\mapsto\mathbb{R}^{m}$ is well defined. Then there are Lipschitz constants $\beta_{F}$ and $\gamma_{F}$ such that we have

1. (i)

Lipschitz continuity of function, tangent and secant models:

[TABLE]

The constant $\beta_{F}$ can be defined by the recurrences $\beta_{v}=\beta_{u}+\beta_{w}$ if $v=u+v$, $\beta_{v}=\beta_{u}$ if $v=|u|$ and

[TABLE] 2. (ii)

Error between function and secant or tangent model:

[TABLE]

where $x\in K$. The constant $\gamma_{F}$ can be defined using the recurrences $\gamma_{v}=\gamma_{u}+\gamma_{w}$ if $v=u+w$, $\gamma_{v}=\gamma_{u}$ if $v=|u|$ and

[TABLE] 3. (iii)

*Lipschitz continuity of secant and tangent model: *

[TABLE]

where $x\in\mathbb{R}^{n}$. 4. (iv)

Lipschitz continuity of the incremental part: Let $x\in K$. Abbreviating $\Delta y=\hat{y}-\check{y}$ and $\Delta z=\hat{z}-\check{z}$ we obtain in the secant case

[TABLE]

which reduces in the tangent case to

[TABLE]

*Here the second bounds applying to the tangent model are always specializations of the previous ones for the secant model. *

Proof.

Since otherwise the bounds can be applied componentwise we may assume without loss of generality that $F$ is a scalar function $f$ and the norm in the range is simply the absolute value $|\cdot|$. The proof proceeds by induction on the intermediate quantities $v$ in the computational graph of $f$. We will define the constants $\beta_{v}$ and $\gamma_{v}$ recursively on the basis of the Lipschitz constants of the elemental functions and their derivatives:

Variable initialization: The initialization of independant variables represent the minimal nodes and all assertions are tivially true with the constants $\beta_{x_{i}}=1$ and $\gamma_{x_{i}}=0$.

Smooth univariate operation: Let $v=\varphi(u)$, $\varphi\in\tilde{\Phi}$, be some elemental function in the computational graph of $f$:

[TABLE]

The first inequality holds due to the Lipschitz continuity of $\varphi$ and $\lozenge_{\check{u}}^{\hat{u}}\varphi$. The latter inequality is the induction hypothesis and since $L_{K}(\lozenge_{\check{u}}^{\hat{u}}\varphi)\leq L_{K}(\varphi)$ holds by the mean value theorem, we may set $\beta_{v}\equiv\beta_{u}L_{K}(\varphi)$.

Absolute value function and sum: The absolute value function $v=\mathbf{abs}(u)$ naturally maintains the Lipschitz constant and the addition summates them.

Thus we have established $(i)$ for all cases.

Now let us consider the approximation property $(ii)$. For additions $v=u+w$ we may set $\gamma_{v}=\gamma_{u}+\gamma_{w}$ and then have by the triangle inequality

[TABLE]

For the absolute value function $v=\mathbf{abs}(u)$ we may also set $\gamma_{v}=\gamma_{u}$, since

[TABLE]

For the univariate functions $v=\varphi(u)$ we have with $\tilde{u}\equiv\lozenge_{\check{x}}^{\hat{x}}u(x)$

[TABLE]

By the mean value theorem and the induction hypothesis, the first term is bounded by

[TABLE]

The second term represents the error in the Hermite interpolation of $\varphi$ between $\check{u}$ and $\hat{u}$. With $L_{K}(\varphi^{\prime})$ a Lipschitz constant of $\varphi^{\prime}$ on $u(K)$ it is bounded by

[TABLE]

where the last bound follows from the fact that according to $(i)$ the approximation $\lozenge_{\check{x}}^{\hat{x}}u(x)$ has the Lipschitz constant $\beta_{u}$ and takes on at $\check{x}$ and $\hat{x}$ the values $\check{u}$ and $\hat{u}$. Hence, we have shown that (ii) holds indeed with

[TABLE]

Next we want to prove Lipschitz continuity of the model as stated in $(iii)$. Again, we find for additions and the abs function that the assertion is almost trivial with the constants $\gamma$ either being summated or just passed on. The challenge is once more the induction through the nonlinear univariates $v=\varphi(u)$. To limit the notational complexity we will connect the two point pairs at the $u$ level by straight lines setting

[TABLE]

Here $x$ is fixed and we assume as induction hypothesis that

[TABLE]

To connect the piecewise linearizations of $v$ define

[TABLE]

The quantity we want to find a bound for is $\Delta\tilde{v}\equiv\tilde{v}(1)-\tilde{v}(0)=\lozenge_{\check{z}}^{\hat{z}}v(x)-\lozenge_{\check{y}}^{\hat{y}}v(x).$ By the mean value theorem we find some $\bar{t}\in[0,1]$ where

[TABLE]

The functions $\check{u},\hat{u},\tilde{u}$ here and following are to be read as evaluated at $t=\bar{t}$. Now introduce $\bar{u}$ as the mean value of $\check{u}$ and $\hat{u}$ at which the difference quotient of $\varphi$ over the intervening interval is equal to its derivative. This yields:

[TABLE]

The two quotients are bounded according to

[TABLE]

where the last equality follows from $\bar{u}$ being between $\check{u}$ and $\hat{u}$. Hence, we find that

[TABLE]

The factors in the middle are easily bounded by

[TABLE]

That leaves us with the second factors, which are linear in $t$ such that

[TABLE]

where the last inequality follows from (i). Analogously, we can derive

[TABLE]

Substituting this into (12) we get

[TABLE]

with $\gamma_{v}\equiv L_{K}(\varphi)\gamma_{u}+L_{K}(\varphi^{\prime})\beta_{u}^{2}$, which completes the proof of $(iii)$. Finally, we have to prove $(iv)$, which gives a bound on the increment part only. Setting $\xi\equiv\tfrac{1}{2}(\mathring{z}+\mathring{y})+\Delta x$ one gets, with the results already proved and a few triangle inequalities, that

[TABLE]

Now, since for example

[TABLE]

both inner maxima can be bounded by the same expression, namely

[TABLE]

so that we obtain the upper bound

[TABLE]

Finally, we can also bound

[TABLE]

which yields the assertion after some elementary modifications. ¿From the secant result we can easily get the bound for the tangent model by setting $\check{z}=\hat{z}=\mathring{z}$ and $\check{y}=\hat{y}=\mathring{y}$. ∎

As one can see by setting $\check{y}=x=\check{z}$ the assertion $(ii)$ almost follows from $(iii)$, except that a factor of 2 is lost in the constants. The proposition also states that the values of $F$ at $\check{x}$ and $\hat{x}$ are reproduced exactly by our approximation as one would expect from a secant approximation. This property clearly nails down the piecewise linearization rules (5) and (4) with (7) for all univariate functions. Also, there is no doubt that addition and subtraction should be linearized according to (2) and that multiplications $v_{i}=c\,v_{j}$ by constants $c$ should yield the differentiated version $\Delta v_{i}=c\,\Delta v_{j}$, which is a special case of (3). For general multiplications $v_{i}=v_{j}\ast v_{k}$ the two values $\check{v}_{i}=\check{v}_{j}\ast\check{v}_{k}$ and $\hat{v}_{i}=\hat{v}_{j}\ast\hat{v}_{k}$ could also be interpolated by linear functions other than the one defined by (3). However, we currently see no possible gain in that flexibility, and maintaing the usual product rule form seems rather attractive.

Stable Surjectivity

Recall that a continuous function $F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$ is called proper if the preimage of every compact set is compact. One can easily see that any piecewise linear function is proper if and only if it maps no affine ray $\{a+\lambda b\;:\;\lambda\geq 0\},(a,b\in\mathbb{R}^{n}),$ to a point. This is trivially the case if $n=m$ and the Jacobians of all selection functions of $F$ are invertible. In [Rad17] it was proved that if a piecewise linear function $F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is proper, there exists a $d\in\mathbb{Z}$ such that for all regular values $y$ of $F$, i.e. all values $y$ such that the Jacobian at all preimages of $y$ exists and is invertible [which is trivially the case if $F^{-1}(y)=\emptyset$], it holds

[TABLE]

Hence, any regular value $y$ of a proper piecewise linear function $F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ has at least one preimage if $\operatorname{deg}(F)\neq 0$. This implies surjectivity of $F$ by the closedness of piecewise linear functions (cf. [Sch12]) and the well known fact that regular values lie dense in the range. We call $\operatorname{deg}(F)$ the degree of $F$.

Lemma 3.3**.**

([Rad17, Cor. 5.2]) Assume a piecewise linear function $F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is composed of $k$ affine functions $F_{i}$ with invertible linear parts $A_{i}$ and define

[TABLE]

Moreover, let $0<\varepsilon<\rho_{F}$ and $l\geq 0$. Then all piecewise linear functions $G:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ with

[TABLE]

are proper and their degree equals that of $F$.

The latter statement enables us to prove another stability result:

Proposition 3.4**.**

Let $F\in\text{span}(\Phi)$ where $F:\mathbb{R}^{n}\to\mathbb{R}^{n}$ and assume a piecewise linearization $\lozenge_{\mathring{x}}F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is composed of $k$ affine functions $F_{i}$ with invertible linear parts $A_{i}$. Moreover, define $\rho_{F}$ as in Lemma 3.3. Then for all tangent and secant mode piecewise linearizations with development points in the ball $B(\mathring{x},\rho_{F}/\gamma_{F})$, where $\gamma_{F}$ is defined as in Proposition 3.2, the mapping degree is well defined and equals that of $\lozenge_{\mathring{x}}F$.

Proof.

Just plug Proposition 3.2. $(iii)$ into Lemma 3.3. ∎

The fact that the mapping degree of a coherently oriented piecewise linear function cannot be zero yields a rather pointed statement:

Corollary 3.5**.**

In the situation of Proposition 3.4, if $\lozenge_{\mathring{x}}F$ is coherently oriented, then all tangent and secant mode piecewise linearizations with development points in the ball $B(\mathring{x},\rho_{F}/\gamma_{F})$ are surjective.

4 Generalized Newton Methods by Piecewise Linearization

We will proceed by proposing and analyzing a possible application of the tangent and secant approximations developed in the previous sections. For this we present generalized Newton’s methods for composite piecewise smooth functions $F:\mathbb{R}^{n}\to\mathbb{R}^{n}$ based on the piecewise linear approximations, both for the tangent and the secant mode. The merit of these methods is the fact that they impose no strong differentiability requirements but require only piecewise differentiability at the root in question.

Definition 1** (Newton operator).**

Let $F\in{\rm span}(\Phi)$ and $x^{\ast}$ be an isolated root of $F$ in an open neighborhood ${\cal D}$. The Newton step for $F$ is definable on ${\cal D}$ in tangent or secant mode if the piecewise linear equation $\lozenge_{\mathring{x}}F(x)=0$ resp. $\lozenge_{\check{x}}^{\hat{x}}F(x)=0$ has at least one root for all $\mathring{x}\in\mathcal{D}$ resp. $\check{x},\hat{x}\in{\cal D}$.

Then the Newton operator is defined in tangent mode as

[TABLE]

and in secant mode as

[TABLE]

Definition 1 is the minimal assumption under which we can conclude, using the approximation and Lipschitz continuity results of Section 3, that the iteration

[TABLE]

converges locally to an isolated root $x^{\ast}$. Assuming that the piecewise linearization $\lozenge_{x^{\ast}}F$ is bijective on a ball ${\cal D}=B(x^{\ast},\bar{R})$ we will proceed to show that close to $x^{\ast}$

• •

the Newton step can be defined

• •

the Newton step stays close to $x^{\ast}$

• •

the tangent mode Newton method converges quadratically and the secant mode Newton method converges with the golden mean as order.

We know that the assumption of bijectivity on a ball can be relaxed through the use of degree theory, but this requires a significant technical effort for which we refer to subsequent works. Note that the minimal assumption for the local convergence of semismooth Newton [Hin10] is similar, namely, that there always exists a generalized Jacobian $J(x)$ that has a uniformly bounded inverse over all $x\in B(x^{\ast},\bar{\rho})$.

For most of the following results up to the last it is sufficient to consider the secant mode piecewise linear approximation of $F$, as results for the tangent mode can be obtained by setting $\check{x}=\hat{x}=\mathring{x}$.

The general bound for the difference of piecewise linearizations with different basis points $\check{z},\hat{z}$ and $\check{y}=\hat{y}=\mathring{y}$ can be refined to a closer bound in the case that one is to express that bound only in distances to $\mathring{y}$.

Lemma 4.1**.**

Let the second order constant $\gamma_{F}$ be valid on some convex set $U$ and $x,\mathring{y},\check{z},\hat{z}\in U$. Then

[TABLE]

Proof.

Select some $N\in\mathbb{N}$ and consider the subdivision of the segments $[\check{z},\mathring{y}]$ and $[\hat{z},\mathring{y}]$ by $\check{u}_{k}=\mathring{y}+\frac{k}{N}(\check{z}-\mathring{y})$ and $\hat{u}_{k}=\mathring{y}+\frac{k}{N}(\hat{z}-\mathring{y})$. Then by Prop. 3.2.$(iii)$

[TABLE]

Summation and limit $N\to\infty$ results in the claim. ∎

We assume hereafter that $\lozenge_{x^{\ast}}F$ is bijective on a ball ${\cal D}=B(x^{\ast},\bar{R})$, where $x^{\ast}$ is an isolated root of $F$, which implies that it is also coherently oriented and metrically regular on the latter. This means that there exists a constant $c>0$ such that for all $x\in{\cal D}$, $y\in\lozenge_{x^{\ast}}F({\cal D})$

[TABLE]

Lemma 4.2** (existence of roots of the PL approximation).**

If $R\leq\min(\bar{R},\frac{1}{c\gamma_{F}})$ and $\rho=\frac{R}{3}$ then for any $\check{x},\hat{x}\in B(x^{\ast},\rho)$ the piecewise linearization $\lozenge_{\check{x}}^{\hat{x}}F$ has a root in $B(x^{\ast},R)$.

Proof.

We intend to apply the Brouwer theorem to the fixed point operator

[TABLE]

Any fixed point is then a root of $\lozenge_{\check{x}}^{\hat{x}}F$. As $T$ is obviously continuous, we only need to show that $T$ maps $B(x^{\ast},R)$ into itself. Using metric regularity and Lemma  4.1 we find

[TABLE]

∎

Lemma 4.3** (containment of the roots).**

Under the same assumptions any root of $\lozenge_{\check{x}}^{\hat{x}}F$ in $B(x^{\ast},R)$ is actually contained in $B(x^{\ast},\rho)$.

Proof.

Let $x$ be a root of $\lozenge_{\check{x}}^{\hat{x}}F$. Then again using metric regularity and Lemma 4.1 the distance to $x^{*}$ has the bound

[TABLE]

so that

[TABLE]

∎

Corollary 4.4** (Newton iteration).**

Let $\lozenge_{x^{\ast}}F$ be bijective on a ball ${\cal D}=B(x^{\ast},\bar{R})$, where $x^{\ast}$ is an isolated root of $F\in{\rm span}(\Phi)$. Then on $B(x^{\ast},\rho)$, where $\rho$ is defined as in Lemma 4.2, the Newton step is defined in both the tangent and secant mode and maps back into $B(x^{\ast},\rho)$. The thus definable Newton iteration converges at least linearly.

Proof.

As the next Newton iterate is among the roots of minimal distance to the basis point(s) of the linearization, one needs to ensure that any root that $\lozenge_{\check{x}}^{\hat{x}}F$ may have outside the ball $B(x^{\ast},R)$ has a larger distance to $\{\check{x},\hat{x}\}$ than the root that is known to exist inside $B(x^{*},\frac{1}{4}\rho)$. The distance from the basis points inside $B(x^{*},\rho)$ to the outside of $B(x^{*},R)$ is at least $R-\rho=\frac{2}{3}R$. The distance from the root of $\lozenge_{\check{x}}^{\hat{x}}F$ inside $B(x^{*},\frac{1}{4}R)$ to the basis points is at most $\frac{1}{4}\rho+\rho=\frac{5}{12}R$ and thus the smaller distance.

By equation (15) of the last lemma we also see that the distance of the root of $\lozenge_{\check{x}}^{\hat{x}}F$ to $x^{*}$ is at most $\frac{1}{4}$ the distance of the basis points to $x^{*}$, which implies linear convergence. ∎

Corollary 4.5** (convergence rates).**

Let $R=\min(\bar{R},\frac{1}{c\gamma_{F}})$ and $\rho=\frac{R}{3}$. Then for all initial points $x_{0}\,(,x_{1})\in B(x^{\ast},\rho)$ the Newton iteration in tangent mode

[TABLE]

converges quadratically resp. in secant mode

[TABLE]

converges with order $\frac{1+\sqrt{5}}{2}$ towards $x^{\ast}$.

Proof.

With the choice of the initial points and by Lemma 4.3, the full iteration sequence stays inside $B(x^{*},\rho)$.

Tangent mode: By replacing $x,\check{x},\hat{x}$ with $x_{j+1},x_{j},x_{j}$ equation (15) in Lemma 4.3 we get

[TABLE]

which implies the quadratic convergence of the sequence $(x_{j})$ towards $x^{*}$,

[TABLE]

Secant mode: Replacing $x,\check{x},\hat{x}$ with $x_{j+2},x_{j+1},x_{j}$ in equation (15) in Lemma 4.3 one finds

[TABLE]

As in the scalar secant method, this implies convergence with rate $\phi=\frac{1+\sqrt{5}}{2}$ or more precisely

[TABLE]

where $(F_{j})$ is the Fibonacci sequence with $F_{0}=0$, $F_{1}=1$ and $F_{j+1}\leq\phi^{j}$. ∎

The inner iterations of (13) require the solution of piecewise linear equations. However, these may possess several solutions $\Delta x_{j}$, of which we must find one of minimal norm. So far, the solvers that we know and surveyed in [Gri13, GBRS15, SGRB14, Rad16] require at least coherent orientation to guarantee convergence. Hence, the actual implementation of successive piecewise linearization needs further study.

5 Singularity Free Implementation

In contrast to the tangent mode of piecewise linearization, the secant mode involves two points of evaluation $\check{v}_{i},\hat{v}_{i}\in\mathbb{R}^{n}$, which means that its computational cost of the primal values are roughly twice that of the tangent mode. These define a line segment $[\check{v}_{i},\hat{v}_{i}]\equiv\{\lambda\check{v}_{i}+(1-\lambda)\hat{v}_{i}\mid\lambda\in[0,1]\}$ for any intermediate operation. The formal definition of the secant slope given in equation (7) may cause numerically unstable divisions when the denominator gets small during the transition from secant to tangent mode, e.g. when the secant mode Newton iteration scheme converges.

However, in this section we will provide singularity free closed form expressions. Therefore an exception handling at $\check{v}_{i}=\hat{v}_{i}$ will no longer be necessary. To that end we move from the line segment representation to a midpoint-radius based representation. Now let $v=\varphi(u)$, where $\varphi\in\{\sin,\exp,\dots\}$ is some elementary operation and

[TABLE]

We adopted the concept of representing intervals via midpoint and radius from interval arithmetic calculus (described in detail e.g. by Siegfrid Rump in [Rum99], or by Götz Alefeld and Jürgen Herzberger in [AH12]). We remark though, that in the present setting the radius $\delta\!v_{i}\in\mathbb{R}$ is allowed to become negative as well.

Now one can rewrite the secant slope of differentiable functions to the new representation

[TABLE]

Using some algebraic manipulations one can find individual formulas for the aforementioned propagation rules of the secant mode:

[TABLE]

Alternatively, we could represent $v_{i}=v_{j}/v_{k}$ as an application of a multiplication on the reciprocal $1/v_{k}$. Furthermore, we can represent the multiplication by the Appolonius identity as above. Moreover, for unary operations we get:

[TABLE]

Note that by, e.g., [Wei16] $\operatorname{\mathbin{\operatorname{sinc}}}(x)$ and $\operatorname{\mathbin{\operatorname{sinhc}}}(x)$ (hyperbolic $\operatorname{\mathbin{\operatorname{sinc}}}(x)$) have regular Taylor expansions

[TABLE]

We want to define $\operatorname{\mathbin{\operatorname{artanhc}}}(x)$ similar to $\operatorname{\mathbin{\operatorname{sinhc}}}(x)$ via its Taylor expansion:

[TABLE]

and it can be implemented in a similar fashion as $\operatorname{\mathbin{\operatorname{sinhc}}}(x)$ from the boost c++ libraries (see [lib]). For Root functions ($v_{i}=\sqrt[c]{v_{j}}$), general Powers ($v_{i}=v_{j}^{c}$, or in a binary fashion $v_{i}=v_{j}^{v_{k}}$) and monomials ($v_{i}=v_{j}^{n}$) one can use the identity

[TABLE]

and apply the rules above. Of course, the base $v_{j}>0$ has to be positive, but there is a less restrictive alternative for monomials:

[TABLE]

Remark 1* (General approximation for unary operations).*

Of course one can find a lot more singularity free formulas for the secant slopes of other operations. Using a Taylor expansion approach one can provide general approximation formulas for the triplet $\mathring{v}_{i},\delta\!v_{i}$ and $c_{ij}$ by:

[TABLE]

6 Numerical Example

Consider the function

[TABLE]

where $c=[1.001,10.01]^{\top}$, and $\angle\!x\in{[0,2\pi[}$ is the angle of $x=(x_{1},x_{2})^{\top}$ in polar coordinate representation. Moreover, $\varphi$, which is defined by

[TABLE]

maps ${[0,2\pi[}$ strictly monotonically onto itself. Hence, $F(x)$ is bijective. Furthermore, the function is differentiable everywhere except at the origin. There it is, just as the Euclidean norm, locally Lipschitz continuous and not even piecewise differentiable in the sense of [Sch12].

We investigated the behavior of the tangent and secant mode Newton on $F$ both with and without noise. That is, we investigated $F$ and $\tilde{F}$, where

[TABLE]

The secant mode was started with the initial values $\check{x}=[-3.7,-2.05]^{\top}$ and $\hat{x}=[7.0,8.0]^{\top}$. The tangent mode was started with the mean value of the latter points.

We recall the well known formula for approximating the convergence rate numerically:

[TABLE]

The first table shows both modes’ residuals in the iterations without noise.

[TABLE]

The next table shows the residuals of the iterations with noise.

[TABLE]

Both methods attain their theoretical convergence rates in either iteration. However, we observe that the secant mode fares better with the problem with noise as it cuts through the latter, while the tangent Newton is thrown off for several iteration steps.

7 Final Remarks

The framework developed in the present paper, as well as in [Gri13, GBRS15, SGRB14, Rad16, GHRS16], aims at presenting viable alternatives to current approaches to piecewise differentiability as it may occur, e.g., in nonsmooth nonlinear systems or ODEs with nonsmooth right hand side. The piecewise linearizations, which were first introduced in [Gri13] can be obtained in an automated fashion by an adaptation of AD tools such as ADOL-C [WG12].

The generalized Newton iterations introduced in Section 5 are intended as an alternative to semismooth Newton [Hin10]. The quadratic convergence rate of the tangent version is in line with that of semismooth Newton. While for one step of semismooth Newton an appropriate element of the generalized derivative has to be calculated, the piecewise linear Newton’s methods solve a piecewise linear system in each step. Both problems are NP-hard in general, but may be solvable in practice with essentially the same effort as an ordinary linear system. Hence, it is likely highly situation dependent, which approach yields better performance.

It is our hope that, in combination with the formulas for numerically stable implementation of the secant linearization, the generalized Newton’s methods can be developed into robust and stable workhorse algorithms, which might even outperform the quadratic methods on selected problems. For example on problems with oscillating noise we observed that the secant method required fewer Newton steps to converge, as it cuts through the oscillations.

Acknowledgements

The work for the article has been partially conducted within the Research Campus MODAL funded by the German Federal Ministry of Education and Research (BMBF) (fund number $05\text{M}14\text{ZAM}$).