Spectral Analysis of Saddle-point Matrices from Optimization problems with Elliptic PDE Constraints

TL;DR

This paper characterizes the spectral properties of saddle-point matrices from PDE-constrained optimization problems, revealing a GLT structure that leads to improved preconditioning strategies for iterative solvers.

Contribution

It uncovers the GLT structure in these matrices, enabling sharper spectral analysis and the development of optimal preconditioners for iterative methods.

Findings

Identification of the GLT structure in saddle-point matrices

Sharper spectral characterization of the matrices

Development of optimal preconditioners for GMRES and Flexible-GMRES

Abstract

The main focus of this paper is the characterization and exploitation of the asymptotic spectrum of the saddle--point matrix sequences arising from the discretization of optimization problems constrained by elliptic partial differential equations. We uncover the existence of a hidden structure in these matrix sequences, namely, we show that these are indeed an example of Generalized Locally Toeplitz (GLT) sequences. We show that this enables a sharper characterization of the spectral properties of such sequences than the one that is available by using only the fact that we deal with saddle--point matrices. Finally, we exploit it to propose an optimal preconditioner strategy for the GMRES, and Flexible-GMRES methods.

Tables

Table 1. Table 1: Comparison of the effective number of eigenvalues of B N subscript 𝐵 𝑁 B_{N} contained in the second interval ( m 2 , M 2 ] subscript 𝑚 2 subscript 𝑀 2 (m_{2},M_{2}] with the expected number n 2 superscript 𝑛 2 n^{2}

$n$	$\#\left\{\lambda \in (m_2,M_2]\right\}$	$n^2$	$\#\left\{\lambda \notin (m_2,M_2]\right\}$	$\#\left\{\lambda \notin (m_2,M_2]\right\}/\sqrt{3n^2}$
10	74	100	26	$0.086$
20	353	400	47	$0.039$
40	1421	1600	179	$0.037$
80	5694	6400	706	$0.036$

Table 2. Table 2: Poisson Control Problem. We compare both the number of iterations, and the solution time for the various preconditioners. Best timings are highlighted in bold face. When the method fails to converge, i.e., the method reaches the maximum number of iterations, a † † \dagger is reported. The inner tolerance for the PCG is set to 1e-8

		GMRES						FGMRES+PCG+IC
		$I_N$		${𝒫}_N$		${𝒫}_{\text{BCT}}$		${𝒫}_N$		${𝒫}_{\text{BCT}}$
$\alpha$	N	IT	T(s)	IT	T(s)	IT	T(s)	IT	T(s)	IT	T(s)
1.0e-03	147	$†$	-	3	3.0e-03	3	2.5e-03	3	4.4e-03	3	4.5e-03
	675	$†$	-	3	6.4e-03	3	3.7e-03	3	4.7e-03	3	4.6e-03
	2883	$†$	-	3	1.0e-02	3	9.9e-03	3	7.3e-03	3	7.3e-03
	11907	$†$	-	2	3.1e-02	2	3.0e-02	2	2.3e-02	2	2.3e-02
	48387	$†$	-	2	2.1e-01	2	1.7e-01	2	1.5e-01	2	1.5e-01
	195075	$†$	-	2	9.4e-01	2	9.0e-01	2	7.3e-01	2	7.4e-01
	783363	$†$	-	1	2.1e+00	1	2.0e+00	1	2.9e+00	1	2.9e+00
1.0e-06	147	$†$	-	15	4.3e-03	15	4.1e-03	15	1.5e-03	15	1.5e-03
	675	$†$	-	14	1.3e-02	14	1.2e-02	14	1.7e-02	14	1.7e-02
	2883	$†$	-	9	3.0e-02	9	3.0e-02	10	2.4e-02	10	2.4e-02
	11907	$†$	-	6	1.0e-01	6	9.3e-02	6	7.0e-02	6	7.4e-02
	48387	$†$	-	4	3.1e-01	4	3.0e-01	4	2.5e-01	4	2.1e-01
	195075	86	3.8e+00	2	8.7e-01	2	8.4e-01	2	7.8e-01	2	7.6e-01
	783363	80	3.0e+01	2	4.3e+00	2	4.3e+00	2	4.5e+00	2	4.6e+00
1.0e-09	147	$†$	-	27	9.6e-03	27	8.7e-03	27	3.2e-03	27	3.4e-03
	675	$†$	-	54	6.3e-02	54	5.8e-02	54	7.9e-02	54	7.9e-02
	2883	$†$	-	52	2.0e-01	52	2.1e-01	52	1.6e-01	52	1.6e-01
	11907	$†$	-	33	5.8e-01	33	6.1e-01	33	5.2e-01	33	5.5e-01
	48387	$†$	-	20	1.5e+00	20	1.5e+00	43	4.8e+00	42	4.8e+00
	195075	86	2.8e+00	33	1.3e+01	33	1.3e+01	37	3.0e+01	36	2.9e+01
	783363	80	3.0e+01	33	2.5e+01	33	2.5e+01	$†$	-	$†$	-

Table 3. Table 3: Poisson Control Problem. We report both the number of iterations, and the solution time for the 𝒫 D subscript 𝒫 𝐷 \mathcal{P}_{D} preconditioner in ( 5.34 ), compare these entries with the last block of rows of Table  2

	GMRES preconditioned by ${𝒫}_D$
$\alpha =$ 1.0e-09	N	147	675	2883	11907	48387	195075	783363
	IT	4	5	6	$†$	$†$	$†$	$†$
	T(s)	1.0e-02	5.4e-03	1.7e-02	-	-	-	-

Table 4. Table 4: Diffusion–Convection–Reaction Control Problem. We compare both the number of iterations, and the solution time for the various preconditioners. Best timings are highlighted in bold face. When the method fails to converge, i.e., the method reaches the maximum number of iterations, a † † \dagger is reported. The tolerances for the inner solvers are set to 1e-8

								FGMRES
		GMRES						PCG/BiCGstab+IC/ILU
		$I_N$		${𝒫}_N$		${𝒫}_{\text{BCT}}$		${𝒫}_N$		${𝒫}_{\text{BCT}}$
$\alpha$	N	IT	T(s)	IT	T(s)	IT	T(s)	IT	T(s)	IT	T(s)
1.0e-03	147	$†$	-	5	8.7e-01	5	4.7e-03	5	7.7e-03	5	7.0e-03
	675	$†$	-	5	9.6e-03	5	8.7e-03	5	7.1e-03	5	6.4e-03
	2883	$†$	-	4	7.2e-02	4	2.7e-02	4	1.4e-01	4	9.5e-03
	11907	$†$	-	3	8.6e-01	3	8.9e-02	3	4.8e-02	3	3.4e-02
	48387	$†$	-	3	1.1e+00	3	4.4e-01	3	3.9e-01	3	2.5e-01
	195075	$†$	-	2	1.7e+00	2	1.7e+00	2	1.7e+00	2	1.1e+00
	783363	$†$	-	2	8.5e+00	2	8.9e+00	2	7.1e+00	2	7.4e+00
1.0e-06	147	$†$	-	24	1.5e-02	24	1.4e-02	24	2.9e-02	24	2.4e-02
	675	$†$	-	26	4.7e-02	26	4.6e-02	26	3.6e-02	27	3.3e-02
	2883	$†$	-	24	1.7e-01	24	1.6e-01	24	7.2e-02	25	6.0e-02
	11907	$†$	-	22	6.9e-01	22	7.0e-01	22	4.0e-01	24	2.9e-01
	48387	$†$	-	19	2.8e+00	19	2.8e+00	19	2.5e+00	22	1.9e+00
	195075	$†$	-	17	1.4e+01	17	1.4e+01	17	1.4e+01	18	1.2e+01
	783363	$†$	-	14	5.9e+01	14	6.1e+01	14	7.9e+01	14	5.9e+01
1.0e-09	147	$†$	-	38	3.8e-02	38	3.5e-02	38	4.2e-02	38	4.4e-02
	675	$†$	-	73	1.5e-01	73	1.6e-01	73	1.2e-01	87	1.4e-01
	2883	$†$	-	84	6.5e-01	73	6.5e-01	86	3.3e-01	73	3.7e-01
	11907	$†$	-	94	3.5e+00	94	3.4e+00	97	2.1e+00	97	1.9e+00
	48387	$†$	-	87	1.4e+01	87	1.4e+01	87	1.2e+01	87	1.1e+01
	195075	$†$	-	77	6.8e+01	77	6.8e+01	$†$	-	$†$	-
	783363	$†$	-	66	5.2e+02	66	5.4e+02	$†$	-	$†$	-

Equations

$${\mathcal{A}}_N=\left[\begin{array}{cc}A& B_1^T\\ B_2& -C\end{array}\right],A\in {\mathbb{R}}^{q\times q},B_1,B_2\in {\mathbb{R}}^{p\times q},C\in {\mathbb{R}}^{p\times p}.$$

$$\left\{\begin{array}[]{rl}\displaystyle\min_{y,u}J(y,u)=&\displaystyle\frac{1}{2}\|y-y_{d}\|_{L^{2}(\Omega)}^{2}+\frac{\alpha}{2}\|u\|_{L^{2}(\Omega)}^{2},\\ \text{ such that }&\begin{array}[]{ll}e(y,u)=0,&\text{ in }\Omega,\\ y=f,&\text{ on }\partial\Omega_{D},\\ \frac{\partial y}{\partial\mathbf{n}}=g,&\text{ on }\partial\Omega_{N},\\ \end{array}\end{array}\right.$$

$$\mathcal{L}(y,u,p)=J(y,u)-⟨p,e(y,u){⟩}_{W^{\ast },W},$$

$${\mathcal{L}}_y^{\prime }(\hat{y},\hat{u},\hat{p})\mathfrak{h}=$$

$${\mathcal{L}}_u^{\prime }(\hat{y},\hat{u},\hat{p})\mathfrak{w}=$$

$${\mathcal{L}}_p^{\prime }(\hat{y},\hat{u},\hat{p})=$$

$$\left\{\begin{array}[]{rl}\displaystyle\min_{y,u}J(y,u)=&\displaystyle\frac{1}{2}\|y-y_{d}\|_{L^{2}(\Omega)}^{2}+\frac{\alpha}{2}\|u\|_{L^{2}(\Omega)}^{2},\\ \text{ such that }&\begin{array}[]{ll}-\nabla^{2}y=u+z,&\text{ in }\Omega,\\ y=f,&\text{ on }\partial\Omega_{D},\\ \frac{\partial y}{\partial\mathbf{n}}=g,&\text{ on }\partial\Omega_{N},\\ \end{array}\end{array}\right.$$

$$\begin{array}[]{ll}\displaystyle\left\{\begin{array}[]{ll}-\nabla^{2}y=u+z,&\text{ in }\Omega,\\ y=f,&\text{ on }\partial\Omega_{D},\\ \frac{\partial y}{\partial\mathbf{n}}=g,&\text{ on }\partial\Omega_{N}.\\ \end{array}\right.&\qquad\text{(State equation)}\\ \\ \displaystyle\left\{\begin{array}[]{ll}-\nabla^{2}p=y-y_{d},&\text{ in }\Omega,\\ y=0,&\text{ on }\partial\Omega_{D},\\ \frac{\partial y}{\partial\mathbf{n}}=0,&\text{ on }\partial\Omega_{N}.\\ \end{array}\right.&\qquad\text{(Adjoint equation)}\\ \\ \displaystyle\alpha u+p=0.&\qquad\text{(Gradient condition)}\end{array}$$

$${\int }_{\Omega }\nabla u\cdot \nabla vdx=$$

$${\int }_{\Omega }\nabla \hat{p}\cdot \nabla vdx=$$

$$\alpha {\int }_{\Omega }uvdx-{\int }_{\Omega }\hat{p}vdx=$$

$$\bar{\mathcal{A}}_{N}\mathbf{x}\equiv\left[\begin{array}[]{cc|c}\bar{M}&O&\bar{K}^{T}\\ &&\\ O&\alpha\bar{M}&-\bar{M}\\ &&\\[-1.00006pt] \hline\cr&&\\ \bar{K}&-\bar{M}&O\end{array}\right]\begin{bmatrix}\mathbf{y}\\ \\ \mathbf{u}\\ \\ \\ \mathbf{p}\end{bmatrix}=\begin{bmatrix}M\mathbf{y}_{d}\\ \\ \mathbf{0}\\ \\ \mathbf{z}\end{bmatrix}\equiv\bar{\mathbf{b}},$$

$$(\bar{M}{)}_{i,j}={\int }_{{\tau }_h}{\varphi }_i{\varphi }_{j}d\mathbf{x},(\bar{K}{)}_{i,j}={\int }_{{\tau }_h}\nabla {\varphi }_i\cdot \nabla {\varphi }_{j}d\mathbf{x},$$

$${\mathbb{P}}_p=\left\{q(x_1,x_2)=\sum _{0\le i+j\le p}{c}_{i,j}{x}_1^i{x}_2^j,c_{i,j}\in \mathbb{R}\right\}.$$

$$V_{\mathbf{n}}^p=\{v\in {\mathcal{C}}^0(\Omega )v{|}_{{\tau }_h}\in {\mathbb{P}}_p,{\tau }_h\in {\Omega }_{N(\mathbf{n})}\}\subset H^1,$$

$$V_{0,\mathbf{n}}^p=\{v\in V_{\mathbf{n}}^p,v=0\text{ on }\partial \Omega \}\subset H_0^1.$$

$$I^{-}\cup I^{+},$$

$$I^{-}=\left[\frac{1}{2}\left({\mu }_n-\sqrt{{\mu }_n^2+4{\sigma }_1^2}\right);\frac{1}{2}\left({\mu }_1-\sqrt{{\mu }_1^2+4{\sigma }_m^2}\right)\right],I^{+}=\left[{\mu }_n;\frac{1}{2}\left({\mu }_1+\sqrt{{\mu }_1^2+4{\sigma }_1^2}\right)\right].$$

$$\widehat{{\text{f}}_{\mathbf{j}}}:=\frac{1}{(2\pi {)}^d}{\int }_{{\mathcal{I}}_d}\text{f}(\mathbf{\theta }){\mathrm{e}}^{-\iota ⟨\mathbf{j},\mathbf{\theta }⟩}\mathrm{d}\mathbf{\theta }\in {\mathbb{C}}^{s\times s},\mathbf{j}=(j_1,\dots ,j_d)\in {\mathbb{Z}}^d,{\iota }^2=-1,$$

$$T_{\mathbf{n}}(\text{f})=\sum _{\mathbf{j}=-(\mathbf{n}-\mathbf{e})}^{\mathbf{n}-\mathbf{e}}{J}_{n_1}^{j_1}\otimes \cdots \otimes J_{n_d}^{j_d}\otimes \widehat{{\text{f}}_{\mathbf{j}}}.$$

$$\text{f}(\mathbf{\theta })=\sum _{\text{j}=-\infty }^{\infty }{\hat{\text{f}}}_{\text{j}}{e}^{\iota ⟨\mathbf{j},\mathbf{\theta }⟩}.$$

$$\{{\mathcal{A}}_N{\}}_{\mathbf{n}\in {\mathbb{N}}^v}{\sim }_{\lambda }(\text{f},G),$$

$$\displaystyle\lim_{{\bf n}\to\infty}\frac{1}{N}\sum_{j=1}^{N}F(\lambda_{j}(\mathcal{A}_{N}))=\frac{1}{\mu_{\ell}(G)}\int_{G}\frac{\displaystyle\sum_{i=1}^{s}F\bigg{(}\left(\lambda^{(i)}(\textbf{f})\right)(\boldsymbol{\theta})\bigg{)}}{s}{\rm d}\boldsymbol{\theta}.$$

$$\{{\mathcal{A}}_N{\}}_{\mathbf{n}\in {\mathbb{N}}^v}{\sim }_{\sigma }(\text{f},G),$$

$$\displaystyle\lim_{{\bf n}\to\infty}\frac{1}{N}\sum_{j=1}^{N}F(\sigma_{j}(\mathcal{A}_{N}))=\frac{1}{\mu_{\ell}(G)}\int_{G}\frac{\displaystyle\sum_{i=1}^{s}F\bigg{(}\left(\sigma^{(i)}(\textbf{f})\right)(\boldsymbol{\theta})\bigg{)}}{s}{\rm d}\boldsymbol{\theta}.$$

$${\left\{T_{\mathbf{n}}(f)\right\}}_{\mathbf{n}\in {\mathbb{N}}^d}{\sim }_{\lambda }(f,{\mathcal{I}}_d).$$

$$\{T_{\mathbf{n}}(\text{f}){\}}_{\mathbf{n}\in {\mathbb{N}}^d}{\sim }_{\lambda }\text{leavevmode}(\text{f},{\mathcal{I}}_d).$$

$$\text{f}(\pm {\theta }_1,\dots ,\pm {\theta }_d)\equiv \text{f}({\theta }_1,\dots ,{\theta }_d),\forall ({\theta }_1,\dots ,{\theta }_d)\in {\mathcal{I}}_d^{+}=[0,\pi {]}^d,$$

$$\{T_{\mathbf{n}}(\text{f}){\}}_{\mathbf{n}\in {\mathbb{N}}^d}{\sim }_{\lambda }\text{leavevmode}(\text{f},{\mathcal{I}}_d^{+}).$$

$$\{P_N^{-1}{A}_N{\}}_N{\sim }_{\mathrm{GLT}}{\xi }^{-1}\kappa ,\{P_N^{-1}{A}_N{\}}_N{\sim }_{\sigma ,\lambda }({\xi }^{-1}\kappa ,{\mathcal{I}}^d).$$

$${\mathcal{A}}_N={\mathcal{D}}_N^{(1)}{\bar{\mathcal{A}}}_N{\mathcal{D}}_N^{(2)}=\left[\begin{array}{ccc}{h}^{4}M& O& K^T\\ O& \alpha M& -M\\ K& -M& O\end{array}\right],h=\frac{1}{n+1},$$

$${\mathcal{D}}_N^{(1)}=\left[\begin{array}{ccc}{h}^2{I}_{n^2}& O& O\\ O& I_{n^2}& O\\ O& O& I_{n^2}\end{array}\right],{\mathcal{D}}_N^{(2)}=\left[\begin{array}{ccc}{I}_{n^2}& O& O\\ O& \frac{1}{h^2}{I}_{n^2}& O\\ O& O& \frac{1}{h^2}{I}_{n^2}\end{array}\right].$$

Full text

{AMS}

62M15, 65F08, 15B05

1 Introduction

Linear systems with saddle–point matrices arises in a wide context of applications and have attracted a great deal of attention [5, 2]. In general form, they can be simply stated as the family of linear systems where the left–hand side is given by a block–matrices of the form

[TABLE]

We are interested here in the analysis of their spectral properties in the very specific context of the discretized version of optimal constraint problems [33]

[TABLE]

where, $\alpha>0$ is a fixed constant that acts as a Tikhonov regularization parameter, $J$ is a cost functional, $\Omega\subset\mathbb{R}^{d}$ is the domain of both the state $y$ and the control $u$, and $\partial\Omega_{D}$ and $\partial\Omega_{N}$ are two disjoint sets that represent the Dirichlet and Neumann boundary respectively and have the whole boundary as union.

Spectral properties of the general case (1.1) have been indeed thoroughly analyzed [26, 22, 6, 3, 23, 18, 31, 7] under several hypotheses on the blocks of $\mathcal{A}_{N}$, e.g., $B_{1}=B_{2}=B$, $C$ semipositive definite, $A$ symmetric and positive definite, and so on. The goal of the latter works has been to provide a sharp localization bounds for their spectrum, and exploit them to devise efficient iterative solvers for such problems. Here we focus on a less general objective, i.e., we intend to exploit finer information on the structure of the blocks of (1.1), a knowledge coming from the coupling of the source problem (1.2) and its discretization, to give an asymptotic description of the spectrum of the matrices $\{\mathcal{A}_{N}\}_{N}$. Specifically, we show that the saddle–point form of $\mathcal{A}_{N}$ obtained from (1.1) hides inside another structure, namely, that the sequence of matrices $\{\mathcal{A}_{N}\}_{N}$ is a Generalized Locally Toeplitz (GLT) sequence [28, 16]. This enables us to obtain a sharper localization of its asymptotic spectrum. Furthermore, we use this characterization to suggest an effective preconditioning strategy for such problems. We stress that an approach of this type has already been exploited for both the saddle–point matrices obtained from a two–dimensional linear elasticity–type problem in [11], and partially explored in [10, 12] for a constrained optimization problem where the constraints $e(y,u)$ were Fractional Differential Equations.

The paper is therefore divided as follows, in Section 2 we describe the discrete form of (1.2) fully specifying the sequence of matrices $\{\mathcal{A}_{N}\}_{N}$. In Section 3 we recall the essential tools needed for working with GLT sequences and apply them to our problem, while in Section 4 we exploit them to devise an efficient preconditioning strategy. In Section 5 we substantiate our claims with some numerical examples, and give conclusions in Section 6.

2 From the Continuous Problem to the Saddle–point sequence $\{\mathcal{A}_{N}\}_{N}$

The first point we need to answer is how we obtain the sequence of saddle–point matrices from (1.2), indeed a way of doing so is going through its Langrangian formulation. Thus, we find the Lagrangian of (1.2) as

[TABLE]

where $e(y,u)$ represents the PDE constraint as an operator between the Banach spaces $Y\times U$ and $W$, and $p$ is the Adjoint status between the space $W$ and its dual $W^{*}$ acting as Lagrange multiplier. Indeed, a solution for the original constrained optimization problem (1.2) is a stationary point for the Lagrangian (2.3). To obtain such stationary point $(\hat{y},\hat{u},\hat{p})\in Y\times U\times W^{*}$ we require that the Gâteaux derivative with respect to each of the variables of (2.3) is zero, i.e.,

[TABLE]

These are called, in general, the first order optimality conditions or the Karush-Kuhn-Tucker conditions (KKT-conditions) for Problem (1.2). Finally, for obtaining such characterization we have to fully specify the operator $e(y,u)$, and consequently all the functional spaces $Y,U$, and $W$. The prototypical elliptic problem in this class is represented by the Poisson distributed control

[TABLE]

where $z$ represents the forcing term.

The KKT conditions for problem (2.4) are expressed as

[TABLE]

By posing $\hat{p}=-p$ and choosing $v\in H^{1}_{0}(\Omega)$ we can rewrite conditions (2.5) in weak form as:

[TABLE]

Finally, the sequence $\{\mathcal{A}_{N}\}$ is obtained by fixing a Finite Element (FEM) approximation of the optimality system (2.6). This means fixing a space $V_{0,\mathbf{n}}(\Omega_{\mathbf{n}})$ with $V_{0,\mathbf{n}}=\operatorname{Span}\{\phi_{1},\ldots,\phi_{N\mathbf{(n)}}\}\subset H^{1}_{0}(\Omega)$ over a mesh $\Omega_{\mathbf{n}}$ on the domain $\Omega$ thus obtaining the linear system

[TABLE]

where

[TABLE]

are the usual (scaled) mass and stiffness matrices, and $O$ is the zero matrix of order $N(\textbf{n})=n_{1}n_{2}\ldots n_{d}$.

2.1 Triangular Lagrangian Elements

To completely specify the linear system (2.7) we need to precise both the mesh $\Omega_{N(\mathbf{n})}$ and the basis functions $\{\phi_{j}\}_{j=1}^{N(\mathbf{n})}$, i.e., chose the element defining our discretization. We focus here on nodal Lagrangian elements [9, Chapter 5] of degree $p$. These are built starting from $\mathbb{P}_{p}$, the vector space of polynomials $q(x_{1},x_{2})$ with scalar coefficients of $\mathbb{R}^{2}$ in $\mathbb{R}$ of degree less than or equal to $p$,

[TABLE]

That is indeed a vector space of dimension $\dim\mathbb{P}_{p}=\frac{1}{2}(p+1)(p+2)$. Then an homogeneous triangulation $\Omega_{N(\mathbf{n})}$ of the unit square domain $\Omega=[0,1]^{2}$ is considered, i.e., a mesh consisting in 2D triangular cells $\tau_{h}$ with straight sides, and a lattice $\Sigma_{p}$ of nodes $\{\mathbf{N}_{i}\}_{i=1}^{\dim\mathbb{P}_{p}}$ on each triangle; see Figure 1.

By this construction, every polynomial $q\in\mathbb{P}_{p}$ is uniquely determined by its values at the points $\{\mathbf{N}_{i}\}_{i=1}^{\dim\mathbb{P}_{p}}$. The finite element method for triangular Lagrange $\mathbb{P}_{p}$ elements is then built on the discrete finite dimensional space

[TABLE]

and its subspace

[TABLE]

We call degrees of freedom of a function $v\in V_{\mathbf{n}}^{p}$ the set of the values of $v$ at the nodes $\mathbf{N}_{j}$ on the entire mesh, then the space $V_{0,\mathbf{n}}^{p}$ has exactly the dimension corresponding to the number of internal degrees of freedom, i.e., excluding the nodes on $\partial\Omega$. For our model grid we find that the degrees of freedom are $N(\mathbf{n})=n_{1}n_{2}=(pn_{x}+1)(pn_{y}+1)$, where $n_{x}$ and $n_{y}$ are the number of elements in the $x$ and $y$ direction, respectively. Thus the dimension $N$ of the matrix in (2.8) will be equal to $3N(\mathbf{n})$. The matrices (2.8) are then constructed by means of the opportune Gauss quadrature formulas, and in terms of the Lagrange basis functions $\{\phi_{i}\}_{i=1}^{N(\mathbf{n})}$. For all the discussion, and computation in the paper we deal with the matrices generated for such elements by the FEniCS library (v.2018.1.0) [1, 21].

3 Spectral analysis of the resulting sequence of saddle point matrices

This section is devoted to the attainment of a characterization of the spectra of a suitable scaling $\{{\mathcal{A}}_{N}\}_{N}$ of the sequence of matrices $\{\bar{\mathcal{A}}_{N}\}_{N}$ in (2.7). Specifically, we are going to answer to the following questions,

1. Q1

can we individuate some (possibly sharp) intervals containing the spectrum with respect to $N$? 2. Q2

For a given $N$ how many eigenvalues are in each interval? 3. Q3

What is the relation between the condition number of a suitably preconditioned matrix sequence and the value of the regularization parameter $\alpha$?

As we mentioned in the introduction, there exist classical localization results for the eigenvalues of a symmetric saddle–point matrix, like the $\mathcal{A}_{N}$ in (1.1).

Theorem 1** (Rusten and Winther [26]).**

Given $\mathcal{A}_{N}$ in (1.1), assume $A$ is symmetric and positive definite, $B_{1}=B_{2}=B$ has full rank, and $C=0$. Let $\mu_{1}$ and $\mu_{n}$ denote the largest and smallest eigenvalues of $A$, and let $\sigma_{1}$ and $\sigma_{m}$ denote the largest and smallest singular values of $B$. Then the spectrum of $\mathcal{A}_{N}$ is contained in

[TABLE]

where

[TABLE]

This bound is indeed very general and versatile, since it requires only information on the symmetry/definiteness of the diagonal blocks, and on the rank of the extradiagonal ones. It can be used to obtain an estimate of the condition number of $\mathcal{A}_{N}$ as function of $N$ in a straightforward way. To this end, an even sharper result can be obtained by means of [3, Theorem 1(c)] that permits to characterize exactly the eigenvalues with the largest and the smallest module. Nevertheless, by exploiting further information on the blocks, we show that finer answers to our question are indeed possible. Specifically, we are going to individuate three disjoint intervals $I_{0}^{-}$, $I_{1}^{+}$, and $I_{2}^{+}$ containing the spectrum of the scaled version of $\bar{\mathcal{A}}_{N}$ , we show that this choice is not arbitrary, and that it stems directly from the structure of the problem, and the selection of the discretization scheme.

In Section 3.1, we start recalling the tools we use, and then we deploy them to achieve these results in Section 3.2.

3.1 Background and definitions

Throughout this paper, we use the following notation. Let $\mathbb{C}^{s\times s}$ be the linear space of the complex $s\times s$ matrices and let $\textbf{f}:G\to\mathbb{C}^{s\times s}$, with $G\subseteq\mathbb{R}^{\ell}$, $\ell\geq 1$, measurable set. We say that f belongs to $L^{1}(G)$ (resp. is measurable) if all its components $\textit{f}_{ij}:G\to\mathbb{C},\ i,j=1,\ldots,s,$ belong to $L^{1}(G)$ (resp. are measurable). We denote by $\mathcal{I}_{d}$ the $d$-dimensional cube $(-\pi,\pi)^{d}$ and define $L^{1}(d,s)$ as the linear space of $d$-variate functions $\textbf{f}:\mathcal{I}_{d}\to\mathbb{C}^{s\times s}$, $\textbf{f}\in L^{1}(\mathcal{I}_{d})$.

Moreover we indicate by $\{\mathcal{A}_{N}\}_{{\bf n}\in\mathbb{N}^{d}}$, or simply $\{\mathcal{A}_{N}\}_{{\bf n}}$, the matrix sequence whose elements are the matrices $\mathcal{A}_{N}$ of dimensions $N\times N=N(s,\textbf{n})\times N(s,\textbf{n})$, with $N(s,\textbf{n})=sN(\textbf{n})=sn_{1}n_{2}\ldots n_{d}$, $\textbf{n}=(n_{1},n_{2},\ldots,n_{d})$.

Definition 1**.**

Let the Fourier coefficients of a given function $\textbf{f}\in L^{1}(d,s)$ be defined as

[TABLE]

where $\left\langle{\bf j},\boldsymbol{\theta}\right\rangle=\sum_{t=1}^{d}j_{t}\theta_{t}$ and the integrals in (3.9) are computed componentwise.

Then, the ${\bf n}$th Toeplitz matrix associated with f is the matrix of order $N(s,\textbf{n})$ given by

[TABLE]

where ${\bf e}=(1,\ldots,1)\in\mathbb{N}^{d},\,{\bf j}=(j_{1},\ldots,j_{d})\in\mathbb{N}^{d}$ and $J^{j_{\xi}}_{n_{\xi}}$ is the $n_{\xi}\times n_{\xi}$ matrix whose $(i,l)$th entry equals 1 if $(i-l)=j_{\xi}$ and [math] otherwise.

The set $\{T_{\bf n}(\textbf{f})\}_{{\bf n}}$ (with ${\bf n}\in\mathbb{N}^{d}$) is called the family of $d$-level Toeplitz matrices generated by f, that in turn is referred to as the generating function or the symbol of $\{T_{\bf n}(\textbf{f})\}_{{\bf n}}$.

Moreover from (3.9) the symbol can be expressed via the Fourier series

[TABLE]

In order to deal with low–rank/small–norm perturbations and to show that they do not affect the symbol of a Toeplitz sequence, we introduce the definition of spectral distribution in the sense of the eigenvalues and of the singular values for a generic matrix-sequence $\{\mathcal{A}_{N}\}_{{\bf n}\in\mathbb{N}^{v}}$, $v\geq 1$, and then the notion of GLT algebra.

Definition 2**.**

Let $\textbf{f}:G\to\mathbb{C}^{s\times s}$ be a measurable function, defined on a measurable set $G\subset\mathbb{R}^{\ell}$ with $\ell\geq 1$, $0<\mu_{\ell}(G)<\infty$. Let $\mathcal{C}_{0}(\mathbb{K})$ be the set of continuous functions with compact support over $\mathbb{K}\in\{\mathbb{C},\mathbb{R}_{0}^{+}\}$ and let $\{\mathcal{A}_{N}\}_{{\bf n}\in\mathbb{N}^{v}}$, $v\geq 1$, be a sequence of matrices with eigenvalues $\lambda_{j}(\mathcal{A}_{N})$, $j=1,\ldots,N$ and singular values $\sigma_{j}(\mathcal{A}_{N})$, $j=1,\ldots,N$.

• •

$\{\mathcal{A}_{N}\}_{{\bf n}\in\mathbb{N}^{v}}$ is distributed as the pair $(\textbf{f},G)$ in the sense of the eigenvalues, in symbols

[TABLE]

if the following limit relation holds for all $F\in\mathcal{C}_{0}(\mathbb{C})$:

[TABLE]

• •

$\{\mathcal{A}_{N}\}_{{\bf n}\in\mathbb{N}^{v}}$ is distributed as the pair $(\textbf{f},G)$ in the sense of the singular values, in symbols

[TABLE]

if the following limit relation holds for all $F\in\mathcal{C}_{0}(\mathbb{R}_{0}^{+})$:

[TABLE]

In this setting the expression ${{\bf n}\to\infty}$ means that every component of the vector ${\bf n}$ tends to infinity, that is, $\displaystyle\min_{i=1,\ldots,v}n_{i}\to\infty$.

Remark 1**.**

We denote by $\lambda^{(1)}(\textbf{f}),\ldots,\lambda^{(s)}(\textbf{f})$ and by $\sigma^{(1)}(\textbf{f}),\ldots,\sigma^{(s)}(\textbf{f})$ the eigenvalues and the singular values of a $s\times s$ matrix-valued function f, respectively. If f is smooth enough, an informal interpretation of the limit relation (3.12) (resp. (3.13)) is that when the matrix-size of $\mathcal{A}_{N}$ is sufficiently large, then $N/s$ eigenvalues (resp. singular values) of $\mathcal{A}_{N}$ can be approximated by a sampling of $\lambda^{(1)}(\textbf{f})$ (resp. $\sigma^{(1)}(\textbf{f})$) on a uniform equispaced grid of the domain $G$. Analogously each following $N/s$ eigenvalues (resp. singular values) can be approximated by an equispaced sampling of the relative $\lambda^{(j)}(\textbf{f})$ (resp. $\sigma^{(j)}(\textbf{f})$), $j=2,\ldots,s$, in the domain.

Remark 2**.**

To perform the sampling in Remark 1 computing a closed analytical expression of any of the eigenvalue functions of f is not the most effective procedure. It is costly and, essentially, useless since for $q=1,\ldots s$ we can provide an “exact” evaluation of $\lambda^{(q)}(\mathbf{f})$ at the grid points $\{\boldsymbol{\theta}_{\bf n}=(\theta_{1}^{(j)},\theta_{2}^{(k)})\}_{j,k=0}^{n-1}$ without actually computing the analytical expression. Indeed the “exact” evaluation for $d=2$ case is achieved by

sampling $\mathbf{f}$ at $\boldsymbol{\theta}_{\bf n-e}=(\theta_{n-1}^{(j)},\theta_{n-1}^{(k)})$, $j,k=0,\ldots,n-1$, and thus obtain $n^{2}$ $s\times s$ matrices, $A_{j,k},$ $j,k=0,\ldots,n-1$; 2. 2.

for each $j,k=0,\ldots,n-1$, compute the $s$ eigenvalues of $A_{j,k}$, $\lambda_{q}(A_{j,k})$, $q=1,\ldots,s$; 3. 3.

for a fixed $q=1,\ldots,s$, the evaluation of $\lambda^{(q)}(\mathbf{f})$ at $\boldsymbol{\theta}_{\bf n-e}$, $j,k=0,\ldots,n-1,$ is given by $\lambda_{q}(A_{j,k})$, $j,k=0,\ldots,n-1$.

3.1.1 Spectral analysis of Hermitian (block) Toeplitz sequences: distribution results

We collect here some classical results concerning the distribution of Hermitian (block) Toeplitz sequences from [19, 32], that we will use extensively in the following.

Theorem 2** (Grenander and Szegő [19]).**

Let $f\in L^{1}(d,1)$ be a real-valued function with $d\geq 1$. Then,

[TABLE]

In the case where f is a Hermitian matrix-valued function, according to Tilli [32], the previous theorem can be extended as follows:

Theorem 3** (Tilli [32]).**

Let $\textbf{f}\in L^{1}(d,s)$ be a Hermitian matrix-valued function with $d\geq 1,s\geq 2$. Then,

[TABLE]

Remark 3**.**

If $\{T_{\bf n}(\textbf{f})\}_{{\bf n}\in\mathbb{N}^{d}}$ is such that each $T_{\bf n}(\textbf{f})$ is symmetric with real symmetric blocks, then the symbol has the additional property that

[TABLE]

and therefore Theorem 3 can be restated as

[TABLE]

3.1.2 GLT sequences: operative features

We list here some properties and operative features from the theory of GLT sequences in their block form; refer to [29, 15, 17] for a full account of the GLT theory.

GLT1

Each GLT sequence has a singular value symbol $\textbf{f}(\textbf{x},\boldsymbol{\theta})$ for $(\textbf{x},\boldsymbol{\theta})\in[0,1]^{d}\times[-\pi,\pi]^{d}$ according to the second Item in Definition 2 with $\ell=2d$. If the sequence is Hermitian, then the distribution also holds in the eigenvalue sense. If $\{\mathcal{A}_{N}\}_{N}$ has a GLT symbol $\textbf{f}(\textbf{x},\boldsymbol{\theta})$ we will write $\{\mathcal{A}_{N}\}_{N}\sim_{\textsc{glt}}\textbf{f}(\textbf{x},\boldsymbol{\theta})$.

GLT2

The set of GLT sequences form a $*$-algebra, i.e., it is closed under linear combinations, products, inversion (whenever the symbol is singular, at most, in a set of zero Lebesgue measure), and conjugation. Hence, the sequence obtained via algebraic operations on a finite set of given GLT sequences is still a GLT sequence and its symbol is obtained by performing the same algebraic manipulations on the corresponding symbols of the input GLT sequences.

GLT3

Every Toeplitz sequence generated by an $L^{1}(d,s)$ function $\textbf{f}=\textbf{f}(\boldsymbol{\theta})$ is a GLT sequence and its symbol is f, with the specifications reported in item GLT1. We note that the function f does not depend on the space variables $\textbf{x}\in[0,1]^{d}$.

GLT4

Every sequence which is distributed as the constant zero in the singular value sense is a GLT sequence with symbol [math]. In particular:

• •

every sequence in which the rank divided by the size tends to zero, as the matrix size tends to infinity;

• •

every sequence in which the trace-norm (i.e., sum of the singular values) divided by the size tends to zero, as the matrix size tends to infinity.

GLT5

If $\{\mathcal{A}_{N}\}_{N}\sim_{\rm GLT}\kappa$ and the matrices $\mathcal{A}_{N}$ are such that $\mathcal{A}_{N}=\mathcal{X}_{N}+\mathcal{Y}_{n}$, where

• •

every $\mathcal{X}_{N}$ is Hermitian,

• •

the spectral norms of $\mathcal{X}_{N}$ and $\mathcal{Y}_{N}$ are uniformly bounded with respect to $N$,

• •

the trace-norm of $\mathcal{Y}_{N}$ divided by the matrix size $N$ converges to 0,

then the distribution holds in the eigenvalue sense.

We highlight that from the previous properties follows that a sequence of Toeplitz matrices is, up to low-rank corrections, a GLT sequence whose symbol is not affected by the low-rank perturbation.

Theorem 4**.**

[16, Section 8.4]** Let $\{A_{N}\}_{N}$ be a sequence of Hermitian matrices such that $\{A_{N}\}_{N}\sim_{GLT}\kappa$, and let $\{P_{N}\}_{N}$ be a sequence of Hermitian positive definite matrices such that $\{P_{N}\}_{N}\sim_{GLT}\xi$ and $\xi\neq 0$ a.e. Then

[TABLE]

3.2 Spectral Analysis of the Sequence $\{\mathcal{A}_{N}\}_{N}$

We can now use the introduced tools to perform the spectral analysis of the matrix sequence $\{\bar{\mathcal{A}}_{N}\}_{N}$ , assuming that $n=n_{1}=n_{2}$, $p=1$. For studying it is easier to consider the equivalent distribution given by the following symmetric diagonal scaling

[TABLE]

with

[TABLE]

From the discretization of the Section 2, the elements of the matrix $\bar{M}$ depend on $n$ as $1/(n+1)^{2}$. Hence, the effect of the proposed scaling permits to eliminate the dependence of $h^{2}$ of the elements in $\bar{M}$, which, for $n$ large, would make the matrix $\mathcal{A}_{N}$ ill-conditioned.

In particular the matrices ${M}=\frac{1}{h^{2}}\bar{M}=T_{\mathbf{n}}(m)$, ${K}=\bar{K}=T_{\mathbf{n}}(\kappa)$ are $n^{2}\times n^{2}$ bi-level Toeplitz matrices with generating functions

[TABLE]

and

[TABLE]

We stress that in this case the matrices $M$ and $K$ are real and symmetric. A property that we will exploit the theoretical analysis, nevertheless we keep the notation $K^{T}$ for the (1,3) block of the matrix $\mathcal{A}_{N}$ for two reasons. On one side, for being consistent with the continuous setting, in which the adjoint is usually explicitly expressed. On the other, to keep the analogy with Section 3.3 in which we will discuss the usage of the advection-diffusion equation as constraint.

Theorem 5**.**

The matrix sequence $\{\mathcal{A}_{N}\}_{N}$ in (3.14) is distributed in the sense of the Eigenvalues as

[TABLE]

i.e., $\{\mathcal{A}_{N}\}_{N}\sim_{\lambda}(\textbf{f},[0,\pi]^{2})$, where

[TABLE]

Proof.

Let $\mathbf{e}_{i}$, $i=1,\ldots,N$ be the $i$th column of the identity matrix of size $N$, we can define a proper $N\times N$ permutation matrix, $\Pi=[P_{1}|P_{2}|P_{3}]$, $P_{l}\in\mathbb{R}^{N\times n^{2}},\,l=1,2,3$, such that the $k$th column of $P_{l}\,l=1,2,3$, is $e_{l+3(k-1)}$. The matrix $\Pi$ transforms $\mathcal{A}_{N}$ as

[TABLE]

where

• •

$T_{\textbf{n}}(\textbf{f})$ is the bi-level $3\times 3$ block Toeplitz $T_{\bf{n}}(\textbf{f})=\left[\hat{\textbf{f}}_{\textbf{i}-\textbf{j}}\right]_{\textbf{i},\textbf{j}=\textbf{e}}^{\bf n}\in\mathbb{C}^{N\times N}$ generated by $\textbf{f}:[-\pi,\pi]^{2}\rightarrow\mathbb{C}^{3\times 3}$ as in (3.11),

• •

$E_{\textbf{n}}$ is a small-norm matrix, with $||E_{\textbf{n}}||<C,$ $C$ constant depending on the bandwidths of $B_{N}$ and $N^{-1}\|E_{\textbf{n}}\|_{1}\to 0$.

This is a congruence transformation, thus if we find the distribution of the sequence $\{B_{N}\}_{N}$ we found also the distribution for the sequence $\{\mathcal{A}_{N}\}_{N}$. Let us observe that the nonzero entries of $T_{\bf n}(\textbf{f})=[\hat{\textbf{f}}_{\textbf{i}-\textbf{j}}]_{\textbf{i},\textbf{j}=\textbf{e}}^{\textbf{n}}$ correspond to the indexes $\textbf{i}=(i_{1},i_{2}),\textbf{j}=(j_{1},j_{2})$ satisfying

[TABLE]

as shown in equation (3.20), for $\textbf{n}=(3,3)$ we find $T_{\bf n}(\textbf{f})$

[TABLE]

Therefore, from (3.11), the generating function f is given by the finite sum

[TABLE]

where $\hat{\textbf{f}}_{(0,0)},\hat{\textbf{f}}_{(-1,0)},\hat{\textbf{f}}_{(0,-1)},\hat{\textbf{f}}_{(1,0)},\hat{\textbf{f}}_{(0,1)},\hat{\textbf{f}}_{(1,1)},\hat{\textbf{f}}_{(-1,-1)}\in\mathbb{R}^{3\times 3}$, that is f is a linear trigonometric polynomial in the variables $\theta_{1}$ and $\theta_{2}$ with matrix coefficients from (3.18). Moreover, using the equalities in (3.18), the symbol in (3.21) can be readily simplified as

[TABLE]

Note, from the latter, that

[TABLE]

thus f is a symmetric matrix-valued function which implies that $T_{\bf n}(\textbf{f})$ is a symmetric matrix. By Theorem 3, we conclude that

[TABLE]

While, from GLT3, we know that $\{T_{\bf n}(\textbf{f})\}_{{\bf n}}$ is a GLT sequence with symbol f. Moreover, let us observe that $\{E_{\bf n}\}$ is a zero–distributed sequence hence $\{E_{\bf n}\}_{{\bf n}}\sim_{\sigma}(\textbf{0},\mathcal{I}_{2}^{+})$. Indeed, $E_{\bf n}$ is the permutation of a matrix that in block position (1,1) collects all the terms that contains the scaling $h^{4}$, deriving from the (1,1) block of $\mathcal{A}_{N}$, and [math] anywhere else. Then it can be written as $E_{\bf n}=h^{4}\tilde{E}_{\bf n}.$

Since the trace norm $\|\cdot\|_{1}$ of $\tilde{E}_{\bf n}$ is equal to a constant $C$ independent on ${\bf n}$, we have

[TABLE]

and hence the zero–distribution follows from GLT4. In addition, from GLT1 and the fact that $E_{\bf n}$ is Hermitian, $\{E_{\bf n}\}_{{\bf n}}\sim_{\lambda}(\textbf{0},\mathcal{I}_{2}^{+})$.

The conclusion of the Theorem is then achieved by applying GLT2 and (3.22), since this proves that $\{T_{\bf n}(\textbf{f})+E_{\bf n}\}_{{\bf n}\in\mathbb{N}^{2}}$ is a GLT sequence with symbol $\mathbf{f}$, i.e., $\{\mathcal{A}_{N}\}_{N}\sim_{\rm GLT}\textbf{f}$. Consequently, by recalling that $T_{\bf n}(\textbf{f})+E_{\bf n}$ is real symmetric for every $\bf n$ and using GLT1, we deduce that the distribution result holds in the sense of the eigenvalues

[TABLE]

Furthermore, since each $B_{N}$ is symmetric and its blocks are symmetric and real, then f is such that $\textbf{f}(\pm\theta_{1},\pm\theta_{2})\equiv\textbf{f}(\theta_{1},\theta_{2})$, $\forall(\theta_{1},\theta_{2})\in[0,\pi]^{2}$ and therefore (3.23) can be rephrased as

[TABLE]

We can now find a first answer to the questions Q1 and Q2. For $N$ sufficiently large, let

[TABLE]

be the eigenvalues of $B_{N}$ from (3.19), i.e., of $\mathcal{A}_{N}$. By Remark 1, with $s=3$, and equation (3.24), we discover that $N/3=n^{2}$ eigenvalues of $B_{N}$, up to a number of outliers infinitesimal in the dimension, can be approximated by a sampling of $\lambda^{(1)}(\textbf{f})$ on an opportune grid (see the following discussion). The next $N/3$ on the second one and the last $n^{2}$ on the sampling of $\lambda^{(3)}(\textbf{f})$. Moreover, obtaining the following proposition, as a specialized version of Theorem 1, is straightforward.

Proposition 1**.**

Let $m_{i}=\operatorname*{ess\,inf}_{\mathcal{I}^{+}_{2}}\lambda^{(i)}(\textbf{f}(\boldsymbol{\theta}))$ and $M_{i}=\operatorname*{ess\,sup}_{\mathcal{I}^{+}_{2}}\lambda^{(i)}(\textbf{f}(\boldsymbol{\theta}))$ be the essential infimum and essential supremum of $\lambda^{(i)}(\textbf{f}(\boldsymbol{\theta}))$ respectively, for $i=1,2,3$. Then, for $N$ sufficiently large, the spectrum $\lambda(\mathcal{A}_{N})$ of the matrix sequence $\{\mathcal{A}_{N}\}_{N}$ is contained in three intervals

[TABLE]

*for $\mathcal{I}^{+}_{2}=[0,\pi]^{2}$. *

Proof.

From the definition of f in (3.17), $\forall\,(\theta_{1},\theta_{2})\in[0,\pi]^{2}$, and matching with the classical analysis for saddle–point matrices in Theorem 1, we find

[TABLE]

i.e.,

[TABLE]

and

[TABLE]

From  [27, Theorem 2.3], we know that the thesis holds true for $T_{\bf n}(\textbf{f})$ and, from the relation $\{\mathcal{A}_{N}\}_{N}\sim_{\lambda}(\textbf{f},[0,\pi]^{2})$ of Theorem 5, we have that asymptotically the inclusion in (3.27) is valid, also involving the small norm correction.

To deliver an actual numerical estimate for these bounds what we need is a reasonable approximation of the eigenvalue functions $\lambda^{(l)}({\bf f})$, $l=1,2,3$, following the procedure from Remark  2 and exploiting Theorem 5, we define the following equispaced grid on $\mathcal{I}^{+}_{2}$

[TABLE]

and consider the following $n^{2}$ Hermitian matrices of size $3\times 3$

[TABLE]

Ordering in ascending way the eigenvalues of $A_{j,k}$

[TABLE]

for any $l=1,2,3$, an evaluation of $\lambda^{(l)}(\mathbf{f})$ at $(\theta_{1}^{(j)},\theta_{2}^{(k)})$ is given by $\lambda_{l}(A_{j,k})$, $j,k=1,\ldots,n$. For a fixed $l$, we denote the vector of all eigenvalues $\lambda_{l}(A_{j,k})$, $j,k=0,\ldots,n-1$ as $\mathbf{P}^{(n)}_{l}$ , i.e.,

[TABLE]

and by $\mathbf{P}^{(n)}$ the vector of all eigenvalues $\lambda_{l}(A_{j,k})$, $j,k=0,\ldots,n-1$ varying $l$, i.e.,

[TABLE]

Note that, refining the grid by increasing $n$, we can provide the evaluation of the eigenvalue functions of f in a larger number of grid points: numerical evidences of this fact are reported in Figure 2,

in which we compare the approximation of $\lambda^{(l)}(\textbf{f})$ on $\boldsymbol{\theta}_{\bf n}$, $n=5,6$ contained in $\mathbf{P}^{(n)}_{l}$ (ordered in ascending way) with the approximation of the same eigenvalue function on a grid that is twice as fine $\boldsymbol{\theta}_{\bf 2n-e}$, $n=5,6$ contained in $\mathbf{P}^{(2n)}_{l}$ (ordered in ascending way as well) for every $l=1,2,3$.

Then, for $n$ sufficiently large, if we order in ascending way $\mathbf{P}^{(n)}_{l}$, its extremes satisfy the following relations

[TABLE]

and we can can compute a satisfactory approximation of the $\{m_{l},M_{l}\}_{l=1}^{3}$ from Proposition 1, e.g., by setting $n=3\cdot 10^{3}$, and $\alpha=\texttt{1.0e-04}$, we obtain the following approximations

[TABLE]

This clearly matches with the fact that the matrix–valued symbol is analytically singular in $(0,0)$, i.e.,

[TABLE]

hence $m_{2}=0$, nevertheless we stress again that this is not in contradiction with the fact that $\mathcal{A}_{N}$ is non singular.

In conclusion, we can exploit Remark 1, to provide an answer to Q2 determining how many eigenvalues are asymptotically contained in each of the three blocks. According to the relations (3.24), (3.26) we expect the eigenvalues of $B_{N}$ to verify

[TABLE]

and then to identify $3$ blocks

[TABLE]

Correspondingly, we can split the vector $\mathbf{P}^{(n)}$ containing the sampling of the eigenvalue functions on $\boldsymbol{\theta}_{\bf n-e}$ as follows

[TABLE]

We stress again that (3.29) allows for a number of outliers that is infinitesimal in the dimension $N$.

For example, for ${\bf n}=(n,n)=(40,40)$ ($N=4800$), approximately $\frac{3n^{2}}{3}=1600$ eigenvalues should be in each block, by a straightforward numerical check one obtains

[TABLE]

Therefore, we expect from that a certain number of eigenvalues of $B_{N}$ are in none of the blocks; in the example the effective $1421$ eigenvalues against the expected $1600$ in the second block. This is confirmed again by Figure 3 in which we highlight represent in blue the whole spectrum of $B_{N}$ and highlight in black the outliers not belonging to the blocks.

On the other hand, such a phenomenon is in line with (3.29), since the order of what is missing/exceeding is infinitesimal in the dimension $N$. As an example, in Table 1 we compare the actual number of eigenvalues of $B_{N}$ contained in the second interval $(m_{2},M_{2}]$ with the expected number $n^{2}$. In such way, we succeed in counting the outliers of $B_{N}$ in $(m_{2},M_{2}]$, whose cardinality behaves as $O(\sqrt{3n^{2}})$.

A further and more natural evidence of relation (3.24) can be obtained by comparing block by block the eigenvalues of $B_{N}$ with the sampling of the eigenvalue functions of f, that is comparing Bl1, Bl2, Bl3, with Eval1, Eval2, Eval3, respectively. Indeed we want to compare the eigenvalues of $B_{N}$ (properly ordered) with the evaluation of $\lambda^{(l)}(\mathbf{f})$ $l=1,2,3$ at $\boldsymbol{\theta}_{\bf n-e}$, using the values that are present in the blocks of $\mathbf{P}^{(n)}$.

More precisely, we compare the elements of Evalt with the elements of Blt by means of the following matching algorithm:

• •

save the couples $(\theta_{n-1}^{(j_{t})},\theta_{n-1}^{(k_{t})})$ of $\boldsymbol{\theta}_{\bf n-e}$ to which the elements of Evalt are associated with;

• •

for a fixed $\lambda\in{\rm Bl}_{t}$ find $\tilde{\eta}\in{\rm Eval}_{t}$ such that

[TABLE]

• •

associate $\lambda$ to the couple $(\theta_{n-1}^{(j_{t})},\theta_{n-1}^{(k_{t})})$ corresponding to $\tilde{\eta}$.

Making use of the previous algorithm, in Figure 4, we compare the eigenvalues of $B_{N}$ with $\lambda^{(l)}(\mathbf{f})$, $l=1,2,3$ displayed as a mesh on $\boldsymbol{\theta}_{\bf n-e}$, for $n=40$. The eigenvalues of $B_{N}$ mimic, up to some outliers shown in the Figure 4b, the sampling of the eigenvalue functions, numerically confirming the result given in Theorem 5.

3.3 From Poisson to advection-diffusion equations

We have built the whole construction using as constraint the Poisson differential equation, this is not restrictive since the analysis can be transparently extended to encompass constraints given by a generic elliptic differential equations, i.e.,

[TABLE]

The matrix sequence (2.7) maintains the same $3\times 3$ block structure, but with a different (1,3) and (3,1) block $\bar{Z}$. The latter, whenever $\mathbf{c}=(c_{1},c_{2})\neq 0$, is no more symmetric since the new constraint is no more self–adjoint. Specifically, the new block $\bar{Z}$ can be decomposed into the sum of three terms,

[TABLE]

with $V\neq V^{T}$. Therefore, the relative scaled version is given by

[TABLE]

By means of a GLT perturbation argument from Section 3.1, and exploiting the analysis in [17, Section 7.4] for the presence of lower order differential terms, we can obtain again a characterization of the eigenvalues of $\mathcal{S}_{N}$ in (3.32) that is analogous to the one we gave in Theorem 5.

Proposition 2**.**

*The matrix sequence $\{\mathcal{S}_{N}\}_{N}$ from (3.32) is distributed in the eigenvalue sense as the matrix–valued function $\mathbf{f}$ from Theorem 5. *

Proof.

Follows from Theorem 5, the techniques adopted in its proof, and from GLT5 applied to $\mathcal{S}_{N}=\mathcal{A}_{N}+\mathcal{Y}_{N}$, where

[TABLE]

4 An optimal preconditioning strategy

In this section we analyze an effective procedure to precondition the GMRES method for the solution of the systems (3.14), and (3.32). There exist indeed many preconditioners for the linear systems of saddle–point type exploiting their block structure, see, e.g, the review [5] the comparisons in [2], and, more specifically, the approaches described in [4, 24, 25, 20]. What we present here belongs to this class, and is built with the objective of obtaining algorithmic scalability, i.e., independence of the number of iteration from $h$, and optimality with respect to the parameter $\alpha$, i.e., independence of the number of iteration also with respect to it. To achieve this kind of results the classical techniques can be broadly divided into three classes, the case of definite Hermitian preconditioners for which it is possible to retrieve a cluster of the eigenvalue sense from a cluster of the singular values [30, 24, 4], that allows also for the use of the MINRES method; the case of the indefinite Hermitian preconditioners, and non Hermitian preconditioner [25, 20]. We focus here on the last approach, while benefiting both from the spectral distribution of the sequence $\{T_{\bf n}(m)\}_{\bf n}$ and $\{T_{\bf n}(\kappa)\}_{\bf n}$ of the Sections 3.2, 3.3, and from the block form of the matrices $\mathcal{A}_{N}$ and $\mathcal{S}_{N}$. Specifically, we propose the following preconditioner

[TABLE]

This is clearly an indefinite, and non Hermitian matrix, nevertheless, the linear systems involving it can be easily solved by the following back–substitution procedure:

Solve $\alpha K^{T}{\mathbf{z}}_{2}=\mathbf{r}_{1}$; 2. 2.

Solve $M{\mathbf{z}}_{3}=\alpha M{\mathbf{z}}_{2}-{\mathbf{r}}_{2}$; 3. 3.

Solve $K{\mathbf{z}}_{1}={\mathbf{r}}_{3}+M{\mathbf{z}}_{2}$.

We stress that this does not require the approximation of any of the possible Schur complements of $\mathcal{A}_{N}$ ($\mathcal{S}_{N}$), thus greatly simplifying the construction of the preconditioner. Moreover, we are going to prove now that this choice provides a strong cluster at 1 for the eigenvalues of the preconditioned linear system while obtaining also the independence of $\alpha$. We obtain this result in two steps by means of the GLT theory showing that the matrix sequence $\{\mathcal{P}_{N}^{-1}\mathcal{A}_{N}\}_{N}$ is distributed in the sense of the eigenvalues as 1. First, in Proposition 3, we show that the eigenvalues of the preconditioned matrix $\mathcal{P}_{N}{{}^{-1}}\mathcal{A}_{N}$ are either $1$, or the generalized eigenvalues of an auxiliary problem, then, in Lemma 1, we prove that the matrix sequence associated to the latter is indeed distributed in the eigenvalue sense as the function 1, thus obtaining that the eigenvalues of the preconditioned system are strictly clustered at $1$.

Proposition 3**.**

Let $\mathcal{A}_{N}$ ($\mathcal{S}_{N}$) be the coefficient matrix in (3.14) (respectively in (3.32)), and let $\mathcal{P}_{N}$ be the associated preconditioner from (4.33). Then, the eigenvalues of the preconditioned matrix $\mathcal{P}_{N}^{-1}\mathcal{A}_{N}$ are

• •

$\lambda_{j}=1$* for $j=1,\ldots,2N(\mathbf{n})$,*

• •

$\lambda_{j}$* for $j=2N(\mathbf{n})+1,\ldots,N(3,\mathbf{n})$ given by the solution of the generalized eigenvalue problem*

[TABLE]

with $\textbf{x}_{1}\neq\textbf{0}\in{\mathbb{R}^{N(\mathbf{n})}}.$

Proof.

For each $n$, $\lambda$ is an eigenvalue of the matrix $\mathcal{P}^{-1}_{N}\mathcal{A}_{N}$ if $(\lambda,\bf{x})$ is an eigenpair of the eigenvalue problem

[TABLE]

with

[TABLE]

That is $(\lambda,\bf{x})$ is solution of

[TABLE]

It is clear from the second and the third “block” equations that $(1,{\bf x})$ is an eigenpair for the latter problem for all the vectors in the $N(2,\mathbf{n})$ subspace of $\mathbb{R}^{N(3,\mathbf{n})}$

[TABLE]

Otherwise, if $\lambda\neq 1$, from the third “block” equation

[TABLE]

follows

[TABLE]

And thus, by substitution, we easily find

[TABLE]

and thus the remaining eigenpairs are given by the solution of

[TABLE]

Lemma 1**.**

The matrix sequence

[TABLE]

associated to the generalized eigenvalue problem

[TABLE]

*is distributed in the eigenvalue sense as ${\bf 1}$ over $\mathcal{I}_{2}^{+}$. *

Proof.

The statement is equivalent to

[TABLE]

since, from (3.15) and (3.16), we have that $M$ and $K$ are the symmetric and positive definite matrices $T_{\bf n}(m)$ and $T_{\bf n}(\kappa)$, respectively.

Moreover the sequence \biggl{\{}\frac{h^{4}}{\alpha}T_{\bf n}(m)\biggr{\}}_{\bf n} is distribuited in the singular value sense as [math] over $\mathcal{I}_{2}^{+}$. Hence from property GLT4 plus properties GLT2-GLT3 we have that the following GLT results hold:

[TABLE]

and

[TABLE]

Exploiting again GLT 2–GLT4 we obtain that

[TABLE]

and

[TABLE]

Since the matrix $T_{\bf n}(\kappa)T_{\bf n}^{-1}(m)T_{\bf n}(\kappa)$ is positive definite, then Theorem 4 implies

[TABLE]

and, hence, the thesis.

Remark 4**.**

Let us stress that the conclusion in Lemma 1 is again an asymptotic result for $h\rightarrow 0$ that is then valid for a fixed value of the parameter $\alpha$. Furthermore, it permits also an answer to Q3 characterizing the condition number of the preconditioned matrix sequence. Specifically, if we let $X$ be the matrix of the generalized eigenvectors for the pencil $(K,M)$, i.e., if $X$ is an invertible matrix such that

[TABLE]

then we find

[TABLE]

It is then straightforward to use (3.15) and (3.16) to estimate the maximum eigenvalues of the generalized eigenvalue problem in Proposition 3 as an $O(\alpha^{-1})$. This means that the asymptotic regime described in Lemma 1 is evident whenever $h^{4}$ becomes smaller than the fixed value of $\alpha$ of the given problem.

We can now answer to question Q3 for both the matrix sequences $\{\mathcal{P}_{N}^{-1}\mathcal{A}_{N}\}_{N}$, and $\{\mathcal{P}_{N}^{-1}S_{N}\}_{N}$ of the Subsection 3.3, where in the definition of the preconditioner (4.33) $Z$ plays the same role of $K$.

Theorem 6**.**

*The matrix sequences $\{\mathcal{P}_{N}^{-1}\mathcal{A}_{N}\}_{N}\sim_{\lambda}(\mathbf{1},\mathcal{I}_{2}^{+})$, $\{\mathcal{P}_{N}^{-1}\mathcal{S}_{N}\}_{N}\sim_{\lambda}(\mathbf{1},\mathcal{I}_{2}^{+})$ independently of $\alpha$. *

Moreover, an analogous spectral result to Theorem 6 can be given for the sequence $\{\mathcal{P}_{\text{BCT}}^{-1}\mathcal{A}_{N}\}_{N}$ (respectively, $\{\mathcal{P}_{\text{BCT}}^{-1}\mathcal{S}_{N}\}_{N}$), for

[TABLE]

Theorem 7**.**

*The matrix sequences $\{\mathcal{P}_{\text{BCT}}^{-1}\mathcal{A}_{N}\}_{N}\sim_{\lambda}(\mathbf{1},\mathcal{I}_{2}^{+})$, $\{\mathcal{P}_{\text{BCT}}^{-1}\mathcal{S}_{N}\}_{N}\sim_{\lambda}(\mathbf{1},\mathcal{I}_{2}^{+})$ independently of $\alpha$. *

Proof.

The proof follows the proofs of the Proposition 3 and Lemma 1, replacing the expression of $\mathcal{P}_{N}$ with that of $\mathcal{P}_{\text{BCT}}$.

This is indeed an example of a block–counter–triangular preconditioner in the style of [4].

Remark 5**.**

The preconditioner proposed in [4] takes the lower anti–triangular part of a different permutation of the system matrix $\mathcal{A}_{N}$, and considers also a different scaling. By this approach, the term that is dropped out in the preconditioner is not a correction of “small” norm, and this makes a substantial difference in the performances of the two approaches. Specifically, comparing the results of Proposition 3, with [4, Theorem 3.1], it is straightforward to observe that in the latter case it is not possible to infer a cluster of the eigenvalues of the preconditioned system, specifically, for the rearranged system

[TABLE]

The non-unit eigenvalues are the one of the matrix sequence $\{I+\alpha h^{-4}M^{-1}KM^{-1}K^{T}\}_{N}$, for which the clustering at one cannot be concluded. Similar observation can be made also for the null–space based block anti–triangular preconditioners [24] arising from the block anti–triangular factorization of the saddle–point matrix. Furthermore, one could consider the preconditioner which neglects the (3,2) block of $\mathcal{\bar{A}}_{N}$, avoiding the reordering and the scaling. This would bring to the case where the non-unit eigenvalues are the solution of the following generalized eigenvalue problem

[TABLE]

and, then, we have a behavior analogous to the case with preconditioner $\tilde{\mathcal{P}}_{\text{BCT}}$ (i.e., the absence of a provable cluster of the preconditioned sequence). Precisely, the non-unit eigenvalues are of the form $\lambda_{i}=1+\mu_{i}$, where, $\mu_{i}$ are the reciprocal of the eigenvalues of the matrix sequence $\{\frac{\alpha}{h^{4}}M^{-1}KM^{-1}K^{T}\}_{N}$.

4.1 Approximate iterative solution of the auxiliary linear systems

The application of the proposed preconditioners requires the solution of auxiliary linear systems with the matrices $K$ ,$K^{T}$, and $M$ or, respectively, $Z$, $Z^{T}$, and $M$ obtained from (4.33). In both cases we are dealing with very common linear systems for which there exist highly efficient and specific solvers, e.g., fast Poisson solvers, multigrid methods of geometric, and algebraic type, inner–outer Krylov solver with incomplete factorization preconditioner, and several combinations of all the previous. Potentially, any optimal preconditioner for these matrices could be included in the present framework without spoiling the overall construction, the actual choice is indeed a matter of computational framework; see, e.g., [8, Chapter 3.8]. For the solution of the systems involving the mass matrix $M$ a straightforward solution is using the unpreconditioned CG method or its preconditioned version. In the latter case, we use either a modified incomplete Cholesky factorization with drop–tolerance 1e-2 or a standard algebraic multigrid. We stress that the solution of the system involving the stiffness matrix can be machine-dependent; see, e.g., Figure 5. We easily observe that the fastest solution with the required accuracy for the system involving the $K=T_{\mathbf{n}}(k)$ is obtained by using the PCG with a standard AMG preconditioner. On the other hand, for the non symmetric case we can use the BiCGstab method together with a modified incomplete LU factorization of Crout type.

Nevertheless, as we discuss in the next Section 5, the time–efficiency in the auxiliary solve it is not so crucial, observe that already the direct method gives acceptable results under this aspect. What really matters is the combination of the achieved accuracy of the auxiliary solve with the presence, and the possible accumulation, of the $\alpha$ factor in the right–hand side of the auxiliary linear systems. This will cause for their solution by a direct method to return better performances for the lowest value of $\alpha$.

5 Numerical Examples

In this section we test the application of the preconditioners analyzed in Section 4 on some test problems. All the numerical tests are made on a laptop running Linux with 8 Gb memory and CPU Intel® Core™ i7–4710HQ CPU with clock 2.50 GHz and MATLAB version 9.4.0.813654 (R2018a). We recall again that all the relevant matrices and right–hand sides are generated by means of the FEniCS library (v.2018.1.0) [1, 21]; see again Section 2 for the details.

We test the solution procedure with the un–restarted GMRES method set to achieve a tolerance on the residual of tol = 1e-6, and a maximum number of iteration maxit = 100, and measure the number of iterations, and the timings in second. As test problem we consider an instance of a Poisson control problem (2.4), and one with the diffusion–advection–reaction constraint from Section 3.3.

Poisson

The first test problem is an instance of the Poisson control problem (2.4), in which we want to obtain the desired state,

[TABLE]

while using the forcing term

[TABLE]

We test the solution for regularization parameter $\alpha=\texttt{1.0e-03},\texttt{1.0e-06},\texttt{1.0e-09}$, and collect the results in Table 2. The approximate preconditioners are applied inside the Flexible–GMRES method as discussed in Section 4.1.

What we observe is that the approximate solution are at an advantage for the higher value of $\alpha$, while perform poorly for the smallest $\alpha=\texttt{1.0e-09}$. We stress that this effect is more connected to the behavior of the accuracy in the computation of the Krylov vectors inside the FGMRES method, than to the optimal behavior of the auxiliary problems. Secondarily, what we observe is indeed the optimal behavior with respect to the iteration discussed in Theorem 7. Indeed, the preconditioning routine becomes asymptotically better with the size of the problem, i.e., we get fewer iteration for bigger problems. Moreover, the decreasing of the $\alpha$ introduces just a latency effect in the solution, i.e., the asymptotic regimes kicks in for slightly bigger problems when $\alpha$ is smaller, we stress that this is exactly the phenomenon described in Remark 4 regarding the asymptotic relation between the value of $h$ going to zero, and the value of $\alpha$ being fixed independently of $h$. To overcome this limitation, one could decouple the system by neglecting the matrix $\alpha\bar{M}^{-1}$, i.e., the $(2,2)$ block in (2.7), thus obtaining the preconditioner

[TABLE]

By computation analogous to the one in Remark 5, we find that the non-unit eigenvalues for this preconditioner are the ones of the matrix sequence $\left\{I+\frac{\alpha}{h^{4}}M^{-1}KM^{-1}K^{T}\right\}_{N}$. The non-unit eigenvalues tend to cluster at one whenever $\alpha h^{-4}\propto\alpha N^{4}$ goes to zero. This means that $\mathcal{P}_{D}$ is efficient for small values of $\alpha$ and moderate values of $N$ and worsen for diverging values of $N$ (keeping fixed $\alpha$), indeed this is confirmed by the numerical test in Table 3.

Diffusion–Convection–Reaction

The second case we consider is the problem (1.2) in which the costraint $e(y,u)$ is given by the Equation (3.31), with coefficients $r=1$, and $\mathbf{c}=(2,3)$. The desired state is given by the sum of the two impulses

[TABLE]

while the forcing term is given by

[TABLE]

We test the solution for regularization parameter $\alpha=\texttt{1.0e-03},\texttt{1.0e-06},\texttt{1.0e-09}$, and collect the results in Table 4.

The results are completely analogous to the one for the Poisson case. We observe a higher number of iteration that is due to the fact that we are using an asymptotic argument both for the sequence $\mathcal{S}_{N}$, and for its block; see Proposition 2, and the discussion in Remark 4 for the asymptotic relationship between $h$, and $\alpha$.

6 Conclusions and future developments

In this paper we have produced a characterization for the saddle–point matrices arising from the application of the discretize–then–optimize approach to quadratic optimization problems with elliptic PDE constraints highlighting the presence of an hidden Generalized Locally Toeplitz structure, i.e., we have proposed an analysis that is sharper and more informative than the one that can be obtained by looking only at the saddle–point structure. We have produced a localization of the spectrum in three intervals, up to a number of outliers infinitesimal in the dimension of the problem, and used this characterization to produce an asymptotically optimal preconditioner, i.e., a preconditioner that is independent of the value of the regularization parameter $\alpha$, and whose performance increases for finer grids.

We plan to extend this analysis in order that it can cover more general constraints, i.e., we would like to discuss also the case of sparse optimization, and bounded controls. Moreover, the GLT spectral analysis techniques we are using have been recently extended for becoming tools for the fast and reliable computation of generalized eigenvalues see, e.g., [13, 14], since we have analyzed the structure of the eigenvectors of our preconditioned problems (Proposition 3), we plan to investigate the possible application of deflation techniques to further accelerate our iterative methods.

Acknowledgment. We are thankful to Prof. S. Serra–Capizzano for the insightful discussions on the spectral distribution results, and to the referee whose suggestion have been extremely helpful in improving the presentation of the material.