Fast Solution Methods for Convex Quadratic Optimization of Fractional Differential Equations

TL;DR

This paper introduces an efficient ADMM-based numerical framework with specialized preconditioners for solving large-scale convex quadratic fractional differential equation constrained optimization problems, enabling solutions on complex space-time domains.

Contribution

The paper develops a novel ADMM approach combined with Krylov solvers and multilevel circulant preconditioners tailored for FDE-constrained optimization, preserving matrix structures and achieving linear storage and N log N complexity.

Findings

Method demonstrates scalability and efficiency in various FDE problem setups.

Preconditioners effectively handle large dense linear systems.

Approach extends to PDE optimization problems previously out of reach.

Abstract

In this paper, we present numerical methods suitable for solving convex quadratic Fractional Differential Equation (FDE) constrained optimization problems, with box constraints on the state and/or control variables. We develop an Alternating Direction Method of Multipliers (ADMM) framework, which uses preconditioned Krylov subspace solvers for the resulting sub-problems. The latter allows us to tackle a range of Partial Differential Equation (PDE) optimization problems with box constraints, posed on space-time domains, that were previously out of the reach of state-of-the-art preconditioners. In particular, by making use of the powerful Generalized Locally Toeplitz (GLT) sequences theory, we show that any existing GLT structure present in the problem matrices is preserved by ADMM, and we propose some preconditioning methodologies that could be used within the solver, to demonstrate the generality of the approach. Focusing on convex quadratic programs with time-dependent 2-dimensional FDE constraints, we derive multilevel circulant preconditioners, which may be embedded within Krylov subspace methods, for solving the ADMM sub-problems. Discretized versions of FDEs involve large dense linear systems. In order to overcome this difficulty, we design a recursive linear algebra, which is based on the Fast Fourier Transform (FFT). We manage to keep the storage requirements linear, with respect to the grid size $N$, while ensuring an order $N \log N$ computational complexity per iteration of the Krylov solver. We implement the proposed method, and demonstrate its scalability, generality, and efficiency, through a series of experiments over different setups of the FDE optimization problem.

Tables

Table 1. Table 1 : Summary of computational and storage complexity

Structure	Operation	Computations	Storage
$d$-level circulant	$C_{𝒏}x$	$O\left(N\left(𝒏\right)\log N\left(𝒏\right)\right)$	$O\left(N\left(𝒏\right)\right)$
$d$-level circulant	$C_{𝒏}^{-1}x$	$O\left(N\left(𝒏\right)\log N\left(𝒏\right)\right)$	$O\left(N\left(𝒏\right)\right)$
$d$-level circulant	$C_{𝒏}^{\left(1\right)}{C}_{𝒏}^{\left(2\right)}$	$O\left(N\left(𝒏\right)\right)$	$O\left(N\left(𝒏\right)\right)$
$d$-level circulant	$C_{𝒏}^{\left(1\right)}+C_{𝒏}^{\left(2\right)}$	$O\left(N\left(𝒏\right)\right)$	$O\left(N\left(𝒏\right)\right)$
$d$-level Toeplitz	$T_{𝒏}x$	$O\left(2^{d}N\left(𝒏\right)\log N\left(𝒏\right)\right)$	$O\left(2^{d}N\left(𝒏\right)\right)$
$d$-level circulant	Construct $C_d\left(T_{𝒏}\right)$	$O\left(N\left(𝒏\right)\log N\left(𝒏\right)\right)$	$O\left(N\left(𝒏\right)\right)$

Table 2. Table 2 : Inequalities on the state: Varying restriction bounds (with N = 50 3 𝑁 superscript 50 3 N=50^{3} , β x = β y = 1.3 subscript 𝛽 𝑥 subscript 𝛽 𝑦 1.3 \beta_{x}=\beta_{y}=1.3 , α = 0.7 𝛼 0.7 \alpha=0.7 , γ = 10 − 4 𝛾 superscript 10 4 \gamma=10^{-4} , δ = 0.1 𝛿 0.1 \delta=0.1 ).

${𝒚}_{𝒂}$	${𝓔}_{{𝑳}^{\mathrm{𝟐}}}\left(𝒚-\overline{𝒚}\right)$	Dual Inf.	Iterations		Time (s)
missingmissingmissing  4-5
	PCG	ADMM
$-7$	5.60 $\times {10}^{-1}$(*)111(*) means that the solution coincides with the equality constrained solution; all the variables lie strictly within the restriction bounds.	2.13 $\times {10}^{-3}$	9	75	142.31
$-5$	5.80 $\times {10}^{-1}$	3.14 $\times {10}^{-3}$	10	105	206.73
$-3$	7.88 $\times {10}^{-1}$	4.11 $\times {10}^{-3}$	10	100	193.77
$-1$	1.38 $\times {10}^0$	8.74 $\times {10}^{-4}$	10	86	173.49

Table 3. Table 3 : Inequalities on the control: Varying restriction bounds (with N = 50 3 𝑁 superscript 50 3 N=50^{3} , β x = β y = 1.3 subscript 𝛽 𝑥 subscript 𝛽 𝑦 1.3 \beta_{x}=\beta_{y}=1.3 , α = 0.7 𝛼 0.7 \alpha=0.7 , γ = 10 − 4 𝛾 superscript 10 4 \gamma=10^{-4} , δ = 0.4 𝛿 0.4 \delta=0.4 ).

${𝒖}_{𝒂}$	${𝓔}_{{𝑳}^{\mathrm{𝟐}}}\left(𝒚-\overline{𝒚}\right)$	Dual Inf.	Iterations		Time (s)
missingmissingmissing  4-5
	PCG	ADMM
$-400$	5.60 $\times {10}^{-1}$(*)	8.43 $\times {10}^{-4}$	16	30	84.46
$-300$	5.65 $\times {10}^{-1}$	8.79 $\times {10}^{-4}$	19	22	72.77
$-200$	6.25 $\times {10}^{-1}$	4.39 $\times {10}^{-4}$	17	28	85.55
$-100$	8.90 $\times {10}^{-1}$	1.49 $\times {10}^{-4}$	18	65	205.69

Table 4. Table 4 : Inequalities on both variables: Varying restriction bounds (with β x = β y = 1.3 subscript 𝛽 𝑥 subscript 𝛽 𝑦 1.3 \beta_{x}=\beta_{y}=1.3 , α = 0.7 𝛼 0.7 \alpha=0.7 , γ = 10 − 4 𝛾 superscript 10 4 \gamma=10^{-4} , δ = 0.4 𝛿 0.4 \delta=0.4 , N = 50 3 𝑁 superscript 50 3 N=50^{3} ).

${𝒚}_{𝒂}$	${𝒖}_{𝒂}$	${𝓔}_{{𝑳}^{\mathrm{𝟐}}}\left(𝒚-\overline{𝒚}\right)$	Dual Inf.	Iter.		Time (s)
missingmissingmissing  5-6
	PCG	ADMM
$-7$	$-400$	5.60 $\times {10}^{-1}$(*)	4.88 $\times {10}^{-3}$	10	36	74.20
$-7$	$-200$	5.94 $\times {10}^{-1}$	2.35 $\times {10}^{-3}$	11	38	80.14
$-4$	$-350$	6.45 $\times {10}^{-1}$	1.99 $\times {10}^{-4}$	18	126	412.66
$-1$	$-400$	1.38 $\times {10}^0$	2.56 $\times {10}^{-4}$	19	109	377.86

Table 5. Table 5 : Inequalities on both variables: Varying grid size (with β x = β y = 1.3 subscript 𝛽 𝑥 subscript 𝛽 𝑦 1.3 \beta_{x}=\beta_{y}=1.3 , α = 0.7 𝛼 0.7 \alpha=0.7 , γ = 10 − 4 𝛾 superscript 10 4 \gamma=10^{-4} , y a = − 4 subscript 𝑦 𝑎 4 y_{a}=-4 , u a = − 350 subscript 𝑢 𝑎 350 u_{a}=-350 ).

$𝑵$	$𝜹$	${𝓔}_{{𝑳}^{\mathrm{𝟐}}}\left(𝒚-\overline{𝒚}\right)$	Dual Inf.	Iter.		Time (s)
missingmissingmissing  5-6
	PCG	ADMM
$8^3$	$2$	3.87 $\times {10}^{-1}$	5.23 $\times {10}^{-4}$	12	86	1.89
${16}^3$	$2$	5.02 $\times {10}^{-1}$	8.68 $\times {10}^{-5}$	13	58	6.06
${32}^3$	$0.4$	6.09 $\times {10}^{-1}$	2.94 $\times {10}^{-4}$	16	62	75.04
${50}^3$	$0.4$	6.45 $\times {10}^{-1}$	1.99 $\times {10}^{-4}$	18	126	412.66
${64}^3$	$0.1$	6.58 $\times {10}^{-1}$	3.34 $\times {10}^{-3}$	17	97	987.12
${80}^3$	$0.1$	6.65 $\times {10}^{-1}$	4.31 $\times {10}^{-3}$	17	102	1,135.83
${100}^3$	$0.1$	6.70 $\times {10}^{-1}$	4.91 $\times {10}^{-3}$	17	119	2,436.17
${128}^3$	$0.1$	6.73 $\times {10}^{-1}$	3.49 $\times {10}^{-3}$	17	169	9,077.08

Table 6. Table 6 : Inequalities on both variables: Varying fractional derivative orders (with N = 32 4 𝑁 superscript 32 4 N=32^{4} , y a = − 4 subscript 𝑦 𝑎 4 y_{a}=-4 , u a = − 350 subscript 𝑢 𝑎 350 u_{a}=-350 , γ = 10 − 4 𝛾 superscript 10 4 \gamma=10^{-4} ).

$𝜶$	$𝜷$	$𝜹$	${𝓔}_{{𝑳}^{\mathrm{𝟐}}}\left(𝒚-\overline{𝒚}\right)$	Dual Inf.	Iter.		Time (s)
missingmissingmissing  6-7
	PCG	ADMM
$0.1$	$1.3$	0.4	6.46 $\times {10}^{-1}$	1.60 $\times {10}^{-4}$	17	126	380.75
$0.3$	$1.3$	0.4	6.46 $\times {10}^{-1}$	2.56 $\times {10}^{-4}$	17	126	385.96
$0.5$	$1.3$	0.4	5.12 $\times {10}^{-1}$	2.83 $\times {10}^{-4}$	18	126	408.47
$0.9$	$1.3$	0.4	6.44 $\times {10}^{-1}$	3.10 $\times {10}^{-4}$	19	125	419.46
$0.7$	$1.1$	0.4	6.48 $\times {10}^{-1}$	1.21 $\times {10}^{-3}$	30	100	508.59
$0.7$	$1.5$	0.1	7.79 $\times {10}^{-1}$	4.23 $\times {10}^{-4}$	15	96	275.03
$0.7$	$1.7$	0.4	1.04 $\times {10}^0$	2.46 $\times {10}^{-4}$	13	113	275.21
$0.7$	$1.9$	0.1	1.36 $\times {10}^0$	1.36 $\times {10}^{-3}$	8	108	180.85

Table 7. Table 7 : Inequalities on both variables: Varying regularization (with N = 50 3 𝑁 superscript 50 3 N=50^{3} , α = 0.7 𝛼 0.7 \alpha=0.7 , β = 1.3 𝛽 1.3 \beta=1.3 ).

$𝜸$	${𝒚}_{𝒂}$	${𝒖}_{𝒂}$	$𝜹$	${𝓔}_{{𝑳}^{\mathrm{𝟐}}}\left(𝒚-\overline{𝒚}\right)$	Dual Inf.	Iter.		Time (s)
missingmissingmissing  7-8
	PCG	ADMM
${10}^{-2}$	$-2$	$-100$	0.1	1.77 $\times {10}^0$(*)	1.38 $\times {10}^{-3}$	11	87	190.99
${10}^{-4}$	$-7$	$-400$	0.4	5.60 $\times {10}^{-1}$(*)	4.88 $\times {10}^{-3}$	10	36	74.20
${10}^{-6}$	$-9$	$-2,800$	10	1.28 $\times {10}^{-1}$(*)	6.03 $\times {10}^{-4}$	8	47	71.13
${10}^{-8}$	$-9$	$-4,000$	100	1.13 $\times {10}^{-1}$(*)	2.18 $\times {10}^{-4}$	6	32	44.70
${10}^{-10}$	$-9$	$-4,000$	100	1.13 $\times {10}^{-1}$(*)	2.18 $\times {10}^{-4}$	5	32	40.77

Equations

$$\bigg{[}\ldots\big{[}[(j_{1},\ldots,j_{d})]_{j_{d}=h_{d},\ldots,k_{d}}\big{]}_{j_{d-1}=h_{d-1},\ldots,k_{d-1}}\ldots\bigg{]}_{j_{1}=h_{1},\ldots,k_{1}}.$$

$$\underset{M\to \infty }{\lim }\underset{n\to \infty }{\mathrm{limsup}}\frac{\#\{i\in \{1,\dots ,N(\mathbf{n})\}:{\sigma }_i(A_{\mathbf{n}})>M\}}{N(\mathbf{n})}=0,$$

$$\underset{M\to \infty }{\lim }\underset{n\to \infty }{\mathrm{limsup}}\frac{\#\{i\in \{1,\dots ,N(\mathbf{n})\}:{\sigma }_i(A_{\mathbf{n}})<1/M\}}{N(\mathbf{n})}=0,$$

$$\#\{j\in \{1,\dots ,d_n\}:{\lambda }_j(A_n)\notin D(S,ϵ)\}=O(1),$$

$$\#\{j\in \{1,\dots ,d_n\}:{\lambda }_j(A_n)\notin D(S,ϵ)\}=o(d_n).$$

$${\varphi }_g:C_c(\mathbb{K})\mapsto \mathbb{C},{\varphi }_g(F)=\frac{1}{{\mu }_d(D)}{\int }_{D}F(g(\mathbf{x}))\mathrm{d}\mathbf{x}.$$

$$\underset{n\to \infty }{\lim }\frac{1}{d_n}\sum _{j=1}^{d_n}F({\lambda }_j(A_n))=\varphi (F),\forall F\in C_c(\mathbb{C}).$$

$$D_{\bm{n}}(v)=\underset{\bm{i}=\bm{1},\ldots,\bm{n}}{\textnormal{diag}}v\bigg{(}\frac{\bm{i}}{\bm{n}}\bigg{)}.$$

$$f_{\mathbf{k}}=\frac{1}{(2\pi {)}^d}{\int }_{[-\pi ,\pi {]}^d}f(\mathbf{\theta })e^{-i⟨\mathbf{k},\mathbf{\theta }⟩}\mathrm{d}\mathbf{\theta },\mathbf{k}\in {\mathbb{Z}}^d,$$

$$T_{\mathbf{n}}(f)=[f_{\mathbf{i}-\mathbf{j}}{]}_{\mathbf{i},\mathbf{j}=1}^{\mathbf{n}}=\sum _{\mathbf{k}=-(\mathbf{n}-1)}^{\mathbf{n}-1}{f}_{\mathbf{k}}{J}_{\mathbf{n}}^{(\mathbf{k})}.$$

$$\big{[}a_{(\bm{i}-\bm{j})\bmod{\bm{n}}}\big{]}_{\bm{i},\bm{j}=\bm{1}}^{\bm{n}}=\sum_{\bm{k}=\bm{0}}^{\bm{n}-\bm{1}}a_{\bm{k}}C_{\bm{n}}^{\bm{k}},$$

$$\sum _{\mathbf{k}=-\mathbf{r}}^{\mathbf{r}}{c}_{\mathbf{k}}{C}_{\mathbf{n}}^{\mathbf{k}}=F_{\mathbf{n}}^{\ast }\left(\underset{\mathbf{j}=0,\dots ,\mathbf{n}-1}{\mathrm{diag}}c\left(\frac{2\pi \mathbf{j}}{\mathbf{n}}\right)\right)F_{\mathbf{n}},$$

$$\lambda\bigg{(}\sum_{\bm{k}=-\bm{r}}^{\bm{r}}c_{\bm{k}}C_{\bm{n}}^{\bm{k}}\bigg{)}=\left\{c\left(\frac{2\pi\bm{j}}{\bm{n}}\right):\bm{j}=\bm{0},\ldots,\bm{n}-\bm{1}\right\}.$$

$$A_n=B_{n,m}+R_{n,m}+N_{n,m},\text{rank}(R_{n,m})\le c(m)d_n,\Vert N_{n,m}\Vert \le \omega (m),$$

$$\underset{m\to \infty }{\lim }c(m)=\underset{m\to \infty }{\lim }{\omega }_m=0.$$

$$\begin{split}LT_{n}^{m}(v,f)=&\ \bigg{(}D_{m}(v)\otimes T_{\lfloor n/m\rfloor}(f)\bigg{)}\oplus O_{n\bmod{m}},\end{split}$$

$$\begin{split}LT_{\bm{n}}^{\bm{m}}(v,f_{1}\otimes\ldots\otimes f_{d})=LT_{n_{1},\ldots,n_{d}}^{m_{1},\ldots,m_{d}}\big{(}v(x_{1},\ldots,x_{d}),f_{1}\otimes\ldots\otimes f_{d}\big{)}.\end{split}$$

$$\{L{T}_{\mathbf{n}}^{\mathbf{m}}(v,f){\}}_n\stackrel{\text{a.c.s.}}{\to }\{A_{\mathbf{n}}{\}}_n,\text{as }\mathbf{m}\to \infty .$$

$$\sum_{i=1}^{N_{\epsilon}}v_{i,\epsilon}\otimes f_{i,\epsilon}\rightarrow\kappa\textnormal{ in measure, and }\big{\{}\sum_{i=1}^{N_{\epsilon}}A_{\bm{n}}^{(i,\epsilon)}\big{\}}_{n}\xrightarrow{\textnormal{a.c.s.}}\{A_{\bm{n}}\}_{n}.$$

$$\mathrm{Dy}(\mathbf{x},t)=\mathrm{g}(\mathbf{x},t),$$

$$\{D_{\mathbf{n}}{y}_{\mathbf{n}}{\}}_n=\{g_{\mathbf{n}}{\}}_n,$$

$$\mathrm{Dy}(\mathbf{x},t)+\mathrm{u}(\mathbf{x},t)=\mathrm{g}(\mathbf{x},t).$$

$$\begin{array}{rl}\underset{\mathrm{y},\mathrm{u}}{\min }& \mathrm{J}(\mathrm{y}(\mathbf{x},t),\mathrm{u}(\mathbf{x},t))={\mathrm{J}}_1(\mathrm{y}(\mathbf{x},t))+{\mathrm{J}}_2(\mathrm{u}(\mathbf{x},t))\\ \text{s.t.}& \mathrm{Dy}(\mathbf{x},t)+\mathrm{u}(\mathbf{x},t)=\mathrm{g}(\mathbf{x},t),\\ & {\mathrm{y}}_{\mathrm{a}}(\mathbf{x},t)\le \mathrm{y}(\mathbf{x},t)\le {\mathrm{y}}_{\mathrm{b}}(\mathbf{x},t),{\mathrm{u}}_{\mathrm{a}}(\mathbf{x},t)\le \mathrm{u}(\mathbf{x},t)\le {\mathrm{u}}_{\mathrm{b}}(\mathbf{x},t).\end{array}$$

$$\begin{split}\min_{y_{\bm{n}},u_{\bm{n}}}\ &\ \bigg{(}\frac{1}{2}(y_{\bm{n}}-\bar{y}_{\bm{n}})^{*}J_{1_{\bm{n}}}(y_{\bm{n}}-\bar{y}_{\bm{n}})+\frac{1}{2}u_{\bm{n}}^{*}J_{2_{\bm{n}}}u_{\bm{n}}\bigg{)}\\ \text{s.t.}\ &\ D_{\bm{n}}y_{\bm{n}}+u_{\bm{n}}=g_{\bm{n}},\\ \ &\ y_{a_{\bm{n}}}\leq y_{\bm{n}}\leq y_{b_{\bm{n}}},\quad u_{a_{\bm{n}}}\leq u_{\bm{n}}\leq u_{b_{\bm{n}}},\end{split}$$

$$\begin{split}\min_{y_{\bm{n}},u_{\bm{n}},z_{y_{\bm{n}}},z_{u_{\bm{n}}}}&\ \bigg{(}\frac{1}{2}(y_{\bm{n}}-\bar{y}_{\bm{n}})^{*}J_{1_{\bm{n}}}(y_{\bm{n}}-\bar{y}_{\bm{n}})+\frac{1}{2}u_{\bm{n}}^{*}J_{2_{\bm{n}}}u_{\bm{n}}\bigg{)}\\ \text{s.t.}\ \ &\ D_{\bm{n}}y_{\bm{n}}+u_{\bm{n}}=g_{\bm{n}},\\ &\ y_{\bm{n}}=z_{y_{\bm{n}}},\ u_{\bm{n}}=z_{u_{\bm{n}}},\\ &\ y_{a_{\bm{n}}}\leq z_{y_{\bm{n}}}\leq y_{b_{\bm{n}}},\ u_{a_{\bm{n}}}\leq z_{u_{\bm{n}}}\leq u_{b_{\bm{n}}}.\end{split}$$

$$\begin{array}{rl}{\mathcal{L}}_{\delta }(y_{\mathbf{n}},u_{\mathbf{n}},z_{y_{\mathbf{n}}},& z_{u_{\mathbf{n}}},p_{{}_{\mathbf{n}}},w_{y_{\mathbf{n}}},w_{u_{\mathbf{n}}})=\frac{1}{2}(y_{\mathbf{n}}-{\bar{y}}_{\mathbf{n}}{)}^{\ast }{J}_{1_{\mathbf{n}}}(y_{\mathbf{n}}-{\bar{y}}_{\mathbf{n}})+\frac{1}{2}{u}_{\mathbf{n}}^{\ast }{J}_{2_{\mathbf{n}}}{u}_{\mathbf{n}}\\ & +p_{\mathbf{n}}^{\ast }(D_{\mathbf{n}}{y}_{\mathbf{n}}+u_{\mathbf{n}}-g_{\mathbf{n}})+w_{y_{\mathbf{n}}}^{\ast }(y_{\mathbf{n}}-z_{y_{\mathbf{n}}})+w_{u_{\mathbf{n}}}^{\ast }(u_{\mathbf{n}}-z_{u_{\mathbf{n}}})\\ & +\frac{1}{2\delta }\left(\Vert D_{\mathbf{n}}{y}_{\mathbf{n}}+u_{\mathbf{n}}-g_{\mathbf{n}}{\Vert }_2^2+\Vert y_{\mathbf{n}}-z_{y_{\mathbf{n}}}{\Vert }_2^2+\Vert u_{\mathbf{n}}-z_{u_{\mathbf{n}}}{\Vert }_2^2\right),\end{array}$$

$$(y_{\mathbf{n}}^{j+1},u_{\mathbf{n}}^{j+1})$$

$$(z_{y_{\mathbf{n}}}^{j+1},z_{u_{\mathbf{n}}}^{j+1})$$

$$p_{\mathbf{n}}^{j+1}$$

$$(w_{y_{\mathbf{n}}}^{j+1},w_{u_{\mathbf{n}}}^{j+1})$$

$$\begin{array}{r}\left[\begin{array}{cc}{J}_{1_{\mathbf{n}}}+\frac{1}{\delta }(D_{\mathbf{n}}^{\ast }{D}_{\mathbf{n}}+I_{\mathbf{n}})& \frac{1}{\delta }{D}_{\mathbf{n}}^{\ast }\\ \frac{1}{\delta }{D}_{\mathbf{n}}& J_{2_{\mathbf{n}}}+\frac{2}{\delta }{I}_{\mathbf{n}}\end{array}\right]\left[\begin{array}{c}{y}_{\mathbf{n}}\\ u_{\mathbf{n}}\end{array}\right]=\left[\begin{array}{c}{\eta }_1\\ {\eta }_2\end{array}\right],\end{array}$$

$${\eta }_1=J_{1_{\mathbf{n}}}{\bar{y}}_{\mathbf{n}}-D_{\mathbf{n}}^{\ast }{p}_{\mathbf{n}}^j-w_{y_{\mathbf{n}}}^j+\frac{1}{\delta }(D_{\mathbf{n}}^{\ast }{g}_{\mathbf{n}}+z_{y_{\mathbf{n}}}^j),{\eta }_2=-p_{\mathbf{n}}^j-w_{u_{\mathbf{n}}}^j+\frac{1}{\delta }(g_{\mathbf{n}}+z_{u_{\mathbf{n}}}^j).$$

Full text

1 Introduction

Optimization problems with Differential Equations (Partial (PDEs) or Ordinary (ODEs)) as constraints have received a great deal of attention within the applied mathematics and engineering communities, due in particular to their wide applicability across many fields of science. In addition to classical differential equation constraints, one may also use Fractional Differential Equations (FDEs) in order to model processes that could not otherwise be modeled using integer derivatives. In fact, there is a wide and increasing use of FDEs in the literature. Among other processes, FDEs have been used to model viscoelasticity (e.g. [43]), chaotic systems (e.g. [76]), turbulent flow, or anomalous diffusion (e.g. [7]). In particular, since the fractional operator is non-local, problems with non-local properties can frequently be modeled accurately using FDEs (see [62] for an extended review).

Availability of closed form solutions for FDEs is rare, and hence various numerical schemes for solving them have been developed and analyzed in the literature (see [50, 51, 52] for finite difference, and [19, 31] for finite element methods). Such numerical schemes produce dense matrices, making the solution or even the storage of FDE-constrained optimization problems extremely difficult for fine grids. Naturally, this behavior is even more severe in the case of multidimensional FDEs. In light of the previous, employing standard (black-box) direct approaches for solving such problems requires $O(N^{3})$ operations and $O(N^{2})$ storage, where $N$ is the number of grid points. Iterative methods with general purpose preconditioners also suffer from similar issues.

Various specialized solution methods have been proposed in the literature, aiming at lowering the computational and storage cost of solving such problems (see for instance [17, 23, 24, 40, 53, 77]). One popular and effective approach is to employ tensor product solvers. Such specialized methods have been proposed for solving high-dimensional FDE-constrained inverse problems with great success, even for very fine discretizations (see for example [23, 40] and the references therein). While these solvers are highly scalable (with respect to the grid size), to date they have been tailored solely to problems with specific cost functionals and without additional algebraic constraints. Another popular approach is based on the observation that multidimensional FDEs possess a multilevel Toeplitz-like structure. It is well known that such matrices can be very well approximated by banded multilevel Toeplitz (see for example [25, 53]) or multilevel circulant matrices (see [14, 15, 32, 44, 45]). The former are usually sparse and can be inverted using specially designed multigrid or factorization methods, while the latter can be inverted or applied to a vector in only $O(N\log N)$ operations using the Fast Fourier Transform (FFT) (e.g. [75]). The idea is to apply a Krylov subspace solver, supported by a banded Toeplitz or circulant preconditioner, in order to solve the optimality conditions of the problem. One is able to redesign the underlying linear algebra, in order to achieve an $O(N\log N)$ iteration complexity for the Krylov solver, with $O(N)$ overall storage requirements (see for example [44, 45, 46]). While such solution methods are certainly more general (although usually slower), when compared to tensor product solvers, they remain rather sensitive in terms of the underlying structure. In particular, to the authors’ knowledge, no such method has been proposed for the solution of more general FDE optimization problems, for instance those which include box constraints on the state and control variables. We highlight that a time-independent problem, with box constraints on the control, is studied in [29], and the authors attempt to solve it using a Limited-memory Broyden–Fletcher–Goldfarb–Shanno (L–BFGS) method.

In this paper, we present an optimization method suitable for solving convex quadratic PDE-constrained optimization problems with box constraints on the state and control variables. In particular, we assume that we are given an arbitrary PDE constrained inverse problem, an associated discretization method, and that the resulting sequences of discrete matrices belong to the class of Generalized Locally Toeplitz (GLT) sequences (we refer the reader to [34, 35, 70] for a comprehensive overview of the powerful GLT theory). Then, we propose the use of an Alternating Direction Method of Multipliers (ADMM) for solving the discretized optimization problems. We employ ADMM in order to separate the equality from the inequality constraints. As a consequence of this choice, we show that the linear systems required to be solved during the iterations of ADMM preserve the GLT structure of the initial problem matrices. Using this structure, we present and analyze some general methodologies for efficiently preconditioning such linear systems, and solving them using an appropriate Krylov subspace method. The Krylov subspace method is in turn, under certain mild assumptions, expected to converge in a number of iterations independent of the grid size. Subsequently, we focus on a certain class of convex quadratic optimization problems with FDE constraints. In particular, we consider time-dependent 2-dimensional FDEs, and we precondition the associated discretized matrices using multilevel circulant preconditioners. We manage to keep the storage requirements linear, with respect to the grid size $N$, while ensuring an order $N\log N$ computational complexity per iteration of the Krylov solver inside ADMM. We implement the proposed method, and demonstrate its robustness and efficiency, through a series of experiments over different setups of the FDE optimization problem.

This paper is structured as follows. In Section 2, we provide the relevant theoretical background as well as the notation used throughout the paper. Subsequently, in Section 3 we present the proposed ADMM framework, as well as possible preconditioning strategies that could be used to accelerate the solution of the ADMM sub-problems, given the assumption that the associated matrices possess a GLT structure. In Section 4 we present the FDE-constrained optimization problem under consideration. Then, in Section 5, we propose the use of a multilevel circulant preconditioner for approximating multilevel Toeplitz matrices arising from the discretization of the FDE under consideration, while demonstrating that such a preconditioner is effective for the problem at hand. In Section 6, we discuss the implementation details of the proposed approach and present some numerical results. Finally, in Section 7, we state our conclusions.

2 Notation and Theoretical Background

In this section, we introduce some notation and provide the theoretical background that will be used in the rest of this manuscript. Firstly, we introduce the notion of $d$-indices which will allow us to compactly represent multilevel matrices. For brevity of presentation, we only discuss the crucial notions that will be used in this paper. A more complete presentation of the notation and theory of this section can be found in [34, 35]. A reader familiar with the theory of GLT sequences can skip directly to Section 3.

Definition 2.1**.**

A multi-index $\bm{i}$ of size $d$ ($d$-index) is a row vector in $\mathbb{Z}^{d}$ with components $i_{1},\ldots,i_{d}$. Using this notation, we define the following notions:

• •

$\bm{0,1,2,}\ldots$, are the row vectors of all zeros, ones, twos, etc.

• •

$N(\bm{i})=\prod_{j=1}^{d}i_{j}$* and we write $\bm{i}\rightarrow\infty$ to indicate that $\min(\bm{i})\rightarrow\infty$.*

• •

Given two d-indices $\bm{h,k}$, we write $\bm{h}\leq\bm{k}$ to express that $h_{j}\leq k_{j}$, for all $j\in\{1,\ldots,d\}$. The d-index range $\bm{h},\ldots,\bm{k}$ is a set of cardinality $N(\bm{k-h+1})$ given by $\{\bm{j}\in\mathbb{Z}^{d}:\bm{h}\leq\bm{j}\leq\bm{k}\}$. The latter set is assumed to be ordered under the lexicographical ordering, that is:

[TABLE]

• •

Let a d-index $\bm{m}\in\mathbb{N}^{d}$, and define $\bm{x}=[x_{\bm{i}}]_{\bm{i}=\bm{1}}^{\bm{m}}$ ($\bm{X}=[x_{\bm{i},\bm{j}}]_{\bm{i},\bm{j}=1}^{\bm{m}}$, respectively). Then $\bm{x}$ ($\bm{X}$, respectively) is a vector of size $N(\bm{m})$ (a matrix of size $N(\bm{m})\times N(\bm{m})$, respectively).

• •

Any operation involving d-indices that has no meaning in the vector space $\mathbb{Z}^{d}$ will be interpreted in a componentwise sense.

A matrix $A$ of size $N$ is a $d$-level matrix with level orders $n_{1},\ldots,n_{d}$ if $N=n_{1}n_{2}\cdots n_{d}$ and it is partitioned into $n_{1}^{2}$ square blocks of size $\frac{N}{n_{1}}$, each of which is partitioned into $n_{2}^{2}$ blocks of size $\frac{N}{n_{1}n_{2}}$, and so on until the last $n_{d}^{2}$ blocks of size $1$. Then, $A$ can be written as $A=[A_{\bm{ij}}]_{\bm{i},\bm{j}=\bm{1}}^{\bm{n}}$, where $A_{\bm{ij}}=A_{i_{1}j_{1};\ldots;i_{d}j_{d}}$, for $\bm{i},\bm{j}=\bm{1},\ldots,\bm{n}$.

Next, we define the notion of a matrix-sequence, which is a fundamental element for studying the asymptotic spectral behavior of structured matrices arising from some discretization of a physical process. In the rest of this manuscript, given an arbitrary matrix $A$, $\sigma(A)$ denotes the set of singular values of the matrix, while $\lambda(A)$ denotes the set of eigenvalues of the matrix $A$ (given that they exist).

Definition 2.2**.**

A matrix-sequence is a sequence of the form $\{A_{n}\}_{n}$, where $n$ varies over some infinite subset of $\mathbb{N}$, $A_{n}$ is a square matrix of size $d_{n}$, and $d_{n}\rightarrow\infty$ as $n\rightarrow\infty$. In particular, a $d$-level matrix sequence is a sequence of the form $\{A_{\bm{n}}\}_{n}$, where $A_{\bm{n}}$ is a matrix of size $N(\bm{n})\times N(\bm{n})$, $n$ varies over some infinite subset of $\mathbb{N}$, and $\bm{n}=\bm{n}(n)\in\mathbb{N}^{d}$ is such that $\bm{n}\rightarrow\infty$, as $n\rightarrow\infty$. Given a d-level matrix-sequence $\{A_{\bm{n}}\}_{n}$, we say that it is sparsely unbounded (and denote that as s.u.) if:

[TABLE]

where $\#S$ denotes the cardinality of a set $S$. Similarly, we say that $\{A_{\bm{n}}\}_{n}$ is sparsely vanishing (and denote that as s.v.) if:

[TABLE]

where we assume that $1/\infty=0$.

An important notion is that of clustering. In order to define it we let, for every $z\in\mathbb{C}$ and any $\epsilon>0$, $D(z,\epsilon)$ represent the disk with center $z$ and radius $\epsilon$. If $S\subseteq\mathbb{C}$ and $\epsilon>0$, $D(S,\epsilon)$ denotes the $\epsilon$-expansion of $S$, defined as $D(S,\epsilon)=\bigcup_{z\in S}D(z,\epsilon)$.

Definition 2.3**.**

Let $\{A_{n}\}_{n}$ be a sequence of matrices, with $A_{n}$ of size $d_{n}\times d_{n}$, and let $S\subseteq\mathbb{C}$ be a non-empty subset of $\mathbb{C}$. We say that $\{A_{n}\}_{n}$ is strongly clustered at $S$ (in the sense of eigenvalues) if $\forall\ \epsilon>0$ we have:

[TABLE]

and weakly clustered at $S$ if $\forall\ \epsilon>0$,

[TABLE]

Clustering in the sense of singular values is defined analogously.

Let $f_{m},\ f:D\subseteq\mathbb{R}^{d}\mapsto\mathbb{C}$ be measurable functions, with respect to the Lebesgue measure $\mu_{d}$ in $\mathbb{R}^{d}$. We say that $f_{m}\rightarrow f$ in measure if, for every $\epsilon>0$, \underset{m\rightarrow\infty}{\lim}\mu_{d}\big{(}\{|f_{m}-f|>\epsilon\}\big{)}=0. Furthermore, $f_{m}\rightarrow f$ a.e. (almost everywhere) if \mu_{d}\big{(}\{f_{m}\nrightarrow f\}\big{)}=0.

Lemma 2.4**.**

Let $f_{m},\ g_{m},\ f,\ g:D\subseteq\mathbb{R}^{d}\mapsto\mathbb{C}$ be measurable functions.

If $f_{m}\rightarrow f$ in measure, then $|f_{m}|\rightarrow|f|$ in measure. 2. 2.

If $f_{m}\rightarrow f$ in measure and $g_{m}\rightarrow g$ in measure, then $\alpha f_{m}+\beta g_{m}\rightarrow\alpha f+\beta g$ in measure for all $\alpha,\ \beta\in\mathbb{C}$. 3. 3.

If $f_{m}\rightarrow f$ in measure, $g_{m}\rightarrow g$ in measure, and $\mu_{d}(D)<\infty$, then $f_{m}g_{m}\rightarrow fg$ in measure.

Proof.

This is stated in [34, Lemma 2.3] and proved in [4, Corollary 2.2.6]. ∎

Let $C_{c}(\mathbb{C})$ ($C_{c}(\mathbb{R})$, respectively) be the space of complex-(real-)valued continuous functions defined on $\mathbb{C}$ (or $\mathbb{R}$) with compact support. Given a field $\mathbb{K}$ ($=\mathbb{C}$ or $\mathbb{R}$) and a measurable function $g:D\subset\mathbb{R}^{d}\mapsto\mathbb{K}$, with $0<\mu_{d}(D)<\infty$, define the functional:

[TABLE]

Definition 2.5**.**

Let $\{A_{n}\}_{n}$ be a matrix-sequence, with $A_{n}$ of size $d_{n}\times d_{n}$. We say that $\{A_{n}\}_{n}$ has an asymptotic eigenvalue (spectral) distribution described by a functional $\phi:C_{c}(\mathbb{C})\mapsto\mathbb{C}$, and we write $\{A_{n}\}_{n}\sim_{\lambda}\phi$, if:

[TABLE]

If $\phi=\phi_{f}$ for some measurable function $f:D\subset\mathbb{R}^{d}\mapsto\mathbb{C}$, where $0<\mu_{d}(D)<\infty$, we say that $\{A_{n}\}_{n}$ has an asymptotic spectral distribution described by $f$ and we write $\{A_{n}\}_{n}\sim_{\lambda}f$. Then, $f$ is referred to as the eigenvalue (spectral) symbol of $\{A_{n}\}_{n}$.

We can define the asymptotic singular value distribution of a matrix sequence similarly to Definition 2.5 (see [33, Definition 2.1]). In that case, we write $\{A_{n}\}_{n}\sim_{\sigma}f$.

Below we define two important classes of matrix sequences, namely diagonal sampling and Toeplitz matrix-sequences.

Definition 2.6**.**

Let a function $v:[0,1]^{d}\mapsto\mathbb{C}$ be given. The $\bm{n}$-th diagonal sampling matrix generated by $v$ is denoted by $D_{\bm{n}}(v)$ and is defined by the following $N(\bm{n})\times N(\bm{n})$ diagonal matrix:

[TABLE]

Definition 2.7**.**

Given a $d$-index $\bm{n}\in\mathbb{N}^{d}$, a matrix of the form $[a_{\bm{i}-\bm{j}}]_{\bm{i},\bm{j}=\bm{1}}^{\bm{n}}\in\mathbb{C}^{N(\bm{n})\times N(\bm{n})}$ is called a d-level Toeplitz matrix. Unilevel Toeplitz matrices ($d=1$) are also defined as matrices that are constant along all of their diagonals.

A characterization of Toeplitz matrix-sequences is given by the following Theorem, the proof of which can be found in [35, Sections 3.1, 3.5]. Before that, let us define a useful matrix. Given an arbitrary $n\in\mathbb{N}$ and $k\in\mathbb{Z}$, define the $n\times n$ matrix $J_{n}^{(k)}$ such that $[J_{n}^{(k)}]_{ij}=1$ if $i-j=k$ and $[J_{n}^{(k)}]_{ij}=0$ otherwise. Given two $d$-indices $\bm{n}\in\mathbb{N}^{d}$ and $\bm{k}\in\mathbb{Z}^{d}$, we define $J_{\bm{n}}^{(\bm{k})}=J_{n_{1}}^{(k_{1})}\otimes J_{n_{2}}^{(k_{2})}\otimes\cdots\otimes J_{n_{d}}^{(k_{d})},$ where $\otimes$ denotes the Kronecker product between two matrices.

Theorem 2.8**.**

Let a function $f:[-\pi,\pi]^{d}\mapsto\mathbb{C}$ belonging to $L^{1}([-\pi,\pi]^{d})$ be given, with Fourier coefficients denoted by:

[TABLE]

where $\langle\bm{k},\bm{\theta}\rangle=\sum_{i=1}^{d}k_{i}\theta_{i}$. The $\bm{n}$-th ($d$-level) Toeplitz matrix associated with $f$ is defined as:

[TABLE]

Every $d$-level matrix sequence of the form $\{T_{\bm{n}}(f)\}_{n}$, with $\{\bm{n}=\bm{n}(n)\}_{n}\subseteq\mathbb{N}^{d}$ such that $\bm{n}\rightarrow\infty$ as $n\rightarrow\infty$, is called a (d-level) Toeplitz sequence generated by $f$, which in turn is referred to as the generating function of $\{T_{\bm{n}}(f)\}_{n}$. Furthermore, $\{T_{\bm{n}}(f)\}_{n}\sim_{\sigma}f.$ If moreover $f$ is real, then $\{T_{\bm{n}}(f)\}_{n}\sim_{\lambda}f.$

A special type of Toeplitz matrices are the circulant matrices, as defined below.

Definition 2.9**.**

A matrix of the form \big{[}a_{(\bm{i}-\bm{j})\bmod{\bm{n}}}\big{]}_{\bm{i},\bm{j}=\bm{1}}^{\bm{n}}\in\mathbb{C}^{N(\bm{n})\times N(\bm{n})}, for some $d$-index $\bm{n}\in\mathbb{N}^{d}$, is called a multilevel ($d$-level) circulant matrix.

Given an arbitrary $n\in\mathbb{N}$, define the $n\times n$ matrix $C_{n}$ such that $[C_{n}]_{ij}=1$ if $(i-j)\bmod{n}=1$ and $[C_{n}]_{ij}=0$ otherwise. Then, for $\bm{n}\in\mathbb{N}^{d}$ and $\bm{k}\in\mathbb{Z}^{d}$, let $C_{\bm{n}}^{\bm{k}}=C_{n_{1}}^{k_{1}}\otimes C_{n_{2}}^{k_{2}}\otimes\ldots\otimes C_{n_{d}}^{k_{d}}$, where $C_{n_{i}}^{k_{i}}$ is the previously defined matrix $C_{n}$ raised to the power $k_{i}$. Let $F_{n}$ denote the unitary discrete Fourier transform of order $n$. For any $\bm{n}\in\mathbb{N}^{d}$ let $F_{\bm{n}}=F_{n_{1}}\otimes\ldots\otimes F_{n_{d}}$. Below we provide a Theorem characterizing multilevel circulant matrices; its proof can be found in [35, Section 3.4].

Theorem 2.10**.**

The $d$-level circulant matrix admits the following expression:

[TABLE]

where $C_{\bm{n}}^{\bm{k}}$ is as defined earlier. Furthermore, letting any $\bm{r}\in\mathbb{N}^{d}$ and $c_{-\bm{r}},\ldots,c_{\bm{r}}\in\mathbb{C}$, we have that any linear combination of the form $\sum_{\bm{d}=-\bm{r}}^{\bm{r}}c_{\bm{k}}C_{\bm{n}}^{\bm{k}}$ is a $d$-level circulant matrix. Then,

[TABLE]

where $c(\bm{\theta})=\sum_{\bm{k}=-\bm{r}}^{\bm{r}}c_{\bm{k}}e^{i\langle\bm{k},\bm{\theta}\rangle}$, and $F_{\bm{n}}$ is the multilevel discrete Fourier transform. Moreover, $\sum_{\bm{k}=-\bm{r}}^{\bm{r}}c_{\bm{k}}C_{\bm{n}}^{\bm{k}}$ is a normal matrix the spectrum of which is given by:

[TABLE]

Let $\mathcal{C}_{\bm{n}}$ be the set of all $d$-level circulant matrices of size $N(\bm{n})\times N(\bm{n})$. In light of Theorem 2.10 we can see that the set $\mathcal{C}_{\bm{n}}$, together with matrix addition and multiplication, is a commutative ring. For more about circulant matrices see [18].

A very important notion of the theory of GLT sequences is that of the approximating class of sequences, which will be denoted as a.c.s.. In particular, it is very common in practice to approximate a “difficult” matrix-sequence by an “easier” sequence of matrix-sequences, which has the same asymptotic singular value or eigenvalue distribution. For example, such an “easier” sequence can be used to construct effective preconditioners inside a suitable Krylov subspace method. For the rest of this manuscript, given a matrix $X$, we denote its spectral norm by $\|X\|$.

Definition 2.11**.**

Let $\{A_{n}\}_{n}$ be a matrix-sequence, with $A_{n}$ of size $d_{n}\times d_{n}$, and let $\{\{B_{n,m}\}_{n}\}_{m}$ be a sequence of matrix-sequences, with $B_{n,m}$ of size $d_{n}\times d_{n}$. We say that $\{\{B_{n,m}\}_{n}\}_{m}$ is an approximating class of sequences (a.c.s.) for $\{A_{n}\}_{n}$ if for every $m$, there exists $n_{m}$ such that, for all $n\geq n_{m}$, we can write:

[TABLE]

where $n_{m},\ c(m)$, and $\omega(m)$ depend only on $m$, and are such that:

[TABLE]

In that case, we write $\{B_{n,m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{A_{n}\}_{n}$.

Below, we provide a result, as reported in [35, Theorem 2.9], which will be very useful when constructing suitable preconditioners later in this paper.

Theorem 2.12**.**

Let two matrix-sequences $\{A_{n}\}_{n},\ \{A^{\prime}_{n}\}_{n}$ be given, with $A_{n},\ A_{n}^{\prime}$ of size $d_{n}\times d_{n}$, and suppose that $\{B_{n,m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{A_{n}\}_{n}$ and $\{B^{\prime}_{n,m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{A^{\prime}_{n}\}_{n}$. The following properties hold:

$\{B_{n,m}^{*}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{A_{n}^{*}\}_{n}$. 2. 2.

$\{c_{1}B_{n,m}+c_{2}B^{\prime}_{n,m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{c_{1}A_{n}+c_{2}A^{\prime}_{n}\}_{n}$, for all $c_{1},\ c_{2}\in\mathbb{C}$. 3. 3.

If $\{A_{n}\}_{n}$ and $\{A^{\prime}_{n}\}_{n}$ are s.u., then $\{B_{n,m}B^{\prime}_{n,m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{A_{n}A^{\prime}_{n}\}_{n}$. 4. 4.

Suppose $\{A_{n}\}_{n}$ is s.v.. If $\{B_{n,m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{A_{n}\}_{n}$ then $\{B^{\dagger}_{n,m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{A^{\dagger}_{n}\}_{n}$.

All the previous definitions are used to define the notion of Locally Toeplitz (LT) sequences, which in turn are generalized to define the notion of GLT sequences. We briefly define this class of matrix-sequences here, and refer the reader to [34, 35] for a complete derivation of this class, and a vast amount of results concerning sequences belonging in the GLT class. This theory was originally developed in [70].

Definition 2.13**.**

Let $m,\ n\in\mathbb{N}$, $v:[0,1]\mapsto\mathbb{C}$, and $f\in L^{1}([-\pi,\pi])$. The 1-level locally Toeplitz operator is defined as the following $n\times n$ matrix:

[TABLE]

where $D_{m}(v)$ is a diagonal sampling matrix generated by $v$, $T_{\lfloor n/m\rfloor}(f)$ a Toeplitz matrix generated by $f$, and $O_{n\bmod{m}}$ a zero matrix. Let also $\bm{m},\ \bm{n}\in\mathbb{N}^{d}$, $v:[0,1]^{d}\mapsto\mathbb{C}$, and $f\in L^{1}([-\pi,\pi]^{d})$. The $d$-level locally Toeplitz operator is recursively defined as the following $N(\bm{n})\times N(\bm{n})$ matrix:

[TABLE]

Definition 2.13 allows us to recall the notion of a multilevel locally Toeplitz sequence.

Definition 2.14**.**

Let $\{A_{\bm{n}}\}_{n}$ be a $d$-level matrix-sequence, let $v:[0,1]^{d}\mapsto\mathbb{C}$ be Riemann-integrable and let $f\in L^{1}([-\pi,\pi]^{d})$. We say that $\{A_{\bm{n}}\}_{n}$ is a ($d$-level) locally Toeplitz sequence with symbol $v\otimes f$, and we write $\{A_{\bm{n}}\}_{n}\sim_{LT}v\otimes f$, if:

[TABLE]

We are now able to define generalized locally Toeplitz sequences.

Definition 2.15**.**

Let a $d$-level matrix-sequence $\{A_{\bm{n}}\}_{n}$, and a measurable function $\kappa:[0,1]^{d}\times[-\pi,\pi]^{d}\mapsto\mathbb{C}$ be given. Suppose that $\forall\ \epsilon>0$ there exists a finite number of $d$-level LT sequences $\{A_{\bm{n}}^{(i,\epsilon)}\}_{n}\sim_{LT}v_{i,\epsilon}\otimes f_{i,\epsilon}$, $i=1,\ldots,N_{\epsilon}$, such that as $\epsilon\rightarrow 0$:

[TABLE]

Then $\{A_{\bm{n}}\}_{n}$ is a $d$-level GLT sequence with symbol $\kappa$, and we write $\{A_{\bm{n}}\}_{n}\sim_{GLT}\kappa$.

The GLT class contains a wide range of matrix-sequences arising from various discretization methods of numerous differential equations. In the following Theorem we present some important properties of GLT sequences that will be used later in this paper. This is only a subset of the properties of multilevel GLT sequences, and the reader is referred to [34, 35] for a detailed derivation of all the results presented in this section. Given a measurable function $\kappa$, we denote its complex conjugate by $\bar{\kappa}$.

Theorem 2.16**.**

Let $\{A_{\bm{n}}\}_{n}$ and $\{B_{\bm{n}}\}_{n}$ be two d-level matrix-sequences and $\kappa,\ \xi:[0,1]^{d}\times[-\pi,\pi]^{d}\mapsto\mathbb{C}$ two measurable functions. Assume that $\{A_{\bm{n}}\}_{n}$ is a GLT sequence with symbol $\kappa$, while $\{B_{\bm{n}}\}_{n}$ a GLT sequence with symbol $\xi$. Then:

If $A_{\bm{n}}$ are Hermitian then $\{A_{\bm{n}}\}_{n}\sim_{\lambda}\kappa$. 2. 2.

$\{A^{*}_{\bm{n}}\}_{n}\sim_{GLT}\bar{\kappa}$. 3. 3.

$\{c_{1}A_{\bm{n}}+c_{2}B_{\bm{n}}\}_{n}\sim_{GLT}c_{1}\kappa+c_{2}\xi,\ \textnormal{for all }c_{1},\ c_{2}\in\mathbb{C}$. 4. 4.

$\{A_{\bm{n}}B_{\bm{n}}\}_{n}\sim_{GLT}\kappa\xi$. 5. 5.

If $\kappa\neq 0$ almost everywhere, then $\{A^{\dagger}_{\bm{n}}\}_{n}\sim_{GLT}\kappa^{-1}$. 6. 6.

Let a sequence of d-level matrix-sequences $\{B_{\bm{n},m}\}_{n}\sim_{GLT}\kappa_{m}$. Then, we have that $\{B_{\bm{n},m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{A_{\bm{n}}\}_{n}$ if and only if $\kappa_{m}\rightarrow\kappa$ in measure. 7. 7.

If $\{A_{\bm{n}}\}_{n}\sim_{GLT}\kappa$ and each $A_{\bm{n}}$ is Hermitian, then $\{f(A_{\bm{n}})\}_{n}\sim_{GLT}f(\kappa)$ for every continuous function $f:\mathbb{C}\mapsto\mathbb{C}$.

3 A Structure Preserving Method

In this section, we will derive an optimization method suitable for solving convex quadratic optimization problems, with linear constraints arising from the discretization of some continuous process. The assumption on the constraints is that the generated (multilevel) matrix-sequence is a GLT sequence. Let us consider the following generic Differential Equation (DE):

[TABLE]

where $\mathrm{D}$ denotes some linear differential operator associated with the DE, $\bm{x}$ is a $(d-1)$-dimensional spatial variable and $t\geq 0$ is the time variable. Since analytical solutions are not readily available for various differential operators, we discretize the previous equation given an arbitrary numerical method, and instead solve a sequence of linear systems of the form:

[TABLE]

with size $d_{n}\times d_{n}$ and $d_{n}=N(\bm{n})$, such that $\bm{n}\rightarrow\infty$ as $n\rightarrow\infty$.

Concerning the objective of the studied model, we assume that it may be summarized by a convex functional $\mathrm{J_{1}}(\mathrm{y}(\bm{x},t))$. Usually, such a functional measures the misfit between the state $\mathrm{y}(\bm{x},t)$ and some desired state $\mathrm{\bar{y}}(\bm{x},t)$, and we will focus our attention on this class of (inverse) problems. In other words, we expect that the discretized version of this functional will be of the form $\frac{1}{2}(y-\bar{y})^{*}J_{1}(y-\bar{y})$, with $J_{1}$ a symmetric positive (semi-)definite matrix. As is common in such problems, the linear systems in (3.1) usually admit more than one solution and hence a regularization functional is usually employed to guarantee that the chosen solution will have some desired properties, depending on the initial DE under consideration. In other words, we introduce a control variable $\mathrm{u}(\bm{x},t)$ which is linked to the state variable as follows:

[TABLE]

The size of the control is measured using some convex functional $\mathrm{J_{2}}(\mathrm{u}(\bm{x},t))$.

Finally, we allow further restrictions on the state and control variables in the form of inequality constraints (which depend on the problem under consideration). By combining all the previous, we obtain the following generic model that is studied in this paper:

[TABLE]

The problem is considered on a given compact space-time domain $\Omega\times(0,T)$, for some $T>0$, where $\Omega\subset\mathbb{R}^{d-1}$ has boundary $\partial\Omega$. The algebraic inequality constraints are assumed to hold a.e. on $\Omega\times(0,T)$. We further note that the restrictions ${\rm y_{a}},\ {\rm y_{b}},\ {\rm u_{a}}$, and ${\rm u_{b}}$ may take the form of constants, or functions in spatial and/or temporal variables. The boundary conditions are not specified since they do not affect the analysis in this section. Notice that problem (3.2) includes the case of equality-constrained optimization, by allowing unbounded restriction functions.

We discretize problem (3.2), using an arbitrary numerical method, to find an approximate solution by solving a sequence of optimization problems of the form:

[TABLE]

in which the associated matrices are of size $d_{n}\times d_{n}$, where $d_{n}=N(\bm{n})$ and $d_{n}\rightarrow\infty$ as $n\rightarrow\infty$. Notice that we only assume $J_{1_{\bm{n}}}$ and $J_{2_{\bm{n}}}$ to be symmetric positive semi-definite. Hence, the presented methodology is applicable to a wide range of convex quadratic programming problems. An entry $n_{j}$ of the multi-index $\bm{n}$ corresponds to the number of discretization points along dimension $j$, with $j\in\{1,\ldots,d\}$, where $n_{d}$ corresponds to the time dimension. Below, we summarize our assumptions for the associated matrices in problem (3.3).

Assumption 1**.**

Given the sequence of problems in (3.3), we assume that:

• •

The sequence $\{D_{\bm{n}}\}_{n}$ is a $d$-level matrix sequence with spectral norm uniformly bounded with respect to $n$, i.e. there exists a constant $C_{D}$ such that $\|D_{\bm{n}}\|\leq C_{D}$ for all $n$. Furthermore, there exists a measurable function $\kappa:[0,1]^{d}\times[-\pi,\pi]^{d}\mapsto\mathbb{C}$, which is the symbol of $\{D_{\bm{n}}\}_{n}$, so that $\{D_{\bm{n}}\}_{n}\sim_{GLT}\kappa$.

• •

The sequences $\{J_{1_{\bm{n}}}\}_{n}$ and $\{J_{2_{\bm{n}}}\}_{n}$ are two $d$-level matrix sequences, with uniformly bounded spectral norms with respect to $n$. Furthermore, there exist two measurable functions $\xi_{1},\ \xi_{2}:[0,1]^{d}\times[-\pi,\pi]^{d}\mapsto\mathbb{R}$, such that $\xi_{1}\geq 0$, $\xi_{2}\geq 0$, $\{J_{1_{\bm{n}}}\}_{n}\sim_{GLT}\xi_{1}$, and $\{J_{2_{\bm{n}}}\}_{n}\sim_{GLT}\xi_{2}$.

We note that a wide range of numerical discretizations of DEs satisfy this assumption (see [34, 35] for a plethora of applications). Notice also that the requirement that $\xi_{1}$ and $\xi_{2}$ are real and non-negative follows from the positive semi-definiteness of $J_{1_{\bm{n}}}$ and $J_{2_{\bm{n}}}$. Towards the end of this section we discuss how one could still apply the presented methodology successfully without requiring the GLT structure of the discretized objective function (i.e. by requiring only boundedness and convexity).

Before presenting the proposed optimization method for solving problems of the form (3.3), we note a negative result concerning a large class of optimization methods. More specifically, problems like (3.3) are often solved using an Interior Point Method (IPM), or some Active-Set (AS) type of method. However, such problems are usually highly structured, and this structure must be exploited, given that the problem size increases indefinitely as one refines the discretization. Obviously, any AS method would fail in maintaining the structure, as only a subset of the constraints of (3.3) is considered at each AS iteration and hence the structure of the AS sub-problems will be unknown. In fact, any optimization method whose sub-problems arise by projecting the variables of the problem in a subspace would face this issue.

On the other hand, IPMs deal with the inequality constraints by introducing logarithmic barriers in the objective (see for example [38]). Then, at every IPM iteration, one forms the optimality conditions of the barrier sub-problem, and approximately solves them using Newton’s method. If we assume that there exists a symbol $f$ which describes the asymptotic eigenvalue distribution of the sequence of Hessian matrices of the logarithmic barrier, we arrive at a contradiction. Indeed, the sequence of Hessian matrices arising from the logarithmic barriers introduced by IPM are not s.u.. This in turns contradicts the assumption that $f$ is the symbol of this matrix sequence, since if an arbitrary matrix sequence is such that $\{L_{n}\}_{n}\sim_{\sigma}f$, for some measurable function $f$, then $\{L_{n}\}_{n}$ must be s.u. (see [34, Section 9–S1]). In particular, any GLT sequence is s.u., and hence the sequence of Hessian matrices of the logarithmic barrier functions cannot be a GLT sequence. As a consequence, the system matrix of the optimality conditions of each barrier sub-problem, within the IPM, will not be in the GLT class.

3.1 Alternating Direction Method of Multipliers

In order to overcome the previous issues, we propose the use of an alternating direction method of multipliers (see [5, Section 5] and the references therein), which separates the equality from the inequality constraints, thus allowing us to preserve the structure found in the matrices associated with (3.3). We should mention here that while ADMM allows us to retain the underlying structure of the problem, it comes at a cost. It is well-known (see e.g. [5]) that ADMM leads to relatively slow convergence and hence is not suitable for finding very accurate solutions. Nevertheless, a 4-digit accurate solution can generally be found in reasonable CPU time. Furthermore, the linear system solved at each ADMM iteration does not change, and hence, if a suitable preconditioner exploiting the problem structure is found, it only needs to be computed once. Finally, linear convergence can also be shown, under certain assumptions on the problem under consideration (such as strong convexity, see [20]).

We begin by rewriting problem (3.3), after introducing some auxiliary variables $z_{y_{\bm{n}}},z_{u_{\bm{n}}}$ of size $N(\bm{n})$:

[TABLE]

Next, we define the augmented Lagrangian function corresponding to (3.4):

[TABLE]

where $p_{\bm{n}},\ w_{y_{\bm{n}}}$, and $w_{u_{\bm{n}}}$ are the dual variables corresponding to each of the equality constraints of (3.4). An ADMM applied to model (3.4) is given in Algorithm 1. We omit specific details of the algorithm. The reader is referred to [5] for a basic proof of convergence of Algorithm 1, as well as a detailed overview of ADMM. For a convergence proof for the case where complex variables and matrices are allowed, the reader is referred to [49]. Linear convergence of a generalization of this algorithm, under certain assumptions, can be found in [20] and the references therein. We should note that the step-length $\rho$ in (3.6c) and (3.6d) plays an important role in the convergence behavior of ADMM, and in fact, convergence of Algorithm 1 is guaranteed for any $\rho\in(0,\frac{\sqrt{5}+1}{2})$ (see [36]).

One can easily observe that the most challenging step of Algorithm 1, is that of solving (3.6a). The optimality conditions of (3.6a), at iteration $j$, read as follows:

[TABLE]

where

[TABLE]

Solving (3.7) directly is not a good idea in our case, since its coefficient matrix is not expected to be cheap or convenient to work with. Instead, we can merge steps (3.6a) and (3.6c) to obtain a more flexible saddle point system. More specifically, to take (3.6c) into account, we substitute $p_{\bm{n}}=p_{\bm{n}}^{j}+\frac{\rho}{\delta}(D_{\bm{n}}y_{\bm{n}}+u_{\bm{n}}-g_{\bm{n}})$ into (3.7), and the optimality conditions of (3.6a) and (3.6c) can then be written as:

[TABLE]

At this point, we have to decide how to solve (3.8). For simplicity of exposition, we present here only one way of solving system (3.8), by forming the normal equations and then employing the Preconditioned Conjugate Gradient method (PCG) to solve the resulting positive definite system, assuming that its $(2,2)$ block will be easily invertible. We note that the developments in this section hold for any Schur complement of the matrix in (3.8) (the choice of which Schur complement to use heavily depends on the problem under consideration). The case where neither the $(1,1)$ nor the $(2,2)$ block is easily invertible, will be treated at the end of this section. Pivoting the second and then the third block equation of this system, yields:

[TABLE]

and the resulting normal equations read as follows:

[TABLE]

Finally, we should mention that problem (3.6b) of Algorithm 1 is trivial, as it admits a closed form solution. More specifically, we perform the optimization by ignoring the box constraints and then projecting the solution onto the box.

In what follows, using Assumption 1, we present some results concerning the asymptotic behavior of the matrix sequence $\{S_{\bm{n}}\}_{n}$, by making use of the Theorems presented in the previous section. The latter is produced by refining an arbitrary discretization applied to (3.2) (assuming it satisfies Assumption 1), employing ADMM to the discretized problem, and forming a certain Schur complement of the joint optimality conditions of (3.6a) and (3.6c). The solution of (3.9) delivers the solution to (3.8), and the remaining ADMM sub-problems can be trivially solved in $O(N(\bm{n}))$ operations. Following practical applications, we assume $\delta$ and $\rho$ to be $\Theta(1)$ and constant along the iterations of ADMM (usually $\delta\in[0.01,100]$ and $\rho\in[1,1.618]$) .

Theorem 3.1**.**

Given Assumption 1, and the sequence $\{S_{\bm{n}}\}_{n}$, with $S_{\bm{n}}$ given in (3.9), we have that there exists a measurable function $\tau:[0,1]^{d}\times[-\pi,\pi]^{d}\mapsto\mathbb{R}$ such that $\tau\geq 0$, $\tau\neq 0$ a.e., and $\{S_{\bm{n}}\}_{n}\sim_{GLT}\tau$. Moreover, $S_{\bm{n}}$ are Hermitian positive definite, $\{S_{\bm{n}}\}_{n}\sim_{\lambda}\tau$, and $\{S^{-1}_{\bm{n}}\}_{n}\sim_{\lambda}\tau^{-1}$.

Proof.

Let Assumption 1 hold. Then, we have that there exist three measurable functions $\kappa,\ \xi_{1},\ \xi_{2}:[0,1]^{d}\times[-\pi,\pi]^{d}\mapsto\mathbb{C}$, such that $\xi_{1}\geq 0$, $\xi_{2}\geq 0$, $\{J_{1_{\bm{n}}}\}_{n}\sim_{GLT}\xi_{1}$, $\{J_{2_{\bm{n}}}\}_{n}\sim_{GLT}\xi_{2}$, and $\{D_{\bm{n}}\}_{n}\sim_{GLT}\kappa$. Furthermore, we can notice that, for any constant $C>0$, $\{CI_{\bm{n}}\}_{n}\sim_{GLT}C$, where $C$ can be considered as a positive constant function (e.g. as a constant on the domain $[-\pi,\pi]^{d}$, generating a diagonal Toeplitz matrix). This, combined with Theorem 2.16 (conditions (2.)–(4.)), yields that:

[TABLE]

Similarly, from Theorem 2.16 (conditions (2.)–(5.)), we have that:

[TABLE]

where we used that $\bar{\kappa}\kappa=|\kappa|^{2}$. Setting \tau=\rho(\xi_{1}+\delta^{-1})+|\kappa|^{2}\bigg{(}\big{(}\rho(\xi_{2}+\delta^{-1})\big{)}^{-1}+\frac{\delta}{\rho}\bigg{)}^{-1} and noticing that $\tau>0$ completes the proof. ∎

Subsequently we present some possible approaches that could allow one to take advantage of the structure preserving property of ADMM. In particular, three possible ways of exploiting the preserved structure are discussed here. However, other approaches could be possible. For this analysis, we will make use of the following proposition:

Proposition 3.2**.**

Let Assumption 1 hold. Then, there exist sequences of $d$-level matrix-sequences $\{\{\tilde{D}_{\bm{n},m}\}_{n}\}_{m},$ $\{\{\tilde{J}_{1_{\bm{n},m}}\}_{n}\}_{m},$ $\{\{\tilde{J}_{2_{\bm{n},m}}\}_{n}\}_{m}$, with uniformly bounded spectral norms with respect to $n$ and $m$, and sequences of measurable functions $\{\kappa_{m}\}_{m},$ $\{\xi_{1_{m}}\}_{m}$, and $\{\xi_{2_{m}}\}_{m}$ such that $\kappa_{m},\ \xi_{1_{m}},\ \xi_{2_{m}}:[0,1]^{d}\times[-\pi,\pi]^{d}\mapsto\mathbb{C}$, $\xi_{1_{m}},\ \xi_{2_{m}}$ are real a.e., non-negative, and:

• •

$\{\tilde{D}_{\bm{n},m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{D_{\bm{n}}\}_{n}$, $\{\tilde{D}_{\bm{n},m}\}_{n}\sim_{GLT}\kappa_{m}$, with $\kappa_{m}\rightarrow\kappa$ in measure,

• •

$\{\tilde{J}_{1_{\bm{n},m}}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{J_{1_{\bm{n}}}\}_{n}$, $\{\tilde{J}_{1_{\bm{n},m}}\}_{n}\sim_{GLT}\xi_{1_{m}}$, with $\xi_{1_{m}}\rightarrow\xi_{1}$ in measure,

• •

$\{\tilde{J}_{2_{\bm{n},m}}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{J_{2_{\bm{n}}}\}_{n}$, $\{\tilde{J}_{2_{\bm{n},m}}\}_{n}\sim_{GLT}\xi_{2_{m}}$, with $\xi_{2_{m}}\rightarrow\xi_{2}$ in measure.

Proof.

The proof can be found in [34, Theorem 8.6] for the unilevel case and [35, Theorem 5.6] for the multilevel case. ∎

For the rest of this section we will assume that we have such sequences of $d$-level GLT sequences available, satisfying the conditions stated in Proposition 3.2. We further assume that these approximate sequences are comprised of matrices that are easy to compute and invert (whenever possible).

3.2 Schur complement approximations

In what follows we present various Schur complement approximations that could potentially serve as preconditioners inside PCG, for solving systems of the form of (3.9) (or any other Schur complement of system (3.8)). The viability of each of the following approaches depends on the structure of the problem, as well as the choice of the discretization. We note that the different approaches are presented for completeness, as well as an indicator of the generality of the proposed methodology. In particular, as the convergence behavior of ADMM does not depend on the choice of preconditioner (assuming that the PCG converges to a desired accuracy), we will only make use of one of the following Schur complement approximations when presenting computational results.

3.2.1 A Schur complement block approximation

Given three sequences of GLT sequences $\{\{\tilde{D}_{\bm{n},m}\}_{n}\}_{m}$, $\{\{\tilde{J}_{1_{\bm{n},m}}\}_{n}\}_{m},$ $\{\{\tilde{J}_{2_{\bm{n},m}}\}_{n}\}_{m}$, satisfying the conditions of Proposition 3.2, we define the following approximation for the matrix in (3.9):

[TABLE]

Theorem 3.3**.**

Let Assumption 1 hold, and assume that we have available the sequences $\{\{\tilde{D}_{\bm{n},m}\}_{n}\}_{m},$ $\{\{\tilde{J}_{1_{\bm{n},m}}\}_{n}\}_{m},$ and $\{\{\tilde{J}_{2_{\bm{n},m}}\}_{n}\}_{m}$, satisfying the conditions of Proposition 3.2. By defining $\tilde{S}_{\bm{n},m}$ as in (3.10), we have:

• •

$\{\tilde{S}_{\bm{n},m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{S_{\bm{n}}\}_{n}$, $\{\tilde{S}_{\bm{n},m}\}_{n}\sim_{GLT}\tau_{m}$, and $\tau_{m}\rightarrow\tau$ in measure, where $\tau$ is given in Theorem 3.1 and:

[TABLE]

• •

The sequence $\{\tilde{S}^{-1}_{\bm{n},m}S_{\bm{n}}\}_{n}$ is weakly clustered at $1$.

• •

For any $n,\ m$, the eigenvalues of $\tilde{S}^{-1}_{\bm{n},m}S_{\bm{n}}$ lie in the interval \bigg{[}\frac{1}{C_{s}},C_{s}\bigg{]}, where $C_{s}$ is a positive constant uniformly bounded with respect to $n$ and $m$.

Proof.

For the first condition, by Proposition 3.2 as well as Theorem 2.12 (conditions (1.)–(4.)), we get that $\{\tilde{S}_{\bm{n},m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{S_{\bm{n}}\}_{n}$. Using Proposition 3.2 again, this time combined with Theorem 2.16 (conditions (2.)–(5.)) and Lemma 2.4, yields that $\tau_{m}$ is given by (3.11), and $\tau_{m}\rightarrow\tau$ in measure.

For the second condition, we firstly note that the sequence under consideration is Hermitian and positive definite by Assumption 1. Then, using that $\{\tilde{S}_{\bm{n},m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{S_{\bm{n}}\}_{n}$ implies that, for every $m$, there exists $n_{m}$ such that for all $n\geq n_{m}$:

[TABLE]

where $n_{m}$, $c(m)$ and $\omega(m)$ depend only on $m$ and are such that:

[TABLE]

By assumption, it is easy to see that any $\tilde{S}_{\bm{n},m}$ has a spectral norm uniformly bounded in $n$ and $m$. Furthermore, since $\delta=\Theta(1)$ and $\rho=\Theta(1)$, we observe that $\|\tilde{S}^{-1}_{\bm{n},m}\|$ is also uniformly bounded in $n$. Hence, by multiplying both sides of (3.12) by $\tilde{S}^{-1}_{\bm{n},m}$:

[TABLE]

where $\tilde{R}_{\bm{n},m}=\tilde{S}^{-1}_{\bm{n},m}R_{\bm{n},m}$ (thus $\textnormal{rank}(\tilde{R}_{\bm{n},m})\leq\textnormal{rank}(R_{\bm{n},m})\leq c(m)d_{n}$) and $\tilde{N}_{\bm{n},m}=\tilde{S}^{-1}_{\bm{n},m}N_{\bm{n},m}$ (and hence $\|\tilde{N}_{\bm{n},m}\|\leq C\omega(m)$, for some constant $C>0$, independent of $m$ and $n$). This, along with the definition of a weak cluster in Definition 2.3, proves the second condition.

For the third condition, let us take some constant $C_{\dagger}$ of $O(1)$, such that:

[TABLE]

We know that such a constant exists by Assumption 1. Then, we have that $\lambda_{\min}(S_{\bm{n}})\geq\frac{\rho}{\delta}$ and $\lambda_{\max}(S_{\bm{n}})\leq\rho C_{\dagger}+\frac{\rho}{\delta}+\frac{\rho}{\delta}C_{\dagger}^{2}$, for any $n$. The exact same bounds hold also for $\tilde{S}_{\bm{n},m}$, for every $n$ and $m$. Using these bounds, we can easily show that:

[TABLE]

for all $n,\ m$. Upon noticing that $\delta=\Theta(1)$ and $\rho=\Theta(1)$, there exists a constant $C_{s}=C_{\dagger}^{2}+\delta C_{\dagger}+1$ uniformly bounded with respect to $n$, satisfying the third condition of the Theorem. ∎

Remark 3.1**.**

Notice that in order to obtain a strong clustering at 1, we would have to employ some extra assumptions. In particular, we would have to require that the sequences given in Assumption 1 are strongly clustered in the essential range of their symbols, which in turn are required to be different from zero a.e.. Furthermore, we would have to assume that the condition in (3.12) is such that $c(m)d_{n}=O(1)$.

Remark 3.2**.**

While Assumption 1 holds for a wide range of problems, and one is able to find easily computable sequences satisfying the conditions in Proposition 3.2, it is not often the case that the preconditioner in (3.10) is easy to compute or invert. If this is the case, then Theorem 3.3 guarantees that such a preconditioner will provide a weak cluster of the eigenvalues of the preconditioned matrix at $1$. We note that while this is not optimal, it is expected to be good enough. This is because of the penalty parameter introduced by ADMM (i.e. $\delta=\Theta(1)$), which (along with the assumption that the involved matrices are uniformly bounded in $n$) guarantees that the normal equations matrix will be relatively well-conditioned, and hence PCG will converge in a number of iterations independent of the grid size (however, possibly depending on the conditioning of the problem matrix as well as the problem parameters).

The use of the preconditioner in (3.10) becomes more obvious in the following example. If the $d$-level approximating matrix sequences $\{\{\tilde{D}_{\bm{n},m}\}_{n}\}_{m},$ $\{\{\tilde{J}_{1_{\bm{n},m}}\}_{n}\}_{m},$ and $\{\{\tilde{J}_{2_{\bm{n},m}}\}_{n}\}_{m}$ satisfy the conditions of Proposition 3.2, and belong to the set of $d$-level circulant matrices of size $N(\bm{n})\times N(\bm{n})$, i.e. $\mathcal{C}_{\bm{n}}$, then the preconditioner in (3.10) will be cheap to form and store, using the fast Fourier transform (requiring O\big{(}N(\bm{n})\log(N(\bm{n}))\big{)} operations and O\big{(}N(\bm{n})\big{)} memory). This is because $\mathcal{C}_{\bm{n}}$ is a commutative ring under standard matrix addition and multiplication (see Theorem 2.10).

3.2.2 A matching Schur complement approximation

As mentioned earlier, many approximating sequences based on the GLT theory would not allow for an easy computation or storage of the preconditioner in (3.10). While the numerical results of this paper will not focus on this case, we present an alternative to the preconditioner in (3.10), which could allow one to use various GLT approximations for the blocks of the matrix in (3.8), and form an easily computable Schur complement approximation for a matrix of the form of (3.9).

In what follows, we define a Schur complement approximation based on the matching strategy, which was proposed in [59] and has been applied in a wide range of applications (e.g. [23, 58, 60]). While this approach can be very general, it is particularly effective under some additional assumptions imposed on problem (3.3). More specifically, we study the properties of this approximation using the GLT theory, and give certain assumptions under which such an approach would be optimal.

Given three sequences $\{\{\tilde{D}_{\bm{n},m}\}_{n}\}_{m}$, $\{\{\tilde{J}_{1_{\bm{n},m}}\}_{n}\}_{m}$, $\{\{\tilde{J}_{2_{\bm{n},m}}\}_{n}\}_{m}$, satisfying the conditions of Proposition 3.2, we define the following matrix:

[TABLE]

using which we can define an approximation for the matrix in (3.9) as

[TABLE]

For simplicity of exposition let us define the following matrices:

[TABLE]

[TABLE]

Under Assumption 1, we have that $\{M_{1_{\bm{n}}}\}_{n}\sim_{GLT}\rho(\xi_{1}+\delta^{-1})$, $\{M_{2_{\bm{n}}}\}_{n}\sim_{GLT}(\rho(\xi_{2}+\delta^{-1}))^{-1}+\frac{\delta}{\rho})$, and $\tilde{M}_{1_{\bm{n}}}\sim_{GLT}\rho(\xi_{1_{m}}+\delta^{-1})$ with $\xi_{1_{m}}\rightarrow\xi_{1}$ in measure, while $\tilde{M}_{2_{\bm{n}}}\sim_{GLT}(\rho(\xi_{2_{m}}+\delta^{-1}))^{-1}+\frac{\delta}{\rho})$, where $\xi_{2_{m}}\rightarrow\xi_{2}$ in measure. Further, notice that all four previous matrix-sequences are comprised of Hermitian and positive definite matrices, each of which admits a square root.

Lemma 3.4**.**

Let $\bm{n}\in\mathbb{N}^{d}$ be a $d$-index and $\{A_{\bm{n}}\}_{n}$ be a multilevel matrix-sequence with $A_{\bm{n}}$ being Hermitian positive definite of size $N(\bm{n})\times N(\bm{n})$ and $\{A_{\bm{n}}\}_{n}\sim_{GLT}\chi$, where $\chi$ is a measurable function $\chi:[0,1]^{d}\times[-\pi,\pi]^{d}\mapsto\mathbb{R}$ such that $\chi\geq 0$ and $\chi\neq 0$ a.e.. Then, the matrices $A_{\bm{n}}$ ($A^{-1}_{\bm{n}}$, respectively) admit a square root $A_{\bm{n}}^{\frac{1}{2}}$ ($A^{-\frac{1}{2}}_{\bm{n}}$, respectively), such that $\{A^{\frac{1}{2}}_{\bm{n}}\}_{n}\sim_{GLT}\chi^{\frac{1}{2}}$ ($\{A^{-\frac{1}{2}}_{\bm{n}}\}_{n}\sim_{GLT}\chi^{-\frac{1}{2}}$, respectively).

Proof.

Let the function $f:(0,\infty)\mapsto(0,\infty)$, be defined as $f(x)=x^{\frac{1}{2}}$. Then, from Theorem 2.16 (condition (7.)), we have that $\{f(A_{\bm{n}})\}_{n}\sim_{GLT}f(\chi)$, where $f(A_{\bm{n}})$ is interpreted as a matrix function applied to the eigenvalues of matrix $A_{\bm{n}}$. ∎

Theorem 3.5**.**

Let Assumption 1 hold, and assume that we have available the sequences $\{\{\tilde{D}_{\bm{n},m}\}_{n}\}_{m}$, $\{\{\tilde{J}_{1_{\bm{n},m}}\}_{n}\}_{m},$ and $\{\{\tilde{J}_{2_{\bm{n},m}}\}_{n}\}_{m}$, satisfying the conditions of Proposition 3.2. By defining $\hat{S}_{\bm{n},m}$ as in (3.13), we have:

• •

$\{\hat{S}_{\bm{n},m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{S_{\bm{n}}+E_{\bm{n}}\}_{n}$, where

[TABLE]

with

[TABLE]

Furthermore, $\{\hat{S}_{\bm{n},m}\}_{n}\sim_{GLT}\tau_{m}+\epsilon_{m}$ and $\tau_{m}+\epsilon_{m}\rightarrow\tau+\epsilon$ in measure, where $\tau$ is defined in Theorem 3.1, and $\tau_{m},\ \epsilon_{m}$, $\epsilon$ are measurable functions having the same domain as $\tau$. If $\tilde{E}_{\bm{n},m}\coloneqq\tilde{M}^{\frac{1}{2}}_{1_{\bm{n},m}}\tilde{M}^{-\frac{1}{2}}_{2_{\bm{n},m}}\tilde{D}_{\bm{n},m}+\tilde{D}^{*}_{\bm{n},m}\tilde{M}^{-\frac{1}{2}}_{2_{\bm{n},m}}\tilde{M}^{\frac{1}{2}}_{1_{\bm{n},m}}$ is positive semi-definite for all $m$ and $n$, then the sequence of preconditioned normal equations’ matrices is such that:

[TABLE]

and there exist positive constants $C_{\dagger_{1}}$, $C_{\dagger_{2}}$, independent of $n,\ m$, such that $\lambda(\hat{S}_{\bm{n},m}^{-1}S_{\bm{n}})\in[C_{\dagger_{1}},C_{\dagger_{2}}]$, for all $n,\ m$.

• •

If the matrix sequences $\{J_{1_{\bm{n}}}\}_{n}$, $\{J_{2_{\bm{n}}}\}_{n}$, and $\{D_{\bm{n}}\}_{n}$ are such that $J_{1_{\bm{n}}}$ and $J_{2_{\bm{n}}}$ are scaled identities or zero matrices, while $D_{\bm{n}}+D^{*}_{\bm{n}}$ is Hermitian positive semi-definite, then the matrix-sequence $\{\hat{S}_{\bm{n},m}^{-1}S_{\bm{n}}\}_{n}$ is weakly clustered at $[\frac{1}{2},1]$. If furthermore the matrix-sequence $\{\tilde{D}_{\bm{n},m}^{-1}D_{\bm{n}}\}_{n}$ is strongly clustered at $1$, then the matrix-sequence $\{\hat{S}_{\bm{n},m}^{-1}S_{\bm{n}}\}_{n}$ is strongly clustered at $[\frac{1}{2},1]$.

Proof.

Firstly, notice that from (3.13) we obtain the following expression:

[TABLE]

where $\tilde{S}_{\bm{n},m}$ is defined as in (3.10). Then, the first part of the Theorem can be proved by employing Lemma 3.4 and by performing a similar analysis to that of the proof of the first and third conditions of Theorem 3.3. For brevity, the latter is omitted.

We proceed by proving the second condition. Notice that if $J_{1_{\bm{n}}}$ and $J_{2_{\bm{n}}}$ are scaled identities or zero matrices (the latter being mostly of theoretical interest), then we can represent them exactly, that is $\tilde{J}_{1_{\bm{n},m}}=J_{1_{\bm{n}}}$ and $\tilde{J}_{2_{\bm{n},m}}=J_{2_{\bm{n}}}$, for all $m$ and $n$. The latter implies that $M_{1_{\bm{n}}}$, $M_{2_{\bm{n}}}$ are scaled identities and we can write $M_{\bm{n}}=M_{1_{\bm{n}}}=\frac{1}{c_{s}}M_{2_{\bm{n}}}$, for some positive constant $c_{s}$. We define the following matrix:

[TABLE]

Following exactly the developments in [60, Theorem 4.1] (since $D+D^{*}\succeq 0$), we can consider the generalized eigenproblem $\bar{S}^{-1}_{\bm{n}}S_{\bm{n}}x=\mu x$, and show that $\lambda(\bar{S}^{-1}_{\bm{n}}S_{\bm{n}})\in[\frac{1}{2},1]$, where $S_{\bm{n}}$ is defined as in (3.9), $\mu$ is an arbitrary eigenvalue of the preconditioned matrix $\bar{S}^{-1}_{\bm{n}}S_{\bm{n}}$ and $x$ the corresponding eigenvector.

Let us now notice that by Assumption 1, the matrix-sequence $\{\bar{S}_{\bm{n}}\}_{n}$ is a GLT sequence. In particular, it is easy to see that $\{S_{\bm{n}}+E_{\bm{n}}\}_{n}\equiv\{\bar{S}_{\bm{n}}\}_{n}$ and hence its symbol is $\tau+\epsilon$, where $\epsilon$ is defined in the first condition of this Theorem. Again, from the first condition of this Theorem, we have that the preconditioner defined in (3.13) is such that $\{\hat{S}_{\bm{n},m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{S_{\bm{n}}+E_{\bm{n}}\}_{n}\equiv\{\bar{S}_{\bm{n}}\}_{n}$, and $\{\hat{S}_{\bm{n},m}\}_{n}\sim_{GLT}\tau_{m}+\epsilon_{m}$ with $\tau_{m}+\epsilon_{m}\rightarrow\tau+\epsilon$ in measure. Then, from Theorem 2.12 we know that $\{\hat{S}^{-1}_{\bm{n},m}\}_{n}\xrightarrow{\textnormal{a.c.s.}}\{\bar{S}^{-1}_{\bm{n}}\}_{n}$. Using Definition 2.11, we have that for all $n\geq n_{m}$, we can write:

[TABLE]

where $n_{m}$, $c(m)$ and $\omega(m)$ depend only on $m$ and are such that:

[TABLE]

In view of the previous, we can analyze the sequence $\{\bar{S}^{-1}_{\bm{n}}S_{\bm{n}}-\hat{S}^{-1}_{\bm{n},m}S_{\bm{n}}\}_{n}$ as follows:

[TABLE]

where $\textnormal{rank}(R_{\bm{n},m}S_{\bm{n}})\leq\textnormal{rank}(R_{\bm{n},m})\leq c(m)N(\bm{n})$ and $\|N_{\bm{n},m}S_{\bm{n}}\|\leq\omega(m)\|S_{\bm{n}}\|=\Theta(\omega(m))$. In other words, we have that $\{\bar{S}^{-1}_{\bm{n}}S_{\bm{n}}-\hat{S}^{-1}_{\bm{n},m}S_{\bm{n}}\}_{n}$ is weakly clustered at zero. Furthermore, as $\lambda(\bar{S}^{-1}_{\bm{n}}S_{\bm{n}})\in[\frac{1}{2},1]$, we conclude that $\{\hat{S}^{-1}_{\bm{n},m}S_{\bm{n}}\}_{n}$ is weakly clustered at $[\frac{1}{2},1]$.

Finally, if we assume that $\{\tilde{D}^{-1}_{\bm{n},m}D_{\bm{n}}\}_{n}$ is strongly clustered at $1$, and by noting that $M_{\bm{n}}$ is a scaled identity (and hence $\tilde{M}_{\bm{n},m}=M_{\bm{n}},$ for all $n,\ m$), we can conclude that (3.14) holds for $R_{\bm{n},m}$ such that $\textnormal{rank}(R_{\bm{n},m})=O(1)$. By employing a similar methodology as before, this yields that $\{\hat{S}^{-1}_{\bm{n},m}S_{\bm{n}}\}_{n}$ is strongly clustered at $[\frac{1}{2},1]$. ∎

Remark 3.3**.**

Let us now briefly discuss the applicability of the preconditioner in (3.13). Firstly, it is important to notice that such a preconditioner is generally only sensible when the approximating matrices $\tilde{J}_{1_{\bm{n},m}}$ and $\tilde{J}_{2_{\bm{n},m}}$ are diagonal, circulant or zero. If this is not the case, we discuss a remedy in Section 3.3.2. In many applications of interest, the preconditioner $\tilde{D}_{\bm{n},m}$ has a diagonal times a multilevel banded Toeplitz structure (e.g. [17, 53]). The application of the preconditioner in (3.13) consists in a single (LU or, if applicable, Cholesky) factorization of $\hat{D}_{\bm{n},m}$ at the beginning of the optimization, and subsequently two backward solves for every ADMM iteration. Such an approach should be feasible, in terms of memory and computational requirements, as long as the problem dimensions are not very large and the bandwidth of the matrix $\tilde{D}_{\bm{n},m}$ is small. For high-dimensional problems, one could employ an incomplete or specialized factorization (e.g. [17]), possibly assisted by suitable low-rank updates, if necessary. In some cases, replacing the factorization with a specialized iterative solver (such as a multigrid method as in [53]) could be beneficial. However, it is important to note that factorization (complete or incomplete) needs to be computed only once. An alternative employing low-rank approximations of the associated matrices is discussed in the following Remark. The suitability of each of the aforementioned approaches depends heavily on the problem under consideration.

Remark 3.4**.**

As mentioned earlier, the proposed preconditioner in (3.13) allows one to use a variety of approximations for the blocks of the matrix in (3.8), based on the GLT theory, under certain conditions (which hold for a wide class of problems, such as the problem considered in Section 4). In fact, this preconditioner can be seen as an approximation of the preconditioner in (3.10), which in turn has limited applicability unless the approximating blocks have a multilevel circulant structure.

The limitations of preconditioner (3.13) depend on the problem under consideration. In particular, if the assumptions of the second condition of Theorem 3.5 hold, then it can serve as a basis for constructing easily computable optimal preconditioners. Furthermore, if the aforementioned assumptions hold for the problem under consideration, one might be able to use a tensor product approach with low-rank approximations of the matrices in (3.8) to solve problem (3.3). Such solvers can be extremely effective, allowing one to solve high-dimensional problems, however they tend to require that various features of the problem (e.g. initial conditions, desired state, boundary conditions, discrete solutions) are approximated in a low-rank format, which is not always the case. The proposed approach would allow one to create a rather general tensor product solver for inverse problems measuring the discrepancy of the state variable $\rm{y}$ from a desired state $\rm{\bar{y}}$ as well as the size of the control $\rm{u}$ in the $L^{2}$-norm, where the structure of the problem allows this. Such solvers have been proposed in [23, 40] for the equality constrained case, and hence the proposed methodology could allow one to further generalize these approaches. Additionally, many low-rank tensor product solvers in the literature require that the objective function has a scaled identity Hessian. This can be alleviated here, by making use of the generalized ADMM presented in Section 3.3.2, alongside the preconditioner in (3.13).

3.2.3 Element-wise Schur complement approximation

In the context of finite element methods, a popular preconditioner is the so-called element-wise (also known as additive or element-by-element) Schur complement approximation. As this approach has been analyzed multiple times, we only mention it here as a viable alternative for preconditioning the normal equations in (3.9) and refer the interested reader to the available literature. In particular, such preconditioners have been analyzed using the GLT theory in [27, 28]. An analysis for general problems can be found in [55] and the references therein. These preconditioners can be very effective (in fact optimal under reasonable and general assumptions). Furthermore, they can efficiently be implemented in a parallel environment, allowing one to solve huge-scale problem instances (see e.g. [26, 28]).

3.3 General quadratic objective function

As we stressed earlier, it could be the case that both $J_{1_{\bm{n}}}$ and $J_{2_{\bm{n}}}$ are general positive semi-definite matrices, whose inverses (if they exist) are expensive to compute. As a consequence, the normal equations could be prohibitively expensive to form. In order to tackle such problems, we propose two alternatives. The former simply avoids forming the normal equations and solves (3.8) instead, using an appropriate Krylov subspace method. The latter approach generalizes the algorithmic framework in Algorithm 1, allowing us to simplify the resulting sub-problems. Then, the simplified sub-problems can be solved using PCG alongside any of the previously presented Schur complement approximations.

3.3.1 A saddle point approximation

In many applications, forming a Schur complement of system (3.8) would be very costly. Instead, one could solve system (3.8), which can be seen as a regularized saddle point system. Among many other iterative methods, one could employ preconditioned MINRES to solve systems of this form. The aforementioned method allows only the use of a positive definite preconditioner, hence, many block preconditioners for (3.8) are not applicable. For instance, block-triangular preconditioners, motivated by the work in [41, 54], would generally require a non-symmetric solver such as GMRES [64]. However, block-diagonal preconditioners have been shown to be very effective and efficient in practice for systems of the form of (3.8) (see for example [2, 57, 69]). To that end, we can define the following positive definite block-diagonal preconditioner:

[TABLE]

where $\tilde{S}_{\bm{n},m}$ can be defined as in (3.10) or as in (3.13) (and indeed any other suitable Schur complement approximation), assuming that we have available two sparse sequences of $d$-level GLT sequences $\{\{\tilde{J}_{1_{\bm{n},m}}\}_{n}\}_{m}$ and $\{\{\tilde{J}_{2_{\bm{n},m}}\}_{n}\}_{m}$, satisfying the conditions of Proposition 3.2. We note that preconditioners similar to (3.15) have been analyzed multiple times in the literature and hence such an analysis is omitted here (see for example [2, 57, 61, 68, 69]). It is important to notice that the quality of the preconditioner in (3.15) depends heavily on the quality of the Schur complement approximation, as well as on the approximations of the $(1,1)$ and $(2,2)$ blocks of the matrix in (3.8), which can be computed by making use of the GLT theory.

3.3.2 Generalized ADMM

Instead of solving the saddle point system in (3.8), one could derive the following generalized ADMM algorithm, as described in Algorithm 2. The following methodology is presented for completeness and is focused on the case where all the associated matrices as well as state and control variables are real. One could apply it to the complex case, however, in that case the theory derived in [20] to support such methods, would no longer hold.

There are two major differences between Algorithm 1 and Algorithm 2. In the latter method, we have added an extra proximal term in problem (3.6a), which belongs to the class of Bregman distances, and indeed is produced by the Bregman function $\|\cdot\|_{R}$, where $R=R_{y}\oplus R_{u}$, $R_{y}\succ 0$, and $R_{u}\succ 0$. For a detailed derivation of proximal methods using Bregman distances, the reader is referred to [30] and the references therein. The second difference is that, in the general case, Algorithm 2 requires that the step-size $\rho$ lies in a smaller interval than that allowed in Algorithm 1. In fact, the allowed values for $\rho$ depend on the choice of $R_{y}$ and $R_{u}$. We refer the reader to [20] for a more general derivation of methods similar to Algorithm 2, in which a precise condition is given for the maximum allowed values of $\rho$, so that the method converges globally. Furthermore, the authors in [20] prove linear convergence of the method under different sets of conditions, one of which requires that $J_{1_{\bm{n}}}\succ 0$ and $J_{2_{\bm{n}}}\succ 0$.

In light of the previous discussion, we can choose:

[TABLE]

where $c_{y},\ c_{u}>0$ are such that $R_{y_{\bm{n}}}\succ 0$, $R_{u_{\bm{n}}}\succ 0$. With these choices of $R_{y_{\bm{n}}}$ and $R_{u_{\bm{n}}}$, the optimality conditions of (3.16a) and (3.16c) involve the coefficient matrix:

[TABLE]

As one can easily observe, the normal equations operator of (3.17) can be efficiently applied to a vector, and all the previously presented Schur complement approximations can be used within PCG to accelerate the solution of the new simplified sub-problems. Furthermore, notice that this way we can ensure that the $(1,1)$ and $(2,2)$ blocks of the matrix in (3.17) are scaled identities, and hence the preconditioner in (3.13) can be particularly effective (see Theorem 3.5).

We should note at this point that similar methodologies can be employed to enforce certain structure on the associated matrices of problem (3.4). While this can be very effective in some cases, by making the ADMM sub-problems easy to solve, it should be used with caution. On the one hand $\rho$ is required, in general, to take values in the interval $(0,1]$. In practice, the larger the value of $\rho$, the faster the convergence of ADMM. More importantly, if the constants $c_{y}$ and $c_{u}$ are large, we essentially regularize the problem strongly (i.e. we force a large $\delta$, in the case of Algorithm 1). This means that tuning $\delta$ in Algorithm 2 will not allow us to accelerate the algorithm significantly (which is not the case for Algorithm 1).

4 The FDE-Constrained Optimization Model

In this section, we present the FDE-constrained optimization problem studied hereon and provide details as to the FDE discretization used. We then highlight some important properties of the resulting discretized matrices.

We define the Caputo derivative of a function ${\rm f}(t)$ defined on $t\in[t_{0},t_{1}]$, of real order $\alpha$ such that $n-1<\alpha<n$ with $n\in\mathbb{N}$, as follows:

[TABLE]

assuming convergence of the above [21, 51, 62]. We also define the left-sided and right-sided Riemann–Liouville derivatives of a function ${\rm f}(x)$ defined on $x\in[x_{0},x_{1}]$, of real order $\beta$ such that $n-1<\beta<n$ with $n\in\mathbb{N}$, as

[TABLE]

respectively. From this, we define the symmetric Riesz derivative as follows [62, 65]:

[TABLE]

We highlight that Caputo derivatives are frequently used for discretization of FDEs in time, given initial conditions, with Riemann–Liouville derivatives correspondingly considered for spatial derivatives, given boundary conditions. We consider the minimization problem:

[TABLE]

where the fractional differential equation and additional algebraic constraints are given on the space-time domain $\Omega\times(0,T)$, where $\Omega\subset\mathbb{R}^{2}$ has boundary $\partial\Omega$, and the spatial coordinates are given by $\bm{x}=[x_{1},x_{2}]^{T}$. We impose the initial condition ${\rm y}(\bm{x},0)=0$ at $t=0$, and the Dirichlet condition ${\rm y}=0$ on $\partial\Omega\times(0,T)$. We assume that the orders of differentiation satisfy $0<\alpha<1$, $1<\beta_{1}<2$, $1<\beta_{2}<2$.

The cost functional $J({\rm y},{\rm u})$ measures the misfit between the state variable ${\rm y}$ and a given desired state $\bar{\rm y}$ in some given norm, and also measures the ‘size’ of the control variable $\rm{u}$. In this paper we consider the cost functional $J({\rm y},{\rm u})$ corresponding to $L^{2}$-norms measuring both terms:

[TABLE]

Here $\gamma>0$ denotes a regularization parameter on the control variable. We note that other variants for $\mathrm{J}({\rm y},{\rm u})$ are possible, including measuring the state misfit and/or the control variable in other norms, as well as alternative weightings within the cost functionals. We also emphasize that it is perfectly reasonable to consider such problems involving FDEs in one or three spatial dimensions (or indeed higher dimensions), rather than in two dimensions as in (4.2), and the methodology in this paper could be readily tailored to such problems.

Upon discretization, we consider the non-shifted Grünwald–Letnikov formula [23, 62, 65, 66] to approximate the Caputo derivative in time:

[TABLE]

where $h_{t}$ is the step-size in time, and $g_{k}^{\alpha}=\frac{\Gamma(k-\alpha)}{\Gamma(-\alpha)\Gamma(k+1)}$ may be computed recursively via $g_{k}^{\alpha}=(1-\frac{\alpha+1}{k})g_{k-1}^{\alpha}$, $k=1,2,...,\nu$, with $g_{0}^{\alpha}=1$ and $\nu\in\mathbb{N}$. This leads to the Caputo derivative matrix for all grid points in the time variable:

[TABLE]

For the (left-sided) spatial derivative we use the $p$-shifted Grünwald–Letnikov formula [3, 50, 52, 62], with shift parameter $p=1$, to minimize the local truncation error:

[TABLE]

where $h_{x}$ is the step-size in space, leading to the matrix

[TABLE]

whereby using the formula (4.1) leads to the following Riemann–Liouville derivative matrix for the symmetrized Riesz derivative:

[TABLE]

Using all the previous definitions, we can write the discretized version of the FDE constraint within (4.2) as

[TABLE]

where $y_{\bm{n}}$, $u_{\bm{n}}$, $g_{\bm{n}}$ represent the discretized variants of ${\rm y}$, ${\rm u}$, ${\rm g}$, $\bm{n}=[n_{x},n_{y},n_{t}]$ is a $3$-index containing the grid sizes along each dimension, and

[TABLE]

For simplicity of exposition, in the rest of the paper we assume that $h_{x_{1}}=h_{x_{2}}=h_{x}$, where $h_{x_{i}}$ is the discretization step in the respective spatial direction, noting that the method readily generalizes to problems where this is not the case.

By using the trapezoidal rule we approximate the two terms in the objective functional (4.3), by

[TABLE]

which is applied to vectors arising from every time-step, apart from the initial time $t=0$. Notice that matrix $J_{\bm{n}}$ is diagonal with only two different values on the diagonal, and hence can be almost exactly approximated by a scaled identity.

We should mention that we assume constant diffusion coefficients in the FDE constraints. In turn, this yields that the discretized constraint matrix has a multilevel Toeplitz structure. As we discuss in the following section, such matrices can be approximated by circulant preconditioners, which in turn allow us to use the preconditioner in (3.10) for accelerating the solution of the resulting ADMM sub-problems. In the presence of non-constant diffusion coefficients, the discretized constraint matrices would belong to the class of multilevel GLT sequences. In this case, circulant preconditioners would no longer be effective and we would have to approximate such matrices using diagonal times multilevel banded Toeplitz matrices (see for example [24, 53]). In light of the discussion in Section 3, we can observe that one could extend the results presented in this paper to the non-constant diffusion coefficient case, by making use of the preconditioner in (3.13) (upon noting that the discretization of the functional in (4.3) yields a diagonal matrix). For brevity of presentation this is left to a future study.

In the following proposition, we summarize some well-known properties of the fractional binomial coefficients that arise above when constructing the matrices $\mathscr{C}_{\alpha}$ and $\mathscr{L}_{\beta}$ (see for example [39, page 397], or [51, 72]):

Proposition 4.1**.**

Let $0<\alpha<1$ and $1<\beta<2$, with $g_{k}^{\alpha},\ g_{k}^{\beta}$ as in (4.4), (4.6). Then, we have that:

[TABLE]

[TABLE]

5 Toeplitz Matrices and Circulant Preconditioners

In this section, we propose a multilevel circulant preconditioner, suitable for approximating multilevel Toeplitz matrices, and then examine the quality of such a preconditioner for the problem at hand, showing that the preconditioner is in fact a.c.s. for a scaled version of the coefficient matrix in (4.9).

Toeplitz and multilevel Toeplitz matrices appear often when (numerically) solving partial, integral, or fractional differential equations, problems in time series analysis, as well as in signal processing (see for example [1, 46, 63, 74], and the references therein). An active area of research is that of solving a huge-scale systems of linear equations, $Ax=b$, where the matrix $A$ has some specific structure, such as Toeplitz, multilevel Toeplitz, or it can be written as a combination of Toeplitz and other structured matrices. There are two major approaches for solving such systems. One alternative is to solve them directly by exploiting the matrix structure (see for example [6, 16, 47, 73]). A more popular approach is to employ some iterative method to solve the system, assisted by an appropriately designed preconditioner, to ensure that the iterative method achieves fast convergence (as in [8, 9, 10, 11, 12, 13, 14, 15, 42, 44, 45, 46, 56, 71]). An equally rich literature exists for preconditioning Toeplitz-like linear systems arising specifically from the discretization of fractional diffusion equations (see [17, 24, 25, 29, 32, 45, 46, 48, 53], among others).

In this paper, we follow the simplest possible approach: that of approximating multilevel Toeplitz matrices using multilevel circulant preconditioners. To do so, we first have to derive a unilevel circulant approximation of an arbitrary unilevel Toeplitz matrix. Given a unilevel Toeplitz matrix $T_{n}\in\mathbb{R}^{n\times n}$, we employ the circulant approximation proposed for the first time in [14] (also called the T. Chan preconditioner for $T_{n}$). More specifically, we define the optimal circulant approximation of $T_{n}$, as the solution of the following optimization problem:

[TABLE]

where $\mathcal{C}_{n}$ is the set of all $n\times n$ circulant matrices, and $\|\cdot\|_{F}$ the Frobenius norm. It turns out that (5.1) admits the following closed form solution:

[TABLE]

Then, we can write $C_{1}(T_{n})=F_{n}^{*}\Lambda_{n}F_{n}$, where $F_{n}$ is the discrete Fourier transform of size $n$ and $\Lambda_{n}$ is a diagonal matrix containing the eigenvalues of $C_{1}(T_{n})$, which can be computed as $\Lambda_{n}=\textnormal{diag}(F_{n}c_{1})$, where $c_{1}$ is the first column of $C_{1}(T_{n})$. Other unilevel circulant approximations are possible, such as those proposed in [12, 13, 71], however, the T. Chan preconditioner seems (empirically) to behave better for the problem under consideration.

We now focus on the discretized FDE given in (4.8). By multiplying this equation on both sides by $\psi=\min\{h_{t}^{\alpha},h_{x}^{\beta_{1}},h_{x}^{\beta_{2}}\}$, we have:

[TABLE]

where $y_{\bm{n}}$, $u_{\bm{n}}$, $g_{\bm{n}}$ represent the discretized variants of ${\rm y}$, ${\rm u}$, ${\rm g}$, $B_{\bm{n}}=\psi D_{\bm{n}}$, with $D_{\bm{n}}$ defined as in (4.9), $h_{t},\ h_{x}$ the time and spatial mesh-sizes, and $\bm{n}=[n_{x_{1}},n_{x_{2}},n_{t}]$. We observe that the matrix $D_{\bm{n}}$ (and hence $B_{\bm{n}}$) enjoys a 3-level Toeplitz structure. In particular, each block of $D_{\bm{n}}$ ($B_{\bm{n}}$) enjoys a quadrantally symmetric block Toeplitz structure (such matrices are analyzed for example in [10]). Given the matrix $B_{\bm{n}}$, we can define its T. Chan-based $3$-level circulant preconditioner as:

[TABLE]

where $\Lambda_{\bm{n}}$ is the diagonal eigenvalue matrix of the preconditioner, computed as:

[TABLE]

with $\Lambda_{\alpha}$, $\Lambda_{\beta_{1}}$, $\Lambda_{\beta_{2}}$ being the diagonal matrices containing the eigenvalues of the T. Chan approximations of the matrices $\mathscr{C}^{\alpha}_{n_{t}}$, $\mathscr{L}^{\beta_{1}}_{n_{x_{1}}}$, and $\mathscr{L}^{\beta_{2}}_{n_{x_{2}}}$, respectively.

The preconditioner in (5.2) can be computed efficiently in $O(N(\bm{n})\log N(\bm{n}))$ operations, using the fast Fourier transform. The storage requirements are $O(N(\bm{n}))$ since we only need to store the eigenvalue matrix, that is $\Lambda_{\bm{n}}$. Clearly, the preconditioner in (5.2) can be defined similarly for FDEs of arbitrary dimension, say $d$. Given a $d$-index $\bm{n}$, containing the level sizes of an arbitrary $d$-level Toeplitz $T_{\bm{n}}$ or circulant matrix $C_{\bm{n}}$, we summarize the computational and storage costs of various recursive linear algebra operations in Table 1.

Using the definition of the matrices used to construct matrix $B_{\bm{n}}$ (see (4.5) and (4.7)), we are now able to derive the generating function of this $3$-level Toeplitz matrix. To that end, let us define the following scalars:

[TABLE]

which are obviously bounded above by $1$, from the definition of $\psi$. Of course in order for these to be theoretically meaningful, we have to assume that $h_{t}^{\alpha}\propto h_{x}^{\beta_{1}}\propto h_{x}^{\beta_{2}}$.

Lemma 5.1**.**

Let $\bm{n}=[n_{x_{1}},n_{x_{2}},n_{t}]$ be a $3$-index and define the matrix $D_{\bm{n}}$ as in (4.9). Then, the symbol generating the matrix-sequence $\{B_{\bm{n}}\}_{n}=\{\psi D_{\bm{n}}\}_{n}$, can be expressed as:

[TABLE]

where $\bm{\theta}=[\theta_{1},\theta_{2},\theta_{3}]$, and $g_{k}^{c}$ is the fractional binomial coefficient, for some $c\in(0,1)\cup(1,2)$ and an arbitrary $k\geq 0$. Thus, we can write $B_{\bm{n}}=T_{\bm{n}}(\phi_{\beta_{1},\beta_{2},\alpha})$.

Proof.

We omit the proof, which follows easily from the definition of the matrices within $B_{\bm{n}}$, that is using the definition of $\mathscr{C}^{\alpha}_{n}$ in (4.5) and $\mathscr{L}^{\beta}_{n}$ in (4.7). The reader is referred to [24, 46, 53], among others, for derivations of similar results. The alternative representation of matrix $B_{\bm{n}}$ follows directly from Theorem 2.8. ∎

To analyze the effectiveness of the proposed $3$-level circulant preconditioner for $B_{\bm{n}}$, we prove that the trigonometric polynomial generating function of matrix $B_{\bm{n}}$ belongs to the Wiener class (that is, it has absolutely summable coefficients).

Lemma 5.2**.**

Assume that $\beta_{1}$ and $\beta_{2}$ are bounded away from $1$. Then, the generating function $\phi_{\beta_{1},\beta_{2},\alpha}(\bm{\theta})$ defined in (5.4) belongs to the Wiener class, that is:

[TABLE]

Proof.

For brevity of presentation, we provide an outline of the proof. Firstly, one has to transform (5.4) to the form $\phi_{\beta_{1},\beta_{2},\alpha}(\bm{\theta})=\sum_{\bm{k}\in\mathbb{Z}^{3}}\phi_{\bm{k}}e^{i\langle\bm{k},\bm{\theta}\rangle}$, by matching the coefficients of the associated trigonometric polynomials. By taking the absolute values of the matched coefficients, applying the triangle inequality, and using the properties of the fractional binomial coefficients, summarized in Proposition 4.1, we obtain:

[TABLE]

The latter completes the proof. ∎

Using the results presented in [10, 44, 45], we can derive the following Theorem, which in fact shows that the $3$-level circulant approximation of matrix $B_{\bm{n}}$ defined in (5.2) is an a.c.s. for it.

Theorem 5.3**.**

Let $B_{\bm{n}}=\psi D_{\bm{n}}$ where $\bm{n}=[n_{x_{1}},n_{x_{2}},n_{t}]$, and $C_{3}(B_{\bm{n}})$ its circulant approximation defined in (5.2). For every $\epsilon(m)>0$, such that $\epsilon(m)\rightarrow 0$ as $m\rightarrow\infty$, there exist constants $N_{x_{1}},N_{x_{2}},N_{t}$, such that for all $n_{x_{1}}>N_{x_{1}}$, $n_{x_{2}}>N_{x_{2}}$, $n_{t}>N_{t}$:

[TABLE]

where

[TABLE]

Proof.

The proof is omitted since it follows exactly the developments in [45, Theorem 3.2 and Theorem 4.1], with the only difference being that the Strang unilevel circulant approximation is used there (for example see [12]) instead of the T. Chan approximation. Notice that the authors in [45] assume invertibility of $C_{3}(B_{\bm{n}})$, using which they prove a weak clustering result. Hence, to prove the result stated here, one needs to follow only part of the proof outlined in [45, Theorem 3.2]. ∎

Remark 5.1**.**

Following [45, Remark 4.1], assuming that $d=O(1)$, we can recursively extend the result of Theorem 5.3 to the $d$-level case, using induction. In other words, the developments discussed in this paper can be extended trivially to higher dimensional FDEs. As expected, the circulant approximation becomes weaker as the dimension of the associated FDE is increased. In particular, the result in [67] shows that in the general case, any multilevel circulant preconditioner for multilevel Toeplitz matrices is not a superlinear preconditioner. Superlinear preconditioners are important, in that they allow preconditioned Conjugate Gradient-like methods to converge in a constant number of iterations, independently of the size of the problem. In general, one could not hope of achieving a strong clustering when preconditioning multilevel Toeplitz matrices using multilevel circulant preconditioners. In light of that, it comes as no surprise that a preconditioner like the one in (5.2) does not asymptotically capture all of the eigenvalues of the approximated multilevel Toeplitz matrix.

Remark 5.2**.**

Let us now notice that a scaled identity approximation for the discretized objective Hessian matrix in (4.10) yields (trivially) a GLT sequence. Similarly, the approximation $C_{3}(B_{\bm{n}})$ in (5.2) for the matrix $B_{\bm{n}}=\psi D_{\bm{n}}$, where $D_{\bm{n}}$ is defined in (4.9), is also a GLT sequence (since it can be considered as a multilevel Toeplitz matrix). In view of the previous, as well as Theorem 2.16 (condition (6.)), we can see that the proposed approximations for the matrices associated to the discretized version of (4.2) satisfy the conditions of Proposition 3.2. Hence, we are able to invoke Theorem 3.3 for the preconditioner in (3.10), which is constructed by using the aforementioned multilevel circulant approximations. Thus, we are able to show that the resulting preconditioned ADMM system matrix, corresponding to the normal equations in (3.9), is weakly clustered at $1$. Furthermore, by the same Theorem, we expect convergence of PCG in a number of iterations independent of the grid-size.

6 Implementation Details and Numerical Results

In this section we discuss specific implementation details and present the numerical results obtained by running the implementation of the proposed method over a variety of settings of the FDE optimization problem.

6.1 Test problem and implementation details

We assess the performance of the proposed method on the following test problem. We attempt to numerically solve problem (4.2). The state and the control are defined on the domain $\Omega\times(0,T)=(0,1)^{2}\times(0,1)$. For some $n\in\mathbb{N}$, the discretized grid contains $n\times n\times n$ uniform points, in space and time (i.e. we make use of the $3$-index $\bm{n}=[n,n,n]$), which yields:

[TABLE]

It is worth mentioning that the choice of the number of discretization points in time should depend on the value of the fractional derivative orders. In particular, in the theory we had to assume that $h_{t}^{\alpha}\propto h_{x}^{\beta_{1}}\propto h_{x}^{\beta_{2}}$. Of course, this could be difficult to satisfy for certain values of $\alpha,\ \beta_{1}$, and $\beta_{2}$. In terms of discretization error, $n_{t}=n$ suffices, as we employ first-order numerical schemes for the space and time fractional derivatives. In what follows, we choose to use $n_{t}=n$ throughout all the experiments, noting that for very large values of $n$, this should be adjusted to take into consideration the values of the fractional derivative orders. Such an increase in the number of discretization points in time, could potentially be tackled by the use of higher-order numerical methods for the space fractional derivatives (see [22] and the references therein for higher-order approximations for the Riemann–Liouville and Riesz fractional derivatives). This is omitted for a future study.

As a desired state function, we set $\bar{{\rm y}}(x_{1},x_{2},t)=10\cos(10x_{1})\sin(x_{1}x_{2})(1-e^{-5t}),$ as in [23, Section 5.1], with homogeneous boundary and initial conditions. Throughout this section, we employ the convention that $n_{x_{1}}=n_{x_{2}}=n_{t}$, and we only present the overall size of the discretized state vector, that is $N(\bm{n})=n_{x_{1}}\cdot n_{x_{2}}\cdot n_{t}=n_{x_{1}}^{3}$. As an indicator of convergence of the numerical method, we apply the trapezoidal rule to roughly approximate the discrepancy between the solution for the state and the desired state on the discrete level, i.e.:

[TABLE]

We should note that the previous measure approximates the misfit between the state and the desired state of the continuous problem, and hence it is not expected to converge to zero. Due to the Dirichlet boundary conditions, there is a mismatch between $y$ and $\bar{y}$ on the boundary. Hence a refinement in the grid size is expected to result in slight increase in the approximate discrepancy measure.

We implement a standard 2-Block ADMM for solving problems of the form of (3.4). The implementation follows exactly the developments in Section 3. We solve system (3.9) using the MATLAB function pcg. We note that while various potential acceleration strategies for ADMMs have been studied in the literature (see for example [5, 37]), the focus of the paper is to illustrate the viability of the proposed approach, and hence the simplest possible ADMM scheme is adopted. The step-size of ADMM is chosen to be close to the maximum allowed one in all computations, that is $\rho=1.618$. The termination criteria of the ADMM are summarized as follows:

[TABLE]

In order to avoid unnecessary computations, we do not require a specific tolerance for the dual infeasibility. Instead, we report the dual infeasibility at the accepted optimal point. The Krylov solver tolerance is set dynamically, based on the accuracy attained at the respective ADMM iteration. In particular, the required tolerance for the Krylov solver is set to:

[TABLE]

at every iteration $j$. Hence, we present the average number of inner iterations in the results to follow. Furthermore, we employ the convention that the discretized restricting functions are of the form $y_{b_{\bm{n}}}=-y_{a_{\bm{n}}}=c\cdot e_{\bm{n}}$ (or $u_{b_{\bm{n}}}=-u_{a_{\bm{n}}}=c\cdot e_{\bm{n}}$), where $e_{\bm{n}}$ is the $N(\bm{n})$-dimensional vector of ones and $c>0$. Thus, we present only the value of the entries of $y_{a_{\bm{n}}}$ ($u_{a_{\bm{n}}}$, respectively).

As we discussed earlier, the FDE constraints were scaled by the constant $\psi$, since this was required from the theory (see Theorem 5.3). By doing this, we ensure that the elements of the matrix $B_{\bm{n}}$ are of order 1 (assuming that $h_{t}^{\alpha}\propto h_{x_{i}}^{\beta_{i}},$ for $i=1,2$). As a result, the discretized control in the FDE constraints is multiplied by $\psi$. In ADMM such a scaling translates to a scaled step of the dual variables corresponding to the FDE constraints. In order to improve the balance of the algorithm, we multiply by $\psi$ the constraints linking $u_{\bm{n}}$ with its copy variables $z_{u_{\bm{n}}}$, thus scaling all the dual multipliers corresponding to these constraints.

The penalty parameter of ADMM, $\delta$, is chosen from a pool of five values which deliver reasonably good behavior of the method. More specifically, for the experiments to follow we choose $\delta\in\{0.1,0.4,2,10,100\}$. We note here that one could tune this parameter for each problem instance and obtain significantly better results. However, as this is not practical, we restrict ourselves to a small set of possible values.

The experiments were conducted on a PC with a 2.2 GHz Intel (hexa-) core i7 processor, run under the Windows 10 operating system. The code is written in MATLAB R2019a.

6.2 Numerical Results

We distinguish three types of problems:

• •

Problems with box constraints on the state $\rm{y}$,

• •

problems with box constraints on the control $\rm{u}$, and

• •

problems with box constraints on both variables.

As expected and verified in practice, the third type of problem is the most difficult one. Hence, we will focus our attention on problems with box constraints on both variables, while presenting a few experiments on problems of the other two types.

Box constraints on the state $\bm{{\rm y}}$

Let us briefly focus on the case where the state variable is required to stay in a box, while the control is free, that is ${\rm y_{a}}\leq{\rm y}\leq{\rm y_{b}},\ -\infty\leq{\rm u}\leq\infty$. Using similar arguments as in [23, 29], we can see that an optimal solution in this case is guaranteed to exist. We run the method for different inequality bounds on the state ${\rm y}$. The results are summarized in Table 6.2. All fixed parameters are provided at the title of the respective Table.