A breakdown-free block conjugate gradient method for large-scale discriminant analysis

1.   Introduction

Dimensionality reduction is a common technique for extracting effective information from high-dimensional data, which can validly reduce the computational complexity and storage demands. Therefore, dimensionality reduction has become a research hotspot in recent years ^[1,2,3]. Linear discriminant analysis (LDA) ^[4,5,6,7,8] is a well-known method for dimensionality reduction, which has been widely used in many applications such as face recognition ^[1,9,10], machine learning ^[11,12], image classification ^[13,14], image reconstruction ^[15], and information retrieval ^[2,16].

Suppose $n$ data points are divided into $c$ classes $X = [{\bf x}_{1}, {\bf x}_{2}, \ldots, {\bf x}_{n}] = [X_1, \ldots, X_c]\in\mathbb{R}^{d\times n}$, where $X_j\in\mathbb{R}^{d\times n_j}$ is the $j$-th class with $n_{j}$ being the number of samples, and $\sum_{j = 1}^{c} n_j = n$. Denote by ${\bf c}_{j}$ the centroid of class $X_j$, and by ${\bf c}$ the global centroid. The within-class scatter, between-class scatter and total scatter matrices are defined as

respectively. LDA looks for a classification that projects high-dimensional data onto a low-dimensional space and achieves the maximum class separability of data ^[17,18]. In order to make LDA more suitable for dimensionality reduction, various extensions of LDA are proposed ^[19,20]. Among those extensions, one popular LDA criterion is

As $S_t$ is a symmetric positive semidefinite matrix, the regularization parameter is often used to overcome singularity problems ^[19]. $S_t$ in (1.1) is replaced by $S = S_t+\alpha I_d$, which in results the regularized linear discriminant analysis (RLDA) problem. The literature ^[21] proposes a branching and binding method for solving problems in the case of reducing to one-dimensional space. The literature ^[22] proposes a new dual parameter regularization BLDA (RBLDA) for MTS data classification and develops an efficient model selection algorithm.

Generally, an LDA problem can be solved equivalently by a generalized eigenvalue problem ^[4,17,18]

However, straightforward implementation of generalized eigenvalue decomposition is very time-consuming and prohibitive for high-dimensional data ^[18,23]. To cure this drawback, researchers use QR decomposition, singular value decomposition, the Krylov subspace method and other techniques to transform the large-scale generalized eigenvalue decomposition problem into a small size matrix decomposition problem, resulting in a series of methods for high-dimensional data, such as QRLDA ^[24], GSVDLDA ^[25], RRDA ^[18], and so on ^[3,25,26,27]. Recently, random sampling ^[17] and randomized SVD (RSVD) ^[26,28] techniques have been used to accelerate LDA algorithms, and it shows that RSVD is very efficient for high-dimensional and large-scale dense data. Among these algorithms, RRDA ^[18] solved the generalized eigenvalue problem using a gradient-like method using the Rayleigh-Ritz framework, showng its superiority to some popular LDA algorithms both on the optimal target value and classification accuracy. However, RRDA is not concise in the implementation framework, hence we are committed to studying its convenient framework in this paper.

In addition to solving LDA by eigenvalue problems, LDA can also be solved by multivariate linear regression with a specific class indicator matrix ^[5,8,29]. In ^[8], the author showed that the solution of multivariate linear regression $\min\limits_{M}\|\frac{1}{\sqrt{n}}\bar{X}^TM-Y\|^2_F$ can solve the corresponding LDA problem (1.1), where

is a class indicator matrix and $\bar{X} = X(I_n-\frac{1}{n}{\bf 1}_n{\bf 1}^T_n)$ is the centralized data matrix. Similarly, ^[27] established an equivalence between RLDA and a ridge regression, and ^[29] investigated the relationship between LDA and the minimum squared error. Chen et al. ^[30] proposed a bidirectional linear discriminant analysis method for image data based on low rank approximation. SRDA ^[23] cast LDA into a regression framework using spectral graph analysis is a representative algorithm for solving high-dimensional and large-scale sparse data. LDADL ^[20] showed an equivalence solution space between LDA and a multivariate linear regression, and proposed an incremental method for large-scale data. Most least square (LS) problems for LDA are based on the data matrix, not the data itself, which will increase computational complexity. In this paper, we consider establishing the relationship between LDA and multivariate linear regression through the original matrix itself.

The structure of this paper is as follows. In Section 2, we first introduce the RRDA method and analyze its advantages and disadvantages, then put forward some settlements for the defects and optimize RRDA. In Section 3, we establish the equivalence of LDA and a new linear symmetric and positive equations system, which contains only data matrix and class indicator matrix. In Section 4, we conduct some numerical experiments by a collection of large-scale and high-dimensional databases from real face images, speech corpus, and video. Numerical results show the effectiveness of our proposed method, and some concluding remarks are given in Section 5.

Here are some notations in our paper. ${\bf 1}_i$ stands for the vector of all ones with dimension $i$, $I_i$ is the identity matrix whose dimension is $i$, ${\bf 0}$ means a zero matrix or vector. $A^T$, $A^{-1}$ and $A^{\dagger}$ are the transposition, inverse, and Moore-Penrose inverse of $A$, respectively. $\|A\|_2$, $\|A\|_F$ represent the 2-norm and Frobenius norm of $A$. $span\{A\}$ consists of the space spanned by the columns of $A$, $rank(A)$ stands for the rank of $A$, $orth(A)$ is an orthonormal basis for the range of $A$, and $tr(A)$ indicates the trace of $A$.

2.   The RRDA method and its improved version

In this section, we briefly review the (RRDA) method ^[18] and optimize it.

2.1.   RRDA

RRDA is a conjugate-gradient-like method based on the Rayleigh-Ritz framework for the generalized eigenvalue problem $S_{b}{\bf w} = \lambda S{\bf w}$ of RLDA, where

For the sake of the Rayleigh-Ritz framework, RRDA rewrote $S_t$ and $S_b$ as ^[18]

where

Therefore, RRDA solves the following problem for dimensionality reduction:

RRDA designs efficient subspace expansion and extraction strategy on the Rayleigh-Ritz framework for generalized eigenvalue problems (2.4). In the expansion phase, it chooses a subspace basis $P_k$ satisfying

such that $V_k = [V_{k-1}, P_k]$ in iteration $k$, where $M_{k-1}$ from the prior extraction phase with the initial $M_0 = {\bf 0}$. $P_k$ in the above formula can be computed by the following theorem ^[18].

Theorem 1. ^[18] Define $\{P_i\}$ and $\{R_i\}$ as

then

In the extraction phase, it computes the eigenpair $(\lambda_i, \eta_i)$ of $(V_k^THH^TV_k, V_k^TSV_k)$ and computes the Ritz vectors ${\bf w}_i = V_k\eta_i$ until convergence. Let $W_k = [{\bf w}_1, {\bf w}_2, \ldots, {\bf w}_c]$, where $c$ is the reduced dimension, then $W_k$ is the approximation of the eigenvectors of (2.4) associated with the first $c$ eigenvalues.

As $S$ is the positive definite matrix, computing the eigenpairs of $(V_k^TSV_k, V_k^THH^TV_k)$ can transform into computing the eigenpairs of $\lambda \eta = (V_k^TSV_k)^{\dagger}(V_k^TH)(V_k^TH)^T\eta$ in the following problem ^[31]. It is known that if $(\lambda, \eta)$ is the nonzero eigenpair of $B^TAB$, then $ABB^T$ has the eigenvector $AB\eta$ with the same eigenvalue. Hence, let $T_k = H^TM_k$, $M_k = V_k(V_k^TSV_k)^{\dagger}V_k^TH$, and $U_k$ be the eigenvectors with the positive eigenvalues of $T_k$. Then, the Ritz pair of $(HH^T, S)$ can be recovered by $V_k\cdot (V_k^TSV_k)^{\dagger}V_k^TH\cdot U_k = M_kU_k$. On the other hand, according to Theorem 1, $M_k = V_k(V_k^TSV_k)^{\dagger}V_k^TH$ can be computed equivalently by

Thus, the specific algorithm of RRDA is shown in Algorithm 1.

2.2.   An improved algorithm of RRDA

Throughout the implementation of RRDA, there are some shortcomings and suggestions here.

First, the large-scale matrix $S = \bar{X}\bar{X}^T+\alpha I_d$ is formed in the initial step, which results in a waste of time and memory. Notice that $S$ is only used in $Q_k = SP_k$. We suggest computing $Q_k$ by its equivalent form $Q_k = \bar{X}(\bar{X}^TP_k)+\alpha P_k$ without forming $S$, which will save a lot of computing memory.

Second, $M_k$ can be updated by $M_k = M_k+P_kE_kL_k$ as (2.6), and (2.6) holds because of the orthogonality conditions of (2.5). Therefore, the premise for the effective implementation of RRDA is the orthogonality of $P_k$ and $R_k$, while their orthogonality may lose gradually during its running due to computational error. Hence, we need to re-orthogonalize $P_k$ in Step 3.1, i.e., set

Fortunately, ^[32] has proved that $R_k^TP_{k-1} = {\bf 0}$, $R_k^TR_{j} = {\bf 0}$ would still be established both on theoretical and computational results with $P_k = orth(P_k)$. Further, $P_k^TQ_k$ is symmetric positive definite matrix under the orthogonality of $P_k$, thus one can compute $E_k = (P_k^TQ_k)^{-1}$ instead of computing the Moore-Penrose inverse $E_k = (P_k^TQ_k)^{\dagger}$ in Step 3.2.1.

Third, we suggest placing the term $T_k = H^TM_k$ in Steps 3.2.2 and the 3.2.3 outside the loop, as the update of $T_k$ and $W_k$ is only related to $M_k$. This is our essential improvement proposal for RRDA. In this way, we can reduce the CPU time of RRDA since it only computes the matrix-matrix products $T_k = H^TM_k$, $W_k = M_kU_k$ and the eigenvalue decomposition of $T_k$ once, while RRDA needs to compute $k$ times. In essence, ignoring the calculation of $T_k$ and $W_k$, the expansion and extraction of RRDA (Algorithm 1) is the process of solving positive definite equations

by the block conjugate gradient method, where $P_k$ is the descent direction, $R_k$ is the residual term, and $M_k$ is the iterative solution. Therefore, the stopping criterion $\|W_k-W_{k-1}\|_F < \epsilon$ in Step 3.3 is improper. In fact, due to the gradual loss of orthogonality of $P_k$, the dimensions of $W_k$ and $W_{k-1}$ may be different during running, resulting in the failure of the RRDA. Thus, we suggest using $\|R_k\|_F < \epsilon$ instead of $\|W_k-W_{k-1}\|_F < \epsilon$ as the stopping criterion.

Remark 1. The solving process of Algorithms 2 and 1 is similar, so convergence analysis can refer to the convergence analysis of Algorithm 1. Let us briefly discuss the computational complexities of this algorithm. The main computational complexity is in Step 3, which costs about $\mathcal{O}(dnc+d^2c+)$ flops. In Case 2, the main computational complexity is in Step 10, which costs about $\mathcal{O}(dnc+c^3+c^2d+d^2c)$ flops. Generally speaking, the class number $c$ is far less than the characteristic dimension $d$ and the sample number $n$, thus the total costs are about $\mathcal{O}(dnc+c^2d+d^2c)$ flops.

In conclusion, we summarize the improved algorithm for RRDA; see Algorithm 2. The essential difference between RRDA and IRRDA is the handing of $T_k$ and $W_k$. RRDA updates $T_k$ and $W_k$ by $M_k$ iteratively in the inner loop until $W_k$ converges, while IRRDA calculates $T_k$ and $W_k$ outside the loop. Throughout RRDA, $T_k$ converges as long as $M_k$ converges because of the continuity of matrix elements. Thus, the iterations of RRDA and IRRDA may be close to each other when they converge. However, RRDA calculates eigenvalue decompositions more than IRRDA, so IRRDA has better computational efficiency.

According to the implementation of IRRDA, IRRDA essentially computes ${\bf w}$ and satisfies the following problems:

Recall that $S$ is symmetric positive definite. We have the following result relating RRDA problem (2.4) and IRRDA problem (2.8).

Theorem 2. Let $M$ be the solution of $SM = H$, then RRDA problem (2.4) and IRRDA problem (2.8) have the same solution.

Proof. As $SM = H$ and $S$ is symmetric positive definite, then $M = S^{-1}H$ and RRDA problem (2.4) can be rewritten as the following eigenvalue problem:

Now, we prove that the solutions of RRDA problem (2.9) and IRRDA problem (2.8) are the same.

First, we prove that ${\bf w}$ in (2.8) solves RRDA problem (2.4). Left multiplying $M$ from both sides of $H^TM{\bf u} = \lambda{\bf u}$, $\lambda > 0$ yields $MH^T(M{\bf u}) = \lambda(M{\bf u})$, which implies that if ${\bf u}$ is an eigenvector of $H^TM$ corresponding to $\lambda$, then ${\bf w} = M{\bf u}$ is an eigenvector of $MH^T$ corresponding to $\lambda$.

Second, we prove that the solution of RRDA satisfies the IRRDA problem (2.8). Left multiplying $H^T$ from both sides of (2.9) yields

which indicates that if ${\bf w}$ is an eigenvector of $MH^T$ corresponding to $\lambda$, then $H^T{\bf w}$ is an eigenvector of $H^TM$ corresponding to $\lambda$. In other words, the eigenvectors of $H^TM$ corresponding to nonnegative eigenvalues can be computed by ${\bf u} = H^T{\bf w}$. Left multiplying $M$ from both sides of (2.10) yields

which indicates that if ${\bf u} = H^T{\bf w}$ is an eigenvector of $H^TM$ corresponding to $\lambda$, then $M{\bf u}$ is an eigenvector of $MH^T$ corresponding to $\lambda$. Note ${\bf w}$ is also the eigenvector of $MH^T$ corresponding to $\lambda$, thus ${\bf w} = M{\bf u}$. □

Theorem 2 depicts that the solution spaces of RRDA and IRRDA are equivalent. Fortunately, algorithm IRRDA has less computational complexity, which shows the effectiveness of the improved algorithm IRRDA. Moreover, let $U$ be the eigenvectors of $H^TM$ associated with the positive eigenvalues. IRRDA and Theorem 2 indicate $W = MU$ and solve the RLDA problem. More importantly, it is clear that the solution space of RLDA is a subset of $span\{M\}$.

3.   A new linear system solving for RLDA

Notice from ^[8] that the solution of LDA and $S_t^{\dagger}H$ are equivalent, thus $M = S^{-1}H = (S_t+\alpha I_d)^{-1}H$ can be used to approximate the RLDA. Computing $S$ and $H$ requires forming $\bar{X} = X(I_n-\frac{1}{n}{\bf 1}_n{\bf 1}_n^T)$, which is time-consuming for large-scale high-dimensional datasets. In large-scale data, the dimension $d$ and sample number $n$ of data are usually huge, and $\frac{1}{n}$ will be very small and even tend to 0. In this section, we consider replacing $\bar{X}$ with $X$, and characterize the rationality of solving RLDA by the following linear system:

where

It is clear that

According to the feature of the data and the number of samples, we compute $\check{M}$ in two different schemes:

(1) When $d\leq n$, we compute $\check{M}$ by (3.3).

(2) When $d > n$, we first compute the solution $\check{M}$ by (3.4).

Here, we apply the breakdown-free block conjugate gradient (BCG) algorithm ^[32] to compute the linear equations for the symmetric semi-positive definite matrix. The specific algorithm of our new model for RLDA is shown in Algorithm 3.

In Algorithm 3, the convergence of Cases 1 and 2 is similar. This article explains the convergence of Case 2 and provides a convergence theorem below.

Theorem 3. Suppose $\check{H}$ is rank deficient with rank $r_0$ ($r_0 < c$). Let $\kappa = \lambda_{max}/\lambda_{min}$, where $\lambda_{max}$ and $\lambda_{min}$ are the maximum and nonzero minimum eigenvalues of matrix $\check{H}$, respectively. Let $E_k = \check{M}^*-\check{M}_k = [e_k^{(1)}, e_k^{(2)}, ..., e_k^{(c)}]$. The minimum error square norm $\|e_k^{(i)}\|^2$ is bounded as

where $t$ is a constant only related to $E_1$.

Proof. Because Algorithm 3 is a special case in the literature ^[32], the proof is trivial and we refer to ^[32][Theorem 4.1]. □

Remark 2. Let us briefly discuss the computational complexities of Algorithm 3. Algorithm 3 is divided into two cases, and their complexity is calculated separately. In Case 1, the main computational complexity is in Step 4, which costs about $\mathcal{O}(dnc+nc^2+cn^2+c^3)$ flops. In Case 2, the main computational complexity is in Step 10, which costs about $\mathcal{O}(dnc+nc^2+cn^2+c^3)$ flops. Generally speaking, the class number $c$ is far less than the characteristic dimension $d$ and the sample number $n$, thus the total costs are about $\mathcal{O}(dnc+cn^2)$ flops.

Algorithm 3 could perfectly avoid the disadvantages of RRDA. For example, $\check{R}_k^T\check{P}_{k-1} = {\bf 0}$, $\check{R}_k^T\check{R}_{j} = {\bf 0}$, $j < k$ are always held, and neither $T_k$, $W_k$ nor eigenvalue decomposition need to be computed in Algorithm 3, which saves a lot of workloads than RRDA.

LDADL ^[20] is an effective RLDA method for large-scale data. It also avoids forming $\bar{X}$, but adds a new component '1' to each datum ${\bf x}_j$. Denote

then LDADL solves the following problem ^[20]:

Let $\underline{M}^* = \left[\begin{array}{c} M^* \\ {\bf m}^T \\ \end{array}\right]$, Zhang et al. ^[20] showed that $M^*$ is a solution of LDA problem (1.1) when $\alpha = 0$. Moreover, it is known that (3.7) has a unique explicit solution given by

Hence,

can solve the LDA problem.

Considering that both (3.1) and (3.7) are linear systems for solving LDA without forming the centralized data $\bar{X}$, we will focus on the connection between our model and LDADL in the following part.

Denote by

Then, we have the following theorem.

Theorem 4. Let $\underline{M}^* = \left[\begin{array}{c} M^* \\ {\bf m}^T \\ \end{array}\right]$ be the solution of (3.7) and $\check{M}$ be the solution of (3.1), and assume that $\beta$, $T$ and $D$ satisfy the formulas (3.9) and (3.6), respectively. Then, $M^* = \check{M}\big(I_c-\frac{1}{1+\beta}T\big)D$ and $span\{\check{M}\} = span\{M^*\}$.

Proof. Formula (2.3) indicates that ${\bf 1}_n = Y{\bf t}$. According to (3.8), we have

Equation (3.10) implies $span\{M^*\}\subseteq span\{\check{M}\}$. Note that

and

thus $I_c-\frac{1}{1+\beta}T$ is full of rank. Recall that $D = diag(\sqrt{n_1}, \sqrt{n_2}, \ldots, \sqrt{n_c})$ is full of rank, and there holds $rank\{\check{M}\} = rank\{M^*\}$. Hence, $span\{\check{M}\} = span\{M^*\}$. □

Since we are only interested in an orthogonal basis matrix of the solution space for the LDA problem, and we know from Theorem 4 that the transformation matrix of LDADL for RLDA and the transformation matrix of our paper for RLDA are the same, we can expect that our method shall achieve performance similar to LDADL. Moreover, as LDADL and SRDA have the similar performance ^[20], the performance of LDADL, SRDA, and our method BBCGLDA are similar. Compared BBCGLDA with LDADL and SRDA, BBCGLDA has two advantages over LDADL and SRDA: (1) The size of the BBCGLDA problem is smaller than that of LDADL and SRDA. (2) Our proposed model (3.1) is a symmetric positive definite linear equation, and the SRDA and LDADL problem are least squares problems. Since linear equations is a special case of the least squares problem, the methods for SRDA and LDADL can also be used to solve our model. In addition, our proposed problem for LDA is similar with the problem proposed in ^[28], which aims to compute $(X)^TY$ by RSVD under the condition that the data matrix $X$ is full of column rank. Therefore, our algorithm is more applicable because it has no restrictions.

4.   Experiments

In this section, we perform some numerical experiments on some real-world datasets, and show the numerical behavior of IRRDA and our new framework BBCGLDA. All the experiments are run on a Hp workstation with 16 cores double Intel(R)Xeon(R) Platinum 8253 processors, and with CPU 2.20 GHz and RAM 256 GB. The operation system is 64-bit Windows 10. All the numerical results are obtained from running the MATLAB R2018b software.

Five real-world databases, including video data, text documents, face images, and audio data, are used in our experiment. The details of these databases, such as dimensionality, and the number of total samples are summarized in Table 1.

$\bullet$ The ${\tt AR}$ database^* contains over 4000 color images corresponding to 126 people's faces (70 men and 56 women). Images feature frontal view these faces with different facial expressions, illumination conditions, and occlusions (e.g., sunglasses and scarf). The pictures were taken at the CVC under strictly controlled conditions. No restrictions on wear (clothes, glasses, etc.), make-up, and hairstyle were imposed to participants. A subset of 100 with 26 images per person, i.e., 2600 images are used in our experiment. We crop and scale the images to $120 \times 165$ pixels.

^*http://www2.ece.ohio-state.edu/~aleix/ARdatabase.html

$\bullet$ The ${\tt Extended YaleB}$^† dataset contains 5760 single light source images of 10 subjects, each seen under 576 viewing conditions (9 different poses and 64 illumination conditions of each person). The images have normal, sleepy, sad, and surprising expressions. A subset of 38 persons with 64 images per person, i.e., 2432 images are used in the example. We crop and scale the images to $100 \times 100$ pixels in our experiment.

^†http://cvc.yale.edu/projects/yalefacesB/yalefacesB.html

$\bullet$ The ${\tt Youtube}$ dataset ^[33] is a large-scale video classification database. It contains 80 million YouTube video links, which are labeled as 4800 knowledge graph entities. Here, we provide a link of 5020 videos with 1595 persons' labels. In our experiments we have employed up to 90 samples of each class, resulting to a dataset of 124819 feature samples data and 1595 classes.

$\bullet$ The ${\tt 20Newsgroups}$ database^‡ is a collection of approximately 20000 newsgroup documents, partitioned (nearly) evenly across 20 different newsgroups. The data contains 18846 documents that are organized into 20 different newsgroups, each corresponding to a different topic.

^‡http://www.cad.zju.edu.cn/home/dengcai/

$\bullet$ The original ${\tt TDT2}$ (Nist Topic Detection and Tracking) corpus^§ collects during the first half of 1998 and is taken from 6 sources. It consists of 11201 on-topic documents which are classified into 96 semantic categories. In this subset, those documents appearing in two or more categories were removed, and only the largest 30 categories were kept, thus leaving us with 9394 documents in total.

^§http://www.cad.zju.edu.cn/home/dengcai/

To show the efficiency of IRRDA and BBCGLDA, we compare them with some state-of-the-art algorithms for large-scale discriminant analysis including RRDA ^[18], SRDA ^[23], LDADL ^[20], and the randomized algorithm in ^[28]. For simplicity, we call the randomized algorithm in ^[28] RandLDA. In RandLDA, we set the iteration power $q = 2$ and the over-sampling parameter $p = 20$. The target rank $r$ is chosen in the following way: $r = 100$ if the number of the training sample is less than 1000, and $r = 200$ if the number of the training sample is in the interval $(1000, \; 3500]$, otherwise, $r = 400$. For RRDA, we find that the dimension of $W_k$ and $W_{k-1}$ may not be consistent in running. In our experiments, we let $W_0 = {\bf 0}_{d\times (c-1)}$, and let $U_k$ be the eigenvectors of $T_k$ associated with the first $c-1$ largest eigenvalues, then the size of $W_k$ is always $d\times (c-1)$. In this case, the stopping criterion $\|W_k-W_{k-1}\|_F < \epsilon$ can be valid. As for SRDA and LDADL, we exploit the lsrq package provided by Cai et al ^[23], to solve the least squares problems.

For all iterative algorithms RRDA, SRDA, LDADL, and BBCGLDA, we set the maximum number of iterations to 20, the tolerance $\epsilon\leq10^{-4}$, and the regularization parameter is $\alpha = 0.01$. In all the experiments, we randomly pick $n = 30\%N, 50\%N$, and $70\%N$ samples as the training set, and the remaining samples are used as the testing set, where $N$ is the total number of samples. We make use of the nearest neighbor classifier (NN) ^[34] for classification. Each experiment will be repeated 10 times, and all the numerical results, i.e., the CPU time for training in seconds, the recognition rate, and the standard deviation (Std-Dev), are the mean from the 10 runs.

Example 4.1. In this example, the AR database was employed to demonstrate the superiority of IRRDA over RRDA and to show the efficiency of BBCGLDA both in CPU time and recognition rates. We set $\alpha = 0.01$, and Table 2 lists the numerical results of the algorithms.

First, Table 2 shows that RRDA is far inferior to other algorithms both in recognition rate and CPU time, which shows the effectiveness and rationality of our improved algorithm IRRDA of RRDA. Second, in terms of CPU time, RandLDA runs the fastest, followed by BBCGLDA. They are about three times faster than IRRDA, tens times faster than RRDA, and about thirty times faster than SRDA and LDADL. Although RandLDA is the fastest, its recognition rate is much lower than other algorithms, even about 7% lower than other algorithms when the training sets are small, expect for RRDA. In contrast, our algorithm BBCGLDA has good performance both in CPU time and recognition rate, as the CPU time of BBCGLDA is only $\mathcal{O}(0.1)$s slower than RandLDA, and the recognition rate is only $\mathcal{O}(0.001)$ lower than the highest recognition rate. In practical application, these numerical gaps can be ignored due to calculation errors and machine performance. Therefore, considering both the recognition rate and CPU time, BBCGLDA has advantages.

On the other hand, the recognition rates of IRRDA, SRDA, LDADL, and BBCGLDA are similar, due to the one to one relationship being established between their explicit solutions, refer to Theorem 3.2 in this paper, Lemma 4.3 in ^[20], and Theorem 3.6 in ^[28]. Moreover, the CPU time of SRDA and LDADL are the slowest. In fact, the AR database is a dense matrix, while SRDA and LDADL are aimed at large-scale sparse data.

Furthermore, we find that the more training samples, the higher recognition rates. That is, the recognition rates for 50%N training samples are 5% higher than those for 30%N training samples, and the recognition rates for 70%N training samples are 2% higher than those for 50%N training samples.

Example 4.2. In this example, we do experiments with the Youtube databse to show that our algorithm has an absolute advantage in high-dimensional dense large sample database. We do the numerical experiments with $\alpha = 0.01$. The numerical results are listed in Table 3.

Table 3 depicts the absolute computational dominant of RandLDA for high-dimensional large-scale dense data in terms of CPU time. When both the dimension $d$ and the number of training samples $n$ are very large, say, $n = 50\%N$, $70\%N$, all the algorithms do not work, except for BBCGLDA and RandLDA. RandLDA is 10 times faster than BBCGLDA, while the recognition rates of RandLDA are about 3% lower than that of BBCGLDA. Even with the increase of the training set, the increase of recognition rates is very weak, and the highest recognition rates of RandLDA only reach 93%. In fact, the recognition rate of RandLDA is closely related to the target rank $r$ for RSVD of training matrix $X$, but the optimal target rank $r$ is difficult to choose, which is a disadvantage of algorithm RandLDA. Next, we will explain the dependence of the algorithm RandLDA on the sampling number, and BBCGLDA is not sensitive to parameter. Therefore, our proposed method BBCGLDA is very powerful for high-dimensional and large-sample dense databases.

Example 4.3. In this example, we try to show the weak influence of the regularization parameter for our proposed algorithm BBCGLDA. Note that except for RandLDA, other algorithms including BBCGLDA, RRDA, SRDA, and LDADL have the regularization parameter. We set the regularization parameter $\alpha = \{10^{-4}, 10^{-3}, 10^{-2}, 10^{-1}, 1\}$ and analyze the numerical performance of these algorithms with the change of regularization parameters. Since RandLDA also needs to determine the target rank $r$ for RSVD of training matrix $X$, we set the target rank $r = \{50, 80,120,150,180\}$ and then focus on its numerical results. The test data is Extended YaleB. In order to observe the numerical results more clearly, we plot the figures of recognition rates VS regularization parameters, and CPU time VS regularization parameters for $n = \{30\%N, \ 50\%N, \ 70\%\}$ in Figure 1.

According to Figure 1, the recognition rate of RRDA is the worst, whose recognition rates of RRDA are about 10% lower than BBCGLDA, IRRDA, SRDA and LDADL, and the CPU time of RRDA is about 3 times slower than IRRDA, which both again emphasize the improved method for RRDA proposed in our paper. Next, we focus on the five algorithms except RRDA, namely, BBCGLDA, IRRDA, SRDA, LDADL, and RandLDA.

The first line of Figure 1 shows that the recognition rates of BBCGLDA, IRRDA, SRDA, and LDADL are comparable, which again emphasizes the validity and rationality of the Theorem 3.2. It is found that when $\alpha\leq0.1$, i.e., $\alpha = \{10^{-4}, 10^{-3}, 10^{-2}, 10^{-1}\}$, the recognition rate of the four algorithms is relatively stable and almost unchanged. The recognition rates decreased significantly when $\alpha = 1$, which is caused by the normalization of data. In order to avoid overflow during running algorithms, we first normalize the columns of data, that is, the norm of each column is 1, which means the maximum element of the dataset is far less than 1. Thus, 1 is too large for solving the multivariate linear regression for RLDA as it drowns out the information in the dataset, and one could result in the dynamic of recognition rate. Similarly, 0.1 may be a little large for some high-dimensional normalized data. The second part of Figure 1 shows that the CPU time of SRDA and LDADL is much slower than IRRDA and BBCGLDA, and BBCGLDA is the fastest algorithms of the four LDA algorithms BBCGLDA, IRRDA, SRDA and LDADL with the regularization parameter. This numerical performance is consistent with that of Example 4.1, which again shows the fast and effectiveness of our proposed algorithm. With the increase of the regularization parameter, the CPU time of IRRDA and BBCGLDA decrease gradually. Indeed, IRRDA and BBCGLDA are based on the BCG algorithm. Theorem 4.1 in ^[32] indicates that the smaller condition number of the coefficient matrix, the faster convergence of BCG, and the condition number decreases with the increase of regularization parameter. Therefore, considering the convergence speed and recognition rate, it is more appropriate to set $\alpha = 0.01$ for IRRDA and BBCGLDA.

In terms of RandLDA, its recognition rates increase with the increase of target rank. For example, when $n = 30\%N$, its recognition rate increases from about 75% to 92%; when $n = 50\%N$, its recognition rate increases about from about 79% to 96%; when $n = 70\%N$, its recognition rate increases about from about 82% to 97%, in which the maximum increase of recognition rate is about 17%. These numerical results indicate the recognition rate of RandLDA is seriously affected by the target rank $r$. However, how to choose the optimal target rank is a difficulty of RandLDA. In addition, it can be observed from the first line of Figure 1 that the highest recognition rates of RandLDA with target rank $r = 180$ are still lower than that of the other four algorithms when $\alpha\leq0.1$. If one wants to have a higher recognition rate of RandLDA, it needs to increase the target rank $r$. However, according to the second line of Figure 1, its CPU time will increase with the increase of the target rank $r$. In other words, RandLDA cannot have both the higher recognition rate and faster time, which is another disadvantage of the algorithm. Note that the recognition rate of BBCGLDA is hardly affected by the regularization parameters, and the convergence speed increases with the increase of the regularization parameters. Therefore, compared with RandLDA, BBCGLDA is a good choice for large-scale discriminant analysis.

Example 4.4. In this example, we aim to show the effectiveness of BBCGLDA for high-dimensional and sparse samples data. The TDT2 and 20Newsgroups databases are used in this example. The numerical results are listed in Tables 4 and 5.

We can see from Tables 4 and 5 that the recognition rates of RRDA and RandLDA can be comparable with others on the TDT2 dataset, while those that are lowest are on the 20Newsgroups dataset, which implies that the applicability of RRDA and RandLDA are not as wide as others. In addition, the standard deviations and recognition rates of the other four algorithms including IRRDA, BBCGLDA, SRDA and LDADL can be comparable, and these results are consistent with the those of the previous examples, which again verifies and determines the correctness of the Theorem 3.2.

In terms of CPU time, BBCGLDA, SRDA, and LDADL are the fastest; they are about 2–4 times faster than RandLDA, about 6 times faster than IRRDA, and about 30–40 times faster than RRDA. In fact, RRDA computes some eigenvalue decompositions and results in a waste of time. In contrast, IRRDA only needs to copmute one eigenvalue decomposition, and its speed is faster than RRDA. However, it also involves matrix decomposition, and runs lower than other algorithms. Moreover, RandLDA computes an RSVD which will destroy the sparse structure of data. We can find from Tables 4 and 5 that the computing speed of RandLDA for sparse data is not as fast as those of dense data, which is much slower than the other three algorithms BBCGLDA, SRDA, and LDADL, which are based on linear systems. This is because RandLDA needs to compute a randomized singular value decomposition, which will destroy the sparse structure of data matrix.

In general, BBCGLDA can effectively solve not only dense data, but also sparse data. Compared with dense data, it is dozens times faster than SRDA and LDADL, and its recognition rate is higher than that of RandLDA. Compared with sparse data, its CPU time is comparable with SRDA and LDADL, and it is about 3–5 times faster than RandLDA. Considering comprehensively, our algorithm BBCGLDA is a good choice for high-dimensional and large-scale discriminant analysis.

It can be seen from Tables 2–5 that the higher the sample size of the training set, the higher the accuracy rate. Tables 2 and 5 indicate that the accuracy of our algorithm in this paper is higher when $d\gg N$. On the contrary, when $d$ and $N$ are similar in size, the efficiency is lower. It shows that the algorithm in this paper is slow for very large-scale data.

5.   Conclusions

RRDA is an efficient algorithm for solving LDA, while its implementation framework has some defects. In this paper, we improve RRDA by re-orthogonalizing the descent direction and putting the eigenvalue decomposition outside the loop. Experiments show the numerical performance of our improved approach IRRDA is better than RRDA in terms of the CPU time, recognition rates, and standard deviations.

To reduce the storage requirement and computational complexity, we replace the centralized matrix $\bar{X}$ with the feature matrix $X$, and propose a linear system of symmetric and positive equations for solving the RLDA problem. The equivalent of our new model and RLDA is established. Numerical experiments show our proposed algorithm has an absolute advantage in computational efficiency and recognition rates both for high-dimensional large-scale dense data and high-dimensional large-scale sparse data.

Compared with other algorithms, our algorithm is the fastest, but it is slower than the RandLDA algorithm. It can be seen from the experiment that our algorithm cannot process the super large-scale data quickly, so the next idea is to use the random gradient of blocks for experimental research to enhance the processing speed of super large-scale data.

Author contributions

Wenya Shi: Material preparation, Data collection, Analysis, Study conception, Design, Write original draft; Zhixiang Chen: Material preparation, Data collection, Analysis, Study conception, Design. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The first author is supported by [National Natural Science Foundation of China] (Grant numbers [12201075]).

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.