# A homogeneous Rayleigh quotient with applications in gradient methods

Giulia Ferrandi  Michiel E. Hochstenbach

###### Abstract

Given an approximate eigenvector, its (standard) Rayleigh quotient and
harmonic Rayleigh quotient are two well-known approximations of the
corresponding eigenvalue. We propose a new type of Rayleigh quotient, the
homogeneous Rayleigh quotient, and analyze its sensitivity with respect to
perturbations in the eigenvector. Furthermore, we study the inverse of this
homogeneous Rayleigh quotient as stepsize for the gradient method for
unconstrained optimization. The notion and basic properties are also extended
to the generalized eigenvalue problem.

###### keywords:

Homogeneous Rayleigh quotient , secant condition , eigenvalue problem ,
projective coordinates , unconstrained optimization , generalized eigenvalue
problem

###### MSC:

65F15 , 65F10 , 65K05 , 90C20 , 90C30 , 65F50

††journal: -

[inst1]organization=Department of Mathematics and Computer Science, TU
Eindhoven,addressline=
PO Box 513, city=Eindhoven, postcode=5600MB, country=The Netherlands

We define a new homogeneous Rayleigh quotient as the minimizer of a residual
for eigenvalue problems.

A nonlinear Galerkin condition for this homogeneous Rayleigh quotient is
derived.

Asymptotic bounds on the relative error to an eigenvalue are obtained.

The quotient is compared with standard and harmonic Rayleigh quotients.

We study the inverse of this quotient as a stepsize for gradient methods.

##  1 Introduction

Let $A$ be an $n\times n$ symmetric matrix with eigenvalues
$\lambda_{1}\leq\cdots\leq\lambda_{n}$. We are interested in the eigenproblem
$A\mathbf{x}=\lambda\mathbf{x}$. Let $\mathbf{u}\in\mathbb{R}^{n}$ be an
approximate eigenvector to $\mathbf{x}$ with unit 2-norm. Traditionally, the
Rayleigh quotient

$\theta=\frac{\mathbf{u}^{T}\\!A\mathbf{u}}{\mathbf{u}^{T}\mathbf{u}}\qquad(=\mathbf{u}^{T}\\!A\mathbf{u})$
(1)

is the standard approach to determine an eigenvalue approximation
corresponding to $\mathbf{u}$, which yields the corresponding eigenvalue when
$\mathbf{u}$ is an exact eigenvector. This $\theta$ satisfies two (related)
optimality properties, a Galerkin condition (cf., e.g., [1, pp. 13–14])

$A\mathbf{u}-\theta\mathbf{u}\perp\mathbf{u},$ (2)

and a minimum residual condition

$\theta=\operatorname*{argmin}_{\gamma}\,\|A\mathbf{u}-\gamma\,\mathbf{u}\|.$
(3)

A second well-known quantity is the harmonic Rayleigh quotient, which is, for
instance, sometimes considered when one is interested in an eigenvalue near a
given target $\tau\in\mathbb{R}$. This is defined as (cf., e.g., [2, p. 294])

$\widetilde{\theta}_{\tau}=\frac{\mathbf{u}^{T}(A-\tau
I)\,A\mathbf{u}}{\mathbf{u}^{T}(A-\tau I)\,\mathbf{u}},$ (4)

provided that the denominator is nonzero. This quotient is also equal to the
eigenvalue when $\mathbf{u}$ is an eigenvector. It satisfies two optimality
conditions: a Galerkin condition
$(A-\widetilde{\theta}_{\tau}I)\,\mathbf{u}\perp(A-\tau I)\,\mathbf{u}$, and a
minimum residual condition (or secant condition)
$\widetilde{\theta}_{\tau}=\operatorname*{argmin}_{\gamma}\|(\gamma-\tau)^{-1}(A-\tau
I)\,\mathbf{u}-\mathbf{u}\|$ (cf. [3, Prop. 5]).

In this paper, we mainly focus on the case $\tau=0$; for nonzero
$\mathbf{u}^{T}\\!A\mathbf{u}$, the harmonic Rayleigh quotient is given by
(cf. [4])

$\widetilde{\theta}=\frac{\mathbf{u}^{T}\\!A^{2}\,\mathbf{u}}{\mathbf{u}^{T}\\!A\mathbf{u}}.$
(5)

For $\tau=0$, the Galerkin condition becomes

$A\mathbf{u}-\widetilde{\theta}\mathbf{u}\perp A\mathbf{u},$ (6)

and the minimum residual condition is

$\widetilde{\theta}=\operatorname*{argmin}_{\gamma}\,\|\gamma^{-1}A\mathbf{u}-\mathbf{u}\|.$
(7)

Condition (7) is easily verified by setting the derivative of
$\|\gamma^{-1}A\mathbf{u}-\mathbf{u}\|$ with respect to $\gamma$ to zero.
Property (7) is perhaps less widely known. Both Properties (3) and (7) have
been exploited as secant conditions to determine a stepsize in gradient-type
optimization methods; for instance for the well-known Barzilai and Borwein
steplengths [5]. We return to this topic in Section 5; see also Vömel [6].

As an alternative to (3) and (7), we propose and study a new homogeneous type
of Rayleigh quotient in this paper. We work in the projective space
$\mathbb{P}=\mathbb{P}_{1}(\mathbb{R})$, also called the _projective line_.
This is a quotient space of $\mathbb{R}^{2}\setminus\\{\bf 0\\}$, with the
following equivalence relation: two points $[\alpha_{1},\alpha_{2}]^{T}$,
$[\widetilde{\alpha}_{1},\widetilde{\alpha}_{2}]^{T}\in\mathbb{R}^{2}\setminus\\{\bf
0\\}$ are equivalent if there exists a nonzero $\zeta$ such that
$[\alpha_{1},\alpha_{2}]^{T}=\zeta\,[\widetilde{\alpha}_{1},\widetilde{\alpha}_{2}]^{T}$.
The coordinates of the representative of a class $(\alpha_{1},\alpha_{2})$ are
called _homogeneous coordinates_. The pair $(1,0)$ corresponds to the point at
infinity. All other equivalence classes $(\alpha_{2}\neq 0)$ can be identified
with a real number $\alpha=\alpha_{1}/\alpha_{2}$.

Homogeneous techniques for eigenvalue problems have already been developed by
various authors. Stewart and Sun [7, p. 283] exploit homogeneous coordinates
for the generalized eigenvalue problem $A\mathbf{x}=\lambda B\,\mathbf{x}$,
for $B\in\mathbb{R}^{n\times n}$, since this is “especially convenient for
treating infinite eigenvalues” [8]. While the standard eigenvalue problem
$A\mathbf{x}=\lambda\mathbf{x}$ does not involve infinite eigenvalues, the
harmonic Rayleigh quotient (5) may be infinite. Dedieu and Tisseur [9, 10]
exploit homogeneous techniques to study perturbation theory for generalized
and polynomial eigenproblems. Inspired by this, a homogeneous Jacobi–Davidson
method for subspace expansion has been proposed in [11]; this may be viewed as
an inexact Newton method.

The contributions and outline of this paper are the following. We define a new
_homogeneous Rayleigh quotient_ as the minimizer of the residual quantity

$\min_{(\alpha_{1},\,\alpha_{2})\,\in\,\mathbb{P}}\
\frac{\|\alpha_{1}\,\mathbf{u}-\alpha_{2}\,A\mathbf{u}\|}{{\color[rgb]{0,0,0}\sqrt{\alpha_{1}^{2}+\alpha_{2}^{2}}}}.$

In Section 2, we show that there is a closed-form solution to express this
homogeneous Rayleigh quotient, and that it can also be obtained from a certain
nonlinear Galerkin condition. To quantify the quality of the homogeneous
Rayleigh quotient as an approximate eigenvalue, a bound for the (chordal)
distance between the homogeneous Rayleigh quotient and the spectrum of $A$ is
derived.

We compare the homogeneous Rayleigh quotient with the standard and harmonic
Rayleigh quotients. Given the relation between the homogeneous residual
quantity and the two residual quantities (3) and (7) (see also Section 2.3),
the homogeneous Rayleigh quotient may be seen as the mediator between the two
other Rayleigh quotients. In Section 3, we highlight this fact in the
theoretical and experimental study on the relative errors of the three
quotients. Asymptotic bounds on the accuracy of the three quotients are also
derived and compared. Interestingly, the homogeneous Rayleigh quotient may be
more accurate than its two competitors in some situations.

In Section 4, we extend the definition of the homogeneous Rayleigh quotient to
the generalized eigenvalue problem $(A,B)$, with $A$ symmetric and $B$
symmetric positive definite (SPD). Finally, in Section 5, we propose a new
stepsize for gradient methods based on the homogeneous Rayleigh quotient.
Conclusions are drawn in Section 6.

We would like to stress that the stepsize (31) presented in Section 5 has been
recently proposed first by Li, Zhang, and Xia [12] (with preliminary version
[13] in 2021). The results in our Section 5 have been developed independently
and the derivation is different: whereas the approach by Li, Zhang, and Xia is
via a total least squares technique, ours is based on exploiting homogeneous
coordinates. We will discuss some more details in Section 5.

##  2 A homogeneous Rayleigh quotient

For a given approximate eigenvector $\mathbf{u}$, we consider the minimization
of the _homogeneous_ residual (as mediator of the two residual quantities in
(3) and (7)):

$\min_{(\alpha_{1},\,\alpha_{2})\,\in\,\mathbb{P}}\
\frac{\|\alpha_{1}\,\mathbf{u}-\alpha_{2}\,A\mathbf{u}\|}{{\color[rgb]{0,0,0}\sqrt{\alpha_{1}^{2}+\alpha_{2}^{2}}}}=\min_{\alpha_{1}^{2}+\alpha_{2}^{2}=1}\
\|[\mathbf{u}\ \
-\\!\\!A\mathbf{u}]{\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{c}\alpha_{1}\\\
\alpha_{2}\end{array}\\!\\!\right]$}}\|.$ (8)

Note that the expression on the left-hand side is well defined on
$\mathbb{P}$. The equality between the two optimization problems can be seen
by taking
$(\alpha_{1},\alpha_{2})\leftarrow(\alpha_{1}^{2}+\alpha_{2}^{2})^{-1/2}\
(\alpha_{1},\alpha_{2})$ as representative of the equivalence class. Then we
minimize the objective function on
$\\{[\alpha_{1},\alpha_{2}]^{T}\in\mathbb{R}^{2}\mid\alpha_{1}^{2}+\alpha_{2}^{2}=1\\}$.
For nonzero $\alpha_{2}$ we call the ratio between the coordinates of the
solution $\alpha=\alpha_{1}/\alpha_{2}$ the _homogeneous Rayleigh quotient_.

###  2.1 Key properties of the homogeneous Rayleigh quotient

We discuss some properties of the solution to (8). First, let us introduce the
$n\times 2$ matrix $C=[\mathbf{u}\ \ \,-\\!\\!A\mathbf{u}]$ and the associated
$2\times 2$ matrix

$C^{T}C=\left[\begin{array}[]{cc}\phantom{-}\mathbf{u}^{T}\mathbf{u}&-\mathbf{u}^{T}\\!A\mathbf{u}\\\
-\mathbf{u}^{T}\\!A\mathbf{u}&\phantom{-}\mathbf{u}^{T}\\!A^{2}\mathbf{u}\end{array}\right]=\mathbf{u}^{T}\\!A\mathbf{u}\cdot\left[\begin{array}[]{cc}\phantom{-}\theta^{-1}&-1\\\
-1&\phantom{-}\widetilde{\theta}\end{array}\right],$ (9)

where the second equality holds when $\mathbf{u}^{T}\\!A\mathbf{u}$ is
nonzero. Although we usually impose $\|\mathbf{u}\|=1$, we will not exploit
this simplification in view of more general use in Section 5.

In the next proposition, assuming that $\mathbf{u}^{T}\\!A\mathbf{u}\neq 0$,
we can provide an explicit formula for the homogeneous Rayleigh quotient
$\alpha$. In addition, we show that $\alpha$ is an eigenvalue of $A$ when
$\mathbf{u}$ is an eigenvector, and that it is located between the Rayleigh
quotient and the harmonic Rayleigh quotient, depending on the sign of
$\mathbf{u}^{T}\\!A\mathbf{u}$.

######  Proposition 1.

Suppose that $\mathbf{u}^{T}\\!A\mathbf{u}\neq 0$ and denote the smallest
eigenvalue of $C^{T}C$ by $\mu$. Then the following properties hold.

  * (i)

The value $\mu$ is a simple eigenvalue of $C^{T}C$, or, equivalently,
$\sqrt{\mu}$ is a simple singular value of $C$; its corresponding eigenvector
$[\alpha_{1},\ \alpha_{2}]^{T}$ is the unique minimizer of (8).

  * (ii)

The value of $\mu$ is

$\displaystyle\mu=\tfrac{1}{2}\,\big{[}\mathbf{u}^{T}\mathbf{u}+\mathbf{u}^{T}\\!A^{2}\,\mathbf{u}-\sqrt{\smash[b]{(\mathbf{u}^{T}\mathbf{u}-\mathbf{u}^{T}\\!A^{2}\,\mathbf{u})^{2}+4\,(\mathbf{u}^{T}\\!A\mathbf{u})^{2}}}\
\big{]}.$ (10)

  * (iii)

For the homogeneous Rayleigh quotient $\alpha$ we have

$\displaystyle\alpha$
$\displaystyle=\frac{\alpha_{1}}{\alpha_{2}}=\frac{\mathbf{u}^{T}\\!A\mathbf{u}}{\mathbf{u}^{T}\mathbf{u}-\mu}=\frac{\mathbf{u}^{T}\\!A^{2}\,\mathbf{u}-\mu}{\mathbf{u}^{T}\\!A\mathbf{u}}$
$\displaystyle=\frac{1}{2\ \mathbf{u}^{T}\\!A\mathbf{u}}\
\big{[}\,\mathbf{u}^{T}\\!A^{2}\,\mathbf{u}-\mathbf{u}^{T}\mathbf{u}+\sqrt{\smash[b]{(\mathbf{u}^{T}\\!A^{2}\,\mathbf{u}-\mathbf{u}^{T}\mathbf{u})^{2}+4\,(\mathbf{u}^{T}\\!A\mathbf{u})^{2}}}\,\big{]}.$
(11)

In terms of the standard and harmonic Rayleigh quotient,

$\alpha=\tfrac{1}{2}\,\big{[}\,\widetilde{\theta}-\theta^{-1}+\sqrt{\smash[b]{(\widetilde{\theta}-\theta^{-1})^{2}+4}}\,\big{]}=\big{(}\tfrac{1}{2}\,\big{[}\,\theta^{-1}-\widetilde{\theta}+\sqrt{\smash[b]{(\theta^{-1}-\widetilde{\theta})^{2}+4}}\,\big{]}\big{)}^{-1}.$

  * (iv)

$\mu=0$ if and only if $\mathbf{u}$ is an eigenvector of $A$. If $\mathbf{u}$
is an eigenvector, then the corresponding homogeneous Rayleigh quotient
$\alpha$ is the corresponding eigenvalue.

  * (v)

$0\leq\mu<\min(\mathbf{u}^{T}\mathbf{u},\ \mathbf{u}^{T}\\!A^{2}\mathbf{u})$.

  * (vi)

We have the following bounds:

$\begin{cases}\theta\leq\alpha\leq\widetilde{\theta}&\text{if}\quad\mathbf{u}^{T}\\!A\mathbf{u}>0,\\\\[1.42262pt]
\widetilde{\theta}\leq\alpha\leq\theta&\text{if}\quad\mathbf{u}^{T}\\!A\mathbf{u}<0.\end{cases}$

###### Proof.

It is well known that the minimizer of (8) is the eigenvector corresponding to
the smallest eigenvalue of the matrix $C^{T}C$; or, equivalently, to the right
singular vector of $C$ corresponding to the smallest singular value. We start
with the eigenvalue problem

$C^{T}C\,{\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{c}\alpha_{1}\\\
\alpha_{2}\end{array}\\!\\!\right]$}}=\mu{\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{c}\alpha_{1}\\\
\alpha_{2}\end{array}\\!\\!\right]$}}.$

The eigenvalue $\mu$ of (9) is simple because $\mathbf{u}^{T}\\!A\mathbf{u}$
is assumed nonzero, so that the discriminant of the quadratic characteristic
polynomial is positive. In particular, this implies that the smallest
eigenvector $[\alpha_{1},\ \alpha_{2}]^{T}$ is well defined. The result in
part (ii) is an explicit expression for the eigenvalue $\mu$. Its
corresponding eigenvector $[\alpha_{1},\ \alpha_{2}]^{T}$ is used in item
(iii) to give an expression of the ratio $\alpha=\alpha_{1}/\alpha_{2}$, which
gives the homogeneous Rayleigh quotient, that can also be written in terms of
the standard and the harmonic Rayleigh quotients.

If $\mu=0$, this means that $C$ has rank 1. This happens precisely when
$\mathbf{u}$ is an eigenvector. Moreover, if $\mathbf{u}$ is an eigenvector
corresponding to some eigenvalue $\lambda$, equations (10) and ((iii)) give
$\mu=0$ and $\alpha=\lambda$ respectively. This concludes the proof of part
(iv).

The nonnegativity of $\mu$ directly follows from the fact that $C^{T}C$ is
positive semidefinite. Since $\mathbf{u}^{T}\\!A\mathbf{u}\neq 0$ we have
$\sqrt{\smash[b]{(\mathbf{u}^{T}\mathbf{u}-\mathbf{u}^{T}\\!A^{2}\,\mathbf{u})^{2}+4\,(\mathbf{u}^{T}\\!A\mathbf{u})^{2}}}>|\mathbf{u}^{T}\mathbf{u}-\mathbf{u}^{T}\\!A^{2}\,\mathbf{u}|$,
and therefore $\mathbf{u}^{T}\mathbf{u}-\mu,\
\mathbf{u}^{T}\\!A^{2}\mathbf{u}-\mu>0$. It follows that
$0\leq\mu<\min(\mathbf{u}^{T}\mathbf{u},\,\mathbf{u}^{T}\\!A^{2}\mathbf{u})$.
In particular, this shows that the sign of $\alpha$ only depends on the
quantity $\mathbf{u}^{T}\\!A\mathbf{u}$. Given this fact, property (vi) is a
direct consequence of part (v). ∎

In the exceptional case of $\mathbf{u}^{T}\\!A\mathbf{u}=0$, the matrix
$C^{T}C$ is diagonal, with zero Rayleigh quotient and infinite harmonic
Rayleigh quotient. The homogeneous Rayleigh quotient
$\alpha=(\alpha_{1},\alpha_{2})$ is either zero (that is, $(0,1)$), or the
point at infinity $(1,0)$. The situation $\mathbf{u}^{T}\\!A\mathbf{u}=0$
cannot happen if $A$ is definite, such as for quadratic optimization problems
with an SPD Hessian; see also Section 5. Moreover,
$\mathbf{u}^{T}\\!A\mathbf{u}=0$ may not be a very common situation in the
context of eigenproblems; this usually means that we are approximating a zero
eigenvalue or that $\mathbf{u}$ is a poor approximate eigenvector.

Rayleigh quotients have some natural invariance properties. It is easy to
check that the standard Rayleigh quotient is invariant under nonzero scalings
of $\mathbf{u}$; it is invariant under nonzero scalings of $A$, and it shifts
naturally with shifts of $A$. The harmonic Rayleigh quotient (with target
$\tau$) only satisfies the first property, but modified properties still hold,
as follows.

######  Proposition 2.

For any $\tau\in\mathbb{R}$, the following properties hold:

  * (i)

$\widetilde{\theta}_{\tau}$ is invariant under scaling of $\mathbf{u}$:
$\widetilde{\theta}_{\tau}(\zeta\,\mathbf{u})=\widetilde{\theta}_{\tau}(\mathbf{u})$
for $\zeta\neq 0$.

  * (ii)

For $\widetilde{\theta}_{\tau}$ as function of scaling of $A$ we have
$\widetilde{\theta}_{\tau}(\zeta\,A)=\zeta\,\widetilde{\theta}_{\zeta^{-1}\,\tau}(A)$
for $\zeta\neq 0$.

  * (iii)

For $\widetilde{\theta}_{\tau}$ as function of shifts of $A$ we have
$\widetilde{\theta}_{\tau}(A-\eta I)=\widetilde{\theta}_{\tau+\eta}(A)-\eta$
for $\zeta\in\mathbb{R}$.

###### Proof.

This is easy to check from (4). ∎

As for the homogeneous Rayleigh quotient, it is easy to see (for instance,
from ((iii))) that the solution to the minimization problem (8) is invariant
under nonzero scaling of $\mathbf{u}$: as function of $\mathbf{u}$, we have
$\alpha(\zeta\,\mathbf{u})=\alpha(\mathbf{u})$ for $\zeta\neq 0$. Second,
while the standard and harmonic Rayleigh quotients are invariant under
multiplications $A\to\zeta A$, there is no easy relation between $\alpha(\zeta
A)$ and $\alpha(A)$. However, in view of Proposition 1(vi), we establish these
bounds:

$\begin{cases}\zeta\,\theta(A)\leq\alpha(\zeta
A)\leq\zeta\,\widetilde{\theta}(A)&\text{if}\quad\mathbf{u}^{T}\\!A\mathbf{u}>0,\\\\[1.42262pt]
\zeta\,\widetilde{\theta}(A)\leq\alpha(\zeta
A)\leq\zeta\,\theta(A)&\text{if}\quad\mathbf{u}^{T}\\!A\mathbf{u}<0.\end{cases}$

The same upper and lower bounds hold for $\zeta\alpha(A)$. Therefore, if the
interval $[\zeta\theta(A),\,\zeta\widetilde{\theta}(A)]$ (or
$[\widetilde{\zeta}\theta(A),\,\zeta\theta(A)]$) is small, we conclude that
$\alpha(\zeta A)\approx\zeta\alpha(A)$.

On the size of this interval, we note that when $\theta\leq\widetilde{\theta}$
(which for instance holds in the SPD case), the “relative size” of the
interval $[\theta,\,\widetilde{\theta}]$ can be expressed in terms of the
angle between $\mathbf{u}$ and $A\mathbf{u}$, since
$(\widetilde{\theta}-\theta)\,/\,\theta=\cos^{-2}(A\mathbf{u},\mathbf{u})-1=\tan^{2}(A\mathbf{u},\mathbf{u})$.
This means that if $\angle(A\mathbf{u},\mathbf{u})$ is small, the interval is
small, and all three Rayleigh quotients are close. Still, Example 1 will show
that one quotient can still be considerably more accurate than the other two.

Moreover, we have

$\frac{\alpha-\theta}{\theta}=\tfrac{1}{2}\,\big{[}\,1+\tan^{2}(A\mathbf{u},\mathbf{u})-\theta^{-2}+\sqrt{\smash[b]{(1+\tan^{2}(A\mathbf{u},\mathbf{u})-\theta^{-2})^{2}+4\,\theta^{-2}}}\,-2\big{]}.$

For $\theta^{-1}\ll\tan(A\mathbf{u},\mathbf{u})$ (that is, for large
eigenvalues), this expression is close to $\tan^{2}(A\mathbf{u},\mathbf{u})$,
which means that $\alpha\approx\widetilde{\theta}$. In contrast, if
$\theta^{-1}\gg\max(\tan(A\mathbf{u},\mathbf{u}),\,1)$ (i.e., for small
eigenvalues), one may check that
$\frac{\alpha-\theta}{\theta}=\mathcal{O}(\theta^{2})\approx 0$, which implies
that $\alpha\approx\theta$. In conclusion, for a small
$\angle(A\mathbf{u},\mathbf{u})$, the homogeneous Rayleigh quotient is close
to the Rayleigh quotient for small eigenvalues, and close to the harmonic
Rayleigh quotient for large ones. Similar remarks can be made when
$\widetilde{\theta}<\theta$.

Simple manipulations of Proposition 1(vi) also lead to some relations between
the relative errors of the three Rayleigh quotients, when the estimated
eigenvalues are either $\lambda_{1}$ or $\lambda_{n}$. If $A$ is SPD, then
both the standard and harmonic Rayleigh quotients lie in the spectrum of $A$,
i.e., $\theta,\,\widetilde{\theta}\in[\lambda_{1},\lambda_{n}]$. Thus the
following inequalities hold:

$\theta-\lambda_{1}\leq\alpha-\lambda_{1}\leq\widetilde{\theta}-\lambda_{1},\qquad\lambda_{n}-\widetilde{\theta}\leq\lambda_{n}-\alpha\leq\lambda_{n}-\theta.$
(12)

We conclude that, in the SPD case, the Rayleigh quotient is more accurate when
we estimate the smallest eigenvalue of $A$, while the harmonic Rayleigh
quotient is more accurate for the largest eigenvalue. The homogeneous Rayleigh
quotient lies in between. This may suggest the perhaps not so well-known
result that to approximate large eigenvalues, the harmonic Rayleigh quotient
may be preferable over the standard Rayleigh quotient.

For indefinite matrices, while $\theta\in[\lambda_{1},\lambda_{n}]$, the
harmonic and the homogeneous Rayleigh quotients can lie outside the spectrum.
In particular, (12) does not hold if $A$ is not SPD. Example 1 shows that the
harmonic and homogeneous Rayleigh quotients may overestimate the largest
eigenvalues, and the sensitivity of the homogeneous Rayleigh quotient may be
smaller than the other two.

######  Example 1.

Let $A=\text{diag}(-\frac{2}{3},\ \frac{1}{3},\ 2)$, and consider the
eigenvector $\mathbf{x}=[0,\,0,\,1]^{T}$ corresponding to $\lambda_{3}=2$.
Given $\mathbf{u}=[-0.02,\ 0.01,\ 1]^{T}$ as approximation to $\mathbf{x}$,
with $\angle(\mathbf{u},\mathbf{x})\approx 2.2\cdot 10^{-2}$, we have that

$|\theta(\mathbf{u})-\lambda_{3}|\approx 1.2\cdot
10^{-3},\qquad|\widetilde{\theta}(\mathbf{u})-\lambda_{3}|\approx 3.3\cdot
10^{-4},\qquad|\alpha(\mathbf{u})-\lambda_{3}|\approx 1.6\cdot 10^{-5}.$

The quotients are $\theta(\mathbf{u})\approx 1.9988<\lambda_{3}$,
$\widetilde{\theta}(\mathbf{u})\approx 2.0003>\lambda_{3}$ and
$\alpha(\mathbf{u})\approx 2.00002>\lambda_{3}$, so the harmonic and the
homogeneous Rayleigh quotients overestimate $\lambda_{3}$. Nevertheless, the
homogeneous Rayleigh quotient is an order of magnitude more accurate than the
harmonic Rayleigh quotient, and two orders more accurate than the standard
Rayleigh quotient. This example highlights that the homogeneous Rayleigh
quotient might be more accurate, especially for exterior eigenvalues of
indefinite matrices. We will show another example in Figure 1 in Section 3.

When the wanted eigenvalue is in the interior of the spectrum, or if $A$ is
indefinite, it is more difficult to derive general conclusions about the
ordering of the relative errors. We will discuss this in more detail in
Section 3.

###  2.2 A Galerkin condition for the homogeneous Rayleigh quotient

In view of (2) and (6), the question arises whether there exists a Galerkin
(orthogonality) condition based on the span of $\mathbf{u}$ and $A\mathbf{u}$
for the homogeneous approach. Let us first express the homogeneous Rayleigh
quotient as a solution to a quadratic equation. As a direct consequence, this
property connects the homogeneous Rayleigh quotient and the harmonic Rayleigh
quotient with target.

######  Proposition 3.

Let $\mathbf{u}^{T}\\!A\mathbf{u}\neq 0$. The homogeneous Rayleigh quotient
((iii)) solves

$(\mathbf{u}^{T}\\!A\mathbf{u})\,\alpha^{2}+(\mathbf{u}^{T}\mathbf{u}-\mathbf{u}^{T}\\!A^{2}\mathbf{u})\,\alpha-(\mathbf{u}^{T}\\!A\mathbf{u})=0$
(13)

and satisfies $\alpha\ \mathbf{u}^{T}\\!A\mathbf{u}>0$. In addition, the
homogeneous Rayleigh quotient is a harmonic Rayleigh quotient with target (4),
where $\tau=-\alpha^{-1}$.

###### Proof.

It is straightforward to check that ((iii)) is a solution to (13). The
solutions to (13) have opposite signs since their product is $-1$. Given that
$0\leq\mu<\min(\mathbf{u}^{T}\mathbf{u},\,\mathbf{u}^{T}\\!A^{2}\mathbf{u})$
from Proposition 1, the homogeneous Rayleigh quotient has the same sign as
$\mathbf{u}^{T}\\!A\mathbf{u}$, thus it corresponds to the solution to (13)
for which $\alpha\ \mathbf{u}^{T}\\!A\mathbf{u}>0$. The second part follows
immediately from the fact that imposing
$\widetilde{\theta}_{-\alpha^{-1}}=\alpha$ is equivalent to solving (13). ∎

We remark that the second solution of (13) is related to the smallest
eigenvector of $C^{T}C$, which maximizes (rather than minimizes) the
homogeneous residual quantity (8).

Another equivalent point of view on this proposition is the following.
Equation (13) can also be expressed as

$\mathbf{u}^{T}(A+\alpha^{-1}I)\,(A-\alpha I)\,\mathbf{u}=0.$ (14)

Therefore, we have the nonlinear Galerkin condition $(A-\alpha
I)\,\mathbf{u}\perp(A+\alpha^{-1}I)\,\mathbf{u}$, or, in homogeneous
coordinates,

$(\alpha_{2}A-\alpha_{1}I)\,\mathbf{u}\perp(\alpha_{1}A+\alpha_{2}I)\,\mathbf{u}.$

This again highlights that $\alpha$ can be viewed as a harmonic Rayleigh
quotient with $-\alpha^{-1}$ as target.

###  2.3 Properties of the homogeneous residual quantity

There is a connection between the two minimal residual conditions of the
standard and harmonic Rayleigh quotients ((3) and (7)) and the quadratic
equation (13). The stationarity conditions for (3) and (7), respectively, can
be stated as

$\mathbf{u}^{T}\mathbf{u}-(\mathbf{u}^{T}\\!A\mathbf{u})\,\gamma^{-1}=0,\qquad(\mathbf{u}^{T}\\!A\mathbf{u})\,\gamma-\mathbf{u}^{T}\\!A^{2}\mathbf{u}=0.$

This means that (14) is a linear combination of the two stationarity
conditions of the standard and the harmonic Rayleigh quotients. This fact,
together with the bounds derived in Proposition 1(vi), lets us interpret the
homogeneous Rayleigh quotient as a “mediator” between the standard and the
harmonic Rayleigh quotient. Experiments in Section 3.3 will also highlight
this mediating behavior.

It is an open question whether it is possible to relate the objective function
(8) to a combination of the two residuals associated with the standard and the
harmonic Rayleigh quotients. Nevertheless, it is possible to show that, when
$\mathbf{u}^{T}\\!A\mathbf{u}\neq 0$, the homogeneous residual quantity is
smaller than the other two residual quantities (see (3) and (7)). As remarked
after Proposition 1, this hypothesis excludes the points $(0,1)$ and $(1,0)$
from the admissible minimizers. Therefore, (8) is restricted to
$(\alpha_{1},\alpha_{2})\in\mathbb{P}$, with $\alpha_{1},\alpha_{2}\neq 0$.
Then we use either $(\gamma,1)$, or $(1,\gamma)$, for a nonzero $\gamma$, as
class representatives. We derive the following minima, which are equivalent to
(8) under $\mathbf{u}^{T}\\!A\mathbf{u}\neq 0$:

$\min_{\gamma\neq 0}\
\frac{1}{\sqrt{\gamma^{2}+1}}\,\|A\mathbf{u}-\gamma\mathbf{u}\|=\min_{\gamma\neq
0}\
\frac{|\gamma|}{\sqrt{\gamma^{2}+1}}\,\|\gamma^{-1}A\mathbf{u}-\mathbf{u}\|.$
(15)

This means that the homogeneous residual quantity corresponds to the standard
and harmonic residual quantities multiplied by a factor smaller than one.

From (15), we may also relate the homogeneous residual to a measure of
distance between the homogeneous Rayleigh quotient and the spectrum of $A$. As
stated in [1, Thm. 4.5.1] and [6, Thm. 3.7], for the standard Rayleigh
quotient and the harmonic Rayleigh quotient there exists certain eigenvalues
of $A$, say $\lambda$ and $\widetilde{\lambda}$, such that

$|\lambda-\theta|\leq\frac{\|(A-\theta\,I)\,\mathbf{u}\|}{\|\mathbf{u}\|},\qquad|\widetilde{\lambda}^{-1}-\widetilde{\theta}^{-1}|\leq\frac{\|(\widetilde{\theta}^{-1}I-A)\,\mathbf{u}\|}{\|A\mathbf{u}\|},$

provided that $A\mathbf{u}\neq\mathbf{0}$ in the second case. These
inequalities hold for any $\theta$ and $\widetilde{\theta}$, but the Rayleigh
quotients $\theta$ and $\widetilde{\theta}$ have the advantage of minimizing
the corresponding residual quantities. A similar result can be stated for the
homogeneous Rayleigh quotient in terms of the chordal metric (see, e.g., [7,
Ch. 6, Def. 1.20]). By applying [1, Thm. 4.5.1], we get the following result:
for any $\alpha$, there exists an eigenvalue $\lambda$ of $A$ such that

$\frac{|\lambda-\alpha|}{\sqrt{1+\lambda^{2}}\
\sqrt{1+\alpha^{2}}}\leq\frac{1}{\sqrt{1+\lambda^{2}}\
\|\mathbf{u}\|}\cdot\frac{\|(A-\alpha\,I)\,\mathbf{u}\|}{\sqrt{1+\alpha^{2}}},$
(16)

where the upper bound is minimized by the choice of the homogeneous Rayleigh
quotient for $\alpha$.

Note that, as explained in [7, Ch. 6], the chordal metric behaves
counterintuitively for large eigenvalues: while the chordal metric is small,
the relative error can still be large. For a fair comparison of the accuracy
of the three Rayleigh quotients, we consider asymptotic bounds on their
relative error when the vector $\mathbf{u}$ approximates an eigenvector in the
next section.

##  3 Sensitivity analysis

We now study the sensitivity of the various Rayleigh quotients with respect to
perturbations in the approximate eigenvector $\mathbf{u}$. Note that this
sensitivity is different from that expressed by the condition number
$\kappa(\lambda)$ of an eigenvalue, which is related to its perturbation as a
function of changes in the matrix $A$. We note that for a simple eigenvalue
$\lambda$ of a symmetric $A$, it holds that $\kappa(\lambda)=1$; this means
that the eigenvalue is perfectly conditioned (see, e.g., [1, p. 16]).

Studying the sensitivity with respect to the approximate eigenvector for the
standard and harmonic Rayleigh quotient is certainly not new (cf., e.g., [14,
15] and the references therein), but we derive a new result for the
homogeneous Rayleigh quotient, and obtain expressions that allow an easy
comparison of the three Rayleigh quotients. We will also comment on the
differences between existing results and ours.

###  3.1 Rayleigh quotient and harmonic Rayleigh quotient

As ansatz, suppose $\mathbf{u}=\mathbf{x}+\mathbf{e}$ is an approximate
eigenvector, where $\|\mathbf{x}\|=1$, $\mathbf{e}\perp\mathbf{x}$, and
$\varepsilon:=\|\mathbf{e}\|$ is small. Then the perturbation of the standard
Rayleigh quotient is

$\theta(\mathbf{u})=\frac{\mathbf{u}^{T}\\!A\mathbf{u}}{\mathbf{u}^{T}\mathbf{u}}=\frac{\lambda+\mathbf{e}^{T}\\!A\mathbf{e}}{1+\mathbf{e}^{T}\mathbf{e}}=\lambda\,(1+\mathbf{e}^{T}(\lambda^{-1}A-I)\,\mathbf{e})+\mathcal{O}(\varepsilon^{4}),$
(17)

where we use the fact that $(1+t)^{-1}=1-t+\mathcal{O}(t^{2})$ for $t\to 0$.
For the harmonic Rayleigh quotient we have

$\displaystyle\widetilde{\theta}(\mathbf{u})=\frac{\mathbf{u}^{T}\\!A^{2}\,\mathbf{u}}{\mathbf{u}^{T}\\!A\mathbf{u}}$
$\displaystyle=\frac{\lambda^{2}+\mathbf{e}^{T}\\!A^{2}\mathbf{e}}{\lambda+\mathbf{e}^{T}\\!A\mathbf{e}}=\lambda\,(1+\lambda^{-2}\,\mathbf{e}^{T}\\!A^{2}\mathbf{e})\,(1-\lambda^{-1}\,\mathbf{e}^{T}\\!A\mathbf{e})+\mathcal{O}(\varepsilon^{4})$
$\displaystyle=\lambda\,(1+\mathbf{e}^{T}\lambda^{-1}A\,(\lambda^{-1}A-I)\,\mathbf{e})+\mathcal{O}(\varepsilon^{4}).$
(18)

For the harmonic Rayleigh quotient with target $\tau$ the expression becomes

$\widetilde{\theta}_{\tau}(\mathbf{u})=\lambda\,(1+\mathbf{e}^{T}\,(\lambda-\tau)^{-1}(A-\tau
I)\,(\lambda^{-1}A-I)\,\mathbf{e})+\mathcal{O}(\varepsilon^{4}).$ (19)

In particular, we note that
$\theta(\mathbf{x})=\widetilde{\theta}(\mathbf{x})=\widetilde{\theta}_{\tau}(\mathbf{x})=\lambda$.
We also point out that we have studied similar expressions for inverse
Rayleigh quotients (as stepsizes for gradient methods) in [3].

To derive the lower and upper bounds on the sensitivity of the approximate
eigenvalues, we make use of this standard result for symmetric operators.

######  Lemma 2.

Suppose $\mathbf{u}=\mathbf{x}+\mathbf{e}$ is an approximate eigenvector
corresponding to a simple eigenvalue $\lambda$, with $\|\mathbf{x}\|=1$,
$\mathbf{e}\perp\mathbf{x}$, and $\varepsilon=\|\mathbf{e}\|$, of a symmetric
$A$. Let $p$ be a polynomial. Then

$\varepsilon^{2}\,\min_{\lambda_{i}\neq\lambda}\ |p(\lambda_{i})|\ \leq\
|\mathbf{e}^{T}p(A)\,\mathbf{e}|\ \leq\
\varepsilon^{2}\,\max_{\lambda_{i}\neq\lambda}\ |p(\lambda_{i})|.$

###### Proof.

Since $A$ is symmetric, it is diagonalizable by an orthogonal transformation,
and we may assume that $A=\text{diag}(\lambda_{1},\dots,\lambda_{n})$. With
the notation $\lambda=\lambda_{j}$, it follows that
$|\mathbf{e}^{T}\\!p(A)\mathbf{e}|=|\sum_{i\neq j}e_{i}^{2}\,p(\lambda_{i})|$,
from which the result follows easily. ∎

The above expressions imply the following lower and upper bounds on the
sensitivity of the approximate eigenvalues as function of the approximate
eigenvector. The assumption $\lambda\neq 0$ in the next propositions may be
viewed as non-restrictive: for zero eigenvalues we can consider the absolute
error $|\theta(\mathbf{u})-\lambda|$ instead.

######  Proposition 4.

Suppose $\mathbf{u}=\mathbf{x}+\mathbf{e}$ is an approximate eigenvector
corresponding to a simple eigenvalue $\lambda\neq 0$ of a symmetric $A$, with
$\|\mathbf{x}\|=1$, $\mathbf{e}\perp\mathbf{x}$, and
$\varepsilon=\|\mathbf{e}\|$. Then, up to
$\mathcal{O}(\varepsilon^{4})$-terms, for the sensitivity of the Rayleigh
quotient (as function of $\mathbf{u}$) it holds that

$\min_{\lambda_{i}\neq\lambda}\frac{|\lambda_{i}-\lambda|}{|\lambda|}\
\varepsilon^{2}\ \lesssim\ \frac{|\theta(\mathbf{u})-\lambda|}{|\lambda|}\
\lesssim\
\max_{\lambda_{i}\neq\lambda}\frac{|\lambda_{i}-\lambda|}{|\lambda|}\
\varepsilon^{2}.$ (20)

If $\tau$ is not equal to an eigenvalue, then for the harmonic Rayleigh
quotient with target $\widetilde{\theta}_{\tau}$:

$\min_{\lambda_{i}\neq\lambda}\frac{|(\lambda_{i}-\tau)(\lambda_{i}-\lambda)|}{|\lambda(\lambda-\tau)|}\,\varepsilon^{2}\
\lesssim\ \frac{|\widetilde{\theta}_{\tau}(\mathbf{u})-\lambda|}{|\lambda|}\
\lesssim\
\max_{\lambda_{i}\neq\lambda}\frac{|(\lambda_{i}-\tau)(\lambda_{i}-\lambda)|}{|\lambda(\lambda-\tau)|}\,\varepsilon^{2}.$
(21)

In particular, for the sensitivity of the harmonic Rayleigh quotient (i.e.,
$\tau=0$) we have

$\min_{\lambda_{i}\neq\lambda}\frac{|\lambda_{i}(\lambda_{i}-\lambda)|}{\lambda^{2}}\,\varepsilon^{2}\
\lesssim\ \frac{|\widetilde{\theta}(\mathbf{u})-\lambda|}{|\lambda|}\
\lesssim\
\max_{\lambda_{i}\neq\lambda}\frac{|\lambda_{i}(\lambda_{i}-\lambda)|}{\lambda^{2}}\,\varepsilon^{2}.$

###### Proof.

This follows from (17), (18), and (19), using Lemma 2 with $p(t)=t-\lambda$
for (17) and $p(t)=t\,(t-\lambda)$ for (19). The bound for (17) is a specific
case of (18), when $\tau=0$. ∎

We mention that the bound (20) has a slightly improved version, as follows
(cf., e.g., [15, Thm. 2.1] for the smallest eigenvalue). For this context we
introduce $\widetilde{\mathbf{e}}=\frac{\mathbf{e}}{\|\mathbf{e}\|}$ of unit
length, $\mathbf{u}=\mathbf{x}+\mathbf{e}$ as before, and
$\widetilde{\mathbf{u}}=\frac{\mathbf{u}}{\|\mathbf{u}\|}$, where
$\|\mathbf{u}\|^{2}=1+\varepsilon^{2}$. We can decompose
$\widetilde{\mathbf{u}}=\cos(\widetilde{\mathbf{u}},\mathbf{x})\,\mathbf{x}+\sin(\widetilde{\mathbf{u}},\mathbf{x})\,\widetilde{\mathbf{e}}$.
Then an easy computation gives
$\widetilde{\mathbf{u}}^{T}\\!A\widetilde{\mathbf{u}}-\lambda=\sin^{2}(\widetilde{\mathbf{u}},\mathbf{x})\
\widetilde{\mathbf{e}}^{T}\\!(A-\lambda I)\,\widetilde{\mathbf{e}}$ and
therefore

$|\widetilde{\mathbf{u}}^{T}\\!A\widetilde{\mathbf{u}}-\lambda|\leq\max_{\lambda_{i}\neq\lambda}|\lambda_{i}-\lambda|\,\sin^{2}(\widetilde{\mathbf{u}},\mathbf{x}).$
(22)

To connect the approximate bound (20) to (22), we note that
$\sin^{2}(\widetilde{\mathbf{u}},\mathbf{x})=\frac{\varepsilon^{2}}{1+\varepsilon^{2}}$,
which is asymptotically equal to $\varepsilon^{2}$. The bound (22) is exact
(i.e., not an asymptotic bound as in Proposition 4) and sharp, but it is
asymptotically equal to our expression.

An exact bound for the error of the harmonic Rayleigh quotient with target is
provided in [14, Thm. 5.2]. Under certain assumptions mentioned in [14] for
the target $\tau$ and $\varepsilon$, with
$\mathbf{x}=\mathbf{x}+\varepsilon\,\widetilde{\mathbf{e}}$, it is shown that

$\frac{\mathbf{u}^{T}A(A-\tau I)\,\mathbf{u}}{\mathbf{u}^{T}(A-\tau
I)\,\mathbf{u}}-\lambda=\varepsilon^{2}\frac{\widetilde{\mathbf{e}}^{T}(A-\lambda
I)(A-\tau
I)\,\widetilde{\mathbf{e}}}{(\lambda-\tau)+\varepsilon^{2}\,(\widetilde{\mathbf{e}}^{T}A\widetilde{\mathbf{e}}-\tau)}\leq\varepsilon^{2}\,\max_{\lambda_{i}\neq\lambda}\frac{(\lambda_{i}-\lambda)(\lambda_{i}-\tau)}{(\lambda-\tau)+\varepsilon^{2}(\lambda_{i}-\tau)}.$

Since we do not make any assumptions on the target or the value of
$\|\mathbf{e}\|$, our upper bound (21) includes absolute values. Discarding
the $\varepsilon^{2}$-term in the denominator yields an approximate upper
bound accurate to $\mathcal{O}(\varepsilon^{4})$-terms. In conclusion, (21) is
less sharp in some situations but more general, and asymptotically the same as
the result of [14]. Another good reason to consider (20) and (21) is that
approximate bounds for the error of the homogeneous Rayleigh quotient can also
be obtained via the second-order approximation of the quotient, as we have
done in Proposition 4. At the end of this section, we will also point out that
the three approximate bounds have the same form.

###  3.2 Homogeneous Rayleigh quotient

We now determine approximate bounds for the sensitivity of the homogeneous
Rayleigh quotient.

######  Proposition 5.

Suppose $\mathbf{u}=\mathbf{x}+\mathbf{e}$ is an approximate eigenvector
corresponding to a simple eigenvalue $\lambda\neq 0$ of a symmetric $A$, with
$\|\mathbf{x}\|=1$, $\mathbf{e}\perp\mathbf{x}$, and
$\varepsilon=\|\mathbf{e}\|$. Then, up to
$\mathcal{O}(\varepsilon^{4})$-terms, for the sensitivity of the homogeneous
Rayleigh quotient (as function of $\mathbf{u}$) it holds that:

$\frac{\varepsilon^{2}}{\lambda^{2}+1}\,\min_{\lambda_{i}\neq\lambda}|(\lambda_{i}-\lambda)(\lambda_{i}+\lambda^{-1})|\,\lesssim\,\frac{|\alpha(\mathbf{u})-\lambda|}{|\lambda|}\,\lesssim\,\frac{\varepsilon^{2}}{\lambda^{2}+1}\,\max_{\lambda_{i}\neq\lambda}|(\lambda_{i}-\lambda)(\lambda_{i}+\lambda^{-1})|.$

###### Proof.

From Proposition 1(iii), the homogeneous Rayleigh quotient can be written as

$\alpha=\sqrt{\omega^{2}+1}-\omega,\qquad\omega=(2\
\mathbf{u}^{T}\\!A\mathbf{u})^{-1}\,(\mathbf{u}^{T}\mathbf{u}-\mathbf{u}^{T}\\!A^{2}\mathbf{u}).$

We can express $\omega$ in terms of $\lambda$ and $\mathbf{e}$ as

$\omega=\frac{1-\lambda^{2}+\|\mathbf{e}\|^{2}-\|A\mathbf{e}\|^{2}}{2\,(\lambda+\mathbf{e}^{T}A\mathbf{e})}.$

Furthermore, when $\|\mathbf{e}\|=\varepsilon\to 0$,

$\omega=\tfrac{1}{2}\,\lambda^{-1}\,(1-\lambda^{-1}\,\mathbf{e}^{T}\\!A\mathbf{e})\,(1-\lambda^{2}+\|\mathbf{e}\|^{2}-\|A\mathbf{e}\|^{2})+\mathcal{O}(\varepsilon^{4}).$

For the square root, we use the first-order approximation
$\sqrt{1+t}=1+\frac{1}{2}t+\mathcal{O}(t^{2})$. Discarding
$\mathcal{O}(\varepsilon^{4})$-terms yields

$\displaystyle\sqrt{\omega^{2}+1}$
$\displaystyle=\tfrac{1}{2}\,\lambda^{-1}\,(1+\lambda^{-1}\,\mathbf{e}^{T}\\!A\mathbf{e})^{-1}\sqrt{(1-\lambda^{2}+\|\mathbf{e}\|^{2}-\|A\mathbf{e}\|^{2})^{2}+4\,(\lambda+\mathbf{e}^{T}A\mathbf{e})^{2}}$
$\displaystyle=\tfrac{1}{2}\lambda^{-1}(1-\lambda^{-1}\mathbf{e}^{T}A\mathbf{e})\big{(}1+\lambda^{2}+\tfrac{1-\lambda^{2}}{1+\lambda^{2}}(\|\mathbf{e}\|^{2}-\|A\mathbf{e}\|^{2})+\tfrac{4\lambda}{1+\lambda^{2}}\,\mathbf{e}^{T}\\!A\mathbf{e}\big{)}.$

Simple computations show that

$\displaystyle\alpha$
$\displaystyle=(1-\lambda^{-1}\mathbf{e}^{T}A\mathbf{e})\Big{(}\lambda-\frac{\lambda}{\lambda^{2}+1}(\|\mathbf{e}\|^{2}-\|A\mathbf{e}\|^{2})+\frac{2}{\lambda^{2}+1}\mathbf{e}^{T}\\!A\mathbf{e}\Big{)}+\mathcal{O}(\varepsilon^{4})$
$\displaystyle=\lambda+\frac{\lambda}{\lambda^{2}+1}\,\mathbf{e}^{T}(A+\lambda^{-1}I)(A-\lambda
I)\,\mathbf{e}+\mathcal{O}(\varepsilon^{4}).$ (23)

The thesis follows from Lemma 2, with the polynomial
$p(t)=(t+\lambda^{-1})(t-\lambda)$.

An alternative method to derive this result is via the Implicit Function
Theorem. From (14), for the pair $(\alpha(\mathbf{u}),\mathbf{u})$, we have
the implicit equation
$F(\alpha(\mathbf{u}),\mathbf{u}):=\mathbf{u}^{T}(A-\alpha
I)(A+\alpha^{-1}I)\,\mathbf{u}=0$. Then we know that

$\nabla\alpha(\mathbf{u})=-\Big{(}\frac{\partial
F}{\partial\alpha}\Big{)}^{-1}\,\frac{\partial F}{\partial\mathbf{u}}\\\ =2\
[\mathbf{u}^{T}(A+\alpha^{-2}A)\,\mathbf{u}]^{-1}\ (A-\alpha
I)(A+\alpha^{-1}I)\,\mathbf{u}.$

Since $\nabla\alpha(\mathbf{x})=\mathbf{0}$, we consider the second-order
approximation of the homogeneous Rayleigh quotient
$\alpha(\mathbf{x}+\mathbf{e})\approx\alpha(\mathbf{x})+\frac{1}{2}\,\mathbf{e}^{T}\nabla^{2}\alpha(\mathbf{x})\,\mathbf{e}$.
The action of $\nabla^{2}\alpha(\mathbf{x})$ can be obtained from the first-
order approximation of the gradient in $\mathbf{x}$:

$\displaystyle\nabla^{2}\alpha(\mathbf{x})\,\mathbf{e}$
$\displaystyle=\nabla\alpha(\mathbf{x}+\mathbf{e})-\nabla\alpha(\mathbf{x})+\mathcal{O}(\varepsilon^{2})$
$\displaystyle=\frac{2\,\lambda}{1+\lambda^{2}}\ (A-\lambda
I)(A+\lambda^{-1}I)\,\mathbf{e}+\mathcal{O}(\varepsilon^{2}).$

This also leads to the expression in (3.2). ∎

Interestingly, by replacing $\tau=-\lambda^{-1}$ in the approximation of the
harmonic Rayleigh quotient with target (19), we get the approximation (3.2)
for the homogeneous Rayleigh quotient.

We now compare the different asymptotic bounds. With a little abuse of
notation, let $\theta(\mathbf{u})$ be any of the three Rayleigh quotients in
$\mathbf{u}$. Then the previous results can be summarized by

$\frac{\varepsilon^{2}}{|\lambda|}\
\min_{\lambda_{i}\neq\lambda}|p_{\lambda}(\lambda_{i})|\ \ \lesssim\ \
\frac{|\theta(\mathbf{u})-\lambda|}{|\lambda|}\ \ \lesssim\ \
\frac{\varepsilon^{2}}{|\lambda|}\
\max_{\lambda_{i}\neq\lambda}|p_{\lambda}(\lambda_{i})|,$

where the different polynomials $p_{\lambda}(\cdot)$ are reported in Table 1,
along with their pointwise asymptotic behavior for $\lambda\to 0$ and
$\lambda\to\pm\infty$. This table provides a naive explanation of what we will
observe in the experiments: the homogeneous Rayleigh quotient tends to have a
similar behavior to the Rayleigh quotient for small eigenvalues, and it
follows the harmonic Rayleigh quotient for large eigenvalues.

Table 1: Polynomials $p_{\lambda}(\cdot)$ appearing in the asymptotic bounds for the relative error of the various Rayleigh quotients. | $p_{\lambda}(t)$ | $\lambda\to 0$ | $\lambda\to\pm\infty$
---|---|---|---
RQ | $t-\lambda$ | $\sim t$ | $\sim-\lambda$
Harmonic RQ | $\frac{1}{\lambda}\,t\,(t-\lambda)$ | $\sim\frac{t^{2}}{\lambda}$ | $\sim-t$
Homogeneous RQ | $\frac{1}{\lambda^{2}+1}(t-\lambda)(\lambda t+1)$ | $\sim t$ | $\sim-t$

We remark that, for the homogeneous Rayleigh quotient,
$p_{\lambda}=p_{-\lambda^{-1}}$. In the special case that both $\lambda$ and
its anti-reciprocal $-\lambda^{-1}$ lie in the spectrum of $A$, the asymptotic
lower bound for the homogeneous Rayleigh quotient will be zero, since in
general $\lambda\neq-\lambda^{-1}$ and both eigenvalues are roots of
$p_{\lambda}$.

For the Rayleigh quotient, we have that
$\max_{\lambda_{i}\neq\lambda}|\lambda_{i}-\lambda|\in\\{|\lambda-\lambda_{1}|,\,|\lambda-\lambda_{n}|\\}$.
In contrast, the upper bounds for the harmonic Rayleigh quotient or the
homogeneous Rayleigh quotient are more complicated. Viewing $p_{\lambda}$ as a
continuous function on $[\lambda_{1},\lambda_{n}]$, the maximum of the
corresponding two upper bounds may be attained at the vertex of the parabola
$p_{\lambda}$, or at the boundary points $\lambda_{1}$ and $\lambda_{n}$.
Thus, in view of the discrete spectrum, we know that

$\max_{\lambda_{i}\neq\lambda}|p_{\lambda}(\lambda_{i})|\leq\max_{\lambda_{i}\in\\{\lambda_{1},\
\lambda_{n},\ \overline{\lambda}\\}}|p_{\lambda}(\lambda_{i})|,$ (24)

where $\overline{\lambda}=\frac{1}{2}\lambda$ for the harmonic Rayleigh
quotient, while $\overline{\lambda}=\frac{1}{2}(\lambda-\lambda^{-1})$ for the
homogeneous Rayleigh quotient. Although $p_{\lambda}(t)$ always has a factor
of $t-\lambda$, due to the quadratic nature of the upper bound for the
harmonic and the homogeneous Rayleigh quotient, we cannot make any a priori
comparison with the bound for the Rayleigh quotient; without additional
information about the spectrum, we cannot improve on these considerations. The
behavior of the relative errors and upper bounds for the various Rayleigh
quotients is shown in the next section.

###  3.3 Comparison of the Rayleigh quotients

Without loss of generality, we may assume that
$A=\text{diag}(\lambda_{1},\dots,\lambda_{n})$. We consider a family of
$100\times 100$ positive definite diagonal matrices where the eigenvalues have
uniform distribution on the interval $[0,2\sigma]$, i.e., $\lambda_{i}\sim
U(0,2\sigma)$, $\sigma>0$. The second family contains $100\times 100$
indefinite diagonal matrices where the eigenvalues have Gaussian distribution
with mean $0$ and variance $\sigma^{2}$, which we indicate by
$\lambda_{i}\sim\mathcal{N}(0,\sigma)$, $\sigma>0$. Since the eigenvalues are
drawn from continuous probability distributions, there is zero probability
that two eigenvalues are equal, so the eigenvectors are well defined with
probability one.

An eigenpair of a generic diagonal $A$ is $(\lambda,\mathbf{x})$, where
$\mathbf{x}$ is one of the vectors of the canonical basis of
$\mathbb{R}^{100}$. To compute a random perturbation
$\mathbf{u}=\mathbf{x}+\varepsilon\,\mathbf{e}$, such that $\|\mathbf{e}\|=1$
and $\mathbf{x}\perp\mathbf{e}$, we start from a vector with random Gaussian
components, project it onto the orthogonal complement of $\mathbf{x}$ and
normalize to get $\mathbf{e}$. Finally we normalize $\mathbf{u}$.

First, we study the behavior of the different Rayleigh quotients when
perturbed, with $\sigma\in\\{0.5,1,5,10\\}$ and $\varepsilon=0.001$. For each
eigenvalue $\mathbf{x}$ we draw $100$ random perturbations $\mathbf{u}_{j}$,
and take the maximum relative error. The quantity of interest for, e.g., the
standard Rayleigh quotient is

$\max_{j}\frac{|\theta(\mathbf{u}_{j})-\lambda|}{|\lambda|}.$

Figure 1 shows the maximum relative error as a function of the estimated
eigenvalue $\lambda$.

Figure 1: Sensitivity of Rayleigh quotients (RQ) for different spectra. The
label of each graph indicates the distribution of the $100$ eigenvalues. Each
eigenvalue is plotted against its corresponding maximum relative error, in
logarithmic scale.

In these examples, there are peaks close to zero since we consider relative
errors; therefore, the Rayleigh quotients are more sensitive when small
eigenvalues are approximated. In particular, the harmonic Rayleigh quotient is
always more sensitive than the other two for small eigenvalues. By looking at
the two families separately, Gaussian and uniformly distributed, we notice
that as $\sigma$ increases, the curves of the Rayleigh quotient and the
harmonic Rayleigh quotient remain almost unchanged. In the uniformly
distributed case, when $\sigma=0.5$, the homogeneous Rayleigh quotient closely
follows the Rayleigh quotient in the first half of the spectrum and slightly
departs in the second half. Starting from $\sigma=1$, the curve gets
increasingly closer to the harmonic Rayleigh quotient until the two curves are
very similar, except for the smaller eigenvalues, where the relative error of
the homogeneous Rayleigh quotient is lower than the one of the harmonic
Rayleigh quotient. This trend can be partially explained by the asymptotic
behavior of the polynomials in Table 1: the homogeneous and the harmonic
Rayleigh quotients have the same behavior for large eigenvalues, while the
homogeneous Rayleigh quotient tends to be similar to the standard Rayleigh
quotient for small eigenvalues.

We observe the same characteristics in the Gaussian family, with the
interesting addition that when $\sigma=1$, the homogeneous Rayleigh quotient
is even less sensitive than the harmonic Rayleigh quotient. We have also
considered different magnitudes of perturbation $\varepsilon$, where we
observe that the relative errors increase as $\varepsilon$ increases, but the
relative positions of the Rayleigh quotients remain the same.

Finally, we look at the asymptotic bounds of Section 3, when
$\lambda_{i}\sim\mathcal{N}(0,1)$ and $\lambda_{i}\sim U(0,2)$. Instead of
just showing the maximum relative error, we draw $10$ random perturbations for
each eigenvector, and plot the corresponding relative errors. The results are
presented in Figure 2.

(a) Eigenvalues drawn from $\mathcal{N}(0,1)$ (b) Eigenvalues drawn from
$U(0,2)$

Figure 2: Sensitivity of the Rayleigh quotients for two different spectra.
Each eigenvalue is plotted against their $10$ corresponding relative errors,
one for each perturbation of size $\varepsilon=0.001$. The asymptotic upper
and lower bounds are shown as red lines. Points that fall outside these bounds
are plotted in red.

We observe that the lower bounds are more erratic than the upper bounds. This
can be explained as follows. In Section 3.2 we have remarked that all bounds
depend on a continuous nonnegative function $|p_{\lambda}(t)|$, which has a
zero in $t=\lambda$ (cf. Table 1). Therefore, given a $\lambda_{j}$, the
minimizer
$\operatorname*{argmin}_{\lambda_{i}\neq\lambda}|p_{\lambda}(\lambda_{i})|$ is
typically either $\lambda_{j-1}$ or $\lambda_{j+1}$; the maximizer is
generally either $\lambda_{1}$ or $\lambda_{n}$ (cf. Table 1). For this
reason, the lower bound shows more variability. As already mentioned, while
this statement holds for the Rayleigh quotient, this is not always true for
the harmonic and the homogeneous Rayleigh quotients.

Regarding the behavior of the relative errors, the plots in Figure 2
reasonably reflect the ones in Figure 1. For the uniformly distributed
eigenvalues, the upper bound seems to be sharp at the extremes of the
spectrum, while for the Gaussian family it is sharper only close to the
smallest eigenvalues (in magnitude). The upper bound for the Rayleigh quotient
is also sharp for the largest eigenvalues. While the lower bound of the
uniformly distributed family poorly reflects the behavior of the relative
error, in the Gaussian family it is tighter. In particular, it is very
accurate for the homogeneous Rayleigh quotient. In Figure 1, we have already
remarked that the maximum relative error in the homogeneous Rayleigh quotient
is lower than for the other Rayleigh quotients. Now we see that for some
perturbations, the sensitivity can be even much lower.

So far we have discussed the properties of the homogeneous Rayleigh quotient
and compared it to the well-known standard and harmonic Rayleigh quotient. Now
we extend the homogeneous Rayleigh quotient for the generalized eigenvalue
problem and present an application of Rayleigh quotients to gradient methods
for unconstrained optimization problems.

##  4 The generalized eigenvalue problem

Let us consider the generalized eigenvalue problem (GEP) $A\mathbf{x}=\lambda
B\,\mathbf{x}$, where $A$ is symmetric (definite or indefinite), and $B$ is
SPD. This problem has $n$ real eigenvalues
$\lambda_{1}\leq\cdots\leq\lambda_{n}$. Although a common generalized Rayleigh
quotient is
$\frac{\mathbf{u}^{T}\\!A\mathbf{u}}{\mathbf{u}^{T}\\!B\,\mathbf{u}}$ (see,
e.g., [1, Ch. 15]), it turns out that, in our context, the most relevant
quantities are

$\theta=\frac{\mathbf{u}^{T}\\!AB\,\mathbf{u}}{\mathbf{u}^{T}\\!B^{2}\mathbf{u}},\qquad\widetilde{\theta}=\frac{\mathbf{u}^{T}\\!A^{2}\mathbf{u}}{\mathbf{u}^{T}\\!AB\,\mathbf{u}},$
(25)

for nonzero $\mathbf{u}^{T}\\!AB\,\mathbf{u}$. The Rayleigh quotient $\theta$
satisfies $A\mathbf{u}-\theta B\,\mathbf{u}\perp B\,\mathbf{u}$, and,
equivalently, it is the solution to $\min_{\gamma}\|A\mathbf{u}-\gamma
B\,\mathbf{u}\|$. The harmonic Rayleigh quotient $\widetilde{\theta}$
satisfies $A\mathbf{u}-\theta B\,\mathbf{u}\perp A\mathbf{u}$ and solves
$\min_{\gamma}\|\gamma^{-1}A\mathbf{u}-B\,\mathbf{u}\|$. The extension of the
homogeneous Rayleigh quotient for the GEP is rather straightforward. It is the
solution to

$\min_{(\alpha_{1},\,\alpha_{2})\,\in\,\mathbb{P}}\
\frac{\|\alpha_{1}B\,\mathbf{u}-\alpha_{2}\,A\mathbf{u}\|}{{\color[rgb]{0,0,0}\sqrt{\alpha_{1}^{2}+\alpha_{2}^{2}}}}=\min_{\alpha_{1}^{2}+\alpha_{2}^{2}=1}\
\|[B\,\mathbf{u}\ \
-\\!\\!A\mathbf{u}]{\mbox{\scriptsize$\left[\\!\\!\begin{array}[]{c}\alpha_{1}\\\
\alpha_{2}\end{array}\\!\\!\right]$}}\|.$ (26)

As for the standard eigenvalue problem, this amounts to solve a reduced SVD of
an $n\times 2$ matrix, or an eigenvalue problem involving a $2\times 2$
matrix.

The next result is a generalization of Proposition 1. Let $C=[B\,\mathbf{u}\ \
-\\!\\!A\mathbf{u}]$.

######  Proposition 6.

Suppose that $\mathbf{u}^{T}\\!AB\,\mathbf{u}\neq 0$ and denote by $\mu$ the
smallest eigenvalue of $C^{T}\\!C$. Then the following properties hold.

  * (i)

$\mu$ is a simple eigenvalue of $C^{T}\\!C$, or equivalently $\sqrt{\mu}$ is a
simple singular value of $C$; its corresponding eigenvector
$[\alpha_{1},\,\alpha_{2}]^{T}$ is the unique minimizer of (26), up to
orthogonal transformations.

  * (ii)

$\mu=\tfrac{1}{2}\,\big{[}\mathbf{u}^{T}\\!A^{2}\mathbf{u}+\mathbf{u}^{T}\\!B^{2}\mathbf{u}-\sqrt{\smash[b]{(\mathbf{u}^{T}\\!A^{2}\mathbf{u}-\mathbf{u}^{T}\\!B^{2}\mathbf{u})^{2}+4\,(\mathbf{u}^{T}\\!AB\,\mathbf{u})^{2}}}\
\big{]}.$

  * (iii)

For the generalized homogeneous Rayleigh quotient $\alpha$ we have

$\displaystyle\alpha$ $\displaystyle=(2\
\mathbf{u}^{T}\\!AB\,\mathbf{u})^{-1}\,\big{[}\,\mathbf{u}^{T}\\!A^{2}\,\mathbf{u}-\mathbf{u}^{T}\\!B^{2}\mathbf{u}+\sqrt{(\mathbf{u}^{T}\\!A^{2}\,\mathbf{u}-\mathbf{u}^{T}\\!B^{2}\mathbf{u})^{2}+4\,(\mathbf{u}^{T}\\!AB\,\mathbf{u})^{2}}\,\big{]}.$
(27)

  * (iv)

$\mu=0$ if and only if $\mathbf{u}$ is an eigenvector. If $\mathbf{u}$ is an
eigenvector, then the corresponding homogeneous Rayleigh quotient $\alpha$ is
the corresponding eigenvalue.

  * (v)

$0\leq\mu<\min(\mathbf{u}^{T}\\!A^{2}\mathbf{u},\
\mathbf{u}^{T}\\!B^{2}\mathbf{u})$.

  * (vi)

We have the following inequalities

$\begin{cases}\theta\leq\alpha\leq\widetilde{\theta}&\text{if}\quad\mathbf{u}^{T}\\!AB\,\mathbf{u}>0,\\\
\widetilde{\theta}\leq\alpha\leq\theta&\text{if}\quad\mathbf{u}^{T}\\!AB\,\mathbf{u}<0.\end{cases}$

###### Proof.

This proof follows the exact same lines as those of Proposition 1. ∎

As in Section 2, we can show that the generalized homogeneous Rayleigh
quotient is a solution to a quadratic equation.

######  Proposition 7.

Let $\mathbf{u}^{T}\\!AB\,\mathbf{u}\neq 0$. The generalized homogeneous
Rayleigh quotient ((iii)) is the solution to

$(\mathbf{u}^{T}\\!AB\,\mathbf{u})\,\alpha^{2}+(\mathbf{u}^{T}\\!B^{2}\mathbf{u}-\mathbf{u}^{T}\\!A^{2}\mathbf{u})\,\alpha-(\mathbf{u}^{T}\\!AB\,\mathbf{u})=0,$
(28)

which satisfies $\alpha\ \mathbf{u}^{T}\\!AB\,\mathbf{u}>0$.

###### Proof.

The proof is similar to that of Proposition 3. ∎

We derive a Galerkin condition from Proposition 7 as follows. Since
$\mathbf{u}^{T}\\!AB\,\mathbf{u}=\mathbf{u}^{T}\\!B\,A\mathbf{u}$, (28) is
equivalent to $\mathbf{u}^{T}\\!(\alpha A+B)(A-\alpha B)\,\mathbf{u}=0$.
Therefore, we can write the corresponding nonlinear Galerkin condition in
homogeneous coordinates as

$(\alpha_{2}A-\alpha_{1}B)\,\mathbf{u}\perp(\alpha_{1}A+\alpha_{2}B)\,\mathbf{u}.$

Finally, we discuss a bound for the chordal metric, analogous to (16). For an
SPD $B$, as a generalization to [1, Thm. 4.5.1], one can show that there
exists a generalized eigenvalue $\lambda=\lambda(A,B)$ of $(A,B)$ such that
$|\lambda-\alpha|\leq(\lambda_{1}(B)\|\mathbf{u}\|)^{-1}\|(A-\alpha\,B)\,\mathbf{u}\|$,
where $\lambda_{1}(B)$ is the smallest eigenvalue of $B$. Indeed, if $\alpha$
is not an eigenvalue of $(A,B)$,

$\displaystyle\|\mathbf{u}\|$
$\displaystyle=\|(A-\alpha\,B)^{-1}\,(A-\alpha\,B)\,\mathbf{u}\|\geq\lambda_{1}(B)\,\|(B^{-1}A-\alpha\,I)^{-1}\|\,\|(A-\alpha\,B)\,\mathbf{u}\|$
$\displaystyle\geq\lambda_{1}(B)\,\min_{i}|\lambda_{i}(A,B)-\alpha|\cdot\|(A-\alpha\,B)\,\mathbf{u}\|.$

Therefore, for the chordal distance between $\lambda$ and $\alpha$ we have

$\frac{|\lambda-\alpha|}{\sqrt{1+\lambda^{2}}\,\sqrt{1+\alpha^{2}}}\leq\frac{1}{\lambda_{1}(B)}\cdot\frac{1}{\sqrt{1+\lambda^{2}}\,\|\mathbf{u}\|}\cdot\frac{\|(A-\alpha\,B)\,\mathbf{u}\|}{\sqrt{1+\alpha^{2}}}.$

This upper bound is minimized for $\alpha$ equal to the generalized
homogeneous Rayleigh quotient. We notice that this inequality differs from
(16) by the factor $(\lambda_{1}(B))^{-1}$.

In the next section, we return to the standard homogeneous Rayleigh quotient,
and study its use as stepsize in gradient methods.

##  5 A homogeneous stepsize for gradient methods

In gradient methods for nonlinear optimization, inverse Rayleigh quotients are
popular choices for stepsizes. We refer to [16] for a nice recent review about
steplength selection. We exploit the inverse of the homogeneous Rayleigh
quotient as a new _homogeneous stepsize_ for gradient-type methods. Consider
the unconstrained minimization of a smooth function
$f:\mathbb{R}^{n}\to\mathbb{R}$,

$\min_{\mathbf{x}\in\mathbb{R}^{n}}\,f(\mathbf{x}).$

Gradient methods are of the form

$\mathbf{x}_{k+1}=\mathbf{x}_{k}-\beta_{k}\,\mathbf{g}_{k}=\mathbf{x}_{k}-\alpha_{k}^{-1}\,\mathbf{g}_{k},$

where $\mathbf{g}_{k}=\nabla f(\mathbf{x}_{k})$, $\beta_{k}>0$ is the
stepsize, and $\alpha_{k}$ the inverse stepsize. As usual we write
$\mathbf{s}_{k-1}=\mathbf{x}_{k}-\mathbf{x}_{k-1}$ and
$\mathbf{y}_{k-1}=\mathbf{g}_{k}-\mathbf{g}_{k-1}$.

In [5], the Barzilai–Borwein stepsizes

$\beta_{k}^{\rm{BB1}}=\frac{\mathbf{s}_{k-1}^{T}\,\mathbf{s}_{k-1}}{\mathbf{s}_{k-1}^{T}\,\mathbf{y}_{k-1}},\qquad\beta_{k}^{\rm{BB2}}=\frac{\mathbf{s}_{k-1}^{T}\,\mathbf{y}_{k-1}}{\mathbf{y}_{k-1}^{T}\,\mathbf{y}_{k-1}}$
(29)

have been introduced. The motivation is that, when we approximate the Hessian
in $\mathbf{x}_{k-1}$ by a scalar multiple of the identity, i.e.,
$\nabla^{2}f(\mathbf{x}_{k-1})\approx\gamma\,I$, the corresponding inverse
stepsizes satisfy the following least squares secant conditions:

$\alpha_{k}^{\rm{BB1}}=\operatorname*{argmin}_{\gamma}\|\mathbf{y}_{k-1}-\gamma\,\mathbf{s}_{k-1}\|,\qquad\alpha_{k}^{\rm{BB2}}=\operatorname*{argmin}_{\gamma}\|\gamma^{-1}\,\mathbf{y}_{k-1}-\mathbf{s}_{k-1}\|,$
(30)

with $\alpha_{k}^{\rm{BB1}}=(\beta_{k}^{\rm{BB1}})^{-1}$ and
$\alpha_{k}^{\rm{BB2}}=(\beta_{k}^{\rm{BB2}})^{-1}$. Moreover, both stepsizes
(29) can be seen as inverse Rayleigh quotients of a certain matrix $H_{k}$ at
$\mathbf{s}_{k-1}$. In fact, for any $f$ it holds

$\mathbf{y}_{k-1}=\nabla f(\mathbf{x}_{k})-\nabla
f(\mathbf{x}_{k-1})=H_{k}\,(\mathbf{x}_{k}-\mathbf{x}_{k-1})=H_{k}\,\mathbf{s}_{k-1},$

where
$H_{k}:=\int_{0}^{1}\nabla^{2}f((1-t)\,\mathbf{x}_{k-1}+t\,\mathbf{x}_{k})\,dt$
is an average Hessian on the line piece between $\mathbf{x}_{k-1}$ and
$\mathbf{x}_{k}$. From this relation, it is easy to see that the minimum
residual conditions (3) and (7) are equivalent to the secant conditions (30)
for $\mathbf{u}=\mathbf{s}_{k-1}$ and $A=H_{k}$.

We introduce the homogeneous BB stepsize (HBB) as the inverse of the
homogeneous Rayleigh quotient of $H_{k}$ in $\mathbf{s}_{k-1}$. The HBB
stepsize is given by the quotient $\beta_{k}^{\rm
HBB}=\alpha_{2,k}/\alpha_{1,k}$, where the pair $(\alpha_{1,k}$,
$\alpha_{2,k})$ solves

$\operatorname*{argmin}_{\alpha_{1}^{2}+\alpha^{2}_{2}=1}\
\|\alpha_{1}\,\mathbf{s}_{k-1}-\alpha_{2}\,\mathbf{y}_{k-1}\|.$

As in (30), this is equivalent to the minimum residual condition (8) for
$\mathbf{u}=\mathbf{s}_{k-1}$ and $A=H_{k}$. Therefore (cf. Proposition 1)

$\beta_{k}^{\rm
HBB}=\frac{\|\mathbf{s}_{k-1}\|^{2}-\|\mathbf{y}_{k-1}\|^{2}+\sqrt{\smash[b]{(\|\mathbf{s}_{k-1}\|^{2}-\|\mathbf{y}_{k-1}\|^{2})^{2}+4\,(\mathbf{s}_{k-1}^{T}\mathbf{y}_{k-1})^{2}}}}{2\
\mathbf{s}_{k-1}^{T}\mathbf{y}_{k-1}}.$ (31)

As mentioned in Section 1, this stepsize has been proposed first by Li, Zhang,
and Xia [12, 13], obtained from a total least squares secant condition. In
this section, we will also introduce an alternating variant of this HBB step.

In Section 5.1 we will carry out some experiments to test the behavior of the
HBB stepsize when plugged into a gradient method for general differentiable
functions. A pseudocode of this method is provided in Algorithm 1.

Algorithm 1 A homogeneous gradient method for minimization of general
functions

Input: Continuous differentiable function $f$, initial guess $\mathbf{x}_{0}$,
initial stepsize $\beta_{0}>0$, tolerance tol; safeguarding parameters
$\beta_{\max}>\beta_{\min}>0$; line search parameters $c_{\rm{ls}}$,
$\sigma_{\rm{ls}}\in(0,1)$; memory integer $M>0$
Output: Approximation to minimizer
$\operatorname*{argmin}_{\mathbf{x}}f(\mathbf{x})$
1: Set $\mathbf{g}_{0}=\nabla f(\mathbf{x}_{0})$ for $k=0,1,\dots$ 2: M
$\nu_{k}=\beta_{k}$, $f_{\text{ref}}=\max\,\\{\,f(\mathbf{x}_{k-j})\,:\,0\leq
j\leq\min(k,M-1)\,\\}$ 3: M while
$f(\mathbf{x}_{k}-\nu_{k}\,\mathbf{g}_{k})>f_{\text{ref}}-c_{\rm{ls}}\,\nu_{k}\,\|\mathbf{g}_{k}\|^{2}$
do $\nu_{k}=\sigma_{\rm{ls}}\,\nu_{k}$ end 4: M Set
$\mathbf{s}_{k}=-\nu_{k}\,\mathbf{g}_{k}$ and update
$\mathbf{x}_{k+1}=\mathbf{x}_{k}+\mathbf{s}_{k}$ 5: M Compute the gradient
$\mathbf{g}_{k+1}=\nabla f(\mathbf{x}_{k+1})$ 6: M if
$\|\mathbf{g}_{k+1}\|\leq{\sf tol}\cdot\|\mathbf{g}_{0}\|$, return, end 7: M
$\mathbf{y}_{k}=\mathbf{g}_{k+1}-\mathbf{g}_{k}$ 8: M if
$\mathbf{s}_{k}^{T}\mathbf{y}_{k}<0$, set
$\beta_{k+1}=\max(\min(\|\mathbf{g}_{k+1}\|_{2}^{-1},\,10^{5}),\ 1)$ 9: M else
compute the homogeneous stepsize $\beta_{k+1}$ (31) end 10: M Set
$\beta_{k+1}=\min(\max(\beta_{k+1},\,\beta_{\min}),\ \beta_{\max})$

This algorithm is similar to the one proposed in [16, 17], with the
homogeneous stepsize (31) as the key difference in Line 9. In presence of an
uphill direction, i.e., when $\mathbf{s}_{k}^{T}\mathbf{y}_{k}<0$, the
homogeneous stepsize is negative (cf. Proposition 1 (vi)), and therefore must
be replaced with a positive quantity. We choose the new stepsize as in [17]
(cf. Line 8 of Algorithm 1). The convergence of the method is not affected by
these choices, since $\beta_{k}$ stays uniformly bounded, i.e.,
$\beta_{k}\in[\beta_{\min},\beta_{\max}]$ for all $k$. Therefore, the proof of
global convergence of Algorithm 1 can be easily adapted from [17, Thm. 2.1].
While the convergence is not affected, choosing the homogeneous stepsize as
starting steplength in the nonmonotone line search might lead to a smaller
number of backtracking steps, compared to classical BB stepsizes.

We finally remark that, as for the BB stepsizes, no line search is required
for HBB steps when the function $f$ is strictly convex quadratic, i.e.,
$f(\mathbf{x})=\tfrac{1}{2}\,\mathbf{x}^{T}A\mathbf{x}-\mathbf{b}^{T}\mathbf{x}$,
with $A$ SPD (see the results in [18, 19]). In fact, since $H_{k}\equiv A$ is
positive definite, from Proposition 1 it holds that
$\alpha_{k}\in[\alpha_{k}^{\rm{BB1}},\,\alpha_{k}^{\rm{BB2}}]$. This is
sufficient to guarantee the R-linear convergence of the corresponding gradient
method (cf. [20, Thm. 4.1] or [3, Thm. 16]).

###  5.1 Numerical experiments

This subsection is devoted to testing the use of the homogeneous stepsize HBB,
and an adaptive variant in Algorithm 1, on a set of unconstrained optimization
problems. We take some general continuously differentiable functions and the
suggested starting points therein from the collection in [17, 21, 22], as
listed in Table 2.

Table 2: Unconstrained optimization test problems. Name | Reference | Name | Reference
---|---|---|---
Diagonal 1 | [21] | Full Hessian FH1 | [21]
Diagonal 2 | [21] | Full Hessian FH2 | [21]
Diagonal 3 | [21] | Generalized Rosenbrock | [21]
Extended Beale | [21] | Generalized White and Holst | [21]
Extended Powell | [22] | Hager | [21]
Extended Rosenbrock | [22] | Perturbed quadratic | [21]
Extended White and Holst | [21] | Strictly Convex 2 | [17]

For all the test functions, the size can be $n\in\\{10^{2},10^{3},10^{4}\\}$.
The generalized Rosenbrock, generalized White and Holst and extended Powell
objective functions have been scaled by the Euclidean norm of the first
gradient.

As for the parameters of Algorithm 1, we set the choices
$\beta_{\min}=10^{-30}$, $\beta_{\max}=10^{30}$, $c_{\rm{ls}}=10^{-4}$,
$\sigma_{\rm{ls}}=\frac{1}{2}$, $M=10$, and $\beta_{0}=1$. The algorithm stops
when $\|\mathbf{g}_{k}\|\leq{\sf tol}\cdot\|\mathbf{g}_{0}\|$, with ${\sf
tol}=10^{-6}$, or when $5\cdot 10^{4}$ iterations are reached. All different
steps in Table 3 are tested, meaning that we change the stepsize choice in
Line 9 of Algorithm 1 and obtain several gradient methods. Along with the
homogeneous stepsize (HBB), BB1 and BB2, we also consider the adaptive ABB
method [23]:

$\beta_{k}^{\rm ABB}=\left\\{\begin{array}[]{ll}\min\\{\beta_{j}^{\rm BB2}\mid
j=\max\\{1,k-m\\},\dots,k\\},&\quad{\rm if}\;\beta_{k}^{\rm
BB2}<\eta\,\beta_{k}^{\rm BB1},\\\\[5.69054pt] \beta_{k}^{\rm
BB1}&\quad\text{otherwise}.\end{array}\right.$ (32)

Default parameters values are $\eta=0.8$ and $m=5$; see, e.g., [16]. This
stepsize is mostly known as $\text{ABB}_{\text{min}}$, while ABB indicates the
case for $m=0$; we choose to indicate it as ABB for any $m$ to ease the
notation. Inspired by this, we also propose a straightforward generalization
of our HBB step: the adaptive HBB method (AHBB), which takes the stepsize

$\beta_{k}^{\rm AHBB}=\left\\{\begin{array}[]{ll}\min\\{\beta_{j}^{\rm
HBB}\mid j=\max\\{1,k-m\\},\dots,k\\},&\ {\rm if}\;\beta_{k}^{\rm
BB2}<\eta\,\beta_{k}^{\rm BB1},\\\\[5.69054pt] \beta_{k}^{\rm BB1}&\
\text{otherwise}.\end{array}\right.$ (33)

Table 3: Tested stepsizes. Method | Description | Reference |
---|---|---|---
BB1 | Inverse Rayleigh quotient |  Cf. [5] |  Eq. (29)
BB2 | Inverse harmonic Rayleigh quotient |  Cf. [5] |  Eq. (29)
HBB | Inverse homogeneous Rayleigh quotient |  |  Eq. (31)
ABB |  Adaptive: BB1 or (past $m=5$) BB2 step |  Cf. [23] |  Eq. (32)
AHBB |  Adaptive: BB1 or (past $m=5$) HBB step |  |  Eq. (33)

Since the stopping criterion is based on the gradient norm, and
$\mathbf{s}_{k}=-\nu_{k}\mathbf{g}_{k}$ (cf. Line 4 in Algorithm 1), the
homogeneous stepsize requires the computation of the two inner products
$\mathbf{y}_{k}^{T}\mathbf{y}_{k}$ and $\mathbf{g}_{k}^{T}\mathbf{y}_{k}$, and
therefore it has the same cost as BB2. In general, all the studied stepsizes
in Table 3 have very similar cost; BB1 is slightly cheaper, since it does not
require the inner product $\mathbf{y}_{k}^{T}\mathbf{y}_{k}$.

Table 4 reports the number of function evaluations and the number of
iterations for each problem and stepsize. We remark that, for each problem,
all methods converged to the same stationary point. It seems that either HBB
or AHBB can be competitive or better than ABB, which is well known for its
generally good behavior. For example, HBB behaves nicely in problems FH1 and
FH2; AHBB is the best method for the generalized White and Holst function and
perturbed quadratic ($n=10^{4}$).

Table 4: Number of function evaluations (NFE) and iterations for each stepsize
and problem. Boldface items are minimal.

Problem and size | NFE | Iterations
---|---|---
|  | BB1 | BB2 | ABB | HBB | AHBB | BB1 | BB2 | ABB | HBB | AHBB
Diagonal 1 |  ($10^{2}$) | 65 | 68 | 65 | 69 | 67 | 57 | 63 | 60 | 63 | 62
Diagonal 1 |  ($10^{3}$) | 305 | 194 | 160 | 159 | 165 | 222 | 180 | 149 | 145 | 156
Diagonal 1 |  ($10^{4}$) | 761 | 433 | 302 | 346 | 310 | 491 | 409 | 289 | 323 | 297
Diagonal 2 |  ($10^{2}$) | 75 | 68 | 59 | 73 | 66 | 68 | 67 | 58 | 68 | 63
Diagonal 2 |  ($10^{3}$) | 447 | 234 | 165 | 281 | 212 | 286 | 219 | 157 | 185 | 175
Diagonal 2 |  ($10^{4}$) | 689 | 702 | 322 | 1221 | 542 | 405 | 663 | 307 | 729 | 407
Diagonal 3 |  ($10^{2}$) | 76 | 83 | 73 | 67 | 67 | 62 | 73 | 65 | 60 | 59
Diagonal 3 |  ($10^{3}$) | 311 | 200 | 152 | 151 | 148 | 206 | 184 | 138 | 137 | 137
Diagonal 3 |  ($10^{4}$) | 349 | 375 | 238 | 283 | 278 | 224 | 356 | 224 | 266 | 263
Ext Beale |  ($10^{2}$) | 50 | 34 | 33 | 33 | 33 | 45 | 30 | 29 | 27 | 29
Ext Powell |  ($10^{2}$) | 138 | 80 | 117 | 136 | 156 | 103 | 74 | 116 | 100 | 131
Ext Powell |  ($10^{3}$) | 117 | 97 | 114 | 132 | 122 | 97 | 84 | 113 | 103 | 104
Ext Powell |  ($10^{4}$) | 164 | 125 | 144 | 161 | 155 | 127 | 112 | 143 | 114 | 124
Ext Rosen |  ($10^{2}$) | 104 | 65 | 88 | 69 | 87 | 54 | 53 | 72 | 55 | 74
Ext WH |  ($10^{2}$) | 99 | 48 | 39 | 37 | 39 | 62 | 37 | 29 | 27 | 29
FH1 |  ($10^{2}$) | 623 | 304 | 337 | 256 | 376 | 408 | 279 | 310 | 221 | 347
FH2 |  ($10^{2}$) | 987 | 574 | 563 | 447 | 499 | 635 | 530 | 537 | 407 | 470
Gen Rosen |  ($10^{2}$) | 4121 | 3162 | 2996 | 4808 | 3567 | 2623 | 2932 | 2736 | 3085 | 2840
Gen Rosen |  ($10^{3}$) | 37769 | 27324 | 28624 | 37829 | 31516 | 24020 | 24972 | 26078 | 24082 | 24700
Gen WH |  ($10^{2}$) | 11040 | 8403 | 8866 | 11227 | 7216 | 6949 | 8107 | 8660 | 7123 | 5602
Hager |  ($10^{2}$) | 24 | 27 | 27 | 25 | 27 | 21 | 24 | 24 | 22 | 24
Hager |  ($10^{3}$) | 41 | 46 | 44 | 44 | 45 | 37 | 42 | 40 | 40 | 41
Hager |  ($10^{4}$) | 84 | 62 | 64 | 66 | 64 | 74 | 55 | 57 | 59 | 57
Pert quad |  ($10^{2}$) | 76 | 98 | 72 | 76 | 73 | 63 | 90 | 64 | 66 | 65
Pert quad |  ($10^{3}$) | 289 | 220 | 169 | 263 | 156 | 194 | 198 | 158 | 243 | 145
Pert quad |  ($10^{4}$) | 618 | 350 | 316 | 270 | 260 | 401 | 330 | 296 | 252 | 242
S Conv 2 |  ($10^{2}$) | 82 | 62 | 61 | 66 | 71 | 72 | 58 | 56 | 57 | 65
S Conv 2 |  ($10^{3}$) | 282 | 156 | 120 | 156 | 141 | 197 | 145 | 111 | 137 | 132
S Conv 2 |  ($10^{4}$) | 467 | 234 | 191 | 213 | 199 | 289 | 221 | 179 | 195 | 185

##  6 Conclusions

We have introduced a new Rayleigh quotient, the homogeneous Rayleigh quotient
$\alpha$, which minimizes the homogeneous residual quantity (8). We have shown
that the value of $\alpha$ lies between the standard and the harmonic Rayleigh
quotient, and so does its sensitivity when approximating extreme eigenvalues
of an SPD matrix $A$. Regarding the other eigenvalues and the indefinite case,
we have provided examples where the homogeneous Rayleigh quotient may be more
accurate than the others. In addition, in all our experiments, we have noticed
that the homogeneous Rayleigh quotient is less sensitive than one of the other
two Rayleigh quotients. To some extent, the homogeneous Rayleigh quotient
seems to leverage the two Rayleigh quotients: it is less sensitive than the
harmonic Rayleigh quotient when estimating the smallest eigenvalues, and less
sensitive than the standard Rayleigh quotient when estimating the largest
eigenvalues (in magnitude).

We have derived a nonlinear Galerkin condition for the homogeneous Rayleigh
quotient, in contrast to the linear Galerkin conditions for the standard and
the harmonic Rayleigh quotient. All the results extend to the homogeneous
Rayleigh quotient for the generalized eigenvalue problem $(A,B)$, when $B$ is
SPD and $A$ is symmetric.

Finally, we have considered the homogeneous Rayleigh quotient as inverse
stepsize (HBB) in gradient methods for unconstrained optimization problems.
This stepsize has first and independently been obtained recently from a
different angle by Li, Zhang, and Xia [12, 13]. We have also proposed the AHBB
steplength as an alternative to the ABB stepsize, based on the homogeneous
stepsize. Experiments show that this variant sometimes performs better than
the classical BB steplengths and ABB.

######  Acknowledgement.

We are grateful to the referees for their very useful suggestions which
considerably improved the quality of the paper. This work has received funding
from the European Union’s Horizon 2020 research and innovation programme under
the Marie Sklodowska-Curie grant agreement No 812912.

## References

  * Parlett [1998] B. N. Parlett, The Symmetric Eigenvalue Problem, SIAM, Philadelphia, PA, 1998\.
  * Stewart [2001] G. W. Stewart, Matrix Algorithms. Vol. II, SIAM, Philadelphia, PA, 2001\.
  * Ferrandi et al. [2023] G. Ferrandi, M. Hochstenbach, N. Krejic,  A harmonic framework for stepsize selection in gradient methods,  Comput. Opt. Appl. (2023).
  * Morgan [1991] R. Morgan,  Computing interior eigenvalues of large matrices,  Linear Algebra Appl. 154 (1991) 289–309.
  * Barzilai and Borwein [1988] J. Barzilai, J. Borwein,  Two-point step size gradient methods,  IMA J. Numer. Anal. 8 (1988) 141–148.
  * Vömel [2010] C. Vömel,  A note on harmonic Ritz values and their reciprocals,  Numer. Linear Algebra Appl. 17 (2010) 97–108.
  * Stewart and Sun [1990] G. W. Stewart, J. G. Sun, Matrix Perturbation Theory, Academic Press Inc., Boston, MA, 1990.
  * Stewart [1975] G. W. Stewart,  Gershgorin theory for the generalized eigenvalue problem $Ax=\lambda Bx$,  Math. Comput. (1975) 600–606.
  * Dedieu [1997] J.-P. Dedieu,  Condition operators, condition numbers, and condition number theorem for the generalized eigenvalue problem,  Linear Algebra Appl. 263 (1997) 1–24.
  * Dedieu and Tisseur [2003] J.-P. Dedieu, F. Tisseur,  Perturbation theory for homogeneous polynomial eigenvalue problems,  Linear Algebra Appl. 358 (2003) 71–94.
  * Hochstenbach and Notay [0708] M. E. Hochstenbach, Y. Notay,  Homogeneous Jacobi–Davidson,  Electron. Trans. Numer. Anal. 29 (2007/08) 19–30.
  * Li et al. [2022] S. Li, T. Zhang, Y. Xia,  A family of Barzilai–Borwein steplengths from the viewpoint of scaled total least squares,  Preprint arXiv:2212.05658 (2022).
  * Li and Xia [2021] S. Li, Y. Xia,  A new Barzilai–Borwein steplength from the viewpoint of total least squares,  Preprint arXiv:2107.06687 (2021).
  * Sleijpen and van den Eshof [2003] G. Sleijpen, J. van den Eshof,  On the use of harmonic Ritz pairs in approximating internal eigenpairs,  Linear Algebra Appl. 358 (2003) 115–137.
  * Sleijpen et al. [2003] G. Sleijpen, J. van den Eshof, P. Smit,  Optimal a priori error bounds for the Rayleigh–Ritz method,  Math. Comput. 72 (2003) 677–684.
  * Di Serafino et al. [2018] D. Di Serafino, V. Ruggiero, G. Toraldo, L. Zanni,  On the steplength selection in gradient methods for unconstrained optimization,  Appl. Math. Comput. 318 (2018) 176–195.
  * Raydan [1997] M. Raydan,  The Barzilai and Borwein gradient method for the large scale unconstrained minimization problem,  SIAM J. Optim. 7 (1997) 26–33.
  * Dai and Liao [2002] Y. Dai, L. Liao,  R-linear convergence of the Barzilai and Borwein gradient method,  IMA J. Numer. Anal. 22 (2002) 1–10.
  * Raydan [1993] M. Raydan,  On the Barzilai and Borwein choice of steplength for the gradient method,  IMA J. Numer. Anal. 13 (1993) 321–326.
  * Dai [2003] Y. Dai,  Alternate step gradient method,  Optimization 52 (2003) 395–415.
  * Andrei [2008] N. Andrei,  An unconstrained optimization test functions collection,  Adv. Model. Optim. 10 (2008) 147–161.
  * Moré et al. [1981] J. Moré, B. Garbow, K. Hillstrom,  Testing unconstrained optimization software,  ACM Trans. Math. Softw. 7 (1981) 17–41.
  * Frassoldati et al. [2008] G. Frassoldati, L. Zanni, G. Zanghirati,  New adaptive stepsize selections in gradient methods,  J. Ind. Manag. 4 (2008) 299.