text:
\section{Introduction}

	Large-scale dynamical systems, including electric circuit design, structural mechanics, and microelectromechanical systems, are often described by considerably high-dimensional mathematical models and hence may require formidable computational demands for numerical simulations, control design, and optimization.
	This spurs the development of various model reduction techniques to generate lower-dimensional approximates, or reduced-order models, that allow us to capture the pertinent characteristics of large-scale systems in a computationally manageable fashion. These reduced-order models can further replace their complex counterparts in the control and optimization process, resulting in a shorter design cycle time, see \cite{obinata2012model,antoulas2005approximation,morAntSG01,Gugercin2004Survey} for an overview.

	Reduced-order models are often constructed to closely approximate the input-output mapping of original complex systems while preserving the essential structural features of large-scale systems presents a fundamental problem.
	Many physical systems possess certain properties or structural properties such as stability, passivity, or dissipativity, which are critical in determining the system behaviour and play a crucial role in the control and optimization of the systems. Hence, it is highly desirable that the reduced-order models also preserve the  structural properties of the system and thereby provide an adequate and lean basis for control and optimization, retaining the physical characteristics, see e.g., \cite{wu2018reduced,zheng2015finite,sadamoto2020fast,cheng2021review}.

	Structure-preserving model reduction has attracted extensive attention, and in the existing literature, different methods have been developed for reducing large-scale dynamical systems with specific structural properties. The balanced truncation approach \cite{moore1981principal} naturally retains system stability and minimality through a reduction process, which is extended to retain the positive real and bounded real properties of a linear system
	\cite{Gugercin2004Survey,antoulas2005approximation,reis2010positive}. Further developments in e.g., \cite{phillips2003guaranteed,guiver2013error,kawano2018structure,wu2018reduced,salehi2021passivity,breiten2021error} extend the balanced truncation approach to preserve dissipativity or the port-Hamiltonian structure. The other model reduction methods, including moment matching or Krylov subspace approaches, generally do not guarantee to preserve stability, while the work in \cite{gugercin2008iterative} shows that by combining the features of the singular value decomposition and Krylov-based approximation, the stability and moment matching can be achieved simultaneously. A generalization of this method is then presented in \cite{selga2012stability} by introducing the notions of contractivity and matrix measure. Furthermore, in e.g., \cite{antoulas2005passive,sorensen2005passivity,ionutiu2008passivity,gugercin2012pH,polyuga2010pH,wolf2010passivity,Ionescu2013,hauschild2019pH}, the Krylov subspace approaches have been applied to the mode reduction of passive systems or port-Hamiltonian systems, where reduced-order models are produced with the preservation of the passivity and hence the stability.  Besides, a reduction method for dissipative Hamiltonian systems is presented in \cite{afkham2019structure}, which adopts a symplectic time integrator to conserve the total energy of a Hamiltonian system. Although the aforementioned methods have demonstrated their effectiveness in structure-preserving model reduction, these methods barely guarantee any optimality with respect to a certain criterion. Furthermore, the methods based on the iterative rational Krylov algorithm (IRKA) do not guarantee convergence to local or global optima \cite{flagg2012convergence}. Recently, \cite{Schwerdtner2020structure,breiten2021passivity} present new developments in passivity-preserving model reduction via techniques such as parameter optimization and spectral factorization.

	In contrast to above methods, the works in e.g., \cite{Yu1999approximate,vanDooren2008H2,beattie2007krylov,beattie2009trust,gugercin2008H2,sato2017structure,sato2018pHsys,sato2021book,Jiang2019model,moser2020new} formulate model reduction problem in the context of nonconvex $\mathcal{H}_2$ optimization, which aims to minimize the $\mathcal{H}_2$ error between original and reduced-order models.
	In \cite{beattie2007krylov,beattie2009trust,gugercin2008H2}, the first-order optimality conditions are evaluated on the basis of poles and residues of the transfer matrices of original and reduced-order models, and thereby the rational Krylov algorithm \cite{beattie2007krylov,gugercin2008H2} or the trust region method \cite{beattie2009trust} can be applied to converge to a reduced-order model satisfying first-order optimality conditions with respect to the $\mathcal{H}_2$ error criteria. Nevertheless, structure preservation is not the primary concern of these methods.
	Within the pole-residue framework, \cite{moser2020new} provides an optimal model reduction method that is capable to preserve the port-Hamiltonian structure, while this method is restricted to single-input and single-output systems.

	The other nonconvex optimization framework for model reduction is formulated in e.g., \cite{Yu1999approximate,sato2017structure,sato2018pHsys,sato2021book,Jiang2019model} and the references therein. In this framework, the first-order conditions for (local) $\mathcal{H}_2$ optimality is established on coupled Lyapunov equations, and then gradient-based algorithms are devised to search for reduced-order models on the Stiefel manifold and Riemannian manifold, respectively, leading to structure-preserving reduced-order models. However, these methods in \cite{Yu1999approximate,sato2017structure,sato2018pHsys,sato2021book} require the system state matrix or its symmetric part to be negative definite. The current paper follows this Lyapunov framework but requires only the stability of state matrices and hence can handle more general systems with different structural properties.

	In this paper, we focus on the Petrov-Galerkin projection framework striving to find an optimal reduced-order model that not only admits the $\mathcal{H}_2$ optimality in reduction error but also preserves structural properties of interest. To this end, the $(Q,S,R)$-\textit{dissipativity} of a linear state-space model is considered and characterized by a linear matrix inequality, whose solution is called a structure matrix and incorporated into constructing the project matrix. It can lead to certain degrees of freedom by altering the structure matrix to maintain desired structural properties  including dissipativity, passivity, and bounded realness. Then, the model reduction problem is formulated as a nonconvex optimization on the basis of a noncompact Stiefel manifold. By using the controllability and observability Gramian matrices of the error system, we derive the explicit expression of the gradient of its objective function on the matrix manifold and provide two gradient descent algorithms to find a (local) optimal projection yielding the smallest reduction error.  Unlike the previous works which only preserve one structural feature, this method allows retaining multiple structural properties in one unified framework by simply altering the structure matrix. Furthermore, compared to \cite{Yu1999approximate}, we formulate the optimal projection problem on a
	noncompact Stiefel manifold that enables the algorithm to search the solutions in a larger topological space.

	The paper is organized as follows. In Section~\ref{sec:Preliminaries} we recapitulate some preliminaries on the $\mathcal{H}_2$ norm and introduce the problem setting for the structure-preserving model reduction. The main results are presented in Section~\ref{sec:projection} and Section~\ref{sec:gradient}, where we formulate the model reduction problem in an optimization framework and solve it using gradient-based algorithms. The proposed methods are illustrated numerically in Section~\ref{sec:example}, and some concluding remarks are made in Section~\ref{sec:conclusion}.

	\section{Preliminaries \& Problem Setting}
	\label{sec:Preliminaries}

	Some important notations and fundamental concepts are defined in this section, and the structure-preserving model reduction problem is then introduced.

	\subsection{Notations}

	The symbol $\mathbb{R}$ denotes the set of real numbers. For a matrix $A \in \mathbb{R}^{n \times m}$, the induced Frobenius norm of $A$ is denoted by $\|A\|^2_F = \langle A, A \rangle = \mathrm{Tr}(A^\top A)$. A symbol $\star$ in a square symmetric matrix stands for the symmetric counterpart, namely, we denote a symmetric matrix
	$
	\begin{bmatrix}
		A & B \\ B^\top & C
	\end{bmatrix}
	$
	simply	by
	$
	\begin{bmatrix}
		A & B \\ \star & C
	\end{bmatrix}
	$.

	Let $\mathbb{S}^{n \times r}_*$ denote the set of all $n \times r$ matrices whose columns are linearly independent, i.e.
	\begin{equation} \label{eq:Stiefel}
		\mathbb{S}^{n \times r}_* = \{ V \in \mathbb{R}^{n \times r} \mid \mathrm{Rank} (V) = r  \}.
	\end{equation}
	This is referred to
	as the manifold $\mathbb{S}^{n \times r}_*$, or the \textit{noncompact Stiefel manifold} of full-rank $n \times r$
	matrices. The set $\mathbb{S}^{n \times r}_*$ is an open subset of $\mathbb{R}^{n \times r}$ since its complement $\{X \in \mathbb{R}^{n \times r} : \det(X^\top X) = 0\}$ is closed. Therefore, $\mathbb{S}^{n \times r}_*$ admits a structure
	of an open submanifold of $\mathbb{R}^{n \times r}$. We refer to  \cite{Absil2009optimization} for more details.

	\subsection{Gramians and $\mathcal{H}_2$ Norm of Linear Stable Systems}

	We recapitulate the concept of the $\mathcal{H}_2$ norm of a stable linear time-invariant system with the following state-space representation.
	\begin{equation}\label{sys:orig}
		\mathbf{\Sigma}: \
		\begin{cases}
			\dot{x}(t) = A x(t) + B u(t), \\
			y(t) = C x(t) + D u(t),
		\end{cases}
	\end{equation}
	where $x(t) \in \mathbb{R}^n$, $u(t) \in \mathbb{R}^p$, and $y(t) \in \mathbb{R}^{q}$ are the vectors of the states, inputs, and outputs of the system $\mathbf{\Sigma}$. $A \in \mathbb{R}^{n \times n}$, $B \in \mathbb{R}^{n \times p}$, $C \in \mathbb{R}^{q \times n}$,  and $D \in \mathbb{R}^{q \times p}$ are constant matrices. If $A$ is Hurwitz, the system $\mathbf{\Sigma}$ is (asymptotically) stable. Then,
	the \textit{controllability and observability Gramians} are properly defined as
	\begin{equation*}
		\label{defn:Gramians}
		\mathcal{P} = \int_{0}^{\infty}
		e^{A t} BB^\top e^{A^\top t} \mathrm{d}{t},
		\quad
		\mathcal{Q} = \int_{0}^{\infty}
		e^{A^\top t} C^\top C e^{A t} \mathrm{d}{t},
	\end{equation*}
	respectively \cite{antoulas2005approximation,gugercin2008H2}. Furthermore, the two positive semidefinite Gramian matrices are solved as unique solutions of the following Lyapunov equations:
	\begin{align} \label{eq:gramins}
		A  \mathcal{P} + \mathcal{P} A^\top + B  B^\top =0,  \quad A^\top  \mathcal{Q} +  \mathcal{Q} A + C^\top C =0.
	\end{align}

	Let the system $\mathbf{\Sigma}$ in \eqref{sys:orig} be strictly causal, i.e. $D = 0$, and  $G(s):=C(sI-A)^{-1} B$ be the transfer matrix of $\mathbf{\Sigma}$. Then, the $\mathcal{H}_2$ norm of $\mathbf{\Sigma}$ is defined and computed as
	\begin{align*}
		\|G(s)\|_{\mathcal{H}_2}^2
		& = \int_{0}^{\infty}  \mathrm{Tr} \left(  C e^{At} B B^\top e^{A^\top t} C^\top \right)  \mathrm{d}{t}
		\\
		& = \mathrm{Tr}(C \mathcal{P}  C^\top) = \mathrm{Tr}(B^\top \mathcal{Q}  B).
	\end{align*}
	In this paper, the $\mathcal{H}_2$ norm is chosen to be a metric for evaluating the approximation accuracy of a reduced-order model such that an $\mathcal{H}_2$ optimal model reduction problem is formulated, where the objective function and its gradient can be formulated on the matrix manifold \eqref{eq:Stiefel}.

	\subsection{Structure-Preserving Model Reduction}
	\label{sec:problem}
	The overall aim of this paper is to develop a general model reduction methodology that is applicable to a wide range of systems with specific structural properties, including passivity (or positive realness) and finite-gain $L_2$ stability (or bounded realness). To this end, we consider the concept of the $(Q, S, R)$-dissipativity, and the aforementioned structural properties can be considered as special cases of the $(Q, S, R)$-dissipativity.
	\begin{definition}[$(Q, S, R)$-Dissipativity]\cite{willems2007dissipative,brogliato2007dissipative,xie1998robust,xia2016passivity}
		For given matrices $Q = Q^\top$, $R = R^\top$
		and $S$, a system $\mathbf{\Sigma}$ in \eqref{sys:orig} is called $(Q,S,R)$-
		dissipative if
		\begin{equation}
			\int_{0}^{\tau} y(t)^\top Q y(t) + u(t)^\top R u(t) + 2 u(t)^\top S y(t) \mathrm{d}{t} \geq 0
		\end{equation}
		holds for all admissible $u(t)$ and all $\tau \geq 0$.
	\end{definition}
	It is further shown in \cite{willems2007dissipative,xia2016passivity,kottenstette2014relationships} that the $(Q, S, R)$-dissipativity of a linear state-space model \eqref{sys:orig} can be characterized by a linear matrix inequality.
	\begin{lemma} \label{lem:QSRdisip}
		A linear system \eqref{sys:orig} is $(Q, S, R)$-dissipative if and only if there exists
		a matrix $K > 0$ such that
		\begin{equation} \label{eq:LMI-QSR}
			\Gamma = \begin{bmatrix}
				A^\top K + K A & KB - C^\top S & - C^\top Q^{\frac{1}{2}} \\
				\star & - R - D^\top S - S^\top D & -D^\top Q^{\frac{1}{2}}  \\
				\star & \star & - I
			\end{bmatrix} \leq 0.
		\end{equation}
		If $\Gamma$ is strictly negative definite, then the system \eqref{sys:orig} is asymptotically stable and $(Q, S, R)$-dissipative.

	\end{lemma}

	This work focuses on the notion of  $(Q, S, R)$-dissipativity as it is a rather general property that can be reduced, with particular weighting matrices $Q$, $S$ and $R$, to passivity and finite-gain $L_2$ stability.
	\begin{enumerate}

		\item Let $Q = 0$, $S = I$, $R = 0$. Then, the $(Q, S, R)$-dissipativity reduces to the \textit{passivity} (or \textit{positive realness}) as $
		\int_{0}^{\tau} u(t)^\top y(t) \mathrm{d}{t} \geq 0
		$
		for all admissible $u(t)$ and all $\tau \geq 0$. The inequality \eqref{eq:LMI-QSR} is simplified as
		\begin{equation} \label{eq:Kpassiv}
			\begin{bmatrix}
				A^\top K + K A &  KB - C^\top
				\\
				\star   &  {-D^\top - D}
			\end{bmatrix} \leq 0,
		\end{equation}

		\item Let $Q = -\gamma^{-1} I$, $S = 0$, $R = \gamma I$, then
		the $(Q, S, R)$-dissipativity becomes the \textit{finite-gain $L_2$ stability} (or \textit{bounded realness}) as
		$$
		\int_{0}^{\tau} y(t)^\top y(t) \mathrm{d}{t} \leq \gamma^2 \int_{0}^{\tau} u(t)^\top u(t) \mathrm{d}{t}
		$$ for all admissible $u(t)$ and all $\tau \geq 0$, and meanwhile \eqref{eq:LMI-QSR} is simplified by {  the Schur complement lemma} as
		\begin{equation} \label{eq:KL2}

			\begin{bmatrix}
				A^\top K + K A &  KB
				\\
				\star   &  - \gamma I
			\end{bmatrix} {  +} \frac{1}{\gamma}
			\begin{bmatrix}
				C^\top
				\\
				D^\top
			\end{bmatrix}
			\begin{bmatrix}
				C & D
			\end{bmatrix}
			\leq 0,
		\end{equation}

	\end{enumerate}

	\begin{remark}
		We emphasize that the inequality in \eqref{eq:LMI-QSR} only provides a characterization for the $(Q,S,R)$-dissipativity, which is a general concept including passivity and bounded realness. However, it is not necessary that the LMI in \eqref{eq:LMI-QSR} has always to be solved numerically to obtain a solution of $K$. For instance, to reduce port-Hamiltonian systems as in e.g., \cite{gugercin2012pH,polyuga2010pH,wolf2010passivity,Ionescu2013,sato2018pHsys,Jiang2019model}, we can simply use positive definite Hamiltonian matrices as $K$. To preserve the stability of {  an} observable system, $K$ can be chosen as the observability Gramian of the system as in \cite{gugercin2008iterative}.
	\end{remark}

	In this paper, we consider a full-order model $\mathbf{\Sigma}$ in \eqref{sys:orig}, which is asymptotically stable and $(Q, S, R)$-dissipative, i.e. $A$, $B$, and $C$ satisfy the inequality \eqref{eq:LMI-QSR}.
	Then, the problem of structure-preserving model reduction in this paper is formulated as follows.
	\begin{problem}
		Given a fixed positive integer $r$ ($ < n$), we aim  to construct a low-dimensional model
		\begin{equation}\label{sys:red0}
			\mathbf{\hat{\Sigma}}: \
			\begin{cases}
				\dot{\hat{x}}(t)  = \hat{A} \hat{x}(t) +  \hat{B} u(t), \\
				\hat{y}(t) = \hat{C}   \hat{x}(t) + \hat{D} u(t),
			\end{cases}
		\end{equation}
		where $\hat{A} \in \mathbb{R}^{r \times r}$, $\hat{B}\in \mathbb{R}^{r \times p}$, $\hat{C} \in \mathbb{R}^{q \times r}$, and $\hat{D} \in \mathbb{R}^{q \times p}$ are reduced matrices such that
		\begin{enumerate}
			\item the reduction procedure guarantees the retention of the stability and the $(Q, S, R)$-dissipativity in the reduced-order model $\mathbf{\hat{\Sigma}}$, and

			\item the reduced-order model $\mathbf{\hat{\Sigma}}$ best approximates the original model $\mathbf{\Sigma}$ in the sense that the $\mathcal{H}_2$ norm of the error between $\mathbf{\Sigma}$ and $\mathbf{\hat{\Sigma}}$ is minimized.
		\end{enumerate}
	\end{problem}
	The objective (1) is {  to} guarantee structural properties of interest retained in the reduced-order model. Under this constraint, the objective in (2) aims to minimize the approximation error in the $\mathcal{H}_2$ norm.

	In this work, a unified framework will be presented to solve the structure-preserving model reduction problem, where the reduced-order model is built based on the Petrov-Galerkin projection. The model reduction problem thereby is reformulated as a nonconvex optimization problem, which is solved efficiently by a gradient descent method.

	\section{Projection-Based Model Reduction With Structure Preservation}
	\label{sec:projection}

	Consider a full-order model in \eqref{sys:orig} with any  structural property discussed in Section~\ref{sec:problem}. To preserve the desired structure in the reduced-order model \eqref{sys:red0}, a novel Petrov-Galerkin projection is adopted to generate a reduced-order model.
	Let $X = X^\top > 0$, $V \in \mathbb{S}^{n \times r}_*$, and denote \begin{equation}\label{eq:pinv}
		V^\dagger : = (V^\top X V)^{-1} V^\top X,
	\end{equation}
	in which $V^\top X V$ is invertible since $V$ has full column rank, and $X \in \mathbb{R}^{n \times n}$ is a positive definite matrix. $V^\dagger$ is called a \textit{reflexive generalized inverse} of $V$. The projection matrix is thereby defined as $\Pi = V V^\dagger $, by which a reduced-order model is constructed as follows.
	\begin{equation}\label{sys:red}
		\mathbf{\hat{\Sigma}}: \
		\begin{cases}
			\dot{\hat{x}}(t)  = \underbrace{V^\dagger A  V}_{\hat{A}} \hat{x}(t) + \underbrace{V^\dagger B}_{\hat{B}} u(t), \\
			\hat{y}(t) = \underbrace{C V}_{\hat C} \hat{x}(t) + D u(t),
		\end{cases}
	\end{equation}
	where $\hat{x} \in \mathbb{R}^{r}$ is the reduced state vector. The matrix $X$ is called a \textit{structure matrix}, which allows us to preserve certain desired structural properties. The following result shows that this particular construction of a reduced-order model can lead to structure preservation.
	\begin{lemma}[Structure Preservation]
		Consider the original system $\mathbf{\Sigma}$ \eqref{sys:orig} and its reduced-order model $\mathbf{\hat{\Sigma}}$ defined in \eqref{sys:red}. If $\mathbf{\Sigma}$ is $(Q, S, R)$-dissipative and $X = K  > 0$ in \eqref{sys:red} is a solution of \eqref{eq:LMI-QSR}, then $\mathbf{\hat{\Sigma}}$ in \eqref{sys:red} is also $(Q, S, R)$-dissipative.
	\end{lemma}
	\begin{proof}
		To verify the $(Q, S, R)$-dissipativity of the reduced-order model $\mathbf{\hat{\Sigma}}$ in \eqref{sys:red}, we let $\hat{K} = V^\top X V  = V^\top K V$. Then, the following equations hold.
		\begin{align*}
			\hat{A}^\top \hat{K} + \hat{K} \hat{A}
			= & V^\top A^\top X V + V^\top X A V \\
			= & V^\top (A^\top K + K A) V
			\\
			\hat{K} \hat{B} - \hat{C}^\top S = &  V^\top (KB - C^\top S),
			\\ \hat{C}^\top Q^{\frac{1}{2}} = &  V^\top C^\top Q^{\frac{1}{2}}.
		\end{align*}
		Based on that, we have
		\begin{align*}
			& \begin{bmatrix}
				\hat{A}^\top \hat{K} + \hat{K} \hat{A} & \hat{K}\hat{B} - \hat{C}^\top S & - \hat{C}^\top Q^{\frac{1}{2}} \\
				\star & - R - D^\top S - S^\top D & -D^\top Q^{\frac{1}{2}}  \\
				\star & \star & - I
			\end{bmatrix} \\
			& =  \begin{bmatrix}
				V^\top & & \\ & I & \\ & & I
			\end{bmatrix}
			\Gamma
			\begin{bmatrix}
				V & & \\ & I & \\ & & I
			\end{bmatrix} \leq 0,
		\end{align*}
		where $\Gamma$ defined in \eqref{eq:LMI-QSR} is negative semidefinite due to the $(Q, S, R)$-dissipativity of the original model $\mathbf{\Sigma}$.
		Therefore, it follows from Lemma~\ref{lem:QSRdisip} that the reduced-order model $\mathbf{\hat{\Sigma}}$ remains $(Q, S, R)$-dissipative.
	\end{proof}

	Different structural properties can be preserved in the reduced-order model  \eqref{sys:red} by altering the structure matrix $X$. More preciously, let $X$ be the solution of \eqref{eq:Kpassiv}, \eqref{eq:KL2} or the {  inequality} $AX + X A^\top < 0$, we can retain passivity, finite-gain $L_2$ stability, or stability, respectively. Therefore, a unified framework is presented here to address the model reduction problems with different structural constraints.   Particularly, this work provides a novel method for the passivity-preserving model reduction problems for port-Hamiltonian systems, which have been extensively studied in the literature, see e.g.,
	\cite{antoulas2005passive,ionutiu2008passivity,gugercin2012pH,polyuga2010pH,wolf2010passivity,Ionescu2013,hauschild2019pH,phillips2003guaranteed,salehi2021passivity,sato2018pHsys,moser2020new}.
	In our approach, the structure matrix $X$ in \eqref{eq:pinv} can be chosen as the energy matrix $K$ of the system i.e. $\frac{1}{2}x(t)^\top K x(t)$ represents the total energy (Hamiltonian) of the system. Consequently, the reduced-order model has the reduced energy matrix as $V^\top K V$.

	\begin{remark}[Multiple Structural
		Properties]
		The proposed model reduction framework can also be applied to preserve multiple structural properties, while $X$ in \eqref{eq:pinv}  is required to be a common solution of the corresponding matrix inequalities or equations. For instance, if there exists $K$ satisfying \eqref{eq:Kpassiv} and \eqref{eq:KL2} simultaneously, then we let $X = K$ in \eqref{eq:pinv} that preserves the passivity and finite-gain $L_2$ stability in the reduced-order model \eqref{sys:red}.
	\end{remark}

	\begin{remark}

		Instead of the $(Q, S, R)$-dissipativity, we may consider retaining controllability or observability of a given system $\mathbf{\Sigma}$. In this regard, we let the structure matrix $X$ be the solution of
		\begin{equation*}
			A^\top X + X A + X BB^\top X = 0,
		\end{equation*}
		or
		\begin{equation*}
			A^\top X + X A + C^\top C = 0,
		\end{equation*}
		using which to construct a reduced-order model \eqref{sys:red} will allow the preservation of controllability or observability. Furthermore, if $X$ is chosen as the observability Gramian of the system $\mathbf{\Sigma}$, then the stability is simultaneously preserved in the reduced-order model \eqref{sys:red} \cite{gugercin2008iterative}.
	\end{remark}

	The reduced-order model in \eqref{sys:red} guarantees the preservation of a desirable structural property by selecting a specific matrix $X$. Next, we aim to minimize the approximation error between the original and reduced models by tuning $V$. To this end, we define the error system as
	\[
	G_e(s) = C(sI_n - A)^{-1} B - C V (sI_r - V^\dagger A V)^{-1} V^\dagger B,
	\]
	where the feedthrough gain $D$ in the two systems has been canceled out. The system $G_e(s)$ has the following state-space representation
	\begin{align} \label{sys:err}
		\begin{cases}
			\dot{z}(t) = {A}_e z(t) + {B}_eu(t )\\
			y_e (t) = {C}_e z(t)
		\end{cases}
	\end{align}
	where $z(t) = [x(t)^\top \ \hat{x}(t)^\top]^\top \in \mathbb{R}^{n+r}$, and
	\begin{align*}
		{A}_e = \begin{bmatrix}
			A & 0 \\ 0 & V^\dagger A  V
		\end{bmatrix}, \
		{B}_e = \begin{bmatrix}
			B  \\  V^\dagger B
		\end{bmatrix}, \
		{C}_e = \begin{bmatrix}
			C  & - CV
		\end{bmatrix}.
	\end{align*}
	To guarantee the stability of the error system, we consider the original system $\mathbf{\Sigma}$ asymptotically stable, i.e. $A$ in \eqref{sys:orig} is Hurwitz. Thereby, let the structure matrix $X$ additionally be a solution of $\Gamma < 0$, instead of $\Gamma \leq 0$ in  \eqref{eq:LMI-QSR}. Note that $\Gamma < 0$ implies
	\begin{equation*}
		A^\top X + X A  < 0,
	\end{equation*}
	which will impose the asymptotic stability of the reduced-order model $\mathbf{\hat{\Sigma}}$.  Then, the square $\mathcal{H}_2$-norm of the approximation error is computed as
	\begin{equation}
		\| G_e(s) \|_{\mathcal{H}_2}^2 =
		\mathrm{Tr}(C_e \mathcal{P} C_e^\top) = \mathrm{Tr}(B_e^\top \mathcal{Q} B_e),
	\end{equation}
	where $\mathcal{P}$ and $\mathcal{Q}$ are the controllability and observability Gramians of the error system \eqref{sys:err}. Algebraically, $\mathcal{P}$ and $\mathcal{Q}$ are solved as the unique solutions of the following Lyapunov equations:
	\begin{subequations}\label{eq_Gramians_err}
		\begin{align}
			\label{W_Lyap}
			A_e \mathcal{P}+\mathcal{P} A_e^\top+ B_e B_e^\top &=0, \\
			\label{M_Lyap}
			A_e^\top\mathcal{Q}+\mathcal{Q} A_e + C_e^\top C_e &=0.
		\end{align}
	\end{subequations}

	Observe that given a structure matrix $X$, the reduced-order model \eqref{sys:red} is parameterized in the matrix $V$.
	We formulate an optimization problem in terms of minimizing the $\mathcal{H}_2$ norm of the error system $G_e(s)$:
	\begin{align}
		\label{pro:optimization}
		& \min_{V,\mathcal{Q}} \; \mathrm{Tr}(B_e^\top \mathcal{Q} B_e)  \\
		& \;\; \text{s.t.:}  \;\; V^\top V > 0, \quad  V \in \mathbb{R}^{n \times r},
		\nonumber
		\\ & \;\; \quad \quad A_e^\top \mathcal{Q} + \mathcal{Q} A_e + C_e^\top C_e =0.
		\nonumber
	\end{align}
	This optimization problem is highly {nonconvex} and has a continuous variable $\mathcal{Q} \in {\mathbb R}^{(n+r) \times (n+r)}$ and a full column rank matrix $V \in \mathbb{R}_*^{n \times r}$, where $\mathcal{Q} =\mathcal{Q}(V)$ is an implicit function of $V$ due to the unique solution of the Lyapunov equation \eqref{M_Lyap}. Hence, for any $V \in \mathbb{R}_*^{n \times r}$, the partial minimization in $\mathcal{Q}$ leads to the optimization problem
	\begin{align}
		\label{eq:optimization1}
		\min_{V \in \mathbb{S}^{n\times r}_*} \;  J(V)  \; : & =  \mathrm{Tr} \left( B_e(V)^\top  \mathcal{Q}(V)  B_e(V) \right)  \\
		& = \mathrm{Tr} \left( \mathcal{Q}(V)  B_e(V) B_e(V)^\top  \right),\nonumber
	\end{align}
	where the set $\mathbb{S}^{n\times r}_*$ is defined in \eqref{eq:Stiefel}.

	\section{Gradient-Based Methods for Optimal Projection}
	\label{sec:gradient}
	This section presents gradient-based methods for solving the model reduction problem as a nonconvex optimization problem. The explicit expression for the gradient of the objective function $J(V)$ is derived, and furthermore, gradient descent algorithms are proposed to minimize $J(V)$ over both noncompact and compact Stiefel manifolds.

	\subsection{Gradient Analysis}

	To solve the optimization problem \eqref{eq:optimization1}, we analyze the key properties of the objective function $J(V)$ such as smoothness and differentiability, based on which we can apply gradient-based methods for solving \eqref{eq:optimization1}  over the manifold $\mathbb{S}^{n \times r}_*$. In particular, we show in the next theorem that $J$ is a differentiable function and derive an explicit expression for its gradient.
	Before proceeding, we partition the controllability and observability Gramians $\mathcal{P}$ and $\mathcal{Q}$   according to the block structure of the error system as:
	\begin{equation} \label{eq:blockMW}
		\mathcal{P} = \begin{bmatrix}
			\mathcal{P}_{11} & \mathcal{P}_{12} \\ \mathcal{P}_{12}^\top & \mathcal{P}_{22}
		\end{bmatrix}, \quad
		\mathcal{Q} = \begin{bmatrix}
			\mathcal{Q}_{11} & \mathcal{Q}_{12} \\ \mathcal{Q}_{12}^\top & \mathcal{Q}_{22}
		\end{bmatrix},
	\end{equation}
	where $P_{11}$ and $Q_{11}$ are Gramians of the high-order system \eqref{sys:orig}, and $P_{22}$ and $Q_{22}$ are the Gramians of the reduced-order model \eqref{sys:red}. Here, the submatrices $\mathcal{P}_{22}$, $\mathcal{Q}_{22}$, $\mathcal{P}_{12}$, and $\mathcal{Q}_{12}$ are characterized by the following Sylvester equations:
	\begin{subequations} \label{eq:Sylv}
		\begin{align}
			\hat{A}	\mathcal{P}_{22} + \mathcal{P}_{22}\hat{A}^\top + \hat{B} \hat{B}^\top = 0,
			\\
			\hat{A}^\top	\mathcal{Q}_{22} + \mathcal{Q}_{22}\hat{A}  + \hat{C}^\top \hat{C}  = 0,
			\\
			{A}	\mathcal{P}_{12} + \mathcal{P}_{12}  \hat{A}^\top +  {B} \hat{B}^\top = 0,
			\label{eq:Sylv_c}
			\\
			{A}^\top	\mathcal{Q}_{12} + \mathcal{Q}_{12}\hat{A}  -  {C}^\top \hat{C}  = 0,
			\label{eq:Sylv_d}
		\end{align}
	\end{subequations}
	where $\hat{A}$, $\hat{B}$, and $\hat{C}$ are computed as in \eqref{sys:red}. Then, the following theorem holds.
	\begin{theorem}[Gradient on Matrix Manifold]
		\label{thm:gradJ}
		The objective function $J(V)$ in  \eqref{eq:optimization1} is differentiable, and the gradient has the following expression:
		\begin{align}  \label{eq:nablaJ}
			\nabla J(V)
			= & 2 X (I - V V^\dagger)  \left[ \mathcal{F}_A (V) + \mathcal{F}_B (V)^\top \right]   (V^\top X V)^{-1}
			\nonumber \\
			& - 2(V^\dagger)^\top  \left[ \mathcal{F}_A (V)^\top + \mathcal{F}_B (V) \right]  (V^\dagger)^\top
			\nonumber
			\\& 	 +  2 A^\top (V^\dagger)^\top \Sigma_{22}^\top   + 2 \mathcal{F}_C (V),
		\end{align}
		where $\Sigma_{22}: = \mathcal{P}_{12}^\top\mathcal{Q}_{12} +  \mathcal{P}_{22}\mathcal{Q}_{22}$, and
		\begin{subequations}\label{eq:FAFBFC}
			\begin{align}
				\mathcal{F}_A (V) : & = A V    \Sigma_{22}, \label{eq:F_A}\\
				\mathcal{F}_B (V) : & = \mathcal{Q}_{12}^\top BB^\top    +  \mathcal{Q}_{22} V^\dagger BB^\top, \label{eq:F_B}\\
				\mathcal{F}_C (V) : & =  C^\top C V \mathcal{P}_{22} - C^\top C\mathcal{P}_{12}.
			\end{align}
		\end{subequations}
	\end{theorem}

	%

	\begin{proof}
		Observe that
		$
		V^\top X = (V^\top X V) V^\dagger.
		$
		We take the differential on both sides, leading to
		\begin{align*}
			\mathrm{d} V^\top X = &\mathrm{d}((V^\top X V) V^\dagger)
			\\
			=&  \mathrm{d} (V^\top X V) V^\dagger + (V^\top X V)  \mathrm{d}  V^\dagger
			\\
			= & (\mathrm{d} V^\top X V + V^\top X \mathrm{d} V) V^\dagger + (V^\top X V)  \mathrm{d} (V^\dagger),
		\end{align*}
		from which, the differential $ \mathrm{d} (V^\dagger)$ is solved as
		\begin{equation}\label{eq:differential0}
			\mathrm{d}(V^\dagger) = (V^\top X V)^{-1}  \mathrm{d} V^\top X (I - V V^\dagger)  - V^\dagger \mathrm{d} V V^\dagger.
		\end{equation}

		To compute the gradient $\nabla J(V)$, we write the derivative  $J'(V) \mathrm{d}{V}$ for some $\mathrm{d}{V} \in  \mathbb R^{n \times r }$ in a gradient form in terms of trace as follows.
		\begin{align*}
			J'(V)  \mathrm{d}{V}  = \mathrm{Tr} \left( \nabla J(V)^\top \mathrm{d}{V}  \right).
		\end{align*}
		Considering \eqref{eq:optimization1}, $J'(V)  \mathrm{d}{V}$ can be specified as
		\begin{align} \label{diff_JV}
			J'(V)  \mathrm{d}{V}
			& = \mathrm{Tr} \left(  \mathcal{Q}(V) \mathrm{d}{\mathcal{B}}(V) + \mathrm{d}  \mathcal{Q}(V)  {\mathcal{B}}(V)      \right)
			\nonumber\\
			& =  \mathrm{Tr}\left(\mathcal{Q}(V) \mathrm{d}{\mathcal{B}}(V)   \right)
			+ \mathrm{Tr} \left(   \mathrm{d}  \mathcal{Q}(V)  {\mathcal{B}}(V) \right) ,
		\end{align}
		where ${\mathcal{B}}(V)$ is denoted for notational simplicity as
		\begin{align}\label{eq:BeBe}
			\mathcal{B}(V) = B_e B_e^\top = \begin{bmatrix}
				B B^\top & B B^\top (V^\dagger)^\top \\
				V^\dagger B B^\top & V^\dagger B B^\top (V^\dagger)^\top
			\end{bmatrix}.
		\end{align}
		In the following, the two terms in \eqref{diff_JV} are computed separately.

		First, we analyze the first term in \eqref{diff_JV}. From \eqref{eq:BeBe}, it follows that
		\begin{align*}
			\mathrm{d} \mathcal{B} (V)
			= \begin{bmatrix}
				0 & BB^\top \mathrm{d} (V^\dagger)^\top
				\\ \mathrm{d} (V^\dagger) BB^\top&  \mathrm{d} (\hat{B} \hat{B}^\dagger)
			\end{bmatrix},
		\end{align*}
		where $\mathrm{d} (\hat{B} \hat{B}^\dagger) = \mathrm{d} (V^\dagger) BB^\top (V^\dagger)^\top +  V^\dagger BB^\top \mathrm{d} (V^\dagger)^\top$.
		Thus, the first term in \eqref{diff_JV} can be specified as
		\begin{align}\label{eq:trace_MdiffB0}
			& \mathrm{Tr} \left( \mathcal{Q}(V) \mathrm{d} \mathcal{B} (V) \right)
			\\ = & \mathrm{Tr} \left(\mathcal{Q}_{12} \mathrm{d}(V^\dagger) BB^\top  + \mathcal{Q}_{12}^\top BB^\top \mathrm{d}(V^\dagger)^\top \right.
			\nonumber\\ & \left. \quad+  \mathcal{Q}_{22} \mathrm{d}(V^\dagger) BB^\top (V^\dagger)^\top + \mathcal{Q}_{22} V^\dagger BB^\top \mathrm{d}(V^\dagger)^\top \right).
			\nonumber
		\end{align}
		Using the standard properties of trace, we have
		\begin{align*}
			\mathrm{Tr} \left(\mathcal{Q}_{12} \mathrm{d}(V^\dagger) BB^\top\right) & = \mathrm{Tr} \left(BB^\top \mathrm{d}(V^\dagger)^\top \mathcal{Q}_{12}^\top  \right) \\
			& = \mathrm{Tr} \left(\mathcal{Q}_{12}^\top BB^\top \mathrm{d}(V^\dagger)^\top   \right),
		\end{align*}
		and
		\begin{align*}
			\mathrm{Tr} \left(\mathcal{Q}_{22} \mathrm{d}(V^\dagger) BB^\top (V^\dagger)^\top \right) & = \mathrm{Tr} \left( V^\dagger BB^\top   \mathrm{d}(V^\dagger)^\top \mathcal{Q}_{22} \right) \\
			& = \mathrm{Tr} \left( \mathcal{Q}_{22} V^\dagger BB^\top   \mathrm{d}(V^\dagger)^\top  \right).
		\end{align*}
		Consequently, \eqref{eq:trace_MdiffB0} is simplified as
		\begin{align*}
			\mathrm{Tr} \left( \mathcal{Q}(V) \mathrm{d} \mathcal{B} (V) \right)
			=  2\mathrm{Tr}   \left( \mathcal{F}_B(V) \mathrm{d}(V^\dagger)^\top\right).
			\nonumber
		\end{align*}
		with $\mathcal{F}_B(V)$ defined in \eqref{eq:F_B}.
		Substituting the differential of $\mathrm{d}(V^\dagger)$ in \eqref{eq:differential0} into the above equation then yields
		\begin{align} \label{trace_MdiffB}
			& \mathrm{Tr} \left( \mathcal{Q}(V) \mathrm{d} \mathcal{B} (V) \right)
			\\ = & 2 \mathrm{Tr} \left[ \mathcal{F}_B(V)(I - (VV^\dagger)^\top) X \mathrm{d} V (V^\top X V)^{-1} \right. \nonumber
			\\ & \left. \quad - \mathcal{F}_B(V) (V^\dagger)^\top \mathrm{d}(V)^\top (V^\dagger)^\top
			\right]
			\nonumber\\
			= & 2\mathrm{Tr} \left[ (V^\top X V)^{-1} \mathcal{F}_B(V)(I - (VV^\dagger)^\top) X \mathrm{d} V  \right. \nonumber
			\\ & \left. \quad - V^\dagger  \mathcal{F}_B(V)^\top  V^\dagger  \mathrm{d} V
			\right].
			\nonumber
		\end{align}

		Next, we derive an explicit expression for the second term in \eqref{diff_JV}, i.e. $\mathrm{Tr} \left(   \mathrm{d}  \mathcal{Q}(V)  {\mathcal{B}}(V) \right)$. To this end, we define the function:
		\begin{equation} \label{eq:PhiVM}
			\Phi(V,\mathcal{Q}) = A_e^\top \mathcal{Q}  + \mathcal{Q}  A_e +  \mathcal{C}(V) = 0,
		\end{equation}
		whose partial derivatives with respect to $V$ and $\mathcal{Q}$ are given as
		\begin{equation*}
			\Phi_V'(V, \mathcal{Q}) \mathrm{d}{V} = (\mathrm{d} A_e(V))^\top {\mathcal{Q}} + {\mathcal{Q}} \mathrm{d} A_e(V) +  \mathrm{d} \mathcal{C}(V),
		\end{equation*}
		and
		\begin{equation*}
			\Phi_\mathcal{Q}' (V, \mathcal{Q}) \mathrm{d}{\mathcal{Q}} = (A_e(V))^\top \mathrm{d}{\mathcal{Q}} + \mathrm{d}{\mathcal{Q}} A_e(V),
		\end{equation*}
		where
		\begin{align} \label{eq:CeCe}
			\mathcal{C}(V) = C_e^\top C_e = \begin{bmatrix}
				C^\top C & - C^\top C V \\
				-V^\top  C^\top C & V^\top  C^\top CV
			\end{bmatrix}.
		\end{align}
		Since $\Phi(V,\mathcal{Q})$ a multivariable function of both $V$ and $\mathcal{Q}$ and $\Phi(V,\mathcal{Q}) = 0$, applying the Implicit Function Theorem yields $ \Phi_V'(V, \mathcal{Q}) \mathrm{d}{V} +  \Phi_\mathcal{Q}' (V, \mathcal{Q}) \mathrm{d}{\mathcal{Q}} = 0$, i.e.
		\begin{multline*}
			(\mathrm{d} A_e)^\top  \mathcal{Q}(V) + \mathcal{Q}(V) \mathrm{d} A_e(V) + \mathrm{d} \mathcal{C}(V)       \\
			+ (A_e(V))^\top \mathrm{d}{\mathcal{Q}} (V) + \mathrm{d}{\mathcal{Q}}(V) A_e(V) = 0.
		\end{multline*}
		Multiplying $\mathcal{P}(V)$ from the left then gives
		\begin{multline} \label{diff_Lyap}
			\mathcal{P}(\mathrm{d} A_e)^\top  \mathcal{Q}  + \mathcal{P} \mathcal{Q}  \mathrm{d} A_e  + \mathcal{P} \mathrm{d} \mathcal{C}    \\
			+ \mathcal{P} A_e^\top \mathrm{d}{\mathcal{Q}}   + \mathcal{P} \mathrm{d}{\mathcal{Q}}  A_e  = 0,
		\end{multline}
		where we suppress the dependence on $V$ for ease of notation. Furthermore, notice that the controllability Gramian $\mathcal{P}(V)$ is the unique solution of the Lyapunov equation \eqref{W_Lyap}, from which we obtain
		\begin{equation}
			\label{W_Lyap-P}
			A_e  \mathcal{P} \mathrm{d}{\mathcal{Q}}+ \mathcal{P} A_e^\top \mathrm{d}{\mathcal{Q}}   = -\mathcal{B} \mathrm{d}{\mathcal{Q}}.
		\end{equation}
		Subtracting \eqref{W_Lyap-P} from \eqref{diff_Lyap} leads to
		\begin{multline*}
			\mathcal{P} A_e^\top \mathrm{d}{\mathcal{Q}} + \mathcal{P} \mathrm{d}{\mathcal{Q}} A_e
			- A_e \mathcal{P} \mathrm{d}{\mathcal{Q}} - \mathcal{P} A_e^\top \mathrm{d}{\mathcal{Q}} \nonumber
			\\ +
			\mathcal{P}(\mathrm{d} A_e)^\top  \mathcal{Q} + \mathcal{P} \mathcal{Q} \mathrm{d} A_e
			- \mathcal{B} \mathrm{d}{\mathcal{Q}}
			+ \mathcal{P} \mathrm{d} \mathcal{C} = 0.
		\end{multline*}
		Then, we take the trace in the above equation that cancels the first four terms and yields an expression for the first term in \eqref{diff_JV} as
		\begin{equation}
			\label{trace_diffMB}
			\mathrm{Tr} \left( \mathrm{d} \mathcal{Q}  \mathcal{B} \right)
			=
			2\mathrm{Tr} \left( \mathcal{P}  \mathcal{Q}  \mathrm{d} A_e \right)
			+ \mathrm{Tr} \left(  \mathcal{P}  \mathrm{d} \mathcal{C}
			\right),
		\end{equation}
		where
		\begin{align*}
			\mathrm{d} A_e(V) =   \begin{bmatrix}
				0 & 0 \\
				0 & \mathrm{d} \hat{A}(V)
			\end{bmatrix},
		\end{align*}
		with $\mathrm{d} \hat{A}(V) = \mathrm{d} ( V^\dagger A V) = \mathrm{d}(V^\dagger) A V +  V^\dagger A \mathrm{d} V$.

		Applying \eqref{eq:differential0}, we therefore obtain
		\begin{align} \label{eq:trace_diffMB1}
			& \mathrm{Tr} \left( \mathcal{P} \mathcal{Q} \mathrm{d} A_e  \right)
			\\
			= & \mathrm{Tr} \left[ \Sigma_{22} \mathrm{d}(V^\dagger) A V + \Sigma_{22}  V^\dagger  A \mathrm{d} V  \right]\nonumber\\
			= &
			\mathrm{Tr} \left[ \Sigma_{22} (V^\top X V)^{-1}  \mathrm{d} V^\top X (I - V V^\dagger) A V   \right. \nonumber\\
			& \quad \left.  - \Sigma_{22} V^\dagger \mathrm{d} V  V^\dagger  AV  + \Sigma_{22} V^\dagger A \mathrm{d} V   \right]
			\nonumber\\
			= &
			\mathrm{Tr} \left[ (V^\top X V)^{-1} \mathcal{F}_A(V)^\top (I - (V V^\dagger)^\top) X  \mathrm{d} V   \right. \nonumber\\
			& \quad \left. - V^\dagger  \mathcal{F}_A(V) V^\dagger \mathrm{d} V   + \Sigma_{22} V^\dagger A \mathrm{d} V  \right],\nonumber
		\end{align}
		with $\mathcal{F}_A(V)$ defined in \eqref{eq:F_A}.
		Moreover, taking into account the structure of $\mathcal{C}(V)$ in \eqref{eq:CeCe}, we expand the second term in \eqref{trace_diffMB} as
		\begin{align}\label{eq:trace_diffMB2}
			&\mathrm{Tr}(\mathcal{P} \mathrm{d} \mathcal{C})
			\\= &
			- \mathrm{Tr}(\mathcal{P}_{12} \mathrm{d}{V}^\top C^\top C) - \mathrm{Tr}(\mathcal{P}_{12}^\top  C^\top C \mathrm{d}{V}) \nonumber \\
			& + \mathrm{Tr}\left[ \mathcal{P}_{22} (\mathrm{d}{V}^\top C^\top C V + {V}^\top C^\top C \mathrm{d} V) \right]
			\nonumber\\
			=  & 2 \mathrm{Tr}( -\mathcal{P}_{12}^\top  C^\top C \mathrm{d}{V} + \mathcal{P}_{22} {V}^\top C^\top C \mathrm{d} V) \nonumber
			\\
			= & 2 \mathrm{Tr}(\mathcal{F}_C(V)^\top \mathrm{d}{V}).
			\nonumber
		\end{align}

		Finally, we combine the three equations in \eqref{trace_MdiffB}, \eqref{eq:trace_diffMB1}, and \eqref{eq:trace_diffMB2} to obtain the gradient in \eqref{eq:nablaJ}
	\end{proof}

	\begin{remark}
		Different from methods in from e.g., \cite{vanDooren2008H2,sato2018pHsys}, which provide three gradient expressions over the variables $\hat{A}$, $\hat{B}$, and $\hat{C}$, respectively, this paper aims to find the optimal reduced basis matrix $V$ for projection, which is the only optimization variable in the optimization problem \eqref{eq:optimization1}. Theorem~\ref{thm:gradJ} thereby provides the expression of the gradient over $V$, which is clearly different from the aforementioned methods.
	\end{remark}

	\begin{remark}
		Note that the computation of the gradient $\nabla J$ does not need the Gramian matrices $\mathcal{P}_{11} \in \mathbb{R}^{n \times n}$ and $\mathcal{Q}_{11} \in \mathbb{R}^{n \times n}$ of the original system. It only requires to compute the Gramians for the reduced-order model $\mathcal{P}_{22} \in \mathbb{R}^{r \times r}$, $\mathcal{Q}_{22} \in \mathbb{R}^{r \times r}$ and the cross terms $\mathcal{P}_{12} \in \mathbb{R}^{n \times r}$, $\mathcal{Q}_{12} \in \mathbb{R}^{n \times r}$. The computation of $\mathcal{P}_{22}$ and $\mathcal{Q}_{22}$ can be efficient as the reduced order $r$ is typically small in model reduction problems. Nevertheless, solving the Sylvester equations in \eqref{eq:Sylv} can still be computationally demanding due to \eqref{eq:Sylv_c} and \eqref{eq:Sylv_d} that contain the matrices from the original system. To make the computation of $\mathcal{P}_{12}$ and $\mathcal{Q}_{12}$ more tractable, we follow \cite{sorensen2002sylvester} and derive explicit expressions for  $\mathcal{P}_{12}$ and $\mathcal{Q}_{12}$ as the sums of rank-one matrices:
		\begin{subequations}
			\label{eq:expSylv}
			\begin{align}
				\mathcal{P}_{12} & = \sum_{j=1}^{r} (\mu_j I_r + {A})^{-1} B \hat{B}^\top  w_j \hat{w}_j^*, \\
				\mathcal{Q}_{12} & = \sum_{j=1}^{r} (\mu_j I_r + {A}^\top)^{-1} {C} \hat{C}^\top  \hat{w}_j w_j^*,
			\end{align}
		\end{subequations}
		where the vectors $w_j$ and $\hat{w}_j$ are the right and left eigenvectors corresponding to the eigenvalue $\mu_j$ of $\hat{A}$, respectively. Computing \eqref{eq:expSylv} now requires  the eigendecomposition of the reduced matrix $\hat{A}$ followed by $r$ inversions of $n \times n$ matrices $\mu_j I  + A$, which may be less expensive than directly applying standard Sylvester equation solvers, particularly when $A$ is sparse, and $r$ is small.
		To further accelerate the method, it is advisable to adopt a low-rank method, see e.g., \cite{jbilou2006low-rank,benner2014low-rank}, to achieve approximate solutions with cheap computations for extremely large-scale systems.
	\end{remark}

	Having the formula for the gradient $\nabla J$ in \eqref{eq:nablaJ}, we further define a gradient flow of $J(V)$ as
	\begin{align} \label{eq:gradflow}
		\frac{\mathrm{d}}{\mathrm{d} t}{V}(t)
		= & -2 X (I - V V^\dagger)  \left[ \mathcal{F}_A (V) + \mathcal{F}_B (V)^\top \right]   (V^\top X V)^{-1}
		\nonumber \\
		& + 2(V^\dagger)^\top  \left[ \mathcal{F}_A (V)^\top + \mathcal{F}_B (V) \right]  (V^\dagger)^\top
		\nonumber
		\\& -  2 A^\top (V^\dagger)^\top \Sigma_{22}^\top   - 2 \mathcal{F}_C (V),
	\end{align}
	with $\mathcal{F}_A (V)$, $\mathcal{F}_B (V)$, and $\mathcal{F}_C (V)$ defined in \eqref{eq:FAFBFC}, and then we analyze the key properties of the gradient flow.
	\begin{theorem} \label{thm:gradflow}
		Let the initial condition of \eqref{eq:gradflow} be given by $V(0) = V_0 \in \mathbb{S}^{n \times r}_*$. Then, the following statements hold.
		\begin{enumerate}
			\item $V(t)^\top \frac{\mathrm{d}}{\mathrm{d} t}{V(t)} = 0$, and  $V(t)^\top V(t)$ is invariant in the gradient flow  \eqref{eq:gradflow}.
			\item $V(t) \in \mathbb{S}^{n \times r}_*$ for all $t \geq 0$.
			\item The objective function is nonincreasing along $V(t)$ with
			\begin{align} \label{eq:nonincreasing}
				J(V(t_1))  \geq J(V(t_2)),
			\end{align}
			for all $0 \leq t_1 \leq t_2$.
			\item Let $V$ be any local minimizer of problem \eqref{eq:optimization1}, the following implication holds
			\begin{equation} \label{eq:JV0}
				\nabla J(V) = 0   \ \Leftrightarrow \ \nabla J(V) V^\top -  V \nabla J(V)^\top = 0.
			\end{equation}
		\end{enumerate}
	\end{theorem}
	\begin{proof}
		We first prove $V^\top \frac{\mathrm{d}}{\mathrm{d} t}{V} = 0$ for any $V \in \mathbb{S}^{n \times r}_*$. We rewrite \eqref{eq:gradflow} as
		$\frac{\mathrm{d}}{\mathrm{d} t}{V} = 2(\mathcal{J}_1 + \mathcal{J}_2)$, where
		\begin{align*}
			\mathcal{J}_1
			= & - X (I - V V^\dagger)  \left[ \mathcal{F}_A (V) + \mathcal{F}_B (V)^\top \right]   (V^\top X V)^{-1},
			\\
			\mathcal{J}_2 = &  (V^\dagger)^\top  \left[ \mathcal{F}_A (V)^\top + \mathcal{F}_B (V) \right]  (V^\dagger)^\top
			\nonumber
			\\& -   A^\top (V^\dagger)^\top \Sigma_{22}^\top   - \mathcal{F}_C (V),
		\end{align*}
		Note that
		$
		V^\top X (I - V V^\dagger) = V^\top X - V^\top X V V^\dagger = 0,
		$
		which leads to
		\begin{equation} \label{eq:VJ1}
			V^\top \mathcal{J}_1 = 0.
		\end{equation}
		Furthermore, we have
		\begin{align} \label{eq:VJ2}
			V^\top \mathcal{J}_2 = &  \Sigma_{22}^\top \hat{A}  + \mathcal{Q}_{12}^\top B \hat{B}^\top +  \mathcal{Q}_{22} \hat{B} \hat{B}^\top
			\nonumber \\
			& - \hat{A}^\top  \Sigma_{22} - \hat{C}^\top \hat{C} \mathcal{P}_{22} + \hat{C}^\top C \mathcal{P}_{12}
		\end{align}
		and it is obtained from the Sylvester equations in \eqref{eq:Sylv} that
		\begin{align*}
			- \mathcal{Q}_{12}^\top B \hat{B}^\top & = \mathcal{Q}_{12}^\top A \mathcal{P}_{12} + \mathcal{Q}_{12}^\top \mathcal{P}_{12} \hat{A}^\top,
			\\
			-  \mathcal{Q}_{22} \hat{B} \hat{B}^\top
			& = \mathcal{Q}_{22} \hat{A} \mathcal{P}_{22} + \mathcal{Q}_{22} \mathcal{P}_{22} \hat{A}^\top,
			\\
			\hat{C}^\top \hat{C} \mathcal{P}_{22}
			& = - \hat{A}^\top \mathcal{Q}_{22} \mathcal{P}_{22} - \mathcal{Q}_{22} \hat{A} \mathcal{P}_{22},
			\\
			- \hat{C}^\top C \mathcal{P}_{12}
			& = - \mathcal{Q}_{12}^\top A \mathcal{P}_{12} - \hat{A}^\top \mathcal{Q}_{12}^\top \mathcal{P}_{12}.
		\end{align*}
		Substituting the above equations into \eqref{eq:VJ2} then gives
		\begin{align} \label{eq:VJ2-2}
			V^\top \mathcal{J}_2 = &   \Sigma_{22}^\top \hat{A} - \mathcal{Q}_{12}^\top \mathcal{P}_{12} \hat{A}^\top - \mathcal{Q}_{22} \mathcal{P}_{22} \hat{A}^\top
			\nonumber \\
			& - \hat{A}^\top \Sigma_{22}  + \hat{A}^\top \mathcal{Q}_{22}\mathcal{P}_{22} + \hat{A}^\top \mathcal{Q}_{12}^\top \mathcal{P}_{12}
			= 0.
		\end{align}
		Therefore, it follows from \eqref{eq:VJ1} and \eqref{eq:VJ2-2} that
		\begin{equation} \label{eq:VdV0}
			V^\top \frac{\mathrm{d}}{\mathrm{d} t}{V} = 2 V^\top \mathcal{J}_1 + 2 V^\top \mathcal{J}_2 = 0,
		\end{equation}
		Furthermore,
		$$\frac{\mathrm{d}}{\mathrm{d} t} (V(t)^\top V(t)) = \frac{\mathrm{d}}{\mathrm{d} t}{V}(t)^\top V(t) + V(t)^\top \frac{\mathrm{d}}{\mathrm{d} t}{V}(t) = 0,$$ that proves the first statement.

		The second statement can be shown with the invariance of $V(t)^\top V(t)$ in the gradient flow.
		It is clear that
		$
		V(t)^\top V(t) = V(0)^\top V(0) = V_0^\top V_0,
		$
		which gives
		\begin{equation*}
			\mathrm{Rank} (V(t)) =

			\mathrm{Rank} (V_0).
		\end{equation*}
		For arbitrary initial condition $V_0 \in \mathbb{S}^{n \times r}_*$, the matrix $V(t)$ at any $t \geq 0$ stays in $\mathbb{S}^{n \times r}_*$.

		Next, we prove the third statement by writing the derivative of $J(V)$ as follows.
		\begin{equation*}
			\frac{\mathrm{d}}{\mathrm{d} t}{J} (V(t)) = \langle \nabla J(V), \frac{\mathrm{d}}{\mathrm{d} t}{V} \rangle =  - 2  \| \nabla J(V) \|^2_F \leq 0.
		\end{equation*}
		Therefore, $J(V)$ is nonincreasing, i.e. \eqref{eq:nonincreasing} holds.

		Finally, denote $W: = \nabla J(V) V^\top -  V \nabla J(V)^\top$, and it is clear that $W = 0$ if $\nabla J(V) = 0$. Moreover, we have $\nabla J(V)^\top V = 0$, $\forall V \in \mathbb{S}^{n \times r}_*$ due to \eqref{eq:VdV0}, then
		$$ W V (V^\top V)^{-1} =  J(V),
		$$
		and thus $J(V) = 0$ if $W = 0$. The implication \eqref{eq:JV0} is obtained.
	\end{proof}

	The features of the gradient flow $\frac{\mathrm{d}}{\mathrm{d} t} V(t)$ reflects that the solution of the ordinary differential equation \eqref{eq:gradflow} always lies on the manifold $\mathbb{S}^{n \times r}_*$ over all the time for
	any given initial condition $V_0$, and moreover, the solution will converge to a critical point of $J$ on $\mathbb{S}^{n \times r}_*$ as the objective function $J(V)$ is always non-negative. It is guaranteed that if the objective function $J(V)$ has only isolated minimum points, the solution $V(t)$ should converge to one of them.

	With the matrix differential equation \eqref{eq:gradflow}, we can apply numerical algorithms that search along a geodesic to solve this gradient
	flow problem as in e.g., \cite{Yu1999approximate,edelman1998geometry,Jiang2015h2}.
	However, this procedure may require a considerable computational resource and hence is not discussed in this paper. Alternatively, we provide iterative algorithms based on gradient descent to solve the nonconvex optimization problem \eqref{eq:optimization1} in a more efficient manner.

	\subsection{Gradient-Based Algorithms}

	With the explicit expression of the gradient $\nabla J(V)$,
	the gradient descent method can be directly used, which has the following simple iteration:

	\begin{algorithm}
		\caption{Gradient Decent Algorithm}
		\begin{algorithmic}[1]
			\State Let $X$ be chosen based on the structural property of the original system \eqref{sys:orig} and $V_0 \in \mathbb{S}^{n \times r}_*$ be an initial solution.
			\Repeat

			\State Solve the Sylvester equations in \eqref{eq:Sylv}.

			\State Compute the gradient $\nabla J(V_k)$ in \eqref{eq:nablaJ}.

			\State Update $V_{k+1}  = V_k - \alpha_k \nabla J(V_k)$.

			\Until{$\| V_{k+1} - V_k \|_{F}^2 \leq \epsilon$}
		\end{algorithmic}
		\label{alg}
	\end{algorithm}
	Algorithm~\ref{alg} terminates when the difference between $V_{k+1}$ and $V_k$ is small enough, i.e. their error norm is less than a small positive scalar $\epsilon$. More detailed discussions on how to choose the stepsize $\alpha_k > 0$ in each iteration $k$ and an initial solution $V_0$ are provided in Remark~\ref{rem:stepsize} and Remark~\ref{rem:initialization}, respectively. Generally, the stepsize $\alpha_k$ should be chosen such that the iteration in Algorithm~\ref{alg} converges. However, the convergence of Algorithm~\ref{alg} is hard to be analyzed theoretically since $\mathbb{S}^{n \times r}_*$ is a noncompact set.

	\begin{coro}
		\label{coro:iter}
		Let $V_0 \in \mathbb{S}^{n \times r}_*$ be the initial condition and $\{V_k\}$ be the generated sequence in Algorithm~\ref{alg}. Then,  $V_k \in \mathbb{S}^{n \times r}_*$ for each iteration $k$.
	\end{coro}
	\begin{proof}
		It is equivalent to show that  $V_{k+1}^\top V_{k+1} > 0$ if $V_{k}^\top V_{k} > 0$.
		Note that
		\begin{align*}
			V_{k+1}^\top V_{k+1} = & V_{k}^\top V_{k} + \alpha_k^2 \nabla J(V_k)^\top \nabla J(V_k) \\
			&- \alpha_k V_{k}^\top   \nabla J(V_k) - \alpha_k \nabla J(V_k)^\top  V_{k}.
		\end{align*}
		We can follow the proof of Theorem~\ref{thm:gradflow} to show that
		\begin{equation}
			V_{k}^\top   \nabla J(V_k) = 0, \ \forall k,
		\end{equation}
		which gives
		\begin{equation*}
			V_{k+1}^\top V_{k+1} =  V_{k}^\top V_{k} + \alpha_k^2 \nabla J(V_k)^\top \nabla J(V_k) > 0.
		\end{equation*}
		Therefore, $V_{k}$ stays in the manifold $\mathbb{S}^{n \times r}_*$ given $V_0 \in \mathbb{S}^{n \times r}_*$.
	\end{proof}

	\begin{remark}
		\label{rem:stepsize}
		The stepsize $\alpha_k$ in Algorithm~\ref{alg} can be fixed. Alternatively, we apply the steepest descent method \cite{Meza2010steepest} to determine the stepsize via an exact line search as follows:
		\begin{equation}\label{eq:steepest}
			\alpha_k = \argmin\limits_{\alpha > 0} J(V_k - \alpha \nabla J(V_k)),
		\end{equation}
		which aims to find the best stepsize along the direction of the gradient at each iteration. We refer to \cite{Meza2010steepest} for more details.
	\end{remark}

	It is noted that the noncompact Stiefel manifold $\mathbb{S}^{n \times r}_*$ is unbounded, and thus it is difficult to show the Lipschitz continuity of the gradient $\nabla J(V)$.
	To establish a convergence guarantee that can be theoretically proven,
	we discuss in the following an iterative optimization algorithm based on a curvilinear search that can generate a sequence of feasible points on an invariant submanifold of $\mathbb{S}^{n \times r}_*$:
	\begin{equation} \label{eq:Stiefel_Comp}
		\mathcal{M} = \{V \in  \mathbb{S}^{n \times r}_* \mid V^\top V = M \},
	\end{equation}
	where $M \in \mathbb{R}^{r \times r}$ is positive definite. In particular, if $M = I_r$, then $\mathcal{M}$ is the compact \textit{Stiefel manifold}, denoted by $\mathrm{St}(n, r)$.
	To this end, we employ the so-called Cayley transform that computes a parametric curve at the point $V \in \mathcal{M}$, with the following closed form:
	\begin{equation} \label{eq:Cayley}
		\mathcal{V}  (\tau) = \left(I - \frac{\tau}{2} W\right)^{-1} \left(I + \frac{\tau}{2} W\right) V ,
	\end{equation}
	where $\tau$ is
	a parameter that represents the length on the curve, and  $W$ is a skew-symmetric matrix given as
	\begin{equation} \label{eq:W}
		W  = \nabla J(V) V^\top -  V \nabla J(V)^\top.
	\end{equation}
	We refer to e.g., \cite{wen2013feasible,Lanlin2020h2} for more details. It is observed that the Cayley transform \eqref{eq:Cayley} leads to interesting properties: $\mathcal{V}  (\tau)$ is a smooth function of $\tau$ with $\mathcal{V}  (0) = V$, and moreover,
	$
	\mathcal{V}(\tau)^\top \mathcal{V} (\tau) = V^\top V,
	$
	for all $\tau \in \mathbb{R}$.

	However, the Cayley transform \eqref{eq:Cayley} involves the computation of a matrix inversion of dimension $n$, and hence the iterative algorithm can be computationally expensive when the original system \eqref{sys:orig} is large-scale. To reduce the computational effort, we make use of the property $\nabla J(V)^\top V = 0$ proven in Theorem~\ref{thm:gradflow} and provide a new update scheme as follows. For simplicity, $\nabla J(V)$ is briefly denoted by $\nabla J$ here.
	\begin{lemma}
		Consider $\mathcal{V}(t)$ in the Cayley transform  \eqref{eq:Cayley} at any point $V \in \mathbb{S}^{n \times r}_*$ with $\tau \in \mathbb{R}$ and the gradient $\nabla J$ in  \eqref{eq:nablaJ}. Then, $\mathcal{V}(t)$ can  be equivalently computed as
		\begin{equation} \label{eq:Vt_new}
			\mathcal{V}(\tau) = V - \tau \left(I+\frac{\tau}{2} V \nabla J^\top \right) \nabla J   \mathcal{U}(\tau)^{-1} V^\top V,
		\end{equation}
		where
		$ \mathcal{U}(\tau) = I +  \frac{\tau^2}{4} V^\top V \nabla J^\top \nabla J.
		$
		If $V \in \mathrm{St}(n, r)$,  the expression \eqref{eq:Vt_new} is simplified as
		\begin{equation} \label{eq:Vt_new1}
			\mathcal{V}(\tau) = V - \tau \left(I+\frac{\tau}{2} V \nabla J^\top \right) \nabla J    \left(I +  \frac{\tau^2}{4}  \nabla J^\top \nabla J\right)^{-1}.
		\end{equation}
	\end{lemma}
	\begin{proof}
		Observe that the matrix $W$ in \eqref{eq:W} can be expressed as the outer product of two low-rank matrices as $W = W_L W_R^\top$ with
		\begin{equation} \label{eq:WLWR}
			W_L = [
			\nabla J (V), \ -V
			], \ \text{and} \ W_R = [
			V, \ \nabla J (V)].
		\end{equation}
		Then, it follows from the matrix inversion lemma that
		\begin{equation*}
			\left(I - \frac{\tau}{2} W\right)^{-1} = I - \frac{\tau}{2} W_L \left(I + \frac{\tau}{2} W_R^\top W_L \right)^{-1} W_R^\top.
		\end{equation*}
		Therefore, we have
		\begin{align} \label{eq:Vt}
			\mathcal{V}(t) & = \left(I - \frac{\tau}{2} W\right)^{-1} \left(I + \frac{\tau}{2} W\right) V
			\nonumber\\
			& = \left[  2\left(I - \frac{\tau}{2} W \right)^{-1} - I \right] V
			\nonumber\\
			& = \left[I - \tau W_L \left(I + \frac{\tau}{2} W_R^\top W_L \right)^{-1} W_R^\top \right] V.
		\end{align}
		Denote $Z: = W_L (I + \frac{\tau}{2} W_R^\top W_L )^{-1} W_R^\top  V$. With the expressions in \eqref{eq:WLWR} into \eqref{eq:Vt}, $Z$ can be rewritten as
		\begin{align*}
			Z & =
			\begin{bmatrix}
				\nabla J & - V
			\end{bmatrix}
			\begin{bmatrix}
				I & \frac{\tau}{2} V^\top V \\
				\frac{\tau}{2} \nabla J^\top \nabla J & I
			\end{bmatrix}^{-1}
			\begin{bmatrix}
				V^\top V \\ 0
			\end{bmatrix}
			\\

			& =  \begin{bmatrix}
				\nabla J  & - V
			\end{bmatrix}
			\begin{bmatrix}
				\mathcal{U}(\tau)^{-1} V^\top V  \\
				\frac{\tau}{2}\nabla J^\top \nabla J \mathcal{U}(\tau)^{-1} V^\top V
			\end{bmatrix}
			\\
			& = \left(I + \frac{\tau}{2} V \nabla J^\top \right) \nabla J \mathcal{U}(\tau)^{-1} V^\top V,
		\end{align*}
		which is substituted into \eqref{eq:Vt} and then yields \eqref{eq:Vt_new}.
	\end{proof}

	The merit of this result is that, rather than computing the inverse of the $n \times n$ matrix $I - \frac{\tau}{2}W$ in \eqref{eq:Cayley}, it presents an alternative to compute  $\mathcal{V}(\tau)$ which only requires the inverse of the $r \times r$ matrix $\mathcal{U}(\tau)$. As in our model reduction problem, the reduced order $r$ is typically much smaller than $n$, and the new update scheme in \eqref{eq:Vt_new} can significantly improve the computational efficiency.

	With the new curve expression \eqref{eq:Vt_new}, the following curvilinear search approach is thereby provided to solve the nonconvex optimization problem \eqref{eq:optimization1}, and the solution $V_k$ in each iteration $k$ remains in the manifold $\mathcal{M}$.
	\begin{algorithm}
		\caption{Iterative Algorithm Based on Curvilinear Search}
		\begin{algorithmic}[1]
			\State Let $X$ be chosen based on the structural property of the original system \eqref{sys:orig} and $V_0 \in \mathbb{S}^{n \times r}_*$ be an initial solution.
			\State $M = V_0^\top V_0$.
			\Repeat

			\State Solve the Sylvester equations in \eqref{eq:Sylv}.

			\State Compute the gradient $\nabla J(V_k)$ in \eqref{eq:nablaJ}.

			\State Update
			$$V_{k+1} = V_k - \tau_k \left(I+\frac{\tau_k}{2} V_k \nabla J(V_k)^\top \right) \nabla J(V_k)   \mathcal{U}_k^{-1}M.$$
			where $\mathcal{U}_k: =  I +  \frac{\tau_k^2}{4} M \nabla J(V_k)^\top \nabla J(V_k)$

			\Until{$\| \nabla J(V_k) \|_{F}^2 \leq \epsilon$}
		\end{algorithmic}
		\label{alg2}
	\end{algorithm}

	Note that Algorithm~\ref{alg2} can treat the optimization problem on the Stiefel manifold $\mathrm{St}(n, r)$ as a special case. It is simply initialized by any $V_0 \in \mathrm{St}(n,r)$ then applies \eqref{eq:Vt_new1} for the update, and the orthogonality constraint $V_k^\top V_k = I_r$ will be retained for all $k$.

	\begin{remark}
		\label{rem:stepsize2}
		At iteration $k$, $\mathcal{V}_k(\tau)$ is the function \eqref{eq:Vt_new} obtained at the point $V_k$, and we can select the stepsize $\tau_k \geq 0$ by minimizing $J(\mathcal{V}_k(\tau))$ along the curve $\mathcal{V}_k(\tau)$ with respect to $\tau$. However, it is computationally demanding to
		obtain the minimizer. Therefore, we can also select the stepsize $\tau_k$ to fulfill  the following Armijo-Wolfe conditions \cite{jorge2006numerical,wen2013feasible}:
		\begin{align}
			J(\mathcal{V}_k(\tau_k)) & \leq J(\mathcal{V}_k(0)) + \rho_1 \tau_k J_\tau'(\mathcal{V}_k(0)),
			\label{eq:Wolfe}\\
			J_\tau'(\mathcal{V}_k(\tau_k)) & \geq \rho_2 J_\tau'(\mathcal{V}_k(0)),
			\nonumber
		\end{align}
		where $0 <\rho_1 < \rho_2 < 1$ are two parameters, and $J_\tau'(\mathcal{V}(\tau))$ denotes the derivative of $J(\mathcal{V}(\tau))$ with respect to $\tau$.
		It is worth emphasizing that there always exists a stepsize $\tau_k$
		satisfying the above
		conditions since $J(\mathcal{V}(\tau))$ is continuously differentiable and bounded, see \cite{jorge2006numerical} for more details.
	\end{remark}

	Finally, the following corollary is presented to formalize the convergence of Algorithm~\ref{alg2}.
	\begin{coro}
			Let $V_0 \in \mathbb{S}^{n \times r}_*$ be the initial condition and $\{V_k\}$ be the generated sequence in Algorithm~\ref{alg2}. Denote $M: = V_0^\top V_0$. If the stepsize $\tau_k$ satisfies the Armijo-Wolfe condition \eqref{eq:Wolfe}, then the sequence $\{V_k\}$ has the following properties.
		\begin{enumerate}
			\item $V_k \in \mathbb{S}^{n \times r}_*$, and $V_k^\top V_k = M$ for each iteration $k$.
			\item the objective function $J(V)$ is non-increasing along the iterations as
			\begin{align} \label{eq:convergence}
				J(V_{k+1}) & \leq J(V_k) - \rho_1 \tau_k \| V_k \nabla J(V_k)^\top \|_F^2 \\
				 & {  \leq J(V_k) - \rho_1  \tau_k \sigma_{\mathrm{min}}(M) \|\nabla J(V_k)^\top \|_F^2} \nonumber
			\end{align}
			for all $k \geq 0$, {  where $\sigma_{\mathrm{m}}(M)$ is the smallest singular value of the constant matrix $M$.}

			\item {  The sequence $\{V_k\}$ converges to a critical point, i.e. $\lim\limits_{k \rightarrow \infty} \| \nabla J(V_k) \|_F = \lim\limits_{k \rightarrow \infty} \| V_{k+1}  - V_k \|_F = 0$.}

						\item {  It takes at most
				\[
				\kappa = \left\lceil  \frac{J(V_0) - J(V_*)}{ \rho_1  \nu \sigma_{\mathrm{min}}(M) \epsilon} \right\rceil
				\] iterations for Algorithm~\ref{alg2} to terminate, where $\nu: = \inf\{\tau_0, \tau_1, ..., \tau_{\kappa-1} \}$ is strictly positive, and $\lceil  \cdot \rceil$ returns the ceiling of a number.}
		\end{enumerate}
	\end{coro}
	\begin{proof}
		The first statement is a direct result of the invariance property of the Cayley transform \eqref{eq:Cayley}. The convergence of Algorithm~\ref{alg2} is guaranteed by the Wolfe condition \eqref{eq:Wolfe}. Using the chain rule,
		\begin{equation}
			J_\tau'(\mathcal{V}(t)) = \mathrm{Tr} \left( \frac{\partial J(\mathcal{V}(\tau))}{\partial \mathcal{V}(\tau)}^\top
			\cdot \frac{\partial \mathcal{V}(\tau)}{\partial \tau} \right),
		\end{equation}
		where
		\begin{equation*}
			\frac{\partial \mathcal{V}(\tau)}{\partial \tau} = - (\nabla J V^\top - V \nabla J^\top ) V = - \nabla J V^\top V,
		\end{equation*}
		and $\frac{\partial J(\mathcal{V}(0))}{\partial \mathcal{V}(0)} = \nabla J(V)$ as $\mathcal{V}(0) = V$.
	Hence,
 $$J_\tau'(\mathcal{V}(0)) = - \mathrm{Tr} (\nabla J^\top \nabla J V^\top V) =
		- \mathrm{Tr} (V \nabla J^\top \nabla J V^\top    )
		.$$  Then the first inequality in \eqref{eq:convergence} is obtained, followed by the condition \eqref{eq:Wolfe}. 		{  The second inequality in \eqref{eq:convergence} is induced from the fact that
			$\mathrm{Tr} (\nabla J^\top \nabla J V^\top V) = \mathrm{Tr} ( \nabla J M \nabla J^\top) \geq \sigma_{\mathrm{m}}(M) \mathrm{Tr} (\nabla J^\top \nabla J)$.

		Since $\rho_1$, $\tau_k$, and $\sigma_{\mathrm{m}}(M)$ are all positive scalars, $J(V_k)$ is monotone decreasing along $k$. Besides, $J(V_{k})$ is always nonnegative by its definition \eqref{eq:optimization1}. Therefore, the sequence $\{J(V_{k})\}$ converges, which implies $\lim_{k \rightarrow \infty}\left[ J(V_k) - J(V_{k+1})\right] = 0$ and hence $
		\lim_{k \to \infty}  \| \nabla J(V_k) \|_F = 0
		$ from \eqref{eq:convergence}. As a result, $\mathcal{U}_k$ in Algorithm~\ref{alg2} converges to the identity matrix as $k \rightarrow \infty$, which leads to $\lim_{k \rightarrow \infty} \| V_{k+1}  - V_k \|_F = 0$. We would also like to refer to \cite{goldfarb2009curvilinear} for the convergence analysis of general curvilinear search methods. The analysis of the convergence rate in item (4)  can be proven following the similar reasoning as in Corollary~\ref{coro:iter} and hence is omitted here.}
	\end{proof}

	\begin{remark}

		Algorithm~\ref{alg2} can also be viewed as a Riemannian optimization method due to the property $\nabla J(V)^\top V = 0$, for all $V \in \mathbb{S}^{n \times r}_*$. Following \cite{Absil2009optimization}, the Riemannian gradient can be written as
		\begin{equation*}
			\mathrm{grad} J(V) : = \nabla J(V) - \frac{1}{2} V (V^\top \nabla J(V) + \nabla J(V)^\top V),
		\end{equation*}
		which is identical to $\nabla J(V)$ whenever $\nabla J(V)^\top V = 0$. Furthermore, the Cayley transformation in \eqref{eq:Cayley} or \eqref{eq:Vt_new} is adopted to preserve the constraint $V_k^\top V_k = M$. This is actually equivalent to a retraction map used in Riemannian gradient descent.
	\end{remark}

	\begin{remark}
		\label{rem:initialization}
		The initialization of $V_0$ in Algorithms~\ref{alg} and \ref{alg2} is discussed. Generally, we can choose $V_0$ as a random $n \times r$ matrix with full column rank. Then, $V_k$ in each iteration of Algorithm~\ref{alg} will have full column rank, and furthermore, $V_k^\top V_k$ will remain constant in Algorithm~\ref{alg2}.
		However, a random initial condition might result in a slow convergence and produce a local optimum with the objective function value much larger than the global minimum value. In contrast, it is preferable to  choose the initial point $V_0$ as an output of the balanced truncation or the Krylov subspace methods. In this manner, the presented approach can be viewed as a posterior procedure of the balanced truncation or the Krylov subspace methods for improving the accuracy of reduced-order models. We will demonstrate this point in the following numerical example.
	\end{remark}

	%

	\section{Numerical Experiments}
	\label{sec:example}

	In this section, the preceding theoretical results of the proposed optimal projection-based model reduction method are illustrated. In the first part, we implement the proposed framework for the stability-preserving model reduction of a building model, which is compared with the standard balanced truncation \cite{morAntSG01} and the projection-based method
	\cite{selga2012stability}. The second part, a model reduction problem for port-Hamiltonian systems is addressed, and comparisons are made among the proposed method and the other existing approaches in \cite{polyuga2012effort,gugercin2012pH,sato2018pHsys}.

	\subsection{Stability-Preserving Model Reduction of Building Model}

	In this example, we aim to preserve the asymptotic stability of a system. Consider a model reduction benchmark example in \cite{morAntSG01} that describes the displacement of a multi-storey building. The full-order model is asymptotically stable and has $48$ state variables, $1$ input, and $1$ output. The objective in this example is to build reduced-order models that remain asymptotically stable.  Note that the optimization-based approaches in  \cite{Yu1999approximate,sato2017structure} are not applicable to this example, since the symmetric part of the system state matrix, i.e. $\frac{1}{2}(A + A^\top)$ with $A$ in \eqref{sys:orig}, is not negative definite.

	For comparison, three projection-based methods are applied: (i) the standard balanced truncation (BT) that naturally preserves stability \cite{morAntSG01}, (ii) the iterative SVD-rational Krylov-based model reduction approach (ISRK) in \cite{gugercin2008iterative}, and (iii) the proposed method based on structure-preserving optimal projection (SPOP) in Algorithm~\ref{alg}, where $X$ in \eqref{eq:pinv} is equal to the observability Gramian of the original system, and the stepsize of each iteration is determined by \eqref{eq:steepest}. Meanwhile, to demonstrate the performance of Algorithm~\ref{alg}, we adopt three different strategies for its initialization. First, we use five randomly generated matrices as initial points to run the iterations, and the one that produces the smallest error is marked as SPOP-R. Besides, we also choose the initial points $V_0$ of  Algorithm~\ref{alg} as the outputs of BT and ISRK, and the corresponding schemes are denoted by SPOP-B and SPOP-K, respectively.

	For each reduced order $r = 3, 6, 9, 12, 15$, we compute and show in Table~\ref{tab:building} the relative errors in terms of the $\mathcal{H}_2$ norm, i.e. $\| {G}(s) - \hat{G}(s) \|_{\mathcal{H}_2}/ \| {G}(s)  \|_{\mathcal{H}_2}$ with ${G}(s)$, $\hat{G}(s)$ the transfer matrices of the full-order and reduced-order models, respectively.
	\begin{table}[t]
		\centering
		\caption{Relative $\mathcal{H}_2$ errors as the reduced order $r$ varies}
		\begin{tabular}[t]{lccccc}
			\toprule
			Order & $r=3$ & $r=6$ & $r=9$ & $r=12$ & $r=15$\\
			\midrule
			BT     & 0.7170 & 0.2905 & 0.2217 & 0.1650 & 0.1644\\
			ISRK   & 0.7945 & 0.2486 & 0.2010 & 0.1399 & 0.1132\\
			SPOP-R & 0.7153 & 0.2607 & 0.2063 & 0.1730 & 0.1583\\
			SPOP-B & \textbf{0.7145} & \textbf{0.2460} & \textbf{0.1761} & 0.1367 & 0.0971\\
			SPOP-K & 0.7809 & \textbf{0.2460} & 0.1780 & \textbf{0.1364} & \textbf{0.0921}\\
			\bottomrule
		\end{tabular}
		\label{tab:building}
	\end{table}%
	From Table~\ref{tab:building}, we see that Algorithm~\ref{alg} is sensitive to the choices of its initial point $V_0$, as with different initial points, the algorithm can converge to distinct local optimal solutions that give different relative errors. When $r = 12$ and $r = 15$, the random selection of initial points of $V_0$ in Algorithm~\ref{alg} produce reduced-order models with greater relative $\mathcal{H}_2$ errors than the ones generated by BT and ISRK. Therefore, the selection of the initial point $V_0$ in Algorithm~\ref{alg} is crucial for achieving a smaller reduction error.

	\begin{figure}[t]
		\centering \includegraphics[width=0.4\textwidth]{fig/bode.pdf}
		\caption{Comparison of the relative error magnitudes versus frequency in different methods when $r = 15$.}
		\label{fig:bode}
	\end{figure}

	Furthermore, when the outputs of BT or ISRK are taken as the initial points, Algorithm~\ref{alg} presents better results than both BT and ISRK. This is because the projections generated by the BT and ISRK methods are generally not optimal in minimizing the approximation error, while the proposed scheme can take their outputs as initial points and leads to a local optimal projection via the iterations on the noncompact Stiefel manifold. Particularly, when $r = 15$, SPOP-B and SPOP-K further reduce the relative errors by $41.0\%$ and $18.6\%$, compared to BT and ISRK, respectively. The detailed comparison of the approximation error magnitudes at different frequencies is presented in Fig.~\ref{fig:bode}, which shows that the SPOP-B and SPOP-K methods outperform the BT or ISRK approaches for most frequencies, particularly in the lower frequency range. This shows the potential of the presented method as a posterior procedure of the balanced truncation or the Krylov subspace methods, as it is able to further improve the quality of the projection and hence gives a better reduced-order model.

	\subsection{Passivity-Preserving Model Reduction of Mass-Spring-Damper System}

	We further show the effectiveness of the proposed framework in the passivity-preserving model reduction problem for port-Hamiltonian systems, which has received considerable attention from the model reduction literature. We adopt the benchmark example of a mass-spring-damper system as shown in Fig.~\ref{fig:msd}, with masses $m_i$, spring  coefficients $k_i$ and damping constants $c_i$, $i = 1,...,\ell$. The displacement of the mass $m_i$ is indicated by $q_i$. The same model parameters are used as in \cite{gugercin2012pH}, which considers $\ell = 50$ masses and hence $100$ states in the full-order model. Furthermore, the system inputs are two external forces $u_1$, $u_2$ applied to the masses $m_1$, $m_2$, whose velocities $\dot{q}_1$ and $\dot{q}_2$ are measured as the outputs $y_1$, $y_2$ of the system. We refer to \cite{gugercin2012pH} for the detailed state-space formulation of this mass-spring-damper system.

	\begin{figure}
		\centering
		\begin{tikzpicture}[scale=0.9, transform shape, every node/.style={outer sep=0pt},thick,
			mass/.style={draw,thick},
			spring/.style={thick,decorate,decoration={zigzag,pre length=0.3cm,post
					length=0.3cm,segment length=6}},
			ground/.style={fill,pattern=north east lines,draw=none,minimum
				width=0.75cm,minimum height=0.3cm},
			dampic/.pic={\fill[white] (-0.1,-0.2) rectangle (0.2,0.2);
				\draw (-0.3,0.2) -| (0.2,-0.2) -- (-0.3,-0.2);
				\draw[line width=1mm] (-0.1,-0.2) -- (-0.1,0.2);}]

			\node[mass,minimum width=1.1cm,minimum height=0.9cm,fill={rgb:orange,1;yellow,2;pink,5}] (m1) {$m_1$};
			\node[mass,minimum width=1.1cm,minimum height=0.9cm,fill={rgb:orange,1;yellow,2;pink,5},right=1.4cm of
			m1] (m2) {$m_2$};
			\node[mass,minimum width=1.1cm,minimum height=0.9cm,fill={rgb:orange,1;yellow,2;pink,5},right=2.1cm of
			m2] (m3) {$m_{\ell}$};
			\node[right=1.4cm of m2] (dots) {$\cdots$};
			\node[right=1.4cm of m3,ground,minimum width=3mm,minimum height=2cm] (g1){};

			\draw[thick] (g1.south west) -- (g1.north west);

			\draw[spring] ([yshift=3mm]g1.west) coordinate(aux)
			-- (m3.east|-aux)
			node[midway,above=1mm]{$k_{\ell}$};
			\draw[spring]  (m1.east|-aux) -- (m2.west|-aux) node[midway,above=1mm]{$k_1$};

			\draw[spring]  (m2.east|-aux) -- (dots.west|-aux) node[midway,above=1mm]{$k_2$};

			\draw ([yshift=-3mm]g1.west) coordinate(aux')
			-- (m3.east|-aux') pic[midway]{dampic} node[midway,below=3mm]{$c_\ell$}
			(m1.east|-aux') -- (m2.west|-aux') pic[midway]{dampic} node[midway,below=3mm]{$c_1$}
			(m2.east|-aux') -- (dots.west|-aux') pic[midway]{dampic} node[midway,below=3mm]{$c_2$};

			\foreach \X in {1,2}
			{
				\draw (m\X.north) -- ++ (0,0.5) coordinate[midway](aux\X);
				\draw[latex-] (aux\X) -- ++ (-0.5,0) node[above]{$u_\X$};
				\draw[thin,dashed] (m\X.south) -- ++ (0,-0.8) coordinate[pos=0.7](aux'\X);
				\draw[-latex] (aux'\X) -- ++ (1,0) node[midway,below]{$q_\X$};
			}

			\draw[thin,dashed] (m3.south) -- ++ (0,-0.8) coordinate[pos=0.7](aux'3);
			\draw[-latex] (aux'3) -- ++ (1,0) node[midway,below]{$q_{\ell}$};

		\end{tikzpicture}
		\caption{A mass-spring-damper system with $\ell$ masses has the dimension of $n = 2 \ell$.}
		\label{fig:msd}
	\end{figure}

	To preserve the passivity of the full-order model, we let $X$ in \eqref{eq:pinv} be chosen as the Hamiltonian matrix of the system, which satisfies the LMI \eqref{eq:LMI-QSR}. Our optimal projection approach in Algorithm~\ref{alg}, shortened as SPOP, is then compared with the other two projection-based methods, namely, the balanced truncation method for port-Hamiltonian systems (BT-PH) in \cite{kawano2018structure}, and the iterative tangential rational interpolation method (IRKA-PH) for port-Hamiltonian systems in \cite{gugercin2012pH}. In Algorithm~\ref{alg}, we fix the number iteration at $100$ and adopt the steepest descent strategy in \eqref{eq:steepest} to compute a stepsize at each iteration. Furthermore, the starting point of the SPOP algorithm is initialized by one-step interpolation. Specifically,
	we choose sets of initial interpolation points $\{s_1, ..., s_r \}$ that are logarithmically spaced between $10^{-3}$ and $10^{-1}$ and tangential directions $\{b_1 , ... , b_r \}$ as the dominant right singular vectors of the transfer matrix $G(s)$ at each interpolation point. Then an initial point is constructed as follows:
	\begin{equation}
		V_0  = \begin{bmatrix}
			(s_1 I - A)^{-1}B b_1, ... ,(s_r I - A)^{-1}  B b_r
		\end{bmatrix},
	\end{equation}
	which is replaced with a real basis if the columns of $V_0$ occur in conjugate pairs, see \cite{gugercin2012pH} for more details.
	The same interpolation points and directions are used to initialize the IRKA-PH algorithm for comparison.  We adopt each of the methods to reduce the full-order model to reduced ones with the order from $r = 2$ to $r = 20$. The resulting relative $\mathcal{H}_2$ errors for different reduced orders are illustrated in Fig.~\ref{fig:errcomp}, from which we observe that the proposed method achieves superior performance to both BT-PH and IRKA-PH, especially when the reduced model dimension between $6$ and $16$. For instance, when $r = 10$, our method achieves a relative error in the $\mathcal{H}_2$ norm that is $49.55\%$ and $46.96\%$ smaller errors than BT-PH and IRKA-PH, respectively.
	The convergence feature of the presented iterative method is illustrated by Fig.~\ref{fig:converg}, in which we see a monotone decay of the objective function in each case when increasing the number of iterations. Eventually, the algorithm can gradually converge to a local optimum as shown in Fig.~\ref{fig:converg}.

	\begin{figure}[t]
		\centering \includegraphics[width=0.35\textwidth]{fig/errorcomp.pdf}
		\caption{The magnitudes of the relative $\mathcal{H}_2$ errors  as the reduced dimension $r$ varies.}
		\label{fig:errcomp}
	\end{figure}
	\begin{figure}[t]
		\centering \includegraphics[width=0.4\textwidth]{fig/converg.pdf}
		\caption{Comparison of the error magnitudes versus frequency in different methods when $r = 15$.}
		\label{fig:converg}
	\end{figure}

	Next, we illustrate both Algorithm~\ref{alg} and Algorithm~\ref{alg2} in reducing a model with $\ell = 1000$, i.e. $n = 2000$ to a reduced-order model with $r = 10$. In both algorithms, fixed stepsizes are imposed, i.e. $\alpha_k = 0.21$ in Algorithm~\ref{alg}, and $\tau_k = 0.21$ in Algorithm~\ref{alg2} at each iteration. Furthermore, we make a comparison between the two algorithms and the {  Riemannian trust-region method} \cite{sato2018pHsys}, where each algorithm is implemented for 100 iterations with the initial point as the output of IRKA-PH. Fig.~\ref{fig:bode2} depicts the amplitude plot of the reduction errors, and Table~\ref{tab:comp} shows the relative $\mathcal{H}_2$ error and total running time on an Intel Core i7-7500U (2.70GHz) CPU.

	From Fig.~\ref{fig:bode2} and Table~\ref{tab:comp}, we see that Algorithm~\ref{alg} and Algorithm~\ref{alg2} achieve nearly the same approximation error, but  Algorithm~\ref{alg2} requires more computation time since it requires the matrix inversion in \eqref{eq:Vt_new} at each iteration. The Riemannian optimization approach does not impose the Petrov-Galerkin projection and thus achieves a better approximation compared to  Algorithm~\ref{alg} and Algorithm~\ref{alg2}. However, it consumes more than three times the computation time of the proposed algorithms.

	\begin{figure}[t]
		\centering \includegraphics[width=0.4\textwidth]{fig/errorcomp2.pdf}
		\caption{Comparison of the relative error magnitudes versus frequency in the three algorithms.}
		\label{fig:bode2}
	\end{figure}

	\begin{table}[t]
		\centering
		\caption{Comparison of three algorithms with 100 iterations}
		\footnotesize	\begin{tabular}[t]{lcc}
			\toprule
			Algorithms & Relative Error & Running Time \\
			\midrule
			Algorithm~\ref{alg} & 0.1616 & 99.46s \\
			Algorithm~\ref{alg2} & 0.1618 & 109.34s\\
			Method in \cite{sato2018pHsys}  & 0.1088 & 344.81s\\
			\bottomrule
		\end{tabular}
		\label{tab:comp}
	\end{table}%

	\section{Conclusions}
	\label{sec:conclusion}
	In this paper, we have presented a novel model reduction method that incorporates nonconvex optimization with the Petrov-Galerkin projection, which provides a unified framework for reducing linear systems with different properties, including stability, passivity, and bounded realness. An optimization problem has been
	formulated on a noncompact Stiefel manifold, aiming to find an optimal projection such that the resulting approximation error is minimized in terms of the $\mathcal{H}_2$ norm.
	The gradient expression for the objective function is derived, which leads to a gradient descent algorithm that produces a locally optimal solution. The convergence property of the algorithm has also been analyzed. Finally, the feasibility and
	performance of the proposed method has been illustrated by two numerical examples, which shows that the method can be applied to both stability-preserving and passivity-preserving model reduction problems and can achieve smaller approximation errors compared to the existing balanced truncation and Krylov subspace approaches.

	\bibliographystyle{abbrv}

redpajama_set_name: RedPajamaArXiv


text:
Богда́н Евге́ньевич Бутко́ (; ) — украинский футболист, защитник луганской «Зари». Бывший игрок сборной Украины.

Карьера

Клубная
Воспитанник донецкого «Шахтёра». Весной 2008 года был переведён в «Шахтёр-3». Во Второй лиге дебютировал 3 апреля 2008 года в матче против луганского «Коммунальника» (1:2), в этом матче Богдан забил один из голов. В сезон 2008/09 вместе с дублем «Шахтёра» стал победителем молодёжного чемпионата Украины. Летом 2010 года на правах аренды перешёл в луцкую «Волынь», за которую в сезоне 2010/11 провёл 27 матчей и забил 3 мяча в украинской Премьер-лиге, а также 1 матч в Кубке Украины. В июне 2011 года также на правах аренды Бутко перешёл в мариупольский «Ильичёвец», в составе которого также провёл 27 матчей в сезоне 2011/12 (1 гол) и 1 игру в Кубке Украины. Подготовку к сезону 2012/13 после окончания Евро 2012, в котором защитник сыграл в одном из матчей, Бутко начал в составе мариупольского клуба. В январе 2015 года перешёл в пермский «Амкар» на правах аренды до конца сезона. Летом 2015 года «Амкар» договорился с «Шахтером» о продлении аренды Бутко на один сезон. В 2016 году вернулся в родной «Шахтер». После дисквалификации капитана «горняков» Дарио Срны, играл основного правого защитника. В феврале 2020 года на правах аренды перешёл в «Лех» из Познани.

В сборной
В юношеской сборной Украины до 17 лет дебютировал 14 ноября 2006 года в матче против Нидерландов (0:0). Всего за сборную до 17 лет провёл 21 матч.

В юношеской сборной до 19 лет дебютировал 5 января 2009 года в матче против Бельгии (0:0). Бутко был включён Юрием Калитвинцевым в состав сборной на чемпионат Европы 2009 года в Донецке и Мариуполе. Сыграл всего 2 матча на групповой стадии против Словении (0:0) и Англии (2:2), из-за полученной травмы. В финале Украина обыграла Англию (0:2).

В национальной сборной первый матч сыграл в 20 лет, 2 сентября 2011 года, отыграв все 90 минут в спарринге с Уругваем, пятой командой мира по рейтингу на тот момент. Его игру прокомментировал главный тренер сборной Олег Блохин: «Очень понравился дебют Бутко, он, по-моему, меньше меня волновался. Вот на таких парнях и будем строить игру».

В мае 2016 года главный тренер сборной Украины Михаил Фоменко включил Бутко в резервную заявку на чемпионат Европы 2016.

Всего за сборную Украины сыграл 33 матча.

Достижения
 Чемпион Европы среди юношей (до 19): 2009
 Чемпион Украины (3): 2016/17, 2017/18, 2018/19
 Обладатель Кубка Украины (3): 2016/17, 2017/18, 2018/19

Примечания

Ссылки

Футболисты Украины
Игроки сборной Украины по футболу (до 21 года)
Игроки сборной Украины по футболу
Чемпионы Европы по футболу среди игроков до 19 лет
Игроки ФК «Шахтёр-3»
Игроки ФК «Шахтёр» Донецк
Игроки ФК «Волынь»
Игроки ФК «Мариуполь»
Игроки ФК «Амкар»
Игроки ФК «Лех»
Игроки ФК «ББ Эрзурумспор»
Игроки ФК «Заря» Луганск

redpajama_set_name: RedPajamaWikipedia


NAN(Lagos): The Joint Admissions and Matriculation Board (JAMB) says the recent lists forwarded to various tertiary institutions was for candidates that qualified for screening based on institutions' capacity.
It added that if, however, the number of candidates on such lists was not sufficient for existing vacancies, the institutions could source from the omnibus printout earlier sent to them by the board.
"The public and all tertiary institutions should note that admission will only be approved by the board after appropriate screening of the candidates by the institutions.
"The list which is made up of candidates who met the national cut-off point within the set criteria is then forwarded to the institutions,'' it said.
The statement noted that the intention of the board was to ensure that available spaces were adequately utilised.


For ages, local rappers tore each other down, as though the city could produce only one star at a time—but in 2016 the whole scene seemed to grasp the value of community.
"I sometimes crave being anonymous, this guy, just playing jazz in some dingy club somewhere."
The Hood Internet look to dominate the Real Internet: their Album Tacos blog is up for an Urlie.
The GZA, the Hood Internet, Dan Deacon, and a funk musician walk into a bar . . .


text:
Thirty On Tap
Tag: movies
15 Fall Activities For Every Mood
October 22, 2017 September 25, 2019 2 Comments
By De and Kate
While there's not a whole lot we can all agree on in 2017, one thing is for certain: Fall is the best season. Between its crisp air, changing colors, delicious flavors, and plentiful activities, it's got a little something for everyone. If you ask us, it's singular biggest downfall is how fast it always seems to go. We hardly feel like we have the time to soak in the glory of the season before it's over and winter has arrived. To express our appreciation for this world in which there are Octobers, we've compiled a list of all the ways to make the most of fall this year. Continue reading →
5 Lady-Boss Inspirations to Take Away From "Girls Trip"
July 28, 2017 July 26, 2017 2 Comments
By Eliza David
Like most women in America, I looked forward to the premiere weekend of the comedy blockbuster, Girls Trip. This was not only a movie centering around black women with a message of friendship, but it was filled with all the naughty comedy a romance writer like me can appreciate. I went to see the film on its opening Friday night and walked away with a headache from laughing so much. Girls Trip boasts an ensemble of talented ladies, great writing, raucous humor, and beautiful men. It was beyond awesome, everything I wanted in a comedy – period. I hadn't laughed that hard at a comedy since Amy Poehler/Tina Fey's Sisters, but the fabulous foursome of Queen Latifah, Jada Pinkett-Smith, Regina Hall, and Tiffany Haddish are on higher plane than any comedy ensemble I've seen in years.
But, beyond the laughs, your resident #GirlBoss Guru walked away with some reminders that soothed both my business mind and busy soul. Continue reading →
The Problem With Netflix's New Movie, 'Death Note'
March 27, 2017 March 27, 2017 1 Comment
By Megan Kramer
I am 100% what present-day society would refer to as a "binge watcher." I have been known to watch full seasons of TV shows in one day, or one weekend if the season is particularly long. And while my binge watching isn't limited to Netflix, most of it occurs there. The streaming service not only has some of my favorite older TV shows to binge, but also new, original programs that have been wonderful: Marvel's Jessica Jones, A Series of Unfortunate Events, Stranger Things and Master of None, to name a few (because trust that I could go on).
But I kind of feel like Netflix is totally killing me.
You may be thinking, "Well, binge watching definitely contributes to a sedentary lifestyle, but it's probably illegal to say Netflix itself is killing you." And you'd be right. But that's not where I'm going with this. Continue reading →
10 Clichés That Remind Us of Disney Movies
November 3, 2016 October 30, 2016
Normally, we tend to avoid clichés like the plague. While they were once clever and meaningful, these overused phrases have lost their punch, and end up making our writing or speech sound cheesy. But in honor of Cliché Day, let us all reflect on the fact that many clichés ring true, and can be found in none other than our favorite Disney movies. Continue reading →
16 Happy Things That Happened in 2016
It's easy to get caught up in a downward spiral of thinking that 2016 has been one terrible, horrible, no good, very bad year. And for a lot of reasons, it has. It's understandable that we'd feel exhausted, both mentally and physically, when we look back at all of the depressing things that have taken place, so many of which has caused us to truly lose faith in humanity.
However, it's important to keep in mind that it's not all doom and gloom out there – even though it sure feels that way sometimes. In an effort to remind everyone – and ourselves – of the good in the world, here are some of the best things that happened this year, amid all the other darkness of 2k16. Continue reading →
Read This, Watch That: Our Favorite Books To Hit The Big Screen
October 6, 2016 October 5, 2016 1 Comment
By De & Kate
Last year, we were quickly captivated by The Girl On The Train and flew through its suspenseful, borderline frustrating, storyline. We often found ourselves wishing that we could hop into the pages of the novel, take a seat next to Rachel, and solve the damn mystery already. Continue reading →
5 Reasons To Be Pumped For "Suicide Squad"
If you go based off of the trailer you can tell that Suicide Squad is a movie not only meant to have the villains portrayed as the heroes but to have them become just as popular as the heroes. DC might be able to pull such a feat off thanks to the characters they've chosen to highlight as well as the actors chosen to portray them. If you aren't all ready hyped here are some reasons you should be. Continue reading →
5 Reasons "Bad" Moms Are Winning
Women are bombarded with images of motherhood and it all starts during childhood. From the very first stuffed animal tucked into our beds at night to watching our own moms make the impossible happen, women are taught that mothers "do it all" with little to no effort. It's easy, flawless, and natural, we're told. As we become adults, some of us become conditioned to this ideal. Continue reading →
#TBT: Don't Tell Mom – 90's Shows & Movies We Watched On The DL
Picture this – you're in your bedroom, or at a friend's house at a sleepover, watching that show you know your mom would kill you for watching. You've got the volume down low, and you're ready to jump into action at the slightest indication of footsteps in the hall. You know the drill. We dedicate this piece to all the 30-year-old women who can still remember the thrill – and slight panic – of watching these shows and movies with remote in hand, hovering the 'last' button at all times. Continue reading →
On His 40th Birthday: 7 Reasons To Love Colin Farrell
May 31, 2016 May 30, 2016 1 Comment
By Jaclyn Lindsay
Today marks Colin Farrell's 40th birthday. As a fellow Irish person – and self-confessed fan of his – I thought it was only fair to pay tribute to him.
The Dublin-born actor has been no stranger to the headlines in the past, with his various addictions to drugs, booze, and women. But in recent years, he has nearly disappeared from gossip columns completely, opting instead for a more private life, all the while, still maintaining his reputation as a respected actor.
View thirtyontap's profile on Facebook
View thirtyontap's profile on Twitter
View thirtyontap's profile on Instagram
View thirtyontap's profile on Pinterest
Subscribe to Thirty on Tap via email
Enter your email address to follow Thirty On Tap and receive notifications of new articles by email.
Crying Over Cake Fails
I'm A Mom. And I'm More.
How Do I Teach My Daughter To Love Her Body When I'm Struggling To Do The Same?
Fill Up. Empty Out.
Only And Already
The Art and Sadness of "Getting Through"
You're Allowed to Be Sad
ABOUT * CONTACT * SUBMISSIONS
© 30ontap.com, 2019. Unauthorized use and/or duplication of this material without express and written permission from this site's author and/or owner is strictly prohibited. Excerpts and links may be used, provided that full and clear credit is given to 30ontap.com with appropriate and specific direction to the original content.

redpajama_set_name: RedPajamaCommonCrawl


Quit calling this an illustration.
Thanks for the strong reaction to "Playing around online" Anonymous but I'm sticking to my nomenclature. I obviously adhere to a less strict definition of the word "illustration" than you.


Smoke, heat and flames represent major fire risks to building users. How materials behave in a fire - from ignition to flashover and a fully developed fire - therefore needs to be optimised at the design stage with correct HVAC insulation materials.
Whether for schools, hospitals, shopping malls, offices or residentials buildings, ISOVER has developed a range of HVAC insulation solutions that provide efficient insulation for your pipes and ductwork. Using ISOVER's HVAC insulation materials will help ensure your building is thermally protected with minimum noise disruption and maximum fire protection.
ULTIMATE™ U PROTECT is a new-generation mineral wool insulation that takes HVAC fire protection to a new level. ISOVER's ULTIMATE™ HVAC insulation uses breakthrough technology to combine the performance and handling benefits of glass wool with the outstanding fire resistance of stone wool.
ULTIMATE™U PROTECT HVAC insulation materials are non-combustible, do not contribute to ignition of fire and do not emit any smoke and burning droplets. They can resist fire up to 2 hours securing your buildings' fire safety in maximum level.
Designed for passive fire protection solutions in HVAC systems, ULTIMATE™ U PROTECT HVAC insulation materials cover different types of pipes, rectangular and circular ventilation and smoke extraction ducts for a complete solution offering excellent, lightweight thermal and acoustic properties, as well as low maintenance costs.
The ULTIMATE™ U PROTECT range provides all the advantages of conventional thermal and acoustic insulation with increased fire protection.
Guaranteeing superior material stability and thermal resistance at high temperatures, the HVAC insulation materials comprise slabs, wired mats and pipe sections that ensure non-combustibility and up to 2 hours of fire protection for pipes and both rectangular and circular ductwork systems.
More information about Ultimate™ U Protect HVAC insulation products and other materials? Contact us!


Looking for Sleep Apnea in Bullhead City? For an appointment at BullheadCityDentist.com call our office today at 928-758-0008 or make an appointment today.
Your dentist will perform an intraoral exam, which involves specialized x-rays, intraoral photos, gag reflex check, orthopedic exams, and the Epworth Sleepiness Scale. These tests will help us determine the most effective position for your lower jaw to open your airway. BullheadCityDentist.com has had success treating sleep apnea with a nighttime oral appliance that is more comfortable than a CPAP mask. This appliance will be customized especially for you. Depending on the severity of your condition, other repositioning devices may be considered.
Call us at 928-758-0008 for a sleep apnea assessment.


Our scientific research aims to understand the development of biodiversity and the interaction between species. The natural history collection of Naturalis Biodiversity Center is the 5th largest in the world with 37 million objects. The arachnid collection at Naturalis contains 13,000 spider specimens and over 120,000 mites. It includes important mite collections Oudemans, Lukoschus, and van der Hammen, and is strong in spiders of Southeast Asia, particularly the Deeleman collection. Please contact Karen van Dorp to arrange loans or visits to the collection.
The fascination of Nature – at the State Museum of Natural History Karlsruhe Nature becomes a great experience for the entire family! Permanent and special exhibitions, the vivarium with its exotic animals, as well as a rich choice of events for both adults and kids make the State Museum of Natural History Karlsruhe a place to discover the diversity and beauty of Nature.