∎

11institutetext: 1 Fujian Key Laboratory of Financial Information Processing,
Putian University, Putian, 351100, China.
2 Faculty of Mathematics, Technische Universit$\ddot{a}$t Chemnitz, Chemnnitz,
D-09107, Germany.

# A splitting-based KPIK method for eddy current optimal control problems in
an all-at-once approach††thanks: The work is supported by the National Natural
Science Foundation of China (Grant Nos. 12126344, 12126337 and 11901324) and
the China Scholarship Council (No. 202308350044).

Min-Li Zeng 1,2,∗ Martin Stoll2 ∗ Corresponding author<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

In this paper, we focus on efficient methods to fast solve discretized linear
systems obtained from eddy current optimal control problems in an all-at-once
approach. We construct a new low-rank matrix equation method based on a
special splitting of the coefficient matrix and the Krylov-plus-inverted-
Krylov (KPIK) algorithm. Firstly, we rewrite the resulting discretized linear
system in a matrix-equation form. Then using the KPIK algorithm, we can obtain
the low-rank approximation solution. The new method is named the splitting-
based Krylov-plus-inverted-Krylov (SKPIK) method. The SKPIK method can not
only solve the large and sparse discretized systems fast but also overcomes
the storage problem. Theoretical results about the existence of the low-rank
solutions are given. Numerical experiments are used to illustrate the
performance of the new low-rank matrix equation method by compared with some
existing classical efficient methods.

###### Keywords:

Eddy Current ProblemsPDE-constrained OptimizationLow-rank ApproximationAll-at-
once ApproachMatrix Splitting

###### MSC:

65F10 65F50 93C20

##  1 Introduction

In this work, we are interested in the numerical solution of the distributed
optimal control problem: Find the state $y$ and the control $u$ that minimizes
the cost functional

$\displaystyle
J(y,u)=\frac{1}{2}\int_{0}^{T}\int_{\Omega}|y-y_{d}|^{2}dxdt+\frac{\beta}{2}\int_{0}^{T}\int_{\Omega}|u|^{2}dxdt,$
(1)

subject to

$\displaystyle\left\\{\begin{array}[]{rlcl}\sigma\frac{\partial}{\partial
t}y+\textbf{curl}(\nu\textbf{curl}y)&=u,&\text{in}&\Omega\times(0,T),\\\
y\times n&=0,&\text{on}&\partial\Omega\times(0,T),\\\
y&=y_{0},&\text{on}&\Omega\times\\{0\\},\\\ \end{array}\right.$ (5)

where $y_{d}$ is the desired state, $\sigma\in L^{\infty}(\Omega)$ denotes the
conductivity, $\nu\in L^{\infty}(\Omega)$ is the reluctivity parameter and
$\beta>0$ is a cost parameter. We assume that $\Omega\in\mathbb{R}^{2}$ or
$\mathbb{R}^{3}$ is a bounded Lipschitz domain.

The solution of eddy current optimal control problems is one of the most
interesting and demanding problems in applied mathematics and scientific
computing. Eddy current optimal control problems are modeled by Maxwell’s
equations kolmbauer2013efficient ; kolmbauer2011frequency . Typically, this
category of issues encompasses a set of interrelated variable functions, such
as the electromagnetic force, electric current, and an associated adjoint
variable function. Over the last decades, the numerical solution of the eddy
current optimal control problems has received lots of attention
axelsson2019note ; axelsson2017preconditioners ; axelsson2019preconditioning ;
bachinger2006efficient .

When the eddy current optimal control problem is time-harmonic$\setminus$time-
periodic, Kolmbauer and Langer used the multi-harmonic finite element method
(FEM) to discretize the problem in kolmbauer2011frequency ;
kolmbauer2012robust  and then constructed a robust preconditioned minimum
residual (MINRES) solver for the corresponding disretized linear system. In
contrast to previous approaches, Kolmbauer further derived the problem setting
in a mixed variational formulation under the cases of different control and
observation domains, observation at a final time, constraints to the control
and the state and observation in certain energy spaces. More details can be
found in kolmbauer2013robust ; kolmbauer2013efficient ; bachinger2006efficient
. During the past two decades, more and more research concerning the
discretization systems and the corresponding preconditioned methods has been
carried out. For example, in wolfmayr2023posteriori , Wolfmayr presents the
multi-harmonic analysis and derivation of functional type a posteriori
estimates of a distributed eddy current optimal control problem. Axelsson and
co-authors presented the efficient preconditioners for the time-harmonic
optimal control eddy-current problems in axelsson2017preconditioners  and for
eddy-current optimally controlled time-harmonic electromagnetic problem in
axelsson2019preconditioning , respectively. We refer to buffa2000justification
; axelsson2019note ; zeng2021respectively  and references therein for more
efficient preconditioning techniques for the time-harmonic$\setminus$time-
periodic eddy current optimal control problems.

When the eddy current problems are not time-harmonic$\setminus$time-periodic,
then the first-order temporal derivative is complicated. Because the classical
time-stepping method, i.e., solving the partial differential equations (PDEs)
one-time step after one-time step, would be time-consuming if the number of
time steps is large. Therefore, the so-called all-at-once approach technique
hinze2008hierarchical ; hinze2012space ; benzi2011preconditioning ;
stoll2013all  becomes popular. The all-at-once approach involves discretizing
the relevant issue within the space-time continuum and resolving it for all
temporal intervals simultaneously. A key benefit of this method is the
deferment of meeting the optimal control problem criteria until the entire
system converges. Nonetheless, a significant drawback is the resultant large-
scale systems inherent to the all-at-once approach. Hence, developing
efficient algorithms for fast solving large-scale systems has always been a
research focus.

The first motivation for fast solving large-scale systems is the development
of parallel-in-time (PinT) methods for evolutionary PDEs lions2001resolution .
Among these, we mention the parareal algorithm maday2002parareal  and a
closely related algorithm multigrid-reduction-in-time (MGRiT) algorithm
falgout2014parallel , which attracted considerable attention in recent years.
The convergence properties of the parareal algorithm and MGRiT are discussed
in maday2002parareal ; gander2007analysis . Many efforts are devoted to
improving these two PinT algorithms, and in particular the authors
wu2018toward ; wu2019acceleration  proposed a novel coarse grid correction,
which shows great potential for increasing the speedup according to the
numerical results in kwok2019schwarz . There are also many other PinT
algorithms with completely different mechanisms from the parareal algorithm
and MGRiT, such as the space-time multigrid algorithms gander2016analysis ;
horton1995space  and the diagonalization-based all-at-once algorithms
gander2016direct ; mcdonald2018preconditioning ; goddard2018note . For an
overview, we refer the interested reader to gander201550 . However, in
contrast to non-PinT techniques, a limitation of PinT methods is their
substantial demand for experimental resources, such as a large array of
processors, leading to increased experimental expenses. Additionally, PinT
methods continue to face storage challenges due to the demands of space-time
discretization.

To solve large-scale systems based on an all-at-once approach efficiently,
people usually consider using iteration methods, for example, Krylov subspace
methods combined with efficient preconditioners. Much progress has been made
in recent years concerning efficient preconditioners for the all-at-once-based
discretized systems for different applications. In rees2010all , Rees and co-
authors discussed implementation details on the all-at-once approach with
different boundary conditions. McDonald, Pestana, and Wathen established a
block circulant preconditioner for the all-at-once evolutionary PDE system in
mcdonald2018preconditioning . More applications can be found in
haber2001preconditioned ; zhao2021preconditioning  and references therein.

However, it is easy to see that the models based on the all-at-once approach
result in millions of variables. Therefore, it needs huge memory to store the
large and dense solution computed by standard algorithms. Nevertheless,
neither PinT methods nor the previously described techniques can overcome the
existing memory limitations. Ultimately, due to their size and complexity,
solving such huge models is still a challenging task.

Therefore, efficient solvers based on low memory requirements are desirable.
Recently, Bünger and co-authors proposed a new algorithm framework
bunger2021low , which can be used to solve PDE-constrained optimization
problems efficiently. To overcome the obstacles brought out by the storage
requirement of vast dimensionality, one can reformulate the Karush-Kuhn-Tucker
(KKT) system into a Sylvester-like matrix equation stoll2015low , which would
provide an efficient representation with a minimal amount of storage. Then
based on this reformulation, we can search for the low-rank approximation
solution by using the idea of low-rank in time technique, which can be found
in markovsky2012low . The main idea of the low-rank solver is to compute a
projection of the solution onto a small subspace. Hence, the low-rank solver
is similar to model order reduction (MOR) approaches. However, in contrast to
classical MOR schemes, the low-rank solvers do not need to compress the full
solution at the end of the algorithm but start with low-rank data and maintain
this form throughout the iteration. Hence, more and more attention has been
paid to the low-rank solvers. For example, the low-rank approximation solution
algorithm has been used for solving the optimal control problem constrained
with a forward Navier-Stokes equation by Dolgov in dolgov2017low , the optimal
control problem constrained with unsteady Stokes-Brinkman equations by Benner
and co-authors in benner2015low  and two-dimensional time-dependent Navier-
Stokes equations with a stochastic by Benner and co-authors in benner2016low .
More low-rank approximation algorithms for optimization of large-scale systems
have been studied in existing papers, and we refer the readers to
benner2013low ; benner2016block ; cichocki2017tensor ; dolgov2016fast ;
dolgov2017low ; benner2023low  for recent accounts.

In this work, we focus on the efficient algorithms based on the low-rank
matrix equation method for solving the eddy current optimal control problems
by an all-at-once approach. Our approach follows the work presented in
bunger2021low . Before rewriting the all-at-once discretized system as a
matrix equation, we make a matrix-splitting based on the special structure of
the coefficient matrix. Following this splitting, we can obtain an equivalent
matrix equation. We reduce the system of the matrix equation into a low-rank
reduced matrix approximation equation under an appropriate approximation
space. For the resulting small-size matrix equation, we use the KPIK algorithm
and then obtain the approximate solution.

The organization of the paper is as follows. In Section 2, we present the
discretized linear system of the eddy current optimization problems (1)-(5)
based on an all-at-once approach. In Section 3, we describe the MOR system and
propose the new low-rank matrix equation method. Theoretical results are
investigated in detail in Section 4. Numerical experiments are carried out in
Section 5 to show the effectiveness of the new method. Finally, we draw some
brief concluding remarks in Section 6.

The following notation is used throughout. Range($V$) denotes the space
spanned by the columns of $V$. We use the notation $\|\cdot\|$ and
$\|\cdot\|_{F}$ to indicate the 2-norm and the Frobenius norm for vectors and
the induced norm for matrices. The dimension of the identity matrix $I$ and
the zero matrix $O$ will be omitted when it is easy to distinguish.

##  2 Discretization and Optimization

In this section, we make a detailed description of the eddy current
optimization problem and then present the discretization of the all-at-once
approach of the model.

We assume the control $u$ is weakly divergence-free. The reluctivity $\nu$ is
supposed to be uniformly positive and independent of $|\textbf{curl}y|$, i.e.,
we assume that the eddy current problem (5) is linear. Eddy current problems
are essentially different for conducting ($\sigma>0$) and nonconducting
regions ($\sigma=0$).

To gain uniqueness in the nonconducting regions, we need to regularize the
state equation (5) by introducing formal regularization operators ($i=1,2,3$)
as:

$\mathcal{R}_{i}(\sigma):=\left\\{\begin{array}[]{lr}\sigma,&i=0,\\\
\sigma,&i=1,\\\ \max(\sigma,\varepsilon),&i=2,\\\
\sigma,&i=3,\end{array}\right.\qquad\qquad\mathcal{Q}_{i}(y):=\left\\{\begin{array}[]{lrl}0,&i=0&(\text{no\
regularization}),\\\ Q(y),&i=1&(\text{exact}),\\\
0,&i=2&(\text{conductivity}),\\\ \varepsilon
y,&i=3&(\text{elliptic}),\end{array}\right.$

where $i=0$ refers to without any regularization. Hence, the regularized
problem can be stated as

$J(y,u)=\min_{(y,u)}\frac{1}{2}\int_{0}^{T}\int_{\Omega}|y-y_{d}|^{2}dxdt+\frac{\beta}{2}\int_{0}^{T}\int_{\Omega}|u|^{2}dxdt,$

subject to

$\displaystyle\left\\{\begin{array}[]{rlcl}\mathcal{R}_{i}(\sigma)\frac{\partial
y}{\partial
t}+\textbf{curl}(\nu\textbf{curl}y)+\mathcal{Q}_{i}(y)&=u,&\text{in}&\Omega\times(0,T),\\\
y\times n&=0,&\text{on}&\partial\Omega\times(0,T),\\\
y&=y_{0},&\text{on}&\Omega\times\\{0\\}.\\\ \end{array}\right.$ (9)

Construct the Lagrangian functional

$\mathcal{L}(y,u,p)=J(y,u)+\frac{1}{2}\int_{0}^{T}\int_{\Omega}(\mathcal{R}_{i}(\sigma)\frac{\partial
y}{\partial t}+\textbf{curl}(\nu\textbf{curl}y)+\mathcal{Q}_{i}(y)-u)pdxdt$

and obtain the necessary optimality conditions

$\left\\{\begin{array}[]{rl}\nabla_{p}\mathcal{L}(y,u,p)&=0,\\\
\nabla_{y}\mathcal{L}(y,u,p)&=0,\\\ \nabla_{u}\mathcal{L}(y,u,p)&=0.\\\
\end{array}\right.$

It follows the system of PDEs: Find the state $y$, the costate $p$, and the
control $u$ such that

$\left\\{\begin{array}[]{rlcl}\mathcal{R}_{i}(\sigma)\frac{\partial
y}{\partial
t}+\textbf{curl}(\nu\textbf{curl}y)+\mathcal{Q}_{i}(y)-u&=0,&\text{in}&\Omega\times(0,T),\\\
-\mathcal{R}_{i}(\sigma)\frac{\partial p}{\partial
t}+\textbf{curl}(\nu\textbf{curl}p)+\mathcal{Q}_{i}(p)+y-y_{d}&=0,&\text{in}&\Omega\times(0,T),\\\
u&=\frac{1}{\beta}p,&\text{in}&\Omega\times(0,T),\\\ y\times
n&=0,&\text{on}&\partial\Omega\times(0,T),\\\ p\times
n&=0,&\text{on}&\partial\Omega\times(0,T),\\\
y&=y_{0},&\text{on}&\Omega\times\\{0\\},\\\
p&=0,&\text{on}&\Omega\times\\{0\\}.\\\ \end{array}\right.$

Next, we discretize space and time using an all-at-once approach. We use the
equal-order finite elements with spatial grid numbers being $N+1$ and to
discretize the space. The temporal discretization is done with the time
interval being split into $m_{T}$ intervals of length $\tau=\frac{T}{m_{T}}$.
Using a rectangle rule, the discretization of (1) and (9) lead to

$\sum_{m=1}^{m_{T}}\frac{\tau}{2}(y_{m}-y_{d,m})^{T}M(y_{m}-y_{d,m})+\frac{\tau\beta}{2}u_{m}^{T}Mu_{m},$
(10)

where $y_{m},y_{d,m}$ and $u_{m}$ are spatial discretizations of $y,y_{d}$ and
$u$ of size $n$ for each time step $m=1,\ldots,m_{T}$. Using an implicit
Euler-scheme the discrete formulation of (5) will lead to

$\displaystyle\frac{M_{\sigma}(y_{m}-y_{m-1})}{\tau}+Ky_{m}$
$\displaystyle=Mu_{m},\mbox{ for }m=1,\ldots,m_{T},$ $\displaystyle
y_{m}\times n$ $\displaystyle=0,$

where the mass matrix $M\in\mathbb{R}^{n\times n}$, the conductivity matrix
$M_{\sigma}\in\mathbb{R}^{n\times n}$ and stiffness matrix
$K\in\mathbb{R}^{n\times n}$ arise from the finite element discretization of
the following bilinear forms elman2005 :

$M:\int_{\Omega}u\cdot vdx,\qquad
M_{\sigma}:\int_{\Omega}\mathcal{R}_{i}(\sigma)u\cdot vdx,$
$K:\int_{\Omega}\nu\textbf{curl}u\cdot\textbf{curl}vdx+\int_{\Omega}\mathcal{Q}_{i}(u)vdx.$

It then follows, for $m=1,2,\ldots,m_{T}$,

$M_{\sigma}(y_{m}-y_{m-1})+\tau Ky_{m}-\tau Mu_{m}=0,$

i.e.,

$(M_{\sigma}+\tau K)y_{m}-M_{\sigma}y_{m-1}-\tau Mu_{m}=0.$

By collecting the discretizations of the variables in matrices
$Y=[y_{1},\ldots,y_{m_{T}}]\in\mathbb{R}^{n\times m_{T}}$ and denote their
vectorization by $\underline{Y}=\mbox{vec}(Y)$. The denotations
$\underline{Y_{d}}$ and $\underline{U}$ are used respectively for $y_{d}$ and
$u$. Therefore, we can rewrite the optimization problem in compact form as

$\min_{Y,U}\frac{\tau}{2}(\underline{Y}-\underline{Y_{d}})^{T}\mathcal{M}(\underline{Y}-\underline{Y_{d}})+\frac{\tau\beta}{2}\underline{U}^{T}\mathcal{M}\underline{U},$

such that

$\mathcal{N}_{\sigma}\underline{Y}-\tau\mathcal{M}\underline{U}=0,$

where

$\displaystyle\mathcal{M}=\begin{bmatrix}M&&&\\\ &M&&\\\ &&\ddots&\\\
&&&M\end{bmatrix},\qquad\mathcal{N}_{\sigma}=\begin{bmatrix}M_{\sigma}+\tau
K&&&\\\ -M_{\sigma}&M_{\sigma}+\tau K&&\\\ &\ddots&\ddots\\\
&&-M_{\sigma}&M_{\sigma}+\tau K\end{bmatrix},$

with the mass matrix $M\in\mathbb{R}^{n\times n}$, the conductivity matrix
$M_{\sigma}\in\mathbb{R}^{n\times n}$ and stiffness matrix
$K\in\mathbb{R}^{n\times n}$ repeating $m_{T}$ times each.

Next, we need to address the system of equations derived from the first-order
optimality conditions. In benzi2005numerical , Benzi and co-authors stated
that an optimal solution must be a saddle point of the Lagrangian of the
problem,

$\nabla\mathcal{L}(\underline{Y}^{*},\underline{U}^{*},\underline{P}^{*})=0.$

The Lagrangian of this problem reads

$\mathcal{L}(\underline{Y},\underline{U},\underline{P})=\frac{\tau}{2}(\underline{Y}-\underline{Y_{d}})^{T}\mathcal{M}(\underline{Y}-\underline{Y_{d}})+\frac{\tau\beta}{2}\underline{U}^{T}\mathcal{M}\underline{U}+\underline{P}^{T}(\mathcal{N}_{\sigma}\underline{Y}-\tau\mathcal{M}\underline{U}).$

Thus, the optimal solution solves the following set of linear equations,

$\left\\{\begin{array}[]{rcl}\nabla_{Y}\mathcal{L}(\underline{Y},\underline{U},\underline{P)}&=&\tau\mathcal{M}(\underline{Y}-\underline{Y_{d}})+\mathcal{N}_{\sigma}^{T}\underline{P}=0,\\\
\nabla_{U}\mathcal{L}(\underline{Y},\underline{U},\underline{P)}&=&\tau\beta\mathcal{M}\underline{U}-\tau\mathcal{M}^{T}\underline{P}=0,\\\
\nabla_{P}\mathcal{L}(\underline{Y},\underline{U},\underline{P})&=&\mathcal{N}_{\sigma}\underline{Y}-\tau\mathcal{M}\underline{U}=0.\end{array}\right.$
(11)

By introducing the auxiliary matrices

$\displaystyle\mathcal{I}=\begin{bmatrix}1&&\\\ &\ddots&\\\
&&1\end{bmatrix}\in\mathbb{R}^{m_{T}\times m_{T}}\mbox{ and
}C=\begin{bmatrix}1&&&\\\ -1&1&&\\\ &\ddots&\ddots&\\\
&&-1&1\end{bmatrix}\in\mathbb{R}^{m_{T}\times m_{T}},$

we can rewrite the above system matrices as Kronecker products
$\mathcal{M}=\mathcal{I}\otimes M$ and
$\mathcal{N}_{\sigma}=\mathcal{I}\otimes\tau K+C\otimes M_{\sigma}$. With
this, our system of KKT conditions (11) becomes

$\left\\{\begin{array}[]{rcc}\tau(\mathcal{I}\otimes
M)\underline{Y}+(\mathcal{I}\otimes\tau K^{T}+C^{T}\otimes
M_{\sigma}^{T})\underline{\Lambda}-\tau(\mathcal{I}\otimes
M)\underline{Y_{d}}&=&0,\\\ \tau\beta(\mathcal{I}\otimes
M)\underline{U}-\tau(\mathcal{I}\otimes M^{T})\underline{\Lambda}&=&0,\\\
(\mathcal{I}\otimes\tau K+C\otimes
M_{\sigma})\underline{Y}-\tau(\mathcal{I}\otimes
M)\underline{U}&=&0.\end{array}\right.$ (12)

Because the mass matrix $M$ arising from a standard Galerkin method is always
symmetric and positive definite, i.e., $M=M^{T}$ and $M^{-1}$ exists.
Therefore, we can eliminate the second equation in Equation (12) by setting
$U=\frac{1}{\beta}\Lambda$ in the remaining two equations. It results in the
following equations:

$\begin{pmatrix}\tau\mathcal{M}&\sqrt{\beta}\mathcal{N}_{\sigma}^{T}\\\
\sqrt{\beta}\mathcal{N}_{\sigma}&-\tau\mathcal{M}\end{pmatrix}\begin{pmatrix}\underline{Y}\\\
\frac{1}{\sqrt{\beta}}\underline{\Lambda}\end{pmatrix}=\begin{pmatrix}\tau(\mathcal{I}\otimes
M)\underline{Y_{d}}\\\ 0\end{pmatrix}.$ (13)

##  3 The splitting-based KPIK method

In this section, we present the splitting-based KPIK method for the linear
system (13). By substituting $\mathcal{M}=I\otimes M$,
$\mathcal{N}_{\sigma}=\mathcal{I}\otimes\tau K+C\otimes M_{\sigma}$ and
$M_{\sigma}=\sigma M$, and $K=K^{T}$, we can rewrite the linear system (13) as

$\begin{pmatrix}\tau\mathcal{I}\otimes M&\sqrt{\beta}(\mathcal{I}\otimes\tau
K+C\otimes\sigma M)^{T}\\\ \sqrt{\beta}(\mathcal{I}\otimes\tau
K+C\otimes\sigma M)&-\tau\mathcal{I}\otimes
M\end{pmatrix}\begin{pmatrix}\underline{Y}\\\
\frac{1}{\sqrt{\beta}}\underline{\Lambda}\end{pmatrix}=\begin{pmatrix}\tau(\mathcal{I}\otimes
M)\underline{Y_{d}}\\\ 0\end{pmatrix},$

i.e.,

$\Big{(}\begin{pmatrix}\tau\mathcal{I}\otimes M&\sqrt{\beta}C^{T}\otimes\sigma
M\\\ \sqrt{\beta}C\otimes\sigma M&-\tau\mathcal{I}\otimes
M\end{pmatrix}+\begin{pmatrix}0&\sqrt{\beta}\mathcal{I}\otimes\tau K\\\
\sqrt{\beta}\mathcal{I}\otimes\tau
K&0\end{pmatrix}\Big{)}\begin{pmatrix}\underline{Y}\\\
\frac{1}{\sqrt{\beta}}\underline{\Lambda}\end{pmatrix}=\begin{pmatrix}\tau(\mathcal{I}\otimes
M)\underline{Y_{d}}\\\ 0\end{pmatrix}.$

We can further rewrite the above matrix equation as

$\Big{(}\begin{pmatrix}\tau\mathcal{I}&\sigma\sqrt{\beta}C^{T}\\\
\sigma\sqrt{\beta}C&-\tau\mathcal{I}\end{pmatrix}\otimes
M+\begin{pmatrix}0&\tau\sqrt{\beta}\mathcal{I}\\\
\tau\sqrt{\beta}\mathcal{I}&0\end{pmatrix}\otimes
K\Big{)}\begin{pmatrix}\underline{Y}\\\
\frac{1}{\sqrt{\beta}}\underline{\Lambda}\end{pmatrix}=\begin{pmatrix}\tau(\mathcal{I}\otimes
M)\underline{Y_{d}}\\\ 0\end{pmatrix}.$ (14)

Using the relation

$(W^{T}\otimes V)\mbox{vec}(X)=\mbox{vec}(VXW),$

we rewrite the linear system (14) in a matrix-equation form as

$MX\begin{pmatrix}\tau\mathcal{I}&\sigma\sqrt{\beta}C^{T}\\\
\sigma\sqrt{\beta}C&-\tau\mathcal{I}\end{pmatrix}+KX\begin{pmatrix}0&\tau\sqrt{\beta}\mathcal{I}\\\
\tau\sqrt{\beta}\mathcal{I}&0\end{pmatrix}=[\tau MY_{d}\quad
0]\in\mathbb{R}^{n\times 2m_{T}}.$ (15)

Multiplying $M^{-1}$ by the left and

$\begin{pmatrix}0&\tau\sqrt{\beta}\mathcal{I}\\\
\tau\sqrt{\beta}\mathcal{I}&0\end{pmatrix}^{-1}$

by the right of (15) on both sides , respectively, we can further obtain

$M^{-1}KX+X\begin{pmatrix}\frac{\sigma}{\tau}C^{T}&\frac{1}{\sqrt{\beta}}\mathcal{I}\\\
-\frac{1}{\sqrt{\beta}}\mathcal{I}&\frac{\sigma}{\tau}C\end{pmatrix}=M^{-1}[\tau
MY_{d}\quad 0]\begin{pmatrix}0&\frac{1}{\tau\sqrt{\beta}}\mathcal{I}\\\
\frac{1}{\tau\sqrt{\beta}}\mathcal{I}&0\end{pmatrix},$

i.e.,

$M^{-1}KX+X\begin{pmatrix}\frac{\sigma}{\tau}C^{T}&\frac{1}{\sqrt{\beta}}\mathcal{I}\\\
-\frac{1}{\sqrt{\beta}}\mathcal{I}&\frac{\sigma}{\tau}C\end{pmatrix}=[0\quad\frac{1}{\sqrt{\beta}}Y_{d}]:=R\in\mathbb{R}^{n\times
2m_{T}},$ (16)

where $X=[Y\quad\frac{1}{\sqrt{\beta}}\Lambda]\in\mathbb{R}^{n\times 2m_{T}}$.

Next, we focus on the low-rank approximation solution of the Sylvester matrix
equation

$AX+XB=R,$

where the left-hand coefficient matrix $A=M^{-1}K$ has size $n\times n$ and
the coefficient matrix

$B=\begin{pmatrix}\frac{\sigma}{\tau}C^{T}&\frac{1}{\sqrt{\beta}}\mathcal{I}\\\
-\frac{1}{\sqrt{\beta}}\mathcal{I}&\frac{\sigma}{\tau}C\end{pmatrix}$ (17)

has size $2m_{T}\times 2m_{T}$, while the right-hand one has size $n\times
2m_{T}$, so that $X\in{\mathbb{R}}^{n\times 2m_{T}}$.

Suppose that there exists a low-rank approximation of the desired state as

$Y_{d}\approx Y_{1}Y_{2}^{T},$

with $Y_{1}\in\mathbb{R}^{n\times r}$, $Y_{2}\in\mathbb{R}^{m_{T}\times r}$
and $r<m_{T}$ of low column and row rank, respectively. Then the approximate
low-rank decomposition of the right-hand side $R\approx R_{1}R_{2}^{T}$ is
given by

$R_{1}=\frac{1}{\sqrt{\beta}}Y_{1}\in{\mathbb{R}}^{n\times
r}\quad\text{and}\quad R_{2}=\begin{pmatrix}O_{m_{T}\times r}\\\
Y_{2}\end{pmatrix}\in{\mathbb{R}}^{2m_{T}\times r},$ (18)

with $R_{1}$ and $R_{2}$ being low column and row rank, respectively. Here,
$O_{m_{T}\times r}$ denotes a zero matrix of size $m_{T}\times r$. We will
omit the subscript when it is easy to distinguish in the following of this
paper.

Because the solution matrix $X\in\mathbb{R}^{n\times 2m_{T}}$ would be dense
and potentially very large, then we need to find an appropriate approximation
space to exploit the new setting of the solution matrix. Suppose that there
exists a low-rank reduced matrix approximation $Z\in\mathbb{R}^{p\times
2m_{T}}$ such that $X\approx V_{p}Z$, where the orthonormal columns of
$V_{p}\in\mathbb{R}^{n\times p}$ generate the approximation space. Hence, we
can construct a reduced version of the matrix equation (17).

Denote the reduced $p\times p$ coefficient matrices as
$A_{r}:=V_{p}^{T}AV_{p}$ and set $R_{1,r}=V_{p}^{T}R_{1}\in\mathbb{R}^{p\times
r}$. Then the resulting reduced equation can be written as

$A_{r}Z+ZB=R_{1,r}R_{2}^{T},$ (19)

which has the same structure as the original matrix equation but with its size
reduced to $p\times 2m_{T}$. Using the relation in (17), we get the small
linear system of equations

$\big{(}(\mathcal{I}_{2m_{T}}\otimes
A_{r})+(B^{T}\otimes\mathcal{I}_{p}))\underline{Z}=\underline{R_{1,r}}\underline{R_{2}}^{T},$
(20)

with $\underline{Z}=\mbox{vec}(Z)$, $\underline{R_{1,r}}=\mbox{vec}(R_{1,r})$,
and $\underline{R_{2}}=\mbox{vec}(R_{2})$.

Because of the small subspace size $p\ll n$, we can either use a direct method
or an iterative method to solve this system of equations, which is
significantly easier to solve than the original system of equations. If the
obtained approximate solution $V_{p}Z$ is not sufficiently good, then the
space can be expanded. Hence, a new approximation can be constructed, giving
rise to an iterative method. More details can be found in stoll2015low ;
breiten2016fast .

At the end of this section, we give a detailed implementation of the SKPIK
method in Algorithm 1.

Algorithm 1 The splitting-based KPIK method

1: Given a tolerance $\epsilon$, a maximum number of iteration
$\text{IT}_{max}$, a spatial grid points number $n$, a temporal grid points
number $m_{T}$ and a tolerance for the truncated singular value decomposition
(SVD) $\varepsilon$. The step-size parameter $\tau=\frac{1}{m_{T}}$.

2: Inputs: $M\in\mathbb{R}^{n\times n}$ a mass matrix,
$K\in\mathbb{R}^{n\times n}$ a stiffness matrix, $Y_{d}\in\mathbb{R}^{n\times
m_{T}}$ a given desired state vector, $\sigma>0$ a conducting constant and
$\beta>0$ a cost parameter.

3: Set $A\approx M^{-1}K$, $B$ given by (17), $R\approx R_{1}R_{2}^{T}$ within
truncated SVD $\varepsilon$, i.e., dropping the singular values of $R$ smaller
than $\varepsilon$ stoll2015low . $U_{1}=gram\\_sh(R,A^{-1}R)$,
$W_{1}=gram\\_sh(R,B^{-1}R)$, $U_{0}=\varnothing$ and $W_{0}=\varnothing$.

4: for $m=1,2,\ldots,\text{IT}_{max}$ do

5: $U_{m}=[U_{m-1},U_{m}]$, $W_{m}=[W_{m-1},W_{m}]$.

6: Set $T^{A}_{m}=W_{m}^{T}AU_{m}$, $T^{B}_{m}=W_{m}^{T}BU_{m}$ and
$R_{1}^{(m)}=U_{m}R_{1}$, $R_{2}^{(m)}=W_{m}^{T}R_{2}$.

7: Solve $T_{m}^{A}Y+Y(T_{m}^{B})^{T}+R_{1}R_{2}^{T}=0$ and set $Y_{m}=Y$.

8: If converged, the $X_{m}=U_{m}Y_{m}W_{m}^{T}$ and stop.

9: Set $U_{m}^{(1)}$: first $s$ columns of $U_{m}$, $U_{m}^{(2)}$: second $s$
columns of $U_{m}$; $W_{m}^{(1)}$: first $s$ columns of $W_{m}$,
$W_{m}^{(2)}$: second $s$ columns of $W_{m}$.

10: $U_{m+1}^{\prime}=[AU_{m}^{(1)},A^{-1}U_{m}^{(2)}]$ and
$W_{m+1}^{\prime}=[BW_{m}^{(1)},B^{-1}W_{m}^{(2)}]$.

11: $\hat{U}_{m+1}\leftarrow$ orthogonalize $U_{m+1}^{\prime}$ w.r.t. $U_{m}$
and $\hat{W}_{m+1}\leftarrow$ orthogonalize $W_{m+1}^{\prime}$ w.r.t. $W_{m}$.

12: $U_{m+1}=gram\\_sh(\hat{U}_{m+1})$ and $W_{m+1}=gram\\_sh(\hat{W}_{m+1})$.

13: if
$\frac{\|AU_{m}Y_{m}W_{m}^{T}+U_{m}Y_{m}W_{m}^{T}B^{T}+R_{1}R_{2}^{T}\|}{\|R_{1}R_{2}^{T}\|}<\epsilon$
then

14: $X_{1}=U_{m}Y_{m}$ and $X_{2}=W_{m}$. Break.

15: end if

16: end for

17: The low-rank approximation solution $X$ is given by $X\approx
X_{1}X_{2}^{T}$.{$\rhd$ The function “gram_sh” performs the modified Gram-
Schmidt orthogonalization, more details can be found in golub2013matrix . }

###### Remark 3.1

Compared with the low-rank MINRES method proposed in stoll2015low , the SKPIK
method in this work needs one to get the solution of size $n\times 2m_{T}$
from the Sylvester equation (16) once, while the low-rank MINRES method,
preconditioned by the block diagonal preconditioner with the approximate Schur
complement proposed in stoll2015low , needs one to obtain the solutions of
size $n\times m_{T}$ from two different Sylvester equations. However, the low-
rank MINRES method also needs one to solve sublinear systems with the
coefficient matrices being $\mathcal{M}$ twice. Hence, the SKPIK method has a
smaller workload than the low-rank MINRES method.

###### Remark 3.2

As the stiffness matrix $K$ is SPSD and the mass matrix $M$ is SPD, then
following the idea in shank2013krylov , we can also shift the matrix equation
$AX+XB=R_{1}R_{2}^{T}$ of (16) as a equivalent matrix equation
$(A+sI)X+X(B-sI)=R_{1}R_{2}^{T}$, where $s>0$ is a given shifted parameter,
such that $A+sI=M^{-1}(K+sM)$ is better conditioned than the matrix $M^{-1}K$.

##  4 The existence of the low-rank solution

In this section, we give the existence of the low-rank solution. The existence
of the low-rank approximant to a Sylverster equation is given in the following
lemma.

###### Lemma 4.1

(Existence of a low rank approximant, Corollary 2 in grasedyck2004existence )
Let $A\in\mathbb{C}^{n\times n}$ and $B\in\mathbb{C}^{n\times n}$ be matrices
with spectrum $\sigma(A)$ and $\sigma(B)$ as in the rectangular case or the
triangular case with constant $\mu$ ($\mu\geq 1$ in the triangular case),
$\Gamma_{A}$ and $\Gamma_{B}$ are paths of index 1 around the spectrum of $A$
and $B$ with $\Lambda$, $\lambda$, paths $\Gamma_{A}$, $\Gamma_{B}$ and a
partitioning $\Gamma_{B}=\dot{\bigcup}_{j=0}^{k_{\sigma}-1}\Gamma_{B,j}$. For
all $j=0,1,2,\ldots,k_{\sigma}-1$, let $\eta_{j}$ be the center of a part
$\Gamma_{B,j}$, $\xi\in\Gamma_{A}$ and $\eta\in\Gamma_{B,j}$. We define

$g_{i,j}(\xi):=(\xi-\eta_{j})^{-1-i}\quad\text{and}\quad
h_{i,j}(\eta):=(\eta-\eta_{j})^{i},\quad 0\leq i<k,\quad 0\leq j<k_{\sigma},$

and

$\kappa_{A}:=\frac{1}{2\pi}\oint_{\Gamma_{A}}\|(\xi
I-A)^{-1}\|_{F}d\xi,\quad\kappa_{B}:=\frac{1}{2\pi}\oint_{\Gamma_{B}}\|(\eta
I-B)^{-1}\|_{F}d\eta,$
$k_{\varepsilon}:=\Big{\lceil}\log_{2}\big{(}1+\frac{6(\|A\|_{F}+\|B\|_{F})\kappa_{A}\kappa_{B}}{\lambda\varepsilon}\big{)}\Big{\rceil},$

where $\varepsilon\in(0,1)$ is given. Then for each right-hand side matrix
$R\in\mathbb{C}^{n\times m}$ of rank $r_{R}$, the matrix

$\tilde{X}:=\frac{1}{4\pi^{2}}\sum_{j=0}^{k_{\sigma}-1}\sum_{j=0}^{k_{\varepsilon}-1}\big{(}\oint_{\Gamma_{A}}(\xi
I-A)^{-1}g_{i,j}(\xi)d\xi\big{)}R\big{(}\oint_{\Gamma_{B}}(\eta
I-B)^{-1}h_{i,j}(\eta)d\eta\big{)}$

approximates the solution $X$ to (16) by

$\|X-\tilde{X}\|_{F}\leq\varepsilon\|X\|_{F}.$

Besides, the rank of $\tilde{X}$ is bounded by
$r_{R}k_{\sigma}k_{\varepsilon}$ with

$k_{\sigma}=\left\\{\begin{array}[]{lll}&\mathcal{O}(\frac{\mu}{\lambda}+\log_{2}(2+\frac{\Lambda}{\lambda}))&(rectangular\
case),\\\
&\mathcal{O}(\mu\log_{2}(2+\frac{\Lambda}{\lambda}))&(triangularcase).\end{array}\right.$

Lemma 4.1 gives the existence of the low-rank solution to the equation (16).
However, due to the complexity of the problems, it is difficult to estimate
the exact rank of the discretization system from the eddy current optimal
problem. The bound of the rank can be obtained from Lemma 4.1 and it is
present in the following theorem.

###### Theorem 4.1

Given the eddy current optimization problems (1)-(5). Let $h$ and $\tau$ be
the spatial step and temporal step, respectively. By making use of the all-at-
once based approach, we have the discretized system (13), which is equivalent
to the matrix equation (16). Suppose that the right-hand side matrix $R$ is
given in the low-rank form (18) with $r_{R}$ being the rank, and the solution
is approximated in the form $\tilde{X}=V_{p}Z(\approx X)$ up to an accuracy
$\varepsilon$. Let $M$ and $K$ be the mass matrix and the stiffness matrix
described previously. Then the rank $r_{X}$ of the solution $\tilde{X}$ is
bounded by

$\bar{r}_{X}=\mathcal{O}((\log\frac{1}{\varepsilon}+\log(\frac{1}{h^{2}}+\frac{1}{\tau}))^{2}r_{R}).$

###### Proof

We can rewrite the resulting matrix equation (16) into a equivalent linear
equation form as

$(I_{n}\otimes M^{-1}K+B\otimes I_{2m_{T}})\underline{X}=\underline{R},$

where the underline of a matrix denotes its vectorization. Denote by
$\mathcal{A}=I_{n}\otimes M^{-1}K+B\otimes I_{2m_{T}}$, then by using the
results in grasedyck2004existence ; dolgov2017low , we can approximate the
inverse of $\mathcal{A}$ in the low-rank form by the exponential quadrature,
i.e., for given $N$ and $k$, denote by $t_{k}=exp(\frac{k\pi}{\sqrt{N}})$,
$c_{x}=\frac{t_{k}\pi}{\sqrt{N}}$, then

$\mathcal{A}^{-1}=(I_{n}\otimes M^{-1}K+B\otimes
I_{2m_{T}})^{-1}\approx\Sigma_{k=-N}^{N}c_{k}\exp(-t_{k}B)\otimes\exp(-t_{k}M^{-1}K),$

where the accuracy is estimated by $\mathcal{O}(\|\mathcal{A}\|_{2}\cdot
e^{-\pi\sqrt{2N}})$, provided that $\|\mathcal{A}^{-1}\|=\mathcal{O}(1)$.

Therefore, the rank of $\mathcal{A}^{-1}$ is estimated by
$\mathcal{O}((\log\frac{1}{\varepsilon}+\log\text{cond}\mathcal{A})^{2})$.
Moreover, as $\text{cond}(\mathcal{A})=\mathcal{O}(h^{-2}+\tau^{-1})$, then
the rank of $\tilde{X}$ is bounded by

$\mathcal{O}((\log\frac{1}{\varepsilon}+\log(\frac{1}{h^{2}}+\frac{1}{\tau}))^{2}r_{R}).$

##  5 Numerical Experiments

We now present the performance and flexibility of the splitting-based KPIK
method (denoted by “ SKPIK ”) on two examples for the eddy current constrained
optimization problems. All experiments were run on a desktop computer with an
Intel(R) Core(TM) i7-10710U CPU @ 1.10GHz 1.61 GHz with 16 GB of RAM. The
lowest order linear Nédélec edge element, which is defined on 2D triangles
nedelec1980mixed  and 3D tetrahedra nedelec1986new , are used to discretize
the state variable, the control variable, and the adjoint variable in our
experiments. The MATLAB package of anjam2015fast  is used to construct the
relevant matrices. To show numerically the feasibility and effectiveness of
the new method, we report the results of all the proposed methods in the sense
of iteration step (denoted as “IT”), elapsed CPU time in seconds (denoted as
“CPU”), and relative residual error (denoted as “RES”). Regarding the low-rank
approximation methods, we report the rank (denoted as “r”) of the final
solutions additionally.

We test all the methods for two examples, i.e., a 2D problem and a 3D problem.
In both examples, the mass matrix $M$ is symmetric positive definite.
Therefore, the coefficient matrix $\mathcal{A}$ of matrix-vector form of the
linear equations (11):

$\mathcal{A}x:=\begin{pmatrix}\tau\mathcal{M}&0&\mathcal{N}_{\sigma}^{T}\\\
0&\tau\beta\mathcal{M}&-\tau\mathcal{M}\\\
\mathcal{N}_{\sigma}&-\tau\mathcal{M}&0\end{pmatrix}\begin{pmatrix}\underline{Y}\\\
\underline{U}\\\
\underline{P}\end{pmatrix}=\begin{pmatrix}\tau\mathcal{M}\underline{Y}_{d}\\\
0\\\ 0\end{pmatrix}=:b$ (21)

is symmetric, then we can consider the low-rank approximate combining with the
MINRES method stoll2015low . To implement the low-rank MINRES method
efficiently, one always apply an appropriate preconditioner, e.g., the block-
diagonal preconditioner with the Schur complement as

$\mathcal{P}=\begin{pmatrix}\tau\mathcal{M}&&\\\ &\tau\beta\mathcal{M}&\\\
&&\mathcal{S}\end{pmatrix},$

where the Schur complement
$\mathcal{S}=\frac{1}{\tau}\mathcal{N}_{\sigma}^{T}\mathcal{M}^{-1}\mathcal{N}_{\sigma}+\frac{\tau}{\beta}\mathcal{M}$.
However, the main difficulty to use the above preconditioner is the heavy
workload of solving the sub-linear system with the coefficient matrix being
the Schur complement $\mathcal{S}$. Hence, we will compare the preconditioner
using the approximation of the Schur-complement
$\hat{\mathcal{S}}\approx\mathcal{S}$, which is presented by Pearson and co-
authors in pearson2012new ; pearson2012regularization , as

$\widehat{\mathcal{P}}=\begin{pmatrix}\tau\mathcal{M}&&\\\
&\tau\beta\mathcal{M}&\\\ &&\hat{\mathcal{S}}\end{pmatrix},$ (22)

with
$\hat{\mathcal{S}}=\frac{1}{\tau}(\mathcal{N}_{\sigma}+\frac{\tau}{\sqrt{\beta}}\mathcal{M})\mathcal{M}^{-1}(\mathcal{N}_{\sigma}+\frac{\tau}{\sqrt{\beta}}\mathcal{M})^{T}$.
Therefore, we compare our new method with the low-rank MINRES method with the
preconditioner $\widehat{\mathcal{P}}$ in (22) for solving the discretizetion
system (21) by all-at-once approach, namely “ LRMINRES ”. Details about the
LRMINRES algorithm can be found in stoll2015low .

We want to point out that when we need to solve the sublinear systems with the
coefficient matrix being sparse and symmetric positive definite, e.g., $M$,
$K$, and so on, we use the sparse Cholesky factorization incorporated with the
symmetric approximation minimum degree reordering. To do so, we use symamd.m
command in the MATLAB toolbox. The same way to use the above command can be
found in bai2003hermitian  and it has been proven to be very effective in
solving large sparse symmetric positive definite linear equations.

During the comparison, we also used the non-all-at-once approach to discretize
the eddy current optimization problem. Then the full-rank MINRES algorithm
paige1975solution  is used to solve the resulting discretization systems step-
by-step. Or equivalently, at each time step, we use the MINRES algorithm with
the block-diagonal preconditioner $P$, which is defined by

$P=\begin{pmatrix}\tau M&&\\\ &\tau\beta M&\\\ &&\hat{S}\end{pmatrix},$

where
$\hat{S}=\frac{1}{\tau}(K+\frac{\tau}{\sqrt{\beta}}M)M^{-1}(K+\frac{\tau}{\sqrt{\beta}}M)$
is approximation to the Schur complement
$S=\frac{1}{\tau}KM^{-1}K+\frac{\tau}{\beta}M$. The corresponding algorithm is
denoted as “ FMINRES ” in our experiments.

The tolerance for all methods is set to be $10^{-6}$. That is to say, all
iteration processes are terminated when the current relative residuals satisfy

$\text{RES}:=\frac{\|b-\mathcal{A}x^{(k)}\|_{2}}{\|b\|_{2}}\leq 10^{-6},$

or the methods do not reach convergence when the numbers of iteration steps
reach the maximum number $\text{IT}_{\max}=500$, or the computing time of the
corresponding method exceeds 1000 in seconds. $x^{(0)}=0$ is the initial guess
and $x^{(k)}$ is the $k$th iterates of the corresponding iteration processes,
respectively. If the current iterates can not reach the above tolerance within
the given maximum iteration number or the above limited time, we denote the
results as ‘-’ in the following tables.

To find the low-rank representation $X=V_{p}Z$ of the solution to the matrix
equation (16), we perform the skinny QR factorization of both matrices, i.e.,
$V_{p}$ and $Z$. The rank of the low-rank approximation solution is produced
by dropping small singular values (depending on some tolerance, denoted as
“truncation tolerance” in our experiments) and then it leads to a low-rank
approximation. The MATLAB function svds is used directly to compute the
singular value decomposition (SVD) of $V_{p}Z$. According to stoll2015low , we
use the truncation tolerance $10^{-10}$ throughout our numerical experiments.
More details about the way to compute the truncated SVD can be found in
stoll2015low .

It is worth mentioning that the iteration counts of the SKPIK method we report
in all the tables refer to the number of iteration counts of the extended
Krylov subspace method, i.e., the KPIK algorithm simoncini2007new . While the
number of iteration counts of the LRMINRES method refers to the number of
iteration counts of the MINRES algorithm. The number of the iteration counts
of the FMINRES method refers to the average of the numbers of iteration
counts, i.e., the sum of the required iterations on all-time steps divided by
the number of time steps $m_{T}$. Besides, when we solve the systems regarding
the coefficient matrices being the Schur-complement approximation in the
LRMINRES method, we employ the inexact KPIK with a fixed number of steps,
i.e., 6 steps according to stoll2015low .

###### Example 1

We consider $\Omega=[0,1]^{2}$ and $\nu=1$. We split the domain into two parts
across the diagonal with $\Omega_{1}=\\{x\in\Omega|x_{1}>x_{2}\\}$ and
$\Omega_{2}=\Omega\setminus\Omega_{1}$. The desired state is given by

$y_{d}(x,t)|_{\Omega_{1}}=\begin{pmatrix}\sin(2\pi x_{1})+2\pi\cos(2\pi
x_{1})(x_{1}-x_{2})\\\ \sin((x_{1}-x_{2})^{2}(x_{1}-1)^{2}x_{2}-\sin(2\pi
x_{1}))\end{pmatrix}\quad\text{and}\quad y_{d}(x,t)|_{\Omega_{2}}=0.$

In this example, our experiments are performed for the final time $T=1$ by
varying the numbers of time steps as $m_{T}=100,200,400,800,1600$ and $3200$.
We first test all the methods regarding the numbers of spatial discretization
nodes being $n=3136$, as well as a large range of problem parameters
concerning $\beta=10^{-2}$, $10^{-4}$, $10^{-6}$ and $10^{-8}$,
$\sigma=10^{-5}$, $10^{-3}$, $10^{-1}$, $10^{1}$, $10^{3}$ and $10^{5}$.

Firstly, we test a small number of time steps and a small size of space
discretization, i.e., $m_{T}=100$ and $n=3136$. The results are listed in
Figure 1, concerning different choices of the parameters $\beta$ and $\sigma$,
where on the left is the SKPIK method, in the middle is the LRMINRES method,
and on the right is the FMINRES method.

Figure 1: The computing time against different values of $\sigma$ and $\beta$
for all the proposed methods ($m_{T}=100$).

We further test all the methods by varying the numbers of time steps as
$m_{T}=$100, 200, 400, 800, 1600, 3200 for the size of space discretization
$n=3136$. Then we plot the computing time against the number of time steps in
Fig. 2 by fixing $\sigma=10$, where again on the left is the SKPIK method, in
the middle is the LRMINRES method, and on the right is the FMINRES method.

Figure 2: The computing time against different time steps $m_{T}$ for all the
proposed methods.

From Figures 1-2, we see the FMINRES method can solve the optimization problem
successfully for small dimensional problems because of the efficiency of the
diagonal-block preconditioning technique and the MINRES algorithm. However, we
also see that the FMINRES method has always consumed lots of time to obtain
the solution, because the FMINRES method needs to obtain the solution at each
time step before solving the system of the next time step, and also needs to
prepare for the coefficient matrices and the right-hand side vectors for the
next time step. Besides, from Fig. 2, we also find that the computing time
needed by the FMINRES method increases rapidly as the number of time steps
increases.

Meanwhile, the LRMINRES method seems to be more efficient than the FMINRES
method for the small number of time steps and the small size of space
discretization in Figure 1. Particularly, when $m_{T}=100$ and 200, and the
parameters $\beta$ and $\sigma$ are small enough, the LRMINRES method performs
very well. But when $m_{T}=400$, 800, 1600, and 3200, the iteration numbers to
get the solution needed by the LRMINRES method vary greatly corresponding to
different problem parameters $\beta$ and $\sigma$. These results are
consistent with the results obtained by Stoll and co-authors stoll2015low .

From all the results plotted in these figures, we see the SKPIK method always
performs the best. It needs the least computing time. Besides, as the number
of time steps increases, the computing time of the SKPIK method increases very
gently. To illustrate this, we further plot in Fig. 3 the computing time of
the SKPIK method against the number of time steps for the number of spatial
discretization nodes being $n=12416$.

Figure 3: The computing time of the SKPIK method against different time steps
$m_{T}$ ($\sigma=10^{-1}$(left), $\sigma=10^{1}$(right)).

Therefore, as seen in Figures 1-3, increasing the discretization size barely
impacts the resulting subspace sizes. The computing time needed by the FMINRES
method increases rapidly as the number of time steps increases. Compared with
the FMINRES method, the LRMINRES method seems to be less impacted concerning
the number of time steps. However, the LRMINRES method depends greatly on the
choices of the parameters $\beta$, $\sigma$ as well as the number of time
steps.

From all the above results, we can see that the SKPIK method shows great
robustness concerning the discretization sizes $n$ and $m_{T}$ as well as the
control parameter $\beta$ and the conductivity $\sigma$. Additionally, the
time needed by the SKPIK method to solve the optimization problems increases
considerably slowly.

To show more details about the computing results of the SKPIK method, at the
end of this example, we test higher dimensional problems, i.e., $n=49408$,
$n=197120$ and $m_{T}=800,1600,3200$, which roughly resembles up to a total of
78 million degrees of freedom. The numerical results are listed in Table 1 and
Table 2, respectively.

Table 1: Numerical results of the SKPIK method for Example 1 (2D). $n=49408$ |  $m_{T}$=800 |  $m_{T}$=1600 |  $m_{T}$=3200
---|---|---|---
$\sigma,\beta$ | r | IT | CPU | RES | r | IT | CPU | RES | r | IT | CPU | RES
$\sigma=10^{-4}$ |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | 6 | 49 | 2.68 | 8.67e-7 | 6 | 51 | 2.92 | 1.63e-7 | 6 | 51 | 3.93 | 1.63e-7
$10^{-4}$ | 6 | 93 | 7.66 | 5.46e-7 | 6 | 93 | 10.09 | 5.46e-7 | 6 | 93 | 14.57 | 5.46e-7
$10^{-6}$ | 4 | 113 | 13.36 | 9.59e-7 | 4 | 113 | 20.78 | 9.59e-7 | 4 | 113 | 21.46 | 9.59e-7
$10^{-8}$ | 4 | 51 | 2.71 | 6.91e-7 | 4 | 51 | 2.92 | 6.91e-7 | 4 | 51 | 3.93 | 6.91e-7
$\sigma=1$ |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | 6 | 87 | 7.16 | 8.29e-7 | 6 | 93 | 11.21 | 5.08e-7 | 6 | 93 | 14.59 | 6.96e-7
$10^{-4}$ | 6 | 93 | 7.66 | 6.87e-7 | 6 | 93 | 11.26 | 8.41e-7 | 6 | 93 | 14.47 | 9.94e-7
$10^{-6}$ | 6 | 113 | 13.39 | 9.56e-7 | 6 | 113 | 20.06 | 9.56e-7 | 6 | 113 | 21.47 | 9.56e-7
$10^{-8}$ | 6 | 51 | 2.71 | 6.91e-7 | 6 | 51 | 2.92 | 6.91e-7 | 6 | 51 | 3.67 | 6.91e-7
$\sigma=10^{4}$ |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | 6 | 21 | 0.63 | 8.57e-7 | 6 | 22 | 0.87 | 9.23e-7 | 6 | 23 | 1.10 | 7.29e-7
$10^{-4}$ | 6 | 22 | 0.64 | 5.99e-7 | 6 | 23 | 1.00 | 6.10e-7 | 6 | 24 | 1.23 | 5.38e-7
$10^{-6}$ | 6 | 22 | 0.64 | 9.32e-7 | 6 | 23 | 1.01 | 9.92e-7 | 6 | 24 | 1.22 | 8.40e-7
$10^{-8}$ | 6 | 24 | 0.75 | 4.93e-7 | 6 | 24 | 1.09 | 8.10e-7 | 6 | 25 | 1.70 | 6.54e-7
Table 2: Numerical results of the SKPIK method for Example 1 (2D). $n=197120$ |  $m_{T}$=800 |  $m_{T}$=1600 |  $m_{T}$=3200
---|---|---|---
$\sigma,\beta$ | r | IT | CPU | RES | r | IT | CPU | RES | r | IT | CPU | RES
$\sigma=10^{-4}$ |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | 6 | 51 | 15.53 | 3.12e-7 | 6 | 51 | 17.47 | 3.13e-7 | 6 | 51 | 22.03 | 3.13e-7
$10^{-4}$ | 6 | 99 | 66.35 | 6.78e-7 | 6 | 99 | 68.28 | 9.39e-7 | 6 | 99 | 103.79 | 9.72e-7
$10^{-6}$ | 4 | 161 | 163.50 | 6.29e-7 | 4 | 161 | 170.83 | 6.29e-7 | 4 | 161 | 203.03 | 6.29e-7
$10^{-8}$ | 4 | 99 | 66.01 | 7.03e-7 | 4 | 99 | 68.90 | 7.03e-7 | 4 | 99 | 76.90 | 7.03e-7
$\sigma=1$ |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | 6 | 111 | 75.39 | 5.94e-7 | 6 | 117 | 86.98 | 5.49e-7 | 6 | 117 | 96.46 | 9.25e-7
$10^{-4}$ | 6 | 111 | 75.41 | 6.78e-7 | 6 | 117 | 86.42 | 5.91e-7 | 6 | 117 | 96.79 | 9.72e-7
$10^{-6}$ | 6 | 161 | 163.88 | 6.27e-7 | 6 | 161 | 170.03 | 6.28e-7 | 6 | 161 | 204.79 | 6.28e-7
$10^{-8}$ | 6 | 99 | 66.75 | 7.03e-7 | 6 | 99 | 69.90 | 7.03e-7 | 6 | 99 | 72.39 | 7.02e-7
$\sigma=10^{4}$ |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | 6 | 34 | 9.65 | 6.97e-7 | 6 | 36 | 10.34 | 7.72e-7 | 6 | 40 | 12.97 | 8.73e-7
$10^{-4}$ | 6 | 34 | 9.01 | 7.72e-7 | 6 | 36 | 12.16 | 9.69e-7 | 6 | 41 | 14.81 | 7.45e-7
$10^{-6}$ | 6 | 38 | 11.56 | 7.21e-7 | 6 | 38 | 12.99 | 8.72e-7 | 6 | 41 | 14.73 | 9.77e-7
$10^{-8}$ | 6 | 47 | 14.62 | 8.23e-7 | 6 | 47 | 17.64 | 8.23e-7 | 6 | 47 | 18.83 | 8.26e-7

The results in Tables 1-2 once again illustrate that the SKPIK method is very
efficient and shows great robustness concerning different problem sizes and
problem parameters.

From the numerical results in Example 1, we see that the SKPIK method
outperforms the other methods especially when the number of time steps is
large. Hence, in the following example, we further test all the methods by a
3D problem to show the robustness of the SKPIK method.

###### Example 2

We consider $\Omega=[0,1]^{3}$ and $\nu=1$. The desired state is given by

$y_{d}(x,t)=\begin{pmatrix}0\\\ 0\\\ \sin(\pi x_{1})\sin(\pi x_{2})\sin(\pi
x_{3})\end{pmatrix}.$

In this example, we test all the methods by varying the number of time steps
as $m_{T}=800,1600,3200$, the cost parameters
$\beta=10^{-2},10^{-4},10^{-6},10^{-8}$ and the conductivity
$\sigma=10^{-4},1,10^{4}$. The number of spatial discretization nodes varies
from $n=1854$, $n=13428$ to $n=102024$. The numerical results corresponding to
all the methods are listed in Table 3 ($m_{T}=800$), Table 4 ($m_{T}=1600$)
and Table 5 ($m_{T}=3200$).

The numerical results in Tables 3-5 show much analogous computational
phenomenon to our previous observations. Therefore, we can conclude that our
new method is very robust and effective for solving large-scale systems from
an all-at-once approach to the discretized eddy current optimization problem.

Table 3: $m_{T}=800$ for Example 2 (3D).  $n$ |  | n=1854 | n=13428 | n=102024
---|---|---|---|---
$\sigma,\beta$ | method | r | IT | CPU | RES | r | IT | CPU | RES | r | IT | CPU | RES
$\sigma=10^{-4}$ |  |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | SKPIK | 6 | 42 | 0.35 | 5.12e-7 | 6 | 68 | 2.06 | 7.66e-7 | 6 | 96 | 52.37 | 9.89e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 23 | 352.17 | 3.18e-7 |  | - | - | - |  | - | - | -
$10^{-4}$ | SKPIK | 6 | 30 | 0.17 | 6.17e-7 | 6 | 58 | 1.49 | 8.48e-7 | 6 | 107 | 55.81 | 8.30e-7
| LRMINRES | 6 | 175 | 41.15 | 9.32e-7 | - | - | - | - | - | - | - | -
| FMINRES |  | 22 | 340.65 | 9.51e-7 |  | - | - | - |  | - | - | -
$10^{-6}$ | SKPIK | 4 | 10 | 0.04 | 4.13e-7 | 4 | 21 | 0.31 | 6.93e-7 | 4 | 43 | 11.80 | 8.30e-7
| LRMINRES | 6 | 133 | 44.95 | 9.72e-7 | 6 | 271 | 1614 | 9.70e-7 | - | - | - | -
| FMINRES |  | 23 | 669.99 | 3.49e-7 |  | - | - | - |  | - | - | -
$10^{-8}$ | SKPIK | 3 | 4 | 0.02 | 4.44e-7 | 4 | 7 | 0.12 | 1.73e-7 | 4 | 13 | 2.55 | 5.99e-7
| LRMINRES | 6 | 27 | 7.81 | 6.66e-7 | 6 | 49 | 326.95 | 6.30e-7 | 6 | 59 | 22639 | 9.65e-7
| FMINRES |  | 15 | 461.66 | 9.87e-8 |  | - | - | - |  | - | - | -
$\sigma=1$ |  |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | SKPIK | 6 | 42 | 0.35 | 5.10e-7 | 6 | 68 | 2.06 | 9.18e-7 | 6 | 104 | 53.88 | 9.92e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 7 | 118.37 | 2.27e-7 |  | - | - | - |  | - | - | -
$10^{-4}$ | SKPIK | 6 | 30 | 0.17 | 5.98e-7 | 6 | 58 | 1.42 | 8.24e-7 | 6 | 107 | 55.64 | 8.09e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 12.59 | 203.14 | 1.32e-7 |  | - | - | - |  | - | - | -
$10^{-6}$ | SKPIK | 4 | 10 | 0.04 | 4.12e-7 | 6 | 21 | 0.31 | 6.91e-7 | 6 | 43 | 11.98 | 8.27e-7
| LRMINRES | 6 | 243 | 100.67 | 9.87e-7 | 6 | 335 | 2416 | 9.83e-7 | - | - | - | -
| FMINRES |  | 9 | 279.38 | 4.28e-7 |  | - | - | - |  | - | - | -
$10^{-8}$ | SKPIK | 4 | 4 | 0.02 | 4.43e-7 | 6 | 7 | 0.11 | 1.73e-7 | 6 | 13 | 3.55 | 5.98e-7
| LRMINRES | 6 | 29 | 8.41 | 5.85e-7 | 6 | 63 | 469.00 | 9.38e-7 | 6 | 67 | 27664 | 8.61e-7
| FMINRES |  | 11 | 358.81 | 1.28e-7 |  | - | - | - |  | - | - | -
$\sigma=10^{4}$ |  |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | SKPIK | 4 | 4 | 0.02 | 2.38e-8 | 5 | 5 | 0.08 | 9.51e-8 | 6 | 7 | 2.08 | 5.95e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 2 | 83.79 | 2.25e-11 |  | - | - | - |  | - | - | -
$10^{-4}$ | SKPIK | 4 | 4 | 0.02 | 2.79e-8 | 5 | 5 | 0.08 | 2.46e-7 | 6 | 7 | 2.06 | 8.52e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 2 | 88.11 | 2.17e-8 |  | - | - | - |  | - | - | -
$10^{-6}$ | SKPIK | 4 | 4 | 0.02 | 9.13e-8 | 5 | 5 | 0.08 | 4.28e-7 | 6 | 8 | 2.11 | 1.68e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 2 | 87.57 | 5.86e-8 |  | - | - | - |  | - | - | -
$10^{-8}$ | SKPIK | 4 | 4 | 0.02 | 1.10e-7 | 5 | 5 | 0.08 | 8.60e-7 | 6 | 8 | 2.10 | 3.61e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 3 | 115.24 | 3.77e-8 |  | - | - | - |  | - | - | -
Table 4: $m_{T}=1600$ for Example 2 (3D).  $n$ |  | n=1854 | n=13428 | n=102024
---|---|---|---|---
$\sigma,\beta$ | method | r | IT | CPU | RES | r | IT | CPU | RES | r | IT | CPU | RES
$\sigma=10^{-4}$ |  |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | SKPIK | 6 | 42 | 0.45 | 5.12e-7 | 6 | 68 | 2.25 | 7.66e-7 | 6 | 96 | 54.99 | 9.90e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 23 | 773.37 | 3.18e-8 |  | - | - | - |  | - | - | -
$10^{-4}$ | SKPIK | 6 | 30 | 0.24 | 6.17e-7 | 6 | 58 | 1.68 | 8.48e-7 | 6 | 107 | 58.82 | 8.03e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 23 | 792.08 | 1.81e-7 |  | - | - | - |  | - | - | -
$10^{-6}$ | SKPIK | 4 | 10 | 0.07 | 4.13e-7 | 4 | 21 | 0.40 | 6.93e-7 | 4 | 43 | 11.17 | 8.30e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 23 | - | 3.46e-7 |  | - | - | - |  | - | - | -
$10^{-8}$ | SKPIK | 3 | 4 | 0.03 | 4.44e-7 | 4 | 7 | 0.12 | 1.73e-7 | 4 | 13 | 3.88 | 5.99e-7
| LRMINRES | 6 | 105 | 67.74 | 9.55e-7 | 6 | 89 | 600.37 | 8.88e-7 | - | - | - | -
| FMINRES |  | 15 | - | 9.80e-8 |  | - | - | - |  | - | - | -
$\sigma=1$ |  |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | SKPIK | 6 | 45 | 0.50 | 8.14e-7 | 6 | 68 | 2.25 | 9.27e-7 | 6 | 105 | 56.35 | 9.39e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 7 | 266.45 | 7.22e-9 |  | - | - | - |  | - | - | -
$10^{-4}$ | SKPIK | 6 | 30 | 0.24 | 5.99e-7 | 6 | 58 | 1.66 | 8.24e-7 | 6 | 107 | 58.52 | 8.09e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 10.10 | 374.82 | 1.71e-7 |  | - | - | - |  | - | - | -
$10^{-6}$ | SKPIK | 6 | 10 | 0.07 | 4.12e-7 | 6 | 21 | 0.40 | 6.91e-7 | 6 | 43 | 11.17 | 8.28e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 11 | 698.08 | 1.98e-7 |  | - | - | - |  | - | - | -
$10^{-8}$ | SKPIK | 4 | 4 | 0.03 | 4.43e-7 | 6 | 7 | 0.12 | 1.73e-7 | 6 | 13 | 3.88 | 5.98e-7
| LRMINRES | 6 | 492 | 243.56 | 9.83e-7 | - | - | - | - | - | - | - | -
| FMINRES |  | 10 | 705.80 | 7.47e-7 |  | - | - | - |  | - | - | -
$\sigma=10^{4}$ |  |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | SKPIK | 4 | 4 | 0.04 | 3.34e-8 | 5 | 5 | 0.09 | 1.32e-7 | 6 | 7 | 1.99 | 8.62e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 2 | 170.62 | 5.51e-8 |  | - | - | - |  | - | - | -
$10^{-4}$ | SKPIK | 4 | 4 | 0.04 | 3.93e-8 | 5 | 5 | 0.09 | 3.49e-7 | 6 | 8 | 2.53 | 1.36e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 2 | 174.64 | 2.00e-8 |  | - | - | - |  | - | - | -
$10^{-6}$ | SKPIK | 4 | 4 | 0.03 | 1.31e-7 | 5 | 5 | 0.09 | 6.12e-7 | 6 | 8 | 2.53 | 2.45e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 2 | 176.18 | 2.61e-8 |  | - | - | - |  | - | - | -
$10^{-8}$ | SKPIK | 4 | 4 | 0.03 | 1.57e-7 | 6 | 6 | 0.11 | 3.07e-8 | 6 | 8 | 2.53 | 5.21e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 2 | 196.44 | 4.44e-7 |  | - | - | - |  | - | - | -
Table 5: $m_{T}=3200$ for Example 2 (3D).  $n$ |  | n=1854 | n=13428 | n=102024
---|---|---|---|---
$\sigma,\beta$ | method | r | IT | CPU | RES | r | IT | CPU | RES | r | IT | CPU | RES
$\sigma=10^{-4}$ |  |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | SKPIK | 6 | 42 | 0.83 | 5.12e-7 | 6 | 68 | 2.94 | 7.66e-7 | 6 | 97 | 57.26 | 9.34e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | - | - | - |  | - | - | - |  | - | - | -
$10^{-4}$ | SKPIK | 6 | 30 | 0.46 | 6.17e-7 | 6 | 58 | 2.63 | 8.48e-7 | 6 | 107 | 63.37 | 8.30e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | - | - | - |  | - | - | - |  | - | - | -
$10^{-6}$ | SKPIK | 4 | 10 | 0.13 | 4.13e-7 | 4 | 21 | 0.54 | 4.95e-7 | 4 | 43 | 14.37 | 8.30e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | - | - | - |  | - | - | - |  | - | - | -
$10^{-8}$ | SKPIK | 3 | 4 | 0.05 | 4.44e-7 | 4 | 7 | 0.16 | 1.73e-7 | 4 | 13 | 4.13 | 5.99e-7
| LRMINRES | 6 | 69 | 36.95 | 8.79e-7 | 6 | 127 | 852.71 | 9.52e-7 | - | - | - | -
| FMINRES |  | - | - | - |  | - | - | - |  | - | - | -
$\sigma=1$ |  |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | SKPIK | 6 | 56 | 1.24 | 8.44e-7 | 6 | 68 | 2.94 | 9.32e-7 | 6 | 106 | 62.52 | 6.25e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 5 | 407.86 | 4.90e-7 |  | - | - | - |  | - | - | -
$10^{-4}$ | SKPIK | 6 | 31 | 0.48 | 6.82e-7 | 6 | 58 | 2.63 | 8.24e-7 | 6 | 107 | 63.35 | 8.10e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 7.84 | 625.07 | 2.45e-7 |  | - | - | - |  | - | - | -
$10^{-6}$ | SKPIK | 6 | 10 | 0.13 | 4.17e-7 | 6 | 21 | 0.54 | 6.91e-7 | 6 | 43 | 18.43 | 8.28e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | - | - | - |  | - | - | - |  | - | - | -
$10^{-8}$ | SKPIK | 5 | 4 | 0.05 | 4.43e-7 | 6 | 7 | 0.16 | 1.73e-7 | 6 | 13 | 4.13 | 5.98e-7
| LRMINRES | 6 | 75 | 42.05 | 9.49e-7 | 6 | 177 | 1205.91 | 9.40e-7 | - | - | - | -
| FMINRES |  | - | - | - |  | - | - | - |  | - | - | -
$\sigma=10^{4}$ |  |  |  |  |  |  |  |  |  |  |  |  |
$10^{-2}$ | SKPIK | 4 | 4 | 0.05 | 4.72e-8 | 5 | 5 | 0.13 | 1.86e-7 | 6 | 8 | 2.47 | 7.24e-8
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 2 | 370.89 | 1.36e-12 |  | - | - | - |  | - | - | -
$10^{-4}$ | SKPIK | 4 | 4 | 0.05 | 5.55e-8 | 5 | 5 | 0.13 | 4.95e-7 | 6 | 8 | 2.47 | 1.96e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 2 | 350.95 | 1.54e-9 |  | - | - | - |  | - | - | -
$10^{-6}$ | SKPIK | 4 | 4 | 0.05 | 1.84e-7 | 5 | 5 | 0.13 | 8.71e-7 | 6 | 8 | 2.47 | 3.52e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 2 | 374.63 | 2.15e-8 |  | - | - | - |  | - | - | -
$10^{-8}$ | SKPIK | 4 | 4 | 0.05 | 2.22e-7 | 6 | 6 | 0.15 | 4.37e-8 | 6 | 8 | 2.47 | 7.45e-7
| LRMINRES | - | - | - | - | - | - | - | - | - | - | - | -
| FMINRES |  | 2 | 362.29 | 1.70e-7 |  | - | - | - |  | - | - | -

##  6 Conclusion and Outlook

We construct a splitting-based KPIK method to solve the discretization from
the eddy current optimal control problems in an all-at-once-based approach.
The SKPIK method relies on the reformulation of a special splitting of the
coefficient matrix from the order-reduced KKT system into a matrix equation
and the KPIK algorithm for solving the matrix equation. Some properties of the
low-rank solution are proposed. Besides, the SKPIK method has been illustrated
to be robust concerning different spatial and temporal discretizations and
parameters by experiments.

However, the LRMINRES method in our experiments performs poorly because of the
bad Schur-complement approximation in the preconditioner. Hence, improving the
Schur-complement approximation would be one of the future work. Another future
work will focus on the low-rank approximation method for solving all-at-once
discretized eddy current optimal control problems with much more complicated
boundary conditions.

## Acknowledgments

This work was supported by the National Natural Science Foundation of China
(Nos. 12126344, 12126337 and 11901324) and the China Scholarship Council (No.
202308350044).

## Conflict of Interests

The authors declare that there is no conflict of interests regarding the
publication of this article.

## References

  * [1] I. Anjam and J. Valdman.  Fast MATLAB assembly of FEM matrices in 2D and 3D: Edge elements.  Applied Mathematics and Computation, 267:252–263, 2015.
  * [2] O. Axelsson and Z.-Z. Liang.  A note on preconditioning methods for time-periodic eddy current optimal control problems.  Journal of Computational and Applied Mathematics, 352:262–277, 2019\.
  * [3] O. Axelsson and D. Lukáš.  Preconditioners for time-harmonic optimal control eddy–current problems.  In International Conference on Large-Scale Scientific Computing, pages 47–54. Springer, 2017.
  * [4] O. Axelsson and D. Lukáš.  Preconditioning methods for eddy-current optimally controlled time-harmonic electromagnetic problems.  Journal of Numerical Mathematics, 27(1):1–21, 2019.
  * [5] F. Bachinger, U. Langer, and J. Schöberl.  Efficient solvers for nonlinear time-periodic eddy current problems.  Computing and Visualization in Science, 9(4):197–207, 2006.
  * [6] Z.-Z. Bai, G. H. Golub, and M. K. Ng.  Hermitian and skew-Hermitian splitting methods for non-Hermitian positive definite linear systems.  SIAM Journal on Matrix Analysis and Applications, 24(3):603–626, 2003.
  * [7] P. Benner and T. Breiten.  Low rank methods for a class of generalized Lyapunov equations and related issues.  Numerische Mathematik, 124(3):441–470, 2013.
  * [8] P. Benner, S. Dolgov, A. Onwunta, and M. Stoll.  Low-rank solvers for unsteady Stokes–Brinkman optimal control problem with random data.  Computer Methods in Applied Mechanics and Engineering, 304:26–54, 2016.
  * [9] P. Benner, J. Heiland, and S. W. Werner.  A low-rank solution method for Riccati equations with indefinite quadratic terms.  Numerical Algorithms, 92(2):1083–1103, 2023.
  * [10] P. Benner, A. Onwunta, and M. Stoll.  Low-rank solution of unsteady diffusion equations with stochastic coefficients.  SIAM/ASA Journal on Uncertainty Quantification, 3(1):622–649, 2015\.
  * [11] P. Benner, A. Onwunta, and M. Stoll.  Block-diagonal preconditioning for optimal control problems constrained by PDEs with uncertain inputs.  SIAM Journal on Matrix Analysis and Applications, 37(2):491–518, 2016.
  * [12] M. Benzi, G. H. Golub, and J. Liesen.  Numerical solution of saddle point problems.  Acta Numerica, 14:1–137, 2005.
  * [13] M. Benzi, E. Haber, and L. Taralli.  A preconditioning technique for a class of PDE–constrained optimization problems.  Advances in Computational Mathematics, 35:149–173, 2011.
  * [14] T. Breiten, V. Simoncini, and M. Stoll.  Fast iterative solvers for fractional differential equations.  Electronic Transactions on Numerical Analysis, 45:751–776, 2016\.
  * [15] A. Buffa, H. Ammari, and J.-C. Nédélec.  A justification of eddy currents model for the Maxwell equations.  SIAM Journal on Applied Mathematics, 60(5):1805–1823, 2000.
  * [16] A. Bünger, V. Simoncini, and M. Stoll.  A low-rank matrix equation method for solving PDE–constrained optimization problems.  SIAM Journal on Scientific Computing, 43(5):S637–S654, 2021.
  * [17] A. Cichocki, A.-H. Phan, Q. Zhao, N. Lee, I. Oseledets, M. Sugiyama, D. P. Mandic, et al.  Tensor networks for dimensionality reduction and large-scale optimization: Part 2 applications and future perspectives.  Foundations and Trends® in Machine Learning, 9(6):431–673, 2017.
  * [18] S. Dolgov, J. W. Pearson, D. V. Savostyanov, and M. Stoll.  Fast tensor product solvers for optimization problems with fractional differential equations as constraints.  Applied Mathematics and Computation, 273:604–623, 2016.
  * [19] S. Dolgov and M. Stoll.  Low-rank solution to an optimization problem constrained by the Navier–Stokes equations.  SIAM Journal on Scientific Computing, 39(1):A255–A280, 2017.
  * [20] R. D. Falgout, S. Friedhoff, T. V. Kolev, S. P. MacLachlan, and J. B. Schroder.  Parallel time integration with multigrid.  SIAM Journal on Scientific Computing, 36(6):C635–C661, 2014.
  * [21] M. J. Gander.  50 years of time parallel time integration.  In Multiple Shooting and Time Domain Decomposition Methods: MuS-TDD, Heidelberg, May 6-8, 2013, pages 69–113. Springer, 2015.
  * [22] M. J. Gander, L. Halpern, J. Ryan, and T. T. B. Tran.  A direct solver for time parallelization.  In Domain decomposition methods in science and engineering XXII, pages 491–499. Springer, 2016.
  * [23] M. J. Gander and M. Neumüller.  Analysis of a new space-time parallel multigrid algorithm for parabolic problems.  SIAM Journal on Scientific Computing, 38(4):A2173–A2208, 2016.
  * [24] M. J. Gander and S. Vandewalle.  Analysis of the parareal time-parallel time-integration method.  SIAM Journal on Scientific Computing, 29(2):556–578, 2007.
  * [25] A. Goddard and A. Wathen.  A note on parallel preconditioning for all-at-once evolutionary PDEs.  Electronic Transactions on Numerical Analysis, 51:135–150, 2019\.
  * [26] G. H. Golub and C. F. Van Loan.  Matrix Computations.  The Johns Hopkins UniversityPress,, Baltimore, MD.
  * [27] L. Grasedyck.  Existence of a low rank or $\mathcal{H}$-matrix approximant to the solution of a Sylvester equation.  Numerical Linear Algebra with Applications, 11(4):371–389, 2004\.
  * [28] E. Haber and U. M. Ascher.  Preconditioned all-at-once methods for large, sparse parameter estimation problems.  Inverse Problems, 17(6):1847, 2001.
  * [29] M. Hinze.  A hierarchical space-time solver for distributed control of the Stokes equation.  Univ., Department Mathematik, 2008.
  * [30] M. Hinze, M. Köster, and S. Turek.  A space–time multigrid method for optimal flow control.  Constrained Optimization and Optimal Control for Partial Differential Equations, pages 147–170, 2012.
  * [31] G. Horton and S. Vandewalle.  A space-time multigrid method for parabolic partial differential equations.  SIAM Journal on Scientific Computing, 16(4):848–864, 1995.
  * [32] M. Kolmbauer.  Efficient solvers for multiharmonic eddy current optimal control problems with various constraints and their analysis.  IMA Journal of Numerical Analysis, 33(3):1063–1094, 2013.
  * [33] M. Kolmbauer and U. Langer.  A frequency-robust solver for the time-harmonic eddy current problem.  In Scientific Computing in Electrical Engineering SCEE 2010, pages 97–105. Springer, 2011.
  * [34] M. Kolmbauer and U. Langer.  A robust preconditioned MinRes solver for distributed time-periodic eddy current optimal control problems.  SIAM Journal on Scientific Computing, 34(6):B785–B809, 2012.
  * [35] M. Kolmbauer and U. Langer.  A robust preconditioned MinRes solver for time-periodic eddy current problems.  Computational Methods in Applied Mathematics, 13(1):1–20, 2013\.
  * [36] F. Kwok and B. W. Ong.  Schwarz waveform relaxation with adaptive pipelining.  SIAM Journal on Scientific Computing, 41(1):A339–A364, 2019.
  * [37] J.-L. Lions, Y. Maday, and G. Turinici.  Résolution d’edp par un schéma en temps $<$pararéel$>$.  Comptes Rendus de l’Académie des Sciences-Series I-Mathematics, 332(7):661–668, 2001.
  * [38] Y. Maday and G. Turinici.  A parareal in time procedure for the control of partial differential equations.  Comptes Rendus Mathematique, 335(4):387–392, 2002.
  * [39] I. Markovsky.  Low rank approximation: algorithms, implementation, applications, volume 906.  Springer, 2012.
  * [40] E. McDonald, J. Pestana, and A. Wathen.  Preconditioning and iterative solution of all-at-once systems for evolutionary partial differential equations.  SIAM Journal on Scientific Computing, 40(2):A1012–A1033, 2018.
  * [41] J.-C. Nédélec.  Mixed finite elements in $\mathbb{R}^{3}$.  Numerische Mathematik, 35:315–341, 1980.
  * [42] J.-C. Nédélec.  A new family of mixed finite elements in $\mathbb{R}^{3}$.  Numerische Mathematik, 50:57–81, 1986.
  * [43] C. C. Paige and M. A. Saunders.  Solution of sparse indefinite systems of linear equations.  SIAM journal on numerical analysis, 12(4):617–629, 1975.
  * [44] J. W. Pearson, M. Stoll, and A. J. Wathen.  Regularization-robust preconditioners for time-dependent PDE-constrained optimization problems.  SIAM Journal on Matrix Analysis and Applications, 33(4):1126–1152, 2012.
  * [45] J. W. Pearson and A. J. Wathen.  A new approximation of the Schur complement in preconditioners for PDE-constrained optimization.  Numerical Linear Algebra with Applications, 19(5):816–829, 2012\.
  * [46] T. Rees, M. Stoll, and A. Wathen.  All-at-once preconditioning in PDE–constrained optimization.  Kybernetika, 46(2):341–360, 2010.
  * [47] S. D. Shank and V. Simoncini.  Krylov subspace methods for large-scale constrained Sylvester equations.  SIAM Journal on Matrix Analysis and Applications, 34(4):1448–1463, 2013.
  * [48] D. Silvester, H. Elman, and A. Wathen.  Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics.  Oxford University Press, United Kingdom, first edition, 2005.
  * [49] V. Simoncini.  A new iterative method for solving large-scale Lyapunov matrix equations.  SIAM Journal on Scientific Computing, 29(3):1268–1288, 2007.
  * [50] M. Stoll and T. Breiten.  A low–rank in time approach to PDE–constrained optimization.  SIAM Journal on Scientific Computing, 37(1):B1–B29, 2015.
  * [51] M. Stoll and A. Wathen.  All–at–once solution of time–dependent Stokes control.  Journal of Computational Physics, 232(1):498–515, 2013.
  * [52] M. Wolfmayr.  A posteriori error estimation for the optimal control of time-periodic eddy current problems.  Computational Methods in Applied Mathematics, 2023.
  * [53] S.-L. Wu.  Toward parallel coarse grid correction for the parareal algorithm.  SIAM Journal on Scientific Computing, 40(3):A1446–A1472, 2018.
  * [54] S.-L. Wu and T. Zhou.  Acceleration of the two-level MGRIT algorithm via the diagonalization technique.  SIAM Journal on Scientific Computing, 41(5):A3421–A3448, 2019.
  * [55] M.-L. Zeng.  Respectively scaled splitting iteration method for a class of block 4-by-4 linear systems from eddy current electromagnetic problems.  Japan Journal of Industrial and Applied Mathematics, 38:489–501, 2021.
  * [56] Y.-L. Zhao, X.-M. Gu, and A. Ostermann.  A preconditioning technique for an all-at-once system from Volterra subdiffusion equations with graded time steps.  Journal of Scientific Computing, 88(1):11, 2021.