1xx462022

# Conditioning of linear systems arising from penalty methods††thanks:
Partially supported by NSF Grant DMS 2110379.

William Layton222 Department of Mathematics, University of Pittsburgh,
Pittsburgh, PA 15260<EMAIL_ADDRESS>Shuxian Xu333Department of
Mathematics, University of Pittsburgh, Pittsburgh, PA 15260,
<EMAIL_ADDRESS>
###### Abstract

Penalizing incompressibility in the Stokes problem leads, under mild
assumptions, to matrices with condition numbers
$\kappa=\mathcal{O}(\varepsilon^{-1}h^{-2})$, $\varepsilon=$ penalty parameter
$<<1$, and $h=$ meshwidth $<1$. Although
$\kappa=\mathcal{O}(\varepsilon^{-1}h^{-2})$ is large, practical tests seldom
report difficulty in solving these systems. In the SPD case, using the
conjugate gradient method, this is usually explained by spectral gaps
occurring in the penalized coefficient matrix. Herein we point out a second
contributing factor. Since the solution is approximately incompressible,
solution components in the eigenspaces associated with the penalty terms can
be small. As a result, the effective condition number can be much smaller than
the standard condition number.

###### keywords:

penalty method, effective condition number

65F35, 15A12

##  1 Estimating conditioning

We dedicate this paper to Professor Owe Axelsson.

Penalty methods have advantages and disadvantages. Disadvantages include the
need to pick a specific value of the penalty parameter $\varepsilon$ and ill-
conditioning of the associated linear system. We show herein that, measured by
the system’s effective condition number, this ill-conditioning is not severe,
Theorem 2.1 below. This observation is developed herein for the standard
penalty approximation to the Stokes problem

$-\triangle u+\nabla p=f(x)\ \ \text{and}\ \ \nabla\cdot u=0,$

in a polygonal domain $\Omega$ with boundary condition $u=0$ on
$\partial\Omega$ and normalization $\int_{\Omega}pdx=0$. Let
$(\cdot,\cdot),||\cdot||$ denote the $L^{2}(\Omega)$ inner product and norm
and $|\cdot|$ the usual euclidean norm of a vector and matrix. $X^{h}$ denotes
a standard, conforming, velocity finite element space of continuous, piecewise
polynomials that vanish on $\partial\Omega$. The penalty approximation results
from replacing $\nabla\cdot u=0$ by $\nabla\cdot u^{\varepsilon}=-\varepsilon
p^{\varepsilon}$ and eliminating $p^{\varepsilon}$. It is: find
$u^{\varepsilon}\in X^{h}$ satisfying

(1) $(\nabla u^{\varepsilon},\nabla v)+\varepsilon^{-1}(\nabla\cdot
u^{\varepsilon},\nabla\cdot v)=(f,v)\text{ for all }v\in X^{h}.$

Picking a basis $\left\\{\phi_{1},\cdot\cdot\cdot,\phi_{N}\right\\}$ for
$X^{h}$ leads to a linear system with coefficient matrix

$A_{ij}=(\nabla\phi_{i},\nabla\phi_{j})+\varepsilon^{-1}(\nabla\cdot\phi_{i},\nabla\cdot\phi_{j}),i,j=1,\cdot\cdot\cdot,N.$

Standard condition number estimates for this system yield bounds like
$\kappa\leq C\varepsilon^{-1}h^{-2}$.

Recall that the standard condition number, [2, 24], introduced in [20, 23] and
still of interest, e.g., [22], measuring the correlation between relative
error and relative residual, the distance to the nearest singular matrix and
the difficulty, cost and worse-case accuracy in solving a linear system with
$A$, is $\kappa:=|A||A^{-1}|$. However, estimates of the above with $\kappa$
are often (but not always, [5]) too rough because they do not consider the
size of solution components in the matrix eigenspaces. Thus, generalizations
exist, such as the Kaporin condition number [16] and extensions based on
generalized inverses [12], that can be used to obtain better predictions. One
important, easy to compute and interpret generalization is the condition
number at the solution, also called the effective condition number.

######  Definition 1.1.

Let $Ac=b$. Then, $\kappa(c)$, the condition number at the solution $c$, is

$\kappa(c):=|A^{-1}|\frac{|Ac|}{|c|}.$

Clearly $\kappa(c)\leq\kappa$ and $\kappa(c)$ takes into account both the
spectrum and the magnitude of solution components across eigenspaces. It is
also known, e.g. [2], that the relative error in approximations is bounded by
$\kappa(c)$ times the relative residual. To our knowledge, this extension is
due to Chan and Foulser [9]. It was soon thereafter developed by Axelsson [1,
2] and has been further developed in important work in [3, 4, 7, 18, 19, 11].

Section 2 gives the proof that $\kappa(u^{\varepsilon})<<\kappa$. Section 3
presents consistent numerical tests.

##  2 Analysis of $\kappa(u^{\varepsilon})$

We assume $X^{h}$ satisfies the following two assumptions typical, e.g. [6,
10], of finite element spaces on quasi-uniform meshes.

A1: [Inverse estimate] For all $v\in X^{h},$ $||\nabla v||\leq Ch^{-1}||v||.$

A2: [Norm equivalence] Let
$v=\sum_{i=1}^{N}a_{i}\phi_{i}(x),N=Ch^{-d},d=\dim(\Omega)=2$ or $3$. Then
$||v||$ and $\sqrt{\frac{1}{N}\sum_{i=1}^{N}a_{i}^{2}}$ are uniform-in-$h$
equivalent norms.

######  Theorem 2.1.

Let A1 and A2 hold. Select $|\cdot|$ to be the euclidean norm. Let $f^{h}$ be
the projection of $f$ into the finite element space and $u^{\varepsilon}$ the
solution of (1). Then,

$\displaystyle\max_{f^{h}}\kappa(u^{\varepsilon})$ $\displaystyle=$
$\displaystyle\kappa\leq C(h^{-2}+\varepsilon^{-1}h^{-2})$
$\displaystyle\kappa(u^{\varepsilon})$ $\displaystyle\leq$ $\displaystyle
C\frac{||f^{h}||}{||u^{\varepsilon}||}\text{, and}$
$\displaystyle\kappa(u^{\varepsilon})$ $\displaystyle\leq$ $\displaystyle
Ch^{-2}\left(1+\frac{h}{\varepsilon}\frac{||\nabla\cdot
u^{\varepsilon}||}{||u^{\varepsilon}||}\right).$

######  Proof 2.2.

We first estimate $|A^{-1}|^{2}=\max_{b}|A^{-1}b|^{2}/|b|^{2}.$ Given an
arbitrary right-hand side $\overrightarrow{b}$, let the linear system for the
undetermined coefficients
$(c_{1},c_{2},\cdot\cdot\cdot,c_{N})^{T}=\overrightarrow{c}$ be denoted
$A\overrightarrow{c}=\overrightarrow{b}.$ We convert in a standard way this
linear system to an equivalent formulation similar to (1). Recall that the
mass matrix associated with this basis is

$M_{ij}=(\phi_{i},\phi_{j}),i,j=1,\cdot\cdot\cdot,N.$

The matrix $M$ is symmetric and positive definite. Under A1 and A2, its
eigenvalues satisfy $0<C_{1}N\leq\lambda(M)\leq C_{2}N$. Given the vector
$\overrightarrow{b}$ let $\overrightarrow{a}=M^{-1}\overrightarrow{b}$. By
construction

$g(x)=\sum_{i=1}^{N}a_{i}\phi_{i}(x)\text{ satisfies
}(g(x),\phi_{j})=b_{j},j=1,\cdot\cdot\cdot,N.$

By A2, $||g||$ and $(N^{-1}\sum a_{i}^{2})^{1/2}$ are uniformly equivalent
norms. The bounds on $\lambda(M)$ (also from A2) imply $(N^{-1}\sum
a_{i}^{2})^{1/2}$ and $(N^{-1}\sum b_{i}^{2})^{1/2}$ are also uniformly
equivalent norms. Next define

$w=\sum_{i=1}^{N}c_{i}\phi_{i}(x).$

By A2 again $||w||$, $(N^{-1}\sum c_{i}^{2})^{1/2}$are equivalent norms. Thus,

$|A^{-1}|^{2}=\max_{b}\frac{|A^{-1}b|^{2}}{|b|^{2}}\leq C\max_{g\in
X^{h}}\frac{||w||^{2}}{||g||^{2}}.$

By construction $w,g\in X^{h}$ satisfy

$(\nabla w,\nabla v)+\varepsilon^{-1}(\nabla\cdot w,\nabla\cdot v)=(g,v)\text{
for all }v\in X^{h}.$

Setting $v=w$ and using simple inequalities gives $||w||^{2}\leq C||g||^{2}$.
Indeed,

$C||w||^{2}\leq||\nabla w||^{2}+\varepsilon^{-1}||\nabla\cdot
w||^{2}=(g,w)\text{ }\leq\frac{C}{2}||w||^{2}+C||g||^{2}.$

Thus $||w||^{2}\leq C||\nabla w||^{2}\leq C||g||^{2}$ and
$||w||^{2}/||g||^{2}\leq C$. This implies $|A^{-1}|\leq C$, uniformly of $h,$
$\varepsilon$.

Next we estimate $|A\overrightarrow{c}|/|\overrightarrow{c}|$ where
$u^{\varepsilon}=\sum_{i=1}^{N}c_{i}\phi_{i}(x)$ is the solution of (1). By
norm equivalence, as above, and (1) this is equivalent to
$||f^{h}||/||u^{\varepsilon}||$ where $f^{h}$ is the $L^{2}$ projection of
$f(x)$ into the finite element space. This gives the first estimate
$\kappa(u^{\varepsilon})\leq C||f^{h}||/||u^{\varepsilon}||$. For the second
estimate, we have

$\displaystyle||f^{h}||$ $\displaystyle=$ $\displaystyle\max_{v\in
X^{h}}\frac{(f^{h},v)}{||v||}=\max_{v\in X^{h}}\frac{(\nabla
u^{\varepsilon},\nabla v)+\varepsilon^{-1}(\nabla\cdot
u^{\varepsilon},\nabla\cdot v)}{||v||}$ $\displaystyle\leq$
$\displaystyle\max_{v\in X^{h}}\frac{||\nabla u^{\varepsilon}||||\nabla
v||+\varepsilon^{-1}||\nabla\cdot u^{\varepsilon}||||\nabla\cdot v||}{||v||}$
(using A1) $\displaystyle\leq$ $\displaystyle\max_{v\in
X^{h}}\frac{Ch^{-2}||u^{\varepsilon}||||v||+\varepsilon^{-1}Ch^{-1}||\nabla\cdot
u^{\varepsilon}||||v||}{||v||}$ $\displaystyle\leq$ $\displaystyle
Ch^{-2}||u^{\varepsilon}||+\varepsilon^{-1}Ch^{-1}||\nabla\cdot
u^{\varepsilon}||\text{, which implies}$
$\displaystyle\frac{||f^{h}||}{||u^{\varepsilon}||}$ $\displaystyle\leq$
$\displaystyle Ch^{-2}+C\varepsilon^{-1}h^{-1}\frac{||\nabla\cdot
u^{\varepsilon}||}{||u^{\varepsilon}||}.$

This implies

$\kappa(u^{\varepsilon})\leq C\frac{||f^{h}||}{||u^{\varepsilon}||}\leq
C\left(h^{-2}+\varepsilon^{-1}h^{-1}\frac{||\nabla\cdot
u^{\varepsilon}||}{||u^{\varepsilon}||}\right).$

Using $||\nabla\cdot u^{\varepsilon}||\leq C||\nabla u||\leq
Ch^{-1}|||u^{\varepsilon}||$ and A1 yields the standard estimate $\kappa\leq
C(h^{-2}+\varepsilon^{-1}h^{-2}).$

##  3 An Illustration

The result in Section 2 gives 3 estimates of conditioning. We explore if the 3
estimates

$\displaystyle\kappa\leq C(h^{-2}+\varepsilon^{-1}h^{-2}),\text{ \ \
}\kappa(u^{\varepsilon})\leq C\frac{||f^{h}||}{||u^{\varepsilon}||},\text{
and}$ $\displaystyle\kappa(u^{\varepsilon})\leq
C\left(h^{-2}+\varepsilon^{-1}h^{-1}\frac{||\nabla\cdot
u^{\varepsilon}||}{||u^{\varepsilon}||}\right)$

yield significantly different predictions. We investigate this question for 2
test problems and 2 finite element spaces. The first test problem has a smooth
solution and the second has an $f(x,y)$ oscillating as fast as possible on the
given mesh. The domain is the unit square and the mesh is a uniform mesh of
squares each divided into 2 right triangles. The first finite element space is
P1, conforming linear elements. The second is P2, conforming quadratics. The
space of conforming linears does not contain a divergence-free subspace. Thus
the coefficient matrix $A$ is expected to show greater ill-conditioning as
$\varepsilon\rightarrow 0$ than with conforming quadratics. The above
constants "$C$" are $\mathcal{O}(1)$ constants, independent of $\varepsilon$
and $h$. Thus, we compute and compare the RHS’ below

$\displaystyle Est.1$ $\displaystyle\simeq$ $\displaystyle
h^{-2}+\varepsilon^{-1}h^{-2},$ $\displaystyle Est.2$ $\displaystyle\simeq$
$\displaystyle\text{ }\frac{||f^{h}||}{||u^{\varepsilon}||}\text{ and }$
$\displaystyle Est.3$ $\displaystyle\simeq$
$\displaystyle\varepsilon^{-1}h^{-1}\frac{||\nabla\cdot
u^{\varepsilon}||}{||u^{\varepsilon}||}.$

We computed these values starting with $\varepsilon=1,h=1$ and halving each
until $\varepsilon=5.96E-08$ and $h=0.00195312.$ We present the results below
selected from equi-spaced $\varepsilon$-values in this data.

Test problem 1: We solve this test problem with P1, conforming linear
elements, and then with P2, conforming, quadratic elements. Choose
$f(x,y)=\left(\sin\left(x+y\right),\cos\left(x+y\right)\right)^{T}$. For P1
elements we have the following results for $Est.2$.

$\displaystyle\begin{array}[]{c|cccccc}{\small\varepsilon\Downarrow\setminus
h\Rightarrow}\text{ \ }&{\small 0.125}&{\small 6.2E-2}&{\small 3.1E-2}&\text{
}{\small 1.6E-2}&{\small 7.8E-3}&{\small 3.9E-3}\\\ \hline\cr{\small
1}&{\small 39}&{\small 38}&{\small 37}&{\small 37}&{\small 37}&{\small 37}\\\
{\small 6.3E-2}&{\small 1.9E2}&{\small 1.7E2}&{\small 1.7E2}&{\small
1.7E2}&{\small 1.7E2}&{\small 1.7E2}\\\ {\small 3.9E-3}&{\small 1.8E3}&{\small
7.8E2}&{\small 4.3E2}&{\small 3.3E2}&{\small 3.0E2}&{\small 2.9E2}\\\ {\small
2.4E-4}&{\small 2.5E4}&{\small 7.5E3}&{\small 2.2E3}&{\small 8.0E2}&{\small
4.3E2}&{\small 3.3E2}\\\ {\small 1.5E-5}&{\small 4.0E5}&{\small 1.1E5}&{\small
2.9E4}&{\small 7.5E3}&{\small 2.1E3}&{\small 8.0E2}\\\ {\small 9.5E-7}&{\small
6.4E6}&{\small 1.8E6}&{\small 4.5E5}&{\small 1.1E5}&{\small 2.9E4}&{\small
7.4E3}\\\ {\small 6.0E-8}&{\small 1.0E8}&{\small 2.9E7}&{\small 7.2E6}&{\small
1.8E6}&{\small 4.5E5}&{\small 1.1E5}\\\ \end{array}$ $\displaystyle\text{
Table 1: Values of }Est.2\simeq\text{
}\frac{||f^{h}||}{||u^{\varepsilon}||}\text{, P1 elements}$

Down the columns (fixing $h$ and decreasing $\varepsilon$), the data for
$Est.2$ shows that this quantity grows as $\varepsilon\downarrow 0$, roughly
like $\varepsilon^{-1}$ for fixed $h$. Across the rows, for fixed
$\varepsilon$ and $h\downarrow 0$, $Est.2$ decreases. Diagonals (necessary if
we choose $\varepsilon\sim h$) show mild growth. Comparing the last row of the
table with the row of $Est.1$ values shows that $Est.2$ consistently yields a
lower estimate of ill-conditioning than $Est.1$. Concerning the growth of
$Est.2$ as $\varepsilon\downarrow 0$ for fixed $h$, since $h,||f||$ do not
change as $\varepsilon$ varies, this growth is due to
$||u^{\varepsilon}||\rightarrow 0$ as $\varepsilon\downarrow 0$. The P1 finite
element space does not have a non-trivial divergence free subspace, [15].
Thus, $\varepsilon\downarrow 0$ forces $||\nabla\cdot u||\rightarrow 0$ which
forces $||u||\rightarrow 0$. This indicates that ill-conditioning of penalty
methods with P1 elements as $\varepsilon\downarrow 0$ is an essential feature
caused by the lack of a divergence-free subspace. For comparison with the last
row in Table 1, the values of $Est.1\simeq h^{-2}+\varepsilon^{-1}h^{-2}$ for
$\varepsilon=6.0E-8$ are

$\begin{array}[]{cccccc}\text{ }{\small 1.0E9}&\text{ }{\small 4.3E9}\text{
}&{\small 1.7E10}&{\small 6.9E10}&{\small 2.7E11}&{\small 1.1E12}\end{array}.$

Clearly $Est.2$ provides a smaller estimate than $Est.1$. Another
interpretation is that a poor choice of finite element space (made to
accentuate ill-conditioning) makes the problem of selecting $\varepsilon$
acute.

Table 2 presents the analogous results for $Est.3$ and P1 elements.

$\displaystyle\begin{array}[]{c|cccccc}{\small\varepsilon\Downarrow\setminus
h\Rightarrow}&{\small 0.125}&{\small 6.2E-2}&{\small 3.1E-2}&\text{ }{\small
1.6E-2}&{\small 7.8E-3}&{\small 3.9E-3}\\\ \hline\cr{\small 1}&{\small
25}&{\small 50}&{\small 1.0E2}&{\small 2.0E2}&{\small 4.0E2}&{\small 8.0E2}\\\
{\small 6.3E-2}&{\small 3.5E2}&{\small 6.6E2}&{\small 1.3E3}&{\small
2.5E3}&{\small 5.0E3}&{\small 1.0E4}\\\ {\small 3.9E-3}&{\small 4.2E3}&{\small
4.1E3}&{\small 4.5E3}&{\small 6.3E3}&{\small 1.1E4}&{\small 2.0E4}\\\ {\small
2.4E-4}&{\small 6.2E4}&{\small 4.9E4}&{\small 4.2E4}&{\small 4.5E4}&{\small
3.5E2}&{\small 5.1E4}\\\ {\small 1.5E-5}&{\small 9.8E5}&{\small 7.6E5}&{\small
6.3E5}&{\small 4.2E4}&{\small 6.0E5}&{\small 6.4E5}\\\ {\small 9.5E-7}&{\small
1.5E7}&{\small 1.2E7}&{\small 1.0E7}&{\small 9.5E6}&{\small 9.4E6}&{\small
9.4E6}\\\ {\small 6.0E-8}&{\small 2.5E8}&{\small 1.9E8}&{\small 1.6E8}&{\small
1.5E8}&{\small 1.5E8}&{\small 1.5E8}\end{array}$ $\displaystyle\text{Table 2:
Values of }Est.3\simeq\varepsilon^{-1}h^{-1}\frac{||\nabla\cdot
u^{\varepsilon}||}{||u^{\varepsilon}||}\text{ , P1 elements}$

For $Est.3$ the data shows similar behavior to $Est.2$ except for fixed
$\varepsilon$ and $h\downarrow 0$. In this case, ill-conditioning for small
$\varepsilon$ and $h\downarrow 0$ is over-estimated compared with Table 1.
This is expected because the data comes from P1 elements and $||\nabla\cdot
u^{\varepsilon}||$ occurs in $Est.3$. Again, $Est.3$ still provides a smaller
estimate of ill-conditioning than $Est.1$.

We have attributed the ill-conditioning observed above as
$\varepsilon\rightarrow 0$ to the use of P1 elements. To test if this
explanation is plausible we repeated this test with P2 elements. Since this
finite element space contains a non-trivial divergence free subspace, [15],
our intuition is that ill-conditioning would be reduced. Table 3 below
presents $Est.3\simeq\varepsilon^{-1}h^{-1}\frac{||\nabla\cdot
u^{\varepsilon}||}{||u^{\varepsilon}||}$ values.
$Est.2\simeq||f^{h}||/||u^{\varepsilon}||$ values, not presented, were
significantly smaller than $Est.3$ and converged to $1.13$ as
$\varepsilon,h\rightarrow 0$.

$\displaystyle\begin{array}[]{c|cccccc}{\small\varepsilon\Downarrow\setminus
h\Rightarrow}&{\small 0.125}&{\small 6.2E-2}&{\small 3.1E-2}&\text{ }{\small
1.6E-2}&{\small 7.8E-3}&{\small 3.9E-3}\\\ \hline\cr{\small 1}&{\small
25}&{\small 50}&{\small 1.0E2}&{\small 2.0E2}&{\small 4.0E2}&{\small 8.0E2}\\\
{\small 6.3E-2}&{\small 3.1E2}&{\small 6.2E2}&{\small 1.2E3}&{\small
2.5E3}&{\small 5.0E3}&{\small 1.0E4}\\\ {\small 3.9E-3}&{\small 6.7E2}&{\small
1.3E3}&{\small 2.5E3}&{\small 5.0E3}&{\small 1.0E4}&{\small 2.0E4}\\\ {\small
2.4E-4}&{\small 7.1E2}&{\small 1.3E3}&{\small 2.6E3}&{\small 5.1E3}&{\small
1.0E4}&{\small 2.0E4}\\\ {\small 1.5E-5}&{\small 7.1E2}&{\small 1.3E3}&{\small
2.6E3}&{\small 5.2E3}&{\small 1.0E4}&{\small 2.0E4}\\\ {\small 9.5E-7}&{\small
7.1E2}&{\small 1.3E3}&{\small 2.7E3}&{\small 5.3E3}&{\small 1.1E4}&{\small
2.1E4}\\\ {\small 6.0E-8}&{\small 7.1E2}&{\small 1.3E3}&{\small 2.7E3}&{\small
5.3E3}&{\small 1.1E4}&{\small 2.1E4}\end{array}$ $\displaystyle\text{Table 3:
Values of }Est.3\simeq\varepsilon^{-1}h^{-1}\frac{||\nabla\cdot
u^{\varepsilon}||}{||u^{\varepsilon}||}\text{, P2 elements}$

The above data indicates that with P2 elements and a problem with a smooth
solution, the contribution of the penalty term to ill-conditioning is small.

Test problem 2: We take

$f(x,y)=\left(10+50\sin\left(\frac{2\pi x}{h}+\frac{2\pi
y}{h}\right),10+50\sin\left(\frac{2\pi x}{h}+\frac{2\pi
y}{h}\right)\right)^{T}.$

Each component oscillates as rapidly as the given mesh allows, see Figure 1
below for $h=0.125$.

Figure 1: Plot of $10+50\sin(\frac{2\pi x}{0.125}+\frac{2\pi y}{0.125})$.

The forcing function and thus the solution change with the mesh size.

$\displaystyle\begin{array}[]{c|cccccc}{\small\varepsilon\Downarrow\setminus
h\Rightarrow}&{\small 0.125}&{\small 6.2E-2}&{\small 3.1E-2}&\text{ }{\small
1.6E-2}&{\small 7.8E-3}&{\small 3.9E-3}\\\ \hline\cr{\small 1}&{\small
1.1E2}&{\small 1.0E2}&{\small 1.0E2}&{\small 1.0E2}&{\small 1.0E2}&{\small
1.0E2}\\\ {\small 6.3E-2}&{\small 5.6E2}&{\small 5.6E2}&{\small 5.6E2}&{\small
5.6E2}&{\small 5.6E2}&{\small 5.6E2}\\\ {\small 3.9E-3}&{\small 7.4E3}&{\small
7.5E3}&{\small 7.6E3}&{\small 7.6E3}&{\small 7.7E3}&{\small 7.7E3}\\\ {\small
2.4E-4}&{\small 1.1E5}&{\small 1.1E5}&{\small 1.2E5}&{\small 1.2E5}&{\small
1.2E5}&{\small 1.2E5}\\\ {\small 1.5E-5}&{\small 1.8E6}&{\small 1.8E6}&{\small
1.8E6}&{\small 1.8E6}&{\small 1.8E6}&{\small 1.9E6}\\\ {\small 9.5E-7}&{\small
2.9E7}&{\small 2.9E7}&{\small 2.9E7}&{\small 2.9E7}&{\small 2.9E7}&{\small
2.9E7}\\\ {\small 6.0E-8}&{\small 4.7E8}&{\small 4.6E8}&{\small 4.6E8}&{\small
4.6E8}&{\small 4.6E8}&{\small 4.6E8}\end{array}$ $\displaystyle\text{Table 3:
Values of }Est.2\simeq\frac{||f^{h}||}{||u^{\varepsilon}||}\text{, P1
elements}$

For comparison with the last row in Table 3, the values of $Est.1\simeq
h^{-2}+\varepsilon^{-1}h^{-2}$ for $\varepsilon=6.0E-8$ are

$\begin{array}[]{cccccc}{\small 1.1E9}&{\small 4.3E9}&{\small 1.7E10}&{\small
6.9E10}&{\small 2.7E11}&{\small 1.1E12}\end{array}.$

In Table 3, as $\varepsilon\downarrow 0,$ $Est.2$ grows roughly like
$\varepsilon^{-1}$. The stable behavior of $Est.2$ as $h\downarrow 0$ is
unexpected. The behavior of
$Est.3\simeq\varepsilon^{-1}h^{-1}\frac{||\nabla\cdot
u^{\varepsilon}||}{||u^{\varepsilon}||}$ for P1 elements was similar, so not
presented.

We now present tests of P2 elements for this problem. In this test $Est.3$ was
larger than $Est.2,$ so we present only $Est.3$ data.

$\displaystyle\begin{array}[]{c|cccccc}{\small\varepsilon\Downarrow\setminus
h\Rightarrow}&{\small 0.125}&{\small 6.2E-2}&{\small 3.1E-2}&\text{ }{\small
1.6E-2}&{\small 7.8E-3}&{\small 3.9E-3}\\\ \hline\cr{\small 1}&{\small
25}&{\small 49}&{\small 1.0E2}&{\small 2.0E2}&{\small 4.0E2}&{\small 7.9E2}\\\
{\small 6.3E-2}&{\small 3.8E2}&{\small 7.5E2}&{\small 1.5E3}&{\small
3.0E3}&{\small 6.0E3}&{\small 1.2E4}\\\ {\small 3.9E-3}&{\small 6.0E3}&{\small
1.2E4}&{\small 2.4E4}&{\small 4.8E4}&{\small 9.6E4}&{\small 1.9E5}\\\ {\small
2.4E-4}&{\small 9.6E4}&{\small 1.9E5}&{\small 3.8E5}&{\small 7.7E5}&{\small
1.5E6}&{\small 3.1E6}\\\ {\small 1.5E-5}&{\small 1.5E6}&{\small 3.1E6}&{\small
6.2E6}&{\small 1.2E7}&{\small 2.5E7}&{\small 4.9E7}\\\ {\small 9.5E-7}&{\small
2.5E7}&{\small 4.9E7}&{\small 9.8E7}&{\small 2.0E8}&{\small 3.9E8}&{\small
7.9E8}\\\ {\small 6.0E-8}&{\small 3.9E8}&{\small 7.9E8}&{\small 1.6E9}&{\small
3.2E9}&{\small 6.3E9}&{\small 1.3E10}\end{array}$ $\displaystyle\text{Table 4:
Values of }Est.3\simeq\varepsilon^{-1}h^{-1}\frac{||\nabla\cdot
u^{\varepsilon}||}{||u^{\varepsilon}||}\text{ },\text{ P2 elements}$

For comparison with the last row in Table 4, $Est.1$ values for
$\varepsilon=6.0E-8$ are

$\begin{array}[]{cccccc}{\small 1.1E9}&{\small 4.3E9}&{\small 1.7E10}&{\small
6.9E10}&{\small 2.7E11}&{\small 1.1E12}\end{array}.$

For this problem, down the columns (fixed $h$, $\varepsilon\downarrow 0$ ) the
computed estimate of $\kappa(u^{\varepsilon})$ grows roughly like
$\varepsilon^{-1}$. Reading across the columns, the values in each row grow
like $h^{-1}$ (not $h^{-2}$). In all cases, $Est.3$ was smaller than $Est.1$.

##  4 Conclusions

The classical view is that two significant issues are impediments to penalty
methods giving low cost and highly accurate velocity approximates. The first
is ill-conditioning of the resulting system matrix. We have shown that the
effective condition number $\kappa(u^{\varepsilon})$ is much smaller than the
usual condition number due to the magnitude of the components in the penalized
eigenspaces being small in a precise sense. Motivated by this theoretical
result, we then compared the derived estimates of ill-conditioning on two test
problems and for two elements. With P1 elements, ill-conditioning was not as
severe as $\mathcal{O}(\varepsilon^{-1}h^{-2})$ but followed
$\varepsilon^{-1}$ growth as $\varepsilon\rightarrow 0$. For P2 elements and
approximating a smooth solution, $\kappa(u^{\varepsilon})$ was smaller the
expected $\kappa\sim h^{-2}$ for the discrete Laplacian. With both P1 and P2
elements, and an academic test problem with data oscillating as fast as the
mesh allows, ill-conditioning was not as severe as the expected
$\mathcal{O}(\varepsilon^{-1}h^{-2})$ but also followed the $\varepsilon^{-1}$
pattern as $\varepsilon\rightarrow 0$.

The second significant issue is the difficulty in the selection of an
effective value of $\varepsilon$. While the most commonly recommended, e.g.
[8, 14], choices are $\varepsilon=$ time step, mesh width, and (machine
epsilon)1/2, none have proven reliably effective. Recent work of [17, 25] may
resolve this impediment by an algorithmic, self-adaptive selection of
$\varepsilon$ based on some indicator of violation of incompressibility. The
estimates in Theorem 2.1 give insight into the resulting conditioning when
this is done. Let $TOL$ denote a specified tolerance. If $\varepsilon$ is
adapted so that the penalized solution $u^{\varepsilon}$ satisfies
$\frac{||\nabla\cdot u^{\varepsilon}||}{||u^{\varepsilon}||}\leq TOL,$ then,
Theorem 2.1 immediately implies $\kappa(u^{\varepsilon})$ satisfies

(2) $\kappa(u^{\varepsilon})\leq
Ch^{-2}\left(1+h\frac{TOL}{\varepsilon}\right).$

If $\varepsilon$ is adapted so that the penalized solution $u^{\varepsilon}$
satisfies $\frac{||\nabla\cdot u^{\varepsilon}||}{||\nabla
u^{\varepsilon}||}\leq TOL,$ then, similarly, Theorem 2.1 implies
$\kappa(u^{\varepsilon})$ satisfies

(3) $\kappa(u^{\varepsilon})\leq
Ch^{-2}\left(1+\frac{TOL}{\varepsilon}\right).$

For (3), rewrite $||\nabla\cdot u^{\varepsilon}||/||u^{\varepsilon}||$ as
$\left(||\nabla\cdot u^{\varepsilon}||/||\nabla
u^{\varepsilon}||\right)\left(||\nabla
u^{\varepsilon}||/||u^{\varepsilon}||\right)$. The inverse estimate, A1,
implies $||\nabla u^{\varepsilon}||/||u^{\varepsilon}||\leq Ch^{-1}$, yielding
(3).

The numerical tests suggest that two factors not considered in Theorem 2.1 are
significant. The first is whether the finite element space has a nontrivial,
divergence-free subspace. The second is the influence of smoothness of the
sought solution or its problem data on effective conditioning. In addition,
highly refined meshes are used in practical flow simulations and the linear
system often has large skew symmetric part. The extension of the analysis
herein to include these effects is an open problem.

## References

  * [1] A.O.H. Axelsson, _Condition Numbers for the Study of the Rate of Convergence of the Conjugate gradient Method_ , p. 3-33 in: Proc. of the 2nd IMACS Internat. Symposium on Iterative Methods in Linear Algebra, Blagoevgrad, Bulgaria, June 17-20, 1995.
  * [2] A.O.H. Axelsson,_Iterative solution methods_ , Cambridge University Press, 1996..
  * [3] A.O.H. Axelsson and I.E. Kaporin,_On the sublinear and superlinear rate of convergence of conjugate gradient methods_. Numerical Algorithms, 25 (2000) 1-22.
  * [4] A.O.H. Axelsson and I.E. Kaporin, _Error norm estimation and stopping criteria in preconditioned conjugate gradient iterations_. Numerical Linear Algebra with Applications, 8 (2001) 265–286.
  * [5] J.M. Banoczi, Nan-Chieh Chiu, Grace E. Cho and I.C.F. Ipsen, _The lack of influence of right-hand side on accuracy of linear system solution_ , SIAM Journal on Scientific Computing, 20 (1998) 203–227.
  * [6] S.C. Brenner, and L.R. Scott, _The mathematical theory of finite element methods_ , New York: Springer, 2008.
  * [7] C. Brezinski, _Error estimates for the solution of linear systems_ , SIAM Journal on Scientific Computing, 21 (1999) 764–781.
  * [8] G.F. Carey and R. Krishnan, _Penalty approximation of Stokes flow_ , Comp. Meth. Appl. Mech. Eng., 35 (1982) 169–206.
  * [9] T.F. Chan and D.E. Foulser,_Effectively well-conditioned linear systems_. SIAM Journal on Scientific and Statistical Computing, 9 (1988) 963–969.
  * [10] P. G. Ciarlet, _The finite element method for elliptic problems_ , vol. 40 of Classics in Applied Mathematics. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2002. (Reprint of the 1978 original).
  * [11] S. Christiansen and P.C. Hansen, _The effective condition number applied to error analysis of certain boundary collocation methods,_ Journal of Computational and Applied Mathematics, 54(1994) 15–36.
  * [12] S.G. Demko, _Condition numbers of rectangular systems and bounds for generalized inverses_ , Linear Algebra and its Applications, 78 (1986) 199-206.
  * [13] Qiang Du, Desheng Wang and Liyong Zhu, _On Mesh Geometry and Stiffness Matrix Conditioning for General Finite Element Spaces_ , SIAM Journal on Numerical Analysis, 47 (2009) 1421-1444.
  * [14] J.C. Heinrich and C.A. Vionnet, _The penalty method for the Navier-Stokes equations_. Archives of Computational Methods in Engineering, 2 (1995) 51-65.
  * [15] V. John, A. Linke, C. Merdon, M. Neilan, and L.G. Rebholz,_On the divergence constraint in mixed finite element methods for incompressible flows_. SIAM Review, 59 (2017), 492-544.
  * [16] I.E. Kaporin, _New convergence results and preconditioning strategies for the conjugate gradient method_. Numerical Linear Algebra with Applications, 2 (1994) 179-210.
  * [17] K. Kean, Xihui Xie and Shuxian Xu, _A Doubly Adaptive Penalty Method for the Navier Stokes Equation_ , preprint, 2021, https://arxiv.org/abs/2201.03978.
  * [18] Z.C. Li and H.T. Huang, _Effective condition number for numerical partial differential equations_. Numerical Linear Algebra with Applications, 7 (2008) 575-94.
  * [19] Z.C. Li, H.T. Huang and J. Huang, _Superconvergence and stability for boundary penalty techniques of finite difference methods_. Numerical Methods for Partial Differential Equations, 24 (2008) 972-90.
  * [20] J. von Neumann and H.H. Goldstine, _Numerical inverting of matrices of high order_ , Bulletin of the American Mathematical Society, 53 (1947) 1021-1099.
  * [21] J.R. Pragua, _New condition number for matrices and linear systems_ , Computing 41 (1989) 211–213.
  * [22] T. Tao and V. Vu, _Smooth analysis of the condition number and the least singular value_ , Mathematics of Computation, 79 (2010) 2333-2352.
  * [23] A.M. Turing, _Rounding-off errors in matrix processes_ , Quarterly Journal of Mechanics and Applied Mathematics, 1 (1948): 287-308.
  * [24] J.H. Wilkinson, _Rounding Errors in Algebraic Processes_ , Volume 32 of National Physical Laboratory Notes on Applied Science, (HMSO), London, 1963.
  * [25] Xihui Xie, _On Adaptive Grad-Div Parameter Selection_ , technical report, https://doi.org/10.48550/arXiv.2108.01766, 2021\.


# TenSEAL: A Library for Encrypted Tensor Operations Using Homomorphic
Encryption

Ayoub Benaissa
École Supérieure en Informatique, Sidi Bel Abbès / OpenMined
Algeria
<EMAIL_ADDRESS>
Bilal Retiat
École Supérieure en Informatique, Sidi Bel Abbès / OpenMined
Algeria
<EMAIL_ADDRESS>
Bogdan Cebere
OpenMined
<EMAIL_ADDRESS>
Alaa Eddine Belfedhal
École Supérieure en Informatique, Sidi Bel Abbès
Algeria
<EMAIL_ADDRESS>
###### Abstract

Machine learning algorithms have achieved remarkable results and are widely
applied in a variety of domains. These algorithms often rely on sensitive and
private data such as medical and financial records. Therefore, it is vital to
draw further attention regarding privacy threats and corresponding defensive
techniques applied to machine learning models. In this paper, we present
TenSEAL, an open-source library for Privacy-Preserving Machine Learning using
Homomorphic Encryption that can be easily integrated within popular machine
learning frameworks. We benchmark our implementation using MNIST and show that
an encrypted convolutional neural network can be evaluated in less than a
second, using less than half a megabyte of communication.

##  1 Introduction

In recent years, we have witnessed the evolution of machine learning as a
service (MLaaS). In a typical scenario, the users need to send their input to
the service provider, which will execute some algorithms on the data and send
back the result.

This new way of making inferences has two critical issues. First, the users
may not want to send their data to the service provider due to privacy
concerns. Second, if we do not send users’ data to the service provider, we
cannot give the users the model due to intellectual property concerns. Using
homomorphic encryption, we can follow the same method, except that users’ data
will always be encrypted. This way, neither the input nor the output will be
visible to the service provider, and the evaluation can still happen on this
encrypted data.

However, the adoption of homomorphic encryption in machine learning is slow.
One reason is that while the available libraries provide an excellent API for
cryptographers, they might be challenging to use for data scientists. The
other blocker is also the cost for evaluation, both in terms of communication
and computation.

###  1.1 Contributions

  * •

We present a flexible open-source library for doing encrypted tensor
computation using homomorphic encryption. The library can directly convert
tensors from popular machine learning frameworks (like PyTorch or Tensorflow)
to their encrypted versions.

  * •

We evaluate a convolution neural network on encrypted data in less than a
second, with less than half a megabyte of communication during inference.

For the rest of the paper, we describe the library’s architecture in Section
2. Then, in Section 3, we detail the algorithms needed for evaluating a
convolutional neural network in the encrypted space. In Section 5, we provide
an experimental evaluation of our library, and conclude with some limitations
of our work in section 6.

##  2 Architecture

TenSEAL  is a library that bridges classical machine learning frameworks to
homomorphic encryption capabilities. It manages all the complexities of
implementing tensor operations on encrypted data. TenSEAL relies on the
implementation of the CKKS (Cheon et al. (2017)) scheme in Microsoft SEAL .
The clients can work with plain or encrypted tensors using one of the
supported frontend languages (C++ or Python). In a client-server scenario, the
message exchange is done using Protocol buffers . The core API is built around
three main components: the context, the plain tensors, and the encrypted
tensors.

###  2.1 The TenSEAL context

The TenSEAL context is the central component of the library. It generates and
stores the necessary keys required by an encrypted computation. The context
generates the secret-key used for decryption, the public-key used for
encryption, the Galois-keys used for rotation, and the relinearization-keys
used for relinearization of ciphertexts. This same object will also handle the
thread-pool, which controls how many jobs should be run in parallel when
performing parallelizable operations. The context can also be configured to do
automatic ciphertext relinearization and rescaling during computation.

###  2.2 The PlainTensor

The PlainTensor is a class that connects unencrypted tensors to the encrypted
implementations. Figure 2 in the appendix describes the process of converting
the tensors.

###  2.3 Encrypted tensors

The EncryptedTensor interface offers an API that needs to be implemented by
every tensor exposed by the library. The interface has a TenSEALContext
object, necessary to make any homomorphic encryption operation. The derived
classes expose different tensor flavors, such as:

  * •

CKKSVector derives the EncryptedVector interface and can hold a vector of real
values by encrypting them into a single ciphertext.

  * •

CKKSTensor follows the same strategy as Boemer et al. (2019b) and holds an
N-dimensional tensor of real values by encrypting them into N-dimensional
tensor of ciphertexts. However, it can batch an axis along with the slots
available in every ciphertext, thus requiring only an (N-1)-dimensional tensor
of ciphertexts.

Figures 2 and 2 in the appendix describe how an encrypted tensor is
constructed. From this point onward, we will focus on CKKSVector, as it will
be the type used to evaluate the MNIST dataset.

##  3 Method

When building a tensor over a Homomorphic Encryption scheme, there are two
significant concerns to tackle: 1.How to encode the tensor before encryption?
and 2.What operations can be performed when using a particular encoding? The
batching feature of the CKKS scheme allows a $N\times N$ matrix to be
encrypted into $N$ ciphertexts, with each row or column as a ciphertext.
Another possibility is to encrypt the whole tensor into a single ciphertext
(Jiang et al. (2018)). Depending on how we put our plain tensor into
ciphertexts, we can perform different operations with varying complexities.
The goal is to use the minimum number of ciphertexts and have a maximum depth
with a minimum runtime, thus optimizing memory and computation. Seeking this
ideal goal, we found out that we can use a single ciphertext to encrypt an
input image and evaluate it on a convolutional neural network. This requires a
pre-processing step on the client-side to encode the image as a matrix,
composed of convolution windows as rows, then flatten it as a vector via a
vertical scan. In TenSEAL, all these functionalities are implemented around
CKKSVector. A CKKSVector holds $N/2$ real values, where $N$ is the polynomial
modulus degree. We can perform element-wise operations with other encrypted or
plain vectors (addition, subtraction, and multiplication). We have a method
for computing the power of an encrypted vector (element-wise) that uses an
optimal circuit, thus using a minimum multiplicative depth. Also, because we
need a polynomial approximation for different activation functions, we built a
method for evaluating a polynomial with the encrypted vector as a variable,
making sure to use a minimal circuit. Apart from the element-wise operations,
we also need matrix operations to perform machine learning tasks. We
implemented a variant of the encrypted vector-plain matrix multiplication
proposed by Halevi & Shoup (2014) that can use multiple threads to run faster.
You can check Table 3 in the appendix for a list of supported operations by
the library’s encrypted tensors.

###  3.1 Dot Product

The library provides a similar algorithm for dot product as in Halevi & Shoup
(2014), but supports vectors of a size that is not a power of two or does not
fill all the slots of a ciphertext. This limitation in the previous method was
due to the right rotation that expects the final element to be the first,
which was not true with regard to the cases we addressed. Our method is
limited to a specific number of dot products if the vector size is not a power
of two. However, this limitation is generally not reached, as the number of
multiplication allowed by the scheme might be lower. We do this by replicating
the input vector as many times as possible into ciphertext slots and only left
ciphertext rotations during computation. Our method has the same algorithmic
complexity as in Halevi & Shoup (2014). We implemented it using CKKS (Cheon et
al. (2017)), for the dot product operation between an encrypted vector and a
plain matrix. Thus, it can be extended to support a dot product between an
encrypted matrix with a plain matrix (matrix multiplication). Figure 3 in the
appendix shows how to perform a dot product between an encrypted vector and a
plain matrix.

###  3.2 2-D Convolution

TenSEAL also supports evaluating convolutions, with a similar implementation
of how modern machine learning frameworks (e.g., PyTorch) are computing them.
We applied the Image Block to Columns (im2col) (Johnson et al. (2016))
technique, which turns a convolution layer into a single matrix multiplication
operation. This technique requires encrypted matrix-plain vector
multiplication, which we implemented by performing element-wise multiplication
of the matrix transpose with replicated plain vector. Finally, it rotates and
accumulates the result into a single vector. This operation uses only one
multiplication operation and $log_{2}(N)$ rotations and additions, where $N$
represents the rows’ number in the matrix. Section A.4 in the appendix
explains how the ”image block to column” algorithm can be applied to encrypted
inputs. It is important to note that the transformation happens in plain data,
and the transformed input image will be encoded and encrypted in a single
ciphertext. This directly implies that stacking two convolutions is not
possible, as reorganizing the slots of a ciphertext is not trivial.

##  4 Related work

In recent years, several research works have made homomorphic encryption
schemes practical for machine learning. Gilad-Bachrach et al. (2016)
implemented CryptoNets, a neural network for making inference on encrypted
data using the YASHE (Bos et al. (2013)) leveled homomorphic encryption
scheme, which has efficient plain addition and multiplication algorithms,
useful for unencrypted models. However, the framework requires large batches
for achieving a good amortized performance, making it less practical for use
cases that evaluate a single instance. Boemer et al. (2019b) used a similar
tensor structure as CryptoNets while implementing different optimization
layers. They used graph level optimizations specific for HE applications,
which reduced the multiplicative depth needed for evaluating operations such
as batch-norm and average pooling. In subsequent work, Boemer et al. (2019a)
evaluated MobileNetV2 (Sandler et al. (2018)) and reported empirical results
on HE encrypted inputs, which was way deeper than the previously used models.
Jiang et al. (2018) used a smaller neural network and only convolution, linear
layers, and the square activation function. Their framework E2DM made
predictions on encrypted data using the CKKS (Cheon et al. (2017)) scheme, but
compared to previous works, the model’s parameters were also encrypted.
Juvekar et al. (2018) proposed the Gazelle framework, which mixes HE with
Garbled Circuits (GC) Yao (1986). They switched between both methods during
computation, choosing the most efficient at a certain point based on the next
operation. Using GC makes it possible to compute the ReLu activation function,
compared to previous works (Gilad-Bachrach et al. (2016); Jiang et al. (2018);
Boemer et al. (2019b; a)) which have mainly used polynomial functions. Even
though the protocol achieves relatively fast run-time, it requires interaction
between participants, resulting in high bandwidth usage. Badawi et al. (2020)
used a GPU-accelerated implementation of BFV(Fan & Vercauteren (2012)) based
on the work from Badawi et al. (2018). Their HCNN of 5 layers could evaluate
in 5 seconds, but could batch more than 8000 images without extra overhead.

All the works are benchmarked using the MNIST dataset (LeCun et al. (2010)),
but with different hardware configurations. We summarize empirical results
reported in each of the corresponding papers in Table 1.

Framework | Method | Batch size | Message Size | Evaluation Time | Accuracy
---|---|---|---|---|---
CryptoNets | HE | 8192 | 595.5 MB | 570 s | 99%
Gazelle | HE, MPC | 1 | 0.5 MB | 0.03 s | -
E2DM | HE | 64 | 23.93 MB | 1.69 s | 98.1%
HCNN-GPU | HE | 8192 | - | 5.16 s | 99%
Table 1: Comparing frameworks and their evaluation results on MNIST.

##  5 Evaluation

To evaluate our library and technique, we implemented a neural network
composed of: a convolutional layer (4 kernels of 7x7, with a stride of 3x3), a
linear layer (input: 256, output: 64), and a final linear layer (input: 64,
output: 10). We used the square activation function after every layer except
for the last. The convolution was done using our image to column
implementation, while the linear layers use the dot product implementation.
The accuracy on the plain test-set was 97.7% in contrast to 97.4% for the
encrypted test-set. We used the CKKSVector implementation, which uses the CKKS
scheme. Knowing that we need 6 multiplications to perform the evaluation and a
security level of 128-bits, we set the polynomial modulus degree to 8192, with
a coefficient modulus of 206-bits, and a scale of 21-bits. The evaluation was
done on Ubuntu Server 20.04 and Python 3.8, using AWS c4.2xlarge (8 vCPUs) and
AWS c4.4xlarge (16 vCPUs) configurations. The measured durations are the
average of 5 rounds of testing, with 10 iterations each. Table 2 contains a
full breakdown for evaluating the neural network over encrypted images sampled
from the MNIST dataset.

Operation | Description | Duration(ms)
---|---|---
AWS c4.2xlarge(8 vCPUs) | AWS c4.4xlarge(16 vCPUs)
Key generation | Generate the context and the encryption keys | 940.01 | 921.04
Input preparation | im2col encoding | 9.8 | 9.8
Convolutional layer evaluation | Input $28\times 28$, kernel $7\times 7,$ stride 3, 4 channels | 236.9 | 237.98
First activation(square) | Square $256$ input values | 8.47 | 8.42
FC1 | Fully connected layer with $256$ inputs and $64$ outputs | 1084.65 | 575.34
Second activation(square) | Square $64$ input values | 4.29 | 4.2
FC2 | Fully connected layer with $64$ inputs and $10$ outputs | 121.36 | 70.36
Full forward step | All the steps above | 1456.29 | 887.06

Table 2: A complete illustration of the encrypted MNIST evaluation, with
durations expressed in milliseconds. We evaluate the methods using two setups:
Amazon c4.2xlarge instance (8 vCPUs, 15 GiB memory) and Amazon c4.4xlarge
instance (16 vCPUs, 30 GiB memory), to underline how the library takes
advantage of the available parallelism.

The results show that the library makes heavy use of the available
parallelism, and it is highly competitive in terms of network communication,
requiring only 427KB of communication to send the encrypted input and receive
the encrypted output. At the same time, TenSEAL does not enforce a specific
batch size for the inference, making it quite practical. The complete
operations benchmarks are open-source, and the results are included in Section
A.5 in the appendix. Table 4 shows the average performance for different
arithmetic operations. Table 5 shows the average performance for matrix
multiplications.

##  6 Limitations and conclusion

The encryption part relies on CKKS Cheon et al. (2017), which is known to be a
leveled homomorphic encryption scheme. This means that depending on our
parameter selection, there is a limit on how many multiplications we can
perform on encrypted data, and this directly impacts the machine learning
model we can use or its depth. Different machine learning models also use non-
linear activation functions, which will need to be approximated using
polynomials in the case of CKKS. A recent work Chillotti et al. (2020b) have
been trying to solve this issue tightly related to machine learning by using
the TFHE scheme Chillotti et al. (2020a), which allows the evaluation of
deeper models, as well as non-linear activation functions.

In conclusion, our results show that it can be practical to do tensorial
operations using the CKKS scheme. Depending on the use case, users can choose
advanced tensor operations (like slicing or broadcasting) or use more
computation-communication optimized implementations. TenSEAL can accommodate
both scenarios while offering a smooth transition from the traditional machine
learning frameworks. Finally, we seek to extend the tensor operations catalog
and to improve the overall performance even further.

## References

  * Badawi et al. (2018) Ahmad Al Badawi, Bharadwaj Veeravalli, Chan Fook Mun, and Khin Mi Mi Aung.  High-performance fv somewhat homomorphic encryption on gpus: An implementation using cuda.  In _IACR Transactions on Cryptographic Hardware and Embedded Systems_ , pp. 70–95, 2018.
  * Badawi et al. (2020) Ahmad Al Badawi, Jin Chao, Jie Lin, Chan Fook Mun, Jun Jie Sim, Benjamin Hong Meng Tan, Xiao Nan, Khin Mi Mi Aung, and Vijay Ramaseshan Chandrasekhar.  Towards the alexnet moment for homomorphic encryption: Hcnn, the first homomorphic cnn on encrypted data with gpus.  In _IEEE Transactions on Emerging Topics in Computing_ , 2020.  doi: 10.1109/TETC.2020.3014636.
  * Boemer et al. (2019a) Fabian Boemer, Anamaria Costache, Rosario Cammarota, and Casimir Wierzynski.  ngraph-he2: A high-throughput framework for neural network inference on encrypted data.  In _Proceedings of the 7th ACM Workshop on Encrypted Computing & Applied Homomorphic Cryptography_, pp. 45–56, 2019a.
  * Boemer et al. (2019b) Fabian Boemer, Yixing Lao, Rosario Cammarota, and Casimir Wierzynski.  ngraph-he: a graph compiler for deep learning on homomorphically encrypted data.  In _Proceedings of the 16th ACM International Conference on Computing Frontiers_ , pp. 3–13, 2019b.
  * Bos et al. (2013) Joppe W Bos, Kristin Lauter, Jake Loftus, and Michael Naehrig.  Improved security for a ring-based fully homomorphic encryption scheme.  In _IMA International Conference on Cryptography and Coding_ , pp. 45–64. Springer, 2013.
  * Cheon et al. (2017) Jung Hee Cheon, Andrey Kim, Miran Kim, and Yongsoo Song.  Homomorphic encryption for arithmetic of approximate numbers.  In _International Conference on the Theory and Application of Cryptology and Information Security_ , pp. 409–437. Springer, 2017.
  * Chillotti et al. (2020a) Ilaria Chillotti, Nicolas Gama, Mariya Georgieva, and Malika Izabachène.  Tfhe: fast fully homomorphic encryption over the torus.  _Journal of Cryptology_ , 33(1):34–91, 2020a.
  * Chillotti et al. (2020b) Ilaria Chillotti, Marc Joye, and Pascal Paillier.  New challenges for fully homomorphic encryption.  _Privacy-preserving Machine Learning (PPML-PriML 2020) NeurIPS 2020 workshop_ , 2020b.
  * Fan & Vercauteren (2012) Junfeng Fan and Frederik Vercauteren.  Somewhat practical fully homomorphic encryption.  Cryptology ePrint Archive, Report 2012/144, 2012.  https://eprint.iacr.org/2012/144.
  * Gilad-Bachrach et al. (2016) Ran Gilad-Bachrach, Nathan Dowlin, Kim Laine, Kristin Lauter, Michael Naehrig, and John Wernsing.  Cryptonets: Applying neural networks to encrypted data with high throughput and accuracy.  In _International Conference on Machine Learning_ , pp. 201–210, 2016.
  * Halevi & Shoup (2014) Shai Halevi and Victor Shoup.  Algorithms in helib.  pp. 554–571, 2014.
  * Jiang et al. (2018) Xiaoqian Jiang, Miran Kim, Kristin Lauter, and Yongsoo Song.  Secure outsourced matrix computation and application to neural networks.  In _Proceedings of the 2018 ACM SIGSAC Conference on Computer and Communications Security_ , pp. 1209–1222, 2018.
  * Johnson et al. (2016) Justin Johnson, Fei-Fei Li, and Andrej Karpathy.  Cnns in practice. convolutional neural networks for visual recognitio.  2016\.
  * Juvekar et al. (2018) Chiraag Juvekar, Vinod Vaikuntanathan, and Anantha Chandrakasan.  $\\{$GAZELLE$\\}$: A low latency framework for secure neural network inference.  In _27th $\\{$USENIX$\\}$ Security Symposium ($\\{$USENIX$\\}$ Security 18)_, pp. 1651–1669, 2018.
  * LeCun et al. (2010) Yann LeCun, Corinna Cortes, and CJ Burges.  Mnist handwritten digit database.  _ATT Labs [Online]. Available: http://yann.lecun.com/exdb/mnist_ , 2, 2010.
  * (16) Microsoft SEAL.  Microsoft SEAL (release 3.6).  https://github.com/Microsoft/SEAL, November 2020.  Microsoft Research, Redmond, WA.
  * (17) Protocol buffers.  Protocol buffers – Google’s data interchange format.  URL https://github.com/protocolbuffers/protobuf.
  * Sandler et al. (2018) Mark Sandler, Andrew Howard, Menglong Zhu, Andrey Zhmoginov, and Liang-Chieh Chen.  Mobilenetv2: Inverted residuals and linear bottlenecks.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pp. 4510–4520, 2018.
  * (19) TenSEAL.  TenSEAL (release 0.3.0).  https://github.com/OpenMined/TenSEAL, February 2021.
  * Yao (1986) Andrew Chi-Chih Yao.  How to generate and exchange secrets.  In _27th Annual Symposium on Foundations of Computer Science (sfcs 1986)_ , pp. 162–167. IEEE, 1986.

##  Appendix A Appendix

###  A.1 Encrypted Tensor classes

Figure 1: A high-level overview of the encrypted tensor construction. The
PlainTensor wraps the tensor representation from popular frameworks, and it is
used as input for the EncryptedTensor interface.

Figure 2: Encrypted tensors relation. The EncryptedTensor interface is derived
into BFVVector, CKKVector, or CKKSTensor classes.

###  A.2 Encrypted Tensor Operations

Operation | Description
---|---
negate | Negate an encrypted tensor
square | Compute the square of an encrypted tensor
power | Compute the power of an encrypted tensor
add | Addition between an encrypted tensor and an encrypted/plain tensor
sub | Subtraction between an encrypted tensor and an encrypted/plain tensor
mul | Multiplication between an encrypted tensor and an encrypted/plain tensor
dot_product | Dot product between an encrypted tensor and an encrypted/plain tensor
polyval | Polynomial evaluation with an encrypted tensor as variable
matmul_plain | Multiplication between an encrypted tensor and an encrypted/plain matrix
conv2d_im2col | Image Block to Columns
Table 3: Supported operations for encrypted tensors

###  A.3 Dot Product

Figure 3 shows how an encrypted vector (in gray) can be multiplied with a
plain matrix using the method from Halevi & Shoup (2014).

Figure 3: Vector-Matrix Multiplication

###  A.4 2D Convolution

A 2D convolution can be performed using a single matrix multiplication,
instead of repeating multiplication on every window. This method is referred
to as image block to column convolution, or image to column convolution.
Figure 4 shows how a convolution can be performed using this method. It first
reorganizes the input matrix into rows representing convolution windows, then
performs a dot product with the flattened kernel.

Figure 4: Image to column convolution

Applying this technique to an encrypted matrix, which is encrypted into a
single ciphertext, is not trivial, as reorganizing slots is not that simple.
We will need to reorganize the matrix as a pre-processing step before
encryption to be ready for convolution. The encrypted-matrix (input image)
with plain-vector (kernel) can be performed with a single element-wise
multiplication and a series of rotations and accumulations. Figure 5 and 6
show the steps for doing that. The first shows how the encrypted matrix
(colored) is encoded and multiplied with the plain kernel. The second step is
to sum different versions of the output that are rotated differently to the
left.

Figure 5: Image to column convolution with CKKS - step 1 Figure 6: Image to
column convolution with CKKS - step 2

###  A.5 Complete benchmarks

In this section, we present the full evaluation of TenSEAL’s operations, ran
on the CKKSVector implementation. The benchmarks are measured on an Amazon EC2
c4.2xlarge instance, with 8 vCPUs at 2.9 GHz(Intel Xeon E5-2666 v3 Processor)
and 15 GiB memory. They are executed using Ubuntu Server 20.04 and Python 3.8.
The measured durations are the average of 5 rounds of testing, with 10
iterations each.

Tables 4 and 5 show the average performance for different arithmetic
operations.

Operation | Tensor shape
---|---
[256] | [1024] | [4096] | [8192] | [16384]
negate | 0.07 | 0.07 | 0.07 | 0.13 | 0.26
square | 4.29 | 4.29 | 4.29 | 8.49 | 17.16
polyval | 10.55 | 10.46 | 10.51 | 21.32 | 42.68
Table 4: Duration in milliseconds for unary operations. The CKKS context is created for polynomial modulus $8192$ and coefficient modulus of 200-bits. The polyval benchmark is executed for $2X^{2}+X$ polynoms. Operation | Tensor shape
---|---
[256] | [1024] | [4096] | [8192] | [16384]
add | 0.08 | 0.08 | 0.08 | 0.16 | 0.31
multiply | 4.45 | 4.34 | 4.43 | 8.84 | 17.75
sub | 0.08 | 0.08 | 0.08 | 0.15 | 0.3
dot | 20.15 | 23.96 | 28.11 | 55.94 | 112.36
add_plain | 0.8 | 0.86 | 1.07 | 2.13 | 4.19
multiply_plain | 1.75 | 1.81 | 2.03 | 4.02 | 7.97
sub_plain | 0.8 | 0.86 | 1.08 | 2.14 | 4.21
dot_plain | 17.37 | 21.36 | 25.63 | 51.14 | 101.82
Table 5: Duration in milliseconds for binary operations. The CKKS context is
created for polynomial modulus $8192$ and coefficient modulus of 200-bits. For
the ”_plain” operations, the operand is a PlainTensor of the same shape. For
the rest, the operand is an encrypted tensor of the same shape.


11institutetext:  Klinikum Rechts der Isar, Technical University of Munich,
Germany 22institutetext: Department of Computing, Imperial College London,
United Kingdom 33institutetext: Medical Image And Data Analysis (MIDAS.lab),
University Hospital of Tübingen, Germany

# Learning-based and unrolled motion-compensated reconstruction for cardiac MR
CINE imaging

Jiazhen Pan Daniel Rueckert Thomas Küstner Kerstin Hammernik11footnotemark: 1
contributed equally  111122331122

###### Abstract

Motion-compensated MR reconstruction (MCMR) is a powerful concept with
considerable potential, consisting of two coupled sub-problems: Motion
estimation, assuming a known image, and image reconstruction, assuming known
motion. In this work, we propose a learning-based self-supervised framework
for MCMR, to efficiently deal with non-rigid motion corruption in cardiac MR
imaging. Contrary to conventional MCMR methods in which the motion is
estimated prior to reconstruction and remains unchanged during the iterative
optimization process, we introduce a dynamic motion estimation process and
embed it into the unrolled optimization. We establish a cardiac motion
estimation network that leverages temporal information via a group-wise
registration approach, and carry out a joint optimization between the motion
estimation and reconstruction. Experiments on 40 acquired 2D cardiac MR CINE
datasets demonstrate that the proposed unrolled MCMR framework can reconstruct
high quality MR images at high acceleration rates where other state-of-the-art
methods fail. We also show that the joint optimization mechanism is mutually
beneficial for both sub-tasks, i.e., motion estimation and image
reconstruction, especially when the MR image is highly undersampled.

##  1 Introduction

Cardiac magnetic resonance imaging (CMR) plays an essential role in evidence-
based diagnostic of cardiovascular disease [16] and serves as the gold-
standard for assessment of cardiac morphology and function [18]. High-quality
cardiac image reconstruction with high spatial and temporal resolution is an
inevitable prerequisite for this assessment. Shorter scan times with higher
spatio-temporal resolution are desirable in this scenario. However, this
requires high acceleration rates which in turn is only achievable if
sufficient spatio-temporal information linked by the cardiac motion is shared
during reconstruction.

A large variety of CMR reconstruction methods have been dedicated to cope with
cardiac motion during the reconstruction. They can be categorized into two
sections: Implicit and explicit motion correction. Implicit motion correction
during the reconstruction sidesteps the non-rigid cardiac motion estimation,
which still remains one of the most challenging problems in CMR. Most research
in this section focused either on exploiting spatio-temporal redundancies in
complementary domains [14, 30], enforcing sparseness/low-rankness along these
dimensions [13, 19, 26], or improving the spatio-temporal image regularization
in an unrolled optimization model [17, 31]. Yet in this sense, motion is only
implicitly corrected without knowing or estimating the true underlying motion.

Motion can also be explicitly corrected during the reconstruction by applying
motion estimation/registration models. The work of Batchelor et al. [6] has
pioneered the field of motion-compensated MR reconstruction (MCMR). It
proposed the idea of embedding the motion information as an explicit general
matrix model into the MR reconstruction process and demonstrated some
successful applications in leg and brain reconstruction. However, in CMR, the
respiratory and cardiac motion is much more complex and therefore a more
sophisticated motion model is required. Some endeavours were made first in the
field of coronary magnetic resonance angiography (CMRA) to compensate
respiratory motion [7, 9] in which non-rigid motion models based on cubic
B-Splines were employed [15, 21]. However, these conventional registrations
require substantial computation times in the order of hours making the
practical implementation of image reconstruction infeasible. More recently,
learning-based registration methods [4, 5, 28] have been proposed, which
leverage trained neural networks to accelerate the estimation in inference. In
[29] these learned registrations are embedded into a CMRA reconstruction for
an unrolled CG-SENSE optimization [27] with a denoiser regularizer [1].
However, MCMR has been rarely studied in the context of cardiac CINE imaging
[8].

MCMR can be recast as two codependent and intertwined sub-optimization
problems: Image reconstruction and motion estimation. This stands in contrast
to methods in which only the MR image is optimized [9, 7, 29, 8], whereas the
motion is pre-computed and assumed to be fixed during the whole optimization
process. A reliable reconstruction relies on precise motion estimation, while
the accuracy of the motion prediction is subject to the quality of images.
Odille et al. [23, 24] introduced an iterative, alternating approach to solve
this problem. Concurrently, [3, 32] proposed to solve this joint optimization
problem by leveraging variational methods. However, to the best of our
knowledge, this problem has not been formulated and explored using deep
learning-based approaches. Furthermore, as algorithm unrolling has been
successfully applied in modern reconstruction approaches [1, 11] by learning
image regularization from data, the possibility to unroll the MCMR
optimization with motion estimation networks has not been explored.

In this work, we introduce a deep learning-based and unrolled MCMR framework
to reconstruct CMR images in the presence of motion. The highlights of the
work can be summarized as follows: First, we propose an unrolled MCMR
optimization process with an embedded group-wise motion estimation network
[12]. A joint optimization between the image reconstruction and motion
estimation is proposed, while the motion estimation is updated iteratively
with the reconstruction progress. Second, in contrast to [7, 8, 29] our
proposed approach does not require a pre-computed motion field, which is
usually obtained from initial non-motion compensated reconstruction. Finally,
instead of reconstructing only one frame by compensating motion from all
frames to one target frame, our proposed method provides a motion-resolved
image sequence in which all frames are motion-corrected. The proposed
framework was trained and tested on in-house acquired 2D CINE data of 40
subjects and compared to state-of-the-art MCMR methods. The conducted
experiments demonstrate that the proposed method outperforms the competing
methods especially in high acceleration situations, concluding that embedding
the dynamic motion into an unrolled MCMR can drastically improve the
reconstruction performance.

##  2 Methods

Assume that $x^{(t)}\in\mathbb{C}^{M}$ is the $t$-th frame of the image
sequence $x=[x^{(1)},\ldots,x^{(N)}]$ with $M$ pixels and
$y^{(t)}\in\mathbb{C}^{MQ}$ denotes its corresponding $t$-th undersampled
k-space data with $Q$ coils. In the context of CMR, the MCMR of the whole MR
sequence with $N$ temporal frames can be formulated as

$\displaystyle\min\limits_{x,\mathbf{U}}\sum_{t_{1}=1}^{N}\sum_{t_{2}=1}^{N}\left\|\mathbf{A}^{(t_{2})}\mathbf{U}^{(t_{1}\rightarrow
t_{2})}x^{(t_{1})}-y^{(t_{2})}\right\|^{2}_{2},$ (1)

where $\mathbf{A}^{(t)}$ denotes the MR multi-coil encoding operator with
$\mathbf{A}^{(t)}=\mathbf{D}^{(t)}\mathbf{F}\mathbf{C}$ for coil-sensitivity
maps $\mathbf{C}$, Fourier transform $\mathbf{F}$ and the undersampling matrix
$\mathbf{D}$. Based on the idea of [6], we build a motion matrix
$\mathbf{U}^{(t_{1}\rightarrow t_{2})}$ representing the warping from the
temporal frame $t_{1}$ to $t_{2}$, which is obtained from the estimated motion
fields $u^{(t_{1}\rightarrow t_{2})}$ (introduced in Section 2.2).

In traditional MCMR [7, 8, 29], only the images $x$ are optimized, whereas the
motion is pre-computed from initial (non-motion compensated) reconstructed
images just for once and assumed to be constant during the whole optimization
process. However, motion estimation on the initial motion-corrupted and
artifact-degraded images can be inaccurate, incurring error-propagation during
the optimization. Therefore, a joint optimization as proposed in [24] for
image $x$ and motion $\mathbf{U}$ as described in Eq. 1 is desired.

###  2.1 Motion-compensated image reconstruction framework

Figure 1: The proposed unrolled MCMR framework. A dynamic joint optimization
is performed between the image reconstruction in the CG-SENSE block and the
motion correction from the motion estimator GRAFT. GRAFT is a group-wise
motion estimation network leveraging the temporal redundancy to conduct a
1-to-$N$ motion estimation. This estimation is carried out $N$ times to
accomplish an $N$-to-$N$ motion correction for all temporal frames.

In this work we unroll the MCMR in Eq. 1 with a motion estimation network and
apply the joint optimization between the image reconstruction and motion
estimation. The framework is illustrated in Fig. 1. First, the zero-filled
undersampled image sequence $x_{0}$ is fed to the conjugate-gradient SENSE
block [27] along with $y$, $\mathbf{C}$, and $\mathbf{D}$. The first SENSE
block conducts a reconstruction without motion embedding. Afterwards, the
reconstructed image sequence $x_{1}$ is passed to a motion estimation network
(GRAFT, introduced in Section 2.2) as inputs and the first motion sequence
$u_{1}$ comprising of all pairs of motion frames is predicted. This motion
sequence is then applied to the input together with $y$, $\mathbf{C}$ and
$\mathbf{D}$ in the next CG-SENSE block, while the previous reconstruction is
used as a regularization. These CG-SENSE blocks solve Eq. 1 with an additional
$\ell_{2}$ regularization by the image from the previous step and output a
motion-corrected image sequence with temporarily freezing $\mathbf{U}_{i-1}$
and $x_{i-1}$, following

$\displaystyle
x_{i}=\arg\min\limits_{x}\sum_{t_{1}=1}^{N}\sum_{t_{2}=1}^{N}\left\|\mathbf{A}^{(t_{2})}\mathbf{U}_{i-1}^{(t_{1}\rightarrow
t_{2})}x^{(t_{1})}-y^{(t_{2})}\right\|^{2}_{2}+\frac{1}{2\lambda}\left\|x-x_{i-1}\right\|^{2}_{2}.$
(2)

This alternating scheme is carried out for $I$ iterations. The motion
estimation difficulty of GRAFT is alleviated gradually with the progress of
this alternating scheme along with the image-quality improvement.

###  2.2 Motion Estimation Network

There exists various motion estimation approaches in the field of medical
imaging. In this work, we select GRAFT [12] for motion estimation, which takes
the full image sequence $x$ together with the target frame $x^{(t)}$ as inputs
to predict the motion sequence $u^{(t)}=\left[u^{(t\rightarrow
1)},\ldots,u^{(t\rightarrow N)}\right]$ for all $N$ frames. GRAFT is suitable
to be embedded into our framework due to its accuracy, speed and efficiency.

Accuracy An accurate motion estimation is the essential requirement for the
unrolled MCMR framework. GRAFT contains a 4D correlation layer and iterative
Gated Recurrent Unit (GRU) blocks to conduct precise motion predictions. It
has shown superior results compared to conventional registration methods [12].
Furthermore, it is a group-wise estimation network considering all temporal
frames and has a dedicated frameblock leveraging temporal redundancies to
mitigate the impact of through-plane-motion. Because of this group-wise
attribute, temporal coherence can be instilled during the training by
appending temporal regularization in the loss function.

Speed Since the motion registration is invoked at each iteration step and for
all pairs of motion frames, a fast motion estimation method is required. The
reconstruction of $N$ frames involves $N^{2}$ motion predictions. In this
context, traditional registration-based methods such as [15, 34] require hours
to perform $N^{2}$ calculations, rendering them impractical for our framework.
In contrast, GRAFT only requires a few seconds to estimate the motion of all
pairs.

Efficiency For end-to-end learning over all iterative stages, the losses and
gradients need to be accumulated. In the context that $N^{2}$ calculations
(and for multiple iterations) are entailed, the training of large network
architectures with more than 20 million trainable parameters [5, 28] becomes
infeasible due to their vast GPU memory footprint (>48GB). In contrast, GRAFT
has only 5 million trainable parameters circumventing the GPU overcharge
problems.

In order to use GRAFT for motion estimation and integrate it in our unrolled
framework we define the following loss function to train GRAFT:

$\displaystyle\mathcal{L}_{i}=\sum_{t_{1}=1}^{N}\sum_{t_{2}=1}^{N}$
$\displaystyle\left\|\rho\left(\mathbf{U}_{i}^{(t_{1}\rightarrow
t_{2})}x_{gt}^{(t_{1})}-x_{gt}^{(t_{2})}\right)\right\|_{1}$ (3)
$\displaystyle+$
$\displaystyle\alpha\sum_{t_{1}=1}^{N}\sum_{t_{2}=1}^{N}\sum_{d\in
x,y}\left\|\nabla_{d}u_{i}^{(t_{1}\rightarrow
t_{2})}\right\|_{1}+\beta\sum_{t_{1}=1}^{N}\left\|\nabla_{t}u_{i}^{(t_{1})}\right\|_{1}.$

The first term ensures data fidelity during the training. We warp the ground-
truth frame $x_{gt}^{(t_{1})}$ to the target ground-truth frame
$x_{gt}^{(t_{2})}$ using the motion $u_{i}^{(t_{1}\rightarrow t_{2})}$, which
is predicted from GRAFT based on the previously reconstructed image $x_{i-1}$.
It should be noted that we apply the ground-truth frame in the loss function
to train GRAFT, so that GRAFT is forced to learn and extract genuine
correlations from the undersampled images and that any motion falsely
originating from image artifacts are not rewarded. The Charbonnier function
$\rho(x)=(x^{2}+10^{-12})^{0.45}$ [33, 25] is employed as the penalty
function. Additional regularization terms for the spatial plane weighted by
$\alpha$ and along the temporal axis weighted by $\beta$ are included in the
loss to mitigate estimation singularities and to ensure motion coherence in
spatial and temporal domain. This loss is calculated after every unrolled
motion estimation step with progressive exponential weight decay
$\mathcal{L}_{total}=\sum_{i=1}^{I}\gamma^{I-i}\mathcal{L}_{i}$.

##  3 Experiments

Training was carried out on 35 subjects (a mix of patients and healthy
subjects) of in-house acquired short-axis 2D CINE CMR, whereas testing was
performed on 5 subjects. Data were acquired with 30/32/34 multiple receiver
coils and 2D balanced steady-state free precession sequence on a 1.5T MR
(Siemens Aera with TE=1.06 ms, TR=2.12 ms,
resolution=$1.9{\times}1.9\text{mm}^{2}$ with 8mm slice thickness, 8 breath-
holds of 15s duration). Image sequence size varies from $176{\times}132$
(smallest) to $192{\times}180$ (largest) with 25 temporal cardiac phases. A
stack of 12 slices along the long axis was collected, resulting in 304/15
image sequence (2D+t) for training/test, while the apical slices without clear
cardiac contour were excluded.

The proposed framework with GRAFT as the embedded motion estimator was trained
on an NVIDIA A40 GPU with AdamW [20] (batch size of 1, learning rate 1e-4 and
weight decay of 1e-5). The number of unrolled iterations is set to 3 during
training, but this can be flexibly adapted during inference and iterations are
stopped when the peak signal to noise ratio (PSNR) converges (PSNR increment
$<0.1$). The trainable weights of GRAFT are shared during the iterative
optimization. The hyperparameters $\alpha$, $\beta$, $\gamma$ and $\lambda$
were set to 10, 10, 0.6 and 2 respectively. Three trainings and tests were
conducted separately on retrospectively undersampled images with VISTA [2] of
R = 8, 12 and 16 acceleration without any prior reconstruction or prior motion
correction as conducted in [7, 8, 29]. During the training and testing, raw
multi-coil k-space data were used for the reconstruction. We compare the
proposed approach with non-motion compensated CG-SENSE (N-CG-SENSE)
reconstruction and with non-iterative MCMR methods. In the latter case, the
motion is pre-computed from GRAFT [12] and Elastix [15] and set to be constant
during the whole optimization process.

##  4 Results and Discussion

Figure 2: Quantitative and qualitative results of proposed MCMR method. A)
Performance measured by peak signal to noise ratio (PSNR) for increasing
unrolled iteration numbers over all test subjects and temporal frames for
retrospective undersampling with R=8, 12 and 16. B) Qualitative image quality
of the proposed MCMR at R=16 over the first three iterations (i=1,2,3) in
comparison to the non-motion compensated CG-SENSE (N-CG-SENSE) at R=16 and the
fully-sampled reference image (R=1).

Fig. 2A reveals the relation between the estimation performance (indicated by
averaged PSNR over all test subjects and temporal frames) and the unrolled
iteration number of the proposed framework during inference for R=8, 12 and
16. The reconstruction accuracy is improved with increasing iteration numbers
for all three cases. The reconstruction for R=12 and R=16 benefit most from
our framework and obtain a significant performance lifting for the first three
iterations. The advantage of unrolling the motion estimator in the joint MCMR
optimization process is qualitatively illustrated in Fig. 2B, in which the MR
image is highly undersampled with R=16. Not only is the image quality improved
with the course of the iteration, but the motion estimator can also deliver
more precise and detailed estimation by virtue of the higher quality image. A
full qualitative analysis with R=8 and 12 is shown in Supplementary Fig. S1.
Furthermore, the motion estimation based on unprocessed and artifact-affected
images (from high undersampling rates) is challenging. We also trained GRAFT
without the proposed unrolled optimization framework [12] to generalize for
motion estimation from images with different undersampling rates. However,
training started to fail for accelerations R=12 and beyond. In contrast,
training difficulty is reduced if we embed GRAFT into the proposed joint
optimization framework while the image quality is restored with the progress
of the optimization. This makes GRAFT capable of estimating meaningful motion
as showcased in Supplementary Fig. S1 for R=12 and even for R=16. Moreover,
this mechanism also introduces a data augmentation process while more high-
quality data is generated during training.

Figure 3: Qualitative reconstruction performance of CINE CMR in spatial (x-y)
and spatio-temporal (y-t) plane of non-motion compensated CG-SENSE (N-CG-
SENSE), proposed MCMR at third iteration, CG-SENSE with pre-estimated motion
from GRAFT and Elastix, for different acceleration rates R=8, 12 and 16 in
comparison to the fully-sampled (R=1) reference. The selected superior-
inferior (y-axis) is marked with a blue line in the reference image. The
obtained PSNR in comparison to the fully-sampled reference is shown at the
bottom left. The highly accelerated cases for Elastix/GRAFT + CG-SENSE did not
converge therefore results are not shown.

A comparison study of the proposed MCMR to N-CG-SENSE and CG-SENSE
reconstruction with pre-computed motion based on GRAFT and Elastix is shown in
Fig. 3. The training of GRAFT alone failed for R=12 and 16, whereas motion
estimation with Elastix for R=16 was not conducted due to the poor quality
results collected at R=12. The results are shown after the third iteration for
the proposed MCMR, providing an optimal trade-off between time and
performance. Our approach restores the undersampled MR sequence with high
quality without artifacts in spatial and temporal domain. A quantitative
analysis over all test subjects including reconstruction of all 25 temporal
frames is shown in Table 1. Although the execution time of our approach is
longer than the original GRAFT if the iteration is set $>1$, an overall more
precise reconstruction is achieved. Furthermore, the overall reconstruction is
within acceptable clinical durations. The proposed MCMR enables to reconstruct
an undersampled image with a much higher acceleration rate and offers the high
flexibility to perform the reconstruction depending on different
time/performance requirement.

Table 1: Quantitative analysis of reconstruction for accelerated CINE CMR (R=8, 12 and 16) using the proposed MCMR method, non-motion compensated CG-SENSE (N-CG-SENSE), GRAFT+CG-SENSE and Elastix+CG-SENSE. Peak signal-to-noise ratio (PSNR) and structural similarity index (SSIM) [35] are used to evaluate all test subjects. Their mean value, standard deviations are shown next to the respective methods execution times. The best results are marked in bold. The failed or inferior experiments are marked with ’N.A.’. Acc R | Methods | SSIM | PSNR | Time (s)
---|---|---|---|---
8 | Proposed MCMR |  0.943 (0.018) |  36.26 (2.22) | 18.81s
GRAFT [12] \+ CG-SENSE | 0.913 (0.019) | 34.93 (1.80) | 6.27s
Elastix [15] \+ CG-SENSE | 0.645 (0.057) | 25.04 (2.10) | 4281s
N-CG-SENSE | 0.821 (0.038) | 30.80 (2.15) | 1.37s
12 | Proposed MCMR |  0.932 (0.018) |  35.45 (2.00) | 18.81s
GRAFT + CG-SENSE | N.A. | N.A. | 6.27s
Elastix + CG-SENSE | 0.568 (0.072) | 23.51 (2.20) | 4281s
N-CG-SENSE | 0.637 (0.062) | 24.40 (2.39) | 1.37s
16 | Proposed MCMR |  0.927 (0.019) |  34.78 (1.86) | 18.81s
GRAFT + CG-SENSE | N.A. | N.A. | 6.27s
Elastix + CG-SENSE | N.A. | N.A. | 4281s
N-CG-SENSE | 0.531 (0.08) | 21.736 (2.45) | 1.37s

The proposed MCMR framework also has some limitations: Currently, the proposed
MCMR cannot guarantee that the estimated motion is diffeomorphic. The
possibility to integrate the scaling and squaring layer [10] to ensure
diffeomorphic motion estiamtion will be investigated in future work.
Furthermore, we have not studied the interaction between our explicit motion
correction framework with implicit motion correction (denoising) networks,
which is also subject to future work. Moreover, we plan to evaluate the
performance of our method on prospectively undersampled data in the future.
Finally, the hyper-parameters applied in this work are estimated empirically
and they might not be the optimal combination. Integration of the learnable
hyper-parameters tuning [22] will be considered in our next step.

##  5 Conclusion

In this work, we proposed a deep learning-based MCMR framework and studied it
in CINE CMR imaging. We explored the possibility to unroll the MCMR
optimization with a motion estimation network while a dynamic joint
optimization between the reconstruction and the motion estimation is carried
out. Although this idea is still at a nascent stage, its potential in high
quality reconstruction of highly accelerated data can be appreciated. The
conducted experiments against baseline methods showcased a rapid, more robust
and more precise reconstruction of our proposed framework.

##  6 Acknowledgements

This work was supported in part by the European Research Council (Grant
Agreement no. 884622).

## References

  * [1] Aggarwal, H.K., Mani, M.P., Jacob, M.: Model based image reconstruction using deep learned priors (modl). In: IEEE International Symposium on Biomedical Imaging (ISBI). pp. 671–674 (2018)
  * [2] Ahmad, R., Xue, H., Giri, S., et al.: Variable density incoherent spatiotemporal acquisition (vista) for highly accelerated cardiac mri. Magnetic Resonance in Medicine 74(5), 1266–1278 (2015)
  * [3] Aviles-Rivero, A.I., Debroux, N., Williams, G., et al.: Compressed sensing plus motion (cs + m): A new perspective for improving undersampled mr image reconstruction. Medical Image Analysis 68, 101933 (2021)
  * [4] Balakrishnan, G., Zhao, A., Sabuncu, M.R., et al.: An unsupervised learning model for deformable medical image registration. In: 2018 IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 9252–9260 (2018)
  * [5] Balakrishnan, G., Zhao, A., Sabuncu, M.R., et al.: Voxelmorph: A learning framework for deformable medical image registration. IEEE Transactions on Medical Imaging 38(8), 1788–1800 (2019)
  * [6] Batchelor, P., Atkinson, D., Irarrazaval, P., Hill, D., et al.: Matrix description of general motion correction applied to multishot images. Magnetic Resonance in Medicine 54, 1273–1280 (2005)
  * [7] Bustin, A., Rashid, I., Cruz, G., et al.: 3d whole-heart isotropic sub-millimeter resolution coronary magnetic resonance angiography with non-rigid motion-compensated prost. Journal of Cardiovascular Magnetic Resonance 22(1) (2020)
  * [8] Cruz, G., Hammernik, K., Kuestner, T., et al.: One-heartbeat cardiac cine imaging via jointly regularized non-rigid motion corrected reconstruction. In: Proc. International Society for Magnetic Resonance in Medicine (ISMRM). p. 0070 (2021)
  * [9] Cruz, G., Atkinson, D., Henningsson, M., et al.: Highly efficient nonrigid motion-corrected 3d whole-heart coronary vessel wall imaging. Magnetic Resonance in Medicine 77(5), 1894–1908 (2017)
  * [10] Dalca, A.V., Balakrishnan, G., Guttag, J., Sabuncu, M.R.: Unsupervised learning of probabilistic diffeomorphic registration for images and surfaces. Medical Image Analysis 57, 226–236 (2019)
  * [11] Hammernik, K., Klatzer, T., Kobler, E., et al.: Learning a variational network for reconstruction of accelerated mri data. Magnetic Resonance in Medicine 79(6), 3055–3071 (2018)
  * [12] Hammernik, K., Pan, J., Rueckert, D., Küstner, T.: Motion-Guided Physics-Based Learning for Cardiac MRI Reconstruction. Asilomar Conference on Signals, Systems, and Computers (2021)
  * [13] Huang, W., Ke, Z., Cui, Z.X., et al.: Deep low-rank plus sparse network for dynamic mr imaging. Medical Image Analysis 73, 102190 (2021)
  * [14] Jung, H., Sung, K., Nayak, K.S., et al.: k-t focuss: A general compressed sensing framework for high resolution dynamic mri. Magnetic Resonance in Medicine 61(1), 103–116 (2009)
  * [15] Klein, S., Staring, M., Murphy, K., et al.: Elastix: a toolbox for intensity-based medical image registration. IEEE Trans Med Imag 29(1), 196–205 (2009)
  * [16] von Knobelsdorff-Brenkenhoff, F., Pilz, G., Schulz-Menger, J.: Representation of cardiovascular magnetic resonance in the aha/acc guidelines. Journal of Cardiovascular Magnetic Resonance 19(1), 1–21 (2017)
  * [17] Küstner, T., Fuin, N., Hammernik, K., et al.: CINENet: deep learning-based 3D cardiac CINE MRI reconstruction with multi-coil complex-valued 4D spatio-temporal convolutions. Scientific Reports 10(1), 1–13 (2020)
  * [18] Lee, D., Markl, M., Dall’Armellina, E., et al.: The growth and evolution of cardiovascular magnetic resonance: A 20-year history of the society for cardiovascular magnetic resonance (scmr) annual scientific sessions. Journal of Cardiovascular Magnetic Resonance 20(1) (Jan 2018)
  * [19] Liu, F., Li, D., Jin, X., et al.: Dynamic cardiac mri reconstruction using motion aligned locally low rank tensor (mallrt). Magnetic Resonance Imaging 66, 104–115 (2020)
  * [20] Loshchilov, I., Hutter, F.: Decoupled weight decay regularization. arXiv preprint arXiv:1711.05101 (19) (2017)
  * [21] Modat, M., Ridgway, G.R., Taylor, Z.A., et al.: Fast free-form deformation using graphics processing units. Computer Methods and Programs in Biomedicine 98(3), 278–284 (2010)
  * [22] Mok, T.C., Chung, A.: Conditional deformable image registration with convolutional neural network. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 35–45. Springer (2021)
  * [23] Odille, F., Vuissoz, P., Marie, P., Felblinger, J.: Generalized reconstruction by inversion of coupled systems (grics) applied to free-breathing mri. Magnetic Resonance in Medicine 60, 146–157 (2008)
  * [24] Odille, F., Menini, A., Escanyé, J.M., et al.: Joint reconstruction of multiple images and motion in mri: Application to free-breathing myocardial ${\rm t}_{2}$ quantification. IEEE Transactions on Medical Imaging 35(1), 197–207 (2016)
  * [25] Pan, J., Rueckert, D., Küstner, T., Hammernik, K.: Efficient image registration network for non-rigid cardiac motion estimation. In: Haq, N., Johnson, P., Maier, A., Würfl, T., Yoo, J. (eds.) Machine Learning for Medical Image Reconstruction. pp. 14–24 (2021)
  * [26] Poddar, S., Jacob, M.: Dynamic mri using smoothness regularization on manifolds (storm). IEEE Transactions on Medical Imaging 35(4), 1106–1115 (2016)
  * [27] Pruessmann, K.P., Weiger, M., Börnert, P., Boesiger, P.: Advances in sensitivity encoding with arbitrary k-space trajectories. Magn Reson Med (46), 638–651 (2001)
  * [28] Qi, H., Fuin, N., Cruz, G., et al.: Non-rigid respiratory motion estimation of whole-heart coronary mr images using unsupervised deep learning. IEEE Transactions on Medical Imaging 40(1), 444–454 (2021)
  * [29] Qi, H., Hajhosseiny, R., Cruz, G., et al.: End-to-end deep learning nonrigid motion-corrected reconstruction for highly accelerated free-breathing coronary mra. Magnetic Resonance in Medicine 86(1), 1983–1996 (2021)
  * [30] Qin, C., Duan, J., Hammernik, K., et al.: Complementary time-frequency domain networks for dynamic parallel mr image reconstruction. Magnetic Resonance in Medicine 86(6), 3274–3291 (2021)
  * [31] Sandino, C.M., Lai, P., Vasanawala, S.S., Cheng, J.Y.: Accelerating cardiac cine mri using a deep learning-based espirit reconstruction. Magnetic Resonance in Medicine 85(1) (2021)
  * [32] Schmoderer, T., Aviles-Rivero, A.I., Corona, V., et al.: Learning optical flow for fast MRI reconstruction. Inverse Problems 37(9), 095007 (aug 2021)
  * [33] Sun, D., Yang, X., Liu, M., Kautz, J.: Pwc-net: Cnns for optical flow using pyramid, warping, and cost volume. In: IEEE Conference on Computer Vision and Pattern Recognition (CVPR). pp. 8934–8943 (2018)
  * [34] Vercauteren, T., Pennec, X., Perchant, A., Ayache, N.: Diffeomorphic demons: Efficient non-parametric image registration. Neuro Image 45(1), S61–S72 (2009)
  * [35] Wang, Z., Bovik, A.C., Sheikh, H.R., Simoncelli, E.P.: Image quality assessment: from error visibility to structural similarity. IEEE Trans Image Process 13(4), 600–612 (2004)