Efﬁcient Evaluation of Matrix Polynomials beyond the Paterson–Stockmeyer Method

: Recently, two general methods for evaluating matrix polynomials requiring one matrix product less than the Paterson–Stockmeyer method were proposed, where the cost of evaluating a matrix polynomial is given asymptotically by the total number of matrix product evaluations. An analysis of the stability of those methods was given and the methods have been applied to Taylor-based implementations for computing the exponential, the cosine and the hyperbolic tangent matrix functions. Moreover, a particular example for the evaluation of the matrix exponential Taylor approximation of degree 15 requiring four matrix products was given, whereas the maximum polynomial degree available using Paterson–Stockmeyer method with four matrix products is 9. Based on this example, a new family of methods for evaluating matrix polynomials more efﬁciently than the Paterson–Stockmeyer method was proposed, having the potential to achieve a much higher efﬁciency, i.e., requiring less matrix products for evaluating a matrix polynomial of certain degree, or increasing the available degree for the same cost. However, the difﬁculty of these family of methods lies in the calculation of the coefﬁcients involved for the evaluation of general matrix polynomials and approximations. In this paper, we provide a general matrix polynomial evaluation method for evaluating matrix polynomials requiring two matrix products less than the Paterson-Stockmeyer method for degrees higher than 30. Moreover, we provide general methods for evaluating matrix polynomial approximations of degrees 15 and 21 with four and ﬁve matrix product evaluations, respectively, whereas the maximum available degrees for the same cost with the Paterson–Stockmeyer method are 9 and 12, respectively. Finally, practical examples for evaluating Taylor approximations of the matrix cosine and the matrix logarithm accurately and efﬁciently with these new methods are given.


Introduction
The authors of [1] presented a new family of methods for the evaluation of matrix polynomials more efficiently than the state-of-the-art method from [2] by Paterson and Stockmeyer (see [3], Section 4.2). These methods are based on the multiplication of matrix polynomials to get a new matrix polynomial with degree given by the sum of the degrees of the original matrix polynomials. The main difficulty in these methods lies in obtaining the coefficients involved for the evaluation of general matrix polynomials. In this sense, the authors of [1] (Section 3) gave two concrete general methods for evaluating matrix polynomials requiring one less matrix product than the Paterson-Stockmeyer method. Regarding the cost of evaluating matrix polynomials, since the cost of a matrix product, denoted by M, is O(n 3 ) for n × n matrices, and both the cost of the sum of two matrices and the cost of a product of a matrix by a scalar are O(n 2 ), similarly to [3]  Evaluating matrix polynomials of degrees greater than 30 with two matrix products less than the Paterson-Stockmeyer method.
Finally, examples for computing Taylor approximations of the matrix cosine and the matrix logarithm efficiently and accurately using those evaluation formulae are given.
Regarding Taylor approximations, if is the Taylor series of the matrix function f (·), where X ∈ C n×n , then is its Taylor approximation of order m (for the convergence of matrix Taylor series, see Theorem 4.7 of [3], p. 76). From [11] (Section 1), a matrix X ∈ C n×n is a logarithm of B ∈ C n×n if e X = B. Therefore, any nonsingular matrix has infinitely many logarithms and we will focus on the principal logarithm, denoted by log(B). For a matrix B ∈ C n×n with no eigenvalues on R − the principal logarithm is the unique logarithm whose eigenvalues have imaginary parts lying in the interval (−π, π). Therefore, in the given examples, we will assume that B has no eigenvalues on R − and we will take the logarithm Taylor series (1) The exponential matrix has been studied in numerous papers (see [3] (Chap. 10), and [5,6,8,12] and the references therein). This matrix function can be defined by The matrix cosine has received attention recently (see [4,7] and the references therein). This matrix function can be defined by Note that if we truncate the Taylor series on the right-hand side of (3) by the term i = m, then the order of the corresponding cosine Taylor approximation is 2m.
Regarding the cost in matrix rational approximations, note that the multiplication by the corresponding matrix inverse is calculated by solving a multiple right-hand side linear system. From [13] (Appendix C), it follows that the cost of the solution of multiple right-hand side linear systems AX = B, where matrices A and B are n × n, denoted by D (see [14], p. 11940) is Therefore, using (4), the cost of computing rational approximations will be also given in terms of M.
In this article, the following notation will be used: x denotes the smallest integer greater than or equal to x, and x the largest integer less than or equal to x. u denotes the unit roundoff in IEEE double precision arithmetic (see [15], Section 2.2). The set of positive integers is denoted as N. The set of real and complex matrices of size n × n are denoted, respectively, by R n×n and C n×n . The identity matrix for both sets is denoted as I. The dependence of a variable y on the variables x 1 , x 2 , . . . , x n is denoted by y = y(x 1 , x 2 , . . . , x n ).
In Section 2, we recall some results for computing matrix polynomials using the Paterson-Stockmeyer method and summarize the matrix polynomial evaluation methods from [1]. In Section 3, we describe the general methods for computing polynomial approximations of degrees 15, 21, and 6s with s = 3, 4, . . . and give examples for the Taylor approximation of the cosine and logarithm matrix functions. Finally, in Section 4, we give some conclusions. In this paper, we provide a method to evaluate matrix polynomials with two matrix products less than the Paterson-Stockmeyer method and one matrix product less than the methods from [1] (Section 3). Moreover, in this paper, we provide methods to evaluate polynomial approximations of matrix functions of degrees 15 and 21 with cost 3M and 4M. These methods are interesting because the maximum available degrees using the other method proposed in this paper are 12 and 18, respectively. All of the new methods proposed can be used in the applications for computing approximations of matrix functions or evaluating matrix polynomials more efficiently than using the state-of-the-art methods.

Paterson-Stockmeyer Method
The Paterson-Stockmeyer method [2] for computing a matrix polynomial consists of calculating P m (A) as where PS m (A) denotes the Paterson-Stockmeyer evaluation Formula (6) and s > 0 is an integer that divides m. Given a number of matrix products, the maximum degrees of P m (A) that are available using the Paterson-Stockmeyer method are the following: m = s 2 , and m = s(s + 1), where s ∈ N, denoted by m * , m * = {1, 2, 4, 6, 9, 12, . . .} [14] (Section 2.1). The cost C PS for computing (6) for the values of m * are given by [14] (Equation (5)), which appear in [14] ( Table 1). In [16], the optimality of the rule m * = (C PS − s + 2)s, where s = C PS /2 + 1, was demonstrated. This rule gives the same results as (7), since if C PS is even then C PS = 2s − 1, and in that case m * = s(s + 1), and if C PS is odd then C PS = 2s, and then m * = s 2 . Note that, for positive integers m / ∈ m * , P m (A) = PS m 0 (A) can be evaluated using (6) taking m 0 = min{m 1 ∈ m * , m 1 > m} and setting some coefficients as zero [1] (Section 2.1).

General Polynomial Evaluation Methods beyond the Paterson-Stockmeyer Method
The authors of [1] (Example 3.1) give a method to compute P 8 (A) from (5) with a cost of 3M with the following evaluation formulae where q 4 , q 3 , r 2 , r 1 , s 2 , s 0 , t 2 , t 1 , and t 0 are complex numbers. In order to compute (5) with m = 8, if we equate y 12 (A) = P m (A) from (5), then the system of eight equations with eight coefficients from (16)-(24) from [1] arises. In this system, some coefficients can be obtained directly from the polynomial P m (A) coefficients as and the remaining equations can be reduced by variable substitution to a quadratic equation on s 2 . This equation gives two solutions for q 4 = √ c 8 and two more solutions for q 4 = − √ c 8 . The remaining coefficients can be obtained from s 2 , q 4 , and q 3 . From (11), one gets for coefficient c 8 in P m (A) from (5).
In order to check the stability of the solutions of q i , r i , and s i rounded to IEEE double precision arithmetic, the authors of [1] (Example 3.1) proposed to compute the relative error for each coefficient c i , for i = 3, 4, . . . , 8 substituting those solutions into the original system of Equations (16)-(24) from [1]. For instance, from (10), it follows that the relative error for c 8 using q 4 rounded to IEEE double precision arithmetic is  Table 4), one of the solutions rounded to IEEE double precision arithmetic for evaluating the Taylor polynomial of the exponential and cosine functions is shown. These solutions were substituted into the original system of equations to calculate the relative error for c i , for i = 3, 4, . . . , 8 (see [1], Example 3.1), giving a relative error of order u, turning out to be stable solutions. Moreover, the numerical tests from [1] (Example 3.2) and [4,5] also show that if the relative error for each coefficient is O(u), then the polynomial evaluation formulae are accurate, and if the relative errors are O(10u) or greater, then the polynomial evaluation formulae are not so accurate.
The authors of [1] (Section 3) also provided a more general method for computing matrix polynomials P m (A) from (5) of degree m = 4s based on the evaluation formulae where s ≥ 2, q s+i , r i , s i and t i are scalar coefficients, q 2s = ± √ c 4s = 0 and then c 4s = 0 for coefficient c 4s from P m (A). Note that A i , i = 2, 3, . . . , s are computed only once. The degree and computing cost of y 1s (A) are given by (36) of [1], i.e., d y 1s = 4s and C y 1s = s + 1, s = 2, 3, . . ., respectively. A general solution for the coefficients in (16) and (17) is given in [1] (Section 3), with the condition c 4s = 0.
Given a cost C(M), the maximum orders that can be reached when using the Formulae (16) and (17) and the Paterson-Stockmeyer method are shown in [1] ( Table 5).
where k = 1, p is a multiple of s and y ks (x) = y 1s (x) is evaluated using (16) and (17). This allows one to increase the degree of the polynomial to be evaluated. The degree of z 1ps (A) and its computational cost are given by (53) of [1], i.e., d z 1ps = 4s + p, C z 1ps = (1 + s + p/s)M, respectively. Ref.
[1] ( Table 6) shows that evaluating a matrix polynomial using (19) requires one less product than using the Paterson-Stockmeyer Formula (6). Proposition 2 from [1] (Section 5) gives general formulae more efficient than the formulae of the previous methods, whenever at least one solution for the coefficients in (62)-(65) from [1] (Prop. 2) exists so that y ks (x) is equal to the polynomial P m to evaluate. The maximum polynomial degree and the computing cost if x = A, A ∈ C n×n , are given by (66) of [1], i.e., d y ks = 2 k+1 s, C y ks = (s + k) where d y ks increases exponentially while C y ks increases linearly. (17) is a particular case of (65) from [1] where k = 1.

Three General Expressions for y 2s (A)
This section gives general procedures to obtain the coefficients of y 2s (A) from (65) from [1] with k = 2, generalizing the results from [5] (Section 3) for the evaluation of the matrix exponential Taylor approximations of degrees 15, 21, 24, and 30, also giving formulae for evaluating matrix polynomials of orders 6s, where s = 2 , 3, . . . The following proposition allows to compute polynomial approximations of order 15 with cost 4M. Note that from [1] (Table 8), the maximum available order with cost 4M is 9 for the Paterson-Stockmeyer method and 12 for the method given by (16) and (17).
and let P 15 (A) be a polynomial of degree 15 with coefficients b i Then, where coefficients a i are functions of the following variables a i = a i (c 8 , c 7 , . . . , c 2 , d 2 , d 1 , e 1 , e 0 , f 0 , g 0 , h 2 , h 1 , h 0 ), i = 0, 1, . . . , 16, and there exist at least one set of values of the 16 coefficients c 8 , c 7 , . . ., and provided the following conditions are fulfilled: Proof of Proposition 1. Note that y 12 (A) from (20) is a matrix polynomial of degree 8. Therefore, if condition (26) holds, then Example [1] (Example 3.1) gives four possible solutions for evaluating y 12 (A) using the evaluation Formulas (8) and (9) with cost 3M. Similarly to [5] (Section 3.2), we will denote these four solutions as nested solutions.
All of the coefficients c 7 , c 6 , . . ., c 3 , d 2 , d 1 , e 1 , e 0 , f 0 , g 0 , can be obtained from c 2 , c 8 and b i , i = 0, 1, . . . 15, and then h i can be obtained from the three last equations of system (30) as Finally, using (20) and (21) For each of those solutions, coefficient a 16 from y 22 (A) in (23) is given by (25). For the particular case of the matrix exponential Taylor approximation from [5] (p. 209), there were two real solutions of c 8 giving Therefore, we selected the first solution (34) since both solutions were stable according to the stability study from Section 2.2 (see [1], p. 243), but (34) had a lower error for a 16 with respect to the corresponding Taylor coefficient 1/16!. Then, considering exact arithmetic, one gets that the matrix exponential approximation from y 22 (A) in evaluation Formulas (10)-(12) from [5] (p. 209) with the coefficients from [5] (Table 3) is more accurate than the exponential Taylor approximation of order 15. For that reason, the corresponding Taylor approximation order was denoted by m = 15+ in [5] (Section 4).
Recently, in [17], an evaluation formula of the type given in Proposition 1 was used to evaluate a Taylor polynomial approximation of degree 15+ of the hyperbolic tangent. However, in this case, all the solutions obtained were complExample We tried different configurations of the evaluation formulae giving degree 15+, but all of them gave complex solutions. Then, we proposed the similar evaluation Formula (11) from [17] (p. 6) with degree 14+ that did give real solutions. Similarly to (34), in the case of the hyperbolic tangent, the relative error of the coefficients a i , i = 15, and 16 was also lower than 1concretely, 0.38 and 0.85, respectively (see [17], p. 6). This method was compared to the Paterson-Stockmeyer method being noticeably more efficient without affecting the accuracy (see [17], Section 3) for details. Proposition 1 allows us to evaluate polynomial approximations of degree 15 not only for the matrix exponential or the hyperbolic tangent but also for other matrix functions. If all the given solutions were complex, we can modify the formula to evaluate approximation formulae with a lower degree, such as 14+, to check if they give real solutions.
Example 1. In [4] (Section 2), we showed that the solutions for the coefficients of the polynomial evaluation method similar to [5] (Section 3.2) of the matrix cosine Taylor approximation of order 2m = 30+ were not stable, giving poor accuracy results. Using Proposition 1, this example gives a stable solution for calculating a Taylor-based approximation of the matrix cosine with a combination of formula (21) with the Paterson-Stockmeyer method from (19). Setting k = p = s = 2 in (19) and y ks = y 22 from (21), one gets The real solutions of system (30) rounded to IEEE double precision arithmetic explored in [4] (Section 2) gave errors of order ≥ 10 −14 , greater than the unit roundoff in IEEE double precision arithmetic u = 2 −53 ≈ 1.11 × 10 −16 . Using MATLAB code fragments 4.1 and 4.2, we checked that there is no solution with a lower error. Then, according to the stability check from Section 2.2, the solutions are unstable, and we checked in [4] that they gave poor accuracy results. However, using Proposition 1, for 2m = 34+, we could find two real solutions of system (30) giving a maximum error of order u. For those two solutions, a 18 gave respectively. Therefore, the solution (37) giving the lowest error was selected. Table 1 gives the corresponding coefficients in IEEE double precision arithmetic from (8) and (9) for computing (20) with three matrix products, and the rest of the needed coefficients for computing y 22 (B) from (21) with s = 2, given finally by Using (39)-(41) with the coefficients from Table 1 and (36), a matrix cosine Taylor approximation of order 2m = 34+ can be computed in IEEE double precision arithmetic with a cost of six matrix products, i.e., B = A 2 , B 2 , three for evaluating (39)-(41), and one more for evaluating (36). The maximum available and stable order given in [4] (Section 2) with six matrix products was 2m = 30. The coefficients from Table 1 were computed with variable precision arithmetic with a precision of 32 and 250 decimal digits to check its correctness, giving the same results.
Taking into account (3) and the selection of the solution in (37), in exact arithmetic, one gets where, using (39)-(41), one gets a 18 = q 4 4 .  (Table 1) corresponding to the backward relative error analysis of the Taylor approximation of the matrix cosine, denoted by E b . Then, if ||B|| = ||A 2 || ≤ Θ, then ||E b || ≤ u for the corresponding Taylor approximations. In [15] (Table 1), Θ for Taylor approximation of order 16 was 9.97 and Θ 20 = 10.18, showing two decimal digits. Then, for our test with order 2m = 34+, we used a set of 48 8 × 8 matrices from the Matrix Computation Toolbox [18] divided by random numbers to give B between 9 and 10. We compared the forward error E f of both functions where function f (A) was cosm and the function using z 222 (B). The "exact value" of cos(A) was computed using the method in [19]. The total cost of the new matrix cosine computation function z 222 summing up the number of matrix products over all the test matrices is denoted by Cost z 222 . Taking into account (4), the cost for the cosm Padé function summing up the number of matrix products and inversions over all the test matrices is denoted by Cost cosm . Then, the following cost comparison was obtained for that set of test matrices i.e., the cost of z 222 is 40.78% lower than the cost of cosm. Moreover, the results were more accurate in 76.60% of the matrices. Therefore, the new formulae are efficient and accurate.

Evaluation of Matrix Polynomial Approximations of Order 21
In this section, we generalize the results from [5] (Section 3.3) for evaluating polynomial approximations of order m = 21 with cost 5M. Note that for that cost, from [1] (Table 8), the maximum available orders using the Paterson-Stockmeyer method and the evaluation Formulas (16) and (17) are 12 and 16, respectively. Applying a similar procedure to that in Section 3.1 to obtain the coefficients for evaluating a matrix polynomial approximation of order 21, in this case, a system of 22 equations with 22 unknown variables arises. This system can be reduced to three equations with three unknowns using variable substitution with the MATLAB Symbolic Toolbox, provided that two of the variables are not zero. The following proposition summarizes the results Proposition 2. Let y 13 (A) and y 23 (A) be and let P 21 (A) be a polynomial of degree 21 with coefficients b i Then, where coefficients a i = a i (c 12 , c 11 , . . . , c 2 , d 3 , d 2 , d 1 , e 1 , e 0 , f 0 , g 0 , h 3 , h 2 , h 1 , h 0 ), i = 0, 1, . . . , 24, and the system of equations arising when equating can be reduced to a system of three equations of variables c 12 , c 11 and c 10 , provided and then variables a i , i = 22, 23 and 24 are a 24 = c 2 12 , a 23 = 2c 11 c 12 , a 22 = c 2 11 + 2c 10 c 12 .
Proof of Proposition 2. The proof of Proposition 2 is similar to the proof of Proposition 1. Analogously, if condition (18) is fulfilled with s = 3, i.e., c 12 = 0, then polynomial y 13 (A) can be evaluated using (16) and (17) with s = 3 and cost 4M, where y 03 is given by (21) of [5] (Section 3.3), i.e., If we apply (48), we obtain a similar system to (30). Using variable substitution with the MATLAB Symbolic Toolbox, the MATLAB code coeffspolm21plus.m (http:// personales.upv.es/jorsasma/Software/coeffspolm21plus.m (accessed on 24 June 2021)) similar to MATLAB code fragments 4.1 and 4.2 is able to reduce the whole nonlinear system of 22 equations to a nonlinear system of three equations with three variables c 10 , c 11 , and c 12 . The MATLAB code coeffspolm21plus.m returns conditions (49) (see the actual code for details.) If there is at least one solution for c 10 , c 11 , and c 12 fulfilling condition (49), all of the other coefficients can be obtained using the values of c 10 , c 11 , c 12 . Then, y 13 (A) from (44) can be evaluated using (16) and (17) giving several possible solutions. Finally, the solutions are rounded to the required precision. Then, according to the stability study from Section 2.2 (see [1], p. 243), the solution giving the least error should be selected.
Similarly to (34) and (35), the degree of y 23 (A) of (45) is 24, but with the proposed method, we can only set the polynomial approximation coefficients of (46) up to order m = 21. The coefficients of a i of the power A i , i =22, 23, and 24 are given by (50). The authors of [5] (Section 3.3) give one particular example of this method for calculating a matrix Taylor approximation of the exponential function, where in exact arithmetic where T 21 is the Taylor approximation of order m = 21 of the exponential function and showing three decimal digits. Again, in exact arithmetic, the approximation y 23 (A) is more accurate than T 21 (A). Therefore, the order of that approximation was denoted as m = 21+ in [5] (Section 4). The experimental results from [5] showed that this method was more accurate and efficient than the Padé method from [6]. Recently, in [17], an evaluation formula similar to (45) was used to evaluate a Taylor polynomial approximation of the hyperbolic tangent. Similarly to (53), in the case of the hyperbolic tangent, the relative error of the coefficients a i , i = 22, 23, and 24 was also lower than 1-concretely, 0.69, 0.69, and 0.70, respectively (see [17], p. 7). This method was compared to the Paterson-Stockmeyer method being noticeably more efficient without affecting the accuracy (see [17], Section 3 for details).
Proposition 2 allows us to evaluate polynomial approximations of degree 21 not only for the matrix exponential or the hyperbolic tangent but also for other matrix functions. In the following example, we show an application for the evaluation of the Taylor approximation of the matrix logarithm.

Example 2.
In this example, we give real coefficients for computing a Taylor-based approximation of the matrix logarithm of order m = 21+ in a stable manner based on the previous results. Evaluating (44) using (16) and (17) with s = 3, and using (45), the following formulae can be used to compute the approximation of order m = 21+ of the principal logarithm log(B) for a square matrix B = I − A with no eigenvalues on R − where the coefficients are numbered correlatively, and using (1), we take The coefficients can be obtained solving first the system of equations arising from (48) with b i = 1/i for i = 1, 2, . . . , 21, b 0 = 0. We used vpasolve (https://es.mathworks.com/help/ symbolic/vpasolve.html (accessed on 24 June 2021)) function from the MATLAB Symbolic Computation Toolbox to solve those equations with variable precision arithmetic. We used the Random option of vpasolve, which allows to obtain different solutions for the coefficients, running it 100 times. The majority of the solutions were complex, but there were two real stable solutions. Then, we obtained the nested solutions for the coefficients of (16) and (17)  Again, we selected the real stable solution given in Table 2. This solution avoids complex arithmetic if the matrix A is real. The relative errors of the coefficients of A 22 , A 23 and A 24 of y 23 (A) with respect to the corresponding Taylor approximation of order 24 of − log(I − A) function are: where a 22 , a 23 , and a 24 are rounded to double precision arithmetic. Then, considering exact arithmetic, one gets which is more accurate than the corresponding Taylor approximation of log(B) of order m = 21. Therefore, similarly to [5] (Section 4), the approximation order of (63) is denoted by m = 21+.   The θ values such that the relative backward errors for the Padé approximations are lower than u are shown in [11] (Table 2.1). The corresponding θ value for the Taylor approximation of log(I − A) of order m = 21+, denoted by θ 21+ , can be computed similarly (see [11] for details), giving θ 21+ = 0.211084493690929, where the value is rounded to IEEE double precision arithmetic.
We compared the results of using (56)-(58) with the coefficient values from Table 2, with the results given by function logm_iss_full from [20]. For that comparison, we used a matrix test set of 43 8 × 8 matrices of the Matrix Computation Toolbox [18]. We reduced their norms so that they are random with a uniform distribution in [0.2, θ 21+ ] in order to compare the Padé approximations of logm_iss_full with the Taylor-based evaluation Formulas (56)-(58) using no inverse scaling in none of the approximations (see [11]).
The "exact" matrix logarithm was computed using the method from [19]. The error of the implementation using Formula (58) was lower than logm_iss_full in 100% of the matrices with a 19.61% lower relative cost in flops. Therefore, evaluation Formulas (56)-(58) are efficient and accurate for a future Taylor-based implementation for computing the matrix logarithm.

Evaluation of Matrix Polynomials of Degree m = 6s
The following proposition generalizes the particular cases of the evaluation of the matrix exponential Taylor approximation with degrees m = 24 and 30 from [5] (Section 3.4) for evaluating general matrix polynomials of degree m = 6s, s = 2, 3, . . .
Corollary 1. If condition (69) holds, then the system of 6s + 1 equations with 7s + 1 variables arising from (70) can be reduced using variable substitution to a system of s equations with s variables, and if there exist at least one solution for that system, then all the coefficients from (64)-(66) can be calculated using the solution of the system.
If there exist at least one solution of system (89), then the values c 3s−k , c 2s−k and c s−k can be calculated for k = 0, . . . , s − 1 (see (72)-(74)), and coefficients f s−k can be calculated for k = 0, . . . , s, using (90), allowing one to obtain all the coefficients from (64)-(66). Table 3, we present the maximum available order for a cost C(M) in the following cases:

Using [1] (Table 6), in
• The Paterson-Stockmeyer evaluation formula. • z kps from (19) with k = 1, denoting the combination of (17) with the Paterson-Stockmeyer formula proposed in [1] (Section 3.1). • z kps from (19) with k = 2, denoting the combination of (66) with the Paterson-Stockmeyer formula, whenever a solution for the coefficients of z 2ps exist.  Table 3 also shows the values of p and s for z 2ps (A) such that s is minimum to obtain the required order, giving the minimum size of the system (89) to solve, i.e., s equations with s unknown variables. Note that it makes no sense to use (66) for s = 1 and cost C = 3M since the order obtained is m = 6s = 6 and for that cost the Paterson-Stockmeyer method obtains the same order. Table 3 shows that evaluation formula z 2ps obtains a greater order than z 1ps for d z 1s > 12. Concretely, for s z 2s ≥ 5, where the available order with z 2s is d z 2s = 30, 36, 42, . . ., z 2ps allows increments 10, 11, 12 . . . of the available order with respect to using the Paterson-Stockmeyer method, and increments of s z2s = 5, 6, 6 . . . with respect to using z 1ps .
In [5], real stable solutions were found for the coefficients of (64)-(66) for the exponential Taylor approximation with degrees 6s with s = 4 and 5, i.e., 24, and 30. The following example deals with the matrix logarithm Taylor approximation.
Note that using the evaluation formulae from Sections 1 and 2 with cost 4M and 5M, one can get an order of approximation 15+ and 21+, respectively, whereas using z 2ps from (19) combining (66) with the Paterson-Stockmeyer method, the orders that can be obtained are lower, i.e., 12 and 18, respectively (see Table 3). Note that for the approximation 15+ where s = 2 (see Section 3.1), one gets order 15+ = (6s + 3)+ and the total degree of the polynomial obtained is 8s = 16. For the approximation 21+ where s = 3, one gets order 21+ = (6s + 3)+ and the total degree of the polynomial degree is 8s = 24. The next step in our research is to extend the evaluation formulae from Propositions 1 and 2 to evaluate polynomial approximations of order (6s + 3)+ of the type Those formulae correspond to a particular case of Formulas (62)-(65) of [1] (Prop. 2) where k = 2. It is easy to show that the degree of y 2s (A) is 8s and the total number of coefficients of y 2s is 6s + 4, i.e., 3s coefficients a i , s coefficients g i , s coefficients h i , s + 1 coefficients l i , and coefficients f 0 , j 0 and k 0 . Using vpasolve in a similar way as in Example 2, we could find solutions for the coefficients of (94)-(96) and (19) so that y 2s (A) and z 2ps allows to evaluate matrix logarithm Taylor-based approximations of orders from 15+ up to 75+. Similarly, we could also find the coefficients for Formulas (94)-(96) to evaluate matrix hyperbolic tangent Taylor approximations of orders higher than 21. Then, our next research step is to show that evaluation Formulas (94)-(96) and its combination with the Paterson-Stockmeyer method from (19) can be used for the general polynomial approximations of matrix functions.

Conclusions
In this paper, we extend the family of methods for evaluating matrix polynomials from [1], obtaining general solutions for new cases of the general matrix polynomial evaluation Formulas (62)-(65) from Proposition 2 from [1] (Section 5). These cases allow to compute matrix polynomial approximations of orders 15 and 21 with a cost of 4M and 5M, respectively, whenever a stable solution for the coefficients exist. Moreover, a general method for computing matrix polynomials of order m = 6s, for s = 3, 4... more efficiently than the methods provided in [1] was provided. Combining this method with the Paterson-Stockmeyer method, polynomials or degree greater than 30 can be evaluated with two matrix products less than using Paterson-Stockmeyer method as shown in Table 3.
Examples for evaluating Taylor approximations of the matrix cosine and the matrix logarithm were given. The accuracy and efficiency results of the proposed evaluation formulae were compared to state-of-the-art Padé algorithms, being competitive for future implementations for computing both functions.
Future work will deal with the generalization of more efficient evaluation formulae based on the evaluation Formulas (62)-(65) from Proposition 2 from [1] (Section 5), its combinations with Paterson-Stockmeyer method (19), and in general, evaluation formulae based on products of matrix polynomials.