Low-rank Linear Fluid-structure Interaction Discretizations

TL;DR

This paper introduces low-rank iterative methods for efficient parameter-dependent fluid-structure interaction discretizations, significantly reducing computational time while maintaining high accuracy in simulations.

Contribution

It develops low-rank GMRES and Chebyshev iteration variants with error bounds, enabling fast, accurate solutions for parameter-dependent FSI models.

Findings

Truncated methods achieve residual norms below 10^{-8}

Methods are twenty times faster than standard approaches

Error bounds for truncation are rigorously derived

Abstract

Fluid-structure interaction models involve parameters that describe the solid and the fluid behavior. In simulations, there often is a need to vary these parameters to examine the behavior of a fluid-structure interaction model for different solids and different fluids. For instance, a shipping company wants to know how the material, a ship's hull is made of, interacts with fluids at different Reynolds and Strouhal numbers before the building process takes place. Also, the behavior of such models for solids with different properties is considered before the prototype phase. A parameter-dependent linear fluid-structure interaction discretization provides approximations for a bundle of different parameters at one step. Such a discretization with respect to different material parameters leads to a big block-diagonal system matrix that is equivalent to a matrix equation as discussed in [KressnerTobler 2011]. The unknown is then a matrix which can be approximated using a low-rank approach that represents the iterate by a tensor. This paper discusses a low-rank GMRES variant and a truncated variant of the Chebyshev iteration. Bounds for the error resulting from the truncation operations are derived. Numerical experiments show that such truncated methods applied to parameter-dependent discretizations provide approximations with relative residual norms smaller than $10^{-8}$ within a twentieth of the time used by individual standard approaches.

Tables

Table 1. Table 1: GMRESTR Compared With Standard GMRES

Method	Approx.	Computation Times
	Storage	(in Minutes)
		Precon.	Comp.	Total
GMRESTR	$O\left[\left(M+m+R\right)R\right]$
($R=200$)	$\approx 306.12$MB	1.24	179.88	181.12
Standard GMRES	$O\left(Mm\right)$
($8000$ times)	$\approx 11.47$GB	6.63	5426.23	5432.86

Table 2. Table 2: ChebyshevT Compared With Standard Chebyshev

Method	Approx.	Computation Times
	Storage	(in Minutes)
		Est.	Precon.	Comp.	Total
ChebyshevT	$O\left[\left(M+m+R\right)R\right]$
($R=200$)	$\approx 306.12$MB	0.013	1.24	177.99	179.243
Standard Chebyshev	$O\left(Mm\right)$
($8000$ times)	$\approx 11.47$GB	$0.013$	$6.63$	$5490.85$	$5497.493$

Table 3. Table 3: Computation Time Comparison of the Truncated Approaches

Method	Computation Times
(R=200)	(in Minutes)
	Est.	Precon.	Comp.	Total
ChebyshevT	$0.013$	$1.24$	$177.99$	$179.24$
GMREST	-	$1.24$	$179.88$	$181.12$
Truncated Bi-CGstab	-	$1.24$	$302.94$	$304.18$

Equations

$$\displaystyle\big{(}A_{0}+\mu_{s}^{i}A_{1}+\lambda_{s}^{i}A_{2}+\rho_{f}^{i}A_{3}\big{)}x_{i}=b_{D}\quad\text{for}\quad i\in\{1,...,m\}\text{,}$$

$$\displaystyle\mathcal{A}\left(\begin{array}[]{c}x_{1}\\ \vdots\\ x_{m}\end{array}\right)=\left(\begin{array}[]{c}b_{D}\\ \vdots\\ b_{D}\end{array}\right)\text{,}$$

$$A_{0}X+A_{1}X{D}_1+A_{2}X{D}_2+A_{3}X{D}_3=B$$

$$\displaystyle X=\left[\begin{array}[]{c|c|c}x_{1}&\cdots&x_{m}\end{array}\right]\in\mathbb{R}^{M\times m}$$

$$\hat{X}=\sum _{j=1}^R{u}_j\otimes v_j^T\text{, }{u}_j\in {\mathbb{R}}^M\text{and}{v}_j\in {\mathbb{R}}^M\forall j\in \{1,...,R\}$$

$$⟨\nabla \cdot v,\xi {⟩}_F$$

$${\mu }_s⟨\nabla u+\nabla u^T,\nabla \phi {⟩}_S+{\lambda }_s⟨\mathrm{tr}(\nabla u)I,\nabla \phi {⟩}_S+{\nu }_f{\rho }_f⟨\nabla v+\nabla v^T,\nabla \phi {⟩}_F$$

$$-⟨p,\nabla \cdot \phi {⟩}_F$$

$$⟨\nabla u,\nabla \psi {⟩}_F$$

$${\nu }_s^p=\frac{{\lambda }_s}{2({\lambda }_s+{\mu }_s)}\text{.}$$

$$\{{\mu }_s^{i_1}{\}}_{i_1\in \{1,...,m_1\}}\subset {\mathbb{R}}^{+}$$

$$\{{\lambda }_s^{i_2}{\}}_{i_2\in \{1,...,m_2\}}\subset {\mathbb{R}}^{+}$$

$$\{{\rho }_f^{i_3}{\}}_{i_3\in \{1,...,m_3\}}\subset {\mathbb{R}}^{+}$$

$$A_1\text{ discretizes }$$

$$A_2\text{ discretizes }$$

$$A_3\text{ discretizes }$$

$$\displaystyle\underbrace{\big{(}A_{0}+(\mu_{s}^{i_{1}}-\mu_{s})A_{1}+(\lambda_{s}^{i_{2}}-\lambda_{s})A_{2}+\nu_{f}(\rho_{f}^{i_{3}}-\rho_{f})A_{3}\big{)}}_{=:A(\mu_{s}^{i_{1}},\lambda_{s}^{i_{2}},\rho_{f}^{i_{3}})}x_{i_{1}i_{2}i_{3}}=b_{D}$$

$$\text{for}(i_1,i_2,i_3)\in \{1,...,m_1\}\times \{1,...,m_2\}\times \{1,...,m_3\}$$

$$\displaystyle\big{(}A_{0}+\mu_{s}^{i_{1}}A_{1}+\lambda_{s}^{i_{2}}A_{2}+$$

$$(i_1,i_2,i_3)\in \{1,...,m_1\}\times \{1,...,m_2\}\times \{1,...,m_3\}\text{.}$$

$$\displaystyle\overbrace{\operatorname*{diag}\limits_{\makebox[0.0pt]{$\begin{subarray}{c}i_{1}\in\{1,...,m_{1}\}\\ i_{2}\in\{1,...,m_{2}\}\\ i_{3}\in\{1,...,m_{3}\}\end{subarray}$}}\big{(}A(\mu_{s}^{i_{1}},\lambda_{s}^{i_{2}},\rho_{f}^{i_{3}})\big{)}}^{=:\mathcal{A}}\left(\begin{array}[]{c}x_{1}\\ \vdots\\ x_{m_{1}m_{2}m_{3}}\end{array}\right)=\left(\begin{array}[]{c}b_{D}\\ \vdots\\ b_{D}\end{array}\right)\text{,}$$

$$D_1$$

$$\text{and}{D}_3$$

$$\underset{=:F(X)}{\underbrace{A_0+A_{1}X{D}_1+A_{2}X{D}_2+{\nu }_f{A}_{3}X{D}_3}}=B:=[b_D|\cdots |b_D]$$

$$\mathcal{A}=I_m\otimes A_0+D_1\otimes A_1+D_2\otimes A_2+{\nu }_f{D}_3\otimes A_3\text{.}$$

$${\mathcal{P}}_{A_0}:=I_m\otimes A_0$$

$$\displaystyle\mathcal{P}_{T}:=I_{m}\otimes\underbrace{\big{(}A_{0}+\bar{\mu}_{s}A_{1}+\bar{\lambda}_{s}A_{2}+\nu_{f}\bar{\rho}_{f}A_{3}\big{)}}_{=:P_{T}}$$

$${\bar{\mu }}_s$$

$${\bar{\lambda }}_s$$

$${\bar{\rho }}_f$$

$$\displaystyle\mathcal{A}\left(\begin{array}[]{c}x_{1}\\ \vdots\\ x_{m_{1}m_{2}m_{3}}\end{array}\right)$$

$$\displaystyle\left(\begin{array}[]{c}x_{1}\\ \vdots\\ x_{m_{1}m_{2}m_{3}}\end{array}\right)$$

$$\hat{X}=\sum _{j=1}^R=(u_j\otimes v_j^T)\text{with}{u}_j\in {\mathbb{R}}^M\text{, }{v}_j\in {\mathbb{R}}^m\forall j\in \{1,...,R\},$$

$$T_{(r_1,r_2)}(V):=\{v\in V_1\otimes V_2:V_1\subset {\mathbb{R}}^M\text{, }\dim (V_1)=r_1\text{, }{V}_2\subset {\mathbb{R}}^m\text{, }\dim (V_2)=r_2\}\text{.}$$

$$\mathcal{T}:{\mathbb{R}}^M\otimes {\mathbb{R}}^m\to T_R$$

$${\mathbb{R}}^M\otimes {\mathbb{R}}^m\cong {\mathbb{R}}^{M\times m}\text{,}$$

Full text

{AMS}

65M22, 65F10, 65M15, 15A69

1 Introduction

A parameter-dependent linear fluid-structure interaction problem as described in Section 2 discretized using bilinear finite elements with a total number of $M\in\mathbb{N}$ degrees of freedom (see Section 3 for details) and $m\in\mathbb{N}$ parameter combinations leads to equations of the form

[TABLE]

where the discretization matrices $A_{0},A_{1},A_{2},A_{3}\in\mathbb{R}^{M\times M}$ and the right hand side $b_{D}\in\mathbb{R}^{M}$ depends on the Dirichlet data and the $i$th finite element solution $x_{i}\in\mathbb{R}^{M}$. The samples of interest are given by the shear moduli $\mu_{s}^{i}\in\mathbb{R}$, the first Lamé parameters $\lambda_{s}^{i}\in\mathbb{R}$ and the fluid densities $\rho_{f}^{i}$ for $i\in\{1,...,m\}$.

Equation Eq. 1 can directly be written as the linear system

[TABLE]

where $\mathcal{A}\in\mathbb{R}^{Mm\times Mm}$ is a block diagonal matrix. Following [9], equation Eq. 8 can then be translated into the matrix equation

[TABLE]

with B:=\left[\begin{array}[]{c|c|c}b_{D}&\cdots&b_{D}\end{array}\right] and diagonal matrices $D_{1},D_{2},D_{3}\in\mathbb{R}^{M\times M}$, where the $i$th diagonal entry of these diagonal matrices is given by $\mu_{s}^{i},\lambda_{s}^{i}$ and $\rho_{f}^{i}$, respectively. In Eq. 9, the unknown

[TABLE]

is a matrix. Now an iterative method to solve linear systems can be modified such that it uses an iterate that is a matrix. It is applied to the big system Eq. 8 but computation is kept in the matrix notation Eq. 9 by representing the iterate as a matrix instead of a vector. The methods used in this paper fix a rank $R\in\mathbb{N}$, $R\ll M,m$ and represent this iterate as a tensor. The goal is to find a low-rank approximation $\hat{X}$ of rank $R$

[TABLE]

that approximates the full matrix $X$ from Eq. 9 and therefore provides (parameter-dependent) finite element approximations for all equations in Eq. 1.

Fluid-structure interaction problems yield non-symmetric system matrices. Hence, the system matrix $\mathcal{A}$ in Eq. 8 is not symmetric. The methods examined in this paper are based on the GMRES method as introduced in a truncated variant in [1] and the Chebyshev method from [11]. These methods will then be compared to a truncated method based on the Bi-CGstab method from [18] similar to Algorithm 3 of [9] and Algorithm 2 of [3].

In Section 2 and Section 3, we derive the matrix equations that appear when dealing with parameter-dependent fluid-structure interaction discretizations. The low-rank framework and related methods are introduced in Section 4 for stationary problems and generalized to non-stationary problems in Section 5. In Section 6, theoretical error bounds for the GMREST and the ChebyshevT method are derived and numerically evaluated in Section 7. The convergences of the truncated approaches presented are compared in numerical experiments in Section 8.

2 The Stationary Linear Fluid-structure Interaction Problem

Let $d\in\{2,3\}$, $\Omega\subset\mathbb{R}^{d}$, $F$, $S\subset\Omega$ such that $\bar{F}\cup\bar{S}=\bar{\Omega}$ and $F\cap S=\emptyset$, where $F$ represents the fluid and $S$ the solid part. Let $\Gamma_{\text{int}}=\partial F\cap\partial S$ and $\Gamma_{f}^{\text{out}}\subset\partial F\setminus\partial S$ denote the boundary part where Neumann outflow conditions hold. $\Gamma_{f}^{D}=\partial F\setminus(\Gamma_{f}^{\text{out}}\cup\Gamma_{\text{int}})$ denotes the boundary part where Dirichlet conditions hold. Consider the Stokes fluid equations from Section 2.4.4 of [13] as a model for the fluid part and the Navier-Lamé equations discussed as Problem 2.23 of [13] as a model for the solid part. Both equations are assumed to have a vanishing right hand side. If these two equations are coupled with the kinematic and the dynamic coupling conditions discussed in Section 3.1 of [13], the weak formulation of the stationary coupled linear fluid-structure interaction problem reads

[TABLE]

with the trial functions $v\in v_{\text{in}}+H_{0}^{1}(\Omega,\Gamma_{f}^{D}\cup\Gamma_{\text{int}})^{d}$ (velocity), where $v_{\text{in}}\in H^{1}(\Omega)^{d}$ is an extension of the Dirichlet data on $\Gamma_{f}^{D}$, $u\in H^{1}_{0}(\Omega)^{d}$ (deformation) and $p\in L^{2}(F)$ (pressure) and the test functions $\xi\in L^{2}(F)$ (divergence equation), $\varphi\in H^{1}_{0}(\Omega,\partial\Omega\setminus\Gamma_{f}^{\text{out}})^{d}$ (momentum equation) and $\psi\in H^{1}_{0}(F)^{d}$ (deformation equation). $\langle\cdot,\cdot\rangle_{S}$ and $\langle\cdot,\cdot\rangle_{F}$ denote the $\mathcal{L}^{2}$ scalar product on $S$ and $F$, respectively. $\nu_{f}\in\mathbb{R}$ denotes the kinematic fluid viscosity and $\rho_{f}\in\mathbb{R}$ the fluid density. The shear modulus $\mu_{s}\in\mathbb{R}$ and the first Lamé parameter $\lambda_{s}\in\mathbb{R}$ determine the Poisson ratio of the solid.

Definition 2.1** (The Poisson Ratio - Definition 2.18 of [13]).**

The Poisson ratio of a solid is given by the number

[TABLE]

*It describes the compressibility of a solid. *

3 Parameter-dependent Discretization

Assume the behavior of a linear fluid-structure interaction model for $m_{1}\in\mathbb{N}$ shear moduli, $m_{2}\in\mathbb{N}$ first Lamé parameters and $m_{3}\in\mathbb{N}$ fluid densities is of interest. The kinematic fluid viscosity $\nu_{f}\in\mathbb{R}$ is assumed to be fixed. Let the samples of interest be given by the following sets.

[TABLE]

In a bilinear finite element discretization of Eq. 10 with a mesh grid size of $N\in\mathbb{N}$, every mesh grid point corresponds to a pressure, a velocity and a deformation variable. In two dimensions, the velocity and deformation are two dimensional vectors, in three dimensions they correspond to a three dimensional vector each. The total number of degrees of freedom is therefore $M=5N$ in two dimensions and $M=7N$ in three dimensions.

Let $\Omega_{h}$ be a matching mesh of the domain $\Omega$ as defined in Definition 5.9 of [13] with $N$ mesh grid points. $A_{0}\in\mathbb{R}^{M\times M}$ is a discrete differential operator restricted to the finite element space with dimension $M$. It discretizes all operators involved in Eq. 10 with a fixed shear modulus $\mu_{s}\in\mathbb{R}$, a fixed first Lamé parameter $\lambda_{s}\in\mathbb{R}$ and a fixed fluid density $\rho_{f}\in\mathbb{R}$. In this paper, $Q_{1}$ finite elements as discussed in Section 4.2.1 of [13] are used and we will denote the discrete differential operators by discretization matrices. Moreover, let $A_{1},A_{2},A_{3}\in\mathbb{R}^{M\times M}$ be the discretization matrices of the following operators:

[TABLE]

The parameter-dependent equation

[TABLE]

is the finite element discretization of Eq. 10 related to a shear modulus $\mu_{s}^{i_{1}}$, a first Lamé parameter $\lambda_{s}^{i_{2}}$ and a fluid density $\rho_{f}^{i_{3}}$. The finite element solution is $x_{i_{1}i_{2}i_{3}}\in\mathbb{R}^{M}$ and the right hand side $b_{D}\in\mathbb{R}^{M}$ depends on the Dirichlet data.

Remark 3.1**.**

If the fixed parameters vanish, namely $\mu_{s}=\lambda_{s}=\rho_{f}=0$, Eq. 11 translates to

[TABLE]

*At first sight, this presentation seems to be more convenient. But choosing, for instance, the parameters $\mu_{s}=\mu_{s}^{1}$, $\lambda_{s}=\lambda_{s}^{1}$ and $\rho_{f}=\rho_{f}^{1}$ minimizes the number of nonzero entries in the diagonal matrices $D_{1}$, $D_{2}$, $D_{3}\in\mathbb{R}^{m\times m}$ that will be introduced in Eq. 19. From a numerical point of view, this is an advantage. Furthermore, vanishing fixed parameters would lead to a singular matrix $A_{0}$. This can become a problem if the preconditioner $\mathcal{P}_{A_{0}}$ from Subsection 3.2 is used. *

Combining all sample combinations in Eq. 11 leads to a total of $m=m_{1}m_{2}m_{3}$ equations. Written as a linear system, these equations translate to

[TABLE]

where $\operatorname{diag}(\cdot)$ denotes the operator introduced in Section 1.2.6 of [5] extended to block diagonalization. Even though $\mathcal{A}\in\mathbb{R}^{Mm\times Mm}$ is of block diagonal structure, solving the blocks on the diagonal $m$ times (potentially in parallel) is often not feasible. If $100$ samples per parameter are considered, one would have to face $100^{3}=10^{6}$ blocks already in such a direct approach. This would lead to huge storage requirements for the solution vectors.

3.1 The Matrix Equation

The diagonal entries of $D_{1}$, $D_{2}$ and $D_{3}\in\mathbb{R}^{m\times m}$ are $\mu_{s}^{i_{1}}-\mu_{s}$, $\lambda_{s}^{i_{2}}-\lambda_{s}$ and $\rho_{f}^{i_{3}}-\rho_{f}$, respectively. The order of the diagonal entries has to be chosen such that every parameter combination occurs only once. If $I_{m_{1}}\in\mathbb{R}^{m_{1}\times m_{1}}$ denotes the $m_{1}\times m_{1}$ identity matrix, a possible sample order would lead to matrices

[TABLE]

As discussed in [9], equation Eq. 8 can then be written as the matrix equation

[TABLE]

where the unknown is the matrix $X=[x_{1}|\cdots|x_{m_{1}m_{2}m_{3}}]\in\mathbb{R}^{M\times m}$ whose $i$th column corresponds to the finite element approximation of Eq. 10 related to the $i$th sample combination.

Definition 3.2** (Vector and Matrix Notation).**

*We refer to the representation (18) and (20) as the vector and the matrix notation, respectively. In (18), the unknown is a vector whereas in (20), the unknown is a matrix. Even though both equations express the same, for the theoretical proofs in Section 6, the vector notation is more suitable since considering spaces that are spanned by matrices is rather uncommon. On the other hand, software implementations exploit the low-rank structure of the matrix $X$ in (20). This is why the matrix notation fits better in these cases. *

In (11), the right hand side $b_{D}$ does not depend on any parameter and $A(\mu_{s}^{i_{1}},\lambda_{s}^{i_{2}},\rho_{f}^{i_{3}})$ depends linearly on each parameter. Assume that $A(\mu_{s}^{i_{1}},\lambda_{s}^{i_{2}},\rho_{f}^{i_{3}})$ is invertible for all parameters. In this case, Theorem 3.6 of [9] is applicable and provides existence of a low-rank approximation of $X$ in (20) with an error bound that implies a stronger error decay than any polynomial in the rank $R$. However, the constant $C$ in Theorem 3.6 of [9] can become very big but we do not want to go into detail here and refer the interested reader to [9].

3.2 Preconditioners

The system matrix $\mathcal{A}$ has the structure

[TABLE]

Promising choices of preconditioners that were used already in [9] are

[TABLE]

or

[TABLE]

with the means

[TABLE]

The preconditioner $\mathcal{P}_{T}$ usually provides faster convergence than $\mathcal{P}_{A_{0}}$, especially if the means $\bar{\mu}_{s},\bar{\lambda}_{s}$ and $\bar{\rho}_{f}$ are big. Left multiplication of $\mathcal{P}_{T}^{-1}$ with

[TABLE]

is equivalent to application of $P_{T}^{-1}$ to $F(X)$ from the left using the matrix notation from Eq. 20.

4 The Low-rank Methods

Now, we discuss iterative methods that can be applied to solve the big system Eq. 18. The iterate is then a vector $x\in\mathbb{R}^{Mm}$. But if the iterate is represented as a matrix instead of a vector, computation can be kept in the matrix notation from Eq. 20. For instance, the matrix-vector multiplication in such a global approach corresponds to the evaluation of the function $F(\cdot)$ from Eq. 20. The Euclidean norm of the vector

[TABLE]

from Eq. 18 then equals the Frobenius norm of the matrix $X$ in Eq. 20, $\|X\|_{F}$. Low-rank methods that use this approach can be based on many methods such as the Richardson iteration or the conjugate gradient method as discussed in Algorithm 1 and Algorithm 2 of [9]. But since for fluid-structure interaction problems, the matrix $\mathcal{A}$ is not symmetric, the focus in this paper lies on methods that base on the GMRES and the Chebyshev method. As proved in Theorem 35.2 of [17], the GMRES method converges in this case, and so does the Chebyshev method, if all eigenvalues of the system matrix lie in an ellipse that does not touch the imaginary axis as proved in [11]. Also, the Bi-CGSTAB method from [18] is considered for a numerical comparison.

As mentioned, the low-rank methods discussed in this paper use an iterate that is, instead of a matrix, a tensor of order two. The iterate is then given by

[TABLE]

where the tensor rank $R\in\mathbb{N}$ is kept small such that $R\ll M,m$. The goal of the method is to find a low-rank approximation $\hat{X}$ that approximates the matrix $X$ in Eq. 20. The methods GMREST (also mentioned in [1]) and ChebyshevT are such methods and will be explained in the following. They are not just faster than the standard methods applied to $m$ individual equations of the form Eq. 11, they also need a smaller amount of storage to store the approximation. If $M$ and $m$ are very big, this plays an important role since the storage amount to store $\hat{X}$ is in O\big{(}(M+m)R\big{)} while the storage amount to store the full matrix $X$ is in $O(Mm)$.

4.1 Tensor Format and Truncation

There are several formats available to represent the tensor $\hat{X}$. For $d=2$, the hierarchical Tucker format (Definition 11.11 of [7]) is equivalent to the Tucker format. It is based on so called minimal subspaces that are explained in Chapter 6 of [7].

Definition 4.1** (Tucker Format - Definition 8.1 of [7] for $d=2$).**

Let $V:=\mathbb{R}^{M}\otimes\mathbb{R}^{m}$, $(r_{1},r_{2})\in\mathbb{N}^{2}$. For $d=2$, the Tucker tensors of Tucker rank $(r_{1},r_{2})$ are given by the set

[TABLE]

*From now on, the set $T_{(R,R)}(V)$ will be denoted by $T_{R}$. By a tensor of rank $R$, a Tucker tensor in $T_{(R,R)}$ is addressed in the following. *

As explained in Section 13.1.4 of [7], summation of two arbitrary Tucker tensors of rank $R$, in general, results in a Tucker tensor of rank $2R$. But to keep a low-rank method fast, the rank of the iterate has to be kept small. This induces the need for a truncation operator.

Definition 4.2** (Truncation Operator).**

The truncation operator

[TABLE]

*maps a Tucker or a full tensor into the set of Tucker tensors of rank $R$. The truncation operator is, in the case of tensors of order $2$, based on the singular value decomposition and projects its arguments to $T_{R}$. For further reading, we refer to Definition 2.5 of [6]. *

Remark 4.3**.**

As proved in Section 3.2.3 of [7], it holds

[TABLE]

where the relation $\cdot\cong\cdot$ denotes spaces that are isomorphic to each other (see Section 3.2.5 of [7]). Since, for our purposes, we consider a matrix that is represented by a tensor, we assume

[TABLE]

and if

[TABLE]

*by $\hat{x}$, the full representation of the tensor in $\mathbb{R}^{Mm}$ in vector notation is addressed. *

Before we proceed, one more definition is needed.

Definition 4.4** (Vectorization restricted to $\mathbb{R}^{M\times m}$).**

The vectorization operator

[TABLE]

stacks matrix entries column wise into a vector. Its inverse maps to an $M\times m$ matrix:

[TABLE]

Remark 4.5**.**

*The argument of the function $F(\cdot)$ from Eq. 20 is tacitly assumed to be a matrix so $F(\hat{x})$ addresses F\big{(}\operatorname{vec}^{-1}(\hat{x})\big{)} for $\hat{x}\in T_{R}$. *

Since truncation is an operation that is applied after nearly every addition of tensors and multiple times in every iteration, the format that provides the truncation with the least complexity is often the preferred one. According to Algorithm 6 of [10], the htucker toolbox [10] for MATLAB® provides truncation with complexity $(2\max(M,m)R^{2}+2R^{4})$ if the input format is in hierarchical Tucker format. The truncation complexity of the TT toolbox [12] for MATLAB that uses the Tensor Train format is in $O(2\max(M,m)R^{3})$ as stated in Algorithm 2 of [12].

4.2 The GMREST and the GMRESTR Method

Consider $\mathcal{A}$ from Eq. 18, a suitable preconditioner $\mathcal{P}=I_{m}\otimes P\in\mathbb{R}^{Mm\times Mm}$, a start vector $x_{0}\in\mathbb{R}^{Mm}$ and

[TABLE]

$l$ GMRES iterations with the preconditioner $\mathcal{P}$ applied to the system Eq. 18 minimize $\|r_{0}-\mathcal{P}^{-1}\mathcal{A}z\|_{2}$ for $z\in\mathbb{R}^{Mm}$ over the Krylov subspace (compare Section 6.2 of [14])

[TABLE]

As mentioned before, from the theoretical point of view, this classical GMRES method is equivalent to the global GMRES method that uses an iterate that is a matrix instead of a vector. But if the iterate is represented by a tensor of a fixed rank $R$, the truncation operator $\mathcal{T}$ generates an additional error every time it is applied to truncate the iterate or tensors involved back to rank $R$. With an initial guess $\hat{x}_{0}:=\mathcal{T}(x_{0})$,

[TABLE]

$l$ iterations of the truncated GMRES method GMREST that is coded in Algorithm 1 minimize \|\operatorname{vec}\Big{(}\mathcal{T}\big{(}\hat{r}_{0}-P^{-1}F(\hat{z})\big{)}\Big{)}\|_{2} for $\hat{z}\in T_{R}$ over the truncated Krylov subspace

[TABLE]

Algorithm 1 is a translation of a preconditioned variant of Algorithm 6.9 of [14] to the low-rank framework. The Arnoldi iteration is used to compute an orthogonal basis of $\mathcal{K}_{l}^{\mathcal{T}}$. The operations involved are translated from the vector notation to the matrix notation. Therefore, the Euclidean norm of a vector translates to the Frobenius norm of a matrix. The scalar product of two vectors corresponds to the Frobenius scalar product of two matrices,

[TABLE]

But even the standard GMRES method can stagnate due to machine precision. This means that at the $l$th iteration, the dimension of the numerical approximation of $\mathcal{K}_{l}$ is smaller than $l$. As we will see later, the truncation operator brings, in addition to the finite precision error (round-off), a truncation error into play. As a result, the GMREST method can stagnate much earlier than the non truncated full approach. As in the full approach, restarting the method with the actual iterate as initial guess can be a remedy. This restarted variant of the GMREST method, called GMRESTR here, is coded in Algorithm 2.

4.3 The ChebyshevT Method

The Chebyshev method converges for non-symmetric system matrices if, in the complex plane, the eigenvalues can be encircled by an ellipse that does not touch the imaginary axis.

The diagonal blocks of the preconditioned system matrix $\mathcal{P}_{T}^{-1}\mathcal{A}$ are

[TABLE]

Moreover, the parameter-dependent matrices Eq. 11 are assumed to be invertible. The eigenvalues of $\mathcal{P}_{T}^{-1}\mathcal{A}$ denoted by $\Lambda(\mathcal{P}_{T}^{-1}\mathcal{A})$ therefore coincide with the set

[TABLE]

In numerical tests, it turned out that $R_{\Lambda}:=\{\operatorname{Re}(\alpha):\alpha\in\Lambda(\mathcal{P}_{T}^{-1}\mathcal{A})\}\subset(0,\infty)$ for the linear fluid-structure interaction problems considered and the maximum and the minimum of $R_{\Lambda}$ do not depend on the number of degrees of freedom, where the operator $\operatorname{Re}(\alpha)$ returns the real part of a complex number $\alpha\in\mathbb{C}$. Unfortunately, so far it is not clear how to derive a useful bound for $R_{\Lambda}$ away from [math] and from above. For a discretization, the quantities

[TABLE]

can therefore be computed using the representation Eq. 22 of $\Lambda(\mathcal{P}_{T}^{-1}\mathcal{A})$ for a small number of degrees of freedom. Since we are using the mean-based preconditioner, the elements in $\Lambda(\mathcal{P}_{T}^{-1}\mathcal{A})$ lie symmetrically around $x=1$ in the complex plane \big{(}compare (21)\big{)}. Consider the ellipse with center

[TABLE]

The imaginary parts of the the elements in $\Lambda(\mathcal{P}_{T}^{-1}\mathcal{A})$ are so small such that this ellipse encircles all eigenvalues of $\mathcal{P}^{-1}\mathcal{A}$. Moreover, it does not touch the imaginary axis since $\Lambda_{\text{min}}>0$. The Chebyshev method from [11] can therefore be generalized in the same manner as the GMRES method in Subsection 4.2 and used to find a low-rank approximation $\hat{X}$ of $X$ in Eq. 20. The resulting truncated Chebyshev variant ChebyshevT is coded in Algorithm 3.

5 Time Discretization

5.1 The Linear Fluid-structure Interaction Problem

Let $[0,T]$ be a time interval for $T\in\mathbb{R}^{+}$ and $t\in[0,T]$ be the time variable. The deformation $u$ and the velocity $v$ now depend, in addition, on the time variable $t$ so we write $u(t,x)$ and $v(t,x)$. With the solid density $\rho_{s}\in\mathbb{R}$, the non stationary Navier-Lamé equations discussed in Section 2.3.1.2 of [13] fulfill

[TABLE]

The time term $\rho_{f}\partial_{t}v$ coming from the Stokes fluid equations as mentioned in (2.42) of [13] is added to the left side of the momentum equation. The weak formulation of the non-stationary coupled linear fluid-structure interaction problem is given by

[TABLE]

with regularity conditions v\in L^{2}\big{(}[0,T];v_{\text{in}}+H_{0}^{1}(\Omega,\Gamma_{f}^{D}\cup\Gamma_{\text{int}})^{d}\big{)}, $\partial_{t}v\in$ L^{2}\big{(}[0,T];H^{-1}(\Omega)^{d}\big{)} for all $(t,x)\in[0,T]\times\Omega$. We use the notation from Eq. 10.

5.2 Time Discretization With the $\boldsymbol{\theta}$-Scheme

Let $A_{t}^{f},A_{t}^{s}\in\mathbb{R}^{M\times M}$ be discretization matrices:

[TABLE]

Now consider a discretization that splits the time interval $[0,T]$ into $s+1\in\mathbb{N}$ equidistant time steps. Let the distance between two time steps be $\Delta_{t}$. The starting time is $t_{0}=0$ and the following times are thus given by $t_{i}:=i\Delta_{t}$ for $i\in\{1,...,s\}$. Let $X^{i}$ be the approximate solution at time $t_{i}$, $X^{0}$ is given as the initial value. The given Dirichlet data $b_{D}^{i}$ at time $t_{i}$ for all $i\in\{0,...,s\}$ yield the time dependent right hand side

[TABLE]

Consider the one-step $\theta$-scheme explained in Section 4.1 of [13]. Using the notation from Eq. 20, at time $t_{i}$, the following equation is to be solved for $X^{i}$:

[TABLE]

where $\theta\in[0,1]$. $F^{i}(\cdot)$ contains only two sum terms more than $F(\cdot)$ from Eq. 20. At time $t_{i}$, both Algorithm 1 and Algorithm 3 can be applied to the quasi stationary problem Eq. 25 with $F^{i}(\cdot)$ instead of $F(\cdot)$ and the right hand side $B^{i}(X^{i-1})$.

5.3 Preconditioner

At all time steps, the full matrix is given by

[TABLE]

The mean-based preconditioner, similar to $\mathcal{P}_{T}$ from Subsection 3.2, is

[TABLE]

Even though the right hand side $B^{i}(X^{i-1})$ changes with every time step, the system matrix does not.

6 Theoretical Error Bounds

The convergence proofs of the GMRES method from Theorem 35.2 of [17] and Section 3.4 of [15] base on the fact that the residual of the $l$th GMRES iterate can be represented as a product of a polynomial in $\mathcal{A}$ and the initial residual since the $l$th GMRES iterate is a linear combination of the start vector $x_{0}$ and the generating elements of $\mathcal{K}_{l}$. Also, the error bound of the Chebyshev method in [4] relies on the fact that the residual of the $l$th Chebyshev iterate is such a product. But even if one considers Algorithm 1 and Algorithm 3 in a non preconditioned version, multiplication with the system matrix $\mathcal{A}$ is always disturbed due to the error induced by the truncation operator. The GMREST method minimizes over $\mathcal{K}_{l}^{\mathcal{T}}$, the truncated Krylov subspace, instead of $\mathcal{K}_{l}$. In Subsection 6.2, the basis elements of $\mathcal{K}_{l}^{\mathcal{T}}$ are represented explicitly taking the truncation accuracy into consideration. Let $x_{l}$ be the $l$th GMRES iterate, $\hat{x}_{l}$ be the $l$th GMREST iterate. An upper bound of

[TABLE]

is derived from the accuracy of the basis elements of $\mathcal{K}_{l}^{\mathcal{T}}$. In relation to Krylov subspace methods, inaccuracies induced by matrix-vector multiplication result in so called inexact Krylov methods and have been discussed in [16]. Iterative processes that involve truncation have been discussed in a general way in [8]. For the $l$th Chebyshev iterate $x_{l}$ and the $l$th ChebyshevT iterate $\hat{x}_{l}$, the error

[TABLE]

is bounded in the same way in Subsection 6.3. These bounds show how the truncation error is propagated iteratively in Algorithm 1 and Algorithm 3 if the machine precision error is neglected.

Remark 6.1**.**

If $v\in\mathbb{R}^{Mm}$, $\mathcal{T}(v)$ addresses $\mathcal{T}(\operatorname{vec}^{-1}(v))$. Thus for the ease of notation, the truncation operator $\mathcal{T}$ from Definition 4.2 is regarded as a map

[TABLE]

*and for $v\in\mathbb{R}^{Mm}$, $\mathcal{T}(v)$ addresses the full representation of the tensor in vector notation, a vector in $\mathbb{R}^{Mm}$. *

Definition 6.2** (Truncation accuracy).**

The truncation operator $\mathcal{T}$ from Definition 4.2 is said to have accuracy $\epsilon>0$ if for any $x\in\mathbb{R}^{Mm}$

[TABLE]

*holds. $\mathcal{E}_{\hat{x}}$ is the error induced by $\mathcal{T}$ when $x$ is truncated. *

6.1 Matrix-Vector Product Evaluation Accuracy

If a tensor is multiplied with a scalar or a matrix, there is no truncation needed since the tensor rank does not grow. But the evaluation of $F(\cdot)$ from Eq. 20 involves 4 sum terms. After an evaluation of $F(\cdot)$ with a tensor as argument, the result has to be truncated. To keep complexity low for $\hat{X}\in T_{R}$, the sum \mathcal{T}\big{(}F(\hat{X})\big{)}, in practice, is truncated consecutively

[TABLE]

$\mathcal{E}_{\hat{F}_{s_{i}}}$ denotes the truncation error induced by the truncation of the $i$th sum term for $i\in\{1,2,3\}$. By Definition 4.2, $\|\mathcal{E}_{\hat{F}_{s_{i}}}\|_{2}\leq\epsilon$ for all $i\in\{1,2,3\}$. In $\mathcal{T}(F(\cdot))$ are, if the number of summands in $F(\cdot)$ is $K\in\mathbb{N}$, a total of $K-1$ truncations hidden. For a truncation accuracy of $\epsilon>0$ we have

[TABLE]

Since $K$ is a small number, usually not bigger than $4$, we will neglect this detail and simply assume

[TABLE]

in the following. To make sure that the stated error bounds are still valid, the truncation accuracy would be asked to, to be exact, less than $\frac{\epsilon}{K-1}$.

6.2 GMREST Error Bounds

Let $x_{l}$ be the $l$th standard GMRES iterate, $\hat{x}_{l}$ be the $l$th GMREST iterate. How big is the difference between the truncated Krylov subspace $\mathcal{K}_{l}^{\mathcal{T}}$ from Subsection 4.2 and the Krylov subspace $\mathcal{K}_{l}$? First we derive explicit representations of the non-normalized basis elements of $\mathcal{K}_{l}^{\mathcal{T}}$. For the following Lemma, we need the truncation errors $\mathcal{E}_{K^{\mathcal{T}_{k}}}\in\mathbb{R}^{Mm}$ for $k\in\mathbb{N}$. They are induced by the truncation operator when the $k$th basis element of the truncated Krylov subspace $\mathcal{K}_{l}^{\mathcal{T}}$ is computed. For a truncation operator with accuracy $\epsilon>0$, it holds $\|\mathcal{E}_{\mathcal{K}^{\mathcal{T}_{k}}}\|_{2}\leq\epsilon$ for all $k\in\mathbb{N}$.

Lemma 6.3** (Basis Representation of $\mathcal{K}_{l}^{\mathcal{T}}$).**

Assume $\operatorname{dim}(\mathcal{K}_{l}^{\mathcal{T}})=l$ and

[TABLE]

Let the truncation operator $\mathcal{T}(\cdot)$ have accuracy $\epsilon>0$. The non-normalized basis elements of $\mathcal{K}_{l}^{\mathcal{T}}$ are given by

[TABLE]

Proof 6.4**.**

*(by induction)

For $k=1$*

[TABLE]

and “$k-1\Rightarrow k$” since

[TABLE]

Remark 6.5** ((Truncation Error of $\hat{r}_{0}$)).**

Consider the line

[TABLE]

of Algorithm 1. In vector notation,

[TABLE]

is solved for $\hat{r}_{0}$. Usually the initial vector $x_{0}$ is chosen such that it can be represented by a tensor of low rank. So we assume $\hat{x}_{0}=\mathcal{T}(x_{0})=x_{0}$. If the linear system Eq. 26 is solved before truncation we have

[TABLE]

But this is rarely implemented this way. In practice, the right-hand side of Eq. 26 is truncated before the linear system is solved for $\hat{r}_{0}$. In this case,

[TABLE]

*holds. The following statements refer to the latter, more practical case. The error bounds that result in the first case can be found in Subsection A.1. *

Remark 6.6** ((Truncation of $\hat{W}$ and Orthogonality)).**

Consider the line

[TABLE]

in Algorithm 1. In the lemma above, this truncation is neglected. When the $k$th basis element $\hat{V}_{k}$ is set up, there are $k$ extra additions involved due to this line. Let $\mathcal{E}_{\hat{W}}$ be the truncation error that occurs when this line is executed. For the sake of readability, we neglect that they differ from loop iteration to loop iteration. As a consequence, we do not add another index to $\mathcal{E}_{\hat{W}}$. The basis elements are then given by

[TABLE]

Furthermore, we neglect round-off errors incurred from finite precision arithmetic as these are assumed to be much smaller than the truncation errors. The reason why the basis elements of $\mathcal{K}_{l}$ and $\mathcal{K}_{l}^{\mathcal{T}}$ obtained from the Arnoldi iteration differ from each other is the truncation error. Even though the elements

[TABLE]

*span $\mathcal{K}_{l}$, they are not orthogonal. But for the sake of notation, we incorporate the error made at the orthogonalization of the basis elements, in the truncated case, into $\mathcal{E}_{\hat{W}}$ and address by Eq. 27 the normalized basis elements that result from the Arnoldi iteration. In other words, we tacitly assume that the basis elements Eq. 27 of $\mathcal{K}_{l}$ are orthonormal, write them in the representation Eq. 27 and incorporate the error we made at orthogonalization into $\mathcal{E}_{\hat{W}}$. This is just one result of the assumption that we use exact precision. *

Lemma 6.7** (Error Bound for Truncated Basis Elements).**

Let $\sigma_{\mathcal{P}}:=\|\mathcal{P}^{-1}\mathcal{A}\|_{2}$. Under the assumptions of Lemma 6.3, for

[TABLE]

it holds that

[TABLE]

Proof 6.8**.**

For $k=0$, we have

[TABLE]

with the convention

[TABLE]

For $k\geq 1$, we use Lemma 6.3.

[TABLE]

*The truncation error coming from the orthogonalization process mentioned in Remark 6.6 adds the term $\epsilon k$ to the error bound. *

The standard GMRES minimizes over the Krylov subspace $\mathcal{K}_{l}$. In terms of Remark 6.6, the standard GMRES method finds coefficients $c_{i}\in\mathbb{R}$ for $i\in\{1,...,l\}$ such that

[TABLE]

In the same way we can write

[TABLE]

where the coefficients $d_{i}$ for $i\in\{1,...,l\}$ refer to the coefficients found by the Arnoldi iteration in the GMREST method. This allows to state the following theorem.

Theorem 6.9** (Approximation Error of GMREST).**

Let $x_{l}$ be the $l$th iterate of the standard GMRES method, $\hat{x}_{l}$ be the $l$th iterate of the GMREST method. It holds

[TABLE]

Proof 6.10**.**

[TABLE]

We assume that the standard GMRES method does an accurate orthogonalization of the Krylov subspace $\mathcal{K}_{l}$ (see Remark 6.6). By the elements $(\star)$ we address the orthonormal basis elements of $\mathcal{K}_{l}$. They all have an Euclidean norm of $1$. Therefore,

[TABLE]

*holds. The additional sum term $\epsilon l$ comes from the last successive sum in the method where the approximation is built. *

6.3 ChebyshevT Error Bounds

Similar to Subsection 6.2, we derive an error bound for the ChebyshevT method coded in Algorithm 3. Let $x_{l}$ denote the $l$th iterate of the standard Chebyshev method and $\hat{x}_{l}$ denote the $l$th iterate of the ChebyshevT method.

Remark 6.11**.**

The $i$th residual is given by the solution $r_{i}$ to

[TABLE]

The truncation of $r_{i}$ yields

[TABLE]

*if the truncation operator $\mathcal{T}$ is assumed to have accuracy $\epsilon$. In analogy to Remark 6.5, the two cases $\|\mathcal{E}_{\hat{r}_{i}}\|_{2}\leq\epsilon$ and $\|\mathcal{E}_{\hat{r}_{i}}\|_{2}\leq\epsilon\|\mathcal{P}^{-1}\|_{2}$ have to be distinguished. In this section, we consider the latter case. The error bounds of Theorem 6.14 for the case $\|\mathcal{E}_{\hat{r}_{i}}\|\leq\epsilon$ can be found in Subsection A.2. *

The start vector and the right hand side are assumed to be of low rank, namely

[TABLE]

In the same way as in Subsection 6.2, the norm $\|\hat{x}_{l}-x_{l}\|_{2}$ is to be estimated. $\mathcal{E}_{\hat{x}_{l}^{a}}$ denotes the total error

[TABLE]

not to be confused with $\mathcal{E}_{\hat{x}_{l}}$, the truncation error with norm $\epsilon$ that occurs when truncating $\hat{x}_{l}$. The iterative Chebyshev method is a three term recursion. Thus, the Chebyshev iterates itself can be represented by a recursive formula.

Lemma 6.12** (Representation of the ChebyshevT Iterates).**

Let the scalars

[TABLE]

be given as defined in Algorithm 3. $\Phi_{i}$ denote the non truncated full matrices corresponding to $\hat{\Phi}_{i}$ if Algorithm 3 is applied and any truncation is neglected. If

[TABLE]

it holds that

[TABLE]

where $\mathcal{E}_{\hat{x}_{j}^{a}}:=\hat{x}_{j}-x_{j}$ for $j\in\{0,...,l\}$. We use the convention

[TABLE]

*If a truncation operator of accuracy $\epsilon>0$ is used, then certainly $\|\mathcal{E}_{\hat{x}_{i}}\|_{2}\leq\epsilon$ but not necessarily $\|\mathcal{E}_{\hat{x}_{i}^{a}}\|_{2}\leq\epsilon$ holds for $i\in\{0,...,l\}$. The error induced by the truncation operator that truncates $\hat{\Phi}_{i}$ is denoted by $\mathcal{E}_{\hat{\Phi}_{i}}$ for $i\in\{0,...,l\}$. *

Proof 6.13**.**

$l=1$*:

Provided that $\hat{x}_{0}=x_{0}=x_{0}+\mathcal{E}_{\hat{x}_{0}^{a}}\Rightarrow\mathcal{E}_{\hat{x}_{0}^{a}}=0$.*

[TABLE]

$l=2$:

[TABLE]

since

[TABLE]

Thus,

[TABLE]

For the proof for $l\geq 3$, we go by induction. For the initial step $l=3$, we need

[TABLE]

and

[TABLE]

Therefore,

[TABLE]

To conclude $l-1\rightarrow l$ we first prove that

[TABLE]

under the assumption that this equation holds for $\hat{\Phi}_{l-2}$. For $\hat{\Phi}_{2}$, this is true since from Eq. 28, we have that

[TABLE]

The induction step for Eq. 29 is as follows.

[TABLE]

With this, it follows that

[TABLE]

Theorem 6.14** (Approximation Error of ChebyshevT).**

Let $\sigma_{\mathcal{P}}:=\|\mathcal{P}^{-1}\mathcal{A}\|_{2}$. Under the assumptions of Lemma 6.12, the following error bounds hold for a truncation operator of accuracy $\epsilon>0$.

[TABLE]

Proof 6.15**.**

*Let $\epsilon_{R}>0$ such that $\|\mathcal{E}_{\hat{r}_{i}}\|_{2}\leq\epsilon_{R}$ for all $i\in\{1,...,l\}$.

$l=1$:*

[TABLE]

$l=2$:

[TABLE]

$l\geq 3$:

[TABLE]

*The estimation $\epsilon_{R}\leq\epsilon\|\mathcal{P}^{-1}\|_{2}$ leads to the claimed error bounds. *

7 Numerical Evaluation of the Error Bounds

In algorithm and software implementations, the accuracy of a truncation operator depends on the truncation rank. If one chooses a rank $R$, the iterate of the GMREST or the ChebyshevT method is truncated to, the accuracy of the truncation operator is still unknown. Most truncation techniques like the HOSVD for hierarchical Tucker tensors (Section 8.3 and Section 10.1.1 of [7]) or the TT-rounding for TT tensors (Algorithm 1 and 2 of [12]) provide quasi optimality for tensors of order $d>2$. For tensors of order $d=2$ they even provide optimality in the sense that the result of the truncation of a matrix to rank $R$ is indeed the best rank $R$ approximation of the matrix. Nonetheless, since the singular value decay of the argument to be truncated is, in general, not known, the truncation operator will be simulated for a numerical evaluation of the error bounds of Theorem 6.9 and Theorem 6.14. Using the MATLAB routine rand(), a vector

[TABLE]

with entries that are uniformly distributed in the interval $(0,1)$ is constructed first. The argument $x\in\mathbb{R}^{Mm}$ is then truncated using the truncation simulator

[TABLE]

Of course, $\tilde{z}$ is computed anew every time $\mathcal{T}_{s}(\cdot)$ is applied. For this subsection, all the computations are therefore made in the full format and whenever a truncation operator is applied, the truncation simulator $\mathcal{T}_{s}(\cdot)$ is evaluated. The main advantage of this strategy is that

[TABLE]

A truncation operator based on the singular value decomposition does not provide such a reliable behavior. Let

[TABLE]

be the singular values of $\operatorname{vec}^{-1}(x)$ and $\exists k\in\{1,...,\min\{M,m\}-1\}$ such that, e.g.,

[TABLE]

A truncation operator with accuracy $\epsilon=10^{-5}$ based on the singular value decomposition would provide an approximation of $x$ with accuracy $10^{-10}$ in this example. So in this sense, the truncation simulator $\mathcal{T}_{s}$ yields the worst case error every time it is applied.

7.1 GMREST Error Bound

We consider the error bound from Theorem 6.9 that reads

[TABLE]

This theoretical error bound is compared with

[TABLE]

where $x_{l}$ denotes the $l$th GMRES iterate and $\hat{x}_{l}$ the $l$th GMREST iterate. As just explained, everything is computed in the full format and every time a truncation is involved (which affects the GMREST iterate $\hat{x}_{l}$ only), $\mathcal{T}_{s}$ from Eq. 30 is evaluated. The $3$d jetty from Subsection 8.1 is considered with size $M=4095$ and a three parameter discretization with a total of $m=8000$ parameter combinations as used in Subsection 8.3. We use the estimate $\sigma_{\mathcal{P}}\approx d+c$ with $c,d$ from Subsection 8.5. In addition, the basis element error bound from Lemma 6.7 is plotted for a truncation accuracy of $\epsilon=10^{-12}$. If one starts with a matrix whose entries are all set to $1$, the error bound Eq. 31 states that $\|x_{10}-\hat{x}_{10}\|_{2}$ is not bigger than $\approx 10^{-2}$, which can be seen in Fig. 1a. The reason for such a tolerant bound is that the first coefficients $d_{1},d_{2},...$ are very big if the initial guess is bad. But if both methods are restarted with $\hat{x}_{6}$ as start matrix, these coefficients become smaller as shown in Fig. 1b. Also, the relative residual norm of the GMRES iterate, $\frac{\|B-F(x)\|_{F}}{\|B\|_{F}}$, and the one of the GMREST iterate, $\frac{\|B-F(\hat{x})\|_{F}}{\|B\|_{F}}$, are plotted. So even though $\|x_{l}-\hat{x}_{l}\|_{2}$ stagnates, the residual of the truncated approach still decreases.

The dominating terms are

[TABLE]

As pointed out above, $d_{1},d_{2},...$ are big for a bad initial guess. Then, in addition, $\epsilon$ can not compensate the (exponential) growth of $\sigma_{\mathcal{P}}^{j-1}$ for $j\in\{1,...,l\}$. Notice that $\sigma_{\mathcal{P}}\approx 1.6$ in this example. Since $\|\mathcal{P}^{-1}\|_{2}$ is rather big, namely $\approx 3\cdot 10^{4}$ in this example, the former of these two terms is bigger for the first $10$ iterations. Furthermore, $1.6^{10}\approx 10^{2}$ and the moduli of the coefficients $d_{j}$ become smaller the bigger the iteration count is. This is why we do not see an exponential growth of the bound in Fig. 1a and Fig. 1b.

7.2 ChebyshevT Error Bound

In this subsection, the approximation error from Theorem 6.14 is numerically examined. Let $x_{l}$ be the $l$th Chebyshev iterate and $\hat{x}_{l}$ be the $l$th ChebyshevT iterate. For the ChebyshevT method, similar to the GMRES comparison, $\mathcal{T}_{s}$ is used as truncation operator.

Remark 7.1**.**

If preconditioned with $\mathcal{P}_{T}$, in theory,

[TABLE]

holds as mentioned in Remark 6.11. But, due to the bad condition of the preconditioner and machine precision, in practice, $\epsilon_{R}$ is often bigger. The line

[TABLE]

is executed in every iteration of Algorithm 3 which sometimes leads to real errors that are higher than the error bound from Theorem 6.14. In contrast to this, the line

[TABLE]

*in Algorithm 1 is less vulnerable. To circumvent this problem, for the error bound of Theorem 6.14, the error $\epsilon_{R}$ is computed explicitly. This value instead of $\epsilon\|\mathcal{P}_{T}^{-1}\|_{2}$ is then used to compute the theoretical error bounds that are compared with the real errors. *

We use the same configuration as used in Subsection 7.1. Even though the theoretical error bound literally explodes, for $\epsilon=10^{-12}$, the truncated method converges roughly as good as the non truncated method until iteration $10$ for the $3$d jetty model from Subsection 8.1 as shown in Fig. 2a. But the convergence of the ChebyshevT method deteriorates remarkably after $4$ iterations for $\epsilon=10^{-6}$ if compared to the full approach (see Fig. 2b).

The two terms in the error bound that are not multiplied with $\epsilon$ are

[TABLE]

The coefficients $\beta_{i}$ in the Chebyshev method have norms that are smaller than $1$. This becomes clear if one considers the recursive computation formula for the Chebyshev polynomials (see (2.4) in [11]) evaluated at $\frac{d}{c}$ with $|c|<|d|$. The coefficients $\beta_{i}$ are then given as a fraction where the numerator has a norm that is smaller than the denominator. The product $(\star)$ becomes smaller the higher the iteration number is and therefore, the term Eq. 33 becomes negligibly small, at least if it is compared with the term Eq. 32. For our configuration, $c=0.6$. Hence, the coefficients $\alpha_{i}$ are bigger than $1$ on the other hand (see (2.24) in [11]). For $\sigma_{\mathcal{P}}\approx c+d=1.6$

[TABLE]

This explains why the first term in Theorem 6.14, the term Eq. 32, dominates the error bound and makes it explode. Thus, when using ChebyshevT, one has to be very careful about the choice of the truncation tolerance.

8 Numerical Examples

8.1 A Three Dimensional Jetty in a Channel

The geometric configuration of a $3$d jetty in a channel is given by

[TABLE]

With the velocity

[TABLE]

the Dirichlet inflow on the left boundary is given by

[TABLE]

The geometric configuration is illustrated in Fig. 3. On the right, for $x=8$, do nothing boundary outflow conditions as discussed in Section 2.4.2 of [13] hold. The surface is at $y=8$. There, $v_{2}$ and $u_{2}$ vanish. Everywhere else on $\partial(\Omega)$, the velocity and the deformation fulfill zero Dirichlet boundary conditions.

8.2 Stabilization of Equal-order Finite Elements

Because we use $Q_{1}$ finite elements for the pressure, velocity and the deformation, we have to stabilize the Stokes fluid equations on the fluid part. For this, stabilized Stokes elements as explained in Lemma 4.47 of [13] are used.

8.3 Three Parameter Discretization

Problem Eq. 10 is discretized with respect to

[TABLE]

The kinematic fluid viscosity is fixed to $\nu_{f}=0.01$. The shear modulus and first Lamé parameter ranges cover solids with Poisson ratios between $\frac{1}{3}$ (e.g. concrete) and $0.43478$ (e.g. clay). The total number of equations is $m=20^{3}=8000$ and the number of degrees of freedom is $M=192423$.

In the following computations, MATLAB 2017b on a CentOS 7.6.1810 64bit with 2 AMD EPYC 7501 and 512GB of RAM is used. The htucker MATLAB toolbox [10] is used to realize the Tucker format $T_{R}$. The preconditioners are decomposed into a permuted LU decomposition using the MATLAB builtin command lu(). All methods start with a start matrix whose entries are all set to $1$.

8.4 GMREST

A standard GMRES approach is compared with the GMRESTR method from Algorithm 2.

By “standard GMRES approach”, the standard GMRES method applied to $m=8000$ different equations of the form Eq. 11 is meant. It is once restarted after $8$ iterations so it uses a total of $16$ iterations per equation. For all $8000$ separate standard GMRES methods, $5$ preconditioners given by

[TABLE]

are set up where the diagonal matrices $\{D_{j}\}_{j\in\{1,2,3\}}$ are the ones from Eq. 19.

The GMRESTR method uses $6$ iterations per restart and is restarted $3$ times. The mean based preconditioner $\mathcal{P}_{T}$ is used. The times to compute the preconditioners (one in the case of GMRESTR, $5$ in the case of standard GMRES) can be found in the column Precon. in Table 1. The method itself took the time that is listed in the column Comp. and the column Total is then the sum of these times. Both methods result in $8000$ approximations. Each of these approximations (x axis) then provides a certain accuracy (y axis) that is plotted in Fig. 4. The standard GMRES method applied to $8000$ equations in this way provides accuracies that are plotted in red within about $90$ hours and $32$ minutes. The approximations the GMRESTR method provides have accuracies that are plotted in blue. The GMRESTR method took only about $181$ minutes to compute these approximations as one can see in Table 1. Also, the storage that is needed to store the approximation varies significantly. The rank $200$ approximation, in the Tucker format, requires only about $306$MB whereas the full matrix requires almost $12$GB.

8.5 ChebyshevT

Before the Chebyshev method can be applied, the extreme eigenvalues of $\mathcal{P}_{T}^{-1}\mathcal{A}$ have to be estimated as explained in Subsection 4.3. An estimation of $\Lambda_{\text{max}}$ and $\Lambda_{\text{min}}$ from Eq. 23 involves the estimation of extreme eigenvalues for $m$ different matrices if the representation Eq. 22 is considered. But we restrict to an estimation of

[TABLE]

With the mean based preconditioner $\mathcal{P}_{T}$, this leads to $d=1$ and $c=0.6$ in this configuration. The time needed to compute $\bar{\Lambda}_{\text{max}}$ and $\bar{\Lambda}_{\text{min}}$ on a coarse grid with $M=735$ degrees of freedom is listed in the column “Est.” in Table 2.

In the same manner as in the preceding subsection, for comparison, a standard Cheyshev approach is applied to $8000$ equations of the form Eq. 11. The standard Chebyshev method uses $20$ iterations at each equation and, in total, the same $5$ preconditioners Eq. 34 as the standard GMRES uses. The ChebyshevT method iterates, in total, $24$ times and uses $\mathcal{P}_{T}$, the mean based preconditioner. The ChebyshevT method is restarted $3$ times with $6$ iterations per restart. Compared to this, $24$ iterations without restart took about the same time but provide approximation accuracies that are slightly worse.

8.6 Comparison With the Bi-CGstab Method

Another method that also works for non-symmetric matrices is the Bi-CGstab method [18]. It was not considered in the first place because it can break down under some circumstances as explained in Section 2.3.8 of [2]. The preconditioned truncated variant similar to Algorithm 3 of [9] but strictly based on [18] is compared with the GMRESTR and the ChebyshevT method. The truncated Bi-CGstab method is applied with $6$ iterations per restart. If once restarted, in total, the method iterates $12$ times. The resulting approximation accuracy is indeed better than the one obtained when iterating $12$ times directly without any restart.

To avoid early stagnation, the residual at step $i$ is computed directly

[TABLE]

The main reason why the Bi-CGstab takes more time, is the truncation. If $M$ and $m$ are big, truncation after every addition is indispensable. In this way, in every for-loop in the preconditioned Bi-CGstab Algorithm[18], a total of $12$ truncations occur ($3$ times, $F(\cdot)$ is evaluated). In a for-loop of Algorithm 3, we have $6$ truncations only. In the outer for-loop of Algorithm 1, $3+i$ truncations occur. It turned out that, in the implementation of the preconditioned Bi-CGStab method of [18], the truncation of the sum

[TABLE]

can be left out. Leaving out further truncations did not lead to a better performance of the truncated Bi-CGstab method.

9 Conclusions

The truncated methods discussed in this paper provide approximations with relative residual norms smaller than $10^{-8}$ within less than a twentieth of the time needed by the correspondent standard approaches that solve the $m$ equations individually. This raises the question how these methods perform when applied to nonlinear problems.

Since the truncation error affects, in addition to the machine precision error, the accuracy of the Arnoldi orthogonalization, the GMREST method should preferably be applied in a restarted version. Mostly, the ChebyshevT method is a bit faster and a bit more accurate than the GMREST method. But the main disadvantage of the ChebyshevT method is that the ellipse that contains the eigenvalues of $\mathcal{P}^{-1}\mathcal{A}$ described by the foci $d\pm c$ has to be approximated newly every time the parameter configuration changes. In this matter, the GMREST method can be seen as a method that is a bit more flexible if compared to the ChebyshevT method.

Also, the ChebyshevT and the truncated BiCGstab methods can and preferably should be applied in a restarted manner. If not restarted, the methods stagnate after a few iterations already. The reason is a numerical issue initiated by the bad condition of the mean-based preconditioner.

There is still further investigation needed regarding the error bounds. If the GMREST method is applied, the coefficients $c_{j}$ are not known. The ChebyshevT bound is rather pessimistic and merely of theoretical nature. The method seems to converge too fast such that the truncation error does not really play a role in the cases examined.

Acknowledgments

We gratefully acknowledge funding received by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - 314838170, GRK 2297 MathCoRe.

Appendix A Error Bounds

A.1 GMREST Error Bounds

Lemma A.1** (Error Bound for Truncated Basis Elements).**

Under the assumptions of Lemma 6.7, with $\|\mathcal{E}_{\hat{r}_{0}}\|_{2}\leq\epsilon$, it holds

[TABLE]

Proof A.2**.**

The case $k=0$ is clear. For $k\geq 1$, we have

[TABLE]

*The term $\epsilon k$ is to be added to the error bound as explained in Remark 6.6. *

Theorem A.3** (Approximation Error of GMREST).**

Under the assumptions of Theorem 6.9 and $\|\mathcal{E}_{\hat{r}_{0}}\|_{2}\leq\epsilon$, it holds

[TABLE]

Proof A.4**.**

We use the proof of Theorem 6.9 and the error bound of Lemma A.1 to estimate

[TABLE]

*Truncation in the last line of Algorithm 1 leads to the additional sum term $\epsilon l$. *

A.2 ChebyshevT Error Bounds

Theorem A.5** (Approximation Error of ChebyshevT).**

Under the assumptions of Theorem 6.14 and $\|\mathcal{E}_{\hat{r}_{i}}\|_{2}\leq\epsilon$ for all $i\in\{1,...,l\}$, the error bounds of the ChebyshevT method translate to

[TABLE]

Proof A.6**.**

*For all cases, we can use the proof of Theorem 6.14 and estimate$\epsilon_{R}\leq\epsilon$. *