Nesterov Accelerated Gradient Descent for Optimizing Fast Harmonic Mean Linear Discriminant Analysis

The availability of HDD in machine learning has surged as technology advances. HDD encompasses high-resolution images, databases, genes, and more. However, HDD poses major challenges for classification algorithms in terms of reliability and efficiency [1]. One core issue that arises is the "curse of dimensionality". This refers to the difficulty of selecting from a vast number of features when building a classifier [2]. As the number of dimensions or features grows, the amount of data needed to populate the space grows exponentially. This makes it incredibly difficult to find meaningful patterns with limited training data. As a result, classification performance suffers. Some algorithms are less susceptible as they rarely emphasize feature selection. Instead, they apply a classification algorithm to every available feature [3]. However, even if the prediction accuracy does not degrade, the computational cost can become prohibitive. This prevents the use of such exhaustive techniques with HDD in many real-world contexts [4]. To address this, dimensionality reduction methods were developed to extract the most useful data and features to facilitate HDD processing. Dimension reduction is often necessary for neural networks, object recognition, and other applications dealing with high-dimensional spaces [5].

Several learning approaches exist for handling HDD classification. These primarily aimed to reduce HDD problems via attribute selection and prediction [6]. Such techniques were also used to ensure diversity of classifications, using sample selection strategies like AdaBoost [7]. However, traditional techniques have some drawbacks: (i) many employ quick dimensionality reduction rather than expanding features, and (ii) refine specific objective functions rather than unsupervised methods like density-based clustering [8].

Several dimensionality reduction techniques have been developed, such as Principal Component Analysis (PCA), LDA, and Independent Component Analysis (ICA) [9-11]. These techniques aim to project data into a lower dimensional space, assuming linear relationships in the data. LDA, for example, maximizes class separability by maximizing between-class variance and minimizing within-class variance. However, it has limitations when dealing with high-dimensional complex data, as it may allow overlap between classes. In many real-world scenarios, the structure of HDD is complex and non-linear, such as in high-resolution image datasets or genomic data. Linear techniques like PCA and LDA fail to capture the intrinsic dimensionality and complex data relationships, resulting in poor representation of the data in the reduced space.Nonlinear dimensionality reduction methods, such as kernel PCA or manifold learning, address this issue by modeling nonlinear structures in data. However, they come with challenges such as higher computational costs. Determining the most suitable dimensionality reduction method requires assessing the inherent structure of the HDD, as linear techniques are insufficient for datasets with clear nonlinear relationships between dimensions. To address this shortcoming, algorithms like Harmonic mean-based LDA (HLDA) and HLDA-pairwise (HLDAp) were proposed [12]. These take into account harmonic average between-class gaps to improve class separation. However, HLDA and HLDAp have a computationally expensive initialization process involving matrix decomposition and inversion. This led to the development of Fast versions – FHLDA and FHLDAp – which use joint diagonalization and Taylor expansion [13] to reduce initialization iterations for discriminant formation. This avoids sweeping procedure, calculating all eigenvector components simultaneously. A first-order approximation of the inverse eigenvector matrix and full matrix were iteratively refined. However, class overlap can still cause misclassification. This was addressed by leveraging between-class scatter matrices to find the optimal discriminant vector. Still, the Stiefel manifold gradient has $O\left(p^2\right)$ complexity, where $p$ is the data dimension, because it requires many matrix multiplications.

Hence this paper proposes AOFHLDA and AOFHLDAp algorithms to reduce the complexity of the Stiefel manifold gradient schemes for FHLDA and FHLDAp. A Nesterov accelerated gradient descent scheme is introduced to handle orthogonality constraints and optimize functions on the Stiefel manifold. This achieves asymptotically optimal error rates for L-smooth convex and L-smooth, $\mu$-strongly convex functions. Thus, it ensures accelerated convergence once the solution is nearer to the local. After applying this scheme, joint diagonalization via Taylor expansion recovers the discriminant vector. This efficiently reduces the complexity of the Stiefel manifold gradient $(O(\sqrt{p}))$.

The proposed algorithms, AOFHLDA and AOFHLDAp, are well-suited for a wider range of HDD tasks compared to prior algorithms. They enable efficient processing of HDD such as hyperspectral images, genomics data, and multivariate time series.

In contrast, HLDA and FHLDA are mostly limited to lower-dimensional data like basic image classification tasks. By reducing the computational burden, the proposed algorithms can scale to much higher dimensionality HDD across various domains. This allows them to tackle more complex learning challenges such as identifying rare cell populations in single-cell genomic data or detecting anomalies in massive industrial sensor datasets. Additionally, these algorithms can resist outliers by estimating the pairwise within-class covariance matrix appropriately for the classification task. Their broad applicability to diverse HDD types is a key advantage over existing algorithms.

The remaining sections of this article are as follows:

Section 2 discusses similar efforts on dimensionality reduction in multiple systems. Section 3 describes the proposed algorithm, while Section 4 demonstrates its effectiveness. Section 5 summarizes this study and discusses future improvements.

1.png

Figure 1. Overall pipeline of the proposed study

The proposed scheme is to alter the function restart method and a below restart criterion is adopted, which forces an adequate reduction in the objective function.

$f\left(g_{t+1}\right)>f\left(g_t\right)-l_R \gamma_t\left\|\nabla f\left(y_t\right)\right\|^2$             (3)

In Eq. (3), $l_R$ denotes a variable, which is small constant, and $\gamma_t$ denotes the step size at step $t$.

This criterion is applied to the optimization issues on the Stiefel manifold, and its convergence for non-convex functions on Euclidean space is analyzed.

Consider $f$ is a differentiable, $L$-smooth function on $\mathbb{R}^n$, i.e., $\nabla f$ is Lipschitz with constant $L$. As well, consider that $f$ is bounded below.

Assume the iteration Eq. (4) with step size $\gamma_t$ selected to satisfy $l /{ }_I \leq \gamma_t \leq \gamma_{t-1}$ for some $l \leq 1$ and $f\left(g_{t+1}\right) \leq$ $f\left(y_t\right)-\left(\gamma_t / 2\right)\left\|\nabla f\left(y_t\right)\right\|_2^2$.

$\begin{gathered}g_0=y_0, g_{t+1}=y_t-\gamma_t \nabla f\left(y_t\right), y_{t+1}=g_{t+1}+ \\ \alpha_t\left(g_{t+1}-g_t\right)\end{gathered}$                     (4)

When this iteration is restarted whenever (3) holds, obtain

$\lim _{t \rightarrow \infty}\left\|\nabla f\left(g_t\right)\right\| \rightarrow 0$             (5)

3.2.2 Exploration and interpolation on the Stiefel manifold

Now, it is understood that how to obtain around knowing the robust convexity and smoothness variables and how to cope with non-convex functions in the process. Consider the difficulty of generalizing the momentum step of Eq. (4) to the manifold as:

$\mathcal{Y}_{t+1}=\mathcal{G}_{t+1}+\alpha_t\left(\mathcal{G}_{t+1}-\mathcal{G}_t\right)$            (6)

Consider the difficulty of effectively extrapolating and interpolating on the Stiefel manifold, i.e., given two points $\mathcal{G}, \mathcal{Y} \in S_{p, c}$ and $\alpha \in \mathbb{R}$, it is essential to determine points $(1-\alpha) \mathcal{G}+\alpha \mathcal{Y}$ on a curve via $\mathcal{G}$ and $\mathcal{Y}$. By assigning $\alpha \in$ $(0,1)$, this gives a way of averaging points on the manifold and by assigning $\alpha>1$ or $\alpha<0$, it can extrapolate as in Eq. (6).

A promising method would be to execute the extrapolation or interpolation in Euclidean space and then project back onto the Stiefel manifold. But this estimation process that comprises the orthogonal part from the joint diagonalization-based decomposition of the matrix is highly expensive compared to the proposed algorithm, particularly for large matrices. One could also substitute the estimation by a reorthogonalization process like Gram-Schmidt. But this is fairly improper when k (the number of vectors) is large and is also not as low-cost as this algorithm, which only comprises the matrix products and inversions.

The method for generalizing $(1-\alpha) \mathcal{G}+\alpha \mathcal{Y}$ is to solve for a $V \in\left(T_G S_{p, c}\right)^*$, which satisfies

$\mathcal{Y}=R\left(\mathcal{G}, \phi_q(V)\right), \mathcal{Y}=\mathcal{G}+V$               (7)

In Eq. (7), R is a retraction, which is predetermined. Then, extrapolate or average by assigning

$\begin{gathered}(1-\alpha) \mathcal{G}+\alpha \mathcal{Y}=R\left(\mathcal{G}, \phi_q(\alpha V)\right) \\ (1-\alpha) \mathcal{G}+\alpha \mathcal{Y}=\mathcal{G}+\alpha V\end{gathered}$                         (8)

Observe that the utilization of $\phi_q$ enables the algorithm to execute in the dual tangent space. The problem with this solving Eq. (7) for $V$, i.e., finding a $V$ such that $R\left(\mathcal{G}, \phi_q(V)\right)=\mathcal{Y}$ for some given $\mathcal{G}$ and $Y$. But when retraction is considered to be $R_1$, then this solves Eq. (9) for $V$.

$\left(I+\frac{1}{2}\left(V \mathcal{G}^T-\mathcal{G} V^T\right)\right) \mathcal{G}=\left(I-\frac{1}{2}\left(V \mathcal{G}^T-\mathcal{G} V^T\right)\right) \mathcal{Y}$                     (9)

Because $\mathcal{G}^T \mathcal{G}=\mathcal{Y}^T \mathcal{Y}=I$, one can simply verify that $V=$ $2 \mathcal{Y}\left(I+\mathcal{G}^T \mathcal{Y}\right)^{-1}$ solves Eq. (8). Here, observing $V \in \mathbb{R}^{p, c}$ as an element of the dual tangent space via the Frobenius inner product. Also, verify that when $V$ is substituted by $V^{\prime}=V+$ $\mathcal{G} S$ for any symmetric $c \times c$ matrix $S$, then $V \mathcal{G}^T-\mathcal{G} V^T=$ $V^{\prime} \mathcal{G}^T-\mathcal{G}\left(\mathcal{G}^{\prime}\right)^T$. This defines that $V^{\prime}$ also satisfies Eq. (9). Particularly, $V$ is substituted by its orthogonal projection onto the dual tangent space $\left(T_G\right)^*$, which then provides the target vector.

One possible problem with this method is the probability that the matrix $\left(I+\mathcal{G}^T \mathcal{Y}\right)$ could be singular or ill-conditioned. To solve this problem, the matrix $\left(I+\mathcal{G}^T \mathcal{Y}\right)$ is wellconditioned providing $\mathcal{G}$ and $\mathcal{Y}$ are not too far apart on $S_{p, c}$. Consider $\mathcal{G}, \mathcal{Y} \in S_{p, c}$. The geodesic distance between $\mathcal{G}$ and $\mathcal{Y}$ w.r.t. the quotient metric as:

$\begin{gathered}d_{S_{p, c}^Q}(\mathcal{G}, \mathcal{Y})^2=\inf _{C(t):[0,1] \rightarrow S_{p, c}} \int_0^1 \operatorname{Tr}\left(C^{\prime}(t)(I-\right.  \left.\left.\frac{1}{2} C(t) C(t)^T\right) C^{\prime}(t)^T\right) d t\end{gathered}$                        (10)

In Eq. (10), the infimum is taken over each path $C(t)$ that connect $\mathcal{G}$ and $\mathcal{Y}$, i.e., for which $C(0)=\mathcal{G}$ and $C(1)=\mathcal{Y}$. Likewise, the geodesic distance is written w.r.t. the embedding metric as:

$d_{S_{p, c}^E}(\mathcal{G}, \mathcal{Y})^2=\inf _{C(t):[0,1] \rightarrow S_{p, c}} \int_0^1 \operatorname{Tr}\left(C^{\prime}(t) C^{\prime}(t)^T\right) d t$                        (11)

In Eq. (11), the infimum is taken over the same set where $C(0)=\mathcal{G}$ and $C(1)=\mathcal{Y}$.

Thus, this method provides computationally efficient algorithm for averaging and extrapolating on the Stiefel manifold.

3.2.3 Gradient restart method

Remember that the gradient restart method reinitiates iteration (4) whenever $\nabla f\left(y_{t-1}\right) \cdot\left(\mathcal{G}_t-\mathcal{G}_{t-1}\right)>0$. Observing that $\mathcal{G}_t=\mathcal{Y}_{t-1}-\gamma_{t-1} \nabla f\left(\mathcal{Y}_{t-1}\right)$ and this criterion is rewritten as:

$-\gamma_{t-1}\left\|\nabla f\left(\mathcal{Y}_{t-1}\right)\right\|_2^2+\nabla f\left(\mathcal{Y}_{t-1}\right) \cdot\left(\mathcal{Y}_{t-1}-\mathcal{G}_{t-1}\right)>0$              (12)

It is observed that on the manifold $\left\|\nabla f\left(\mathcal{Y}_{t-1}\right)\right\|_2^2$ must become $\left\|\nabla f\left(\mathcal{Y}_{t-1}\right)\right\|_{q^*}^2$. Here, seeing $\nabla f\left(\mathcal{Y}_{t-1}\right)$ as an element of the dual tangent space. The difficult part is generalizing $+\nabla f\left(\mathcal{Y}_{t-1}\right) \cdot\left(\mathcal{Y}_{t-1}-\mathcal{G}_{t-1}\right)$. To solve for $V \in\left(S_{y_{t-1}}\right)^*$, the method is presented such that

$\mathcal{G}_{t-1}=R\left(\mathcal{Y}_{t-1}, \phi_q(V)\right)$             (13)

Then, this element $V$ acts as $\mathcal{G}_{t-1}-\mathcal{Y}_{t-1}$  and the correspondent of the gradient restart criterion becomes

$-\gamma_{t-1}\left\|\nabla f\left(\mathcal{Y}_{t-1}\right)\right\|_{q^*}^2-\left\langle\nabla f\left(\mathcal{Y}_{t-1}\right), V\right\rangle_{q^*}>0$                   (14)

Eq. (13) can be effectively resolved for V when the retraction is used as $R_1$.

3.2.4 Accelerated gradient descent on the Stiefel manifold

The above-studied concepts are combined to design an accelerated gradient descent scheme with the gradient restart on the Stiefel manifold as in Algorithm 1.

Algorithm 1 Accelerated gradient descent with gradient restart method

Input: Smooth function $f$, tolerance $\epsilon$, initial step size $\gamma_0$, variables needed for line search $\lambda_d, c_L$

Output: A point $\mathcal{G}_t$ such that $\left\|\nabla f\left(\mathcal{G}_t\right)\right\|_{q^*}<\epsilon$

1. Initialize
2. $\mathcal{G}_0 \leftarrow$ initial point;
3. $\mathcal{Y}_0 \leftarrow \mathcal{G}_0, t \leftarrow 0, c \leftarrow 0$;
4. while $\left(\left\|\nabla f\left(\mathcal{G}_t\right)\right\|_{q^*} \geq \epsilon\right)$
5. $\mathcal{G}_{t+1} \leftarrow R_1\left(\mathcal{Y}_t, \phi_q\left(-\gamma_t \nabla f\left(\mathcal{Y}_t\right)\right)\right)$, execute a line search to guarantee the Armijo criterion $f\left(\mathcal{G}_{t+1}\right) \leq f\left(\mathcal{Y}_t\right)-\frac{1}{2}\left\|\nabla f\left(\mathcal{Y}_t\right)\right\|_{q^*}^2$ is satisfied;
6. $\boldsymbol{w h i l e}\left(f\left(\mathcal{G}_{t+1}\right) \leq f\left(\mathcal{Y}_t\right)-l_L \gamma_t\left\|\nabla f\left(\mathcal{Y}_t\right)\right\|_{q^*}^2\right)$
7. $\gamma_t \leftarrow \lambda_d \gamma_t$
8. $\mathcal{G}_{t+1} \leftarrow R_1\left(\mathcal{Y}_t, \phi_q\left(-\gamma_t \nabla f\left(\mathcal{Y}_t\right)\right)\right)$;
9. end while
10. while $\left(f\left(\mathcal{G}_{t+1}\right)>f\left(\mathcal{Y}_t\right)-\frac{1}{2} \gamma_t\left\|\nabla f\left(\mathcal{Y}_t\right)\right\|_{q^*}^2\right)$
11. $\gamma_t \leftarrow{ }^{\gamma_t} / \lambda_d$;
12. $\mathcal{G}_{t+1} \leftarrow R_1\left(\mathcal{Y}_t, \phi_q\left(-\gamma_t \nabla f\left(\mathcal{Y}_t\right)\right)\right)$;
13. end while
14. $W_t \leftarrow 2 \mathcal{G}_t\left(I+\mathcal{G}_t^T \mathcal{Y}_t\right)^{-1}$;
15. $W_t \leftarrow W_t-\frac{1}{2} \mathcal{G}_t\left(W_t^T \mathcal{G}_t+\mathcal{G}_t^T W_t\right)$;
16. $\boldsymbol{i} \boldsymbol{f}\left(\left(\nabla f\left(\mathcal{Y}_t\right), W_t\right)_{q^*}<-\gamma_t\left\|\nabla f\left(\mathcal{Y}_t\right)\right\|_{q^*}^2\right)$
17. $\mathcal{G}_{t+1} \leftarrow \mathcal{G}_t, \mathcal{Y}_t \leftarrow \mathcal{G}_{t+1}, c \leftarrow 0$;
18. else
19. $V_t \leftarrow 2 \mathcal{G}_{t+1}\left(I+\mathcal{G}_{t+1}^T \mathcal{G}_t\right)^{-1}$;
20. $V_t \leftarrow V_t-\frac{1}{2} \mathcal{G}_t\left(V_t^T \mathcal{G}_t+\mathcal{G}_t^T V_t\right)$;
21. $\mathcal{Y}_{t+1} \leftarrow R_1\left(\mathcal{G}_t,\left(1+\frac{c}{c+3}\right) \phi_q\left(V_t\right)\right)$;
22. $c \leftarrow c+1$;
23. end if
24. $t \leftarrow t+1$;
25. $\gamma_{t+1}=\gamma_t$;
26. end while
27. End

After that, $\mathcal{G}$ is retrieved to the manifold based on the joint diagonalization via Taylor expansion with less computational complexity [13]. Accordingly, the proposed accelerated optimization can achieve faster convergence and reduce the complexity of Stiefel gradient manifolds effectively.

3.3 Accelerated optimization with FHLDAp

For FHLDAp, a pairwise within-class scatter matrix $\left(\mathcal{W}_{k l}\right)$ is also determined along with $\mathcal{B}_{k \ell}$. The objective function for FHLDAp is rewritten as:

$\min _G \mathcal{J}_2(\mathcal{G})=\sum_{k<l} n_k n_{\ell} \frac{\operatorname{Tr}\left(\mathcal{G}^T \mathcal{W}_{k \ell} \mathcal{G}\right)}{\operatorname{Tr}\left(\mathcal{G}^T \mathcal{B}_{k \ell \mathcal{G}}\right)}$, s.t. $\mathcal{G}^T \mathcal{G}=I$                 (15)

Similar to the FHLDA, the accelerated gradient descent on the Stiefel manifold is applied to solve (15). Algorithm 2 presents the overall steps in solving Eq. (1) or Eq. (15) iteratively.

Algorithm 2 Accelerated optimization on FHLDA or FHLDAp algorithm

Input: $X \in \mathbb{R}^{p \times n}$, class matrix $\mathcal{Y} \in \mathbb{R}^{n \times K}$ and $c$

putput: Projection matrix $\mathcal{G} \in \mathbb{R}^{p \times c}$

1. Initialize $\mathcal{G}$ using standard LDA if $c \leq K-1$, or usino trace ratio $I D A$ if $c>K-1$;
2. Calculate Calculate $\mathcal{S}_w$ and $\mathcal{B}_{k l}$;
3. Calculate the discriminant value of all classes denoted as $\left(\alpha_k^{\prime}, \beta_k^{\prime}\right), i=1, \ldots, K$;
4. while $\left(\begin{array}{c}\text { objective function Eq. (1) or } \\ \text { Eq.(16) not converge }\end{array}\right)$
5. Determine the accelerated gradient descent on Stiefel manifold based on Algorithm 1;
6. Retrieve $\mathcal{G}$ to the manifold utilizing the joint diagonalization;
7. end while

The algorithm initializes the projection matrix $\mathcal{G}$ using standard LDA or trace ratio LDA based on the subspace dimension. It then performs an iterative optimization process using accelerated gradient descent on the Stiefel manifold (Algorithm 1). The optimization is guided by the specified objective function Eq. (1) or Eq. (16), and the resulting $\mathcal{G}$ is retrieved to the manifold using joint diagonalization. The algorithm continues iterating until convergence of the objective function is achieved, aiming to accelerate the optimization process for FHLDA or FHLDAp and potentially enhance convergence speed and overall performance.

Figure 2 illustrates the convergence behavior of various algorithms on the Steifel manifold. It compares the number of iterations to the condition number on a logarithmic scale. The accelerated approach shows superior convergence, with the number of iterations scaling slightly better than the square root of the condition number. Therefore, the accelerated gradient descent with gradient restart method outperforms others, even for small condition numbers.

2.png

Figure 2. Convergence behavior of different algorithms on Steifel manifold

The subsequent section will focus on the empirical validation of the proposed AOFHLDA and AOFHLDAp algorithms. Extensive experiments are conducted on diverse benchmark datasets including facial images, object images, and video datasets. The performance of AOFHLDA and AOFHLDAp is compared against existing state-of-the-art LDA variants. Evaluation metrics and computational efficiency are reported to fully assess the effectiveness and practical utility of the proposed algorithms. Details including dataset statistics, experimental setup, and parameter tuning, will be provided to ensure transparency and reproducibility of the results.