∎

11institutetext: Xin-long Luo 22institutetext: Corresponding author. School of
Artificial Intelligence,
Beijing University of Posts and Telecommunications, P. O. Box 101,
Xitucheng Road No. 10, Haidian District, 100876, Beijing China
22email<EMAIL_ADDRESS>33institutetext: Jia-hui Lv 44institutetext:
School of Artificial Intelligence,
Beijing University of Posts and Telecommunications, P. O. Box 101,
Xitucheng Road No. 10, Haidian District, 100876, Beijing China
44email<EMAIL_ADDRESS>55institutetext: Hang Xiao 66institutetext: School
of Artificial Intelligence,
Beijing University of Posts and Telecommunications, P. O. Box 101,
Xitucheng Road No. 10, Haidian District, 100876, Beijing China
66email<EMAIL_ADDRESS>
# Explicit continuation methods with L-BFGS updating formulas for linearly
constrained optimization problems

Xin-long Luo$\ast$ Jia-hui Lv  Hang Xiao

(Received: date / Accepted: date)

###### Abstract

This paper considers an explicit continuation method with the trusty time-
stepping scheme and the limited-memory BFGS (L-BFGS) updating formula (Eptctr)
for the linearly constrained optimization problem. At every iteration, Eptctr
only involves three pairs of the inner product of vector and one matrix-vector
product, other than the traditional and representative optimization method
such as the sequential quadratic programming (SQP) or the latest continuation
method such as Ptctr LLS2020 , which needs to solve a quadratic programming
subproblem (SQP) or a linear system of equations (Ptctr). Thus, Eptctr can
save much more computational time than SQP or Ptctr. Numerical results also
show that the consumed time of EPtctr is about one tenth of that of Ptctr or
one fifteenth to 0.4 percent of that of SQP. Furthermore, Eptctr can save the
storage space of an $(n+m)\times(n+m)$ large-scale matrix, in comparison to
SQP. The required memory of Eptctr is about one fifth of that of SQP. Finally,
we also give the global convergence analysis of the new method under the
standard assumptions.

###### Keywords:

continuation method trust-region method SQP structure-preserving algorithm
generalized projected gradient flow large-scale optimization

###### MSC:

65J15 65K05 65L05

††journal: Journal of XXX

##  1 Introduction

In this article, we consider the following linearly equality-constrained
optimization problem

$\displaystyle\min_{x\in\Re^{n}}\;f(x)$ $\displaystyle\text{subject
to}\;\;Ax=b,$ (1)

where matrix $A\in\Re^{m\times n}$ and vector $b\in\Re^{m}$ may have random
noise. This problem has many applications in engineering fields such as the
visual-inertial navigation of an unmanned aerial vehicle maintaining the
horizontal flight CMFO2009 ; LLS2020 , and there are many practical methods to
solve it such as the sequential quadratic programming (SQP) method
Bertsekas2018 ; NW1999  or the penalty function method FM1990 .

For the constrained optimization problem (1), the continuation method AG2003 ;
CKK2003 ; Goh2011 ; KLQCRW2008 ; Pan1992 ; Tanabe1980  is another method other
than the traditional optimization method such as SQP or the penalty function
method. The advantage of the continuation method over the SQP method is that
the continuation method is capable of finding many local optimal points of the
non-convex optimization problem by tracking its trajectory, and it is even
possible to find the global optimal solution BB1989 ; Schropp2000 ;
Yamashita1980 . However, the computational efficiency of the continuation
method may be higher than that of SQP. Recently, Luo, Lv and Sun LLS2020  give
a continuation method with the trusty time-stepping scheme and its consumed
time is about one fifth of that of SQP for the linearly constrained
optimization problem (1). Their method only needs to solve a linear system of
equations with an $n\times n$ symmetric definite coefficient matrix at every
iteration, which involves about $\frac{1}{3}n^{3}$ flops. SQP needs to solve a
linear system of equations with an $(m+n)\times(m+n)$ coefficient matrix,
which involves about $\frac{2}{3}(m+n)^{3}$ flops. In order to improve the
computational efficiency further and save the storage of the continuation
method LLS2020  for the large-scale optimization problem, we consider a
special limited-memory BFGS updating formula and the trusty time-stepping
scheme in this article.

The rest of the paper is organized as follows. In section 2, we give a new
continuation method with the trusty time-stepping scheme and the L-BFGS
updating formula for the linearly equality-constrained optimization problem
(1). In section 3, we analyze the global convergence of this new method. In
section 4, we report some promising numerical results of the new method, in
comparison to the traditional optimization method (SQP) and the latest
continuation method (Ptctr) for some large-scale problems. Finally, we give
some discussions and conclusions in section 5.

##  2 The explicit continuation method with L-BFGS updating formulas

In this section, we construct an explicit continuation method with the
adaptive time-stepping scheme based on the trust-region updating strategy
Yuan2015  for the linearly equality-constrained optimization problem (1).
Firstly, we construct a generalized projected gradient flow based on the KKT
conditions of linearly constrained optimization problem. Then, in order to
efficiently follow the generalized gradient flow, we construct an explicit
continuation method with an adaptive time-stepping scheme for this special
ordinary differential equations (ODEs). Furthermore, we give a preprocessing
method for the infeasible initial point.

###  2.1 The generalized projected gradient flow

For the linearly constrained optimization problem (1), it is well known that
its optimal solution $x^{\ast}$ needs to satisfy the Karush-Kuhn-Tucker
conditions (p. 328, NW1999 ) as follows:

$\displaystyle\nabla_{x}L(x,\,\lambda)$ $\displaystyle=\nabla
f(x)+A^{T}\lambda=0,$ (2) $\displaystyle Ax-b$ $\displaystyle=0,$ (3)

where the Lagrangian function $L(x,\,\lambda)$ is defined by

$\displaystyle L(x,\,\lambda)=f(x)+\lambda^{T}(Ax-b).$ (4)

Similarly to the method of the negative gradient flow for the unconstrained
optimization problem LKLT2009 , from the first-order necessary conditions
(2)-(3), we can construct a dynamical system of differential-algebraic
equations for problem (1) LL2010 ; Luo2012 ; LLW2013 ; Schropp2003  as
follows:

$\displaystyle\frac{dx}{dt}=-\nabla L_{x}(x,\,\lambda)=-\left(\nabla
f(x)+A^{T}\lambda\right),$ (5) $\displaystyle Ax-b=0.$ (6)

By differentiating the algebraic constraint (6) with respect to $t$ and
replacing it into the differential equation (5), we obtain

$\displaystyle A\frac{dx}{dt}=-A\left(\nabla f(x)+A^{T}\lambda\right)=-A\nabla
f(x)-AA^{T}\lambda=0.$ (7)

If we assume that matrix $A$ has full row rank further, from equation (7), we
obtain

$\displaystyle\lambda=-\left(AA^{T}\right)^{-1}A\nabla f(x).$ (8)

By replacing $\lambda$ of equation (8) into equation (5), we obtain the
dynamical system of ordinary differential equations (ODEs) as follows:

$\displaystyle\frac{dx}{dt}=-\left(I-A^{T}\left(AA^{T}\right)^{-1}A\right)\nabla
f(x).$ (9)

Thus, we also obtain the projected gradient flow for the constrained
optimization problem Tanabe1980 .

For convenience, we denote the projection matrix $P$ as

$\displaystyle P=I-A^{T}\left(AA^{T}\right)^{-1}A.$ (10)

It is not difficult to verify $P^{2}=P$ and $AP=0$. That is to say, $P$ is a
symmetric projection matrix and its eigenvalues are 0 or 1. From Theorem 2.3.1
in p. 73 of GV2013 , we know that its matrix 2-norm is

$\displaystyle\|P\|=1.$ (11)

We denote $P^{+}$ as the Moore-Penrose generalized inverse of $P$ (p. 11,
SY2006 ). Since $P$ is symmetric and its eigenvalues are 0 or 1, it is not
difficult to verify

$\displaystyle P^{+}=P.$ (12)

Thus, for a full rank matrix $B\in\Re^{n\times n}$, we obtain the generalized
inverse $(PB)^{+}$ of $PB$ as follows:

$\displaystyle(PB)^{+}=B^{+}P^{+}=B^{-1}P.$ (13)

Similarly to the generalized gradient flow for an unconstrained optimization
problem (p. 361, HM1996 ), from the projected gradient flow (9), we can
construct the generalized projected gradient flow for the constrained
optimization problem (1) as follows:

$\displaystyle\frac{dx}{dt}=-(PH(x)P)P\nabla f(x),\;x(0)=x_{0},$ (14)

where $H(x)$ is a symmetric positive definite matrix for any $x\in\Re^{n}$.
Here, $H(x)$ may be selected as the inverse of the Hessian matrix
$\nabla^{2}f(x)$ of $f(x)$ and $PH(x)P$ can be regarded as a pre-conditioner
of $P\nabla f(x)$ to mitigate the stiffness of the ODEs (14). Consequently, we
can adopt the explicit numerical method to compute the trajectory of the ODEs
(14) efficiently LXL2020 ; LY2021 .

###### Remark 1

If $x(t)$ is the solution of the ODEs (14), it is not difficult to verify that
$x(t)$ satisfies $A(dx/dt)=0$. That is to say, if the initial point $x_{0}$
satisfies $Ax_{0}=b$, the solution $x(t)$ of the generalized projected
gradient flow (14) also satisfies $Ax(t)=b,\;\forall t\geq 0$. This property
is very useful when we construct a structure-preserving algorithm HLW2006 ;
Simos2013  to follow the trajectory of the ODEs (14) and obtain its
equilibrium point $x^{\ast}$.

###### Remark 2

If we assume that $x(t)$ is the solution of the ODEs (14), from equations
(10)-(11) and the positive definite property of $H(x)$, we obtain

$\displaystyle\frac{df(x)}{dt}$ $\displaystyle=\left(\nabla
f(x)\right)^{T}\frac{dx}{dt}=-(\nabla f(x))^{T}PH(x)P\nabla f(x)$
$\displaystyle=-(P\nabla f(x))^{T}H(x)(P\nabla f(x))\leq 0.$

That is to say, $f(x)$ is monotonically decreasing along the solution curve
$x(t)$ of the dynamical system (14). Furthermore, the solution $x(t)$
converges to $x^{\ast}$ when $f(x)$ is lower bounded and $t$ tends to infinity
HM1996 ; Schropp2000 ; Tanabe1980 , where $x^{\ast}$ satisfies the first-order
Karush-Kuhn-Tucker conditions (2)-(3). Thus, we can follow the trajectory
$x(t)$ of the ODEs (14) to obtain its equilibrium point $x^{\ast}$, which is
also one saddle point of the original optimization problem (1).

###  2.2 The explicit continuation method

The solution curve of the degenerate ordinary differential equations is not
efficiently followed on an infinite interval by the traditional ODE method
AP1998 ; BCP1996 ; LF2000 , so one needs to construct the particular method
for this problem (14). We apply the first-order implicit Euler method SGT2003
to the ODEs (14), then we obtain

$\displaystyle x_{k+1}=x_{k}-\Delta t_{k}(PH(x_{k}))(P\nabla f(x_{k})),$ (15)

where $\Delta t_{k}$ is the time-stepping size.

Since the system of equations (15) is a nonlinear system which is not directly
solved, we seek for its explicit approximation formula. We denote
$s_{k}=x_{k+1}-x_{k}$. By using the first-order Taylor expansion, we have the
linear approximation $\nabla f(x_{k})+\nabla^{2}f(x_{k})s_{k}$ of $\nabla
f(x_{k+1})$. By substituting it into equation (15) and using the zero-order
approximation $H(x_{k})$ of $H(x_{k+1})$, we have

$\displaystyle s_{k}$ $\displaystyle\approx-\Delta
t_{k}(PH(x_{k}))\left(P\left(\nabla
f(x_{k})+\nabla^{2}f(x_{k})s_{k}\right)\right)$ $\displaystyle=-\Delta
t_{k}(PH(x_{k}))(P\nabla f(x_{k}))-\Delta
t_{k}P(H(x_{k})P)(P\nabla^{2}f(x_{k}))s_{k}.$ (16)

From equation (15) and $P^{2}=P$, we have $Ps_{k}=s_{k}$. Let
$H(x_{k})=(\nabla^{2}f(x_{k}))^{-1}$. Then, we have
$H(x_{k})P=(P\nabla^{2}f(x_{k}))^{+}$. Thus, we regard

$\displaystyle
P(H(x_{k})P)(P\nabla^{2}f(x_{k})Ps_{k})=P(P\nabla^{2}f(x_{k}))^{+}(P\nabla^{2}f(x_{k}))Ps_{k}\approx
Ps_{k}=s_{k}.$ (17)

By substituting it into equation (16), we obtain the explicit continuation
method as follows:

$\displaystyle s_{k}=-\frac{\Delta t_{k}}{1+\Delta t_{k}}(PH_{k})(Pg_{k}),$
(18) $\displaystyle x_{k+1}=x_{k}+s_{k},$ (19)

where $g_{k}=\nabla f(x_{k})$ and $H_{k}=(\nabla^{2}f(x_{k}))^{-1}$ or its
quasi-Newton approximation in the projective space
$S_{p}^{k}=\\{x:\;x=x_{k}+Pd,\,d\in\Re^{n}\\}$.

If we let the projection matrix $P=I$, the formula (18) is equivalent to the
explicit continuation method given by Luo, Xiao and Lv LXL2020  for nonlinear
equations. The explicit continuation method (18)-(19) is similar to the
projected damped Newton method if we let $\alpha_{k}=\Delta t_{k}/(1+\Delta
t_{k})$ in equation (18). However, from the view of the ODE method, they are
different. The projected damped Newton method is obtained by the explicit
Euler scheme applied to the generalized projected gradient flow (14), and its
time-stepping size $\alpha_{k}$ is restricted by the numerical stability
SGT2003 . That is to say, the large time-stepping size $\alpha_{k}$ can not be
adopted in the steady-state phase.

The explicit continuation method (18)-(19) is obtained by the implicit Euler
approximation method applied to the generalized projected gradient flow (14),
and its time-stepping size $\Delta t_{k}$ is not restricted by the numerical
stability. Therefore, the large time-stepping size can be adopted in the
steady-state phase for the explicit continuation method (18)-(19), and it
mimics the Newton method near the equilibrium solution $x^{\ast}$ such that it
has the fast local convergence rate. The most of all, the new step size
$\alpha_{k}=\Delta t_{k}/(\Delta t_{k}+1)$ is favourable to adopt the trust-
region updating technique for adaptively adjusting the time-stepping size
$\Delta t_{k}$ such that the explicit continuation method (18)-(19) accurately
tracks the trajectory of the generalized projected gradient flow in the
transient-state phase and achieves the fast convergence rate near the
equilibrium point $x^{\ast}$.

###### Remark 3

From equation (18) and the property $AP=0$ of the projected matrix $P$, it is
not difficult to verify $As_{k}=0$. Thus, if the initial point $x_{0}$
satisfies the linear constraint $Ax_{0}=b$, the point $x_{k}$ also satisfies
the linear constraint $Ax_{k}=b$. That is to say, the explicit continuation
method (18)-(19) is a structure-preserving method.

###  2.3 The L-BFGS quasi-Newton updating formula

For the large-scale problem, the numerical evaluation of the Hessian matrix
$\nabla^{2}f(x_{k})$ consumes much time and stores an $n\times n$ matrix. In
order to overcome these two shortcomings, we use the L-BFGS quasi-Newton
formula (BNY1987 ; Goldfarb1970  or pp. 222-230, NW1999 ) to approximate the
generalized inverse $H(x_{k})P$ of $P\nabla^{2}f(x_{k})$. Recently, Ullah,
Sabi and Shah USS2020  give an efficient L-BFGS updating formula for the
system of monotone nonlinear equations. Here, in order to suit the generalized
projected gradient flow (14), we revise their L-BFGS updating formula as

$\displaystyle
H_{k+1}=\begin{cases}I-\frac{y_{k}s_{k}^{T}+s_{k}y_{k}^{T}}{y_{k}^{T}s_{k}}+2\frac{y_{k}^{T}y_{k}}{(y_{k}^{T}s_{k})^{2}}s_{k}s_{k}^{T},\;\text{if}\;|s_{k}^{T}y_{k}|>\theta\|s_{k}\|^{2},\\\
I,\;\text{otherwise}.\end{cases}$ (20)

where $s_{k}=x_{k+1}-x_{k},\;y_{k}=P\nabla f(x_{k+1})-P\nabla f(x_{k})$ and
$\theta$ is a small positive constant such as $\theta=10^{-6}$. The initial
matrix $H_{0}$ can be simply selected by the identity matrix. When
$|s_{k}^{T}y_{k}|\geq\theta\|s_{k}\|^{2}$, from equation (20), it is not
difficult to verify

$\displaystyle H_{k+1}y_{k}=\frac{y_{k}^{T}y_{k}}{y_{k}^{T}s_{k}}s_{k}.$

That is to say, $H_{k+1}$ satisfies the scaling quasi-Newton property. By
using the Sherman-Morrison-Woodburg formula, from equation (20), when
$|s_{k}^{T}y_{k}|\geq\theta\|s_{k}\|^{2}$, we have

$\displaystyle
B_{k+1}=H_{k+1}^{-1}=I-\frac{s_{k}s_{k}^{T}}{s_{k}^{T}s_{k}}+\frac{y_{k}y_{k}^{T}}{y_{k}^{T}y_{k}}.$

The L-BFGS updating formula (20) has some nice properties such as the
symmetric positive definite property and the positive lower bound of its
eigenvalues.

###### Lemma 1

Matrix $H_{k+1}$ defined by equation (20) is symmetric positive definite and
its eigenvalues are greater than $1/2$.

Proof. (i) For any nonzero vector $z\in\Re^{n}$, from equation (20), when
$|s_{k}^{T}y_{k}|>\theta\|s_{k}\|^{2}$, we have

$\displaystyle
z^{T}H_{k+1}z=\|z\|^{2}-2{(z^{T}y_{k})(z^{T}s_{k})}/{y_{k}^{T}s_{k}}+2(z^{T}s_{k})^{2}{\|y_{k}\|^{2}}/{(y_{k}^{T}s_{k})^{2}}$
$\displaystyle\quad=\left(\|z\|-\left|{z^{T}s_{k}}/{y_{k}^{T}s_{k}}\right|\|y_{k}\|\right)^{2}+2\|z\|\left|{z^{T}s_{k}}/{y_{k}^{T}s_{k}}\right|\|y_{k}\|$
$\displaystyle\quad\quad-2{(z^{T}y_{k})(z^{T}s_{k})}/{y_{k}^{T}s_{k}}+\|y_{k}\|^{2}(z^{T}s_{k}/y_{k}^{T}s_{k})^{2}\geq
0.$ (21)

In the last inequality of equation (21), we use the Cauchy-Schwartz inequality
$\|z^{T}y\|\leq\|z\|\|y_{k}\|$ and its equality holds if only if $z=ty_{k}$.
When $z=ty_{k}$, from equation (21), we have
$z^{T}H_{k+1}z=t^{2}\|y_{k}\|^{2}=\|z\|^{2}>0$. When $z^{T}s_{k}=0$, from
equation (21), we also have $z^{T}H_{k+1}z=\|z\|^{2}>0$. Therefore, we
conclude that $H_{k+1}$ is a symmetric positive definite matrix when
$|s_{k}^{T}y_{k}|>\theta\|s_{k}\|^{2}$. From equation (20), We apparently
conclude that $H_{k+1}$ is a symmetric positive definite matrix since
$H_{k+1}=I$ when $|s_{k}^{T}y_{k}|\leq\theta\|s_{k}\|^{2}$.

(ii) It is not difficult to know that it exists at least $n-2$ linearly
independent vectors $z_{1},\,z_{2},\,\ldots,\,z_{n-2}$ such that
$z_{i}^{T}s_{k}=0,\,z_{i}^{T}y_{k}=0\,(i=1:(n-2))$ hold. That is to say,
matrix $H_{k+1}$ defined by equation (20) has at least $(n-2)$ linearly
independent eigenvectors whose corresponding eigenvalues are 1. We denote the
other two eigenvalues of $H_{k+1}$ as $\mu_{i}^{k+1}\,(i=1:2)$ and their
corresponding eigenvalues as $p_{1}$ and $p_{2}$, respectively. Then, from
equation (20), we know that the eigenvectors $p_{i}\,(i=1:2)$ can be
represented as $p_{i}=y_{k}+\beta_{i}s_{k}$ when $\mu_{i}^{k+1}\neq
1\,(i=1:2)$. From equation (20) and
$H_{k+1}p_{i}=\mu_{i}^{k+1}p_{i}\,(i=1:2)$, we have

$\displaystyle-\left(\mu_{i}^{k+1}+\beta_{i}\frac{s_{k}^{T}s_{k}}{s_{k}^{T}y_{k}}\right)y_{k}+\left(\frac{y_{k}^{T}y_{k}}{y_{k}^{T}s_{k}}+2\beta_{i}\frac{(y_{k}^{T}y_{k})(s_{k}^{T}s_{k})}{(y_{k}^{T}s_{k})^{2}}-\mu_{i}^{k+1}\beta_{i}\right)s_{k}=0,\;i=1:n.$
(22)

When $y_{k}=ts_{k}$, from equation (20), we have $H_{k+1}=I$. In this case, we
conclude that the eigenvalues of $H_{k+1}$ are greater than $1/2$. When
vectors $y_{k}$ and $s_{k}$ are linearly independent, from equation (22), we
have

$\displaystyle\mu_{i}^{k+1}+\beta_{i}{s_{k}^{T}s_{k}}/{s_{k}^{T}y_{k}}=0,$
$\displaystyle{y_{k}^{T}y_{k}}/{y_{k}^{T}s_{k}}+2\beta_{i}{(y_{k}^{T}y_{k})(s_{k}^{T}s_{k})}/{(y_{k}^{T}s_{k})^{2}}-\mu_{i}^{k+1}\beta_{i}=0,\;i=1:n.$

That is to say, $\mu_{i}^{k+1}\,(i=1:2)$ are the two solutions of the
following equation:

$\displaystyle\mu^{2}-2\mu(y_{k}^{T}y_{k})(s_{k}^{T}s_{k})/(s_{k}^{T}y_{k})^{2}+(y_{k}^{T}y_{k})(s_{k}^{T}s_{k})/(s_{k}^{T}y_{k})^{2}=0.$
(23)

Consequently, from equation (23), we obtain

$\displaystyle\mu_{1}^{k+1}+\mu_{2}^{k+1}=2(y_{k}^{T}y_{k})(s_{k}^{T}s_{k})/(s_{k}^{T}y_{k})^{2},\;\mu_{1}^{k+1}\mu_{2}^{k+1}=(y_{k}^{T}y_{k})(s_{k}^{T}s_{k})/(s_{k}^{T}y_{k})^{2}.$
(24)

From equation (24), it is not difficult to obtain

$\displaystyle{1}/{\mu_{1}^{k+1}}+{1}/{\mu_{2}^{k+1}}=2,\;\mu_{i}^{k+1}>0,\;i=1:2.$
(25)

Therefore, from equation (25), we conclude that
$\mu_{i}^{k+1}>\frac{1}{2}\,(i=1:2)$. Consequently, the eigenvalues of
$H_{k+1}$ are greater than 1/2. ∎

If $s_{k-1}$ is obtained from the explicit continuation method (18), we have
$Ps_{k-1}=s_{k-1}$ since $P^{2}=P$. By combining it with the L-BFGS updating
formula (20), the explicit continuation method (18)-(19) can be simplified by

$\displaystyle
d_{k}=\begin{cases}-p_{g_{k}},\;\text{if}\;|s_{k-1}^{T}y_{k-1}|\leq\theta\|s_{k-1}\|^{2},\\\
-p_{g_{k}}+\frac{y_{k-1}(s_{k-1}^{T}p_{g_{k}})+s_{k-1}(y_{k-1}^{T}p_{g_{k}})}{y_{k-1}^{T}s_{k-1}}-2\frac{\|y_{k-1}\|^{2}(s_{k-1}^{T}p_{g_{k}})}{(y_{k-1}^{T}s_{k-1})^{2}}s_{k-1},\,\text{otherwise},\end{cases}$
(26) $\displaystyle s_{k}=\frac{\Delta t_{k}}{1+\Delta
t_{k}}d_{k},\;x_{k+1}=x_{k}+s_{k},$ (27)

where $p_{g_{k}}=Pg_{k}=P\nabla f(x_{k})$ and $y_{k-1}=P\nabla
f(x_{k})-P\nabla f(x_{k-1})$. Thus, it does not need to store the matrix
$H_{k}$ in practical computation. Furthermore, it only requires three pairs of
the inner product of vector and one matrix-vector product ($p_{g_{k}}=Pg_{k}$)
to obtain the trial step $s_{k}$ and involves in $O((n-m)n)$ flops when we use
the QR decomposition or the singular value decomposition to obtain the
projection matrix $P$ in subsection 2.4.

###  2.4 The treatments of infeasible initial points and projection matrices

We need to compute the projected gradient $Pg_{k}$ at every iteration in the
updating formula (26). In order to reduce the computational complexity, we use
the QR factorization (pp.276-278, GV2013 ) to factor $A^{T}$ into a product of
an orthogonal matrix $Q\in\Re^{n\times n}$ and an upper triangular matrix
$R\in\Re^{n\times m}$:

$\displaystyle
A^{T}=QR=\begin{bmatrix}Q_{1}|Q_{2}\end{bmatrix}\begin{bmatrix}R_{1}\\\
0\end{bmatrix},$ (28)

where $Q_{1}=Q(1:n,\,1:m),\;Q_{2}=Q(1:n,\,m+1:n)$, $R_{1}=R(1:r,\,1:m)$ is
upper triangular and nonsingular. Then, from equations (10), (28), we simplify
the projection matrix $P$ as

$\displaystyle P=I-Q_{1}Q_{1}^{T}=Q_{2}Q_{2}^{T}.$ (29)

In practical computation, we adopt the different formulas of the projection
$P$ according to $m\leq n/2$ or $m>n/2$. Thus, we give the computational
formula of the projected gradient $Pg_{k}$ as follows:

$\displaystyle
Pg_{k}=\begin{cases}g_{k}-Q_{1}\left(Q_{1}^{T}g_{k}\right),\;\text{if}\;m\leq\frac{1}{2}n,\\\
Q_{2}\left(Q_{2}^{T}g_{k}\right),\;\text{otherwise}.\end{cases}$ (30)

For a real-world optimization problem (1), we probably meet the infeasible
initial point $x_{0}$. That is to say, the initial point can not satisfy the
constraint $Ax=b$. We handle this problem by solving the following projection
problem:

$\displaystyle\min_{x\in\Re^{n}}\;\left\|x-x_{0}\right\|^{2}\;\text{subject
to}\hskip 5.69054ptQ_{1}^{T}x=b_{r},$ (31)

where $b_{r}=\left(R_{1}R_{1}^{T}\right)^{-1}\left(R_{1}b\right)$. By using
the Lagrangian multiplier method and the QR factorization (28) of matrix
$A^{T}$ to solve problem (31), we obtain the initial feasible point
$x_{0}^{F}$ of problem (1) as follows:

$\displaystyle
x_{0}^{F}=x_{0}-Q_{1}\left(Q_{1}^{T}Q_{1}\right)^{-1}\left(Q_{1}^{T}x_{0}-b_{r}\right)=x_{0}-Q_{1}\left(Q_{1}^{T}x_{0}-b_{r}\right).$
(32)

For convenience, we set $x_{0}=x_{0}^{F}$ in line 4, Algorithm 1.

###  2.5 The trusty time-stepping scheme

Another issue is how to adaptively adjust the time-stepping size $\Delta
t_{k}$ at every iteration. We borrow the adjustment method of the trust-region
radius from the trust-region method due to its robust convergence and fast
local convergence CGT2000 . After the preprocess of the initial point $x_{0}$,
it is feasible. According to the structure-preserving property of the explicit
continuation method (18)-(20), $x_{k+1}$ will preserve the feasibility. That
is to say, $x_{k+1}$ satisfies $Ax_{k+1}=b$. Therefore, we use the objective
function $f(x)$ instead of the nonsmooth penalty function
$f(x)+\sigma\|Ax-b\|$ as the cost function.

When we use the trust-region updating strategy to adaptively adjust time-
stepping size $\Delta t_{k}$ Higham1999 , we need to construct a local
approximation model of the objective $f(x)$ around $x_{k}$. Here, we adopt the
following quadratic function as its approximation model:

$\displaystyle
q_{k}(x_{k}+s)=f(x_{k})+s^{T}g_{k}+\frac{1}{2}s^{T}H_{k}^{-1}s.$ (33)

In practical computation, we do not store the matrix $H_{k}$. Thus, we use the
explicit continuation method (18)-(20) and regard $(H_{k}P)(H_{k}P)^{+}\approx
I$ to simplify the quadratic model $q_{k}(x_{k}+s_{k})-q(x_{k})$ as follows:

$\displaystyle m_{k}(s_{k})=g_{k}^{T}s_{k}-\frac{0.5\Delta t_{k}}{1+\Delta
t_{k}}g_{k}^{T}s_{k}=\frac{1+0.5\Delta t_{k}}{1+\Delta
t_{k}}g_{k}^{T}s_{k}\approx q_{k}(x_{k}+s_{k})-q_{k}(x_{k}).$ (34)

where $g_{k}=\nabla f(x_{k})$. We enlarge or reduce the time-stepping size
$\Delta t_{k}$ at every iteration according to the following ratio:

$\displaystyle\rho_{k}=\frac{f(x_{k})-f(x_{k+1})}{m_{k}(0)-m_{k}(s_{k})}.$
(35)

A particular adjustment strategy is given as follows:

$\displaystyle\Delta t_{k+1}=\begin{cases}\gamma_{1}\Delta t_{k},&{if\hskip
2.84526pt0\leq\left|1-\rho_{k}\right|\leq\eta_{1},}\\\ \Delta t_{k},&{if\hskip
2.84526pt\eta_{1}<\left|1-\rho_{k}\right|<\eta_{2},}\\\ \gamma_{2}\Delta
t_{k},&{if\hskip 2.84526pt\left|1-\rho_{k}\right|\geq\eta_{2},}\end{cases}$
(36)

where the constants are selected as
$\eta_{1}=0.25,\;\gamma_{1}=2,\;\eta_{2}=0.75,\;\gamma_{2}=0.5$ according to
numerical experiments. When $\rho_{k}\geq\eta_{a}$, we accept the trial step
$s_{k}$ and let $x_{k+1}=x_{k}+s_{k}$, where $\eta_{a}$ is a small positive
number such as $\eta_{a}=1.0\times 10^{-6}$. Otherwise, we discard it and let
$x_{k+1}=x_{k}$.

According to the above discussions, we give the detailed implementation of the
explicit continuation method with the trusty time-stepping scheme for the
linearly equality-constrained optimization problem (1) in Algorithm 1.

Algorithm 1 The explicit continuation method with the trusty time-stepping
scheme for linearly constrained optimization (the Eptctr method)

0: the objective function $f(x)$, the linear constraint $Ax=b$, the initial
point $x_{0}$ (optional), the terminated parameter $\epsilon$ (optional).

0: the optimal approximation solution $x^{\ast}$.

1: Set $x_{0}=\text{ones}(n,\,1)$ and $\epsilon=10^{-6}$ as the default
values.

2: Initialize the parameters:
$\eta_{a}=10^{-6},\;\eta_{1}=0.25,\;\gamma_{1}=2,\;\eta_{2}=0.75,\;\gamma_{2}=0.5,\;\theta=10^{-6}$.

3: Factorize matrix $A^{T}$ with the QR factorization (28).

4: Compute  $x_{0}\leftarrow x_{0}-Q_{1}\left(Q_{1}^{T}x_{0}-b_{r}\right),$
such that $x_{0}$ satisfies the linear system of constraints $Ax=b$.

5: Set $k=0$. Evaluate $f_{0}=f(x_{0})$ and $g_{0}=\nabla f(x_{0})$.

6: Compute the projected gradient $p_{g_{0}}=Pg_{0}$ according to the formula
(30). Set $y_{-1}=0$ and $s_{-1}=0$.

7: Set $\Delta t_{0}=10^{-2}$.

8: while $\|p_{g_{k}}\|>\epsilon$ do

9: if $\left(|s_{k-1}^{T}y_{k-1}|>\theta s_{k-1}^{T}s_{k-1}\right)$ then

10: $s_{k}=-\frac{\Delta t_{k}}{1+\Delta
t_{k}}\left(p_{g_{k}}-\frac{y_{k-1}(s_{k-1}^{T}p_{g_{k}})+s_{k-1}(y_{k-1}^{T}p_{g_{k}})}{y_{k-1}^{T}s_{k-1}}+2\frac{\|y_{k-1}\|^{2}(s_{k-1}^{T}p_{g_{k}})}{(y_{k-1}^{T}s_{k-1})^{2}}s_{k-1}\right)$.

11: else

12: $s_{k}=-\frac{\Delta t_{k}}{1+\Delta t_{k}}p_{g_{k}}$.

13: end if

14: Compute $x_{k+1}=x_{k}+s_{k}$.

15: Evaluate $f_{k+1}=f(x_{k+1})$ and compute the ratio $\rho_{k}$ from
equations (34)-(35).

16: if $\rho_{k}\leq\eta_{a}$ then

17: Set
$x_{k+1}=x_{k},\;f_{k+1}=f_{k},\;p_{g_{k+1}}=p_{g_{k}},\;g_{k+1}=g_{k},\;y_{k}=y_{k-1}.$

18: else

19: Evaluate $g_{k+1}=\nabla f(x_{k+1})$.

20: Compute $p_{g_{k+1}}=Pg_{k+1}$ according to the formula (30). Set
$y_{k}=p_{g_{k+1}}-p_{g_{k}}$ and $s_{k}=x_{k+1}-x_{k}$.

21: end if

22: Adjust the time-stepping size $\Delta t_{k+1}$ based on the trust-region
updating scheme (36).

23: Set $k\leftarrow k+1$.

24: end while

##  3 Algorithm Analysis

In this section, we analyze the global convergence of the explicit
continuation method (18)-(19) with the trusty time-stepping scheme and the
L-BFGS updating formula (20) for the linearly equality-constrained
optimization problem (i.e. Algorithm 1). Firstly, we give a lower-bounded
estimate of $m_{k}(0)-m_{k}(s_{k})$ $(k=1,\,2,\,\ldots)$. This result is
similar to that of the trust-region method for the unconstrained optimization
problem Powell1975 . For simplicity, we assume that the rank of matrix $A$ is
full.

###### Lemma 2

Assume that the quadratic model $q_{k}(x)$ is defined by equation (34) and
$s_{k}$ is computed by the explicit continuation method (18)-(20). Then, we
have

$\displaystyle m_{k}(0)-m_{k}(s_{k})\geq\frac{\Delta t_{k}}{4(1+\Delta
t_{k})}\left\|p_{g_{k}}\right\|^{2},$ (37)

where $p_{g_{k}}=Pg_{k}=P\nabla f(x_{k})$ and the projection matrix $P$ is
defined by equation (10).

Proof. From equation (20) and Lemma 1, we know that $H_{k}$ is symmetric
positive definite and its eigenvalues are greater than 1/2. According to the
eigenvalue decomposition of $H_{k}$, we know that it exists an orthogonal
matrix $Q_{k}$ such that $H_{k}=Q_{k}^{T}S_{k}Q_{k}$ holds, where
$S_{k}=\text{diag}(\mu_{1}^{k},\,\ldots,\,\mu_{n}^{k})$ and
$\mu_{i}^{k}\,(i=1:n)$ are the eigenvalues of $H_{k}$. We denote the smallest
eigenvalue of $H_{k}$ is $\mu_{min}^{k}$. From the explicit continuation
method (18) and $P^{2}=P$, we know that $s_{k}=Ps_{k}$. By combining it with
the explicit continuation method (18) and the quadratic model (34), we have

$\displaystyle
m_{k}(0)-m_{k}(s_{k})\geq-\frac{1}{2}g_{k}^{T}s_{k}=-\frac{1}{2}g_{k}^{T}Ps_{k}=-\frac{1}{2}(Pg_{k})^{T}s_{k}=\frac{\Delta
t_{k}}{2(1+\Delta t_{k})}(Pg_{k})^{T}H_{k}(Pg_{k})$
$\displaystyle\quad=\frac{\Delta t_{k}}{2(1+\Delta
t_{k})}(Pg_{k})^{T}(Q_{k}^{T}S_{k}Q_{k})(Pg_{k})=\frac{\Delta
t_{k}}{2(1+\Delta t_{k})}(QPg_{k})^{T}S_{k}(QPg_{k})$
$\displaystyle\quad\geq\mu_{min}^{k}\frac{\Delta t_{k}}{2(1+\Delta
t_{k})}\|QPg_{k}\|^{2}\geq\frac{\Delta t_{k}}{4(1+\Delta
t_{k})}\|QPg_{k}\|^{2}=\frac{\Delta t_{k}}{4(1+\Delta t_{k})}\|Pg_{k}\|^{2}.$
(38)

In the first inequality in equation (38), we use the property $(1+0.5\Delta
t_{k})/(1+\Delta t_{k})\geq 0.5$ when $\Delta t_{k}\geq 0$. Consequently, we
prove the result (37). ∎

In order to prove that $p_{g_{k}}$ converges to zero when $k$ tends to
infinity, we need to estimate the lower bound of time-stepping sizes $\Delta
t_{k}\,(k=1,\,2,\,\ldots)$. We denote the constrained level set $S_{f}$ as

$\displaystyle S_{f}=\\{x:\;f(x)\leq f(x_{0}),\;Ax=b\\}.$ (39)

###### Lemma 3

Assume that $f:\;\Re^{n}\rightarrow\Re$ is continuously differentiable and its
gradient $g(x)$ satisfies the following Lipschitz continuity:

$\displaystyle\|g(x)-g(y)\|\leq L_{c}\|x-y\|,\;\forall x,\,y\in S_{f}.$ (40)

where $L_{c}$ is the Lipschitz constant. We suppose that the sequence
$\\{x_{k}\\}$ is generated by Algorithm 1. Then, there exists a positive
constant $\delta_{\Delta t}$ such that

$\displaystyle\Delta t_{k}\geq\gamma_{2}\delta_{\Delta t}$ (41)

holds for all $k=1,\,2,\,\dots$, where $\Delta t_{k}$ is adaptively adjusted
by the trust-region updating scheme (34)-(36).

Proof. From Lemma 1, we know that the eigenvalues of $H_{k}$ is greater than
1/2 and it has at least $n-2$ eigenvalues which equal 1. When
$|s_{k-1}^{T}y_{k-1}|\geq\theta\|s_{k-1}\|^{2}$, we denote the other two
eigenvalues of $H_{k}$ as $\mu_{1}^{k}$ and $\mu_{2}^{k}$. By substituting it
into equation (24), we obtain

$\displaystyle\mu_{1}^{k}\mu_{2}^{k}=\frac{\|y_{k-1}\|^{2}\|s_{k-1}\|^{2}}{(s_{k-1}^{T}y_{k-1})^{2}}\leq\frac{\|y_{k-1}\|^{2}\|s_{k-1}\|^{2}}{\theta^{2}\|s_{k-1}\|^{4}}=\frac{1}{\theta^{2}}\frac{\|y_{k-1}\|^{2}}{\|s_{k-1}\|^{2}}.$
(42)

From Lemma 2 and Algorithm 1, we know $f(x_{k})\leq
f(x_{0})\,(k=1,\,2,\,\ldots)$. From the explicit continuation method (18)-(20)
and Remark 3, we know that $Ax_{k}=Ax_{0}=b\,(k=1,\,2,\,\ldots)$. Thus, from
the Lipschitz continuity (40) of $g(x)$, we have

$\displaystyle\|y_{k-1}\|\leq\|P\|\|g(x_{k})-g(x_{k-1})\|\leq
L_{C}\|x_{k}-x_{k-1}\|=L_{C}\|s_{k-1}\|.$ (43)

By substituting it into equation (42) and using
$\mu_{i}^{k}>\frac{1}{2}\,(i=1,\,2)$, we obtain

$\displaystyle\frac{1}{2}<\mu_{i}^{k}<\frac{2L_{C}^{2}}{\theta^{2}},\;i=1,\,2.$
(44)

That is to say, the eigenvalues of $H_{k}$ are less than or equal to $M_{H}$,
where $M_{H}=\max\\{1,\,{2L_{C}^{2}}/{\theta^{2}}\\}$. According to the
eigenvalue decomposition theorem, we know that there exists an orthogonal
matrix $Q_{k}$ such that $H_{k}=Q_{k}^{T}S_{k}Q_{k}$ holds, where
$S_{k}=\text{diag}(\mu_{1}^{k},\,\ldots,\,\mu_{n}^{k})$ and
$\mu_{i}^{k}\,(i=1:n)$ are the eigenvalues of $H_{k}$. Consequently, we have

$\displaystyle\|H_{k}(Pg_{k})\|=\|(Q_{k}^{T}S_{k}Q_{k})Pg_{k}\|=\|S_{k}(Q_{k}Pg_{k})\|\leq
M_{H}\|Pg_{k}\|.$ (45)

From the first-order Taylor expansion, we have

$\displaystyle
f(x_{k}+s_{k})=f(x_{k})+\int_{0}^{1}s_{k}^{T}g(x_{k}+ts_{k})dt.$ (46)

Thus, from equations (34)-(37), (46) and the Lipschitz continuity (40) of
$g(x)$, we have

$\displaystyle\left|\rho_{k}-1\right|=\left|\frac{(f(x_{k})-f(x_{k}+s_{k}))-(m_{k}(0)-m_{k}(s_{k}))}{m_{k}(0)-m_{k}(s_{k})}\right|$
$\displaystyle\quad\leq\frac{1+\Delta t_{k}}{1+0.5\Delta
t_{k}}\frac{\left|\int_{0}^{1}s_{k}^{T}(g(x_{k}+ts_{k})-g(x_{k}))dt\right|}{\left|s_{k}^{T}g_{k}\right|}+\frac{0.5\Delta
t_{k}}{1+0.5\Delta t_{k}}$ $\displaystyle\quad\leq\frac{2L_{C}(1+\Delta
t_{k})}{\Delta t_{k}}\frac{\|s_{k}\|^{2}}{\|Pg_{k}\|^{2}}+\frac{0.5\Delta
t_{k}}{1+0.5\Delta t_{k}}.$ (47)

By substituting equation (18) and equation (45) into equation (47), we have

$\displaystyle\left|\rho_{k}-1\right|\leq\frac{2L_{C}\Delta t_{k}}{1+\Delta
t_{k}}\frac{\|PH_{k}(Pg_{k})\|^{2}}{\|Pg_{k}\|^{2}}+\frac{0.5\Delta
t_{k}}{1+0.5\Delta t_{k}}$ $\displaystyle\leq\frac{2L_{C}\Delta
t_{k}}{1+\Delta
t_{k}}\frac{\|P\|^{2}\|H_{k}(Pg_{k})\|^{2}}{\|Pg_{k}\|^{2}}+\frac{0.5\Delta
t_{k}}{1+0.5\Delta t_{k}}\leq\frac{(2L_{C}M_{H}^{2}+0.5)\Delta
t_{k}}{1+0.5\Delta t_{k}}.$ (48)

In the last inequality of equation (48), we use the property $\|P\|=1$. We
denote

$\displaystyle\delta_{\Delta
t}\triangleq\frac{\eta_{1}}{2L_{C}M_{H}^{2}+0.5}.$ (49)

Then, from equation (48)-(49), when $\Delta t_{k}\leq\delta_{\Delta t}$, it is
not difficult to verify

$\displaystyle\left|\rho_{k}-1\right|\leq(2L_{C}M_{H}^{2}+0.5)\Delta
t_{k}\leq\eta_{1}.$ (50)

We assume that $K$ is the first index such that $\Delta
t_{K}\leq\delta_{\Delta t}$ where $\delta_{\Delta t}$ is defined by equation
(49). Then, from equations (49)-(50), we know that $|\rho_{K}-1|\leq\eta_{1}$.
According to the time-stepping adjustment formula (36), $x_{K}+s_{K}$ will be
accepted and the time-stepping size $\Delta t_{K+1}$ will be enlarged.
Consequently, the time-stepping size $\Delta t_{k}$ holds $\Delta
t_{k}\geq\gamma_{2}\delta_{\Delta t}$ for all $k=1,\,2,\ldots$. ∎

By using the results of Lemma 2 and Lemma 3, we prove the global convergence
of Algorithm 1 for the linearly constrained optimization problem (1) as
follows.

###### Theorem 3.1

Assume that $f:\;\Re^{n}\rightarrow\Re$ is continuously differentiable and its
gradient $\nabla f(x)$ satisfies the Lipschitz continuity (40). Furthermore,
we suppose that $f(x)$ is lower bounded when $x\in S_{f}$, where the
constrained level set $S_{f}$ is defined by equation (39). The sequence
$\\{x_{k}\\}$ is generated by Algorithm 1. Then, we have

$\displaystyle\lim_{k\to\infty}\inf\|Pg_{k}\|=0,$ (51)

where $g_{k}=\nabla f(x_{k})$ and matrix $P$ is defined by equation (10).

Proof. According to Lemma 3 and Algorithm 1, we know that there exists an
infinite subsequence $\\{x_{k_{i}}\\}$ such that trial steps $s_{k_{i}}$ are
accepted, i.e., $\rho_{k_{i}}\geq\eta_{a},\,i=1,\,2,\ldots$. Otherwise, all
steps are rejected after a given iteration index, then the time-stepping size
will keep decreasing, which contradicts (41). Therefore, from equations (35)
and (37), we have

$\displaystyle
f(x_{0})-\lim_{k\to\infty}f(x_{k})=\sum_{k=0}^{\infty}(f(x_{k})-f(x_{k+1}))$
$\displaystyle\geq\eta_{a}\sum_{i=0}^{\infty}\left(m_{k_{i}}(0)-m_{k_{i}}(s_{k_{i}})\right)\geq\eta_{a}\sum_{i=0}^{\infty}\frac{\Delta
t_{k_{i}}}{4(\Delta t_{k_{i}}+1)}\|Pg_{k_{i}}\|.$ (52)

From the result (41) of Lemma 3, we know that $\Delta
t_{k}\geq\gamma_{2}\delta_{\Delta t}\,(k=1,\,2,\,\dots)$. By substituting it
into equation (52), we have

$\displaystyle
f(x_{0})-\lim_{k\to\infty}f(x_{k})\geq\eta_{a}\sum_{i=0}^{\infty}\frac{\gamma_{2}\delta_{\Delta
t}}{4(\gamma_{2}\delta_{\Delta t}+1)}\|Pg_{k_{i}}\|.$ (53)

Since $f(x)$ is lower bounded when $x\in S_{f}$ and the sequence
$\\{f(x_{k})\\}$ is monotonically decreasing, we have
$\lim_{k\to\infty}f(x_{k})=f^{\ast}$. By substituting it into equation (53),
we obtain the result (51). ∎

##  4 Numerical Experiments

In this section, some numerical experiments are executed to test the
performance of Algorithm 1 (the Eptctr method). The codes are executed by a
Dell G3 notebook with the Intel quad-core CPU and 20Gb memory. We compare
Eptctr with SQP (the built-in subroutine fmincon.m of the MATLAB2018a
environment) FP1963 ; Goldfarb1970 ; MATLAB ; NW1999 ; Wilson1963  Ptctr
LLS2020 ) for some large-scale linearly constrained-equality optimization
problems which are listed in Appendix A. SQP is the traditional and
representative optimization for the constrained optimization problem. Ptctr is
significantly better than SQP for linearly constrained optimization problems
according to the numerical results in LLS2020 . Therefore, we select these two
typical methods as the basis for comparison.

The termination conditions of the three compared methods are all set by

$\displaystyle\|\nabla_{x}L(x_{k},\,\lambda_{k})\|_{\infty}\leq 1.0\times
10^{-6},$ (54) $\displaystyle\|Ax_{k}-b\|_{\infty}\leq 1.0\times
10^{-6},\;k=1,\,2,\,\ldots,$ (55)

where the Lagrange function $L(x,\,\lambda)$ is defined by equation (4) and
$\lambda$ is defined by equation (8).

We test those ten problems with $n\approx 5000$. The numerical results are
arranged in Table 1 and illustrated in Figure 1. From Table 1, we find that
three methods can correctly solve those ten test problems and the consumed
time of Eptctr is significantly less than those of the other two methods for
every test problem, respectively. The consumed time of Eptctr is about one
tenth of that Ptctr or one fifteenth to 0.4 percent of that of SQP for the
test problem.

From those test data, we find that Eptctr works significantly better than the
other two methods, respectively. One of the reasons is that Eptctr only
involves three pairs of the inner product of two vectors and one matrix-vector
product ($p_{g_{k}}=Pg_{k}$) to obtain the trial step $s_{k}$ and involves
about $(n-m)n$ flops at every iteration. However, Ptctr needs to solve a
linear system of equations with an $n\times n$ symmetric definite coefficient
matrix and involves about $\frac{1}{3}n^{3}$ flops (p. 169, GV2013 ) at every
iteration. SQP needs to solve a linear system of equations with dimension
$(m+n)$ when it solves a quadratic programming subproblem at every iteration
(pp. 531-532, NW1999 ) and involves about $\frac{2}{3}(m+n)^{3}$ flops (p.
116, GV2013 ). Furthermore, Eptctr can save the storage space of an
$(n+m)\times(n+m)$ large-scale matrix, in comparison to SQP. The required
memory of Eptctr is about one fifth of that of SQP.

Table 1: Numerical results of test problems with $n\approx 5000$. Problems | Ptctr | Eptctr | SQP
---|---|---|---
| steps
---
(time)
$f(x^{\star})$ |  | Mem
---
/Gb
| steps
---
(time)
$f(x^{\star})$ |  | Mem
---
/Gb
| steps
---
(time)
$f(x^{\star})$ |  | Mem
---
/Gb
| Exam. 1
---
(n = 5000,
m = n/2)
| 11
---
(15.46)
3.64E+04 | 3.41 |  | 11
---
(1.94)
3.64E+04 | 0.51 |  | 2
---
(34.93)
3.64E+04 | 2.43
| Exam. 2
---
(n = 4800,
m = n/2)
| 17
---
(16.06)
5.78E+03 | 4.09 |  | 15
---
(1.58)
5.78E+03 | 0.31 |  | 14
---
(116.16)
5.78E+03 | 1.53
| Exam. 3
---
(n = 4800,
m = 2/3n)
| 12
---
(23.51)
2.86E+03 | 3.40 |  | 12
---
(2.11)
2.86E+03 | 0.54 |  | 3
---
(56.06)
2.86E+03 | 3.08
| Exam. 4
---
(n = 5000,
m = n/2)
| 11
---
(17.07)
493.79 | 3.41 |  | 11
---
(2.02)
493.79 | 0.51 |  | 8
---
(123.04)
493.79 | 2.43
| Exam. 5
---
(n = 5000,
m = n/2)
| 14
---
(16.79)
432.15 | 3.97 |  | 13
---
(2.04)
432.15 | 0.51 |  | 11
---
(178.41)
432.15 | 2.43
| Exam. 6
---
(n = 4800,
m = 2/3n)
| 13
---
(23.58)
2.06E+03 | 3.57 |  | 15
---
(2.17)
2.06E+03 | 0.54 |  | 11
---
(211.35)
2.06E+03 | 3.10
| Exam. 7
---
(n = 5000,
m = n/2)
| 10
---
(15.02)
5.94E+04 | 3.22 |  | 13
---
(2.07)
5.94E+04 | 0.51 |  | 20
---
(344.37)
5.94E+04 | 2.43
| Exam. 8
---
(n = 4800,
m = n/3)
| 38
---
(38.29)
776.88 | 7.36 |  | 133
---
(3.21)
-1.21E+04 | 0.31 |  | 28
---
(243.09)
784.94 | 1.53
| Exam. 9
---
(n = 5000,
m = n/2)
| 12
---
(22.59)
2.21E+05 | 3.59 |  | 8
---
(1.94)
2.21E+05 | 0.51 |  | 29
---
(490.36)
2.21E+05 | 2.43
| Exam. 10
---
(n = 4800,
m = n/3)
| 16
---
(13.51)
2.00 | 3.92 |  | 20
---
(1.39)
2.00 | 0.31 |  | 14
---
(109.59)
2.00 | 1.53
Figure 1: The consumed CPU time (s) of Ptctr, Eptctr and SQP for test problems
with $n\approx 5000$.

##  5 Conclusion and Future Work

In this paper, we give an explicit continuation method with the trusty time-
stepping scheme and the L-BFGS updating formula (Eptctr) for linearly
equality-constrained optimization problems. This method only involves three
pairs of the inner product of vector and one matrix-vector product
($p_{g_{k}}=Pg_{k}$) at every iteration, other than the traditional
optimization method such as SQP or the latest continuation method such as
Ptctr LLS2020 , which needs to solve a quadratic programming subproblem (SQP)
or a linear system of equations (Ptctr). Thus, Eptctr involves about $(n-m)n$
flops, Ptctr involves about $\frac{1}{3}n^{3}$ flops, and SQP involves about
$\frac{2}{3}(m+n)^{3}$ flops at every iteration. This means that Eptctr can
save much more computational time than SQP or Ptctr. Numerical results also
show that the consumed time of EPtctr is about one tenth of that Ptctr or one
fifteenth to 0.4 percent of that of SQP for the test problem with $n\approx
5000$. Furthermore, Eptctr can save the storage space of an $(n+m)\times(n+m)$
large-scale matrix, in comparison to SQP. The required memory of Eptctr is
about one fifth of that of SQP. Therefore, Eptctr is worth exploring further,
and we will extend it to the general nonlinear optimization problem in the
future.

## Acknowledgments

This work was supported in part by Grant 61876199 from National Natural
Science Foundation of China, Grant YBWL2011085 from Huawei Technologies Co.,
Ltd., and Grant YJCB2011003HI from the Innovation Research Program of Huawei
Technologies Co., Ltd.. The first author is grateful to Prof. Ya-xiang Yuan
and Prof. Li-zhi Liao for their suggestions.

##  Appendix A Test Problems

Example 1.

$\displaystyle\quad m$ $\displaystyle=n/2$
$\displaystyle\min_{x\in\Re^{n}}\;f(x)$
$\displaystyle=\sum_{k=1}^{n/2}\;\left(x_{2k-1}^{2}+10x_{2k}^{2}\right),\;\text{subject
to}\;x_{2i-1}+x_{2i}=4,\;i=1,\,2,\ldots,\,m.$

This problem is extended from the problem of Kim2010 . We assume that the
feasible initial point is $(2,\,2,\,\ldots,\,2,\,2)$.

Example 2.

$\displaystyle\quad m$ $\displaystyle=n/3$
$\displaystyle\min_{x\in\Re^{n}}\;f(x)$
$\displaystyle=\sum_{k=1}^{n/2}\;\left(\left(x_{2k-1}-2\right)^{2}+2\left(x_{2k}-1\right)^{4}\right)-5,\;\text{subject
to}\;x_{3i-2}+4x_{3i-1}+2x_{3i}=3,\;i=1,\,2,\ldots,\,n/3.$

We assume that the infeasible initial point is
$(-0.5,\,1.5,\,1,\,0,\,\ldots,\,0,\,0)$.

Example 3.

$\displaystyle\quad m$ $\displaystyle=(2/3)n$
$\displaystyle\min_{x\in\Re^{n}}\;f(x)$
$\displaystyle=\sum_{k=1}^{n}\;x_{k}^{2},\;\text{subject
to}\;x_{3i-2}+2x_{3i-1}+x_{3i}=1,\;2x_{3i-2}-x_{3i-1}-3x_{3i}=4,\;i=1,\,2,\ldots,\,n/3.$

This problem is extended from the problem of Osborne2016 . The infeasible
initial point is $(1,\,0.5,\,-1,\,\ldots,\,1,\,0.5,\,-1)$.

Example 4.

$\displaystyle\quad m$ $\displaystyle=n/2$
$\displaystyle\min_{x\in\Re^{n}}\;f(x)$
$\displaystyle=\sum_{k=1}^{n/2}\;\left(x_{2k-1}^{2}+x_{2k}^{6}\right)-1,\;\text{subject
to}\;x_{2i-1}+x_{2i}=1,\;i=1,\,2,\,\ldots,\,n/2.$

This problem is modified from the problem of MAK2019 . We assume that the
infeasible initial point is $(1,\,1,\,\ldots,\,1)$.

Example 5.

$\displaystyle\quad m$ $\displaystyle=n/2$
$\displaystyle\min_{x\in\Re^{n}}\;f(x)$
$\displaystyle=\sum_{k=1}^{n/2}\;\left(\left(x_{2k-1}-2\right)^{4}+2\left(x_{2k}-1\right)^{6}\right)-5,\;\text{subject
to}\;x_{2i-1}+4x_{2i}=3,\;i=1,\,2,\,\ldots,\,m.$

We assume that the feasible initial point is
$(-1,\,1,\,-1,\,1,\,\ldots,\,-1,\,1)$.

Example 6.

$\displaystyle\quad m$ $\displaystyle=(2/3)n$
$\displaystyle\min_{x\in\Re^{n}}\;f(x)$
$\displaystyle=\sum_{k=1}^{n/3}\;\left(x_{3k-2}^{2}+x_{3k-1}^{4}+x_{3k}^{6}\right),$
subject to
$\displaystyle{x_{3i-2}}+2{x_{3i-1}}+{x_{3i}}=1,\;2{x_{3i-2}}-{x_{3i-1}}-3{x_{3i}}=4,\;i=1,\,2,\,\ldots,\,m/2.$

This problem is extended from the problem of Osborne2016 . We assume that the
infeasible initial point is $(2,\,0,\,\ldots,\,0)$.

Example 7.

$\displaystyle\quad m$ $\displaystyle=n/2$
$\displaystyle\min_{x\in\Re^{n}}\;f(x)$
$\displaystyle=\sum_{k=1}^{n/2}\;\left(x_{2k-1}^{4}+3x_{2k}^{2}\right),\;\text{subject
to}\;x_{2i-1}+x_{2i}=4,\;i=1,\,2,\,\ldots,\,n/2.$

This problem is extended from the problem of Carlberg2009 . We assume that the
infeasible initial point is $(2,\,2,\,0,\,\ldots,\,0,\,0)$.

Example 8.

$\displaystyle\quad m$ $\displaystyle=n/3$
$\displaystyle\min_{x\in\Re{{}^{n}}}\;f(x)$
$\displaystyle=\sum_{k=1}^{n/3}\;\left(x_{3k-2}^{2}+x_{3k-2}^{2}\,x_{3k}^{2}+2x_{3k-2}\,x_{3k-1}+x_{3k-1}^{4}+8x_{3k-1}\right),$
subject to $\displaystyle
2x_{3i-2}+5x_{3i-1}+x_{3i}=3,\;i=1,\,2,\,\ldots,\,m.$

We assume that the infeasible initial point is $(1.5,\,0,\,0,\,\ldots,\,0)$.

Example 9.

$\displaystyle\quad m$ $\displaystyle=n/2$
$\displaystyle\min_{x\in\Re^{n}}\;f(x)$
$\displaystyle=\sum_{k=1}^{n/2}\;\left(x_{2k-1}^{4}+10x_{2k}^{6}\right),\;\text{subject
to}\;x_{2i-1}+x_{2i}=4,\;i=1,\,2,\,\ldots,\,m.$

This problem is extended from the problem of Kim2010 . We assume that the
feasible initial point is $(2,\,2,\,\ldots,\,2,\,2)$.

Example 10.

$\displaystyle\quad m$ $\displaystyle=n/3$
$\displaystyle\min_{x\in\Re^{n}}\;f(x)$
$\displaystyle=\sum_{k=1}^{n/3}\;\left(x_{3k-2}^{8}+x_{3k-1}^{6}+x_{3k}^{2}\,\right),\;\text{subject
to}\;x_{3i-2}+2x_{3i-1}+2x_{3i}=1,\;i=1,\,2,\,\ldots,\,m.$

This problem is modified from the problem of Yamashita1980 . The feasible
initial point is $(1,\,0,\,0,\,\ldots,\,1,\,0,\,0)$.

## References

  * (1) E. L. Allgower and K. Georg, _Introduction to Numerical Continuation Methods_ , SIAM, Philadelphia, PA, 2003.
  * (2) U. M. Ascher and L. R. Petzold, _Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations_ , SIAM, Philadelphia, PA, 1998.
  * (3) D. P. Bertsekas, _Nonlinear Programming (3rd Edition)_ , Tsinghua University Press, 2018.
  * (4) A. A. Brown and M. C. Bartholomew-Biggs, _ODE versus SQP methods for constrained optimization_ , Journal of Optimization and Theory Applications, 62 (3): 371-386, 1989.
  * (5) K. E. Brenan, S. L. Campbell and L. R. Petzold, _Numerical solution of initial-value problems in differential-algebraic equations_ , SIAM, Philadelphia, PA, 1996.
  * (6) R. Byrd, J. Nocedal and Y. X. Yuan, _Global convergence of a class of quasi-Newton methods on convex problems_ , SIAM Journal of Numerical Analysis, 24: 1171-1189, 1987.
  * (7) K. Carlberg, _Lecture notes of constrained optimization_ , https://www.sandia.gov/~ktcarlb/opt_class/OPT_Lecture3.pdf, 2009.
  * (8) F. Caballero, L. Merino, J. Ferruz and A. Ollero, _Vision-based odometry and SLAM for medium and high altitude flying UAVs_ , Journal of Intelligent and Robotic Systems, 54 (1-3): 137-161, 2009.
  * (9) T. S. Coffey, C. T. Kelley and D. E. Keyes, _Pseudotransient continuation and differential-algebraic equations_ , SIAM Journal on Scientific Computing, 25: 553-569, 2003.
  * (10) A. R. Conn, N. Gould and Ph. L. Toint, _Trust-Region Methods_ , SIAM, Philadelphia, USA, 2000.
  * (11) A. V. Fiacco and G. P. McCormick, _Nonlinear programming: Sequential Unconstrained Minimization Techniques_ , SIAM, 1990.
  * (12) R. Fletcher and M. J. D. Powell, _A rapidly convergent descent method for minimization_ , Computer Journal, 6: 163-168, 1963.
  * (13) B. S. Goh, _Approximate greatest descent methods for optimization with equality constraints_ , Journal of Optimization Theory and Applications 148 (3): 505-527, 2011.
  * (14) D. Goldfarb, _A family of variable metric updates derived by variational means_ , Mathematics of Computing, 24: 23-26, 1970.
  * (15) G. H. Golub and C. F. Van Loan, _Matrix Computations_ , 4th ed., The Johns Hopkins University Press, 2013.
  * (16) E. Hairer, C. Lubich and G. Wanner, _Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations_ , 2nd ed., Springer, Berlin, 2006.
  * (17) U. Helmke and J. B. Moore, _Optimization and Dynamical Systems_ , 2nd ed., Springer-Verlag, London, 1996.
  * (18) D. J. Higham, _Trust region algorithms and timestep selection_ , SIAM Journal on Numerical Analysis, 37: 194-210, 1999.
  * (19) C. T. Kelley, L.-Z. Liao, L. Qi, M. T. Chu, J. P. Reese and C. Winton, _Projected Pseudotransient Continuation_ , SIAM Journal on Numerical Analysis, 46: 3071-3083, 2008.
  * (20) D. G. Liu and J. G. Fei, _Digital Simulation Algorithms for Dynamic Systems_ (in Chinese), Science Press, Beijing, 2000.
  * (21) S.-T. Liu and X.-L. Luo, _A method based on Rayleigh quotient gradient flow for extreme and interior eigenvalue problems_ , Linear Algebra and its Applications, 432 (7): 1851-1863, 2010.
  * (22) X.-L. Luo, _A dynamical method of DAEs for the smallest eigenvalue problem_ , Journal of Computational Science, 3 (3): 113-119, 2012.
  * (23) X.-L. Luo, C. T. Kelley, L.-Z. Liao and H.-W. Tam, _Combining trust-region techniques and Rosenbrock methods to compute stationary points_ , Journal of Optimization Theory and Applications, 140 (2): 265-286, 2009.
  * (24) X.-L. Luo, J.-R. Lin and W.-L. Wu, _A prediction-correction dynamic method for large-scale generalized eigenvalue problems_ , Abstract and Applied Analysis, Article ID 845459, 1-8, http://dx.doi.org/10.1155/2013/845459, 2013.
  * (25) X.-L. Luo, J.-H. Lv and G. Sun, _Continuation method with the trusty time-stepping scheme for linearly constrained optimization with noisy data_ , published online in http://arxiv.org/abs/2005.05965 or http://doi.org/10.1007/s11081-020-09590-z, Optimization and Engineering, Accepted, January 10, 2021.
  * (26) X.-L. Luo, H. Xiao and J.-H. Lv, _Continuation Newton methods with the residual trust-region time-stepping scheme for nonlinear equations_ , June 2020, arXiv preprint, http://arxiv.org/abs/2006.02634.
  * (27) X.-L. Luo and Y.Y. Yao, _Primal-dual path-following methods and the trust-region updating strategy for linear programming with noisy data_ , June 2020, arXiv preprint available at http://arxiv.org/abs/2006.07568, minor revision resubmitted to Journal of Computational Mathematics, January 16, 2021.
  * (28) M.-W. Mak, _Lecture notes of constrained optimization and support vector machines_ , http://www.eie.polyu.edu.hk/~mwmak/EIE6207/ContOpt-SVM-beamer.pdf, 2019.
  * (29) MATLAB 9.4.0 (R2018a), The MathWorks Inc., http://www.mathworks.com, 2018.
  * (30) J. Nocedal and S. J. Wright, _Numerical Optimization_ , Springer-Verlag, 1999.
  * (31) N. H. Kim, _Leture notes of constrained optimization_ , https://mae.ufl.edu/nkim/eas6939/ConstrainedOpt.pdf, 2010.
  * (32) M. J. Obsborne, _Mathematical methods for economic theory_ , https://mjo.osborne.economics.utoronto.ca/index.php/tutorial/index/1/mem, 2016.
  * (33) P.-Q. Pan, _New ODE methods for equality constrained optimization (2): algorithms_ , Journal of Computational Mathematics, 10 (2): 129-146, 1992.
  * (34) M. J. D. Powell, _Convergence properties of a class of minimization algorithms_ , in: O.L. Mangasarian, R. R. Meyer and S. M. Robinson, eds., _Nonlinear Programming 2_ , Academic Press, New York, 1-27, 1975.
  * (35) J. Schropp, _A dynamical systems approach to constrained minimization_ , Numerical Functional Analysis and Optimization, 21 (3-4): 537-551, 2000.
  * (36) J. Schropp, _One- and multistep discretizations of index 2 differential algebraic systems and their use in optimization_ , Journal of Computational and Applied Mathematics, 150: 375-396, 2003.
  * (37) L. F. Shampine, I. Gladwell and S. Thompson, _Solving ODEs with MATLAB_ , Cambridge University Press, Cambridge, 2003.
  * (38) W. Y. Sun and Y. X. Yuan, _Optimization Theory and Methods: Nonlinear Programming_ , Springer, New York, 2006.
  * (39) K. Tanabe, _A geometric method in nonlinear programming_ , Journal of Optimization Theory and Applications, 30 (2): 181-210, 1980.
  * (40) T. E. Simos, _New open modified Newton Cotes type formulae as multilayer symplectic integrators_ , Applied Mathematical Modelling 37: 1983-1991, 2013\.
  * (41) N. Ullah, J. Sabi and A. Shah, _A derivative-free scaled memoryless BFGS method for solving a system of monotone nonlinear equations_ , Submitted to Numerical Linear Algebra with Applications, October 2020.
  * (42) R. B. Wilson, _A Simplicial Method for Convex Programming_ , Phd thesis, Harvard University, 1963.
  * (43) H. Yamashita, _A differential equation approach to nonlinear programming_ , Mathematical Programming, 18: 155-168, https://doi.org/10.1007/BF01588311, 1980\.
  * (44) Y. Yuan, _Recent advances in trust region algorithms_ , Mathematical Programming, 151: 249-281, 2015.