Interpolatory Model Reduction of Parameterized Bilinear
Dynamical Systems ∗

arXiv:1712.07308v1 [math.NA] 20 Dec 2017

Andrea Carracedo Rodriguez†

Serkan Gugercin†

Jeff Borggaard†

Originally submitted for publication in July 2017
Abstract
Interpolatory projection methods for model reduction of nonparametric linear dynamical systems
have been successfully extended to nonparametric bilinear dynamical systems. However, this is not the
case for parametric bilinear systems. In this work, we aim to close this gap by providing a natural
extension of interpolatory projections to model reduction of parametric bilinear dynamical systems. We
introduce necessary conditions that the projection subspaces must satisfy to obtain parametric tangential
interpolation of each subsystem transfer function. These conditions also guarantee that the parameter
sensitivities (Jacobian) of each subsystem transfer function is matched tangentially by those of the
corresponding reduced order model transfer function. Similarly, we obtain conditions for interpolating
the parameter Hessian of the transfer function by including extra vectors in the projection subspaces. As
in the parametric linear case, the basis construction for two-sided projections does not require computing
the Jacobian or the Hessian.

Keywords: Model reduction, Parametric, Bilinear, Interpolation.

1

Introduction

Simulation of dynamical systems has become an essential part in the development of science to study complex
physical phenomena. However, as the ever increasing need for accuracy has lead to ever larger dimensional
dynamical systems, this increased dimension often makes the desired numerical simulations prohibitively
expensive to perform. Model order reduction (MOR) is one remedy for this predicament. MOR tackles this
issue by constructing a much lower dimensional representation of the corresponding full-order dynamical
system, which is cheap to simulate, yet provides high-fidelity, i.e., it provides a good approximation to
the original quantity of interest. In many applications such as optimization, design, control, uncertainty
quantification, and inverse problems, the dynamics of the system are defined by a set of parameters that
describe initial conditions, material properties, etc. Since carrying out model reduction for every parameter
value is not computationally feasible, the goal in the parameterized setting is to construct a parametric
reduced model that can approximate one or more quantities of interest well for the whole parameter range
of interest. This lead to the parametric model reduction framework. For more specific details on both
parametric and nonparametric model reduction, we refer the reader to [4, 5, 10, 13, 15, 33] and the references
therein.
In this paper, we will focus on large-scale bilinear systems parametrized with the parameter vector p ∈ Rν
and represented in state-space form

m
X

 E(p)ẋ(t; p) = A(p)x(t; p) +
Nj (p)x(t)uj (t) + B(p)u(t),
(1)
j=1


y(t; p) = C(p)x(t; p),
∗ Supported in part by the National Science Foundation under contract DMS–1522616 and the National Institute for Occupational Safety and Health under contract 200-2014-59669. The work of Gugercin was also supported in part by the Alexander
von Humboldt Foundation.
† Department of Mathematics, Virginia Tech, Blacksburg, VA, USA.
{crandrea, gugercin, jborggaard}@vt.edu

1

where x(t; p) ∈ Rn , y(t; p) ∈ R` , and u(t) = [u1 (t), u2 (t), . . . , um (t)]> ∈ Rm denote the states, outputs
(measurements/quantities of interest), and inputs (excitation/forcing) of the bilinear dynamical system,
respectively. Thus, the corresponding state-matrices have the dimensions E(p), A(p), Nj (p) ∈ Rn×n , for
j = 1, . . . m, B(p) ∈ Rn×m , and C(p) ∈ R`×n . In this paper, we assume that the matrix E(p) is nonsingular
for every parameter value p ∈ Rν . Bilinear systems of the form (1) appear in a variety of applications such as
the study of biological species and nuclear fission, are used in the context of stochastic control problems, and
frequently appear in modeling nonlinear phenomena of small magnitude, for instance, [17,18,31,36,37,42,44].
We are interested in large-scale settings where simulating/solving (1) for a wide variety of inputs u(t) and
parameters p solely to determine the output y(t; p) is too expensive. Therefore, our goal is to construct a
reduced parametric bilinear system of order r  n in state-space form

m
X

 E(p)
e
e
e j (p)e
e
e˙ (t; p) = A(p)e
x
x(t; p) +
N
x(t; p)uj (t) + B(p)u(t),
e:
Σ
(2)
j=1


e
e (t; p) = C(p)e
y
x(t; p),
e
e
e j (p) ∈ Rr×r , for j = 1, . . . m, B(p)
e
e
where E(p),
A(p),
N
∈ Rr×m , and C(p)
∈ R`×r such that the reduced
e (t; p) provides a good approximation to the original output y(t; p) for a variety of inputs u(t) and
output y
a range of parameters p.
Non-parametric bilinear systems where the state-space matrices E, A, Nj for j = 1, . . . , m, B and C are
constant, have been studied thoroughly, and input-independent/optimal model reduction techniques from the
linear case (Nj = 0 for j = 1, . . . , m) have been successfully generalized to non-parametric bilinear systems.
For example, [1, 8, 16, 20, 41] have extended model reduction via rational interpolation [6, 11] from linear to
non-parametric bilinear systems. The optimal model reduction of linear dynamical systems in the H2 norm
via the iterative rational Krylov algorithm (IRKA) [29] has been generalized to bilinear systems via bilinear
IRKA (B-IRKA) [12]. Later, [26] showed that, as with IRKA and H2 model reduction in the linear case,
the reduced model via B-IRKA also yields a Hermite interpolation in this case, but in the sense of Volterra
series interpolation. Similarly, gramians and balanced truncation (BT) for linear dynamical systems [38, 39]
have also been generalized to nonparametric bilinear systems [2, 17, 27, 31]. Moreover, [7] has applied the
Loewner framework [35] to bilinear systems.
A plethora of work exists on model reduction of parametrized linear dynamical systems, i.e., Nj = 0 for
j = 1, . . . , m in (1); see, for example, [3,9,14,15,21,23,24,30,40] and the references therein. In this paper, we
are interested in input-independent (transfer function-based) model reduction of parametric bilinear systems
where only the state-space matrices enter into the model reduction process and there is no need to choose a
specific input u(t) nor to simulate the full model (1). More specifically, we are focused on parametric model
reduction that uses the concept of (parametric) rational interpolation. These methods, also referred to as
interpolatory parametric model reduction, have been successfully applied to parametric linear dynamical
systems; see, e.g., [9, 14, 23]. However, unlike the extensions of interpolation theory and IRKA to nonparametric bilinear systems, interpolatory methods have not yet been generalized to parametric bilinear
systems. In this paper, we close this gap and provide a natural extension of interpolatory projections to
parametric bilinear dynamical systems. Our framework yields a reduced parametric bilinear model whose
subsystem transfer functions will (tangentially) interpolate the original subsystem transfer functions together
with the parameter sensitivities and Hessians at the sampled frequencies and parameter values along chosen
directions. Note that we are not focusing on the problem of selecting parameter samples, but rather on
ensuring tangential interpolation of the full and reduced models at the chosen points and directions. One
can then couple this interpolatory model reduction algorithm to a desired sampling strategy.
The remainder of the paper is organized as follows: Section 2 introduces the problem description and
presents the main theoretical results. Section 3 illustrates the theory using two numerical examples. This is
followed by the conclusions and future directions in Section 4.

2

Problem Description

In this section, we introduce the ingredients of the model reduction problem for parametric bilinear systems
such as projection, subsystem transfer function, and tangential interpolation. We then present the main
results of the paper.
2

2.1

Projection-based model reduction of parametric bilinear systems via global
basis

We construct the reduced parametric bilinear system (2) via projection. We follow the global basis approach
(as opposed to using a local basis and performing extrapolation [32, 46] or interpolation [3, 24, 45]). Thus
we construct two constant global model reduction bases, namely V ∈ Cn×r and W ∈ Cn×r , that capture
the parametric dependence of the underlying system using the information from various sampling points.
We refer the reader to [15] for detailed explanations regarding global and local bases, and different sampling
options. The subspaces V and W are computed to enforce specific interpolation conditions as discussed in
Section 2.4.
Once the model reduction bases V and W are constructed, the reduced model quantities in (2) are
obtained via Petrov-Galerkin projection:
e
E(p)
= W> E(p)V,
e
C(p) = C(p)V,
and

e = W> A(p)V,
A
e j = W> Nj (p)V
N

e
B(p)
= W> B(p),
for j = 1, . . . , m.

(3)

b ∈ Rν . Consider the
Now consider reevaluating reduced model quantities in (3) for a new parameter value p
>
e p). This will require re-evaluating the projection E(b
e p) = W E(b
case of E(b
p)V where the operations depend
on the original system dimension n. In practice, many problems exhibit an affine parametric structure, which
makes the projection step numerically efficient. For simplicity, continue to consider the matrix E(p) only.
Assume that E(p) has the following affine parametric form
E(p) = E0 +

N
X

fi (p)Ei ,

(4)

i=1

where fi are scalar (nonlinear) functions reflecting the parametric dependency, and Ei ∈ Rn×n for i =
0, . . . , N are constant matrices. Then, the reduced matrix Er (p) in (3) is given by
Er (p) = W> E0 V +

N
X

fi (p)W> Ei V,

(5)

i=1

where W> Ei V, for i = 0, . . . , N have to be computed once in an offline phase then can be recombined
for efficient computation of Er (b
p) in any online phase. The same discussion applies to other matrices in
(3) as well. When E(p) does not admit such an affine parametrization as in (4), one usually performs
an affine approximation of E(p) first, usually via a matrix version of (Discrete) Empirical Interpolation
Method [22,28]; see [15] for details. We will revisit this issue in the second numerical example in Section 3.2.

2.2

Interpolatory projections for parametric linear systems

A powerful framework in the case of linear dynamical systems
E(p)ẋ(t; p) = A(p)x(t; p) + B(p)u(t),

y(t; p) = C(p)x(t; p),

(6)

is to transform the problem into the frequency domain via Laplace transform. To do so, let Y(s; p) and
U(s) denote the Laplace transforms of y(t; p) and u(t), respectively. Then, applying the Laplace transform
to (6) leads to
Y(s; p) = H(s; p)U(s),

where

H(s; p) = C(p) (s E(p) − A(p))

−1

B(p)

is the transfer function of (6). Then, the goal is to construct a reduced parametric linear model
e
e
e
e˙ (t; p) = A(p)e
E(p)
x
x(t; p) + B(p)u(t),

e
e (t; p) = C(p)e
y
x(t; p),

(7)


−1
e p) = C(p)
e
e
e
e
whose reduced parametric transfer function H(s;
s E(p)
− A(p)
B(p)
approximates H(s; p)
e p) = (H(s; p) − H(s;
e p))U(s).
e (t; p) ≈ y(t; p) since Y(s; p) − Y(s;
well, which would in turn imply y
3

e p) ≈ H(s; p) is via rational interpolation: Given the frequency interpolation
One way to enforce H(s;
points {σ1 , . . . , σqs } ⊂ C, the right tangential directions {b1 , . . . , bqs } ⊂ Cm , the left tangential directions
b qp } ⊂ Rν , find a reduced model (7)
{c1 , . . . , cqs } ⊂ C` , and the parameter interpolation samples {b
p1 , . . . , p
e p) is a Hermite tangential interpolant to H(s; p) at the selected samples, i.e.,
such that H(s;
∂
∂s

e i; p
b j )bi = H(σ
H(σi ; p
 b j )bi ,


∂
>
e
b j )bi ,
b j )bi = ∂s c>
ci H(σi ; p
i H(σi ; p

and ∇p

e i; p
b j ) = c>
b j ),
c>
i H(σi ; p
i H(σ


>
e
b j )bi = ∇p c>
b j )bi ,
ci H(σi ; p
i H1 (σi ; p

for i = 1, . . . , qs and j = 1, . . . , qp . In other words, the reduced model tangentially matches transfer function
values in addition to its frequency and parametric derivatives at the sampled points. One can impose higher
order interpolation conditions at the frequency and parameter samples as well, such as the parameter Hessian.
We omit it for brevity here. The following result from [9] shows how to construct model reduction bases V
and W that satisfy the desired interpolation conditions.

−1
−1
e p) = C(p)
e
e
e
e
Theorem 1. Given H(s; p) = C(p) (s E(p) − A(p)) B(p), let H(s;
s E(p)
− A(p)
B(p)
b ∈ Rν , b ∈ Cm \ {0}, and
be obtained via Petrov-Galerkin projection using the bases V and W. Let σ ∈ C, p
c ∈ C` \ {0}. Define
A(s; p) = sE(p) − A(p)

e
e
e p) = sE(p)
A(s;
− A(p).

and

(8)

b )−1 B(b
(a) If A(σ; p
p)b ∈ Ran(V), then
e
b )b = H(σ;
b )b;
H(σ; p
p
b )−> C(b
(b) If A(σ; p
p)> c ∈ Ran(W), then
e
b ) = c> H(σ;
b );
c> H(σ; p
p
(c) If both (a) and (b) hold simultaneously, then


∂  >e
∂
b )b =
b )b
c> H(σ; p
c H(σ; p
∂s
∂s

and




e
b )b = ∇p c> H(σ;
b )b ,
∇p c> H(σ; p
p

e p
b ) and A(σ;
b ) are invertible.
provided A(σ; p
Theorem 1 shows how to construct V and W to fulfill the required interpolation conditions. All one has
b )−1 B(b
to do is to compute the vectors, e.g., the vector A(σ; p
p)b, for the desired frequency interpolation
b , and use these vectors as columns of V. We refer the reader to
points σ and parameter interpolation point p
the original source [9] for more details. The goal of this paper is to extend this result to parametric bilinear
systems.

2.3

Interpolatory parametric bilinear model reduction problem

Re-consider the full-order parametric bilinear system in (1):

m
X

 E(p)ẋ(t; p) = A(p)x(t; p) +
Nj (p)x(t)uj (t) + B(p)u(t),
j=1


y(t; p) = C(p)x(t; p).

(1)

Even though this system is nonlinear, due to the terms involving Nj (p), the concept of transfer function
can still be applied via Volterra series representation [42]. Given the bilinear system (1), we first introduce
some notation to make the presentation of the Volterra series representation more compact:


N1 (p)
 N2 (p) 
k


N(p) = 
I⊗
⊗ · · · ⊗ Im ,
(9)
N(p) = [N1 (p) N2 (p) · · · Nm (p)],
 , and
..
m =I
|m {z
}


.
Nm (p)
4

k times

where ⊗ denotes the Kronecker product.
The output y(t; p) of (1) can be represented as a Volterra series
y(t; p) =

∞ Z
X
k=1

0

t1

Z

t2

Z

tk

···
0

hk (t1 , t2 , . . . , tk ; p) u(t −
0

k
X

!
ti ) ⊗ · · · ⊗ u(t − tk ) dtk · · · dt1 ,

(10)

i=1

where hk (t1 , t2 , . . . , tk ; p)’s are the regular Volterra kernels, also called subsystem kernels. Then, taking
the multivariable Laplace transform of the degree k regular kernel hk leads to the k th subsystem transfer
function:
k−2

Hk (s1 , . . . , sk ; p) = C(p)A(sk ; p)−1 × N(p)[Im ⊗ A(sk−1 ; p)−1 N(p)] · · · [I⊗
m

⊗A(s2 ; p)−1 N(p)]
×

k−1
[I⊗
m

(11)
−1

⊗ A(s1 ; p)

B(p)],

k

where A(s; p) is as defined in (8), and N(p) and I⊗
m are as defined in (9). For details of this analysis, we
refer the reader to [42, 43]. The Volterra series representation of bilinear systems has been successfully used
for interpolation-based input-independent, optimal model reduction of non-parametric bilinear systems; see,
e.g., [12, 26].
Similarly, for the reduced bilinear system (2), define
e
e 1 (p) N
e 2 (p) · · · N
e m (p)].
N(p)
= [N

(12)

Then, the k th subsystem transfer function of the reduced model (2) is given by
k−2
−1 e
e k (s1 , . . . , sk ; p) = C(p)
e
e
e
e k ; p)−1 ×N(p)[I
e
e 2 ; p)−1 N(p)]
H
A(s
N(p)] · · · [I⊗
⊗A(s
m ⊗ A(sk−1 ; p)
m
k−1

⊗
× [Im

(13)

e
e 1 ; p)−1 B(p)].
⊗ A(s

This allows us to formulate the parametric interpolatory model reduction problem in our setting: Given
interpolation frequencies {σ1 , . . . , σq } ⊂ C, nontrivial right direction b ∈ Cm , nontrivial left direction c ∈ C` ,
b ∈ Rν , find V, W such that the reduced model (2) constructed via
and interpolation parameter sample p
projection as in (3) satisfies the following interpolation conditions for any k ∈ {1, . . . , q}:
k−1

k−1
e k (σ1 , . . . , σk ; p
b )(I⊗
⊗ b) = H
⊗ b),
m
>
>e
b ) = c Hk (σ1 , . . . , σk ; p
b ),
c Hk (σ1 , . . . , σk ; p




k−1
k−1
∂
∂
e k (σ1 , . . . , σk ; p
b )(I⊗
b )(I⊗
c> Hk (σ1 , . . . , σk ; p
⊗ b) =
c> H
⊗ b) ,
m
m
∂si
∂si




>
⊗k−1
⊗k−1
e k (σ1 , . . . , σk ; p
b )(Im
b )(Im
Jp c Hk (σ1 , . . . , σk ; p
⊗ b) = Jp c> H
⊗ b) ,




k−1
>e
⊗k−1
b )(I⊗
b
Hp c> Hk (σ1 , . . . , σk ; p
⊗
b)
=
H
c
H
(σ
,
.
.
.
,
σ
;
p
)(I
⊗
b)
,
p
k
1
k
m
m

b )(I⊗
Hk (σ1 , . . . , σk ; p
m

(14)
(15)
i ∈ {1, . . . , k}

(16)
(17)
(18)

where Jp (·) denotes the matrix of sensitivities (Jacobian) and Hp (·) denotes the Hessian (tensor) with
respect to p. In other words, we would like to construct a reduced parametric bilinear system whose leading
subsystems interpolate (both in frequency and parameter space) the corresponding leading subsystems of
the full order parametric model. Note that we are not only enforcing Lagrange interpolation. We require
the reduced model to match the parameter sensitivities and Hessians as well, which is important, especially,
in the setting of reduced models in optimization. Moreover, these conditions can then be generalized for
different reordering of the frequencies, multiple tangential directions, and several parameter values.

2.4

Subspace conditions for parametric bilinear interpolation

In this section, we establish the subspace conditions to enforce the desired interpolation conditions (14)–
(18) for parametric bilinear systems. Note that even for the parametric bilinear system (1) we consider here,
some of these interpolation conditions, e.g, (14), do not involve parameter gradient and/or parameter Hessian
interpolation; and thus can be interpreted as regular tangential bilinear subsystem interpolation for a fixed
5

b as considered in [16]. However, even though our subspace conditions for (14)-(15) will look
parameter p = p
similar to those in [16], we include the corresponding theorem (Theorem 2 below) and its complete proof for
the following reasons. Although [16] considers tangential interpolation for non-parametric bilinear systems,
the tangential interpolation conditions appear differently. In our formulation, tangential directions appear
in Kronecker product form due to the structure of Hk (s1 , s2 , . . . , sk ; p) as defined in (11) and illustrates
that Hk (s1 , s2 , . . . , sk ; p) can be considered to have mk inputs. Our conditions result in regular tangential
interpolation in the subblocks of Hk (s1 , s2 , . . . , sk ; p); details will be given below. Moreover, we provide
different proofs for (14) and (15), which we later use in the proof of Theorem 3. Finally, we include the E(p)
term in the full model. Clearly, the subspaces conditions (17) for matching the parameter gradient and (18)
for matching the parameter Hessian are new and will be fully discussed.
Theorem 2. Let q be the number of subsystems we wish to interpolate. Consider {σ1 , . . . , σq } ⊂ C and
b ∈ Rν such that A(σi ; p
b ) is invertible for all i ∈ {1, . . . , q}. Consider also the nontrivial vectors b ∈ Cm
p
`
and c ∈ C . Define
b )−1 B(b
V1 = A(σ1 ; p
p)b,
b ))
W1 = (A(σq ; p

−>

b )−1 N(b
Vk = A(σk ; p
p)(Im ⊗ Vk−1 ),
>

C(b
p) c,

b ))
Wk = (A(σq+1−k ; p

−>

>

N(b
p) (Im ⊗ Wk−1 ),

for k = 2, . . . , q,

(19)

for k = 2, . . . , q.

(20)

If
q
[

Vk ⊆ Ran(V),

(21)

k=1

then, for k = 1, . . . , q,
k−1

b )(I⊗
Hk (σ1 , . . . , σk ; p
m

k−1
e k (σ1 , . . . , σk ; p
b )(I⊗
⊗ b) = H
⊗ b).
m

(22)

If
q
[

Wk ⊆ Ran(W),

(23)

k=1

then, for k = 1, . . . , q,
e k (σq+1−k , . . . , σq ; p
b ) = c> H
b ).
c> Hk (σq+1−k , . . . , σq ; p

(24)

Proof. Define
e p)−1 W> A(s; p),
P(s; p) = VA(s;
e p)−1 W> ,
Q(s; p) = A(s; p)VA(s;
k−2

fk (s1 , . . . , sk ; p) = A(sk ; p)−1 N(p)(Im ⊗ A(sk−1 ; p)−1 N(p)) · · · (I⊗
m
×

k−1

(I⊗
m

−1

⊗ A(s1 ; p)

⊗ A(s2 ; p)−1 N(p))

B(p)b), and
k−2

gk> (s1 , . . . , sk ; p) = c> C(p)A(sk ; p)−1 N(p)(Im ⊗ A(sk−1 ; p)−1 N(p)) · · · (I⊗
m
k−1

× (I⊗
m

(25)

⊗ A(s2 ; p)−1 N(p))

⊗ A(s1 ; p)−1 ).

Note that P(s; p) is a skew projector onto Ran(V) while Q(s; p) is a skew projector along Ker(W> ). Also
k−1
⊗k−1
note that Hk (s1 , . . . , sk ; p)(I⊗
⊗b) = C(p)fk (s1 , . . . , sk ; p) and c> Hk (s1 , . . . , sk ; p) = gk> (s1 , . . . , sk ; p)(Im
⊗
m
B(p)).
First we prove (22). Suppose (21) holds. We know that (22) is true for k = 1 by Theorem 1. Assume
that the result is true for k − 1; recall k < q. Then using the definitions of P(s; p) and fk (s1 , . . . , sk ; p) from
(25), we obtain
k−1

b )(I⊗
Hk (σ1 , . . . , σk ; p
m

k−1
e k (σ1 , . . . , σk ; p
b )(I⊗
b ))fk (σ1 , . . . , σk ; p
b ), (26)
⊗ b) − H
⊗ b) = C(b
p)(In − P(σk ; p
m

b ) using the right interpolation of Hk−1 (σ1 , . . . , σk−1 ; p
b ) due to
where we factor out the term fk (σ1 , . . . , σk ; p
the induction assumption. Then, what is left to show is that (26) is zero: By the construction of V in (21),
b ) ∈ Ran(V). Hence fk (σ1 , . . . , σk ; p
b ) ∈ Ran(P(σk ; p
b )), which implies
we obtain fk (σ1 , . . . , σk ; p
b ))fk (σ1 , . . . , σk ; p
b ) = 0,
(In − P(σk ; p
6

(27)

since P(s; p) is a skew projector onto Ran(V). The proof of (24) follows similarly. Suppose (23) holds. Once
again, the result is true for k = 1. Assume that it holds for k − 1. Similar to (26), we obtain
k−1

e k (σq+1−k , . . . , σq ; p
b ) − c> H
b ) = gk> (σq+1−k , . . . , σq ; p
b )[I⊗
c> Hk (σq+1−k , . . . , σq ; p
m

b ))B(b
⊗ (In − Q(σq ; p
p)].
(28)
Again we show that this expression is zero. Note that by the definition of gk and the construction of W in
k−1
b ) ⊥ Ker(I⊗
b )); thus
(23), we have gk (σq+1−k , . . . , σq ; p
⊗ Q(σq ; p
m
k−1

b )[I⊗
gk> (σq+1−k , . . . , σq ; p
m

b ))] = 0.
⊗ (In − Q(σq ; p

(29)


Theorem 2 provides tangential interpolation of Hk (s1 , . . . , sk ; p) in a specific order of the frequencies,
namely in the order {σ1 , . . . , σk }. However one might also consider enforcing interpolation at the frequency
samples {σ1 , . . . , σk } in any order, including repetitions, as is considered in [16]. Indeed, as we will show
in Theorem 3, interpolation of the transfer function sensitivities will require this. The result is a direct
extension of Theorem 2; thus we skip the details. It simply requires the subspaces to contain all possible
combinations:
Corollary 1. Let q be the number of subsystems we wish to interpolate. Consider {σ1 , . . . , σq } ⊂ C and
b ∈ Rν such that A(σi ; p
b ) is invertible for all i ∈ {1, . . . , q}. Consider also the nontrivial vectors b ∈ Cm
p
`
and c ∈ C . Define
b )−1 B(b
b )−1 B(b
V1 = [A(σ1 ; p
p)b, · · · , A(σq ; p
p)b],
b )−1 N(b
b )−1 N(b
Vk = [A(σ1 ; p
p)(Im ⊗ Vk−1 ), · · · , A(σq ; p
p)(Im ⊗ Vk−1 )],
b ))
W1 = [(A(σ1 ; p

−>

b ))
C(b
p)> c, · · · , (A(σq ; p

−>

b ))
Wk = [(A(σ1 ; p

−>

b ))
N(b
p)> (Im ⊗ Wk−1 ), · · · , (A(σq ; p

k = 2, . . . , q,
(30)

C(b
p)> c]
−>

N(b
p)> (Im ⊗ Wk−1 )],

k = 2, . . . , q.

If
q
[

Vk ⊆ Ran(V),

(31)

k=1

then, for k = 1, . . . , q, and for any i1 , . . . , ik ∈ {1, . . . , q},
k−1

b )(I⊗
Hk (σi1 , . . . , σik ; p
m

k−1
e k (σi , . . . , σi ; p
b )(I⊗
⊗ b) = H
⊗ b).
m
1
k

If
q
[

Wk ⊆ Ran(W),

(32)

k=1

then, for k = 1, . . . , q, and for any i1 , . . . , ik ∈ {1, . . . , q},
e k (σi , . . . , σi ; p
b ) = cT H
b ).
cT Hk (σi1 , . . . , σik ; p
1
k
So far, we proved the interpolation conditions using either only V or only W; i.e., we assumed interpolation information only in one of the subspaces, considering one-sided projection. The next theorem shows that
when both subspaces are considered one automatically matches the sensitivities (derivatives) with respect
to the frequencies and parameter, indeed without computing the sensitivities to be matched.
Theorem 3. Assume the hypotheses of Corollary 1. Let Vk and Wk be constructed as in (30) for k =
1, 2, . . . , q. If both
q
q
[
[
Vk ⊆ Ran(V)
and
Wk ⊆ Ran(W),
k=1

k=1

then for k = 1, . . . , q and for i = 1, . . . , k:


k−1
∂  >e
∂  >
⊗k−1
b )(I⊗
b )(Im
c Hk (σ1 , . . . , σk ; p
⊗ b) =
c Hk (σ1 , . . . , σk ; p
⊗ b) ,
m
∂si
∂si
7

and





k−1
>e
⊗k−1
b )(I⊗
b
Jp c> Hk (σ1 , . . . , σk ; p
⊗
b)
=
J
c
H
(σ
,
.
.
.
,
σ
;
p
)(I
⊗
b)
.
k
1
k
p
m
m

where Jp (·) denotes the matrix of sensitivities (Jacobian) with respect to p.
Proof. Recall the definitions P(s; p), Q(s; p), fk (s1 , . . . , sk ; p), and gk> (s1 , . . . , sk ; p) from (25). Similarly,
define
k−2
−1 e
e
e
e
e k ; p)−1 N(p)(I
e
e 2 ; p)−1 N(p))
fk (s1 , . . . , sk ; p) = A(s
N(p)) · · · (I⊗
⊗ A(s
m ⊗ A(sk−1 ; p)
m
k−1

× (I⊗
m

e
e 1 ; p)−1 B(p)b),
⊗ A(s
and

k−2
−1 e
e
e
e
e k ; p)−1 N(p)(I
e
e 2 ; p)−1 N(p))
ek> (s1 , . . . , sk ; p) = c> C(p)
g
A(s
N(p)) · · · (I⊗
⊗ A(s
m ⊗ A(sk−1 ; p)
m
k−1

× (I⊗
m

e 1 ; p)−1 ).
⊗ A(s

Let k ∈ {1, . . . , q} and i ∈ {1, . . . , k}. Under the assumptions of the theorem, we know (27) and (29) are
satisfied for any choice of frequencies due to Corollary 1; in particular,
b ))fκ (σ1 , . . . , σκ ; p
b ) = 0,
(In − P(σκ ; p

and

κ−ι

>
b )[I⊗
(σι , . . . , σκ ; p
gκ−ι+1
m

b ))] = 0,
⊗ (In − Q(σκ ; p

for any κ ∈ {1, . . . , q} and ι < κ. Recalling the definitions of the full and reduced transfer functions (11)
and (13), together with (26) and (28), and our definitions in this proof, we can rewrite these two terms as
κ−ι
>
>
b ) − Ve
b ) = 0, gκ−ι+1
b) − g
eκ−ι+1
b )(I⊗
fκ (σ1 , . . . , σκ ; p
fκ (σ1 , . . . , σκ ; p
(σι , . . . , σκ ; p
(σι , . . . , σκ ; p
⊗W> ) = 0,
m

or equivalently
κ−ι

>
>
b ) = Ve
b ) and gκ−ι+1
b) = g
eκ−ι+1
b )(I⊗
fκ (σ1 , . . . , σκ ; p
fκ (σ1 , . . . , σκ ; p
(σι , . . . , σκ ; p
(σι , . . . , σκ ; p
m

⊗ W> ).
(33)

e
Now fix k ∈ {1, . . . , q} and i ∈ {1, . . . , k}, and recall that E(p)
= W> E(p)V. If i 6= 1, then

k−1
k−i
∂  >e
>
b )(I⊗
ek−i+1
b )(I⊗
b)
c Hk (σ1 , . . . , σk ; p
⊗ b) = g
(σi , . . . , σk ; p
⊗ W> )E(b
p)Ve
fi (σ1 , . . . , σi ; p
m
m
∂si
>
b )E(b
b)
= gk−i+1
(σi , . . . , σk ; p
p)fi (σ1 , . . . , σi ; p
(by (33))

k−1
∂  >
b )(I⊗
c Hk (σ1 , . . . , σk ; p
⊗ b) .
=
m
∂si
If i = 1, then
k−1
k−1
∂ >e
b )(I⊗
ek> (σ1 , . . . , σk ; p
b )(I⊗
b)
c Hk (σ1 , . . . , σk ; p
⊗ b) = g
⊗ W> )E(b
p)Ve
f1 (σ1 ; p
m
m
∂si
b )E(b
b)
= gk> (σ1 , . . . , σk ; p
p)f1 (σ1 ; p
k−1
∂ >
b )(I⊗
c Hk (σ1 , . . . , σk ; p
⊗ b).
=
m
∂si

(by (33))

Similarly, we can justify interpolation of the parameter gradients. Since the expression of the parameter
gradient for a general subsystem transfer function Hk (s1 , . . . , sk ; p) becomes too involved to properly present
in a single page, we provide the proof for the second subsystem (the result for the first subsystem follows
Theorem 1) and sketch the proof for a general subsystem. Let pj refer to any entry of the parameter vector
p ∈ Rν . Consider

∂  >
e 2 (σ1 , σ2 ; p
b )(Im ⊗ b) − c> H
b )(Im ⊗ b) .
c H2 (σ1 , σ2 ; p
(34)
∂pj
Let Mpj (p) denote the partial derivative of M(p) with respect to pj . Then, by taking the partial derivae p (p) = Cp (p)V,
tives in (34), using interpolation of the first subsystem, rearranging terms and using C
j
j

8

e p (p) = W> Np (p)(Im ⊗ V), and B
e p (p) = W> Bp (p), one can write
Aepj (s; p) = W> Apj (s; p)V, N
j
j
j
j

∂  >
e 2 (σ1 , σ2 ; p
b )(Im ⊗ b) − c> H
b )(Im ⊗ b)
c H2 (σ1 , σ2 ; p
∂pj
b )Apj (σ2 ; p
b ))(I − P(σ2 ; p
b ))f2 (σ1 , σ2 ; p
b)
= (c> Cpj (b
p) − g1 (σ2 ; p

b ) Im ⊗ (I − Q(σ1 ; p
b ))(Bpj (b
b )f1 ) ,
+ g2> (σ1 , σ2 ; p
p)b − Apj (σ1 ; p
which can be justified by multiplying out the right-hand side and re-grouping. We know that in the secondb ))f2 (σ1 , σ2 ; p
b ) = 0 and in the third-line of this expression, we
line of this expression, we obtain (I − P(σ2 ; p
b ) (Im ⊗ (I − Q(σ1 ; p
b ))) = 0 using k = 2 in the proof of Theorem 2. Thus, we have
obtain g2> (σ1 , σ2 ; p


∂  >e
∂
b )(Im ⊗ b) =
b )(Im ⊗ b) .
c> H2 (σ1 , σ2 ; p
c H2 (σ1 , σ2 ; p
∂pj
∂pj



e 2 (σ1 , σ2 ; p
b )(Im ⊗ b) = Jp c> H
b )(Im ⊗ b)
Since pj was an arbitrary entry of p, this yields Jp c> H2 (σ1 , σ2 ; p
as desired. Now, for the general case, consider

k−1
∂  >e
b )(I⊗
b)
c Hk (σ1 , . . . , σk ; p
⊗
b)
= c> Cpj (b
p)Ve
fk (σ1 , . . . , σk ; p
(35)
m
∂pj
e> (σk ; p
b )W> Ap (σk ; p
b )Ve
b)
+g
fk (σ1 , . . . , σk ; p
(36)
1

+

j

e1> (σk ; p
b )W> Npj (b
g
p)(Im

b)
⊗ V)e
fk−1 (σ1 , . . . , σk−1 ; p
k−1

ek> (σ1 , . . . , σk ; p
b )(I⊗
+ ··· + g
m

⊗ W> Bpj (b
p)b).

(37)
(38)

b ) with fk (σ1 , . . . , σk ; p
b ).
Consider the right-hand side of (35). Using (33), one can replace Ve
fk (σ1 , . . . , σk ; p
e1> (σk ; p
b )W> with g1> (σk ; p
b ). Continuing in this
Similarly, in (36), once again using (33), one replaces g
fashion, we obtain

k−1
∂  >e
b )(I⊗
b)
c Hk (σ1 , . . . , σk ; p
⊗ b) = c> Cpj (b
p)fk (σ1 , . . . , σk ; p
m
∂pj
b )Apj (σk ; p
b )fk (σ1 , . . . , σk ; p
b)
+ g1> (σk ; p
b )Npj (b
b)
+ g1> (σk ; p
p)fk−1 (σ1 , . . . , σk−1 ; p
k−1

⊗
b )(Im
+ · · · + gk> (σ1 , . . . , σk ; p
⊗ Bpj (b
p)b)


k−1
∂
b )(I⊗
c> Hk (σ1 , . . . , σk ; p
⊗ b)
=
m
∂pj

as desired.



We now present the final theoretical result, showing the interpolation of the parameter Hessian. As the
expressions become too involved for a general subsystem transfer function, we write and proof the conditions
for the first and second subsystems only, but the results can be generalized similarly.
Theorem 4. Assume the hypotheses of Corollary 1 for q = 2. Define


V1 (p) = A(σ1 ; p)−1 B(p)b, A(σ2 ; p)−1 B(p)b ,


V2 (p) = A(σ1 ; p)−1 N(p)(Im ⊗ V1 (p)), A(σ2 ; p)−1 N(p)(Im ⊗ V1 (p)) ,
h
i
−>
−>
W1 (p) = (A(σ1 ; p)) C(p)> c, (A(σ2 ; p)) C(p)> c ,
h
i
−>
−>
W2 (p) = (A(σ1 ; p)) N(p)> (Im ⊗ W1 (p)), (A(σ2 ; p)) N(p)> (Im ⊗ W1 (p)) .
Assume

2
[

Vk (b
p) ⊆ Ran(V),

and

k=1

2
[
k=1

9

Wk (b
p) ⊆ Ran(W).

(39)

If either
2 [
ν
[
∂
Vk (b
p) ⊆ Ran(V)
∂p
j
j=1

or

k=1

then

2 [
ν
[
∂
Wk (b
p) ⊆ Ran(W),
∂p
j
j=1

(40)

k=1




e 2 (σ1 , σ2 ; p
b )(Im ⊗ b) = Hp c> H
b )(Im ⊗ b) .
Hp c> H2 (σ1 , σ2 ; p

where Hp (·) denotes the Hessian with respect to p.

b )b
Proof. Assume we have the extra conditions on V. First note that we have interpolation of Hp c> H1 (σi ; p
for i ∈ {1, 2} since this is the linear case (see [9] for details). Let pi and pj refer to any entries in the parameter
vector p. Recall the definition of the second transfer function of the full and reduced model, respectively,
H2 (s1 , s2 ; p) = C(p)A(s1 ; p)−1 N(p)(Im ⊗ A(s1 ; p)−1 B(p)),
−1 e
e 2 (s1 , s2 ; p) = C(p)
e
e
e 1 ; p)−1 N(p)(I
e
H
A(s
B(p)).
m ⊗ A(s1 ; p)
e 2 (s1 , s2 ; p) with respect to pj and pi , apply the definition of the
We take the second partial derivative of H
reduced order matrices, rearrange the terms and use the notation in previous proofs to obtain

∂2  > e
b )(Im ⊗ b)
c H2 (σ1 , σ2 ; p
∂pj ∂pi
b)
= c> Cp p (b
p)Ve
f2 (σ1 , σ2 ; p

(41)

j i

+

e2> (σ1 , σ2 ; p
b )W> Bpj pi (b
g
p)b

(42)

e1> (σ2 ; p
b )W> Npj pi (b
g
p)(Im

(43)

b ))
⊗ Ve
f1 (σ1 ; p
e1> (σ2 ; p
b )W> Apj pi (σ2 ; p
b )Ve
b)
+g
f2 (σ1 , σ2 ; p
>
>
e2 (σ1 , σ2 ; p
b )(Im ⊗ W Apj pi (σ1 ; p
b )Ve
b ))
+g
f1 (σ1 ; p

+

b )W> Npj (b
{e
g1> (σ2 ; p
p)

e2> (σ1 , σ2 ; p
b )W> Apj (σ1 ; p
b )}(Im
g

b ))
⊗ V[e
f1 ]pi (σ1 ; p
b )W> Npi (b
e2> (σ1 , σ2 ; p
b )W> Api (σ1 ; p
b )}(Im ⊗ V[e
b ))
+ {e
g1> (σ2 ; p
p) + g
f1 ]pj (σ1 ; p
>
>
>
e1 (σ2 ; p
b )W Apj (σ2 ; p
b )}V[e
b)
+ {c Cpj (b
p) + g
f2 ]pi (σ1 , σ2 ; p
e1> (σ2 ; p
b )W> Api (σ2 ; p
b )}V[e
b ),
+ {c> Cpi (b
p) + g
f2 ]pj (σ1 , σ2 ; p

+

+

(44)
(45)
(46)
(47)
(48)
(49)

where Mpj pi (p) denotes the second partial derivative of M(p) with respect to pj and pi , and [fk ]pi denotes
the partial derivative of fk with respect pi . Then, we follow the similar manipulations used in the proof of
Theorem 3 for (35)-(38): Equations (41)-(45) contain the same terms, and thus we follow the same reasonings
that we used for (35)-(38).
b ) and [f2 ]pj (σ1 , σ2 ; p
b ), the same manipuThen, even though (46)-(49) contain the new terms [f1 ]pj (σ1 ; p
b ) and [f2 ]pj (σ1 , σ2 ; p
b ) are
lations still apply here due to the construction of V in (39) and (40), [f1 ]pj (σ1 ; p
now also spanned by Ran(V) for any pj . Therefore, we obtain


∂2
∂2  > e
b )(Im ⊗ b) =
b )(Im ⊗ b) .
c H2 (σ1 , σ2 ; p
c> H2 (σ1 , σ2 ; p
∂pj ∂pi
∂pj ∂pi
Since pj and pi were arbitrary, we obtain the Hessian matching as desired. The proof would be analogous if
we assumed the extra conditions on W instead. Only the rearrangement of the terms would change so that
b ) and [g2 ]pj (σ1 , σ2 ; p
b ) instead.
the expression depends on [g1 ]pj (σ2 ; p

Remark 1. As we stated above, one can write the conditions for matching the parameter Hessian of the
higher index subsystems. Let q be the number of subsystems we wish to interpolate. To obtain the parameter
Hessian matching for the general k th order subsystem, i.e., to satisfy




k−1
>e
⊗k−1
b )(I⊗
b
Hp c> Hk (σ1 , . . . , σk ; p
⊗
b)
=
H
c
H
(σ
,
.
.
.
,
σ
;
p
)(I
⊗
b)
,
k = 1, . . . , q
p
k 1
k
m
m
10

we would need


V1 (p) = A(σ1 ; p)−1 B(p)b, · · · , A(σq ; p)−1 B(p)b ,


Vk (p) = A(σ1 ; p)−1 N(p)(Im ⊗ Vk−1 (p)), · · · , A(σq ; p)−1 N(p)(Im ⊗ Vk−1 (p)) ,
h
i
−>
−>
W1 (p) = (A(σ1 ; p)) C(p)> c, · · · , (A(σq ; p)) C(p)> c ,
i
h
−>
−>
Wk (p) = (A(σ1 ; p)) N(p)> (Im ⊗ W1 (p)), · · · , (A(σq ; p)) N(p)> (Im ⊗ W1 (p))

k = 1, . . . , q

k = 1, . . . , q.

b ) to be contained in the ranges of the basis V and W, respectively, together with either the
(evaluated at p
partial derivatives of the Vk ’s or the partial derivatives of the Wk ’s with respect to the parameter entries
b ).
(evaluated at p

3

Numerical Examples

In this section, we illustrate the theoretical discussion from Section 2 using two examples: A nonlinear RC
b (i) (or
circuit in Section 3.1 and an advection-diffusion equation in Section 3.2. Throughout this section, p
(i)
p̂ when the parameter is a scalar) denotes the parameter sampling points we used in constructing the
model reduction bases V and W, and p(i) (or p(i) ) denotes the parameter points (which are not sampled)
at which we evaluate both reduced and full models to investigate the accuracy of the reduced model.

3.1

A nonlinear RC circuit

We begin with a modified version of a standard benchmark problem for bilinear systems, namely a nonlinear
RC circuit [8, 41]. The original benchmark problem leads to a non-parametric bilinear system. We have
revised the problem to add parametric dependence. To clearly motivate this parametric dependence, we
include details of the model derivation.
Consider the following SISO parametric nonlinear system

v̇(t; p) = f (v(t); p) + bu(t)
(50)
y(t; p) = c> v(t; p),
>

where v(t; p) ∈ Rn , b = c = [1 0 · · · 0] ∈ Rn ,

−g(v1 ; p) − g(v1 − v2 ; p)

g(v
1 − v2 ; p) − g(v2 − v3 ; p)


..

.
f (v; p) = 
 g(vk−1 − vk ; p) − g(vk − vk+1 ; p)


..

.
g(vN −1 − vN ; p)






,





(51)

and
g(v; p) = epv + v − 1.

(52)

System (50) models a nonlinear RC circuit with N resistors where the state variable v(t; p) is the voltage
at each node, u(t) is the input signal to the current source, the ouput y(t; p) is the voltage between node
1 and ground, and g(ν; p) gives the current-voltage dependency at each resistor. We have introduced a
parameter dependency p ∈ R in the exponential term of this current-voltage dependency, which models the
influence of the operating temperature on the current. Following [8, 41], we apply Carleman bilinearization
to f (v; p) ≈ A1 (p)v + A2 (p)(v ⊗ v) and a second-order approximation of g(v; p) ≈ (p + 1)v + 21 p2 v 2 leading
to an approximation of the nonlinear dynamics (50) by the following parametric bilinear system:

Eẋ(t; p) = A(p)x(t; p) + Nx(t; p)u(t) + bu(t)
(53)
yb (t; p) = c> x(t; p),
11

where

x(t; p) =

A(p) =

v(t; p)
v(t; p) ⊗ v(t; p)
A1 (p)
0


,

E = In ,

A2 (p)
,
A1 (p) ⊗ I + I ⊗ A1 (p)



1
c=b=
,
0

0
N=
b⊗I+I⊗b

(54)
0
0


,

(55)

and we use yb (t; p) to denote the output of the full-order parametric bilinear system. Note that the dimension
2
of the bilinear system is n = N + N 2 and the matrices A1 (p) ∈ RN ×N and A2 (p) ∈ RN ×N are given by


−2 1
 1 −2 1





.
.
.
..
..
..
A1 (p) = (1 + p) 
,



1 −2 1 
1 −1
and, for k = 2, . . . , N − 1,
[A2 (p)](1,1) = −p2 ,
[A2 (p)](1,2) = [A2 (p)](1,N +1) = [A2 (p)](k,(k−2)N +k−1) = [A2 (p)](k,(k−1)N +k+1) =
= [A2 (p)](k,kN +k) = [A2 (p)](N,(N −2)N +N −1) = [A2 (p)](N,(N −1)N +N ) =

p2
,
2

[A2 (p)](1,N +2) = [A2 (p)](k,(k−2)N +k) = [A2 (p)](k,(k−1)N +k−1) = [A2 (p)](k,kN +k+1) =
= [A2 (p)](N,(N −2)N +N ) = [A2 (p)](N,(N −1)N +N −1) = −

p2
,
2

where [A2 (p)](i,j) denotes the (i, j)th entry of A2 (p). Note that both A1 (p) and A2 (p) have the desired affine
structure (4) with the nonlinear scalar parametric functions −p2 , p12 , and − p12 . As in the original benchmark
problem, we choose N = 200, and thus obtain a parametric bilinear system of dimension n = 40, 200. We
are interested in the parameter range p ∈ [0, 70], and choose two parameter sampling points, p̂(1) = 1 and
p̂(2) = 50. For each sampling point, we focus on the leading q = 2 subsystems. We choose {σ1 , σ2 } by running
IRKA on the linearized model (by setting N = 0); i.e, {σ1 , σ2 } correspond to optimal sampling points for the
linear model. With these frequencies, we construct the basis V1 and W1 (using Theorem 4) that guarantees
interpolation of H1 (s; p), H2 (s1 , s2 ; p), and their sensitives for p = p̂(1) = 1, and at {σ1 , σ2 }. Similarly, we
construct V2 and W2 for p̂(2) = 50. We then construct the global bases V = [V1 V2 ] and W = [W1 W2 ]
and obtain a reduced parametric bilinear model of dimension r = 12 using the projection described in (3);
thus we are approximating a parametric bilinear system of dimension n = 40, 200 by a reduced parametric
bilinear model of dimension r = 12. To test the accuracy of the parametric reduced model, we simulate
and compare the outputs of the original nonlinear model (50), the full bilinear model (53), and the reduced
bilinear model for two different inputs, u(t) = e−t and u(t) = 21 (cos(5πt) + 1), and three different parameter
values, p(1) = 18, p(2) = 40, and p(3) = 62. Note that these parameter values are not the sampled values;
indeed p3 = 62 is even outside the sampling range [1, 50]. Moreover, note that the inputs u(t) = e−t and
u(t) = 21 (cos(5πt) + 1) were not used in the model reduction step, i.e., the reduced model is not informed
by these choices of excitation. As Figure 1 shows, the parametric reduced bilinear system provides a very
accurate approximation to the full bilinear model; their responses are almost indistinguishable. Relative L2
errors in the outputs for our three parameter values p(1) , p(2) , p(3) are listed in Table 1, showing a relative
error on the order of 10−3 . We also emphasize that the only deviations visible in Figure 1 are deviations
from the original nonlinear system, due to Carleman bilinearization, and not due to the model reduction
step. We also note that even though the responses might look similar for different parameter values, the
scales of the outputs are different. To make the numerical investigations more detailed, we performed a
parameter sweep using 103 linearly sampled points in the interval [0, 70], then identified the performance of
the reduced model measured in terms of the relative L2 error in the output yb (t; p) for each of the inputs.
We found that for the first input u(t) = e−t , in the worst-case scenario, the reduced model led to a relative
L2 output error of 6.31 × 10−3 . For the second input u(t) = 12 (cos(5πt) + 1), the worst performance yielded
12

Input

p = p(1)

p = p(2)

p = p(3)

u(t) = e−t

2.54 × 10−3

2.91 × 10−3

1.53 × 10−3

u(t) = 12 (cos(5πt) + 1)

2.54 × 10−3

4.33 × 10−3

4.57 × 10−3

Table 1: Relative L2 output error
a relative L2 output error of 5.96 × 10−3 . These numbers further illustrate the ability of the reduced model
to accurately approximate the full order parametric bilinear system.

3.2

Advection-diffusion equation

For our second example consider a model of the transport and diffusion of a passive scalar field T (representing
a chemical concentration, temperature, etc.) on the domain Ω = [−1, 1] × [−1, 1]. The transport of T is
controlled using a background velocity field described with two input parameters u1 and u2 , and two velocity
fields v1 and v2 . Thus the background velocity field is v(x, y) = u1 (t)v1 (x, y) + u2 (t)v2 (x, y). The value of
the passive scalar on the boundary of Ω (∂Ω) is controlled by an input u3 . We model the diffusion using the
viscosity parameter p1 , and include a source term centered at (p2 , p3 ) ∈ Ω with an area of affect described
by p4 given by


(x − p2 )2 + (y − p3 )2
f (x, y; p2 , p3 , p4 ) = exp −
.
p4
The strength of the source term is controlled by an input u4 .
Our passive scalar field T then satisfies
Ṫ (x, y, t) = p1 ∆T (x, y, t) − v · ∇T (x, y, t) + u4 (t)f (x, y; p2 , p3 , p4 ),

(x, y, t) ∈ Ω × (0, ∞),
(x, y) ∈ Ω,

T (x, y, 0) = T0 (x, y),

(x, y, t) ∈ ∂Ω × (0, ∞).

T (x, y, t) = u3 (t),

(56)
(57)
(58)

Thus our model depends on the parameter vector



p1
 p2 

p=
 p3  .
p4
We will consider the following as our parameter range
−3 ≤ ln p1 ≤ 1,

(p2 , p3 ) ∈ Ω,

1 ≤ p4 ≤ 10.

PN
We approximate solutions to (56) using a finite element discretization TN (x, y, t) = j=1 xj (t)ϕj (x, y), where
the {ϕj }N
j=1 arise from quadratic (P2) triangular elements. For convenience, we will split the summation
above into two disjoint parts, one with indices corresponding to boundary nodes (B) and the remainder
corresponding to interior nodes (I). Thus {1, 2, . . . , N } = B ∪ I. Upon substituting this into the weak form
of (56) and suppressing function arguments, we arrive at

 



 



˙
N
N
N
X
X
X


xj ϕj , ϕi  = − p1 ∇ 
xj ϕj  , ∇ϕi  − v · ∇ 
xj ϕj  , ϕi  + (u4 f, ϕi ) ,
∀i ∈ I,

j=1

j=1

j=1

where the boundary integrals vanish since ϕi are zero on the boundary when i ∈ I. Interchanging integration
in the L2 -inner products with the summation leads to
X
X
X
(ϕj , ϕi ) ẋj = −p1
(∇ϕj , ∇ϕi ) xj − p1 u3
(∇ϕj , ∇ϕi )
j∈I

j∈I

− u1

X

j∈B

(v1 · ∇ϕj , ϕi ) xj − u2

j∈I

X
j∈I

13

(v2 · ∇ϕj , ϕi ) xj + u4 (f, ϕi ) ,

p = 18

0.03

p = 18

0.035
original
full
reduced

0.025

original
full
reduced

0.03

0.025
0.02

y(t)

y(t)

0.02
0.015

0.015
0.01
0.01
0.005

0.005

0

0
0

0.5

1

1.5

2

0

t (seconds)
p = 40

0.014

0.5

0.016
original
full
reduced

0.012

1

1.5

2

t (seconds)
p = 40
original
full
reduced

0.014
0.012

0.01

0.01

y(t)

y(t)

0.008
0.008

0.006
0.006
0.004

0.004

0.002

0.002

0

0
0

0.5

1

1.5

2

0

t (seconds)
p = 62

0.01

0.5

0.012
original
full
reduced

0.009

1

1.5

original
full
reduced

0.01

0.008
0.007

2

t (seconds)
p = 62

0.008

y(t)

y(t)

0.006
0.005

0.006

0.004
0.004

0.003
0.002

0.002

0.001
0

0
0

0.5

1

1.5

2

0

t (seconds)

0.5

1

1.5

2

t (seconds)

Figure 1: Solution to (50) (denoted by “original”), (53) (denoted by “full”), and reduced order model
(denoted by “reduced”) for different inputs and parameter values. Left column: u(t) = e−t . Right column:
u(t) = 12 (cos(5πt) + 1).

14

for each i ∈ I. Letting x(t) = [x1 (t) x2 (t) . . . x|I| (t)]> and
[N1 ]ij = − (v1 · ∇ϕj , ϕi ) ,

[E]ij = (ϕi , ϕj ) ,

[b3 ]i = −p1

X

(∇ϕi , ∇ϕk ) ,

k∈B

[A]ij = −p1 (∇ϕi , ∇ϕj ) ,

[N2 ]ij = − (v2 · ∇ϕj , ϕi ) ,

[b4 ]i = (f (·; p2 , p3 , p4 ), ϕi (·)) ,

for i, j ∈ I. We can write our discrete problem as
Eẋ(t; p) = A(p)x(t; p) + N1 x(t; p)u1 (t) + N2 x(t; p)u2 (t) + b3 (p)u3 (t) + b4 (p)u4 (t),
or
Eẋ(t; p) = A(p)x(t; p) +

4
X

Ni x(t; p)ui (t) + B(p)u(t),

i=1

where
N = [N1 N2 N3 N4 ] = [N1 N2 0 0],
B(p) = [0 0 b3 (p) b4 (p)], and
u(t) = [u1 (t) u2 (t) u3 (t) u4 (t)]> .
We also include an output y(t; p) = c> x(t; p) that represents the average of our scalar field over [0.5, 1] ×
[0.5, 1]. In summary, we have a bilinear parametric multi-input/single-output system

P4
Eẋ(t; p) = A(p)x(t; p) + i=1 Ni x(t; p)ui (t) + B(p)u(t)
(59)
y(t; p) = c> x(t; p),
which can be reduced using the strategy presented in the previous example.
For our simulations, we chose a 21-by-21 FEM mesh (which results in a FOM of dimension n = |I| = 361)
and velocity with v1 (x, y) = [−y, x]> and v2 (x, y) = 21 (cos(π(x − y)) + 1)[1, 1]> .
Note that the full order matrices E, N1 , N2 , c are constant, [A(p)]ij = −p1 [A]ij , and [b3 (p)]i = −p1 [b3 ]i ,
hence the ROM is given by

e 3 (p)u3 (t) + b
e 4 (p)u4 (t)
ex
e
e 1x
e 2x
e˙ (t; p) = A(p)e
e(t; p)u1 (t) + N
e(t; p)u2 (t) + b
E
x(t; p) + N
(60)
>
e(t; p)
ye(t; p) = e
c x
where
e = W> EV
E
e j = W> Nj V
N
>

e
A(p)
= −p1 W> AV
e 3 (p) = −p1 W> b3
b

>

>

e
c =c V

e 4 (p) = W b4 (p).
b

(61)
(62)
(63)

Even though the dimension in (60) is lower, the reduction of the vector b4 (p) can not be done offline
(as for the rest of the system matrices), hence we aim to reduce the cost of computing b4 (p) by means of
DEIM approximation as we discussed in Section 2.1. In this case, since b4 (p) is a vector, there is no need
for a matrix-version and the original DEIM formulation suffices. Applying DEIM, we want to find a basis
U ∈ Rn×M where M  n and a row selector S so that
b4 (p) ≈ U(S> U)−1 S> b4 (p)

e 4 (p) ≈ W> U(S> U)−1 S> b4 (p)
and b

are good approximations. This way W> U(S> U)−1 can be precomputed offline, while the online computation
of S> b4 (p) will now only require us to compute the entries in b4 (p) indicated by S. Clearly the accuracy of
this approximation depends on U and S.
We first find U by Proper Orthogonal Decomposition (POD) [19, 34]. That is, we generate a matrix of
snapshots of the vector b4 (p) and select the leading M left singular vectors to be the columns of U. By
15

viscosity

source center

source reach

0.1
1
e−3
e−3
0.0529
0.2392
0.1261
0.0754

(0.25, 0.8)
(0, 0)
(1, 1)
(−0.5, −1)
(0.975, 0.9275)
(0.6914, 0.3149)
(−0.7224, −0.7623)
(−0.3214, 0.4988)

1
9
4
1
1.6636
3.6730
5.1100
2.816

b (1)
p
b (2)
p
b (3)
p
b (4)
p
p(1)
p(2)
p(3)
p(4)

Table 2: Advection-diffusion model. Parameter values.
taking enough snapshots and singular vectors, we expect the range of our basis U to represent the values of
b4 (p) over the parameter domain. We chose a tolerance of 10−5 to truncate the singular values in the POD
basis, resulting in a DEIM approximation of order M = 33. To chose the interpolation indices (row selector)
in S, we use the Q-DEIM algorithm [25] which determines S using a pivoted QR factorization of U> .
In Figure 2 we show the relative error of the Q-DEIM approximation of b4 (p) over 104 random parameter
values in the entire parameter domain. Note that the maximum relative error is on the order of 10−4 , showing
e 4 (p) ≈ W> U(S> U)−1 S> b4 (p) in
the accuracy of the DEIM approximation. Thus, we can confidently use b
our reduced model.
×10-4

4

3.5

10
3

p4

2.5

5

2
1.5

0
1

1

1
0

p3

0
-1

-1

0.5
0

p2

Figure 2: Relative error in the DEIM approximation of b4 (p).

b (i) for i = 1, 2, 3, 4 (see the leading four
To construct our ROM, we sample at four parameter values p
rows in Table 2). We calculate the corresponding projection matrices V1 , V2 , V3 , V4 , W1 , W2 , W3 , and W4
using Theorem 4 that guarantee interpolation and sensitivity matching at the frequency interpolation points
(generated via IRKA once again) and tangential directions corresponding to each parameter value sampled.
e and A,
e we concatenate all of the projection matrices and consider a one-sided
To maintain symmetry in E
projection, i.e., V = [V1 V2 V3 V4 W1 W2 W3 W4 ] and W = V. We truncate the basis (using SVD) and
obtain a ROM with dimension r = 20.
To illustrate the accuracy of the reduced model, we test it for two different inputs sets (see Table 3)
and for four different parameter samples p(i) for i = 1, 2, 3, 4 (see Table 2, rows 5–9) that were not part
of the sampling set. We show the results, the full-order and reduced-order outputs in Figures 3 and 4 for
two different inputs (see Table 3). Both figures show that for each input selection (neither of which entered
into our transfer function-based model reduction process), the parametric reduced bilinear model provides a
high-quality approximation, only showing slight variations at the parameter values that were not sampled.
16

Input 1
Input 2

u1

u2

u3

u4

sin t
0.5

cos t
0.25

-1
1

0.5
-1

Table 3: Advection-diffusion model. Input values.
p (1)

0

FOM
ROM

-0.02

FOM
ROM

-0.05

temperature

-0.04

temperature

p (2)

0

-0.06
-0.08
-0.1

-0.1

-0.15

-0.12
-0.14

-0.2
0

5

10

15

20

0

5

time
p (3)

0

10

15

20

time
p (4)

0
FOM
ROM

FOM
ROM

temperature

temperature

-0.05

-0.1

-0.05

-0.1

-0.15

-0.2

-0.15
0

5

10

15

20

0

time

5

10

15

20

time

Figure 3: Solution of the full (FOM) and reduced order model (ROM) corresponding to non-sampled parameter values (see values in Table 2) for Input 1 with entires u1 (t) = sin t, u2 (t) = cos t, u3 (t) = −1,
u4 (t) = 0.5.
As in the previous example, to ensure a fair comparison, we performed an exhaustive search via 104
uniform random samples in our full parameter domain (except that we fixed the fourth parameter entry at
p4 = 5 so that we can present the results with a 3-dimensional plot) for both input selections from Table 3.
Out of these 104 parameter selections, in Figure 5, we display in the left-plot the relative errors at every
sampling point and in the right-plot, the outputs for the worst performance of the reduced model for Input
1. Note that even for the worst parameter sample, the parametric reduced model still provides an accurate
approximation with a relative L2 error of 4.79 × 10−2 . We repeat the procedure for Input 2 in Figure 6 and
obtain similar results.

17

p (1)

0.04

p (2)

0.15
FOM
ROM

FOM
ROM

temperature

temperature

0.03

0.02

0.1

0.05

0.01

0

0
0

5

10

15

20

0

5

time
p (3)

0.14

15

20

p (4)

0.08
FOM
ROM

0.12

FOM
ROM

0.06

temperature

0.1

temperature

10

time

0.08
0.06
0.04

0.04

0.02

0.02
0

0
0

5

10

15

20

0

5

time

10

15

20

time

Figure 4: Solution of the full (FOM) and reduced order model (ROM) corresponding to non-sampled parameter values (see values in Table 2) for for Input 1 with entires u1 (t) = 0.5, u2 (t) = 0.25, u3 (t) = 1,
u4 (t) = −1.
worst performance of ROM out of 10000 samples

0.1

0.03

0
0.02

-0.5
-1
1

0.01

0.05

temperature

0.5

p3

FOM
ROM

0.04

1

0

-0.05

-0.1

3
0

2
1

p2

-1

0

0

-0.15
0

p1

5

10

15

20

time

Figure 5: For Input 1: on the left, the relative L2 output error; on the right, solution of the full order model
and the reduced order model corresponding to the highest relative L2 error.

4

Conclusions and Future Work

In this paper, we presented conditions that ensure Hermite interpolation conditions for parametric bilinear
systems. These conditions also ensure that parametric directional derivatives of the reduced-order transfer
functions match the full-order transfer function at given interpolation points and directions. We demonstrate
18

worst performance of ROM out of 10000 samples

0.08

0.035

1

0.025

0

0.02

-0.5

0.015

-1
1

0.06

temperature

0.5

p3

FOM
ROM

0.03

0.04

0.02

0.01

3
0

2
1

p2

-1

0

0

0.005
0

-0.02
0

p1

5

10

15

20

time

Figure 6: For Input 2: on the left, the relative L2 output error; on the right, solution of the full order model
and the reduced order model corresponding to the highest relative L2 error.
the quality of our model reduction using two examples, one a well-known benchmark and the other an
interesting advection-diffusion equation. The performance was very good, and we emphasize that no effort
was made to select the sample points in parameter space. In fact, this approach is agnostic to the parameter
choices and can be easily embedded in well-known parameter selection schemes. The next natural steps are
to test this algorithm with different schemes and more challenging problems. We also intend to extend this
approach to parametric quadratic nonlinear systems.

References
[1] Mian Ilyas Ahmad, Peter Benner, and Pawan Goyal. Krylov subspace-based model reduction for a class
of bilinear descriptor systems. Journal of Computational and Applied Mathematics, 315:303–318, 2017.
[2] SAMIR A Al-Baiyat and MAAMAR Bettayeb. A new model reduction scheme for k-power bilinear
systems. In Decision and Control, 1993., Proceedings of the 32nd IEEE Conference on, pages 22–27.
IEEE, 1993.
[3] D. Amsallem and C. Farhat. An online method for interpolating linear parametric reduced-order models.
SIAM Journal on Scientific Computing, 33(5):2169–2198, 2011.
[4] A. C. Antoulas, D. C. Sorensen, and S. Gugercin. A survey of model reduction methods for large-scale
systems. Contemporary Mathematics, 280:193–219, 2001.
[5] A.C. Antoulas. Approximation of Large-Scale Dynamical Systems. SIAM Publications, Philadelphia,
PA, 2005.
[6] A.C. Antoulas, C.A. Beattie, and S. Gugercin. Interpolatory model reduction of large-scale dynamical
systems. In Efficient Modeling and Control of Large-Scale Systems, pages 3–58. Springer, 2010.
[7] Athanasios C Antoulas, Ion Victor Gosea, and Antonio Cosmin Ionita. Model reduction of bilinear
systems in the Loewner framework. SIAM Journal on Scientific Computing, 38(5):B889–B916, 2016.
[8] Z. Bai and D. Skoogh. A projection method for model reduction of bilinear dynamical systems. Linear
Algebra and its Applications, 415:406–425, 2006.
[9] U. Baur, C. Beattie, P. Benner, and S. Gugercin. Interpolatory projection methods for parameterized
model reduction. SIAM Journal on Scientific Computing, 33(5):2489–2518, 2011.
[10] U. Baur, P. Benner, and L. Feng. Model order reduction for linear and nonlinear systems: A systemtheoretic perspective. Arch. Comput. Methods Eng., 21(4):331–358, 2014.
19

[11] C.A. Beattie and S. Gugercin. Model reduction by rational interpolation. In P. Benner, A. Cohen,
M. Ohlberger, and K. Willcox, editors, Model Reduction and Approximation: Theory and Algorithms.
Available as http://arxiv.org/abs/1409.2140. SIAM, Philadelphia, PA, USA, 2017.
[12] P. Benner and T. Breiten. Interpolation-based H2 -model reduction of bilinear control systems. SIAM
Journal on Matrix Analysis & Applications, 33(3):859–885, 2012.
[13] P. Benner, A. Cohen, M. Ohlberger, and K. Willcox. Model Reduction and Approximation: Theory and
Algorithms. SIAM, 2017.
[14] P. Benner and L. Feng. A robust algorithm for parametric model order reduction based on implicit
moment matching. In A. Quarteroni and G. Rozza, editors, Reduced Order Methods for Modeling and
Computational Reduction, volume 9 of MS&A – Modeling, Simulation and Applications, pages 159–185.
Springer, 2014.
[15] P. Benner, S. Gugercin, and K. Willcox. A survey of projection-based model reduction methods for
parametric dynamical systems. SIAM Review, 57(4):483–531, 2015.
[16] Peter Benner, Tobias Breiten, and Tobias Damm. Generalised tangential interpolation for model reduction of discrete-time MIMO bilinear systems. International Journal of Control, 84(8):1398–1407,
2011.
[17] Peter Benner and Tobias Damm. Lyapunov equations, energy functionals, and model order reduction
of bilinear and stochastic systems. SIAM journal on control and optimization, 49(2):686–711, 2011.
[18] Peter Benner, Tobias Damm, and Yolanda R. Rodriguez Cruz. Dual pairs of generalized Lyapunov
inequalities and balanced truncation of stochastic linear systems. IEEE Transactions on Automatic
Control, 62(2):782–791, 2017.
[19] G. Berkooz, P. Holmes, and J.L. Lumley. The proper orthogonal decomposition in the analysis of
turbulent flows. Annual Review of Fluid Mechanics, 25:539–575, 1993.
[20] T. Breiten and T. Damm. Krylov subspace methods for model order reduction of bilinear control
systems. Systems & Control Letters, 59:443–450, 2010.
[21] T. Bui-Thanh, K. Willcox, and O. Ghattas. Model reduction for large-scale systems with highdimensional parametric input space. SIAM Journal on Scientific Computing, 30(6):3270–3288, 2008.
[22] S. Chaturantabut and D. Sorensen. Nonlinear model reduction via discrete empirical interpolation.
SIAM Journal on Scientific Computing, 32(5):2737–2764, 2010.
[23] L. Daniel, O.C. Siong, S.C. Low, K.H. Lee, and J. White. A multiparameter moment matching model
reduction approach for generating geometrically parameterized interconnect performance models. IEEE
Transactions on Computer-Aided Design of Integrated Circuits and Systems, 23(5):678–693, 2004.
[24] J. Degroote, J. Vierendeels, and K. Willcox. Interpolation among reduced-order matrices to obtain
parameterized models for design, optimization and probabilistic analysis. International Journal for
Numerical Methods in Fluids, 63(2):207–230, 2010.
[25] Z. Drmac and S. Gugercin. A new selection operator for the discrete empirical interpolation method
– improved a priori error bound and extensions. SIAM Journal on Scientific Computing, 38(2):A631–
A648, 2016.
[26] G. Flagg and S. Gugercin. Multipoint Volterra series interpolation and H2 optimal model reduction of
bilinear systems. SIAM Journal on Matrix Analysis and Applications, 2015.
[27] W Steven Gray and Joseph Mesko. Energy functions and algebraic gramians for bilinear systems. IFAC
Proceedings Volumes, 31(17):101–106, 1998.

20

[28] M.A. Grepl, Y. Maday, N.C. Nguyen, and A.T. Patera. Efficient reduced-basis treatment of nonaffine
and nonlinear partial differential equations. Mathematical Modelling and Numerical Analysis (M2AN),
41(3):575–605, 2007.
[29] S. Gugercin, A. C. Antoulas, and C. Beattie. H2 model reduction for large-scale linear dynamical
systems. SIAM Journal on Matrix Analysis and Applications, 30(2):609–638, 2008.
[30] P.K. Gunupudi, R. Khazaka, M.S. Nakhla, T. Smy, and D. Celo. Passive parameterized time-domain
macromodels for high-speed transmission-line networks. IEEE Transactions on Microwave Theory and
Techniques, 51(12):2347–2354, 2003.
[31] Carsten Hartmann, Boris Schäfer-Bung, and Anastasia Thöns-Zueva. Balanced averaging of bilinear
systems with applications to stochastic control. SIAM Journal on Control and Optimization, 51(3):2356–
2378, 2013.
[32] A. Hay, J. Borggaard, and D. Pelletier. Local improvements to reduced-order models using sensitivity
analysis of the proper orthogonal decomposition. Journal of Fluid Mechanics, 629:41–72, 2009.
[33] J. S Hesthaven, G. Rozza, and B. Stamm. Certified reduced basis methods for parametrized partial
differential equations. Springer, 2016.
[34] J.L. Lumley. The Structures of Inhomogeneous Turbulent Flow. Atmospheric Turbulence and Radio
Wave Propagation, pages 166–178, 1967.
[35] AJ Mayo and AC Antoulas. A framework for the solution of the generalized realization problem. Linear
algebra and its applications, 425(2):634–662, 2007.
[36] R. R. Mohler. Natural bilinear control processes. IEEE Transactions on Systems Science and Cybernetics, 6(3):192–197, 1970.
[37] R. R. Mohler. Nonlinear systems (vol. 2): applications to bilinear control. Prentice-Hall, Inc. Upper
Saddle River, NJ, USA, 1991.
[38] B. Moore. Principal component analysis in linear systems: Controllability, observability, and model
reduction. IEEE Transactions on Automatic Control, 26(1):17–32, 1981.
[39] C. Mullis and R. Roberts. Synthesis of minimum roundoff noise fixed point digital filters. IEEE
Transactions on Circuits and Systems, 23(9):551–562, 1976.
[40] H. Panzer, J. Mohring, R. Eid, and B. Lohmann. Parametric model order reduction by matrix interpolation. at–Automatisierungstechnik, 58(8):475–484, 2010.
[41] J. R. Phillips. Projection-based approaches for model reduction of weakly nonlinear, time-varying
systems. IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 22:171–
187, 2003.
[42] W.J. Rugh. Nonlinear system theory. Johns Hopkins University Press Baltimore, MD, 1981.
[43] T. Siu and M. Schetzen. Convergence of Volterra series representation and BIBO stability of bilinear
systems. International Journal of Systems Science, 22(12):2679–2684, 1991.
[44] D. D. Weiner and J. F. Spina. Sinusoidal Analysis and Modeling of Weakly Nonlinear Circuits: With
Application to Nonlinear Interference Effects. Van Nostrand Reinhold, 1980.
[45] R. Zimmermann. A locally parametrized reduced-order model for the linear frequency domain approach
to time-accurate computational fluid dynamics. SIAM Journal on Scientific Computing, 36(3):B508–
B537, 2014.
[46] Ralf Zimmermann. Local parametrization of subspaces on matrix manifolds via derivative information.
SIAM Journal on Scientific Computing, 36(3):B508–B537, 2014.

21

