# Stochastic Hessian Fitting on Lie Group

Xi-Lin Li

###### Abstract

This paper studies the fitting of Hessian or its inverse with stochastic
Hessian-vector products. A Hessian fitting criterion, which can be used to
derive most of the commonly used methods, e.g., BFGS, Gaussian-Newton,
AdaGrad, etc., is used for the analysis. Our studies reveal different
convergence rates for different Hessian fitting methods, e.g., sublinear rates
for gradient descent in the Euclidean space and a commonly used closed-form
solution, linear rates for gradient descent on the manifold of symmetric
positive definite (SPL) matrices and certain Lie groups. The Hessian fitting
problem is further shown to be strongly convex under mild conditions on a
specific yet general enough Lie group. To confirm our analysis, these methods
are tested under different settings like noisy Hessian-vector products, time
varying Hessians, and low precision arithmetic. These findings are useful for
stochastic second order optimizations that rely on fast, robust and accurate
Hessian estimations.

##  1 Introduction

This paper considers the problem of Hessian fitting for second order
stochastic optimizations where the objective functions are typically defined
as an expectation, say $f(\theta)=E_{z}[\ell(\theta,z)]$ with
$\theta\in\mathbb{R}^{m}$ being the parameters to be optimized, and $z$ a
random vector that can be sampled to evaluate an unbiased estimate of
$f(\theta)$. We always assume that the considered problem is second-order
differentiable such that a preconditioner $P\succ 0$, say the inverse of the
Hessian of $f(\theta)$ when it is convex, can be used to accelerate the
convergence of the following iteration,

$\theta_{t+1}=\theta_{t}-\mu_{t}P_{t}\,{\partial
f_{t}(\theta)}/{\partial\theta}\left.\right|_{\theta=\theta_{t}}$ (1)

to an optimal solution for $\theta$, where subscript $t$ is the time index,
and $f_{t}(\theta)$ a sampled noisy estimation of $f(\theta)$ at step $t$. To
simplify the notations, let us ignore the time index $t$ when doing so causing
no misunderstanding, and use $g={\partial f(\theta)}/{\partial\theta}$ and
$h={\partial(v^{T}g)}/{\partial\theta}=Hv$ to denote the stochastic gradient
and Hessian-vector product associated with a vector $v$ at a certain iteration
step, respectively, where
$H={\partial^{2}f(\theta)}/{(\partial\theta^{T}\partial\theta)}$ is the
Hessian and generally, we do not assume its definiteness. Then, with the
following preconditioner estimation criterion proposed in [Li18],

$c(P;v,h)=h^{T}Ph+v^{T}P^{-1}v$ (2)

we can derive most of the commonly used preconditioners as listed in Table 1
[Li19]. For example, solution $Phh^{T}P=vv^{T}$ reduces to the secant equation
$Ph=v$ when $v^{T}h>0$, which is used for deriving the BFGS like methods
[BV04]. Solution $P=(H^{2})^{-0.5}$ reduces to the Newton’s method when
$H\succ 0$, which uses preconditioner $P=H^{-1}$. Solution $P^{-2}=E[gg^{T}]$
reduces to the Fisher information matrix if $g$ is the sample-wise gradient,
which is used in the Gaussian-Newton and natural gradient methods [Ama98]. The
Fisher information matrix and Hessian are closely related when $\ell(\cdot)$
is a negative logarithm likelihood (NLL) loss from the exponential familiy
[Mur23]. Solution $P_{t}=\left(\sum_{i=1}^{t}g_{i}g_{i}^{T}\right)^{-0.5}$
reduces to the preconditioner used in AdaGrad [DHS11], a large family of
optimizers including many popular examples like the Adam(W), Shampoo, etc.
Thus, (2) is a reasonable choice for analyzing the performance of a wide range
of preconditioner fitting methods. Previously, the preconditioned stochastic
gradient descent (PSGD) related works always fit $P$ on certain Lie groups
[Li18, Li19, Li22], and recent empirical results have confirmed their
excellent performance for a handful of large scale machine learning related
stochastic optimization problems [PL24]. Practically, we know that these Lie
group preconditioner fitting methods can avoid a few problematic operations
inherent to other popular methods, e.g., line search required by the BFGS like
quasi-Newton methods in a stochastic setting, matrix inverse or matrix square
root inverse in the closed-form solutions listed in Table 1. However, these
works do not answer a few theoretically important questions like the
efficiency and necessity of fitting $P$ on those Lie groups, advantages over
other simpler methods like BFGS and the closed-form solutions listed in Table
1, disregarding their certain practical difficulties. One extraordinary
finding of this work is that criterion (2) is strongly convex on a certain
general enough Lie group under mild conditions, see Proposition 5 in Section
3.1 for details.

Notations in the paper are fairly standard. For example, $A\succ 0$ means that
$A$ is a symmetric positive definite (SPD) matrix, $E_{v}[\cdot]$ takes
expectation with respect to a certain random variable $v$, $\mathcal{N}(0,I)$
denotes the standard normal distribution, $A^{T}$ is the transpose of $A$,
$A^{0.5}$ denotes the principal square root of $A$, $A^{-2T}$ is the shorthand
for $(A^{-2})^{T}$, $\|A\|_{2}$ is the spectral norm of $A$, ${\rm tr}(A)$
takes the trace of $A$, $\lambda_{\min}(A)$ takes the minimum eigenvalue of
$A$, ${\rm vec}(A)$ vectorizes matrix $A$ by stacking its columns
successively, GL$(n,\mathbb{R})$ stands for the general linear group of
$n\times n$ invertible matrices, and $\otimes$ denotes Kronecker product. We
say that $A$ and $B$ commutes if $[A,B]=AB-BA=0$. For a symmetric matrix
$A\in\mathbb{R}^{n\times n}$, we only need to vectorize either its upper or
lower triangular part, and denote this operation as ${\rm vech}(A)$. Clearly,
${\rm vec}(A)=S\,{\rm vech}(A)$, where
$S\in\mathbb{R}^{n^{2}\times[0.5n(n+1)]}$ is a sparse full column rank matrix
with elements being $0$ or $1$.

Table 1: Variations of (2) for preconditioner fitting. We assume that $v\sim\mathcal{N}(0,I)$ in $E_{v}[\cdot]$.  Criterion | Solution for $P$ | Notes
---|---|---
$c(P;v,h)$ | $Phh^{T}P=vv^{T}$ |  $Ph=v$ if $v^{T}h>0$
$E_{v}[c(P;v,h)]$ | $P=(H^{2})^{-0.5}$ |  $P=H^{-1}$ if $H\succ 0$
$E_{v}[c(P;v,g)]$ | $P=\left(E[gg^{T}]\right)^{-0.5}$ |  $g$ is sample-wise gradient
$\sum_{i=1}^{t}E_{v_{i}}\left[c(P;v_{i},g_{i})\right]$ | $P_{t}=\left(\sum_{i=1}^{t}g_{i}g_{i}^{T}\right)^{-0.5}$ |  $g_{i}$ is stochastic gradient

##  2 Candidates for stochastic Hessian fitting

###  2.1 The problem setups

We are to fit the Hessian or its inverse sequentially, i.e., updating the
Hessian estimate once given a pair of vector and its associated Hessian-vector
product, $(v,h)$, with model,

$h=Hv+\epsilon$ (3)

where $H\in\mathbb{R}^{n\times n}$ is the unknown Hessian matrix to be
estimated, and $\epsilon$ denotes any possible modeling errors, e.g.,
truncation and quantization errors when the Hessian-vector product $h$ is
approximated with a finite-difference method, or the randomness error when a
stochastic yet unbiased proxy of the target function
$f(\theta)=E_{z}[\ell(\theta,z)]$ is adopted for Hessian-vector product
evaluation, or the structural errors caused by using a sparse matrix form to
fit a dense Hessian. We assume that $v$ and $\epsilon$ are independent.

We choose to update $P$ by minimizing the preconditioner fitting criterion in
(2). With $v$ drawn from the standard normal distribution, i.e.,
$v\sim\mathcal{N}(0,I)$, from (2), it is ready to show that

$\displaystyle c(P)=$ $\displaystyle E_{v\sim\mathcal{N}(0,I)}[c(P;v,h)]$
$\displaystyle=$ $\displaystyle{\rm
tr}\left(PE[hh^{T}]+P^{-1}E[vv^{T}]\right)$ $\displaystyle=$
$\displaystyle{\rm
tr}\left(P\left(H^{2}+E[\epsilon\epsilon^{T}]\right)+P^{-1}\right)$ (4)

Then, the optimal positive definite $P$ for minimizing $c(P)$ can be shown to
be

$P=\left(H^{2}+E[\epsilon\epsilon^{T}]\right)^{-0.5}$ (5)

by solving equation ${\partial c(P)}/{\partial P}=0$. Thus, $P$ is called the
preconditioner. One can replace the $P$ in (2) with $P^{-1}$ to obtain another
criterion for fitting the Hessian directly. Here, we only consider the fitting
criterion in (2) since the intention is to use $P$ as a preconditioner for
preconditioning the stochastic gradient descent (SGD) in (1).

As revealed in (5), the error term $\epsilon$ in (3) helps to damp the
preconditioner estimation and has been shown to regularize the optimizer in a
perfectly balanced way [Li18]. Use of $(H^{2})^{-0.5}$ as a preconditioner
directly could severely amplify the gradient noise [Li18]. Here, we will
ignore the term $\epsilon$ most of the time to simplify our performance
studies. Definiteness of the Hessian $H$ is important for optimization of the
target function $f(\theta)$, but not interested to our setups since (4) and
(5) suggest that signs of $\lambda(H)$ do not matter. Hence, without the loss
of generality, we always assume that $H$ is an invariant positive definite
matrix for the study. Just bear in mind that except for the BFGS, our
preconditioners derived from (2) always work regardless of the definiteness of
$H$.

Criterion (2) is a well defined function for any $P$ with $\det(P)\neq 0$,
i.e., the group GL$(n,\mathbb{R})$. For our purpose, we further require
$P\succ 0$ such that it is a valid preconditioner for the iterations given in
(1). This symmetry constraint, i.e., $P=P^{T}$, will slightly complicate the
calculations of the derivatives of $c(P;v,h)$ with respect to $P$. We will
bring up this subtleness explicitly when it matters. From now on, by default,
the gradient and Hessian always refer to the ones of $c(P;v,h)$ or $c(P)$ with
respect to $P$. Gradient descent (GD) or SGD also refers to the iterations for
updating $P$, not $\theta$, if not clarified.

###  2.2 Expansion of criterion (2)

Let $dP$ denote a small perturbation of $P$. We are to expand the criterion in
(2) up to its second order terms of $dP$. To begin with, we expand
$(P+dP)^{-1}$ around $P^{-1}$ as below

$\displaystyle dP^{-1}=$ $\displaystyle(P+dP)^{-1}-P^{-1}$ $\displaystyle=$
$\displaystyle(I-P^{-1}dP+P^{-1}dP\,P^{-1}dP+\mathcal{O}((dP)^{3}))P^{-1}-P^{-1}$
$\displaystyle=$
$\displaystyle-P^{-1}dP\,P^{-1}+P^{-1}dP\,P^{-1}dP\,P^{-1}+\mathcal{O}((dP)^{3})$

where $\mathcal{O}((dP)^{3})$ includes the third and higher order terms of
$dP$ that are negligible when $dP$ is small enough. Now, with the above
expression, we can expand (2) as,

$\displaystyle dc(P;v,h)=$ $\displaystyle h^{T}dP\,h+v^{T}dP^{-1}\,v$
$\displaystyle=$ $\displaystyle
h^{T}dP\,h-v^{T}P^{-1}dP\,P^{-1}v+v^{T}P^{-1}dP\,P^{-1}dP\,P^{-1}v+\mathcal{O}((dP)^{3})$
(6)

To maximize the reduction of $c(P;v,h)$, we truncate the $dc(P;v,h)$ in (6) up
to its quadratic term of $dP$, and solve for $dP$ from the following equation,

$hh^{T}-P^{-T}vv^{T}P^{-T}+P^{-T}vv^{T}P^{-T}dP^{T}\,P^{-T}+P^{-T}dP^{T}\,P^{-T}vv^{T}P^{-T}=0$
(7)

From (7), it is ready to identify the gradient of $c(P;v,h)$ with respect to
$P$ as below,

$\frac{\partial c(P;v,h)}{\partial P}=hh^{T}-{P}^{-T}vv^{T}{P}^{-T}$ (8)

Note that (8) gives the steepest ascent direction in the space
$\mathbb{R}^{n\times n}$. With constraint $P=P^{T}$, we should double the off-
diagonal elements on the left side of (8) to obtain the steepest ascent
direction in the space $\mathbb{R}^{0.5n(n+1)}$. We can choose either gradient
for the update of $P$. But, we do need to calculate the Hessian with respect
to $P$ in space $\mathbb{R}^{0.5n(n+1)}$. We typically assume that
$v\sim\mathcal{N}(0,I)$, and consider the expectation of criterion (2), i.e.,
the $c(P)=E_{v}[c(P;v,h]$ given in (4), for performance analysis. We have the
following convex property of $c(P)$.

_Proposition 1:_ $c(P)$ is a convex function of $P$ on set $\\{P|P\succ 0\\}$.

_Proof:_ From (6), we can write the second order term of $dc(P)$ as

$\displaystyle\quad
E_{v\sim\mathcal{N}(0,I)}\left[v^{T}P^{-1}dP\,P^{-1}dP\,P^{-1}v\right]$
$\displaystyle={\rm
tr}\left(E_{v}\left[P^{-1}dP\,P^{-1}dP\,P^{-1}vv^{T}\right]\right)$
$\displaystyle={\rm tr}\left(P^{-1}dP\,P^{-1}dP\,P^{-1}\right)$
$\displaystyle=\left({\rm
vec}\left(\left(P^{-1}dP\right)^{T}\right)\right)^{T}{\rm
vec}\left(P^{-1}dP\,P^{-1}\right)$ $\displaystyle=\left({\rm
vec}\left(dP\,P^{-1}\right)\right)^{T}{\rm vec}\left(P^{-1}dP\,P^{-1}\right)$
$\displaystyle=\left(\left(P^{-1}\otimes I\right){\rm
vec}(dP)\right)^{T}\left(\left(P^{-1}\otimes P^{-1}\right){\rm
vec}(dP)\right)$ $\displaystyle=\left({\rm
vec}(dP)\right)^{T}\left(P^{-1}\otimes I\right)\left(P^{-1}\otimes
P^{-1}\right){\rm vec}(dP)$ $\displaystyle=\left({\rm
vec}(dP)\right)^{T}\left(P^{-2}\otimes P^{-1}\right){\rm vec}(dP)$
$\displaystyle=\left({\rm vech}(dP)\right)^{T}S^{T}\left(P^{-2}\otimes
P^{-1}\right)S\,{\rm vech}(dP)$ (9)

where we have used the relation ${\rm tr}(A^{T}B)=\left({\rm
vec}(A)\right)^{T}{\rm vec}(B)$ to arrive the fourth line from the third one,
${\rm vec}(ABC)=(C^{T}\otimes A){\rm vec}(B)$ from the fifth to the sixth
line, and $(AB)\otimes(CD)=(A\otimes C)(B\otimes D)$ from the third to the
last to the next line. Since $P^{-2}\otimes P^{-1}\succ 0$ and $S$ has full
column rank, the Hessian with respect to $P$, i.e., $S^{T}\left(P^{-2}\otimes
P^{-1}\right)S$ from (9), is positive definite. $\square$

Note that following the proof of Proposition 1, we also can show that
$c(P;v,h)$ is weakly convex. In general, the quadratic term in (6) can be
negative without constraint $dP=dP^{T}$ even when $P\succ 0$. Thus, the
symmetry constraint is essential for the convexity of $c(P)$.

###  2.3 Hessian fitting in the Euclidean space

####  2.3.1 Gradient descent

One obvious choice is to update $P$ with the gradient given in (8). Let us
start from an initial guess $P_{0}\succ 0$, and update $P$ at step $t$ as

$P_{t+1}=P_{t}-\mu(h_{t}h_{t}^{T}-{P}_{t}^{-1}v_{t}v_{t}^{T}{P}_{t}^{-1})$
(10)

where $\mu$ is a small enough positive step size,
$(v_{t},h_{t})=(v_{t},Hv_{t})$ with $v_{t}\sim\mathcal{N}(0,I)$ is a pair of
vector and its associated Hessian-vector product at step $t$. Note that
$P_{t}$ is always symmetric as long as $P_{0}$ is.

From the proof of Proposition 1, we see that the Hessian of $c(P)$ with
respect to $P$ is neither lower nor upper tightly bounded. Thus, criterion
$c(P)$ satisfies neither the Polyak-Łojasiewicz nor the Lipschitz continuity
condition. Still, following the regret analysis [DHS11], we can show that with
a proper constant step size $\mu$, the sequence given by (10) converges to
$H^{-1}$ at least sublinearly. We only list the key steps here. Assuming that
$\mu$ is small enough such that we always have $P_{t}\succ 0$ as long as
$P_{0}\succ 0$. Then, via the convexity of $c(P;v,h)$, we have

$\displaystyle\|P_{t+1}-H^{-1}\|_{F}^{2}=$
$\displaystyle\left\|P_{t}-\mu\frac{\partial c(P_{t};v_{t},h_{t})}{\partial
P_{t}}-H^{-1}\right\|_{F}^{2}$ $\displaystyle=$
$\displaystyle\left\|P_{t}-H^{-1}\right\|_{F}^{2}+\mu^{2}\left\|\frac{\partial
c(P_{t};v_{t},h_{t})}{\partial P_{t}}\right\|_{F}^{2}+2\mu\,{\rm
tr}\left((H^{-1}-P_{t})\frac{\partial c(P_{t};v_{t},h_{t})}{\partial
P_{t}}\right)$ $\displaystyle\leq$
$\displaystyle\left\|P_{t}-H^{-1}\right\|_{F}^{2}+\mu^{2}\left\|\frac{\partial
c(P_{t};v_{t},h_{t})}{\partial
P_{t}}\right\|_{F}^{2}+2\mu\left(c(H^{-1};v_{t},h_{t})-c(P_{t};v_{t},h_{t})\right)$

where the last line is due to the convexity of $c(P;v,h)$. Now, let us
rearrange the above inequality to have

$c(P_{t};v_{t},h_{t})-c(H^{-1};v_{t},h_{t})\leq\frac{\left\|P_{t}-H^{-1}\right\|_{F}^{2}-\left\|P_{t+1}-H^{-1}\right\|_{F}^{2}}{2\mu}+\frac{\mu}{2}\left\|\frac{\partial
c(P_{t};v_{t},h_{t})}{\partial P_{t}}\right\|_{F}^{2}$

Lastly, we take expectation on both sides of the above equation with respect
to $v$, and sum over index $t$ to have

$c\left(\sum_{i=1}^{t}P_{i}/t\right)-c(H^{-1})=\mathcal{O}(1/\sqrt{t})$

by setting $\mu\propto 1/\left(\sqrt{t}\,\max_{1\leq i\leq
t}\left\|\frac{\partial c(P_{i})}{\partial P_{i}}\right\|_{F}\right)$. Still,
the above analysis does not exclude possible better convergence results with
variant step sizes.

####  2.3.2 The Newton’s method

Another choice is to take the expectation on the left side of (7) to have

$H^{2}-P^{-2}+P^{-2}dP\,P^{-1}+P^{-1}dP\,P^{-2}=0$ (11)

and then solve for the optimal $dP$ directly as

$\displaystyle{\rm vec}(dP)=$ $\displaystyle\left(P^{-1}\otimes
P^{-2}+P^{-2}\otimes P^{-1}\right)^{-1}{\rm vec}\left(P^{-2}-H^{2}\right)$
$\displaystyle=$ $\displaystyle\left(I\otimes P^{-1}+P^{-1}\otimes
I\right)^{-1}(P\otimes P){\rm vec}\left(P^{-2}-H^{2}\right)$ $\displaystyle=$
$\displaystyle\left(I\otimes P^{-1}+P^{-1}\otimes I\right)^{-1}{\rm
vec}\left(I-PH^{2}P\right)$ (12)

The above equation gives a numerically stable but costly solution for $dP$.
But, if we assume that $P$ and $dP$ commutes such that
$P^{-2}dP\,P^{-1}=P^{-1}dP\,P^{-2}$, then (11) immediately gives the following
solution,

$dP=-0.5P(H^{2}-P^{-2})P^{2}=0.5(P-PH^{2}P^{2})$ (13)

Now, (13) suggests the following Newton’s method for updating $P$,

$P_{t+1}=1.5P_{t}-0.5P_{t}H^{2}P_{t}^{2}$ (14)

_Proposition 2:_ Starting from an initial guess $P_{0}\succ 0$ that commutes
with $H$ and satisfies condition $0<\lambda(HP)<(\sqrt{17}-1)/2$, the Newton’s
method in (14) converges to $H^{-1}$ quadratically.

_Proof:_ First, let us show that $P_{t+1}$ commutes with $H$ if $P_{t}$ does.
With (14), we have

$\displaystyle P_{t+1}H-HP_{t+1}=$ $\displaystyle
1.5(P_{t}H-HP_{t})-0.5(P_{t}H^{2}P_{t}^{2}H-HP_{t}H^{2}P_{t}^{2})$
$\displaystyle=$ $\displaystyle-0.5(P_{t}H^{3}P_{t}^{2}-P_{t}H^{3}P_{t}^{2})$
$\displaystyle=$ $\displaystyle 0$

Hence, $P_{t}$ always commutes with $H$ as long as $P_{0}$ does. Also, we can
easily show that $P_{t}$ is symmetric as long as $P_{0}$ is.

Now, let us define $R_{t}=HP_{t}-I$. Similarly, we can show that $P_{t}$, $H$
and $R_{t}$ all commute with each other and are symmetric. Then, by (14), we
have

$\displaystyle R_{t+1}=$ $\displaystyle 1.5HP_{t}-0.5HP_{t}H^{2}P_{t}^{2}-I$
$\displaystyle=$ $\displaystyle 1.5HP_{t}-0.5(HP_{t})^{3}-I$ $\displaystyle=$
$\displaystyle 1.5(R_{t}+I)-0.5(R_{t}+I)^{3}-I$ $\displaystyle=$
$\displaystyle-0.5(R_{t}+3I)R_{t}^{2}$ (15)

We require

$\|0.5(R_{t}+3I)R_{t}\|_{2}<1$

to make (15) become a contraction mapping such that
$\|R_{t+1}\|_{2}<\|R_{t}\|_{2}$, which suggests that the eigenvalues of
$P_{t}H$ should be bounded by

$0<\lambda(P_{t}H)<\frac{\sqrt{17}-1}{2}$ (16)

Then, with the bounds given in (16) and the last line of (15), we can show the
following quadratic convergence,

$\frac{\|R_{t+1}\|_{2}}{\|R_{t}\|_{2}^{2}}=0.5\|R_{t}+3I\|_{2}<\frac{\sqrt{17}+3}{4}$

Note that the lower bound of $\lambda(P_{t}H)$ in (16) already implies that
$H$ cannot be singular such that $H^{-1}$ exists. $\square$

It is trivial to choose a starting point $P_{0}$ that commutes with $H$, e.g.,
a $P_{0}\propto I$. However, in practice, the Newton’s method given in (14) is
vulnerable to round-off errors that eventually break the commuting property
between $P_{t}$ and $H$. Also, the Newton’s method require perfect knowledge
of $H^{2}$. We present it here mainly for theoretical interests.

####  2.3.3 A few closed-form solutions

For the sake of completeness, let us also briefly discuss the convergence rate
of a commonly used closed-form solution,
$P_{t}=\left(\left(P_{0}^{-2}+\sum_{i=1}^{t}h_{i}h_{i}^{T}\right)/(t+1)\right)^{-0.5}$,
which can be put in a more fancier form as below

$P_{t+1}=\left(\frac{t+1}{t+2}P_{t}^{-2}+\frac{1}{t+2}h_{t}h_{t}^{T}\right)^{-0.5}$
(17)

As $h_{t}=Hv_{t}$ and $v_{t}\sim\mathcal{N}(0,I)$, we have

$\displaystyle P_{t}=$
$\displaystyle\left(H\left(P_{0}^{-2}/(t+1)+\sum_{i=1}^{t}v_{i}v_{i}^{T}/(t+1)\right)H\right)^{-0.5}$
$\displaystyle=$
$\displaystyle\left(H\left(I+\mathcal{O}(1/t)\right)H\right)^{-0.5}$
$\displaystyle=$ $\displaystyle H^{-1}+\mathcal{O}(1/t)$

Hence, the closed-form solution given in (17) converges to $H^{-1}$
sublinearly.

Another closed-formed solution can be obtained from (7) by solving for $P$
directly after setting $dP=0$. Specifically, (7) suggests that a stationary
point $P$ satisfies

$P\left(\sum_{i=1}^{t}h_{i}h_{i}^{T}\right)P-\sum_{i=1}^{t}v_{i}v_{i}^{T}=0$
(18)

Then, we can solve the Riccati equation in (18) for $P$ directly when both
$\sum_{i=1}^{t}v_{i}v_{i}^{T}$ and $\sum_{i=1}^{t}h_{i}h_{i}^{T}$ are not
singular, see the proof of Proposition 3 in [Li18]. Similar to the Newton’s
method, this closed-form solution is mainly for theoretical interests

One more closed-form solution worth to mention is the BFGS like methods that
update $P_{t}$ to $P_{t+1}$ with a rank-2 update such that
$P_{t+1}h_{t}=v_{t}$ [BV04]. These methods have two limitations: the pair
$(v_{t},h_{t})$ needs to satisfy the curvature condition $v_{t}^{T}h_{t}>0$;
they diverge easily with noisy Hessian-vector products.

###  2.4 Fitting on the manifold of SPD

One possible way to define a local coordinate for $P$ on the manifold of SPD
is to let $dP$ be

$dP=P\mathcal{E}+\mathcal{E}P^{T}$ (19)

where $\mathcal{E}$ is a sufficiently small matrix, but not necessarily
symmetric in general. Clearly, $dP$ is symmetric if $\mathcal{E}$ is,
regardless of the symmetry property of $P$. On the other hand, we can show
that $\mathcal{E}$ is symmetric if $P\succ 0$.

_Proposition 3:_ Eq. (19) induces a new coordinate system for manifold $P$
with $P\succ 0$.

_Proof:_ We are to show that (19) defines a bijective mapping between $dP$ and
$\mathcal{E}$ when $P\succ 0$, which is partially done after we vectorize both
sides of (19) to solve for $\mathcal{E}$ as,

${\rm vec}(\mathcal{E})=(I\otimes P+P\otimes I)^{-1}{\rm vec}(dP)$ (20)

Still, we need to show that (20) preserves the symmetry properties of $dP$ and
$\mathcal{E}$. Clearly, a symmetric $\mathcal{E}$ maps to a symmetric $dP$ by
(19). What left is to show that (20) also maps a symmetric $dP$ to a symmetric
$\mathcal{E}$. To begin with, let us introduce a permutation matrix $\Pi$ such
that ${\rm vec}(A^{T})=\Pi\,{\rm vec}(A)$. Then, for any matrix $X$ with the
same size as $P$, we have

$\displaystyle(P\otimes I+I\otimes P)\Pi\,{\rm vec}(X)=$
$\displaystyle(P\otimes I+I\otimes P){\rm vec}(X^{T})$ $\displaystyle=$
$\displaystyle{\rm vec}(X^{T}P+PX^{T})$ $\displaystyle=$ $\displaystyle{\rm
vec}\left((PX+XP)^{T}\right)$ $\displaystyle=$ $\displaystyle\Pi{\rm
vec}(PX+XP)$ $\displaystyle=$ $\displaystyle\Pi(I\otimes P+P\otimes I){\rm
vec}(X)$

Since $X$ is arbitrary, the above equation suggests that

$\Pi^{-1}(P\otimes I+I\otimes P)=(I\otimes P+P\otimes I)\Pi^{-1}$ (21)

Then, from (20), we have

$\displaystyle{\rm vec}(\mathcal{E}^{T})=$ $\displaystyle\Pi\,{\rm
vec}(\mathcal{E})$ $\displaystyle=$ $\displaystyle\Pi(I\otimes P+P\otimes
I)^{-1}{\rm vec}(dP)$ $\displaystyle=$ $\displaystyle\left((I\otimes
P+P\otimes I)\Pi^{-1}\right)^{-1}{\rm vec}(dP)$ $\displaystyle=$
$\displaystyle\left(\Pi^{-1}(I\otimes P+P\otimes I)\right)^{-1}{\rm vec}(dP)$
$\displaystyle=$ $\displaystyle(I\otimes P+P\otimes I)^{-1}\Pi{\rm vec}(dP)$
$\displaystyle=$ $\displaystyle(I\otimes P+P\otimes I)^{-1}{\rm vec}(dP^{T})$
$\displaystyle=$ $\displaystyle(I\otimes P+P\otimes I)^{-1}{\rm vec}(dP)$
$\displaystyle=$ $\displaystyle{\rm vec}(\mathcal{E})$

where we have used (21) from the third to the fourth line. Now, it is clear
from the above equation that $\mathcal{E}=\mathcal{E}^{T}$ when $P\succ 0$.
$\square$

Hence, (19) induces a new coordinate system for $P\succ 0$. For a positive
scalar $P$, integral on the right side of (20) gives us the closed-form
transform from coordinate $P$ to the new one, i.e., $0.5\log P+{\rm
constant}$. In general, we have no need to write down this transform
explicitly for a matrix $P$, and (19) is enough for our purpose.

Now, we can rewrite the $dc(P;v,h)$ in (6) in the new coordinate as below,

$\displaystyle dc(P;v,h)=$ $\displaystyle
h^{T}(P\mathcal{E}+\mathcal{E}P^{T})h-v^{T}P^{-1}(P\mathcal{E}+\mathcal{E}P^{T})P^{-1}v$
$\displaystyle+v^{T}P^{-1}(P\mathcal{E}+\mathcal{E}P^{T})P^{-1}(P\mathcal{E}+\mathcal{E}P^{T})P^{-1}v+\mathcal{O}(\mathcal{E}^{3})$
$\displaystyle=$ $\displaystyle
h^{T}P\mathcal{E}h+h^{T}\mathcal{E}Ph-v^{T}\mathcal{E}P^{-1}v-v^{T}P^{-1}\mathcal{E}v$
$\displaystyle+v^{T}\mathcal{E}^{2}P^{-1}v+v^{T}\mathcal{E}P^{-1}\mathcal{E}v+v^{T}P^{-1}\mathcal{E}P\mathcal{E}P^{-1}v+v^{T}P^{-1}\mathcal{E}^{2}v+\mathcal{O}(\mathcal{E}^{3})$
(22)

where in the last line we have assumed that $P=P^{T}$. After ignoring the
higher order terms of $\mathcal{E}$ in (22), we can solve for the
$\mathcal{E}$ that reduces $dc(P;v,h)$ the most from the following equation,

$\displaystyle Phh^{T}+hh^{T}P-vv^{T}P^{-1}-P^{-1}vv^{T}+$ $\displaystyle
vv^{T}P^{-1}\mathcal{E}^{T}+\mathcal{E}^{T}vv^{T}P^{-1}+vv^{T}\mathcal{E}^{T}P^{-1}+P^{-1}\mathcal{E}^{T}vv^{T}+$
$\displaystyle
P^{-1}vv^{T}P^{-1}\mathcal{E}^{T}P+P\mathcal{E}^{T}P^{-1}vv^{T}P^{-1}+P^{-1}vv^{T}\mathcal{E}^{T}+\mathcal{E}^{T}P^{-1}vv^{T}=0$
(23)

Now, given a pair $(v_{t},h_{t})$ at step $t$, (23) suggests the following SGD
in the new coordinate,

$\displaystyle\mathcal{E}_{t}=$
$\displaystyle-\mu(P_{t}h_{t}h_{t}^{T}+h_{t}h_{t}^{T}P_{t}-v_{t}v_{t}^{T}P_{t}^{-1}-P_{t}^{-1}v_{t}v_{t}^{T})$
$\displaystyle P_{t+1}=$ $\displaystyle
P_{t}+(P_{t}\mathcal{E}_{t}+\mathcal{E}_{t}P_{t})$ (24)

Clearly, with (24), $P_{t}$ always is symmetric as long as $P_{0}$ is.

_Proposition 4:_ Starting from an initial guess $P_{0}\succ 0$ that commutes
with $H$ and with a constant step size
$0<\mu<0.5/\max_{t}{\lambda_{\max}(H+H^{2}P_{t})}$, the sequence $P_{t}$ given
by (24) converges to $H^{-1}$ linearly with rate $1-8\mu\,\lambda_{\min}(H)$.

_Proof:_ We take expectation on both sides of (24) to have

$\displaystyle E[P_{t+1}]=$ $\displaystyle
E[P_{t}]-\mu\left(E[P_{t}^{2}]H^{2}+2E[P_{t}H^{2}P_{t}]+H^{2}E[P_{t}^{2}]-4I\right)$
$\displaystyle=$ $\displaystyle
E[P_{t}]-\mu\left(\left(E[P_{t}]\right)^{2}H^{2}+2E[P_{t}]H^{2}E[P_{t}]+H^{2}\left(E[P_{t}]\right)^{2}-4I+\mathcal{O}(\left(P_{t}-E[P_{t}]\right)^{2})\right)$
(25)

where $E[\cdot]$ is a shorthand for
$E_{v_{1}\sim\mathcal{N}(0,I),\ldots,v_{t}\sim\mathcal{N}(0,I)}[\cdot]$. The
term $\mathcal{O}(\left(P_{t}-E[P_{t}]\right)^{2})$ in the second line of (25)
is negligible with a small $\mu$. We further drop the notation $E[\cdot]$ in
(25), and just use $P_{t}$ to denote its expectation to simplify our writings.

First, we are to show that $P_{t+1}$ commutes with $H$ if $P_{t}$ does.
Assuming $P_{t}H=HP_{t}$, we can simplify (25) to

$P_{t+1}=P_{t}-4\mu(H^{2}P_{t}^{2}-I)$ (26)

Then, it is not difficult to show that $P_{t+1}$ also commutes with $H$,
similar to the proof of Proposition 2. Thus, $P_{t}$ always commutes with $H$
as long as $P_{0}$ does. Second, let us define fitting error $R_{t}=HP_{t}-I$
and show that $\|R_{t}\|_{2}\rightarrow 0$ for large $t$. Clearly, $H$,
$P_{t}$ and $R_{t}$ all are symmetric and commutes with each other. Now, from
(26), we have

$\displaystyle R_{t+1}=$ $\displaystyle HP_{t+1}-I$ $\displaystyle=$
$\displaystyle HP_{t}-4\mu H\left((HP_{t})^{2}-I\right)-I$ $\displaystyle=$
$\displaystyle(R_{t}+I)-4\mu H\left((R_{t}+I)^{2}-I\right)-I$ $\displaystyle=$
$\displaystyle(I-8\mu H-4\mu HR_{t})R_{t}$ $\displaystyle=$
$\displaystyle(I-4\mu H-4\mu H^{2}P_{t})R_{t}$ (27)

From the last line of (27), we see that $\mu$ needs to satisfy

$0<\mu<\frac{1}{2\lambda_{\max}(H+H^{2}P_{t})}$

to ensure the monotonic convergence of $P_{t}$. The second to the last line of
(27) gives the asymptotic convergence rate as

$1>\frac{\|R_{t+1}\|_{2}}{\|R_{t}\|_{2}}=\|I-8\mu H\|_{2}\geq
1-8\mu\,\lambda_{\rm min}(H)\quad{\rm for}\;\|R_{t}\|_{2}\rightarrow 0$

Thus, $P_{t}$ converges linearly to $H^{-1}$ when $\lambda_{\rm min}(H)>0$.
$\square$

Again, similar to the proof of Proposition 2, the condition that $P$ and $H$
commutes simplifies the proof, but could be too strong to be satisfied in
practice. Yet, the linear convergence rate $1-8\mu\,\lambda_{\min}(H)$ matches
well with observations. Following a similar way for the proof of Proposition
1, from (22), we can show that the Hessian in the new coordinates is

$S^{T}\left[3(P^{-1}\otimes I)+P^{-2}\otimes P\right]S$

which again is neither lower nor upper bounded. Still, unlike the SGD in the
Euclidean space in (10), Proposition 4 shows that SGD on the manifold of SPL
could converge linearly to $H^{-1}$ with a proper constant step size.
Intuitively, from (19), we see that $dP$ is roughly proportional to $P$ to
mimic a Euclidean space varying step size adaptation.

It is interesting to check the Newton algorithm in the new coordinates as
well. Following the developments of the Newton’s method in the Euclidean
space, we can take expectation of the left side of (23) and solve for the
optimal $\mathcal{E}$ directly. Generally, the two Newton’s methods are
different as the coordinate transform defined by (19) are not linear. However,
if we assume that $\mathcal{E}$ and $P$ commutes, we can readily show that
$dP$ and $P$ commutes as well due to (19), and eventually the two Newton’s
methods are the same in this case.

A connected Lie group is a special manifold equipped with extra properties
that can further facilitate the development of Hessian fitting methods. We
leave this topic to the next section.

###  2.5 Summary

We summarize all the above candidates for Hessian fitting into Table 2. Fig. 1
shows a set of typical convergence curves when $H$ is a small Hilbert matrix.
The empirical convergence rates match well with the predictions in Table 2.

Table 2: Methods and convergence rates for Hessian fitting. Note that $0<\rho<1$ is not necessarily the same for different methods.  Method | Equation | Space | Err. at $t$
---|---|---|---
SGD | (10) | Euclidean | $\mathcal{O}(1/\sqrt{t})$
Closed-form | (17) | Euclidean | $\mathcal{O}(1/t)$
SGD | (24) | SPD | $\mathcal{O}(\rho^{t})$
SGD | (36) | Lie group | $\mathcal{O}(\rho^{t})$
Newton | (14) | Euclidean/SPD | $\mathcal{O}(\rho^{2^{t}})$
Figure 1: Typical convergence curves for the five Hessian fitting methods in
Table 2 when $H$ is a $3\times 3$ Hilbert matrix, i.e., $H_{i,j}=1/(i+j-1)$
with $1\leq i,j\leq n$. Except for the Newton method that starts from
$P_{0}=0.02I$, all the other methods start from $P_{0}=I$. The step sizes of
SGD in the Euclidean space and SPD manifold are tuned to achieve the maximum
rate of convergence. Except for the SGD (10) in the Euclidean space and the
closed-form solution (17) that converges sublinearly, all the other methods
can converge to $H^{-1}$ within machine precision after reasonable numbers of
iterations.

##  3 Hessian fitting on Lie groups

###  3.1 Strong convexity on Lie group

It is still not handy to fit the Hessian on the manifold of SPL as we need to
constantly guard the positive definiteness of $P$. Also, knowledge of $H$ is
required to determine the step size according to Proposition 4. In this
section, we let

$P=Q^{T}Q$ (28)

where $Q$ is on a certain Lie group. When $Q$ is a triangular matrix, (28) is
known as the Cholesky decomposition of $P$. But, in general, $Q$ can take many
possible forms. It is natural to come up with (28) since with coordinate
change $\vartheta=Q^{-T}\theta$, the preconditioned SGD update in (1) becomes
a standard SGD in the new coordinate of $\vartheta$. The group structure helps
to preserve certain invariance after transform $\theta\rightarrow\vartheta$,
e.g., orientation with group GL${}^{+}(n,\mathbb{R})$, vector length with
group O$(n)$, etc. Such a constraint actually helps to regularize the update
of $Q$ as shown in our following analysis.

With the Lie group constraint, we let $I+\mathcal{E}=e^{\varepsilon G}$, where
$\varepsilon\rightarrow 0$ and $G$ is the group generator. Then, we can have
the following two choices for $dQ$

$\displaystyle{\rm choice\,I:}\quad$ $\displaystyle dQ=e^{\varepsilon
G}Q-Q=\mathcal{E}Q$ $\displaystyle{\rm choice\,II:}\quad$ $\displaystyle
dQ=Qe^{\varepsilon G}-Q=Q\mathcal{E}$ (29)

Correspondingly, the $dP$ in (28) will have the following two possible forms,

$\displaystyle{\rm choice\,I:}\quad dP=$
$\displaystyle(Q+\mathcal{E}Q)^{T}(Q+\mathcal{E}Q)-P$ $\displaystyle=$
$\displaystyle Q^{T}(\mathcal{E}+\mathcal{E}^{T}+\mathcal{E}^{T}\mathcal{E})Q$
$\displaystyle{\rm choice\,II:}\quad dP=$
$\displaystyle(Q+Q\mathcal{E})^{T}(Q+Q\mathcal{E})-P$ $\displaystyle=$
$\displaystyle\mathcal{E}^{T}P+P\mathcal{E}+\mathcal{E}^{T}P\mathcal{E}$ (30)

The second choice of $dP$ in (30) looks almost the same as that in (19) after
ignoring the second order term of $\mathcal{E}$. But, now $\mathcal{E}$ is a
small matrix proportional to group generator, and no longer required to be
symmetric. Substituting (30) back to (6) and noting that
$P^{-1}=Q^{-1}Q^{-T}$, we arrive at the following two possible choices for
$dc(P;v,h)$,

$\displaystyle{\rm choice\,I:}\quad dc(P;v,h)=$ $\displaystyle
h^{T}Q^{T}(\mathcal{E}+\mathcal{E}^{T})Qh-v^{T}Q^{-1}(\mathcal{E}+\mathcal{E}^{T})Q^{-T}v$
$\displaystyle+h^{T}Q^{T}\mathcal{E}^{T}\mathcal{E}Qh+v^{T}Q^{-1}(\mathcal{E}^{2}+\mathcal{E}\mathcal{E}^{T}+\mathcal{E}^{2T})Q^{-T}v+\mathcal{O}(\mathcal{E}^{3})$
$\displaystyle{\rm choice\,II:}\quad dc(P;v,h)=$ $\displaystyle
h^{T}(\mathcal{E}^{T}P+P\mathcal{E})h-v^{T}(P^{-1}\mathcal{E}^{T}+\mathcal{E}P^{-1})v$
$\displaystyle+h^{T}\mathcal{E}^{T}P\mathcal{E}h+v^{T}(\mathcal{E}^{2}P^{-1}+\mathcal{E}P^{-1}\mathcal{E}^{T}+P^{-1}\mathcal{E}^{2T})v+\mathcal{O}(\mathcal{E}^{3})$
(31)

The first choice in (31) not only looks more balanced and symmetric, but also
have nicer properties than the second one as will be revealed in our following
developments. Thus, we mainly focus on the first choice.

One striking difference between fitting Hessian on the Lie group and the
manifold of SPL is that criterion $c(P)$ could be strongly convex on the Lie
group under certain mild conditions like $\lambda_{\min}(H)>0$. Clearly,
$c(P)$ is not convex on GL$(n,\mathbb{R})$ as $Q$ and $UQ$ give the same value
for $c(P)$ for any orthogonal matrix $U$. But, $c(P)$ could behave better
after properly aggregating out this rotation ambiguity. To this end, we need
to use the polar decomposition $A=U\Gamma$ to define a proper Lie group that
is irrelevant of such rotations, where $A\in\mathbb{R}^{n\times n}$ is a
square matrix, $U$ is an orthogonal matrix, and $\Gamma\succeq 0$ is semi-
positive definite. With the polar decomposition, let us introduce the
following equivalence relation.

_Equivalence relation $R_{\rm polar}$: _For $Q_{1}$ and $Q_{2}$ in
$\mathbb{R}^{n\times n}$, we say that $Q_{1}\sim Q_{2}$ if and only if their
polar decompositions share the same $\Gamma$, i.e., $Q_{1}=U_{1}\Gamma$ and
$Q_{2}=U_{2}\Gamma$ for certain orthogonal matrices $U_{1}$ and $U_{2}$.

Then we can define the quotient set GL$(n,\mathbb{R})/R_{\rm polar}$. It is
not difficult to show that all the elements of GL$(n,\mathbb{R})/R_{\rm
polar}$ forms a Lie group. We skip the proof. Note that since group $O(n)$ has
dimension $0.5n(n-1)$, elements in this group have dimension
$n^{2}-0.5n(n-1)=0.5n(n+1)$. However, we cannot represent elements from this
group as positive definite symmetric matrices since set $\\{A|A\succ 0\\}$ is
not closed with matrix multiplication operation. Abstract mathematical
properties of this group is not interested here.

_Proposition 5:_ $c(P)$ is strongly convex on group GL$(n,\mathbb{R})/R_{\rm
polar}$ with the coordinate system induced by choice I in (29) when
$\lambda_{\min}(H)\geq\lambda_{0}>0$.

_Proof:_ First, let us briefly show that with $A\succ 0$ and eigenvalue
decomposition $A=UDU^{T}$,
$aA+bA^{-1}-2\sqrt{ab}I=U(aD+bD^{-1}-2\sqrt{ab})U^{T}\succeq 0$, where $a$ and
$b$ are two positive numbers.

Second, we show that it is sufficient to consider a symmetric $\mathcal{E}$ on
group GL$(n,\mathbb{R})/R_{\rm polar}$. Starting from $Q$, let
$Q_{1}=Q+\mathcal{E}_{1}Q$ be a neighbor of $Q$ on GL$(n,\mathbb{R})/R_{\rm
polar}$ with an arbitrary $\mathcal{E}_{1}$ small enough such that
$\|\mathcal{E}_{1}\|_{2}<1$. Let the polar decomposition of $Q_{1}$ and
$I+\mathcal{E}_{1}$ be

$Q_{1}=U_{1}\Gamma_{1},\qquad I+\mathcal{E}_{1}=U_{2}\Gamma_{2}$

respectively. Then, we have

$\Gamma_{2}Q=U_{2}^{T}(I+\mathcal{E}_{1})Q=U_{2}^{T}Q_{1}=(U_{2}^{T}U_{1})\Gamma_{1}$

Thus, by letting $\mathcal{E}_{2}=\Gamma_{2}-I$, we have

$(I+\mathcal{E}_{2})Q=(U_{2}^{T}U_{1})\Gamma_{1}\sim
U_{1}\Gamma_{1}=(I+\mathcal{E}_{1})Q$ (32)

Since $\Gamma_{2}$ is symmetric, $\mathcal{E}_{2}=\Gamma_{2}-I$ is symmetric
as well. Thus, (32) suggests that it is enough to only consider a symmetric
$\mathcal{E}$ on group GL$(n,\mathbb{R})/R_{\rm polar}$. This is not a
surprise as a symmetric $\mathcal{E}$ just have dimension $0.5n(n+1)$.

The last step is to check the lower bound of the second order terms of
$dc(P)$. Now, without the loss of generality, let us assume that
$\mathcal{E}=\mathcal{E}^{T}$ due to (32). Starting from (31), the expectation
of the second order terms of $dc(P;v,h)$ with choice $dQ=\mathcal{E}Q$ is
given by

$\displaystyle
E_{v}\left[h^{T}Q^{T}\mathcal{E}^{T}\mathcal{E}Qh+v^{T}Q^{-1}(\mathcal{E}^{2}+\mathcal{E}\mathcal{E}^{T}+\mathcal{E}^{2T})Q^{-T}v\right]$
$\displaystyle=$ $\displaystyle{\rm
tr}\left(\mathcal{E}^{2}QH^{2}Q^{T}+3\mathcal{E}^{2}Q^{-T}Q^{-1}\right)$
$\displaystyle\geq$ $\displaystyle{\rm
tr}\left(\left(\lambda_{\min}(H)\right)^{2}\mathcal{E}^{2}QQ^{T}+3\mathcal{E}^{2}\left(QQ^{T}\right)^{-1}\right)$
$\displaystyle\geq$ $\displaystyle{\rm
tr}\left(\mathcal{E}^{2}\left(\lambda_{0}^{2}QQ^{T}+3(QQ^{T})^{-1}\right)\right)$
$\displaystyle\geq$ $\displaystyle 2\sqrt{3}\lambda_{0}\,{\rm
tr}\left(\mathcal{E}^{2}\right)$ $\displaystyle=$ $\displaystyle
2\sqrt{3}\lambda_{0}\left({\rm vec}(\mathcal{E})\right)^{T}{\rm
vec}(\mathcal{E})$ $\displaystyle=$ $\displaystyle
2\sqrt{3}\lambda_{0}\left({\rm vech}(\mathcal{E})\right)^{T}S^{T}S\,{\rm
vech}(\mathcal{E})$ $\displaystyle=$ $\displaystyle\left({\rm
vech}(\mathcal{E})\right)^{T}\left(2\sqrt{3}\lambda_{0}\,S^{T}S\right){\rm
vech}(\mathcal{E})$

where we have used expectations $E_{v}[hh^{T}]=H^{2}$ and $E_{v}[vv^{T}]=I$
from the first to the second line, $aA+bA^{-1}\succ 2\sqrt{ab}I$ for $A\succ
0$ to arrive the fifth line from the fourth one, and definition ${\rm
vec}(A)=S{\rm vech}(A)$ for a symmetric matrix $A$ to arrive the last line.
Thus the Hessian of $c(P)$ in the new coordinate is lower bounded by
$2\sqrt{3}\lambda_{0}\,S^{T}S\succ 2\sqrt{3}\lambda_{0}I$ since $S^{T}S\succ
I$. $\square$

Unfortunately, Hessian fitting on the Lie group does not always make life
easier. For example, $c(P)$ is generally not convex with the coordinate system
induced by (29), choice II, even with constraint
$\mathcal{E}=\mathcal{E}^{T}$. Empirical results also have confirmed that
Hessian fitting with choice II converges significantly slower than that with
choice I. Also, $c(P)$ is generally not convex on the group of triangular
matrices with the coordinate system induced by (29), choice I. Yet, it is not
difficult to solve the triangular group Hessian fitting problem through
Proposition 5, see details in the next subsection.

###  3.2 The “standard” Lie group Hessian fitting method and its variants

As our Hessian fitting criterion is strongly convex on group
GL$(n,\mathbb{R})/R_{\rm polar}$, we can readily follow the conventional
convex optimization theory to design a proper Hessian fitting method. If we
could find an upper bound for the eigenvalues of the Hessian of $c(P)$, say
$L$, we know that gradient descent with step size $1/L$ converges linearly to
the optimal solution $H^{-1}$. From (31), choice I, we know that sum of the
second order terms of $\mathcal{E}$ is upper bounded by

$\displaystyle\left|h^{T}Q^{T}\mathcal{E}^{T}\mathcal{E}Qh+v^{T}Q^{-1}(\mathcal{E}^{2}+\mathcal{E}\mathcal{E}^{T}+\mathcal{E}^{2T})Q^{-T}v\right|$
$\displaystyle\leq$
$\displaystyle\left|h^{T}Q^{T}\mathcal{E}^{T}\mathcal{E}Qh\right|+\left|v^{T}Q^{-1}\mathcal{E}^{2}Q^{-T}v\right|+\left|v^{T}Q^{-1}\mathcal{E}\mathcal{E}^{T}Q^{-T}v\right|+\left|v^{T}Q^{-1}\mathcal{E}^{2T}Q^{-T}v\right|$
$\displaystyle\leq$
$\displaystyle\left(\|Qh\|^{2}+3\|Q^{-T}v\|^{2}\right)\|\mathcal{E}\|_{2}^{2}$
(33)

Thus, from (33), we see that an extremely rough estimate of $L$ could be

$0.5L_{t}\approx\max_{1\leq i\leq
t}\left(\|Q_{i}h_{i}\|^{2}+3\|Q_{i}^{-T}v_{i}\|^{2}\right)$ (34)

At the same time, from the linear terms in (31), choice I, we see that the
gradient on group GL$(n,\mathbb{R})$ or GL$(n,\mathbb{R})/R_{\rm polar}$ is

$0.5\nabla_{Q}=Qhh^{T}Q^{T}-Q^{-T}vv^{T}Q^{-1}$ (35)

Combining (33) and (35) gives the following “standard” Lie group SGD Hessian
fitting method,

$Q_{t+1}=Q_{t}-\mu\frac{Q_{t}h_{t}h_{t}^{T}Q_{t}^{T}-Q_{t}^{-T}v_{t}v_{t}^{T}Q_{t}^{-1}}{\|Q_{t}h_{t}\|^{2}+\|Q_{t}^{-T}v_{t}\|^{2}}Q_{t}$
(36)

where $0<\mu\leq 1$ is a normalized step size, and $(v_{t},h_{t})$ is the pair
of vector and its associated Hessian-vector product at step $t$. Note that in
(36), we do not use the more conservative $L$ estimate given in (34).
Numerical results suggest that (36) is stable enough.

####  3.2.1 Fitting on the group of triangular matrices

We can simply take the triangular part of (35) to get the gradient for
updating $Q$ on the group of triangular matrices as having done in [Li18].
But, in general $c(P)$ is not convex on this group. Still, we can factorize
out the rotation part of the updated $Q$ from (36) to have the following
update formula for $Q$ on the group of upper triangular matrix,

$Q_{t+1}=\left[I-\mu\frac{Q_{t}h_{t}h_{t}^{T}Q_{t}^{T}-Q_{t}^{-T}v_{t}v_{t}^{T}Q_{t}^{-1}}{\|Q_{t}h_{t}\|^{2}+\|Q_{t}^{-T}v_{t}\|^{2}}\right]_{R}Q_{t}$
(37)

where $[A]_{R}$ keeps the triangular part of the QR decomposition of a square
matrix $A$. In this way, $Q_{t}$ is still updated on group
GL$(n,\mathbb{R})/R_{\rm polar}$, but keeps to be a triangular matrix as long
as $Q_{0}$ is triangular. Noting that only two rank-1 QR decomposition updates
[GL96] are required to obtain the QR decomposition of the matrix inside
$[\cdot]_{R}$ in (37).

Furthermore, when $\mu$ is small enough, we can cheaply approximate the QR
decomposition needed in (37). Let $I+\delta=(I+\delta_{1})(I+\delta_{2})$,
where $\|\delta\|_{2}<1$ is a small symmetric matrix, $I+\delta_{1}$ is an
orthogonal matrix such that $(I+\delta_{1})^{T}(I+\delta_{1})=I$, and
$\delta_{2}$ is an upper triangular matrix. After ignoring the second order
terms $\delta_{1}^{T}\delta_{1}$ and $\delta_{1}\delta_{2}$, we have
$I+\delta\approx I+\delta_{1}+\delta_{2}$ and
$\delta_{1}+\delta_{1}^{T}\approx 0$. Thus, we have $\delta_{1}\approx{\rm
tril}(\delta,-1)-{\rm triu}(\delta,1)$ and $\delta_{2}={\rm triu}(\delta)+{\rm
triu}(\delta,1)$ as the approximate solutions, where ${\rm tril}(A,-1)$ and
${\rm triu}(A,1)$ exclude the diagonal elements of ${\rm tril}(A)$ and ${\rm
triu}(A)$, respectively. In this way, (37) resembles the SGD method on the
group of triangular matrices proposed in [Li18].

####  3.2.2 Step size annealing, numerical stability and Newton’s methods

Clearly, the $Q_{t}$ given in (36) and (37) will converge to $H^{-1}$ with a
noise-free Hessian-vector product model, i.e., $\epsilon=0$ in (3), due to the
strong convexity of $c(P)$. However, with a noisy Hessian-vector product
model, i.e., $\epsilon\neq 0$, the gradient given in (35) cannot settle down
to zero. Hence, we may need to anneal down $\mu$ to a small enough positive
number to encourage convergence.

From (35), we observe that we need to solve the linear system $Q^{T}x=v$ for
gradient calculation. From (36), we see that $Q_{t+1}-Q_{t}$ is a matrix with
rank two, and thus we can recursively update the inverse of $Q$ with the
Woodbury matrix identity as in BFGS. However, this may not be a reliable
procedure with single or half precision arithmetic. With (37), $Q$ is
triangular, and we can solve this linear system with backward substitution.
This procedure is empirically shown to be numerically stable even with half
precision arithmetic.

Lastly, we want to mention that the Newton’s methods on the Lie groups can be
derived from (31) by solving for the optimal $\mathcal{E}$ that reduces $c(P)$
the most. But, they turn out to be way more complicated than the form in (14),
and seem to be of little practical value. Nevertheless, the Hessian for the
case considered in Proposition 5 still has the following neat form,

$S^{T}\left[I\otimes(QH^{2}Q^{T}+3Q^{-T}Q^{-1})\right]S$

###  3.3 Sparse Lie groups for large scale Hessian fitting

The methods in (36) and (37) are only good for small to middle scale problems
with roughly $n\leq 10^{4}$ parameters to optimize. For problems with millions
to billions parameters, we need significantly sparse Hessian estimators. Table
3 lists a few examples of groups that are previously shown to work well for
large scale stochastic machine learning problems, e.g., the Kronecker product
preconditioner that can be decomposed into the direct sum and products of
smaller ones [Li19], the low-rank approximation (LRA) preconditioner that is
good for general purpose optimizations [Li22]. The theories developed here may
help to further improve the designs of these preconditioners in the future
work.

Table 3: Lie group examples for Hessian fitting Group | Matrix form | Storage for $\theta\in\mathbb{R}^{m\times m}$ | Notes
---|---|---|---
GL$(n,\mathbb{R})$ | $Q\in\mathbb{R}^{n\times n}$ | $m^{4}$ | invertible
Triangular matrices | Triangular matrices | $0.5m^{2}(m^{2}+1)$ | positive diagonals
Diagonal matrices | $Q={\rm diag}(q_{1},\ldots,q_{n})$ | $m^{2}$ | positive diagonals
Decomposable group | $Q=\oplus_{i}\left(\otimes_{j}Q_{i,j}\right)$ | $\leq 2m^{2}$ | $Q_{i,j}\in{\rm GL}(n_{i,j},\mathbb{R})$
Low-rank approx. | $Q=UV^{T}+{\rm diag}(d)$ | $(2r+1)m^{2}$ | $U,V\in R^{n\times r},r\ll n$

###  3.4 Numerical results for empirical performance studies

Fig. 2 and 3 show a few typical Hessian fitting convergence curves from
different methods with different settings. The BFGS Hessian update formula
diverges easily with a noisy Hessian-vector product model. The closed-form
solution (17) needs either double precision or a carefully tuned damping
schedule to stabilize it. Fitting on group GL$(n,\mathbb{R})/R_{\rm polar}$ is
well understood through our above analysis, and we know that it is bounded to
be stable as long as $\lambda_{\min}(H)>0$.

Figure 2: Typical convergence curves when $H$ is a $50\times 50$ matrix with
$H_{i,i}=1$ and $H_{i,j}=0.5$ for $|i-j|=1$. We have
$\lambda_{\min}(H)\rightarrow 0$ and $\lambda_{\max}(H)\rightarrow 2$ when
$n\rightarrow\infty$ for this type of matrix. For Hessian fitting on Lie
group, one candidate is (36) with a constant step size $1$, and another one is
from [Li18] which fits $Q$ on the group of triangular matrices, normalizes the
gradient in (35) with $\|\nabla_{Q}\|_{2}$, and anneals $\mu$ exponentially
from $0.1$ to $10^{-7}$. All methods are tested with single precision. On the
left panel, we test a noiseless model of (3). The closed-form solution (17)
diverges, while others converges to $H^{-1}$ within machine precision. On the
right panel, we test a noisy model of (3) with
$E[\epsilon\epsilon]=\sigma_{\epsilon}^{2}I$ and
$H^{\prime}=(H^{2}+\sigma_{\epsilon}^{2}I)^{0.5}$. Both the closed-form
solution (17) and BFGS diverge, while the Lie group Hessian fittings still
work as expected.  Figure 3: Typical convergence curves when $H$ is a time
varying process defined by $H_{t+1}=0.999H_{t}+0.001uu^{T}$, where
$u_{i}\sim\mathcal{U}(0,1)$ for $1\leq i\leq n=50$ and $H_{0}$ is a matrix
with all elements being $1/4$. For Hessian fitting on Lie group, one candidate
is (36) with $\mu=1$, and another one is the method proposed in [Li18], which
uses $\|\nabla_{Q}\|_{2}$ to normalize the gradient. Expectation in the
closed-form solution (17) is replaced by moving average with factor $0.999$.

##  4 More experimental results

The focus here is on the performance analysis for a wide range of stochastic
Hessian fitting methods. One may refer to [PL24] for the comparison results on
a large set of stochastic machine learning benchmarks. The empirical results
reported there suggest that it is generally beneficial to incorporate the
curvature information in the training for faster convergences and better
generalizations.

The results in Fig. 2 and 3 also imply that the Lie group Hessian fitting
method could be useful to the traditional convex optimizations as well. Fig. 4
reports a set of results comparing BFGS and the preconditioned SGD (PSGD)
method equipped with different Lie group preconditioners listed in Table 3 on
a tensor rank decomposition benchmark. Indeed, PSGD takes a clear advantage in
running time as it does not rely on accurate line search.

Scripts to generate all the figures in the paper are put at here.

Figure 4: A set of typical convergence curves on a tensor rank decomposition
benchmark defined by
$\min_{X,Y,Z}\sum_{i=1,j=1,k=1}^{i=20,j=50,k=100}\|T_{i,j,k}-\sum_{r=1}^{10}X_{r,i}Y_{r,j}Z_{r,k}\|^{2}$.
We have compared gradient descent (GD), LM-BFGS, and preconditioned SGD (PSGD)
[Li18] with variant Lie groups for $Q$ introduced in [Li19] and [Li22]. Left
panel shows fitting loss versus number of iterations, and the right one shows
losses versus wall time on a RTX 2080 GPU. Note that what not shown here are
that GD has a good chance to solve this problem successfully, while LM-BFGS
also has a good chance to fail. LM-BFGS per iteration is consistently slower
than other methods due to its need of line search.

## References

  * [Ama98] S. Amari.  Natural gradient works efficiently in learning.  Neural Computation, 10(2):251–276, 1998.
  * [BV04] S. Boyd and L. Vandenberghe.  Convex Optimization.  Cambridge University Press, 2004.
  * [DHS11] J. Duchi, E. Hazan, and Y. Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  Journal of Machine Learning Research, 12(7), 2011.
  * [GL96] Gene H. Golub and Charles Van Loan.  Matrix Computations.  Johns Hopkins University Press, 1996.
  * [Li18] X. L. Li.  Preconditioned stochastic gradient descent.  IEEE Trans. Neural Networks and Learning Systems, 29(5):1454–1466, 2018.
  * [Li19] X. L. Li.  Preconditioner on matrix lie group for SGD.  In 7th International Conference on Learning Representations, ICLR 2019, New Orleans, LA, USA, May 6-9, 2019. OpenReview.net, 2019.
  * [Li22] X. Li.  Black box lie group preconditioners for SGD.  In HOOML, 2022.
  * [Mur23] Kevin P. Murphy.  Probabilistic Machine Learning: Advanced Topics.  MIT Press, 2023.
  * [PL24] O. Pooladzandi and X. L. Li.  Curvature-informed SGD via general purpose lie group preconditioners, 2024.