# A new block diagonal preconditioner for a class of $3\times 3$ block saddle
point problems

Maryam Abdolmaleki  Faculty of Intelligent Systems Engineering and Data
Science, Persian Gulf University, Bushehr, Iran<EMAIL_ADDRESS>Saeed
Karimi∗ Faculty of Intelligent Systems Engineering and Data Science, Persian
Gulf University, Bushehr, Iran<EMAIL_ADDRESS>Davod Khojasteh Salkuyeh
Faculty of Mathematical Sciences, University of Guilan, Rasht, Iran
Center of Excellence for Mathematical Modelling, Optimization and
Combinational Computing (MMOCC), University of Guilan, Rasht, Iran
<EMAIL_ADDRESS>
(Date: January 1, 2020)

###### Abstract.

We study the performance of a new block preconditioner for a class of $3\times
3$ block saddle point problems which arise from finite element methods for
solving time-dependent Maxwell equations and some other practical problems. We
also estimate the lower and upper bounds of eigenvalues of the preconditioned
matrix.  Finally, we examine our new preconditioner to accelerate the
convergence speed of the GMRES method which shows the effectiveness of the
preconditioner.

###### Key words and phrases:

$3\times 3$ block saddle point, GMRES, block diagonal, preconditioner,
eigenvalue

###### 1991 Mathematics Subject Classification:

65F10, 65F50, 65F08.

##  1\. Introduction

Consider the following $3\times 3$ block saddle point problems:

$\mathcal{B}u:\equiv\left(\begin{array}[]{ccc}A&B^{T}&0\\\ B&0&C^{T}\\\
0&C&0\\\ \end{array}\right)\left(\begin{array}[]{c}x\\\ y\\\ z\\\
\end{array}\right)=\left(\begin{array}[]{c}f\\\ g\\\ h\\\ \end{array}\right),$
(1.1)

where $A\in\mathbb{R}^{n\times n}$ is a  symmetric positive definite (SPD)
matrix, and $B\in\mathbb{R}^{m\times n}$ and $C\in\mathbb{R}^{l\times m}$ have
full row rank. In addition, $f\in\mathbb{R}^{n}$, $g\in\mathbb{R}^{m}$ and
$h\in\mathbb{R}^{l}$ are given vectors, and $u$ is an unknown vector which is
to be determined. Furthermore, we presume that the matrices $A$ and $B$ are
large and sparse. It is not difficult to check that under the above conditions
the coefficient matrix $\mathcal{B}$ is nonsingular and as a result system
(1.1) has a unique solution [27]. In this work, we focus on preconditioned
Krylov-subspace methods, especially,  the preconditioned GMRES method (see
[23, 24]).

Linear systems of the form (1.1) appear in a variety of scientific and
engineering problems, for instance, full discrete finite element methods for
solving the time-dependent Maxwell equations with discontinuous coefficients
[1, 8, 13], the following quadratic program [15]:

$\min\\{\frac{1}{2}x^{T}Ax+r^{T}x+q^{T}y\\},\quad$ $s.t.\quad
Bx+c^{T}y=b,\quad x\in\mathbb{R}^{n},\quad y\in\mathbb{R}^{l},$

where $r\in\mathbb{R}^{n}$ and $q\in\mathbb{R}^{l}$ are given vectors, Picard
iteration schemes for variational formulation of the stationary incompressible
magnetohydrodynamics system [21] and least squares problems [28].

Linear systems of $2\times 2$ block form,

$\left(\begin{array}[]{cc}\tilde{A}&{\tilde{B}}^{T}\\\ \tilde{B}&-\tilde{C}\\\
\end{array}\right)\left(\begin{array}[]{c}x\\\ y\\\
\end{array}\right)=\left(\begin{array}[]{c}f\\\ g\\\ \end{array}\right),$
(1.2)

known as traditional saddle point problems, have been extensively studied for
decades, where $\tilde{A}$ and $\tilde{C}$ are positive and positive semi-
definite matrices respectively, and $\tilde{B}$ is a full row rank matrix.
Constructing preconditioners to improve the convergence speed of Krylov-
subspace methods for solving (1.2) has produced a considerable amount of
literature, e.g., shift-splitting preconditioners [5, 11, 12, 29], block
triangular preconditioners [10, 14, 20, 25], inexact constraint
preconditioners[2, 6, 30] and so on. It is obvious that the $3\times 3$ block
linear system (1.1) can be seen as a special case of the traditional $2\times
2$ form (1.2) using the following partitioning strategies,

$\begin{pmatrix}A&B^{T}&\vdots&0\\\ B&0&\vdots&C^{T}\\\
\cdots&\cdots&\vdots&\cdots\\\ 0&C&\vdots&0\end{pmatrix},\qquad$ (1.3)

or

$\begin{pmatrix}A&\vdots&B^{T}&0\\\ \cdots&\cdots&\cdots&\cdots\\\
B&\vdots&0&C^{T}\\\ 0&\vdots&C&0\end{pmatrix}.$ (1.4)

However,  most of the methods for the latter matrices cannot be directly
applied to (1.1). This is because the properties of the sub-matrix in (1.3)
and (1.4) are different from the traditional form (1.2). Indeed, the $(1,1)$
leading block matrix in (1.3) is symmetric indefinite and a standard saddle
point matrix. Thus, the linear system with the coefficient matrix (1.3) can be
considered as a double saddle point problem studied recently [3, 4]. In
addition, the $(2,2)$ block in (1.4) is a symmetric matrix and $(1,2)$ block
is rank deficient. Therefore, it is extremely interesting to find new
efficient preconditioners for the $3\times 3$ block saddle point problems
(1.1).

Some preconditioners have been studied to accelerate convergence rate of the
Krylov-subspace methods for solving this class of $3\times 3$ systems, for
example [9, 16, 17]. Recently, Huang and Ma [18] proposed two block diagonal
(BD) preconditioners,

$\mathcal{P}_{BD1}=\begin{pmatrix}A&0&0\\\ 0&S&0\\\
0&0&CS^{-1}C^{T}\end{pmatrix}\qquad\text{and}\qquad\mathcal{P}_{BD2}=\begin{pmatrix}\hat{A}&0&0\\\
0&\hat{S}&0\\\ 0&0&CS^{-1}C^{T}\end{pmatrix},$ (1.5)

for solving (1.1) in which $S=BA^{-1}B^{T}$, $\hat{A}$ and $\hat{S}$ are SPD
approximations of $A$ and $S$, respectively. They also derive all the
eigenpairs of preconditioned matrix. Subsequently, Xie and Li [27] proposed
three new preconditioners for solving the linear system (1.1) that can be seen
as follows

$\displaystyle\mathcal{P}_{1}=\begin{pmatrix}A&0&0\\\ B&-S&C^{T}\\\
0&0&-CS^{-1}C^{T}\end{pmatrix},$
$\displaystyle\mathcal{P}_{2}=\begin{pmatrix}A&0&0\\\ B&-S&C^{T}\\\
0&0&CS^{-1}C^{T}\end{pmatrix},$ (1.6)
$\displaystyle\mathcal{P}_{3}=\begin{pmatrix}A&B^{T}&0\\\ B&-S&C^{T}\\\
0&0&-CS^{-1}C^{T}\end{pmatrix},$

where $S=BA^{-1}B^{T}$. They analyzed spectral properties of corresponding
preconditioned matrices and showed that the proposed preconditioners
significantly accelerate the convergence rate of GMRES method.  However, when
the inner systems are solved inexactly, the elapsed CPU time is increased
drastically and often give unacceptable solutions.

Here, we consider the following equivalent form of (1.1):

${\mathcal{A}}u:\equiv\left(\begin{array}[]{ccc}A&B^{T}&0\\\ -B&0&-C^{T}\\\
0&C&0\\\ \end{array}\right)\left(\begin{array}[]{c}x\\\ y\\\ z\\\
\end{array}\right)=\left(\begin{array}[]{c}f\\\ -g\\\ h\\\
\end{array}\right)\equiv b.$ (1.7)

Although the coefficient matrix of the system (1.7) is not symmetric, it has
some desirable properties. For instance, the matrix $\mathcal{A}$ is positive
semi-definite, i.e., $\mathcal{A}+{\mathcal{A}}^{T}$ is symmetric positive
semi-definite. This is a signification for the GMRES method. In fact, the
restarted version of GMRES($m$) converges for all $m\geq 1$.  It is noticeable
that setting up the exact preconditioners mentioned above are very time-
consuming and the inexact preconditioners need approximations of $A$ and $S$.
In this paper, we establish a new block preconditioner for solving the linear
system (1.7) which is easy to implement and has much better computing
efficiently than the preconditioners studied recently.

It is noteworthy that, few analytical results on spectral bounds are available
for a $3\times 3$ block matrix of the form $\mathcal{A}$. In fact, contrary to
the case for $2\times 2$ block matrix, we get a cubic equation from the eigen-
system of a $3\times 3$ block matrix. Therefore, estimating the bounds of
eigenvalues brings some difficulties. Consider the following monic polynomial
of degree $n\geq 2$:

$p(z)=z^{n}+\sum_{k=0}^{n-1}a_{k}z^{k},\qquad a_{i}\in\mathbb{C},\quad
i=0,1,\ldots,n-1.$ (1.8)

Finding approximately the roots of $p(z)$ through simple operations with its
coefficients has led to publish a plenty of studies, for example, see the
comprehensive surveys [22, 26] and references therein for further details. We
will derive the spectral bounds of corresponding preconditioned system by
using one of the classical and sharp bounds that has been used to obtain
simple lower and upper bounds on the absolute value of the roots of $p(z)$.

This paper is divided into three sections, schemed as follows. In Section 2,
the new preconditioner is presented and the clustering properties of
corresponding preconditioned system are discussed. Numerical results are given
in Section 3 to demonstrate the effectiveness of the new preconditioner. The
paper is ended by some concluding remarks in Section 4.

We end this section with an introduction of some notation that will be used in
the subsequent sections. The symbol $x^{*}$ is used for the conjugate
transpose of the vector $x$. For any square matrix $A$ with real eigenvalues,
the minimum and maximum eigenvalues of $A$ are indicated by
$\lambda_{\min}(A)$ and $\lambda_{\max}(A)$, respectively. The norm $\|.\|$
indicates the  Euclidean norm.  Moreover, we use Matlab notation $(x;y;z)$ to
denote the vector $(x^{T},y^{T},z^{T})^{T}$.

##  2\. The new block diagonal preconditioner

We propose the following block preconditioner for solving the linear system
(1.7)

$\mathcal{M}=\begin{pmatrix}A&0&0\\\ 0&\alpha I+\beta BB^{T}&0\\\ 0&0&\alpha
I+\beta CC^{T}\end{pmatrix},$ (2.1)

where $\alpha,\beta>0$.  The main advantage of the preconditioner
$\mathcal{M}$ over the preconditioners mentioned in Section 1 is that it is
free of the Schur complement matrix $S=BA^{-1}B^{T}$ and easy to implement.
Concerning the clustering properties of the eigenvalues of the preconditioned
matrix $\mathcal{M}^{-1}\mathcal{A}$, we have the following theorems.

######  Theorem 2.1.

[19] Let $p(z)=z^{n}+\sum_{k=0}^{n-1}a_{k}z^{k}$ be a monic polynomials with
complex coefficients and $\lambda$ be any root of $p(z)$. Then $|\lambda|$
satisfies the following inequalities

  * •

Cauchy’s lower and upper bounds

$\displaystyle\frac{\left|a_{0}\right|}{\max\\{1,\left|a_{0}\right|+\left|a_{1}\right|,\left|a_{0}\right|+\left|a_{2}\right|,\ldots,\left|a_{0}\right|+\left|a_{n-1}\right|\\}}$
$\displaystyle\leq\left|\lambda\right|\leq$
$\displaystyle\max\\{\left|a_{0}\right|,1+\left|a_{1}\right|,\ldots,1+\left|a_{n-1}\right|\\}.$

  * •

Montel’s lower and upper bounds

$\frac{\left|a_{0}\right|}{\max\\{\left|a_{0}\right|,1+\left|a_{1}\right|+\left|a_{2}\right|+\cdots+\left|a_{n-1}\right|\\}}\leq\left|\lambda\right|\leq\max\\{1,\left|a_{0}\right|+\left|a_{1}\right|+\cdots+\left|a_{n-1}\right|\\}.$

  * •

Carmichael-Mason’s lower and upper bounds

$\frac{\left|a_{0}\right|}{\sqrt{1+\left|a_{0}\right|^{2}+\left|a_{1}\right|^{2}+\cdots+\left|a_{n-1}\right|^{2}}}\leq\left|\lambda\right|\leq\sqrt{1+\left|a_{0}\right|^{2}+\left|a_{1}\right|^{2}+\cdots+\left|a_{n-1}\right|^{2}}.$

  * •

Frobenius’ lower and upper bounds

$\frac{\left|a_{0}\right|}{\sqrt{1+(n-1)\left|a_{0}\right|^{2}+\left|a_{1}\right|^{2}+\cdots+\left|a_{n-1}\right|^{2}}}\leq\left|\lambda\right|\leq\sqrt{(n-1)+\left|a_{0}\right|^{2}+\cdots+\left|a_{n-1}\right|^{2}}.$

It is noted that, Cauchy’s bounds are essentially the most sharpest ones in
Theorem 2.1 and have been used for testing the sharpness of other new bounds.
Hence, we will utilize these bounds to locate the eigenvalues of corresponding
preconditioned matrix.

######  Theorem 2.2.

Suppose that $A\in\mathbb{R}^{n\times n}$ is SPD, and $B\in\mathbb{R}^{m\times
n}$ and $C\in\mathbb{R}^{l\times m}$ have full row rank. Then the
preconditioned matrix $\mathcal{M}^{-1}\mathcal{A}$ has an eigenvalue $1$ of
algebraic multiplicity $n-m$, i.e. $\lambda_{i}^{(1)}=1$, $i=1,2,\ldots,n-m$,
and the corresponding eigenvectors are of the form $(x_{i};0;0),$ where
$\\{x_{i}\\}_{i=1}^{n-m}$ is a basis for the null space of $B$. Moreover, for
the remaining eigenvalues we have the following statements:

  * •

If $p+q>1$, $\qquad\frac{p}{1+2p+q}\leq\left|\lambda-1\right|<2+p+q$,

  * •

If $p+q\leq 1$, $\qquad\quad\frac{p}{2+p}\leq\left|\lambda-1\right|\leq 3$,

where

$p=\frac{y^{*}BA^{-1}B^{T}y}{y^{*}(\alpha I+\beta
BB^{T})y}\qquad\text{and}\qquad q=\frac{y^{*}C^{T}(\alpha I+\beta
CC^{T})^{-1}Cy}{y^{*}(\alpha I+\beta BB^{T})y},$

for some $0\neq y\in\mathbb{R}^{m}$ and $\alpha,\beta>0$.

###### Proof.

Let $(\lambda,(x;y;z))$ be an eigenpair of $\mathcal{M}^{-1}\mathcal{A}$.
Then, we have

$\left(\begin{array}[]{ccc}A&B^{T}&0\\\ -B&0&-C^{T}\\\ 0&C&0\\\
\end{array}\right)\left(\begin{array}[]{c}x\\\ y\\\ z\\\
\end{array}\right)=\lambda\left(\begin{array}[]{ccc}A&0&0\\\ 0&\alpha I+\beta
BB^{T}&0\\\ 0&0&\alpha I+\beta CC^{T}\\\
\end{array}\right)\left(\begin{array}[]{c}x\\\ y\\\ z\\\ \end{array}\right),$

which can be rewritten as

$\displaystyle Ax+B^{T}y=\lambda Ax,$ (2.2a) $\displaystyle
Bx+C^{T}z=-\lambda(\alpha I+\beta BB^{T})y,$ (2.2b) $\displaystyle
Cy=\lambda(\alpha I+\beta CC^{T})z.$ (2.2c)

If $x=0$, then from (2.2a) we obtain $B^{T}y=0$, which shows that $y=0$. This
along with (2.2c) yield that $z=0$, which is a contradiction with the fact
that $(x;y;z)$ is an eigenvector. Next, we complete the proof in the following
cases.

Firstly, we consider the case that $y=0$. Since $\lambda\neq 0$ and $\alpha
I+\beta CC^{T}$ is nonsingular, we deduce that $z=0$. Then (2.2a) and (2.2b)
are reduced to

$Ax=\lambda Ax,\qquad Bx=0,$

respectively. This shows that $\lambda=1$ and the corresponding eigenvectors
are of the form $(x;0;0)$ with $x\in null(B).$

Next, we consider the case that $y\neq 0$. Then from (2.2a), we know that
$\lambda\neq 1$. Noticing that both of matrices $\mathcal{A}$ and $\alpha
I+\beta CC^{T}$ are nonsingular, it follows from (2.2a) and (2.2c) that

$x=\frac{1}{\lambda-1}A^{-1}B^{T}y,\qquad z=\frac{1}{\lambda}(\alpha I+\beta
CC^{T})^{-1}Cy.$

Substituting the preceding equalities into (2.2b), gives

$\frac{1}{\lambda-1}BA^{-1}B^{T}y+\frac{1}{\lambda}C^{T}(\alpha I+\beta
CC^{T})^{-1}Cy=-\lambda(\alpha I+\beta BB^{T})y.$

By premultiplying the above equation by $\lambda(\lambda-1)$ and $y^{*}$, we
obtain the following cubic equation:

$\lambda^{3}-\lambda^{2}+(p+q)\lambda-q=0,$ (2.3)

where

$p=\frac{y^{*}BA^{-1}B^{T}y}{y^{*}(\alpha I+\beta
BB^{T})y}\qquad\text{and}\qquad q=\frac{y^{*}C^{T}(\alpha I+\beta
CC^{T})^{-1}Cy}{y^{*}(\alpha I+\beta BB^{T})y}.$ (2.4)

To obtain the clustering properties of the eigenvalues around the point
$(1,0)$ we set $\mu=\lambda-1$. By substituting $\mu$ in (2.3) we get the
following cubic equation

$\mu^{3}+2\mu^{2}+(1+p+q)\mu+p=0,$ (2.5)

and find upper and lower bounds for $\mu$. From Theorem 2.1, we deduce that
the Cauchy’s lower and upper bounds for $\mu$ are

$\displaystyle\frac{p}{1+2p+q}\leq\left|\mu\right|<$ $\displaystyle 2+p+q,$
$\displaystyle\text{if}\quad p+q>1,$ (2.6a)
$\displaystyle\frac{p}{2+p}\leq\left|\mu\right|\leq$ $\displaystyle 3,$
$\displaystyle\text{if}\quad p+q\leq 1,$ (2.6b)

which completes the proof. ∎

As we see in Theorem 2.2, the lower and upper bounds of some eigenvalues of
the preconditioned matrix $\mathcal{M}^{-1}\mathcal{A}$ depend on the
parameters $\alpha$ and $\beta$. Using (2.4) and the Schur decomposition of
matrices $A$ and $\alpha I+\beta CC^{T}$, and straightforward computations the
following upper and lower bounds for $p$ and $q$ are obtained as

$\frac{\lambda_{\min}(A^{-1})\lambda_{\min}(BB^{T})}{\alpha+\beta\lambda_{\max}(BB^{T})}\leq
p\leq\frac{\lambda_{\max}(A^{-1})\lambda_{\max}(BB^{T})}{\alpha+\beta\lambda_{\min}(BB^{T})},$

and

$\frac{{\lambda_{\min}(C^{T}C)}(\alpha+\beta\lambda_{\max}(CC^{T}))^{-1}}{\alpha+\beta\lambda_{\max}(BB^{T})}\leq
q\leq\frac{\lambda_{\max}(CC^{T})(\alpha+\beta\lambda_{\min}(CC^{T}))^{-1}}{\alpha+\beta\lambda_{\min}(BB^{T})}.$

Hence, we conclude that

$\frac{\lambda_{\min}(A^{-1})\lambda_{\min}(BB^{T})+{\lambda_{\min}(C^{T}C)}(\alpha+\beta\lambda_{\max}(CC^{T}))^{-1}}{\alpha+\beta\lambda_{\max}(BB^{T})}\\\
\leq
p+q\leq\frac{\lambda_{\max}(A^{-1})\lambda_{\max}(BB^{T})+\lambda_{\max}(CC^{T})(\alpha+\beta\lambda_{\min}(CC^{T}))^{-1}}{\alpha+\beta\lambda_{\min}(BB^{T})}.$
(2.7)

Now, we consider the following two cases:

case I: $p+q>1$: From (2.7) we deduce that all we need is to seek the
parameters $\alpha$ and $\beta$ such that

$1<\frac{\lambda_{\min}(A^{-1})\lambda_{\min}(BB^{T})+{\lambda_{\min}(C^{T}C)}(\alpha+\beta\lambda_{\max}(CC^{T}))^{-1}}{\alpha+\beta\lambda_{\max}(BB^{T})}.$

For the sake of the simplicity, let

$\eta=\lambda_{\max}(BB^{T})+\lambda_{\max}(CC^{T}),\qquad\theta=\lambda_{\max}(BB^{T})\lambda_{\max}(CC^{T}),$
$\kappa=\lambda_{\min}(A^{-1})\lambda_{\min}(BB^{T}),\qquad\gamma=\lambda_{\min}(A^{-1})\lambda_{\min}(BB^{T})\lambda_{\max}(CC^{T}),$
${\zeta=\lambda_{\min}(C^{T}C)}.$

Therefore, we have the following inequality

$\alpha^{2}+(\eta\beta-\kappa)\alpha+\theta\beta^{2}-\gamma\beta-\zeta<0.$

It is clear that matrices $CC^{T}$ and $BB^{T}$ are both SPD under our
assumptions. Hence, we infer that $\eta\neq 0$. Therefore, we can set

$\beta=\frac{\kappa}{\eta}\equiv\frac{\lambda_{\min}(A^{-1})\lambda_{\min}(BB^{T})}{\lambda_{\max}(BB^{T})+\lambda_{\max}(CC^{T})}$

and derive that

$0<\alpha<\sqrt{\frac{\gamma\kappa}{\eta}-\frac{\theta\kappa^{2}}{\eta^{2}}+\zeta}.$

Note that $\theta\kappa<\eta\gamma$, and then it follows that
$\frac{\kappa}{\eta}(\gamma-\frac{\theta\kappa}{\eta})>0$. Therefore, the
above inequality can be easily deduced.

case II: $p+q\leq 1$: In this case, it is enough to find $\alpha$ and $\beta$
such that

$\frac{\lambda_{\max}(A^{-1})\lambda_{\max}(BB^{T})+\lambda_{\max}(CC^{T})(\alpha+\beta\lambda_{\min}(CC^{T}))^{-1}}{\alpha+\beta\lambda_{\min}(BB^{T})}\leq
1.$

Hence, we have

$\alpha^{2}+(\eta^{\prime}\beta-\kappa^{\prime})\alpha+\beta^{2}\theta^{\prime}-\beta\gamma^{\prime}-\zeta^{\prime}\geq
0,$

where

$\eta^{\prime}=\lambda_{\min}(BB^{T})+\lambda_{\min}(CC^{T}),\quad\theta^{\prime}=\lambda_{\min}(CC^{T})\lambda_{\min}(BB^{T}),$
$\kappa^{\prime}=\lambda_{\max}(A^{-1})\lambda_{\max}(BB^{T}),\quad\gamma^{\prime}=\lambda_{\max}(A^{-1})\lambda_{\max}(BB^{T})\lambda_{\min}(CC^{T}),$
$\zeta^{\prime}=\lambda_{\max}(CC^{T}).$

Similar to the argument in the previous case we set

$\beta=\frac{\kappa^{\prime}}{\eta^{\prime}}\equiv\frac{\lambda_{\max}(A^{-1})\lambda_{\max}(BB^{T})}{\lambda_{\min}(BB^{T})+\lambda_{\min}(CC^{T})},$

and finally conclude that

$\alpha\geq\sqrt{\frac{\gamma^{\prime}\kappa^{\prime}}{\eta^{\prime}}-\frac{\theta^{\prime}\kappa^{\prime^{2}}}{\eta^{\prime^{2}}}+\zeta^{\prime}}.$

According to the above discussion we state the following theorem that gives
conditions on $\alpha$ and $\beta$ under which $p+q>1$ or $p+q\leq 1$.

######  Theorem 2.3.

We have the following cases:

case I: If

$\beta=\frac{\kappa}{\eta}\qquad\text{and}\qquad
0<\alpha<\sqrt{\frac{\gamma\kappa}{\eta}-\frac{\theta\kappa^{2}}{\eta^{2}}+\zeta},$

then, $p+q>1$.

case II: If

$\beta=\frac{\kappa^{\prime}}{\eta^{\prime}}\qquad\text{and}\qquad\alpha\geq\sqrt{\frac{\gamma^{\prime}\kappa^{\prime}}{\eta^{\prime}}-\frac{\theta^{\prime}\kappa^{\prime^{2}}}{\eta^{\prime^{2}}}+\zeta^{\prime}},$

then, $p+q\leq 1$.

In each iteration of a Krylov subspace method for solving the preconditioned
system $\mathcal{M}^{-1}{\mathcal{A}}u=\mathcal{M}^{-1}b$ we need to compute a
vector of the form $(z_{1};z_{2};z_{3})=\mathcal{M}^{-1}(r_{1};r_{2};r_{3})$
where $r_{1}\in\mathbb{R}^{n}$, $r_{2}\in\mathbb{R}^{m}$ and
$r_{3}\in\mathbb{R}^{\ell}$. To do so, it is enough to solve the system
$\mathcal{M}(z_{1};z_{2};z_{3})=(r_{1};r_{2};r_{3})$ for
$(z_{1};z_{2};z_{3})$. Since, $\mathcal{M}$ is a block diagonal matrix,
solution of the system is reduced to the solution of three systems with the
coefficient matrices $A$, $\alpha I+\beta BB^{T}$ and $\alpha I+\beta CC^{T}$,
which are all SPD. By summarizing the above notes we can state Algorithm 1.

Algorithm 1 Computation of
$(z_{1};z_{2};z_{3})=\mathcal{M}^{-1}(r_{1};r_{2};r_{3})$

  1. (1)

Solve $\mathcal{M}_{1}z_{1}\equiv Az_{1}=r_{1}$ for $z_{1}$.

  2. (2)

Solve $\mathcal{M}_{2}z_{2}\equiv(\alpha I+\beta BB^{T})z_{2}=r_{2}$ for
$z_{2}$.

  3. (3)

Solve $\mathcal{M}_{3}z_{3}\equiv(\alpha I+\beta CC^{T})z_{3}=r_{3}$ for
$z_{3}$.

Since the coefficient matrices of the subsystems in Algorithm 1 are SPD, they
can be solved exactly using either the Cholesky factorization or inexactly
using the conjugate gradient (CG) iteration method. It is noted that, the
shift matrix $\alpha I$ in the steps 3 and 4 of Algorithm 1 increases the
convergence speed of the CG method considerably.

##  3\. Numerical experiments

In this section, we experimentally compare the effectiveness our proposed
preconditioner with two preconditioners $\mathcal{P}_{BD2}$ and
$\mathcal{P}_{2}$ defined, respectively, in Eqs. (1.5) and (1.6) for solving
the saddle point linear system (1.7). All the numerical experiments were
computed in double precision using some Matlab codes on a Laptop with Intel
Core i7 CPU 2.9 GHz, 16GB RAM.

We apply the preconditioners to accelerate the convergence of the GMRES
(flexible and full versions) method [23, 24]. In the implementation of the
preconditioners $\mathcal{M}$, $\mathcal{P}_{BD2}$ and $\mathcal{P}_{2}$
within the GMRES method three subsystems with SPD coefficient matrices need to
be solved. When the subsystems are solved inexactly we apply the flexible
GMRES (FGMRES) method to solve the preconditioned system and if the subsystems
are solved exactly we employ the full version of GMRES (Full-GMRES) method.

In the subsequent presented numerical results, the right-hand side vector $b$
is set to $b=\mathcal{A}u^{\ast}$, where $u^{\ast}$ is a vector of all ones.
We use a null vector as an initial guess for the GMRES method and the
iteration is stopped once the relative residual 2-norm satisfies

$R_{k}=\frac{\|b-\mathcal{A}u^{(k)}\|_{2}}{\|b\|_{2}}<10^{-6},$ (3.1)

where $u^{(k)}$ is the $k$th computed approximate solution. The maximum number
of iterations is set to be $maxit=1000$. To show the accuracy of the methods
we report the values

$E_{k}=\frac{\|u^{(k)}-u^{*}\|_{2}}{\|u^{*}\|_{2}}.$ (3.2)

In the implementation of the preconditioners $\mathcal{M}$,
$\mathcal{P}_{BD2}$ and $\mathcal{P}_{2}$ three subsystems with SPD
coefficient matrices need to be solved in each iteration of the GMRES method.
These systems are solved using the Cholesky factorization of matrices, when we
apply the Full-GMRES method for solving the preconditioned system. On the
other hand, the subsystems are solved using the CG method, when we employ the
FGMRES method for solving the preconditioned system. In this case, for solving
the subsystems the initial guess is set to be a zero vector and the iterations
is stopped as soon as the residual norm is reduced by a factor of $10^{3}$ or
the number of iterations exceeds $500$.

Numerical results are presented in the tables. To show the effectiveness of
the preconditioners we also report the numerical results of Full-GMRES without
any preconditioner.

######  Example 1.

[18, 9] Consider the saddle point problem (1.7) for which

$A=\begin{pmatrix}I\otimes T+T\otimes I&0\\\ 0&I\otimes T+T\otimes
I\end{pmatrix}\in\mathbb{R}^{2p^{2}\times 2p^{2}},$ $B=\begin{pmatrix}I\otimes
F&F\otimes I\end{pmatrix}\in\mathbb{R}^{p^{2}\times 2p^{2}},\quad C=E\otimes
F\in\mathbb{R}^{p^{2}\times p^{2}}$

and

$T=\frac{1}{h^{2}}\text{tridiag}(-1,2,-1)\in\mathbb{R}^{p\times p},\quad
F=\frac{1}{h}\text{tridiag}(0,1,-1)\in\mathbb{R}^{p\times p}$
$E=\text{diag}(1,p+1,\ldots,p^{2}-p+1),$

with $\otimes$ being the Kroneker product symbol and $h={1}/{(p+1)}$ the
discretization meshsize. Here, the total number of unknowns is $4p^{2}$ and
$m=l$.

Numerical results are presented in Table 1 for different values of $p$. This
table shows the number of iterations (denoted by “Iters”) and CPU time
(denoted by “CPU”) of the FGMRES without preconditioning, and with the
preconditioners $\mathcal{M}$, $\mathcal{P}_{2}$ and $\mathcal{P}_{BD2}$. We
set $\hat{S}=B{diag(A)}^{-1}B^{T}$. Moreover, we provide the elapsed CPU time
(denoted by “Prec.CPU”) for setting up the preconditioners $\mathcal{P}_{2}$
and $\mathcal{P}_{BD2}$ and the total CPU time (denoted by “Total.CPU”).
Furthermore, for the preconditioner $\mathcal{M}$ we set $\alpha=10^{-3}$ and
$\beta=1$. In the all tables, the symbols $\dagger$ and ${\ddagger}$ are used
to indicate that the method has not converged in 1000 seconds and $maxit$,
respectively. As we see, the preconditioner $\mathcal{M}$ outperforms the
other examined preconditioners from the iteration counts, CPU time and the
accuracy of computed solution points of view. It should be mentioned that for
$p=64,128,256$ the FGMRES method without preconditioning, and the FGMRES with
the preconditioners $\mathcal{P}_{2}$ and $\mathcal{P}_{BD2}$ fail to converge
in 1000 iterations. Therefore, our preconditioner is more effective and
practical than the preconditioners $\mathcal{P}_{2}$ and $\mathcal{P}_{BD2}$
for solving saddle point problems of the form (1.7).

Numerical results of the Full-GMRES method in conjunction with the three
preconditioners $\mathcal{M}$, $\mathcal{P}_{2}$ and $\mathcal{P}_{BD2}$ are
shown in Table 2. We observe that the preconditioner $\mathcal{M}$ has
provided quite suitable results and this results are in good agreement with
what we claimed above for the preconditioned FGMRES method. In addition, when
$p$ is large the preconditioned GMRES method fails to converge for the
preconditioners $\mathcal{P}_{2}$ and $\mathcal{P}_{BD2}$ in 1000 seconds.
However, the FGMRES method with the new preconditioner requires less CPU time.
It is noted that we could not set up the preconditioners $\mathcal{P}_{2}$ and
$\mathcal{P}_{BD2}$ for $p=512$, because of memory limitation.

Table 1. Numerical results of FGMRES for Example 1.  $p$ | 16 | 32 | 64 | 128 | 256
---|---|---|---|---|---
No Preconditioning
Iters | 425 | 949 | ${\ddagger}$ | ${\ddagger}$ | ${\ddagger}$
CPU | 0.24 | 5.02 | $\dagger$ | $\dagger$ | $\dagger$
$R_{k}$ | $8.6\times 10^{-7}$ | $9.9\times 10^{-7}$ | - | - | -
$E_{k}$ | $2.6\times 10^{-6}$ | $2.4\times 10^{-5}$ | - | - | -
$\mathcal{M}$
Iters | 109 | 80 | 65 | 71 | 78
CPU | 0.15 | 0.25 | 1.17 | 6.66 | 36.99
$R_{k}$ | $7.0\times 10^{-7}$ | $8.4\times 10^{-7}$ | $7.1\times 10^{-7}$ | $9.0\times 10^{-7}$ | $1.0\times 10^{-6}$
$E_{k}$ | $3.0\times 10^{-7}$ | $6.0\times 10^{-7}$ | $1.5\times 10^{-6}$ | $8.3\times 10^{-6}$ | $1.4\times 10^{-5}$
$\mathcal{P}_{2}$
Iters | 124 | 616 | ${\ddagger}$ | $-$ | $-$
Prec.CPU | 0.0006 | 0.0015 | 0.0082 | $-$ | $-$
CPU | 0.81 | 25.53 | 288.93 | $\dagger$ | $\dagger$
Total.CPU | 0.81 | 25.53 | 288.93 | $-$ | $-$
$R_{k}$ | $9.9\times 10^{-7}$ | $8.3\times 10^{-7}$ | $1.8\times 10^{-1}$ | $-$ | $-$
$E_{k}$ | $1.7\times 10^{-6}$ | $1.7\times 10^{-6}$ | $6.8\times 10^{-1}$ | $-$ | $-$
$\mathcal{P}_{BD2}$
Iters | 163 | 785 | ${\ddagger}$ | $-$ | $-$
Prec.CPU | 0.0006 | 0.0015 | 0.0082 | $-$ | $-$
CPU | 0.99 | 34.66 | 286.02 | $\dagger$ | $\dagger$
Total.CPU | 0.99 | 34.36 | 286.03 | $-$ | $-$
$R_{k}$ | $8.1\times 10^{-7}$ | $7.1\times 10^{-7}$ | $1.2\times 10^{-1}$ | $-$ | $-$
$E_{k}$ | $1.0\times 10^{-6}$ | $1.1\times 10^{-6}$ | $6.5\times 10^{-1}$ | $-$ | $-$
Table 2. Numerical results of Full-GMRES for Example 1. $p$ | 16 | 32 | 64 | 128 | 256
---|---|---|---|---|---
No Preconditioning
Iters | 425 | 949 | ${\ddagger}$ | ${\ddagger}$ | ${\ddagger}$
CPU | 0.21 | 5.35 | 14.27 | 45.34 | 182.99
$R_{k}$ | $8.6\times 10^{-7}$ | $9.9\times 10^{-7}$ | $2.7\times 10^{-3}$ | $6.7\times 10^{-3}$ | $4.9\times 10^{-2}$
$E_{k}$ | $2.6\times 10^{-6}$ | $2.4\times 10^{-5}$ | $1.8\times 10^{-1}$ | $5.5\times 10^{-1}$ | $7.8\times 10^{-1}$
$\mathcal{M}$
Iters | 109 | 75 | 54 | 60 | 74
Prec.CPU | 0.0016 | 0.0037 | 0.0389 | 0.2334 | 1.8452
CPU | 0.05 | 0.08 | 0.24 | 1.67 | 14.08
Total.CPU | 0.05 | 0.08 | 0.28 | 1.91 | 15.92
$R_{k}$ | $7.3\times 10^{-7}$ | $9.9\times 10^{-7}$ | $7.8\times 10^{-7}$ | $8.8\times 10^{-7}$ | $6.0\times 10^{-7}$
$E_{k}$ | $2.9\times 10^{-7}$ | $9.3\times 10^{-7}$ | $2.5\times 10^{-6}$ | $5.5\times 10^{-6}$ | $7.8\times 10^{-6}$
$\mathcal{P}_{2}$
Iters | 115 | 402 | ${\ddagger}$ | ${\ddagger}$ | -
Prec.CPU | 0.0286 | 0.2170 | 6.0734 | 261.1566 | -
CPU | 0.05 | 1.77 | 42.95 | 454.22 | $\dagger$
Total.CPU | 0.08 | 1.98 | 49.02 | 715.38 | -
$R_{k}$ | $8.4\times 10^{-7}$ | $9.5\times 10^{-7}$ | $1.6\times 10^{-3}$ | $7.0\times 10^{-1}$ | -
$E_{k}$ | $1.3\times 10^{-6}$ | $7.9\times 10^{-6}$ | $2.2\times 10^{-4}$ | $1.7\times 10^{0}$ | -
$\mathcal{P}_{BD2}$
Iters | 156 | 569 | ${\ddagger}$ | ${\ddagger}$ | -
Prec.CPU | 0.0286 | 0.2170 | 6.0734 | 261.1566 | -
CPU | 0.05 | 3.15 | 41.72 | 451.88 | $\dagger$
Total.CPU | 0.08 | 3.37 | 47.80 | 713.04 | -
$R_{k}$ | $8.7\times 10^{-7}$ | $8.8\times 10^{-7}$ | $1.9\times 10^{-3}$ | $4.9\times 10^{-1}$ | -
$E_{k}$ | $1.1\times 10^{-6}$ | $1.3\times 10^{-5}$ | $5.5\times 10^{-4}$ | $4.7\times 10^{-1}$ | -

######  Example 2.

[18] Consider the saddle point problem (1.1), for which

$A=\text{diag}(2W^{T}W+D_{1},D_{2},D_{3})\in\mathbb{R}^{n\times n}$

is a block-diagonal matrix,

$B=[E,-I_{2\widetilde{p}},I_{2\widetilde{p}}]\in\mathbb{R}^{m\times
n}\quad\text{and}\quad C=E^{T}\in\mathbb{R}^{\ell\times m}$

are both full row-rank matrices, where $\widetilde{p}=p^{2}$,
$\widehat{p}=p(p+1)$; $W=(w_{ij})\in\mathbb{R}^{\widehat{p}\times\widehat{p}}$
with $w_{ij}=e^{-2((i/3)^{2}+(j/3)^{2})}$; $D_{1}=I_{\widehat{p}}$ is an
identity matrix;
$D_{i}=\text{diag}(d_{j}^{(i)})\in\mathbb{R}^{2\widetilde{p}\times
2\widetilde{p}}$, $i=2,3,$ are diagonal matrix with

$d_{j}^{(2)}=\begin{cases}1,&\text{for}\quad 1\leq j\leq\widetilde{p},\\\
10^{-5}(j-\widetilde{p})^{2},&\text{for}\quad\widetilde{p}+1\leq j\leq
2\widetilde{p},\end{cases}$
$d_{j}^{(3)}=10^{-5}(j+\widetilde{p})^{2},\qquad\text{for}\quad 1\leq j\leq
2\widetilde{p},$

and

$E=\begin{pmatrix}\widehat{E}\otimes I_{p}\\\
I_{p}\otimes\widehat{E}\end{pmatrix},\quad\widehat{E}=\begin{pmatrix}2&-1&&&\\\
&2&1&&\\\ &&\ddots&\ddots&\\\
&&&2&-1\end{pmatrix}\in\mathbb{R}^{p\times(p+1)}.$

Similar to Example 1, the Matrix $\hat{S}$ is set to be
$\hat{S}=B{diag(A)}^{-1}B^{T}$. In addition, we choose the parameters
$\alpha=10^{-1}$ and $\beta=1$. Numerical results for different values of $p$
are listed in Table 3. As we observe, the numerical results illustrate that
the preconditioner $\mathcal{M}$ considerably reduces the CPU time of the
FGMRES method without preconditioning.

We have also applied the Full-GMRES iteration method incorporated with
preconditioners $\mathcal{M}$, $\mathcal{P}_{2}$ and $\mathcal{P}_{BD2}$. In
this case, the subsystems were solved exactly using the Cholesky
factorization. Our numerical results show that, the Full-GMRES method
outperforms the Full-GMRES method in conjunction with the preconditioners (all
the three preconditioners) from the elapsed CPU time point of view for large
values of $p$. Hence, we have not reported the numerical results.

Table 3. Numerical results of FGMRES for Example 2. $p$ | 16 | 32 | 64 | 128 | 256 | 512
---|---|---|---|---|---|---
No Preconditioning
Iters | 186 | 190 | 187 | 180 | 172 | 163
CPU | 0.08 | 0.37 | 0.97 | 2.86 | 17.78 | 96.70
$R_{k}$ | $1.0\times 10^{-6}$ | $9.9\times 10^{-7}$ | $1.0\times 10^{-6}$ | $9.8\times 10^{-7}$ | $9.6\times 10^{-7}$ | $1.0\times 10^{-6}$
$E_{k}$ | $1.3\times 10^{-5}$ | $1.4\times 10^{-5}$ | $1.4\times 10^{-5}$ | $1.4\times 10^{-5}$ | $1.4\times 10^{-5}$ | $1.4\times 10^{-5}$
$\mathcal{M}$
Iters | 70 | 69 | 68 | 65 | 63 | 60
CPU | 0.09 | 0.29 | 0.85 | 2.18 | 10.30 | 45.97
$R_{k}$ | $1.0\times 10^{-6}$ | $9.5\times 10^{-7}$ | $8.8\times 10^{-7}$ | $9.3\times 10^{-7}$ | $8.5\times 10^{-7}$ | $9.5\times 10^{-7}$
$E_{k}$ | $5.6\times 10^{-6}$ | $5.7\times 10^{-6}$ | $5.0\times 10^{-6}$ | $5.5\times 10^{-6}$ | $4.9\times 10^{-6}$ | $5.2\times 10^{-6}$
$\mathcal{P}_{2}$
Iters | 13 | 13 | 13 | 13 | $-$ | $-$
Prec.CPU | 0.0015 | 0.01 | 0.37 | 13.75 | $-$ | $-$
CPU | 0.02 | 0.18 | 2.59 | 34.55 | $\dagger$ | $\dagger$
Total.CPU | 0.03 | 0.20 | 2.96 | 48.31 | $-$ | $-$
$R_{k}$ | $5.1\times 10^{-7}$ | $3.3\times 10^{-7}$ | $3.4\times 10^{-7}$ | $2.5\times 10^{-7}$ | $-$ | $-$
$E_{k}$ | $2.1\times 10^{-6}$ | $1.5\times 10^{-6}$ | $1.1\times 10^{-6}$ | $8.0\times 10^{-7}$ | $-$ | $-$
$\mathcal{P}_{BD2}$
Iters | 19 | 18 | 18 | 18 | $-$ | $-$
Prec.CPU | 0.0015 | 0.01 | 0.37 | 13.75 | $-$ | $-$
CPU | 0.03 | 0.25 | 3.30 | 49.54 | $\dagger$ | $\dagger$
Total.CPU | 0.03 | 0.26 | 3.67 | 63.29 | $-$ | $-$
$R_{k}$ | $1.2\times 10^{-7}$ | $9.6\times 10^{-7}$ | $8.4\times 10^{-7}$ | $9.8\times 10^{-7}$ | $-$ | $-$
$E_{k}$ | $6.1\times 10^{-7}$ | $1.9\times 10^{-6}$ | $1.8\times 10^{-6}$ | $3.2\times 10^{-6}$ | $-$ | $-$

##  4\. Conclusion

A new block diagonal preconditioner has been presented for a class of $3\times
3$ block saddle point problems. This preconditioner is based on augmentation
and performs well in practice. Also, it is easy to implement and has much
better efficiency than the recently existing preconditioners. We have further
estimated the lower and upper bounds of eigenvalues of the preconditioned
matrix. We have examined the new preconditioner to accelerate the convergence
speed of the FGMRES method as well as Full-GMRES. Our numerical tests
illustrate that the proposed preconditioner is quite suitable and is superior
to the other tested preconditioners in the literature.

## Acknowledgments

The authors would like to thank the anonymous referees for their useful
comments and suggestions.

## References

  * [1] F. Assous, P. Degond, E. Heintze, P.A. Raviart, J. Segre, On a finite-element method for solving the three-dimentional Maxwell equations. J. Comput. Phys. 109 (1993), 222-237.
  * [2] Z.-Z. Bai, M.K. Ng, Z.-Q. Wang, Constraint preconditioners for symmetric indefinite matrices. SIAM J. Matrix Anal. Appl. 31 (2009), 410-433.
  * [3] F.P.A. Beik, M. Benzi, Iterative methods for double saddle point systems. SIAM J. Matrix Anal. Appl. 39 (2018), 902-921.
  * [4] F.P.A. Beik, M. Benzi, Block Preconditioners for Saddle Point Systems Arising from Liquid Crystal Directors Modeling. Calcolo, 55 (2018), Article 29.
  * [5] M. Benzi, G.H. Golub, A preconditioner for generalized saddle point problems. SIAM J. Matrix Anal. Appl. 26 (2004), 20-41.
  * [6] L. Bergamaschi, On eigenvalue distribution of constraint-preconditioned symmetric saddle point matrices. Numer. Linear Algebra Appl. 19 (2012), 754-772.
  * [7] E.P. Bertsekas, Nonlinear Programming. 2nd edn, Athena Scientific, (1999).
  * [8] Z.-M. Chen, Q. DU, J. Zou, Finite element methods with matching and nonmatching meshes of Maxwell equation with discontinues coefficients. SIAM J.Numer. Anal. 37 (2000), 1542-1570.
  * [9] Y. Cao, Shift-splitting preconditioners for a class of block three-by-three saddle point problems. Appl. Math. Lett. 96 (2019), 40-46.
  * [10] Z.H. Cao, Positive stable block triangular preconditioners for symmetric saddle point problems. Appl. Numer. Math. 57 (2007), 899-910.
  * [11] Y. Cao, J. Du, Q. Niu, Shift-splitting preconditioners for saddle point problems. J. Comput. Appl. Math. 272 (2014), 239-250.
  * [12] C. Chen, C.-F Ma, A generalized shift-splitting preconditioner for saddle point problems. Appl. Math. Lett. 43 (2015), 49-55.
  * [13] P. Ciarlet, J. Zou, Finite element convergence for Darwin model to Maxwell’s equations. RAIRO Math. Modeling Numer. Anal. 31 (1997), 213-249.
  * [14] H.C. Elman, D.J. Silvester, A.J. Wathen, Performance and analysis of saddle point preconditioners for the discrete steady-state Navier-Stokes equations. Numer. Math. 90 (2002), 665-688.
  * [15] D. Han , X. Yuan, Local linear convergence of the alternating direction method of multipliers for quadratic programs. SIAM J.Numer. Anal. 51 (2013), 3446-3457.
  * [16] N. Huang, Variable parameter Uzawa method for solving a class of block three-by-three saddle point problems. Numerical Algorithms, 85 (2020), 1233-1254.
  * [17] N. Huang, Y.H. Dai, Q. Hu, Uzawa methods for a class of block three-by-three saddle point problems. Numer Linear Algebra Appl. 26 (2019), e2265.
  * [18] N. Huang, C.-F. Ma, Spectral analysis of the preconditioned system for the $3\times 3$ block saddle point problem. Numer. Algorithms 81 (2019), 421-444.
  * [19] R.A. Horn, C.R. Johnson, Matrix Analysis. 2nd Edition, Cambridge University Press, Cambridge, (2012).
  * [20] M.-Q. Jiang, Y. Cao, L.-Q. Yao, On parameterized block triangular preconditioners for generalized saddle point problems. Appl. Math. Comput. 216 (2010), 1777-1789.
  * [21] K.B. Hu, J.C. XU, Structure-preserving finite element methods for stationary MHD models. Math. Comp. 88 (2019), 553-581.
  * [22] M. Madern, Geometry of Polynomials. Math. Surveys Monger., vol. 3, AMS, Providence, RI, USA, (1966).
  * [23] Y. Saad, A flexible inner-outer preconditioned GMRES algorithm. SIAM J. Sci. Comput. 14 (1993), 461-469.
  * [24] Y. Saad, Iterative methods for sparse linear systems. Second edition PWS, New York, 1995.
  * [25] D.K. Salkuyeh, M. Abdolmaleki, S. Karimi. On a splitting preconditioner for saddle point problems. J. Appl. Math. & Informatics. 36 (2018), 459-474.
  * [26] Bl. Sendov, A. Andreev, N. Kjurkchiev, Numerical solution of polynomial equations. in: P.G. Ciarlet, J. L. Lions (Eds.), Handbook of Numerical Analysis, vol. III: Solution of Equations in $\mathbb{R}^{n}$ (Part2), Elsevier, Amsterdam, (1994).
  * [27] X. Xie, H.-B. Li, A note on preconditioning for the $3\times 3$ block saddle point problem. J. Comput. Math. Appl. 79 (2020), 3289-3296.
  * [28] J.-Y. Yuan, Numerical methods for generalized least squares problems. J. Compute. Appl. Math. 66 (1996), 571-584.
  * [29] J.-L. Zhang, C.-Q. Gu, K. Zhang, A relaxed positive-definite and skew-Hermitian splitting preconditioner for saddle point problems. Appl. Math Comput. 249 (2014), 468-479.
  * [30] N.-M. Zhang, P. Shen, Constraint preconditioners for solving singular saddle point problems. J. Comput. Appl. Math. 238 (2013), 116-125.


††thanks: More precisely, the usefulness of linear algebra for the same.

# Simulating the DFT Algorithm for Audio Processing

Omkar Deshpande<EMAIL_ADDRESS>Kharanshu Solanki
<EMAIL_ADDRESS>Sree Pujitha Suribhatla<EMAIL_ADDRESS>Department of Physics, School of Arts and Sciences, Ahmedabad University,
Ahmedabad - 380009
Sanya Zaveri<EMAIL_ADDRESS>Luv Ghodasara<EMAIL_ADDRESS>Department of Comptuter Science, School of Arts and Sciences, Ahmedabad
University, Ahmedabad - 380009

###### Abstract

Since the evolution of digital computers, the storage of data has always been
in terms of discrete bits that can store values of either 1 or 0. Hence, all
computer programs (such as MATLAB), convert any input continuous signal into a
discrete dataset. Applying this to oscillating signals, such as audio, opens a
domain for processing as well as editing. The Fourier transform, which is an
integral over infinite limits, for the use of signal processing is discrete.
The essential feature of the Fourier transform is to decompose any signal into
a combination of multiple sinusoidal waves that are easy to deal with. The
discrete Fourier transform (DFT) can be represented as a matrix, with each
data point acting as an orthogonal point, allowing one to perform complicated
transformations on individual frequencies. Due to this formulation, all the
concepts of linear algebra and linear transforms prove to be extremely useful
here. In this paper, we first explain the theoretical basis of audio
processing using linear algebra, and then focus on a simulation coded in
MATLAB, to process and edit various audio samples. The code is open ended and
easily expandable by just defining newer matrices which can transform over the
original audio signal. Finally, this paper attempts to highlight and briefly
explain the results that emerge from the simulation.

††preprint: APS/123-QED

##  I Introduction

The Fourier transform, in the most basic sense, is used as a tool to break
down a waveform into its composite amplitudes and frequencies. It has the
ability to construct any repetitive function of time or a waveform as a
combination of individual sinusoidal vibrations, which are individually
functions of time.

Although Joseph Fourier first used the notion of a Fourier transform to
explain certain scenarios pertaining to heat transfer [1], it now finds its
uses in a huge number of fields. These include practical applications like
signal processing, optics, image processing, as well as theoretical
applications like calculating eigenvalues of unitary operators in quantum
algorithms.

Fourier transform finds its most crucial application in the field of signal
processing, i.e., to study a given signal which is a waveform (a function of
time), and then obtain useful and important information from these waveforms.
Waveforms can be the range of any function changing with time - from
electromagnetic waves to sound waves and waves produced in strings to
electrical signals.

Since signals and waveforms are so essential and are present everywhere in the
world around us, it is necessary that one understands signal processing and
the role of the Fourier transform as an important tool to analyze it. More
precisely computers are equipped to implement the discrete Fourier transform,
which deals with discrete numbers of equally spaced samples, rather than
continuous systems.

This paper aims to analyze an audio form by converting it into a waveform and
then applying the DFT algorithm to it, via a coded simulation. Once the
waveform is analyzed, the code obtains the number of highest amplitudes or
frequencies that constitute the signal. Further, we try to demonstrate how and
why these composite amplitudes are important in various fields. These
objectives are to be carried out with the aid of the MATLAB programming
languge, which provides simple and helpful functions to carry out the same.

##  II Theory and Significance

###  II.1 Discrete Data and Digital Signals

Data can be broadly classified into discrete and continuous sets of data. If a
piece of data is defined to be between two particular real values, and
further, if the data can be anything between these two points, then it can be
called as continuous data. On the contrary, if the data can only be defined
such that there is a finite gap on either side of it and it cannot be taken as
anything else, it would be defined as discrete.

An illustration of this is the binary of analog signal and digital signal.
Digital signals are the discrete quantities which convey information regarding
something with respect to a predefined set of values. For example, the binary
codes used in computers can only read “on” (true or 1) and “off” (false or 0).
On the other hand, an analog signal is a continuous signal with one of the
quantities changing with time (like temperature or humidity). Theoretically,
analog signals are meant to have infinitely many values and the data can take
any value within the infinitesimal range of another value. For example, a
light bulb can have not only “on” or “off” but it can be described as “dimmer”
or “brighter” on a relative scale [2].

Signals or waveforms are usually continuous data, and although Fourier
transform is an important tool to analyze them, our computers are not equipped
to perform these functions which are continuous in nature. Our computers work
in the “on/off” system, which means that they can only deal with discrete data
points. So then, how do we deal with continuous data? Well, digital signal
processing takes care of it. In digital signal processing, a physical or a
natural signal which is continuous, gets sampled and quantized [3]. The
signal’s values are picked out at equal intervals of time and only those data
points are considered essential to study the signal. The digital signal is an
ordered code obtained from a finite set of values.

Figure 1:  A digital signal picks out values which are spaced out in equal
intervals of time.

If one understands this, then one can understand how the Fourier transform and
discrete Fourier transform (DFT) would work. The Fourier transform is
essentially an integral transform; it deals with continuous data from a
waveform. However, a DFT will only pick out a specific required number of
sample points from the available data, and then analyze them and give out
crucial information, regarding the signal that we originally had [4].
Therefore, the DFT is a vital tool in signal processing and we shall explore
the use of it in later sections, once we understand the working of Fourier
transform itself.

###  II.2 Integral Transforms

The Fourier transform is basically a special case of a more general kind of
representation of functions. This general class of representations is known as
an integral transform. What is an integral transform? Well, the moment we see
the term transform, it implies some kind of a change. When we say that we are
applying a Fourier transform (or in general any integral transform) on a given
function, we mean that we are changing from one particular kind of function to
another.

Let’s say we have a function $f(t)$, which is dependent on the variable $t$,
and another function $g(x)$, which is dependent on the variable $x$. An
integral transform then gives us a way to interchange (or transform) between
these two functions, by means of an integral over a specific interval. In
terms of linear algebra, we could say the function $f(t)$ is a vector (or a
set of vectors) in the vector space $V$ with basis $t$, and the function
$g(x)$ is a vector (or a set of vectors) in the vector space $W$ with basis
$x$. In this case, and integral transform would be analogous to a linear
transformation $T:V\rightarrow W$ (or rather a change of basis from t to x),
that would tell us how to transform from $f(t)$ to $g(x)$, and vice-versa. The
vice-versa case is usually called an inverse transformation, or an inverse
integral transform. The functions $f(t)$ and $g(x)$ are related through an
integral transform as,

$g(x)=\int_{a}^{b}f(t)K(x,t)dt$ (1)

So, we get the function $g(x)$ by multiplying the function $f(t)$ by another
function $K(x,t)$ (which is dependent on both $x$ and $t$), and then
integrating this product w.r.t. $t$, over the interval $[a,b]$.Here, $a$ and
$b$ are constants, and $K$ is also called the kernel function. We can also say
that $g(x)$ itself is the integral transform of the given function $f(t)$. For
the sake of notation, it is useful to introduce a linear operator,

$g(x)=\int_{a}^{b}K(x,t)dt=T[f(t)]$ (2)

The values of $a$ and $b$, and the form of $K(x,t)$ is different for different
functions. This gives us a compact form of the integral transform as,

$T[f(t)]=\int_{a}^{b}f(t)K(x,t)dt\mid$ (3)

It turns out that the integral transform operator $(T)$ is in fact a linear
operator, as it has the following properties of linearity:

  * •

Closure under addition: The integral transform of a sum of two
functions($f_{1}$ and $f_{2}$) is equal to the sum of the integral transforms
of the two functions individually.

$T\left[f_{1}(t)+f_{2}(t)\right]=T\left[f_{1}(t)\right]+T\left[f_{2}(t)\right]$
(4)

  * •

Closure under scalar multiplication: The integral transform of a scalar
multiple of a function $f(t)$ is equal to the scalar $c$ times the integral
transform of the function.

$\begin{array}[]{c}T[cf(t)]=cT[f(t)]\end{array}$ (5)

Now, the kind of integral transforms that are useful and have varied
applications (like the Fourier transform) are generally invertible. In terms
of the linear operator notation, this means that $T^{-}1$ exists, and usually
there’s a convenient method to evaluate it. So, given that,

$g(x)=T[f(t)]$ (6)

there exists a $T^{-}1$ such that,

$T^{-1}[g(x)]=f(t)$ (7)

Apart from the Fourier transform, other kinds of integral transforms include
the Laplace transform, the Hankel transform and the Mellin transform.

Why is it that the Fourier transform (or any kind of integral transform) is
used a lot for solving both theoretical, and real-world practical problems?
Well, any kind of mathematical tool is generally introduced when solving a
problem becomes an exceedingly difficult task, or when solving the problem
becomes much easier while using the mathematical tool in a certain context.
For example, it is often easier to calculate spherical volumes (via volume
integrals) using the spherical coordinate system rather than the Cartesian
one. It is easier to use the Lagrangian or Hamiltonian formulation to analyze
certain complex dynamical systems rather than the Newtonian one [5]. It is
easier to decouple Maxwell’s inhomogeneous equations using the Lorenz gauge
than using the Coulomb gauge [6]. Symmetry principles are so often used in
physics to simplify massive calculations. The same is the case with the
Fourier transform. It is quite a sophisticated mathematical tool, which makes
seemingly difficult problems, easy to solve.

Figure 2:  The integral transform algorithm [7]

The integral transform algorithm is very well explained by Arfken, Weber and
Harris [7] using the following diagram (Figure 2), which is self-explanatory.

###  II.3 The Fourier Transform

The Fourier transform can also be viewed as an extension of the idea of a
Fourier series. The Fourier transform is to a Fourier series, what a tensor is
to a vector. In other words, the Fourier transform is a sort of generalization
of the concept of the Fourier series.

The Fourier series is basically a way to represent functions that have some
periodicity, i.e., it is a way to represent periodic functions. The important
thing to note about the Fourier series is that the periodic functions are
defined over a finite period. For a given function $f(t)$, its Fourier series
expansion is given as the sum of a constant and a series of sinusoids (sines
and cosines).

$f(t)=\frac{a_{0}}{2}+\sum_{n=1}^{\infty}\left[a_{n}\cos\left(\frac{n\pi
t}{L}\right)+b_{n}\sin\left(\frac{n\pi t}{L}\right)\right]$ (8)

So, what if we want to represent functions with an infinite period (which
basically implies an aperiodic function)? In that case, we use the Fourier
transform. We saw how an integral transform takes a function in one kind of
space (or domain) to another. The Fourier transform specifically takes a
function from the time domain to another function in the frequency domain. How
is this important or useful? Well, as we will see, this has an interesting
application in signal processing and frequency analysis of signals. For a
given function $f(t)$, its Fourier transform is given as,

$F(\omega)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^{\infty}f(t)e^{i\omega t}dt$
(9)

The only requirement or constraint on the function $f(t)$ is that the integral
of the absolute value of $f(t)$ must be a finite quantity, i.e.,
$\int_{-\infty}^{\infty}|f(t)|dt$ is finite.

###  II.4 Application of the Fourier Transform to Audio Decomposition

So far, we have discussed the basic notion of a Fourier transform and some of
its properties. However, we haven’t yet touched upon how linear algebra is
used along with the Fourier transform to perform audio decomposition
algorithms or frequency analysis of signals.

Let’s start with thinking about what an audio waveform or a signal is.
Generally, an audio file will contain a single wave with varying amplitudes
and intensities with time. Therefore, we can think of any audio file (such as
a .mp3 file or a .wav file) as intensity fluctuations in time, or simply a
function of time.

Figure 3:  An audio waveform as a combination of different frequencies.

Now, what use is a Fourier transform here? Well, a majority of the times, it
is the case that a single audio waveform is a combination of different waves
of different frequencies, rather than being a single sinusoidal wave of a
specific frequency. What if we want to find out what different kinds of
frequencies have been combined to give the final audio? The Fourier transform
comes in at this point. By applying Fourier transform on the audio waveform,
we can basically know which individual frequencies the audio is made up of.
But what is the use of this? Why would someone want to know what kind of
frequencies an audio waveform is made up of? Audio artists, or video editors
(and other people as well) find it useful to do so. Why? Well, consider that a
music composer is given an audio waveform, which has a disturbing high pitch
frequency to it, which is annoying to listen to. If the composer wants to
remove this high frequency disturbance in the audio, then he/she would first
want to know what different kinds of frequencies are individually contributing
to the final audio waveform. To do so, he/she can apply a Fourier transform to
the audio waveform, which then tells him/her what is that high frequency which
is causing an annoying disturbance; and then he/she can remove that high
frequency contribution from the audio waveform!

We still don’t know where linear algebra comes into play. Well consider this.
An average human (or even a superhuman probably) would not be able to listen
to an audio and create its waveform by hand! Computing a Fourier transform for
that waveform would be another mammoth task to do by mind. Therefore,
computational methods are required to convert an audio into its waveform, and
then to apply a Fourier transform on the waveform. The thing with computers
though, is that they are not continuous systems (like the Fourier transform),
but they work in discrete units called bits (0s and 1s). Therefore, it would
be impossible to perform the ordinary continuous Fourier transform algorithm
using a computer. This is where we introduce an incredibly important algorithm
known as the Discrete Fourier Transform (DFT) algorithm. And linear algebra
comes in handy in the DFT algorithm [8]. The fundamental process of the DFT
algorithm is that of approximating a Fourier series approximation of the audio
waveform on a finite period where the waveform is periodic. So, what we are
actually approximating is the Fourier series expansion, and not the infinite
Fourier transform. The DFT actually leads to one of the most important
algorithms - the Fast Fourier Transform (FFT) algorithm - which has
applications in image compression, audio compression, and other highly
advanced computer applications. For now, let us focus on the DFT. The reason
why linear algebra is an important ingredient of the DFT algorithm, is because
the DFT is a mathematical transformation that can be written in terms of a
huge matrix multiplication (the matrix is called the DFT matrix).

So, we talked about approximating functions that are periodic using infinite
sums of sines and cosines, using the Fourier series. But in most cases, we
don’t necessarily have an analytic (or continuous) function to deal with. In
computer related algorithms, we usually have to deal with some sort of
measurement data or with discrete data sets from a simulation or an
experiment. So, what we have is a function $f$ defined at discrete locations
$x_{0},x_{1},x_{2}...x_{n}$, as shown below in the figure.

Figure 4:  A function $f$ defined at discrete locations
$x_{0},x_{1},x_{2}...x_{n}$

These discrete data or function values $f_{0},f_{1},f_{2},...,f_{n}$ at the
respective locations $x_{1},x_{2},x_{3}...x_{n}$, can be represented as a
vector of data as follows.

$f=\left[\begin{array}[]{c}f_{0}\\\ f_{1}\\\ f_{2}\\\ \vdots\\\
f_{n}\end{array}\right]$ (10)

Now, what we want to do is compute the DFT of this vector of data, i.e., we
want to break the data up into a sum of sines and cosines. And this comes in
handy while analyzing audio. If the data is an audio data, then we can use the
DFT to figure out what sort of tones add up to make that data. Further, based
on the frequency content of the data, we can diagnose whether something is
wrong with the audio (like some unwanted deviations or spikes in the data).
So, in a nutshell, turning a data vector into its sinusoidal components using
the DFT algorithm, can be very useful. The important point to press is that a
DFT is just a Fourier series expansion of discrete data, rather than a
continuous, analytic function. Recall how while expanding a continuous
function using the Fourier series, we took the function and turned them into
the Fourier coefficients $a_{n}$ and $b_{n}$, that were multiplied to the
sines and cosines. With the DFT, we are essentially doing the same thing. We
take our data $\hat{f}$, and obtain its Fourier transform vector, , which is
given as,

$\hat{f}=\left[\begin{array}[]{c}\hat{f}_{0}\\\ \hat{f}_{1}\\\ \hat{f}_{2}\\\
\vdots\\\ \hat{f}_{n}\end{array}\right]$ (11)

The Fourier transform vector of the data, consists of frequency components of
the data (because recall that a Fourier transform takes a function from the
time domain to the frequency domain). So, in the $\hat{f}$ vector, the first
frequency component $\hat{f}_{0}$ tells us how much of that first low-
frequency is in the data $f$, the second component $\hat{f_{1}}$ tells us how
much of the second (slightly higher) frequency is there in the data, and so on
till the $n$th frequency (which is the highest possible frequency). If we let
$k=0,1,2,3,...,n$, then we can calculate the $k$th Fourier coefficient, or the
$k$th frequency component as,

$\hat{f}_{k}=\sum_{j=0}^{n-1}f_{j}e^{-\frac{i2\pi jk}{n}}$ (12)

We will get to what the exponential term means in a while, but this is
essentially how we calculate the Fourier coefficients or frequency components.
For now, it is important to keep in mind that if we have data, we can get its
Fourier transform, but if we have a Fourier transform, we should also be able
to reconstruct the original data. This sort-of inverse DFT is given as,

$f_{k}=\frac{1}{n}\sum_{j=0}^{n-1}\widehat{f}_{j}e^{\frac{i2\pi jk}{n}}$ (13)

All this basically means that if we have the data, and we apply the DFT to it,
and get the Fourier coefficients. Each of these Fourier coefficients literally
just tell us how much of each specific frequency is required to reconstruct
the original data from the coefficients. An interesting thing to note is that
all of the terms in a DFT are multiplied in the form of an exponential, which
is an integral multiple of $e^{-2\pi i/n}$. This term actually defines a
fundamental frequency, given as,

$\omega_{n}=e^{-2\pi i/n}$ (14)

This fundamental frequency is related to what kind of sines and cosines can be
approximated with $n$ discrete values in the domain $x$ (refer figure 3). In
other words, $\omega_{n}$ is the kind of fundamental frequency that we need to
work with, if we have an interval of n data points. Now, all of these Fourier
transforms are adding up integer multiples of that fundamental frequency times
the data values (and the same with the inverse Fourier transform). So, we can
use this fundamental frequency to compute a matrix, which when multiplied by
the data, would give the corresponding Fourier transform (or Fourier
coefficients).

Computing the DFT for each data point can be a tedious task by hand, leave
alone creating an entire computer program that does that! Instead, what we do
with the DFT, is that we try to write what the sum would be for each $k$, in
terms of a matrix operation. This matrix, called the DFT matrix, has
components which are in terms of the fundamental frequency $\omega_{n}$. It
can be shown that this is exactly analogous to the infinite series of sines
and cosines that we get with the ordinary Fourier series. Now, for the zeroth
frequency $(f_{0})$, then all of the exponents in the DFT would be $e^{0}$,
i.e., 1. So, all the data points will have a unit coefficient. This gives us a
linear equation of the form,

$f_{0}+f_{1}+f_{2}+\ldots+f_{n}=\hat{f}_{0}$ (15)

Similarly, for $k=1$, we get the equation,

$1f_{0}+\omega_{n}f_{1}+\omega_{n}^{2}f_{2}+\ldots+\omega_{n}^{n-1}f_{n}=\hat{f}_{1}$
(16)

If we find the equation for each particular $k$-value, then we find that we
can write these linear equations in the matrix representation as,

$\left[\begin{array}[]{c}\hat{f}_{0}\\\ \hat{f}_{1}\\\ \hat{f}_{2}\\\
\vdots\\\
\hat{f}_{n}\end{array}\right]=\left[\begin{array}[]{ccccc}1&1&1&\cdots&1\\\
1&\omega_{n}&\omega_{n}^{2}&\cdots&\omega_{n}^{n-1}\\\
1&\omega_{n}^{2}&\omega_{n}^{4}&\cdots&\omega_{n}^{2(n-1)}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
1&\omega_{n}^{n-1}&\omega_{n}^{2(n-1)}&\cdots&\omega_{n}^{(n-1)^{2}}\end{array}\right]\left[\begin{array}[]{c}f_{0}\\\
f_{1}\\\ f_{2}\\\ \vdots\\\ f_{n}\end{array}\right]$ (17)

Here, note that the fundamental frequencies are all complex numbers (which can
be easily shown using Euler’s formula). So, the DFT matrix is a complex
matrix, and all the Fourier coefficients on the LHS of the matrix equation are
complex valued.

###  II.5  Applications of the DFT matrix

Now that we have understood how Fourier Transforms and DFT work, we can delve
deeper into the applications of DFT and the DFT matrix. The Fourier Transform
vector,

$\hat{f}=\left[\begin{array}[]{c}\hat{f}_{0}\\\ \hat{f}_{1}\\\ \hat{f}_{2}\\\
\vdots\\\ \hat{f}_{n}\end{array}\right]$ (18)

which we obtained in the previous section is known as the frequency spectrum;
it tells us about the contribution of the individual frequencies of the
original waveform.

Vital information is stored in the sinusoids which form a signal. Our universe
is an encompassment of many oscillating objects, and understanding these
signals is important. For example, all the sounds that we hear are a result of
various vibrations and combinations of different sound waves; we are able to
speak because of the vibration produced in our vocal cords; there is a
periodic movement in the working of a water turbine, etc. All of these things
give out signals which are waveforms and DFT is essential to understand the
frequency, amplitude, and phase of the constituent sinusoidal waves [9].

This DFT algorithm is used vastly in our project to filter out the ingredient
frequencies of our considered audio data. This can further be used to de-noise
signals or modify them in any way. However, note that the FFT (fast Fourier
transform) would be better to use in case the data set is extremely huge,
because that’s what an FFT basically is. The FFT is just an efficient method
to calculate the DFT for large datasets.

##  III Audio Processing Samples

###  III.1 The Audio Processing Code Implemented in MATLAB

The MATLAB code file can be found here:
https://drive.google.com/file/d/10ndp1fgEiaBCc_0osNp_-
ihh_AL6ZyRt/view?usp=sharing.

The MATLAB code is designed with a user interface which is easily expandable
as one needs it. It will first take in an audio file name as an input, import
the audio file and create its waveform along with its Fourier transform. Then,
one can perform a list of tasks with it. The working of each task is explained
as comments in the program.

The code has incorporated certain ways to calculate multiplication and
transforms of matrices such that the RAM usage as well as the performance
optimized for regular use. After performing various operations on the matrix,
the user can then choose to do the inverse transform which will give back the
transformed matrix in time-dependent domain. The task ‘plots’, creates the
plot of the edited waveform as well as its Fourier counterpart that can be
observed as graphs. Additional functions to save the edited sample as well as
to listen to it (multiple times) are added for secondary observations.

###  III.2 Results from the Simulation

####  III.2.1 DFT of a Simple Chord Progression

The program takes in a sound waveform which is a superposition of different
frequencies and looks as shown below in figure 5.

Figure 5:  The amplitudes of all the frequencies which constitute the audio.

When the DFT function is applied onto it, we obtain the following graph:

Figure 6:  Frequencies of individual waves which constitute the audio wave

Once the DFT is applied on the original sound, it processes and decomposes it;
thus, yielding individual frequencies of the different waves which combine to
form the original wave. In Figure 6, we see that we obtain three different
frequencies which can be identified with their respective musical notes.

The standard note to frequency table relevant to our obtained frequencies is
as shown below.

Table 1:  Standard notes and their frequencies [10]. Note | Frequency
---|---
$C_{4}$ | 261.63
$E_{4}$ | 329.63
$G_{4}$ | 392.00

Therefore, we can obtain a rough correlation between our obtained frequencies
and the standard frequencies. We have 261.49 which pertain to the $C_{4}$
note, then 329.49 pertaining to $E_{4}$ and 391.98 pertain to $G_{4}$. As the
basic laws of music theory says, it corresponds to a $C$-major chord which
could be heard from the original signal.

####  III.2.2 DFT of Complex Music Signals

The sound wave in the previous section of composed of two different audios
whose waveforms are shown below:

  * •

Sample: Every Tear Drop is a Waterfall by Coldplay.

Figure 7:  The amplitudes of all the frequencies which constitute the audio.
Figure 8:  Frequencies of individual waves which constitute the audio wave

  * •

Sample- Master of Puppets by Metallica:

Figure 9:  The amplitudes of all the frequencies which constitute the audio.
Figure 10:  Frequencies of individual waves which constitute the audio wave

In both of the instances above, initially we have the Amplitude vs. Time graph
of the audio which shows just how many waves of different amplitudes have
superposed to create the audio that we’re analyzing. The amplitudes combine
and change as a whole in different patterns with time and create a beautiful
combination of individual waves.

Second, we have the frequency graph which is obtained once we apply the DFT on
the original sound wave. The DFT picks out the peak frequency of individual
waves within the original audio wave and plots them. Toward the left side of
the graph, the lower frequencies pertain to the bass in the music while as we
move toward the right, the higher frequencies pertain to the guitar and
contributions of other musical instruments.

This is again very interesting because it precisely shows what the main
frequencies which, in the form of sinusoids, act as building blocks to the
audio are. It also shows how many frequencies are really important to us if we
wanted to study the properties of the original audio. This operation of DFT
helps us to modify the original audio as well because if know what frequencies
are constitute the audio, we can also remove unwanted frequencies or modify
and add others (this is explained further in the next section).

###  III.3 Equalization using Eigen-matrices

Equalization is a technique in music where amplitude of specific range of
frequencies can be modified individually. For instance, in hip-hop songs, the
low frequency bass is often amplified to give a bass boost that marks the
signature music. Similarly, in Jazz and Orchestras, the middle as well as the
high frequencies (also called treble) are signified more.

Doing such changes in a time dependent mixed audio signal (especially when all
the instruments are recorded together) will be a tedious task by the process
of simple addition. One would have to add the waves of the same exact shapes
at those exact times, which would require a special solution in each case.
However, since Fourier transform is a matrix which converts the basis from
time to frequency domain, we can matrix-multiply it to more matrix and easily
define a transform that can switch frequencies. Mathematically,

$X=FAF^{-1}$ (19)

Here $A$ is a diagonal amplification matrix. Why a diagonal matrix? If we
break the above three step transform, firstly we do the Fourier transform and
get to the basis of frequencies which we simply want to scale (not shear or
rotate). Once the appropriate frequencies are scaled (we can add new
frequencies by scaling up and delete by scaling down). Then all we have to do
is do the inverse transform which will give us back the transformed signal in
time domain with the modified contribution of each frequency.

This form is the same as the Eigen matrix representation of any diagonalizable
matrix. All the complications of different shapes and quality of sound waves
(Timbre) is reduced down to sinusoids which are very easy to deal with. Each
sinusoid is scaled by the respective element in the diagonal, giving us the
equalized output.

Figure 11:  Original Sample in frequency domain. Figure 12:  Original Sample
in time domain Figure 13:  High treble transformed sample (frequency) Figure
14:  High treble transformed sample (time) Figure 15:  Bass boost transformed
sample (frequency) Figure 16:  Bass boost transformed sample (time)

As one can see from figures 11 to 16, we have chosen a sample (1.5 sec sound
wave from ‘Every Teardrop is a Waterfall’ by Coldplay) and then used the fast
Fourier transform algorithm to perform the DFT. Then all we did was define a
square diagonal matrix where each pivot represented the amplitude of the
frequency corresponding in that row. For high treble (figures 13 and 14), the
amplitudes were scaled by 0.1, 0.25, 0.5, 1, and 1, for ranges of frequencies
from low to high respectively. As one can see, the initial frequencies (around
0-160hz, in our case) are scaled down by a huge amount which is clearly
evident in the frequency domain, but less so in time domain. However, it can
be observed that the time-domain wave is denser as the resultant frequency is
higher (number of oscillations per sec will increase, making the wave closely
packed together). Therefore, desired result is achieved.

In the case of bass boost, the opposite was done, which one can clearly see in
the frequency domain where the relative peaks of low frequencies are higher
than the high frequencies compared to that of the original wave. As a result,
the time-dependent wave is less dense but has the same shape of wave hence not
changing the audio itself.

We can see as the result that the shape of the wave is maintained throughout
the transform and we do not observe any change of frequencies, yet we are able
to amplify and de-amplify specific frequencies by such a transform. The
Fourier transform did the tedious task of designing the wave while we just had
to input the amplitude for individual frequencies. Due to this the timbre and
the range of active frequencies (notes of musical instruments) is preserved,
yet their amplitudes are changed independently.

We can see that this simple algorithm can be exploited further if we introduce
shear and rotational transforms to it, therefore, giving us a lot of potential
in non-trivial signal processing.

##  IV Concluding Remarks

This paper began with the understanding of discrete and continuous data and
how continuous data can be converted to discrete data. It further explains the
Fourier transform theoretically and the working of discrete Fourier transform
in the field of Signal Processing with the help of computing. A code is built
and implemented in MATLAB programming software which can process audio samples
and yield interesting results, which opens a domain for huge variety of
functions which can be applied on an audio sample.

Hopefully, the significance of the Fourier transform in signal processing, and
more importantly, the beauty of linear algebra which aids the transform so
well, was conveyed to the reader. We can take the advantage of converting
continuous data into discrete, tangible vectors and apply these transforms as
a consequence of linear algebra.

## References

  * Fourier [1878] J. Fourier,  _The Analytical Theory of Heat_ (Cambridge University Press, 1878).
  * Horowitz [1980] P. Horowitz,  _The Art of Electronics_ (Cambridge University Press, 1980).
  * Khanna [2009] V. Khanna,  _Digital Signal Processing_ (S. Chand and Company Ltd., 2009).
  * III [2007] J. S. III,  _Mathematics of the Discrete Fourier Transform (DFT), with Audio Applications - Second Edition_ (W3K Publishing, 2007).
  * Schwichtenberg [2015] J. Schwichtenberg,  _Physics from Symmetry - Second Edition_ (Springer, 2015).
  * Griffiths [2018] D. Griffiths,  _Introduction to Electrodynamics_ (Cambridge University Press, 2018).
  * G. Arfken and Harris [2013] H. W. G. Arfken and F. Harris,  _Mathematical Methods for Physicists - Seventh Edition_ (Elsevier, 2013).
  * Brunton and Kutz [2019] S. Brunton and J. Kutz,  _Data-driven Science and Engineering_ (Cambridge Univeristy Press, 2019).
  * Bozkurt  _et al._ [2014] A. Bozkurt, M. T. Arslan, R. A. Sevimli, C. E. Akbas, and A. E. Çetin, Approximate computation of dft without performing any multiplications: Applications to radar signal processing (2014), arXiv:1402.0532 [cs.CC] .
  * Suits [1998] B. H. Suits,  _Frequencies for Equal-tempered Scale, A4 = 440 Hz_ (Online: https://pages.mtu.edu/ suits/PofMCopyRight.html, 1998).