# On the Frank-Wolfe algorithm for non-compact constrained optimization
problems ††thanks: The first author was supported in part by CNPq grants
305158/2014-7 and 302473/2017-3, FAPEG/PRONEM- 201710267000532 and CAPES. The
second author was supported in part by Fundação de Apoio à Pesquisa do
Distrito Federal (FAP-DF) by the grant 0193.001695/2017 and PDE 05/2018

O. P. Ferreira Instituto de Matemática e Estatística, Universidade Federal de
Goiás, CEP 74001-970 - Goiânia-GO, Brazil<EMAIL_ADDRESS>W. S. Sosa
Programa de Pós-Graduação em Economia, Universidade Católica de Brasília, CEP
70790-160 - Brasília-DF, Brazil, E-mail<EMAIL_ADDRESS>
###### Abstract

This paper is concerned with the Frank–Wolfe algorithm for a special class of
non-compact constrained optimization problems. The notion of asymptotic cone
is used to introduce this class of problems as well as to establish that the
algorithm is well defined. These problems, with closed and convex constraint
set, are characterized by two conditions on the gradient of the objective
function. The first establishes that the gradient of the objective function is
Lipschitz continuous, which is quite usual in the analysis of this algorithm.
The second, which is new in this subject, establishes that the gradient
belongs to the interior of the dual asymptotic cone of the constraint set.
Classical results on the asymptotic behavior and iteration-complexity bounds
for the sequence generated by the Frank–Wolfe algorithm are extended to this
new class of problems. Examples of problems with non-compact constraints and
objective functions satisfying the aforementioned conditions are also
provided.

Keywords: Frank-Wolfe method; constrained optimization problem; non-compact
constraint.

AMS subject classification: 90C25, 90C60, 90C30, 65K05.

##  1 Introduction

The Frank–Wolfe algorithm or conditional gradient method is one of the oldest
methods for finding minimizers of differentiable functions on compact convex
sets. It was initially proposed in 1956 [1] for solving quadratic programming
problems with linear constraints (see also [2, 3]), and it has since attracted
considerable attention owing to its simplicity and ease of implementation, as
it only requires access to a linear minimization oracle over the constraint
set. In particular, it allows low storage cost and readily exploits
separability and sparsity; therefore, it can be effectively applied to large-
scale problems. With the emergence of machine-learning applications, this
method has recently gained increasing popularity [4, 5, 6]. Accordingly, it
has been extensively studied, and several variants thereof have been developed
[7, 8, 9, 10, 11, 12, 13, 14] and references therein.

The aim of this study is to extend the Frank–Wolfe algorithm to a special
class of constrained optimization problems ${\rm Minimize}_{x\in{C}}f(x)$,
where $f:{\mathbb{R}}^{n}\to{\mathbb{R}}$ is a continuously differentiable
function, and $C\subset{\mathbb{R}}^{n}$ is a closed and convex but not
necessarily compact set. In addition to the classical assumptions (i.e., the
gradient of $f$ is Lipschitz continuous), we assume that $\nabla
f(x)\in\mbox{int}(C_{\infty})^{*}$ for all $x\in C$, where $\nabla f$ and
$\mbox{int}(C_{\infty})^{*}$ denote the gradient of $f$ and the interior of
the positive dual asymptotic cone of $C$, respectively. For this class of
functions, with classical assumptions, we also extend previous results on the
asymptotic behavior and iteration-complexity bounds for the sequence generated
by the Frank–Wolfe algorithm.

This paper is organized as follows. In Section 2, notations and auxiliary
results are presented. In Section 3, we formulate the Frank–Wolfe algorithm.
In Section 3.1, we establish that the sequence generated by this algorithm is
well defined. Section 3.2 is devoted to the study of the asymptotic-
convergence properties of Algorithm 1, and Section 3.3 to the study of the
iteration-complexity bounds. In Section 4, we present some examples. We
conclude the paper in Section 5.

##  2 Preliminaries

Herein, we present notations, definitions, and auxiliary results. Let
${\mathbb{R}}^{n}$ be the $n$-dimensional Euclidean space with the usual inner
product $\langle\cdot,\cdot\rangle$ and norm $\|\cdot\|$. We denote by
${\mathbb{R}}^{m\times n}$ the set of all $m\times n$ matrices with real
entries (${\mathbb{R}}^{n}\equiv{\mathbb{R}}^{n\times 1}$), by $e^{i}$ the
$i$-th canonical unit vector in ${\mathbb{R}}^{n}$, and by ${\rm I_{n}}$ the
$n\times n$ identity matrix. A set ${\cal K}\subseteq{\mathbb{R}}^{n}$ is
called a cone if for any $\alpha>0$ and $x\in\cal{K}$, we have $\alpha
x\in\cal{K}$. A cone ${\cal K}\subseteq{\mathbb{R}}^{n}$ is called convex if
for any $x,y\in\cal{K}$, we have $x+y\in\cal{K}$. The positive dual cone of a
cone ${\cal{K}}\subseteq{\mathbb{R}}^{n}$ is the cone
${\cal{K}}^{*}\\!\\!:=\\!\\{x\in{\mathbb{R}}^{n}:~{}x^{{\scriptscriptstyle\mathrm{T}}}y\geq\\!0,~{}\forall\,y\\!\in\\!{\cal{K}}\\}$,
and its interior is denoted by
$\mbox{int}{\cal{K}}^{*}\\!\\!:=\\!\\{x\in{\mathbb{R}}^{n}:~{}x^{{\scriptscriptstyle\mathrm{T}}}y>0,~{}\forall\,y\\!\in\\!{\cal{K}}{\setminus\\{0\\}}\\}$.
Let $C\subset{\mathbb{R}}^{n}$ be a closed convex set; then, we define the
asymptotic cone of $C$ by

$C_{\infty}:=\Big{\\{}d\in{\mathbb{R}}^{n}:~{}\exists~{}(t_{k})_{k\in\mathbb{N}}\subset(0,\infty),~{}\exists~{}(x^{k})_{k\in\mathbb{N}}\subset
C;\lim_{k\to\infty}t_{k}=0,\lim_{k\to\infty}t_{k}x_{k}=d\Big{\\}},$

or equivalently, $C_{\infty}:=\\{d\in{\mathbb{R}}^{n}:~{}x+td\in
C,\forall~{}x\in C,~{}\forall~{}t\geq 0\\}$ [15, pp. 39]. Let
$f:{\mathbb{R}}^{n}\to{\mathbb{R}}$ be a continuously differentiable function.
We consider the problem of determining an optimum point of $f$ in a closed
convex set ${C}\subset{\mathbb{R}}^{n}$, that is, a point $x^{*}\in{C}$ such
that $f(x^{*})\leq f(x)$ for all $x\in{C}$. We denote this constrained problem
as

$\displaystyle{\rm Minimize}_{x\in{C}}f(x).$ (1)

The optimal value of $f$ on ${C}$ is denoted by $f^{*}$, that is,
$f^{*}:=\inf_{x\in{C}}f(x)$. The first-order optimality condition  for problem
(1) is stated as

$\nabla f(x^{*})^{T}(x-x^{*})\geq 0,\qquad\forall~{}x\in{C}.$ (2)

In general, the condition (2) is necessary but not sufficient for optimality.
Thus, a point $x^{*}\in C$ satisfying condition (2) is called a stationary
point to problem (1).

######  Definition 2.1.

Let $C\subset{\mathbb{R}}^{n}$ be a convex set. A function
$f:{\mathbb{R}}^{n}\to{\mathbb{R}}$ is called $M$-strongly convex with
parameter $M\geq 0$ on $C$ if the inequality $f(tx+(1-t)y)\leq
tf(x)+(1-t)f(y)-\frac{1}{2}Mt(1-t)\|y-x\|^{2}$ holds for all $t\in[0,1]$ and
$x,y\in{C}$. In particular, for $M=0$, $f$ is called convex rather than
$0$-strongly convex.

The following results provide a useful characterization of convex/strongly
convex differentiable functions, see the proof in [16, Theorem 4.1.1, pp.
183].

######  Proposition 2.2.

Let $C\subset{\mathbb{R}}^{n}$ be a convex set,
$f:{\mathbb{R}}^{n}\to{\mathbb{R}}$ be a continuously differentiable function,
and $M\geq 0$. Then, $f$ is $M$-strongly convex in $C$ if and only if
$f(y)\geq f(x)+\nabla f(x)^{T}(y-x)+M\|x-y\|^{2}/2$ for all $x,y\in C$.

######  Remark 1.

It is well known that if $f$ is $M$-strong convex, then (2) is sufficient for
optimality, that is, any point $x^{*}\in C$ satisfying (2) is a minimizer to
problem (1).

The proof of the next lemma can be found in [17, Lemma 6].

######  Lemma 2.3.

Let $\\{a_{k}\\}$ be a nonnegative sequence of real numbers. If $\Gamma
a_{k}^{2}\leq a_{k}-a_{k+1}$ for some $\Gamma>0$ and for any $k=1,...,\ell$,
then $a_{\ell}\leq a_{0}/(1+\Gamma a_{0}\ell)<1/(\Gamma\ell)$.

##  3 Frank–Wolfe algorithm

Herein, we formulate the Frank–Wolfe algorithm to solve problem (1). To this
end, we henceforth assume that the constraint set $C\subset{\mathbb{R}}^{n}$
is closed and convex (not necessarily compact), the objective function
$f:{\mathbb{R}}^{n}\to{\mathbb{R}}$ of problem (1) is continuously
differentiable, and its gradient satisfies the following condition:

  * (A)

$\|\nabla f(x)-\nabla f(y)\|\leq L\|x-y\|$ for all $x,y\in{C}$ and $L>0$.

######  Remark 2.

In Section 4, we present examples of problem (1) with objective function
satisfying (A).

To formulate the Frank–Wolfe algorithm, we should assume that there exists a
linear-optimization oracle (LO oracle) capable of minimizing linear functions
over $C$.

######  Algorithm 1.

CondGC,f method

0.

Select $x^{0}\in{C}$. Set $k=0$.

1.

Use an “LO oracle” to compute an optimal solution $p^{k}$ and the optimal
value $v_{k}^{*}$ as

$p^{k}\in{\rm argmin}_{p\in{C}}\nabla f(x^{k})^{T}(p-x^{k}),\qquad
v_{k}^{*}:=\nabla f(x^{k})^{T}(p^{k}-x^{k}).$ (3)

2.

If $v^{*}_{k}=0$, then stop; otherwise, compute the step size
$\lambda_{k}\in(0,1]$ as

$\lambda_{k}:=\mbox{min}\left\\{1,\frac{|v_{k}^{*}|}{L\|p^{k}-x^{k}\|^{2}}\right\\}={\rm
argmin}_{\lambda\in(0,1]}\left\\{v_{k}^{*}\lambda+\frac{L}{2}\|p(x^{k})-x^{k}\|^{2}\lambda^{2}\right\\},$
(4)

and set the next iterate $x^{k+1}$ as

$x^{k+1}:=x^{k}+\lambda_{k}(p^{k}-x^{k}).$ (5)

3.

Set $k\leftarrow k+1$, and go to step 1.

We conclude this section by stating a basic inequality for functions
satisfying assumption (A) [18, Lemma 2.4.2].

######  Lemma 3.1.

Let $f:\mathbb{R}^{n}\to\mathbb{R}$ be a continuously differentiable function
satisfying condition (A), $x\in{C}$, and $\lambda\in[0,1]$. Then,

$f(x+\lambda(p-x))\leq f(x)+\nabla
f(x)^{{\scriptscriptstyle\mathrm{T}}}(p-x)\lambda+\frac{L}{2}\|p-x\|^{2}\lambda^{2},\qquad\forall~{}p\in{C}.$
(6)

###  3.1 Well-definedness

Herein, we establish that the sequence $(x^{k})_{k\in\mathbb{N}}$ generated by
Algorithm 1 is well defined. To this end, we assume that the gradient of the
objective function $f:{\mathbb{R}}^{n}\to{\mathbb{R}}$ of problem (1)
satisfies the following condition:

  * (B)

$\nabla f(x)\in\mbox{int}(C_{\infty})^{*}$, for all $x\in C$.

######  Remark 3.

It follows from [15, Proposition 2.2.3] that a closed convex set
$C\subset{\mathbb{R}}^{n}$ is compact if and only if $C_{\infty}=\\{0\\}$.
Then, $\mbox{int}(C_{\infty})^{*}={\mathbb{R}}^{n}$, and therefore (B) holds
trivially if $C$ is compact. In Section 4, we present examples of problem (1)
with $\nabla f$ satisfying (B) with unbounded constraint set $C$.

We now use ${\bf(B})$ to prove a general result that implies the existence of
a solution to the problem (3). It is worth mentioning that ${\bf(A})$ is not
required in the proof.

######  Proposition 3.2.

The following three assertions hold:

(i)

For each $x\in C$, the set $D_{x}:=\left\\{p\in C:~{}\nabla f(x)^{T}(p-x)\leq
0\right\\}$ is compact.

(ii)

For each $x\in C$, the linear problem

$\displaystyle{\rm Minimize}_{p\in{C}}\nabla f(x)^{T}(p-x)$ (7)

has a solution.

(iii)

If $D\subset C$ is a bounded set, then the set

$\bigcup_{x\in D}\Big{\\{}q_{x}\in C:~{}q_{x}\in{\rm argmin}_{p\in{C}}\nabla
f(x)^{T}(p-x)\Big{\\}},$ (8)

is also bounded.

###### Proof.

To prove (i), let $x\in C$. We assume toward a contradiction that $D_{x}$ is
unbounded. Thus, let $(q^{k})_{k\in\mathbb{N}}\subset D_{x}$ such that
$\lim_{k\to\infty}\|q^{k}\|=\infty$, and let
$(t_{k})_{k\in\mathbb{N}}\subset(0,\infty)$ be the sequence defined by
$t_{k}:=1/\|q^{k}\|$ for all $k=0,1\ldots$. Then, $\lim_{k\to\infty}t_{k}=0$.
As $t_{k}q^{k}=q^{k}/\|q^{k}\|$, we conclude that $\|t_{k}p_{k}\|=1$ for all
$k=0,1,\ldots$. Hence, there exist subsequences
$(q^{k_{j}})_{j\in\mathbb{N}}\subset D_{x}$ and
$(t_{k_{j}})_{j\in\mathbb{N}}\subset(0,\infty)$ such that
$\lim_{k_{j}\to\infty}t_{k_{j}}q^{k_{j}}=d\in C_{\infty}$. Thus, the
definition of $D_{x}$ implies

$\nabla f(x)^{T}\left(t_{k_{j}}q^{k_{j}}-t_{k_{j}}x\right)=\nabla
f(x)^{T}\left(\frac{q^{k_{j}}}{\|q^{k_{j}}\|}-\frac{x}{\|q^{k_{j}}\|}\right)\leq
0.$

Taking the limit in the last inequality as $j$ tends to $\infty$, we conclude
that $\nabla f(x)^{T}d\leq 0$, which is absurd. Indeed, assumption (B) implies
that $\nabla f(x)\in\mbox{int}(C_{\infty})^{*}$, and as $d\in C_{\infty}$, we
have $\nabla f(x)^{T}d>0$. Therefore, (i) is proved. To prove (ii), it is
sufficient to note that the problem (7) has $D_{x}$ as a sublevel set, which
by (i) is compact. We now prove (iii). We assume toward a contradiction that
the set in (8) is unbounded. Then, there exist sequences
$(x^{k})_{k\in\mathbb{N}}\subset D$ and $(q_{x^{k}})_{k\in\mathbb{N}}\subset
C$ such that $\lim_{k\to\infty}\|q_{x^{k}}\|=\infty$. Thus, as $D$ is bounded,
we have

$\lim_{k\to\infty}\tau_{k}=0,\qquad\mbox{where}\quad\tau_{k}:=\frac{1}{\|q_{x^{k}}-x^{k}\|},\qquad
k=0,1,\ldots.$ (9)

However, as $C$ is convex, and $(x^{k})_{k\in\mathbb{N}}$ and
$\\{q_{x^{k}}\\}$ belong to $C$, we have

$x^{k}+t\left(q_{x^{k}}-x^{k}\right)\in C,\qquad k=0,1,\ldots$ (10)

for any $t\in(0,1)$. As
$\tau_{k}(q_{x^{k}}-x^{k})=(q_{x^{k}}-x^{k})/\|q_{x^{k}}-x^{k}\|$, we have
$\|\tau_{k}(q_{x^{k}}-x^{k})\|=1$ for all $k=0,1,\ldots$. Thus, there exist
subsequences $(x^{k_{j}})_{j\in\mathbb{N}}\subset D$,
$(q_{x^{k_{j}}})_{j\in\mathbb{N}}\subset C$, and
$(t_{k_{j}})_{j\in\mathbb{N}}\subset(0,\infty)$ such that

$\lim_{k_{j}\to\infty}\tau_{k_{j}}(q_{x^{k_{j}}}-x^{k_{j}})=v.$ (11)

As $(x^{k_{j}})_{j\in\mathbb{N}}\subset D$ and $D$ is bounded, (9), (10), and
(11) yield

$\lim_{k_{j}\to\infty}\tau_{k_{j}}\left[x^{k_{j}}+t\left(q_{x^{k_{j}}}-x^{k_{j}}\right)\right]=\lim_{k_{j}\to\infty}\left[\tau_{k_{j}}x^{k_{j}}+t\tau_{k_{j}}\left(q_{x^{k_{j}}}-x^{k_{j}}\right)\right]=tv\in
C_{\infty}.$ (12)

As $q_{x^{k_{j}}}\in{\rm argmin}_{p\in{C}}\nabla
f(x^{k_{j}})^{T}(p-x^{k_{j}})$ and $x^{k_{j}}\in C$, we have $\nabla
f(x^{k_{j}})^{T}(q_{x^{k_{j}}}-x^{k_{j}})\leq 0$ for all $j=0,1,\ldots$. Then,
as $(t_{k_{j}})_{j\in\mathbb{N}}\subset(0,\infty)$, we have

$\nabla f(x^{k_{j}})^{T}\left(\tau_{k_{j}}(q_{x^{k_{j}}}-x^{k_{j}})\right)\leq
0,\qquad j=0,1,\ldots.$

Considering that $D$ is bounded and $(x^{k_{j}})_{j\in\mathbb{N}}\subset D$,
we can assume without loss of generality that
$\lim_{k_{j}\to\infty}x^{k_{j}}={\bar{x}}.$ Thus, taking the limit in the last
inequality as $j$ tends to $\infty$ and using (11), we obtain $\nabla
f({\bar{x}})^{T}v\leq 0$, which is absurd because, by (12), we have $v\in
C_{\infty}$, and by assumption (B), we have $\nabla
f({\bar{x}})\in\mbox{int}(C_{\infty})^{*}$. Therefore, the proof of (iii) is
complete. ∎

In the next lemma, we establish that the sequence $(x^{k})_{k\in\mathbb{N}}$
generated by Algorithm 1 is well defined. We also obtain results related to
the optimal value $v_{k}^{*}$ defined in (3).

######  Lemma 3.3.

The sequences $\\{p^{k}\\}_{k\in\mathbb{N}}\subset C$ and
$\\{x^{k}\\}_{k\in\mathbb{N}}\subset C$ are well defined. Moreover, the
following assertions hold:

(i)

$v_{k}^{*}\leq 0$ for all $k=0,1,\ldots$;

(ii)

$v_{k}^{*}=0$ if and only if $x^{k}$ is a stationary point of problem (1);

(iii)

$v_{k}^{*}<0$ if and only if $\lambda_{k}>0$ and $p^{k}\neq x^{k}$.

###### Proof.

For each $x^{k}\in C$, it follows from Proposition 3.2 that $p^{k}$ and
$v_{k}^{*}$ in (3) can be computed and $p^{k}\in C$. As $x^{0}\in C$ and
$0\leq\lambda\leq 1$, by using (5) and an inductive argument, we conclude that
$(p^{k})_{k\in\mathbb{N}}$, $(v_{k}^{*})_{k\in\mathbb{N}}$, and
$(x^{k})_{k\in\mathbb{N}}$ are well defined, and that
$(p^{k})_{k\in\mathbb{N}}$ and $(x^{k})_{k\in\mathbb{N}}$ belong to $C$. To
prove (i), it suffices to note that the optimality in (3) implies
$v_{k}^{*}\leq\nabla f(x^{k})^{T}(x^{k}-x^{k})=0$ for all $k=0,1,\ldots$. To
prove(ii), we note that (3) implies that $v_{k}^{*}\leq\nabla
f(x^{k})^{T}(p-x^{k})$ for all $p\in C$. Thus, if $v_{k}^{*}=0$, we conclude
that $0\leq\nabla f(x^{k})^{T}(p-x^{k})$ for all $p\in C$. Hence, $x^{k}$
satisfies (2), that is, $x^{k}$ is a stationary point of problem (1).
Conversely, if $x^{k}$ is a stationary point of problem (1), then (2) implies
$0\leq\nabla f(x^{k})^{T}(p-x^{k})$ for all $p\in C$. As $p^{k}\in C$, we
conclude that $0\leq\nabla f(x^{k})^{T}(p^{k}-x^{k})=v_{k}^{*}$. Therefore, by
(i), we have $v_{k}^{*}=0$. We now prove (iii). It is immediate from (i) and
(4) that $v_{k}^{*}<0$ if and only if $\lambda_{k}>0$ and $p^{k}\neq x^{k}$;
this concludes the proof. ∎

It follows from Lemma 3.3 that Algorithm 1 generates either an infinite
sequence or a finite sequence $\\{x^{k}\\}_{k\in\mathbb{N}}\subset C$, the
last iterate of which is a stationary point of problem (1). Henceforth, let
$\\{p^{k}\\}_{k\in\mathbb{N}}\subset C$ and
$\\{x^{k}\\}_{k\in\mathbb{N}}\subset C$ be sequences generated by Algorithm 1;
we assume that these sequences are infinite.

###  3.2 Asymptotic convergence analysis

Herein, we study the asymptotic convergence of Algorithm 1. We first prove an
important inequality.

######  Lemma 3.4.

The following inequality holds:

$f(x^{k+1})\leq f(x^{k})-\frac{1}{2}|v_{k}^{*}|\lambda_{k},\qquad
k=0,1,\ldots.$ (13)

Consequently, $f(x^{k})>f(x^{k+1})$ for all $k=0,1,\ldots$.

###### Proof.

Let $x^{k}\in{C}$ be defined as in Algorithm 1, and $v_{k}^{*}$ as in (3). We
first recall that we have assumed that $(x^{k})_{k\in\mathbb{N}}$ is infinite.
Thus, Lemma 3.3 implies that $v_{k}^{*}<0$ and $p^{k}\neq x^{k}$. Applying
Lemma 3.1 with $x=x^{k}$, $p=p^{k}$, and $\lambda=\lambda_{k}$, we have

$f(x^{k+1})\leq
f(x^{k})+v_{k}^{*}\lambda_{k}+\frac{L}{2}\|p^{k}-x^{k}\|^{2}\lambda_{k}^{2}.$
(14)

We separately consider two cases:
$\lambda_{k}=|v_{k}^{*}|/{(L\|p^{k}-x^{k}\|^{2})}$ and $\lambda_{k}=1$. In the
former, it follows from (14) that

$f(x^{k+1})\leq f(x^{k})-\frac{1}{2}|v_{k}^{*}|\lambda_{k}.$ (15)

If now $\lambda_{k}=1$, then (14) becomes $f(x^{k+1})\leq
f(x^{k})-|v_{k}^{*}|+L\|p^{k}-x^{k}\|^{2}/2$, and (4) yields
$\lambda_{k}=1\leq|v_{k}^{*}|/(L\|p^{k}-x^{k}\|^{2})$. Thus, we obtain
$f(x^{k+1})\leq f(x^{k})-(|v_{k}^{*}|/2)\lambda_{k}.$ Therefore, combining
this inequality with (15) yields (13). As we assumed that $v_{k}^{*}<0$ for
all $k=0,1,\ldots$, the second part follows, and the proof is complete. ∎

The next result shows a partial asymptotic-convergence property of the
Frank–Wolfe algorithm; it requires neither convexity nor strong convexity on
$f$.

######  Theorem 3.5.

Each limit point ${\bar{x}}\in C$ of $(x^{k})_{k\in\mathbb{N}}$ is stationary
for the problem (1).

###### Proof.

Let ${\bar{x}}\in C$ be a limit point of the sequence
$(x^{k})_{k\in\mathbb{N}}$, and let $(x^{k_{j}})_{j\in\mathbb{N}}$ be a
subsequence of $(x^{k})_{k\in\mathbb{N}}$ such that
$\lim_{j\to\infty}x^{k_{j}}={\bar{x}}$. Hence,
$\lim_{j\to\infty}f(x^{k_{j}})=f({\bar{x}})$. As Lemma 3.4 implies that
$(f(x^{k}))_{k\in\mathbb{N}}$ is a decreasing sequence, we conclude that
$(f(x^{k}))_{k\in\mathbb{N}}$ converges to $f({\bar{x}})$. Thus, in
particular,

$\lim_{k\to\infty}[f(x^{k})-f(x^{k+1})]=0.$ (16)

Moreover, as $(x^{k_{j}})_{j\in\mathbb{N}}$ is bounded, by combining the
inclusion in (3) with (iii) of Proposition 3.2, it follows that
$(p^{k_{j}})_{j\in\mathbb{N}}$ is also bounded. Let
$(p^{k_{\ell}})_{\ell\in\mathbb{N}}$ be a subsequence of
$(p^{k_{j}})_{j\in\mathbb{N}}$ such that
$\lim_{\ell\to\infty}p^{k_{\ell}}={\bar{p}}$. If ${\bar{x}}={\bar{p}}$, then
by (3) and the continuity of $\nabla f$, we have
$\lim_{\ell\to\infty}v_{k_{\ell}}^{*}=\nabla
f(x^{k_{\ell}})^{T}(p^{k_{\ell}}-x^{k_{\ell}})=0$. We now assume that
${\bar{x}}\neq{\bar{p}}$. Hence, combining Lemma 3.4 with (16), we conclude
that $\lim_{\ell\to\infty}|v_{k_{\ell}}^{*}|\lambda_{k_{\ell}}=0$, whereas (4)
yields

$|v_{k_{\ell}}^{*}|\lambda_{k_{\ell}}=\mbox{min}\left\\{|v_{k_{\ell}}^{*}|,\frac{|v_{k_{\ell}}^{*}|^{2}}{L\|p^{k_{\ell}}-x^{k_{\ell}}\|^{2}}\right\\}.$
(17)

As $\lim_{\ell\to\infty}x^{k_{\ell}}={\bar{x}}$,
$\lim_{\ell\to\infty}p^{k_{\ell}}={\bar{p}}$, and ${\bar{x}}\neq{\bar{p}}$, we
obtain that
$\lim_{\ell\to\infty}\|p^{k_{\ell}}-x^{k_{\ell}}\|=\|{\bar{x}}-{\bar{p}}\|\neq
0$. Thus, as $\lim_{\ell\to\infty}|v_{k_{\ell}}^{*}|\lambda_{k_{\ell}}=0$, it
follows from (17) that $\lim_{\ell\to\infty}|v_{k_{\ell}}^{*}|=0$. Moreover,
the optimality of $v_{k_{\ell}}^{*}$ in (3) implies

$v_{k_{\ell}}^{*}\leq\nabla
f(x^{k_{\ell}})^{{\scriptscriptstyle\mathrm{T}}}(p-x^{k_{\ell}}),\qquad\forall~{}p\in
C.$ (18)

As $\lim_{\ell\to\infty}v_{k_{\ell}}^{*}=0$, taking the limit in (18) and
using the continuity of $\nabla f$, we have $\nabla
f(\bar{x})^{{\scriptscriptstyle\mathrm{T}}}(p-\bar{x})\geq 0$ for all $p\in
C$. Therefore, $\bar{x}$ is stationary for the problem (1). ∎

We now show that if $f$ is assumed convex, we can improve the previous result.

######  Theorem 3.6.

The following assertions hold:

(i)

If $f$ is a convex function and $x^{*}$ is a cluster point of the sequence
$(x^{k})_{k\in\mathbb{N}}$, then $x^{*}$ is a solution for the problem (1).

(ii)

If $f$ is an $M$-strongly convex function, then $(x^{k})_{k\in\mathbb{N}}$
converges to a point $x^{*}\in C$ that is a solution of problem (1). Moreover,
$\|x^{k}-x^{*}\|\leq\sqrt{2(f(x^{k})-f(x^{*}))/M}$ for all $k=0,1,\ldots$.

###### Proof.

To prove(i), we assume that $f$ is convex and $x^{*}$ is a cluster point of
$(x^{k})_{k\in\mathbb{N}}$. As $f$ is convex, applying Proposition 2.2 with
$M=0$, we obtain that $f(p)\geq f(x^{*})+\nabla f(x^{*})^{T}(p-x^{*})$ for all
$p\in C$. Therefore, considering that Theorem 3.5 implies that $\nabla
f(x^{*})^{{\scriptscriptstyle\mathrm{T}}}(p-x^{*})\geq 0$ for all $p\in C$, we
conclude that $f(p)\geq f(x^{*})$ for all $p\in C$. Then, $x^{*}$ is a
solution for the problem (1). To prove (ii), we first note that as $f$ is
$M$-strongly convex, the level set ${\cal L}_{f(x^{0})}:=\\{x\in C:~{}f(x)\leq
f(x^{0})\\}$ is bounded. Lemma 3.4 now implies that
$(x^{k})_{k\in\mathbb{N}}\subset{\cal L}_{f(x^{0})}$. Hence,
$(x^{k})_{k\in\mathbb{N}}$ is also bounded. Let $x^{*}\in C$ be a cluster
point of $(x^{k})_{k\in\mathbb{N}}$. It follows from (i) that $x^{*}$ is a
solution of problem (1). Furthermore, combining (2) with Proposition 2.2, we
obtain

$f(x^{k})-f(x^{*})\geq\frac{M}{2}\|x^{*}-x^{k}\|^{2},\qquad k=0,1,\ldots.$
(19)

As Lemma 3.4 implies that $(f(x^{k}))_{k\in\mathbb{N}}$ is a decreasing
sequence, we have $\lim_{k\to+\infty}f(x^{k})=f(x^{*})$. Therefore, (19)
implies that $(x^{k})_{k\in\mathbb{N}}$ converges to $x^{*}$. Finally, we note
that (19) is equivalent to the inequality in (ii), and the proof is complete.
∎

###  3.3 Iteration-complexity analysis

Herein, we derive two iteration-complexity bounds for the sequence
$(x^{k})_{k\in\mathbb{N}}$ generated by Algorithm 1. To this end, we assume
that $\lim_{k\to\infty}x_{k}=x^{*}$. Then, (iii) of Proposition 3.2 implies
that $(p^{k})_{k\in\mathbb{N}}$ is bounded. Therefore, we define

$0<\sigma:=\sup_{k}\\{\|p^{k}-x^{k}\|:~{}k=0,1,\ldots\\}<\infty.$ (20)

We also define the following constants:

$\Gamma:=\min\left\\{\frac{1}{2\gamma\sigma},\frac{1}{2L\sigma^{2}}\right\\}>0,\qquad\gamma:=\max\left\\{\|\nabla
f(x^{k})\|:~{}k=0,1,\ldots\right\\}>0.$ (21)

######  Theorem 3.7.

The following assertions hold:

(i)

If $f$ is a convex function, then $f(x^{k})-f^{*}\leq\Gamma^{-1}/k$ for all
$k=1,2,\ldots$.

(ii)

If $f$ is $M$-strongly convex, then $\|x^{k}-x^{*}\|\leq\sqrt{2/(\Gamma
M)}/\sqrt{k}$ for all $k=1,2,\ldots$.

###### Proof.

To prove (i), we first prove the following inequality:

$\Gamma{v_{k}^{*}}^{2}\leq f(x^{k})-f(x^{k+1}),\qquad\forall~{}k=0,1,\dots,$
(22)

where $\Gamma$ is defined in (21). By using Lemma 3.4, with $\sigma$ defined
in(20), and considering that ${v_{k}^{*}}<0$, we conclude after some algebraic
manipulations that

$\min\left\\{\frac{1}{2|v_{k}^{*}|},\frac{1}{2L\sigma^{2}}\right\\}{v_{k}^{*}}^{2}\leq
f(x^{k})-f(x^{k+1}),\qquad\forall~{}k=0,1,...$ (23)

Further, by combining (3) with (20) and the second equality in (21), we obtain
$0<|v_{k}^{*}|\leq\|\nabla f(x^{k})\|\|x^{k}-p^{k}\|\leq\gamma\sigma$ for all
$k=0,1,\ldots$, which implies

$\frac{1}{\gamma\sigma}\leq\frac{1}{|v_{k}^{*}|},\qquad k=0,1,\ldots.$

Thus, (22) follows from (23), the previous inequality, and (21). It now
follows from (ii) of Theorem 3.6 that $x^{*}$ is a solution of problem (1). As
$f$ is convex, we have $f^{*}=f(x^{*})\geq f(x^{k})+\nabla
f(x^{k})^{T}(x^{*}-x^{k})$ for all $k$. Thus, (3) implies that
$f^{*}-f(x^{k})\geq\nabla f(x^{k})^{T}(x^{*}-x^{k})\geq
f(x^{k})^{T}(p^{k}-x^{k})=v_{k}^{*}$ for all $k$. As $f^{*}\leq f(x^{k})$ for
all $k$, we conclude that $v_{k}^{*}\leq f^{*}-f(x^{k})\leq 0$ for all $k$.
Therefore, we obtain $(f(x^{k})-f^{*})^{2}\leq v_{k}^{*2}$ for all
$k=0,1,\ldots$; this, combined with (22), yields

$(f(x^{k})-f^{*})-(f(x^{k+1})-f^{*})\geq\Gamma(f(x^{k})-f^{*})^{2},\qquad
k=0,1,\ldots.$ (24)

As $\Gamma>0$, if we define $a_{k}:=f(x^{k})-f^{*}$, we conclude from (24)
that $a_{k}-a_{k+1}\geq\Gamma a_{k}^{2}$. Thus, (i) follows by applying Lemma
2.3. To prove (ii), we first note that the $M$-strong convexity of $f$ and (2)
imply that $\|x^{k}-x^{*}\|\leq\sqrt{2(f(x^{k})-f(x^{*}))/M}$ for all
$k=0,1,\ldots$. Therefore, (ii) follows by using the inequality in (i), and
the proof is complete. ∎

##  4 Examples

Herein, we present examples of problem (1), with $f$ satisfying (A) and (B),
and $C$ unbounded. We first present a general class of functions $f$
satisfying (A) and (B). To this end, let ${\cal{K}}\subset\mathbb{R}^{n}$ be a
closed convex cone such that
$\mbox{int}{\cal{K}}\cap\mbox{int}{\cal{K}}^{*}\neq\varnothing$,
$G:\mathbb{R}^{n}\to\mathbb{R}^{n}$ be a differentiable function, and
$G^{\prime}$ be its Jacobian. We assume that for constants $L_{1}>0$ and
$L_{2}>0$, the function $G$ satisfies the following conditions:

(C1)

$\|G(x)-G(y)\|\leq L_{1}\|x-y\|$ for all $x,y\in{\cal{K}}$;

(C2)

$\|G^{\prime}(x)x-G^{\prime}(y)y\|\leq L_{2}\|x-y\|$ for all
$x,y\in{\cal{K}}$;

(C3)

$G(x)+G^{\prime}(x)x\in{\cal{K}}^{*}$ for all $x\in{\cal{K}}$.

Let $a\in\mbox{int}{\cal{K}}^{*}$, and a function
$f:\mathbb{R}^{n}\to\mathbb{R}$ be defined by

$f(x):=a^{T}x+G(x)^{T}x.$ (25)

######  Lemma 4.1.

Let $C\subset{\cal{K}}$ be closed and convex. Then, the gradient $\nabla f$ of
$f$ defined in (25) satisfies conditions (A) and (B) on $C$.

###### Proof.

We first note that the gradient $\nabla f$ of the function $f$ defined in (25)
is given by

$\nabla f(x)=a+G(x)+G^{\prime}(x)x.$ (26)

Thus, using (C1) and (C2), we conclude that $\|\nabla f(x)-\nabla
f(y)\|\leq(L_{1}+L_{2})\|x-y\|$ for all $x,y\in C$. Hence, $\nabla f$
satisfies (A). Finally, as $a\in\mbox{int}{\cal{K}}^{*}$, (C3) implies that
$\nabla f(x)\in\mbox{int}{\cal{K}}^{*}$ for all $x\in C$. Furthermore, as
$C\subset{\cal{K}}$, we obtain that $C_{\infty}\subset{\cal{K}}$. Hence,
${\cal{K}}^{*}\subset C_{\infty}^{*}$, which implies that
$\mbox{int}{\cal{K}}^{*}\subset\mbox{int}(C_{\infty}^{*})$. Thus, we conclude
that $\nabla f(x)\in\mbox{int}(C_{\infty})^{*}$ for all $x\in C$. Therefore,
$\nabla f$ satisfies (B) on $C$. ∎

We now recall a well-known result about Lipschitz functions, which is an
immediate consequence of the mean-value inequality.

######  Theorem 4.2.

Let $C\subset{\mathbb{R}}^{n}$ be a convex set,
$F:{\mathbb{R}}^{n}\to{\mathbb{R}}^{m}$ be a continuously differentiable
function, and $F^{\prime}$ be its Jacobian. We assume that there exists a
constant $L\geq 0$ such that $\|F^{\prime}(x)\|\leq L$ for all $x\in C$. Then,
$F$ is Lipschitz continuous with constant $L$ on $C$, that is,
$\|F(x)-F(y)\|\leq L\|x-y\|$ for all $x,y\in C$.

Moreover, the following characterization of convex functions is required; its
proof can be found in [16, Theorem 4.1.1, pp. 190].

######  Theorem 4.3.

Let $C\subset{\mathbb{R}}^{n}$ be a convex set,
$f:{\mathbb{R}}^{n}\to{\mathbb{R}}$ be a twice continuously differentiable
function, $\nabla^{2}f$ be its Hessian, and $M\geq 0$ be a constant. Then, $f$
is $M$-strongly convex on $C$ if and only if $v^{T}\nabla^{2}f(x)v\geq
M\|v\|^{2}$ for all $x\in C$ and $v\in{\mathbb{R}}^{n}$.

In the following, we present specific examples of functions $G$ satisfying
(C1)–(C3) on the cone ${\cal{K}}=\mathbb{R}^{n}_{+}$ such $f$ in (25) is
convex and $\nabla f$ satisfies (A) and (B) for any closed convex set
$C\subset\mathbb{R}^{n}_{+}$.

######  Example 4.4.

Let $Q=(q_{ij})\in{\mathbb{R}}^{n\times n}$ with entries $q_{ij}\geq 0$ for
all $i,j$ and $a\in\mathbb{R}^{n}_{++}$. Let
$G:\mathbb{R}^{n}\to\mathbb{R}^{n}$ be a linear function defined by $G(x)=Qx$.
Direct calculations show that
$\|G(x)-G(y)\|=\|G^{\prime}(x)x-G^{\prime}(y)y\|\leq\|Q\|\|x-y\|$ for all
$x,y\in\mathbb{R}^{n}$. Thus, $G$ satisfies (C1) and (C2), with
$L_{1}=L_{2}=\|Q\|$, for any cone ${\cal{K}}$. As $q_{ij}\geq 0$ and
$a\in\mathbb{R}^{n}_{++}$, $G$ also satisfies (C3) in
${\cal{K}}=\mathbb{R}^{n}_{+}$. Moreover, if $Q$ satisfies the condition
$v^{T}Qv\geq M\|v\|^{2}$ for all $v\in{\mathbb{R}}^{n}$ and some $M\geq 0$,
then Theorem 4.3 implies that the problem (1) with the associated quadratic
function $f(x):=a^{T}x+x^{T}Qx$ is $M$-strong convex. Moreover, Lemma 4.1
implies that $\nabla f(x)=a+Qx$ satisfies (A) and (B) for any closed and
convex set $C\subset\mathbb{R}^{n}_{+}$. For instance,
$C:=\left\\{x\in\mathbb{R}^{p}:~{}x_{1}\geq x_{2}\geq\cdots\geq x_{p}\geq
0\right\\}$, the monotone nonnegative cone.

######  Example 4.5.

Let $e:=(1,\ldots,n)\in\mathbb{R}^{n}$ be a vector, and $\alpha>0$ and
$\beta>0$ be constants satisfying $2\alpha>3\beta^{3/2}\sqrt{n}$. Then, the
function $G:\mathbb{R}^{n}\to\mathbb{R}^{n}$ defined by

$G(x):=\alpha x+\frac{\beta}{\sqrt{1+\beta x^{T}x}}e,$ (27)

satisfies (C1)–(C3) in ${\cal{K}}=\mathbb{R}^{n}_{+}$. Moreover, the problem
(1) with the associated function

$f(x):=a^{T}x+\alpha x^{T}x+\frac{\beta}{\sqrt{1+\beta x^{T}x}}e^{T}x,$ (28)

with $a\in\mathbb{R}^{n}_{++}$ is $M$-strong convex. Moreover, $\nabla f$
satisfies (A) and (B). Indeed, we first note that $f(x):=a^{T}x+G(x)^{T}x$.
Some calculations show that

$G^{\prime}(x)=\alpha{\rm I_{n}}+\frac{-\beta^{2}}{(1+\beta
x^{T}x)^{3/2}}ex^{T},\qquad G(x)+G^{\prime}(x)x=2\alpha
x+\frac{\beta}{(1+\beta x^{T}x)^{3/2}}e,$ (29)

where ${\rm I_{n}}\in{\mathbb{R}}^{n\times n}$ is the identity matrix. Thus,
$\nabla f(x)=a+G(x)+G^{\prime}(x)x$. Hence, after some calculations, we have

$\nabla^{2}f(x)=2\alpha{\rm I_{n}}+\frac{-3\beta^{2}}{(1+\beta
x^{T}x)^{5/2}}ex^{T}.$ (30)

The first equality in (29) yields
$\|G^{\prime}(x)\|\leq\beta^{3/2}\sqrt{n}+\alpha$ for all $x\in C$, and (30)
implies

$0<\left(2\alpha-3\beta^{3/2}\sqrt{n}\right)|v\|^{2}\leq
v^{T}\nabla^{2}f(x)v\leq\left(2\alpha+3\beta^{3/2}\sqrt{n}\right)\|v\|^{2},$
(31)

for all $x\in C$ and all $v\in{\mathbb{R}}^{n}$. As
$\|G^{\prime}(x)\|\leq\beta^{3/2}\sqrt{n}+\alpha$, it follows from Theorem 4.2
that $G$ also satisfies (C1) with $L_{1}=\beta^{3/2}\sqrt{n}+\alpha$. In
particular, (31) implies that $\|\nabla^{2}f(x)\|\leq
2\alpha+3\beta^{3/2}\sqrt{n}$, for all $x\in C$. Hence, as
$\|G^{\prime}(x)\|\leq\beta^{3/2}\sqrt{n}+\alpha$ for all $x\in C$, Theorem
4.2 also implies that

$\|G^{\prime}(x)x-G^{\prime}(y)y\|\leq\|G(x)-G(y)\|+\|\nabla f(x)-\nabla
f(y)\|\leq\left(4\beta^{3/2}\sqrt{n}+3\alpha\right)\|x-y\|,$

for all $x,y\in C$. Thus, $G$ also satisfies (C2) with
$L_{2}=4\beta^{3/2}\sqrt{n}+3\alpha$. The second inequality in (29) implies
that $G(x)+G^{\prime}(x)x\in{\mathbb{R}^{n}_{++}}$ for all
$x\in\mathbb{R}^{n}_{+}$. As $(\mathbb{R}^{n}_{+})^{*}=\mathbb{R}^{n}_{+}$,
$G$ satisfies (C3). Therefore, Lemma 4.1 implies that $\nabla f$ of $f$ in
(28) satisfies conditions (A) and (B) for any closed and convex set
$C\subset\mathbb{R}^{n}_{+}$. For instance,
$C:=\left\\{x\in\mathbb{R}^{p}:~{}x_{1}\geq x_{2}\geq\cdots\geq x_{p}\geq
0\right\\}$, the monotone nonnegative cone. Finally, using (31), it follows
from Theorem 4.3 that $f$ in (28) is $M$-strong convex with
$M=2\alpha-3\beta^{3/2}\sqrt{n}$.

In the next example, we present directly a convex function $f$ satisfying (A)
and (B) in a closed convex set $C\subset\mathbb{R}^{n}_{+}$

######  Example 4.6.

Let $\beta>0$, $a\in\mathbb{R}^{n}_{++}$, and $f:\mathbb{R}^{n}\to\mathbb{R}$
be defined by

$f(x):=a^{T}x+\sqrt{1+\beta x^{T}x}.$ (32)

Let $C:=\left\\{x\in\mathbb{R}^{p}:~{}x_{1}\geq x_{2}\geq\cdots\geq x_{p}\geq
0\right\\}$ be the monotone nonnegative cone. We note that, in this case, the
gradient and the Hessian of $f$ are given by

$\nabla f(x)=a+\frac{\beta}{\sqrt{1+\beta
x^{T}x}}x,\qquad\nabla^{2}f(x)=\frac{\beta}{\sqrt{1+\beta x^{T}x}}{\rm
I_{n}}-\frac{\beta^{2}}{(1+\beta x^{T}x)^{3/2}}xx^{T},$

respectively, where ${\rm I_{n}}\in{\mathbb{R}}^{n\times n}$ is the identity
matrix. Some calculations show that

$\frac{\beta}{(1+\beta x^{T}x)^{3/2}}v^{T}v\leq v^{T}\nabla^{2}f(x)v\leq\beta
v^{T}v,\qquad\forall~{}v\in\mathbb{R}^{n},$

which implies that $\nabla^{2}f(x)$ is positive definite and
$\|\nabla^{2}f(x)\|\leq\beta$. Thus, using Theorems 4.3 and 4.2, we conclude
that $f$ is convex, and $\nabla f$ is Lipschitz continuous with constant
$\beta$. Moreover, for any closed and convex set $C\subset\mathbb{R}^{n}_{+}$,
we have $C_{\infty}\subset\mathbb{R}^{n}_{+}$. Hence,
$\mathbb{R}^{n}_{+}=(\mathbb{R}^{n}_{+})^{*}\subset C_{\infty}^{*}$, which
implies that $\mathbb{R}^{n}_{++}\subset int(C_{\infty}^{*})$. As $\nabla
f(x)\in\mathbb{R}^{n}_{++}$ for all $x\in C$, we conclude that $\nabla f(x)\in
int(C_{\infty}^{*})$ for all $x\in C$. Finally, the convexity of $g$ implies
that $C$ is convex. Therefore, the problem (1) with the objective function
(32) is convex, and $f$ satisfies conditions (B) and (A).

Let us present two more examples of convex functions and the respective convex
sets satisfying conditions (B) and (A).

(i)

Let $f:\mathbb{R}^{n}\to\mathbb{R}$ be defined by
$f(x)=\ln(e^{x_{1}}+e^{x_{2}}+\ldots+e^{x_{n})}$. This function satisfies (A)
[19, Example 5.15, pp. 115]. Some calculations show that it also satisfies (B)
for any convex set $C\subset\mathbb{R}^{n}_{+}$.

(ii)

Let $C\subset\mathbb{R}^{n}_{++}$ be closed and convex, and
$d_{C}(x):=\min_{y\in C}\|x-y\|$ for $x\in\mathbb{R}^{n}$. We define the
convex function $\psi_{C}:\mathbb{R}^{n}\to\mathbb{R}$ by
$\psi_{C}(x):=\frac{1}{2}\|x\|^{2}-\frac{1}{2}d^{2}_{C}(x)$ [19, Example
2.17.4, pp. 22]. We can prove that $\nabla\psi_{C}(x)=P_{C}(x)$, where $P_{C}$
denotes the orthogonal projection onto $C$ [19, Example 3.49, pp. 61]. As
$\mathbb{R}^{n}_{++}\subset int(C_{\infty}^{*})$ and
$C\subset\mathbb{R}^{n}_{++}$, we conclude that $\nabla\psi_{C}(x)\in
C_{\infty}^{*}$, which implies that $\psi_{C}$ satisfies (B). Moreover, the
nonexpansivity of the projection implies that $\psi_{C}$ also satisfies (A).

We conclude this section by presenting examples of unbounded convex sets that
appear as constraints in optimization problems.

(a)

$C:=\\{x\in\mathbb{R}^{n}_{+}:~{}1\leq x_{1}\ldots x_{n}\\}$;

(b)

$C:=\\{x\in\mathbb{R}^{n}_{+}:~{}1\leq x_{1}+\ldots+x_{n}\\}$;

(c)

$C:=\\{x\in\mathbb{R}^{n}_{+}:~{}b\leq Ax\\}$, where
$A=(a_{ij})\in{\mathbb{R}}^{n\times n}$ with $a_{ij}>0$ and
$b\in\mathbb{R}^{n}_{++}$.

Let $\Omega\subset\mathbb{R}^{n}$ be a closed, convex set, and
$g:\Omega\to\mathbb{R}$ be a convex function. The epigraph of $g$ is defined
by $\mbox{epi}(g):=\\{(x,t)\in\Omega\times\mathbb{R}:~{}g(x)\leq t\\}.$ The
set $C=\mbox{epi}(g)\subset{\mathbb{R}}^{n}\times{\mathbb{R}}$ is convex and
unbounded. We now provide specific examples.

(1)

(Lorentz cone)
$\left\\{(x,t)\in\mathbb{R}^{n}\times\mathbb{R}:~{}\|x\|_{2}\leq t\right\\}$,
where $\|\cdot\|_{2}$ denotes the $2$-norm;

(2)

$\left\\{(x,t)\in\mathbb{R}^{n}\times\mathbb{R}:~{}\|x\|_{1}\leq t\right\\}$,
where $\|\cdot\|_{1}$ denotes the $1$-norm;

(3)

$\left\\{(x_{1},\ldots,x_{n},t)\in\
\mathbb{R}^{n}_{++}\times\mathbb{R}:~{}1/(x_{1}\ldots x_{n})\leq t\right\\}$.

We point out that projecting on the sets $(1)$, $(2)$, and $(3)$ is not a
particularly expensive task [19, Chapter 6]. Finally, we present an example on
the cone of positive semidefinite matrices.

######  Example 4.7.

The cone of positive semidefinite (resp., definite) $n\times n$ symmetric
matrices is denoted by ${\mathbb{S}}^{n}_{+}$ (resp., ${\mathbb{S}}^{n}_{++}$)
and is selfdual, that is, ${{\mathbb{S}}^{n}_{+}}^{*}={\mathbb{S}}^{n}_{+}$.
The trace of $X=(X_{ij})\in{\mathbb{S}}^{n}$ is denoted by
$trX:=\sum_{i=1}^{p}X_{ii}$. Given $X$ and $Y$ in ${\mathbb{S}}^{n}$, their
inner product is defined as $\langle
X,Y\rangle:=trXY=\sum_{i=1,j=1}^{n,m}X_{ij}Y_{ij}$, whereas the norm of $X$ is
defined by $\|X\|:=\langle X,X\rangle^{1/2}$. Let $\beta>0$ and $\gamma>0$ be
such that $2\beta>\gamma$, and let $f:{\mathbb{S}}^{n}_{+}\to\mathbb{R}$ be
defined by $f(X):=\beta\|X||^{2}+\gamma\ln(1+trX).$ We consider the set
$C:=\\{X\in{\mathbb{S}}^{n}_{+}:~{}g(X)\leq 0\\}$, where
$g:{\mathbb{S}}^{n}_{+}\to\mathbb{R}$ is a convex function. The gradient of
$f$ is given by

$\nabla f(X)=2\beta X+\frac{1}{1+trX}{\rm
I_{n}}\in{\mathbb{S}}^{n}_{++},\qquad\forall~{}X\in{\mathbb{S}}^{n}_{+}.$

Calculating the Hessian of $F$, we obtain

$\langle\nabla^{2}f(X)V,V\rangle=2\beta\|V\|^{2}-\frac{\gamma}{\left(1+trX\right)^{2}}(trV)^{2},\quad\forall~{}X\in{\mathbb{S}}^{n}_{+},\quad\forall~{}V\in{\mathbb{S}}^{n}.$

Moreover, $0<(2\beta-\gamma)\|V\|^{2}\langle\nabla f(X)V,V\rangle\leq
2\beta\|V\|^{2}$ for all $X\in{\mathbb{S}}^{n}_{+}$ and $0\neq
V\in{\mathbb{S}}^{n}$. Thus, Theorems 4.2 and 4.3 imply that $\nabla f$ is
Lipschitz continuous with constant $2\beta$, and $f$ is strongly convex with
constant $2\beta-\gamma$, respectively. Furthermore, as
$C\subset{\mathbb{S}}^{n}_{+}$, we conclude that
$C_{\infty}\subset{\mathbb{S}}^{n}_{+}$. Hence,
${\mathbb{S}}^{n}_{+}={{\mathbb{S}}^{n}_{+}}^{*}\subset C_{\infty}^{*}$, which
implies that ${\mathbb{S}}^{n}_{++}\subset int(C_{\infty}^{*})$. As $\nabla
f(x)\in{\mathbb{S}}^{n}_{++}$ for all $x\in C$, we conclude that $\nabla
f(x)\in int(C_{\infty}^{*})$ for all $x\in C$. Moreover the convexity of $g$
implies that $C$ is convex. Therefore, the problem (1) is convex, and $f$
satisfies conditions (A) and (B). In particular, if
$g:{\mathbb{S}}^{n}_{++}\to\mathbb{R}$ given by $g(x)=1/det(X)$ or
$g(x)=1/tr(X)$, then the set $C\subset{\mathbb{S}}^{n}_{++}$ is convex and
unbounded.

##  5 Conclusions

In this study, we considered the classical Frank–Wolf algorithm for nonempty,
closed, convex, not necessarily compact constraints. To study its convergence
properties, we used recession techniques. The examples in Section 4
demonstrated that the Frank–Wolf algorithm can indeed be applied to several
optimization problems (not necessarily convex) with non-compact constraint
sets.

## Funding

The first author was supported in part by CNPq grants 305158/2014-7 and
302473/2017-3, FAPEG/PRONEM- 201710267000532, and CAPES. The second author was
supported in part by Fundação de Apoio à Pesquisa do Distrito Federal (FAP-DF)
through grants 0193.001695/2017 and PDE 05/2018. This research was partly
carried out during a visit of the second author to the Center for Mathematical
Research (CRM) (while Western Catalonia was in a state of alert), in the
framework of the Research-in-pairs call in 2020. The CRM is a paradise for
research, and the second author appreciates its hospitality and support.

## References

  * [1] Frank M, Wolfe P. An algorithm for quadratic programming. Nav Res Log. 1956;:95–110.
  * [2] Demyanov VF, Rubinov AM. Approximate methods in optimization problems. American Elsevier Publishing Co., Inc., New York; 1970.
  * [3] Levitin E, Polyak B. Constrained minimization methods. USSR Computational Mathematics and Mathematical Physics. 1966;6(5):1–50.
  * [4] Jaggi M. Revisiting frank-wolfe: Projection-free sparse convex optimization. Proceedings of the 30th International Conference on International Conference on Machine Learning - Volume 28. 2013;ICML’13:I–427–I–435.
  * [5] Lacoste-Julien S, Jaggi M. On the global linear convergence of frank-wolfe optimization variants. arXiv e-prints. 2015;arXiv:1511.05932.
  * [6] Lan G. The complexity of large-scale convex programming under a linear optimization oracle. arXiv e-prints.
  * [7] Beck A, Teboulle M. A conditional gradient method with linear rate of convergence for solving convex linear systems. Math Methods Oper Res. 2004;59(2):235–247.
  * [8] Harchaoui Z, Juditsky A, Nemirovski A. Conditional gradient algorithms for norm-regularized smooth. convex optimization. Math Program. 2015;152(1-2, Ser. A):75–112.
  * [9] Boyd N, Schiebinger G, Recht B. The alternating descent conditional gradient method for sparse inverse problems. SIAM J Optim. 2017;27(2):616–639.
  * [10] Luss R, Teboulle M. Conditional gradient algorithms for rank-one matrix approximations with a sparsity constraint. SIAM Rev. 2013;55(1):65–98.
  * [11] Freund RM, Grigas P, Mazumder R. An extended Frank-Wolfe method with “in-face” directions, and its application to low-rank matrix completion. SIAM J Optim. 2017;27(1):319–346.
  * [12] Konnov IV. Simplified versions of the conditional gradient method. Optimization. 2018;67(12):2275–2290.
  * [13] Ghadimi S. Conditional gradient type methods for composite nonlinear and stochastic optimization. Math Program. 2019;173(1-2, Ser. A):431–464.
  * [14] Lan G, Zhou Y. Conditional gradient sliding for convex optimization. SIAM J Optim. 2016;26(2):1379–1409.
  * [15] Hiriart-Urruty JB, Lemaréchal C. Fundamentals of convex analysis. Springer-Verlag, Berlin; 2001. Grundlehren Text Editions; abridged version of ıt Convex analysis and minimization algorithms. I [Springer, Berlin, 1993\.
  * [16] Hiriart-Urruty JB, Lemaréchal C. Convex analysis and minimization algorithms. I. (Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences]; Vol. 305). Springer-Verlag, Berlin; 1993\.
  * [17] Polyak BT. Introduction to optimization. Optimization Software, New York; 1987. Translations Series in Mathematics and Engineering.
  * [18] Dennis JE Jr, Schnabel RB. Numerical methods for unconstrained optimization and nonlinear equations. (Classics in Applied Mathematics; Vol. 16). Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; 1996.
  * [19] Beck A. First-order methods in optimization. (MOS-SIAM Series on Optimization; Vol. 25). Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA; Mathematical Optimization Society, Philadelphia, PA; 2017.


# Modelling pulsatile blood flow through surgically coupled microvascular
anastomoses.

M. T. Gallagher  School of Mathematics, University of Birmingham,Birmingham.
B15 2TT. UK

###### Abstract

The significant rates of post-operative arterial thrombus formation following
the sutured anastomosis (surgical connection) of vessels in microvascular
reconstruction has kindled the development of engineered solutions in the form
of surgical coupling devices. These devices reduce thrombosis rates, which
recent numerical studies suggest to be due to decreases in shear strain rate
across the anastomosis site when compared to sutured vessels. In this work we
develop and analytically solve the first mathematical model for the problem of
microvascular anastomosis, leading to the discovery that the rates of
thrombosis using these surgical coupling devices could be further reduced by
decreasing the gradient of the wall deformation at the join.

##  1 Introduction

The ability to create a well structured anastomosis, the surgical connection
between two vessels, is essential in a whole range of medical applications,
particularly in the case of microvascular reconstruction. It is well known,
for example, that the survival of tissue post-breast reconstructive surgery is
highly dependent on the ability of the connected vasculature to flow freely
(Koul et al., 2013), with a leading cause of failure being that of thrombosis
(Kroll et al., 1996). This has led to the recent series of papers by Wain et
al. (Wain et al. (2014), Wain et al. (2016), Wain et al. (2018)) who performed
the first detailed numerical investigations into the flow properties through
microvascular anastomoses.

Currently the preferred method for vascular anastomoses involves the end-to-
end suturing of vessels (Carrel, 1902). However, for smaller vessels in the
microvasculature, even when executed by the most highly trained of surgeons,
this technique can result in thrombosis rates of around $3.3\%$ (Yap et al.,
2006) and can be very time consuming. As a result of this, there has
engineered solutions have been developed in the form of coupling devices such
as the UNILINK coupler (Ostrup & Berggren, 1986), see figure 1a, which has
been shown to result in thrombosis rates as low as $0.6\%$ (Jandali et al.,
2010), with several sources reporting the decrease in both the failure rates
and the time required for anastomoses when compared to the use of sutures de
Bruijn & Marck (1996). The UNILINK device was designed for venous anastomoses
and as such is not recommended for arterial anastomoses, although this has not
prevented the developing use of such devices in breast reconstruction (Spector
et al., 2006), and free tissue transfer (Grewal et al., 2012). This push
towards the use of surgical coupling devices on small vessels for
microvascular arterial anastomoses means that it is more important than ever
that we fully understand the detailed properties of how the flow is affected
by this procedure.

The key influence on the rates of thrombosis is the shear rate on the wall of
the vessel, with increases in such having been shown to lead to increases in
thrombus formation (Sakariassen et al., 2015) due to increased platelet
activation and aggregation at the site of anastomosis (Murray, 1926). In
addition to increases in the magnitude of shear rate, it has been shown by
Nesbitt et al. (2009) that shear gradients play a significant role in the rate
and size of platelet aggregation, and thus the potential formation of
thrombus.

With this in mind, recent studies have been investigating in detail the effect
of both surgical coupling and suturing on the flow through small vessels
approximately $1$ mm in radius, with particular attention being paid to the
shear rate across the anastomosis. These studies, performed numerically with
computational fluid dynamics (CFD) packages such as ANSYS-CFX, have aimed to
develop increasingly realistic models of microvascular anastomoses, beginning
by investigating the flow due to a constant pressure gradient (Wain et al.,
2014), followed by detailed studies into the effect of suture placement (Wain
et al., 2016), and, most recently, considering the effects of including flows
driven by pulsatile pressure gradients (Wain et al., 2018). To complement
these studies, we develop here a mathematical model for the problem, which,
once solved analytically, will allow the more detailed analysis of how the
shear rate is influenced by the inclusion of a microvascular anastomosis.

The use of mathematical models has been fundamental to understanding all
aspects of fluid dynamics, not-least those problems relating to the flow of
fluid through arteries. Of particular interest for the present work is the
first analytical expression for the velocity of a viscous incompressible fluid
through an axisymmetric pipe, derived by Womersley (1955). In this study we
will utilise the perturbation techniques of Van Dyke (1964) to consider the
asymptotic correction terms to the Womersley flow profile due to a
microvascular anastomosis, simplifying the problem by noting that we can
consider sections of vessel which have a radial length scale that is much
smaller than the axial length scale. The use of such simplifications, as well
as other mathematical techniques, have been widely used to great success for
related problems in blood flow; for a detailed explanation of the subject see
Pedley (1980). These techniques will enable significant analytical progress to
be made in understanding the changes in shear rate across the site of
anastomosis, and thus will bring new insight into the ways in which rates of
thrombosis may be further reduced.

In this work, we develop the first mathematical model for coupled
microvascular anastomosis in small vessels. By modelling the anastomosis as a
small axisymmetric perturbation to a pipe with rigid walls, and considering an
extruding bump as a model of wall deformation due to a coupling device, and an
intruding bump as a first approximation to the deformation due to sutures, we
investigate in detail the asymptotic correction to the flow profile driven by
a pulsatile pressure gradient. We begin by deriving the governing equations
for the problem in §2, followed by the construction of the associated
analytical solutions in §3. The solutions are discussed in §4, with particular
attention paid to the physiological importance of the results.

##  2 Equations of flow

We consider the flow of an incompressible Newtonian fluid driven by a
pulsatile pressure gradient through a pipe of length $L^{*}$ and radius
$R^{*}$, where $\epsilon=R^{*}/L^{*}\ll 1$. We model the site of anastomosis
by an $O\left(\epsilon\right)$ axisymmetric deformation to the pipe wall
centred about $z^{*}=0$. Here, $\left(r^{*},z^{*}\right)$ are the radial and
axial components of a cylindrical polar coordinate system, with $z^{*}$
pointing in the direction of flow, $r^{*}$ pointing outwardly from the centre
of the pipe, with time $t^{*}>0$, and ∗ denoting dimensionful quantities. A
definition sketch of the geometry is included in figure 1b. We study the flow
driven by an oscillatory pressure gradient of the form

$p_{0,z^{*}}^{*\pm}\left(r^{*},z^{*},t^{*}\right)=P_{0}^{*}+P_{1}^{*}e^{i2\pi\omega^{*}t^{*}},$
(1)

for some constants $P_{0}^{*}$ and $P_{1}^{*}$, with $\omega^{*}$ being the
frequency of oscillation. Here, and throughout this paper, we use subscript
notation to denote partial differentiation. We note that, when discussing the
results of this model, we will consider only the real components of fluid
pressure and velocity, however for algebraic simplicity we will not include
this in our notation. We denote the axisymmetric fluid velocity by
${\mathbf{u}^{*}\left(r^{*},z^{*},t^{*}\right)=\left(u^{*}\left(r^{*},z^{*},t^{*}\right),0,w^{*}\left(r^{*},z^{*},t^{*}\right)\right)}$,
and introduce the dimensionless variables

$\displaystyle r^{*}$ $\displaystyle=$ $\displaystyle R^{*}r,\enspace
z^{*}=\epsilon^{-1}R^{*}z,\enspace t^{*}=\left(2\pi\omega\right)^{*-1}t,$
$\displaystyle u^{*}$ $\displaystyle=$
$\displaystyle\epsilon\omega^{*}R^{*}u,\enspace
w^{*}=\omega^{*}R^{*}w,\enspace p^{*}=\epsilon^{-1}\omega^{*}\mu^{*}~{}p,$ (2)

where $\mu^{*}$ is the fluid viscosity. The flow is characterised by the
dimensionless Womersley number (Womersley, 1955), a measure of the relative
importance of the pulsatile flow frequency to viscous effects, namely

$\alpha=R^{*}\sqrt{\frac{\rho^{*}\omega^{*}}{\mu^{*}}},$ (3)

with $\rho^{*}$ the fluid density. We model the site of anastomosis by an
axisymmetric Gaussian bump, with the boundary location given by

$r=h^{\pm}\left(z\right)=1\pm\epsilon e^{-z^{2}/\beta^{2}},$ (4)

where $\beta$ is the width parameter for the deformation. Here, the case of a
protruding bump ($h^{+}$ in (4)) is a model for an arterial coupling device,
while the intruding bump ($h^{-}$ in (4)) is a first approximation towards
modelling sutures. We write the velocity of the fluid under axisymmetry as

$\mathbf{u}^{\pm}\left(r,z,t\right)=\epsilon
u^{\pm}\left(r,z,t\right)\hat{\mathbf{r}}+w^{\pm}\left(r,z,t\right)\hat{\mathbf{z}},\quad\left(r,z\right)\in\mathcal{D}^{\pm}\left(t\right),\
t>0,$ (5)

where $\hat{\mathbf{r}}$ and $\hat{\mathbf{z}}$ are unit vectors in the radial
and axial directions respectively, and the domain occupied by the fluid at
time $t$ is then given by

$\mathcal{D}^{\pm}\left(t\right)=\left\\{\left(r,\theta,z\right)\in\mathbb{R}^{3}:r\in\left[0,h^{\pm}\left(z\right)\right),\
\theta\in\left[0,2\pi\right),z\in\left(-\infty,\infty\right)\right\\},\quad
t\in\left[0,\infty\right).$ (6)

In what follows we will drop the superscript ± for ease of notation.

Figure 1: Definition sketches showing (a) the surgical coupler (redrawn from
Wain et al. (2014)) and (b) the modelled wall geometry post-surgical coupling.

The flow is governed by the Navier-Stokes equations, which are given, in terms
of the dimensionless variables (2), in component form, as

$\displaystyle\epsilon^{2}\alpha^{2}\left[u_{t}+\epsilon\left(uu_{r}+wu_{z}\right)\right]$
$\displaystyle=$ $\displaystyle-
p_{r}+\epsilon^{2}\left[\frac{1}{r}\left(ru_{r}\right)_{r}-\frac{1}{r^{2}}u\right]+\epsilon^{4}u_{zz},$
(7)
$\displaystyle\alpha^{2}\left[w_{t}+\epsilon\left(uw_{r}+ww_{z}\right)\right]$
$\displaystyle=$ $\displaystyle-
p_{z}+\frac{1}{r}\left(rw_{r}\right)_{r}+\epsilon^{2}w_{zz},$ (8)

together with conservation of mass,

$\frac{1}{r}\left(ru\right)_{r}+w_{z}=0,$ (9)

each evaluated at all points in the domain
$\left(r,z\right)\in\mathcal{D}\left(t\right)$, $t>0$, and completed by the
no-slip, and no-penetration boundary conditions, namely

$\displaystyle\mathbf{u}\cdot\hat{\mathbf{t}}$
$\displaystyle=0\quad\mbox{on}\quad r=h\left(z\right),\
z\in\left(-\infty,\infty\right),\ t>0;$ (10)
$\displaystyle\mathbf{u}\cdot\hat{\mathbf{n}}$
$\displaystyle=0\quad\mbox{on}\quad r=h\left(z\right),\
z\in\left(-\infty,\infty\right),\ t>0;$ (11)

and the axisymmetry conditions,

$u=0,\quad w_{r}=0\quad\mbox{on}\quad r=0,\ z\in\left(-\infty,\infty\right),\
t>0.$ (12)

Here, the tangent vector in the direction of flow, $\hat{\mathbf{t}}$, and
outwardly pointing normal vector, $\hat{\mathbf{n}}$, are given, after
expanding in powers of $\epsilon$, as

$\displaystyle\hat{\mathbf{t}}$ $\displaystyle=$
$\displaystyle\left(\mp\epsilon\frac{2z}{\beta^{2}}e^{-z^{2}/\beta^{2}}\pm\epsilon^{3}\frac{4z^{3}}{\beta^{6}}e^{-{3}z^{2}/\beta^{2}}+O\left(\epsilon^{5}\right)\right)\hat{\mathbf{r}}$
(13)
$\displaystyle\qquad\qquad+\left(1-\epsilon^{2}\frac{2z^{2}}{\beta^{4}}e^{-{2}z^{2}/\beta^{2}}+O\left(\epsilon^{4}\right)\right)\hat{\mathbf{z}},$
$\displaystyle\hat{\mathbf{n}}$ $\displaystyle=$
$\displaystyle\left(1-\epsilon^{2}\frac{2z^{2}}{\beta^{4}}e^{-{2}z^{2}/\beta^{2}}+O\left(\epsilon^{4}\right)\right)\hat{\mathbf{r}}$
(14)
$\displaystyle\qquad\qquad+\left(\pm\epsilon\frac{2z}{\beta^{2}}e^{-z^{2}/\beta^{2}}\mp\epsilon^{3}\frac{4z^{3}}{\beta^{6}}e^{-{3}z^{2}/\beta^{2}}+O\left(\epsilon^{5}\right)\right)\hat{\mathbf{z}}.$

We expand the boundary conditions on the wall, (10) and (11), about $r=1$,
utilising the $O\left(\epsilon\right)$ nature of the boundary perturbation,
resulting in the simplified the no-slip condition

$w\left(1,z,t\right)\pm\epsilon
e^{-z^{2}/\beta^{2}}w_{r}\left(1,z,t\right)+O\left(\epsilon^{2}\right)=0,$
(15)

and equivalent no-penetration condition

$\displaystyle u\left(1,z,t\right)$ $\displaystyle\pm$
$\displaystyle\frac{2z}{\beta^{2}}e^{-z^{2}/\beta^{2}}w\left(1,z,t\right)$
(16)
$\displaystyle\quad+\epsilon\left(\frac{2z}{\beta^{2}}e^{-{2}z^{2}/\beta^{2}}w_{r}\left(1,z,t\right)\pm
e^{-z^{2}/\beta^{2}}u_{r}\left(1,z,t\right)\right)+O\left(\epsilon^{2}\right)=0,$

each with $z\in\left(-\infty,\infty\right)$ and $t>0$. Finally, the flow is
driven by the oscillatory pressure gradient

$p_{0,z}=P_{0}+P_{1}e^{it},\qquad r\in\left[0,1\right),\
z\in\left(-\infty,\infty\right),\ t>0,$ (17)

with $P_{0}$ and $P_{1}$ dimensionless constants.

We will now solve the boundary value problem (7) - (9), and (15) - (17),
hereafter termed [BVP], through constructing the $O\left(\epsilon\right)$
correction to the flow caused by the anastomosis.

##  3 Analytic solution to the flow across an arterial anastomosis

We begin by constructing the leading order approximation to [BVP]. Writing
${w^{\pm}=w_{0}+O\left(\epsilon\right)}$, we obtain the leading order problem

$\displaystyle\alpha^{2}w_{0,t}-\frac{1}{r}\left(rw_{0,r}\right)_{r}+P_{0}+P_{1}e^{it}$
$\displaystyle=$ $\displaystyle 0,\qquad-\infty<z<\infty,\ 0\leq r<1,\ t>0;$
(18) $\displaystyle w_{0,r}\left(0,z,t\right)$ $\displaystyle=$ $\displaystyle
0,\qquad-\infty<z<\infty,\ t>0;$ (19) $\displaystyle w_{0}\left(1,z,t\right)$
$\displaystyle=$ $\displaystyle 0,\qquad-\infty<z<\infty,\ t>0.$ (20)

The solution to the boundary value problem (18) - (20) is given by Womersley
(1955) for the case $P_{0}=0$, $P_{1}\neq 0$, and, in the case $P_{0}\neq 0$,
$P_{1}=0$, by the solution to axisymmetric flow through a pipe. Thus, we have

$w_{0}\left(r,z,t\right)=\frac{1}{4}P_{0}\left(r^{2}-1\right)+\frac{iP_{1}}{\alpha^{2}}\left(1-\frac{J_{0}\left(i^{3/2}\alpha
r\right)}{J_{0}\left(i^{3/2}\alpha\right)}\right)e^{it},$ (21)

on $-\infty<z<\infty$, and $0\leq r<1$, for $t>0$, where $J_{k}$ is the
$k^{th}$ order Bessel function of the first kind. It is important to note that
the anastomosis has no effect on the leading order contribution to the flow.

We now turn our attention to calculating the $O\left(\epsilon\right)$
correction to the flow caused by the anastomosis. Writing
${u^{\pm}=u_{1}+O\left(\epsilon\right)}$, application of the solution (21) to
the continuity equation (9), together with the no-penetration condition (16),
immediately gives that $u_{1}\equiv 0$. Furthermore, following the form of the
boundary conditions (15) and (16), we expand the velocity and pressure terms
$w$, and $p$ as

$w=w_{0}\left(r,z,t\right)+\epsilon
w_{1}\left(r,z,t\right)+O\left(\epsilon^{2}\right),\quad
p=p_{0}\left(z,t\right)+\epsilon
p_{1}\left(r,z,t\right)+O\left(\epsilon^{2}\right),$ (22)

with $w_{0}$ and $p_{0,z}$ given by (21), and (17) respectively. On
substitution of (22) into [BVP] we obtain the boundary value problem at next
order, namely

$\displaystyle\alpha^{2}w_{1,t}-\frac{1}{r}\left(rw_{1,r}\right)_{r}+p_{1,z}$
$\displaystyle=$ $\displaystyle 0,\qquad-\infty<z<\infty,\ 0\leq r<1,\ t>0;$
(23) $\displaystyle p_{1,r}$ $\displaystyle=$ $\displaystyle
0,\qquad-\infty<z<\infty,\ 0\leq r<1,\ t>0;$ (24) $\displaystyle
w_{1,r}\left(0,z,t\right)$ $\displaystyle=$ $\displaystyle
0,\qquad-\infty<z<\infty,\ t>0;$ (25) $\displaystyle
w_{1}\left(1,z,t\right)\pm e^{-z^{2}/\beta^{2}}w_{0,r}\left(1,z,t\right)$
$\displaystyle=$ $\displaystyle 0,\qquad-\infty<z<\infty,\ t>0;$ (26)

where $w_{0,r}\left(1,z,t\right)$, calculated via (21), is given by

$w_{0,r}\left(1,z,t\right)=\frac{1}{2}P_{0}-\frac{i^{1/2}\bar{\mathrm{J}}P_{1}}{\alpha}e^{it},$
(27)

with the constant
$\bar{\mathrm{J}}=J_{1}\left(i^{3/2}\alpha\right)/J_{0}\left(i^{3/2}\alpha\right)$
having been introduced for algebraic convenience. The boundary value problem
(23) - (26) admits the exact solution

$\displaystyle\\!\\!\\!\\!\\!\\!w_{1}$ $\displaystyle=$
$\displaystyle\pm\left[\frac{P_{0}}{2}\left(1-2r^{2}\right)+\frac{i2P_{1}\bar{\mathrm{J}}}{\alpha\left(\alpha+2\bar{\mathrm{J}}i^{1/2}\right)}\left(\bar{\mathrm{J}}-i^{3/2}\alpha\frac{J_{0}\left(i^{3/2}\alpha
r\right)}{J_{0}\left(i^{3/2}\alpha\right)}\right)e^{it}\right]e^{-z^{2}/\beta^{2}},$
(28) $\displaystyle\\!\\!\\!\\!\\!\\!p_{1,z}$ $\displaystyle=$
$\displaystyle\mp\left[4P_{0}-\frac{2\alpha
P_{1}}{\alpha+2\bar{\mathrm{J}}i^{1/2}}\bar{\mathrm{J}}^{2}e^{it}\right]e^{-z^{2}/\beta^{2}}.$
(29)

For completeness, the first non-zero component of the radial velocity is
given, through expanding $u^{\pm}=\epsilon u_{2}+O\left(\epsilon^{2}\right)$
and solving (9), together with (16) and (28), to give

$u_{2}=\pm\left[\frac{P_{0}}{2\beta^{2}}\left(r-r^{3}\right)z+\frac{i2P_{1}\bar{\mathrm{J}}}{\alpha\beta^{2}\left(\alpha+2\bar{\mathrm{J}}i^{1/2}\right)}\left(\bar{\mathrm{J}}r-\frac{J_{1}\left(i^{3/2}\alpha
r\right)}{J_{0}\left(i^{3/2}\alpha\right)}\right)ze^{it}\right]e^{-z^{2}/\beta^{2}}.$
(30)

Through substitution of the correction terms (28), (29), and (30), together
with the leading order solution (21), into the fluid velocity expansion (5)
and taking the real part we obtain an analytical expression for the flow
through a microvascular anastomosis. This modelled velocity profile enables
the detailed and methodical investigation into the physical effects of
introducing of a coupling device, as well as providing a first comparison with
the case of sutured anastomosis.

##  4 Results and analysis of the shear rate across the anastomosis

In §3 we constructed the $O\left(\epsilon\right)$ correction term to the
Womersley flow profile, which arrises from modelling the case of a
microvascular anastomosis through the inclusion of a small axisymmetric
perturbation to the radial geometry. Where previous papers in this context
(see Wain et al. (2014), Wain et al. (2016), Wain et al. (2018)) have focussed
on using CFD packages such as ANSYS-CFX to solve anastomoses problems with
realistic geometries, the mathematical approach used here provides the basis
for a detailed, rational exploration of the effect of conjoining small
arteries with a radially symmetric coupling device.

Throughout this section we will choose a realistic parameter set for the
calculation and analysis of the results, namely the following. As the height
of the bump due to a coupling device is difficult to accurately measure, in
the present work we take the small parameter $\epsilon=0.04/1.25$ to be equal
to the size of a typical suture $\left(0.04\mathrm{~{}mm}\right)$ relative to
the typical radius of the vessels we are investigating
$\left(1.25\mathrm{~{}mm}\right)$. The width parameter $\beta=0.3$ is chosen
arbitrarily for our calculations. The viscosity and density of blood give
$\mu^{*}=0.0035$ Pa s, and $\rho^{*}=1060$ kg m-3, and choosing the beat
frequency $\omega^{*}=1.5$ Hz then results in an approximate Womersley number
$\alpha\approx 0.84$. Finally the pressure scales in (17) are chosen to give a
maximum dimensional velocity of $0.62$ ms-1 (to match Wain et al. (2018)) and
to ensure no back flow, thus $P_{0}=P_{1}\approx-660$.

In figure 2 we see (a) the leading order and (c) the $O\left(\epsilon\right)$
correction to the axial velocity through a surgically coupled microvascular
anastomosis, with (b) showing, for completeness, the
$O\left(\epsilon^{2}\right)$ correction to the radial velocity, and (d)
showing a selection of streamlines across the coupling. We see that, as
expected, there is a small deformation to the streamlines, following the
vessel wall, when compared to the pristine case. Far from the site of
anastomosis the velocity profile relaxes to the pristine case.

Figure 2: Velocity profile across the coupled vessel at time $t=0$. In figure
(a) we see the leading order contribution to the axial velocity, with the
corrections to the velocity in the radial
$\left(O\left(\epsilon^{2}\right)\right)$ and axial
$\left(O\left(\epsilon\right)\right)$ directions in figures (b) and (c)
respectively. Figure (d) shows a selection of streamlines through the not-to-
scale coupled-vessel geometry and shear rate of figure 3c.

As discussed in §1, key to the rate of thrombus formation is both the shear
rate and the spatial gradient of shear rate across the site of anastomosis.
The dimensionless shear rate $\dot{\gamma}$, non-dimensionalised with respect
to the frequency parameter $\omega^{*-1}$, is given, on substitution of the
real part of the fluid velocity (5), as

$\dot{\gamma}=\left[\left(\frac{\partial w}{\partial
r}\right)^{2}+\epsilon^{2}\left(\frac{2u^{2}}{r^{2}}+\left(\frac{\partial{u}}{\partial{r}}\right)^{2}+\left(\frac{\partial{w}}{\partial{z}}\right)^{2}+2\frac{\partial{u}}{\partial{z}}\frac{\partial{w}}{\partial{r}}\right)+\epsilon^{4}\left(\frac{\partial{u}}{\partial{z}}\right)^{2}\right]^{1/2}.$
(31)

Substitution of the asymptotic forms for $u$ and $w$ into (31), provided
$\partial w_{0}/\partial r\gg\epsilon$, yields the expansion

$\dot{\gamma}\sim\left|\frac{\partial{w_{0}^{R}}}{\partial{r}}\right|+\epsilon\,\mathrm{sgn}\left(\frac{\partial{w_{0}^{R}}}{\partial{r}}\right)\frac{\partial{w_{1}^{R}}}{\partial{r}},$
(32)

where $\mathrm{sgn}$ is the signum function. It is clear from (32), together
with (28), that the correction to the shear rate in the case of an intruding
bump differs from that of a protruding bump only in sign. Figure 3 plots (a)
the leading-order and (b) the $O\left(\epsilon\right)$ correction to the
dimensionless shear rate with the parameter set defined above. Here we see
that, in the vicinity of the anastomosis, the effect of the coupling device is
a decrease in shear rate when compared to the pristine vessel, whereas the
sutures effect a corresponding increase. The maximum change in shear rate
compared to the pristine case is $\sim 15\%$, and is found at the peak of the
anastomosis. While this is much lower in the sutured case than the increase
calculated by Wain et al. (2018), this is to be expected: in the case of a
more realistic model of a suture the wall deformation would be concentrated
over a smaller area, which is likely to effect a greater increase in shear
rate than the initial model given here. It has also been noted by Wain et al.
(2016) that the angle of suture placement can have a significant effect on the
shear stress when compared to those placed in-line with the direction of flow.

Detailed examination of the asymptotic form of the shear rate (32), together
with the axial velocity (28), reveals that there is a linear relationship
between the height of the wall deformation and the corresponding change to the
shear rate, while the width parameter $\beta$ has no effect on the magnitude
of the $O\left(\epsilon\right)$ change in the shear rate, but controls the
spatial gradient of the shear rate across the anastomosis through the spatial
exponential factor in (28).

We regards to the overarching aim of reducing the rates of thrombosis in
microvascular anastomoses we can now draw the following conclusions. The
results we have presented here agree with existing studies, showing a
significant decrease in the shear rates due to the use of a surgical coupling
device when compared to sutures, and thus coupling devices should be less
likely to lead to thrombus growth. However, as discussed in §1, the magnitude
of the shear rate does not tell the whole story: the spatial gradients of
shear rate are also linked to platelet aggregation, and thus to increased
rates of thrombosis. Through the detailed mathematical work presented here we
have shown that the spatial gradients of shear rate, across the site of
anastomosis, can be reduced by decreasing the gradient of the wall at the
join, even if the maximum displacement on the wall remains the same (i.e. by
increasing $\beta$ in the model).

Figure 3: Dimensionless shear rate across the vessel wall with time,
calculated in (32). In figure (a) we see the leading order spatially
homogeneous contribution to the shear rate along the wall. The
$O\left(\epsilon\right)$ correction to the shear rate is shown in (b), with
the combined shear rate shown on the not-to-scale coupled-vessel geometry.

##  5 Discussion

In this paper we have modelled the effect of microvascular anastomosis in
small vessels of approximately $1$ mm in radius on the flow driven by a
pulsatile pressure gradient, with wall perturbations taking the form of both
extruding and intruding axisymmetric deformations relating to the use of a
coupling device as well as (a first analytical approximation to) the use of
sutures. Through this work, the first analytical approach to these problems,
we have been able to replicate the key findings of recent numerical work in
this area, while providing new detailed insight into the ways in which new
devices could be designed to minimise areas of potential biomedical impact.
Key to this is the discovery that the width of the wall perturbation plays the
significant role in the biologically important shear rate gradient, an insight
which would not have been immediately apparent without the detailed analytical
results presented here.

It is important to note that the asymptotic expansion for the shear rate (32)
is only valid providing $\left|\partial w_{0}/\partial r\right|\gg\epsilon$.
When this is not the case for some point in time, i.e. in the case of flow
reversal, a more detailed asymptotic examination is required. An interesting
extension of this work would be to investigate the effect of a more accurate
pressure gradient (for example, that derived from Doppler velocimetry in Wain
et al. (2018)). It is likely that in this case, some flow reversal would be
apparent, requiring the more careful construction of the shear rate expansion.

A potential extension of the present work lies in the assumption that the
fluid is Newtonian i.e. there is a linear relationship between the viscous
stresses in the fluid and the strain rate. In regions of high shear rate this
is likely to be a very good approximation, however when shear rates are more
moderate a more physiologically accurate constitutive viscosity relationship
may be required. In this case it is unlikely that a fully analytic solution
would be forthcoming and as a result some amount of numerical approximation
would be required.

The results from this model of microvascular anastomosis agree with the main
findings from existing numerical studies in this area, while also providing
new insight into potential ways for further reducing the thrombosis rates of
surgical coupling devices through decreasing the change in spatial gradients
of the shear rate. These findings support the growing clinical evidence in
favour of using such devices in preference to the more traditional methods of
suturing.

## Acknowledgements

M.T.G. is supported by the Engineering and Physical Sciences Research Council
award EP/N021096/1. The author would like to thank David Smith and Thomas
Montenegro-Johnson for helpful insights and discussions, Rich Wain and Justin
Whitty for surgical and engineering related discussions, as well as the
members of the Multi-scale Biology Study Group, University of Birmingham
(12-15th December 2016) for valuable discussions regarding this subject.

The Multi-scale Biology Study Group, University of Birmingham (12-15th
December 2016) was jointly funded by POEMS (Predictive modelling for
healthcare technology through maths - EP/L001101/1) and MSB-Net (UK Multi-
Scale Biology Network - BB/M025888/1).

## References

  * de Bruijn & Marck (1996) de Bruijn, H. P. & Marck, K. W. 1996 Coupling the venous anastomosis: safe and simple. Microsurg. 17 (7), 414–416.
  * Carrel (1902) Carrel, A. 1902 La techniue operatoire des anastomoses vasculaires et la transplantation des visceres. Lyon. Med. pp. 859–864.
  * Grewal et al. (2012) Grewal, A. S., Erovic, B., Strumas, N., Enepekides, D. J. & Higgins, K. M. 2012 The utility of the microvascular anastomotic coupler in free tissue transfer. Can. J. Plast. Surg. 20 (2), 98–102.
  * Jandali et al. (2010) Jandali, S., Wu, L. C., Vega, S. J., Kovach, S. J. & Serletti, J. M. 2010 1000 consecutive venous anastomoses using the microvascular anastomotic coupler in breast reconstruction. Plast. Reconstr. Surg. 125 (3), 792–798.
  * Koul et al. (2013) Koul, A. R., Patil, R. K. & Nahar, S. 2013 Unfavourable results in free tissue transfer. Indian J. Plast. Surg. 46 (2), 247.
  * Kroll et al. (1996) Kroll, S. S., Schusterman, M. A., Reece, G. P., Miller, M. J., Evans, G. R., Robb, G. L. & Baldwin, B. J. 1996 Choice of flap and incidence of free flap success. Plast. Reconstr. Surg. 98 (3), 459–463.
  * Murray (1926) Murray, C. D. 1926 The physiological principle of minimum work i. the vascular system and the cost of blood volume. Proc. Natl. Acad. Sci. USA 12 (3), 207–214.
  * Nesbitt et al. (2009) Nesbitt, W. S., Westein, E., Tovar-Lopez, F. J., Tolouei, E., Mitchell, A., Fu, J., Carberry, J., Fouras, A. & Jackson, S. P. 2009 A shear gradient–dependent platelet aggregation mechanism drives thrombus formation. Nat. Med. 15 (6), 665.
  * Ostrup & Berggren (1986) Ostrup, L. T. & Berggren, A. 1986 The UNILINK instrument system for fast and safe microvascular anastomosis. Ann. Plas. Surg. 17 (6), 521–525.
  * Pedley (1980) Pedley, T. J. 1980 The Fluid Mechanics of Large Blood Vessels. Cambridge University Press.
  * Sakariassen et al. (2015) Sakariassen, K. S., Orning, L. & Turitto, V. T. 2015 The impact of blood shear rate on arterial thrombus formation. Future Sci. OA 1 (4).
  * Spector et al. (2006) Spector, J. A., Draper, L. B., Levine, J. P. & Ahn, C. Y. 2006 Routine use of microvascular coupling device for arterial anastomosis in breast reconstruction. Ann. Plas. Surg. 56 (4), 365–368.
  * Van Dyke (1964) Van Dyke, M. 1964 Perturbation Methods in Fluid Mechanics. Applied mathematics and mechanics v. 8. Academic Press.
  * Wain et al. (2016) Wain, R. A. J., Hammond, D., McPhillips, M., Whitty, J. P. M. & Ahmed, W. 2016 Microarterial anastomoses: A parameterised computational study examining the effect of suture position on intravascular blood flow. Microvasc. Res. 105, 141–148.
  * Wain et al. (2018) Wain, R. A. J., Smith, D. J., Hammond, D. R. & Whitty, J. P. M. 2018 Influence of microvascular sutures on shear strain rate in realistic pulsatile flow. Microvasc. Res. .
  * Wain et al. (2014) Wain, R. A. J., Whitty, J. P. M., Dalal, M. D., Holmes, M. C. & Ahmed, W. 2014 Blood flow through sutured and coupled microvascular anastomoses: a comparative computational study. J. Plast. Reconstr. Aes. 67 (7), 951–959.
  * Womersley (1955) Womersley, J.R. 1955 Method for the calculation of velocity, rate of flow and viscous drag in arteries when the pressure gradient is known. J. Physol. 127 (3), 553–563.
  * Yap et al. (2006) Yap, L. H., Constantinides, J. & Butler, C. E. 2006 Venous thrombosis in coupled versus sutured microvascular anastomoses. Ann. Plas. Surg. 57 (6), 666–669.