<repo_name>swischuk/operator_inference<file_sep>/Makefile
.PHONY: test clean

REMOVE = rm -rfv
PYTHON = python3
PYTEST = pytest --cov --cov-report html
TARGET = tests/*.py # test_core.py


test:
	$(PYTHON) check_docs.py
	$(PYTEST) $(TARGET)
	open htmlcov/index.html

clean:
	find . -type d -name "*.egg*" | xargs $(REMOVE)
	find . -type f -name ".coverage*" | xargs $(REMOVE)
	find . -type d -name ".pytest_cache" | xargs $(REMOVE)
	find . -type d -name "__pycache__" | xargs $(REMOVE)
	find . -type d -name ".ipynb_checkpoints" | xargs $(REMOVE)
	find . -type d -name "htmlcov" | xargs $(REMOVE)
<file_sep>/examples/project_and_save.py
# project_and_save.py
"""An old data management script for computing a POD basis, projecting data,
computing velocities, and saving the projected data for later.
"""

import os
import sys
import h5py
import numpy as np
from sklearn.utils.extmath import randomized_svd

import rom_operator_inference as roi


def offdiag_penalizers(reg, r, m):
    """Construct a list of regularization matrices that penalize the
    off-diagonal elements of A and all elements of the remaining operators,
    where the model has the form dx / dt = c + Ax + H(x⊗x) + Bu.
    """
    regI = np.sqrt(reg) * np.eye(1 + r + r*(r+1)//2 + m)
    Gs = [regI] * r
    for i in range(1, r+1):
        Gs[i][i,i] = 0
    return Gs


def compute_randomized_svd(data, savepath, n_components=500):
    """Compute and save randomized SVD following
    https://scikit-learn.org/stable/modules/generated/sklearn.utils.extmath.randomized_svd.html

    Parameters
    ----------
    data : (n_variables, n_samples) ndarray
        Snapshot matrix for which the SVD will be computed.

    savepath : str
        Name of file to save the SVD to.

    n_components : float
        Number of singular values and vectors to compute and save.
    """
    Vr,Sigma,_ = randomized_svd(data, n_components, n_iter=15, random_state=42)

    with h5py.File(savepath, 'w') as svd_h5file:
        svd_h5file.create_dataset("U", data=Vr)
        svd_h5file.create_dataset("S", data=Sigma)

    print(f"SVD computed and saved in {savepath}")
    return Vr,Sigma


def project_data(data,
                 projected_data_folder,
                 projected_xdot_folder,
                 Vr,
                 dt,
                 r_vals=-1):
    """
    Compute xdot and project both data and xdot onto the POD basis.

    Projected data is saved as "<projected_data_folder>/data_reduced_%d.h5"
    Projected xdot data is saved as "<projected_xdot_folder>/xdot_reduced_%d.h5"

    Parameters
    ----------
    data : (n_variables, n_samples) ndarray
        Data to be projected and used to compute, project xdot.

    projected_data_folder : str
        Folder to save projected data to.

    projected_xdot_folder : str
        Folder to save projected xdot data to.

    Vr: (n_variables, r) ndarray
        The POD basis.

    dt: float
        Timestep size between data (for computing xdot)

    r_vals: list of ints <= Vr.shape[1] (the number of columns in Vr).
        Basis sizes to compute projection for.
    """

    if data.shape[1] < 2:
        raise ValueError("At least two snapshots required for x' computation")

    os.makedirs(os.path.dirname(projected_data_folder), exist_ok=True)
    os.makedirs(os.path.dirname(projected_xdot_folder), exist_ok=True)

    if r_vals == -1:
        r_vals = [Vr.shape[1]]

    for r in r_vals:
        print(f"r = {r}", flush=True)
        # select first r basis vectors
        VrT = Vr[:,:r].T

        # project snapshots
        print("\tProjecting snapshots...", end='', flush=True)
        data_reduced = VrT @ data
        print("done.")

        # project velocities
        print("\tComputing projected velocities...", end='', flush=True)
        if data.shape[1] <= 4:  # Too few time steps for 4th order derivatives
            xdot_reduced = np.gradient(data_reduced, dt, edge_order=2, axis=-1)
        else:
            xdot_reduced = roi.pre.xdot_uniform(data_reduced, dt, order=4)
        print("done.")

        # save results
        print("\tSaving files...", end='', flush=True)
        fname1 = os.path.join(projected_data_folder, f"data_reduced_{r:d}.h5")
        with open h5py.File(fname1,'w') as savexf:
            savexf.create_dataset('data', data=data_reduced)

        fname2 = os.path.join(projected_xdot_folder, f"xdot_reduced_{r:d}.h5")
        with open h5py.File(fname2,'w') as savexdotf:
            savexdotf.create_dataset('xdot', data=xdot_reduced)
        print("done.")


# Some custom error evaluators ------------------------------------------------
def mape(truth, pred):
    """Mean absolute prediction error."""
    mask = truth > 1e-12
    return (((np.abs(truth-pred)/truth)[mask]).sum()+pred[~mask].sum())/len(truth)


def rpe(truth,pred):
    """Relative prediction error."""
    mask = abs(truth) > 1e-10
    rpe_mat = np.zeros(mask.shape)
    rpe_mat[mask] = (np.abs(truth-pred)/abs(truth))[mask]
    rpe_mat[~mask] = abs(pred[~mask])
    return rpe_mat


if __name__ == '__main__':

    # make command line arguments for trainsize and whether to compute svd
    trainsize = int(input("How many snapshots?"))
    svdfile = os.path.join("data", f"svd_nt{trainsize}.h5")
    compute_svd = input(f"Compute SVD? (True or False -- Type False if it exists as '{svdfile}')")

    fulldatapath = "data/data_renee.h5"
    with h5py.File(fulldatapath,'r') as hf:
        data = hf['data'][:,:]

    if compute_svd != "False":
        print("Computing svd...", end='')
        V,_ = compute_randomized_svd(data[:,:trainsize], svdfile)
        print("done.")

    else:
        print("Loading svd...", end='')
        with h5py.File(svdfile, 'r') as hfs:
            V = hfs['Und'][:,:]
        print("done.")

    print("Computing xdot and reducing data...")
    target_folder = os.path.join("data", f"data_reduced_minmax_nt{trainsize}")
    project_data(data[:,:trainsize],
                 target_folder,
                 target_folder,
                 V,
                 1e-7,
                 [5,10,15])
<file_sep>/README.md
# Operator Inference

This is a Python implementation of Operator Inference for constructing projection-based reduced-order models of dynamical systems with a polynomial form.
The procedure is **data-driven** and **non-intrusive**, making it a viable candidate for model reduction of black-box or complex systems.
The methodology was introduced in [\[1\]](https://www.sciencedirect.com/science/article/pii/S0045782516301104).
See [**References**](#references) for more papers that use or build on Operator Inference.

**Contributors**: [<NAME>](https://github.com/swischuk), [<NAME>](https://github.com/shanemcq18), [<NAME>](https://github.com/elizqian), [<NAME>](https://github.com/bokramer), [<NAME>](https://kiwi.oden.utexas.edu/).

See [this repository](https://github.com/elizqian/operator-inference) for a MATLAB implementation and [DOCUMENTATION.md](DOCUMENTATION.md) for the code documentation.

## Problem Statement

Consider the (possibly nonlinear) system of _n_ ordinary differential equations with state variable **x**, input (control) variable **u**, and independent variable _t_:

<p align="center"><img src="img/prb/eq1.svg"/></p>

where

<p align="center"><img src="img/prb/eq2.svg"/></p>

This system is called the _full-order model_ (FOM).
If _n_ is large, as it often is in high-consequence engineering applications, it is computationally expensive to numerically solve the FOM.
This package provides tools for constructing a _reduced-order model_ (ROM) that is up to quadratic in the state **x** with optional linear control inputs **u**.
The procedure is data-driven, non-intrusive, and relatively inexpensive.
In the most general case, the code can construct and solve a reduced-order system with the polynomial form

<p align="center"><img src="img/prb/eq3.svg"/></p>

<!-- <p align="center">
  <img src="https://latex.codecogs.com/svg.latex?\dot{\hat{\mathbf{x}}}(t)=\hat{A}\hat{\mathbf{x}}(t)+\hat{H}(\hat{\mathbf{x}}\otimes\hat{\mathbf{x}})(t)+\hat{B}\mathbf{u}(t)+\sum_{i=1}^m\hat{N}_{i}\hat{\mathbf{x}}(t)u_{i}(t)+\hat{\mathbf{c}},"/>
</p> -->

where now

<p align="center"><img src="img/prb/eq4.svg"/></p>
<p align="center"><img src="img/prb/eq5.svg"/></p>

<!-- <p align="center">
  <img src="https://latex.codecogs.com/svg.latex?\hat{A}\in\mathbb{R}^{r\times%20r},\qquad\hat{H}\in\mathbb{R}^{r\times%20r^2},\qquad\hat{B}\in\mathbb{R}^{r\times%20m},\qquad\hat{N}_{i}\in\mathbb{R}^{r\times%20r}."/>
</p> -->

This reduced low-dimensional system approximates the original high-dimensional system, but it is much easier (faster) to solve because of its low dimension _r_ << _n_.

See [DETAILS.md](DETAILS.md) for more mathematical details and an index of notation.


## Quick Start

#### Installation

Install from the command line with the following single command (requires [`pip`](https://pypi.org/project/pip/) and [`git`](https://git-scm.com/)).
```bash
$ pip3 install git+https://github.com/swischuk/rom-operator-inference-Python3.git
```

#### Usage

Given a linear basis `Vr`, snapshot data `X`, and snapshot velocities `Xdot`, the following code learns a reduced model for a problem of the form _d**x**/dt = **c** + A**x**(t)_, then solves the reduced system for _0 ≤ t ≤ 1_.

```python
import numpy as np
import rom_operator_inference as roi

# Define a model of the form  dx / dt = c + Ax(t).
>>> model = roi.InferredContinuousROM(modelform="cA")

# Fit the model to snapshot data X, the snapshot derivative Xdot,
# and the linear basis Vr by solving for the operators c_ and A_.
>>> model.fit(Vr, X, Xdot)

# Simulate the learned model over the time domain [0,1] with 100 timesteps.
>>> t = np.linspace(0, 1, 100)
>>> X_ROM = model.predict(X[:,0], t)
```


## Examples

The [`examples/`](examples/) folder contains scripts and notebooks that set up and run several examples:
- [`examples/tutorial.ipynb`](https://nbviewer.jupyter.org/github/Willcox-Research-Group/rom-operator-inference-Python3/blob/master/examples/tutorial.ipynb): A walkthrough of a very simple heat equation example.
- [`examples/heat_1D.ipynb`](https://nbviewer.jupyter.org/github/Willcox-Research-Group/rom-operator-inference-Python3/blob/master/examples/heat_1D.ipynb): A more complicated one-dimensional heat equation example [\[1\]](https://www.sciencedirect.com/science/article/pii/S0045782516301104).
- [`examples/data_driven_heat.ipynb`](https://nbviewer.jupyter.org/github/Willcox-Research-Group/rom-operator-inference-Python3/blob/master/examples/data_driven_heat.ipynb): A purely data-driven example using data generated from a one-dimensional heat equation \[4\].
<!-- - `examples/TODO.ipynb`: Burgers' equation [\[1\]](https://www.sciencedirect.com/science/article/pii/S0045782516301104). -->
<!-- - `examples/TODO.ipynb`: Euler equation [\[2\]](https://arc.aiaa.org/doi/10.2514/6.2019-3707). -->
<!-- This example uses MATLAB's Curve Fitting Toolbox to generate the random initial conditions. -->

(More examples coming)


## References

- \[1\] <NAME>. and <NAME>.,
[Data-driven operator inference for non-intrusive projection-based model reduction.](https://www.sciencedirect.com/science/article/pii/S0045782516301104)
_Computer Methods in Applied Mechanics and Engineering_, Vol. 306, pp. 196-215, 2016.
([Download](https://kiwi.oden.utexas.edu/papers/Non-intrusive-model-reduction-Peherstorfer-Willcox.pdf))<details><summary>BibTeX</summary><pre>
@article{Peherstorfer16DataDriven,
    title     = {Data-driven operator inference for nonintrusive projection-based model reduction},
    author    = {<NAME>. and <NAME>.},
    journal   = {Computer Methods in Applied Mechanics and Engineering},
    volume    = {306},
    pages     = {196--215},
    year      = {2016},
    publisher = {Elsevier}
}</pre></details>

- \[2\] <NAME>., <NAME>., <NAME>., and <NAME>.,
[Transform & Learn: A data-driven approach to nonlinear model reduction](https://arc.aiaa.org/doi/10.2514/6.2019-3707).
In the AIAA Aviation 2019 Forum & Exhibition, Dallas, TX, June 2019. ([Download](https://kiwi.oden.utexas.edu/papers/learn-data-driven-nonlinear-reduced-model-Qian-Willcox.pdf))<details><summary>BibTeX</summary><pre>
@inbook{QKMW2019aviation,
    title     = {Transform \\& Learn: A data-driven approach to nonlinear model reduction},
    author    = {<NAME>. and <NAME>. and <NAME>. and <NAME>.},
    booktitle = {AIAA Aviation 2019 Forum},
    doi       = {10.2514/6.2019-3707},
    URL       = {https://arc.aiaa.org/doi/abs/10.2514/6.2019-3707},
    eprint    = {https://arc.aiaa.org/doi/pdf/10.2514/6.2019-3707}
}</pre></details>

- \[3\] <NAME>., <NAME>., <NAME>., and <NAME>.,
[Projection-based model reduction: Formulations for physics-based machine learning.](https://www.sciencedirect.com/science/article/pii/S0045793018304250)
_Computers & Fluids_, Vol. 179, pp. 704-717, 2019.
([Download](https://kiwi.oden.utexas.edu/papers/Physics-based-machine-learning-swischuk-willcox.pdf))<details><summary>BibTeX</summary><pre>
@article{swischuk2019projection,
    title     = {Projection-based model reduction: Formulations for physics-based machine learning},
    author    = {<NAME>. and <NAME>. and <NAME>. and <NAME>.},
    journal   = {Computers \\& Fluids},
    volume    = {179},
    pages     = {704--717},
    year      = {2019},
    publisher = {Elsevier}
}</pre></details>

- \[4\] <NAME>., [Physics-based machine learning and data-driven reduced-order modeling](https://dspace.mit.edu/handle/1721.1/122682). Master's thesis, Massachusetts Institute of Technology, 2019. ([Download](https://dspace.mit.edu/bitstream/handle/1721.1/122682/1123218324-MIT.pdf))<details><summary>BibTeX</summary><pre>
@phdthesis{swischuk2019physics,
    title  = {Physics-based machine learning and data-driven reduced-order modeling},
    author = {<NAME>},
    year   = {2019},
    school = {Massachusetts Institute of Technology}
}</pre></details>

- \[5\] <NAME>. [Sampling low-dimensional Markovian dynamics for pre-asymptotically recovering reduced models from data with operator inference](https://arxiv.org/abs/1908.11233). arXiv:1908.11233.
([Download](https://arxiv.org/pdf/1908.11233.pdf))<details><summary>BibTeX</summary><pre>
@article{peherstorfer2019sampling,
    title   = {Sampling low-dimensional Markovian dynamics for pre-asymptotically recovering reduced models from data with operator inference},
    author  = {<NAME>},
    journal = {arXiv preprint arXiv:1908.11233},
    year    = {2019}
}</pre></details>

- \[6\] <NAME>., <NAME>., <NAME>., and <NAME>., [Learning physics-based reduced-order models for a single-injector combustion process](https://arc.aiaa.org/doi/10.2514/1.J058943). _AIAA Journal_, published online March 2020. Also in Proceedings of 2020 AIAA SciTech Forum & Exhibition, Orlando FL, January, 2020. Also Oden Institute Report 19-13. ([Download](https://kiwi.oden.utexas.edu/papers/learning-reduced-model-combustion-Swischuk-Kramer-Huang-Willcox.pdf))<details><summary>BibTeX</summary><pre>
@article{SKHW2019_learning_ROMs_combustor,
    title   = {Learning physics-based reduced-order models for a single-injector combustion process},
    author  = {<NAME>. and <NAME>. and <NAME>. and <NAME>.},
    journal = {AIAA Journal},
    volume  = {},
    pages   = {Published Online: 19 Mar 2020},
    url     = {},
    year    = {2020}
}</pre></details>

- \[7\] <NAME>., <NAME>., <NAME>., and <NAME>. [Lift & Learn: Physics-informed machine learning for large-scale nonlinear dynamical systems](https://www.sciencedirect.com/science/article/abs/pii/S0167278919307651). _Physica D: Nonlinear Phenomena_, Volume 406, May 2020, 132401. ([Download](https://kiwi.oden.utexas.edu/papers/lift-learn-scientific-machine-learning-Qian-Willcox.pdf))<details><summary>BibTeX</summary><pre>
@article{QKPW2020_lift_and_learn,
    title   = {Lift \\& Learn: Physics-informed machine learning for large-scale nonlinear dynamical systems.},
    author  = {<NAME>. and <NAME>. and <NAME>. and <NAME>.},
    journal = {Physica {D}: {N}onlinear {P}henomena},
    volume  = {406},
    pages   = {132401},
    url     = {https://doi.org/10.1016/j.physd.2020.132401},
    year    = {2020}
}</pre></details>
<file_sep>/check_docs.py
# check_docs.py
"""Script for catching errors in README and DETAILS files."""

import os
import re
import difflib


# Global variables ============================================================
README = "README.md"
DETAILS = "DETAILS.md"
DOCS = "DOCUMENTATION.md"
SETUP = "setup.py"
INIT = os.path.join("rom_operator_inference", "__init__.py")
MDFILES = [README, DETAILS, DOCS]
DOCFILES = [DETAILS, DOCS]
VNFILES = [SETUP, INIT]

REDSPACE = '\x1b[41m \x1b[49m'
VERSION = re.compile(r'_{0,2}?version_{0,2}?\s*=\s*"([\d\.]+?)"',
                     re.MULTILINE)
TEX = re.compile(
                r'''<img\ src=                          #  begin image tag
                "https://latex.codecogs.com/svg.latex\? #   LaTeX web prefix
                (.+?)                                   #   LATEX CODE
                "/>                                     #  end image tag
                ''',
                re.VERBOSE | re.MULTILINE)


# Tests =======================================================================
def check_latex_code_for_spaces(filelist=MDFILES):
    """Make sure there are no spaces in the latex code snippets, since that
    keeps them from displaying correctly on Github.
    """
    # Read files and search for errors.
    errors = []
    for filename in filelist:
        with open(filename, 'r') as infile:
            data = infile.read()
        for index,tex in enumerate(TEX.findall(data)):
            if ' ' in tex:
                errors.append((filename, index, tex))

    # Report errors.
    for name,index,tex in errors:
        print(f"bad space in {name}, LaTeX occurrence #{index+1}:\n\t",
              tex.replace(' ', REDSPACE),
              sep='', end='\n')
    if errors:
        nerrs = len(errors)
        raise SyntaxError(f"{nerrs} LaTeX error{'s' if nerrs > 1 else ''}")


def check_references_sections_are_the_same(filelist=MDFILES):
    """Make sure that the "References" sections in each doc file are the same.
    """
    # Read both '## References' sections.
    refsections = []
    for filename in filelist:
        with open(filename, 'r') as infile:
            data = infile.read()
        refsections.append(data[data.index("## References"):].splitlines())

    # Compare sections and report errors.
    errors = False
    for i in range(len(filelist)-1):
        file1, file2 = filelist[i:i+2]
        for line in difflib.unified_diff(refsections[i], refsections[i+1],
                                         fromfile=file1, tofile=file2):
            print(line)
            errors = True
        if errors:
            raise SyntaxError(f"'References' of {file1} and {file2}"
                               " do not match")


def check_index_sections_are_the_same(filelist=DOCFILES):
    """Make sure that the "Index of Notation" sections in each doc file are
    the same.
    """
    idxsections = []
    for filename in filelist:
        with open(filename, 'r') as infile:
            data = infile.read()
        start = data.index("## Index of Notation")
        end = data.index("## References")
        idxsections.append(data[start:end].splitlines())

    # Compare sections and report errors.
    errors = False
    for i in range(len(filelist)-1):
        file1, file2 = filelist[i:i+2]
        for line in difflib.unified_diff(idxsections[i], idxsections[i+1],
                                         fromfile=file1, tofile=file2):
            print(line)
            errors = True
        if errors:
            raise SyntaxError(f"'Index of Notation' of {file1} and {file2} "
                               " do not match")


def check_version_numbers_match(filelist=VNFILES):
    """Make sure that the version number in setup.py and __init__.py match."""
    if len(filelist) != 2:
        raise ValueError("can only compare 2 files at a time")
    file1, file2 = filelist

    # Get the version number listed in each file.
    versions = []
    for filename in filelist:
        with open(filename, 'r') as infile:
            data = infile.read()
        versions.append(VERSION.findall(data)[0])

    if versions[0] != versions[1]:
        raise ValueError(f"Version numbers in {file1} and {file2} do "
                         f"not match ('{versions[0]}' != '{versions[1]}')")


if __name__ == "__main__":
    check_latex_code_for_spaces()
    check_references_sections_are_the_same()
    check_index_sections_are_the_same()
    check_version_numbers_match()
<file_sep>/DETAILS.md
# Summary of Mathematical Details

This document gives a short explanation of the mathematical details behind the package.
For a full treatment, see [\[1\]](https://www.sciencedirect.com/science/article/pii/S0045782516301104).
However, note that some notation has been altered for coding convenience and clarity.

**Contents**
- [**Problem Statement**](#problem-statement)
- [**Projection-based Model Reduction**](#projection-based-model-reduction)
- [**Operator Inference via Least Squares**](#operator-inference-via-least-squares)
- [**Index of Notation**](#index-of-notation)
- [**References**](#references)


## Problem Statement

Consider the (possibly nonlinear) system of _n_ ordinary differential equations with state variable **x**, input (control) variable **u**, and independent variable _t_:

<p align="center"><img src="img/prb/eq1.svg"/></p>

where

<p align="center"><img src="img/prb/eq2.svg"/></p>

This system is called the _full-order model_ (FOM).
If _n_ is large, as it often is in high-consequence engineering applications, it is computationally expensive to numerically solve the FOM.
This package provides tools for constructing a _reduced-order model_ (ROM) that is up to quadratic in the state **x** with optional linear control inputs **u**.
The procedure is data-driven, non-intrusive, and relatively inexpensive.
In the most general case, the code can construct and solve a reduced-order system with the polynomial form

<p align="center"><img src="img/prb/eq3.svg"/></p>

<!-- <p align="center">
  <img src="https://latex.codecogs.com/svg.latex?\dot{\hat{\mathbf{x}}}(t)=\hat{A}\hat{\mathbf{x}}(t)+\hat{H}(\hat{\mathbf{x}}\otimes\hat{\mathbf{x}})(t)+\hat{B}\mathbf{u}(t)+\sum_{i=1}^m\hat{N}_{i}\hat{\mathbf{x}}(t)u_{i}(t)+\hat{\mathbf{c}},"/>
</p> -->

where now

<p align="center"><img src="img/prb/eq4.svg"/></p>
<p align="center"><img src="img/prb/eq5.svg"/></p>

<!-- <p align="center">
  <img src="https://latex.codecogs.com/svg.latex?\hat{A}\in\mathbb{R}^{r\times%20r},\qquad\hat{H}\in\mathbb{R}^{r\times%20r^2},\qquad\hat{B}\in\mathbb{R}^{r\times%20m},\qquad\hat{N}_{i}\in\mathbb{R}^{r\times%20r}."/>
</p> -->

This reduced low-dimensional system approximates the original high-dimensional system, but it is much easier (faster) to solve because of its low dimension _r_ << _n_.


## Projection-based Model Reduction

Model reduction via projection occurs in three steps:
1. **Data Collection**: Gather snapshot data, i.e., solutions to the full-order model (the FOM) at various times / parameters.
2. **Compression**: Compute a low-rank basis (which defines a low-dimensional linear subspace) that captures most of the behavior of the snapshots.
3. **Projection**: Use the low-rank basis to construct a low-dimensional ODE (the ROM) that approximates the FOM.

<!-- These steps comprise what is called the _offline phase_ in the literature, since they can all be done before the resulting ROM is simulated. -->

This package focuses on step 3, constructing the ROM given the snapshot data and the low-rank basis from steps 1 and 2, respectively.

Let _X_ be the _n_ x _k_ matrix whose _k_ columns are each solutions to the FOM of length _n_ (step 1), and let _V_<sub>_r_</sub> be an orthonormal _n_ x _r_ matrix representation for an _r_-dimensional subspace (step 2).
For example, a common choice for _V_<sub>_r_</sub> is the POD Basis of rank _r_, the matrix comprised of the first _r_ singular vectors of _X_.
We call _X_ the _snapshot matrix_ and _V_<sub>_r_</sub> the _reduced basis matrix_.

The classical approach to the projection step is to make the Ansatz

<p align="center"><img src="img/dtl/eq01.svg"/></p>

Inserting this into the FOM and multiplying both sides by the transpose of _V_<sub>_r_</sub> yields

<p align="center"><img src="img/dtl/eq02.svg"/></p>

This new system is _r_-dimensional in the sense that

<p align="center"><img src="img/dtl/eq03.svg"/></p>

If the FOM operator **f** is known and has a nice structure, this reduced system can be solved cheaply by precomputing any involved matrices and then applying a time-stepping scheme.
For example, if **f** is linear in **x** and there is no input **u**, then

<p align="center"><img src="img/dtl/eq04.svg"/></p>

where

<p align="center"><img src="img/dtl/eq05.svg"/></p>

However, _this approach breaks down if the FOM operator **f** is unknown, uncertain, or highly nonlinear_.

## Operator Inference via Least Squares

Instead of directly computing the reduced operators, the Operator Inference framework takes a data-driven approach: assuming a specific structure of the ROM (linear, quadratic, etc.), solve for the involved operators that best fit the data.
For example, suppose that we seek a ROM of the form

<p align="center"><img src="img/dtl/eq06.svg"/></p>

We have only the snapshot matrix _X_, the low-rank basis matrix _V_<sub>_r_</sub> (which was derived from _X_), the inputs _U_, and perhaps the snapshot velocities _X'_ (if not, these must be approximated).
Here the (_ij_)<sup>th</sup> entry of _U_ is the _i_<sup>th</sup> component of **u** at the time corresponding to the _j_<sup>th</sup> snapshot.
To solve for the linear operators on the right-hand side of the preceding equation, we project the snapshot data via the basis matrix,

<p align="center"><img src="img/dtl/eq07.svg"/></p>

then solve the least squares problem

<p align="center"><img src="img/dtl/eq08.svg"/></p>

where **1** is a _k_-vector of 1's and

<p align="center"><img src="img/dtl/eq09.svg"/></p>

For our purposes, the ⊗ operator between matrices denotes a column-wise Kronecker product (also called the [Khatri-Rao product](https://en.wikipedia.org/wiki/Kronecker_product#Khatri%E2%80%93Rao_product)).
The minimization problem given above decouples into _r_ independent ordinary least-squares problems, one for each of the columns of _O<sup>T</sup>_:

<p align="center"><img src="img/dtl/eq10.svg"/></p>

The entire routine is relatively inexpensive to solve.
The code also allows for a Tikhonov regularization matrix or list of matrices (the `P` keyword argument for `predict()` methods), in which case the problem being solved is

<p align="center"><img src="img/dtl/eq11.svg"/></p>

It can be shown [\[1\]](https://www.sciencedirect.com/science/article/pii/S0045782516301104) that, under some idealized assumptions, these inferred operators converge to the operators computed by explicit projection.
The key idea, however, is that _the inferred operators can be cheaply computed without knowing the full-order model_.
This is very convenient in situations where the FOM is given by a "black box," such as a legacy code for complex fluid simulations.

#### The Discrete Case

The framework described above can also be used to construct reduced-order models for approximating _discrete_ dynamical systems.
For instance, consider the full-order model

<p align="center"><img src="img/dtl/eq12.svg"/></p>

Instead of collecting snapshot velocities, we collect _k+1_ snapshots and let _X_ be the _n x k_ matrix whose columns are the first _k_ snapshots and _X'_ be the _n x k_ matrix whose columns are the last _k_ snapshots.
That is, the columns **x**<sub>_k_</sub> of _X_ and **x**<sub>_k_</sub>' satisfy

<p align="center"><img src="img/dtl/eq13.svg"/></p>

Then we set up the same least squares problem as before, but now the right-hand side matrix is

<p align="center"><img src="img/dtl/eq14.svg"/></p>

The resulting reduced-order model has the form

<p align="center"><img src="img/dtl/eq15.svg"><p>


<!-- TODO: #### Re-projection and Recovering Intrusive Models -->


#### Implementation Note: The Kronecker Product

The vector [Kronecker product](https://en.wikipedia.org/wiki/Kronecker_product) ⊗ introduces some redundancies.
For example, the product **x** ⊗ **x** contains both _x_<sub>1</sub>_x_<sub>2</sub> and _x_<sub>2</sub>_x_<sub>1</sub>.
To avoid these redundancies, we introduce a "compact" Kronecker product <img src="img/dtl/eq16.svg" height=10/> which only computes the unique terms of the usual vector Kronecker product:

<p align="center"><img src="img/dtl/eq17.svg"/></p>

When the compact Kronecker product is used, we call the resulting operator _H<sub>c</sub>_ instead of _H_.
Thus, the reduced order model becomes

<p align="center"><img src="img/dtl/eq18.svg"/></p>

and the corresponding Operator Inference least squares problem is

<p align="center"><img src="img/dtl/eq19.svg"/></p>


## Index of Notation

We generally denote scalars in lower case, vectors in bold lower case, matrices in upper case, and indicate low-dimensional quantities with a hat.
In the code, a low-dimensional quantity ends with an underscore, so that the model classes follow some principles from the [scikit-learn](https://scikit-learn.org/stable/index.html) [API](https://scikit-learn.org/stable/developers/contributing.html#apis-of-scikit-learn-objects).

### Dimensions

| Symbol | Code | Description |
| :----: | :--- | :---------- |
| <img src="img/ntn/n.svg"/> | `n`  | Dimension of the full-order system (large) |
| <img src="img/ntn/r.svg"/> | `r`  | Dimension of the reduced-order system (small) |
| <img src="img/ntn/m.svg"/> | `m`  | Dimension of the input **u** |
| <img src="img/ntn/k.svg"/> | `k`  | Number of state snapshots, i.e., the number of training points |
| <img src="img/ntn/s.svg"/> | `s`  | Number of parameter samples for parametric training |
| <img src="img/ntn/p.svg"/> | `p` | Dimension of the parameter space |
| <img src="img/ntn/d.svg"/> | `d` | Number of columns of the data matrix _D_ |
| <img src="img/ntn/nt.svg"/> | `nt`  | Number of time steps in a simulation |

<!-- | <img src="https://latex.codecogs.com/svg.latex?\ell"/> | `l` | Dimension of the output **y** | -->


### Vectors

<!-- \sigma_j\in\text{diag}(\Sigma) &= \textrm{singular value of }X\\
\boldsymbol{\mu}\in\mathcal{P} &= \text{system parameter}\\
\mathcal{P}\subset\mathbb{R}^{p} &= \text{parameter space}\\
\Omega\subset\mathbb{R}^{d} &= \text{spatial domain}\\
% \omega\in\Omega &= \text{spatial point (one dimension)}\\
\boldsymbol{\omega}\in\Omega &= \text{spatial point}\\
t\ge 0 &= \text{time}\\
\hat{} &= \textrm{reduced variable, e.g., }\hat{\mathbf{x}}\textrm{ or }\hat{A}\\
\dot{} = \frac{d}{dt} &= \text{time derivative} -->


| Symbol | Code | Size | Description |
| :----: | :--- | :--: | :---------- |
| <img src="img/ntn/x.svg"/> | `x` | <img src="img/ntn/n.svg"/> | Full-order state vector |
| <img src="img/ntn/xhat.svg"/> | `x_` | <img src="img/ntn/r.svg"/> | Reduced-order state vector |
| <img src="img/ntn/xhatdot.svg"/> | `xdot_` | <img src="img/ntn/r.svg"/> | Reduced-order state velocity vector |
| <img src="img/ntn/xrom.svg"/> | `x_ROM` | <img src="img/ntn/n.svg"/> | Approximation to **x** produced by ROM |
| <img src="img/ntn/chat.svg"/> | `c_` | <img src="img/ntn/m.svg"/> | Learned constant term  |
| <img src="img/ntn/u.svg"/> | `u` | <img src="img/ntn/m.svg"/> | Input vector  |
| <img src="img/ntn/f.svg"/> | `f(t,x,u)` | <img src="img/ntn/n.svg"/>  | Full-order system operator |
| <img src="img/ntn/fhat.svg"/> | `f_(t,x_,u)` | <img src="img/ntn/n.svg"/>  | Reduced-order system operator |
| <img src="img/ntn/kronx.svg"/> | `np.kron(x,x)` | <img src="img/ntn/n2.svg"/> | Kronecker product of full state (quadratic terms) |
| <img src="img/ntn/kronxhat.svg"/> | `np.kron(x_,x_)` | <img src="img/ntn/r2.svg"/>  | Kronecker product of reduced state (quadratic terms) |
| <img src="img/ntn/kronxhatc.svg"/> | `kron2c(x_)` | <img src="img/ntn/r2c.svg"/>  | Compact Kronecker product of reduced state (quadratic terms) |
| <img src="img/ntn/vj.svg"/> | `vj` | <img src="img/ntn/n.svg"/> | _j_<sup>th</sup> subspace basis vector, i.e., column _j_ of _V_<sub>_r_</sub> |

<!-- | **y**  | `y`             | Output vector | -->
<!-- | **y_ROM**, **y~** | `y_ROM`      | Approximation to **y** produced by ROM | -->

### Matrices

| Symbol | Code | Shape | Description |
| :----: | :--- | :---: | :---------- |
| <img src="img/ntn/Vr.svg"/> | `Vr` | <img src="img/ntn/nxr.svg"/> | low-rank basis of rank _r_ (usually the POD basis) |
| <img src="img/ntn/XX.svg"/> | `X` | <img src="img/ntn/nxk.svg"/> | Snapshot matrix |
| <img src="img/ntn/XXdot.svg"/> | `Xdot` | <img src="img/ntn/nxk.svg"/> | Snapshot velocity matrix |
| <img src="img/ntn/UU.svg"/> | `U` | <img src="img/ntn/mxk.svg"/> | Input matrix (inputs corresonding to the snapshots) |
| <img src="img/ntn/XXhat.svg"/> | `X_` | <img src="img/ntn/rxk.svg"/> | Projected snapshot matrix |
| <img src="img/ntn/XXhatdot.svg"/> | `Xdot_` | <img src="img/ntn/rxk.svg"/> | Projected snapshot velocity matrix |
| <img src="img/ntn/DD.svg"/> | `D` | <img src="img/ntn/kxd.svg"/> | Data matrix |
| <img src="img/ntn/OO.svg"/> | `O` | <img src="img/ntn/dxr.svg"/> | Operator matrix |
| <img src="img/ntn/RR.svg"/> | `R` | <img src="img/ntn/kxr.svg"/> | Right-hand side matrix |
| <img src="img/ntn/PP.svg"/> | `P` | <img src="img/ntn/dxd.svg"/> | Tikhonov regularization matrix |
| <img src="img/ntn/AAhat.svg"/> | `A_` | <img src="img/ntn/rxr.svg"/> | Learned state matrix |
| <img src="img/ntn/HHhat.svg"/> | `H_` | <img src="img/ntn/rxr2.svg"/> | Learned matricized quadratic tensor |
| <img src="img/ntn/HHhatc.svg"/> | `Hc_` | <img src="img/ntn/rxr2c.svg"/> | Learned matricized quadratic tensor without redundancy (compact) |
| <img src="img/ntn/BBhat.svg"/> | `B_` | <img src="img/ntn/rxm.svg"/> | Learned input matrix |

<!-- | <img src="https://latex.codecogs.com/svg.latex?\hat{N}_i"/> | `Ni_` | <img src="https://latex.codecogs.com/svg.latex?r\times%20r"/> | Bilinear state-input matrix for _i_th input | -->

<!-- | <img src="https://latex.codecogs.com/svg.latex?\hat{C}"/> | `C_` | <img src="https://latex.codecogs.com/svg.latex?q\times%20r"/> | Learned output matrix | -->

<!-- I_{a\times%20a}\in\mathbb{R}^{a\times a} | | identity matrix\\ -->
<!-- \Sigma \in \mathbb{R}^{\ell\times\ell} &= \text{diagonal singular value matrix}\\ -->


## References

- \[1\] <NAME>. and <NAME>.,
[Data-driven operator inference for non-intrusive projection-based model reduction.](https://www.sciencedirect.com/science/article/pii/S0045782516301104)
_Computer Methods in Applied Mechanics and Engineering_, Vol. 306, pp. 196-215, 2016.
([Download](https://kiwi.oden.utexas.edu/papers/Non-intrusive-model-reduction-Peherstorfer-Willcox.pdf))<details><summary>BibTeX</summary><pre>
@article{Peherstorfer16DataDriven,
    title     = {Data-driven operator inference for nonintrusive projection-based model reduction},
    author    = {<NAME>. and <NAME>.},
    journal   = {Computer Methods in Applied Mechanics and Engineering},
    volume    = {306},
    pages     = {196--215},
    year      = {2016},
    publisher = {Elsevier}
}</pre></details>

- \[2\] <NAME>., <NAME>., <NAME>., and <NAME>.,
[Transform & Learn: A data-driven approach to nonlinear model reduction](https://arc.aiaa.org/doi/10.2514/6.2019-3707).
In the AIAA Aviation 2019 Forum & Exhibition, Dallas, TX, June 2019. ([Download](https://kiwi.oden.utexas.edu/papers/learn-data-driven-nonlinear-reduced-model-Qian-Willcox.pdf))<details><summary>BibTeX</summary><pre>
@inbook{QKMW2019aviation,
    title     = {Transform \\& Learn: A data-driven approach to nonlinear model reduction},
    author    = {<NAME>. and <NAME>. and <NAME>. and <NAME>.},
    booktitle = {AIAA Aviation 2019 Forum},
    doi       = {10.2514/6.2019-3707},
    URL       = {https://arc.aiaa.org/doi/abs/10.2514/6.2019-3707},
    eprint    = {https://arc.aiaa.org/doi/pdf/10.2514/6.2019-3707}
}</pre></details>

- \[3\] <NAME>., <NAME>., <NAME>., and <NAME>.,
[Projection-based model reduction: Formulations for physics-based machine learning.](https://www.sciencedirect.com/science/article/pii/S0045793018304250)
_Computers & Fluids_, Vol. 179, pp. 704-717, 2019.
([Download](https://kiwi.oden.utexas.edu/papers/Physics-based-machine-learning-swischuk-willcox.pdf))<details><summary>BibTeX</summary><pre>
@article{swischuk2019projection,
    title     = {Projection-based model reduction: Formulations for physics-based machine learning},
    author    = {<NAME>. and <NAME>. and <NAME>. and <NAME>.},
    journal   = {Computers \\& Fluids},
    volume    = {179},
    pages     = {704--717},
    year      = {2019},
    publisher = {Elsevier}
}</pre></details>

- \[4\] <NAME>., [Physics-based machine learning and data-driven reduced-order modeling](https://dspace.mit.edu/handle/1721.1/122682). Master's thesis, Massachusetts Institute of Technology, 2019. ([Download](https://dspace.mit.edu/bitstream/handle/1721.1/122682/1123218324-MIT.pdf))<details><summary>BibTeX</summary><pre>
@phdthesis{swischuk2019physics,
    title  = {Physics-based machine learning and data-driven reduced-order modeling},
    author = {<NAME>},
    year   = {2019},
    school = {Massachusetts Institute of Technology}
}</pre></details>

- \[5\] <NAME>. [Sampling low-dimensional Markovian dynamics for pre-asymptotically recovering reduced models from data with operator inference](https://arxiv.org/abs/1908.11233). arXiv:1908.11233.
([Download](https://arxiv.org/pdf/1908.11233.pdf))<details><summary>BibTeX</summary><pre>
@article{peherstorfer2019sampling,
    title   = {Sampling low-dimensional Markovian dynamics for pre-asymptotically recovering reduced models from data with operator inference},
    author  = {<NAME>},
    journal = {arXiv preprint arXiv:1908.11233},
    year    = {2019}
}</pre></details>

- \[6\] <NAME>., <NAME>., <NAME>., and <NAME>., [Learning physics-based reduced-order models for a single-injector combustion process](https://arc.aiaa.org/doi/10.2514/1.J058943). _AIAA Journal_, published online March 2020. Also in Proceedings of 2020 AIAA SciTech Forum & Exhibition, Orlando FL, January, 2020. Also Oden Institute Report 19-13. ([Download](https://kiwi.oden.utexas.edu/papers/learning-reduced-model-combustion-Swischuk-Kramer-Huang-Willcox.pdf))<details><summary>BibTeX</summary><pre>
@article{SKHW2019_learning_ROMs_combustor,
    title   = {Learning physics-based reduced-order models for a single-injector combustion process},
    author  = {<NAME>. and <NAME>. and <NAME>. and <NAME>.},
    journal = {AIAA Journal},
    volume  = {},
    pages   = {Published Online: 19 Mar 2020},
    url     = {},
    year    = {2020}
}</pre></details>

- \[7\] <NAME>., <NAME>., <NAME>., and <NAME>. [Lift & Learn: Physics-informed machine learning for large-scale nonlinear dynamical systems](https://www.sciencedirect.com/science/article/abs/pii/S0167278919307651). _Physica D: Nonlinear Phenomena_, Volume 406, May 2020, 132401. ([Download](https://kiwi.oden.utexas.edu/papers/lift-learn-scientific-machine-learning-Qian-Willcox.pdf))<details><summary>BibTeX</summary><pre>
@article{QKPW2020_lift_and_learn,
    title   = {Lift \\& Learn: Physics-informed machine learning for large-scale nonlinear dynamical systems.},
    author  = {<NAME>. and <NAME>. and <NAME>. and <NAME>.},
    journal = {Physica {D}: {N}onlinear {P}henomena},
    volume  = {406},
    pages   = {132401},
    url     = {https://doi.org/10.1016/j.physd.2020.132401},
    year    = {2020}
}</pre></details>
<file_sep>/img/README.md
The files here are latex svg images embedded into the markdown-based documentation.
