Search

This work deals with the numerical solution of systems of oscillatory second-order differential equations which often arise from the semi-discretization in space of partial differential equations. Since these differential equations exhibit (pronounced or highly) oscillatory behavior, standard numerical methods are known to perform poorly. Our approach consists in directly discretizing the problem by means of Gautschi-type integrators based on sinc matrix functions. The novelty contained here is that of using a suitable rational approximation formula for the sinc matrix function to apply a rational Krylov-like approximation method with suitable choices of poles. In particular, we discuss the application of the whole strategy to a finite element discretization of the wave equation. [ABSTRACT FROM AUTHOR]

ABSTRACT Among randomized numerical linear algebra strategies, so‐called sketching procedures are emerging as effective reduction means to accelerate the computation of Krylov subspace methods for, for example, the solution of linear systems, eigenvalue computations, and the approximation of matrix functions. While there is plenty of experimental evidence showing that sketched Krylov solvers may dramatically improve performance over standard Krylov methods, especially when combined with a truncated orthogonalization step, many features of these schemes are still unexplored. We derive a new sketched Arnoldi‐type relation that allows us to obtain several different new theoretical results. These lead to an improvement of our understanding of sketched Krylov methods, in particular by explaining why the frequently occurring sketched Ritz values far outside the spectral region of A$$ A $$ do not negatively influence the convergence of sketched Krylov methods for f(A)b$$ f(A)b $$. Our findings also help to identify, among several possible equivalent formulations, the most suitable sketched approximations according to their numerical stability properties. These results are also employed to analyze the error of sketched Krylov methods in the approximation of the action of matrix functions, significantly contributing to the theory available in the current literature. [ABSTRACT FROM AUTHOR]

We obtain an expression for the error in the approximation of f(A)b$$ f(A)\boldsymbol{b} $$ and bTf(A)b$$ {\boldsymbol{b}}^Tf(A)\boldsymbol{b} $$ with rational Krylov methods, where A$$ A $$ is a symmetric matrix, b$$ \boldsymbol{b} $$ is a vector and the function f$$ f $$ admits an integral representation. The error expression is obtained by linking the matrix function error with the error in the approximate solution of shifted linear systems using the same rational Krylov subspace, and it can be exploited to derive both a priori and a posteriori error bounds. The error bounds are a generalization of the ones given in Chen et al. for the Lanczos method for matrix functions. A technique that we employ in the rational Krylov context can also be applied to refine the bounds for the Lanczos case. [ABSTRACT FROM AUTHOR]

This article considers the computation of the matrix exponential eA{{\rm{e}}}^{A} with numerical quadrature. Although several quadrature-based algorithms have been proposed, they focus on (near) Hermitian matrices. In order to deal with non-Hermitian matrices, we use another integral representation including an oscillatory term and consider applying the double exponential (DE) formula specialized to Fourier integrals. The DE formula transforms the given integral into another integral whose interval is infinite, and therefore, it is necessary to truncate the infinite interval. In this article, to utilize the DE formula, we analyze the truncation error and propose two algorithms. The first one approximates eA{{\rm{e}}}^{A} with the fixed mesh size, which is a parameter in the DE formula affecting the accuracy. The second one computes eA{{\rm{e}}}^{A} based on the first one with automatic selection of the mesh size depending on the given error tolerance.

If the matrix function f(At) posses the properties of f(At)=gf(tkA, then the recurrence formula fi−1=gfi,i=N,N−1,⋯,1,f(tA)=f0, can be established. Here, fN=f(AN)=∑j=0majANj,AN=tkNA. This provides an algorithm for computing the matrix function f(At). By specifying the calculation accuracy p, a method is presented to determine m and N in a way that minimizes the time of the above algorithm, thus providing a fast algorithm for f(At). It is important to note that m only depends on the calculation accuracy p and is independent of the matrix A and t. Therefore, f(AN) has a fixed calculation format that is easily computed. On the other hand, N depends not only on A, but also on t. This provides a means to select t such that N is equal to 0, a property of significance. In summary, the algorithm proposed in this article enables users to establish a desired level of accuracy and then utilize it to select the appropriate values for m and N to minimize computation time. This approach ensures that both accuracy and efficiency are addressed concurrently. We develop a general algorithm, then apply it to the exponential, trigonometric, and logarithmic matrix functions, and compare the performance with that of the internal system functions of Mathematica and Pade approximation. In the last section, an example is provided to illustrate the rapid computation of numerical solutions for linear differential equations.

If the matrix function f (A t) posses the properties of f (A t) = g f (t k A , then the recurrence formula f i − 1 = g f i , i = N , N − 1 , ⋯ , 1 , f (t A) = f 0 , can be established. Here, f N = f (A N) = ∑ j = 0 m a j A N j , A N = t k N A . This provides an algorithm for computing the matrix function f (A t) . By specifying the calculation accuracy p, a method is presented to determine m and N in a way that minimizes the time of the above algorithm, thus providing a fast algorithm for f (A t) . It is important to note that m only depends on the calculation accuracy p and is independent of the matrix A and t. Therefore, f (A N) has a fixed calculation format that is easily computed. On the other hand, N depends not only on A, but also on t. This provides a means to select t such that N is equal to 0, a property of significance. In summary, the algorithm proposed in this article enables users to establish a desired level of accuracy and then utilize it to select the appropriate values for m and N to minimize computation time. This approach ensures that both accuracy and efficiency are addressed concurrently. We develop a general algorithm, then apply it to the exponential, trigonometric, and logarithmic matrix functions, and compare the performance with that of the internal system functions of Mathematica and Pade approximation. In the last section, an example is provided to illustrate the rapid computation of numerical solutions for linear differential equations. [ABSTRACT FROM AUTHOR]

The area of matrix functions has received growing interest for a long period of time due to their growing applications. Having a numerical algorithm for a matrix function, the ideal situation is to have an estimate or bound for the error returned alongside the solution. Condition numbers serve this purpose; they measure the first-order sensitivity of matrix functions to perturbations in the input data. We have observed that the existing unstructured condition number leads most of the time to inferior bounds of relative forward errors for some matrix functions at triangular and quasi-triangular matrices. We propose a condition number of matrix functions exploiting the structure of triangular and quasi-triangular matrices. We then adapt an existing algorithm for exact computation of the unstructured condition number to an algorithm for exact evaluation of the structured condition number. Although these algorithms are direct rather than iterative and useful for testing the numerical stability of numerical algorithms, they are less practical for relatively large problems. Therefore, we use an implicit power method approach to estimate the structured condition number. Our numerical experiments show that the structured condition number captures the behavior of the numerical algorithms and provides sharp bounds for the relative forward errors. In addition, the experiment indicates that the power method algorithm is reliable to estimate the structured condition number. [ABSTRACT FROM AUTHOR]

The traditional linear discriminant analysis (LDA) is a classical dimensionality reduction method. But there are two problems with LDA. One is the small-sample-size (SSS) problem, and the other is the classification ability of LDA is not very strong. In this paper, by studying several common methods of LDA, an implicit rule of the criterion of LDA is found, and then a general framework of LDA combined with matrix function is presented. Based on this framework, the polynomials are used to reconstruct the LDA criterion, and then a new method called polynomial linear discriminant analysis (PLDA) is proposed. The proposed PLDA method has two contributions: first, it addresses the small-sample-size problem of LDA, and second, it enhances the distance of the between-class sample through polynomial function mapping, improving the classification ability. Experimental results on ORL, FERET, AR, and Coil datasets show the superiority of PLDA over existing variations of LDA. [ABSTRACT FROM AUTHOR]

We extend the notion of nonbacktracking walks from unweighted graphs to graphs whose edges have a nonnegative weight. Here the weight associated with a walk is taken to be the product over the weights along the individual edges. We give two ways to compute the associated generating function, and corresponding node centrality measures. One method works directly on the original graph and one uses a line graph construction followed by a projection. The first method is more efficient, but the second has the advantage of extending naturally to time-evolving graphs. Based on these generating functions, we define and study corresponding centrality measures. Illustrative computational results are also provided. [ABSTRACT FROM AUTHOR]

We present a theory for general partial derivatives of matrix functions of the form f(A(x)), where A(x) is a matrix path of several variables (x = (x1, . . . ,xj)). Building on results by Mathias [SIAM J. Matrix Anal. Appl., 17 (1996), pp. 610--620] for the first order derivative, we develop a block upper triangular form for higher order partial derivatives. This block form is used to derive conditions for existence and a generalized Daleckiī--Kreī formula for higher order derivatives. We show that certain specializations of this formula lead to classical formulas of quantum perturbation theory. We show how our results are related to earlier results for higher order Frechet derivatives. Block forms of complex step approximations are introduced, and we show how those are related to evaluation of derivatives through the upper triangular form. These relations are illustrated with numerical examples. [ABSTRACT FROM AUTHOR]

The paper is concerned with efficient numerical methods for solving a linear system ϕ (A) x = b , where ϕ (z) is a ϕ -function and A ∈ R N × N. In particular in this work we are interested in the computation of ϕ (A) − 1 b for the case where ϕ (z) = ϕ 1 (z) = e z − 1 z and ϕ (z) = ϕ 2 (z) = e z − 1 − z z 2 . Under suitable conditions on the spectrum of A we design fast algorithms for computing both ϕ ℓ (A) − 1 and ϕ ℓ (A) − 1 b based on Newton's iteration and Krylov-type methods, respectively. Adaptations of these schemes for structured matrices are considered. In particular the cases of banded and more generally quasiseparable matrices are investigated. Numerical results are presented to show the effectiveness of our proposed algorithms. [ABSTRACT FROM AUTHOR]

The computation of the sign function of a matrix plays a crucial role in various mathematical applications. It provides a matrix-valued mapping that determines the sign of each eigenvalue of a nonsingular matrix. In this article, we present a novel iterative algorithm designed to efficiently calculate the sign of an invertible matrix, emphasizing the enlargement of attraction basins. The proposed solver exhibits convergence of order four, making it highly efficient for a wide range of matrices. Furthermore, the method demonstrates global convergence properties. We validate the theoretical outcomes through numerical experiments, which confirm the effectiveness and efficiency of our proposed algorithm. [ABSTRACT FROM AUTHOR]

An algorithm based on the Paterson-Stockmeyer (PS) scheme and filtering technology is developed to compute the large matrix functions in dynamic problems accurately and efficiently. With the assistance of analysis on truncation error and error caused by filtering, the error growth law during the computation is studied, based on which an adaptive filtering threshold is proposed to help the proposed algorithm more efficiently achieve similar accuracy as the original PS scheme. Numerical examples, including 30 random matrices with different bandwidths, 10 adjacency matrices in the complex network dynamic problems, and a trampoline vibration problem, are given to verify the efficiency and accuracy of the proposed algorithm. Numerical results suggest that the proposed algorithm can achieve good accuracy and efficiency in computing the matrix function in the considered dynamic problems. [ABSTRACT FROM AUTHOR]

The Poincaré problem is solved for a system of second-order Fuchs-class differential equations with three and four singular points. A method is proposed for finding the matrices of the monodromy group of an equation in terms of elements of its differential substitutions. The theory of functions from matrices is applied. [ABSTRACT FROM AUTHOR]

The computation of f(A)b, the action of a matrix function on a vector, is a task arising in many areas of scientific computing. In many applications, the matrix A is sparse but so large that only a rather small number of Krylov basis vectors can be stored. Here we discuss a new approach to overcome this limitation by randomized sketching combined with an integral representation of f(A)b. Two different approximation methods are introduced, one based on sketched FOM and another based on sketched GMRES. The convergence of the latter method is analyzed for Stieltjes functions of positive real matrices. We also derive a closed-form expression for the sketched FOM approximant and bound its distance to the full FOM approximant. Numerical experiments demonstrate the potential of the presented sketching approaches. [ABSTRACT FROM AUTHOR]

Any defective matrix M in C n × n has a diagonalizable dilation M , but it is not guaranteed that M preserves the eigenvalues or is a power dilation of M. We construct a diagonalizable dilation with a predictable eigenstructure that retains the eigenvalues of M. Although the construction is not a direct power dilation of M , we show that the powers of M can still be easily computed using the powers of M , and furthermore, this method can be used to compute matrix functions. [ABSTRACT FROM AUTHOR]

The main aim of this paper is to introduce the definitions of generalized order and generalized type of the entire function of several complex matrix variables in hyperspherical region and then study some of their properties. By considering the concepts of generalized order and generalized type, we will extend some results of Abul-Ez et al. [1]. [ABSTRACT FROM AUTHOR]

This paper considers the computation of approximations of matrix functionals of form F (A) : = v T f (A) v , where A is a large symmetric positive definite matrix, v is a vector, and f is a Stieltjes function. The functional F (A) is approximated by a rational Gauss quadrature rule with poles on the negative real axis (or part thereof) in the complex plane, and we focus on the allocation of the poles. Specifically, we propose that the poles, when considered positive point charges, be allocated to make the negative real axis (or part thereof) approximate an equipotential curve. This is easily achieved with the aid of conformal mapping. [ABSTRACT FROM AUTHOR]

The target of this article is to propose a novel numerical scheme that can compete its corresponding solvers of the same type for computing the sign matrix. To do this, a three‐step root solver is designed to gain as much as possible of global convergence when finding the matrix sign. It is proved that the novel numerical scheme is of quartical convergence order. The stability of the scheme is discussed in detail as well. Finally, numerical tests are provided to uphold the theoretical findings. [ABSTRACT FROM AUTHOR]

We propose an integral representation for the higher-order Fréchet derivative of analytic matrix functions f (A) which unifies known results for the first-order Fréchet derivative of general analytic matrix functions and for higher-order Fréchet derivatives of A − 1. We highlight two applications of this integral representation: On the one hand, it allows to find the exact value of the level-2 condition number (i.e., the condition number of the condition number) of f (A) for a large class of functions f when A is Hermitian. On the other hand, it also allows to use numerical quadrature methods to approximate higher-order Fréchet derivatives. We demonstrate that in certain situations—in particular when the derivative order k is moderate and the direction terms in the derivative have low-rank structure—the resulting algorithm can outperform established methods from the literature by a large margin. [ABSTRACT FROM AUTHOR]

Marginal fisher analysis (MFA) is a dimensionality reduction method based on a graph embedding framework. In contrast to traditional linear discriminant analysis (LDA), which requires the data to follow a Gaussian distribution, MFA is suitable for non-Gaussian data, and it has better pattern classification ability. However, MFA has the small-sample-size (SSS) problem. This paper aims to solve the small-sample-size problem while increasing the classification performance of MFA. Based on a matrix function dimensionality reduction framework, the criterion of the MFA method is reconstructed by using the polynomials matrix function transformation, and then a new MFA method is proposed, named PMFA (polynomial marginal fisher analysis). The major contributions of the proposed PMFA method are that it solves the small-sample-size problem of MFA, and it can enlarge the distance between marginal sample points of inter-class, so that it can get better pattern classification performance. Experiments on the public face datasets show that PMFA can get a better classification ability than MFA and its improved methods.

We study walk-based centrality measures for time-ordered network sequences. For the case of standard dynamic walk-counting, we show how to derive and compute centrality measures induced by analytic functions. We also prove that dynamic Katz centrality, based on the resolvent function, has the unique advantage of allowing computations to be performed entirely at the node level. We then consider two distinct types of backtracking and develop a framework for capturing dynamic walk combinatorics when either or both is disallowed. [ABSTRACT FROM AUTHOR]

Hölder matrix functions of the second order are considered. We assume that one element is arbitrary, diagonal elements do not vanish on the contour, and the choice of the last element determines the possibility of their effective factorization. [ABSTRACT FROM AUTHOR]

The derivative of a nn-order matrix function on the complex field is usually defined as a n2{n}^{2}-order matrix, which is not suitable for generalizing Cauchy integral formula of matrix functions to its higher-order derivative form. In this paper, a new kind of derivative of matrix functions is defined, and the higher-order derivative form of Cauchy integral formula of matrix functions is also proved to be true under the new kind of definition of derivative. At the same time, some examples about calculating the values of matrix functions by using Cauchy integral formula of matrix functions and its higher-order derivative form are given.

The computation of the sign function of a matrix plays a crucial role in various mathematical applications. It provides a matrix-valued mapping that determines the sign of each eigenvalue of a nonsingular matrix. In this article, we present a novel iterative algorithm designed to efficiently calculate the sign of an invertible matrix, emphasizing the enlargement of attraction basins. The proposed solver exhibits convergence of order four, making it highly efficient for a wide range of matrices. Furthermore, the method demonstrates global convergence properties. We validate the theoretical outcomes through numerical experiments, which confirm the effectiveness and efficiency of our proposed algorithm.

We investigate the problem of approximating the matrix function f(A) by r(A), with f a Markov function, r a rational interpolant of f, and A a symmetric Toeplitz matrix. In a first step, we obtain a new upper bound for the relative interpolation error 1 − r/f on the spectral interval of A. By minimizing this upper bound over all interpolation points, we obtain a new, simple, and sharp a priori bound for the relative interpolation error. We then consider three different approaches of representing and computing the rational interpolant r. Theoretical and numerical evidence is given that any of these methods for a scalar argument allows to achieve high precision, even in the presence of finite precision arithmetic. We finally investigate the problem of efficiently evaluating r(A), where it turns out that the relative error for a matrix argument is only small if we use a partial fraction decomposition for r following Antoulas and Mayo. An important role is played by a new stopping criterion which ensures to automatically find the degree of r leading to a small error, even in presence of finite precision arithmetic. [ABSTRACT FROM AUTHOR]

Identifying the most important nodes according to specific centrality indices is an important issue in network analysis. Node metrics based on the computation of functions of the adjacency matrix of a network were defined by Estrada and his collaborators in various papers. This paper describes a MATLAB toolbox for computing such centrality indices using efficient numerical algorithms based on the connection between the Lanczos method and Gauss-type quadrature rules. [ABSTRACT FROM AUTHOR]

Let T be a square matrix with a real spectrum, and let f be an analytic function. The problem of the approximate calculation of f(T) is discussed. Applying the Schur triangular decomposition and the reordering, one can assume that T is triangular and its diagonal entries tii are arranged in increasing order. To avoid calculations using the differences tii − tjj with close (including equal) tii and tjj, it is proposed to represent T in a block form and calculate the two main block diagonals using interpolating polynomials. The rest of the f(T) entries can be calculated using the Parlett recurrence algorithm. It is also proposed to perform some scalar operations (such as the building of interpolating polynomials) with an enlarged number of significant decimal digits. [ABSTRACT FROM AUTHOR]

Many dimensionality reduction methods in the manifold learning field have the so-called small-sample-size (SSS) problem. Starting from solving the SSS problem, we first summarize the existing dimensionality reduction methods and construct a unified criterion function of these methods. Then, combining the unified criterion with the matrix function, we propose a general matrix function dimensionality reduction framework. This framework is configurable, that is, one can select suitable functions to construct such a matrix transformation framework, and then a series of new dimensionality reduction methods can be derived from this framework. In this article, we discuss how to choose suitable functions from two aspects: 1) solving the SSS problem and 2) improving pattern classification ability. As an extension, with the inverse hyperbolic tangent function and linear function, we propose a new matrix function dimensionality reduction framework. Compared with the existing methods to solve the SSS problem, these new methods can obtain better pattern classification ability and have less computational complexity. The experimental results on handwritten digit, letters databases, and two face databases show the superiority of the new methods. [ABSTRACT FROM AUTHOR]

Existing algorithms for computing the matrix cosine are tightly coupled to a specific precision of floating-point arithmetic for optimal efficiency so they do not conveniently extend to an arbitrary precision environment. We develop an algorithm for computing the matrix cosine that takes the unit roundoff of the working precision as input, and so works in an arbitrary precision. The algorithm employs a Taylor approximation with scaling and recovering and it can be used with a Schur decomposition or in a decomposition-free manner. We also derive a framework for computing the Fr'echet derivative, construct an efficient evaluation scheme for computing the cosine and its Fréchet derivative simultaneously in arbitrary precision, and show how this scheme can be extended to compute the matrix sine, cosine, and their Fréchet derivatives all together. Numerical experiments show that the new algorithms behave in a forward stable way over a wide range of precisions. The transformation-free version of the algorithm for computing the cosine is competitive in accuracy with the state-of-the-art algorithms in double precision and surpasses existing alternatives in both speed and accuracy in working precisions higher than double. [ABSTRACT FROM AUTHOR]

This work is concerned with approximating matrix functions for banded matrices, hierarchically semiseparable matrices, and related structures. We develop a new divide-and-conquer method based on (rational) Krylov subspace methods for performing low-rank updates of matrix functions. Our convergence analysis of the newly proposed method proceeds by establishing relations to best polynomial and rational approximation. When only the trace or the diagonal of the matrix function is of interest, we demonstrate-in practice and in theory-that convergence can be faster. For the special case of a banded matrix, we show that the divide-and-conquer method reduces to a much simpler algorithm, which proceeds by computing matrix functions of small submatrices. Numerical experiments confirm the effectiveness of the newly developed algorithms for computing large-scale matrix functions from a wide variety of applications. [ABSTRACT FROM AUTHOR]

Matrix functions play important roles in various branches of science and engineering. In numerical computations and physical measurements there are several sources of error which significantly affect the main results obtained from solving the problems. This effect also influences the matrix computations. In this paper, we propose some approaches to enclose the matrix functions. We then present some analytical arguments to ensure that the obtained enclosures contain the exact result. Numerical experiments are given to illustrate the performance and effectiveness of the proposed approaches. [ABSTRACT FROM AUTHOR]

This paper considers the computation of approximations of matrix functionals of form F(A):=vTf(A)v, where A is a large symmetric positive definite matrix, v is a vector, and f is a Stieltjes function. The functional F(A) is approximated by a rational Gauss quadrature rule with poles on the negative real axis (or part thereof) in the complex plane, and we focus on the allocation of the poles. Specifically, we propose that the poles, when considered positive point charges, be allocated to make the negative real axis (or part thereof) approximate an equipotential curve. This is easily achieved with the aid of conformal mapping.

The Fréchet derivative Lf(A,E) of the matrix function f(A) plays an important role in many different applications, including condition number estimation and network analysis. We present several different Krylov subspace methods for computing low‐rank approximations of Lf(A,E) when the direction term E is of rank one (which can easily be extended to general low rank). We analyze the convergence of the resulting methods both in the Hermitian and non‐Hermitian case. In a number of numerical tests, both including matrices from benchmark collections and from real‐world applications, we demonstrate and compare the accuracy and efficiency of the proposed methods. [ABSTRACT FROM AUTHOR]

In this work, we introduce and investigate a class of matrix-free regularization techniques for discrete linear ill-posed problems based on the approximate computation of a special matrix-function. In order to produce a regularized solution, the proposed strategy employs a regular approximation of the Heavyside step function computed into a small Krylov subspace. This particular feature allows our proposal to be independent from the structure of the underlying matrix. If on the one hand, the use of the Heavyside step function prevents the amplification of the noise by suitably filtering the responsible components of the spectrum of the discretization matrix, on the other hand, it permits the correct reconstruction of the signal inverting the remaining part of the spectrum. Numerical tests on a gallery of standard benchmark problems are included to prove the efficacy of our approach even for problems affected by a high level of noise. [ABSTRACT FROM AUTHOR]

We describe a complete theory for walk-based centrality indices in complex networks defined in terms of Mittag-Leffler functions. This overarching theory includes as special cases wellknown centrality measures like subgraph centrality and Katz centrality. The indices we introduce are parametrized by two numbers; by letting these vary, we show that Mittag-Leffler centralities interpolate between degree and eigenvector centrality, as well as between resolvent-based and exponential-based indices. We further discuss modelling and computational issues, and provide guidelines on parameter selection. The theory is then extended to the case of networks that evolve over time. Numerical experiments on synthetic and real-world networks are provided. [ABSTRACT FROM AUTHOR]

The need to compute matrix functions of the form f(A)v, where A ∈ ℝ N × N is a large symmetric matrix, f is a function such that f(A) is well defined, and v≠ 0 is a vector, arises in many applications. This paper is concerned with the situation when A is so large that the evaluation of f(A) is prohibitively expensive. Then, an approximation of f(A)v often is computed by applying a few, say 1 ≤ n ≪ N, steps of the symmetric Lanczos process to A with initial vector v to determine a symmetric tridiagonal matrix T n ∈ ℝ n × n and a matrix V n ∈ ℝ N × n , whose orthonormal columns span a Krylov subspace. The expression Vnf(Tn)e1∥v∥ furnishes an approximation of f(A)v. The evaluation of f(Tn) is inexpensive, because the matrix Tn is small. It is important to be able to estimate the error in the computed approximation. This paper describes a novel approach that is based on a technique proposed by Spalević for estimating the error in Gauss quadrature rules. [ABSTRACT FROM AUTHOR]

The k th Fréchet derivative of a matrix function f is a multilinear operator from a cartesian product of k subsets of the space ℂ n × n into itself. We show that the k th Fréchet derivative of a real-valued matrix function f at a real matrix A in real direction matrices E1, E2, ... , Ek can be computed using the complex step approximation. We exploit the algorithm of Higham and Relton (SIAM J. Matrix Anal. Appl.35(3):1019–1037, 2014) with the complex step approximation and mixed derivative of complex step and central finite difference scheme. Comparing with their approach, our cost analysis and numerical experiment reveal that half and seven-eighths of the computational cost can be saved for the complex step and mixed derivative, respectively. When f has an algorithm that computes its action on a vector, the computational cost drops down significantly as the dimension of the problem and k increase. [ABSTRACT FROM AUTHOR]

The Schur--Parlett algorithm, implemented in MATLAB as funm, evaluates an analytic function f at an n X n matrix argument by using the Schur decomposition and a block recurrence of Parlett. The algorithm requires the ability to compute f and its derivatives, and it requires that / have a Taylor series expansion with a suitably large radius of convergence. We develop a version of the Schur--Parlett algorithm that requires only function values and not derivatives. The algorithm requires access to arithmetic of a matrix-dependent precision at least double the working precision, which is used to evaluate f on the diagonal blocks of order greater than 2 (if there are any) of the reordered and blocked Schur form. The key idea is to compute by diagonalization the function of a small random diagonal perturbation of each diagonal block, where the perturbation ensures that diagonalization will succeed. Our algorithm is inspired by Davies's randomized approximate diago-nalization method, but we explain why that is not a reliable numerical method for computing matrix functions. This multiprecision Schur--Parlett algorithm is applicable to arbitrary analytic functions / and, like the original Schur--Parlett algorithm, it generally behaves in a numerically stable fashion. The algorithm is especially useful when the derivatives of f are not readily available or accurately computable. We apply our algorithm to the matrix Mittag--Leffler function and show that it yields results of accuracy similar to, and in some cases much greater than, the state-of-the-art algorithm for this function. [ABSTRACT FROM AUTHOR]

We discuss the general Rodrigues problem and we give explicit determinant formulae for the coefficients when the eigenvalues of the matrix have double multiplicity (Theorem 5). When n = 4 explicit formulae and effective computations for the exponential map on the Lie group SO(4) are presented. [ABSTRACT FROM AUTHOR]

Identifying the most important nodes according to specific centrality indices is an important issue in network analysis. Node metrics based on the computation of functions of the adjacency matrix of a network were defined by Estrada and his collaborators in various papers. This paper describes a MATLAB toolbox for computing such centrality indices using efficient numerical algorithms based on the connection between the Lanczos method and Gauss-type quadrature rules.