New algorithms for approximating $$\varphi$$ -functions and their condition numbers for large sparse matrices

New algorithms for approximating $$\varphi$$ -functions and their condition numbers for large... Abstract Given an $$n$$-by-$$n$$ large sparse matrix $$A$$ with low rank or with fast decaying singular values, we propose new algorithms for consecutively approximating the functions $$\varphi_{\ell}(A)$$ as well as their 2-norm condition numbers for $$\ell=0,1,\ldots,p$$. The key is to reduce the computation of $$\varphi_{\ell}$$-functions of the large matrix to that of $$\varphi_{\ell+1}$$-functions of some smaller matrices of size $$r$$-by-$$r$$, where $$r$$ is the numerical rank of $$A$$. For storage reasons, the new algorithms need to store only two $$n$$-by-$$r$$ sparse matrices and some $$r$$-by-$$r$$ matrices, rather than some $$n$$-by-$$n$$ possibly dense matrices. Some error analysis is given. Numerical experiments illustrate the numerical behavior of our new algorithms and show the practical utility of new strategies. 1. Introduction In recent years, a great deal of attention has been focused on efficient and accurate evaluation of matrix functions closely related to $$\varphi$$-functions (Hochbruck et al., 1998; Lu, 2003; Moler & Van Loan, 2003; Berland et al., 2007; Schmelzer & Trefethen, 2007; Higham, 2008a; Skaflestad & Wright, 2009; Hochbruck & Ostermann, 2010; Niesen & Wright, 2012; Wu et al., 2016). For instance, exponential integrators make use of the matrix exponential and related matrix functions within the formulations of the numerical methods, and the evaluation of matrix functions is crucial for accuracy, stability and efficiency of exponential integrators (Hochbruck & Ostermann, 2010). The $$\varphi$$-functions are defined for scalar arguments by the integral representation as \begin{equation} \varphi_{0}(z)={\rm exp}(z)\quad{\rm and}\quad\varphi_{\ell}(z)=\frac{1}{(\ell-1)!}{\int_0^1 {\rm exp}\big((1-\theta)z\big)\theta^{\ell-1}\,{\rm d}\theta},\quad \ell=1,2,\ldots, \end{equation} (1.1) and they satisfy the recurrence relations \begin{equation} \varphi_{\ell}(z)=z\varphi_{\ell+1}(z)+\frac{1}{\ell!},\quad \ell=0,1,2,\ldots. \end{equation} (1.2) Moreover, $$\varphi_{\ell}(z)$$ can be expressed as the following power series whose radius of convergence is $$+\infty$$: \begin{equation} \varphi_{\ell}(z)=\sum_{i =0}^{\infty}\frac{z^i}{(i+\ell)!}. \end{equation} (1.3) The above definitions can be extended to matrices instead of scalars by using any of the available definitions of matrix functions (Higham, 2008a; Hochbruck & Ostermann, 2010). In a wide range of applications, such as the matrix exponential discriminant analysis (EDA) method for data dimensionality reduction (Yan & Pan, 2010; Zhang et al., 2010; Wang et al., 2011; Ahmed, 2013; Dornaika & Bosaghzadeh, 2013; Wang et al., 2014; Wu et al., 2017) and the complex network analysis method based on matrix functions (Estrada & Rodríguez-Velázquez, 2005; Estrada & Hatano, 2008; Estrada & Higham, 2010; Benzi et al., 2013), it is required to compute the matrix exponential of large-scale and low-rank matrices. Given a large-scale matrix $$A$$ with low rank or with fast decaying singular values, in this article we are interested in approximating several $$\varphi$$-functions consecutively. Let $$\sigma_j$$ be the $$j$$th largest singular value of $$A$$. By ‘fast decaying singular values’, we mean $$\sigma_j=\mathcal{O}(\rho^{-j}),\rho>1$$ or $$\sigma_j=\mathcal{O}(j^{-\alpha}),\ \alpha>1$$ (Hofmannn, 1986). Such matrices appear frequently in diverse application areas such as data dimensionality reduction (Duda et al., 2000), complex network analysis (Estrada & Higham, 2010), discretization of ill-posed operator equations that model many inverse problems (Hansen, 1998), randomized algorithms for matrix approximations (Halko et al., 2011; Mahoney, 2011), finite element discretization (Davis & Hu, 2011) and various other applications. In spite of the high demand for efficient methods to solve the matrix $$\varphi$$-functions in various fields of computational sciences there is no easy way to solve this type of problem. Indeed, when $$A$$ is large, both the computational cost and the storage requirements are demanding; moreover, $$\varphi_{\ell}(A)$$ is generally dense even if $$A$$ is sparse (Higham, 2008a). Some available methods are suitable only for medium-sized matrices. For instance, a MATLAB toolbox called EXPINT was provided by Berland et al. (2007). Kassam and Trefethen proposed to approximate $$\varphi$$-functions with a contour integral, which works well as long as the contour of integration is suitably chosen (Kassam & Trefethen, 2005). However, the contour is generally problem dependent and difficult to determine in advance. Another way is to reduce the computation of $$\varphi_{\ell}(A)$$ to that of a matrix exponential with a larger size (Sidje, 1998; Lu, 2003; Al-Mohy & Higham, 2011), which is unfeasible as $$A$$ is large. A third way is based on a modification of the scaling and squaring technique (Skaflestad & Wright, 2009), the most commonly used approach for computing the matrix exponential (Higham, 2009). The most well-known methods for approximating $$\varphi$$-functions for large sparse matrices are the Krylov subspace methods (Niesen & Wright, 2012; Wu et al., 2016). However, Krylov subspace methods are applicable to the computation of $$\varphi$$-functions on given (block) vectors, while the main aim of this article is to compute approximations of $$\varphi$$-functions of large sparse matrices. In practical calculations, it is important to know how accurate the computed solution is and how small perturbations in input data can affect outputs Higham (2002). Therefore, it is crucial to give some error analysis and to understand the sensitivity of matrix functions to perturbations in the data. Sensitivity is measured by condition numbers. For matrix functions, condition number can be expressed in terms of the norm of the Fréchet derivative, and it is often measured by using the 1-norm Higham (2008a,b). In this work, we focus on the 2-norm, for which there is a tighter relation between the norm of the Fréchet derivative and the norm of the Kronecker matrix, and we consider how to evaluate the Fréchet 2-norm condition numbers of $$\varphi$$-functions efficiently. Given a large scale matrix $$A$$ with low rank or with fast decaying singular values, we propose new algorithms for evaluating several $$\varphi$$-functions and their absolute and relative condition numbers consecutively. Our new algorithms are based on the sparse column–row approximation (SCR) of large sparse matrices (Stewart, 1999, 2005; Berry et al., 2005). An advantage of our algorithms is that there is no need to form and store the $$\varphi$$-functions of large matrices. The overhead is only to compute $$\varphi$$-functions of some $$r$$-by-$$r$$ matrices and to store two $$n$$-by-$$r$$ sparse matrices, where $$r$$ is the (numerical) rank of $$A$$. This article is organized as follows. In Section 2, we present the main algorithm and give some error analysis on this method. In Section 3, we propose a novel strategy for estimating the absolute and the relative 2-norm condition numbers of $$\varphi$$-functions. In Section 4, numerical experiments are given to show the effectiveness of our new strategies. Some concluding remarks are given in Section 5. Throughout this article, we denote by $$\widetilde{A}=XTY^{\rm T}$$ an SCR to $$A$$, and by $$r$$ the (numerical) rank of $$A$$. Let $$\|\cdot\|_2,\|\cdot\|_1$$ be the 2-norm and the 1-norm of a vector or a matrix, and let $$\|\cdot\|_{\rm F}$$ be the Frobenius norm of a matrix. We denote by $$\otimes$$ the Kronecker product, and by $${\rm vec}(\cdot)$$ the ‘vec operator’ that stacks the columns of a matrix into a long vector. Let $$I$$ and $$\mathbf{0}$$ be the identity matrix and the zero matrix, respectively, whose order is clear from context. We focus only on real matrices in this article. However, all the results can be extended to complex matrices in a similar way. 2. A new algorithm for approximating $$\varphi$$-functions consecutively In this section, we will present a new method for approximating $$\varphi$$-functions of large sparse matrices with low rank or with fast decaying singular values and give some error analysis on this method. Given an $$n\times n$$ large sparse matrix $$A$$, we first find a reduced-rank approximation $$XTY^{\rm T}$$ to it, where both $$X$$ and $$Y$$ are full column rank and $$T$$ is nonsingular. This type of problem arises in a number of applications such as information retrieval, computational biology and complex network analysis (Berry et al., 1995, 1999, 2005; Berry & Browne, 1999; Jiang & Berry, 2000; Zhang et al., 2002; Stuart & Berry, 2003). A widely used reduced-rank approximation is the truncated singular value decomposition (TSVD) (Golub & Van Loan, 2013), which is known to be optimal in the sense that the Frobenius norm $$\|A-XTY^{\rm T}\|_F$$ is minimized. Unfortunately, this method computes the full decomposition and is not suitable for very large matrices. An alternative is the randomized singular value decomposition (Halko et al., 2011; Mahoney, 2011), which generally gives results comparable to TSVD. However, for a large and sparse matrix, the situation is not so simple: since the resulting factors $$X, T$$ and $$Y$$ are generally not sparse, the storage requirements and operation counts will become proportional to the number of nonzero elements. In Stewart (1999), Stewart introduced a quasi-Gram–Schmidt algorithm that produces a sparse QR factorization to $$A$$. Based on the quasi-Gram-Schmidt algorithm, an SCR algorithm was proposed (Stewart, 1999). This algorithm first applies the quasi-Gram–Schmidt algorithm to the columns of $$A$$ to get a representative set of columns $$X$$ of $$A$$, and an upper triangular matrix $$R$$. Let the error in the corresponding reduced-rank decomposition be $$\epsilon_{\rm col}$$. It then applies the same algorithm to $$A^{\rm T}$$ to get a representative set $$Y^{\rm T}$$ of rows and another upper triangular matrix $$S$$. Let the error be $$\epsilon_{\rm row}$$. Then the SCR method seeks a matrix $$T$$ such that $$\|A-XTY^{\rm T}\|_F^2=\min$$, and the minimizer turns out to be $$ T=R^{-1}R^{-{\rm T}}(X^{\rm T}AY)S^{-1}S^{-{\rm T}}$$ (see Stewart, 1999; Berry et al., 2005). Moreover, we have \begin{equation} \|A-XTY^{\rm T}\|_F^2\leq \epsilon_{\rm col}^2+\epsilon_{\rm row}^2. \end{equation} (2.1) The matrix $$XTY^{\rm T}$$ is called an SCR to $$A$$, where $$X,Y\in\mathbb{R}^{n\times r}$$ are sparse and of full column rank, $$T\in\mathbb{R}^{r\times r}$$ is nonsingular and $$r$$ is the numerical rank of $$A$$. In this approximation, $$X$$ consists of a selection of the columns of $$A$$, and $$Y$$ consists of a selection of the rows of $$A$$, so that when $$A$$ is sparse then so are both $$X$$ and $$Y$$. Error analysis of the quasi-Gram–Schmidt algorithm was given in Stewart (2005). One is recommended to see Berry et al. (2005) and Stewart (1999) for more details on this algorithm. Given any rank-revealing decomposition of $$A$$, the following theorem indicates that the computation of a matrix function with respect to $$A$$ can be reduced to that of another matrix function with respect to an $$r\times r$$ matrix, where $$r$$ is the (numerical) rank of $$A$$. Theorem 2.1 Suppose that the power series $$f(z)=\sum_{i=0}^{\infty}\alpha_i z^i$$ has radius of convergence $$\rho$$. Let $$\widetilde{A}=XTY^{\rm T}\in\mathbb{R}^{n\times n}$$, where $$X,Y\in\mathbb{R}^{n\times r}$$ and $$T\in\mathbb{R}^{r\times r}$$. Let $$g(z)=\sum_{i=1}^{\infty}\alpha_i z^{i-1}$$ and suppose that $$\|\widetilde{A}\|<\rho$$; then \begin{align} f(\widetilde{A})=f(\mathbf{0})+Xg(Z)TY^{\rm T}, \end{align} (2.2) where $$Z=T(Y^{\rm T}X)\in\mathbb{R}^{r\times r}$$, and $$\mathbf{0}$$ stands for the zero matrix. Proof. As $$k\geq1$$, we see that $$\widetilde{A}^k=XZ^{k-1}TY^{\rm T}$$. Thus, \begin{align*} f(\widetilde{A})&=\alpha_0I+\alpha_1\widetilde{A}+\cdots+\alpha_k\widetilde{A}^k+\cdots\\ &=\alpha_0I+X\big(\alpha_1I+\alpha_2Z+\cdots+\alpha_kZ^{k-1}+\cdots\big)TY^{\rm T}\\ &=f(\mathbf{0})+Xg(Z)TY^{\rm T}. \end{align*} □ Remark 2.2 By (1.3) and Theorem 2.1, the computation of $$\varphi_{\ell}(A)$$ can be reduced to that of $$\varphi_{\ell+1}(Z)$$. More precisely, \begin{equation} \varphi_{\ell}(XTY^{\rm T})=\frac{1}{\ell!}I+X\big(\varphi_{\ell+1}(Z)T\big)Y^{\rm T},\quad \ell=0,1,2,\ldots. \end{equation} (2.3) We mention that for the matrix exponential, the main idea of Theorem 2.1 has already been given in (Celledon & Iserles, 2000, Proposition 3) and both (2.2) and (2.3) can be viewed as generalizations of this proposition. Let $$\widetilde{A}=XTY^{\rm T}$$ be an SCR to $$A$$. We make use of $$\varphi_{\ell}(\widetilde{A})$$ as an approximation to $$\varphi_{\ell}(A)$$. The main reason we choose the SCR is that both $$X$$ and $$Y$$ are sparse if $$A$$ is sparse. Given a large matrix $$A$$ with low rank or with fast decaying singular values, we have the following algorithm for consecutively approximating $$\varphi_{\ell}(\widetilde{A}),\,\ell=0,1,\ldots,p$$. Algorithm 1 An algorithm for approximating $$\varphi$$-functions of large sparse matrices consecutively. 1. Compute a reduced-rank approximation $$XTY^{\rm T}$$ to $$A$$ by using the sparse column–row approximation (SCR) algorithm. 2. Compute $$\varphi$$-functions of small-sized matrices: $$\varphi_{\ell+1}(TY^{\rm T}X),~\ell=0,1,\ldots,p$$. 3. Store $$X,Y,T$$ and $$\varphi_{\ell+1}(TY^{\rm T}X)$$ for $$\varphi_{\ell}(\widetilde{A}),~\ell=0,1,\ldots,p$$. If desired, form $$\varphi_{\ell}(\widetilde{A})$$ in terms of (2.3) and use them as approximations to $$\varphi_{\ell}(A),~\ell=0,1,\ldots,p$$. We stress that the main overhead of Algorithm 1 consists of computing the SCR approximation of $$A$$ provided that $$r$$ is small. This involves performing two quasi-Gram–Schmidt algorithms with respect to $$A$$ and $$A^{\rm T}$$. In the quasi-Gram–Schmidt algorithm, it is required to extract the pivot column from $$A$$ or $$A^{\rm T}$$ and to calculate several (sparse) matrix–vector operations. For a discussion of the computational complexities of the SCR approximation, refer to Berry et al. (2005, pp. 265, 268). In order to allow a fair comparison, the computation of the SCR approximation is also included in the computational comparisons and CPU timings of our new algorithm. Remark 2.3 Obviously, an advantage of the proposed method is its simplicity. In conventional methods, the cost amounts to $$\mathcal{O}\big((p+1) n^3\big)$$ flops for the computation of $$\varphi_{\ell}(A)$$ (Beylkin et al., 1998; Berland et al., 2007; Higham, 2008a; Skaflestad & Wright, 2009). Given a sparse reduced-rank approximation to $$A$$, Theorem 2.1 reduces the computation of $$\varphi_{\ell}(A)$$ to that of $$\varphi_{\ell+1}$$ functions with respect to the $$r$$-by-$$r$$ matrix $$Z=T(Y^{\rm T}X)$$, in $$\mathcal{O}\big((p+1) r^3\big)$$ flops, $$\ell=0,1,\ldots,p$$. For storage reasons, it needs to store only two $$n$$-by-$$r$$ sparse matrices $$X,Y$$, and some smaller matrices of size $$r$$-by-$$r$$, rather than the $$n$$-by-$$n$$ possibly dense matrices $$\varphi_{\ell}(A),~\ell=0,1,\ldots,p$$. Thus, the new method can consecutively approximate $$\varphi_{\ell}(A),\ell=0,1,\ldots,p$$ and reduce the computational complexities significantly as $$r\ll n$$. In practice, however, most data are inexact or uncertain. Indeed, even if the data were exact, the computations would be subject to rounding errors. So it is important to give error analysis and understand the sensitivity of matrix functions to perturbations in the input data. Sensitivity is measured by condition numbers. For a matrix function, the condition number can be expressed in terms of the norm of the Fréchet derivative. The Fréchet derivative $$L_f(A,E)$$ of a matrix function $$f: \mathbb{R}^{n\times n}\rightarrow\mathbb{R}^{n\times n}$$ at a point $$A\in\mathbb{R}^{n\times n}$$ is a linear mapping such that for all $$E\in\mathbb{R}^{n\times n}$$, \begin{equation} f(A+E)-f(A)-L_f(A,E)=o(\|E\|) \end{equation} (2.4) (see Higham, 2008a, p. 56). The absolute and the relative condition numbers of a matrix function $$f(A)$$ are defined as \begin{equation} {\rm cond}_{\rm abs}(f,A)=\|L_f(A)\|=\max_{E\neq \mathbf{0}}\frac{\|L_f(A,E)\|}{\|E\|} \end{equation} (2.5) and \begin{equation} {\rm cond}_{\rm rel}(f,A)=\frac{\|L_f(A)\|\|A\|}{\|f(A)\|}, \end{equation} (2.6) respectively (see Rice, 1966; Higham, 2008a). Theoretically, the Fréchet derivative can be obtained from applying any existing method for computing the matrix function of a $$2n\times 2n$$ matrix and using the formula \begin{equation} f\left(\left[\begin{array}{@{}cc@{}} A & E \\ \mathbf{0} & A \\\end{array} \right]\right) =\left[\begin{array}{@{}cc@{}} f(A) & L_f(A,E) \\ \mathbf{0} & f(A)\\ \end{array} \right] \end{equation} (2.7) (see Higham, 2008a, Algorithm 3.17). However, it requires $$\mathcal{O}(n^5)$$ flops, assuming that the computation of $$L_f(A,E)$$ takes $$\mathcal{O}(n^3)$$ flops (Higham, 2008a). So it is computationally demanding for large matrices. Let $$\widetilde{A}=XTY^{\rm T}$$ be an SCR to $$A$$, and let $$E=A-\widetilde{A}$$; then we see from (2.1) that $$\|E\|_F\leq\sqrt{\epsilon_{\rm col}^2+\epsilon_{\rm row}^2}$$. Thus, it is interesting to combine the existing error analysis for SCR with the theory of matrix functions to obtain error bounds for the new method. We first present the following theorem for Fréchet derivatives of $$\varphi$$-functions. Theorem 2.4 For any matrix $$E\in\mathbb{R}^{n\times n}$$, we have \begin{equation} \varphi_{\ell}(A+E)-\varphi_{\ell}(A)=\int_{0}^{1}\exp\big((1-s)A\big)s^{\ell}E\varphi_{\ell}\big(s(A+E)\big)\,{\rm d}s,\quad \ell=0,1,2,\ldots. \end{equation} (2.8) Proof. For the matrix exponential (Higham, 2008a, p. 238), we have \begin{equation} \exp(A+E)=\exp(A)+\int_0^1\exp\big((1-s)A\big)E\exp\big(s(A+E)\big)\,{\rm d}s, \end{equation} (2.9) and thus (2.8) holds for $$\ell=0$$. When $$\ell>0$$, we know that the function $$X(t)=t^{\ell}\varphi_{\ell}\big(t(A+E)\big)$$ is the solution of the initial value problem (Niesen & Wright, 2012, Lemma 2.1) $$ X'(t)=(A+E)X(t)+\frac{t^{\ell-1}}{(\ell-1)!}I=AX(t)+G(t,X),\quad X(0)={\bf 0}, $$ where $$ G(t,X)=EX(t)+\frac{t^{\ell-1}}{(\ell-1)!}I. $$ By the variation of constants formula, we get \begin{eqnarray*} X(t)&=&\exp(tA)X(0)+\int_0^t \exp\big((t-s)A\big)G\big(s,X(s)\big)\,{\rm d}s\\ &=&\int_0^t \exp\big((t-s)A\big)\left(EX(s)+\frac{s^{\ell-1}}{(\ell-1)!}I\right)\,{\rm d}s\\ &=&\int_0^t \exp\big((t-s)A\big)EX(s)\,{\rm d}s+\int_0^t \exp\big((t-s)A\big)\frac{s^{\ell-1}}{(\ell-1)!}\,{\rm d}s. \end{eqnarray*} Notice that $$X(s)=s^{\ell}\varphi_{\ell}\big(s(A+E)\big)$$ and $$t^{\ell}\varphi_{\ell}(tA)$$ is the solution of the initial value problem $$ U'(t)=AU(t)+\frac{t^{\ell-1}}{(\ell-1)!}I,\quad U(0)={\bf 0}. $$ Similarly, from the variation of constants formula, we have $$ t^{\ell}\varphi_{\ell}(tA)=\int_0^t \exp\big((t-s)A\big)\frac{s^{\ell-1}}{(\ell-1)!}\,{\rm d}s. $$ Therefore, $$X(t)=\int_0^t \exp\big((t-s)A\big)Es^\ell\varphi_{\ell}\big(s(A+E)\big)\,{\rm d}s+t^{\ell}\varphi_{\ell}(tA),$$ from which we get our result (2.8) for $$t=1$$. □ Using (2.8) to substitute for $$\varphi_{\ell}\big(s(A+E)\big)$$ inside the integral, we get \begin{equation} \varphi_{\ell}(A+E)=\varphi_{\ell}(A)+\int_{0}^{1}\exp\big((1-s)A\big)s^{\ell}E\varphi_{\ell}(sA)\,{\rm d}s+\mathcal{O}(\|E\|^2),\quad \ell=0,1,2,\ldots. \end{equation} (2.10) Combining with (2.4), we can present the following result for the Fréchet derivatives of $$\varphi$$-functions. Corollary 2.5 The Fréchet derivatives of $$\varphi$$-functions at $$A$$ in the direction $$E$$ are given by \begin{equation} L_{\varphi_{\ell}}(A,E)=\int_{0}^{1}\exp\big((1-s)A\big)s^{\ell}E\varphi_{\ell}(sA)\,{\rm d}s,\quad \ell=0,1,2,\ldots \end{equation} (2.11) As a result, \begin{equation} \varphi_{\ell}(A+E)=\varphi_{\ell}(A)+L_{\varphi_{\ell}}(A,E)+o(\|E\|). \end{equation} (2.12) The following theorem shows that the values of $$\epsilon_{\rm col}$$ and $$\epsilon_{\rm row}$$ used during the SCR will have a direct impact on the final accuracy of approximating $$\varphi_{\ell}(A)$$. Note that $$XTY^{\rm T}$$ can be any low-rank approximation to $$A$$ in this theorem. Theorem 2.6 Let $$XTY^{\rm T}$$ be an SCR to $$A$$, and $$\|A-XTY^{\rm T}\|=\varepsilon$$. Then we have \begin{equation} \|\varphi_{\ell}(A)-\varphi_{\ell}(XTY^{\rm T})\|\lesssim{\rm cond}_{\rm abs}(\varphi_{\ell},A)~\varepsilon \end{equation} (2.13) and \begin{equation} \frac{\|\varphi_{\ell}(A)-\varphi_{\ell}(XTY^{\rm T})\|}{\|\varphi_{\ell}(A)\|}\lesssim{\rm cond}_{\rm rel}(\varphi_{\ell},A)\frac{\varepsilon}{\|A\|}, \end{equation} (2.14) where $$\lesssim$$ represents omitting the high-order term $$o(\|E\|)$$. Proof. Let $$E=XTY^{\rm T}-A$$; then we get from (2.12) that \begin{eqnarray*} \|\varphi_{\ell}(A)-\varphi_{\ell}(XTY^{\rm T})\|&\lesssim&\|L_{\varphi_{\ell}}(A)\|\|E\| =\|L_{\varphi_{\ell}}(A)\|\varepsilon\nonumber\\ &=&{\rm cond}_{\rm abs}(\varphi_{\ell},A)\varepsilon, \end{eqnarray*} where we omit the high-order term $$o(\|E\|)$$. The upper bound (2.14) for the relative error is derived from (2.6) and (2.13). □ 3. A new strategy for estimating the absolute and the relative 2-norm condition numbers consecutively In terms of Theorem 2.6, it is crucial to consider how to estimate the absolute and the relative condition numbers $${\rm cond}_{\rm abs}(\varphi_{\ell},A)$$ and $${\rm cond}_{\rm rel}(\varphi_{\ell},A)$$ efficiently. Notice that \begin{equation} {\rm vec}(ACB)=(B^{\rm T}\otimes A){\rm vec}(C). \end{equation} (3.1) Since $$L_f$$ is a linear operator, we have \begin{equation} {\rm vec}(L_f(A,E))=K_f(A){\rm vec}(E) \end{equation} (3.2) for some $$K_f(A)\in\mathbb{R}^{n^2\times n^2}$$ that is independent of $$E$$. The matrix $$K_f(A)$$ is referred to as the Kronecker form of the Fréchet derivative (Higham, 2008a, p. 60). Specifically, we have \begin{equation} \|L_f(A)\|_F=\lambda_{\max}\big(K_f(A)^{\rm T} K_f(A)\big)^{1/2}=\|K_f(A)\|_2. \end{equation} (3.3) To estimate $$\|K_f(A)\|_2$$, the power method can be applied (Higham, 2008a, Algorithm 3.20). However, the power method lacks convergence tests, and because of its linear convergence rate the number of iterations required is unpredictable. Alternatively, the condition number can be based on the 1-norm (Higham, 2008a, Algorithm 3.22). One reason for using the 1-norm is its cheapness. Although there is no analogue to the relation (3.3) for the 1-norm, the next result gives a relation between $$\|K_f (A)\|_1$$ and $$\|L_f (A)\|_1$$. Theorem 3.1 (Higham, 2008a, Lemma 3.18). For $$A\in\mathbb{R}^{n\times n}$$ and any function $$f$$, \begin{equation} \frac{\|L_f(A)\|_1}{n}\leq \|K_f(A)\|_1\leq n\|L_f(A)\|_1. \end{equation} (3.4) In this article, we choose to work with the 2-norm, for which there is a tighter relation between the norm of the Fréchet derivative and the norm of the Kronecker matrix. The following theorem establishes a relation between $$\|K_f (A)\|_2$$ and $$\|L_f (A)\|_2$$. Theorem 3.2 For $$A\in\mathbb{R}^{n\times n}$$ and any function $$f$$, \begin{equation} \frac{\|L_f(A)\|_2}{\sqrt{n}}\leq \|K_f(A)\|_2\leq \sqrt{n}\|L_f(A)\|_2. \end{equation} (3.5) Proof. The proof is analogous to that by Higham (2008a, Lemma 3.18). □ Remark 3.3 Notice that $$\sqrt{n}$$ is used in the lower and upper bounds of (3.5), instead of $$n$$ in (3.4). Hence, there is a tighter relation between the 2-norm of the Fréchet derivative and the 2-norm of the Kronecker matrix. We are ready to propose a novel strategy to evaluate the absolute and the relative 2-norm condition numbers. Recall that the key idea of our new strategy in Section 2 is to reduce the computation of $$\varphi_{\ell}(A)$$ to that of $$\varphi_{\ell+1}(Z)$$. As $$ \varphi_{\ell}(\widetilde{A})=\frac{1}{\ell!}I+X\big(\varphi_{\ell+1}(Z)T\big)Y^{\rm T}, $$ if we denote $$ \widetilde{\varphi_{\ell}(\widetilde{A})}=\frac{1}{\ell!}I+X\big(\varphi_{\ell+1}(Z+F)T\big)Y^{\rm T} $$ then it follows from (2.4) that $$ L_{\varphi_{\ell+1}}(Z,F)=\varphi_{\ell+1}(Z+F)-\varphi_{\ell+1}(Z)+o(\|F\|)\quad \forall\,F\in\mathbb{R}^{r\times r}. $$ Hence, \begin{eqnarray} \widetilde{\varphi_{\ell}(\widetilde{A})}-\varphi_{\ell}(A)&=&X\big(\varphi_{\ell+1}(Z+F)-\varphi_{\ell+1}(Z)\big)TY^{\rm T}\nonumber\\ &=&XL_{\varphi_{\ell+1}}(Z,F)TY^{\rm T}+o(\|F\|). \end{eqnarray} (3.6) Let $$X=Q_1R_1,Y=Q_2R_2$$ be the (sparse) QR decompositions of $$X$$ and $$Y$$, respectively, where $$R_1,R_2\in\mathbb{R}^{r\times r}$$. We propose to use \begin{equation} {\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})=\|XL_{\varphi_{\ell+1}}(Z)TY^{\rm T}\|_2=\|R_1L_{\varphi_{\ell+1}}(Z)TR_2^{\rm T}\|_2 \end{equation} (3.7) as an approximation to the absolute 2-norm condition number $${\rm cond}_{\rm abs}(\varphi_{\ell},A)$$. Notice that the key of our strategy is (3.6) and the fact that both $$Q_1$$ and $$Q_2$$ are orthonormal. On the other hand, Theorem 2.1 also provides a cheap way to estimate 2-norms of $$\varphi_{\ell}(A),~\ell=0,1,\ldots,p$$. Indeed, from (2.3) we get \begin{equation} \Big|\|R_1\big(\varphi_{\ell+1}(Z)T\big)R_2^{\rm T}\|_2-1/\ell!\Big|\leq\|\varphi_{\ell}(\widetilde{A})\|_2\leq 1/\ell!+\|R_1\varphi_{\ell+1}(Z)TR_2^{\rm T}\|_2. \end{equation} (3.8) Thus, we exploit \begin{equation} \eta_{\ell}=\big\|R_1\varphi_{\ell+1}(Z)TR_2^{\rm T}\big\|_2,\quad \ell=0,1,\ldots,p \end{equation} (3.9) as approximations to $$\|\varphi_{\ell}({A})\|_2$$. Then the relative 2-norm condition number $${\rm cond}_{\rm rel}(\varphi_{\ell},A)$$ can be approximated by \begin{equation} {\rm cond}_{\rm rel}^{\rm est}(\varphi_{\ell},\widetilde{A})=\frac{\|A\|_2\|R_1L_{\varphi_{\ell+1}}(Z)TR_2^{\rm T}\|_2}{\|R_1\varphi_{\ell+1}(Z)TR_2^{\rm T}\|_2}; \end{equation} (3.10) refer to (2.14). In summary, we have the following algorithm for estimating 2-norm condition numbers of $$\varphi$$-functions. Algorithm 2 An algorithm for estimating the absolute and the relative 2-norm condition numbers of $$\varphi$$-functions consecutively. 1. Compute a reduced-rank approximation $$XTY^{\rm T}$$ to $$A$$ by using the sparse column–row approximation (SCR) algorithm. 2. Compute the R-factors $$R_1$$ and $$R_2$$ by using the (sparse) QR decompositions of $$X$$ and $$Y$$, respectively. 3. Solve $$\varphi_{\ell+1}(Z)$$ and $$L_{\varphi_{\ell+1}}(Z)$$ for the $$r$$-by-$$r$$ matrix $$Z=T(Y^{\rm T}X),~\ell=0,1,\ldots,p$$. 4. Use (3.7) and (3.10) as approximations to the absolute and the relative 2-norm condition numbers of $$\varphi_{\ell}(A),~\ell=0,1,\ldots,p$$. Remark 3.4 Two remarks are given. First, there is no need to form and store the Q-factors $$Q_1$$ and $$Q_2$$ in the QR decompositions of $$X$$ and $$Y$$, which may be dense in practice. Indeed, it is only necessary to evaluate 2-norms of some $$r$$-by-$$r$$ matrices in Algorithm 2. Second, similarly to Algorithm 1, one can evaluate the relative and absolute 2-norm condition numbers consecutively for $$\ell=0,1,\ldots,p$$. Next, we will give another motivation for using $$L_{\varphi_{\ell+1}}(Z)$$ for approximating the norm of $$L_{\varphi_{\ell}}(\widetilde{A})$$, where $$\widetilde{A}=XTY^{\rm T}$$. Indeed, the theorem below motivates computing the norm of the approximating formula (3.7) and (3.10). Theorem 3.5 With the notation of Theorem 2.1, if we define $$\widetilde{A}=XTY^{\rm T}$$ and $$W=YT^{\rm T}$$ then \begin{equation} K_f(\widetilde{A})={\it{\Psi}}_1+{\it{\Psi}}_2+{\it{\Psi}}_3, \end{equation} (3.11) where \begin{align*} {\it{\Psi}}_1 & =\alpha_1 I\otimes I,\\ {\it{\Psi}}_2 & =\left(W\otimes I\right)\sum_{i=2}^{\infty}\alpha_i\left((Z^{\rm T})^{i-2}\otimes I\right)\left(X\otimes I\right)^{\rm T}+ \left(I\otimes X\right)\sum_{i=2}^{\infty}\alpha_i\left(I\otimes Z^{i-2}\right)\left(I\otimes W\right)^{\rm T} \end{align*} and $$ {\it{\Psi}}_3=\left(W\otimes X\right)\left(\sum_{i=3}^{\infty}\alpha_i\sum_{j=2}^{i-1}\big((Z^{\rm T})^{i-j-1}\otimes Z^{j-2}\big)\right) \left(X\otimes W\right)^{\rm T}, $$ where the $$\{\alpha_i\}$$’s are the coefficients from the expression of $$f(x)$$. Proof. For any $$E\in\mathbb{R}^{n\times n}$$, from (2.7) and the expression of $$f(x)$$ we have \begin{eqnarray} L_f(\widetilde{A},E)&=&\sum_{i=1}^{\infty}\alpha_i\sum_{j=1}^{i}\widetilde{A}^{j-1}E\widetilde{A}^{i-j}\nonumber\\ &=&\alpha_1E+\alpha_2(E\widetilde{A}+\widetilde{A}E)+\sum_{i=3}^{\infty}\alpha_i \left(E\widetilde{A}^{i-1}+\sum_{j=2}^{i-1}\widetilde{A}^{j-1}E\widetilde{A}^{i-j}+\widetilde{A}^{i-1}E\right)\nonumber\\ &=&\alpha_1E+\sum_{i=2}^{\infty}\alpha_i(E\widetilde{A}^{i-1}+\widetilde{A}^{i-1}E)+\sum_{i=3}^{\infty}\alpha_i\sum_{j=2}^{i-1}\widetilde{A}^{j-1}E\widetilde{A}^{i-j}. \end{eqnarray} (3.12) As a result, \begin{equation} {\rm vec}(\alpha_1E)=(\alpha_1 I\otimes I){\rm vec}(E)={\it{\Psi}}_1 {\rm vec}(E). \end{equation} (3.13) If $$i\geq 2$$, we obtain from $$\widetilde{A}^{i-1}=XZ^{i-2}TY^{\rm T}=XZ^{i-2}W^{\rm T}~$$ that \begin{eqnarray} {\rm vec}\left(\sum_{i=2}^{\infty}\alpha_i(E\widetilde{A}^{i-1}+\widetilde{A}^{i-1}E)\right)&=&{\rm vec}\left(\sum_{i=2}^{\infty}\alpha_i\big(EXZ^{i-2}W^{\rm T}+XZ^{i-2}W^{\rm T}E\big)\right)\nonumber\\ &=&\sum_{i=2}^{\infty}\alpha_i\big(W(Z^{\rm T})^{i-2}X^{\rm T}\otimes I+I\otimes XZ^{i-2}W^{\rm T}\big){\rm vec}(E)\nonumber\\ &=&\left(\sum_{i=2}^{\infty}\alpha_i(W\otimes I)\big((Z^{\rm T})^{i-2}\otimes I\big)(X\otimes I)^{\rm T}\right.\nonumber\\ &&\left.+\sum_{i=2}^{\infty}\alpha_i(I\otimes X)\big(I\otimes Z^{i-2}\big)(I\otimes W)^{\rm T}\big]\right){\rm vec}(E)\nonumber\\ &=&{\it{\Psi}}_2{\rm vec}(E). \end{eqnarray} (3.14) Moreover, \begin{eqnarray} {\rm vec}\left(\sum_{i=3}^{\infty}\alpha_i\sum_{j=2}^{i-1}\widetilde{A}^{j-1}E\widetilde{A}^{i-j}\right)&=&\sum_{i=3}^{\infty}\alpha_i\sum_{j=2}^{i-1}{\rm vec}\big(XZ^{j-2}W^{\rm T} E XZ^{i-j-1}W^{\rm T}\big)\nonumber\\ &=&\sum_{i=3}^{\infty}\alpha_i\sum_{j=2}^{i-1}\big((W(Z^{\rm T})^{i-j-1}X^{\rm T}) \otimes (XZ^{j-2}W^{\rm T})\big){\rm vec}(E)\nonumber\\ &=&\sum_{i=3}^{\infty}\alpha_i\sum_{j=2}^{i-1}(W\otimes X)\big((Z^{\rm T})^{i-j-1}\otimes Z^{j-2}\big) (X\otimes W)^{\rm T}{\rm vec}(E)\nonumber\\ &=&{\it{\Psi}}_3{\rm vec}(E), \end{eqnarray} (3.15) and (3.11) follows from (3.12)–(3.15). □ Remark 3.6 We note that \begin{equation} K_g(Z)=\sum_{i=2}^{\infty}\alpha_i\sum_{j=1}^{i-1}\big((Z^{\rm T})^{i-j-1}\otimes Z^{j-1}\big) \end{equation} (3.16) (see Higham, 2008a, Problem 3.6). As a result, Theorem 3.5 indicates that $$L_f(\widetilde{A})$$ and $$L_g(Z)$$ are closely related. More precisely, if we define $$ \psi_i(Z)=\sum_{j=1}^{i-1}\big((Z^{\rm T})^{i-j-1}\otimes Z^{j-1}\big), \quad i=2,3,\ldots, $$ then $$K_g(Z)=\sum_{i=2}^{\infty}\alpha_i\psi_i(Z)$$ and $${\it{\Psi}}_3=\big(W\otimes X\big)\big(\sum_{i=3}^{\infty}\alpha_i\psi_{i-1}(Z)\big)\big(X\otimes W\big)^{\rm T}$$. 4. Numerical experiments In this section, we perform some numerical experiments to illustrate the numerical behavior of Algorithm 1 and Algorithm 2. All the numerical experiments were run on a Dell PC with eight cores Intel(R) Core(TM)i7-2600 processor with 3.4 GHz and 16.0 GB of RAM. As an operating system, we used Windows 7 with 64-bit architecture. All the numerical results were obtained from MATLAB R2015b implementations with unit roundoff $$\epsilon\approx 2.22\times 10^{-16}$$ (Higham, 2002, Chapter 2). In all the examples, the SCR of $$A$$ is computed using the MATLAB functions $${\tt scra.m}$$ and $${\tt spqr.m}$$ due to G.W. Stewart,1 where the tolerance $${\tt tol}$$ is taken as $$\epsilon_{\rm col}=\epsilon_{\rm row}=10^{-5}$$. The matrix exponential is calculated by using the MATLAB built-in function $${\tt expm.m}$$, while the $$\varphi_{\ell}\,(\ell\geq 1)$$ functions are computed using the $${\tt phipade.m}$$ function in the MATLAB package EXPINT (Berland et al., 2007). As was pointed out in Remark 2.3, if $$A$$ is large and sparse then $$\varphi_{\ell}(A)$$ can be dense and there is no need to form or store them in practice. In this section, however, we form the matrix functions $$\varphi_{\ell}(\widetilde{A})$$, so that the relative and the absolute errors of $$\varphi_{\ell}(\widetilde{A})$$ with respect to $$\varphi_{\ell}(A)$$ can be explicitly computed for $$\ell=0,1,\ldots,p$$. 4.1 An application to data dimensionality reduction In this example, we show the efficiency of our new method for computing the matrix exponential of large scale-and low-rank matrices. Many data-mining problems involve data sets represented in very high-dimensional spaces. In order to handle high-dimensional data, the dimensionality needs to be reduced. Linear discriminant analysis (LDA) is one of notable subspace transformation methods for dimensionality reduction (Duda et al., 2000). LDA encodes discriminant information by maximizing the between-class scatter while minimizing the within-class scatter in the projected subspace. Let $$X=[{\bf a}_{1},{\bf a}_{2},\ldots,{\bf a}_{m}]$$ be a set of training samples in an $$n$$-dimensional feature space, and assume that the original data is partitioned into $$K$$ classes as $$X=[X_1,X_2,\ldots,X_K]$$. We denote by $$m_{j}$$ the number of samples in the $$j$$th class, and thus $$\sum_{j=1}^{K}m_{j}=m$$. Let $${\bf c}_{j}$$ be the centroid of the $$j$$th class, and $${\bf c}$$ be the global centroid of the training data set. If we define $${\rm 1}\kern-0.24em{\rm I}_{j}=[1,1,\ldots,1]^{\rm T} \in \mathbb{R}^{m_{j}}$$ then the within-class scatter matrix is defined as $$ S_{\rm W}=\sum_{j=1}^{K}\sum_{{\bf a}_{i}\in X_{j}}({\bf a}_{i}-{\bf c}_{j})({\bf a}_{i}-{\bf c}_{j})^{\rm T}=H_WH_W^{\rm T}, $$ where $$H_{\rm W}=[X_{1}-{\bf c}_{1}{\rm 1}\kern-0.24em{\rm I}_{1}^{\rm T},\ldots,X_{K}-{\bf c}_{K}{\rm 1}\kern-0.24em{\rm I}_{K}^{\rm T}]\in\mathbb{R}^{n\times m}$$. The between-class scatter matrix is defined as $$ S_{\rm B}=\sum_{j=1}^{K}n_{j}({\bf c}_{j}-{\bf c})({\bf c}_{j}-{\bf c})^{\rm T}=H_BH_B^{\rm T}, $$ where $$H_B=[\sqrt{n_{1}}({\bf c}_{1}-{\bf c}),\sqrt{n_{2}}({\bf c}_{2}-{\bf c}),\ldots,\sqrt{n_{K}}({\bf c}_{K}-{\bf c})]\in\mathbb{R}^{n \times K}$$. The optimal projection matrix in the LDA method can be obtained by solving the large-scale generalized eigenproblem \begin{equation} S_{\rm B}{\bf x}= \lambda S_{\rm W}{\bf x}. \end{equation} (4.1) In practice, the dimension of real data usually exceeds the number of training samples (i.e., $$n\gg m$$), which is called the small-sample size (SSS) or the undersampled problem (Duda et al. 2000; Park & Park 2008). It is an intrinsic limitation of the classical LDA method and is also a common problem in classification applications (Park & Park, 2008). Note that both $$S_W$$ and $$S_B$$ are singular. Indeed, if $$X$$ is linearly independent then $${\rm rank}(S_B)=K-1$$ and $${\rm rank}(S_W)=m-K$$ (Duda et al., 2000; Wu & Feng, 2015). Thus, the ranks of $$S_B$$ and $$S_W$$ are much smaller than the dimensionality $$n$$. In other words, the SSS problem stems from generalized eigenproblems with singular matrices. To rectify this drawback, a novel method based on the matrix exponential, called the EDA method, was proposed in Zhang et al. (2010). Instead of (4.1), the EDA method solves the following generalized matrix exponential eigenproblem (Zhang et al., 2010; Wu et al., 2017): \begin{equation} \textrm{exp}(S_{\rm B}){\bf x}= \lambda\,\textrm{exp}(S_{\rm W}){\bf x}. \end{equation} (4.2) As both $$\textrm{exp}(S_{\rm W})$$ and $$\exp(S_B)$$ are symmetric positive definite (SPD), this is a definite generalized eigenproblem, and there is a complete set of eigenvectors (Davies et al., 2001; Higham, 2009); moreover, all the eigenvalues are real. The EDA method is described as follows; for more details, refer to Zhang et al. (2010). Algorithm 3 Zhang et al. (2010) The exponential discriminant analysis (EDA). Input: The data matrix $$X=[{\bf a}_{1},{\bf a}_{2},\ldots,{\bf a}_{m}]\in\mathbb{R}^{n\times m}$$, where $${\bf a}_{j}$$ represents the $$j$$th training image. Output: The projection matrix $$V$$. 1. Compute the matrices $$S_{\rm B}$$, $$S_{\rm W}$$, $${\rm exp}(S_{\rm B})$$, and $${\rm exp}(S_{\rm W})$$. 2. Compute the eigenvalues $$\{\lambda_{i}\}$$ and associated eigenvectors $$\{{\bf x}_{i}\}$$ of $${\rm exp}(S_{\rm W})^{-1}{\rm exp}(S_{\rm B})$$. 3. Sort the eigenvectors $$V=\{{\bf x}_{i}\}$$ according to $$\{\lambda_{i}\}$$ in decreasing order. 4. Orthogonalize the columns of the projection matrix $$V$$. The framework of the EDA method for dimensionality reduction has gained wide attention in recent years (Yan & Pan, 2010; Zhang et al., 2010; Wang et al., 2011; Ahmed, 2013; Dornaika & Bosaghzadeh, 2013; Wang et al., 2014). However, the time complexity of Algorithm 3 is dominated by the computation of $${\rm exp}(S_B)$$ and $${\rm exp}(S_W)$$, which is demanding since the data dimension $$n$$ is large (Zhang et al., 2010). Moreover, the problem of (4.2) should not be solved as in Algorithm 3 and not by the inverting transformation shown. One refers to a recent article on inexact Krylov subspace methods for solving this type of generalized eigenproblem (Wu et al., 2017). By Theorem 2.1, we can compute the large matrix exponentials as follows. Corollary 4.1 Under the above notation, we have \begin{equation} \exp(S_B)=I+H_B\big(\varphi_1(H_B^{\rm T}H_B)\big) H_B^{\rm T} \end{equation} (4.3) and \begin{equation} \exp(S_W)=I+H_W\big(\varphi_1(H_W^{\rm T}H_W)\big) H_W^{\rm T}. \end{equation} (4.4) So we have the following algorithm for the matrix EDA. Algorithm 4 An algorithm for approximating $$\exp(S_B)$$ and $$\exp(S_W)$$. 1. Given the data matrix $$X=[{\bf a}_{1},{\bf a}_{2},\ldots,{\bf a}_{m}]\in\mathbb{R}^{n\times m}$$, form $$H_B$$ and $$H_W$$. 2. Compute $$\varphi_1(H_B^{\rm T}H_B)$$ and $$\varphi_1(H_W^{\rm T}H_W)$$. 3. Store $$H_B,H_W$$ and $$\varphi_1(H_B^{\rm T}H_B)$$, $$\varphi_1(H_W^{\rm T}H_W)$$ for $$\exp(S_B)$$ and $$\exp(S_W)$$. If desired, form $$\exp(S_B)$$ and $$\exp(S_W)$$ in terms of (4.3) and (4.4). Note that both $$H_B$$ and $$H_W$$ are already available in step 1, and so there is no need to perform rank-revealing decompositions of $$S_B$$ and $$S_W$$. As a result, the computation of the two $$n\times n$$ matrix exponentials $$\exp(S_B),\exp(S_W)$$ reduces to that of $$\varphi_1(H_B^{\rm T}H_B)\in\mathbb{R}^{K\times K}$$ and $$\varphi_1(H_W^{\rm T}H_W)\in\mathbb{R}^{m\times m}$$, where $$K,m\ll n$$. Next we illustrate the efficiency of Algorithm 4 for matrix EDA. There are three real-world databases in this example. The first one is the $${\tt ORL}$$ database2 that contains 400 face images of $$40$$ individuals, and the original image size is $$92\times 112=10304$$. The second test set is the $${\tt Yale}$$ face database taken from the Yale Center for Computational Vision and Control.3 It contains $$165$$ grayscale images of 15 individuals. The original image size is $$320\times 243=77760$$. The third test set is the $${\tt Extended\,YaleB}$$ database.4 This database contains 5760 single light source images of 10 subjects, each seen under 576 viewing conditions. A subset of $$38$$ classes with 2432 images are used in this example, 64 images per individual with illumination. In the $${\tt ORL}$$ database, the images are aligned based on eye coordinates and are cropped and scaled to $$n=32\times 32$$ and $$64\times 64$$, respectively; the original image size with $$n=92\times 112$$ is also considered. In the $${\tt Yale}$$ and the $${\tt Extended\,YaleB}$$ databases, all images are aligned based on eye coordinates and are cropped and scaled to $$n=32\times 32,~64\times 64$$ and $$100\times 100$$, respectively. In this example, a random subset with $$3$$ images per subject is taken to form the training set, and the rest of the images are used as the testing set. Each column of the data matrices is scaled by its 2-norm. In Algorithm 4, the CPU time consists of computing $$H_B,H_W$$, $$\varphi_1(H_B^{\rm T}H_B)$$ and $$\varphi_1(H_W^{\rm T}H_W)$$ and forming $$\exp(S_B)$$ and $$\exp(S_W)$$ in terms of (4.3) and (4.4). In the original EDA algorithm (Algorithm 3), the CPU time consists of forming $$H_B,H_W$$ and the computation of $$\exp(S_B)$$ and $$\exp(S_W)$$ by using the MATLAB built-in function $${\tt expm.m}$$. Let $$\exp(S_B),\exp(S_W)$$ be the ‘exact solutions’ obtained from running $${\tt expm.m}$$, and let $$\widetilde{\exp(S_B)},\widetilde{\exp(S_W)}$$ be the approximations obtained from (4.3) and (4.4). We define $$ {\bf Err_B}=\frac{\|\exp(S_B)-\widetilde{\exp(S_B)}\|_F}{\|\exp(S_B)\|_F}\quad {\rm and} \quad {\bf Err_W}=\frac{\|\exp(S_W)-\widetilde{{\rm exp}(S_W)}\|_F}{\|\exp(S_W)\|_F} $$ as the relative errors of the approximations $$\widetilde{{\rm exp}(S_B)},\widetilde{{\rm exp}(S_W)}$$, respectively, and we denote by $$ {\tt Rel\_ErrF}=\max({\bf Err_B},{\bf Err_W}) $$ the maximal value of the two relative errors. Table 1 lists the CPU time in seconds of Algorithm 4, $${\tt expm.m}$$ and the values of the maximal relative errors $${\tt Rel\_ErrF}$$. Table 1 Example $$1$$: CPU time in seconds and the relative errors for approximating $$\exp(S_B)$$ and $$\exp(S_W)$$ $${\tt Database}$$ $${\tt n}$$ $${\tt Algorithm\,4}$$ $${\tt expm.m}$$ $${\tt Rel\_ErrF}$$ ORL 1024 0.08 0.26 $$2.20\times 10^{-15}$$ $$4096$$ 0.26 17.3 $$3.22\times 10^{-15}$$ $$10304$$ 1.28 261.1 $$4.49\times 10^{-15}$$ Yale $$1024$$ 0.08 0.25 $$2.11\times 10^{-15}$$ $$4096$$ 0.24 17.3 $$3.10\times 10^{-15}$$ $$10000$$ 1.13 238.5 $$4.35\times 10^{-15}$$ Extended YaleB $$1024$$ 0.08 0.27 $$2.13\times 10^{-15}$$ $$4096$$ 0.27 17.3 $$3.14\times 10^{-15}$$ $$10000$$ 1.22 238.6 $$4.50\times 10^{-15}$$ $${\tt Database}$$ $${\tt n}$$ $${\tt Algorithm\,4}$$ $${\tt expm.m}$$ $${\tt Rel\_ErrF}$$ ORL 1024 0.08 0.26 $$2.20\times 10^{-15}$$ $$4096$$ 0.26 17.3 $$3.22\times 10^{-15}$$ $$10304$$ 1.28 261.1 $$4.49\times 10^{-15}$$ Yale $$1024$$ 0.08 0.25 $$2.11\times 10^{-15}$$ $$4096$$ 0.24 17.3 $$3.10\times 10^{-15}$$ $$10000$$ 1.13 238.5 $$4.35\times 10^{-15}$$ Extended YaleB $$1024$$ 0.08 0.27 $$2.13\times 10^{-15}$$ $$4096$$ 0.27 17.3 $$3.14\times 10^{-15}$$ $$10000$$ 1.22 238.6 $$4.50\times 10^{-15}$$ Table 1 Example $$1$$: CPU time in seconds and the relative errors for approximating $$\exp(S_B)$$ and $$\exp(S_W)$$ $${\tt Database}$$ $${\tt n}$$ $${\tt Algorithm\,4}$$ $${\tt expm.m}$$ $${\tt Rel\_ErrF}$$ ORL 1024 0.08 0.26 $$2.20\times 10^{-15}$$ $$4096$$ 0.26 17.3 $$3.22\times 10^{-15}$$ $$10304$$ 1.28 261.1 $$4.49\times 10^{-15}$$ Yale $$1024$$ 0.08 0.25 $$2.11\times 10^{-15}$$ $$4096$$ 0.24 17.3 $$3.10\times 10^{-15}$$ $$10000$$ 1.13 238.5 $$4.35\times 10^{-15}$$ Extended YaleB $$1024$$ 0.08 0.27 $$2.13\times 10^{-15}$$ $$4096$$ 0.27 17.3 $$3.14\times 10^{-15}$$ $$10000$$ 1.22 238.6 $$4.50\times 10^{-15}$$ $${\tt Database}$$ $${\tt n}$$ $${\tt Algorithm\,4}$$ $${\tt expm.m}$$ $${\tt Rel\_ErrF}$$ ORL 1024 0.08 0.26 $$2.20\times 10^{-15}$$ $$4096$$ 0.26 17.3 $$3.22\times 10^{-15}$$ $$10304$$ 1.28 261.1 $$4.49\times 10^{-15}$$ Yale $$1024$$ 0.08 0.25 $$2.11\times 10^{-15}$$ $$4096$$ 0.24 17.3 $$3.10\times 10^{-15}$$ $$10000$$ 1.13 238.5 $$4.35\times 10^{-15}$$ Extended YaleB $$1024$$ 0.08 0.27 $$2.13\times 10^{-15}$$ $$4096$$ 0.27 17.3 $$3.14\times 10^{-15}$$ $$10000$$ 1.22 238.6 $$4.50\times 10^{-15}$$ We observe from Table 1 that Algorithm 4 runs much faster than $${\tt expm.m}$$, especially when $$n$$ is large. For instance, when the dimensionality of the data sets is around $$10^4$$, $${\tt expm.m}$$ requires about 240 s, while our new method needs only about 1.2 s, which is a significant improvement. Furthermore, the relative errors of our approximations are of the order of $$\mathcal{O}(10^{-15})$$. As a result, our new method is very efficient and reliable for solving the large matrix exponential problems arising in the EDA framework. 4.2 Computing approximations to $$\varphi$$-functions of large matrices with low rank or with fast decaying singular values We show the efficiency of Algorithm 1 for consecutively approximating several $$\varphi$$-functions of $$A$$ with low rank or with fast decaying singular values. The test matrices are available from Davis & Hu (2011) and Tspras, and Table 2 lists the problem characteristics and the numerical rank, as well as the error (measured by the Frobenius norm) obtained from the SCR. Table 2 Example $$2$$: Problem characteristics of the test matrices, $$r$$, and the error from the SCR approximation Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_F$$ Description $${\tt man5976}$$ 5976 2341 $$2.83\times 10^{-9}$$ Structural problem $${\tt Movies}$$ 5757 1053 $$2.01\times 10^{-11}$$ Directed network $${\tt lock3491}$$ 3491 603 $$1.96\times 10^{-8}$$ Structural problem $${\tt cegb3306}$$ 3306 537 $$7.76\times 10^{-11}$$ Finite element framework $${\tt zenios}$$ 2873 250 $$7.75\times 10^{-6}$$ Optimization problem $${\tt watt\_1}$$ 1856 151 $$9.97\times 10^{-6}$$ Computational fluid dynamics $${\tt watt\_2}$$ 1856 151 $$9.98\times 10^{-6}$$ Computational fluid dynamics $${\tt eris1176}$$ 1176 774 $$2.49\times 10^{-10}$$ Power network problem Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_F$$ Description $${\tt man5976}$$ 5976 2341 $$2.83\times 10^{-9}$$ Structural problem $${\tt Movies}$$ 5757 1053 $$2.01\times 10^{-11}$$ Directed network $${\tt lock3491}$$ 3491 603 $$1.96\times 10^{-8}$$ Structural problem $${\tt cegb3306}$$ 3306 537 $$7.76\times 10^{-11}$$ Finite element framework $${\tt zenios}$$ 2873 250 $$7.75\times 10^{-6}$$ Optimization problem $${\tt watt\_1}$$ 1856 151 $$9.97\times 10^{-6}$$ Computational fluid dynamics $${\tt watt\_2}$$ 1856 151 $$9.98\times 10^{-6}$$ Computational fluid dynamics $${\tt eris1176}$$ 1176 774 $$2.49\times 10^{-10}$$ Power network problem Table 2 Example $$2$$: Problem characteristics of the test matrices, $$r$$, and the error from the SCR approximation Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_F$$ Description $${\tt man5976}$$ 5976 2341 $$2.83\times 10^{-9}$$ Structural problem $${\tt Movies}$$ 5757 1053 $$2.01\times 10^{-11}$$ Directed network $${\tt lock3491}$$ 3491 603 $$1.96\times 10^{-8}$$ Structural problem $${\tt cegb3306}$$ 3306 537 $$7.76\times 10^{-11}$$ Finite element framework $${\tt zenios}$$ 2873 250 $$7.75\times 10^{-6}$$ Optimization problem $${\tt watt\_1}$$ 1856 151 $$9.97\times 10^{-6}$$ Computational fluid dynamics $${\tt watt\_2}$$ 1856 151 $$9.98\times 10^{-6}$$ Computational fluid dynamics $${\tt eris1176}$$ 1176 774 $$2.49\times 10^{-10}$$ Power network problem Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_F$$ Description $${\tt man5976}$$ 5976 2341 $$2.83\times 10^{-9}$$ Structural problem $${\tt Movies}$$ 5757 1053 $$2.01\times 10^{-11}$$ Directed network $${\tt lock3491}$$ 3491 603 $$1.96\times 10^{-8}$$ Structural problem $${\tt cegb3306}$$ 3306 537 $$7.76\times 10^{-11}$$ Finite element framework $${\tt zenios}$$ 2873 250 $$7.75\times 10^{-6}$$ Optimization problem $${\tt watt\_1}$$ 1856 151 $$9.97\times 10^{-6}$$ Computational fluid dynamics $${\tt watt\_2}$$ 1856 151 $$9.98\times 10^{-6}$$ Computational fluid dynamics $${\tt eris1176}$$ 1176 774 $$2.49\times 10^{-10}$$ Power network problem In this example, we compare Algorithm 1 with $${\tt expm.m/phipade.m}$$, that is, $${\tt expm.m}$$ for the matrix exponential $$\exp(A)$$ and $${\tt phipade.m}$$ for $$\varphi_{\ell}(A),~\ell=1,2,3,4$$. In Algorithm 1, the CPU time consists of computing the SCR and the evaluation of $$\varphi_{\ell}(Z)~(\ell=1,2,3,4,5)$$ by using $${\tt phipade.m}$$, as well as forming $$\varphi_{\ell}(A)$$ in terms of (2.3), $$\ell=0,1,2,3,4$$. In $${\tt expm.m}$$/$${\tt phipade.m}$$, the CPU time consists of computing $$\varphi_{\ell}(A)$$ by using $${\tt expm.m}$$ ($$\ell=0$$) and $${\tt phipade.m}$$ ($$\ell=1,2,3,4$$). In order to measure the accuracy of the solutions, we define the maximal relative error as $$ {\tt Rel\_ErrF}=\max_{0\leq\ell\leq 4}\frac{\|\varphi_{\ell}(A)-{\varphi_{\ell}(\widetilde{A})}\|_F}{\|\varphi_{\ell}(A)\|_F}, $$ where $$\varphi_{\ell}(A)$$ is the ‘exact solution’ obtained from $${\tt expm.m}$$ for $$\ell=0$$ and $${\tt phipade.m}$$ for $$\ell=1,2,3,4$$, and $${\varphi_{\ell}(\widetilde{A})}$$ is the approximation obtained from running Algorithm 1. Table 3 lists the numerical results. Table 3 Example $$2$$: CPU time in seconds and the maximal relative error for approximating the $$\varphi_{\ell}$$-functions, $$\ell=0,1,2,3,4$$, where ‘$$^{*}$$’ denotes that we compute $$\varphi_{\ell}(-A)$$ instead of $$\varphi_{\ell}(A)$$ $${\tt Test matrix}$$ $${\tt Algorithm\,1}$$ $${\tt expm.m/phipade.m}$$ $${\tt Rel\_ErrF}$$ $${\tt man5976}$$$$^{*}$$ 205.5 10314.0 $$1.10\times 10^{-12}$$ $${\tt Movies}$$ 360.6 456.1 $$4.32\times 10^{-14}$$ $${\tt lock}$$3491$$^{*}$$ 5.49 2469.6 $$5.12\times 10^{-13}$$ $${\tt cegb}$$3306$$^{*}$$ 3.11 1599.7 $$1.24\times 10^{-13}$$ $${\tt zenios}$$ 0.63 3.63 $$1.17\times 10^{-7}$$ $${\tt watt}$$_1 0.19 30.4 $$1.93\times 10^{-7}$$ $${\tt watt}$$_2 0.19 39.1 $$2.02\times 10^{-7}$$ $${\tt eris1176}$$$$^{*}$$ 6.74 85.7 $$2.86\times 10^{-12}$$ $${\tt Test matrix}$$ $${\tt Algorithm\,1}$$ $${\tt expm.m/phipade.m}$$ $${\tt Rel\_ErrF}$$ $${\tt man5976}$$$$^{*}$$ 205.5 10314.0 $$1.10\times 10^{-12}$$ $${\tt Movies}$$ 360.6 456.1 $$4.32\times 10^{-14}$$ $${\tt lock}$$3491$$^{*}$$ 5.49 2469.6 $$5.12\times 10^{-13}$$ $${\tt cegb}$$3306$$^{*}$$ 3.11 1599.7 $$1.24\times 10^{-13}$$ $${\tt zenios}$$ 0.63 3.63 $$1.17\times 10^{-7}$$ $${\tt watt}$$_1 0.19 30.4 $$1.93\times 10^{-7}$$ $${\tt watt}$$_2 0.19 39.1 $$2.02\times 10^{-7}$$ $${\tt eris1176}$$$$^{*}$$ 6.74 85.7 $$2.86\times 10^{-12}$$ Table 3 Example $$2$$: CPU time in seconds and the maximal relative error for approximating the $$\varphi_{\ell}$$-functions, $$\ell=0,1,2,3,4$$, where ‘$$^{*}$$’ denotes that we compute $$\varphi_{\ell}(-A)$$ instead of $$\varphi_{\ell}(A)$$ $${\tt Test matrix}$$ $${\tt Algorithm\,1}$$ $${\tt expm.m/phipade.m}$$ $${\tt Rel\_ErrF}$$ $${\tt man5976}$$$$^{*}$$ 205.5 10314.0 $$1.10\times 10^{-12}$$ $${\tt Movies}$$ 360.6 456.1 $$4.32\times 10^{-14}$$ $${\tt lock}$$3491$$^{*}$$ 5.49 2469.6 $$5.12\times 10^{-13}$$ $${\tt cegb}$$3306$$^{*}$$ 3.11 1599.7 $$1.24\times 10^{-13}$$ $${\tt zenios}$$ 0.63 3.63 $$1.17\times 10^{-7}$$ $${\tt watt}$$_1 0.19 30.4 $$1.93\times 10^{-7}$$ $${\tt watt}$$_2 0.19 39.1 $$2.02\times 10^{-7}$$ $${\tt eris1176}$$$$^{*}$$ 6.74 85.7 $$2.86\times 10^{-12}$$ $${\tt Test matrix}$$ $${\tt Algorithm\,1}$$ $${\tt expm.m/phipade.m}$$ $${\tt Rel\_ErrF}$$ $${\tt man5976}$$$$^{*}$$ 205.5 10314.0 $$1.10\times 10^{-12}$$ $${\tt Movies}$$ 360.6 456.1 $$4.32\times 10^{-14}$$ $${\tt lock}$$3491$$^{*}$$ 5.49 2469.6 $$5.12\times 10^{-13}$$ $${\tt cegb}$$3306$$^{*}$$ 3.11 1599.7 $$1.24\times 10^{-13}$$ $${\tt zenios}$$ 0.63 3.63 $$1.17\times 10^{-7}$$ $${\tt watt}$$_1 0.19 30.4 $$1.93\times 10^{-7}$$ $${\tt watt}$$_2 0.19 39.1 $$2.02\times 10^{-7}$$ $${\tt eris1176}$$$$^{*}$$ 6.74 85.7 $$2.86\times 10^{-12}$$ It is seen from Table 3 that Algorithm 1 works quite well, and we provide a competitive candidate for consecutively approximating several $$\varphi$$-functions of large sparse matrices with low rank or with fast decaying singular values. Two remarks are given. First, Algorithm 1 runs much faster than $${\tt expm.m}$$/$${\tt phipade.m}$$. For instance, 205.5 s vs. 10314.0 s for the $${\tt man5976}$$ matrix, 5.49 s vs. 2469.6 s for the $${\tt lock3491}$$ matrix, and 3.11 s vs. 1599.7 s for the $${\tt cegb3306}$$ matrix. Second, the accuracy of our approximations is satisfactory in most cases. However, for the $${\tt zenios}$$, $${\tt watt\_1}$$ and $${\tt watt\_2}$$ matrices, the relative errors $${\tt Rel\_ErrF}$$ are of the order of $$\mathcal{O}(10^{-7})$$. Let us discuss it in detail. In Fig. 1, we plot the singular values of the $${\tt eris1176}$$ matrix and the $${\tt watt\_1}$$ matrix. It is observed that the $${\tt eris1176}$$ matrix has (relatively) faster decaying singular values than the $${\tt watt\_1}$$ matrix. Indeed, the errors $$\|A-\widetilde{A}\|_F$$ from the SCR approximations, with respect to the three matrices $${\tt zenios}$$, $${\tt watt\_1}$$ and $${\tt watt\_2}$$, are about $$7.75\times 10^{-6}$$, $$9.97\times 10^{-6}$$ and $$9.98\times 10^{-6}$$, respectively, while that of the $${\tt eris1176}$$ matrix is about $$2.49\times 10^{-10}$$. In Table 4, we provide the relative F-norm condition numbers for the $$\varphi_{\ell}$$-functions of the four matrices, which are computed by using the $${\tt funm\_condest\_fro.m}$$ function in the Matrix Function Toolbox (Higham, 2008b). One can see that the maximal condition numbers with respect to the four matrices are 0.98, 0.25, 1.71 and 60.0. By Theorem 2.6, the accuracy of the computed solution of the $${\tt eris1176}$$ matrix is about $$1.87\times 10^{-10}$$, while that of the $${\tt watt\_1}$$ matrix is about $$2.49\times 10^{-6}$$. Hence, the accuracy with respect to the $${\tt eris1176}$$ matrix is higher than that of the other three matrices. Further, it is possible to get a smaller $${\tt Rel\_ErrF}$$ if a larger $$r$$ is used in our new method, since the error $$\|A-\widetilde{A}\|_F$$ will be smaller. Consequently, our new method is suitable for $$\varphi$$-functions of large matrices with low rank or with fast decaying singular values. Fig. 1. View largeDownload slide Example 2: Singular values of the $${\tt eris1176}$$ and the $${\tt watt\_1}$$ matrices. Fig. 1. View largeDownload slide Example 2: Singular values of the $${\tt eris1176}$$ and the $${\tt watt\_1}$$ matrices. Table 4 Example $$2$$: Relative F-norm condition numbers for the $$\varphi_{\ell}$$-functions of the four matrices $${\tt zenios}$$, $${\tt watt\_1}$$, $${\tt watt\_2}$$ and $${\tt eris1176}$$, $$\ell=0,1,2,3,4$$ $$\ell$$ $${\tt zenios}$$ $${\tt watt\_1}$$ $${\tt watt\_2}$$ $${\tt eris1176}$$ 0 $$1.32\times 10^{-2}$$ $$1.66\times 10^{-2}$$ $$2.02\times 10^{-2}$$ 0.12 1 0.98 0.25 1.71 60.0 2 0.40 0.14 0.75 25.2 3 0.19 $$9.34\times 10^{-2}$$ 0.41 11.8 4 0.13 $$6.83\times 10^{-2}$$ 0.26 4.44 $$\ell$$ $${\tt zenios}$$ $${\tt watt\_1}$$ $${\tt watt\_2}$$ $${\tt eris1176}$$ 0 $$1.32\times 10^{-2}$$ $$1.66\times 10^{-2}$$ $$2.02\times 10^{-2}$$ 0.12 1 0.98 0.25 1.71 60.0 2 0.40 0.14 0.75 25.2 3 0.19 $$9.34\times 10^{-2}$$ 0.41 11.8 4 0.13 $$6.83\times 10^{-2}$$ 0.26 4.44 Table 4 Example $$2$$: Relative F-norm condition numbers for the $$\varphi_{\ell}$$-functions of the four matrices $${\tt zenios}$$, $${\tt watt\_1}$$, $${\tt watt\_2}$$ and $${\tt eris1176}$$, $$\ell=0,1,2,3,4$$ $$\ell$$ $${\tt zenios}$$ $${\tt watt\_1}$$ $${\tt watt\_2}$$ $${\tt eris1176}$$ 0 $$1.32\times 10^{-2}$$ $$1.66\times 10^{-2}$$ $$2.02\times 10^{-2}$$ 0.12 1 0.98 0.25 1.71 60.0 2 0.40 0.14 0.75 25.2 3 0.19 $$9.34\times 10^{-2}$$ 0.41 11.8 4 0.13 $$6.83\times 10^{-2}$$ 0.26 4.44 $$\ell$$ $${\tt zenios}$$ $${\tt watt\_1}$$ $${\tt watt\_2}$$ $${\tt eris1176}$$ 0 $$1.32\times 10^{-2}$$ $$1.66\times 10^{-2}$$ $$2.02\times 10^{-2}$$ 0.12 1 0.98 0.25 1.71 60.0 2 0.40 0.14 0.75 25.2 3 0.19 $$9.34\times 10^{-2}$$ 0.41 11.8 4 0.13 $$6.83\times 10^{-2}$$ 0.26 4.44 4.3 Estimating the relative and the absolute condition numbers This example aims to demonstrate the efficiency of our new strategy for estimating the absolute and the relative 2-norm condition numbers of $$\varphi$$-functions. There are six test matrices that are available from Davis & Hu (2011) and Tspras. For the first three matrices, we compute the absolute and the relative condition numbers for a given $$\ell$$, while for the rest of the three matrices, we evaluate the condition numbers consecutively for $$\ell=0,1,2,3,4$$. Table 5 lists the problem characteristics the numerical rank, as well as the error (measured by the 2-norm) obtained from the SCR of these matrices. Table 5 Example $$3$$: Problem characteristics of the test matrices, $$r$$, and the error from the SCR approximation Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_2$$ $${\tt EVA}$$ 8497 1301 $$1.64\times 10^{-13}$$ $${\tt EPA}$$ 4772 951 $$3.57\times 10^{-12}$$ $${\tt eris1176}$$ 1176 774 $$1.85\times 10^{-10}$$ $${\tt Abortion}$$ 2293 633 $$9.41\times 10^{-12}$$ $${\tt Computational\_Geometry}$$ 1226 415 $$1.96\times 10^{-11}$$ $${\tt Computational\_Complexity}$$ 884 254 $$7.84\times 10^{-14}$$ Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_2$$ $${\tt EVA}$$ 8497 1301 $$1.64\times 10^{-13}$$ $${\tt EPA}$$ 4772 951 $$3.57\times 10^{-12}$$ $${\tt eris1176}$$ 1176 774 $$1.85\times 10^{-10}$$ $${\tt Abortion}$$ 2293 633 $$9.41\times 10^{-12}$$ $${\tt Computational\_Geometry}$$ 1226 415 $$1.96\times 10^{-11}$$ $${\tt Computational\_Complexity}$$ 884 254 $$7.84\times 10^{-14}$$ Table 5 Example $$3$$: Problem characteristics of the test matrices, $$r$$, and the error from the SCR approximation Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_2$$ $${\tt EVA}$$ 8497 1301 $$1.64\times 10^{-13}$$ $${\tt EPA}$$ 4772 951 $$3.57\times 10^{-12}$$ $${\tt eris1176}$$ 1176 774 $$1.85\times 10^{-10}$$ $${\tt Abortion}$$ 2293 633 $$9.41\times 10^{-12}$$ $${\tt Computational\_Geometry}$$ 1226 415 $$1.96\times 10^{-11}$$ $${\tt Computational\_Complexity}$$ 884 254 $$7.84\times 10^{-14}$$ Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_2$$ $${\tt EVA}$$ 8497 1301 $$1.64\times 10^{-13}$$ $${\tt EPA}$$ 4772 951 $$3.57\times 10^{-12}$$ $${\tt eris1176}$$ 1176 774 $$1.85\times 10^{-10}$$ $${\tt Abortion}$$ 2293 633 $$9.41\times 10^{-12}$$ $${\tt Computational\_Geometry}$$ 1226 415 $$1.96\times 10^{-11}$$ $${\tt Computational\_Complexity}$$ 884 254 $$7.84\times 10^{-14}$$ We compare Algorithm 2 with the $${\tt funm\_condest1.m}$$ function in the Matrix Function Toolbox (Higham, 2008b). In Algorithm 2, the matrices $$R_1,R_2$$ are obtained from the QR decompositions of $$X$$ and $$Y$$, respectively, by using the MATLAB built-in function $${\tt qr.m}$$. The matrix $$L_{\varphi_{\ell+1}}(Z)$$ in (3.10) is evaluated by using the $${\tt funm\_condest\_fro.m}$$ function, and $$\|A\|_2$$ is evaluated by using the MATLAB built-in function $${\tt svds.m}$$. When $$\ell=0$$, the parameter ‘$${\tt fun}$$’ in $${\tt funm\_condest1.m}$$ is called by the MATLAB built-in function $${\tt expm.m}$$, whereas when $$\ell>0$$ this parameter is called by the $${\tt phipade.m}$$ function in the EXPINT package. The CPU time for Algorithm 2 is composed of computing the SCR and estimating $$\|A\|_2$$, as well as calculating (3.7) and (3.10). Tables 6–11 report the numerical results, where ‘$${\tt Relative\_Est}$$’ and ‘$${\tt Absolute\_Est}$$’ denote estimations to the relative and the absolute condition numbers, respectively. Table 6 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$ and the CPU time in seconds where ‘$$>$$3h’ denotes that the algorithm fails to converge within 3 h. The $${\tt EVA}$$ matrix, $$n=8497, r=1303$$ $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$7.74\times 10^3$$ $$1.83\times 10^{4}$$ 700.0 $${\tt Algorithm\,2}$$ $$1.07\times 10^3$$ $$1.32\times 10^{3}$$ 420.1 1 $${\tt funm\_condest1}$$ $$7.07\times 10^3$$ $$7.22\times 10^{3}$$ $$5846.1$$ $${\tt Algorithm\,2}$$ $$538.7$$ $$269.6$$ 425.2 2 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$321.7$$ 53.5 429.8 3 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$122.5$$ 5.08 440.9 4 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ 73.3 0.61 447.1 $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$7.74\times 10^3$$ $$1.83\times 10^{4}$$ 700.0 $${\tt Algorithm\,2}$$ $$1.07\times 10^3$$ $$1.32\times 10^{3}$$ 420.1 1 $${\tt funm\_condest1}$$ $$7.07\times 10^3$$ $$7.22\times 10^{3}$$ $$5846.1$$ $${\tt Algorithm\,2}$$ $$538.7$$ $$269.6$$ 425.2 2 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$321.7$$ 53.5 429.8 3 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$122.5$$ 5.08 440.9 4 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ 73.3 0.61 447.1 Table 6 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$ and the CPU time in seconds where ‘$$>$$3h’ denotes that the algorithm fails to converge within 3 h. The $${\tt EVA}$$ matrix, $$n=8497, r=1303$$ $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$7.74\times 10^3$$ $$1.83\times 10^{4}$$ 700.0 $${\tt Algorithm\,2}$$ $$1.07\times 10^3$$ $$1.32\times 10^{3}$$ 420.1 1 $${\tt funm\_condest1}$$ $$7.07\times 10^3$$ $$7.22\times 10^{3}$$ $$5846.1$$ $${\tt Algorithm\,2}$$ $$538.7$$ $$269.6$$ 425.2 2 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$321.7$$ 53.5 429.8 3 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$122.5$$ 5.08 440.9 4 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ 73.3 0.61 447.1 $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$7.74\times 10^3$$ $$1.83\times 10^{4}$$ 700.0 $${\tt Algorithm\,2}$$ $$1.07\times 10^3$$ $$1.32\times 10^{3}$$ 420.1 1 $${\tt funm\_condest1}$$ $$7.07\times 10^3$$ $$7.22\times 10^{3}$$ $$5846.1$$ $${\tt Algorithm\,2}$$ $$538.7$$ $$269.6$$ 425.2 2 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$321.7$$ 53.5 429.8 3 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$122.5$$ 5.08 440.9 4 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ 73.3 0.61 447.1 Table 7 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A)$$, $$\ell=0,1,2,3,4$$ and the CPU time in seconds. The $${\tt EPA}$$ matrix, $$n=4772, r=951$$ $$\ell$$ Method Rel_Est Abs_Est CPU 0 $${\tt funm\_condest1}$$ $$9.90\times 10^{3}$$ $$2.60\times 10^{4}$$ 90.1 $${\tt Algorithm\,2}$$ $$3.34\times 10^4$$ $$1.52\times 10^{5}$$ 95.4 1 $${\tt funm\_condest1}$$ $$1.09\times 10^4$$ $$8.05\times 10^{3}$$ 272.1 $${\tt Algorithm\,2}$$ $$2.25\times 10^4$$ $$2.74\times 10^{4}$$ 100.6 2 $${\tt funm\_condest1}$$ $$8.95\times 10^3$$ $$2.03\times 10^{3}$$ 443.7 $${\tt Algorithm\,2}$$ $$1.47\times 10^4$$ $$4.27\times 10^3$$ 105.5 3 $${\tt funm\_condest1}$$ $$7.86\times 10^3$$ $$425.9$$ 773.0 $${\tt Algorithm\,2}$$ $$9.71\times 10^3$$ $$587.4$$ 110.7 4 $${\tt funm\_condest1}$$ $$7.17\times 10^3$$ 75.2 1357.7 $${\tt Algorithm\,2}$$ $$5.29\times 10^3$$ 57.6 116.7 $$\ell$$ Method Rel_Est Abs_Est CPU 0 $${\tt funm\_condest1}$$ $$9.90\times 10^{3}$$ $$2.60\times 10^{4}$$ 90.1 $${\tt Algorithm\,2}$$ $$3.34\times 10^4$$ $$1.52\times 10^{5}$$ 95.4 1 $${\tt funm\_condest1}$$ $$1.09\times 10^4$$ $$8.05\times 10^{3}$$ 272.1 $${\tt Algorithm\,2}$$ $$2.25\times 10^4$$ $$2.74\times 10^{4}$$ 100.6 2 $${\tt funm\_condest1}$$ $$8.95\times 10^3$$ $$2.03\times 10^{3}$$ 443.7 $${\tt Algorithm\,2}$$ $$1.47\times 10^4$$ $$4.27\times 10^3$$ 105.5 3 $${\tt funm\_condest1}$$ $$7.86\times 10^3$$ $$425.9$$ 773.0 $${\tt Algorithm\,2}$$ $$9.71\times 10^3$$ $$587.4$$ 110.7 4 $${\tt funm\_condest1}$$ $$7.17\times 10^3$$ 75.2 1357.7 $${\tt Algorithm\,2}$$ $$5.29\times 10^3$$ 57.6 116.7 Table 7 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A)$$, $$\ell=0,1,2,3,4$$ and the CPU time in seconds. The $${\tt EPA}$$ matrix, $$n=4772, r=951$$ $$\ell$$ Method Rel_Est Abs_Est CPU 0 $${\tt funm\_condest1}$$ $$9.90\times 10^{3}$$ $$2.60\times 10^{4}$$ 90.1 $${\tt Algorithm\,2}$$ $$3.34\times 10^4$$ $$1.52\times 10^{5}$$ 95.4 1 $${\tt funm\_condest1}$$ $$1.09\times 10^4$$ $$8.05\times 10^{3}$$ 272.1 $${\tt Algorithm\,2}$$ $$2.25\times 10^4$$ $$2.74\times 10^{4}$$ 100.6 2 $${\tt funm\_condest1}$$ $$8.95\times 10^3$$ $$2.03\times 10^{3}$$ 443.7 $${\tt Algorithm\,2}$$ $$1.47\times 10^4$$ $$4.27\times 10^3$$ 105.5 3 $${\tt funm\_condest1}$$ $$7.86\times 10^3$$ $$425.9$$ 773.0 $${\tt Algorithm\,2}$$ $$9.71\times 10^3$$ $$587.4$$ 110.7 4 $${\tt funm\_condest1}$$ $$7.17\times 10^3$$ 75.2 1357.7 $${\tt Algorithm\,2}$$ $$5.29\times 10^3$$ 57.6 116.7 $$\ell$$ Method Rel_Est Abs_Est CPU 0 $${\tt funm\_condest1}$$ $$9.90\times 10^{3}$$ $$2.60\times 10^{4}$$ 90.1 $${\tt Algorithm\,2}$$ $$3.34\times 10^4$$ $$1.52\times 10^{5}$$ 95.4 1 $${\tt funm\_condest1}$$ $$1.09\times 10^4$$ $$8.05\times 10^{3}$$ 272.1 $${\tt Algorithm\,2}$$ $$2.25\times 10^4$$ $$2.74\times 10^{4}$$ 100.6 2 $${\tt funm\_condest1}$$ $$8.95\times 10^3$$ $$2.03\times 10^{3}$$ 443.7 $${\tt Algorithm\,2}$$ $$1.47\times 10^4$$ $$4.27\times 10^3$$ 105.5 3 $${\tt funm\_condest1}$$ $$7.86\times 10^3$$ $$425.9$$ 773.0 $${\tt Algorithm\,2}$$ $$9.71\times 10^3$$ $$587.4$$ 110.7 4 $${\tt funm\_condest1}$$ $$7.17\times 10^3$$ 75.2 1357.7 $${\tt Algorithm\,2}$$ $$5.29\times 10^3$$ 57.6 116.7 Table 8 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(-A), \ell=0,1,2,3,4$$, and the CPU time in seconds. The $${\tt eris1176}$$ matrix, $$n=1176$$, $$r=774$$ $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$575.1$$ $$1.45\times 10^{3}$$ 2.47 $${\tt Algorithm\,2}$$ $$2.72\times 10^3$$ $$4.66\times 10^{3}$$ 11.1 1 $${\tt funm\_condest1}$$ $$1.09\times 10^3$$ $$548.7$$ 38.2 $${\tt Algorithm\,2}$$ $$2.66\times 10^3$$ $$890.5$$ 14.5 2 $${\tt funm\_condest1}$$ $$1.75\times 10^3$$ $$171.0$$ 60.2 $${\tt Algorithm\,2}$$ $$2.71\times 10^3$$ $$167.6$$ 17.6 3 $${\tt funm\_condest1}$$ $$2.37\times 10^3$$ $$42.2$$ 82.2 $${\tt Algorithm\,2}$$ $$2.82\times 10^3$$ 29.4 20.5 4 $${\tt funm\_condest1}$$ $$2.85\times 10^3$$ 8.49 103.4 $${\tt Algorithm\,2}$$ $$2.88\times 10^3$$ 4.61 24.2 $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$575.1$$ $$1.45\times 10^{3}$$ 2.47 $${\tt Algorithm\,2}$$ $$2.72\times 10^3$$ $$4.66\times 10^{3}$$ 11.1 1 $${\tt funm\_condest1}$$ $$1.09\times 10^3$$ $$548.7$$ 38.2 $${\tt Algorithm\,2}$$ $$2.66\times 10^3$$ $$890.5$$ 14.5 2 $${\tt funm\_condest1}$$ $$1.75\times 10^3$$ $$171.0$$ 60.2 $${\tt Algorithm\,2}$$ $$2.71\times 10^3$$ $$167.6$$ 17.6 3 $${\tt funm\_condest1}$$ $$2.37\times 10^3$$ $$42.2$$ 82.2 $${\tt Algorithm\,2}$$ $$2.82\times 10^3$$ 29.4 20.5 4 $${\tt funm\_condest1}$$ $$2.85\times 10^3$$ 8.49 103.4 $${\tt Algorithm\,2}$$ $$2.88\times 10^3$$ 4.61 24.2 Table 8 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(-A), \ell=0,1,2,3,4$$, and the CPU time in seconds. The $${\tt eris1176}$$ matrix, $$n=1176$$, $$r=774$$ $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$575.1$$ $$1.45\times 10^{3}$$ 2.47 $${\tt Algorithm\,2}$$ $$2.72\times 10^3$$ $$4.66\times 10^{3}$$ 11.1 1 $${\tt funm\_condest1}$$ $$1.09\times 10^3$$ $$548.7$$ 38.2 $${\tt Algorithm\,2}$$ $$2.66\times 10^3$$ $$890.5$$ 14.5 2 $${\tt funm\_condest1}$$ $$1.75\times 10^3$$ $$171.0$$ 60.2 $${\tt Algorithm\,2}$$ $$2.71\times 10^3$$ $$167.6$$ 17.6 3 $${\tt funm\_condest1}$$ $$2.37\times 10^3$$ $$42.2$$ 82.2 $${\tt Algorithm\,2}$$ $$2.82\times 10^3$$ 29.4 20.5 4 $${\tt funm\_condest1}$$ $$2.85\times 10^3$$ 8.49 103.4 $${\tt Algorithm\,2}$$ $$2.88\times 10^3$$ 4.61 24.2 $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$575.1$$ $$1.45\times 10^{3}$$ 2.47 $${\tt Algorithm\,2}$$ $$2.72\times 10^3$$ $$4.66\times 10^{3}$$ 11.1 1 $${\tt funm\_condest1}$$ $$1.09\times 10^3$$ $$548.7$$ 38.2 $${\tt Algorithm\,2}$$ $$2.66\times 10^3$$ $$890.5$$ 14.5 2 $${\tt funm\_condest1}$$ $$1.75\times 10^3$$ $$171.0$$ 60.2 $${\tt Algorithm\,2}$$ $$2.71\times 10^3$$ $$167.6$$ 17.6 3 $${\tt funm\_condest1}$$ $$2.37\times 10^3$$ $$42.2$$ 82.2 $${\tt Algorithm\,2}$$ $$2.82\times 10^3$$ 29.4 20.5 4 $${\tt funm\_condest1}$$ $$2.85\times 10^3$$ 8.49 103.4 $${\tt Algorithm\,2}$$ $$2.88\times 10^3$$ 4.61 24.2 Table 9 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$; the $${\tt Abortion}$$ matrix, $$n=2293, r=633$$. In the new algorithm, we compute the condition numbers consecutively, while in $${\tt funm\_condest1}$$ we have to compute them one by one. The CPU timings of Algorithm 2 and $${\tt funm\_condest1}$$ are $${104.4}$$ and $$336.9$$ s, respectively $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$1.37\times 10^3$$ $$9.87\times 10^{4}$$ $${\tt Algorithm\,2}$$ 194.4 $$2.90\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$1.93\times 10^3$$ $$1.77\times 10^{4}$$ $${\tt Algorithm\,2}$$ 170.8 $$3.24\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$2.72\times 10^3$$ $$3.14\times 10^3$$ $${\tt Algorithm\,2}$$ 144.0 343.2 3 $${\tt funm\_condest1}$$ $$3.64\times 10^3$$ 522.9 $${\tt Algorithm\,2}$$ 135.2 39.8 4 $${\tt funm\_condest1}$$ $$4.58\times 10^3$$ 79.1 $${\tt Algorithm\,2}$$ 130.7 4.58 $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$1.37\times 10^3$$ $$9.87\times 10^{4}$$ $${\tt Algorithm\,2}$$ 194.4 $$2.90\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$1.93\times 10^3$$ $$1.77\times 10^{4}$$ $${\tt Algorithm\,2}$$ 170.8 $$3.24\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$2.72\times 10^3$$ $$3.14\times 10^3$$ $${\tt Algorithm\,2}$$ 144.0 343.2 3 $${\tt funm\_condest1}$$ $$3.64\times 10^3$$ 522.9 $${\tt Algorithm\,2}$$ 135.2 39.8 4 $${\tt funm\_condest1}$$ $$4.58\times 10^3$$ 79.1 $${\tt Algorithm\,2}$$ 130.7 4.58 Table 9 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$; the $${\tt Abortion}$$ matrix, $$n=2293, r=633$$. In the new algorithm, we compute the condition numbers consecutively, while in $${\tt funm\_condest1}$$ we have to compute them one by one. The CPU timings of Algorithm 2 and $${\tt funm\_condest1}$$ are $${104.4}$$ and $$336.9$$ s, respectively $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$1.37\times 10^3$$ $$9.87\times 10^{4}$$ $${\tt Algorithm\,2}$$ 194.4 $$2.90\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$1.93\times 10^3$$ $$1.77\times 10^{4}$$ $${\tt Algorithm\,2}$$ 170.8 $$3.24\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$2.72\times 10^3$$ $$3.14\times 10^3$$ $${\tt Algorithm\,2}$$ 144.0 343.2 3 $${\tt funm\_condest1}$$ $$3.64\times 10^3$$ 522.9 $${\tt Algorithm\,2}$$ 135.2 39.8 4 $${\tt funm\_condest1}$$ $$4.58\times 10^3$$ 79.1 $${\tt Algorithm\,2}$$ 130.7 4.58 $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$1.37\times 10^3$$ $$9.87\times 10^{4}$$ $${\tt Algorithm\,2}$$ 194.4 $$2.90\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$1.93\times 10^3$$ $$1.77\times 10^{4}$$ $${\tt Algorithm\,2}$$ 170.8 $$3.24\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$2.72\times 10^3$$ $$3.14\times 10^3$$ $${\tt Algorithm\,2}$$ 144.0 343.2 3 $${\tt funm\_condest1}$$ $$3.64\times 10^3$$ 522.9 $${\tt Algorithm\,2}$$ 135.2 39.8 4 $${\tt funm\_condest1}$$ $$4.58\times 10^3$$ 79.1 $${\tt Algorithm\,2}$$ 130.7 4.58 Table 10 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$; the $${\tt Computational\_Geometry}$$ matrix, $$n=1226, r=415$$. In the new algorithm, we compute the condition numbers $${\tt consecutively}$$, while in $${\tt funm\_condest1}$$ we have to compute them one by one. The CPU timings of Algorithm 2 and $${\tt funm\_condest1}$$ are $${22.8}$$ and $$66.2$$ s, respectively $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$5.35\times 10^3$$ $$2.73\times 10^{4}$$ $${\tt Algorithm\,2}$$ $$3.15\times 10^3$$ $$2.35\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$4.00\times 10^3$$ $$6.33\times 10^{3}$$ $${\tt Algorithm\,2}$$ 998.6 $$2.07\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$3.20\times 10^3$$ $$1.29\times 10^3$$ $${\tt Algorithm\,2}$$ 860.8 420.1 3 $${\tt funm\_condest1}$$ $$2.69\times 10^3$$ 229.5 $${\tt Algorithm\,2}$$ 824.9 80.4 4 $${\tt funm\_condest1}$$ $$2.36\times 10^3$$ 36.0 $${\tt Algorithm\,2}$$ 262.1 4.38 $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$5.35\times 10^3$$ $$2.73\times 10^{4}$$ $${\tt Algorithm\,2}$$ $$3.15\times 10^3$$ $$2.35\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$4.00\times 10^3$$ $$6.33\times 10^{3}$$ $${\tt Algorithm\,2}$$ 998.6 $$2.07\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$3.20\times 10^3$$ $$1.29\times 10^3$$ $${\tt Algorithm\,2}$$ 860.8 420.1 3 $${\tt funm\_condest1}$$ $$2.69\times 10^3$$ 229.5 $${\tt Algorithm\,2}$$ 824.9 80.4 4 $${\tt funm\_condest1}$$ $$2.36\times 10^3$$ 36.0 $${\tt Algorithm\,2}$$ 262.1 4.38 Table 10 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$; the $${\tt Computational\_Geometry}$$ matrix, $$n=1226, r=415$$. In the new algorithm, we compute the condition numbers $${\tt consecutively}$$, while in $${\tt funm\_condest1}$$ we have to compute them one by one. The CPU timings of Algorithm 2 and $${\tt funm\_condest1}$$ are $${22.8}$$ and $$66.2$$ s, respectively $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$5.35\times 10^3$$ $$2.73\times 10^{4}$$ $${\tt Algorithm\,2}$$ $$3.15\times 10^3$$ $$2.35\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$4.00\times 10^3$$ $$6.33\times 10^{3}$$ $${\tt Algorithm\,2}$$ 998.6 $$2.07\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$3.20\times 10^3$$ $$1.29\times 10^3$$ $${\tt Algorithm\,2}$$ 860.8 420.1 3 $${\tt funm\_condest1}$$ $$2.69\times 10^3$$ 229.5 $${\tt Algorithm\,2}$$ 824.9 80.4 4 $${\tt funm\_condest1}$$ $$2.36\times 10^3$$ 36.0 $${\tt Algorithm\,2}$$ 262.1 4.38 $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$5.35\times 10^3$$ $$2.73\times 10^{4}$$ $${\tt Algorithm\,2}$$ $$3.15\times 10^3$$ $$2.35\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$4.00\times 10^3$$ $$6.33\times 10^{3}$$ $${\tt Algorithm\,2}$$ 998.6 $$2.07\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$3.20\times 10^3$$ $$1.29\times 10^3$$ $${\tt Algorithm\,2}$$ 860.8 420.1 3 $${\tt funm\_condest1}$$ $$2.69\times 10^3$$ 229.5 $${\tt Algorithm\,2}$$ 824.9 80.4 4 $${\tt funm\_condest1}$$ $$2.36\times 10^3$$ 36.0 $${\tt Algorithm\,2}$$ 262.1 4.38 Table 11 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$; the $${\tt Computational\_Complexity}$$ matrix, $$n=884, r=254$$. In the new algorithm, we compute the condition numbers $${\tt consecutively}$$, while in $${\tt funm\_condest1}$$ we have to compute them one by one. The CPU timings of Algorithm 2 and $${\tt funm\_condest1}$$ are $${6.24}$$ and $$27.5$$ s, respectively $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$2.92\times 10^3$$ $$5.27\times 10^{3}$$ $${\tt Algorithm\,2}$$ 873.9 $$2.17\times 10^{3}$$ 1 $${\tt funm\_condest1}$$ $$2.52\times 10^3$$ $$1.60\times 10^{3}$$ $${\tt Algorithm\,2}$$ 530.1 $$438.4$$ 2 $${\tt funm\_condest1}$$ $$2.04\times 10^3$$ 397.0 $${\tt Algorithm\,2}$$ 331.8 73.8 3 $${\tt funm\_condest1}$$ $$1.70\times 10^3$$ 82.0 $${\tt Algorithm\,2}$$ 211.1 10.8 4 $${\tt funm\_condest1}$$ $$1.49\times 10^3$$ 14.4 $${\tt Algorithm\,2}$$ 140.9 1.37 $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$2.92\times 10^3$$ $$5.27\times 10^{3}$$ $${\tt Algorithm\,2}$$ 873.9 $$2.17\times 10^{3}$$ 1 $${\tt funm\_condest1}$$ $$2.52\times 10^3$$ $$1.60\times 10^{3}$$ $${\tt Algorithm\,2}$$ 530.1 $$438.4$$ 2 $${\tt funm\_condest1}$$ $$2.04\times 10^3$$ 397.0 $${\tt Algorithm\,2}$$ 331.8 73.8 3 $${\tt funm\_condest1}$$ $$1.70\times 10^3$$ 82.0 $${\tt Algorithm\,2}$$ 211.1 10.8 4 $${\tt funm\_condest1}$$ $$1.49\times 10^3$$ 14.4 $${\tt Algorithm\,2}$$ 140.9 1.37 Table 11 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$; the $${\tt Computational\_Complexity}$$ matrix, $$n=884, r=254$$. In the new algorithm, we compute the condition numbers $${\tt consecutively}$$, while in $${\tt funm\_condest1}$$ we have to compute them one by one. The CPU timings of Algorithm 2 and $${\tt funm\_condest1}$$ are $${6.24}$$ and $$27.5$$ s, respectively $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$2.92\times 10^3$$ $$5.27\times 10^{3}$$ $${\tt Algorithm\,2}$$ 873.9 $$2.17\times 10^{3}$$ 1 $${\tt funm\_condest1}$$ $$2.52\times 10^3$$ $$1.60\times 10^{3}$$ $${\tt Algorithm\,2}$$ 530.1 $$438.4$$ 2 $${\tt funm\_condest1}$$ $$2.04\times 10^3$$ 397.0 $${\tt Algorithm\,2}$$ 331.8 73.8 3 $${\tt funm\_condest1}$$ $$1.70\times 10^3$$ 82.0 $${\tt Algorithm\,2}$$ 211.1 10.8 4 $${\tt funm\_condest1}$$ $$1.49\times 10^3$$ 14.4 $${\tt Algorithm\,2}$$ 140.9 1.37 $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$2.92\times 10^3$$ $$5.27\times 10^{3}$$ $${\tt Algorithm\,2}$$ 873.9 $$2.17\times 10^{3}$$ 1 $${\tt funm\_condest1}$$ $$2.52\times 10^3$$ $$1.60\times 10^{3}$$ $${\tt Algorithm\,2}$$ 530.1 $$438.4$$ 2 $${\tt funm\_condest1}$$ $$2.04\times 10^3$$ 397.0 $${\tt Algorithm\,2}$$ 331.8 73.8 3 $${\tt funm\_condest1}$$ $$1.70\times 10^3$$ 82.0 $${\tt Algorithm\,2}$$ 211.1 10.8 4 $${\tt funm\_condest1}$$ $$1.49\times 10^3$$ 14.4 $${\tt Algorithm\,2}$$ 140.9 1.37 As was pointed out by Higham (2008a, p. 64), for the absolute and the relative condition numbers, ‘what is needed is an estimate that is of the correct order of magnitude in practice—more than one correct significant digit is not needed’. Recall that the $${\tt funm\_condest1.m}$$ function estimates the 1-norm relative and absolute condition numbers, while Algorithm 2 estimates the 2-norm condition numbers. Compared with the numerical results of $${\tt funm\_condest1}$$, we see from Tables 6–11 that our new strategy captures the correct order of magnitude of the condition numbers in many cases. The difference from $${\tt funm\_condest1}$$ can partially be explained by the change of norm, and the condition numbers presented differ in (at most) up to one digit. On the other hand, Algorithm 2 often runs much faster than $${\tt funm\_condest1}$$, especially when $$\ell\geq 1$$. For instance, for the $${\tt EVA}$$ matrix, the new algorithm used 425.2 s while $${\tt funm\_condest1}$$ used 5846.1 s for $$\ell=1$$. Specifically, for $$\ell\geq 2$$, $${\tt funm\_condest1}$$ fails to converge within 3 h, and Algorithm 2 runs quite well. Moreover, we see from Tables 9 to 11 that Algorithm 2 can compute the condition numbers consecutively, whereas $${\tt funm\_condest1}$$ has to compute them one by one. So we benefit from our new strategy, and we provide a competitive algorithm for estimating the relative and the absolute condition numbers of $$\varphi$$-functions. 4.4 Sharpness of Theorem 2.6 In this example, we aim to show the sharpness of our error analysis for approximating $$\varphi$$-functions. The test matrix is the $${\tt watt\_1}$$ matrix that was used in Example 2. define $$ {\tt Abs\_Err2}=\|\varphi_{\ell}(A)-{\varphi_{\ell}(\widetilde{A})}\|_2, $$ and $$ {\tt Rel\_Err2}=\frac{\|\varphi_{\ell}(A)-{\varphi_{\ell}(\widetilde{A})}\|_2}{\|\varphi_{\ell}(A)\|_2}, $$ the ‘exact’ absolute and the relative errors of the computed solutions $${\varphi_{\ell}(\widetilde{A})}$$ with respect to $$\varphi_{\ell}(A)$$ in terms of the 2-norm. In our estimations, the values of $${\rm cond}_{\rm abs}(\varphi_{\ell},A)$$ and $${\rm cond}_{\rm rel}(\varphi_{\ell},A)$$ in the upper bounds of (2.13) and (2.14) are approximated by using (3.7) and (3.10). The corresponding estimations are denoted by ‘(2.13)–NewStr’ and ‘(2.14)–NewStr’, respectively. Table 12 lists the numerical results, as well as the estimated absolute and relative condition numbers $${\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})$$, $${\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})$$. Table 12 Example $$4$$: Absolute and relative errors and their estimations in terms of the $$2$$-norm, $$\ell=0,1,2,3,4$$. The $${\tt watt\_1}$$ matrix, with $$\|A-\widetilde{A}\|_2\approx 1.07\times 10^{-6}$$. Here eqrefeq2.16-NewStr, (2.14)-NewStr denote that the absolute and the relative condition numbers $${\rm cond}_{\rm abs}(\varphi_{\ell},A)$$ and $${\rm cond}_{\rm rel}(\varphi_{\ell},A)$$ in the upper bounds of (2.13) and (2.14) are estimated by using (3.7) and (3.10), respectively $$\ell=0$$ $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $${\tt Abs_Err2}$$ $$1.07\times 10^{-6}$$ $$5.35\times 10^{-7}$$ $$1.78\times 10^{-7}$$ $$4.45\times 10^{-8}$$ $$8.91\times 10^{-9}$$ $${\tt Rel_Err2}$$ $$3.93\times 10^{-7}$$ $$3.11\times 10^{-7}$$ $$2.48\times 10^{-7}$$ $$2.04\times 10^{-7}$$ $$1.73\times 10^{-7}$$ (2.13)–NewStr $$1.83\times 10^{-7}$$ $$5.24\times 10^{-8}$$ $$1.19\times 10^{-8}$$ $$2.22\times 10^{-9}$$ $$3.47\times 10^{-10}$$ (2.14)–NewStr $$1.06\times 10^{-7}$$ $$7.30\times 10^{-8}$$ $$5.43\times 10^{-8}$$ $$4.30\times 10^{-8}$$ $$3.49\times 10^{-8}$$ $${\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ 0.17 $$6.73\times 10^{-2}$$ $$5.01\times 10^{-2}$$ $$3.92\times 10^{-2}$$ $$3.23\times 10^{-2}$$ $${\rm cond}_{\rm rel}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ $$9.81\times 10^{-2}$$ $$4.83\times 10^{-2}$$ $$1.09\times 10^{-2}$$ $$2.00\times 10^{-3}$$ $$3.21\times 10^{-4}$$ $$\ell=0$$ $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $${\tt Abs_Err2}$$ $$1.07\times 10^{-6}$$ $$5.35\times 10^{-7}$$ $$1.78\times 10^{-7}$$ $$4.45\times 10^{-8}$$ $$8.91\times 10^{-9}$$ $${\tt Rel_Err2}$$ $$3.93\times 10^{-7}$$ $$3.11\times 10^{-7}$$ $$2.48\times 10^{-7}$$ $$2.04\times 10^{-7}$$ $$1.73\times 10^{-7}$$ (2.13)–NewStr $$1.83\times 10^{-7}$$ $$5.24\times 10^{-8}$$ $$1.19\times 10^{-8}$$ $$2.22\times 10^{-9}$$ $$3.47\times 10^{-10}$$ (2.14)–NewStr $$1.06\times 10^{-7}$$ $$7.30\times 10^{-8}$$ $$5.43\times 10^{-8}$$ $$4.30\times 10^{-8}$$ $$3.49\times 10^{-8}$$ $${\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ 0.17 $$6.73\times 10^{-2}$$ $$5.01\times 10^{-2}$$ $$3.92\times 10^{-2}$$ $$3.23\times 10^{-2}$$ $${\rm cond}_{\rm rel}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ $$9.81\times 10^{-2}$$ $$4.83\times 10^{-2}$$ $$1.09\times 10^{-2}$$ $$2.00\times 10^{-3}$$ $$3.21\times 10^{-4}$$ Table 12 Example $$4$$: Absolute and relative errors and their estimations in terms of the $$2$$-norm, $$\ell=0,1,2,3,4$$. The $${\tt watt\_1}$$ matrix, with $$\|A-\widetilde{A}\|_2\approx 1.07\times 10^{-6}$$. Here eqrefeq2.16-NewStr, (2.14)-NewStr denote that the absolute and the relative condition numbers $${\rm cond}_{\rm abs}(\varphi_{\ell},A)$$ and $${\rm cond}_{\rm rel}(\varphi_{\ell},A)$$ in the upper bounds of (2.13) and (2.14) are estimated by using (3.7) and (3.10), respectively $$\ell=0$$ $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $${\tt Abs_Err2}$$ $$1.07\times 10^{-6}$$ $$5.35\times 10^{-7}$$ $$1.78\times 10^{-7}$$ $$4.45\times 10^{-8}$$ $$8.91\times 10^{-9}$$ $${\tt Rel_Err2}$$ $$3.93\times 10^{-7}$$ $$3.11\times 10^{-7}$$ $$2.48\times 10^{-7}$$ $$2.04\times 10^{-7}$$ $$1.73\times 10^{-7}$$ (2.13)–NewStr $$1.83\times 10^{-7}$$ $$5.24\times 10^{-8}$$ $$1.19\times 10^{-8}$$ $$2.22\times 10^{-9}$$ $$3.47\times 10^{-10}$$ (2.14)–NewStr $$1.06\times 10^{-7}$$ $$7.30\times 10^{-8}$$ $$5.43\times 10^{-8}$$ $$4.30\times 10^{-8}$$ $$3.49\times 10^{-8}$$ $${\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ 0.17 $$6.73\times 10^{-2}$$ $$5.01\times 10^{-2}$$ $$3.92\times 10^{-2}$$ $$3.23\times 10^{-2}$$ $${\rm cond}_{\rm rel}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ $$9.81\times 10^{-2}$$ $$4.83\times 10^{-2}$$ $$1.09\times 10^{-2}$$ $$2.00\times 10^{-3}$$ $$3.21\times 10^{-4}$$ $$\ell=0$$ $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $${\tt Abs_Err2}$$ $$1.07\times 10^{-6}$$ $$5.35\times 10^{-7}$$ $$1.78\times 10^{-7}$$ $$4.45\times 10^{-8}$$ $$8.91\times 10^{-9}$$ $${\tt Rel_Err2}$$ $$3.93\times 10^{-7}$$ $$3.11\times 10^{-7}$$ $$2.48\times 10^{-7}$$ $$2.04\times 10^{-7}$$ $$1.73\times 10^{-7}$$ (2.13)–NewStr $$1.83\times 10^{-7}$$ $$5.24\times 10^{-8}$$ $$1.19\times 10^{-8}$$ $$2.22\times 10^{-9}$$ $$3.47\times 10^{-10}$$ (2.14)–NewStr $$1.06\times 10^{-7}$$ $$7.30\times 10^{-8}$$ $$5.43\times 10^{-8}$$ $$4.30\times 10^{-8}$$ $$3.49\times 10^{-8}$$ $${\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ 0.17 $$6.73\times 10^{-2}$$ $$5.01\times 10^{-2}$$ $$3.92\times 10^{-2}$$ $$3.23\times 10^{-2}$$ $${\rm cond}_{\rm rel}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ $$9.81\times 10^{-2}$$ $$4.83\times 10^{-2}$$ $$1.09\times 10^{-2}$$ $$2.00\times 10^{-3}$$ $$3.21\times 10^{-4}$$ It is seen from Table 12 that both (2.13) and (2.14) are very sharp, which justifies our new strategy for estimating $${\rm cond}_{\rm abs}(\varphi_{\ell},A)$$ and $${\rm cond}_{\rm rel}(\varphi_{\ell},A)$$. However, we find that the values of (2.13)-NewStr and (2.14)-NewStr may be a little smaller than those of $${\tt Abs\_Err2}$$ and $${\tt Rel\_Err2}$$ in some cases. In fact, the new strategy gives only some approximations to the absolute and the relative condition numbers, which are neither upper bounds nor lower bounds theoretically. 5. Concluding remarks In this article, we consider the computations, error analysis, implementations and applications of $$\varphi$$-functions for large sparse matrices with low rank or with fast decaying singular values. Given an SCR of $$A$$, we take into account how to compute the approximations to matrix functions $$\varphi_{\ell}(A)$$ efficiently for $$\ell=0,1,2,\ldots,p$$ and to estimate their 2-norm Fréchet relative and absolute condition numbers effectively. The key is to reduce the computation of $$\varphi_{\ell}$$-functions of a large matrix to that of $$\varphi_{\ell+1}$$-functions of some smaller matrices of size $$r$$-by-$$r$$, where $$r$$ is the numerical rank of the large matrix in question. Moreover, the new algorithms need to store only two $$n$$-by-$$r$$ sparse matrices and some $$r$$-by-$$r$$ matrices, rather than some $$n$$-by-$$n$$ possibly dense matrices. Our theoretical analysis and numerical experiments indicate that the error of the computed $$\varphi$$-function is of the order of the product of the error from SCR approximation with the condition number. Numerical experiments demonstrate the numerical performance of our new algorithms for $$\varphi$$-functions and for the condition numbers and show the effectiveness of our theoretical results. On the other hand, the efficiency of our new approach is closely related to that of reduced-rank approximation of large sparse matrices. Thus, a promising research area is to investigate new technologies to improve the performance of the SCR algorithm on very large matrices. Another interesting topic is to combine other advanced algorithms such as the randomized singular value decomposition algorithm by Halko et al. (2011) and Mahoney (2011) with our new strategies for the approximation of functions of large sparse matrices. Further, our new strategy also applies to other functions that can be expanded into power series. These are interesting topics and deserve further investigation. Acknowledgements We wish to thank the associate editor and the anonymous reviewers for their insightful comments and constructive suggestions that greatly improved the presentation of this article. Specifically, we thank a reviewer for the concise proof of Theorem 2.4. Meanwhile, we are grateful to Juan-juan Tian for helpful discussions. Funding National Science Foundation of China (11371176); Natural Science Foundation of Jiangsu Province; Talent Introduction Program of China University of Mining and Technology to G. Wu. Footnotes 1ftp://ftp.cs.umd.edu/pub/stewart/reports/Contents.html. 2http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html. 3http://vision.ucsd.edu/datasets/yale_face_dataset_original/yalefaces.zip. 4http://vision.ucsd.edu/~leekc/ExtYaleDatabase/Yale%20Face%20Database.html. References Ahmed N. ( 2013 ) Exponential discriminant regularization using nonnegative constraint and image descriptor. IEEE 9th International Conference on Emerging Technologies , pp. 1 – 6 . Al-Mohy A. & Higham N. J. ( 2011 ) Compution the action of the matrix exponential, with an application to exponential integrators. SIAM J. Sci. Comput ., 33 , 488 – 511 . Google Scholar CrossRef Search ADS Benzi M. , Estrada E. & Klymko C. ( 2013 ) Ranking hubs and authorities using matrix functions. Linear Algebra Appl. , 438 , 2447 – 2474 . Google Scholar CrossRef Search ADS Berland H. , Skaflestad B. & Wright W. ( 2007 ) EXPINT–A MATLAB package for exponential integrators. ACM Trans. Math. Software , 33 , Article 4. Berry M. & Browne M. ( 1999 ) Understanding Search Engines: Mathematical Modeling and Text Retrieval . Philadelphia : SIAM. Berry M. , Drmač, Z. & Jessup E. ( 1999 ) Matrices, vector spaces, and information retrieval. SIAM Rev. , 41 , 335 – 362 . Google Scholar CrossRef Search ADS Berry M. , Dumais S. & O’Brien G. ( 1995 ) Using linear algebra for intelligent information retrieval. SIAM Rev. , 37 , 573 – 595 . Google Scholar CrossRef Search ADS Berry M. , Pulatova S. & Stewart G. W. ( 2005 ) Computing sparse reduced-rank approximations to sparse matrices. ACM Trans. Math. Software , 31 , 252 – 269 . Google Scholar CrossRef Search ADS Beylkin G. , Keiser J. & Vozovoi L. ( 1998 ) A new class of time discretization schemes for the solution of nonlinear PDEs. J. Comput. Phys. , 147 , 362 – 387 . Google Scholar CrossRef Search ADS Celledon E. & Iserles A. ( 2000 ) Approximating the exponential from a Lie algebra to a Lie group. Math. Comp. , 69 , 1457 – 1480 . Google Scholar CrossRef Search ADS Davies P. , Higham N. J. & Tisseur F. ( 2001 ) Analysis of the Cholesky method with iterative refinement for solving the symmetric definite generalized eigenproblem. SIAM J. Matrix Anal. Appl. , 23 , 472 – 493 . Google Scholar CrossRef Search ADS Davis T. & Hu Y. ( 2011 ) The University of Florida sparse matrix collection. ACM Trans. Math. Software , 38 , Article 1 . Available at http://www.cise.ufl.edu/research/sparse/matrices/. Dornaika F. & Bosaghzadeh A. ( 2013 ) Exponential local discriminant embedding and its application to face recognition. IEEE Trans. Syst. Man Cybern. B Cybern. , 43 , 921 – 934 . Duda R. , Hart P. & Stork D. ( 2000 ) Pattern Classification , 2nd edn. New York : Wiley. Estrada E. & Hatano N. ( 2008 ) Communicability in complex networks Phys. Rev. E , 77 , 036111 . Google Scholar CrossRef Search ADS Estrada E. & Higham D. J. ( 2010 ) Network properties revealed through matrix functions. SIAM Rev. , 52 , 671 – 696 . Google Scholar CrossRef Search ADS Estrada E. & Rodríguez-Velázquez J. ( 2005 ) Subgraph centrality in complex networks. Phys. Rev. E , 71 , 056103. Google Scholar CrossRef Search ADS Golub G. H. & Van Loan C. F. ( 2013 ) Matrix Computations , 4th edn. Baltimore, MD : Johns Hopkins University Press. Halko N. , Martinsson P. & Tropp J. ( 2011 ) Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Rev. , 53 , 217 – 288 . Google Scholar CrossRef Search ADS Hansen P. ( 1998 ) Rank-Deficient and Discrete Ill-Posed Problems. Philadelphia : SIAM . Google Scholar CrossRef Search ADS Higham N. J. ( 2002 ) Accuracy and Stability of Numerical Algorithms , 2nd edn . Philadelphia : SIAM. Google Scholar CrossRef Search ADS Higham N. J. ( 2008a ) Functions of Matrices: Theory and Computation . Philadelphia : SIAM. Google Scholar CrossRef Search ADS Higham N. J. ( 2008b ) The Matrix Function Toolbox . Available at http://www.maths.manchester.ac.uk/~higham/ mftoolbox. Higham N. J. ( 2009 ) The scaling and squaring method for the matrix exponential revisited. SIAM Rev. , 51 , 747 – 764 . Google Scholar CrossRef Search ADS Higham N. J. , Mackey D. & Tisseur F. ( 2009 ) Definite matrix polynomials and their linearization by definite pencils. SIAM J. Matrix Anal. Appl. , 31 , 478 – 502 . Google Scholar CrossRef Search ADS Hochbruck M. , Lubich C. & Selhofer H. ( 1998 ) Exponential integrators for large systems of differential equations. SIAM J. Sci. Comput. , 19 , 1552 – 1574 . Google Scholar CrossRef Search ADS Hochbruck M. & Ostermann A. ( 2010 ) Exponential integrators. Acta Numer. , 209 – 286 . Hofmannn B. ( 1986 ) Regularization for Applied Inverse and Ill-posed Problems . Stuttgart, German : Teubner. Google Scholar CrossRef Search ADS Jiang P. & Berry M. ( 2000 ) Solving total least squares problems in information retrieval. Linear Algebra Appl. , 316 , 137 – 156 . Google Scholar CrossRef Search ADS Kassam A. & Trefethen L. N. ( 2005 ) Fourth-order time-stepping for stiff PDEs. SIAM J. Sci. Comput. , 26 , 1214 – 1233 . Google Scholar CrossRef Search ADS Lu Y. ( 2003 ) Computing a matrix function for exponential integrators. J. Comput. Appl. Math. , 161 , 203 – 216 . Google Scholar CrossRef Search ADS Mahoney M. ( 2011 ) Randomized algorithms for matrices and data . Foundations and Trends in Machine Learning, vol. 3, Boston : Now publishers inc., pp. 123 – 224 . Moler C. & Van Loan C. F. ( 2003 ) Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM Rev. , 45 , 3 – 49 . Google Scholar CrossRef Search ADS Niesen J. & Wright W. ( 2012 ) Algorithm 919: a Krylov subspace algorithm for evaluating the $$\varphi$$-functions appearing in exponential integrators. ACM Trans. Math. Software , 38 , Article 22. Park C. & Park H. ( 2008 ) A comparision of generalized linear discriminant analysis algorithms. Pattern Recogn. , 41 , 1083 – 1097 . Google Scholar CrossRef Search ADS Rice J. ( 1966 ) A theory of condition. SIAM J. Numer. Anal. , 3 , 287 – 310 . Google Scholar CrossRef Search ADS Schmelzer T. & Trefethen L. N. ( 2007 ) Evaluating matrix functions for exponential integrators via Carathéodory-Fejér approximation and contour integrals. Electron. Trans. Numer. Anal. , 29 , 1 – 18 . Sidje R. ( 1998 ) EXPOKIT: Software package for computing matrix exponentials. ACM Trans. Math. Software , 24 , 130 – 156 . Google Scholar CrossRef Search ADS Skaflestad B. & Wright W. ( 2009 ) The scaling and squaring method for matrix functions related to the exponential. Appl. Numer. Math. , 59 , 783 – 799 . Google Scholar CrossRef Search ADS Stewart G. W. ( 1999 ) Four algorithms for the efficient computation of truncated pivoted QR approximations to a sparse matrix. Numer. Math. , 83 , 313 – 323 . Google Scholar CrossRef Search ADS Stewart G. W. ( 2005 ) Error analysis of the quasi-Gram–Schmidt algorithm. SIAM J. Matrix Anal. Appl ., 27 , 493 – 506 . Google Scholar CrossRef Search ADS Stuart G. & Berry M. ( 2003 ) A comprehensive whole genome bacterial phylogeny using correlated peptide motifs defined in a high dimensional vector space. J. Bioinform. Comput. Bio. , 1 , 475 – 493 . Google Scholar CrossRef Search ADS Tspras P. Datasets for Experiments on Link Analysis Ranking Algorithms, Available at http://www.cs.toronto.edu/~tsap/experiments/datasets/index.html. Wang S. , Chen H. , Peng X. & Zhou C. ( 2011 ) Exponential locality preserving projections for small sample size problem. Neurocomputing , 74 , 36 – 54 . Wang S. , Yan S. , Yang J. , Zhou C. & Fu X. ( 2014 ) A general exponential framework for dimensionality reduction. IEEE Trans. Image Process. , 23 , 920 – 930 . Google Scholar CrossRef Search ADS PubMed Wu G. & Feng T. ( 2015 ) A theoretical contribution to the fast implementation of null linear discriminant analysis with random matrix multiplication. Numer. Linear Algebra Appl. , 22 , 1180 – 1188 . Google Scholar CrossRef Search ADS Wu G. , Feng T. , Zhang L. & Yang M. ( 2017 ) Inexact implementation using Krylov subspace methods for large scale exponential discriminant analysis with applications to high dimensionality reduction problems. Pattern Recogn. , 66 , 328 – 341 . Google Scholar CrossRef Search ADS Wu G. , Zhang L. & Xu T. ( 2016 ) A framework of the harmonic Arnoldi method for evaluating $$\varphi$$-functions with applications to exponential integrators. Adv. Comput. Math. , 42 , 505 – 541 . Google Scholar CrossRef Search ADS Yan L. & Pan J. ( 2010 ) Two-dimensional exponential discriminant analysis and its application to face recognition. International Conference on Computational Aspects of Social Networks (CASoN) , pp. 528 – 531 . Zhang T. , Fang B. , Tang Y. , Shang Z. & Xu B. ( 2010 ) Generalized discriminant analysis: a matrix exponential approach. IEEE Trans. Syst. Man Cybern. B Cybern. , 40 , 186 – 197 . Google Scholar CrossRef Search ADS PubMed Zhang Z. , Zha H. & Simon H. ( 2002 ) Low-rank approximations with sparse factors I: basic algorithms and error analysis. SIAM J. Matrix Anal. Appl. , 23 , 706 – 727 . Google Scholar CrossRef Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

New algorithms for approximating $$\varphi$$ -functions and their condition numbers for large sparse matrices

Loading next page...
 
/lp/ou_press/new-algorithms-for-approximating-varphi-functions-and-their-condition-5UiSRjwZGR
Publisher
Oxford University Press
Copyright
© The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved.
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx035
Publisher site
See Article on Publisher Site

Abstract

Abstract Given an $$n$$-by-$$n$$ large sparse matrix $$A$$ with low rank or with fast decaying singular values, we propose new algorithms for consecutively approximating the functions $$\varphi_{\ell}(A)$$ as well as their 2-norm condition numbers for $$\ell=0,1,\ldots,p$$. The key is to reduce the computation of $$\varphi_{\ell}$$-functions of the large matrix to that of $$\varphi_{\ell+1}$$-functions of some smaller matrices of size $$r$$-by-$$r$$, where $$r$$ is the numerical rank of $$A$$. For storage reasons, the new algorithms need to store only two $$n$$-by-$$r$$ sparse matrices and some $$r$$-by-$$r$$ matrices, rather than some $$n$$-by-$$n$$ possibly dense matrices. Some error analysis is given. Numerical experiments illustrate the numerical behavior of our new algorithms and show the practical utility of new strategies. 1. Introduction In recent years, a great deal of attention has been focused on efficient and accurate evaluation of matrix functions closely related to $$\varphi$$-functions (Hochbruck et al., 1998; Lu, 2003; Moler & Van Loan, 2003; Berland et al., 2007; Schmelzer & Trefethen, 2007; Higham, 2008a; Skaflestad & Wright, 2009; Hochbruck & Ostermann, 2010; Niesen & Wright, 2012; Wu et al., 2016). For instance, exponential integrators make use of the matrix exponential and related matrix functions within the formulations of the numerical methods, and the evaluation of matrix functions is crucial for accuracy, stability and efficiency of exponential integrators (Hochbruck & Ostermann, 2010). The $$\varphi$$-functions are defined for scalar arguments by the integral representation as \begin{equation} \varphi_{0}(z)={\rm exp}(z)\quad{\rm and}\quad\varphi_{\ell}(z)=\frac{1}{(\ell-1)!}{\int_0^1 {\rm exp}\big((1-\theta)z\big)\theta^{\ell-1}\,{\rm d}\theta},\quad \ell=1,2,\ldots, \end{equation} (1.1) and they satisfy the recurrence relations \begin{equation} \varphi_{\ell}(z)=z\varphi_{\ell+1}(z)+\frac{1}{\ell!},\quad \ell=0,1,2,\ldots. \end{equation} (1.2) Moreover, $$\varphi_{\ell}(z)$$ can be expressed as the following power series whose radius of convergence is $$+\infty$$: \begin{equation} \varphi_{\ell}(z)=\sum_{i =0}^{\infty}\frac{z^i}{(i+\ell)!}. \end{equation} (1.3) The above definitions can be extended to matrices instead of scalars by using any of the available definitions of matrix functions (Higham, 2008a; Hochbruck & Ostermann, 2010). In a wide range of applications, such as the matrix exponential discriminant analysis (EDA) method for data dimensionality reduction (Yan & Pan, 2010; Zhang et al., 2010; Wang et al., 2011; Ahmed, 2013; Dornaika & Bosaghzadeh, 2013; Wang et al., 2014; Wu et al., 2017) and the complex network analysis method based on matrix functions (Estrada & Rodríguez-Velázquez, 2005; Estrada & Hatano, 2008; Estrada & Higham, 2010; Benzi et al., 2013), it is required to compute the matrix exponential of large-scale and low-rank matrices. Given a large-scale matrix $$A$$ with low rank or with fast decaying singular values, in this article we are interested in approximating several $$\varphi$$-functions consecutively. Let $$\sigma_j$$ be the $$j$$th largest singular value of $$A$$. By ‘fast decaying singular values’, we mean $$\sigma_j=\mathcal{O}(\rho^{-j}),\rho>1$$ or $$\sigma_j=\mathcal{O}(j^{-\alpha}),\ \alpha>1$$ (Hofmannn, 1986). Such matrices appear frequently in diverse application areas such as data dimensionality reduction (Duda et al., 2000), complex network analysis (Estrada & Higham, 2010), discretization of ill-posed operator equations that model many inverse problems (Hansen, 1998), randomized algorithms for matrix approximations (Halko et al., 2011; Mahoney, 2011), finite element discretization (Davis & Hu, 2011) and various other applications. In spite of the high demand for efficient methods to solve the matrix $$\varphi$$-functions in various fields of computational sciences there is no easy way to solve this type of problem. Indeed, when $$A$$ is large, both the computational cost and the storage requirements are demanding; moreover, $$\varphi_{\ell}(A)$$ is generally dense even if $$A$$ is sparse (Higham, 2008a). Some available methods are suitable only for medium-sized matrices. For instance, a MATLAB toolbox called EXPINT was provided by Berland et al. (2007). Kassam and Trefethen proposed to approximate $$\varphi$$-functions with a contour integral, which works well as long as the contour of integration is suitably chosen (Kassam & Trefethen, 2005). However, the contour is generally problem dependent and difficult to determine in advance. Another way is to reduce the computation of $$\varphi_{\ell}(A)$$ to that of a matrix exponential with a larger size (Sidje, 1998; Lu, 2003; Al-Mohy & Higham, 2011), which is unfeasible as $$A$$ is large. A third way is based on a modification of the scaling and squaring technique (Skaflestad & Wright, 2009), the most commonly used approach for computing the matrix exponential (Higham, 2009). The most well-known methods for approximating $$\varphi$$-functions for large sparse matrices are the Krylov subspace methods (Niesen & Wright, 2012; Wu et al., 2016). However, Krylov subspace methods are applicable to the computation of $$\varphi$$-functions on given (block) vectors, while the main aim of this article is to compute approximations of $$\varphi$$-functions of large sparse matrices. In practical calculations, it is important to know how accurate the computed solution is and how small perturbations in input data can affect outputs Higham (2002). Therefore, it is crucial to give some error analysis and to understand the sensitivity of matrix functions to perturbations in the data. Sensitivity is measured by condition numbers. For matrix functions, condition number can be expressed in terms of the norm of the Fréchet derivative, and it is often measured by using the 1-norm Higham (2008a,b). In this work, we focus on the 2-norm, for which there is a tighter relation between the norm of the Fréchet derivative and the norm of the Kronecker matrix, and we consider how to evaluate the Fréchet 2-norm condition numbers of $$\varphi$$-functions efficiently. Given a large scale matrix $$A$$ with low rank or with fast decaying singular values, we propose new algorithms for evaluating several $$\varphi$$-functions and their absolute and relative condition numbers consecutively. Our new algorithms are based on the sparse column–row approximation (SCR) of large sparse matrices (Stewart, 1999, 2005; Berry et al., 2005). An advantage of our algorithms is that there is no need to form and store the $$\varphi$$-functions of large matrices. The overhead is only to compute $$\varphi$$-functions of some $$r$$-by-$$r$$ matrices and to store two $$n$$-by-$$r$$ sparse matrices, where $$r$$ is the (numerical) rank of $$A$$. This article is organized as follows. In Section 2, we present the main algorithm and give some error analysis on this method. In Section 3, we propose a novel strategy for estimating the absolute and the relative 2-norm condition numbers of $$\varphi$$-functions. In Section 4, numerical experiments are given to show the effectiveness of our new strategies. Some concluding remarks are given in Section 5. Throughout this article, we denote by $$\widetilde{A}=XTY^{\rm T}$$ an SCR to $$A$$, and by $$r$$ the (numerical) rank of $$A$$. Let $$\|\cdot\|_2,\|\cdot\|_1$$ be the 2-norm and the 1-norm of a vector or a matrix, and let $$\|\cdot\|_{\rm F}$$ be the Frobenius norm of a matrix. We denote by $$\otimes$$ the Kronecker product, and by $${\rm vec}(\cdot)$$ the ‘vec operator’ that stacks the columns of a matrix into a long vector. Let $$I$$ and $$\mathbf{0}$$ be the identity matrix and the zero matrix, respectively, whose order is clear from context. We focus only on real matrices in this article. However, all the results can be extended to complex matrices in a similar way. 2. A new algorithm for approximating $$\varphi$$-functions consecutively In this section, we will present a new method for approximating $$\varphi$$-functions of large sparse matrices with low rank or with fast decaying singular values and give some error analysis on this method. Given an $$n\times n$$ large sparse matrix $$A$$, we first find a reduced-rank approximation $$XTY^{\rm T}$$ to it, where both $$X$$ and $$Y$$ are full column rank and $$T$$ is nonsingular. This type of problem arises in a number of applications such as information retrieval, computational biology and complex network analysis (Berry et al., 1995, 1999, 2005; Berry & Browne, 1999; Jiang & Berry, 2000; Zhang et al., 2002; Stuart & Berry, 2003). A widely used reduced-rank approximation is the truncated singular value decomposition (TSVD) (Golub & Van Loan, 2013), which is known to be optimal in the sense that the Frobenius norm $$\|A-XTY^{\rm T}\|_F$$ is minimized. Unfortunately, this method computes the full decomposition and is not suitable for very large matrices. An alternative is the randomized singular value decomposition (Halko et al., 2011; Mahoney, 2011), which generally gives results comparable to TSVD. However, for a large and sparse matrix, the situation is not so simple: since the resulting factors $$X, T$$ and $$Y$$ are generally not sparse, the storage requirements and operation counts will become proportional to the number of nonzero elements. In Stewart (1999), Stewart introduced a quasi-Gram–Schmidt algorithm that produces a sparse QR factorization to $$A$$. Based on the quasi-Gram-Schmidt algorithm, an SCR algorithm was proposed (Stewart, 1999). This algorithm first applies the quasi-Gram–Schmidt algorithm to the columns of $$A$$ to get a representative set of columns $$X$$ of $$A$$, and an upper triangular matrix $$R$$. Let the error in the corresponding reduced-rank decomposition be $$\epsilon_{\rm col}$$. It then applies the same algorithm to $$A^{\rm T}$$ to get a representative set $$Y^{\rm T}$$ of rows and another upper triangular matrix $$S$$. Let the error be $$\epsilon_{\rm row}$$. Then the SCR method seeks a matrix $$T$$ such that $$\|A-XTY^{\rm T}\|_F^2=\min$$, and the minimizer turns out to be $$ T=R^{-1}R^{-{\rm T}}(X^{\rm T}AY)S^{-1}S^{-{\rm T}}$$ (see Stewart, 1999; Berry et al., 2005). Moreover, we have \begin{equation} \|A-XTY^{\rm T}\|_F^2\leq \epsilon_{\rm col}^2+\epsilon_{\rm row}^2. \end{equation} (2.1) The matrix $$XTY^{\rm T}$$ is called an SCR to $$A$$, where $$X,Y\in\mathbb{R}^{n\times r}$$ are sparse and of full column rank, $$T\in\mathbb{R}^{r\times r}$$ is nonsingular and $$r$$ is the numerical rank of $$A$$. In this approximation, $$X$$ consists of a selection of the columns of $$A$$, and $$Y$$ consists of a selection of the rows of $$A$$, so that when $$A$$ is sparse then so are both $$X$$ and $$Y$$. Error analysis of the quasi-Gram–Schmidt algorithm was given in Stewart (2005). One is recommended to see Berry et al. (2005) and Stewart (1999) for more details on this algorithm. Given any rank-revealing decomposition of $$A$$, the following theorem indicates that the computation of a matrix function with respect to $$A$$ can be reduced to that of another matrix function with respect to an $$r\times r$$ matrix, where $$r$$ is the (numerical) rank of $$A$$. Theorem 2.1 Suppose that the power series $$f(z)=\sum_{i=0}^{\infty}\alpha_i z^i$$ has radius of convergence $$\rho$$. Let $$\widetilde{A}=XTY^{\rm T}\in\mathbb{R}^{n\times n}$$, where $$X,Y\in\mathbb{R}^{n\times r}$$ and $$T\in\mathbb{R}^{r\times r}$$. Let $$g(z)=\sum_{i=1}^{\infty}\alpha_i z^{i-1}$$ and suppose that $$\|\widetilde{A}\|<\rho$$; then \begin{align} f(\widetilde{A})=f(\mathbf{0})+Xg(Z)TY^{\rm T}, \end{align} (2.2) where $$Z=T(Y^{\rm T}X)\in\mathbb{R}^{r\times r}$$, and $$\mathbf{0}$$ stands for the zero matrix. Proof. As $$k\geq1$$, we see that $$\widetilde{A}^k=XZ^{k-1}TY^{\rm T}$$. Thus, \begin{align*} f(\widetilde{A})&=\alpha_0I+\alpha_1\widetilde{A}+\cdots+\alpha_k\widetilde{A}^k+\cdots\\ &=\alpha_0I+X\big(\alpha_1I+\alpha_2Z+\cdots+\alpha_kZ^{k-1}+\cdots\big)TY^{\rm T}\\ &=f(\mathbf{0})+Xg(Z)TY^{\rm T}. \end{align*} □ Remark 2.2 By (1.3) and Theorem 2.1, the computation of $$\varphi_{\ell}(A)$$ can be reduced to that of $$\varphi_{\ell+1}(Z)$$. More precisely, \begin{equation} \varphi_{\ell}(XTY^{\rm T})=\frac{1}{\ell!}I+X\big(\varphi_{\ell+1}(Z)T\big)Y^{\rm T},\quad \ell=0,1,2,\ldots. \end{equation} (2.3) We mention that for the matrix exponential, the main idea of Theorem 2.1 has already been given in (Celledon & Iserles, 2000, Proposition 3) and both (2.2) and (2.3) can be viewed as generalizations of this proposition. Let $$\widetilde{A}=XTY^{\rm T}$$ be an SCR to $$A$$. We make use of $$\varphi_{\ell}(\widetilde{A})$$ as an approximation to $$\varphi_{\ell}(A)$$. The main reason we choose the SCR is that both $$X$$ and $$Y$$ are sparse if $$A$$ is sparse. Given a large matrix $$A$$ with low rank or with fast decaying singular values, we have the following algorithm for consecutively approximating $$\varphi_{\ell}(\widetilde{A}),\,\ell=0,1,\ldots,p$$. Algorithm 1 An algorithm for approximating $$\varphi$$-functions of large sparse matrices consecutively. 1. Compute a reduced-rank approximation $$XTY^{\rm T}$$ to $$A$$ by using the sparse column–row approximation (SCR) algorithm. 2. Compute $$\varphi$$-functions of small-sized matrices: $$\varphi_{\ell+1}(TY^{\rm T}X),~\ell=0,1,\ldots,p$$. 3. Store $$X,Y,T$$ and $$\varphi_{\ell+1}(TY^{\rm T}X)$$ for $$\varphi_{\ell}(\widetilde{A}),~\ell=0,1,\ldots,p$$. If desired, form $$\varphi_{\ell}(\widetilde{A})$$ in terms of (2.3) and use them as approximations to $$\varphi_{\ell}(A),~\ell=0,1,\ldots,p$$. We stress that the main overhead of Algorithm 1 consists of computing the SCR approximation of $$A$$ provided that $$r$$ is small. This involves performing two quasi-Gram–Schmidt algorithms with respect to $$A$$ and $$A^{\rm T}$$. In the quasi-Gram–Schmidt algorithm, it is required to extract the pivot column from $$A$$ or $$A^{\rm T}$$ and to calculate several (sparse) matrix–vector operations. For a discussion of the computational complexities of the SCR approximation, refer to Berry et al. (2005, pp. 265, 268). In order to allow a fair comparison, the computation of the SCR approximation is also included in the computational comparisons and CPU timings of our new algorithm. Remark 2.3 Obviously, an advantage of the proposed method is its simplicity. In conventional methods, the cost amounts to $$\mathcal{O}\big((p+1) n^3\big)$$ flops for the computation of $$\varphi_{\ell}(A)$$ (Beylkin et al., 1998; Berland et al., 2007; Higham, 2008a; Skaflestad & Wright, 2009). Given a sparse reduced-rank approximation to $$A$$, Theorem 2.1 reduces the computation of $$\varphi_{\ell}(A)$$ to that of $$\varphi_{\ell+1}$$ functions with respect to the $$r$$-by-$$r$$ matrix $$Z=T(Y^{\rm T}X)$$, in $$\mathcal{O}\big((p+1) r^3\big)$$ flops, $$\ell=0,1,\ldots,p$$. For storage reasons, it needs to store only two $$n$$-by-$$r$$ sparse matrices $$X,Y$$, and some smaller matrices of size $$r$$-by-$$r$$, rather than the $$n$$-by-$$n$$ possibly dense matrices $$\varphi_{\ell}(A),~\ell=0,1,\ldots,p$$. Thus, the new method can consecutively approximate $$\varphi_{\ell}(A),\ell=0,1,\ldots,p$$ and reduce the computational complexities significantly as $$r\ll n$$. In practice, however, most data are inexact or uncertain. Indeed, even if the data were exact, the computations would be subject to rounding errors. So it is important to give error analysis and understand the sensitivity of matrix functions to perturbations in the input data. Sensitivity is measured by condition numbers. For a matrix function, the condition number can be expressed in terms of the norm of the Fréchet derivative. The Fréchet derivative $$L_f(A,E)$$ of a matrix function $$f: \mathbb{R}^{n\times n}\rightarrow\mathbb{R}^{n\times n}$$ at a point $$A\in\mathbb{R}^{n\times n}$$ is a linear mapping such that for all $$E\in\mathbb{R}^{n\times n}$$, \begin{equation} f(A+E)-f(A)-L_f(A,E)=o(\|E\|) \end{equation} (2.4) (see Higham, 2008a, p. 56). The absolute and the relative condition numbers of a matrix function $$f(A)$$ are defined as \begin{equation} {\rm cond}_{\rm abs}(f,A)=\|L_f(A)\|=\max_{E\neq \mathbf{0}}\frac{\|L_f(A,E)\|}{\|E\|} \end{equation} (2.5) and \begin{equation} {\rm cond}_{\rm rel}(f,A)=\frac{\|L_f(A)\|\|A\|}{\|f(A)\|}, \end{equation} (2.6) respectively (see Rice, 1966; Higham, 2008a). Theoretically, the Fréchet derivative can be obtained from applying any existing method for computing the matrix function of a $$2n\times 2n$$ matrix and using the formula \begin{equation} f\left(\left[\begin{array}{@{}cc@{}} A & E \\ \mathbf{0} & A \\\end{array} \right]\right) =\left[\begin{array}{@{}cc@{}} f(A) & L_f(A,E) \\ \mathbf{0} & f(A)\\ \end{array} \right] \end{equation} (2.7) (see Higham, 2008a, Algorithm 3.17). However, it requires $$\mathcal{O}(n^5)$$ flops, assuming that the computation of $$L_f(A,E)$$ takes $$\mathcal{O}(n^3)$$ flops (Higham, 2008a). So it is computationally demanding for large matrices. Let $$\widetilde{A}=XTY^{\rm T}$$ be an SCR to $$A$$, and let $$E=A-\widetilde{A}$$; then we see from (2.1) that $$\|E\|_F\leq\sqrt{\epsilon_{\rm col}^2+\epsilon_{\rm row}^2}$$. Thus, it is interesting to combine the existing error analysis for SCR with the theory of matrix functions to obtain error bounds for the new method. We first present the following theorem for Fréchet derivatives of $$\varphi$$-functions. Theorem 2.4 For any matrix $$E\in\mathbb{R}^{n\times n}$$, we have \begin{equation} \varphi_{\ell}(A+E)-\varphi_{\ell}(A)=\int_{0}^{1}\exp\big((1-s)A\big)s^{\ell}E\varphi_{\ell}\big(s(A+E)\big)\,{\rm d}s,\quad \ell=0,1,2,\ldots. \end{equation} (2.8) Proof. For the matrix exponential (Higham, 2008a, p. 238), we have \begin{equation} \exp(A+E)=\exp(A)+\int_0^1\exp\big((1-s)A\big)E\exp\big(s(A+E)\big)\,{\rm d}s, \end{equation} (2.9) and thus (2.8) holds for $$\ell=0$$. When $$\ell>0$$, we know that the function $$X(t)=t^{\ell}\varphi_{\ell}\big(t(A+E)\big)$$ is the solution of the initial value problem (Niesen & Wright, 2012, Lemma 2.1) $$ X'(t)=(A+E)X(t)+\frac{t^{\ell-1}}{(\ell-1)!}I=AX(t)+G(t,X),\quad X(0)={\bf 0}, $$ where $$ G(t,X)=EX(t)+\frac{t^{\ell-1}}{(\ell-1)!}I. $$ By the variation of constants formula, we get \begin{eqnarray*} X(t)&=&\exp(tA)X(0)+\int_0^t \exp\big((t-s)A\big)G\big(s,X(s)\big)\,{\rm d}s\\ &=&\int_0^t \exp\big((t-s)A\big)\left(EX(s)+\frac{s^{\ell-1}}{(\ell-1)!}I\right)\,{\rm d}s\\ &=&\int_0^t \exp\big((t-s)A\big)EX(s)\,{\rm d}s+\int_0^t \exp\big((t-s)A\big)\frac{s^{\ell-1}}{(\ell-1)!}\,{\rm d}s. \end{eqnarray*} Notice that $$X(s)=s^{\ell}\varphi_{\ell}\big(s(A+E)\big)$$ and $$t^{\ell}\varphi_{\ell}(tA)$$ is the solution of the initial value problem $$ U'(t)=AU(t)+\frac{t^{\ell-1}}{(\ell-1)!}I,\quad U(0)={\bf 0}. $$ Similarly, from the variation of constants formula, we have $$ t^{\ell}\varphi_{\ell}(tA)=\int_0^t \exp\big((t-s)A\big)\frac{s^{\ell-1}}{(\ell-1)!}\,{\rm d}s. $$ Therefore, $$X(t)=\int_0^t \exp\big((t-s)A\big)Es^\ell\varphi_{\ell}\big(s(A+E)\big)\,{\rm d}s+t^{\ell}\varphi_{\ell}(tA),$$ from which we get our result (2.8) for $$t=1$$. □ Using (2.8) to substitute for $$\varphi_{\ell}\big(s(A+E)\big)$$ inside the integral, we get \begin{equation} \varphi_{\ell}(A+E)=\varphi_{\ell}(A)+\int_{0}^{1}\exp\big((1-s)A\big)s^{\ell}E\varphi_{\ell}(sA)\,{\rm d}s+\mathcal{O}(\|E\|^2),\quad \ell=0,1,2,\ldots. \end{equation} (2.10) Combining with (2.4), we can present the following result for the Fréchet derivatives of $$\varphi$$-functions. Corollary 2.5 The Fréchet derivatives of $$\varphi$$-functions at $$A$$ in the direction $$E$$ are given by \begin{equation} L_{\varphi_{\ell}}(A,E)=\int_{0}^{1}\exp\big((1-s)A\big)s^{\ell}E\varphi_{\ell}(sA)\,{\rm d}s,\quad \ell=0,1,2,\ldots \end{equation} (2.11) As a result, \begin{equation} \varphi_{\ell}(A+E)=\varphi_{\ell}(A)+L_{\varphi_{\ell}}(A,E)+o(\|E\|). \end{equation} (2.12) The following theorem shows that the values of $$\epsilon_{\rm col}$$ and $$\epsilon_{\rm row}$$ used during the SCR will have a direct impact on the final accuracy of approximating $$\varphi_{\ell}(A)$$. Note that $$XTY^{\rm T}$$ can be any low-rank approximation to $$A$$ in this theorem. Theorem 2.6 Let $$XTY^{\rm T}$$ be an SCR to $$A$$, and $$\|A-XTY^{\rm T}\|=\varepsilon$$. Then we have \begin{equation} \|\varphi_{\ell}(A)-\varphi_{\ell}(XTY^{\rm T})\|\lesssim{\rm cond}_{\rm abs}(\varphi_{\ell},A)~\varepsilon \end{equation} (2.13) and \begin{equation} \frac{\|\varphi_{\ell}(A)-\varphi_{\ell}(XTY^{\rm T})\|}{\|\varphi_{\ell}(A)\|}\lesssim{\rm cond}_{\rm rel}(\varphi_{\ell},A)\frac{\varepsilon}{\|A\|}, \end{equation} (2.14) where $$\lesssim$$ represents omitting the high-order term $$o(\|E\|)$$. Proof. Let $$E=XTY^{\rm T}-A$$; then we get from (2.12) that \begin{eqnarray*} \|\varphi_{\ell}(A)-\varphi_{\ell}(XTY^{\rm T})\|&\lesssim&\|L_{\varphi_{\ell}}(A)\|\|E\| =\|L_{\varphi_{\ell}}(A)\|\varepsilon\nonumber\\ &=&{\rm cond}_{\rm abs}(\varphi_{\ell},A)\varepsilon, \end{eqnarray*} where we omit the high-order term $$o(\|E\|)$$. The upper bound (2.14) for the relative error is derived from (2.6) and (2.13). □ 3. A new strategy for estimating the absolute and the relative 2-norm condition numbers consecutively In terms of Theorem 2.6, it is crucial to consider how to estimate the absolute and the relative condition numbers $${\rm cond}_{\rm abs}(\varphi_{\ell},A)$$ and $${\rm cond}_{\rm rel}(\varphi_{\ell},A)$$ efficiently. Notice that \begin{equation} {\rm vec}(ACB)=(B^{\rm T}\otimes A){\rm vec}(C). \end{equation} (3.1) Since $$L_f$$ is a linear operator, we have \begin{equation} {\rm vec}(L_f(A,E))=K_f(A){\rm vec}(E) \end{equation} (3.2) for some $$K_f(A)\in\mathbb{R}^{n^2\times n^2}$$ that is independent of $$E$$. The matrix $$K_f(A)$$ is referred to as the Kronecker form of the Fréchet derivative (Higham, 2008a, p. 60). Specifically, we have \begin{equation} \|L_f(A)\|_F=\lambda_{\max}\big(K_f(A)^{\rm T} K_f(A)\big)^{1/2}=\|K_f(A)\|_2. \end{equation} (3.3) To estimate $$\|K_f(A)\|_2$$, the power method can be applied (Higham, 2008a, Algorithm 3.20). However, the power method lacks convergence tests, and because of its linear convergence rate the number of iterations required is unpredictable. Alternatively, the condition number can be based on the 1-norm (Higham, 2008a, Algorithm 3.22). One reason for using the 1-norm is its cheapness. Although there is no analogue to the relation (3.3) for the 1-norm, the next result gives a relation between $$\|K_f (A)\|_1$$ and $$\|L_f (A)\|_1$$. Theorem 3.1 (Higham, 2008a, Lemma 3.18). For $$A\in\mathbb{R}^{n\times n}$$ and any function $$f$$, \begin{equation} \frac{\|L_f(A)\|_1}{n}\leq \|K_f(A)\|_1\leq n\|L_f(A)\|_1. \end{equation} (3.4) In this article, we choose to work with the 2-norm, for which there is a tighter relation between the norm of the Fréchet derivative and the norm of the Kronecker matrix. The following theorem establishes a relation between $$\|K_f (A)\|_2$$ and $$\|L_f (A)\|_2$$. Theorem 3.2 For $$A\in\mathbb{R}^{n\times n}$$ and any function $$f$$, \begin{equation} \frac{\|L_f(A)\|_2}{\sqrt{n}}\leq \|K_f(A)\|_2\leq \sqrt{n}\|L_f(A)\|_2. \end{equation} (3.5) Proof. The proof is analogous to that by Higham (2008a, Lemma 3.18). □ Remark 3.3 Notice that $$\sqrt{n}$$ is used in the lower and upper bounds of (3.5), instead of $$n$$ in (3.4). Hence, there is a tighter relation between the 2-norm of the Fréchet derivative and the 2-norm of the Kronecker matrix. We are ready to propose a novel strategy to evaluate the absolute and the relative 2-norm condition numbers. Recall that the key idea of our new strategy in Section 2 is to reduce the computation of $$\varphi_{\ell}(A)$$ to that of $$\varphi_{\ell+1}(Z)$$. As $$ \varphi_{\ell}(\widetilde{A})=\frac{1}{\ell!}I+X\big(\varphi_{\ell+1}(Z)T\big)Y^{\rm T}, $$ if we denote $$ \widetilde{\varphi_{\ell}(\widetilde{A})}=\frac{1}{\ell!}I+X\big(\varphi_{\ell+1}(Z+F)T\big)Y^{\rm T} $$ then it follows from (2.4) that $$ L_{\varphi_{\ell+1}}(Z,F)=\varphi_{\ell+1}(Z+F)-\varphi_{\ell+1}(Z)+o(\|F\|)\quad \forall\,F\in\mathbb{R}^{r\times r}. $$ Hence, \begin{eqnarray} \widetilde{\varphi_{\ell}(\widetilde{A})}-\varphi_{\ell}(A)&=&X\big(\varphi_{\ell+1}(Z+F)-\varphi_{\ell+1}(Z)\big)TY^{\rm T}\nonumber\\ &=&XL_{\varphi_{\ell+1}}(Z,F)TY^{\rm T}+o(\|F\|). \end{eqnarray} (3.6) Let $$X=Q_1R_1,Y=Q_2R_2$$ be the (sparse) QR decompositions of $$X$$ and $$Y$$, respectively, where $$R_1,R_2\in\mathbb{R}^{r\times r}$$. We propose to use \begin{equation} {\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})=\|XL_{\varphi_{\ell+1}}(Z)TY^{\rm T}\|_2=\|R_1L_{\varphi_{\ell+1}}(Z)TR_2^{\rm T}\|_2 \end{equation} (3.7) as an approximation to the absolute 2-norm condition number $${\rm cond}_{\rm abs}(\varphi_{\ell},A)$$. Notice that the key of our strategy is (3.6) and the fact that both $$Q_1$$ and $$Q_2$$ are orthonormal. On the other hand, Theorem 2.1 also provides a cheap way to estimate 2-norms of $$\varphi_{\ell}(A),~\ell=0,1,\ldots,p$$. Indeed, from (2.3) we get \begin{equation} \Big|\|R_1\big(\varphi_{\ell+1}(Z)T\big)R_2^{\rm T}\|_2-1/\ell!\Big|\leq\|\varphi_{\ell}(\widetilde{A})\|_2\leq 1/\ell!+\|R_1\varphi_{\ell+1}(Z)TR_2^{\rm T}\|_2. \end{equation} (3.8) Thus, we exploit \begin{equation} \eta_{\ell}=\big\|R_1\varphi_{\ell+1}(Z)TR_2^{\rm T}\big\|_2,\quad \ell=0,1,\ldots,p \end{equation} (3.9) as approximations to $$\|\varphi_{\ell}({A})\|_2$$. Then the relative 2-norm condition number $${\rm cond}_{\rm rel}(\varphi_{\ell},A)$$ can be approximated by \begin{equation} {\rm cond}_{\rm rel}^{\rm est}(\varphi_{\ell},\widetilde{A})=\frac{\|A\|_2\|R_1L_{\varphi_{\ell+1}}(Z)TR_2^{\rm T}\|_2}{\|R_1\varphi_{\ell+1}(Z)TR_2^{\rm T}\|_2}; \end{equation} (3.10) refer to (2.14). In summary, we have the following algorithm for estimating 2-norm condition numbers of $$\varphi$$-functions. Algorithm 2 An algorithm for estimating the absolute and the relative 2-norm condition numbers of $$\varphi$$-functions consecutively. 1. Compute a reduced-rank approximation $$XTY^{\rm T}$$ to $$A$$ by using the sparse column–row approximation (SCR) algorithm. 2. Compute the R-factors $$R_1$$ and $$R_2$$ by using the (sparse) QR decompositions of $$X$$ and $$Y$$, respectively. 3. Solve $$\varphi_{\ell+1}(Z)$$ and $$L_{\varphi_{\ell+1}}(Z)$$ for the $$r$$-by-$$r$$ matrix $$Z=T(Y^{\rm T}X),~\ell=0,1,\ldots,p$$. 4. Use (3.7) and (3.10) as approximations to the absolute and the relative 2-norm condition numbers of $$\varphi_{\ell}(A),~\ell=0,1,\ldots,p$$. Remark 3.4 Two remarks are given. First, there is no need to form and store the Q-factors $$Q_1$$ and $$Q_2$$ in the QR decompositions of $$X$$ and $$Y$$, which may be dense in practice. Indeed, it is only necessary to evaluate 2-norms of some $$r$$-by-$$r$$ matrices in Algorithm 2. Second, similarly to Algorithm 1, one can evaluate the relative and absolute 2-norm condition numbers consecutively for $$\ell=0,1,\ldots,p$$. Next, we will give another motivation for using $$L_{\varphi_{\ell+1}}(Z)$$ for approximating the norm of $$L_{\varphi_{\ell}}(\widetilde{A})$$, where $$\widetilde{A}=XTY^{\rm T}$$. Indeed, the theorem below motivates computing the norm of the approximating formula (3.7) and (3.10). Theorem 3.5 With the notation of Theorem 2.1, if we define $$\widetilde{A}=XTY^{\rm T}$$ and $$W=YT^{\rm T}$$ then \begin{equation} K_f(\widetilde{A})={\it{\Psi}}_1+{\it{\Psi}}_2+{\it{\Psi}}_3, \end{equation} (3.11) where \begin{align*} {\it{\Psi}}_1 & =\alpha_1 I\otimes I,\\ {\it{\Psi}}_2 & =\left(W\otimes I\right)\sum_{i=2}^{\infty}\alpha_i\left((Z^{\rm T})^{i-2}\otimes I\right)\left(X\otimes I\right)^{\rm T}+ \left(I\otimes X\right)\sum_{i=2}^{\infty}\alpha_i\left(I\otimes Z^{i-2}\right)\left(I\otimes W\right)^{\rm T} \end{align*} and $$ {\it{\Psi}}_3=\left(W\otimes X\right)\left(\sum_{i=3}^{\infty}\alpha_i\sum_{j=2}^{i-1}\big((Z^{\rm T})^{i-j-1}\otimes Z^{j-2}\big)\right) \left(X\otimes W\right)^{\rm T}, $$ where the $$\{\alpha_i\}$$’s are the coefficients from the expression of $$f(x)$$. Proof. For any $$E\in\mathbb{R}^{n\times n}$$, from (2.7) and the expression of $$f(x)$$ we have \begin{eqnarray} L_f(\widetilde{A},E)&=&\sum_{i=1}^{\infty}\alpha_i\sum_{j=1}^{i}\widetilde{A}^{j-1}E\widetilde{A}^{i-j}\nonumber\\ &=&\alpha_1E+\alpha_2(E\widetilde{A}+\widetilde{A}E)+\sum_{i=3}^{\infty}\alpha_i \left(E\widetilde{A}^{i-1}+\sum_{j=2}^{i-1}\widetilde{A}^{j-1}E\widetilde{A}^{i-j}+\widetilde{A}^{i-1}E\right)\nonumber\\ &=&\alpha_1E+\sum_{i=2}^{\infty}\alpha_i(E\widetilde{A}^{i-1}+\widetilde{A}^{i-1}E)+\sum_{i=3}^{\infty}\alpha_i\sum_{j=2}^{i-1}\widetilde{A}^{j-1}E\widetilde{A}^{i-j}. \end{eqnarray} (3.12) As a result, \begin{equation} {\rm vec}(\alpha_1E)=(\alpha_1 I\otimes I){\rm vec}(E)={\it{\Psi}}_1 {\rm vec}(E). \end{equation} (3.13) If $$i\geq 2$$, we obtain from $$\widetilde{A}^{i-1}=XZ^{i-2}TY^{\rm T}=XZ^{i-2}W^{\rm T}~$$ that \begin{eqnarray} {\rm vec}\left(\sum_{i=2}^{\infty}\alpha_i(E\widetilde{A}^{i-1}+\widetilde{A}^{i-1}E)\right)&=&{\rm vec}\left(\sum_{i=2}^{\infty}\alpha_i\big(EXZ^{i-2}W^{\rm T}+XZ^{i-2}W^{\rm T}E\big)\right)\nonumber\\ &=&\sum_{i=2}^{\infty}\alpha_i\big(W(Z^{\rm T})^{i-2}X^{\rm T}\otimes I+I\otimes XZ^{i-2}W^{\rm T}\big){\rm vec}(E)\nonumber\\ &=&\left(\sum_{i=2}^{\infty}\alpha_i(W\otimes I)\big((Z^{\rm T})^{i-2}\otimes I\big)(X\otimes I)^{\rm T}\right.\nonumber\\ &&\left.+\sum_{i=2}^{\infty}\alpha_i(I\otimes X)\big(I\otimes Z^{i-2}\big)(I\otimes W)^{\rm T}\big]\right){\rm vec}(E)\nonumber\\ &=&{\it{\Psi}}_2{\rm vec}(E). \end{eqnarray} (3.14) Moreover, \begin{eqnarray} {\rm vec}\left(\sum_{i=3}^{\infty}\alpha_i\sum_{j=2}^{i-1}\widetilde{A}^{j-1}E\widetilde{A}^{i-j}\right)&=&\sum_{i=3}^{\infty}\alpha_i\sum_{j=2}^{i-1}{\rm vec}\big(XZ^{j-2}W^{\rm T} E XZ^{i-j-1}W^{\rm T}\big)\nonumber\\ &=&\sum_{i=3}^{\infty}\alpha_i\sum_{j=2}^{i-1}\big((W(Z^{\rm T})^{i-j-1}X^{\rm T}) \otimes (XZ^{j-2}W^{\rm T})\big){\rm vec}(E)\nonumber\\ &=&\sum_{i=3}^{\infty}\alpha_i\sum_{j=2}^{i-1}(W\otimes X)\big((Z^{\rm T})^{i-j-1}\otimes Z^{j-2}\big) (X\otimes W)^{\rm T}{\rm vec}(E)\nonumber\\ &=&{\it{\Psi}}_3{\rm vec}(E), \end{eqnarray} (3.15) and (3.11) follows from (3.12)–(3.15). □ Remark 3.6 We note that \begin{equation} K_g(Z)=\sum_{i=2}^{\infty}\alpha_i\sum_{j=1}^{i-1}\big((Z^{\rm T})^{i-j-1}\otimes Z^{j-1}\big) \end{equation} (3.16) (see Higham, 2008a, Problem 3.6). As a result, Theorem 3.5 indicates that $$L_f(\widetilde{A})$$ and $$L_g(Z)$$ are closely related. More precisely, if we define $$ \psi_i(Z)=\sum_{j=1}^{i-1}\big((Z^{\rm T})^{i-j-1}\otimes Z^{j-1}\big), \quad i=2,3,\ldots, $$ then $$K_g(Z)=\sum_{i=2}^{\infty}\alpha_i\psi_i(Z)$$ and $${\it{\Psi}}_3=\big(W\otimes X\big)\big(\sum_{i=3}^{\infty}\alpha_i\psi_{i-1}(Z)\big)\big(X\otimes W\big)^{\rm T}$$. 4. Numerical experiments In this section, we perform some numerical experiments to illustrate the numerical behavior of Algorithm 1 and Algorithm 2. All the numerical experiments were run on a Dell PC with eight cores Intel(R) Core(TM)i7-2600 processor with 3.4 GHz and 16.0 GB of RAM. As an operating system, we used Windows 7 with 64-bit architecture. All the numerical results were obtained from MATLAB R2015b implementations with unit roundoff $$\epsilon\approx 2.22\times 10^{-16}$$ (Higham, 2002, Chapter 2). In all the examples, the SCR of $$A$$ is computed using the MATLAB functions $${\tt scra.m}$$ and $${\tt spqr.m}$$ due to G.W. Stewart,1 where the tolerance $${\tt tol}$$ is taken as $$\epsilon_{\rm col}=\epsilon_{\rm row}=10^{-5}$$. The matrix exponential is calculated by using the MATLAB built-in function $${\tt expm.m}$$, while the $$\varphi_{\ell}\,(\ell\geq 1)$$ functions are computed using the $${\tt phipade.m}$$ function in the MATLAB package EXPINT (Berland et al., 2007). As was pointed out in Remark 2.3, if $$A$$ is large and sparse then $$\varphi_{\ell}(A)$$ can be dense and there is no need to form or store them in practice. In this section, however, we form the matrix functions $$\varphi_{\ell}(\widetilde{A})$$, so that the relative and the absolute errors of $$\varphi_{\ell}(\widetilde{A})$$ with respect to $$\varphi_{\ell}(A)$$ can be explicitly computed for $$\ell=0,1,\ldots,p$$. 4.1 An application to data dimensionality reduction In this example, we show the efficiency of our new method for computing the matrix exponential of large scale-and low-rank matrices. Many data-mining problems involve data sets represented in very high-dimensional spaces. In order to handle high-dimensional data, the dimensionality needs to be reduced. Linear discriminant analysis (LDA) is one of notable subspace transformation methods for dimensionality reduction (Duda et al., 2000). LDA encodes discriminant information by maximizing the between-class scatter while minimizing the within-class scatter in the projected subspace. Let $$X=[{\bf a}_{1},{\bf a}_{2},\ldots,{\bf a}_{m}]$$ be a set of training samples in an $$n$$-dimensional feature space, and assume that the original data is partitioned into $$K$$ classes as $$X=[X_1,X_2,\ldots,X_K]$$. We denote by $$m_{j}$$ the number of samples in the $$j$$th class, and thus $$\sum_{j=1}^{K}m_{j}=m$$. Let $${\bf c}_{j}$$ be the centroid of the $$j$$th class, and $${\bf c}$$ be the global centroid of the training data set. If we define $${\rm 1}\kern-0.24em{\rm I}_{j}=[1,1,\ldots,1]^{\rm T} \in \mathbb{R}^{m_{j}}$$ then the within-class scatter matrix is defined as $$ S_{\rm W}=\sum_{j=1}^{K}\sum_{{\bf a}_{i}\in X_{j}}({\bf a}_{i}-{\bf c}_{j})({\bf a}_{i}-{\bf c}_{j})^{\rm T}=H_WH_W^{\rm T}, $$ where $$H_{\rm W}=[X_{1}-{\bf c}_{1}{\rm 1}\kern-0.24em{\rm I}_{1}^{\rm T},\ldots,X_{K}-{\bf c}_{K}{\rm 1}\kern-0.24em{\rm I}_{K}^{\rm T}]\in\mathbb{R}^{n\times m}$$. The between-class scatter matrix is defined as $$ S_{\rm B}=\sum_{j=1}^{K}n_{j}({\bf c}_{j}-{\bf c})({\bf c}_{j}-{\bf c})^{\rm T}=H_BH_B^{\rm T}, $$ where $$H_B=[\sqrt{n_{1}}({\bf c}_{1}-{\bf c}),\sqrt{n_{2}}({\bf c}_{2}-{\bf c}),\ldots,\sqrt{n_{K}}({\bf c}_{K}-{\bf c})]\in\mathbb{R}^{n \times K}$$. The optimal projection matrix in the LDA method can be obtained by solving the large-scale generalized eigenproblem \begin{equation} S_{\rm B}{\bf x}= \lambda S_{\rm W}{\bf x}. \end{equation} (4.1) In practice, the dimension of real data usually exceeds the number of training samples (i.e., $$n\gg m$$), which is called the small-sample size (SSS) or the undersampled problem (Duda et al. 2000; Park & Park 2008). It is an intrinsic limitation of the classical LDA method and is also a common problem in classification applications (Park & Park, 2008). Note that both $$S_W$$ and $$S_B$$ are singular. Indeed, if $$X$$ is linearly independent then $${\rm rank}(S_B)=K-1$$ and $${\rm rank}(S_W)=m-K$$ (Duda et al., 2000; Wu & Feng, 2015). Thus, the ranks of $$S_B$$ and $$S_W$$ are much smaller than the dimensionality $$n$$. In other words, the SSS problem stems from generalized eigenproblems with singular matrices. To rectify this drawback, a novel method based on the matrix exponential, called the EDA method, was proposed in Zhang et al. (2010). Instead of (4.1), the EDA method solves the following generalized matrix exponential eigenproblem (Zhang et al., 2010; Wu et al., 2017): \begin{equation} \textrm{exp}(S_{\rm B}){\bf x}= \lambda\,\textrm{exp}(S_{\rm W}){\bf x}. \end{equation} (4.2) As both $$\textrm{exp}(S_{\rm W})$$ and $$\exp(S_B)$$ are symmetric positive definite (SPD), this is a definite generalized eigenproblem, and there is a complete set of eigenvectors (Davies et al., 2001; Higham, 2009); moreover, all the eigenvalues are real. The EDA method is described as follows; for more details, refer to Zhang et al. (2010). Algorithm 3 Zhang et al. (2010) The exponential discriminant analysis (EDA). Input: The data matrix $$X=[{\bf a}_{1},{\bf a}_{2},\ldots,{\bf a}_{m}]\in\mathbb{R}^{n\times m}$$, where $${\bf a}_{j}$$ represents the $$j$$th training image. Output: The projection matrix $$V$$. 1. Compute the matrices $$S_{\rm B}$$, $$S_{\rm W}$$, $${\rm exp}(S_{\rm B})$$, and $${\rm exp}(S_{\rm W})$$. 2. Compute the eigenvalues $$\{\lambda_{i}\}$$ and associated eigenvectors $$\{{\bf x}_{i}\}$$ of $${\rm exp}(S_{\rm W})^{-1}{\rm exp}(S_{\rm B})$$. 3. Sort the eigenvectors $$V=\{{\bf x}_{i}\}$$ according to $$\{\lambda_{i}\}$$ in decreasing order. 4. Orthogonalize the columns of the projection matrix $$V$$. The framework of the EDA method for dimensionality reduction has gained wide attention in recent years (Yan & Pan, 2010; Zhang et al., 2010; Wang et al., 2011; Ahmed, 2013; Dornaika & Bosaghzadeh, 2013; Wang et al., 2014). However, the time complexity of Algorithm 3 is dominated by the computation of $${\rm exp}(S_B)$$ and $${\rm exp}(S_W)$$, which is demanding since the data dimension $$n$$ is large (Zhang et al., 2010). Moreover, the problem of (4.2) should not be solved as in Algorithm 3 and not by the inverting transformation shown. One refers to a recent article on inexact Krylov subspace methods for solving this type of generalized eigenproblem (Wu et al., 2017). By Theorem 2.1, we can compute the large matrix exponentials as follows. Corollary 4.1 Under the above notation, we have \begin{equation} \exp(S_B)=I+H_B\big(\varphi_1(H_B^{\rm T}H_B)\big) H_B^{\rm T} \end{equation} (4.3) and \begin{equation} \exp(S_W)=I+H_W\big(\varphi_1(H_W^{\rm T}H_W)\big) H_W^{\rm T}. \end{equation} (4.4) So we have the following algorithm for the matrix EDA. Algorithm 4 An algorithm for approximating $$\exp(S_B)$$ and $$\exp(S_W)$$. 1. Given the data matrix $$X=[{\bf a}_{1},{\bf a}_{2},\ldots,{\bf a}_{m}]\in\mathbb{R}^{n\times m}$$, form $$H_B$$ and $$H_W$$. 2. Compute $$\varphi_1(H_B^{\rm T}H_B)$$ and $$\varphi_1(H_W^{\rm T}H_W)$$. 3. Store $$H_B,H_W$$ and $$\varphi_1(H_B^{\rm T}H_B)$$, $$\varphi_1(H_W^{\rm T}H_W)$$ for $$\exp(S_B)$$ and $$\exp(S_W)$$. If desired, form $$\exp(S_B)$$ and $$\exp(S_W)$$ in terms of (4.3) and (4.4). Note that both $$H_B$$ and $$H_W$$ are already available in step 1, and so there is no need to perform rank-revealing decompositions of $$S_B$$ and $$S_W$$. As a result, the computation of the two $$n\times n$$ matrix exponentials $$\exp(S_B),\exp(S_W)$$ reduces to that of $$\varphi_1(H_B^{\rm T}H_B)\in\mathbb{R}^{K\times K}$$ and $$\varphi_1(H_W^{\rm T}H_W)\in\mathbb{R}^{m\times m}$$, where $$K,m\ll n$$. Next we illustrate the efficiency of Algorithm 4 for matrix EDA. There are three real-world databases in this example. The first one is the $${\tt ORL}$$ database2 that contains 400 face images of $$40$$ individuals, and the original image size is $$92\times 112=10304$$. The second test set is the $${\tt Yale}$$ face database taken from the Yale Center for Computational Vision and Control.3 It contains $$165$$ grayscale images of 15 individuals. The original image size is $$320\times 243=77760$$. The third test set is the $${\tt Extended\,YaleB}$$ database.4 This database contains 5760 single light source images of 10 subjects, each seen under 576 viewing conditions. A subset of $$38$$ classes with 2432 images are used in this example, 64 images per individual with illumination. In the $${\tt ORL}$$ database, the images are aligned based on eye coordinates and are cropped and scaled to $$n=32\times 32$$ and $$64\times 64$$, respectively; the original image size with $$n=92\times 112$$ is also considered. In the $${\tt Yale}$$ and the $${\tt Extended\,YaleB}$$ databases, all images are aligned based on eye coordinates and are cropped and scaled to $$n=32\times 32,~64\times 64$$ and $$100\times 100$$, respectively. In this example, a random subset with $$3$$ images per subject is taken to form the training set, and the rest of the images are used as the testing set. Each column of the data matrices is scaled by its 2-norm. In Algorithm 4, the CPU time consists of computing $$H_B,H_W$$, $$\varphi_1(H_B^{\rm T}H_B)$$ and $$\varphi_1(H_W^{\rm T}H_W)$$ and forming $$\exp(S_B)$$ and $$\exp(S_W)$$ in terms of (4.3) and (4.4). In the original EDA algorithm (Algorithm 3), the CPU time consists of forming $$H_B,H_W$$ and the computation of $$\exp(S_B)$$ and $$\exp(S_W)$$ by using the MATLAB built-in function $${\tt expm.m}$$. Let $$\exp(S_B),\exp(S_W)$$ be the ‘exact solutions’ obtained from running $${\tt expm.m}$$, and let $$\widetilde{\exp(S_B)},\widetilde{\exp(S_W)}$$ be the approximations obtained from (4.3) and (4.4). We define $$ {\bf Err_B}=\frac{\|\exp(S_B)-\widetilde{\exp(S_B)}\|_F}{\|\exp(S_B)\|_F}\quad {\rm and} \quad {\bf Err_W}=\frac{\|\exp(S_W)-\widetilde{{\rm exp}(S_W)}\|_F}{\|\exp(S_W)\|_F} $$ as the relative errors of the approximations $$\widetilde{{\rm exp}(S_B)},\widetilde{{\rm exp}(S_W)}$$, respectively, and we denote by $$ {\tt Rel\_ErrF}=\max({\bf Err_B},{\bf Err_W}) $$ the maximal value of the two relative errors. Table 1 lists the CPU time in seconds of Algorithm 4, $${\tt expm.m}$$ and the values of the maximal relative errors $${\tt Rel\_ErrF}$$. Table 1 Example $$1$$: CPU time in seconds and the relative errors for approximating $$\exp(S_B)$$ and $$\exp(S_W)$$ $${\tt Database}$$ $${\tt n}$$ $${\tt Algorithm\,4}$$ $${\tt expm.m}$$ $${\tt Rel\_ErrF}$$ ORL 1024 0.08 0.26 $$2.20\times 10^{-15}$$ $$4096$$ 0.26 17.3 $$3.22\times 10^{-15}$$ $$10304$$ 1.28 261.1 $$4.49\times 10^{-15}$$ Yale $$1024$$ 0.08 0.25 $$2.11\times 10^{-15}$$ $$4096$$ 0.24 17.3 $$3.10\times 10^{-15}$$ $$10000$$ 1.13 238.5 $$4.35\times 10^{-15}$$ Extended YaleB $$1024$$ 0.08 0.27 $$2.13\times 10^{-15}$$ $$4096$$ 0.27 17.3 $$3.14\times 10^{-15}$$ $$10000$$ 1.22 238.6 $$4.50\times 10^{-15}$$ $${\tt Database}$$ $${\tt n}$$ $${\tt Algorithm\,4}$$ $${\tt expm.m}$$ $${\tt Rel\_ErrF}$$ ORL 1024 0.08 0.26 $$2.20\times 10^{-15}$$ $$4096$$ 0.26 17.3 $$3.22\times 10^{-15}$$ $$10304$$ 1.28 261.1 $$4.49\times 10^{-15}$$ Yale $$1024$$ 0.08 0.25 $$2.11\times 10^{-15}$$ $$4096$$ 0.24 17.3 $$3.10\times 10^{-15}$$ $$10000$$ 1.13 238.5 $$4.35\times 10^{-15}$$ Extended YaleB $$1024$$ 0.08 0.27 $$2.13\times 10^{-15}$$ $$4096$$ 0.27 17.3 $$3.14\times 10^{-15}$$ $$10000$$ 1.22 238.6 $$4.50\times 10^{-15}$$ Table 1 Example $$1$$: CPU time in seconds and the relative errors for approximating $$\exp(S_B)$$ and $$\exp(S_W)$$ $${\tt Database}$$ $${\tt n}$$ $${\tt Algorithm\,4}$$ $${\tt expm.m}$$ $${\tt Rel\_ErrF}$$ ORL 1024 0.08 0.26 $$2.20\times 10^{-15}$$ $$4096$$ 0.26 17.3 $$3.22\times 10^{-15}$$ $$10304$$ 1.28 261.1 $$4.49\times 10^{-15}$$ Yale $$1024$$ 0.08 0.25 $$2.11\times 10^{-15}$$ $$4096$$ 0.24 17.3 $$3.10\times 10^{-15}$$ $$10000$$ 1.13 238.5 $$4.35\times 10^{-15}$$ Extended YaleB $$1024$$ 0.08 0.27 $$2.13\times 10^{-15}$$ $$4096$$ 0.27 17.3 $$3.14\times 10^{-15}$$ $$10000$$ 1.22 238.6 $$4.50\times 10^{-15}$$ $${\tt Database}$$ $${\tt n}$$ $${\tt Algorithm\,4}$$ $${\tt expm.m}$$ $${\tt Rel\_ErrF}$$ ORL 1024 0.08 0.26 $$2.20\times 10^{-15}$$ $$4096$$ 0.26 17.3 $$3.22\times 10^{-15}$$ $$10304$$ 1.28 261.1 $$4.49\times 10^{-15}$$ Yale $$1024$$ 0.08 0.25 $$2.11\times 10^{-15}$$ $$4096$$ 0.24 17.3 $$3.10\times 10^{-15}$$ $$10000$$ 1.13 238.5 $$4.35\times 10^{-15}$$ Extended YaleB $$1024$$ 0.08 0.27 $$2.13\times 10^{-15}$$ $$4096$$ 0.27 17.3 $$3.14\times 10^{-15}$$ $$10000$$ 1.22 238.6 $$4.50\times 10^{-15}$$ We observe from Table 1 that Algorithm 4 runs much faster than $${\tt expm.m}$$, especially when $$n$$ is large. For instance, when the dimensionality of the data sets is around $$10^4$$, $${\tt expm.m}$$ requires about 240 s, while our new method needs only about 1.2 s, which is a significant improvement. Furthermore, the relative errors of our approximations are of the order of $$\mathcal{O}(10^{-15})$$. As a result, our new method is very efficient and reliable for solving the large matrix exponential problems arising in the EDA framework. 4.2 Computing approximations to $$\varphi$$-functions of large matrices with low rank or with fast decaying singular values We show the efficiency of Algorithm 1 for consecutively approximating several $$\varphi$$-functions of $$A$$ with low rank or with fast decaying singular values. The test matrices are available from Davis & Hu (2011) and Tspras, and Table 2 lists the problem characteristics and the numerical rank, as well as the error (measured by the Frobenius norm) obtained from the SCR. Table 2 Example $$2$$: Problem characteristics of the test matrices, $$r$$, and the error from the SCR approximation Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_F$$ Description $${\tt man5976}$$ 5976 2341 $$2.83\times 10^{-9}$$ Structural problem $${\tt Movies}$$ 5757 1053 $$2.01\times 10^{-11}$$ Directed network $${\tt lock3491}$$ 3491 603 $$1.96\times 10^{-8}$$ Structural problem $${\tt cegb3306}$$ 3306 537 $$7.76\times 10^{-11}$$ Finite element framework $${\tt zenios}$$ 2873 250 $$7.75\times 10^{-6}$$ Optimization problem $${\tt watt\_1}$$ 1856 151 $$9.97\times 10^{-6}$$ Computational fluid dynamics $${\tt watt\_2}$$ 1856 151 $$9.98\times 10^{-6}$$ Computational fluid dynamics $${\tt eris1176}$$ 1176 774 $$2.49\times 10^{-10}$$ Power network problem Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_F$$ Description $${\tt man5976}$$ 5976 2341 $$2.83\times 10^{-9}$$ Structural problem $${\tt Movies}$$ 5757 1053 $$2.01\times 10^{-11}$$ Directed network $${\tt lock3491}$$ 3491 603 $$1.96\times 10^{-8}$$ Structural problem $${\tt cegb3306}$$ 3306 537 $$7.76\times 10^{-11}$$ Finite element framework $${\tt zenios}$$ 2873 250 $$7.75\times 10^{-6}$$ Optimization problem $${\tt watt\_1}$$ 1856 151 $$9.97\times 10^{-6}$$ Computational fluid dynamics $${\tt watt\_2}$$ 1856 151 $$9.98\times 10^{-6}$$ Computational fluid dynamics $${\tt eris1176}$$ 1176 774 $$2.49\times 10^{-10}$$ Power network problem Table 2 Example $$2$$: Problem characteristics of the test matrices, $$r$$, and the error from the SCR approximation Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_F$$ Description $${\tt man5976}$$ 5976 2341 $$2.83\times 10^{-9}$$ Structural problem $${\tt Movies}$$ 5757 1053 $$2.01\times 10^{-11}$$ Directed network $${\tt lock3491}$$ 3491 603 $$1.96\times 10^{-8}$$ Structural problem $${\tt cegb3306}$$ 3306 537 $$7.76\times 10^{-11}$$ Finite element framework $${\tt zenios}$$ 2873 250 $$7.75\times 10^{-6}$$ Optimization problem $${\tt watt\_1}$$ 1856 151 $$9.97\times 10^{-6}$$ Computational fluid dynamics $${\tt watt\_2}$$ 1856 151 $$9.98\times 10^{-6}$$ Computational fluid dynamics $${\tt eris1176}$$ 1176 774 $$2.49\times 10^{-10}$$ Power network problem Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_F$$ Description $${\tt man5976}$$ 5976 2341 $$2.83\times 10^{-9}$$ Structural problem $${\tt Movies}$$ 5757 1053 $$2.01\times 10^{-11}$$ Directed network $${\tt lock3491}$$ 3491 603 $$1.96\times 10^{-8}$$ Structural problem $${\tt cegb3306}$$ 3306 537 $$7.76\times 10^{-11}$$ Finite element framework $${\tt zenios}$$ 2873 250 $$7.75\times 10^{-6}$$ Optimization problem $${\tt watt\_1}$$ 1856 151 $$9.97\times 10^{-6}$$ Computational fluid dynamics $${\tt watt\_2}$$ 1856 151 $$9.98\times 10^{-6}$$ Computational fluid dynamics $${\tt eris1176}$$ 1176 774 $$2.49\times 10^{-10}$$ Power network problem In this example, we compare Algorithm 1 with $${\tt expm.m/phipade.m}$$, that is, $${\tt expm.m}$$ for the matrix exponential $$\exp(A)$$ and $${\tt phipade.m}$$ for $$\varphi_{\ell}(A),~\ell=1,2,3,4$$. In Algorithm 1, the CPU time consists of computing the SCR and the evaluation of $$\varphi_{\ell}(Z)~(\ell=1,2,3,4,5)$$ by using $${\tt phipade.m}$$, as well as forming $$\varphi_{\ell}(A)$$ in terms of (2.3), $$\ell=0,1,2,3,4$$. In $${\tt expm.m}$$/$${\tt phipade.m}$$, the CPU time consists of computing $$\varphi_{\ell}(A)$$ by using $${\tt expm.m}$$ ($$\ell=0$$) and $${\tt phipade.m}$$ ($$\ell=1,2,3,4$$). In order to measure the accuracy of the solutions, we define the maximal relative error as $$ {\tt Rel\_ErrF}=\max_{0\leq\ell\leq 4}\frac{\|\varphi_{\ell}(A)-{\varphi_{\ell}(\widetilde{A})}\|_F}{\|\varphi_{\ell}(A)\|_F}, $$ where $$\varphi_{\ell}(A)$$ is the ‘exact solution’ obtained from $${\tt expm.m}$$ for $$\ell=0$$ and $${\tt phipade.m}$$ for $$\ell=1,2,3,4$$, and $${\varphi_{\ell}(\widetilde{A})}$$ is the approximation obtained from running Algorithm 1. Table 3 lists the numerical results. Table 3 Example $$2$$: CPU time in seconds and the maximal relative error for approximating the $$\varphi_{\ell}$$-functions, $$\ell=0,1,2,3,4$$, where ‘$$^{*}$$’ denotes that we compute $$\varphi_{\ell}(-A)$$ instead of $$\varphi_{\ell}(A)$$ $${\tt Test matrix}$$ $${\tt Algorithm\,1}$$ $${\tt expm.m/phipade.m}$$ $${\tt Rel\_ErrF}$$ $${\tt man5976}$$$$^{*}$$ 205.5 10314.0 $$1.10\times 10^{-12}$$ $${\tt Movies}$$ 360.6 456.1 $$4.32\times 10^{-14}$$ $${\tt lock}$$3491$$^{*}$$ 5.49 2469.6 $$5.12\times 10^{-13}$$ $${\tt cegb}$$3306$$^{*}$$ 3.11 1599.7 $$1.24\times 10^{-13}$$ $${\tt zenios}$$ 0.63 3.63 $$1.17\times 10^{-7}$$ $${\tt watt}$$_1 0.19 30.4 $$1.93\times 10^{-7}$$ $${\tt watt}$$_2 0.19 39.1 $$2.02\times 10^{-7}$$ $${\tt eris1176}$$$$^{*}$$ 6.74 85.7 $$2.86\times 10^{-12}$$ $${\tt Test matrix}$$ $${\tt Algorithm\,1}$$ $${\tt expm.m/phipade.m}$$ $${\tt Rel\_ErrF}$$ $${\tt man5976}$$$$^{*}$$ 205.5 10314.0 $$1.10\times 10^{-12}$$ $${\tt Movies}$$ 360.6 456.1 $$4.32\times 10^{-14}$$ $${\tt lock}$$3491$$^{*}$$ 5.49 2469.6 $$5.12\times 10^{-13}$$ $${\tt cegb}$$3306$$^{*}$$ 3.11 1599.7 $$1.24\times 10^{-13}$$ $${\tt zenios}$$ 0.63 3.63 $$1.17\times 10^{-7}$$ $${\tt watt}$$_1 0.19 30.4 $$1.93\times 10^{-7}$$ $${\tt watt}$$_2 0.19 39.1 $$2.02\times 10^{-7}$$ $${\tt eris1176}$$$$^{*}$$ 6.74 85.7 $$2.86\times 10^{-12}$$ Table 3 Example $$2$$: CPU time in seconds and the maximal relative error for approximating the $$\varphi_{\ell}$$-functions, $$\ell=0,1,2,3,4$$, where ‘$$^{*}$$’ denotes that we compute $$\varphi_{\ell}(-A)$$ instead of $$\varphi_{\ell}(A)$$ $${\tt Test matrix}$$ $${\tt Algorithm\,1}$$ $${\tt expm.m/phipade.m}$$ $${\tt Rel\_ErrF}$$ $${\tt man5976}$$$$^{*}$$ 205.5 10314.0 $$1.10\times 10^{-12}$$ $${\tt Movies}$$ 360.6 456.1 $$4.32\times 10^{-14}$$ $${\tt lock}$$3491$$^{*}$$ 5.49 2469.6 $$5.12\times 10^{-13}$$ $${\tt cegb}$$3306$$^{*}$$ 3.11 1599.7 $$1.24\times 10^{-13}$$ $${\tt zenios}$$ 0.63 3.63 $$1.17\times 10^{-7}$$ $${\tt watt}$$_1 0.19 30.4 $$1.93\times 10^{-7}$$ $${\tt watt}$$_2 0.19 39.1 $$2.02\times 10^{-7}$$ $${\tt eris1176}$$$$^{*}$$ 6.74 85.7 $$2.86\times 10^{-12}$$ $${\tt Test matrix}$$ $${\tt Algorithm\,1}$$ $${\tt expm.m/phipade.m}$$ $${\tt Rel\_ErrF}$$ $${\tt man5976}$$$$^{*}$$ 205.5 10314.0 $$1.10\times 10^{-12}$$ $${\tt Movies}$$ 360.6 456.1 $$4.32\times 10^{-14}$$ $${\tt lock}$$3491$$^{*}$$ 5.49 2469.6 $$5.12\times 10^{-13}$$ $${\tt cegb}$$3306$$^{*}$$ 3.11 1599.7 $$1.24\times 10^{-13}$$ $${\tt zenios}$$ 0.63 3.63 $$1.17\times 10^{-7}$$ $${\tt watt}$$_1 0.19 30.4 $$1.93\times 10^{-7}$$ $${\tt watt}$$_2 0.19 39.1 $$2.02\times 10^{-7}$$ $${\tt eris1176}$$$$^{*}$$ 6.74 85.7 $$2.86\times 10^{-12}$$ It is seen from Table 3 that Algorithm 1 works quite well, and we provide a competitive candidate for consecutively approximating several $$\varphi$$-functions of large sparse matrices with low rank or with fast decaying singular values. Two remarks are given. First, Algorithm 1 runs much faster than $${\tt expm.m}$$/$${\tt phipade.m}$$. For instance, 205.5 s vs. 10314.0 s for the $${\tt man5976}$$ matrix, 5.49 s vs. 2469.6 s for the $${\tt lock3491}$$ matrix, and 3.11 s vs. 1599.7 s for the $${\tt cegb3306}$$ matrix. Second, the accuracy of our approximations is satisfactory in most cases. However, for the $${\tt zenios}$$, $${\tt watt\_1}$$ and $${\tt watt\_2}$$ matrices, the relative errors $${\tt Rel\_ErrF}$$ are of the order of $$\mathcal{O}(10^{-7})$$. Let us discuss it in detail. In Fig. 1, we plot the singular values of the $${\tt eris1176}$$ matrix and the $${\tt watt\_1}$$ matrix. It is observed that the $${\tt eris1176}$$ matrix has (relatively) faster decaying singular values than the $${\tt watt\_1}$$ matrix. Indeed, the errors $$\|A-\widetilde{A}\|_F$$ from the SCR approximations, with respect to the three matrices $${\tt zenios}$$, $${\tt watt\_1}$$ and $${\tt watt\_2}$$, are about $$7.75\times 10^{-6}$$, $$9.97\times 10^{-6}$$ and $$9.98\times 10^{-6}$$, respectively, while that of the $${\tt eris1176}$$ matrix is about $$2.49\times 10^{-10}$$. In Table 4, we provide the relative F-norm condition numbers for the $$\varphi_{\ell}$$-functions of the four matrices, which are computed by using the $${\tt funm\_condest\_fro.m}$$ function in the Matrix Function Toolbox (Higham, 2008b). One can see that the maximal condition numbers with respect to the four matrices are 0.98, 0.25, 1.71 and 60.0. By Theorem 2.6, the accuracy of the computed solution of the $${\tt eris1176}$$ matrix is about $$1.87\times 10^{-10}$$, while that of the $${\tt watt\_1}$$ matrix is about $$2.49\times 10^{-6}$$. Hence, the accuracy with respect to the $${\tt eris1176}$$ matrix is higher than that of the other three matrices. Further, it is possible to get a smaller $${\tt Rel\_ErrF}$$ if a larger $$r$$ is used in our new method, since the error $$\|A-\widetilde{A}\|_F$$ will be smaller. Consequently, our new method is suitable for $$\varphi$$-functions of large matrices with low rank or with fast decaying singular values. Fig. 1. View largeDownload slide Example 2: Singular values of the $${\tt eris1176}$$ and the $${\tt watt\_1}$$ matrices. Fig. 1. View largeDownload slide Example 2: Singular values of the $${\tt eris1176}$$ and the $${\tt watt\_1}$$ matrices. Table 4 Example $$2$$: Relative F-norm condition numbers for the $$\varphi_{\ell}$$-functions of the four matrices $${\tt zenios}$$, $${\tt watt\_1}$$, $${\tt watt\_2}$$ and $${\tt eris1176}$$, $$\ell=0,1,2,3,4$$ $$\ell$$ $${\tt zenios}$$ $${\tt watt\_1}$$ $${\tt watt\_2}$$ $${\tt eris1176}$$ 0 $$1.32\times 10^{-2}$$ $$1.66\times 10^{-2}$$ $$2.02\times 10^{-2}$$ 0.12 1 0.98 0.25 1.71 60.0 2 0.40 0.14 0.75 25.2 3 0.19 $$9.34\times 10^{-2}$$ 0.41 11.8 4 0.13 $$6.83\times 10^{-2}$$ 0.26 4.44 $$\ell$$ $${\tt zenios}$$ $${\tt watt\_1}$$ $${\tt watt\_2}$$ $${\tt eris1176}$$ 0 $$1.32\times 10^{-2}$$ $$1.66\times 10^{-2}$$ $$2.02\times 10^{-2}$$ 0.12 1 0.98 0.25 1.71 60.0 2 0.40 0.14 0.75 25.2 3 0.19 $$9.34\times 10^{-2}$$ 0.41 11.8 4 0.13 $$6.83\times 10^{-2}$$ 0.26 4.44 Table 4 Example $$2$$: Relative F-norm condition numbers for the $$\varphi_{\ell}$$-functions of the four matrices $${\tt zenios}$$, $${\tt watt\_1}$$, $${\tt watt\_2}$$ and $${\tt eris1176}$$, $$\ell=0,1,2,3,4$$ $$\ell$$ $${\tt zenios}$$ $${\tt watt\_1}$$ $${\tt watt\_2}$$ $${\tt eris1176}$$ 0 $$1.32\times 10^{-2}$$ $$1.66\times 10^{-2}$$ $$2.02\times 10^{-2}$$ 0.12 1 0.98 0.25 1.71 60.0 2 0.40 0.14 0.75 25.2 3 0.19 $$9.34\times 10^{-2}$$ 0.41 11.8 4 0.13 $$6.83\times 10^{-2}$$ 0.26 4.44 $$\ell$$ $${\tt zenios}$$ $${\tt watt\_1}$$ $${\tt watt\_2}$$ $${\tt eris1176}$$ 0 $$1.32\times 10^{-2}$$ $$1.66\times 10^{-2}$$ $$2.02\times 10^{-2}$$ 0.12 1 0.98 0.25 1.71 60.0 2 0.40 0.14 0.75 25.2 3 0.19 $$9.34\times 10^{-2}$$ 0.41 11.8 4 0.13 $$6.83\times 10^{-2}$$ 0.26 4.44 4.3 Estimating the relative and the absolute condition numbers This example aims to demonstrate the efficiency of our new strategy for estimating the absolute and the relative 2-norm condition numbers of $$\varphi$$-functions. There are six test matrices that are available from Davis & Hu (2011) and Tspras. For the first three matrices, we compute the absolute and the relative condition numbers for a given $$\ell$$, while for the rest of the three matrices, we evaluate the condition numbers consecutively for $$\ell=0,1,2,3,4$$. Table 5 lists the problem characteristics the numerical rank, as well as the error (measured by the 2-norm) obtained from the SCR of these matrices. Table 5 Example $$3$$: Problem characteristics of the test matrices, $$r$$, and the error from the SCR approximation Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_2$$ $${\tt EVA}$$ 8497 1301 $$1.64\times 10^{-13}$$ $${\tt EPA}$$ 4772 951 $$3.57\times 10^{-12}$$ $${\tt eris1176}$$ 1176 774 $$1.85\times 10^{-10}$$ $${\tt Abortion}$$ 2293 633 $$9.41\times 10^{-12}$$ $${\tt Computational\_Geometry}$$ 1226 415 $$1.96\times 10^{-11}$$ $${\tt Computational\_Complexity}$$ 884 254 $$7.84\times 10^{-14}$$ Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_2$$ $${\tt EVA}$$ 8497 1301 $$1.64\times 10^{-13}$$ $${\tt EPA}$$ 4772 951 $$3.57\times 10^{-12}$$ $${\tt eris1176}$$ 1176 774 $$1.85\times 10^{-10}$$ $${\tt Abortion}$$ 2293 633 $$9.41\times 10^{-12}$$ $${\tt Computational\_Geometry}$$ 1226 415 $$1.96\times 10^{-11}$$ $${\tt Computational\_Complexity}$$ 884 254 $$7.84\times 10^{-14}$$ Table 5 Example $$3$$: Problem characteristics of the test matrices, $$r$$, and the error from the SCR approximation Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_2$$ $${\tt EVA}$$ 8497 1301 $$1.64\times 10^{-13}$$ $${\tt EPA}$$ 4772 951 $$3.57\times 10^{-12}$$ $${\tt eris1176}$$ 1176 774 $$1.85\times 10^{-10}$$ $${\tt Abortion}$$ 2293 633 $$9.41\times 10^{-12}$$ $${\tt Computational\_Geometry}$$ 1226 415 $$1.96\times 10^{-11}$$ $${\tt Computational\_Complexity}$$ 884 254 $$7.84\times 10^{-14}$$ Test matrix $$n$$ $$r$$ $$\|A-\widetilde{A}\|_2$$ $${\tt EVA}$$ 8497 1301 $$1.64\times 10^{-13}$$ $${\tt EPA}$$ 4772 951 $$3.57\times 10^{-12}$$ $${\tt eris1176}$$ 1176 774 $$1.85\times 10^{-10}$$ $${\tt Abortion}$$ 2293 633 $$9.41\times 10^{-12}$$ $${\tt Computational\_Geometry}$$ 1226 415 $$1.96\times 10^{-11}$$ $${\tt Computational\_Complexity}$$ 884 254 $$7.84\times 10^{-14}$$ We compare Algorithm 2 with the $${\tt funm\_condest1.m}$$ function in the Matrix Function Toolbox (Higham, 2008b). In Algorithm 2, the matrices $$R_1,R_2$$ are obtained from the QR decompositions of $$X$$ and $$Y$$, respectively, by using the MATLAB built-in function $${\tt qr.m}$$. The matrix $$L_{\varphi_{\ell+1}}(Z)$$ in (3.10) is evaluated by using the $${\tt funm\_condest\_fro.m}$$ function, and $$\|A\|_2$$ is evaluated by using the MATLAB built-in function $${\tt svds.m}$$. When $$\ell=0$$, the parameter ‘$${\tt fun}$$’ in $${\tt funm\_condest1.m}$$ is called by the MATLAB built-in function $${\tt expm.m}$$, whereas when $$\ell>0$$ this parameter is called by the $${\tt phipade.m}$$ function in the EXPINT package. The CPU time for Algorithm 2 is composed of computing the SCR and estimating $$\|A\|_2$$, as well as calculating (3.7) and (3.10). Tables 6–11 report the numerical results, where ‘$${\tt Relative\_Est}$$’ and ‘$${\tt Absolute\_Est}$$’ denote estimations to the relative and the absolute condition numbers, respectively. Table 6 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$ and the CPU time in seconds where ‘$$>$$3h’ denotes that the algorithm fails to converge within 3 h. The $${\tt EVA}$$ matrix, $$n=8497, r=1303$$ $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$7.74\times 10^3$$ $$1.83\times 10^{4}$$ 700.0 $${\tt Algorithm\,2}$$ $$1.07\times 10^3$$ $$1.32\times 10^{3}$$ 420.1 1 $${\tt funm\_condest1}$$ $$7.07\times 10^3$$ $$7.22\times 10^{3}$$ $$5846.1$$ $${\tt Algorithm\,2}$$ $$538.7$$ $$269.6$$ 425.2 2 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$321.7$$ 53.5 429.8 3 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$122.5$$ 5.08 440.9 4 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ 73.3 0.61 447.1 $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$7.74\times 10^3$$ $$1.83\times 10^{4}$$ 700.0 $${\tt Algorithm\,2}$$ $$1.07\times 10^3$$ $$1.32\times 10^{3}$$ 420.1 1 $${\tt funm\_condest1}$$ $$7.07\times 10^3$$ $$7.22\times 10^{3}$$ $$5846.1$$ $${\tt Algorithm\,2}$$ $$538.7$$ $$269.6$$ 425.2 2 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$321.7$$ 53.5 429.8 3 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$122.5$$ 5.08 440.9 4 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ 73.3 0.61 447.1 Table 6 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$ and the CPU time in seconds where ‘$$>$$3h’ denotes that the algorithm fails to converge within 3 h. The $${\tt EVA}$$ matrix, $$n=8497, r=1303$$ $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$7.74\times 10^3$$ $$1.83\times 10^{4}$$ 700.0 $${\tt Algorithm\,2}$$ $$1.07\times 10^3$$ $$1.32\times 10^{3}$$ 420.1 1 $${\tt funm\_condest1}$$ $$7.07\times 10^3$$ $$7.22\times 10^{3}$$ $$5846.1$$ $${\tt Algorithm\,2}$$ $$538.7$$ $$269.6$$ 425.2 2 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$321.7$$ 53.5 429.8 3 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$122.5$$ 5.08 440.9 4 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ 73.3 0.61 447.1 $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$7.74\times 10^3$$ $$1.83\times 10^{4}$$ 700.0 $${\tt Algorithm\,2}$$ $$1.07\times 10^3$$ $$1.32\times 10^{3}$$ 420.1 1 $${\tt funm\_condest1}$$ $$7.07\times 10^3$$ $$7.22\times 10^{3}$$ $$5846.1$$ $${\tt Algorithm\,2}$$ $$538.7$$ $$269.6$$ 425.2 2 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$321.7$$ 53.5 429.8 3 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ $$122.5$$ 5.08 440.9 4 $${\tt funm\_condest1}$$ — — $$>$$3h $${\tt Algorithm\,2}$$ 73.3 0.61 447.1 Table 7 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A)$$, $$\ell=0,1,2,3,4$$ and the CPU time in seconds. The $${\tt EPA}$$ matrix, $$n=4772, r=951$$ $$\ell$$ Method Rel_Est Abs_Est CPU 0 $${\tt funm\_condest1}$$ $$9.90\times 10^{3}$$ $$2.60\times 10^{4}$$ 90.1 $${\tt Algorithm\,2}$$ $$3.34\times 10^4$$ $$1.52\times 10^{5}$$ 95.4 1 $${\tt funm\_condest1}$$ $$1.09\times 10^4$$ $$8.05\times 10^{3}$$ 272.1 $${\tt Algorithm\,2}$$ $$2.25\times 10^4$$ $$2.74\times 10^{4}$$ 100.6 2 $${\tt funm\_condest1}$$ $$8.95\times 10^3$$ $$2.03\times 10^{3}$$ 443.7 $${\tt Algorithm\,2}$$ $$1.47\times 10^4$$ $$4.27\times 10^3$$ 105.5 3 $${\tt funm\_condest1}$$ $$7.86\times 10^3$$ $$425.9$$ 773.0 $${\tt Algorithm\,2}$$ $$9.71\times 10^3$$ $$587.4$$ 110.7 4 $${\tt funm\_condest1}$$ $$7.17\times 10^3$$ 75.2 1357.7 $${\tt Algorithm\,2}$$ $$5.29\times 10^3$$ 57.6 116.7 $$\ell$$ Method Rel_Est Abs_Est CPU 0 $${\tt funm\_condest1}$$ $$9.90\times 10^{3}$$ $$2.60\times 10^{4}$$ 90.1 $${\tt Algorithm\,2}$$ $$3.34\times 10^4$$ $$1.52\times 10^{5}$$ 95.4 1 $${\tt funm\_condest1}$$ $$1.09\times 10^4$$ $$8.05\times 10^{3}$$ 272.1 $${\tt Algorithm\,2}$$ $$2.25\times 10^4$$ $$2.74\times 10^{4}$$ 100.6 2 $${\tt funm\_condest1}$$ $$8.95\times 10^3$$ $$2.03\times 10^{3}$$ 443.7 $${\tt Algorithm\,2}$$ $$1.47\times 10^4$$ $$4.27\times 10^3$$ 105.5 3 $${\tt funm\_condest1}$$ $$7.86\times 10^3$$ $$425.9$$ 773.0 $${\tt Algorithm\,2}$$ $$9.71\times 10^3$$ $$587.4$$ 110.7 4 $${\tt funm\_condest1}$$ $$7.17\times 10^3$$ 75.2 1357.7 $${\tt Algorithm\,2}$$ $$5.29\times 10^3$$ 57.6 116.7 Table 7 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A)$$, $$\ell=0,1,2,3,4$$ and the CPU time in seconds. The $${\tt EPA}$$ matrix, $$n=4772, r=951$$ $$\ell$$ Method Rel_Est Abs_Est CPU 0 $${\tt funm\_condest1}$$ $$9.90\times 10^{3}$$ $$2.60\times 10^{4}$$ 90.1 $${\tt Algorithm\,2}$$ $$3.34\times 10^4$$ $$1.52\times 10^{5}$$ 95.4 1 $${\tt funm\_condest1}$$ $$1.09\times 10^4$$ $$8.05\times 10^{3}$$ 272.1 $${\tt Algorithm\,2}$$ $$2.25\times 10^4$$ $$2.74\times 10^{4}$$ 100.6 2 $${\tt funm\_condest1}$$ $$8.95\times 10^3$$ $$2.03\times 10^{3}$$ 443.7 $${\tt Algorithm\,2}$$ $$1.47\times 10^4$$ $$4.27\times 10^3$$ 105.5 3 $${\tt funm\_condest1}$$ $$7.86\times 10^3$$ $$425.9$$ 773.0 $${\tt Algorithm\,2}$$ $$9.71\times 10^3$$ $$587.4$$ 110.7 4 $${\tt funm\_condest1}$$ $$7.17\times 10^3$$ 75.2 1357.7 $${\tt Algorithm\,2}$$ $$5.29\times 10^3$$ 57.6 116.7 $$\ell$$ Method Rel_Est Abs_Est CPU 0 $${\tt funm\_condest1}$$ $$9.90\times 10^{3}$$ $$2.60\times 10^{4}$$ 90.1 $${\tt Algorithm\,2}$$ $$3.34\times 10^4$$ $$1.52\times 10^{5}$$ 95.4 1 $${\tt funm\_condest1}$$ $$1.09\times 10^4$$ $$8.05\times 10^{3}$$ 272.1 $${\tt Algorithm\,2}$$ $$2.25\times 10^4$$ $$2.74\times 10^{4}$$ 100.6 2 $${\tt funm\_condest1}$$ $$8.95\times 10^3$$ $$2.03\times 10^{3}$$ 443.7 $${\tt Algorithm\,2}$$ $$1.47\times 10^4$$ $$4.27\times 10^3$$ 105.5 3 $${\tt funm\_condest1}$$ $$7.86\times 10^3$$ $$425.9$$ 773.0 $${\tt Algorithm\,2}$$ $$9.71\times 10^3$$ $$587.4$$ 110.7 4 $${\tt funm\_condest1}$$ $$7.17\times 10^3$$ 75.2 1357.7 $${\tt Algorithm\,2}$$ $$5.29\times 10^3$$ 57.6 116.7 Table 8 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(-A), \ell=0,1,2,3,4$$, and the CPU time in seconds. The $${\tt eris1176}$$ matrix, $$n=1176$$, $$r=774$$ $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$575.1$$ $$1.45\times 10^{3}$$ 2.47 $${\tt Algorithm\,2}$$ $$2.72\times 10^3$$ $$4.66\times 10^{3}$$ 11.1 1 $${\tt funm\_condest1}$$ $$1.09\times 10^3$$ $$548.7$$ 38.2 $${\tt Algorithm\,2}$$ $$2.66\times 10^3$$ $$890.5$$ 14.5 2 $${\tt funm\_condest1}$$ $$1.75\times 10^3$$ $$171.0$$ 60.2 $${\tt Algorithm\,2}$$ $$2.71\times 10^3$$ $$167.6$$ 17.6 3 $${\tt funm\_condest1}$$ $$2.37\times 10^3$$ $$42.2$$ 82.2 $${\tt Algorithm\,2}$$ $$2.82\times 10^3$$ 29.4 20.5 4 $${\tt funm\_condest1}$$ $$2.85\times 10^3$$ 8.49 103.4 $${\tt Algorithm\,2}$$ $$2.88\times 10^3$$ 4.61 24.2 $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$575.1$$ $$1.45\times 10^{3}$$ 2.47 $${\tt Algorithm\,2}$$ $$2.72\times 10^3$$ $$4.66\times 10^{3}$$ 11.1 1 $${\tt funm\_condest1}$$ $$1.09\times 10^3$$ $$548.7$$ 38.2 $${\tt Algorithm\,2}$$ $$2.66\times 10^3$$ $$890.5$$ 14.5 2 $${\tt funm\_condest1}$$ $$1.75\times 10^3$$ $$171.0$$ 60.2 $${\tt Algorithm\,2}$$ $$2.71\times 10^3$$ $$167.6$$ 17.6 3 $${\tt funm\_condest1}$$ $$2.37\times 10^3$$ $$42.2$$ 82.2 $${\tt Algorithm\,2}$$ $$2.82\times 10^3$$ 29.4 20.5 4 $${\tt funm\_condest1}$$ $$2.85\times 10^3$$ 8.49 103.4 $${\tt Algorithm\,2}$$ $$2.88\times 10^3$$ 4.61 24.2 Table 8 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(-A), \ell=0,1,2,3,4$$, and the CPU time in seconds. The $${\tt eris1176}$$ matrix, $$n=1176$$, $$r=774$$ $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$575.1$$ $$1.45\times 10^{3}$$ 2.47 $${\tt Algorithm\,2}$$ $$2.72\times 10^3$$ $$4.66\times 10^{3}$$ 11.1 1 $${\tt funm\_condest1}$$ $$1.09\times 10^3$$ $$548.7$$ 38.2 $${\tt Algorithm\,2}$$ $$2.66\times 10^3$$ $$890.5$$ 14.5 2 $${\tt funm\_condest1}$$ $$1.75\times 10^3$$ $$171.0$$ 60.2 $${\tt Algorithm\,2}$$ $$2.71\times 10^3$$ $$167.6$$ 17.6 3 $${\tt funm\_condest1}$$ $$2.37\times 10^3$$ $$42.2$$ 82.2 $${\tt Algorithm\,2}$$ $$2.82\times 10^3$$ 29.4 20.5 4 $${\tt funm\_condest1}$$ $$2.85\times 10^3$$ 8.49 103.4 $${\tt Algorithm\,2}$$ $$2.88\times 10^3$$ 4.61 24.2 $$\ell$$ Method Relative_Est Absolute_Est CPU 0 $${\tt funm\_condest1}$$ $$575.1$$ $$1.45\times 10^{3}$$ 2.47 $${\tt Algorithm\,2}$$ $$2.72\times 10^3$$ $$4.66\times 10^{3}$$ 11.1 1 $${\tt funm\_condest1}$$ $$1.09\times 10^3$$ $$548.7$$ 38.2 $${\tt Algorithm\,2}$$ $$2.66\times 10^3$$ $$890.5$$ 14.5 2 $${\tt funm\_condest1}$$ $$1.75\times 10^3$$ $$171.0$$ 60.2 $${\tt Algorithm\,2}$$ $$2.71\times 10^3$$ $$167.6$$ 17.6 3 $${\tt funm\_condest1}$$ $$2.37\times 10^3$$ $$42.2$$ 82.2 $${\tt Algorithm\,2}$$ $$2.82\times 10^3$$ 29.4 20.5 4 $${\tt funm\_condest1}$$ $$2.85\times 10^3$$ 8.49 103.4 $${\tt Algorithm\,2}$$ $$2.88\times 10^3$$ 4.61 24.2 Table 9 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$; the $${\tt Abortion}$$ matrix, $$n=2293, r=633$$. In the new algorithm, we compute the condition numbers consecutively, while in $${\tt funm\_condest1}$$ we have to compute them one by one. The CPU timings of Algorithm 2 and $${\tt funm\_condest1}$$ are $${104.4}$$ and $$336.9$$ s, respectively $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$1.37\times 10^3$$ $$9.87\times 10^{4}$$ $${\tt Algorithm\,2}$$ 194.4 $$2.90\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$1.93\times 10^3$$ $$1.77\times 10^{4}$$ $${\tt Algorithm\,2}$$ 170.8 $$3.24\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$2.72\times 10^3$$ $$3.14\times 10^3$$ $${\tt Algorithm\,2}$$ 144.0 343.2 3 $${\tt funm\_condest1}$$ $$3.64\times 10^3$$ 522.9 $${\tt Algorithm\,2}$$ 135.2 39.8 4 $${\tt funm\_condest1}$$ $$4.58\times 10^3$$ 79.1 $${\tt Algorithm\,2}$$ 130.7 4.58 $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$1.37\times 10^3$$ $$9.87\times 10^{4}$$ $${\tt Algorithm\,2}$$ 194.4 $$2.90\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$1.93\times 10^3$$ $$1.77\times 10^{4}$$ $${\tt Algorithm\,2}$$ 170.8 $$3.24\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$2.72\times 10^3$$ $$3.14\times 10^3$$ $${\tt Algorithm\,2}$$ 144.0 343.2 3 $${\tt funm\_condest1}$$ $$3.64\times 10^3$$ 522.9 $${\tt Algorithm\,2}$$ 135.2 39.8 4 $${\tt funm\_condest1}$$ $$4.58\times 10^3$$ 79.1 $${\tt Algorithm\,2}$$ 130.7 4.58 Table 9 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$; the $${\tt Abortion}$$ matrix, $$n=2293, r=633$$. In the new algorithm, we compute the condition numbers consecutively, while in $${\tt funm\_condest1}$$ we have to compute them one by one. The CPU timings of Algorithm 2 and $${\tt funm\_condest1}$$ are $${104.4}$$ and $$336.9$$ s, respectively $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$1.37\times 10^3$$ $$9.87\times 10^{4}$$ $${\tt Algorithm\,2}$$ 194.4 $$2.90\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$1.93\times 10^3$$ $$1.77\times 10^{4}$$ $${\tt Algorithm\,2}$$ 170.8 $$3.24\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$2.72\times 10^3$$ $$3.14\times 10^3$$ $${\tt Algorithm\,2}$$ 144.0 343.2 3 $${\tt funm\_condest1}$$ $$3.64\times 10^3$$ 522.9 $${\tt Algorithm\,2}$$ 135.2 39.8 4 $${\tt funm\_condest1}$$ $$4.58\times 10^3$$ 79.1 $${\tt Algorithm\,2}$$ 130.7 4.58 $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$1.37\times 10^3$$ $$9.87\times 10^{4}$$ $${\tt Algorithm\,2}$$ 194.4 $$2.90\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$1.93\times 10^3$$ $$1.77\times 10^{4}$$ $${\tt Algorithm\,2}$$ 170.8 $$3.24\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$2.72\times 10^3$$ $$3.14\times 10^3$$ $${\tt Algorithm\,2}$$ 144.0 343.2 3 $${\tt funm\_condest1}$$ $$3.64\times 10^3$$ 522.9 $${\tt Algorithm\,2}$$ 135.2 39.8 4 $${\tt funm\_condest1}$$ $$4.58\times 10^3$$ 79.1 $${\tt Algorithm\,2}$$ 130.7 4.58 Table 10 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$; the $${\tt Computational\_Geometry}$$ matrix, $$n=1226, r=415$$. In the new algorithm, we compute the condition numbers $${\tt consecutively}$$, while in $${\tt funm\_condest1}$$ we have to compute them one by one. The CPU timings of Algorithm 2 and $${\tt funm\_condest1}$$ are $${22.8}$$ and $$66.2$$ s, respectively $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$5.35\times 10^3$$ $$2.73\times 10^{4}$$ $${\tt Algorithm\,2}$$ $$3.15\times 10^3$$ $$2.35\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$4.00\times 10^3$$ $$6.33\times 10^{3}$$ $${\tt Algorithm\,2}$$ 998.6 $$2.07\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$3.20\times 10^3$$ $$1.29\times 10^3$$ $${\tt Algorithm\,2}$$ 860.8 420.1 3 $${\tt funm\_condest1}$$ $$2.69\times 10^3$$ 229.5 $${\tt Algorithm\,2}$$ 824.9 80.4 4 $${\tt funm\_condest1}$$ $$2.36\times 10^3$$ 36.0 $${\tt Algorithm\,2}$$ 262.1 4.38 $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$5.35\times 10^3$$ $$2.73\times 10^{4}$$ $${\tt Algorithm\,2}$$ $$3.15\times 10^3$$ $$2.35\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$4.00\times 10^3$$ $$6.33\times 10^{3}$$ $${\tt Algorithm\,2}$$ 998.6 $$2.07\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$3.20\times 10^3$$ $$1.29\times 10^3$$ $${\tt Algorithm\,2}$$ 860.8 420.1 3 $${\tt funm\_condest1}$$ $$2.69\times 10^3$$ 229.5 $${\tt Algorithm\,2}$$ 824.9 80.4 4 $${\tt funm\_condest1}$$ $$2.36\times 10^3$$ 36.0 $${\tt Algorithm\,2}$$ 262.1 4.38 Table 10 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$; the $${\tt Computational\_Geometry}$$ matrix, $$n=1226, r=415$$. In the new algorithm, we compute the condition numbers $${\tt consecutively}$$, while in $${\tt funm\_condest1}$$ we have to compute them one by one. The CPU timings of Algorithm 2 and $${\tt funm\_condest1}$$ are $${22.8}$$ and $$66.2$$ s, respectively $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$5.35\times 10^3$$ $$2.73\times 10^{4}$$ $${\tt Algorithm\,2}$$ $$3.15\times 10^3$$ $$2.35\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$4.00\times 10^3$$ $$6.33\times 10^{3}$$ $${\tt Algorithm\,2}$$ 998.6 $$2.07\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$3.20\times 10^3$$ $$1.29\times 10^3$$ $${\tt Algorithm\,2}$$ 860.8 420.1 3 $${\tt funm\_condest1}$$ $$2.69\times 10^3$$ 229.5 $${\tt Algorithm\,2}$$ 824.9 80.4 4 $${\tt funm\_condest1}$$ $$2.36\times 10^3$$ 36.0 $${\tt Algorithm\,2}$$ 262.1 4.38 $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$5.35\times 10^3$$ $$2.73\times 10^{4}$$ $${\tt Algorithm\,2}$$ $$3.15\times 10^3$$ $$2.35\times 10^{4}$$ 1 $${\tt funm\_condest1}$$ $$4.00\times 10^3$$ $$6.33\times 10^{3}$$ $${\tt Algorithm\,2}$$ 998.6 $$2.07\times 10^{3}$$ 2 $${\tt funm\_condest1}$$ $$3.20\times 10^3$$ $$1.29\times 10^3$$ $${\tt Algorithm\,2}$$ 860.8 420.1 3 $${\tt funm\_condest1}$$ $$2.69\times 10^3$$ 229.5 $${\tt Algorithm\,2}$$ 824.9 80.4 4 $${\tt funm\_condest1}$$ $$2.36\times 10^3$$ 36.0 $${\tt Algorithm\,2}$$ 262.1 4.38 Table 11 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$; the $${\tt Computational\_Complexity}$$ matrix, $$n=884, r=254$$. In the new algorithm, we compute the condition numbers $${\tt consecutively}$$, while in $${\tt funm\_condest1}$$ we have to compute them one by one. The CPU timings of Algorithm 2 and $${\tt funm\_condest1}$$ are $${6.24}$$ and $$27.5$$ s, respectively $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$2.92\times 10^3$$ $$5.27\times 10^{3}$$ $${\tt Algorithm\,2}$$ 873.9 $$2.17\times 10^{3}$$ 1 $${\tt funm\_condest1}$$ $$2.52\times 10^3$$ $$1.60\times 10^{3}$$ $${\tt Algorithm\,2}$$ 530.1 $$438.4$$ 2 $${\tt funm\_condest1}$$ $$2.04\times 10^3$$ 397.0 $${\tt Algorithm\,2}$$ 331.8 73.8 3 $${\tt funm\_condest1}$$ $$1.70\times 10^3$$ 82.0 $${\tt Algorithm\,2}$$ 211.1 10.8 4 $${\tt funm\_condest1}$$ $$1.49\times 10^3$$ 14.4 $${\tt Algorithm\,2}$$ 140.9 1.37 $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$2.92\times 10^3$$ $$5.27\times 10^{3}$$ $${\tt Algorithm\,2}$$ 873.9 $$2.17\times 10^{3}$$ 1 $${\tt funm\_condest1}$$ $$2.52\times 10^3$$ $$1.60\times 10^{3}$$ $${\tt Algorithm\,2}$$ 530.1 $$438.4$$ 2 $${\tt funm\_condest1}$$ $$2.04\times 10^3$$ 397.0 $${\tt Algorithm\,2}$$ 331.8 73.8 3 $${\tt funm\_condest1}$$ $$1.70\times 10^3$$ 82.0 $${\tt Algorithm\,2}$$ 211.1 10.8 4 $${\tt funm\_condest1}$$ $$1.49\times 10^3$$ 14.4 $${\tt Algorithm\,2}$$ 140.9 1.37 Table 11 Example $$3$$: Estimation of the relative and the absolute condition numbers of $$\varphi_{\ell}(A), \ell=0,1,2,3,4$$; the $${\tt Computational\_Complexity}$$ matrix, $$n=884, r=254$$. In the new algorithm, we compute the condition numbers $${\tt consecutively}$$, while in $${\tt funm\_condest1}$$ we have to compute them one by one. The CPU timings of Algorithm 2 and $${\tt funm\_condest1}$$ are $${6.24}$$ and $$27.5$$ s, respectively $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$2.92\times 10^3$$ $$5.27\times 10^{3}$$ $${\tt Algorithm\,2}$$ 873.9 $$2.17\times 10^{3}$$ 1 $${\tt funm\_condest1}$$ $$2.52\times 10^3$$ $$1.60\times 10^{3}$$ $${\tt Algorithm\,2}$$ 530.1 $$438.4$$ 2 $${\tt funm\_condest1}$$ $$2.04\times 10^3$$ 397.0 $${\tt Algorithm\,2}$$ 331.8 73.8 3 $${\tt funm\_condest1}$$ $$1.70\times 10^3$$ 82.0 $${\tt Algorithm\,2}$$ 211.1 10.8 4 $${\tt funm\_condest1}$$ $$1.49\times 10^3$$ 14.4 $${\tt Algorithm\,2}$$ 140.9 1.37 $$\ell$$ Method Relative_Est Absolute_Est 0 $${\tt funm\_condest1}$$ $$2.92\times 10^3$$ $$5.27\times 10^{3}$$ $${\tt Algorithm\,2}$$ 873.9 $$2.17\times 10^{3}$$ 1 $${\tt funm\_condest1}$$ $$2.52\times 10^3$$ $$1.60\times 10^{3}$$ $${\tt Algorithm\,2}$$ 530.1 $$438.4$$ 2 $${\tt funm\_condest1}$$ $$2.04\times 10^3$$ 397.0 $${\tt Algorithm\,2}$$ 331.8 73.8 3 $${\tt funm\_condest1}$$ $$1.70\times 10^3$$ 82.0 $${\tt Algorithm\,2}$$ 211.1 10.8 4 $${\tt funm\_condest1}$$ $$1.49\times 10^3$$ 14.4 $${\tt Algorithm\,2}$$ 140.9 1.37 As was pointed out by Higham (2008a, p. 64), for the absolute and the relative condition numbers, ‘what is needed is an estimate that is of the correct order of magnitude in practice—more than one correct significant digit is not needed’. Recall that the $${\tt funm\_condest1.m}$$ function estimates the 1-norm relative and absolute condition numbers, while Algorithm 2 estimates the 2-norm condition numbers. Compared with the numerical results of $${\tt funm\_condest1}$$, we see from Tables 6–11 that our new strategy captures the correct order of magnitude of the condition numbers in many cases. The difference from $${\tt funm\_condest1}$$ can partially be explained by the change of norm, and the condition numbers presented differ in (at most) up to one digit. On the other hand, Algorithm 2 often runs much faster than $${\tt funm\_condest1}$$, especially when $$\ell\geq 1$$. For instance, for the $${\tt EVA}$$ matrix, the new algorithm used 425.2 s while $${\tt funm\_condest1}$$ used 5846.1 s for $$\ell=1$$. Specifically, for $$\ell\geq 2$$, $${\tt funm\_condest1}$$ fails to converge within 3 h, and Algorithm 2 runs quite well. Moreover, we see from Tables 9 to 11 that Algorithm 2 can compute the condition numbers consecutively, whereas $${\tt funm\_condest1}$$ has to compute them one by one. So we benefit from our new strategy, and we provide a competitive algorithm for estimating the relative and the absolute condition numbers of $$\varphi$$-functions. 4.4 Sharpness of Theorem 2.6 In this example, we aim to show the sharpness of our error analysis for approximating $$\varphi$$-functions. The test matrix is the $${\tt watt\_1}$$ matrix that was used in Example 2. define $$ {\tt Abs\_Err2}=\|\varphi_{\ell}(A)-{\varphi_{\ell}(\widetilde{A})}\|_2, $$ and $$ {\tt Rel\_Err2}=\frac{\|\varphi_{\ell}(A)-{\varphi_{\ell}(\widetilde{A})}\|_2}{\|\varphi_{\ell}(A)\|_2}, $$ the ‘exact’ absolute and the relative errors of the computed solutions $${\varphi_{\ell}(\widetilde{A})}$$ with respect to $$\varphi_{\ell}(A)$$ in terms of the 2-norm. In our estimations, the values of $${\rm cond}_{\rm abs}(\varphi_{\ell},A)$$ and $${\rm cond}_{\rm rel}(\varphi_{\ell},A)$$ in the upper bounds of (2.13) and (2.14) are approximated by using (3.7) and (3.10). The corresponding estimations are denoted by ‘(2.13)–NewStr’ and ‘(2.14)–NewStr’, respectively. Table 12 lists the numerical results, as well as the estimated absolute and relative condition numbers $${\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})$$, $${\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})$$. Table 12 Example $$4$$: Absolute and relative errors and their estimations in terms of the $$2$$-norm, $$\ell=0,1,2,3,4$$. The $${\tt watt\_1}$$ matrix, with $$\|A-\widetilde{A}\|_2\approx 1.07\times 10^{-6}$$. Here eqrefeq2.16-NewStr, (2.14)-NewStr denote that the absolute and the relative condition numbers $${\rm cond}_{\rm abs}(\varphi_{\ell},A)$$ and $${\rm cond}_{\rm rel}(\varphi_{\ell},A)$$ in the upper bounds of (2.13) and (2.14) are estimated by using (3.7) and (3.10), respectively $$\ell=0$$ $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $${\tt Abs_Err2}$$ $$1.07\times 10^{-6}$$ $$5.35\times 10^{-7}$$ $$1.78\times 10^{-7}$$ $$4.45\times 10^{-8}$$ $$8.91\times 10^{-9}$$ $${\tt Rel_Err2}$$ $$3.93\times 10^{-7}$$ $$3.11\times 10^{-7}$$ $$2.48\times 10^{-7}$$ $$2.04\times 10^{-7}$$ $$1.73\times 10^{-7}$$ (2.13)–NewStr $$1.83\times 10^{-7}$$ $$5.24\times 10^{-8}$$ $$1.19\times 10^{-8}$$ $$2.22\times 10^{-9}$$ $$3.47\times 10^{-10}$$ (2.14)–NewStr $$1.06\times 10^{-7}$$ $$7.30\times 10^{-8}$$ $$5.43\times 10^{-8}$$ $$4.30\times 10^{-8}$$ $$3.49\times 10^{-8}$$ $${\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ 0.17 $$6.73\times 10^{-2}$$ $$5.01\times 10^{-2}$$ $$3.92\times 10^{-2}$$ $$3.23\times 10^{-2}$$ $${\rm cond}_{\rm rel}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ $$9.81\times 10^{-2}$$ $$4.83\times 10^{-2}$$ $$1.09\times 10^{-2}$$ $$2.00\times 10^{-3}$$ $$3.21\times 10^{-4}$$ $$\ell=0$$ $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $${\tt Abs_Err2}$$ $$1.07\times 10^{-6}$$ $$5.35\times 10^{-7}$$ $$1.78\times 10^{-7}$$ $$4.45\times 10^{-8}$$ $$8.91\times 10^{-9}$$ $${\tt Rel_Err2}$$ $$3.93\times 10^{-7}$$ $$3.11\times 10^{-7}$$ $$2.48\times 10^{-7}$$ $$2.04\times 10^{-7}$$ $$1.73\times 10^{-7}$$ (2.13)–NewStr $$1.83\times 10^{-7}$$ $$5.24\times 10^{-8}$$ $$1.19\times 10^{-8}$$ $$2.22\times 10^{-9}$$ $$3.47\times 10^{-10}$$ (2.14)–NewStr $$1.06\times 10^{-7}$$ $$7.30\times 10^{-8}$$ $$5.43\times 10^{-8}$$ $$4.30\times 10^{-8}$$ $$3.49\times 10^{-8}$$ $${\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ 0.17 $$6.73\times 10^{-2}$$ $$5.01\times 10^{-2}$$ $$3.92\times 10^{-2}$$ $$3.23\times 10^{-2}$$ $${\rm cond}_{\rm rel}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ $$9.81\times 10^{-2}$$ $$4.83\times 10^{-2}$$ $$1.09\times 10^{-2}$$ $$2.00\times 10^{-3}$$ $$3.21\times 10^{-4}$$ Table 12 Example $$4$$: Absolute and relative errors and their estimations in terms of the $$2$$-norm, $$\ell=0,1,2,3,4$$. The $${\tt watt\_1}$$ matrix, with $$\|A-\widetilde{A}\|_2\approx 1.07\times 10^{-6}$$. Here eqrefeq2.16-NewStr, (2.14)-NewStr denote that the absolute and the relative condition numbers $${\rm cond}_{\rm abs}(\varphi_{\ell},A)$$ and $${\rm cond}_{\rm rel}(\varphi_{\ell},A)$$ in the upper bounds of (2.13) and (2.14) are estimated by using (3.7) and (3.10), respectively $$\ell=0$$ $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $${\tt Abs_Err2}$$ $$1.07\times 10^{-6}$$ $$5.35\times 10^{-7}$$ $$1.78\times 10^{-7}$$ $$4.45\times 10^{-8}$$ $$8.91\times 10^{-9}$$ $${\tt Rel_Err2}$$ $$3.93\times 10^{-7}$$ $$3.11\times 10^{-7}$$ $$2.48\times 10^{-7}$$ $$2.04\times 10^{-7}$$ $$1.73\times 10^{-7}$$ (2.13)–NewStr $$1.83\times 10^{-7}$$ $$5.24\times 10^{-8}$$ $$1.19\times 10^{-8}$$ $$2.22\times 10^{-9}$$ $$3.47\times 10^{-10}$$ (2.14)–NewStr $$1.06\times 10^{-7}$$ $$7.30\times 10^{-8}$$ $$5.43\times 10^{-8}$$ $$4.30\times 10^{-8}$$ $$3.49\times 10^{-8}$$ $${\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ 0.17 $$6.73\times 10^{-2}$$ $$5.01\times 10^{-2}$$ $$3.92\times 10^{-2}$$ $$3.23\times 10^{-2}$$ $${\rm cond}_{\rm rel}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ $$9.81\times 10^{-2}$$ $$4.83\times 10^{-2}$$ $$1.09\times 10^{-2}$$ $$2.00\times 10^{-3}$$ $$3.21\times 10^{-4}$$ $$\ell=0$$ $$\ell=1$$ $$\ell=2$$ $$\ell=3$$ $$\ell=4$$ $${\tt Abs_Err2}$$ $$1.07\times 10^{-6}$$ $$5.35\times 10^{-7}$$ $$1.78\times 10^{-7}$$ $$4.45\times 10^{-8}$$ $$8.91\times 10^{-9}$$ $${\tt Rel_Err2}$$ $$3.93\times 10^{-7}$$ $$3.11\times 10^{-7}$$ $$2.48\times 10^{-7}$$ $$2.04\times 10^{-7}$$ $$1.73\times 10^{-7}$$ (2.13)–NewStr $$1.83\times 10^{-7}$$ $$5.24\times 10^{-8}$$ $$1.19\times 10^{-8}$$ $$2.22\times 10^{-9}$$ $$3.47\times 10^{-10}$$ (2.14)–NewStr $$1.06\times 10^{-7}$$ $$7.30\times 10^{-8}$$ $$5.43\times 10^{-8}$$ $$4.30\times 10^{-8}$$ $$3.49\times 10^{-8}$$ $${\rm cond}_{\rm abs}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ 0.17 $$6.73\times 10^{-2}$$ $$5.01\times 10^{-2}$$ $$3.92\times 10^{-2}$$ $$3.23\times 10^{-2}$$ $${\rm cond}_{\rm rel}^{\rm est}(\varphi_{\ell},\widetilde{A})$$ $$9.81\times 10^{-2}$$ $$4.83\times 10^{-2}$$ $$1.09\times 10^{-2}$$ $$2.00\times 10^{-3}$$ $$3.21\times 10^{-4}$$ It is seen from Table 12 that both (2.13) and (2.14) are very sharp, which justifies our new strategy for estimating $${\rm cond}_{\rm abs}(\varphi_{\ell},A)$$ and $${\rm cond}_{\rm rel}(\varphi_{\ell},A)$$. However, we find that the values of (2.13)-NewStr and (2.14)-NewStr may be a little smaller than those of $${\tt Abs\_Err2}$$ and $${\tt Rel\_Err2}$$ in some cases. In fact, the new strategy gives only some approximations to the absolute and the relative condition numbers, which are neither upper bounds nor lower bounds theoretically. 5. Concluding remarks In this article, we consider the computations, error analysis, implementations and applications of $$\varphi$$-functions for large sparse matrices with low rank or with fast decaying singular values. Given an SCR of $$A$$, we take into account how to compute the approximations to matrix functions $$\varphi_{\ell}(A)$$ efficiently for $$\ell=0,1,2,\ldots,p$$ and to estimate their 2-norm Fréchet relative and absolute condition numbers effectively. The key is to reduce the computation of $$\varphi_{\ell}$$-functions of a large matrix to that of $$\varphi_{\ell+1}$$-functions of some smaller matrices of size $$r$$-by-$$r$$, where $$r$$ is the numerical rank of the large matrix in question. Moreover, the new algorithms need to store only two $$n$$-by-$$r$$ sparse matrices and some $$r$$-by-$$r$$ matrices, rather than some $$n$$-by-$$n$$ possibly dense matrices. Our theoretical analysis and numerical experiments indicate that the error of the computed $$\varphi$$-function is of the order of the product of the error from SCR approximation with the condition number. Numerical experiments demonstrate the numerical performance of our new algorithms for $$\varphi$$-functions and for the condition numbers and show the effectiveness of our theoretical results. On the other hand, the efficiency of our new approach is closely related to that of reduced-rank approximation of large sparse matrices. Thus, a promising research area is to investigate new technologies to improve the performance of the SCR algorithm on very large matrices. Another interesting topic is to combine other advanced algorithms such as the randomized singular value decomposition algorithm by Halko et al. (2011) and Mahoney (2011) with our new strategies for the approximation of functions of large sparse matrices. Further, our new strategy also applies to other functions that can be expanded into power series. These are interesting topics and deserve further investigation. Acknowledgements We wish to thank the associate editor and the anonymous reviewers for their insightful comments and constructive suggestions that greatly improved the presentation of this article. Specifically, we thank a reviewer for the concise proof of Theorem 2.4. Meanwhile, we are grateful to Juan-juan Tian for helpful discussions. Funding National Science Foundation of China (11371176); Natural Science Foundation of Jiangsu Province; Talent Introduction Program of China University of Mining and Technology to G. Wu. Footnotes 1ftp://ftp.cs.umd.edu/pub/stewart/reports/Contents.html. 2http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html. 3http://vision.ucsd.edu/datasets/yale_face_dataset_original/yalefaces.zip. 4http://vision.ucsd.edu/~leekc/ExtYaleDatabase/Yale%20Face%20Database.html. References Ahmed N. ( 2013 ) Exponential discriminant regularization using nonnegative constraint and image descriptor. IEEE 9th International Conference on Emerging Technologies , pp. 1 – 6 . Al-Mohy A. & Higham N. J. ( 2011 ) Compution the action of the matrix exponential, with an application to exponential integrators. SIAM J. Sci. Comput ., 33 , 488 – 511 . Google Scholar CrossRef Search ADS Benzi M. , Estrada E. & Klymko C. ( 2013 ) Ranking hubs and authorities using matrix functions. Linear Algebra Appl. , 438 , 2447 – 2474 . Google Scholar CrossRef Search ADS Berland H. , Skaflestad B. & Wright W. ( 2007 ) EXPINT–A MATLAB package for exponential integrators. ACM Trans. Math. Software , 33 , Article 4. Berry M. & Browne M. ( 1999 ) Understanding Search Engines: Mathematical Modeling and Text Retrieval . Philadelphia : SIAM. Berry M. , Drmač, Z. & Jessup E. ( 1999 ) Matrices, vector spaces, and information retrieval. SIAM Rev. , 41 , 335 – 362 . Google Scholar CrossRef Search ADS Berry M. , Dumais S. & O’Brien G. ( 1995 ) Using linear algebra for intelligent information retrieval. SIAM Rev. , 37 , 573 – 595 . Google Scholar CrossRef Search ADS Berry M. , Pulatova S. & Stewart G. W. ( 2005 ) Computing sparse reduced-rank approximations to sparse matrices. ACM Trans. Math. Software , 31 , 252 – 269 . Google Scholar CrossRef Search ADS Beylkin G. , Keiser J. & Vozovoi L. ( 1998 ) A new class of time discretization schemes for the solution of nonlinear PDEs. J. Comput. Phys. , 147 , 362 – 387 . Google Scholar CrossRef Search ADS Celledon E. & Iserles A. ( 2000 ) Approximating the exponential from a Lie algebra to a Lie group. Math. Comp. , 69 , 1457 – 1480 . Google Scholar CrossRef Search ADS Davies P. , Higham N. J. & Tisseur F. ( 2001 ) Analysis of the Cholesky method with iterative refinement for solving the symmetric definite generalized eigenproblem. SIAM J. Matrix Anal. Appl. , 23 , 472 – 493 . Google Scholar CrossRef Search ADS Davis T. & Hu Y. ( 2011 ) The University of Florida sparse matrix collection. ACM Trans. Math. Software , 38 , Article 1 . Available at http://www.cise.ufl.edu/research/sparse/matrices/. Dornaika F. & Bosaghzadeh A. ( 2013 ) Exponential local discriminant embedding and its application to face recognition. IEEE Trans. Syst. Man Cybern. B Cybern. , 43 , 921 – 934 . Duda R. , Hart P. & Stork D. ( 2000 ) Pattern Classification , 2nd edn. New York : Wiley. Estrada E. & Hatano N. ( 2008 ) Communicability in complex networks Phys. Rev. E , 77 , 036111 . Google Scholar CrossRef Search ADS Estrada E. & Higham D. J. ( 2010 ) Network properties revealed through matrix functions. SIAM Rev. , 52 , 671 – 696 . Google Scholar CrossRef Search ADS Estrada E. & Rodríguez-Velázquez J. ( 2005 ) Subgraph centrality in complex networks. Phys. Rev. E , 71 , 056103. Google Scholar CrossRef Search ADS Golub G. H. & Van Loan C. F. ( 2013 ) Matrix Computations , 4th edn. Baltimore, MD : Johns Hopkins University Press. Halko N. , Martinsson P. & Tropp J. ( 2011 ) Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Rev. , 53 , 217 – 288 . Google Scholar CrossRef Search ADS Hansen P. ( 1998 ) Rank-Deficient and Discrete Ill-Posed Problems. Philadelphia : SIAM . Google Scholar CrossRef Search ADS Higham N. J. ( 2002 ) Accuracy and Stability of Numerical Algorithms , 2nd edn . Philadelphia : SIAM. Google Scholar CrossRef Search ADS Higham N. J. ( 2008a ) Functions of Matrices: Theory and Computation . Philadelphia : SIAM. Google Scholar CrossRef Search ADS Higham N. J. ( 2008b ) The Matrix Function Toolbox . Available at http://www.maths.manchester.ac.uk/~higham/ mftoolbox. Higham N. J. ( 2009 ) The scaling and squaring method for the matrix exponential revisited. SIAM Rev. , 51 , 747 – 764 . Google Scholar CrossRef Search ADS Higham N. J. , Mackey D. & Tisseur F. ( 2009 ) Definite matrix polynomials and their linearization by definite pencils. SIAM J. Matrix Anal. Appl. , 31 , 478 – 502 . Google Scholar CrossRef Search ADS Hochbruck M. , Lubich C. & Selhofer H. ( 1998 ) Exponential integrators for large systems of differential equations. SIAM J. Sci. Comput. , 19 , 1552 – 1574 . Google Scholar CrossRef Search ADS Hochbruck M. & Ostermann A. ( 2010 ) Exponential integrators. Acta Numer. , 209 – 286 . Hofmannn B. ( 1986 ) Regularization for Applied Inverse and Ill-posed Problems . Stuttgart, German : Teubner. Google Scholar CrossRef Search ADS Jiang P. & Berry M. ( 2000 ) Solving total least squares problems in information retrieval. Linear Algebra Appl. , 316 , 137 – 156 . Google Scholar CrossRef Search ADS Kassam A. & Trefethen L. N. ( 2005 ) Fourth-order time-stepping for stiff PDEs. SIAM J. Sci. Comput. , 26 , 1214 – 1233 . Google Scholar CrossRef Search ADS Lu Y. ( 2003 ) Computing a matrix function for exponential integrators. J. Comput. Appl. Math. , 161 , 203 – 216 . Google Scholar CrossRef Search ADS Mahoney M. ( 2011 ) Randomized algorithms for matrices and data . Foundations and Trends in Machine Learning, vol. 3, Boston : Now publishers inc., pp. 123 – 224 . Moler C. & Van Loan C. F. ( 2003 ) Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later. SIAM Rev. , 45 , 3 – 49 . Google Scholar CrossRef Search ADS Niesen J. & Wright W. ( 2012 ) Algorithm 919: a Krylov subspace algorithm for evaluating the $$\varphi$$-functions appearing in exponential integrators. ACM Trans. Math. Software , 38 , Article 22. Park C. & Park H. ( 2008 ) A comparision of generalized linear discriminant analysis algorithms. Pattern Recogn. , 41 , 1083 – 1097 . Google Scholar CrossRef Search ADS Rice J. ( 1966 ) A theory of condition. SIAM J. Numer. Anal. , 3 , 287 – 310 . Google Scholar CrossRef Search ADS Schmelzer T. & Trefethen L. N. ( 2007 ) Evaluating matrix functions for exponential integrators via Carathéodory-Fejér approximation and contour integrals. Electron. Trans. Numer. Anal. , 29 , 1 – 18 . Sidje R. ( 1998 ) EXPOKIT: Software package for computing matrix exponentials. ACM Trans. Math. Software , 24 , 130 – 156 . Google Scholar CrossRef Search ADS Skaflestad B. & Wright W. ( 2009 ) The scaling and squaring method for matrix functions related to the exponential. Appl. Numer. Math. , 59 , 783 – 799 . Google Scholar CrossRef Search ADS Stewart G. W. ( 1999 ) Four algorithms for the efficient computation of truncated pivoted QR approximations to a sparse matrix. Numer. Math. , 83 , 313 – 323 . Google Scholar CrossRef Search ADS Stewart G. W. ( 2005 ) Error analysis of the quasi-Gram–Schmidt algorithm. SIAM J. Matrix Anal. Appl ., 27 , 493 – 506 . Google Scholar CrossRef Search ADS Stuart G. & Berry M. ( 2003 ) A comprehensive whole genome bacterial phylogeny using correlated peptide motifs defined in a high dimensional vector space. J. Bioinform. Comput. Bio. , 1 , 475 – 493 . Google Scholar CrossRef Search ADS Tspras P. Datasets for Experiments on Link Analysis Ranking Algorithms, Available at http://www.cs.toronto.edu/~tsap/experiments/datasets/index.html. Wang S. , Chen H. , Peng X. & Zhou C. ( 2011 ) Exponential locality preserving projections for small sample size problem. Neurocomputing , 74 , 36 – 54 . Wang S. , Yan S. , Yang J. , Zhou C. & Fu X. ( 2014 ) A general exponential framework for dimensionality reduction. IEEE Trans. Image Process. , 23 , 920 – 930 . Google Scholar CrossRef Search ADS PubMed Wu G. & Feng T. ( 2015 ) A theoretical contribution to the fast implementation of null linear discriminant analysis with random matrix multiplication. Numer. Linear Algebra Appl. , 22 , 1180 – 1188 . Google Scholar CrossRef Search ADS Wu G. , Feng T. , Zhang L. & Yang M. ( 2017 ) Inexact implementation using Krylov subspace methods for large scale exponential discriminant analysis with applications to high dimensionality reduction problems. Pattern Recogn. , 66 , 328 – 341 . Google Scholar CrossRef Search ADS Wu G. , Zhang L. & Xu T. ( 2016 ) A framework of the harmonic Arnoldi method for evaluating $$\varphi$$-functions with applications to exponential integrators. Adv. Comput. Math. , 42 , 505 – 541 . Google Scholar CrossRef Search ADS Yan L. & Pan J. ( 2010 ) Two-dimensional exponential discriminant analysis and its application to face recognition. International Conference on Computational Aspects of Social Networks (CASoN) , pp. 528 – 531 . Zhang T. , Fang B. , Tang Y. , Shang Z. & Xu B. ( 2010 ) Generalized discriminant analysis: a matrix exponential approach. IEEE Trans. Syst. Man Cybern. B Cybern. , 40 , 186 – 197 . Google Scholar CrossRef Search ADS PubMed Zhang Z. , Zha H. & Simon H. ( 2002 ) Low-rank approximations with sparse factors I: basic algorithms and error analysis. SIAM J. Matrix Anal. Appl. , 23 , 706 – 727 . Google Scholar CrossRef Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Aug 13, 2017

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off