A Fast Projected Gradient Algorithm for Quaternion Hermitian Eigenvalue Problems

1. Introduction

In 1843, Sir William Rowan Hamilton [9] introduced quaternions in an effort to expand the concept of complex numbers into spaces of higher dimensions. Quaternions and quaternion matrices play a critical role in many applications, such as quantum mechanics, computer graphics, quaternion principal component analysis (QPCA) and image processing [2,11,23,24]. Due to the non-commutative property of quaternion multiplication, the eigenvalues of quaternion matrices are distinguished into left and right types, with the right eigenvalue problem having garnered widespread attention[4,6,14,16,25].

In recent years, a series of numerical methods have been developed to compute the eigenvalues of quaternion matrices, particularly focusing on the eigenvalue problems of Hermitian matrices. These numerical methods can be broadly categorized into three classes: The first class involves direct quaternion arithmetic operations. For instance, Bunse-Gerstner proposed a quaternion QR algorithm for solving the right eigenvalue problem of quaternion matrices [1]. However, due to the complexity of quaternion arithmetic, this algorithm requires significant computational effort. The second class is based on the real or complex counterparts of quaternion matrices. By studying the real or complex counterpart structures and properties of quaternion matrices, and leveraging stable orthogonal transformations, real or complex structure-preserving methods have been developed to solve the right eigenvalue problem of quaternion Hermitian matrices [8,10]. The third class is based on the real counterparts of quaternion matrices, leading to the development of numerous structure-preserving iterative algorithms. Examples include the explicitly restarted quaternion Arnoldi Method (ERQAM) [18], designed to compute standard right eigenpairs of general quaternion matrices, and a novel quaternion power method introduced in [13] for computing the dominant standard right eigenvalue and its corresponding eigenvector. Structure-preserving methods exhibit significant advantages in terms of storage space and computational efficiency.

In the field of quaternion optimization, significant progress has been made with the generalized HR calculus (GHR) [20,21,22]. GHR leverages quaternion rotations within a general orthogonal system, offering a way to compute the derivatives and gradients of functions with quaternion variables, thereby providing a solid theoretical foundation for the development of quaternion optimization methods. Subsequently, based on the generalized HR calculus (GHR), Diao et al. [3] proposed a gradient projection algorithm for maximizing the quaternion Rayleigh quotient under unit constraints. This algorithm demonstrated good performance and contributed to the development of quaternion optimization algorithms.

In this paper, we first equivalently transform the principal eigenvalue problem of quaternion Hermitian matrices into a maximization optimization problem over the quaternion skew field. Leveraging generalized HR calculus, we propose a quaternion Nesterov’s accelerated gradient projection algorithm to solve it. Subsequently, we conduct a convergence analysis of the quaternion Nesterov’s accelerated gradient projection algorithm, proving that a real differentiable function with Lipschitz continuous gradient possesses a quadratic upper bound. Furthermore, we demonstrate that the algorithm possesses a convergence rate of $(O(1/t^2)$. Finally, we compare our algorithm with two other methods, and numerical experiments indicate that our algorithm exhibits superior performance in terms of both accuracy and time efficiency.

The rest of this paper is organized as follows. Section 2 introduces some basic notations and fundamental properties of quaternions, including definitions of quaternion modulus, similarity, and rotation, with a particular emphasis on reviewing the relevant definitions and properties of generalized HR integrals. In Section 3, we design a quaternion Nesterov accelerated gradient projection algorithm to solve the principal eigenvalue and corresponding eigenvector of quaternion Hermitian matrices. Section 4 provides a convergence analysis of the quaternion Nesterov accelerated gradient projection algorithm. In Section 5, we conduct numerical experiments to validate the proposed method. Finally, in Section 6, we summarize this paper.

2. Preliminaries

In this section, some quaternion notations and basic definitions are introduced, which will be used in the rest of the paper.

2.1. Notations

Throughout this paper, to distinguish scalars, vectors, real or complex matrices and quaternion matrices, scalars will be denoted by lower case Greek letters, e.g., $\alpha$, $\beta$, quaternions will be denoted by lowercase letters, e.g., $\mathbf{p},\mathbf{q}$, and quaternion vectors are denoted by $\mathbf{x},\mathbf{y}$, real or complex matrices will be defined by uppercase letters, e.g., $A,B$, and quaternion matrices will be denoted by bold uppercase letters, e.g., $\mathbf{A},\mathbf{B}$. $\mathbf{I}$ denotes the $n\times n$ identity quaternion matrix. The operators ${(·)}^T$ and ${(·)}^H$ represent transpose and conjugate transpose, respectively. The MATLAB function command will be denoted by typewriter letters, e.g., $[V,D]=\mathtt{eig}\left(A\right)$.

2.2. Quaternions and Quaternion Matrices

Denote the set of quaternions as

$$\mathbb{Q}=\mathrm{span}\{1,\mathbf{i},\mathbf{j},\mathbf{k}\}\triangleq \left\{\mathbf{q}=q_1+q_2\mathbf{i}+q_3\mathbf{j}+q_4\mathbf{k}|q_1,q_2,q_3,q_4\in \mathbb{R}\right\},$$

(1)

where $\mathbf{i},\mathbf{j},\mathbf{k}$ are three imaginary units of quaternions, satisfying

$$\begin{array}{cc}\hfill & {\mathbf{i}}^2={\mathbf{j}}^2={\mathbf{k}}^2=\mathbf{ijk}=-1,\hfill \\ \hfill & \mathbf{ij}=-\mathbf{ji}=\mathbf{k},\mathbf{jk}=-\mathbf{kj}=\mathbf{i},\mathbf{ki}=-\mathbf{ik}=\mathbf{j}.\hfill \end{array}$$

The scalar (real) part of $\mathbf{q}$ is denoted by $\Re \left(\mathbf{q}\right)=q_1$. And the vector (imaginary) part of $\mathbf{q}$ is denoted by $\Im \left(\mathbf{q}\right)=q_2\mathbf{i}+q_3\mathbf{j}+q_4\mathbf{k}$. A quaternion is called imaginary when its real part is equal to zero. The multiplication of quaternions adheres to the distributive law but is noncommutative.

The zero element in $\mathbb{Q}$ is $\mathbf{0}=0+0\mathbf{i}+0\mathbf{j}+0\mathbf{k}$ and the unit element is $\mathbf{1}=1+0\mathbf{i}+0\mathbf{j}+0\mathbf{k}$. For any $\mathbf{q}=q_0+q_1\mathbf{i}+q_2\mathbf{j}+q_3\mathbf{k}\in \mathbb{Q}$, the conjugate of a quaternion is defined as

$${\mathbf{q}}^{*}=q_0-q_1\mathbf{i}-q_2\mathbf{j}-q_3\mathbf{k}.$$

The magnitude of $\mathbf{q}$ is $\left|\mathbf{q}\right|=\sqrt{{\mathbf{q}}^{*}\mathbf{q}}=\sqrt{q_0^2+q_1^2+q_2^2+q_3^2}$, it follows that the inverse of a nonzero quaternion $\mathbf{q}$ is given by ${\mathbf{q}}^{-1}={\mathbf{q}}^{*}/{\left|\mathbf{q}\right|}^2$.

Two quaternions $\mathbf{a}$ and $\mathbf{b}$ are said to be similar if there exists a nonzero quaternion $\mathbf{c}$ such that ${\mathbf{c}}^{-1}\mathbf{a}\mathbf{c}=\mathbf{b}$, this is written as $\mathbf{a}\sim \mathbf{b}$. Obviously, $\mathbf{a}$ and $\mathbf{b}$ are similar if and only if there is a unit quaternion $\mathbf{d}$ such that ${\mathbf{d}}^{-1}\mathbf{a}\mathbf{d}=\mathbf{b}$, and two similar quaternions have the same norm. It is routine to check that ∼ is an equivalence relation on the quaternions. We denote by $\left[\mathbf{a}\right]$ the equivalence class containing $\mathbf{a}$. If $\mathbf{q}=q_0+q_1\mathbf{i}+q_2\mathbf{j}+q_3\mathbf{k}\in \mathbb{Q}$, then $\mathbf{q}$ and $q_0+\sqrt{q_1^2+q_2^2+q_3^2}\mathbf{i}$ are similar, namely, $\mathbf{q}\in \left[q_0+\sqrt{q_1^2+q_2^2+q_3^2}\mathbf{i}\right]$.

Quaternions can also be expressed in polar form as $\mathbf{q}=\left|\mathbf{q}\right|(cos\theta +\widehat{\mathbf{q}}sin\theta )$, where $\widehat{\mathbf{q}}=\Im \left(\mathbf{q}\right)/|\Im \left(\mathbf{q}\right)|$ is a pure unit quaternion and $\theta =arccos(\Re (\mathbf{q})/|\mathbf{q}\left|\right)\in \mathbb{R}$ denotes the angle (or argument) of the quaternion. Next, we will introduce the quaternion rotation and involution operators.

Definition 1

(Quaternion rotation[19]). For any quaternion $\mathbf{q}$, the transformation

$${\mathbf{q}}^{\mu }\triangleq \mu \mathbf{q}{\mu }^{-1}$$

(2)

geometrically describes a three-dimensional rotation of the vector part of $\mathbf{q}$ by an angle $2\theta$ about the vector part of μ, where $\mu =\left|\mu \right|(cos\theta +\widehat{\mu }sin\theta )$ is any nonzero quaternion.

Specifically, if $\mu$ in (2) is an imaginary unit, then the quaternion rotation (2) reduces to quaternion involution [5], defined by

$$\begin{array}{cc}\hfill {\mathbf{q}}^{\mathbf{i}}& =-\mathbf{i}\mathbf{q}\mathbf{i}=q_0+\mathbf{i}{q}_1-\mathbf{j}{q}_2-\mathbf{k}{q}_3,\hfill \\ \hfill {\mathbf{q}}^{\mathbf{j}}& =-\mathbf{j}\mathbf{q}\mathbf{j}=q_0-\mathbf{i}{q}_1+\mathbf{j}{q}_2-\mathbf{k}{q}_3,\hfill \\ \hfill {\mathbf{q}}^{\mathbf{k}}& =-\mathbf{k}\mathbf{q}\mathbf{k}=q_0-\mathbf{i}{q}_1-\mathbf{j}{q}_2+\mathbf{k}{q}_3,\hfill \end{array}$$

where $\mathbf{q}=q_0+\mathbf{i}{q}_1+\mathbf{j}{q}_2+\mathbf{k}{q}_3\in \mathbb{Q}$. Below, we will list some properties of quaternion rotation, including

$${\left(\mathbf{p}\mathbf{q}\right)}^{\mu }={\mathbf{p}}^{\mu }{\mathbf{q}}^{\mu },\mathbf{p}\mathbf{q}={\mathbf{q}}^{\mathbf{p}}\mathbf{p}=\mathbf{q}{\mathbf{p}}^{\left({\mathbf{q}}^{*}\right)},\forall \mathbf{p},\mathbf{q}\in \mathbb{Q}$$

and

$${\mathbf{q}}^{\mu \nu }={\left({\mathbf{q}}^{\nu }\right)}^{\mu },{\mathbf{q}}^{\mu *}\triangleq {\left({\mathbf{q}}^{*}\right)}^{\mu }={\left({\mathbf{q}}^{\mu }\right)}^{*}\triangleq {\mathbf{q}}^{*\mu },\forall \nu ,\mu \in \mathbb{Q}.$$

Note that the representation in (1) can be extended to a general orthogonal basis $\left\{1,{\mathbf{i}}^{\mu },{\mathbf{j}}^{\mu },{\mathbf{k}}^{\mu }\right\}$, where the following properties hold [19]:

$${\mathbf{i}}^{\mu }{\mathbf{i}}^{\mu }={\mathbf{j}}^{\mu }{\mathbf{j}}^{\mu }={\mathbf{k}}^{\mu }{\mathbf{k}}^{\mu }={\mathbf{i}}^{\mu }{\mathbf{j}}^{\mu }{\mathbf{k}}^{\mu }=-1.$$

Denote the set of quaternion matrices as

$${\mathbb{Q}}^{m\times n}=\left\{\mathbf{A}=A_0+A_1\mathbf{i}+A_2\mathbf{j}+A_3\mathbf{k}|A_0,A_1,A_2,A_3\in {\mathbb{R}}^{m\times n}\right\}.$$

The conjugate transpose of $\mathbf{A}$ is ${\mathbf{A}}^{*}=A_0-A_1\mathbf{i}-A_2\mathbf{j}-A_3\mathbf{k}$. We say that a square quaternion matrix $\mathbf{A}\in {\mathbb{Q}}^{n\times n}$ is normal if ${\mathbf{A}}^{*}\mathbf{A}=\mathbf{A}{\mathbf{A}}^{*}$; Hermitian if ${\mathbf{A}}^{*}=\mathbf{A}$, i.e. $A_0=A_0^{\top }$ and $A_i=-A_i^{\top },i=1,2,3$; Unitary if ${\mathbf{A}}^{*}\mathbf{A}=\mathbf{I}$, where $\mathbf{I}$ is the identity matrix; Invertible (nonsingular) if there exists a matrix $\mathbf{B}\in {\mathbb{Q}}^{n\times n}$ such that $\mathbf{A}\mathbf{B}=\mathbf{B}\mathbf{A}=\mathbf{I}$. In this case, we denote ${\mathbf{A}}^{-1}=\mathbf{B}$. We have ${\left(\mathbf{A}\mathbf{B}\right)}^{-1}={\mathbf{B}}^{-1}{\mathbf{A}}^{-1}$ if $\mathbf{A}$ and $\mathbf{B}$ are invertible, and ${\left({\mathbf{A}}^{*}\right)}^{-1}={\left({\mathbf{A}}^{-1}\right)}^{*}$ if $\mathbf{A}$ is invertible.

2.3. GHR calculus

We now introduce the generalized HR derivatives which comprise both the product and chain rules, see [20,22] for more details.

Definition 2

(real-differentiability [20]). Let $\mathbf{q}=q_a+q_b\mathbf{i}+q_c\mathbf{j}+q_d\mathbf{k}\in \mathbb{Q}$, then a function $f\left(\mathbf{q}\right)=f_a\left(q_a,q_b,q_c,q_d\right)+f_b\left(q_a,q_b,q_c,q_d\right)\mathbf{i}+f_c\left(q_a,q_b,q_c,q_d\right)\mathbf{j}+f_d\left(q_a,q_b,q_c,q_d\right)\mathbf{k}$ is called real differentiable when $f_a\left(q_a,q_b,q_c,q_d\right),f_b\left(q_a,q_b,q_c,q_d\right),f_c\left(q_a,q_b,q_c,q_d\right)$ and $f_d\left(q_a,q_b,q_c,q_d\right)$ are differentiable with respect to the real variables $q_a,q_b,q_c$ and $q_d$, respectively.

Definition 3

(GHR derivatives [20]). If $f:\mathbb{Q}\to \mathbb{Q}$ is real differentiable, then the GHR derivatives of $f\left(\mathbf{q}\right)$ with respect to ${\mathbf{q}}^{\mu }$ and ${\mathbf{q}}^{\mu *}$$(0\ne \mu \in \mathbb{Q})$ are defined as

$${\displaystyle \frac{\partial f}{\partial {\mathbf{q}}^{\mu }}}={\displaystyle \frac{1}{4}}\left({\displaystyle \frac{\partial f}{\partial q_a}}-{\displaystyle \frac{\partial f}{\partial q_b}}{\mathbf{i}}^{\mu }-{\displaystyle \frac{\partial f}{\partial q_c}}{\mathbf{j}}^{\mu }-{\displaystyle \frac{\partial f}{\partial q_d}}{\mathbf{k}}^{\mu }\right),$$

and

$${\displaystyle \frac{\partial f}{\partial {\mathbf{q}}^{\mu *}}}={\displaystyle \frac{1}{4}}\left({\displaystyle \frac{\partial f}{\partial q_a}}+{\displaystyle \frac{\partial f}{\partial q_b}}{\mathbf{i}}^{\mu }+{\displaystyle \frac{\partial f}{\partial q_c}}{\mathbf{j}}^{\mu }+{\displaystyle \frac{\partial f}{\partial q_d}}{\mathbf{k}}^{\mu }\right),$$

where $\mathbf{q}=q_a+q_b\mathbf{i}+q_c\mathbf{j}+q_d\mathbf{k}\in \mathbb{Q}$, $q_a,q_b,q_c,q_d\in \mathbb{R}$, $\partial f/\partial q_a,\partial f/\partial q_b,\partial f/\partial q_c$ and $\partial f/\partial q_d$ are the partial derivatives of f with respect to $q_a,q_b,q_c$ and $q_d$, while the set $\left\{1,{\mathbf{i}}^{\mu },{\mathbf{j}}^{\mu },{\mathbf{k}}^{\mu }\right\}$ is an orthogonal basis of $\mathbb{Q}$.

Definition 4

(Quaternion gradient[20]). Let $f\left(\tilde{\mathbf{q}}\right):{\mathbb{Q}}^{n\times 1}\to \mathbb{Q}$ and $\tilde{\mathbf{q}}=\left({\mathbf{q}}_1,{\mathbf{q}}_2,\dots ,{\mathbf{q}}_n\right)\in {\mathbb{Q}}^{n\times 1}$, then the two quaternion gradients of f are defined as

$${\nabla }_{\tilde{\mathbf{q}}}f\triangleq {\left({\displaystyle \frac{\partial f}{\partial \tilde{\mathbf{q}}}}\right)}^T={\left({\displaystyle \frac{\partial f}{\partial {\mathbf{q}}_1}},\dots ,{\displaystyle \frac{\partial f}{\partial {\mathbf{q}}_n}}\right)}^T\in {\mathbb{Q}}^{n\times 1}$$

and

$${\nabla }_{{\tilde{\mathbf{q}}}^{*}}f\triangleq {\left({\displaystyle \frac{\partial f}{\partial {\tilde{\mathbf{q}}}^{*}}}\right)}^T={\left({\displaystyle \frac{\partial f}{\partial {\mathbf{q}}_1^{*}}},\dots ,{\displaystyle \frac{\partial f}{\partial {\mathbf{q}}_n^{*}}}\right)}^T\in {\mathbb{Q}}^{n\times 1}.$$

Based on the definitions of GHR provided above, we consider a simple quadratic function $f\left(\mathbf{x}\right)={\mathbf{x}}^H\mathbf{A}\mathbf{x}$, where $\mathbf{x}\in {\mathbb{Q}}^n$ and $\mathbf{A}\in {\mathbb{Q}}^{n\times n}$ is a quaternion Hermitian matrix, then the gradient of this function f is given by

$${\nabla }_{\mathbf{x}}f\left(\mathbf{x}\right)={\displaystyle \frac{1}{2}}{\left(\mathbf{A}\mathbf{x}\right)}^{*},{\nabla }_{{\mathbf{x}}^{*}}f\left(\mathbf{x}\right)={\displaystyle \frac{1}{2}}\mathbf{A}\mathbf{x},$$

in which ${\nabla }_{{\mathbf{x}}^{*}}f$ is the steepest ascent direction [22].

3. Quaternion Nesterov’s Accelerated Gradient (QNAG)

In this section, we will introduce the quaternion Nesterov accelerated projected gradient algorithm. To this end, we first review the definition related to the eigenvalue of quaternion Hermitian matrices.

Definition 5.

Let $\mathbf{A}\in {\mathbb{Q}}^{n\times n}$ be a quaternion Hermitian matrix. Then $\lambda \in \mathbb{R}$ is called an eigenvalue of $\mathbf{A}$ if there exists a nonzero vector $\mathbf{x}\in {\mathbb{Q}}^n$ such that

$$\mathbf{A}\mathbf{x}=\mathbf{x}\lambda .$$

Here $\mathbf{x}$ is called the eigenvector corresponding to the eigenvalue λ.

Due to the non-commutativity of quaternion multiplication, general quaternion square matrices have distinct left and right eigenvalues. However, for a quaternion Hermitian matrix $\mathbf{A}$, if $\mathbf{A}$ has a right eigenvalue $\lambda$ and its corresponding eigenvector $\mathbf{x}$, it is straightforward to show that

$${\mathbf{x}}^H\mathbf{A}\mathbf{x}={\mathbf{x}}^H\mathbf{x}\lambda ,$$

by dividing both sides of the above equation by ${\mathbf{x}}^H\mathbf{x}$, we obtain

$$\lambda ={\displaystyle \frac{{\mathbf{x}}^H\mathbf{A}\mathbf{x}}{{\mathbf{x}}^H\mathbf{x}}}={\displaystyle \frac{{\mathbf{x}}^H\mathbf{A}\mathbf{x}}{{\parallel \mathbf{x}\parallel }^2}}\in \mathbb{R},$$

(3)

where (3) represents the Rayleigh quotient on the quaternion skew field. Therefore, the eigenvalues of quaternion Hermitian matrices are all real numbers, and thus there is no distinction between left and right eigenvalues.

Our goal is to compute the principal eigenvalue of a given quaternion Hermitian matrix. By defining the objective function as $f\left(\mathbf{x}\right)={\mathbf{x}}^H\mathbf{A}\mathbf{x}$ and imposing the normalization constraint ${\parallel \mathbf{x}\parallel }^2=1$, we can equivalently transform the problem of finding the principal eigenvalue of a given quaternion Hermitian matrix into the following maximization optimization problem on the quaternion skew-field

$$\begin{array}{cc}\hfill \underset{\mathbf{x}\in {\mathbb{Q}}^n}{max}& f\left(\mathbf{x}\right)={\mathbf{x}}^H\mathbf{A}\mathbf{x}\hfill \\ \hfill \mathrm{s}.\mathrm{t}.& {\parallel \mathbf{x}\parallel }^2=1.\hfill \end{array}$$

(4)

The above problem (4) can be addressed using the quaternion gradient projection algorithm. Since the introduction of the Nesterov’s Accelerated Gradient (NAG) method [15], the incorporation of momentum has become a conventional approach to overcome the shortsighted issue in gradient algorithms [7,12]. To tackle the problem (4), we propose a quaternion Nesterov’s accelerated gradient projection algorithm (QNAG). Given an initial point ${\mathbf{x}}_0$, and set ${\mathbf{x}}_1={\mathbf{x}}_0$, the QNAG method repeats, for $t\ge 0$,

$$\begin{array}{cc}\hfill & {\mathbf{y}}_t={\mathbf{x}}_t+\beta \left({\mathbf{x}}_t-{\mathbf{x}}_{t-1}\right),\hfill \\ \hfill & {\mathbf{z}}_{k+1}={\mathbf{y}}_t+\alpha {\nabla }_{{\mathbf{y}}^{*}}f\left({\mathbf{y}}_t\right),\hfill \\ \hfill & {\mathbf{x}}_{t+1}={\mathbf{z}}_{t+1}/\parallel {\mathbf{z}}_{t+1}\parallel ,\hfill \end{array}$$

(5)

where $\alpha$ and $\beta$ are the step-size and momentum parameters, respectively. When the momentum parameter $\beta =0$, QNAG simplifies to standard gradient ascent (GA). When $\beta >0$ it is possible to achieve accelerated rates of convergence for certain combinations of $\alpha$ and $\beta$ in the deterministic setting. The framework of the proposed algorithm is detailed below.

Remark 1.

After obtaining the principal eigenvalue and its corresponding eigenvector of the quaternion Hermitian matrix, we can employ a deflation technique by updating $A=A-\lambda \mathbf{x}{\mathbf{x}}^H$ and continue to apply the QNAG algorithm. By repeating this process, all eigenvalues and their corresponding eigenvectors can be obtained.

4. Convergence Analysis of QNAG

In this section, we will theoretically prove the convergence properties of the quaternion Nesterov’s accelerated gradient projection algorithm. We begin our analysis with the following Lemma 1.

Lemma 1.

If $f:{\mathbb{Q}}^n\to \mathbb{R}$ is a real differentiable and gradient Lipschitz continuous function with constant $L>0$,

$$∥{\nabla }_{{\mathbf{x}}^{*}}f(\mathbf{x}+\Delta \mathbf{x})-{\nabla }_{{\mathbf{x}}^{*}}f\left(\mathbf{x}\right)∥\le L\parallel \Delta \mathbf{x}\parallel ,\forall \mathbf{x},\Delta \mathbf{x}\in {\mathbb{Q}}^n.$$

Algorithm 1: Quaternion Nesterov’s Accelerated Gradient (QNAG).

Input:  Given a quaternion Hermitian matrix $\mathbf{A}\in {\mathbb{Q}}^{n\times n}$, step size $0<\alpha \in \mathbb{R}$, momentum coefficient $0<\beta \in \mathbb{R}$ and the maximum iteration number $I_{max}$. Output:  The principle eigenvalue $\lambda$ and its corresponding eigenvector ${\mathbf{x}}_{t+1}$. 1: Initialize: a unit quaternion vector ${\mathbf{x}}_0\in {\mathbb{Q}}^n$. 2: ${\mathbf{x}}_1={\mathbf{x}}_0$. 3: for $t=1$ to $I_{max}$ do 4:    Momentum extrapolation: ${\mathbf{y}}_t={\mathbf{x}}_t+\beta ({\mathbf{x}}_t-{\mathbf{x}}_{t-1})$. 5:    Compute the gradient at the extrapolation point: ${\nabla }_{{\mathbf{y}}^{*}}f\left({\mathbf{y}}_t\right)={\displaystyle \frac{1}{2}}\mathbf{A}{\mathbf{y}}_t$. 6:    Gradient ascent: ${\mathbf{z}}_{k+1}={\mathbf{y}}_t+\alpha {\nabla }_{{\mathbf{y}}^{*}}f\left({\mathbf{y}}_t\right)$. 7:    Normalization: ${\mathbf{x}}_{t+1}={\mathbf{z}}_{t+1}/\parallel {\mathbf{z}}_{t+1}\parallel$. 8:    $\lambda ={\mathbf{x}}_{t+1}^H\mathbf{A}{\mathbf{x}}_{t+1}$. 9:    if $\parallel A{\mathbf{x}}_{t+1}-{\mathbf{x}}_{t+1}\lambda \parallel <{10}^{-10}$ then 10:      Break 11:    end if 12: end for

Then, f has the following quadratic upper bound

$$\left|f\left(\mathbf{y}\right)-f\left(\mathbf{x}\right)-4\Re \left\{{\nabla }_{{\mathbf{x}}^{*}}f{\left(\mathbf{x}\right)}^H(\mathbf{y}-\mathbf{x})\right\}\right|\le 2L{\parallel \mathbf{y}-\mathbf{x}\parallel }^2,\forall \mathbf{x},\mathbf{y}\in {\mathbb{Q}}^n.$$

(6)

Proof. 

For any $\mathbf{x},\mathbf{y}\in {\mathbb{Q}}^n$, let $\Delta \mathbf{x}=\mathbf{y}-\mathbf{x}$. Define the parameterized path $g\left(t\right)=f(\mathbf{x}+t\Delta \mathbf{x})$ for $t\in [0,1]$. Then the difference in the function can be expressed as

$$f\left(\mathbf{y}\right)-f\left(\mathbf{x}\right)={\int }_0^1{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}g\left(t\right)\mathrm{d}t.$$

Using the chain rule, the derivative $g^{\prime }\left(t\right)$ corresponds to the directional derivative of f along $\Delta \mathbf{x}$. In terms of the quaternion gradient [22], this is given by

$$g^{\prime }\left(t\right)=4\Re \left\{{\nabla }_{{\mathbf{x}}^{*}}f{(\mathbf{x}+t\Delta \mathbf{x})}^H\Delta \mathbf{x}\right\}.$$

Thus, the function difference becomes

$$f\left(\mathbf{y}\right)-f\left(\mathbf{x}\right)={\int }_0^{1}4\Re \left\{{\nabla }_{{\mathbf{x}}^{*}}f{(\mathbf{x}+t\Delta \mathbf{x})}^H\Delta \mathbf{x}\right\}\mathrm{d}t.$$

Subtracting the linear term $4\Re \left\{{\nabla }_{{\mathbf{x}}^{*}}f{\left(\mathbf{x}\right)}^H\Delta \mathbf{x}\right\}$ from both sides

$$f\left(\mathbf{y}\right)-f\left(\mathbf{x}\right)-4\Re \left\{{\nabla }_{{\mathbf{x}}^{*}}f{\left(\mathbf{x}\right)}^H\Delta \mathbf{x}\right\}={\int }_0^{1}4\Re \left\{{\left({\nabla }_{{\mathbf{x}}^{*}}f(\mathbf{x}+t\Delta \mathbf{x})-{\nabla }_{{\mathbf{x}}^{*}}f\left(\mathbf{x}\right)\right)}^H\Delta \mathbf{x}\right\}\mathrm{d}t,$$

then using the absolute value inequality and the Cauchy-Schwarz inequality, we have

$$\left|f\left(\mathbf{y}\right)-f\left(\mathbf{x}\right)-4\Re \left\{{\nabla }_{{\mathbf{x}}^{*}}f{\left(\mathbf{x}\right)}^H\Delta \mathbf{x}\right\}\right|\le 4{\int }_0^1∥{\nabla }_{{\mathbf{x}}^{*}}f(\mathbf{x}+t\Delta \mathbf{x})-{\nabla }_{{\mathbf{x}}^{*}}f\left(\mathbf{x}\right)∥·\parallel \Delta \mathbf{x}\parallel \mathrm{d}t.$$

(7)

By the gradient Lipschitz condition, we have

$$∥{\nabla }_{{\mathbf{x}}^{*}}f(\mathbf{x}+t\Delta \mathbf{x})-{\nabla }_{{\mathbf{x}}^{*}}f\left(\mathbf{x}\right)∥\le L\parallel t\Delta \mathbf{x}\parallel =Lt\parallel \Delta \mathbf{x}\parallel ,$$

substituting the above inequality into the integral (7), we get

$$\mathrm{Right}\text{-}\mathrm{hand}\mathrm{side}\le {4L\parallel \Delta \mathbf{x}\parallel }^2{\int }_0^{1}t\mathrm{d}t={4L\parallel \Delta \mathbf{x}\parallel }^2·{\displaystyle \frac{1}{2}}=2L{\parallel \Delta \mathbf{x}\parallel }^2.$$

Thus, for all $\mathbf{x},\mathbf{y}\in {\mathbb{Q}}^n$, we have

$$\left|f\left(\mathbf{y}\right)-f\left(\mathbf{x}\right)-4\Re \left\{{\nabla }_{{\mathbf{x}}^{*}}f{\left(\mathbf{x}\right)}^H(\mathbf{y}-\mathbf{x})\right\}\right|\le 2L{\parallel \mathbf{y}-\mathbf{x}\parallel }^2.$$

This completes the proof of the quadratic upper bound.

   □

Theorem 1.

Let $f:{\mathbb{Q}}^n\to \mathbb{R}$ is a real differentiable and gradient Lipschitz continuous function with constant $L>0$. If the step size $\alpha ={\displaystyle \frac{1}{2L}}$ and momentum parameter ${\beta }_t={\displaystyle \frac{t-1}{t+2}}$, then

$$f\left({\mathbf{x}}_t^{*}\right)-f\left({\mathbf{x}}_t\right)\le C/{(t+1)}^2$$

is true for any $t\ge 0$, where $f\left({\mathbf{x}}_t^{*}\right)={\lambda }_{max}$ and C is a constant related to the initial conditions.

Proof. 

From Lemma 1, for the update ${\mathbf{z}}_{t+1}={\mathbf{y}}_t+\alpha {\nabla }_{{\mathbf{y}}^{*}}f\left({\mathbf{y}}_t\right)$, we have

$$f\left({\mathbf{z}}_{t+1}\right)\ge f\left({\mathbf{y}}_t\right)+4\Re \left\{{\nabla }_{{\mathbf{y}}^{*}}f{\left({\mathbf{y}}_t\right)}^H\left({\mathbf{z}}_{t+1}-{\mathbf{y}}_t\right)\right\}-2L{∥{\mathbf{z}}_{t+1}-{\mathbf{y}}_t∥}^2$$

Substituting ${\mathbf{z}}_{t+1}-{\mathbf{y}}_t=\alpha {\nabla }_{{\mathbf{y}}^{*}}f\left({\mathbf{y}}_t\right)$, we derive

$$f\left({\mathbf{z}}_{t+1}\right)\ge f\left({\mathbf{y}}_t\right)+4\alpha {∥\nabla f\left({\mathbf{y}}_t\right)∥}^2-2L{\alpha }^2{∥\nabla f\left({\mathbf{y}}_t\right)∥}^2.$$

Choosing $\alpha ={\displaystyle \frac{1}{2L}}$, this simplifies to

$$f\left({\mathbf{z}}_{t+1}\right)\ge f\left({\mathbf{y}}_t\right)+{\displaystyle \frac{3}{2L}}{∥\nabla f\left({\mathbf{y}}_t\right)∥}^2.$$

(8)

After projecting ${\mathbf{z}}_{t+1}$ onto the unit sphere, ${\mathbf{x}}_{t+1}={\mathbf{z}}_{t+1}/∥{\mathbf{z}}_{t+1}∥$. Applying Lemma 1 again to ${\mathbf{x}}_{t+1}$ and ${\mathbf{z}}_{t+1}$, we get

$$f\left({\mathbf{x}}_{t+1}\right)\ge f\left({\mathbf{z}}_{t+1}\right)+4\Re \left\{{\nabla }_{{\mathbf{z}}_{t+1}^{*}}f{\left({\mathbf{z}}_{t+1}\right)}^H\left({\mathbf{x}}_{t+1}-{\mathbf{z}}_{t+1}\right)\right\}-2L{∥{\mathbf{x}}_{t+1}-{\mathbf{z}}_{t+1}∥}^2.$$

By the optimality condition of the projection, we have

$$\Re \left\{{\nabla }_{{\mathbf{z}}_{t+1}^{*}}f{\left({\mathbf{z}}_{t+1}\right)}^H\left({\mathbf{x}}_{t+1}-{\mathbf{z}}_{t+1}\right)\right\}\ge 0$$

which implies

$$f\left({\mathbf{x}}_{t+1}\right)\ge f\left({\mathbf{z}}_{t+1}\right)-2L{∥{\mathbf{x}}_{t+1}-{\mathbf{z}}_{t+1}∥}^2.$$

(9)

Next, we consider the bound of the projection error

$$\begin{array}{cc}\hfill \parallel {\mathbf{x}}_{t+1}-{\mathbf{z}}_{t+1}\parallel & =∥{\displaystyle \frac{{\mathbf{z}}_{t+1}}{\parallel {\mathbf{z}}_{t+1}\parallel }}-{\mathbf{z}}_{t+1}∥\hfill \\ \hfill & =∥{\displaystyle \frac{{\mathbf{z}}_{t+1}\left(1-\parallel {\mathbf{z}}_{t+1}\parallel \right)}{\parallel {\mathbf{z}}_{t+1}\parallel }}∥\hfill \\ \hfill & ={\displaystyle \frac{\left|1-\parallel {\mathbf{z}}_{t+1}\parallel \right|}{\parallel {\mathbf{z}}_{t+1}\parallel }}\parallel {\mathbf{z}}_{t+1}\parallel \hfill \\ \hfill & =\left|1-\parallel {\mathbf{z}}_{t+1}\parallel \right|.\hfill \end{array}$$

Combining ${\mathbf{z}}_{t+1}={\mathbf{y}}_t+\alpha {\nabla }_{{\mathbf{y}}^{*}}f\left({\mathbf{y}}_t\right)$ and normalizing such that $\parallel {\mathbf{y}}_t\parallel =1$, we obtain

$$\parallel {\mathbf{z}}_{t+1}\parallel \le \parallel {\mathbf{y}}_t\parallel +\alpha \parallel {\nabla }_{{\mathbf{y}}^{*}}f\left({\mathbf{y}}_t\right)\parallel ,$$

which implies that

$$\parallel {\mathbf{x}}_{t+1}-{\mathbf{z}}_{t+1}\parallel \le \alpha \parallel {\nabla }_{{\mathbf{y}}^{*}}f\left({\mathbf{y}}_t\right)\parallel .$$

Thus, the projection error in (9) is a higher-order term, we get

$$f\left({\mathbf{x}}_{t+1}\right)\ge f\left({\mathbf{z}}_{t+1}\right)-2L{\alpha }^2{∥\nabla f\left({\mathbf{y}}_t\right)∥}^2.$$

The formulas (8) and (9) give a lower bound

$$f\left({\mathbf{x}}_{t+1}\right)\ge f\left({\mathbf{y}}_t\right)+{\displaystyle \frac{3}{2L}}{∥\nabla f\left({\mathbf{y}}_t\right)∥}^2-{\displaystyle \frac{1}{2L}}{∥\nabla f\left({\mathbf{y}}_t\right)∥}^2=f\left({\mathbf{y}}_t\right)+{\displaystyle \frac{1}{L}}{∥\nabla f\left({\mathbf{y}}_t\right)∥}^2.$$

(10)

We then define the Lyapunov function as

$${\mathrm{\Phi }}_t={\theta }_t^2\left(f\left({\mathbf{x}}^{*}\right)-f\left({\mathbf{x}}_t\right)\right)+{\displaystyle \frac{1}{2}}{∥{\mathbf{v}}_t-{\mathbf{x}}^{*}∥}^2,$$

where the auxiliary sequence ${\theta }_t={\displaystyle \frac{t+1}{2}}$ satisfies ${\beta }_t={\displaystyle \frac{{\theta }_t-1}{{\theta }_{t+1}}}$, and ${\mathbf{v}}_t$ follows the update rule

$${\mathbf{v}}_{t+1}={\mathbf{v}}_t+{\theta }_t\alpha \nabla f\left({\mathbf{y}}_t\right).$$

Our goal is to prove that ${\mathrm{\Phi }}_{t+1}\le {\mathrm{\Phi }}_t$. Compute ${\mathrm{\Phi }}_{t+1}-{\mathrm{\Phi }}_t$, we obtain

$${\mathrm{\Phi }}_{t+1}-{\mathrm{\Phi }}_t={\theta }_{t+1}^2\left(f\left({\mathbf{x}}^{*}\right)-f\left({\mathbf{x}}_{t+1}\right)\right)-{\theta }_t^2\left(f\left({\mathbf{x}}^{*}\right)-f\left({\mathbf{x}}_t\right)\right)+{\displaystyle \frac{1}{2}}\left({∥{\mathbf{v}}_{t+1}-{\mathbf{x}}^{*}∥}^2-{∥{\mathbf{v}}_t-{\mathbf{x}}^{*}∥}^2\right).$$

(11)

For the term ${∥{\mathbf{v}}_{t+1}-{\mathbf{x}}^{*}∥}^2$, we have

$$\begin{array}{cc}\hfill {∥{\mathbf{v}}_{t+1}-{\mathbf{x}}^{*}∥}^2& ={∥{\mathbf{v}}_t-{\mathbf{x}}^{*}+{\theta }_t\alpha \nabla f\left({\mathbf{y}}_t\right)∥}^2\hfill \\ \hfill & ={∥{\mathbf{v}}_t-{\mathbf{x}}^{*}∥}^2+2{\theta }_t\alpha \Re \left\{\nabla f{\left({\mathbf{y}}_t\right)}^H\left({\mathbf{v}}_t-{\mathbf{x}}^{*}\right)\right\}+{\theta }_t^2{\alpha }^2{∥\nabla f\left({\mathbf{y}}_t\right)∥}^2.\hfill \end{array}$$

Thus we have

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{t+1}-{\mathrm{\Phi }}_t& =\left({\theta }_{t+1}^2-{\theta }_t^2\right)\left(f\left({\mathbf{x}}^{*}\right)-f\left({\mathbf{x}}_{t+1}\right)\right)+{\theta }_t^2\left(f\left({\mathbf{x}}_{t+1}\right)-f\left({\mathbf{x}}_t\right)\right)\hfill \\ \hfill & {\theta }_t\alpha \Re \left\{\nabla f{\left({\mathbf{y}}_t\right)}^H\left({\mathbf{v}}_t-{\mathbf{x}}^{*}\right)\right\}+{\displaystyle \frac{{\theta }_t^2+{\alpha }^2}{2}}{∥\nabla f\left({\mathbf{y}}_t\right)∥}^2.\hfill \end{array}$$

Suppose that

$$f\left({\mathbf{x}}^{*}\right)\le f\left({\mathbf{y}}_t\right)+\Re \left\{\nabla f{\left({\mathbf{y}}_t\right)}^H\left({\mathbf{x}}^{*}-{\mathbf{y}}_t\right)\right\}.$$

Substituting this into the Lyapunov difference and combining with (10), we arrive at

$$\begin{array}{cc}\hfill {\mathrm{\Phi }}_{t+1}-{\mathrm{\Phi }}_t& \le \left({\theta }_{t+1}^2-{\theta }_t^2-{\theta }_t^2\right)\left(f\left({\mathbf{x}}^{*}\right)-f\left({\mathbf{y}}_t\right)\right)\hfill \\ \hfill & +{\theta }_t\alpha \Re \left\{\nabla f{\left({\mathbf{y}}_t\right)}^H\left({\mathbf{v}}_t-{\mathbf{x}}^{*}\right)\right\}+{\displaystyle \frac{{\theta }_t^2{\alpha }^2}{2}}{∥\nabla f\left({\mathbf{y}}_t\right)∥}^2.\hfill \end{array}$$

Using ${\theta }_{t+1}={\theta }_t+{\displaystyle \frac{1}{2}}$ and $\alpha ={\displaystyle \frac{1}{2L}}$, we simplify

$${\theta }_{t+1}^2-{\theta }_t^2={\theta }_t+{\displaystyle \frac{1}{4}}\mathrm{and}{\displaystyle \frac{{\theta }_t^2}{L}}\ge {\theta }_t^2{\alpha }^{2}L.$$

After algebraic manipulation, we show that

$${\mathrm{\Phi }}_{t+1}\le {\mathrm{\Phi }}_t-{\displaystyle \frac{{\theta }_t^2}{4L}}{∥\nabla f\left(y_t\right)∥}^2\le {\mathrm{\Phi }}_t.$$

From ${\mathrm{\Phi }}_t\le {\mathrm{\Phi }}_0$ and ${\theta }_t={\displaystyle \frac{t+1}{2}}$, we have

$${\displaystyle \frac{{(t+1)}^2}{4}}\left(f\left({\mathbf{x}}^{*}\right)-f\left({\mathbf{x}}_t\right)\right)\le {\mathrm{\Phi }}_0={\displaystyle \frac{1}{4}}\left(f\left({\mathbf{x}}^{*}\right)-f\left(x_0\right)\right)+{\displaystyle \frac{1}{2}}{∥{\mathbf{v}}_0-{\mathbf{x}}^{*}∥}^2$$

Letting ${\mathbf{v}}_0={\mathbf{x}}_0$, we conclude

$$f\left({\mathbf{x}}^{*}\right)-f\left({\mathbf{x}}_t\right)\le {\displaystyle \frac{C}{{(t+1)}^2}},$$

where $C=f\left({\mathbf{x}}^{*}\right)-f\left({\mathbf{x}}_0\right)+2{∥{\mathbf{v}}_0-{\mathbf{x}}^{*}∥}^2$ is a constant related to the initial setup. Algorithm 1 achieves an $O(1/t^2)$ convergence rate. This completes the proof.

   □

5. Numerical Experiments

In this section, we provide numerical examples to demonstrate the feasibility and effectiveness of the quaternion Nesterov’s accelerated gradient projection algorithm for the eigenvalue problem of quaternion Hermitian matrices. In the specific implementation of Algorithm 1, we set the constant step size and momentum parameter to $\alpha =0.05$ and $\beta =0.9$, respectively.

All the experiments are performed under Windows 11 and MATLAB version 23.2.0.2365128 (R2023b) with an AMD Ryzen 7 5800H with Radeon Graphics CPU at 3.20 GHz and 16 GB of memory.

Example 1.

Given quaternion matrix $\mathbf{A}=A_0+A_1\mathbf{i}+A_2\mathbf{j}+A_3\mathbf{k}$ with

$$A_0=\left(\begin{array}{ccc}17.6331& -1.6420& -1.2730\\ -1.6420& 8.3929& -1.7952\\ -1.2730& -1.7952& 15.1089\end{array}\right),A_1=\left(\begin{array}{ccc}0& 1.2315& 1.5751\\ -1.2315& 0& 2.5700\\ -1.5751& -2.5700& 0\end{array}\right)$$

and

$$A_2=\left(\begin{array}{ccc}0& 0.6530& 3.2730\\ -0.6530& 0& 1.2301\\ -3.2730& -1.2301& 0\end{array}\right),A_3=\left(\begin{array}{ccc}0& -4.3909& -9.2817\\ 4.3909& 0& 1.9585\\ 9.2817& -1.9585& 0\end{array}\right).$$

In this experiment, we employ the quaternion Nesterov’s accelerated gradient (QNAG) method (Algorithm 1) to compute all eigenvalues of $\mathbf{A}$, which are ${\lambda }_1=27.0543,{\lambda }_2=12.4577,{\lambda }_3=1.6229$, with their corresponding eigenvectors being

$$\begin{array}{cc}\hfill & {\mathbf{x}}_1=\left(\begin{array}{c}0.4485+0.4371\mathbf{i}-0.1443\mathbf{j}-0.3829\mathbf{k}\\ -0.0365+0.0228\mathbf{i}+0.0147\mathbf{j}+0.1706\mathbf{k}\\ 0.3077+0.1006\mathbf{i}+0.2161\mathbf{j}+0.5076\mathbf{k}\end{array}\right),\hfill \\ \hfill & {\mathbf{x}}_2=\left(\begin{array}{c}0.1208-0.0357\mathbf{i}+0.3112\mathbf{j}-0.2267\mathbf{k}\\ 0.0627-0.7303\mathbf{i}-0.0928\mathbf{j}-0.0756\mathbf{k}\\ -0.2492+0.4640\mathbf{i}+0.0583\mathbf{j}-0.0588\mathbf{k}\end{array}\right),\hfill \\ \hfill & {\mathbf{x}}_3=\left(\begin{array}{c}-0.0011-0.1557\mathbf{i}+0.3809\mathbf{j}+0.3270\mathbf{k}\\ 0.5029+0.3354\mathbf{i}+0.0994\mathbf{j}+0.2046\mathbf{k}\\ 0.1149+0.4714\mathbf{i}+0.1735\mathbf{j}+0.2025\mathbf{k}\end{array}\right).\hfill \end{array}$$

Additionally, we obtain the following three residuals:

$$\begin{array}{cc}\hfill & \parallel \mathbf{A}{\mathbf{x}}_1-{\mathbf{x}}_1{\lambda }_1\parallel \le 6.7432e-12,\hfill \\ & \parallel \mathbf{A}{\mathbf{x}}_2-{\mathbf{x}}_2{\lambda }_2\parallel \le 5.1990e-12,\hfill \\ & \parallel \mathbf{A}{\mathbf{x}}_3-{\mathbf{x}}_3{\lambda }_3\parallel \le 4.3378e-12.\hfill \end{array}$$

It is evident that the residuals are controlled within an ideal range, demonstrating the feasibility and effectiveness of Algorithm 1 in computing the eigenvalues of quaternion Hermitian matrices.

Example 2.

In this experiment, we utilize MATLAB’s built-in functions to randomly generate three quaternion Hermitian matrices of different sizes and compare Algorithm 1 with the QPGA algorithm [3] and theeigfunction in the Quaternion Toolbox for MATLAB (QTFM)[17].

We first test the performance of Algorithm 1, QPGA, and the eig function in computing the principal eigenvalues of three different types and sizes of quaternion Hermitian matrices. The numerical experimental results are presented in Table 1, which includes three evaluation metrics: the number of iterations, residuals, and runtime. Superior results are highlighted in bold. It can be observed that Algorithm 1 outperforms the other algorithms in terms of the number of iterations, problem residuals, and runtime, demonstrating significant advantages when computing large-scale quaternion Hermitian matrices.

Subsequently, we plotted the variation curves of the objective function values for the first 50 iterations of Algorithm 1 and the QNGA algorithm, as shown in Figure 1. It is evident that our algorithm achieves a faster increase in the objective function, demonstrating higher efficiency in obtaining the maximum eigenvalue.

Figure 2 illustrates the residual variation curves generated by Algorithm 1 and the QNGA algorithm with respect to the number of iterations. Across the tested matrix dimensions, our algorithm consistently achieves higher accuracy and efficiency.

6. Conclusions

In this paper, leveraging the innovative generalized Hamilton-real (GHR) calculus, we have introduced a novel quaternion Nesterov’s accelerated projected gradient algorithm designed to compute the dominant eigenvalue and corresponding eigenvector of quaternion Hermitian matrices. The incorporation of momentum terms and look-ahead updates has enabled the algorithm to attain an accelerated convergence rate. Theoretical analysis has confirmed the convergence of the quaternion Nesterov’s accelerated projected gradient algorithm. Empirical results from numerical experiments indicate that the proposed method surpasses both the Quaternion Projected Gradient Ascent (QPGA) and conventional algebraic approaches in terms of both computational precision and efficiency in runtime.