Decompositions of Some Special Block Tridiagonal Matrices

Hsin-Chu Chen

Department of Cyber-Physical Systems, Clark Atlanta University, Atlanta, USA.

DOI: 10.4236/alamt.2021.112005                

Abstract

In this paper, we present a unified approach to decomposing a special class of block tridiagonal matrices K (α ,β ) into block diagonal matrices using similarity transformations. The matrices K (α ,β )∈ R^{pq× pq} are of the form K (α ,β = block-tridiag[β B,A,α B] for three special pairs of (α ,β ): K (1,1), K (1,2) and K (2,2) , where the matrices A and B, A, B∈ R^{p× q} , are general square matrices. The decomposed block diagonal matrices (α ,β ) for the three cases are all of the form: (α ,β ) = D₁ (α ,β ) ⊕ D₂ (α ,β ) ⊕---⊕ D_q (α ,β ) , where D_k (α ,β ) = A+ 2cos ( θ_k (α ,β )) B, in which θ_k (α ,β ) , k = 1,2, --- q , depend on the values of α and β. Our decomposition method is closely related to the classical fast Poisson solver using Fourier analysis. Unlike the fast Poisson solver, our approach decomposes K (α ,β ) into q diagonal blocks, instead of p blocks. Furthermore, our proposed approach does not require matrices A and B to be symmetric and commute, and employs only the eigenvectors of the tridiagonal matrix T (α ,β ) = tridiag[β b, a,αb] in a block form, where a and b are scalars. The transformation matrices, their inverses, and the explicit form of the decomposed block diagonal matrices are derived in this paper. Numerical examples and experiments are also presented to demonstrate the validity and usefulness of the approach. Due to the decoupled nature of the decomposed matrices, this approach lends itself to parallel and distributed computations for solving both linear systems and eigenvalue problems using multiprocessors.

Keywords

Block Tridiagonal Matrices, Block Fourier Decomposition, Linear Systems, Eigenvalue Problems

1. Introduction

In this paper, we present explicit similarity transformations to decompose block tridiagonal matrices $K\left(\alpha \mathrm{,}\beta \right)\in {\mathcal{R}}^{pq\times pq}$ of the following form:

$K\left(\alpha ,\beta \right)=\left[\begin{array}{ccccc}A& \alpha B& & & \\ B& A& B& & \\ & \ddots & \ddots & \ddots & \\ & & B& A& B\\ & & & \beta B& A\end{array}\right]$ (1)

for some special pairs of $\left(\alpha ,\beta \right)$, where $A\mathrm{,}B\in {\mathcal{R}}^{p\times p}$, into block diogonal matrices. The $\left(\alpha ,\beta \right)$ pairs to be considered in this paper include (1,1), (1,2), and (2,2). We shall show that the transformations $K\left(\alpha \mathrm{,}\beta \right)$ for these three $\left(\alpha ,\beta \right)$ pairs all lead to the block diagonal matrices $\tilde{K}\left(\alpha \mathrm{,}\beta \right)$ of the following single unified form:

$\tilde{K}\left(\alpha \mathrm{,}\beta \right)=D_1\left(\alpha \mathrm{,}\beta \right)\oplus D_2\left(\alpha \mathrm{,}\beta \right)\oplus \cdots \oplus D_q\left(\alpha \mathrm{,}\beta \right)$ (2)

where the operation symbol $\oplus$ denotes the matrix direct sum and the diagonal submatrices $D_k\left(\alpha \mathrm{,}\beta \right)=A+2\cos \left({\theta }_k\left(\alpha \mathrm{,}\beta \right)\right)B$ in which ${\theta }_k\left(\alpha \mathrm{,}\beta \right)$, $k=1,2,\cdots ,q$, are explicitly known, although they depend on the values of $\alpha$ and $\beta$. Our decomposition method is closely related to the classical fast Poisson solver [1] [2] using Fourier analysis.

The block decomposition scheme to be addressed has been presented by the author in [3] and formal proof was given for $K\left(1,1\right)$. The decompositions for the other two cases were simply mentioned without derivations. Unfortunately, that paper consists of two errors. First, the eigenvectors used to form the transformation matrix for decomposing $K\left(\mathrm{2,2}\right)$ and the decomposed submatrices are incorrect, which will be addressed in Theorem 2 of Section 2 in this paper. Second, the transformation matrix Q for $K\left(1,2\right)$ is not orthogonal, although the eigenvectors and the decomposed submatrices presented are correct, meaning that the expression $Q^{\text{T}}KQ$ should read $Q^{-1}KQ$. The main purposes of this paper include the following tasks.

1) We show that the transformation matrix for decomposing $K\left(1,1\right)$ is orthogonal in Theorem 1.

2) We take this opportunity to correct the errors made in our previous paper by providing formal proof with the transformation matrix formed by the correct eigenvectors for decomposing $K\left(\mathrm{2,2}\right)$ in Theorem 2.

3) We also present a formal derivation for the decomposition of $K\left(1,2\right)$ into a block diagonal matrix in Theorem 3.

4) Numerical examples and experiments for all three cases will be given in Section 3 to demonstrate the validity and the advantages of the decompositions.

The block decompositions are all based on similarity transformations with known eigenvectors of certain tridiagonal matrices and they all yield a single unified form of block diagonal matrices. Since similarity transformations preserve all eigenvalues, the eigenvalues of the original matrix $K\left(\alpha \mathrm{,}\beta \right)$, which is of size pq by pq, can be obtained from the q diagonal blocks $D_k\left(\alpha \mathrm{,}\beta \right)$, each of size only p by p. This block decomposition scheme provides a much more efficient means for solving eigenvalue problems with this type of coefficient matrices. It can also be employed for solving linear systems with efficiency because the transformation matrices are explicitly known. In addition, the decoupled structure of the transformed matrix lends itself to parallel computation with coarse-grain parallelism.

2. Decompositions

In this following, we present our key observations that lead to the proposed block decomposition method for this class of matrices $K\left(\alpha \mathrm{,}\beta \right)$ using transformation matrices whose entries are inherent in the special block tridiagonal form of $K\left(\alpha \mathrm{,}\beta \right)$. Whenever there is no confusion, we shall simply use K to denote $K\left(\alpha \mathrm{,}\beta \right)$. Throughout the paper, the operation symbols $\oplus$ and $\otimes$ are used to denote the matrix direct sum and the Kronecker product.

Theorem 1. When $\alpha =\beta =1$ , the block tridiagonal matrix $K\left(\alpha \mathrm{,}\beta \right)$ is orthogonally similar to the block diagonal matrix $D=D_1\oplus D_2\oplus \cdots \oplus D_q$ ,

$D_k=A+2\cos \left({\theta }_k\right)B$, where ${\theta }_k=\left(\frac{k}{q+1}\right)\pi$ , $k=1,2,\cdots ,q$.

Proof. Let $v_k^{\text{T}}=\sqrt{\frac{2}{q+1}}\left[\sin \left({\theta }_k\right),\sin \left(2{\theta }_k\right),\cdots ,\sin \left(q{\theta }_k\right)\right]$ and $V_k=v_k\otimes I_p$,

$k=1,2,\cdots ,q$, where $I_p$ is the identity matrix of dimension p and the symbol $\otimes$ denotes the Kronecker product. It has been shown in [3] that

$V_i^{\text{T}}K{V}_j=\left(v_i^{\text{T}}{v}_j\right)\left[A+2\cos \left({\theta }_j\right)B\right]=\{\begin{array}{ll}A+2\cos \left({\theta }_i\right)B\hfill & \text{for}i=j\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$

by stating that $v_i^{\text{T}}{v}_j$ is orthonormal. This paper will skip the proof and just provide the details to show that $v_i^{\text{T}}{v}_j$ is indeed orthonormal. To this end, we need the following formula [4]:

${\displaystyle \underset{k=1}{\overset{q}{\sum }}}\cos \left(kx\right)=\frac{1}{2}\left[S\left(x\right)-1\right],$

$S\left(x\right)=\left[\frac{\sin \left(\left(2q+1\right)\left(x/2\right)\right)}{\sin \left(x/2\right)}\right].$

Let $t_{ij}$ and $u_{ij}$ denote $\left({\theta }_i-{\theta }_j\right)$ and $\left({\theta }_i+{\theta }_j\right)$, respectively. We have

$\begin{array}{c}{v}_i^{\text{T}}{v}_j=\frac{2}{q+1}{\displaystyle \underset{k=1}{\overset{q}{\sum }}}\sin \left(k{\theta }_i\right)\sin \left(k{\theta }_j\right)=\frac{1}{q+1}{\displaystyle \underset{k=1}{\overset{q}{\sum }}}\left[\cos k\left({\theta }_i-{\theta }_j\right)-\cos k\left({\theta }_i+{\theta }_j\right)\right]\\ =\frac{1}{q+1}{\displaystyle \underset{k=1}{\overset{q}{\sum }}}\left(\cos k{t}_{ij}-\cos k{u}_{ij}\right)=\frac{1}{2\left(q+1\right)}\left[S\left(t_{ij}\right)-S\left(u_{ij}\right)\right].\end{array}$

Note that, by L’Hospital’s rule,

$S\left(x\right)=\left(2q+1\right)\left[\frac{\cos \left(\left(2q+1\right)\left(x/2\right)\right)}{\cos \left(x/2\right)}\right]$

if $\sin \left(x/2\right)=0$. This will be the case when $x=2k\pi$ for any integer k. Now since $-\left(q-1\right)\le \left(i-j\right)\le \left(q-1\right)$, we have $-\pi <t_{ij}={\theta }_i-{\theta }_j=\frac{i-j}{q+1}\pi <\pi$. The

denominator of $S\left(t_{ij}\right)$ will be equal to zero only if $i=j$. When $i=j$, we have

$S\left(t_{ij}\right)=S\left(0\right)=2q+1.$

When $i\ne j$, $S\left(t_{ij}\right)$ can be simplified as

$\begin{array}{c}S\left(t_{ij}\right)=\frac{\sin \left(\left(2q+1\right)\left(t_{ij}/2\right)\right)}{\sin \left(t_{ij}/2\right)}\\ =\frac{\sin \left(2\left(q+1\right)\left(t_{ij}/2\right)\right)\cos \left(t_{ij}/2\right)-\cos \left(2\left(q+1\right)\left(t_{ij}/2\right)\right)\sin \left(t_{ij}/2\right)}{\sin \left(t_{ij}/2\right)}\\ =\frac{\sin \left(\left(i-j\right)\pi \right)\cos \left(t_{ij}/2\right)-\cos \left(\left(i-j\right)\pi \right)\sin \left(t_{ij}/2\right)}{\sin \left(t_{ij}/2\right)}\\ =-\cos \left(\left(i-j\right)\pi \right)={\left(-1\right)}^{i-j+1}\end{array}$

since $\sin \left(i-j\right)\pi =0$, where we have uesd the fact that $\sin \left(x-y\right)=\sin x\cos y-\sin y\cos x$. Therefore,

$S\left(t_{ij}\right)=\{\begin{array}{ll}2q+1\hfill & \text{for}i=j\hfill \\ -{\left(-1\right)}^{i-j}\hfill & \text{otherwise}\hfill \end{array}$ (3)

Likewise, since $0<\left(i+j\right)\le 2q$, we have

$0<u_{ij}={\theta }_i+{\theta }_j=\frac{i+j-2}{q+1}\pi <2\pi .$

The denominator of $S\left(u_{ij}\right)$ will never be equal to zero. Accordingly,

$S\left(u_{ij}\right)=-\cos \left(\left(i+j\right)\pi \right)={\left(-1\right)}^{i+j+1}.$ (4)

Finally from (3) and (4), we obtain

$S\left(t_{ij}\right)-S\left(u_{ij}\right)=\{\begin{array}{ll}2\left(q+1\right)\hfill & \text{for}i=j\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$

and, therefore,

$v_i^{\text{T}}{v}_j=\frac{1}{2\left(q+1\right)}\left[S\left(t_{ij}\right)-S\left(u_{ij}\right)\right]=\{\begin{array}{cc}1& \text{for}i=j\\ 0& \text{otherwise}\end{array}$

indicating the vectors $v_i$ are orthonormal. Accordingly, we have

$V_i^{\text{T}}K{V}_j=\{\begin{array}{ll}A+2\cos \left({\theta }_j\right)B\hfill & \text{for}i=j\hfill \\ 0\hfill & \text{otherwise}\hfill \end{array}$

This completes the proof.

Theorem 2. When $\alpha =\beta =2$ , the block tridiagonal matrix $K\left(\alpha \mathrm{,}\beta \right)$ is similar to the block diagonal matrix $D=D_1\oplus D_2\oplus \cdots \oplus D_q$ , $D_k=A+2\cos \left({\theta }_k\right)B$ ,

where ${\theta }_k=\left(\frac{k-1}{q-1}\right)\pi$ , $k=1,2,\cdots ,q$ .

Proof. This block diagonalization was mentioned previously in Corollary 2 in [3] without a proof. Unfortunately, the eigenvectors $v_k^{\text{T}}$ used to form the transformation matrix Q and the decomposed submatrices D_k consist of errors,

in which 1) the vector $\sqrt{\frac{2}{q}}\left[\cos \left({\theta }_k\right),\cos \left(2{\theta }_k\right),\cdots ,\cos \left(q{\theta }_k\right)\right]$ as stated in that paper should be replaced by $\sqrt{\frac{2}{q-1}}\left[1,\cos \left({\theta }_k\right),\cos \left(2{\theta }_k\right),\cdots ,\cos \left(\left(q-1\right){\theta }_k\right)\right]$ with ${\theta }_k=\left(\frac{k-1}{q-1}\right)\pi$ ; and 2) the expression $\sin \left({\theta }_k\right)$ in D_k should read $\cos \left({\theta }_k\right)$.

In this paper, we give a formal proof with the correct eigenvectors and provide the explicit form of the inverse of the transformation matrix Q for $K\left(\mathrm{2,2}\right)$.

Let $v_k^{\text{T}}=\sqrt{\frac{2}{q-1}}\left[1,\cos \left({\theta }_k\right),\cos \left(2{\theta }_k\right),\cdots ,\cos \left(\left(q-1\right){\theta }_k\right)\right]$. Let also

$V_k=v_k\otimes I_p$, $k=1,2,\cdots ,q$ ; and $Q=\left[V_1\mathrm{,}{V}_2\mathrm{,}\cdots \mathrm{,}{V}_q\right]$. We now show that the similarity transformation $Q^{-1}KQ$ block-diagonalizes K into D. It deserves mentioning that Q in this case is not orthogonal. $Q^{-1}$, however, exists and is explicitly known. Therefore, It suffices to show that $K{V}_j=V_j{D}_j$ for all $j=1,2,\cdots ,q$.

$\begin{array}{c}K{V}_j=\sqrt{\frac{2}{q-1}}\left[\begin{array}{ccccc}A& 2B& & & \\ B& A& B& & \\ & \ddots & \ddots & \ddots & \\ & & B& A& B\\ & & & 2B& A\end{array}\right]\left[\begin{array}{c}I\\ \cos \left({\theta }_j\right)I\\ ⋮\\ \cos \left(\left(q-2\right){\theta }_j\right)I\\ \cos \left(\left(q-1\right){\theta }_j\right)I\end{array}\right]\\ =\sqrt{\frac{2}{q-1}}\left[\begin{array}{c}A+2\cos \left({\theta }_j\right)B\\ \cos \left({\theta }_j\right)A+\left[1+\cos \left(2{\theta }_j\right)\right]B\\ \cos \left(2{\theta }_j\right)A+\left[\cos \left({\theta }_j\right)+\cos \left(3{\theta }_j\right)\right]B\\ ⋮\\ \cos \left(\left(q-2\right){\theta }_j\right)A+\left[\cos \left(\left(q-3\right){\theta }_j\right)+\cos \left(\left(q-1\right){\theta }_j\right)\right]B\\ \cos \left(\left(q-1\right){\theta }_j\right)A+2\cos \left(\left(q-2\right){\theta }_j\right)B\end{array}\right]\end{array}$

Applying the identities:

$1+\cos \left(2x\right)=2{\cos }^2\left(x\right)\text{and}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\cos \left(x+y\right)+\cos \left(x-y\right)=2\cos \left(x\right)\cos (\; y\; )$

and noting that:

$\begin{array}{c}\cos \left(\left(q-2\right){\theta }_j\right)=\cos \left(\left(q-1\right){\theta }_j\right)\cos \left({\theta }_j\right)+\sin \left(\left(q-1\right){\theta }_j\right)\sin \left({\theta }_j\right)\\ =\cos \left(\left(q-1\right){\theta }_j\right)\cos (\; \theta \; j\; )\end{array}$

where we have used the fact that $\sin \left(\left(q-1\right){\theta }_j\right)=\sin \left(\left(j-1\right)\pi \right)=0$, we obtain

$K{V}_j=\sqrt{\frac{2}{q-1}}\left[\begin{array}{c}I\\ \cos \left({\theta }_j\right)I\\ \cos \left(2{\theta }_j\right)I\\ ⋮\\ \cos \left(\left(q-2\right){\theta }_j\right)I\\ \cos \left(\left(q-1\right){\theta }_j\right)I\end{array}\right]\left(A+2B\cos \left({\theta }_j\right)\right)=V_j{D}_j.$ (5)

Equation (5) holds for all j, $1\le j\le q$. Accordingly, by arranging all $V_j$ together to form the matrix Q, we obtain $KQ=QD$. In other words, the transformation $Q^{-1}KQ$ block-diagonalizes the matrix $K\left(\mathrm{2,2}\right)$ :

$Q^{-1}KQ=D$ (6)

where $Q=\left[V_1,V_2,\cdots ,V_q\right]$ and $D=D_1\oplus D_2\oplus \cdots \oplus D_q$.

This is a similarity transformation and, therefore, all eigenvalues of $K\left(\mathrm{2,2}\right)$ are preserved in the decomposed matrix D. It is worth mentioning that obtaining all the eigenvalues from D is far more efficient than from the original matrix K since D consists of only q diagonal blocks: $D_j,j=1,\cdots ,q$. When it comes to solving linear systems in the transformed space that involves Q, however, one needs to employ the LU decomposition of Q or to find $Q_{-1}$. Normally, finding the LU decomposition is more efficient and preferred. However, it does not make sense to find the LU decomposition of Q if the inverse of Q is readily available. In the following, we show that $Q^{-1}$ can be obtained explicitly.

Let C be the matrix formed by $v_k$ : $C=\left[v_1\mathrm{,}{v}_2\mathrm{,}\cdots \mathrm{,}{v}_q\right]$, whose explicit form is

$C=\sqrt{\frac{2}{q-1}}\left[\begin{array}{ccccc}1& 1& \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\cdots \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& 1& 1\\ 1& u\left(\mathrm{1,1}\right)& \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\cdots \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& u\left(\mathrm{1,}q-2\right)& -1\\ 1& u\left(\mathrm{2,1}\right)& \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\cdots \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& u\left(\mathrm{2,}q-2\right)& 1\\ & & \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\ddots \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& & \\ 1& u\left(q-\mathrm{2,1}\right)& \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\cdots \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& u\left(q-\mathrm{2,}q-2\right)& {\left(-1\right)}^{q-2}\\ 1& -1& \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\cdots \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& {\left(-1\right)}^{q-2}& {\left(-1\right)}^{q-1}\end{array}\right]$

where the symbol $u\left(i\mathrm{,}j\right)$ denotes $\cos \left(\frac{ij}{q-1}\pi \right)$. It is well known that

$C^{-1}=\sqrt{\frac{2}{q-1}}\left[\begin{array}{ccccc}1/4& 1/2& \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\cdots \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& 1/2& 1/4\\ 1/2& u\left(1,1\right)& \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\cdots \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& u\left(1,q-2\right)& -1/2\\ 1/2& u\left(2,1\right)& \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\cdots \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& u\left(2,q-2\right)& 1/2\\ & & \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\ddots \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& & \\ 1/2& u\left(q-2,1\right)& \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\cdots \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& u\left(q-2,q-2\right)& {\left(-1\right)}^{q-2}/2\\ 1/4& -1/2& \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\cdots \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}& {\left(-1\right)}^{q-2}/2& {\left(-1\right)}^{q-1}/4\end{array}\right]$

Let $S=\text{diag}\left(\sqrt{2}\mathrm{,1,}\cdots \mathrm{,1,}\sqrt{2}\right)$, a diagonal matrix. It can easily be seen that

$C^{-1}=S^{-2}C{S}^{-2}.$

Recall that $V_k=v_k\otimes I_p,k=1,2,\cdots ,q$ ; and $Q=\left[V_1\mathrm{,}{V}_2\mathrm{,}\cdots \mathrm{,}{V}_q\right]$, which can be expressed as:

$\begin{array}{c}Q=\left[v_1\otimes I_p,v_2\otimes I_p,\cdots ,v_q\otimes I_p\right]\\ =\left[v_1,v_2,\cdots ,v_q\right]\otimes I_p=C\otimes I_p.\end{array}$

Therefore,

$\begin{array}{c}{Q}^{-1}={\left(C\otimes I_p\right)}^{-1}=C^{-1}\otimes I_p=\left(S^{-2}C{S}^{-2}\right)\otimes I_p\\ =\left(S^{-2}\otimes I_p\right)\left(C\otimes I_p\right)\left(S^{-2}\otimes I_p\right)\\ =\left(S^{-2}\otimes I_p\right)Q\left(S^{-2}\otimes I_p\right)\end{array}$

where we have used the following two properties of Kronecker products [5]:

${\left(A\otimes B\right)}^{-1}=A^{-1}\otimes B^{-1}\text{and}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\left(AB\right)\otimes I=\left(A\otimes I\right)\left(B\otimes I\right).$

Note that $S^{-2}\otimes I_p$ is a diagonal matrix and is almost an identity matrix except the 1st and the last blocks. This shows that $Q^{-1}$ is almost the same as Q. Computationally, $Q^{-1}$ can be obtained directly from $C^{-1}$, which is explicitly known, since $Q^{-1}$ is simply a block structure of $C^{-1}$. Note that C is symmetric, but not orthogonal.

Theorem 3. When $\alpha =1$ and $\beta =2$ , the block tridiagonal matrix $K\left(\alpha \mathrm{,}\beta \right)$ is similar to the block diagonal matrix $D=D_1\oplus D_2\oplus \cdots \oplus D_q$ ,

$D_k=A+2\cos \left({\theta }_k\right)B$, where ${\theta }_k=\left(\frac{2k-1}{2q}\right)\pi$ , $k=1,2,\cdots ,q$.

Proof. Let $v_k^{\text{T}}=\sqrt{\frac{2}{q}}\left[\sin \left({\theta }_k\right),\sin \left(2{\theta }_k\right),\cdots ,\sin \left(q{\theta }_k\right)\right]$ ; $V_k=v_k\otimes I_p$,

$k=1,2,\cdots ,q$ ; and $Q=\left[V_1\mathrm{,}{V}_2\mathrm{,}\cdots \mathrm{,}{V}_q\right]$. We have

$\begin{array}{c}K{V}_j=\sqrt{\frac{2}{q}}\left[\begin{array}{ccccc}A& B& & & \\ B& A& B& & \\ & \ddots & \ddots & \ddots & \\ & & B& A& B\\ & & & 2B& A\end{array}\right]\left[\begin{array}{c}\sin \left({\theta }_j\right)I\\ \sin \left(2{\theta }_j\right)I\\ ⋮\\ \sin \left(\left(q-1\right){\theta }_j\right)I\\ \sin \left(q{\theta }_j\right)I\end{array}\right]\\ =\sqrt{\frac{2}{q}}\left[\begin{array}{c}\sin \left({\theta }_j\right)A+\sin \left(2{\theta }_j\right)B\\ \sin \left(2{\theta }_j\right)A+\left[\sin \left({\theta }_j\right)+\sin \left(3{\theta }_j\right)\right]B\\ ⋮\\ \sin \left(\left(q-1\right){\theta }_j\right)A+\left[\sin \left(\left(q-2\right){\theta }_j\right)+\sin \left(q{\theta }_j\right)\right]B\\ \sin \left(q{\theta }_j\right)A+2\sin \left(\left(q-1\right){\theta }_j\right)B\end{array}\right]\end{array}$

Using the identities $\sin \left(2x\right)=2\sin \left(x\right)\cos \left(x\right)$ and $\sin \left(x+y\right)+\sin \left(x-y\right)=2\sin \left(x\right)\cos \left(y\right)$ yields

$K{V}_j=\sqrt{\frac{2}{q}}\left[\begin{array}{c}\sin \left({\theta }_j\right)I\\ \sin \left(2{\theta }_j\right)I\\ ⋮\\ \sin \left(\left(q-1\right){\theta }_j\right)I\\ \sin \left(q{\theta }_j\right)I\end{array}\right]\left(A+2\cos \left({\theta }_j\right)B\right)=V_j{D}_j$

where we have used the fact that

$\sin \left(\left(q-1\right){\theta }_j\right)=\sin \left(q{\theta }_j\right)\cos \left({\theta }_j\right)-\cos \left(q{\theta }_j\right)\sin \left({\theta }_j\right)=\sin \left(q{\theta }_j\right)\cos \left({\theta }_j\right)$ since

$\cos \left(q{\theta }_j\right)=\cos \left(\frac{2j-1}{2}\pi \right)=0$. Note that the matrix $Q$ here is not orthogonal

either. It can be shown that $Q^{-1}$ exists. The transformation $Q^{-1}KQ$, therefore, block-diagonalizes K into D.

In the following, we show that $Q^{-1}$ is almost identical to $Q^{\text{T}}$ and, therefore, can be explicitly obtained from $Q^{\text{T}}$ without any difficulty. Again, let C be the matrix formed by $v_k$ : $C=\left[v_1\mathrm{,}{v}_2\mathrm{,}\cdots \mathrm{,}{v}_q\right]$, $V_k=v_k\otimes I_p$, and $Q=\left[V_1\mathrm{,}{V}_2\mathrm{,}\cdots \mathrm{,}{V}_q\right]$ as was done in the previos section. We have in this case:

$C=\sqrt{\frac{2}{q}}\left[\begin{array}{cccc}s\left(1,1\right)& s\left(1,2\right)& \cdots & s\left(1,q\right)\\ s\left(2,1\right)& s\left(2,2\right)& \cdots & s\left(2,q\right)\\ & & \ddots & \\ s\left(\left(q-1\right),1\right)& s\left(\left(q-1\right),2\right)& \cdots & s\left(\left(q-1\right),q\right)\\ 1& -1& \ddots & {\left(-1\right)}^{q-1}\end{array}\right].$

whose inverse is:

$C^{-1}=\sqrt{\frac{2}{q}}\left[\begin{array}{cccccc}s\left(1,1\right)& s\left(2,1\right)& \cdots & s\left(\left(q-1\right),1\right)& |& 1/2\\ s\left(1,2\right)& s\left(2,2\right)& \cdots & s\left(\left(q-1\right),2\right)& |& -1/2\\ & & \ddots & & |& ⋮\\ s\left(1,q\right)& s\left(2,q\right)& \cdots & s\left(\left(q-1\right),q\right)& |& {\left(-1\right)}^{q-1}/2\end{array}\right]$

where the symbol $s\left(i\mathrm{,}j\right)$ denotes $\sin \left(i{\theta }_j\right)$ with ${\theta }_j=\left(\frac{2j-1}{2q}\right)\pi$. If C is partitioned as $\left[\frac{R}{r}\right]$, with R consisting of the first $\left(q-1\right)$ rows and r being the last row of C, we clearly see that $C^{-1}=\left[R^{\text{T}}\mathrm{<mo>\; |\; </mo>}\frac{1}{2}{r}^{\text{T}}\right]$. In other words, $C^{-1}$ is the

same as $C^{\text{T}}$ except the last column. Now let $S=\text{diag}\left(\mathrm{1,}\cdots \mathrm{,1,}\sqrt{2}\right)$, a diagonal matrix. We have $C^{-1}=S^{-2}C{S}^{-2}$. Following the same derivation as we have done in Theorem 2, we conclude that:

$Q^{-1}={\left(C\otimes I_p\right)}^{-1}=C^{-1}\otimes I_p=\left(S^{-2}\otimes I_p\right)Q\left(S^{-2}\otimes I_p\right)$

since $Q=C\otimes I_p$. Note that C in this case is neither symmetric nor orthogonal.

3. Numerical Experiments

To demonstrate the validity and advantage of this block decomposition approach, we present numerical experiments for all three cases of the $\left(\alpha \mathrm{,}\beta \right)$ pair discussed in this paper. Our main task in the experiments is to find all the eigenvalues of the following matrix $K\left(\alpha \mathrm{,}\beta \right)\in {\mathcal{R}}^{pq\times pq}$, with $p=4$ and $q=5$ :

$K\left(\alpha ,\beta \right)=\left[\begin{array}{ccccc}A& \alpha B& & & \\ B& A& B& & \\ & \ddots & \ddots & \ddots & \\ & & B& A& B\\ & & & \beta B& A\end{array}\right]$ (7)

where,

$A=\left[\begin{array}{rrrr}\hfill 4& \hfill -2& \hfill & \hfill \\ \hfill -2& \hfill 8& \hfill -2& \hfill \\ \hfill & \hfill -2& \hfill 8& \hfill -2\\ \hfill & \hfill & \hfill -2& \hfill 4\end{array}\right]\text{and}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}B=\left[\begin{array}{cccc}-1& & & \\ & -2& & \\ & & -2& \\ & & & -1\end{array}\right]$

The unshown entries in A and B are all zeros. The matrix $K\left(\mathrm{1,1}\right)$ can be obtained from a finite element discretization of the heat or membrane equation over a rectangular domain, subject to Dirichlet boundary condition along two opposite sides of the boundary [6]. The matrices $K\left(\mathrm{1,2}\right)$ and $K\left(\mathrm{2,2}\right)$ can be obtained from discretization of the same problem subject to certain Neumann boundary conditions. The values of p and q depend on the number of grid lines in the discretization domain.

We intentionally keep the dimension of A and B small so that one can easily see that $AB\ne BA$ by explicit multiplications. In other words, the matrices A and B do not commute and, therefore, the traditional fast Poisson solver fails to simultaneously decompose A and B. Using the block decomposition approach presented in the previous section, it is not difficult to see that $K\left(\alpha \mathrm{,}\beta \right)$ is similar to a block diagonal matrix $\tilde{K}\left(\alpha \mathrm{,}\beta \right)$ of the form:

$\tilde{K}\left(\alpha ,\beta \right)=D_1\left(\alpha ,\beta \right)\oplus D_2\left(\alpha ,\beta \right)\oplus \cdots \oplus D_q\left(\alpha ,\beta \right)$

for all three cases where,

$D_k\left(\alpha \mathrm{,}\beta \right)=\left[\begin{array}{cccc}4-2\cos \left({\theta }_k\right)& -2& & \\ -2& 8-4\cos \left({\theta }_k\right)& -2& \\ & -2& 8-4\cos \left({\theta }_k\right)& -2\\ & & -2& 4-2\cos \left({\theta }_k\right)\end{array}\right]$ (8)

with ${\theta }_k=\{\begin{array}{l}\left(\frac{k}{q+1}\right)\pi \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{for}\left(\alpha ,\beta \right)=\left(1,1\right)\hfill \\ \left(\frac{k-1}{q-1}\right)\pi \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>}\mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{for}\left(\alpha ,\beta \right)=\left(2,2\right)\hfill \\ \left(\frac{2k-1}{2q}\right)\pi \mathrm{<mtext>\; \hspace{0.17em}\; </mtext>\; <mtext>\; \hspace{0.17em}\; </mtext>}\text{for}\left(\alpha ,\beta \right)=\left(1,2\right)\hfill \end{array}$, $k=1,2,\cdots ,q$.

In the following, we present the results from our experiments, which were performed using the Octave software. Specifically, the Octave built-in function $eig\left(X\right)$ is used to obtain the eigenvalues of a square matrix X. For comparison purposes, we first display the eigenvalues computed from the original matrices $K\left(\alpha \mathrm{,}\beta \right)$ (Equation (7)) without applying the block decomposition for all three $\left(\alpha \mathrm{,}\beta \right)$ cases. Note that each $K\left(\alpha \mathrm{,}\beta \right)$ is a 20 by 20 matrix. The numerical results are listed in Table 1 where all eigenvalues are listed in the order produced by Octave without reordering.

We then present in Table 2, Table 3, and Table 4 the eigenvalues obtained directly from the decomposed diagonal blocks D₁ through D₅ (Equation (8)) of $\tilde{K}\left(\mathrm{1,1}\right)$, $\tilde{K}\left(\mathrm{1,2}\right)$, and $\tilde{K}\left(\mathrm{2,2}\right)$, respectively, where each $D_k$ is a 4 by 4 matrix. As can be seen from these tables, all eigenvalues are preserved after the block decomposition. For example, all the eigenvalues shown in Table 2 are identical to those of $K\left(\mathrm{1,1}\right)$ in Table 1, except for the ordering. This is due to the fact that all of the three block decompositions are similarity transformations. The advantage of using the decomposed matrices to compute the eigenvalues is apparent because the diagonal blocks of the decomposed matrix are explicitly known and need only A, B, and ${\theta }_k$, without the need of forming the block-tridiagonal matrix $K\left(\alpha \mathrm{,}\beta \right)$.

To close this section, it deserves mentioning that the computational complexity of finding all eigenvalues of a square matrix of size $n\times n$ is $\mathcal{O}\left(n^3\right)$ in general. For the matrices $K\left(\alpha \mathrm{,}\beta \right)$ presented in this paper, $n=pq$. Without decompositions, the computational complexity is $\mathcal{O}\left(p^3{q}^3\right)$. With the proposed

Table 1. Eigenvalues computed from the original matrices $K\left(\alpha \mathrm{,}\beta \right)$ without decomposition.

Table 2. Eigenvalues computed directly from D₁ through D₅ of $\tilde{K}\left(\mathrm{1,1}\right)$.

Table 3. Eigenvalues computed directly from D₁ through D₅ of $\tilde{K}\left(\mathrm{1,2}\right)$.

Table 4. Eigenvalues computed directly from D₁ through D₅ of $\tilde{K}\left(\mathrm{2,2}\right)$.

decomposition, the computational complexity reduces to only $\mathcal{O}\left(q{p}^3\right)$, a significant saving in computation. The advantage of the decomposition is obvious, not to mention the additional advantage that can be exploited from the coarse-grain parallelism offered by the block decomposition when the problem is to be solved using multiple processors.

4. Conclusions

In this paper, we have presented a unified block decomposition scheme for three special cases of block tridiagonal matrices of the form $K\left(\alpha \mathrm{,}\beta \right)$, as shown in Equation (1). This class of block tridiagonal matrices arises frequently from the finite difference approximation to solving certain partial differential equations such as the Laplace’s, Poisson’s, or Helmholtz equations using five- or nine-point schemes, over a rectangular or cubic domain [7]. They can also arise from some finite-element discretization of the same equation [8] and from surface fitting with the B-spline functions [9]. The values of $\alpha$ and $\beta$ typically depend on the boundary conditions of the physical problem: $\alpha =\beta =1$ for Dirichlet-Dirichlet conditions, $\alpha =1$ and $\beta =2$ for Dirichlet-Neumann conditions, and $\alpha =\beta =2$ for Neumann-Neumann conditions.

The block decompositions presented are all based on similarity transformations with known eigenvectors of tridiagonal matrices of the same form as $K\left(\alpha \mathrm{,}\beta \right)$, with the submatrices A and B replaced by scalars $a$ and $b$. We have also derived the explicit decomposed block diagonal matrices $\tilde{K}\left(\alpha \mathrm{,}\beta \right)$ for all of the three cases ( $K\left(\mathrm{1,1}\right)$, $K\left(\mathrm{1,2}\right)$, and $K\left(\mathrm{2,2}\right)$ ):

$\tilde{K}\left(\alpha ,\beta \right)=D_1\left(\alpha ,\beta \right)\oplus D_2\left(\alpha ,\beta \right)\oplus \cdots \oplus D_q\left(\alpha ,\beta \right)$

where $D_k\left(\alpha ,\beta \right)=A+2\cos \left({\theta }_k\left(\alpha ,\beta \right)\right)B$ in which ${\theta }_k\left(\alpha \mathrm{,}\beta \right)$, $k=1,2,\cdots ,q$, depend on the values of $\alpha$ and $\beta$, as can be seen in Section 2. The availability of the explicit decomposed matrices offers great computational and programming advantages. Numerical experiments have been conducted using the software Octave to demonstrate its validity and advantages. Although analogous to the classical fast Poisson solver, this approach does not require matrices A and B be symmetric and commute. This approach also exploits large-grain parallelism and lends itself to parallel and distributed computations for finding the solution of both linear systems and eigenvalue problems.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Buzbee, B.L., Golub, G.H. and Nielson, C.W. (1970) On Direct Methods for Solving Poisson’s Equations. SIAM Journal on Numerical Analysis, 7, 627-656. https://doi.org/10.1137/0707049
[2]	Strang, G. (1986) Introduction to Applied Mathematics. Wellesley-Cambridge Press, Cambridge, MA.
[3]	Chen, H.C. (2002) A Block Fourier Decomposition Method. PARA 2002, 2367, 351-358. https://doi.org/10.1007/3-540-48051-X_35
[4]	Tolstov, G.P. and Series, F. (1962) Translated from the Russian by R.A. Silverman. Dover Publications, Inc., New York.
[5]	Marco, T. (2017) Properties of the Kronecker Product. https://www.statlect.com/matrix-algebra/Kronecker-product-properties
[6]	Meirovitch, L. (1980) Computational Methods in Structural Dynamics. Sijthoff & Noordhoff, Maryland.
[7]	Milne, W.E. (1970) Numerical Solution of Differential Equations. Dover Publications, Inc., New York.
[8]	Wait, R. and Mitchell, A.R. (1985) Finite Element Analysis and Applications. John Wiley & Sons, New York.
[9]	Chen, F. and Goshtasby, A. (1989) A Parallel B-Spline Surface Fitting Algorithm. ACM Transactions on Graphics, 8, 41-50. https://doi.org/10.1145/49155.214377