††thanks: Corresponding author<EMAIL_ADDRESS>
# Formation and stability of vortex solitons in nematic liquid crystals

Pawel S. Jung  CREOL, The College of Optics and Photonics, University of
Central Florida, Orlando, Florida 32816, USA  Faculty of Physics, Warsaw
University of Technology, Warsaw, Poland  Yana V. Izdebskaya  Laser Physics
Centre, Research School of Physics, The Australian National University,
Canberra ACT 0200, Australia  Nonlinear Physics Centre, Research School of
Physics, The Australian National University, Canberra ACT 0200, Australia
Vladlen G. Shvedov  Laser Physics Centre, Research School of Physics, The
Australian National University, Canberra ACT 0200, Australia  Demetrios N.
Christodoulides  CREOL, The College of Optics and Photonics, University of
Central Florida, Orlando, Florida 32816, USA  Wieslaw Krolikowski  Laser
Physics Centre, Research School of Physics, The Australian National
University, Canberra ACT 0200, Australia  Science Program, Texas A & M
University at Qatar, Doha, Qatar

###### Abstract

We study the propagation dynamics of bright optical vortex solitons in nematic
liquid crystals with a nonlocal reorientational nonlinear response. We
investigate the role of optical birefringence on the stability of these
solitons. In agreement with recent experimental observations, we show that the
birefringence-induced astigmatism can eventually destabilize these vortex
solitons. However, for low and moderate birefringence, vortex solitons can
propagate stably over experimentally relevant distances.

Vortex solitons represent non-diffracting singular bright rings propagating in
nonlinear media by properly balancing diffraction and self-focusing nonlinear
effects [1, 2, 3]. They are rather elusive objects as they tend to break-up
during evolution. In typical homogenous local nonlinear media such as those
exhibiting Kerr-like and saturable nonlinearities, vortex solitons are subject
to azimuthal instabilities that can break them into individual bright solitons
[4, 5]. A few approaches have been proposed to arrest this instability. They
involve, for instance, use of vector solitons formed by the incoherent
superposition of a vortex component and a fundamental beam [6, 7] or two
vortices with opposite charges [8]. In addition, systems with nonlinear gain
and loss may support stable vortex solitons [9].

As far as single-component vortex solitons in passive nonlinear media are
concerned, their azimuthal instability can be prevented by the simultaneous
action of competing nonlinearities. In this nonlinear model, a substantial
defocusing contribution to the otherwise focusing nonlinearity may stabilize
the vortex soliton as has been predicted for the so-called cubic-quintic model
[10]. This model has been later employed to theoretically demonstrate
stabilization of even higher dimensional vortices [11].

Figure 1:  Schematic representation of beam propagation in a nematic liquid
crystal cell. The uniform molecular orientation defines the ordinary (x-y) and
extraordinary planes (y-z). An arbitrarily polarized input beam splits into an
ordinary (freely diffracting) and an extraordinary (self-focusing) component.

On the other hand, it appears that the spatially nonlocal nonlinearity may
stabilize vortex solitons. Nonlocality creates a broad refractive index
structure that serves as a waveguide (or trapping potential) and may prevent a
soliton from breaking up [12, 13]. The spatial nonlocality of nonlinearity is
commonly associated with certain transport processes such as heat transfer in
thermal media [14, 15], diffusion of atoms or molecules [16]. It is also
natural for media with long-range inter-molecular or atomic interactions, such
as in cold atomic gases [17] or liquid crystals (LC) [18]. The latter is
particularly important because their large noninstantaneous orientational and
spatially nonlocal nonlinearity and strong electro-optical response allows one
to realize various types of solitons, the so-called nematicons.

The nonlinear response of LC is caused by the light-induced reorientation of
their constituent molecules and subsequent increase of the extraordinary
refractive index. To this end, the molecules of liquid crystals must be
externally pre-oriented. It is commonly done by placing the liquid crystal
between two planar boundaries, which enforce initial molecular orientation by
either applying an external AC voltage [19] or rubbing the surfaces, which
pre-orients the molecules. Alternatively, one may use an external magnetic
field, which also orients molecules as they exhibit strong diamagnetism [20].
Until recently, the vortex-like solitons, which were observed in liquid
crystal cells, were vector vortex solitons, i.e., structures formed by a co-
propagating vortex beam with another bell-like beam which served as a guide
for the vortex component [7, 21].

However, formation and stable propagation of real (single beam, or scalar)
vortex solitons have been reported very recently in two experimental studies.
Izdebskaya et al. [20] employed a magnetic field for controlling the pre-
orientation and the degree of nonlocality [22], while Laudyn et al. [23] used
the mechanically conditioned liquid crystal cell’s surfaces to enforce the
proper initial molecular orientation. The latter work demonstrated that the
vortices’ stability critically depends on the optical anisotropy of the liquid
crystal. It was shown that liquid crystals with low birefringence do allow the
formation and propagation of vortex solitons. High birefringence introduces
strong astigmatism, which, subsequently, destabilizes the vortex soliton,
which quickly disintegrates during propagation.

In this work, we address the stability of vortex solitons in nematic liquid
crystals with a molecular orientation-based nonlinearity.In particular, we
show how this stability depends on the birefringence induced astigmatism. Our
results provide a detailed explanation of vortex instability in nematic liquid
crystals and hence complement recent experimental and earlier theoretical
works on this subject.

We consider the propagation of optical beams in the bulk of a nematic liquid
crystal having a uniform molecular initial orientation, as depicted in Fig.1.
The initial orientation can be achieved using an external low-frequency AC
electric or DC magnetic field. In general, the LC represents a uniaxial
optical medium. Therefore an arbitrary polarized optical beam that is incident
onto the LC cell will split into an ordinary and an extraordinary component.
The former propagates linearly, experiencing diffraction. The latter may
induce reorientation of molecules, thus increasing the effective refractive
index if its intensity is high, and, consequently, could experience self-
focusing in which case it will propagate as a spatial optical soliton. For
simplicity and without any loss of generality, in what follows, we will use a
simplified description of beam evolution in a LC, where we treat the medium as
a uniaxial optical crystal. Assuming that the input beam is extraordinarily
polarized, the spatial evolution of the electric field’s amplitude $E(x,z)$,
propagating predominantly along the y-axis (neglecting absorption), is
described by the following wave equation [24, 25]:

$\displaystyle 2ik_{o}n(\theta_{o})\left(\frac{\partial E}{\partial
y}+\tan\delta(\theta)\frac{\partial E}{\partial
z}\right)+\frac{\partial^{2}E}{\partial
x^{2}}+D_{z}\left(\theta\right)\frac{\partial^{2}E}{\partial z^{2}}$ (1)
$\displaystyle+k_{o}^{2}\left(n^{2}(\theta)-n^{2}(\theta_{o})\right)E=0,$

where $k_{o}$=2$\pi/\lambda_{o}$, $\delta\left(\theta\right)$ is the walk-off
angle calculated along the beam axis, $\theta_{o}$ is the initial molecular
orientation (in the absence of light),
$D_{z}=\cos^{2}\theta+\gamma^{2}\sin^{2}\theta$ is the diffraction coefficient
across $z$ and
$n\left(\theta\right)=\left(\cos^{2}\theta/n_{o}^{2}+\sin^{2}\theta/n^{2}_{e}\right)^{-1/2}$
is an effective index of refraction for the y-polarized (i.e., extraordinary)
light. Here $n_{o}$ and $n_{e}$ represents ordinary and extraordinary
refractive indices, respectively and $\gamma^{2}=n_{e}^{2}/n_{o}^{2}$. The
index $n(\theta)$ and walk-off angle $\delta$ [26] depend on the local
orientation of molecules $\theta$ which is governed by the following relation
[27]:

$\frac{\partial^{2}\theta}{\partial z^{2}}+\frac{\partial^{2}\theta}{\partial
x^{2}}+\frac{\Delta\varepsilon\varepsilon_{o}}{2K}|{E}|^{2}\sin 2\theta=0.$
(2)

In Eq.(2) $K$ is an effective elastic constant and
$\Delta\varepsilon=n_{e}^{2}-n_{o}^{2}$ is the dielectric [28] anisotropy of
the liquid crystal.

In what follows, we will neglect the walk-off effect (as defined by $\delta$)
as it does not play a role in the vortex soliton instability. A crucial aspect
associated with Eq.(1) is the presence of anisotropy in the Laplacian,
reflected by an effective an diffraction coefficient $D_{z}(\theta)$. This
term reflects an anisotropy in the diffraction along the extraordinary and
ordinary axes. In effect, it is a source of astigmatism which affects beam
propagation, and as we will show later, it is directly responsible for the
instability of the vortex soliton. The diffraction coefficient is a function
of the LC’s optical anisotropy (if $D_{z}=1$ the medium is isotropic).
Fig.2(a) depicts the diffraction coefficient as a function of the initial
molecular orientation ($\theta_{o}$) and an optical anisotropy ($\Delta
n=n_{e}-n_{o}$). We used values ranging from $\Delta n=0.1$ to $\Delta n=0.4$
, which correspond to widely available types of nematic liquid crystals. The
lack of cylindrical symmetry in Eq.(1) directly affects the optical beam’s
propagation even in the linear regime.

Figure 2:  (a) The effective diffraction coefficient along
$z^{\prime}$-direction for few values of optical anisotropy. (b-c)
Illustrating the effect of optical anisotropy on propagation of vortex beams
with charge $m=1$ and waist $w_{o}=6{\mu}$m in liquid crystal ($\lambda=0.8$
$\mu m$) (b,c) linear regime (low intensity), (d,e) nonlinear regime (high
intensity). $D_{z}=1$ in (b,d) and $D_{z}=1.6$ in (c,e).

This effect is illustrated in Figs.2(b-c), which depict the propagation of a
single charge optical vortex beam with (b) and without (c) anisotropy. In the
latter case, the vortex retains its the structure while gradually diffracting.
However, in the anisotropic medium, ensued astigmatism completely breaks down
the vortex beam. To better understand the effect of optical anisotropy and
initial molecular orientation of LC on the stability of vortex solitons, we
conducted a series of numerical simulations by varying the strength of
anisotropy ($\Delta n$) and initial molecular orientation ($\theta_{o}$). In
all subsequent numerical simulations, we use a Laguerre-Gaussian-like beam as
an initial condition. The beam size was kept at 6$\mu$m. The input power was
chosen such that the beam would constitute an exact vortex soliton solution in
the isotropic case ($D_{z}=1$). Consequently, the input power was adjusted to
account for the reorientation nonlinearity dependence on $\theta_{0}$,
reaching its maximum near $\theta_{0}\approx 45^{\circ}$. In Fig.2 (e), we
show the nonlinear propagation of a single charge vortex soliton in the
presence of strong anisotropy. Clearly, the soliton becomes quickly unstable.
For comparison, in Fig.2(d), we show the same soliton dynamics without an
anisotropy, which as expected is stable.

Figure 3:  Illustrating the effect of optical anisotropy on the vortex
solitons structural instability having charge 1 in nematic liquid crystals
with various degrees of anisotropicity. (a) Beam intensity ($|E|^{2}$) and
phase ($\varphi$) distribution at the output z=0.5cm. (b) The asymmetry
coefficient Q as a function of distance for two values of anisotropy.

In Fig.3, we summarize the results of numerical simulations of a single charge
vortex soliton propagating in the presence of an optical anisotropy. These
plots depict the output intensity distribution of an initially exact vortex
soliton, after 0.5 cm propagation distance in a LC (Fig.3(b)). Each row
corresponds to a fixed initial molecular orientation while each column to a
specific value of optical anisotropy, ranging from low to high. The vortex
soliton instability is quantified in Fig.3(b) which depicts the asymmetry
coefficient
$Q=\left(\iint_{z,x=0}^{z,x=+\infty}|E|^{2}dxdz\right)/\iint_{z,x=-\infty}^{z,x=+\infty}|E|^{2}dxdz$,
for two values of anisotropy. The plots confirm that the larger the initial
molecular orientation angle is, the faster the astigmatism-induced vortex
soliton instability develops, even for a small linear anisotropy. On the other
hand, one can achieve a relatively stable vortex soliton even in highly
anisotropic crystals, provided the initial molecular orientation angle is
sufficiently small. For small orientation angles the nonlinearity is and hance
a higher input power is required to maintain the soliton. However, as the plot
shows, even for $\theta_{0}=40^{\circ}$, the charge one vortex soliton may, in
principle, propagate over 5$mm$ or more.

Figure 4:  Astigmatism-induced structural instability of a high charge vortex
solitons in a nematic liquid crystal with a relatively large anisotropy
($\Delta n=0.4$) . In each case, graphs depict the vortex solitons intensity
($|E|^{2}$) and phase ($\varphi$) profiles at the same output distance
$z=0.5~{}$cm.

The physics behind the instability of these vortex solitons can be explained
by using a waveguide approach. For the vortex soliton to be stable, it must be
a mode of a self-induced (via nonlinearity) waveguide. In the isotropic case,
the waveguide has a cylindrical symmetry, and the vortex beam can be
represented as in the quadrature superposition of two first excited (dipole-
like) modes. Degeneracy ensures that both modes propagate with the same phase
velocities, thus preserving the stationary singular phase structure. In our
situation, the mode degeneracy is lifted because of the difference in the
diffraction coefficients along the principal x and z’ directions.
Consequently, the constituent in-quadrature dipole modes travel with different
phase velocities, thus destroying the initial vortex phase and intensity
structures. The larger the $D_{z}(\theta$), the quicker the vortex
destabilizes and ultimately breaks up. On the other hand, it is well known
that vortices do not have to exhibit circular symmetry. The so-called non-
canonical vortices can display an elliptical intensity distribution and an
angularly nonuniform phase gradient [29]. Since the LC anisotropy breaks the
circular symmetry, one might expect that it could, perhaps, support elliptical
vortex solitons. Such elliptical vortex solitons were shown to exist in Bose-
Einstein condensates, trapped in anisotropic harmonic potentials [30]. In a
liquid crystal, the nonlocal self-induced refractive index change plays the
role of this ”external” potential. However, because of the radial symmetry of
the nonlocal response in LC, which tends to weaken the anisotropic nonlinear
index change, the elliptical vortex soliton cannot maintain its structure and
breaks down, as we confirmed in additional simulations. The nonexistence of
elliptic vortex solitons has been also confirmed in earlier studies of
nonlocal media with anisotropic linear or nonlinear properties [31, 32].
However, a co-propagating vortex beam with a fundamental bell-like beam should
stabilize the former even for a strong anisotropy. In such a case, the
fundamental beam plays the role of an external anisotropic potential. Together
with the nonlinear contribution from the vortex itself, it may restore the
degeneracy of the two principal orthogonal dipole-like modes in the self-
induced waveguide, that supporting the vortex’s stable soliton propagation.

For completeness, we also studied the behavior of high charge vortex solitons.
A summary of numerical simulations for vortices with charges 1-4 is depicted
in Fig.4. It is clear that these solitons quickly lose stability and breakup,
first into their constituent fundamental components followed by a complete
disintegration. This happens even for small initial molecular orientation
angles. The breakup mechanism is similar to that seen for charge 1 vortices -
lack of degeneracy between constituent modes forming high-charge vortices. In
addition, there is a natural tendency for high charge vortices to split into
fundamental vortices under small perturbations, even in the linear regime.
Notice that even though the vortex soliton ultimately breaks down during
propagation, its remnants stay close to the propagation axis. In this case the
strong spatial LC nonlocality induces an overall waveguide-like structure that
traps the original soliton and its remnants during evolution.

In conclusion, we numerically studied the structural stability of vortex
solitons in nematic liquid crystals with a nonlocal reorientation nonlinearity
response. We showed that the anisotropy-induced astigmatism can destabilize
the solitons and can lead to a breakup. However, the process can be slow for
low anisotropies and small initial orientation angles, thus allowing almost
stable propagation of red charge vortex solitons over an experimentally
relevant distance. On the other hand, vortex solitons with higher charges
quickly lose their stability and disintegrate during propagation even for
moderate degrees of astigmatism. Our results explain recent experimental works
on vortex solitons in low and high birefringent nematic liquid crystals.

Founding: The Qatar National Research Fund (grant NPRP12S-0205-190047), the
Polish National Science Center (Contract No. UMO-2013/08/A/ST3/00708) and the
Polish Ministry of Science and Higher Education (1654/MOB/V/2017/0),the ONR
MURI (N00014-20-1-2789),the MPS Simons collaboration (Simons grant 733682).

## References

  * Kruglov and Vlasov [1985] V. Kruglov and R. Vlasov, Spiral self-trapping propagation of optical beams in media with cubic nonlinearity, Physics Letters A 111, 401 (1985).
  * Denschlag  _et al._ [2000] J. Denschlag, J. E. Simsarian, D. L. Feder, C. W. Clark, L. A. Collins, J. Cubizolles, L. Deng, E. W. Hagley, K. Helmerson, W. P. Reinhardt, _et al._ , Generating solitons by phase engineering of a bose-einstein condensate, Science 287, 97 (2000).
  * Desyatnikov  _et al._ [2005] A. S. Desyatnikov, L. Torner, and Y. S. Kivshar, Optical vortices and vortex solitons, Prog. Opt. 47, 291 (2005).
  * Soto-Crespo  _et al._ [1991] J. Soto-Crespo, D. Heatley, E. M. Wright, and N. Akhmediev, Stability of the higher-bound states in a saturable self-focusing medium, Physical Review A 44, 636 (1991).
  * Malomed  _et al._ [2005] B. A. Malomed, D. Mihalache, F. Wise, and L. Torner, Spatiotemporal optical solitons, Journal of Optics B: Quantum and Semiclassical Optics 7, R53 (2005).
  * Yang and Pelinovsky [2003] J. Yang and D. E. Pelinovsky, Stable vortex and dipole vector solitons in a saturable nonlinear medium, Physical Review E 67, 016608 (2003).
  * Izdebskaya  _et al._ [2015] Y. Izdebskaya, G. Assanto, and W. Krolikowski, Observation of stable-vector vortex solitons, Optics Letters 40, 4182 (2015).
  * Izdebskaya  _et al._ [2012] Y. V. Izdebskaya, J. Rebling, A. S. Desyatnikov, and Y. S. Kivshar, Observation of vector solitons with hidden vorticity, Optics Letters 37, 767 (2012).
  * Porras  _et al._ [2016] M. A. Porras, M. Carvalho, H. Leblond, and B. A. Malomed, Stabilization of vortex beams in kerr media by nonlinear absorption, Physical Review A 94, 053810 (2016).
  * Quiroga-Teixeiro and Michinel [1997] M. Quiroga-Teixeiro and H. Michinel, Stable azimuthal stationary state in quintic nonlinear optical media, Journal of the Optical Society of America B 14, 2004 (1997).
  * Mihalache  _et al._ [2002] D. Mihalache, D. Mazilu, L.-C. Crasovan, I. Towers, A. Buryak, B. A. Malomed, L. Torner, J. Torres, and F. Lederer, Stable spinning optical solitons in three dimensions, Physical Review Letters 88, 073902 (2002).
  * Briedis  _et al._ [2005] D. Briedis, D. Petersen, D. Edmundson, W. Krolikowski, and O. Bang, Ring vortex solitons in nonlocal nonlinear media, Optics Express 13, 435 (2005).
  * Ramaniuk  _et al._ [2019] A. Ramaniuk, P. S. Jung, D. N. Christodoulides, W. Krolikowski, and M. Trippenbach, Absorption-mediated stabilization of nonlinear propagation of vortex beams in nematic liquid crystals, Optics Communications 451, 338 (2019).
  * Dabby and Whinnery [1968] F. Dabby and J. Whinnery, Thermal self-focusing of laser beams in lead glasses, Applied Physics Letters 13, 284 (1968).
  * Kwaśny  _et al._ [2020] M. Kwaśny, M. A. Karpierz, G. Assanto, and U. A. Laudyn, Optothermal vortex-solitons in liquid crystals, Optics Letters 45, 2451 (2020).
  * Suter and Blasberg [1993] D. Suter and T. Blasberg, Stabilization of transverse solitary waves by a nonlocal response of the nonlinear medium, Physical Review A 48, 4583 (1993).
  * Dalfovo  _et al._ [1999] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Theory of bose-einstein condensation in trapped gases, Reviews of Modern Physics 71, 463 (1999).
  * Peccianti and Assanto [2012] M. Peccianti and G. Assanto, Nematicons, Physics Reports 516, 147 (2012).
  * Peccianti  _et al._ [2000] M. Peccianti, A. De Rossi, G. Assanto, A. De Luca, C. Umeton, and I. Khoo, Electrically assisted self-confinement and waveguiding in planar nematic liquid crystal cells, Applied Physics Letters 77, 7 (2000).
  * Izdebskaya  _et al._ [2018] Y. V. Izdebskaya, V. G. Shvedov, P. S. Jung, and W. Krolikowski, Stable vortex soliton in nonlocal media with orientational nonlinearity, Optics Letters 43, 66 (2018).
  * Izdebskaya  _et al._ [2016] Y. Izdebskaya, W. Krolikowski, N. F. Smyth, and G. Assanto, Vortex stabilization by means of spatial solitons in nonlocal media, Journal of Optics 18, 054006 (2016).
  * Izdebskaya  _et al._ [2017] Y. Izdebskaya, V. Shvedov, G. Assanto, and W. Krolikowski, Magnetic routing of light-induced waveguides, Nature Communications 8, 1 (2017).
  * Laudyn  _et al._ [2020] U. A. Laudyn, M. Kwaśny, M. A. Karpierz, and G. Assanto, Vortex nematicons in planar cells, Optics Express 28, 8282 (2020).
  * Fleck and Feit [1983] J. A. Fleck and M. D. Feit, J. Opt. Soc. Am. , 920 (1983).
  * Alberucci  _et al._ [2010] A. Alberucci, A. Piccardi, M. Peccianti, M. Kaczmarek, and G. Assanto, Propagation of spatial optical solitons in a dielectric with adjustable nonlinearity, Physical Review A 82, 023806 (2010).
  * Piccardi  _et al._ [2010] A. Piccardi, A. Alberucci, and G. Assanto, Soliton self-deflection via power-dependent walk-off, Applied Physics Letters 96, 061105 (2010).
  * Assanto and Peccianti [2003] G. Assanto and M. Peccianti, Spatial solitons in nematic liquid crystals, IEEE Journal of Quantum Electronics 39, 13 (2003).
  * Khoo  _et al._ [1993] I. C. Khoo, H. Li, and Y. Liang, Optically induced extraordinarily large negative orientational nonlinearity in dye-doped liquid crystal, IEEE Journal of Quantum Electronics 29, 1444 (1993).
  * Molina-Terriza  _et al._ [2001] G. Molina-Terriza, E. M. Wright, and L. Torner, Propagation and control of noncanonical optical vortices, Optics Letters 26, 163 (2001).
  * Ye  _et al._ [2010] F. Ye, L. Dong, B. A. Malomed, D. Mihalache, and B. Hu, Elliptic vortices in optical waveguides and self-attractive bose-einstein condensates, Journal of the Optical Society of America B 27, 757 (2010).
  * Ge  _et al._ [2014] L. Ge, M. Shen, T. Zang, and L. Dai, Stability of elliptic vortex solitons in anisotropic nonlocal media, Chinese Optics Letters 12, 121901 (2014).
  * Gao  _et al._ [2020] M. Gao, Y. Chen, and M. Shen, Elliptic fundamental, dipole and vortex solitons in nonlocal nonlinear media with linear anisotropic diffraction, Journal of Optics 22, 025502 (2020).


# 2nd-order Updates with 1st-order Complexity

Michael F. Zimmer
<EMAIL_ADDRESS>

###### Abstract

It has long been a goal to efficiently compute and use second order
information on a function ($f$) to assist in numerical approximations. Here it
is shown how, using only basic physics and a numerical approximation, such
information can be accurately obtained at a cost of ${\cal O}(N)$ complexity,
where $N$ is the dimensionality of the parameter space of $f$. In this paper,
an algorithm (VA-Flow) is developed to exploit this second order information,
and pseudocode is presented. It is applied to two classes of problems, that of
inverse kinematics (IK) and gradient descent (GD). In the IK application, the
algorithm is fast and robust, and is shown to lead to smooth behavior even
near singularities. For GD the algorithm also works very well, provided the
cost function is locally well-described by a polynomial.

##  1 Introduction

The ability to compute second order information of a function adds obvious
benefits over only using first order information [Bishop, 2006]. Two obvious
applications are in optimization (gradient descent (GD)) and in robotics
(inverse kinematics (IK)). To date, the main challenge has been the complexity
costs (i.e., computations and storage), which are ${\cal O}(N^{2})$ when there
are $N$ parameters; in contrast, first-order computations have ${\cal O}(N)$
complexity. This paper will introduce an algorithm which extracts the benefits
of a second order algorithm, but only incurs ${\cal O}(N)$ complexity.

After discussing relevant research in the section Related Work, the new ideas
and algorithms are introduced in the sections Vector Flows, Usage, and The VA-
Flow Algorithm. The following three sections contain a simple demonstration on
an ellipse, numerical experiments on IK and GD, as well as a concluding
discussion. Finally, the reader will find additional material in the
Supplement and on the author’s Github site [Zimmer, 2021a].

##  2 Related Work

Work that relates to second order corrections (e.g., Hessians) and to
applications that make use of such corrections (e.g., GD and IK) are of
primary interest here. In the Summary sub-section, a general characterization
of the applications is given, which will be used throughout this paper.

### Inverse Kinematics (IK)

This application concerns itself with adjusting the pose of an articulated
figure [e.g., see Sciavicco and Siciliano, 2001, Buss, 2009, Aristidou and
Lasenby, 2009]. For concreteness, if a human wishes to move his hand from
point A to point B, while keeping his torso fixed, the problem becomes one of
determining the changes in the orientation of the shoulder, elbow, and wrist.
More abstractly, the arm is described as links (the bones of the arm)
connected by joints (e.g., the elbow), and the standard IK problem is to
determine how to adjust the joints in order to move the end of the last link
(the end-effector) from point A to B. The position of end-effector may be
formulated as

$\displaystyle{\bm{r}}=f({\bm{\theta}})\,\,\,,$

where ${\bm{\theta}}=(\theta_{1},\dots,\theta_{N})$ is a parameterization of
the $N$ joints, and ${\bm{s}}$ is the position of the end-effector.

Following Aristidou and Lasenby [2009], there are four main classes of
solutions for this problem: Jacobian inverse, Newton, Heuristics, and Data-
driven. (Heuristic and Data-Driven methods are not of direct interest for this
paper, and will not be discussed.) The Jacobian inverse method centers on
making small changes $d{\bm{r}}$, $d{\bm{\theta}}$, and are related by

$\displaystyle d{\bm{r}}={\bm{J}}d{\bm{\theta}}\,\,\,,$ (1)

where ${\bm{J}}=df/d{\bm{\theta}}$ is the Jacobian. (For example, if the goal
is to move the end-effector from ${\bm{r}}_{A}$ to ${\bm{r}}_{B}$, one could
define ${\bm{s}}={\bm{r}}_{B}-{\bm{r}}_{A}$ and
$\hat{{\bm{s}}}={\bm{s}}/\|{\bm{s}}\|$ so that
$d{\bm{r}}=\alpha\hat{{\bm{s}}}$. Here $\alpha$ is the magnitude of the move,
which may be taken as the minimum of $\|{\bm{s}}\|$ and a user-defined
maximum.) Eq. 1 is normally an underdetermined system of equations; it
additionally offers the challenge of singular points in ${\bm{\theta}}$-space,
which are associated with a change in rank of the matrix $J$. Some of the
approaches that have appeared in this class of solutions are: Jacobian
pseudoinverse [Whitney, 1969], damped least-squares [Wampler, 1986, Nakamura
and Hanafusa, 1986], singular value filtering [Colomè and Torras, 2012],
pseudoinverse damped least-squares [Maciejewski, 1990, Press et al., 2007],
selectively damped least-squares [Buss and Kim, 2005], and others. Of these
named methods, the fastest (and least robust) is the Jacobian pseudoinverse,
and will be the baseline for comparisons herein. It obtains a solution in the
form [Zimmer, 2020]

$\displaystyle d{\bm{\theta}}$
$\displaystyle={\bm{G}}d{\bm{r}}+d{\bm{\theta}}_{null}$
$\displaystyle{\bm{P}}$ $\displaystyle=1-{\bm{G}}{\bm{J}}\,\,\,,$

where $d{\bm{\theta}}_{null}$ is in the null space of ${\bm{J}}$, ${\bm{G}}$
is the generalized inverse [Ben-Israel and Greville, 2003] of ${\bm{J}}$, and
${\bm{P}}$ is a projection operator for the null space of ${\bm{J}}$. When
${\bm{G}}$ is of type {1234}, it is known as a Moore-Penrose inverse (often
written as ${\bm{J}}^{\dagger}$). The operator ${\bm{P}}$ can be used for
secondary goals [Liegeois, 1977] of the articulated arm, such as satisfying
joint constraints. Also, Newton methods give an accurate 2nd-order method, but
are computationally very expensive to implement. In summary, the challenge
with the IK problem is that existing algorithms have been either fast or
robust, but not both.

IK algorithms find applications wherever articulated limbs are used, such as
computer graphics and robotic arms (in manufacturing, human-like robots, arms
in space travel, prosthetic limbs, etc.). Thus any opinion of the ethics of
this application is simply a reflection of the ethics of those existing
technologies, which in the author’s opinion has been positive thus far.

### Gradient Descent (GD)

Gradient descent [Cauchy, 1847, Lemaréchal, 2010, Press et al., 2007] is an
iterative optimization algorithm to be used for the minimization of a
differentiable function. With respect to the cost function (CF) $f$, which
depends on the parameter vector ${\bm{\theta}}=(\theta_{1},\theta_{2},...)$,
an initial value ${\bm{\theta}}_{old}$ is updated to ${\bm{\theta}}_{new}$ as

$\displaystyle{\bm{\theta}}_{new}$
$\displaystyle={\bm{\theta}}_{old}-\alpha{\bm{\nabla}}f\,\,\,,$

where ${\bm{\nabla}}f$ is the gradient of $f$ with respect to ${\bm{\theta}}$,
and $\alpha$ is a positive number known as the learning rate.

The topic of iterative approaches to optimization has been well investigated
(e.g., see [Nesterov, 2018, Sra et al., 2011, Nocedal and Wright, 2006]), and
most methods are classified as either line search or trust region approaches.
Other useful techniques include momentum and stochastic gradient descent
[e.g., see Goodfellow et al., 2016].

There is a closely related class of algorithms that involve momentum [Polyak,
1964, Nesterov, 1983] being applied to the gradient and/or gradient-squared,
effectively creating a different learning rate for each component of
${\bm{\theta}}$. In this category are the algorithms AdaGrad [Duchi et al.,
2011], AdaDelta [Zeiler, 2012], RMS Prop [Tieleman and Hinton, 2012], Adam &
AdaMax [Kingma and Ba, 2015], Nadam [Dozat, 2016], AdamNC [Reddi et al.,
2018], and others [see Ruder, 2017]. At this point, Adam seems to be the
preferred choice, and it will be used for comparisons herein. These algorithms
were originally meant to be used with stochastic gradient descent, but may be
used in batch GD as well.

Additionally, the recent Neograd algorithm [Zimmer, 2021b] formulaically
adjusts the learning rate to keep a metric near its target value. It uses the
philosophy that the deviation from the first-order update should be relatively
small compared to the first-order update itself. In batch optimizations
compared to Adam, it was shown to reach lower CF values by factors such as
$10^{8}$ to $10^{12}$.

Finally, GD finds use in myriad applications, one of them being supervised
learning for image detection. That application has ethically positive uses
(e.g., detecting tumors from X-ray images), as well as dystopian ones (e.g.,
pervasive facial recognition technology in society). Thus one will have to
survey all such applications to decide whether the net ethical impact is good
or bad.

### Hessians

In a Taylor series expansion of a function $f$, the Hessian ($H$) appears as
the matrix in the second order term. If $f$ depends on $N$ parameters, then
$H$ has $N^{2}$ terms, and thus a direct computation of $H$ becomes
impracticable for large N. A quasi-Newton method [cf. Nocedal and Wright,
2006] can be used as an approximation toward computing it, and then
determining its impact. Other approximations for computing $H$ include
Broyden’s method, Powell’s method, Siciliano’s method, as well as the BFGS
method and its variants [see Fletcher, 1987, Siciliano, 1990, Chin et al.,
1997]. In general, these methods are known for being accurate and more stable,
but computationally costly.

In two independent efforts, Møller [1993a, d] and Pearlmutter [1994]
determined a way to use $H$ without explicitly calculating and storing its
$N^{2}$ entries. Noting that the Hessian normally appears as $H{\bm{w}}$,
where ${\bm{w}}$ is a suitable vector, they each observed that it was
accessible by a numerical approximation via
$H{\bm{w}}\approx(g({\bm{\theta}}+\epsilon{\bm{w}})-g({\bm{\theta}}))/\epsilon$,
where $g({\bm{\theta}})$ is the gradient of $f$ and $\epsilon$ is a small
number. While they primarily employed this in the context of neural nets, they
both also discussed how it can be used in the line search algorithm to set the
step size ($\alpha$), such as with scaled conjugate gradient [Møller, 1993b,
c]. They both also noted its use in estimating eigenvalues and eigenvectors of
$H$.

### Summary

The applications of IK and GD may be treated side by side by identifying the
vector field ${\bm{v}}({\bm{\theta}})$

${\bm{v}}({\bm{\theta}})=\left\\{\begin{array}[]{@{\quad}l@{\quad}l@{}}-{\bm{\nabla}}f|_{\bm{\theta}}&\text{
(gradient descent) }\\\ G({\bm{\theta}})\hat{{\bm{s}}}&\text{ (inverse
kinematics) },\\\ \end{array}\right.$ (2)

where all quantities are as defined above. In the algorithms of interest the
parameter ${\bm{\theta}}$ is iteratively updated by the amount

$\displaystyle d{\bm{\theta}}=\alpha{\bm{v}}\,\,\,,$

where $\alpha$ is a positive number.

##  3 Vector Flows

Computing updates of ${\bm{\theta}}$ with a single vector ${\bm{v}}$ is
convenient and quick, but a more accurate prediction would involve the rate of
change of ${\bm{v}}$. Here a numerical estimate is made of the rate of change
of ${\bm{v}}$ at the two nearby values ${\bm{\theta}}_{0}$ and
${\bm{\theta}}_{1}$. Instead of making a single update using $\alpha$, a
smaller update using $\epsilon$ is made. In this paper this value is set as

$\displaystyle\epsilon$ $\displaystyle=\alpha/M\ll\alpha\,\,\,,$

where $M$ is a suitably large number (cf. Sec. 5). This allows a calculation
of what is called the acceleration

$\displaystyle{\bm{a}}$
$\displaystyle\equiv({\bm{v}}_{1}-{\bm{v}}_{0})/\epsilon\,\,\,,$ (3)

where ${\bm{v}}_{k}={\bm{v}}({\bm{\theta}}_{k})$ and
${\bm{\theta}}_{1}={\bm{\theta}}_{0}+\epsilon{\bm{v}}_{0}$. Next, the
compounded effect of ${\bm{a}}$ on the scale of $\epsilon$ is computed, to
determine its influence on the scale of $\alpha$. Assuming that ${\bm{a}}$
remains constant for these subsequent updates, the net effect can be computed
as (with $n\geq 1$)

$\displaystyle{\bm{v}}_{n}$ $\displaystyle={\bm{v}}_{n-1}+\epsilon{\bm{a}}$
$\displaystyle={\bm{v}}_{0}+n\epsilon{\bm{a}}$ (4a)
$\displaystyle{\bm{\theta}}_{n}$
$\displaystyle={\bm{\theta}}_{n-1}+\epsilon{\bm{v}}_{n-1}$
$\displaystyle={\bm{\theta}}_{0}+\epsilon({\bm{v}}_{0}+\cdots+{\bm{v}}_{n-2}+{\bm{v}}_{n-1})$
$\displaystyle={\bm{\theta}}_{0}+n\epsilon{\bm{v}}_{0}+\frac{1}{2}n(n-1)\epsilon^{2}{\bm{a}}\,\,\,.$
(4b)

Assuming large $n$, and identifying $\alpha=n\epsilon$, it approximately
follows that

$\displaystyle{\bm{v}}_{\alpha}$ $\displaystyle={\bm{v}}_{0}+\alpha{\bm{a}}$
(5a) $\displaystyle{\bm{\theta}}_{\alpha}$
$\displaystyle={\bm{\theta}}_{0}+\alpha{\bm{v}}_{0}+\frac{\alpha^{2}}{2}{\bm{a}}\,\,\,.$
(5b)

These two update formulas should be immediately recognizable as two equations
from elementary physics: the former for the velocity ${\bm{v}}$ due to an
acceleration ${\bm{a}}$ over a time $\alpha$ given an initial velocity
${\bm{v}}_{0}$, and the latter for the change in position ${\bm{\theta}}$ due
to the same.

###  3.1 Further Aspects

This section collects several incidental results related to the previous
section. Eq. 3 can be written as

$\displaystyle{\bm{a}}$
$\displaystyle=[{\bm{v}}({\bm{\theta}}_{0}+d\hat{{\bm{\theta}}})-{\bm{v}}({\bm{\theta}}_{0})]/\epsilon$
$\displaystyle\approx(d\hat{{\bm{\theta}}}\cdot{\bm{\nabla}}){\bm{v}}/\epsilon$
$\displaystyle\approx({\bm{v}}_{0}\cdot{\bm{\nabla}}){\bm{v}}\,\,\,,$ (6a)

where $d\hat{{\bm{\theta}}}=\epsilon{\bm{v}}_{0}$. Evaluations of the
derivative are understood to be made at ${\bm{\theta}}_{0}$. Note that in the
context of GD where ${\bm{v}}=-{\bm{\nabla}}f$, ${\bm{a}}$ becomes
${\bm{\nabla}}f({\bm{\nabla}}{\bm{\nabla}}f)$, which was studied by Møller
[1993a, d] and Pearlmutter [1994]. Using this expression for ${\bm{a}}$, and
now setting $d{\bm{\theta}}=\alpha{\bm{v}}_{0}+\alpha^{2}{\bm{a}}/2$, a
second-order approximation for the CF $f({\bm{\theta}})$ can be written as

$\displaystyle f$
$\displaystyle=f_{0}+({\bm{\nabla}}f)^{T}d{\bm{\theta}}+\frac{1}{2}d{\bm{\theta}}^{T}({\bm{\nabla}}{\bm{\nabla}}f)d{\bm{\theta}}+{\cal
O}(\|d{\bm{\theta}}\|^{3})$ $\displaystyle=f_{0}-\alpha
v_{0}^{2}-\alpha^{2}{\bm{v}}_{0}\cdot{\bm{a}}+{\cal O}(\alpha^{3})\,\,\,,$
(7a)

where $v_{0}^{2}={\bm{v}}_{0}\cdot{\bm{v}}_{0}$, and
$f_{0}=f({\bm{\theta}}_{0})$. Note the improvement in estimating $f$ due to
the term $-\alpha^{2}{\bm{v}}_{0}\cdot{\bm{a}}$.

Finally, a differential geometric [e.g., see Lee, 2013, Frankel, 2011, Tu,
2011] interpretation of ${\bm{a}}$ may be made, beginning from an expansion of
the vector field ${\bm{v}}$ at ${\bm{\theta}}_{0}$:

$\displaystyle{\bm{v}}={\bm{v}}_{0}+(d{\bm{\theta}}\cdot{\bm{\nabla}}){\bm{v}}+\cdots\,\,\,.$
(8)

Here $d{\bm{\theta}}={\bm{\theta}}-{\bm{\theta}}_{0}$ is the only
${\bm{\theta}}$ dependence on the right hand side, since the gradient is
evaluated at ${\bm{\theta}}_{0}$. In the neighborhood of ${\bm{\theta}}_{0}$,
define the vector fields $X$ and $Y$ based on the contributions from the first
and second order terms of Eq. 8

$\displaystyle X$ $\displaystyle={\bm{v}}_{0}\cdot{\bm{\nabla}}$
$\displaystyle Y$
$\displaystyle={\bm{v}}_{0}\cdot{\bm{\nabla}}+[(d{\bm{\theta}}\cdot{\bm{\nabla}}){\bm{v}}]\cdot{\bm{\nabla}}\,\,\,.$

The vector field $X$ represents basic GD, while $Y$ also includes the
contribution from a changing ${\bm{v}}$. Thus we are interested in the rate of
change of $Y$ due to the flow created by $X$. This naturally motivates the
consideration of the Lie derivative [introduced by Ślebodziński, 1931] of $Y$
with respect to $X$, which when evaluated at ${\bm{\theta}}_{0}$ is

$\displaystyle{\cal L}_{X}Y$ $\displaystyle=[X,Y]$
$\displaystyle=[({\bm{v}}_{0}\cdot{\bm{\nabla}}){\bm{v}}]\cdot{\bm{\nabla}}$
$\displaystyle={\bm{a}}\cdot{\bm{\nabla}}\,\,\,,$

where the last line made use of Eq. 6a. Thus, the vector ${\bm{a}}$ in the
small-$\epsilon$ limit of Eq. 3 may also be understood to be equivalent to the
Lie derivative ${\cal L}_{X}Y$. This entire procedure could be repeated for
higher order corrections in Eq. 8. Also, note that the original expression for
${\bm{a}}$ (i.e., Eq. 3) was valid because the underlying space is Euclidean.
Thus even though ${\bm{v}}_{0}$ and ${\bm{v}}$ live in different tangent
spaces and in general can’t be directly compared, they can be in this case
(cf. p.228 in Lee [2013]).

##  4 Usage

This section introduces several heuristics for using the second order update
formulas in Eqs. 4. They’re based on the principal that the correction to the
first-order update should be smaller than the first-order update itself [cf.
Zimmer, 2021b]. Note that this is different from line search techniques, which
attempt to lower the CF as much as possible in a single update. (Although,
Wolfe conditions [Armijo, 1966, Wolfe, 1969, 1971] are a refinement on this.)

###  4.1 Approach A

The metric $\rho_{A}$ is defined as the ratio of the magnitudes of the second
term to the first from Eq. 4b

$\displaystyle\rho_{A}$
$\displaystyle=\frac{\frac{1}{2}n(n-1)\epsilon^{2}\|{\bm{a}}\|}{n\epsilon\|{\bm{v}}\|}=\frac{(n-1)\epsilon}{2}\frac{\|{\bm{a}}\|}{\|{\bm{v}}\|}\,\,\,.$

This equation provides a way to solve for the value of $n=n^{*}$ that leads to
a specified value for $\rho_{A}$, which in this case is the target value
$\rho_{targ}$

$\displaystyle n^{*}$
$\displaystyle=1+\frac{2\rho_{targ}}{\epsilon}\frac{\|{\bm{v}}\|}{\|{\bm{a}}\|}\,\,\,.$
(9)

This value of $n^{*}$ reveals how many 111The value of $n^{*}$ should be cast
as an integer, but that restriction isn’t necessary for its use here. updates
should be taken to maintain the metric at the value $\rho_{targ}$. This value
should be substituted into Eq. 4b to obtain the new parameter value. Thus,
this approach provides for a single update $\theta_{n^{*}}$, starting from
$\theta_{0}$. Note that for this approach to make sense, $n^{*}\gg 1$; if not
then then it’s essentially just GD. This is the preferred approach in this
paper.

###  4.2 Approach B

The approach here is similar to that for $\rho_{A}$, except now the
contributions from the right-hand side of Eq. 4b are measured along the
direction of $\hat{{\bm{v}}}={\bm{v}}/\|{\bm{v}}\|$. The metric is

$\displaystyle\rho_{B}$
$\displaystyle=\left|\frac{\hat{{\bm{v}}}\cdot[\frac{n(n-1)\epsilon^{2}}{2}{\bm{a}}]}{\hat{{\bm{v}}}\cdot[n\epsilon{\bm{v}}]}\right|=\frac{(n-1)\epsilon}{2}\frac{|{\bm{v}}\cdot{\bm{a}}|}{v^{2}}\,\,\,.$

Note this is the same (up to a factor of 2) as taking the ratio of the second
to the first correction in Eq. 7a. As before, this equation provides for a way
to solve for the value of $n=n^{*}$ that leads to the value
$\rho_{B}=\rho_{targ}$

$\displaystyle n^{*}$
$\displaystyle=1+\frac{2\rho_{targ}}{\epsilon}\frac{v^{2}}{|{\bm{v}}\cdot{\bm{a}}|}\,\,\,.$

This value should be used in the same manner as with Approach A. Finally, in
the context of GD where ${\bm{a}}=({\bm{v}}_{0}\cdot{\bm{\nabla}}){\bm{v}}$
and ${\bm{v}}=-{\bm{\nabla}}f$, it easily follows that in the small
$d{\bm{\theta}}$ limit, it becomes identical to the metric from Zimmer [2021b]

$\displaystyle\rho_{B}$
$\displaystyle=\frac{1}{2}\left|\frac{d{\bm{\theta}}^{T}({\bm{\nabla}}{\bm{\nabla}}f)d{\bm{\theta}}}{{\bm{\nabla}}f^{T}d{\bm{\theta}}}\right|\,\,\,.$

That metric was also based on keeping the deviation from the first-order
correction small.

###  4.3 Approach C

In a similar manner, comparing the second to the first terms of the right hand
side of Eq. 4a leads to the metric

$\displaystyle\rho_{C}=\frac{\|n\epsilon{\bm{a}}\|}{\|{\bm{v}}\|}\,\,\,.$

In this case the number of updates of size $\epsilon$ needed to reach
$\rho_{C}=\rho_{targ}$ is

$\displaystyle
n^{*}=\frac{\rho_{targ}}{\epsilon}\frac{\|{\bm{v}}\|}{\|{\bm{a}}\|}\,\,\,.$

For large $n^{*}$, this becomes equivalent to Eq. 9, except it is missing the
factor of 2.

##  5 The VA-Flow Algorithm

Pseudocode for the ideas discussed thus far is given below in Algo. 1 in what
is called the VA-Flow algorithm. Note that lines 4-7 correspond to the
computations of ${\bm{v}}$ and ${\bm{a}}$ in Sec. 3, with the function get_v
implementing the appropriate choice from Eq. 2. Line 9, which includes the
function get_dtheta, implements Eq. 4b; get_nstar implements Eq. 9. (Approach
B or C could also be used). These are the most essential elements of the
algorithm. The parameter values generally used by the author were
$\rho_{targ}=0.1$ and $M=100$; the value of $num$ depends on the application.
(This value of $M$ was chosen to provide room for $n^{*}$ to adjust while
still keeping $n^{*}\gg 1$.)

Algorithm 1 : VA-Flow, version 0

1: Input: (${\bm{\theta}}$, $\alpha$, $\rho_{targ}$, $M$, $num$ )

2: for i = 1 to num do

3: $\epsilon=\alpha/M$

4: ${\bm{v}}=\text{get\\_v}({\bm{\theta}})$

5: ${\bm{\theta}}_{1}={\bm{\theta}}+\epsilon{\bm{v}}$

6: ${\bm{v}}_{1}=\text{get\\_v}({\bm{\theta}}_{1})$

7: ${\bm{a}}=({\bm{v}}_{1}-{\bm{v}})/\epsilon$

8: $n^{*}=\text{get\\_nstar}({\bm{v}},{\bm{a}},\rho_{targ},\epsilon)$

9:
${\bm{\theta}}={\bm{\theta}}+\text{get\\_dtheta}({\bm{\theta}},{\bm{v}},{\bm{a}},n^{*},\epsilon)$

10: $\alpha=n^{*}\epsilon$

11: end for

12: Return: ${\bm{\theta}}$

The reader should note the interplay between the statements on lines 3 and 10,
which allow $\alpha$ and $\epsilon$ to adjust. On line 10, $\alpha$ is reset
using the most recent $n^{*}$ value, and on line line 3, $\epsilon$ is again
rescaled to be much smaller than $\alpha$. It is important that $n^{*}\gg 1$,
to conform to the approach given in Sec. 3. If $n^{*}\approx 1$, then it is
effectively the same as GD, and no second-order benefits can be expected.
Also, the reader may instead use get_dtheta to simply make a first-order
update of ${\bm{\theta}}$; in this case, this algorithm is only used as a
means to set $\alpha$ in an otherwise first-order GD. Finally, notice that VA-
Flow is an ${\cal O}(N)$ algorithm, since the additional operations beyond GD
in lines 6-9 are all ${\cal O}(N)$. (It is assumed that get_v is ${\cal
O}(N)$.)

Note that the benefit of using ${\bm{a}}$ is perhaps greatest when
$\|{\bm{a}}\|/\|{\bm{v}}\|$ is relatively large. It is indicative of a rapidly
changing ${\bm{v}}$, and thus provides important information regarding the
next update of ${\bm{\theta}}$. In contrast, when that ratio is relatively
small ${\bm{a}}$ becomes less important, and it can lead to large values of
$n^{*}$. In this case it may be necessary to use other information to limit
the size of $\alpha$. (Alternatively, one could simply default to the
underlying first-order algorithm.) It will be shown that for IK and GD on a
polynomial CF it is not an issue, but it can be for GD on a CF with
exponential terms.

Finally, there are modifications that can be made to improve performance and
robustness, and they are discussed in the Supplement. There it is shown how to
keep $n^{*}\gg 1$ in extreme cases, by repeating certain lines of code with a
reduced $\alpha$. In addition, VA-Flow can be modified with the heuristic of
including momentum on the updates of ${\bm{\theta}}$. (Although, it will be
shown that momentum isn’t needed as much as it is with basic GD, since the
updates are more accurate.) Also, there are a number of coding improvements
that should be utilized in a real implementation, but which weren’t shown in
Algo. 1 for simplicity. Examples include checks to avoid dividing by zero; see
the implementation by Zimmer [2021a].

##  6 Simple Demonstration

The purpose of this section is to illustrate the difference between Basic GD
and the new VA-Flow algorithm on a simple example, that of the optimization of
the ellipse cost function $f$

$\displaystyle f$
$\displaystyle=\frac{\theta_{1}^{2}}{c_{1}^{2}}+\frac{\theta_{2}^{2}}{c_{2}^{2}}\,\,\,.$

This comparison will be illustrated using a single update of ${\bm{\theta}}$.
Following Eq. 5b, the update for ${\bm{\theta}}$ for Basic GD (i.e.,
${\bm{\theta}}_{GD}$), and for VA-Flow (i.e., ${\bm{\theta}}_{Flow}$) appear
as the following functions of $\alpha$

$\displaystyle{\bm{\theta}}_{GD}$
$\displaystyle={\bm{\theta}}_{0}+\alpha{\bm{v}}$
$\displaystyle{\bm{\theta}}_{Flow}$
$\displaystyle={\bm{\theta}}_{0}+\alpha{\bm{v}}+\frac{\alpha^{2}}{2}{\bm{a}}\,\,\,,$

where ${\bm{\theta}}_{0}$ is the initial condition for the parameter. Note
that ${\bm{\theta}}_{Flow}$ has an additional quadratic $\alpha$-dependence
that ${\bm{\theta}}_{GD}$ does not have; this allows it to more closely follow
the exact path from infinitesimal updates. Using the parameter choices
$c_{1}=6$, $c_{2}=2$ and ${\bm{\theta}}_{0}^{T}=(4,1.5)$, it follows

$\displaystyle{\bm{v}}^{T}$
$\displaystyle=-2(\theta_{1}/c_{1}^{2},\theta_{2}/c_{2}^{2})=(-0.22...,-0.75)$
$\displaystyle{\bm{a}}^{T}$
$\displaystyle=4(\theta_{1}/c_{1}^{4},\theta_{2}/c_{2}^{4})=(-0.012...,0.375)\,\,\,,$

from which both ${\bm{\theta}}_{GD}$ and ${\bm{\theta}}_{Flow}$ are plotted as
functions of $\alpha$ in Fig. 1.

Figure 1:  A close-up view of the updates for basic GD (orange) and VA-Flow
(green) for a range of $\alpha$ values when seeking the minimum of an ellipse
function. The light gray lines are contours of the CF, and the light red line
is the exact path toward the minimum. The black dot is the starting point for
both, and the red dot is the recommended update from VA-Flow.

The green line represents ${\bm{\theta}}_{Flow}$, while the orange line is
${\bm{\theta}}_{GD}$. Also, ${\bm{\theta}}_{0}$ is depicted with a black dot.
Using Approach A and requiring $\rho_{targ}=0.2$, it follows
$\alpha_{A}\approx 0.834$. Using $\alpha_{A}$, the recommended
${\bm{\theta}}_{Flow}$ for a single update can be computed (following Approach
A) and appears as a red dot in the figure. The dotted red line is the exact
path from infinitesimal GD updates; it is described by
$\theta_{2}\sim\theta_{1}^{c}$ with $c=(c_{1}/c_{2})^{2}$, as shown in the
Supplement. This figure makes clear the advantages VA-Flow has to offer: a
means to determine a best value of $\alpha$, and a more accurate update of
${\bm{\theta}}$.

##  7 Application: Inverse Kinematics

This application is based on a 3-link arm, where the links (each of length
$1.0$) are constrained to move in a 2D plane. The parameter vector
${\bm{\theta}}=(\theta_{1},\theta_{2},\theta_{3})$ describes the angles for
the three joints, indicated by blue dots in the left plot of Fig. 2; there are
no constraints on these angles.

Figure 2:  In the left plot, point B (i.e., the red dot) is the target, while
point A (i.e., the green dot) is the starting location for the end effector.
The three blue dots are the joints of the arm, and are parametrized by
${\bm{\theta}}$. The light blue shading indicates the range of motion of the
arm. The right plot graphs the distance between ${\bm{s}}_{A}$ and
${\bm{s}}_{B}$ for each algorithm, revealing a smooth and efficient approach
for VA-Flow and an erratic approach for the Jacobian Pseudoinverse method.

A simple test suffices to illustrate the difference between an arm controlled
by VA-Flow (using Approach A) and the method of Jacobian Pseudoinverse
discussed earlier. The test is merely to move the end-effector from point A to
B, starting from the initial pose shown in the left plot of Fig. 2; the two
points are ${\bm{s}}_{A}=(1.71,1.71)$ and ${\bm{s}}_{B}=(2.132,2.132)$. Point
B is just beyond the range of motion of the arm. The results for each
algorithm are shown in the right plot of the figure, which measures the
distance $\|{\bm{s}}_{B}-{\bm{s}}_{A}\|$ versus iterations. It reveals that
VA-Flow leads to a very smooth approach to the target, with no oscillations or
erratic behavior. In contrast, the graph for the basic Jacobian Pseudoinverse
approach reveals very erratic updates: as the end effector becomes close to
the target, it suddenly moves away. In general, it was observed that the
Jacobian Pseudoinverse algorithm performed well when the target was within the
range of motion, but poorly outside this range. Also, it performs poorly when
the arm is near a singular configuration (i.e., when two or more rows of the
Jacobian are nearly equal). In contrast, the VA-Flow algorithm performed well
in all circumstances. The instability was not manifested with VA-Flow because
it is a second-order algorithm, and can better navigate such singular
configurations. Also, as pointed out earlier, second-order methods based on
the quasi-Newton method do not exhibit erratic behavior either. Although those
methods are computationally much more costly in comparison to VA-Flow.

Finally, there were some minor adjustments made to the VA-Flow algorithm for
the IK application; these are discussed in the Supplement. Primary among them
is that $\alpha$ should be limited by (1) the maximum desired update size of
the end-effector, and (2) the remaining distance to the target.

##  8 Application: Gradient Descent

The VA-Flow algorithm can largely be applied to a GD example in the same way
as for IK. Here the example is to determine the minimum of Beale’s function,
defined by

$f=(1.5-\theta_{1}+\theta_{1}\theta_{2})^{2}+(2.25-\theta_{1}+\theta_{1}\theta_{2}^{2})^{2}+(2.625-\theta_{1}+\theta_{1}\theta_{2}^{3})^{2}\,\,\,.$

The known global minimum of this function is
$(\theta_{1},\theta_{2})=(3,0.5)$; the initial point used is $(4,3)$, from
which the minimum is reachable. In Fig. 3 a comparison is made between VA-Flow
and Adam on this CF.

Figure 3:  A comparison figure for the VA-Flow and Adam algorithms. In the
left plot, the log of the CF is plotted for each algorithm, revealing that VA-
Flow reaches a value $10^{-4}$ the size of the CF reached by Adam. The right
plot measures the distance from the current ${\bm{\theta}}$ to the known
global minimum, and likewise reveals superior performance by VA-Flow.

For VA-Flow, the learning rate was $\alpha=4.8\times 10^{-6}$, Approach A was
used, and $M=100$. Also, no momentum was used for VA-Flow. For Adam, $\alpha$
was set to $0.15$, and the other parameters had default settings. The most
obvious difference between these two algorithms is that, as shown in the left
plot, VA-Flow is able to reach much lower CF values, which is the desired
result. Also, the right plot displays the distance between the current
${\bm{\theta}}$ and the target, also revealing the superior performance of VA-
Flow.

Another example investigated by the author was using VA-Flow to minimize the
cross entropy CF, which was defined on a neural net (NN) model for classifying
handwritten digits [Pedregosa et al., 2011, Scikit-learn, 2020, Dua and Graff,
2017]. The NN had one hidden layer with 30 nodes; there were 64 input nodes
and 10 output nodes. The CF on the final layer was a cross entropy cost
function using softmax, and the hidden layer used a tanh activation function.
Altogether, the dimensionality of ${\bm{\theta}}$ was 2260. As it turned out,
the VA-Flow algorithm fared poorly in the initial iterations, since
$\|{\bm{a}}\|\ll\|{\bm{v}}\|$. This led to a large $n^{*}$ (and $\alpha$) and
subsequently poor behavior. In retrospect this is perhaps not surprising. Even
though the VA-Flow algorithm is using 2nd-order information, it is purely
local information. Such localized information is of limited value on the cross
entropy CF, which is composed of exponentials, and is characterized by very
flat and very steep regions. A better approach is one that contains
information over a neighborhood of values surrounding the current
${\bm{\theta}}$, such as the Neograd algorithm [Zimmer, 2021b].

##  9 Discussion

A method has been introduced to efficiently calculate the rate of change of a
vector field, and in so doing, compute with ${\cal O}(N)$ complexity a task
which has long been construed to be an ${\cal O}(N^{2})$ task. From this, the
algorithm VA-Flow was created, which was based on the philosophy that the
correction to a first-order update should be smaller than the first-order
update itself. The value of the parameter $\alpha$ determined by this
algorithm can be used within a first or second order version of the algorithm.
It is noted that a numerical estimate of the quantity ${\bm{a}}$ was
previously done by Møller [1993a, d] and Pearlmutter [1994]. However, they
made different applications of it, and in particular did not derive the second
order update formula for ${\bm{\theta}}$ (Eq. 5b), which was essential for the
VA-Flow algorithm.

When applied to inverse kinematics (IK), the algorithm proved to be extremely
successful. It was capable of moving an articulated limb in ${\cal O}(N)$
time, without any erratic movement whatsoever. Before VA-Flow, that was only
capable with algorithms that required ${\cal O}(N^{2})$ computations.

VA-Flow was also applied to a polynomial test function, and performed much
better than Adam, reaching a CF that was about $10^{4}$ times smaller than
that reached with Adam. However, when it was applied to a more realistic
application (a modestly sized NN), it didn’t fare as well. The author reasons
it is because the VA-Flow places too much importance on measurements made
approximately at the same point. While this works well on polynomials, it is
perhaps not well-suited for functions involving steep and flat regions caused
by exponentials. An algorithm such as Neograd [Zimmer, 2021b] seems to be a
better choice.

There are a number of ways to move forward with the results presented here.
First, one might consider an additional measurement of
${\bm{v}}({\bm{\theta}})$ in the neighborhood of ${\bm{\theta}}_{0}$, in order
to assess the rate of change of ${\bm{a}}$ itself. Another approach is to
assume that ${\bm{v}}$ changes in an exponential manner, and model it
appropriately (perhaps using two or three measurements of ${\bm{v}}$). This
could be helpful when working with CFs from NNs.

Finally, the mathematical results presented in this paper have myriad
applications with varying ethical impacts. Ethical considerations were made on
the IK and GD applications in Sec. 2.

## References

  * Aristidou and Lasenby [2009] A. Aristidou and J. Lasenby.  Inverse kinematics: a review of existing techniques and introduction of a new fast iterative solver.  Technical Report CUED/F-INFENG/TR-632, University of Cambridge, 09 2009\.
  * Armijo [1966] L. Armijo.  Minimization of functions having Lipschitz continuous first partial derivatives.  _Pacific J. Math._ , 16(1):1–3, 1966.
  * Ben-Israel and Greville [2003] A. Ben-Israel and T.N.E. Greville.  _Generalized Inverses, Theory and Applications_.  Springer, 2nd edition, 2003.
  * Bishop [2006] C.M. Bishop.  _Pattern Recognition and Machine Learning_.  Springer, 2006.
  * Buss [2009] S.R. Buss.  Introduction to inverse kinematics with jacobian transpose, pseudoinverse and damped least squares methods.  Technical report, University of California, 2009.
  * Buss and Kim [2005] S.R. Buss and J.-S. Kim.  Selectively damped least squares for inverse kinematics.  _Journal of Graphics Tools_ , 10(3):37–49, 2005\.
  * Cauchy [1847] A. Cauchy.  Méthode générale pour la résolution des systèmes d’équations simultanées.  _C. R. Acad. Sci. Paris_ , 25:536–538, 1847.
  * Chin et al. [1997] K.W. Chin, B.R. von Konsky, and A. Marriott.  Closed form and generalized inverse kinematics solutions for the analysis of human motion.  In _In Proceedings of the IEEE Conference of Engineering in Medicine and Biology Society_ , volume 5, page 1911–1914, 1997.
  * Colomè and Torras [2012] A. Colomè and C. Torras.  Redundant inverse kinematics: Experimental comparative review and two enhancements.  In _In Pro- ceedings of IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , page 5333–5340, October 2012.
  * Dozat [2016] T. Dozat.  Incorporating Nesterov momentum into Adam.  _ICLR Workshop_ , 1:2013–206, 2016.
  * Dua and Graff [2017] D. Dua and C. Graff.  UCI machine learning repository, 2017.  URL http://archive.ics.uci.edu/ml.
  * Duchi et al. [2011] J. Duchi, E. Hazan, and Y. Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  _Journal of Machine Learning Research (JMLR)_ , 12:2121–2159, 2011.
  * Fletcher [1987] R. Fletcher.  _Practical Methods of Optimization_.  Wiley-Interscience, 2nd edition, 1987.
  * Frankel [2011] T. Frankel.  _The Geometry of Physics: An Introduction_.  Cambridge, 3rd edition, 2011.
  * Goodfellow et al. [2016] I. Goodfellow, Y. Bengio, and A. Courville.  _Deep Learning_.  The MIT Press, 2016.
  * Kingma and Ba [2015] D.P. Kingma and J.L. Ba.  Adam: A method for stochastic optimization.  In _International Conference on Learning Representations (ICLR)_ , pages 1–13, 2015.
  * Lee [2013] J.M. Lee.  _Introduction to Smooth Manifolds_.  Springer, 2nd edition, 2013.
  * Lemaréchal [2010] C. Lemaréchal.  Cauchy and the gradient method, 2010.  URL math.uni-bielefeld.de/documents/vol-ismp/40_lemarechal-claude.pdf.
  * Liegeois [1977] A. Liegeois.  Automatic supervisory control of the configuration and behavior of multibody mechanisms.  _IEEE Transactions on Systems, Man and Cybernetics_ , 7(12):868–871, 1977.
  * Maciejewski [1990] A.A. Maciejewski.  Dealing with ill-conditioned equations of motion for articulated figures.  _IEEE Computer Graphics and Applications_ , 10(3):63–71, 1990.
  * Møller [1993a] M. Møller.  Exact calculation of the product of the Hessian matrix of feed-forward network error functions and a vector in O(N) time.  Technical Report DAIMI PB432, Computer Science Department, Aarhus University, Denmark, 1993a.
  * Møller [1993b] M. Møller.  A scaled conjugate gradient algorithm for fast supervised learning.  _Neural Networks_ , 6(4):525–533, 1993b.
  * Møller [1993c] M. Møller.  Supervised learning on large redundant training sets.  _International Journal of Neural Systems_ , 4(1), 1993c.
  * Møller [1993d] M. Møller.  _Efficient Training of Feedforward Neural Networks_.  PhD thesis, Århus Universitet, Danmark, 1993d.
  * Nakamura and Hanafusa [1986] Y. Nakamura and H. Hanafusa.  Inverse kinematic solutions with singularity robustness for robot manipulator control.  _Journal of Dynamic Systems, Measurement, and Control_ , 108(3):163–171, 1986.
  * Nesterov [1983] Y. Nesterov.  A method for unconstrained convex minimization problem with the rate of convergence O$(1/k^{2})$.  _Soviet Mathematics Doklady_ , 27:372–376, 1983.
  * Nesterov [2018] Y. Nesterov.  _Lectures on Convex Optimization_.  Springer, 2nd edition, 2018.
  * Nocedal and Wright [2006] J. Nocedal and S.J. Wright.  _Numerical Optimization_.  Springer, 2nd edition, 2006.
  * Pearlmutter [1994] B.A. Pearlmutter.  Fast exact multiplication by the Hessian.  _Neural Computation_ , 6(1):147–160, 1994.
  * Pedregosa et al. [2011] F. Pedregosa, , G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, v. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay.  Scikit-learn: Machine learning in Python.  _Journal of Machine Learning Research_ , 12:2825–2830, 2011\.
  * Polyak [1964] B.T. Polyak.  Some methods of speeding up the convergence of iteration methods.  _USSR Computational Mathematics and Mathematical Physics_ , 4(5):1–17, 1964.
  * Press et al. [2007] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery.  _Numerical Recipes, The Art of Scientific Computing_.  Cambridge University Press, 3rd edition, 2007.
  * Reddi et al. [2018] S.J. Reddi, S. Kale, and S. Kumar.  On the convergence of Adam and beyond.  In _International Conference on Learning Representations (ICLR)_ , 2018.
  * Ruder [2017] S. Ruder.  An overview of gradient descent optimization algorithms.  _arXiv preprint_ , arXiv:1609.04747, 2017.
  * Sciavicco and Siciliano [2001] L. Sciavicco and B. Siciliano.  _Modelling and Control of Robot Manipulators_.  CreateSpace Independent Publishing Platform, 2nd edition, 2001.
  * Scikit-learn [2020] Scikit-learn.  Digits dataset, 2020.  URL https://scikit-learn.org/stable/modules/generated/sklearn.datasets.load_digits.html.
  * Siciliano [1990] B. Siciliano.  A closed-loop inverse kinematic scheme for on-line joint-based robot control.  _Robotica_ , 8(7):231–243, 1990.
  * Ślebodziński [1931] W. Ślebodziński.  Sur les équations de Hamilton.  _Bull. Acad. Roy. d. Belg._ , 17(5):864–870, 1931\.
  * Sra et al. [2011] S. Sra, S. Nowozin, and S.J. Wright.  _Optimization for Machine Learning_.  The MIT Press, 2011.
  * Tieleman and Hinton [2012] T. Tieleman and G. Hinton.  Lecture 6.5 - RMSProp, Coursera: Neural networks for machine learning.  Technical report, 2012.
  * Tu [2011] L.W. Tu.  _An Introduction to Manifolds_.  Springer, 2nd edition, 2011.
  * Wampler [1986] C.W. Wampler.  Manipulator inverse kinematic solutions based on vector formulations and damped least-squares methods.  _IEEE Transactions on Systems, Man and Cybernetics_ , 16(1):93–101, 1986.
  * Whitney [1969] D.E. Whitney.  Resolved motion rate control of manipulators and human prostheses.  _IEEE Transactions on Man-Machine Systems_ , 10:47–53, 1969.
  * Wolfe [1969] P. Wolfe.  Convergence conditions for ascent methods.  _SIAM Review_ , 11(2):226–235, 1969.
  * Wolfe [1971] P. Wolfe.  Convergence conditions for ascent methods. II: Some corrections.  _SIAM Review_ , 13(2):185–188, 1971.
  * Zeiler [2012] M.D. Zeiler.  Adadelta: an adaptive learning rate method.  _arXiv preprint_ , arXiv:1212.5701, 2012.
  * Zimmer [2020] M.F. Zimmer.  Two direct solvers for a system of linear equations, 2020.  URL https://www.arxiv.org/pdf/1611.06633pdf.
  * Zimmer [2021a] M.F. Zimmer.  Github repository, 2021a.  URL www.github.com/mfzimmer.
  * Zimmer [2021b] M.F. Zimmer.  Neograd: Gradient descent with a near-ideal learning rate, 2021b.  URL https://www.arxiv.org/pdf/2010.07873.pdf.