## 1 Introduction  
Many problems from the applied sciences that involve a dynamic process under study, either to be solved, simulated, or even controlled, are described by models that could contain a considerable amount of variables and/or many inputs and outputs. This is particularly true in the context of fluid dynamics. For instance, equations, such as Navier-Stokes, Burgers', or Kuramoto-Sivashinsky, are nonlinear partial differential equations (PDEs), which are complex and challenging to deal with by means of conventional/direct methods. Thus, these equations are often discretized and linearized, resulting in linear ordinary differential equations (ODEs) described by a huge amount of states (internal variables). The large dimension results in high computational complexity, which leads to challenges in terms of simulation, analysis, and design. Consequently, there is a need to obtain a simplified model of smaller dimensions, which can still be used to provide accurate information for the simulation, analysis, and control of such systems. Model reduction is commonly viewed as a methodology used for reducing the computational complexity of large scale complex models in numerical simulations. The goal is to construct a smaller system with the same structure and similar response characteristics as the original. For an overview of conventional model reduction methods, we refer the reader to the books [1; 2; 3]. For a fairly general overview on system-theoretical model order reduction (MOR) approaches, we refer the reader to [4]. In this work, we assume that the nonlinear systems to be modeled contains quadratic nonlinearities or it can be equivalently transformed to such a structure. This class of systems is of interest since most smooth nonlinear systems can be exactly reformulated (lifted) as quadratic or quadratic-bilinear time-invariant systems (QBTIs), provided that the nonlinearities are analytical. Thus, even if the original model does not involve explicitly (only) quadratic nonlinearities, we could assume this structure (as we will see in the next sections). The main lifting approaches involve artificially introducing new variables into the state-space representation, computing derivatives, and rewriting all relations in the required format. MOR methods specifically tailored for reducing QBTI systems have been proposed in [5; 6; 7; 8]. Such approaches are fairly new, and most of them represent extensions of the techniques proposed for model reduction of linear time-invariant (LTI in short) systems and bilinear time-invariant (BTI) systems [9]. The latter class of systems represent a class of mildly nonlinear dynamical systems, which inherits many properties from LTIs; more specifically speaking, BTIs have many systems and control theoretical quantities that are very closely related or can be easily extended from LTIs (infinite Gramians, system norms, stability properties, etc.). Additionally, model reduction of BTIs is a fairly established sub-field of MOR that has been constantly developed in recent decades. We mention here the research endeavors [10; 11; 12; 13; 14] that preceded the more recent works in [15; 16; 17; 18; 19]. With an ever-increasing availability of measured data in many engineering fields, the need for incorporating measurements in the modeling process has steadily grown over recent decades. The main challenge lies in effectively using the available data in order to construct/learn models that can accurately represent the dynamics of the underlying dynamic process. Sometimes, in order to satisfy accuracy requirements, the fitted/learned models might have large dimensions and, hence, not be suitable for fast numerical simulation. Thus, it is of interest to learn reliable reduced-order surrogate models to be used instead. Here, we assume that such models have particular dynamic structure, i.e., the nonlinearities include bilinear or quadratic-bilinear terms. As seen in the sequel, this is by no means restrictive.