- Email: [email protected]
On determinants of cyclic pentadiagonal matrices with Toeplitz structure
Computers and Mathematics with Applications (
)
–
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
On determinants of cyclic pentadiagonal matrices with Toeplitz structure Jiteng Jia a,∗ , Sumei Li b a
School of Mathematics and Statistics, Xidian University, Xi’an, Shaanxi 710071, China
b
School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China
article
abstract
info
Article history: Received 25 September 2016 Received in revised form 13 November 2016 Accepted 23 November 2016 Available online xxxx
Cyclic pentadiagonal matrices with Toeplitz structure have received tremendous attention in recent years. In the current paper, we present a block upper triangular transformation of the cyclic pentadiagonal Toeplitz matrices. By using the transformation, the determinant of an n-by-n cyclic pentadiagonal Toeplitz matrix can be readily evaluated since one just needs to compute the determinant of a 4-by-4 matrix obtained from the transformation. In addition, an efficient numerical algorithm of O(n) is derived for computing nth order cyclic pentadiagonal Toeplitz determinants. Some numerical experiments are given to show the performance of the proposed algorithm. © 2016 Elsevier Ltd. All rights reserved.
Keywords: Cyclic pentadiagonal matrix Toeplitz matrix Block matrix Triangular transformation Computational cost
1. Introduction and objectives In the current paper, we mainly consider the determinant of a cyclic pentadiagonal Toeplitz matrix of the form d b
e A=
c a
a d
c a
c
..
..
..
..
. .
.. ..
e
. . .
.. .. e
c
. . .
.. .. ..
. . .
b e
.. ..
b e
. .
d b
∈ Rn×n , c
(1.1)
a d
where a, b, c, d, e ∈ R, and we assume that c ̸= 0 and n ≥ 5. In fact, similar results can be obtained by the use of the argument given in Section 2 if c = 0. From practical point of view, pentadiagonal matrices and cyclic pentadiagonal matrices frequently arise from boundary value problems (BVPs) involving 4th-order derivatives [1,2], and, computational algorithms for the determinants are linked to the problem of obtaining efficient test for the existence of unique solutions of the PDEs [3,4]. Usually, the determinant of a square matrix can be computed by the well-known Laplace formula [5]. However, algorithms based on such formula for
∗
Corresponding author. E-mail address: [email protected] (J. Jia).
http://dx.doi.org/10.1016/j.camwa.2016.11.031 0898-1221/© 2016 Elsevier Ltd. All rights reserved.
2
J. Jia, S. Li / Computers and Mathematics with Applications (
)
–
large matrices may be time-consuming due to the number of required operations grows quickly. Therefore, more involved techniques have been developed for the determinant evaluation. In order to evaluate the determinant of the cyclic pentadiagonal matrices, some authors have devised numerical (or symbolic) algorithms recently, see [6–8]. The main objective of this paper is to develop an efficient numerical algorithm for the determinant of an n-by-n cyclic pentadiagonal Toeplitz matrix. Recent developments of algorithms for the determinants and permanents of other related matrices, see e.g. [9–16]. The present paper is organized as follows. The main results are given in the next section. In Section 3, we give the results of some numerical experiments to show the performance of our algorithm. Finally we make some concluding remarks in Section 4. 2. Main results In this section, we first present an approach for the determinant of a cyclic pentadiagonal Toeplitz matrix. Based on some elementary optimization methods, we then improve the proposed approach and develop an efficient numerical algorithm with the cost of 7n + O(log n) for evaluating nth order cyclic pentadiagonal Toeplitz determinants. 2.1. An approach based on block upper triangular transformation Let {C1 , C2 , . . . , Cn } be the column vectors of the matrix A as in (1.1). If we successively interchange the first two columns C1 and C2 with other columns to make them as the last two columns, we then obtain a nearly lower triangular Toeplitz matrix T with the column vectors {C3 , C4 , . . . , Cn , C1 , C2 }, thus
c a
d T = b e
e
b e
c
.. .. .. ..
..
.
..
.
..
.
..
.
.
..
.
..
.
..
.
e
. . .
.. ..
b
d b
a d
e
b
e ∈ Rn×n .
. .
c a
d
(2.1)
c
Since we made 2n − 2 column interchanges in the matrix reordering above, we readily obtain the following relation: det(A) = det(T ). We now consider transforming the matrix T (2.1) into a 2-by-2 block upper triangular matrix of the form:
U =
c · In − 4 0
U12 U22
∈ Rn×n ,
(2.2)
where U12 and U22 are matrices of size (n − 4)-by-4 and 4-by-4, respectively. Below, we show that the above matrix can be obtained from the original matrix T and a suitable choice ωi+1,i , ωi+2,i , ωi+3,i , and ωi+4,i of the following matrices:
1 ω21 ω31 Ω1 = ω41 ω51
I
1
,
1 1 1
..
.
i −1
1
Ωi =
ωi+1,i ωi+2,i ωi+3,i ωi+4,i
,
1 1 1 1 In−i−4
1
for i = 2, 3, . . . , n − 4. Here, Ii−1 and In−i−4 are identity matrices of order i − 1 and n − i − 4. For the sake of illustration, we consider the case n = 7 to show the transformation. From (2.1) and setting ω21 = −a/c, ω31 = −d/c, ω41 = −b/c, and ω51 = −e/c, we have
1 −a/c −d/c Ω1 T = −b/c −e/c 0 0
0 1 0 0 0 0 0
0 0 1 0 0 0 0
0 0 0 1 0 0 0
0 0 0 0 1 0 0
0 0 0 0 0 1 0
0 c 0 a 0 d 0 b 0 e 0 0 1 0
0 c a d b e 0
0 0 c a d b e
e 0 0 c a d b
b e 0 0 c a d
d b e 0 0 c a
a d b e 0 0 c
J. Jia, S. Li / Computers and Mathematics with Applications (
c 0 0 = 0 0 0 0
0 c a d b e 0
0 0 c a d b e
e t˜24 t˜34 t˜44 t˜54 d b
b t˜25 t˜35 t˜45 t˜55 a d
d t˜26 t˜36 t˜46 t˜56 c a
)
–
3
a t˜27 t˜37 t˜47 . t˜57 0 c
As we can see from the above, multiplication on the left by the matrix Ω1 only changes a few values in the first column and the last four columns of the matrix T (i.e., ti1 → 0, ti,n−3 → t˜i,n−3 , ti,n−2 → t˜i,n−2 , ti,n−1 → t˜i,n−1 , tin → t˜in , i = 2, 3, 4, 5). Similarly, choosing ω32 = ω43 = −a/c, ω42 = ω53 = −d/c, ω52 = ω63 = −b/c, ω62 = ω73 = −e/c, we finally have the following 7-by-7 matrix c 0 0 Ω3 Ω2 Ω1 T = 0 0 0 0
0 c 0 0 0 0 0
0 0 c 0 0 0 0
e t˜24 t˜34 t˜44 t˜54 t˜64 t˜74
b t˜25 t˜35 t˜45 t˜55 t˜65 t˜75
d t˜26 t˜36 t˜46 t˜56 t˜66 t˜76
a t˜27 t˜37 t˜47 . t˜57 t˜67 t˜77
Generally, the values ωi+1,i , ωi+2,i , ωi+3,i , ωi+4,i are defined as
ωi+1,i = −a/c ,
ωi+2,i = −d/c ,
ωi+3,i = −b/c ,
ωi+4,i = −e/c ,
i = 1, 2, . . . , n − 4.
Denoting Ω = Ωn−4 Ωn−5 · · · Ω1 , it follows that
ΩT = U, where U is the 2-by-2 block upper triangular matrix defined in (2.2). Since we can readily see from the definition of Ωi that det(Ω ) = 1, we directly deduce the following relation: det(A) = det(T ) = det(Ω T ) = det(U ) = c n−4 · det(U22 ). Hence, the determinant of A can be implicitly evaluated by computing the determinant of the 4-by-4 matrix U22 . Now, we are ready for describing the summary of our approach. The summary is given below. Step 1. Reorder the columns in the original cyclic pentadiagonal Toeplitz matrix to form a nearly lower triangular matrix with Toeplitz structure. Step 2. Transform the nearly lower triangular matrix into a 2-by-2 block upper triangular matrix. Step 3. Compute the determinant of the block upper triangular matrix. If we measure the computational cost of the approach in terms of total number of flops, where each flop represents one of the basic arithmetic floating point operations, then, the above approach may give det(A) in 29n + O(1) operations, since costs for the steps 2 and 3 are 28n + O(1) and n + O(1), respectively. It should be mentioned that the costs of our approach are more expensive than those of some existing methods [7,8]. In the next section, we propose a more efficient numerical algorithm with the cost of only 7n + O(log n) for evaluating nth order cyclic pentadiagonal Toeplitz determinants. 2.2. An efficient algorithm for the determinant evaluation From the previous subsection, we note that division by constant c appears multiple times, as it is used in transforming the nearly lower triangular matrix T into a block upper triangular matrix. In order to reduce the number of such operations, we set a′ = a/c ,
b′ = b/c ,
c ′ = 1,
d′ = d/c ,
e′ = e/c .
If we use these new constants in the definition of the matrix A (1.1), we obtain another cyclic pentadiagonal Toeplitz matrix A′ . Obviously, we have det(A) = c n · det(A′ ).
(2.3)
For the matrix A′ , we may follow the similar procedure applied to A so that we obtain a nearly lower triangular Toeplitz matrix T ′ corresponding to the matrix T as in (2.1). In fact, T ′ is nothing but T with new entries a′ , b′ , c ′ , d′ and e′ . Also, we can form the matrix Ωi′ corresponding to the matrix Ωi . Denoting Ω ′ = Ωn′ −4 Ωn′ −5 · · · Ω1′ , it yields
Ω ′T ′ =
In−4 0
′ U12 ′ U22
∈ R n×n ,
(2.4)
4
J. Jia, S. Li / Computers and Mathematics with Applications (
)
–
′ ′ where U12 and U22 are matrices of size (n − 4)-by-4 and 4-by-4, respectively. In addition, based on the structure of Ωi′ and simple matrix multiplication, we can readily obtain
1
(n−4)
w1 (n−4) w 2 (n−4) w3 ′ Ω = . . . . .. (n−4) wn−2
1
w1(n−5) w2(n−5) .. . .. . wn(n−−35) wn(n−−25)
wn(n−−14)
(1)
where w1
(k)
wi
(1)
..
.
..
.
1
..
.
w1(1)
..
.
w2(1)
···
(1)
′
1 1
w3
1
(1)
w4
··· (1)
= −a , w2 = −d , w3 (k−1) , wi ( k − 1 ) − a′ wk(k−)1 , w k = wk(k+−11) − d′ wk(k−)1 , wk(k+−21) − b′ wk(k−)1 , ′ (k) −e wk−1 , ′
,
(2.5)
1 (1)
= −b , w4 = −e , and ′
′
for i = 1, 2, . . . , k − 1, for i = k, for i = k + 1, for i = k + 2, for i = k + 3,
for k = 2, 3, . . . , n − 4. Finally, it follows from (2.3)–(2.5) that we have ′ det(A) = c n · det(A′ ) = c n · det(Ω ′ T ′ ) = c n · det(U22 ),
where
′ U22
(n−4) wn−4 (n−4) wn−3 = (n−4) wn−2
wn(n−−14)
wn(n−−55)
wn(n−−66)
wn(n−−45)
wn(n−−56)
wn(n−−35)
wn(n−−46)
wn(n−−25)
wn(n−−36)
wn(n−−77) ′ e wn(n−−67) 0 wn(n−−57) 0 wn(n−−47)
0
b′ e′ 0 0
d′ b′ e′ 0
1 a′ d ′ a′ + d′ b′ ′ b′ e
0 1 a′ d′
0 0 1 a′
0 0 . 0 1
Now, we state a numerical algorithm for evaluating nth order cyclic pentadiagonal Toeplitz determinants. Algorithm 2.1 Step 1. Input the elements a, b, c, d, e and the order n. (1) (1) (1) (1) Step 2. Set w1 = −a/c, w2 = −d/c, w3 = −b/c, w4 = −e/c. For i = 2, 3, . . . , n − 4 compute
w1(i) = w2(i−1) + w1(1) ∗ w1(i−1) , w2(i) = w3(i−1) + w2(1) ∗ w1(i−1) , w3(i) = w4(i−1) + w3(1) ∗ w1(i−1) , w4(i) = w4(1) ∗ w1(i−1) , End. Step 3. Form the following three 4-by-4 matrices
w1(n−4) w1(n−5) w1(n−6) w1(n−7) ( n − 4 ) ( n − 5 ) ( n − 6 ) ( n − 7 ) w w2 w2 w2 2 W1 = w(n−4) w (n−5) w(n−6) w(n−7) , 3 3 3 3 w4(n−4) w4(n−5) w4(n−6) w4(n−7) −w4(1) −w3(1) −w2(1) −w1(1) ( 1 ) ( 1 ) ( 1 ) 0 −w4 −w3 −w2 , W2 = (1) 0 0 −w4 −w3(1) 0 0 0 −w4(1)
J. Jia, S. Li / Computers and Mathematics with Applications (
1
0 1
−w1(1) W3 = −w(1) 2 −w3(1)
0 0 1
−w1(1) −w2(1)
−w1(1)
)
–
5
0 0 , 0 1
then compute B = W1 ∗ W2 + W3 . Step 4. Output the determinant: det(A) = c n · det(B). We can see from the above algorithm that the total computational cost for calculating the determinant of an n-by-n cyclic pentadiagonal Toeplitz matrix is 8n + O(1), since costs for steps 2, 3, and 4 are 7n + O(1), O(1), and n + O(1), respectively. In fact, the computational cost of Algorithm 2.1 can be reduced to 7n + O(log n) since the computation of c n (in Step 4) can be performed in O(log n) operations, for details, see [17, p. 226]. On the other hand, it should be mentioned that the algorithms given in [7,8] require 15n + O(1) and 8n + O(1), respectively. Remark 2.1. Although we work with real cyclic pentadiagonal Toeplitz matrices in this paper, the results and the proposed algorithm are also valid for cyclic pentadiagonal Toeplitz matrices with complex entries (i.e., a, b, c, d, e ∈ C).
3. Numerical examples In this section, we give the results of some simple numerical experiments to confirm that the proposed algorithm works as well as other already existing algorithms. All tests were performed in MATLAB R2014b using double precision arithmetic. Example 3.1. In this example, we evaluate the determinant of an n-by-n symmetric cyclic pentadiagonal Toeplitz matrix originating from [7] 2.5 1.5
1.5 2.5
1 1.5
.
.
1 A=
..
..
1 1.5
.
..
.. ..
1 1
.. ..
.
..
.
. . .
1 1
.. .. ..
. . .
1.5 1
.. ..
. .
2.5 1.5
1.5 1
. 1 1.5 2.5
The results with DETQPT algorithm [7], DCPT algorithm [8], and Algorithm 2.1 are shown in Table 1. From the results shown in Table 1, we see that our algorithm generated almost the same values as those given by other algorithms. Example 3.2. Now, we consider an ill-conditioned nonsymmetric cyclic pentadiagonal Toeplitz matrix that corresponds to a slight modification of Example 3.3 in [8]
1.3 + h2 A=
1.2 0.4
1 1.5
1.5 1.3 + h2
..
..
. .
1 1.5
..
.. ..
. . .
0.4 1
.. .. ..
. . .
0.4 1
.. .. ..
. . .
1.2 0.4
.. ..
1.2 0.4
. .
1.3 + h2 1.2
1 1.5 1.3 + h2
.
Here, we set h = 1/n. For n = 500, 1000, 2000, 3000, 6000, we also use DETQPT algorithm, DCPT algorithm and Algorithm 2.1 to compute the determinants. The absolute differences between the determinants derived by these three algorithms and those of MATLAB built-in function (i.e., AD1 = |DDETQPT − DMATLAB |, AD2 = |DDCPT − DMATLAB |, AD3 = |DAlgorithm 2.1 − DMATLAB |) are shown in Table 2. The condition numbers of the matrix A and the CPU times in the computation of the determinants are also given.
6
J. Jia, S. Li / Computers and Mathematics with Applications (
)
–
Table 1 Numerical results of the determinant for Example 3.1. n
DETQPT algorithm
DCPT algorithm
Algorithm 2.1
100 500 1 000 5 000 8 000 10 000 20 000
0.054821825453014 1.321092500749511 4.702608204485831 0.680884214220257 7.031833797460169 2.569002419156596 8.076085200081945
0.054821825453016 1.321092500749579 4.702608204485895 0.680884214220271 7.031833797460265 2.569002419156551 8.076085200081494
0.054821825453014 1.321092500749508 4.702608204485832 0.680884214220254 7.031833797460159 2.569002419156595 8.076085200081929
Table 2 Numerical results of the determinant for Example 3.2. Algorithms
n Cond(A)
500 1.4126e+019
1000 2.3253e+034
2000 7.8222e+082
3000 4.2945e+131
6000 6.7214e+276
DETQPT
AD1 CPU times
1.7341e−014 0.0146
1.7603e−014 0.0261
4.8559e−014 0.0398
1.7726e−013 0.0712
2.1512e−013 0.1433
DCPT
AD2 CPU times
1.7446e−014 0.0009
1.6206e−014 0.0017
4.7361e−014 0.0026
1.7615e−013 0.0051
2.1192e−013 0.0069
Algorithm 2.1
AD3 CPU times
1.7210e−014 0.0006
1.7269e−014 0.0017
4.8068e−014 0.0028
1.7590e−013 0.0049
2.1217e−013 0.0072
4. Concluding remarks In this paper, first we have presented an approach to computing the determinant of the cyclic pentadiagonal Toeplitz matrices based on block upper triangular transformation. However, the price we should pay is that the costs of the approach are more expensive than those of some existing methods. Second, we have described a numerical algorithm (Algorithm 2.1) for the determinant evaluation. The proposed algorithm is based on the use of a certain type of matrix reordering and equalities involving the products of special kind of relevant matrices, and it leads to a significant decrease in the number of operations. Moreover, we have shown that the complexity of our algorithm is less than or almost equal to those of some recent algorithms. Acknowledgments The authors would like to thank anonymous referees for useful comments that enhanced the quality of the manuscript. This work was supported by the Natural Science Foundation of China (NSFC) under grant 11601408 and the Fundamental Research Funds for the Central Universities under grant JB160701, JBG160704. References [1] S. Barnett, Matrices: Methods and Applications, Oxford University Press, New York, 1990. [2] J. Monterde, H. Ugail, A general 4th-order PDE method to generate Bézier surfaces from the boundary, Comput. Aided Geom. Design 23 (2006) 208–225. [3] R.B. Marr, G.H. Vineyard, Five-diagonal Toeplitz determinants and their relation to Chebyshev polynomials, SIAM J. Matrix Anal. Appl. 9 (4) (1988) 579–586. [4] M.E. Kanal, Parallel algorithm on inversion for adjacent pentadiagonal matrices with MPI, J. Supercomput. 59 (2012) 1071–1078. [5] G.H. Golub, C.F. Van Loan, Matrix Computations, third ed., The Johns Hopkins University Press, Baltimore and London, 1996. [6] J.T. Jia, T. Sogabe, A novel algorithm for solving quasi penta-diagonal linear systems, J. Math. Chem. 51 (2013) 881–889. [7] Y.L. Jiang, J.T. Jia, On the determinant evaluation of quasi penta-diagonal matrices and quasi penta-diagonal Toeplitz matrices, J. Math. Chem. 51 (2013) 2503–2513. [8] J.T. Jia, S.M. Li, Numerical algorithm for the determinant evaluation of cyclic pentadiagonal matrices with Toeplitz structure, Numer. Algorithms 71 (2016) 337–348. [9] R.A. Sweet, A recursive relation for the determinant of a pentadiagonal matrix, Comm. ACM. 12 (1969) 330–332. [10] T. Sogabe, A fast numerical algorithm for the determinant of a pentadiagonal matrix, Appl. Math. Comput. 196 (2008) 835–841. [11] T. Sogabe, A note on A fast numerical algorithm for the determinant of a pentadiagonal matrix, Appl. Math. Comput. 201 (2008) 561–564. [12] Z. Cinkir, An elementary algorithm for computing the determinant of pentadiagonal Toeplitz matrices, J. Comput. Appl. Math. 236 (2012) 2298–2305. [13] J.T. Jia, S.M. Li, On the inverse and determinant of general bordered tridiagonal matrices, Comput. Math. Appl. 69 (2015) 503–509. [14] J.T. Jia, T. Sogabe, S.M. Li, A generalized symbolic Thomas algorithm for the solution of opposite-bordered tridiagonal linear systems, J. Comput. Appl. Math. 290 (2015) 423–432. [15] J.T. Jia, B.T. Yang, S.M. Li, On a homogeneous recurrence relation for the determinants of general pentadiagonal Toeplitz matrices, Comput. Math. Appl. 71 (2016) 1036–1044. [16] M.E.A. El-Mikkawy, A fast and reliable algorithm for evaluating nth order pentadiagonal determinants, Appl. Math. Comput. 202 (2008) 210–215. [17] K.H. Rosen, Discrete Mathematics and its Applications, McGraw-Hill, New York, 2007.
)
–
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
On determinants of cyclic pentadiagonal matrices with Toeplitz structure Jiteng Jia a,∗ , Sumei Li b a
School of Mathematics and Statistics, Xidian University, Xi’an, Shaanxi 710071, China
b
School of Mathematics and Statistics, Xi’an Jiaotong University, Xi’an, Shaanxi 710049, China
article
abstract
info
Article history: Received 25 September 2016 Received in revised form 13 November 2016 Accepted 23 November 2016 Available online xxxx
Cyclic pentadiagonal matrices with Toeplitz structure have received tremendous attention in recent years. In the current paper, we present a block upper triangular transformation of the cyclic pentadiagonal Toeplitz matrices. By using the transformation, the determinant of an n-by-n cyclic pentadiagonal Toeplitz matrix can be readily evaluated since one just needs to compute the determinant of a 4-by-4 matrix obtained from the transformation. In addition, an efficient numerical algorithm of O(n) is derived for computing nth order cyclic pentadiagonal Toeplitz determinants. Some numerical experiments are given to show the performance of the proposed algorithm. © 2016 Elsevier Ltd. All rights reserved.
Keywords: Cyclic pentadiagonal matrix Toeplitz matrix Block matrix Triangular transformation Computational cost
1. Introduction and objectives In the current paper, we mainly consider the determinant of a cyclic pentadiagonal Toeplitz matrix of the form d b
e A=
c a
a d
c a
c
..
..
..
..
. .
.. ..
e
. . .
.. .. e
c
. . .
.. .. ..
. . .
b e
.. ..
b e
. .
d b
∈ Rn×n , c
(1.1)
a d
where a, b, c, d, e ∈ R, and we assume that c ̸= 0 and n ≥ 5. In fact, similar results can be obtained by the use of the argument given in Section 2 if c = 0. From practical point of view, pentadiagonal matrices and cyclic pentadiagonal matrices frequently arise from boundary value problems (BVPs) involving 4th-order derivatives [1,2], and, computational algorithms for the determinants are linked to the problem of obtaining efficient test for the existence of unique solutions of the PDEs [3,4]. Usually, the determinant of a square matrix can be computed by the well-known Laplace formula [5]. However, algorithms based on such formula for
∗
Corresponding author. E-mail address: [email protected] (J. Jia).
http://dx.doi.org/10.1016/j.camwa.2016.11.031 0898-1221/© 2016 Elsevier Ltd. All rights reserved.
2
J. Jia, S. Li / Computers and Mathematics with Applications (
)
–
large matrices may be time-consuming due to the number of required operations grows quickly. Therefore, more involved techniques have been developed for the determinant evaluation. In order to evaluate the determinant of the cyclic pentadiagonal matrices, some authors have devised numerical (or symbolic) algorithms recently, see [6–8]. The main objective of this paper is to develop an efficient numerical algorithm for the determinant of an n-by-n cyclic pentadiagonal Toeplitz matrix. Recent developments of algorithms for the determinants and permanents of other related matrices, see e.g. [9–16]. The present paper is organized as follows. The main results are given in the next section. In Section 3, we give the results of some numerical experiments to show the performance of our algorithm. Finally we make some concluding remarks in Section 4. 2. Main results In this section, we first present an approach for the determinant of a cyclic pentadiagonal Toeplitz matrix. Based on some elementary optimization methods, we then improve the proposed approach and develop an efficient numerical algorithm with the cost of 7n + O(log n) for evaluating nth order cyclic pentadiagonal Toeplitz determinants. 2.1. An approach based on block upper triangular transformation Let {C1 , C2 , . . . , Cn } be the column vectors of the matrix A as in (1.1). If we successively interchange the first two columns C1 and C2 with other columns to make them as the last two columns, we then obtain a nearly lower triangular Toeplitz matrix T with the column vectors {C3 , C4 , . . . , Cn , C1 , C2 }, thus
c a
d T = b e
e
b e
c
.. .. .. ..
..
.
..
.
..
.
..
.
.
..
.
..
.
..
.
e
. . .
.. ..
b
d b
a d
e
b
e ∈ Rn×n .
. .
c a
d
(2.1)
c
Since we made 2n − 2 column interchanges in the matrix reordering above, we readily obtain the following relation: det(A) = det(T ). We now consider transforming the matrix T (2.1) into a 2-by-2 block upper triangular matrix of the form:
U =
c · In − 4 0
U12 U22
∈ Rn×n ,
(2.2)
where U12 and U22 are matrices of size (n − 4)-by-4 and 4-by-4, respectively. Below, we show that the above matrix can be obtained from the original matrix T and a suitable choice ωi+1,i , ωi+2,i , ωi+3,i , and ωi+4,i of the following matrices:
1 ω21 ω31 Ω1 = ω41 ω51
I
1
,
1 1 1
..
.
i −1
1
Ωi =
ωi+1,i ωi+2,i ωi+3,i ωi+4,i
,
1 1 1 1 In−i−4
1
for i = 2, 3, . . . , n − 4. Here, Ii−1 and In−i−4 are identity matrices of order i − 1 and n − i − 4. For the sake of illustration, we consider the case n = 7 to show the transformation. From (2.1) and setting ω21 = −a/c, ω31 = −d/c, ω41 = −b/c, and ω51 = −e/c, we have
1 −a/c −d/c Ω1 T = −b/c −e/c 0 0
0 1 0 0 0 0 0
0 0 1 0 0 0 0
0 0 0 1 0 0 0
0 0 0 0 1 0 0
0 0 0 0 0 1 0
0 c 0 a 0 d 0 b 0 e 0 0 1 0
0 c a d b e 0
0 0 c a d b e
e 0 0 c a d b
b e 0 0 c a d
d b e 0 0 c a
a d b e 0 0 c
J. Jia, S. Li / Computers and Mathematics with Applications (
c 0 0 = 0 0 0 0
0 c a d b e 0
0 0 c a d b e
e t˜24 t˜34 t˜44 t˜54 d b
b t˜25 t˜35 t˜45 t˜55 a d
d t˜26 t˜36 t˜46 t˜56 c a
)
–
3
a t˜27 t˜37 t˜47 . t˜57 0 c
As we can see from the above, multiplication on the left by the matrix Ω1 only changes a few values in the first column and the last four columns of the matrix T (i.e., ti1 → 0, ti,n−3 → t˜i,n−3 , ti,n−2 → t˜i,n−2 , ti,n−1 → t˜i,n−1 , tin → t˜in , i = 2, 3, 4, 5). Similarly, choosing ω32 = ω43 = −a/c, ω42 = ω53 = −d/c, ω52 = ω63 = −b/c, ω62 = ω73 = −e/c, we finally have the following 7-by-7 matrix c 0 0 Ω3 Ω2 Ω1 T = 0 0 0 0
0 c 0 0 0 0 0
0 0 c 0 0 0 0
e t˜24 t˜34 t˜44 t˜54 t˜64 t˜74
b t˜25 t˜35 t˜45 t˜55 t˜65 t˜75
d t˜26 t˜36 t˜46 t˜56 t˜66 t˜76
a t˜27 t˜37 t˜47 . t˜57 t˜67 t˜77
Generally, the values ωi+1,i , ωi+2,i , ωi+3,i , ωi+4,i are defined as
ωi+1,i = −a/c ,
ωi+2,i = −d/c ,
ωi+3,i = −b/c ,
ωi+4,i = −e/c ,
i = 1, 2, . . . , n − 4.
Denoting Ω = Ωn−4 Ωn−5 · · · Ω1 , it follows that
ΩT = U, where U is the 2-by-2 block upper triangular matrix defined in (2.2). Since we can readily see from the definition of Ωi that det(Ω ) = 1, we directly deduce the following relation: det(A) = det(T ) = det(Ω T ) = det(U ) = c n−4 · det(U22 ). Hence, the determinant of A can be implicitly evaluated by computing the determinant of the 4-by-4 matrix U22 . Now, we are ready for describing the summary of our approach. The summary is given below. Step 1. Reorder the columns in the original cyclic pentadiagonal Toeplitz matrix to form a nearly lower triangular matrix with Toeplitz structure. Step 2. Transform the nearly lower triangular matrix into a 2-by-2 block upper triangular matrix. Step 3. Compute the determinant of the block upper triangular matrix. If we measure the computational cost of the approach in terms of total number of flops, where each flop represents one of the basic arithmetic floating point operations, then, the above approach may give det(A) in 29n + O(1) operations, since costs for the steps 2 and 3 are 28n + O(1) and n + O(1), respectively. It should be mentioned that the costs of our approach are more expensive than those of some existing methods [7,8]. In the next section, we propose a more efficient numerical algorithm with the cost of only 7n + O(log n) for evaluating nth order cyclic pentadiagonal Toeplitz determinants. 2.2. An efficient algorithm for the determinant evaluation From the previous subsection, we note that division by constant c appears multiple times, as it is used in transforming the nearly lower triangular matrix T into a block upper triangular matrix. In order to reduce the number of such operations, we set a′ = a/c ,
b′ = b/c ,
c ′ = 1,
d′ = d/c ,
e′ = e/c .
If we use these new constants in the definition of the matrix A (1.1), we obtain another cyclic pentadiagonal Toeplitz matrix A′ . Obviously, we have det(A) = c n · det(A′ ).
(2.3)
For the matrix A′ , we may follow the similar procedure applied to A so that we obtain a nearly lower triangular Toeplitz matrix T ′ corresponding to the matrix T as in (2.1). In fact, T ′ is nothing but T with new entries a′ , b′ , c ′ , d′ and e′ . Also, we can form the matrix Ωi′ corresponding to the matrix Ωi . Denoting Ω ′ = Ωn′ −4 Ωn′ −5 · · · Ω1′ , it yields
Ω ′T ′ =
In−4 0
′ U12 ′ U22
∈ R n×n ,
(2.4)
4
J. Jia, S. Li / Computers and Mathematics with Applications (
)
–
′ ′ where U12 and U22 are matrices of size (n − 4)-by-4 and 4-by-4, respectively. In addition, based on the structure of Ωi′ and simple matrix multiplication, we can readily obtain
1
(n−4)
w1 (n−4) w 2 (n−4) w3 ′ Ω = . . . . .. (n−4) wn−2
1
w1(n−5) w2(n−5) .. . .. . wn(n−−35) wn(n−−25)
wn(n−−14)
(1)
where w1
(k)
wi
(1)
..
.
..
.
1
..
.
w1(1)
..
.
w2(1)
···
(1)
′
1 1
w3
1
(1)
w4
··· (1)
= −a , w2 = −d , w3 (k−1) , wi ( k − 1 ) − a′ wk(k−)1 , w k = wk(k+−11) − d′ wk(k−)1 , wk(k+−21) − b′ wk(k−)1 , ′ (k) −e wk−1 , ′
,
(2.5)
1 (1)
= −b , w4 = −e , and ′
′
for i = 1, 2, . . . , k − 1, for i = k, for i = k + 1, for i = k + 2, for i = k + 3,
for k = 2, 3, . . . , n − 4. Finally, it follows from (2.3)–(2.5) that we have ′ det(A) = c n · det(A′ ) = c n · det(Ω ′ T ′ ) = c n · det(U22 ),
where
′ U22
(n−4) wn−4 (n−4) wn−3 = (n−4) wn−2
wn(n−−14)
wn(n−−55)
wn(n−−66)
wn(n−−45)
wn(n−−56)
wn(n−−35)
wn(n−−46)
wn(n−−25)
wn(n−−36)
wn(n−−77) ′ e wn(n−−67) 0 wn(n−−57) 0 wn(n−−47)
0
b′ e′ 0 0
d′ b′ e′ 0
1 a′ d ′ a′ + d′ b′ ′ b′ e
0 1 a′ d′
0 0 1 a′
0 0 . 0 1
Now, we state a numerical algorithm for evaluating nth order cyclic pentadiagonal Toeplitz determinants. Algorithm 2.1 Step 1. Input the elements a, b, c, d, e and the order n. (1) (1) (1) (1) Step 2. Set w1 = −a/c, w2 = −d/c, w3 = −b/c, w4 = −e/c. For i = 2, 3, . . . , n − 4 compute
w1(i) = w2(i−1) + w1(1) ∗ w1(i−1) , w2(i) = w3(i−1) + w2(1) ∗ w1(i−1) , w3(i) = w4(i−1) + w3(1) ∗ w1(i−1) , w4(i) = w4(1) ∗ w1(i−1) , End. Step 3. Form the following three 4-by-4 matrices
w1(n−4) w1(n−5) w1(n−6) w1(n−7) ( n − 4 ) ( n − 5 ) ( n − 6 ) ( n − 7 ) w w2 w2 w2 2 W1 = w(n−4) w (n−5) w(n−6) w(n−7) , 3 3 3 3 w4(n−4) w4(n−5) w4(n−6) w4(n−7) −w4(1) −w3(1) −w2(1) −w1(1) ( 1 ) ( 1 ) ( 1 ) 0 −w4 −w3 −w2 , W2 = (1) 0 0 −w4 −w3(1) 0 0 0 −w4(1)
J. Jia, S. Li / Computers and Mathematics with Applications (
1
0 1
−w1(1) W3 = −w(1) 2 −w3(1)
0 0 1
−w1(1) −w2(1)
−w1(1)
)
–
5
0 0 , 0 1
then compute B = W1 ∗ W2 + W3 . Step 4. Output the determinant: det(A) = c n · det(B). We can see from the above algorithm that the total computational cost for calculating the determinant of an n-by-n cyclic pentadiagonal Toeplitz matrix is 8n + O(1), since costs for steps 2, 3, and 4 are 7n + O(1), O(1), and n + O(1), respectively. In fact, the computational cost of Algorithm 2.1 can be reduced to 7n + O(log n) since the computation of c n (in Step 4) can be performed in O(log n) operations, for details, see [17, p. 226]. On the other hand, it should be mentioned that the algorithms given in [7,8] require 15n + O(1) and 8n + O(1), respectively. Remark 2.1. Although we work with real cyclic pentadiagonal Toeplitz matrices in this paper, the results and the proposed algorithm are also valid for cyclic pentadiagonal Toeplitz matrices with complex entries (i.e., a, b, c, d, e ∈ C).
3. Numerical examples In this section, we give the results of some simple numerical experiments to confirm that the proposed algorithm works as well as other already existing algorithms. All tests were performed in MATLAB R2014b using double precision arithmetic. Example 3.1. In this example, we evaluate the determinant of an n-by-n symmetric cyclic pentadiagonal Toeplitz matrix originating from [7] 2.5 1.5
1.5 2.5
1 1.5
.
.
1 A=
..
..
1 1.5
.
..
.. ..
1 1
.. ..
.
..
.
. . .
1 1
.. .. ..
. . .
1.5 1
.. ..
. .
2.5 1.5
1.5 1
. 1 1.5 2.5
The results with DETQPT algorithm [7], DCPT algorithm [8], and Algorithm 2.1 are shown in Table 1. From the results shown in Table 1, we see that our algorithm generated almost the same values as those given by other algorithms. Example 3.2. Now, we consider an ill-conditioned nonsymmetric cyclic pentadiagonal Toeplitz matrix that corresponds to a slight modification of Example 3.3 in [8]
1.3 + h2 A=
1.2 0.4
1 1.5
1.5 1.3 + h2
..
..
. .
1 1.5
..
.. ..
. . .
0.4 1
.. .. ..
. . .
0.4 1
.. .. ..
. . .
1.2 0.4
.. ..
1.2 0.4
. .
1.3 + h2 1.2
1 1.5 1.3 + h2
.
Here, we set h = 1/n. For n = 500, 1000, 2000, 3000, 6000, we also use DETQPT algorithm, DCPT algorithm and Algorithm 2.1 to compute the determinants. The absolute differences between the determinants derived by these three algorithms and those of MATLAB built-in function (i.e., AD1 = |DDETQPT − DMATLAB |, AD2 = |DDCPT − DMATLAB |, AD3 = |DAlgorithm 2.1 − DMATLAB |) are shown in Table 2. The condition numbers of the matrix A and the CPU times in the computation of the determinants are also given.
6
J. Jia, S. Li / Computers and Mathematics with Applications (
)
–
Table 1 Numerical results of the determinant for Example 3.1. n
DETQPT algorithm
DCPT algorithm
Algorithm 2.1
100 500 1 000 5 000 8 000 10 000 20 000
0.054821825453014 1.321092500749511 4.702608204485831 0.680884214220257 7.031833797460169 2.569002419156596 8.076085200081945
0.054821825453016 1.321092500749579 4.702608204485895 0.680884214220271 7.031833797460265 2.569002419156551 8.076085200081494
0.054821825453014 1.321092500749508 4.702608204485832 0.680884214220254 7.031833797460159 2.569002419156595 8.076085200081929
Table 2 Numerical results of the determinant for Example 3.2. Algorithms
n Cond(A)
500 1.4126e+019
1000 2.3253e+034
2000 7.8222e+082
3000 4.2945e+131
6000 6.7214e+276
DETQPT
AD1 CPU times
1.7341e−014 0.0146
1.7603e−014 0.0261
4.8559e−014 0.0398
1.7726e−013 0.0712
2.1512e−013 0.1433
DCPT
AD2 CPU times
1.7446e−014 0.0009
1.6206e−014 0.0017
4.7361e−014 0.0026
1.7615e−013 0.0051
2.1192e−013 0.0069
Algorithm 2.1
AD3 CPU times
1.7210e−014 0.0006
1.7269e−014 0.0017
4.8068e−014 0.0028
1.7590e−013 0.0049
2.1217e−013 0.0072
4. Concluding remarks In this paper, first we have presented an approach to computing the determinant of the cyclic pentadiagonal Toeplitz matrices based on block upper triangular transformation. However, the price we should pay is that the costs of the approach are more expensive than those of some existing methods. Second, we have described a numerical algorithm (Algorithm 2.1) for the determinant evaluation. The proposed algorithm is based on the use of a certain type of matrix reordering and equalities involving the products of special kind of relevant matrices, and it leads to a significant decrease in the number of operations. Moreover, we have shown that the complexity of our algorithm is less than or almost equal to those of some recent algorithms. Acknowledgments The authors would like to thank anonymous referees for useful comments that enhanced the quality of the manuscript. This work was supported by the Natural Science Foundation of China (NSFC) under grant 11601408 and the Fundamental Research Funds for the Central Universities under grant JB160701, JBG160704. References [1] S. Barnett, Matrices: Methods and Applications, Oxford University Press, New York, 1990. [2] J. Monterde, H. Ugail, A general 4th-order PDE method to generate Bézier surfaces from the boundary, Comput. Aided Geom. Design 23 (2006) 208–225. [3] R.B. Marr, G.H. Vineyard, Five-diagonal Toeplitz determinants and their relation to Chebyshev polynomials, SIAM J. Matrix Anal. Appl. 9 (4) (1988) 579–586. [4] M.E. Kanal, Parallel algorithm on inversion for adjacent pentadiagonal matrices with MPI, J. Supercomput. 59 (2012) 1071–1078. [5] G.H. Golub, C.F. Van Loan, Matrix Computations, third ed., The Johns Hopkins University Press, Baltimore and London, 1996. [6] J.T. Jia, T. Sogabe, A novel algorithm for solving quasi penta-diagonal linear systems, J. Math. Chem. 51 (2013) 881–889. [7] Y.L. Jiang, J.T. Jia, On the determinant evaluation of quasi penta-diagonal matrices and quasi penta-diagonal Toeplitz matrices, J. Math. Chem. 51 (2013) 2503–2513. [8] J.T. Jia, S.M. Li, Numerical algorithm for the determinant evaluation of cyclic pentadiagonal matrices with Toeplitz structure, Numer. Algorithms 71 (2016) 337–348. [9] R.A. Sweet, A recursive relation for the determinant of a pentadiagonal matrix, Comm. ACM. 12 (1969) 330–332. [10] T. Sogabe, A fast numerical algorithm for the determinant of a pentadiagonal matrix, Appl. Math. Comput. 196 (2008) 835–841. [11] T. Sogabe, A note on A fast numerical algorithm for the determinant of a pentadiagonal matrix, Appl. Math. Comput. 201 (2008) 561–564. [12] Z. Cinkir, An elementary algorithm for computing the determinant of pentadiagonal Toeplitz matrices, J. Comput. Appl. Math. 236 (2012) 2298–2305. [13] J.T. Jia, S.M. Li, On the inverse and determinant of general bordered tridiagonal matrices, Comput. Math. Appl. 69 (2015) 503–509. [14] J.T. Jia, T. Sogabe, S.M. Li, A generalized symbolic Thomas algorithm for the solution of opposite-bordered tridiagonal linear systems, J. Comput. Appl. Math. 290 (2015) 423–432. [15] J.T. Jia, B.T. Yang, S.M. Li, On a homogeneous recurrence relation for the determinants of general pentadiagonal Toeplitz matrices, Comput. Math. Appl. 71 (2016) 1036–1044. [16] M.E.A. El-Mikkawy, A fast and reliable algorithm for evaluating nth order pentadiagonal determinants, Appl. Math. Comput. 202 (2008) 210–215. [17] K.H. Rosen, Discrete Mathematics and its Applications, McGraw-Hill, New York, 2007.