11institutetext:  Academy of Mathmatics and Systems Science Chinese Academy of
Science, Beijing, China, 11email: zhanganpeng, fengxt,
<EMAIL_ADDRESS>
# Size optimization of CNOT circuits on NISQ

Anpeng Zhang  Xiutao Feng and Shengyuan Xu

###### Abstract

Quantum computers in practice today require strict memory constraints, where
2-qubit operations can only be performed between the qubits closest to each
other in a graph structure. So a quantum circuit must undergo a transformation
to the graph before it can be implemented. In this paper, we study the
optimization of the CNOT circuits on some noisy intermediate-scale
quantum(NISQ) devices. Compared with previous works, we decompose it into two
sub-problems: optimization with a given initial qubit distribution and
optimization without limitations of initial qubit distribution. We find that
most of the previous researches focused on the first sub-problem, and ignored
the influence of different distribution of qubits in the same topology
structure on the optimization results. In this paper, We take both sub-
problems into account and give some new optimization algorithms. In short, our
method is divided into two steps: matrix optimization and routing
optimization. We implement matrix optimization with the algorithm in [XZL+20]
and put forward a new heuristic algorithm with MILP method which can solve the
second step. We implement our algorithm on IBM20 and some other NISQ devices,
the results are better than most other methods in our experiment.

###### keywords:

CNOT circuits topologically-constrained graph structure MILP NISQ

##  1 Introduction

With the rapid development of quantum computing technology, quantum computer
has gradually become a reality. Now noisy intermediate-scale quantum (NISQ)
[Pre18] computers with 10-80 qubits have been made by some laboratories such
as IBM and Google[Nay19].

In current quantum computers, many physical implementations of quantum
computers such as superconducting quantum circuits [Nay19] [ROT+18] and ion
traps[BSK+12] [BHL+16] impose strict memory limits, where 2-qubit operations
(such as CNOT gates) can only be performed between the nearest qubits in a
lattice or graphical structure. For a given quantum circuit with no limitation
of memory operation, how to transform it into a circuit that can be realized
on a real quantum computer has become an urgent problem to be solved. The
traditional solution to this problem is to exchange qubits through the SWAP
gate operation between adjacent qubits to achieve the realization of two
quantum gates between two non-adjacent qubits[ZPW18] [HS18]. However, this
scheme is too costly, which usually increases the depth of the circuit and the
number of gates by 1.5 to 3 times. In fact, finding the exact optimal solution
appears intractable [HND17].

In 2019, Kissinger $et\ al$[KdG19] showed that a pure CNOT circuit can be
compiled directly on an NISQ without significantly increasing the number of
CNOT gates. Inspired by classic SLP problem, they transformed a CNOT circuit
(linear map) into a matrix and then performed Restricted Gauss elimination to
find its optimal implementation under the limit of some topological graphs.
Later Nash $et\ al$ proposed a similar algorithm which can implement any CNOT
circuit with a $4n^{2}$-size equivalent CNOT circuit on any connected
graph[NGM20]. In 2021, Xiaoming Sun $et\ al$ put forward a new algorithm wihch
gives a $2n^{2}$-size equivalent CNOT circuit for any CNOT circuit on any
connected graph $G$ and proved that the lower bound of this problem is
$O(n^{2}/log\delta)$, where $\delta$ denotes the maximum degree of vertices in
$G$[WHY+19].

It has shown that finding the minimum CNOT gates in implementation for a CNOT
circuit without restriction can correspond to a Shortest Linear Program (SLP)
problem, and the latter was proved to be an NP-hard problem [BMP13]. In 2020,
Xiang Zejun $et\ al$ presented a new heuristic algorithm to search the optimal
implementations of reasonable large matrices based on the s-Xor operation
[XZL+20]. Their basic idea is to find various matrix decompositions for a
given matrix and optimize those decompositions to get the best implementation.
In order to optimize matrix decompositions, they present several matrix
multiplication rules over $F_{2}$, which are proved to be very powerful in
reducing the implementation cost.

In this paper we mainly discuss the optimization of the CNOT circuit under
some topologically-constrained quantum memories. First we abstract it into a
concrete mathematical problem, and give a new heuristic algorithm. Our
algorithm includes two steps:
1) find the optimal implementation of the corresponding matrix of the CNOT
circuit by Xiang $et\ al$’s method;
2) optimize the qubit distribution according to the result of Step 1) by means
of the solver Gurobi.
For consistency, we use the same set of randomly-generated CNOT circuits on 9,
16, and 20 qubits, and list the results in Table 1. The specific algorithm and
implementation can be found on Github. It is seen that our results are better
than the QuilC [SCZ16], $t\ket{ket}$ compilers [CDD+19] and stenier tree
[KdG19] in most cases.

The paper is organized as follows. We start in Section 2 by reviewing some
basic knowledge of quantum computing. In Section 3, we abstract the
optimization of CNOT circuit into mathematical problems and detail our method,
including matrix optimization and distribution optimization of qubits. In
Section 4, we compared our experimental results with previous works and apply
our method in the quantum implementation of AES’s Mixcolumns.

##  2 Preliminaries

###  2.1 CNOT circuit

Classical computer circuits consist of wires and logic gates [NC02]. The wires
are used to carry information around the circuit, while the logic gates
perform manipulations of the information, converting it from one form to
another. For example, the logic gate Xor can transform bits $a$ and $b$ to
$a\oplus b$, where $a,b=0\ or\ 1$. But different from classical bit, the qubit
is in superposition:

$\ket{\psi}=\alpha\ket{0}+\beta\ket{1}$

And we can write it in a vector notation as

$\ket{\psi}=\begin{bmatrix}\alpha\\\ \beta\end{bmatrix}$

In order to realize XOR operation on quantum computers, we introduce a multi-
qubit gate, which is called the controlled-NOT or CNOT gate. This gate has two
input qubits, known as the control qubit and the target qubit, respectively.

$\Qcircuit@C=1.5em@R=1em{&\ket{A}\qquad\ctrl{1}\qw\qquad\ket{A}\\\
\ket{B}\qquad\targ\qw\qquad\ket{A\oplus B}\\\ }$ Figure 1: Controlled-NOT gate

Another way of describing the action of the CNOT gate is to give a matrix
representation. In fact, every qubit gate can correspond to a unitary matrix,
we use $U_{CN}$ to denote the CNOT gate.

$U_{CN}=\begin{bmatrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&0&1\\\ 0&0&1&0\end{bmatrix}$

One can easily check that $U_{CN}^{{\dagger}}U_{CN}=I$. Of course, there are
many interesting quantum gates other than the controlled-NOT. But here we only
focus on quantum circuit which is composed of CNOT gates.

Definition 1: A quantum circuit is called a CNOT circuit if it contains only
CNOT gates.

Figure 2: An example of CNOT circuit swapping two qubits, and an equivalent
schematic symbol notation

In fact, since the CNOT circuit is linear, it is easy to map it to an $n\times
n$ matrix over $GF(2)$, and each CNOT gate can correspond to an elementary row
operation of a matrix, such as $CNOT_{i,j}$ is actually adding the row $i$ of
the matrix to the row $j$.

###  2.2 CNOT basis transformations

In some practical quantum computers such as IBM-QX5 in figure 3, there is a
strict limit between the control bit and the target bit, and the target bit
cannot be used as the control bit for controlled operation, but if we only
consider CNOT circuits, we can ignore the limitations by the follow fact.

Figure 3: a quantum circuit interchanging controlled bit and target bit of a
CNOT gate

Thus, with four hamdamard gates, we implement the interchange of a controlled
bit and a target bit. All graphs we will mention in the rest of this paper are
undirected graphs.

Figure 4: The structure of IBM-QX5, where an arrow from qubit A to qubit B
means that A can act as control qubit of a CNOT gate with target qubit B

##  3 Optimization of the CNOT circuit

Noisy intermediate-scale quantum (NISQ) has many different architectures, such
as 3 × 3 and 4 × 4 square lattices, 16-qubit architectures of the IBM QX-5 and
Rigetti Aspen devices, and the 20-qubit IBM Q20 Tokyo architecture. For each
architecture, we can represent it as a graph $G(V,E)$, and CNOT operations are
typically only possible between pairs of adjacent systems. In other words, we
can only perform row operation between specific rows when doing Gauss-Jordan
elimination.

Definition 2: Given a graph $G(V,E)$ and the order $1\sim|V|$ of vertices in
$G$. $S(G)$ is a set of elementary matrices in which each element is an
elementary row operation $E(i+j)$ where $E(i+j)$ is the resulting matrix by
adding the $j$th row to the $i$th row of an identity matrix and labels of two
rows are the ends of an edge in $G$.

Then the optimization of CNOT circuits can be expressed as the following
questions.

Problem 1: Given a matrix $A$ over $GF(2)$ and $S(G)$ of a graph $G(V,E)$, how
to decompose $A$ into the product of elements in $S(G)$ with the least
elements.

One can easily check that different orders of $V$ of a same graph $G$ will
derive different $S(G)$ and lead to different matrix optimization
implementations. Naturally, we consider the following question.

Problem 2: Given matrix $A$ over $GF(2)$ and graph $G(V,E)$, find an order of
$V$ and some matrixes in $S(G)$ such that we can decompose $A$ into the
product of the least elements in $S(G)$.

Problem 2 is a generalization of Problem 1, and Steiner tree method in [KdG19]
is for Problem 1, without considering the distribution of qubits.
Unfortunately, since Problem 2 involves two critical factors including
distribution of qubits and matrix optimization implementation, and these two
factors are mutually restricted, the distribution of qubits will greatly
affect the results of matrix optimization implementation, and the results of
matrix optimization implementation will also affect the selection of qubit
position. Therefore, it is difficult to optimize the two aspects at the same
time. In this case, without considering the limitation of the topology
structure, we firstly optimize the matrix, and then optimize the qubit
distribution through the result of matrix optimization implementation.

###  3.1 Main methods

####  3.1.1 Matrix optimization implementation

For each matrix $A$ over $GF(2)$, we perform Gauss-Jordan elimination on $A$,
and each elementary row operation corresponds to some elementary matrix, we
will obtain $A_{1}A_{2}...A_{k}A=I$, where $A_{i}$ is an elementary matrice.
Reverse the order of $A_{i}$ we can represent $A$ as product of some
elementary matrices.
Since the decomposition of a matrix is not unique, there is a natural problem
that how to find the least elementary row operations to perform Gauss-Jordan
elimination. This problem corresponds to the Shortest Linear Program (SLP)
problem. As for this problem, we can follow some rules of reduction to find
the optimal implementation[XZL+20]. Here is an example.
Let

$A=\begin{bmatrix}1&1&0&1&1\\\ 0&0&1&1&0\\\ 1&0&1&0&1\\\ 1&1&0&1&0\\\
1&1&1&1&0\end{bmatrix}$

First we can represent A as product of some elementary matrices $E(i+j)$ and
$E(i\leftrightarrow j)$, where $E(i+j)$ is the resulting matrix by adding the
$j$th row to the $i$th row of an identity matrix, $E(i\leftrightarrow j)$ is
the resulting matrix by exchanging the $i$th and $j$th row of an identity
matrix.
Then find the optimal implementation according to the following rules.

R1 | $E(k+i)E(k+j)E(i+j)=E(i+j)E(k+i)$
---|---
R2 | $E(i+k)E(k+j)E(i+j)=E(k+j)E(i+k)$
R3 | $E(i+k)E(j+k)E(i+j)=E(i+j)E(j+k)$
R4 | $E(j+k)E(i+k)E(i+j)=E(i+j)E(j+k)$
R5 | $E(k+j)E(k+i)E(i+j)=E(i+j)E(k+i)$
R6 | $E(k+j)E(i+k)E(i+j)=E(i+k)E(k+j)$
R7 | $E(j+i)E(i+j)=E(i\leftrightarrow j)E(j+i)$

We get

$A=A^{\prime}\times E(2+4)\times E(4+3)\times E(2+1)\times E(1+5)\times
E(1+3)\times E(5+2)\times E(3+2)$

Where $A^{\prime}$ is a permutation matrix

$\begin{bmatrix}0&0&0&0&1\\\ 0&0&0&1&0\\\ 1&0&0&0&0\\\ 0&1&0&0&0\\\
0&0&1&0&0\end{bmatrix}$

Then we get a CNOT circuit of A in Figure 5.

$\Qcircuit@C=1.5em@R=1em{&a_{1}\qquad\qw\qw\targ\targ\targ\qw\qw\qw a_{3}\\\
a_{2}\qquad\targ\qw\ctrl{-1}\qw\qw\ctrl{3}\ctrl{1}\qw a_{4}\\\
a_{3}\qquad\qw\targ\qw\qw\ctrl{-2}\qw\targ\qw a_{5}\\\
a_{4}\qquad\ctrl{-2}\ctrl{-1}\qw\qw\qw\qw\qw\qw a_{2}\\\
a_{5}\qquad\qw\qw\qw\ctrl{-4}\qw\targ\qw\qw a_{1}\\\ }$ Figure 5: the CNOT
circuit of $A$ without topological constraints

####  3.1.2 Distribution optimization of qubits

In order to describe the optimization of CNOT circuits as a quantifiable
mathematical problem, we first define the distance on any graph
Defination 3: For two vertices $v_{1}$,$v_{2}$ in a graph $G(V,E)$, define
their distance as the length of the least path connecting the two vertices,
written as $d(u_{1},u_{2})$.

Then we have the following lemma.
Lemma 1 : Given a CNOT gate $CNOT_{1,i}$, and denote $d$ to be distance of $1$
and $i$ in a graph $G(V,E)$, then $CNOT_{1,i}$ can be implemented on $G$ with
$s$ CNOT gates and depth $s$ without changing any other qubits, where

$s=\left\\{\begin{array}[]{rcl}4d-4&&{d>1}\\\ 1&&{d=1}\\\ \end{array}\right.$

$Proof.$ It is trivial for the case of $d=1$. If $d>1$, without loss of
generality denote $i=d+1$ and $\\{1,2,...,d+1\\}$ to be $d$ points connecting
$1$ and $i$, and perform the following procedure.
$(a)$ Perform $CNOT_{j,j+1}$, $j$ from $d$ decrease to $1$. ($d$ depth)
$(b)$ Perform $CNOT_{j,j+1}$, $j$ from $2$ increase to $d$. ($d-1$ depth)
Then the input of vertex $1$ is added to all other vertices.
$(c)$ Perform $CNOT_{j,j+1}$, $j$ from $d-1$ decrease to $1$. ($d-1$ depth)
$(d)$ Perform $CNOT_{j,j+1}$, $j$ from $2$ increase to $d-1$. ($d-2$ depth)
Thus, we implement $CNOT_{1,i}$ without changing any other bits. $\Box$

$\Qcircuit<EMAIL_ADDRESS>\qw\qw\qquad\qquad\ctrl{-1}\targ\qw\qw\qw\targ\ctrl{-1}\targ\qw\qw\qw\targ\qw\\\
\qw\qw=\qquad\qw\ctrl{-1}\targ\qw\targ\ctrl{-1}\qw\ctrl{-1}\targ\qw\targ\ctrl{-1}\qw\\\
\qw\qw\qquad\qquad\qw\qw\ctrl{-1}\targ\ctrl{-1}\qw\qw\qw\ctrl{-1}\targ\ctrl{-1}\qw\qw\\\
\ctrl{-4}\qw\qquad\qquad\qw\qw\qw\ctrl{-1}\qw\qw\qw\qw\qw\ctrl{-1}\qw\qw\qw\\\
}$
Figure 6: an example of CNOT circuit between two qubits of distance 4

Lemma 1 guarantees that for any CNOT circuit, as long as our topological
constraint is a connected graph, the quantum circuit can be transformed into
an equivalent quantum circuit that can be implemented under the topological
constraint.

Definition 4: Given an $n\times n$ matrix $A$ over $GF(2)$ and a
representation of $A$ as product of some elementary matrices $E(i+j)$ and a
permutation matrix $A^{\prime}$, denote $S$ as set of these elementary
matrices. Define a multigraph $G(V,E,A)$ corresponds to $A$, where $|V|=n$,
$E=\\{(v_{i},v_{j})|E(i+j)\ or\ E(j+i)\in S\\}$.

Thus we describe the matrix implementation of a CNOT circuit as a graph $G$,
each vertex in the graph corresponds to a qubit, then the problem becomes how
to match $G$ to $G^{\prime}$ such that the cost of implementing all CNOT gates
in $G^{\prime}$ is minimal, where $G^{\prime}$ is the topological constraint
of quantum memories.

Problem 3: Given $G_{1}(V_{1},E_{1})$, $G_{2}(V_{2},E_{2})$,
$|V_{1}|=|V_{2}|$, finding a mapping from $G_{1}$ to $G_{2}$,
$G_{1}\rightarrow G_{2},v_{i}\rightarrow x_{i}$, such that the sum
$\sum_{(v_{i},v_{j})\in E_{1}}s(x_{i},x_{j})$ is the minimum, where $s$ is
defined in Lemma 1.

Unfortunately, there isn’t an algorithm to obtain the theoretical optimal
solution for this problem at present. Here we put forward an optimization
algorithm to find a solution. Our idea is to transform this problem into a
integer programming problem($IP$), then use the solver Gurobi to solve it. We
further optimize the results of Gurobi and finally get a better result.
First of all, we define distance functions for different topological
constraints. Here we cansider several different architectures, such as
16q-square and IBM 20 in Figure 7.

$1$$2$$3$$4$$8$$7$$6$$5$$9$$10$$11$$12$$16$$15$$14$$13$$1$$2$$3$$4$$5$$6$$7$$8$$9$$10$$11$$12$$13$$14$$15$$16$$17$$18$$19$$20$
Figure 7: 16q-square and IBM 20 Figure 8: Distance function of IBM 20

For vertice in 16q-square, we can easliy check that

$d(x_{i},x_{j})=|(x_{i}-1)\ mod\ 4-(x_{j}-1)\ mod\
4|+|\lfloor(x_{i}-1)/4\lfloor-\lfloor(x_{j}-1)/4\lfloor|$

But in some other architectures such as IBM 20, we have to define the function
as a table in Figure 8.

Now we can write our objective function $\sum_{(v_{i},v_{j})\in
E_{1}}s(x_{i},x_{j})$ as a mathematical expression, where $x_{i}\in[1,n]$ and
$x_{i}\neq x_{j}$ if $i\neq j$. Minimizing this function is a typical MILP
problem. For MILP problems, there are relatively mature solvers such as Gurobi
and SAT. Here we use the former. Puting our objective function into Gurobi
with , we can get the optimal solution under the current running time based on
algorithms in Gurobi with Algorithm1. If the running time was infinite, we
would get the optimal solution. But it is impossible for most practical
problems. We can only run Gurobi for a certain amount of time and then get a
current optimal solution.

Algorithm 1 Build and solve MILP models

Input: Point set $X$
Output: Minimum distance

procedure ConstructMILPModel($X$)

2: Precompute the full linear integer inequality characterization of distance
function;

Initialize an empty Model $\mathcal{M}$;

4: $\mathcal{M}.addVar(x_{i,j})$, $i\in[1,n],j\in[0,t]$, where $x_{i,j}$ are
binary variables;

$\mathcal{M}.addVar(s_{i,j}^{k})$, $(i,j)\in E_{1},k\in[0,t]$, where
$s_{i,j}^{k}$ are binary variables;

6: for $i,j\in[1,n]$, $i\neq j$ do

$x_{i}=\sum_{k=0}^{t}2^{k}{x_{i,k}}$;$\triangleright$ Binary to integer
representation

8: $x_{j}=\sum_{k=0}^{t}2^{k}{x_{j,k}}$;

$\mathcal{M}.addConstr(abs(x_{i}-x_{j})\geq 1)$;$\triangleright$ distinct

10: Adding constrains of $x_{i},x_{j}$ and $s_{i,j}^{k}$ to $\mathcal{M}$
according to the full linear integer inequality characterization of distance
function;

end for

12: $\mathcal{M}.setObjective\leftarrow Minimize\sum_{i,j}\sum_{k}2^{k}\cdot
s_{i,j}^{k}$;

$\mathcal{M}.Optimize()$;

14: $D\leftarrow$ the best solution of $\mathcal{M}$; return $D$;

end procedure

As for the results of Gurobi, it is natural to think whether further
optimization can be done. Our idea is permutation of the points. We exchanged
the grid points on the rectangular lattice to calculate the objective
function. If it is reduced, the current result is replaced with the result
after exchange. The number of points exchanged can be changed based on the
different scale of the problem. Specific details are shown in Algorithm 2.

Algorithm 2 Permutation of points

Input: $\\{x_{i}\\},i\in[1,n],x_{i}\in[1,n],x_{i}\neq x_{j}\ if\ i\neq j$,
$Objectivefunction\ \ S$

procedure Permutation($\\{x_{i}\\}$)

2: Calculate $S=S_{initial}=\sum_{(v_{i},v_{j})\in E}s(x_{i},x_{j})$

for each subset of order k $D_{k}=\\{i_{1},i_{2},...,i_{k}\\}\subseteq[1,n]$
do

4: for each $\sigma\in S_{k}$,where $S_{k}$ is a permutation group on $D_{k}$
do

Calculate $S^{\prime}=\sum_{(v_{i},v_{j})\in
E}s(\sigma(x_{i}),\sigma(x_{j}))$,where $\sigma(x_{i})=x_{i}$ if $x_{i}\notin
D_{k}$

6: if $S^{\prime}<S_{initial}$ then

$x_{i}=\sigma(x_{i})\quad S=S^{\prime}$

8: end if

end for

10: end for

while $S<S_{initial}$ do

12: Repeat $Step\ 2\sim 10$

end while

14: return $\\{x_{i}\\}$

end procedure

Here we take a $5\times 5$ matrix in [WHY+19] as an example. The matrix is
same as $A$ in section 3.1.1 and the topological constraint graph $G$ is shown
in Figure 9. By decomposition of $A$, $G(V,E,A)$ is a graph as in Figure 9.
The number of points is so small that we can get the optimal solution only by
algorithm 2. Then our objective function $\sum_{(v_{i},v_{j})\in
G(V,E,A)}s(x_{i},x_{j})$ goes to a minimum of 16 in Figure 10.

$1$$2$$3$$4$$5$$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$ Figure 9: The topological
constraint $G$ and representation of $A$, $G(V,E,A)$
$v_{1}$$v_{2}$$v_{3}$$v_{4}$$v_{5}$
Figure 10: The map between $G(V,E,A)$ and $G$

Then we can get a $CNOT$ circuit of A under topological constraint of $G$ with
$size=16$ as in Figure 11. It is less 4 than the example of method ROWCOL in
[WHY+19].

$\Qcircuit@C=1.5em@R=1em{&a_{1}\qquad\qw\qw\qw\qw\qw\targ\targ\targ\qw\targ\qw\ctrl{4}\targ\ctrl{4}\targ\qw\qw
a_{3}\\\
a_{2}\qquad\targ\ctrl{2}\targ\ctrl{2}\targ\ctrl{-1}\qw\ctrl{-1}\targ\ctrl{-1}\targ\qw\ctrl{-1}\qw\ctrl{-1}\ctrl{1}\qw
a_{4}\\\
a_{3}\qquad\qw\qw\ctrl{-1}\qw\ctrl{-1}\qw\qw\qw\ctrl{-1}\qw\ctrl{-1}\qw\qw\qw\qw\targ\qw
a_{5}\\\
a_{4}\qquad\ctrl{-2}\targ\qw\targ\qw\qw\qw\qw\qw\qw\qw\qw\qw\qw\qw\qw\qw
a_{2}\\\
a_{5}\qquad\qw\qw\qw\qw\qw\qw\ctrl{-4}\qw\qw\qw\qw\targ\qw\targ\qw\qw\qw
a_{1}\\\ }$ Figure 11: the CNOT circuit of $A$ with topological constraints

##  4 Experiments

###  4.1 Comparison with other methods

We compare our algorithm with the most commonly used Steiner Tree algorithm.
For the sake of justice, we use the same matrixes with [KdG19], the number of
each kind of CNOT circuit is 20, and we use a test set of 380 random circuits
in total. The result is listed in Table 1.

Architecture | # | QuilC | t|ket> | Steiner | our method | <Steiner
---|---|---|---|---|---|---
9q-square | 3 | 3.8 | 3.6 | 3 | 2.95 | 1.6%
9q-square | 5 | 10.82 | 6.4 | 5.2 | 4.6 | 11.5%
9q-square | 10 | 20.08 | 16.95 | 11.6 | 11.7 | -0.9%
9q-square | 20 | 46.24 | 40.75 | 25.85 | 23.8 | 7.9%
9q-square | 30 | 72.89 | 66.15 | 35.55 | 31.3 | 12.0%
16q-square | 4 | 6.14 | 5.8 | 4.44 | 4 | 10.0%
16q-square | 8 | 19.68 | 12.95 | 12.41 | 7.65 | 38.4%
16q-square | 16 | 48.13 | 36.2 | 33.08 | 25 | 24.4%
16q-square | 32 | 106.75 | 94.45 | 82.95 | 68.5 | 17.4%
16q-square | 64 | 225.69 | 203.75 | 147.38 | 138.15 | 6.3%
16q-square | 128 | 457.35 | 436.25 | 168.12 | 150.25 | 10.6%
16q-square | 256 | 925.85 | 922.65 | 169.28 | 153.65 | 9.2%
ibm-q20-tokyo | 4 | 5.5 | 6.05 | 4 | 4 | 0
ibm-q20-tokyo | 8 | 17.3 | 12 | 7.69 | 7.8 | -1.4%
ibm-q20-tokyo | 16 | 43.83 | 29.05 | 20.44 | 14.85 | 27.3%
ibm-q20-tokyo | 32 | 93.58 | 78.15 | 66.93 | 49.35 | 26.3%
ibm-q20-tokyo | 64 | 215.9 | 181.25 | 165.6 | 124.2 | 25.0%
ibm-q20-tokyo | 128 | 432.65 | 391.85 | 237.64 | 217.95 | 8.3%
ibm-q20-tokyo | 256 | 860.74 | 789.3 | 245.84 | 219.5 | 10.7%
Table 1: The second colunmn shows the number of CNOT gates of random circuits.
The remaining columns show the average 2-qubit gate count after mapping 20
random circuits.

###  4.2 Application in AES’s Mixcolumns

The quantum circuit implementation of AES is an important problem[BNPS19]
[GLRS16], but most of the researchers neglects the structure of quantum
computer. They acquiesced that the cost of AES’s Mixcolumns is 92 CNOT gates
without constrain which seems a little unreasonable on real quantum computers.
In this section we apply our method in the MC layer of AES which is written as
a $32\times 32$ matrix. We use a $4\times 8$ lattice as topologically-
constrained quantum memories, and get a result of 159 rather than 92 CNOT
gates of the cost of AES’s Mixcolumns. Here is our result.

$v_{3}$$v_{11}$$v_{12}$$v_{28}$$v_{21}$$v_{5}$$v_{30}$$v_{6}$$v_{27}$$v_{10}$$v_{4}$$v_{20}$$v_{13}$$v_{29}$$v_{22}$$v_{14}$$v_{19}$$v_{18}$$v_{1}$$v_{17}$$v_{24}$$v_{32}$$v_{23}$$v_{7}$$v_{2}$$v_{26}$$v_{9}$$v_{25}$$v_{16}$$v_{8}$$v_{15}$$v_{31}$
Figure 12: The routing of 32 qubits in AES’s Mixcolumns, Size=159

## References

  * [BHL+16] CJ Ballance, TP Harty, NM Linke, MA Sepiol, and DM Lucas.  High-fidelity quantum logic gates using trapped-ion hyperfine qubits.  Physical review letters, 117(6):060504, 2016.
  * [BMP13] Joan Boyar, Philip Matthews, and René Peralta.  Logic minimization techniques with applications to cryptology.  Journal of Cryptology, 26(2):280–312, 2013.
  * [BNPS19] Xavier Bonnetain, María Naya-Plasencia, and André Schrottenloher.  Quantum security analysis of aes.  IACR Transactions on Symmetric Cryptology, 2019(2):55–93, 2019\.
  * [BSK+12] Joseph W Britton, Brian C Sawyer, Adam C Keith, C-C Joseph Wang, James K Freericks, Hermann Uys, Michael J Biercuk, and John J Bollinger.  Engineered two-dimensional ising interactions in a trapped-ion quantum simulator with hundreds of spins.  Nature, 484(7395):489–492, 2012.
  * [CDD+19] Alexander Cowtan, Silas Dilkes, Ross Duncan, Alexandre Krajenbrink, Will Simmons, and Seyon Sivarajah.  On the qubit routing problem.  arXiv preprint arXiv:1902.08091, 2019.
  * [GLRS16] Markus Grassl, Brandon Langenberg, Martin Roetteler, and Rainer Steinwandt.  Applying grover’s algorithm to aes: quantum resource estimates.  In Post-Quantum Cryptography, pages 29–43. Springer, 2016.
  * [HND17] Daniel Herr, Franco Nori, and Simon J Devitt.  Optimization of lattice surgery is np-hard.  Npj quantum information, 3(1):1–5, 2017.
  * [HS18] Steven Herbert and Akash Sengupta.  Using reinforcement learning to find efficient qubit routing policies for deployment in near-term quantum computers.  arXiv preprint arXiv:1812.11619, 2018.
  * [KdG19] Aleks Kissinger and Arianne Meijer-van de Griend.  Cnot circuit extraction for topologically-constrained quantum memories.  arXiv preprint arXiv:1904.00633, 2019.
  * [Nay19] C Nay.  Ibm unveils world’s first integrated quantum computing system for commercial use.  Last Accessed: July 4th, 2019.
  * [NC02] Michael A Nielsen and Isaac Chuang.  Quantum computation and quantum information, 2002.
  * [NGM20] Beatrice Nash, Vlad Gheorghiu, and Michele Mosca.  Quantum circuit optimizations for nisq architectures.  Quantum Science and Technology, 5(2):025010, 2020.
  * [Pre18] John Preskill.  Quantum computing in the nisq era and beyond.  Quantum, 2:79, 2018.
  * [ROT+18] Matthew Reagor, Christopher B Osborn, Nikolas Tezak, Alexa Staley, Guenevere Prawiroatmodjo, Michael Scheer, Nasser Alidoust, Eyob A Sete, Nicolas Didier, Marcus P da Silva, et al.  Demonstration of universal parametric entangling gates on a multi-qubit lattice.  Science advances, 4(2):eaao3603, 2018.
  * [SCZ16] Robert S Smith, Michael J Curtis, and William J Zeng.  A practical quantum instruction set architecture.  arXiv preprint arXiv:1608.03355, 2016.
  * [WHY+19] Bujiao Wu, Xiaoyu He, Shuai Yang, Lifu Shou, Guojing Tian, Jialin Zhang, and Xiaoming Sun.  Optimization of CNOT circuits under topological constraints.  arXiv e-prints, page arXiv:1910.14478, October 2019.
  * [XZL+20] Zejun Xiang, Xiangyoung Zeng, Da Lin, Zhenzhen Bao, and Shasha Zhang.  Optimizing implementations of linear layers.  IACR Transactions on Symmetric Cryptology, pages 120–145, 2020\.
  * [ZPW18] Alwin Zulehner, Alexandru Paler, and Robert Wille.  An efficient methodology for mapping quantum circuits to the ibm qx architectures.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 38(7):1226–1236, 2018.


# A fast multi-fidelity method with uncertainty quantification for complex
data correlations: Application to vortex-induced vibrations of marine risers

Xuhui Meng  Zhicheng Wang  Dixia Fan  Michael Triantafyllou  George Em
Karniadakis 111Corresponding author<EMAIL_ADDRESS>Division of
Applied Mathematics, Brown University, Providence, RI 02906, USA  Department
of Mechanical Engineering, Massachusetts Institute Technology, Cambridge, MA
02139, USA  Sea Grant College Program, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA

###### Abstract

We develop a fast multi-fidelity modeling method for very complex correlations
between high- and low-fidelity data by working in modal space to extract the
proper correlation function. We apply this method to infer the amplitude of
motion of a flexible marine riser in cross-flow, subject to vortex-induced
vibrations (VIV). VIV are driven by an absolute instability in the flow, which
imposes a frequency (Strouhal) law that requires a matching with the impedance
of the structure; this matching is easily achieved because of the rapid
parametric variation of the added mass force. As a result, the wavenumber of
the riser spatial response is within narrow bands of uncertainty. Hence, an
error in wavenumber prediction can cause significant phase-related errors in
the shape of the amplitude of response along the riser, rendering correlation
between low- and high-fidelity data very complex. Working in modal space as
outlined herein, dense data from low-fidelity data, provided by the semi-
empirical computer code VIVA, can correlate in modal space with few high-
fidelity data, obtained from experiments or fully-resolved CFD simulations, to
correct both phase and amplitude and provide predictions that agree very well
overall with the correct shape of the amplitude response. We also quantify the
uncertainty in the prediction using Bayesian modeling and exploit this
uncertainty to formulate an active learning strategy for the best possible
location of the sensors providing the high fidelity measurements.

###### keywords:

elsarticle.cls, LaTeX, Elsevier , template

###### MSC:

[2010] 00-01, 99-00 modal decomposition , multi-fidelity , uncertainty
quantification , active learning , vortex-induced vibrations

††journal: Journal of LaTeX Templates

##  1 Introduction

Vortex-Induced Vibrations. Long slender structures within steady oncoming flow
are subject to vibrations caused by vortical structures forming due to a
distributed flow instability in their wake [1], as shown in Fig. 1, and the
flexible cylinder deflects due to an increase in drag caused by vortex
shedding. The problem has considerable theoretical interest as it constitutes
a fundamental nonlinear flow-structure interaction (FSI) system [2], while it
is very important for the design of offshore industry systems, such as risers,
cables, and hawsers, which are subject to large drag loads and potentially
catastrophic fatigue damage as a consequence of the vortex-induced vibration
(VIV) [3]. Therefore, a large amount of research has focused on better
understanding and predicting the flexible cylinder VIV response, as well as
devising suppression methods to mitigate excessive fatigue damage.

Figure 1: Large eddy simulation of the vortex induced vibration of a flexible
riser in linearly sheared flow. Top figure, instantaneous vortical flow
patterns; bottom figure, the corresponding riser displacement (red curve).
Note that, the magnitude of the displacement has been enlarged in order to
visualize it more clearly.

Similar to rigid cylinder VIV [4, 5, 6], flexible cylinder VIV display FSI
features, including “lock-in” [7], dual resonance of in-line (IL) and cross-
flow (CF) vibration [8, 9], fluid-structure energy exchange that is high
correlated with the form of the IL-and-CF motion trajectory [10, 11, 12], and
large variation of the effective fluid added mass at different flow speeds
[13, 14]. However, flexible cylinder VIV have additional features caused by
the distributed flow instability along their length, which can cause multiple
frequencies of response and variable spatial correlation coherence of the
vortical patterns. A few well-controlled laboratory experiments and large-
scale field test have helped to shed some light on these features. Fan et al.
[9] measured directly the displacement of a flexible cylinder in uniform flow
over a wide range of reduced velocities at moderate Reynolds number using an
underwater optical measurement system. Their work showed that various modes
are excited in the IL and the CF directions, documented jumps in maximum
amplitude as function of oncoming flow velocity, and mapped the distributed
hydrodynamic coefficient changes during a modal switch. Similar results were
reported by Chaplin et al. [15] from a larger-scale experiment of a flexible
cylinder in a stepped incoming flow. The results were compared with
predictions from various numerical and semi-empirical simulations [16]. The
predictions, although capturing some flexible cylinder VIV features, also
showed discrepancies from the experiments. Under the Norwegian Deepwater
Program (NDP), experiments were conducted on a 38-m flexible riser model in
both uniform and linearly sheared flow [17]. It was shown that the riser model
would exhibit chaotic response, switching randomly between steady-state and
chaotic states, with strong high harmonic components in the lift force during
steady-state response [18, 19, 20]. Field tests, such as towing long flexible
cylinders within the Gulf Stream [21, 22, 23], displayed also the presence of
mixed traveling and standing waves, a strong interaction between IL and CF
modal responses, and a response that consisted of multiple frequencies.

Various computational fluid dynamics (CFD) tools have been developed to
predict the structural and wake response of flexible cylinders in cross-flow
[24, 25, 26, 27], but due to the large computation resources required, CFD
research has focused on relatively low Reynolds numbers ($Re\sim
O(10)-O(10^{3})$) and small to medium aspect ratios ($L/d\sim
O(10)-O(10^{3})$) , which are different from full-scale conditions in riser
deployed in offshore platforms ($Re\sim O(10^{5})$ and $L/d\sim O(10^{4})$).
Therefore, the state-of-the-art riser VIV prediction tools widely adopted in
the industry are often based on semi-empirical model predictions [28, 29, 30].
These codes have been developed based on the assumption of strip theory [9,
31] and employing various spatial correlation models. The hydrodynamic model
is hence simplified, adapting the hydrodynamic coefficients of a rigid
cylinder conducting forced vibrations within a cross-stream [32, 33]. A key
factor for riser VIV prediction via the semi-empirical approach is the
availability of accurately measured hydrodynamic coefficients as distributed
along the riser span. However, various studies showed that these coefficients
are sensitive to several parameters, including Reynolds number [34], riser
configuration [35], and surface roughness [36]. Hence, a systematic
development of a single hydrodynamic database is virtually impossible. In
addition, during the lifetime of a riser in the field, long-term effects, such
as equipment aging and biofouling, inevitably alter the hydrodynamic
coefficients, making long-term riser prediction and monitoring even more
challenging.

The objective of this study is to develop a fast vibration inference system by
integrating low fidelity, low cost, semi-empirical prediction tools with high
fidelity, high cost numerical and experimental approaches. Here, the high
fidelity data consists of three sources: First, we use LES simulation results
of VIV of a flexible cylinder in linearly sheared inflow. In the simulation,
the reduced velocity is $U_{r}=\frac{U_{m}}{f_{n1}D}=15.65$, where
$U_{m}=(U_{max}+U_{min})/2$ is the mean inflow velocity, $f_{n1}$ is the first
modal natural frequency, $D$ is diameter of the cylinder. Here $U_{max}=1.4$
and $U_{min}=0.6$ are the inflow velocity at two ends of the current profile,
respectively. The Reynolds number $Re=U_{m}\,D/\nu=800$, where $\nu$ is the
kinematic viscosity. Moreover, in the simulation, the aspect ratio of the
flexible cylinder $L/D=240$. Second, we use optical measurement experiments
that were conducted at the MIT Tow Tank facility, using a flexible cylinder
with $L/D=240$, towed to generate a uniform inflow at various speeds, to
achieve $Re$ from 250 to 2,300, while the $U_{r}$ varies from 4.8 to 36. Note
that, most of the high fidelity data presented in this paper are from the
optical experiment. Third, we use the measurement data from NDP case 2430
[17], to demonstrate our method is capable of predicting the VIV motions
accurately in realistic conditions. In the NDP experiment a bare riser of
$L=38$ m and $D=0.027$ m is towed horizontally at $U_{min}=0$, $U_{max}=1.5$
m/s. Having the high fidelity data, we subsequently performed calculations
using the VIVA program [28] to obtain the low fidelity data, with the same
flow and structural parameters.

Multi-fidelity Modeling. To predict the VIV response of a riser, we combine a
dense set of data from a low fidelity model because of its low cost of use,
with a relatively sparse set of experimental and field measurements that are
expensive to obtain. With abundant simulation data from VIVA and a small set
of accurate measurements, a natural idea is to utilize multi-fidelity
modeling, such as multi-fidelity Gaussian process regression (GPR) [37, 38,
39, 40, 41], or multi-fidelity deep neural networks [42, 43], to assimilate
data with different accuracy. In particular, we will use here the VIVA data as
low-fidelity and experimental measurements as high-fidelity data.

The key idea of multi-fidelity modeling is to discover the cross-correlation
between the low- and high-fidelity data rather than approximating the high-
fidelity data directly. Conventional multi-fidelity modeling approaches work
well for cases in which the correlation is easier to obtain than the high-
fidelity function itself. However, the correlation between the low- and the
high-fidelity data in the present study is much more complicated than the
high-fidelity function, due to phase errors between the multi-fidelity data.
The aforementioned multi-fidelity GPR [38] has been extended to resolve the
complicated correlation caused by the phase error between the low- and high-
fidelity data by adding more dimensions, i.e., shifts of low-fidelity
functions, to the input space. To obtain accurate predictions, additional
dimensions are used for various low- and high-fidelity data in [44], but how
to choose a proper number of additional dimensions remains unclear.
Furthermore, the increasing input dimensions will clearly lead to higher
computational cost. Therefore, we need to develop alternative multi-fidelity
approaches, which can overcome these drawbacks. To this end, we develop a
multi-fidelity approach, which we refer to as VIV-MFnet, to assimilate the VIV
simulation and the experimental data that can handle phase errors effectively
and efficiently. The rest of the paper is organized as follows: we present the
VIV-MFnet method in Sec. 2, and we show our results in Sec. 3. Finally, we
give a summary for this work in Sec. 4.

##  2 Methodology

###  2.1 Fourier decomposition of displacement

Figure 2:  Correlation between low- and high-fidelity data for two
representative cases. (a) Uniform flow past a marine riser at $Ur=12.66$. Left
panel: Low- and high-fidelity data. LF: low-fidelity data that are from VIVA
model, HF: high-fidelity data that are from MIT experiment. Right panel:
Correlation between the low- and high-fidelity data. (b) Uniform flow past a
marine riser at $Ur=14.54$. Left panel: Low- and high-fidelity data. LF: low-
fidelity data from the VIVA model, HF: high-fidelity data from the MIT
experiment. Right panel: Correlation between the low- and high-fidelity data.
$A_{y}$: displacement, $z$: location, $A_{y,l}$: low-fidelity data
(displacements from VIVA), $A_{y,h}$: high-fidelity data (measurements from
experiments or field observations). Note that all the displacements
($A_{y}$/$A_{y,l}/A_{y,h}$) and locations ($z$) in this study are normalized
by the diameter ($d$) and the length of the cylinder ($L$), respectively.

Correlations between the low- and high-fidelity models for two representative
cases are illustrated as examples in Fig. 2. Existing multi-fidelity methods
perform poorly with such complex correlation functions. We observe that there
are differences between the low- and high-fidelity data both in the amplitude
as well as the spatial phase. Hence, for accurate inference we need to be able
to vary the amplitude and phase of the low fidelity data. Therefore, we
transform the displacements in the physical space into the modal space using
the fast Fourier transformation (FFT) as

$\displaystyle A_{m,k}=\sum^{N-1}_{n=0}A_{y,m,n}\exp(-i2\pi
kn/N),\leavevmode\nobreak\ n=0,...,N-1,$ (1)

where $m=l$ or $h$ represents the data at the low- and high-fidelity level,
respectively, $A_{m,k}=(\alpha(0),...,\alpha(k),...)$ denotes the displacement
in the modal space, $k$ is the wavenumber, $\alpha(k)$ is a complex number
denoting the coefficient for the $k_{th}$ wave, and $N$ is the total number of
the modes, which is 1,024 in the present work. As displayed in Fig. 3, the
low- and high-fidelity data are quite similar in modal space. For example, the
module in Fig. 3 first decreases and then increases with the increasing $k$ in
subdomain I, and then it decreases in subdomain II with increasing $k$.
Similar results can also be observed in Fig. 3. Hence, we conduct multi-
fidelity modeling in the modal space.

We now use the case shown in Fig. 3 as an example to explain the multi-
fidelity model employed in the present study. As observed in the low-fidelity
data of Fig. 3, the modules of the wavenumbers for $k>6$ are smaller than
$0.01$. Therefore, we can keep the first 7 modes and set the remaining modules
as equal to zero. Note that the last 6 modes are kept the same as
$\alpha(1)-\alpha(6)$ due to the symmetry property of FFT. We then employ the
inverse FFT (iFFT) to transform them into the physical space to obtain the
low-fidelity function with good accuracy. Due to the similarity of the
functional shape of the low- and high-fidelity data in the modal space, we
assume that we can accurately predict the high-fidelity function using the
first 7 modes. Predicting the high-fidelity function is now simplified to
determining the coefficients for the retained modes in the modal space.

Algorithm 1 VIV-Multi-fidelity network (VIV-MFnet)

Low- (VIVA simulations) and high-fidelity (Experiment/Field observations)
data.

$\bullet$ Perform FFT for low-fidelity data to obtain $A_{l,k}$.

$\bullet$ Obtain $k_{c}$ and $k_{t}$ based on the low-fidelity data in modal
space.

$\bullet$ Set $\alpha(k)$ as zero for $k>k_{t}$.

$\bullet$ Optimize $\mathcal{L}(\bm{\alpha})$ in Eq. (3) using the Adam
optimizer to obtain $\alpha(k)$ for $k\leq k_{t}$.

$\bullet$ Conduct HMC to get samples for $\alpha(k)$ from $P(\bm{\alpha}|D)$
in Eq. (4).

$\bullet$ Get predictions $\hat{A}_{y,h}$ from Step 4 or 5.

Figure 3:  Low- and high-fidelity data in the frequency space for two
representative cases. (a) Uniform flow past a marine riser at $Ur=12.66$.
Upper panel: Low-fidelity data in the frequency space. Lower panel: Low-
fidelity data in the frequency space. (b) Uniform flow past a marine riser at
$Ur=14.54$. Upper panel: Low-fidelity data in the frequency space. Lower
panel: High-fidelity data in the frequency space. $k$: wavenumber,
$|A_{l,k}|$: module for the $k_{th}$ wave of the low-fidelity data, and
$|A_{h,k}|$: module for the $k_{th}$ wave of the high-fidelity data. The red
dashed line divides the frequency domain into two sub-domain I and II. The
maximal amplitude for waves in I is $O(0.1)$, while it is $O(0.01)$ in II.

###  2.2 Maximum a posteriori probability (MAP) estimate

To obtain the coefficients in the modal space, we can minimize the following
loss function:

$\displaystyle\mathcal{L}(\bm{\alpha})=|A_{y,h}^{*}-\hat{A}_{y,h}|^{2}+\lambda|\bm{\alpha}|^{2},$
(2)

where $A_{y,h}^{*}$ represents the high-fidelity measurements, $\hat{A}_{y,h}$
is the high-fidelity prediction obtained from the iFFT, and the second term is
the $\mathcal{L}_{2}$ regularization used to prevent overfitting. In addition,
$\bm{\alpha}$ denotes the coefficients in the modal space, and $\lambda$ is
the weight for $\bm{\alpha}$.

The entire frequency space is divided into two subdomains, I and II, as shown
in Fig. 3. We observe that the maximum modules for the low-fidelity data in
subdomains I and II are of order $O(0.1)$ and $O(0.01)$, respectively. We
assume that this also holds for the high-fidelity data, and then utilize
different weights for the coefficients in these two subdomains. Specifically,
we can rewrite Eq. (2) as follows:

$\displaystyle\mathcal{L}(\bm{\alpha})=|A_{y,h}^{*}-\hat{A}_{y,h}|^{2}+\lambda_{I}|\bm{\alpha}_{I}|^{2}+\lambda_{II}|\bm{\alpha}_{II}|^{2},$
(3)

where $\lambda_{I}(\bm{\alpha}_{I})$ and $\lambda_{II}(\bm{\alpha}_{II})$ are
the weights in subdomain I and II, respectively. In addition, the first term
in Eq. (3) is generally as small as O($10^{-4}$) based on our experience.
Therefore, we can set $\lambda_{I}\thicksim O(0.01)$ and
$\lambda_{II}\thicksim O(0.1)$ or $\lambda_{II}\thicksim O(1)$ to guarantee
that the three terms in Eq. (3) are of the similar order, i.e., O($10^{-4}$).
Note that no regulation is employed for $k=0$.

The choice of $\lambda_{I}$ and $\lambda_{II}$ is crucial to ensure accurate
inference. Generally, we can set $\lambda_{I}=0.01$ and
$\lambda_{II}=\lambda_{I}|\alpha_{I,max}|/|\alpha_{II,max}|$, where
$\alpha_{I,max}$ and $\alpha_{II,max}$ are the maxima in subdomain I and II,
respectively. In particular, we can tune $\lambda_{II}$ based on $\lambda_{I}$
and the ratio between the maximum modules in subdomain I and II. In addition,
the critical wavenumber between the two subdomains is denoted as $k_{c}$, and
the truncation wavenumber is denoted as $k_{t}$. We note that the values of
$\alpha_{I,max}$, $\alpha_{II,max}$, $k_{c}$, and $k_{t}$ are determined from
the low-fidelity data. Specifically, we employ $\lambda_{I}=0.01$ and
$\lambda_{II}=0.25$, which works well for all cases studied in this work.

###  2.3 Bayesian inference (BI) for uncertainty quantification

In the Bayesian framework, each coefficient is a distribution rather than the
point estimate in Sec. 2.2. Specifically, by denoting the observed data as
$\mathcal{D}$, the posterior distribution of $\bm{\alpha}$ can be expressed as
follows, based on Bayes’ rule:

$\displaystyle
P(\bm{\alpha}|\mathcal{D})=\frac{P(\bm{\alpha})P(\mathcal{D}|\bm{\alpha})}{P(\mathcal{D})},$
(4)

where $P({\bm{\alpha}})$ is the prior distribution of the unknown
coefficients, and $P(\mathcal{D}|\bm{\alpha})$ is the likelihood.
$P(\bm{\alpha}|\mathcal{D})$ is basically intractable since $P(\mathcal{D})$
is unknown. However, we can use the Hamiltonian Monte Carlo (HMC) [45, 46] to
sample from $P(\bm{\alpha}|\mathcal{D})$. Similar as the regularization in
Sec. 2.2, we divide the prior $P({\bm{\alpha}})$ into two parts and set
different priors for each subdomain, to prevent overfitting. In particular, we
employ $P({\bm{\alpha}_{I}})\thicksim\mathcal{N}(0,0.1^{2})$ and
$P({\bm{\alpha}_{II}})\thicksim\mathcal{N}(0,0.01^{2})$ in this study. We
would like to point out that MAP is computationally very efficient but does
not predict uncertainty, while BI is relatively more expensive but is capable
of handling noisy data as well as quantifying uncertainty in predictions.

##  3 Results and Discussion

In this section, we first apply the approach described in Sec. 2.2, i.e., MAP,
to predict riser displacements in uniform and sheared flows. In particular, we
conduct a comparison between the present method and the state-of-art multi-
fidelity approach, MF-DNNs [42], which is capable of capturing nonlinear
cross-correlation between the low-and high-fidelity data. We then employ the
method described in Sec. 2.3, i.e., BI, to quantify uncertainty in our
predictions. In MAP, the Adam optimizer with an initial learning rate is set
as $10^{-3}$, and the number of training steps is 10,000. In HMC, we set the
number of burn-in steps as 1,000 with the leapfrog step 20; also, the initial
time step is set as $0.1$. Finally, 1,000 samples are employed to compute the
mean and standard deviation after the burn-in steps.

###  3.1 Results from MAP

####  3.1.1 Lab experiments for uniform and shear flows

We fist test the case with uniform flow at $Ur=12.66$. The critical and
truncated wavenumbers for this case are $k_{c}=2$ and $k_{t}=6$, respectively,
according to Fig. 3. In our multi-fidelity modeling, we assume that we have 4
high-fidelity measurements, which are randomly distributed in $z\in[0,1]$ as
displayed in Fig. 4. The predicted high-fidelity profile with
$\lambda_{I}=0.01$ and $\lambda_{II}=0.25$ is illustrated in Fig. 4, which
agrees well with the experimental results. In particular, predicting the
maximum displacement is significant in the design of riser, because it
determines the stress level: the maximum displacements for the multi-fidelity
modeling in Fig. 4 and the experiment are 0.6968 and 0.6769, respectively. In
addition, we use the recently developed multi-fidelity deep neural networks
because they are highly expressive in capturing nonlinear correlation between
the low- and high-fidelity data [42]. Here we also present the results using
the MF-DNNs in Fig. 4 for comparison. As shown, the present method outperforms
the MF-DNNs, especially for predictions in the interval $z\in[0.5,0.7]$.

Figure 4:  Importance of regularization: Displacements for uniform flow past a
marine riser at $Ur=12.66$. (a) Multi-fidelity predictions with
$\mathcal{L}_{2}$ regularization $\lambda_{I}=0.01$, and $\lambda_{II}=0.25$.
(b) Multi-fidelity predictions with no $\mathcal{L}_{2}$ regularization, i.e.,
$\lambda_{I}=\lambda_{II}=0$. (c) Multi-fidelity predictions with
$\mathcal{L}_{2}$ regularization $\lambda_{I}=\lambda_{II}=0.01$. Blue dot:
low-fidelity training data, Blue circle: high-fidelity training data. HF: MIT
experiment, VIV-MFnet: multi-fidelity modeling in Sec. 2.2, MF-DNN: multi-
fidelity modeling using [42]. The first DNN in MF-DNNs has 2 hidden layers
with 40 neurons per layer, and the second DNN in MF-DNNs has 2 hidden layers
with 20 neurons per layer. The hyperbolic tangent function is employed as the
activation function in both DNNs. The Adam optimizer is used for optimization.
More details for MF-DNNs can be found in [42].

Next, we study the effect of $\mathcal{L}_{2}$ regularization on the accuracy
of the inference. Specifically, two additional scenarios are considered, viz.,
no $\mathcal{L}_{2}$ regularization is employed in either subdomain, and the
same weight ($\bm{\lambda}>0$) is used in both subdomains. We observe that:
(1) it is easy to get overfitting if no $\mathcal{L}_{2}$ regularization is
used (Fig. 4), and (2) the usage of different weights for waves in different
subdomains can enhance the predicted accuracy when comparing the results in
Fig. 4 to those in Fig. 4.

We now consider the case in Fig. 2, in which the critical and truncated
wavenumbers are set as $k_{c}=3$ and $k_{t}=6$, respectively, based on Fig. 3.
Similarly, we perform two tests with different regularizations, i.e., (1)
$\lambda_{I}=0.01$, $\lambda_{II}=0.25$, and (2) $\lambda_{I}=\lambda_{II}=0$.
We assume that we have 6 random samples for training data here. As shown in
Fig. 5, the predictions from multi-fidelity modeling with $\mathcal{L}_{2}$
regularization are in good agreement with the experimental results, while
overfitting is observed in the results without $\mathcal{L}_{2}$
regularization.

Figure 5:  Multi-fidelity predictions of displacements for uniform flow past a
marine riser at $Ur=14.54$. (a) With $\mathcal{L}_{2}$ regularization
$\lambda_{I}=0.01$, and $\lambda_{II}=0.25$. (b) With no $\mathcal{L}_{2}$
regularization, i.e., $\lambda_{I}=\lambda_{II}=0$. HF: MIT experiment. Blue
dot: low-fidelity training data, Blue circle: high-fidelity training data.

Next, we test the performance of the proposed method for predicting
displacements of the riser within uniform flow at different velocities, i.e.,
$Ur=21.64$ and $26.26$; 8 and 10 random samples are employed as the training
data for these two cases, respectively. As shown in Figs. 6-7, the multi-
fidelity predictions agree well with the experimental results.

Figure 6:  Multi-fidelity predictions of displacements for uniform flow past a
marine riser at $Ur=21.64$. (a) Left panel: Low- and high-fidelity data. LF:
low-fidelity data from the VIVA model, HF: high-fidelity data from the MIT
experiment. Right panel: Correlation between the low- and high-fidelity data.
(b) Left panel: Low- and high-fidelity data in the frequency space. First row:
low-fidelity data, Second row: high-fidelity data. The critical and truncated
wavenumbers are $k_{c}=4$ and $k_{t}=9$, respectively. Right panel: With
$\lambda_{I}=0.01$, and $\lambda_{II}=0.25$. Blue dot: low-fidelity training
data, Blue circle: high-fidelity training data.

Figure 7:  Multi-fidelity predictions of displacements for uniform flow past a
marine riser at $Ur=26.26$. (a) Left panel: Low- and high-fidelity data. LF:
low-fidelity data from the VIVA model, HF: high-fidelity data from the MIT
experiment. Right panel: Correlation between the low- and high-fidelity data.
(b) Left panel: Low- and high-fidelity data in the frequency space. First row:
low-fidelity data, Second row: high-fidelity data. The critical and truncated
wavenumbers are $k_{c}=5$ and $k_{t}=12$, respectively. Right panel: With
$\lambda_{I}=0.01$, and $\lambda_{II}=0.25$. Blue dot: low-fidelity training
data, Blue circle: high-fidelity training data.

Note that the number of waves in all the cases tested above is the same for
the data from VIVA simulations and experiments. However, there are some cases
in which VIVA predicts a different number of waves from the experiments.
Specifically, two different cases, i.e., $Ur=17.23$, and 35.28, are
considered. In both cases, the number of half-wavelengths for the low-fidelity
data is one more than that in the high-fidelity data. As for the first test
case in Fig. 8, the critical and truncated wavenumbers here are $k_{c}=4$ and
$k_{t}=10$, respectively. We first employ 6 random samples as the training
data. As illustrated in the last panel of Fig. 8, the phase error between the
low- and high-fidelity data can be corrected. In addition, the predicted
maximum displacement is close to that from the experiment. To further improve
the predicted accuracy, we then randomly add one more training data (blue
cross in the right panel in Fig. 8). We see that the error between the maximum
displacements from the present method and experiment is further reduced. In
the second case (Fig. 9), the number of half-wavelengths in the low-fidelity
data is 7, while we have 6 waves in the high-fidelity data. The critical and
truncated wavenumbers are set as 7 and 22, respectively. The multi-fidelity
predictions for both the phase as well as the amplitude are in good agreement
with the experiments as 13 random training samples are used.

Figure 8:  Multi-fidelity predictions of displacements for uniform flow past a
marine riser at $Ur=17.23$. (a) Left panel: Low- and high-fidelity data. LF:
low-fidelity data that are from VIVA model, HF: high-fidelity data that are
from MIT experiment. Middle panel: Correlation between the low- and high-
fidelity data. Right panel: Low- and high-fidelity data in the frequency
space. First row: low-fidelity data, Second row: high-fidelity data. Note that
the critical wavenumber is different in the low- (i.e., $k_{c}=4$) and high-
fidelity (i.e., $k_{c}=3$) data. (b) Left panel: With 6 random training data.
Right panel: With 7 random training data. In both cases, $\lambda_{I}=0.01$,
and $\lambda_{II}=0.25$. Blue dot: low-fidelity training data, Blue cross:
high-fidelity training data, Blue cross: added high-fidelity training data.

Figure 9:  Multi-fidelity predictions of displacements for uniform flow past a
marine riser at $Ur=35.28$. (a) Left panel: Low- and high-fidelity data. LF:
low-fidelity data from the VIVA model, HF: high-fidelity data from the MIT
experiments. Right panel: Correlation between the low- and high-fidelity data.
(b) Left panel: Low- and high-fidelity data in the frequency space. First row:
low-fidelity data, Second row: high-fidelity data. The critical and truncated
wavenumbers are $k_{c}=7$ and $k_{t}=22$, respectively. Right panel: With
$\lambda_{I}=0.01$, and $\lambda_{II}=0.25$. Blue dot: low-fidelity training
data, Blue circle: high-fidelity training data.

Next, We predict displacements of a riser caused by linear shear flow.
Specifically, the maximum velocity of the shear flow is $Ur=15.65$. The low-
and high-fidelity data are shown in Fig. 10. Here the high-fidelity data are
from a large eddy simulation (LES) rather than from experiments. Similarly, we
assume that we have 6 random samples as the training data. The predictions
from the multi-fidelity modeling agree well with results from the high-
resolution LES.

Figure 10:  Multi-fidelity predictions of displacements for linearly sheared
flow past a marine riser at $Ur=15.65$. (a) Left panel: Low- and high-fidelity
data. LF: low-fidelity data from the VIVA model, HF: high-fidelity data from
LES simulations. Right panel: Correlation between the low- and high-fidelity
data. (b) Left panel: Low- and high-fidelity data in the frequency space.
First row: low-fidelity data, Second row: high-fidelity data. The critical and
truncated wavenumbers are $k_{c}=3$ and $k_{t}=9$, respectively. Right panel:
With $\lambda_{I}=0.01$, and $\lambda_{II}=0.25$. Blue dot: low-fidelity
training data, Blue circle: high-fidelity training data.

###  3.2 Results from UQ and active learning

In this section, we use the Bayesian inference in Sec. 2.3 to quantify the
predicted uncertainty. Moreover, we also conduct active learning to find the
best location for adding more training samples based on the predicted
uncertainty as in [47]. Note that the uncertainty is represented by the
standard deviation.

We first test the performance of the BI using the case in Fig. 2. Note that
both $k_{c}$ and $k_{t}$ for the present case are kept the same as in Fig. 4,
which are obtained from the low-fidelity data. In addition, we randomly
sampled 4 high-fidelity data as the training data. As mentioned in Sec. 2.3,
$P({\bm{\alpha}_{I}})\thicksim\mathcal{N}(0,0.1^{2})$, and
$P({\bm{\alpha}_{II}})\thicksim\mathcal{N}(0,0.01^{2})$ are employed as priors
in the two subdomains. We observe in Fig. 11 that: (1) the predicted means are
in good agreement with the experimental data, and (2) the errors between the
predicted means and the experiments are bounded by the two standard deviations
in most regions. Furthermore, we also present the results with different prior
distributions, i.e., $P({\bm{\alpha}_{I}})\thicksim\mathcal{N}(0,0.1^{2})$,
and $P({\bm{\alpha}_{II}})\thicksim\mathcal{N}(0,0.1^{2})$. We see that: (1)
the predicted means are not as accurate as in Fig. 11, and (2) fluctuations
are observed in the predicted uncertainties, especially for $z\in[0,0.4]$,
which may be attributed to the overfitting caused by the inappropriate prior
distributions. All these results demonstrate that the use of different prior
distributions for different wavenumbers can enhance the predicted accuracy,
which is consistent with the results in Sec. 3.1.

Figure 11:  Multi-fidelity Bayesian modeling of displacements for uniform flow
past a marine riser at $Ur=12.66$: Importance of prior distributions. (a)
$P(\bm{\alpha}_{I})\thicksim\mathcal{N}(0,0.1^{2})$, and
$P(\bm{\alpha}_{II})\thicksim\mathcal{N}(0,0.01^{2})$. (b)
$P(\bm{\alpha}_{I})\thicksim\mathcal{N}(0,0.1^{2})$, and
$P(\bm{\alpha}_{II})\thicksim\mathcal{N}(0,0.1^{2})$. HF: MIT experiment,
Mean: predicted means from multi-fidelity predictions, 2 Stds: two standard
deviations. Blue dot: low-fidelity training data, Blue circle: high-fidelity
training data.

We now implement an active learning strategy based on the case in Fig. 9.
Similarly, we employ $P({\bm{\alpha}_{I}})\thicksim\mathcal{N}(0,0.1^{2})$,
and $P({\bm{\alpha}_{II}})\thicksim\mathcal{N}(0,0.01^{2})$ as the prior
distributions in the two subdomains. We start with 7 training data, which are
randomly distributed in $z\in[0,1]$. As shown in Fig. 12, the phase errors are
well corrected while the predicted means do not fit the experiments quite
well, which is also reflected in the large uncertainty. In addition, the
computational errors between the predicted means and the experimental results
are bounded by the two standard deviations. To reduce the uncertainty as well
as enhance the predicted accuracy, we add more training data using the active
learning strategy. Specifically, we can add one more training data at the
location where the standard deviation is maximum, and then perform Bayesian
inference using the new training dataset. We see in Figs. 12-12 that: (1) the
predicted uncertainty decreases as we add more training data, (2) the
predicted means are in good agreement with the experimental results as we
employ 15 training data as shown in Fig. 12, and (3) the computational errors
are bounded by the two standard deviations for all results in Fig. 12.

Figure 12:  Multi-fidelity Bayesian modeling of displacements for uniform flow
past a marine riser at $Ur=35.28$. (a) With 7 training data. (b) With 10
training data. (c) With 15 training data. Note that the low-fidelity training
data are the same as in Fig. 9, which are not presented here. HF: MIT
experiment, Mean: predicted means from multi-fidelity predictions, 2 Stds: two
standard deviations. Blue plus: added high-fidelity training data based on the
predicted uncertainty.

Finally, we consider a shear flow case in which the low- and high-fidelity
data are from VIVA simulations and the NDP experimental data, respectively
(left panel in Fig. 13). Based on the low-fidelity data in the modal space
(right panel in Fig. 13), modules for all the waves are below 0.1, we then use
the same prior distribution for $\bm{\alpha}$ in the entire domain, i.e.,
$P(\bm{\alpha})=\mathcal{N}(0,0.02^{2})$ (Note that 0.02 is about half of the
maximum $A_{l,k}$). We start with 16 random training samples. As shown in Fig.
13, the predicted means are not in good agreement with the NP data, which is
also reflected in the large estimated uncertainty. Similarly, we can improve
the predicted accuracy by performing active learning to add more training
samples. In particular, the computational errors are bounded by the two
standard deviations in most of the areas of $z\in[0,1]$.

Figure 13:  Multi-fidelity Bayesian modeling of displacements for linearly
sheared flow past a marine riser (NDP case 2430). (a) Left panel: Low- and
high-fidelity data. LF: low-fidelity data that are from VIVA model, HF: high-
fidelity data that are from the NDP data. Middle panel: Correlation between
the low- and high-fidelity data. Right panel: Low- and high-fidelity data in
the frequency space. First row: low-fidelity data, Second row: high-fidelity
data. The truncation wavenumber is $k_{t}=20$. (b) From left to right: 16, 24,
and 28 training data. HF: NDP data, Mean: predicted means from multi-fidelity
predictions, 2 Stds: two standard deviations. Blue dot: low-fidelity training
data, Blue circle: initial high-fidelity training data, Blue plus: added high-
fidelity training data based on the predicted uncertainty.

##  4 Summary

We present a multi-fidelity framework to predict vortex-induced vibrations of
marine risers in modal space. This method, which we call VIV-MFnet, is capable
of assimilating dense inexpensive low-fidelity data, obtained from the semi-
empirical code VIVA, with a small set of high-fidelity data, consisting of
displacement measurements obtained from either fully-resolved 3D LES or
experiments on large-aspect-ratio flexible risers. Our aim is to enhance the
predicted accuracy compared to single-fidelity modeling, which consists of
performing regression with high-fidelity data only. In particular, the low-
fidelity data are used to extract useful information on selecting basis
functions as well as regularization in modal space to prevent overfitting,
while the high-fidelity data are utilized to correct the phase/amplitude
errors between the low- and the high-fidelity data.

We applied two optimization approaches for the VIV-MFnet, the Maximum A
Posteriori probability (MAP) estimate, and the Bayesian inference (BI). The
selection of the two methods provides a balance between computational
efficiency and advanced learning capability. Based on gradient descent, the
first approach of MAP is more computationally efficient, about five times
faster than BI, while the BI provides quantified uncertainty prediction,
enabling active learning to adaptively improve the accuracy of prediction.

We apply the VIV-MFnet on multiple datasets from various sources: small-scale
lab experiments, LES simulations, and large-scale experiments from the
Norwegian Deepwater Program (NP), to test and demonstrate the effectiveness
and robustness of the present method. Our main findings are as follows:

  1. 1.

Using MAP, the vortex-induced vibration of a flexible cylinder in both uniform
and shear flows is studied, and vibration displacements along the entire span
are accurately reconstructed given only a small set of high-fidelity data and
dense samples from low-fidelity data.

  2. 2.

Comparison with the prediction results from state-of-the-art multi-fidelity
deep neural networks shows that our approach outperforms these methods, which
fail to capture the complex correlation between the low- and high-fidelity
data that results in spatial phase errors.

  3. 3.

Using Bayesian inference, we are able to quantify the uncertainty in our
predictions along the riser span, given limited high-fidelity observations.

  4. 4.

Guided by adaptive sampling strategy to minimize the global uncertainty, we
illustrate that active learning can serve to adaptively select the optimal
high-fidelity measurement locations along the span in order to improve the
predicted accuracy.

Finally, it is noteworthy that the current approach using VIV-MFnet is
computationally light, as it takes only about four minutes to learn the NDP
experimental data of riser response on a personal laptop with 2-core CPU
(Intel i7-8650U with 1.9 GHz). Such efficiency is important for future
implementation of a digital twin of marine risers that is capable of providing
real-time prediction and monitoring of production riser VIVs over their
lefetime.

## Acknowledgement

We gratefully acknowledge support from the DigiMaR Consortium, consisting of
ABS, ExxonMobil, Petrobras, SAIPEM, Shell, and Subsea 7. G.E.K., Z.W., and
X.M. would like to acknowledge the support by PhILMS grant DE-SC0019453, and
OSD/AFOSR MURI grant FA9550-20-1-0358.

## References

  * [1] X. Wu, F. Ge, Y. Hong, A review of recent studies on vortex-induced vibrations of long slender cylinders, Journal of Fluids and Structures 28 (2012) 292–308.
  * [2] C. H. K. Williamson, R. N. Govardhan, Vortex-induced vibrations, Annual Review of Fluid Mechanics 36 (2004) 413–455.
  * [3] J. Wang, D. Fan, K. Lin, A review on flow-induced vibration of offshore circular cylinders, Journal of Hydrodynamics 32 (3) (2020) 415–440.
  * [4] R. Gabbai, H. Benaroya, An overview of modeling and experiments of vortex-induced vibration of circular cylinders, Journal of Sound and Vibration 282 (3-5) (2005) 575–616.
  * [5] C. H. K. Williamson, R. N. Govardhan, A brief review of recent results in vortex-induced vibrations, Journal of Wind engineering and industrial Aerodynamics 96 (6-7) (2008) 713–735.
  * [6] P. Bearman, Circular cylinder wakes and vortex-induced vibrations, Journal of Fluids and Structures 27 (5-6) (2011) 648–658.
  * [7] P. W. Bearman, Vortex shedding from oscillating bluff bodies, Annual Review of Fluid Mechanics 16 (1) (1984) 195–222.
  * [8] J. M. Dahl, F. S. Hover, M. S. Triantafyllou, et al., Resonant vibrations of bluff bodies cause multivortex shedding and high frequency forces, Physical Review Letters 99 (14) (2007) 144503.
  * [9] D. Fan, Z. Wang, M. S. Triantafyllou, et al., Mapping the properties of the vortex-induced vibrations of flexible cylinders in uniform oncoming flow, Journal of Fluid Mechanics 881 (2019) 815–858.
  * [10] J. M. Dahl, F. S. Hover, M. S. Triantafyllou, Two-degree-of-freedom vortex-induced vibrations using a force assisted apparatus, Journal of Fluids and Structures 22 (6) (2006) 807–818.
  * [11] R. Bourguet, Y. Modarres-Sadeghi, G. E. Karniadakis, et al., Wake-body resonance of long flexible structures is dominated by counterclockwise orbits, Physical Review Letters 107 (13) (2011) 134502.
  * [12] D. Fan, M. S. Triantafyllou, Vortex induced vibration of riser with low span to diameter ratio buoyancy modules, in: The 27th International Ocean and Polar Engineering Conference, International Society of Offshore and Polar Engineers, 2017.
  * [13] C. H. K. Williamson, Vortex dynamics in the cylinder wake, Annual Review of Fluid Mechanics 28 (1) (1996) 477–539.
  * [14] T. Sarpkaya, A critical review of the intrinsic nature of vortex-induced vibrations, Journal of Fluids and Structures 19 (4) (2004) 389–447.
  * [15] J. Chaplin, P. Bearman, F. H. Huarte, et al., Laboratory measurements of vortex-induced vibrations of a vertical tension riser in a stepped current, Journal of Fluids and Structures 21 (1) (2005) 3–24.
  * [16] J. Chaplin, P. Bearman, Y. Cheng, et al., Blind predictions of laboratory measurements of vortex-induced vibrations of a tension riser, Journal of Fluids and Structures 21 (1) (2005) 25–40.
  * [17] H. Braaten, H. Lie, NDP riser high mode VIV tests, Norwegian Marine Technology Research Institute, Technical Report (512394.00) (2004) 01.
  * [18] Y. Modarres-Sadeghi, H. Mukundan, J. M. Dahl, et al., The effect of higher harmonic forces on fatigue life of marine risers, Journal of Sound and Vibration 329 (1) (2010) 43–55.
  * [19] Y. Modarres-Sadeghi, F. Chasparis, M. S. Triantafyllou, et al., Chaotic response is a generic feature of vortex-induced vibrations of flexible risers, Journal of Sound and Vibration 330 (11) (2011) 2565–2579.
  * [20] H. Zheng, R. E. Price, Y. Modarres-Sadeghi, et al., On fatigue damage of long flexible cylinders due to the higher harmonic force components and chaotic vortex-induced vibrations, Ocean Engineering 88 (2014) 318–329.
  * [21] J. Vandiver, H. Marcollo, S. Swithenbank, V. Jhingran, et al., High mode number vortex-induced vibration field experiments, in: Offshore Technology Conference, Offshore Technology Conference, 2005.
  * [22] J. K. Vandiver, S. B. Swithenbank, V. Jaiswal, V. Jhingran, Fatigue damage from high mode number vortex-induced vibration, in: Proceedings of 25th International Conference on Ocean, Offshore and Arctic Engineering, 2006.
  * [23] J. K. Vandiver, V. Jaiswal, V. Jhingran, Insights on vortex-induced, traveling waves on long risers, Journal of Fluids and Structures 25 (4) (2009) 641–653.
  * [24] D. Lucor, L. Imas, G. E. Karniadakis, Vortex dislocations and force distribution of long flexible cylinders subjected to sheared flows, Journal of Fluids and Structures 15 (3-4) (2001) 641–650.
  * [25] R. Bourguet, G. E. Karniadakis, M. S. Triantafyllou, Vortex-induced vibrations of a long flexible cylinder in shear flow, Journal of Fluid Mechanics 677 (2011) 342–382.
  * [26] R. Bourguet, G. E. Karniadakis, M. S. Triantafyllou, Distributed lock-in drives broadband vortex-induced vibrations of a long flexible cylinder in shear flow, Journal of Fluid Mechanics 717 (2013) 361–375.
  * [27] H. Zhu, D. Zhou, Y. Bao, R. Wang, et al., Wake characteristics of stationary catenary risers with different incoming flow directions, Ocean Engineering 167 (2018) 142–155.
  * [28] M. S. Triantafyllou, G. Triantafyllou, Y. Tein, et al., Pragmatic riser VIV analysis, in: Offshore technology conference, Offshore Technology Conference, 1999\.
  * [29] C. M. Larsen, K. Vikestad, R. Yttervik, et al., VIVANA theory manual, Marintek, Trondheim, Norway (2001).
  * [30] F. E. Roveri, J. K. Vandiver, Slenderex: Using shear7 for assessment of fatigue damage caused by current induced vibrations, in: Proceedings of the 20th International Conference on Offshore Mechanics and Arctic Engineering, 2001, pp. 3–8.
  * [31] Z. Wang, D. Fan, M. S. Triantafyllou, et al., A large-eddy simulation study on the similarity between free vibrations of a flexible cylinder and forced vibrations of a rigid cylinder, Journal of Fluid and Structures (under review).
  * [32] R. Gopalkrishnan, Vortex-induced forces on oscillating bluff cylinders, Ph.D. thesis, Massachusetts Institute of Technology (1993).
  * [33] D. Fan, G. Jodin, T. R. Consi, et al., A robotic intelligent towing tank for learning complex fluid-structure dynamics, Science Robotics 4 (36) (2019).
  * [34] Y. Xu, S. Fu, Y. Chen, et al., Experimental investigation on vortex induced forces of oscillating cylinder at high Reynolds number, Ocean Systems Engineering 3 (3) (2013) 167–180.
  * [35] Y. Chen, S. Fu, Y. Xu, et al., Hydrodynamic characters of a near-wall circular cylinder oscillating in cross flow direction in steady current, Acta Physica Sinica 62 (6) (2013) 064701.
  * [36] C. C. J. Chang, R. A. Kumar, M. M. Bernitsas, VIV and galloping of single circular cylinder with surface roughness at $3.0\times 10^{4}<\mbox{Re}<1.2\times 10^{5}$, Ocean Engineering 38 (16) (2011) 1713–1732.
  * [37] A. I. Forrester, A. Sóbester, A. J. Keane, Multi-fidelity optimization via surrogate modelling, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 463 (2088) (2007) 3251–3269.
  * [38] P. Perdikaris, M. Raissi, A. Damianou, et al., Nonlinear information fusion algorithms for data-efficient multi-fidelity modelling, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473 (2198) (2017) 20160751.
  * [39] L. Bonfiglio, P. Perdikaris, S. Brizzolara, G. Karniadakis, Multi-fidelity optimization of super-cavitating hydrofoils, Computer Methods in Applied Mechanics and Engineering 332 (2018) 63–85.
  * [40] F. S. Costabal, P. Perdikaris, E. Kuhl, D. E. Hurtado, Multi-fidelity classification using Gaussian processes: accelerating the prediction of large-scale computational models, Computer Methods in Applied Mechanics and Engineering 357 (2019) 112602.
  * [41] K. Tian, Z. Li, L. Huang, K. Du, L. Jiang, B. Wang, Enhanced variable-fidelity surrogate-based optimization framework by Gaussian process regression and fuzzy clustering, Computer Methods in Applied Mechanics and Engineering 366 (2020) 113045.
  * [42] X. Meng, G. E. Karniadakis, A composite neural network that learns from multi-fidelity data: Application to function approximation and inverse PDE problems, Journal of Computational Physics 401 (2020) 109020.
  * [43] X. Zhang, F. Xie, T. Ji, Z. Zhu, Y. Zheng, Multi-fidelity deep neural network surrogate model for aerodynamic shape optimization, Computer Methods in Applied Mechanics and Engineering 373 (2020) 113485.
  * [44] S. Lee, F. Dietrich, G. E. Karniadakis, I. G. Kevrekidis, Linking Gaussian process regression with data-driven manifold embeddings for nonlinear data fusion, Interface focus 9 (3) (2019) 20180083.
  * [45] R. M. Neal, et al., MCMC using Hamiltonian dynamics, Handbook of Markov Chain Monte Carlo 2 (11) (2011) 2.
  * [46] M. Betancourt, A conceptual introduction to Hamiltonian Monte Carlo, arXiv preprint arXiv:1701.02434 (2017).
  * [47] L. Yang, X. Meng, G. E. Karniadakis, B-PINNs: Bayesian physics-informed neural networks for forward and inverse PDE problems with noisy data, Journal of Computational Physics 415 (2021) 109913.