opinf.operators

`opinf.operators`#

Operator classes for the individual terms of a model.

Nonparametric Operators

`OperatorTemplate`	Template for general operators \(\Ophat_{\ell}(\qhat,\u).\)
`InputMixin`	Mixin for operators whose `apply()` method acts on the input \(\u\).
`OpInfOperator`	Template for nonparametric operators that can be calibrated through Operator Inference, i.e., \(\Ophat_{\ell}(\qhat, \u) = \Ohat_{\ell}\d(\qhat, \u)\).
`ConstantOperator`	Constant operator \(\Ophat_{\ell}(\qhat,\u) = \chat \in \RR^{r}\).
`LinearOperator`	Linear state operator \(\Ophat_{\ell}(\qhat,\u) = \Ahat\qhat\) where \(\Ahat \in \RR^{r \times r}\).
`QuadraticOperator`	Quadratic state operator \(\Ophat_{\ell}(\qhat,\u) = \Hhat[\qhat\otimes\qhat]\) where \(\Hhat\in\RR^{r \times r^{2}}\).
`CubicOperator`	Cubic state operator \(\Ophat_{\ell}(\qhat,\u) = \Ghat[\qhat\otimes\qhat\otimes\qhat]\) where \(\Ghat\in\RR^{r \times r^{3}}\).
`InputOperator`	Linear input operator \(\Ophat_{\ell}(\qhat,\u) = \Bhat\u\) where \(\Bhat \in \RR^{r \times m}\).
`StateInputOperator`	Linear state / input interaction operator \(\Ophat_{\ell}(\qhat,\u) = \Nhat[\u\otimes\qhat]\) where \(\Nhat \in \RR^{r \times rm}\).

Parametric Operators

`ParametricOperatorTemplate`	Template for operators that depend on external parameters, \(\Ophat_{\ell}(\qhat,\u;\bfmu).\)
`ParametricOpInfOperator`	Template for operators that depend on external parameters, and which can be calibrated through operator inference, i.e., \(\Ophat_\ell(\qhat,\u;\bfmu) = \Ohat_\ell(\bfmu)\d_\ell(\qhat,\u)\).
`AffineConstantOperator`	Affine-parametric constant operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \chat_{\ell}(\bfmu) = \sum_{a=0}^{A_{\ell}-1}\theta_\ell^{(a)}\!(\bfmu)\,\chat_{\ell}^{(a)}.\)
`AffineLinearOperator`	Affine-parametric linear operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Ahat_{\ell}(\bfmu)\qhat = \left( \sum_{a=0}^{A_{\ell}-1}\theta_{\ell}^{(a)}\!(\bfmu)\,\Ahat_{\ell}^{(a)} \right)\qhat.\)
`AffineQuadraticOperator`	Affine-parametric quadratic operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Hhat_{\ell}(\bfmu)[\qhat\otimes\qhat] = \left( \sum_{a=0}^{A_{\ell}-1}\theta_{\ell}^{(a)}\!(\bfmu)\,\Hhat_{\ell}^{(a)} \right)[\qhat\otimes\qhat].\)
`AffineCubicOperator`	Affine-parametric cubic operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Ghat_{\ell}(\bfmu)[\qhat\otimes\qhat\otimes\qhat] = \left( \sum_{a=0}^{A_{\ell}-1}\theta_{\ell}^{(a)}\!(\bfmu)\,\Ghat_{\ell}^{(a)} \right)[\qhat\otimes\qhat\otimes\qhat].\)
`AffineInputOperator`	Affine-parametric input operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Bhat_{\ell}(\bfmu)\u = \left( \sum_{a=0}^{A_{\ell}-1}\theta_{\ell}^{(a)}\!(\bfmu)\,\Bhat_{\ell}^{(a)} \right)\u.\)
`AffineStateInputOperator`	Affine-parametric state-input operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Nhat_{\ell}(\bfmu)\qhat = \left( \sum_{a=0}^{A_{\ell}-1}\theta_{\ell}^{(a)}\!(\bfmu)\,\Nhat_{\ell}^{(a)} \right)[\u\otimes\qhat].\)
`InterpConstantOperator`	Parametric constant operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \chat(\bfmu) \in \RR^r\) where the parametric dependence is handled with elementwise interpolation.
`InterpLinearOperator`	Parametric linear operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Ahat(\bfmu)\qhat\) where \(\Ahat(\bfmu) \in \RR^{r \times r}\) and the parametric dependence is handled with elementwise interpolation.
`InterpQuadraticOperator`	Parametric quadratic operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Hhat(\bfmu)[\qhat\otimes\qhat]\) where \(\Ahat(\bfmu) \in \RR^{r \times r^2}\) and the parametric dependence is handled with elementwise interpolation.
`InterpCubicOperator`	Parametric cubic operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Ghat(\bfmu)[\qhat\otimes\qhat\otimes\qhat]\) where \(\Ghat(\bfmu) \in \RR^{r \times r^3}\) and the parametric dependence is handled with elementwise interpolation.
`InterpInputOperator`	Parametric input operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Bhat(\bfmu)\u\) where \(\Bhat(\bfmu) \in \RR^{r \times m}\) and the parametric dependence is handled with elementwise interpolation.
`InterpStateInputOperator`	Parametric state-input operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Nhat(\bfmu)[\u\otimes\qhat]\) where \(\Nhat(\bfmu) \in \RR^{r \times rm}\) and the parametric dependence is handled with elementwise interpolation.

Overview

opinf.operators classes represent the individual terms in a model. Nonparametric operators are functions of the state and input, while parametric operators are also dependent on one or more external parameters.
Operators that can be written as the product of a matrix and a known vector-valued function can be calibrated through solving a regression problem.
opinf.models objects are constructed with a list of operator objects. A model’s fit() method constructs and solves a regression problem to learn the operator matrices.

import numpy as np
import scipy.linalg as la

import opinf

Operators#

Models based on Operator Inference are systems of ordinary differential equations (or discrete-time difference equations) that can be written as a sum of terms,

(11)#\[ \begin{aligned} \ddt\qhat(t) = \sum_{\ell=1}^{n_\textrm{terms}}\Ophat_{\ell}(\qhat(t),\u(t)), \end{aligned} \]

where each \(\Ophat_{\ell}:\RR^{r}\times\RR^{m}\to\RR^{r}\) is a vector-valued function of the reduced state \(\qhat\in\RR^{r}\) and the input \(\u\in\RR^{m}\). We call these functions operators on this page.

Operator Inference calibrates operators that can be written as the product of a matrix \(\Ohat_{\ell}\in\RR^{r \times d_\ell}\) and a known (possibly nonlinear) vector-valued function \(\d_{\ell}:\RR^{r}\times\RR^{m}\to\RR^{d_ell}\),

\[ \begin{aligned} \Ophat_{\ell}(\qhat,\u) = \Ohat_{\ell}\d_{\ell}(\qhat,\u). \end{aligned} \]

Operators with this structure are called OpInf operators on this page. The goal of Operator Inference is to learn the operator matrix \(\Ohat_\ell\) for each OpInf operator in the model.

This module defines classes representing various operators.

Example

To represent a linear time-invariant (LTI) system

(12)#\[ \begin{align} \ddt\qhat(t) = \Ahat\qhat(t) + \Bhat\u(t), \qquad \Ahat\in\RR^{r \times r}, ~ \Bhat\in\RR^{r \times m}, \end{align} \]

we use the following operator classes.

Class	Definition	Operator matrix	data vector
`LinearOperator`	\(\Ophat_{1}(\qhat,\u) = \Ahat\qhat\)	\(\Ohat_{1} = \Ahat \in \RR^{r\times r}\)	\(\d_{1}(\qhat,\u) = \qhat\in\RR^{r}\)
`InputOperator`	\(\Ophat_{2}(\qhat,\u) = \Bhat\u\)	\(\Ohat_{2} = \Bhat \in \RR^{r\times m}\)	\(\d_{2}(\qhat,\u) = \u\in\RR^{m}\)

An opinf.models.ContinuousModel object can be instantiated with a list of operators objects to represent (12):

LTI_model = opinf.models.ContinuousModel(
    operators=[
        opinf.operators.LinearOperator(),
        opinf.operators.InputOperator(),
    ]
)

The operator matrices \(\Ahat\) and \(\Bhat\) are calibrated by calling LTI_model.fit().

Nonparametric Operators#

A nonparametric operator is a function of the state and input only. For OpInf operators, this means the operator matrix \(\Ohat_\ell\) is constant. See Parametric Operators for operators that also depend on one or more external parameters.

Available nonparametric OpInf operators are listed below.

`ConstantOperator`	Constant operator \(\Ophat_{\ell}(\qhat,\u) = \chat \in \RR^{r}\).
`LinearOperator`	Linear state operator \(\Ophat_{\ell}(\qhat,\u) = \Ahat\qhat\) where \(\Ahat \in \RR^{r \times r}\).
`QuadraticOperator`	Quadratic state operator \(\Ophat_{\ell}(\qhat,\u) = \Hhat[\qhat\otimes\qhat]\) where \(\Hhat\in\RR^{r \times r^{2}}\).
`CubicOperator`	Cubic state operator \(\Ophat_{\ell}(\qhat,\u) = \Ghat[\qhat\otimes\qhat\otimes\qhat]\) where \(\Ghat\in\RR^{r \times r^{3}}\).
`InputOperator`	Linear input operator \(\Ophat_{\ell}(\qhat,\u) = \Bhat\u\) where \(\Bhat \in \RR^{r \times m}\).
`StateInputOperator`	Linear state / input interaction operator \(\Ophat_{\ell}(\qhat,\u) = \Nhat[\u\otimes\qhat]\) where \(\Nhat \in \RR^{r \times rm}\).

Nonparametric OpInf operators can be instantiated without arguments. If the operator matrix is known, it can be passed into the constructor or set later with set_entries(). The operator matrix is stored as the entries attribute and can be accessed with slicing operations [:].

N = opinf.operators.StateInputOperator()
print(N)

StateInputOperator
  state_dimension: None
  input_dimension: None
  entries.shape:   None

r = 4
m = 2
Ohat = np.arange(r * r * m).reshape((r, r * m))

N.set_entries(Ohat)
print(N)

StateInputOperator
  state_dimension: 4
  input_dimension: 2
  entries.shape:   (4, 8)

N.entries  # or N[:]

array([[ 0,  1,  2,  3,  4,  5,  6,  7],
       [ 8,  9, 10, 11, 12, 13, 14, 15],
       [16, 17, 18, 19, 20, 21, 22, 23],
       [24, 25, 26, 27, 28, 29, 30, 31]])

There are two ways to determine operator matrices in the context of model reduction:

Non-intrusive Operator Inference: Learn operator matrices from data.
Intrusive (Petrov–)Galerkin Projection: Compress an existing high-dimensional operator.

Once the entries are set, the following methods are used to compute the action of the operator or its derivatives.

apply(): compute the operator action \(\Ophat_\ell(\qhat, \u)\).
jacobian(): construct the state Jacobian \(\ddqhat\Ophat_\ell(\qhat, \u)\).

qhat = np.zeros(r)
u = np.ones(m)

N.apply(qhat, u)

array([0., 0., 0., 0.])

N.jacobian(qhat, u)

array([[ 4.,  6.,  8., 10.],
       [20., 22., 24., 26.],
       [36., 38., 40., 42.],
       [52., 54., 56., 58.]])

Learning Operators from Data#

Attention

This section describes the crux of Operator Inference.

Operator Inference requires state, input, and “left-hand side” data \(\{(\qhat_j,\u_j,\z_j)\}_{j=0}^{k-1}\) that approximately satisfy the desired model dynamics. For time-continuous models, \(\z_j\) is the time derivative of the state data; for fully discrete models, \(\z_j\) is the ``next state,’’ usually \(\qhat_{j+1}\). For (11), and assuming each operator is an OpInf operator, the data should approximately satisfy

(13)#\[ \begin{aligned} \z_j \approx \Ophat(\qhat_j, \u_j) = \sum_{\ell=1}^{n_\textrm{terms}} \Ophat_{\ell}(\qhat_j, \u_j) = \sum_{\ell=1}^{n_\textrm{terms}} \Ohat_{\ell}\d_{\ell}(\qhat_j, \u_j), \qquad j = 0, \ldots, k-1. \end{aligned} \]

Operator Inference determines the operator matrices \(\Ohat_1,\ldots,\Ohat_{n_\textrm{terms}}\) through a regression problem, written generally as

\[ \begin{aligned} \text{find}\quad\Ohat\quad\text{such that}\quad \Z \approx \Ohat\D\trp \quad\Longleftrightarrow\quad \D\Ohat\trp \approx \Z\trp, \end{aligned} \]

where \(\D\) is the \(k \times d\) data matrix formed from state and input snapshots and \(\Z\) is the \(r \times k\) matrix of left-hand side data. To arrive at this problem, we write

\[\begin{split} \begin{aligned} \z_j \approx \sum_{\ell=1}^{n_\textrm{terms}} \Ohat_{\ell}\d_{\ell}(\qhat_j, \u_j) = [~\Ohat_{1}~~\cdots~~\Ohat_{n_\textrm{terms}}~] \left[\begin{array}{c} \d_{1}(\qhat_j, \u_j) \\ \vdots \\ \d_{n_\textrm{terms}}(\qhat_j, \u_j) \end{array}\right]. \end{aligned} \end{split}\]

Collecting this matrix-vector product for each \(j=0,\ldots,k-1\) results in the matrix-matrix product system

\[\begin{split} \begin{aligned} \left[\begin{array}{c|c|c} & & \\ \z_0 & \cdots & \z_{k-1} \\ & & \end{array}\right] \approx [~\Ohat_{1}~~\cdots~~\Ohat_{n_\textrm{terms}}~] \left[\begin{array}{ccc} \d_{1}(\qhat_0, \u_0) & \cdots & \d_{1}(\qhat_{k-1}, \u_{k-1}) \\ \vdots & & \vdots \\ \d_{n_\textrm{terms}}(\qhat_0, \u_0) & \cdots & \d_{n_\textrm{terms}}(\qhat_{k-1}, \u_{k-1}) \end{array}\right], \end{aligned} \end{split}\]

which is \(\Z \approx \Ohat\D\trp\) with

\[\begin{split} \begin{aligned} \Z &= [~\z_0~~\cdots~~\z_{k-1}~] \in \RR^{r\times k}, \\ & \\ \Ohat &= [~\Ohat_{1}~~\cdots~~\Ohat_{n_\textrm{terms}}~] \in \RR^{r \times d}, \\ & \\ \D\trp &= \left[\begin{array}{ccc} \d_{1}(\qhat_0, \u_0) & \cdots & \d_{1}(\qhat_{k-1}, \u_{k-1}) \\ \vdots & & \vdots \\ \d_{n_\textrm{terms}}(\qhat_0, \u_0) & \cdots & \d_{n_\textrm{terms}}(\qhat_{k-1}, \u_{k-1}) \end{array}\right] \in \RR^{d \times k}, \end{aligned} \end{split}\]

where \(d = \sum_{i=1}^{n_\textrm{terms}}d_\ell\).

Nonparametric OpInf operator classes have two static methods that facilitate constructing the operator regression problem.

operator_dimension(): given the state dimension \(r\) and the input dimension \(r\), return the data dimension \(d_\ell\).
datablock(): given the state-input data pairs \(\{(\qhat_j,\u_j)\}_{j=0}^{k-1}\), form the matrix

\[\begin{split} \begin{aligned} \D_{\ell}\trp = \left[\begin{array}{c|c|c|c} & & & \\ \d_{\ell}(\qhat_0,\u_0) & \d_{\ell}(\qhat_1,\u_1) & \cdots & \d_{\ell}(\qhat_{k-1},\u_{k-1}) \\ & & & \end{array}\right] \in \RR^{d_{\ell} \times k}. \end{aligned} \end{split}\]

The complete data matrix \(\D\) is the concatenation of the data matrices from each operator:

\[\begin{split} \begin{aligned} \D = \left[\begin{array}{ccc} & & \\ \D_1 & \cdots & \D_{n_\textrm{terms}} \\ & & \end{array}\right]. \end{aligned} \end{split}\]

Model Classes Do the Work

Model classes from opinf.models are instantiated with a list of operators. The model’s fit() method calls the datablock() method of each OpInf operator to assemble the full data matrix \(\D\), solves the regression problem for the full data matrix \(\Ohat\) (see opinf.lstsq), and extracts from \(\Ohat\) the individual operator matrix \(\Ohat_{\ell}\) for each \(\ell = 1, \ldots, n_{\textrm{terms}}\) using the operator_dimension().

Example

For the LTI system (12), the Operator Inference regression \(\Z \approx \Ohat\D\trp\) has the following matrices.

\[\begin{split} \begin{aligned} \Z &= [~\z_0~~\cdots~~\z_{k-1}~] \in \RR^{r\times k}, \\ \\ \Ohat &= [~\Ahat~~\Bhat~] \in \RR^{r \times (r + m)}, \\ \\ \D\trp &= \left[\begin{array}{ccc} \qhat_0 & \cdots & \qhat_{k-1} \\ \u_0 & \cdots & \u_{k-1} \end{array}\right] \in \RR^{(r + m) \times k}. \end{aligned} \end{split}\]

In this setting, \(\z_j = \dot{\qhat}_j\), the time derivative of \(\qhat_j\). Collecting the state snapshots in the matrix \(\Qhat = [~\qhat_0~~\cdots~~\qhat_{k-1}~]\in\RR^{r\times k}\) and the inputs in the matrix \(\U = [~\u_0~~\cdots~~\u_{k-1}~]\), the full data matrix can be abbreviated as \(\D = [~\Qhat\trp~~\U\trp~]\).

If the regression \(\Z \approx \Ohat\D\trp\) is treated as an ordinary least-squares problem, the optimization problem to solve is given by

\[ \begin{aligned} \min_{\Ohat}\left\| \D\Ohat\trp - \Z\trp \right\|_F^2 = \min_{\Ahat,\Bhat}\sum_{j=0}^{k-1}\left\| \Ahat\qhat_j + \Bhat\u_j - \z_j \right\|_2^2. \end{aligned} \]

That is, the ordinary least-squares Operator Inference regression minimizes a sum of residuals of the model equation (12) with respect to available data.

Operators with Entries are Not Recalibrated

Only operators whose entries are not initialized (set to None) when a model is constructed are learned with Operator Inference when fit() is called. For example, suppose for the LTI system (12) an appropriate input matrix \(\Bhat\) is known and stored as the variable B_.

LTI_model = opinf.models.ContinuousModel(
    operators=[
        opinf.operators.LinearOperator(),   # No entries specified.
        opinf.operators.InputOperator(B_),  # Entries set to B_.
    ]
)

In this case, LIT_model.fit() only determines the entries of the LinearOperator object using Operator Inference. The known information for \(\Bhat\) is absorbed into the matrix \(\Z\), so the Operator Inference regression \(\Z \approx \Ohat\D\trp\) is now defined with the matrices

\[\begin{split} \begin{aligned} \Z &= [~(\z_0 - \Bhat\u_0)~~\cdots~~(\z_{k-1} - \Bhat\u_{k-1})~] \in \RR^{r\times k}, \\ \Ohat &= \Ahat \in \RR^{r \times r}, \qquad \D\trp = \Qhat \in \RR^{r \times k}. \end{aligned} \end{split}\]

Using ordinary least-squares regression, the optimization problem is given by

\[ \begin{aligned} &\min_{\Ahat}\sum_{j=0}^{k-1}\left\| \Ahat\qhat_j - (\z_j - \Bhat\u_j) \right\|_2^2. \end{aligned} \]

A more likely scenario is that \(\Bhat\in\RR^{r \times m}\) is derived from a known \(\B\in\RR^{n \times m}\), which is the subject of the next section.

Learning Operators via Projection#

The goal of Operator Inference is to learn operator entries from data because the operators of a high-fidelity / full-order model are unknown or computationally inaccessible. However, in some scenarios a subset of the full-order model operators are known, in which case the corresponding reduced-order model operators can be determined through intrusive projection. Consider a full-order operator \(\Op:\RR^{n}\times\RR^{m}\to\RR^{n}\), written \(\Op_{\ell}(\q,\u)\), where

\(\q\in\RR^n\) is the full-order state, and
\(\u\in\RR^m\) is the input.

Given a trial basis \(\Vr\in\RR^{n\times r}\) and a test basis \(\Wr\in\RR^{n\times r}\), the corresponding intrusive projection of \(\Op_{\ell}\) is the operator \(\Ophat_{\ell}:\RR^{r}\times\RR^{m}\to\RR^{r}\) defined by

\[ \begin{aligned} \Ophat_{\ell}(\qhat, \u) = (\Wr\trp\Vr)^{-1}\Wr\trp\Op_{\ell}(\Vr\qhat, \u) \end{aligned} \]

where

\(\qhat\in\RR^{r}\) is the reduced-order state, and
\(\u\in\RR^{m}\) is the input (the same as before).

This approach uses the low-dimensional state approximation \(\q \approx \Vr\qhat\). If \(\Wr = \Vr\), the result is called a Galerkin projection. Note that if \(\Vr\) has orthonormal columns, we have in this case the simplification

\[ \begin{aligned} \Ophat_{\ell}(\qhat, \u) = \Vr\trp\Op_{\ell}(\Vr\qhat, \u). \end{aligned} \]

If \(\Wr \neq \Vr\), the result is called a Petrov–Galerkin projection.

Example

Consider the bilinear operator \(\Op(\q,\u) = \N[\u\otimes\q]\) where \(\N\in\RR^{n \times nm}\). This type of operator can represented as a StateInputOperator. The intrusive Petrov–Galerkin projection of \(\Op\) is the bilinear operator

\[ \begin{aligned} \Ophat(\qhat,\u) = (\Wr\trp\Vr)^{-1}\Wr\trp\N[\u\otimes\Vr\qhat] = \Nhat[\u\otimes\qhat] \end{aligned} \]

where \(\Nhat = (\Wr\trp\Vr)^{-1}\Wr\trp\N(\I_m\otimes\Vr) \in \RR^{r\times rm}\). Hence, \(\Ophat\) can also be represented as a StateInputOperator. Using Galerkin projection (\(\Wr = \Vr\)), \(\Nhat\) simplifies to \(\Nhat = \Vr\trp\N(\I_m\otimes\Vr)\).

Every operator class has a galerkin() method that performs intrusive projection.

Custom Nonparametric Operators#

New nonparametric operators can be defined by inheriting from OperatorTemplate or, for operators that can be calibrated through Operator Inference, OpInfOperator. Once implemented, the verify() method may be used to test for consistency between apply() and the other methods outlined below.

For an arbitrary operator (not calibrated through Operator Inference), use the following inheritance template.

See OperatorTemplate for details on the arguments for each method.

For a new Operator Inference operator, use the following inheritance template.

See OpInfOperator for details on the arguments for each method.

Developer Notes

If the operator depends on the input \(\u\), the class should also inherit from InputMixin and set the input_dimension attribute.
The jacobian() method is optional, but opinf.models objects have a jacobian() method that calls jacobian() for each of its operators. In an opinf.models.ContinuousModel, the Jacobian is required for some time integration strategies used in predict().
The galerkin() method is optional, but opinf.models objects have a galerkin() method that calls galerkin() for each of its operators.
The save() and load() methods should be implemented using opinf.utils.hdf5_savehandle() and opinf.utils.hdf5_loadhandle(), respectively.

Example: Hadamard Product with a Fixed Vector#

Consider the operator \(\Ophat_{\ell}(\qhat, \u) = \qhat \ast \hat{\s}\), where \(\hat{\s}\in\RR^{r}\) is a constant vector and \(\ast\) denotes the Hadamard (elementwise) product (* in NumPy). To implement an operator for this class, we first calculate its state Jacobian and determine the operator produced by (Petrov–)Galerkin projection.

Let \(\qhat = [~\hat{q}_1~~\cdots~~\hat{q}_r~]\trp\) and \(\hat{\s} = [~\hat{s}_1~~\cdots~~\hat{s}_r~]\trp\), i.e., the \(i\)-th entry of \(\Ophat_{\ell}(\qhat, \u)\) is \(\hat{q}_i\hat{s}_i\). Then the \((i,j)\)-th entry of the Jacobian is

\[\begin{split} \begin{aligned} \frac{\partial}{\partial \hat{q}_j}\left[\hat{q}_i\hat{s}_i\right] = \begin{cases} \hat{s}_i & \textrm{if}~i = j, \\ 0 & \textrm{else}. \end{cases} \end{aligned} \end{split}\]

That is, \(\ddqhat\Ophat_{\ell}(\qhat, \u) = \operatorname{diag}(\hat{\s}).\)

Now consider a version of this operator with a large state dimension, \(\Op_{\ell}(\q, \u) = \q \ast \s\) for \(\q,\s\in\mathbb{R}^{n}\). For basis matrices \(\Vr,\Wr\in\mathbb{R}^{n \times r}\), the Petrov–Galerkin projection of \(\Op_{\ell}\) is given by

\[ \begin{aligned} \Ophat_{\ell}(\qhat, \u) = (\Wr\trp\Vr)^{-1}\Wr\trp\Op_{\ell}(\Vr\qhat, \u) = (\Wr\trp\Vr)^{-1}\Wr\trp((\Vr\qhat)\ast\s). \end{aligned} \]

It turns out that this product can be written as a matrix-vector product \(\Ahat\qhat\) where \(\Ahat\) depends on \(\Vr\) and \(\s\). Therefore, galerkin() should return a LinearOperator with entries matrix \(\Ahat\).

The following class inherits from OperatorTemplate, stores \(\hat{\s}\) and sets the state dimension \(r\) in the constructor, and implements the methods outlined the inheritance template.

class HadamardOperator(opinf.operators.OperatorTemplate):
    """Custom non-OpInf nonparametric operator for the Hadamard product."""

    # Constructor -------------------------------------------------------------
    def __init__(self, s):
        """Construct the operator and set the state_dimension."""
        self.svector = np.array(s)
        self._jac = np.diag(self.svector)

    # Required properties and methods -----------------------------------------
    @property
    def state_dimension(self):
        """Dimension of the state that the operator acts on."""
        return self.svector.shape[0]

    def apply(self, state, input_=None):
        """Apply the operator to the given state / input."""
        svec = self.svector
        if state.ndim == 2:
            svec = svec.reshape((self.state_dimension, 1))
        return state * svec

    # Optional methods --------------------------------------------------------
    def jacobian(self, state, input_=None):
        """Construct the state Jacobian of the operator."""
        return self._jac

    def galerkin(self, Vr, Wr=None):
        """Get the (Petrov-)Galerkin projection of this operator."""
        if Wr is None:
            Wr = Vr
        n = self.state_dimension
        r = Vr.shape[1]

        M = la.khatri_rao(Vr.T, np.eye(n)).T
        Ahat = Wr.T @ (M.reshape((n, r, n)) @ self.svector)
        if not np.allclose((WrTVr := Wr.T @ Vr), np.eye(r)):
            Ahat = la.solve(WrTVr, Ahat)
        return opinf.operators.LinearOperator(Ahat)

    def save(self, savefile, overwrite=False):
        """Save the operator to an HDF5 file."""
        with opinf.utils.hdf5_savehandle(savefile, overwrite) as hf:
            hf.create_dataset("svector", data=self.svector)

    @classmethod
    def load(cls, loadfile):
        """Load an operator from an HDF5 file."""
        with opinf.utils.hdf5_loadhandle(loadfile) as hf:
            return cls(hf["svector"][:])

    def copy(self):
        """Make a copy of the operator."""
        return self.__class__(self.svector)

r = 10
svec = np.random.standard_normal(r)
hadamard = HadamardOperator(svec)

print(hadamard)

HadamardOperator
  state_dimension: 10

hadamard.verify()

apply() is consistent with state_dimension
jacobian() is consistent with state_dimension
jacobian() finite difference relative errors
  ------------------------------------------
  h = 1.00e-01	error = 5.7463e-16
  h = 1.00e-02	error = 1.3123e-14
  h = 1.00e-03	error = 7.3735e-14
  h = 1.00e-04	error = 7.7231e-13
  h = 1.00e-05	error = 3.0133e-12
  h = 1.00e-06	error = 9.6158e-11
  h = 1.00e-07	error = 7.1203e-10
  h = 1.00e-08	error = 5.0540e-09
  h = 1.00e-09	error = 9.4283e-08
  h = 1.00e-10	error = 3.7839e-07
galerkin() is consistent with apply()
copy() is consistent with state_dimension
copy() preserves the results of apply()
save()/load() is consistent with state_dimension
save()/load() preserves the results of apply()

Example: Weighted Hadamard Input Operator#

Consider an operator

\[ \begin{aligned} \Ophat_{\ell}(\qhat,\u) = \Ohat_{\ell}(\u\ast\u), \end{aligned} \]

where \(\Ohat_{\ell}\in\RR^{r \times m}\) and \(\ast\) is the Hadamard (elementwise) product. The matrix \(\Ohat_{\ell}\) can be calibrated with Operator Inference. Since \(\Ophat_{\ell}\) does not depend on the state \(\qhat\), the state Jacobian is zero.

class HadamardInputOperator(
    opinf.operators.OpInfOperator,
    opinf.operators.InputMixin,
):
    """Custom nonparametric OpInf operator: Hadamard of inputs."""

    # Required methods --------------------------------------------------------
    @property
    def input_dimension(self):
        """Dimension of the input that the operator acts on."""
        return self.entries.shape[1] if self.entries is not None else None

    @opinf.utils.requires("entries")
    def apply(self, state=None, input_=None):
        """Apply the operator to the given state / input."""
        return self.entries @ (input_**2)

    @staticmethod
    def datablock(states, inputs=None):
        """Return the data matrix block corresponding to the operator."""
        return inputs**2

    @staticmethod
    def operator_dimension(r=None, m=None):
        """Column dimension of the operator entries matrix."""
        return m

    # Optional methods --------------------------------------------------------
    def galerkin(self, Vr, Wr=None):
        """Get the (Petrov-)Galerkin projection of this operator."""
        if Wr is None:
            Wr = Vr
        r = Vr.shape[1]

        Ohat = Wr.T @ self.entries
        if not np.allclose((WrTVr := Wr.T @ Vr), np.eye(r)):
            Ohat = la.solve(WrTVr, Ohat)

        return self.__class__(Ohat)

# Instantiate an operator without entries.
hadamard_input = HadamardInputOperator()

print(hadamard_input)

HadamardInputOperator
  state_dimension: None
  input_dimension: None
  entries.shape:   None

hadamard_input.verify()

operator_dimension() is consistent with datablock()
cannot verify apply() when entries=None

m = 4
Ohat = np.random.random((r, m))
hadamard_input.set_entries(Ohat)

print(hadamard_input)

HadamardInputOperator
  state_dimension: 10
  input_dimension: 4
  entries.shape:   (10, 4)

hadamard_input.verify()

operator_dimension() is consistent with datablock()
apply() is consistent with state_dimension and input_dimension
jacobian() = 0
galerkin() is consistent with apply()
copy() is consistent with state_dimension and input_dimension
copy() preserves the results of apply()
save()/load() is consistent with state_dimension and input_dimension
save()/load() preserves the results of apply()

Parametric Operators#

An operator is called parametric if it depends on an independent parameter vector \(\bfmu\in\RR^{p}\), i.e., \(\Ophat_{\ell} = \Ophat_{\ell}(\qhat,\u;\bfmu)\) When the parameter vector is fixed, a parametric operator becomes nonparametric. In particular, a parametric operator’s evaluate() method accepts a parameter vector \(\bfmu\) and returns an instance of a nonparametric operator.

Parametric OpInf operators have the form \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Ohat_{\ell}(\bfmu)\d_{\ell}(\qhat,\u)\) defined by the matrix-valued function \(\Ohat_{\ell}:\RR^{p}\to\RR^{r\times d_\ell}\) and (as in the nonparametric case) the data vector \(\d_{\ell}:\RR^{r}\times\RR^{m}\to\RR^{d_\ell}\). This module provides two options for the parameterization of \(\Ohat_{\ell}(\bfmu)\): affine expansion and elementwise interpolation. In each case, Operator Inference begins with \(s\) training parameter values \(\bfmu_{0},\ldots,\bfmu_{s-1}\) and corresponding state, input, and left-hand side data \(\{(\qhat_{i,j},\u_{i,j},\z_{i,j})\}_{j=0}^{k_{i}-1}\) for each training parameter value \(\bfmu_{i}\). A regression of the form \(\Z \approx \Ohat\D\trp\) is formed as in the nonparametric case, with the structure of the matrices \(\Ohat\) and \(\D\) depending on the choice of parameterization for each \(\Ohat_{\ell}(\bfmu)\).

Example

Let \(\bfmu = [~\mu_{0}~~\mu_{1}~]\trp\). The operator

\[ \begin{aligned} \Ophat_1(\qhat,\u;\bfmu) = (\mu_{0}\Ahat_{0} + \mu_{1}^{2}\Ahat_{1})\qhat \end{aligned} \]

is a parametric OpInf operator because it can be written as \(\Ophat_1(\qhat,\u;\bfmu) = \Ohat_1(\bfmu)\d_1(\qhat,\u)\) with \(\Ohat_1(\bfmu) = \mu_{0}\Ahat_{0} + \mu_{1}^{2}\Ahat_{1}\) and \(\d_1(\qhat,\u) = \qhat\).

This operator can be represented with an AffineLinearOperator. For a given parameter value, the evaluate() method returns a LinearOperator instance.

Affine Operators#

Affine parametric OpInf operators \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Ohat_{\ell}(\bfmu)\d_{\ell}(\qhat,\u)\) parameterize the operator matrix \(\Ohat_{\ell}(\bfmu)\) as a sum of constant matrices with parameter-dependent scalar coefficients,

\[ \begin{aligned} \Ohat_{\ell}(\bfmu) &= \sum_{a=0}^{A_{\ell}-1}\theta_{\ell}^{(a)}\!(\bfmu)\,\Ohat_{\ell}^{(a)}, \end{aligned} \]

where each \(\theta_{\ell}^{(a)}:\RR^{p}\to\RR\) is a scalar-valued function and each \(\Ohat_{\ell}^{(a)}\in\RR^{r\times d_\ell}\) is constant. Affine expansions are grouped such that the coefficient functions \(\theta_{\ell}^{(0)},\ldots,\theta_{\ell}^{(A_{\ell}-1)}\) are linearly independent.

Affine parametric operators arise in Operator Inference settings because linear projection preserves affine structure.

Available affine parametric operators are listed below.

`AffineConstantOperator`	Affine-parametric constant operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \chat_{\ell}(\bfmu) = \sum_{a=0}^{A_{\ell}-1}\theta_\ell^{(a)}\!(\bfmu)\,\chat_{\ell}^{(a)}.\)
`AffineLinearOperator`	Affine-parametric linear operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Ahat_{\ell}(\bfmu)\qhat = \left( \sum_{a=0}^{A_{\ell}-1}\theta_{\ell}^{(a)}\!(\bfmu)\,\Ahat_{\ell}^{(a)} \right)\qhat.\)
`AffineQuadraticOperator`	Affine-parametric quadratic operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Hhat_{\ell}(\bfmu)[\qhat\otimes\qhat] = \left( \sum_{a=0}^{A_{\ell}-1}\theta_{\ell}^{(a)}\!(\bfmu)\,\Hhat_{\ell}^{(a)} \right)[\qhat\otimes\qhat].\)
`AffineCubicOperator`	Affine-parametric cubic operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Ghat_{\ell}(\bfmu)[\qhat\otimes\qhat\otimes\qhat] = \left( \sum_{a=0}^{A_{\ell}-1}\theta_{\ell}^{(a)}\!(\bfmu)\,\Ghat_{\ell}^{(a)} \right)[\qhat\otimes\qhat\otimes\qhat].\)
`AffineInputOperator`	Affine-parametric input operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Bhat_{\ell}(\bfmu)\u = \left( \sum_{a=0}^{A_{\ell}-1}\theta_{\ell}^{(a)}\!(\bfmu)\,\Bhat_{\ell}^{(a)} \right)\u.\)
`AffineStateInputOperator`	Affine-parametric state-input operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Nhat_{\ell}(\bfmu)\qhat = \left( \sum_{a=0}^{A_{\ell}-1}\theta_{\ell}^{(a)}\!(\bfmu)\,\Nhat_{\ell}^{(a)} \right)[\u\otimes\qhat].\)

Affine parametric operators are instantiated with a function \(\boldsymbol{\theta}_{\ell}(\mu) = [~\theta_{\ell}^{(0)}(\bfmu)~~\cdots~~\theta_{\ell}^{(A_{\ell}-1)}(\bfmu)~]\trp\) for the affine expansion coefficients, the number of terms \(A_{\ell}\) in the expansion, and with or without the operator matrices \(\Ohat_{\ell}^{(1)},\ldots,\Ohat_{\ell}^{(A_{\ell})}\).

thetas = lambda mu: np.array([mu[0], mu[1] ** 2])
A = opinf.operators.AffineLinearOperator(thetas, nterms=2)
print(A)

AffineLinearOperator
  state_dimension:     None
  parameter_dimension: None
  expansion terms:     2
  evaluate(parameter) -> LinearOperator

# Set the constant operator matrices in the affine expansion.
r = 5
Ahats = [np.ones((r, r)), np.eye(r)]
A.set_entries(Ahats, fromblock=False)

print(A)

AffineLinearOperator
  state_dimension:     5
  parameter_dimension: None
  expansion terms:     2
  evaluate(parameter) -> LinearOperator

A.entries

[array([[1., 1., 1., 1., 1.],
        [1., 1., 1., 1., 1.],
        [1., 1., 1., 1., 1.],
        [1., 1., 1., 1., 1.],
        [1., 1., 1., 1., 1.]]),
 array([[1., 0., 0., 0., 0.],
        [0., 1., 0., 0., 0.],
        [0., 0., 1., 0., 0.],
        [0., 0., 0., 1., 0.],
        [0., 0., 0., 0., 1.]])]

# Evaluate the parametric operator at a fixed parameter value,
# resulting in a nonparametric operator.
A_nonparametric = A.evaluate([2, 4])
print(A_nonparametric)

LinearOperator
  state_dimension: 5
  entries.shape:   (5, 5)

A_nonparametric.entries

array([[18.,  2.,  2.,  2.,  2.],
       [ 2., 18.,  2.,  2.,  2.],
       [ 2.,  2., 18.,  2.,  2.],
       [ 2.,  2.,  2., 18.,  2.],
       [ 2.,  2.,  2.,  2., 18.]])

# The parameter dimension p = 2 was recorded when A was evaluated.
print(A)

AffineLinearOperator
  state_dimension:     5
  parameter_dimension: 2
  expansion terms:     2
  evaluate(parameter) -> LinearOperator

Operator Inference for Affine Operators

Consider a model with a single affine operator,

\[ \begin{aligned} \z = \Ophat_{1}(\qhat,\u;\bfmu) &= \left(\sum_{a=0}^{A_{1}-1}\theta_{1}^{(a)}\!(\bfmu)\,\Ohat_{1}^{(a)}\right)\d_1(\qhat,\u). \end{aligned} \]

For the moment we will neglect the \(\ell = 1\) subscript. With data \(\{(\qhat_{i,j},\u_{i,j},\z_{i,j})\}_{j=0}^{k_{i}-1}\) corresponding to \(s\) training parameter values \(\bfmu_0,\ldots,\bfmu_{s-1}\), we seek the \(A\) operator matrices \(\Ohat^{(0)},\ldots,\Ohat^{(A-1)}\) such that for each parameter index \(i\) and time index \(j\), we have

\[\begin{split} \begin{aligned} \z_{i,j} \approx \Ophat(\qhat_{i,j},\u_{i,j};\bfmu_{i}) &= \left(\sum_{a=0}^{A-1}\theta^{(a)}\!(\bfmu_{i})\,\Ohat^{(a)}\right)\d(\qhat_{i,j},\u_{i,j}) \\ &= \left[\begin{array}{ccc} \Ohat^{(0)} & \cdots & \Ohat^{(A-1)} \end{array}\right] \underbrace{\left[\begin{array}{c} \theta^{(0)}\!(\bfmu_{i})\,\d(\qhat_{i,j},\u_{i,j}) \\ \vdots \\ \theta^{(A-1)}\!(\bfmu_{i})\,\d(\qhat_{i,j},\u_{i,j}) \end{array}\right]}_{\d_{i,j}\in\RR^{dA}}, \end{aligned} \end{split}\]

where \(d\) is the output dimension of \(\d\). Collecting these expressions for each time index \(j = 0, \ldots, k_i - 1\) (but keeping the parameter index \(i\) fixed for the moment) results in

\[ \begin{aligned} \underbrace{\left[\begin{array}{ccc} \z_{i,0} & \cdots & \z_{i,k_i-1} \end{array}\right]}_{\Z_i\in\RR^{r\times k_i}} \approx \left[\begin{array}{ccc} \Ohat^{(0)} & \cdots & \Ohat^{(A-1)} \end{array}\right] \underbrace{\left[\begin{array}{ccc} \d_{i,0} & \cdots & \d_{i,k_i-1} \end{array}\right]}_{\D_i\trp\in\RR^{dA \times k_i}}. \end{aligned} \]

Finally, we concatenate each of these expressions for each parameter index \(i = 0,\ldots, s-1\) to arrive at

\[ \begin{aligned} \underbrace{\left[\begin{array}{ccc} \Z_{0} & \cdots & \Z_{s-1} \end{array}\right]}_{\Z\in\RR^{r\times K}} \approx \left[\begin{array}{ccc} \Ohat^{(0)} & \cdots & \Ohat^{(A-1)} \end{array}\right] \underbrace{\left[\begin{array}{ccc} \D_{0}\trp & \cdots & \D_{s-1}\trp \end{array}\right]}_{\D\trp\in\RR^{dA \times K}}, \end{aligned} \]

where \(K = \sum_{i=0}^{s-1}k_i\), the total number of available data instances. This is the familiar \(\Z \approx \Ohat\D\trp\) where \(\Ohat = [~\Ohat^{(1)}~~\cdots~~\Ohat^{(A)}~]\), which can be solved for \(\Ohat\) using opinf.lstsq.

The construction of \(\D\) is taken care of through the datablock() method of the affine

For models with multiple affine operators, the operator matrix \(\Ohat\) is further concatenated horizontally to accommodate the operator matrices from each affine expansion, and the data matrix \(\D\trp\) gains additional block rows.

Interpolatory Operators#

Interpolatory parametric OpInf operators define the parametric dependence of the operator matrix on \(\bfmu\) through elementwise interpolation. That is,

\[ \begin{aligned} \Ophat_{\ell}(\qhat,\u;\bfmu) = \Ohat_{\ell}(\bfmu)\d_{\ell}(\qhat,\u), \end{aligned} \]

where \(\Ohat_{\ell}(\bfmu)\) is determined by interpolating \(s\) matrices \(\Ohat_{\ell}^{(0)},\ldots,\Ohat_{\ell}^{(s-1)}\). In the context of Operator Inference, \(s\) is the number of training parameter values.

Available interpolatory operators are listed below.

`InterpConstantOperator`	Parametric constant operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \chat(\bfmu) \in \RR^r\) where the parametric dependence is handled with elementwise interpolation.
`InterpLinearOperator`	Parametric linear operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Ahat(\bfmu)\qhat\) where \(\Ahat(\bfmu) \in \RR^{r \times r}\) and the parametric dependence is handled with elementwise interpolation.
`InterpQuadraticOperator`	Parametric quadratic operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Hhat(\bfmu)[\qhat\otimes\qhat]\) where \(\Ahat(\bfmu) \in \RR^{r \times r^2}\) and the parametric dependence is handled with elementwise interpolation.
`InterpCubicOperator`	Parametric cubic operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Ghat(\bfmu)[\qhat\otimes\qhat\otimes\qhat]\) where \(\Ghat(\bfmu) \in \RR^{r \times r^3}\) and the parametric dependence is handled with elementwise interpolation.
`InterpInputOperator`	Parametric input operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Bhat(\bfmu)\u\) where \(\Bhat(\bfmu) \in \RR^{r \times m}\) and the parametric dependence is handled with elementwise interpolation.
`InterpStateInputOperator`	Parametric state-input operator \(\Ophat_{\ell}(\qhat,\u;\bfmu) = \Nhat(\bfmu)[\u\otimes\qhat]\) where \(\Nhat(\bfmu) \in \RR^{r \times rm}\) and the parametric dependence is handled with elementwise interpolation.

Interpolatory operators can be instantiated with no arguments.

B = opinf.operators.InterpInputOperator()
print(B)

InterpInputOperator
  state_dimension:     None
  input_dimension:     None
  parameter_dimension: None
  training parameters: None
  type(interpolator):  None
  evaluate(parameter) -> InputOperator

s = 9  # Number of training parameters
p = 1  # Dimension of the training parameters.
r = 4  # Dimension of the states.
m = 2  # Dimension of the inputs.

training_parameters = np.random.standard_normal((s, p))
operator_matrices = [np.random.random((r, m)) for _ in range(s)]

B.set_training_parameters(training_parameters)
print(B)

InterpInputOperator
  state_dimension:     None
  input_dimension:     None
  parameter_dimension: 1
  training parameters: 9
  type(interpolator):  None
  evaluate(parameter) -> InputOperator

B.set_entries(operator_matrices, fromblock=False)
print(B)

InterpInputOperator
  state_dimension:     4
  input_dimension:     2
  parameter_dimension: 1
  training parameters: 9
  type(interpolator):  CubicSpline
  evaluate(parameter) -> InputOperator

B_nonparametric = B.evaluate(np.random.standard_normal(p))
print(B_nonparametric)

InputOperator
  state_dimension: 4
  input_dimension: 2
  entries.shape:   (4, 2)

opinf.operators

Contents

opinf.operators#

Operators#

Nonparametric Operators#

Learning Operators from Data#

Learning Operators via Projection#

Custom Nonparametric Operators#

Example: Hadamard Product with a Fixed Vector#

Example: Weighted Hadamard Input Operator#

Parametric Operators#

Affine Operators#

Interpolatory Operators#

`opinf.operators`#