Stability-preserving model order reduction for linear stochastic Galerkin systems

Mathematical modeling often yields linear dynamical systems in science and engineering. We change physical parameters of the system into random variables to perform an uncertainty quantification. The stochastic Galerkin method yields a larger linear dynamical system, whose solution represents an approximation of random processes. A model order reduction (MOR) of the Galerkin system is advantageous due to the high dimensionality. However, asymptotic stability may be lost in some MOR techniques. In Galerkin-type MOR methods, the stability can be guaranteed by a transformation to a dissipative form. Either the original dynamical system or the stochastic Galerkin system can be transformed. We investigate the two variants of this stability-preserving approach. Both techniques are feasible, while featuring different properties in numerical methods. Results of numerical computations are demonstrated for two test examples modeling a mechanical application and an electric circuit, respectively.


Introduction
Numerical simulation of mathematical models represents the main issue in scientific computing.We consider linear dynamical systems, which play an important role in mechanics and electrical engineering, for example.Furthermore, uncertainty quantification becomes more and more relevant in many fields of applications, see [14], for instance.A common approach is to replace uncertain parameters by random variables, see [32,34].Statistics of the stochastic model can be computed by sampling methods or quadrature rules.Alternatively, the stochastic Galerkin method changes the random-dependent linear dynamical system into a larger deterministic linear dynamical system.
The dimension of the stochastic Galerkin system becomes huge in the case of large numbers of random variables.Methods of model order reduction (MOR) are able to decrease the complexity.Transient solutions of a reduced system allow for an efficient numerical simulation.Several MOR methods are available for general linear dynamical systems, see [1,2,3,29].MOR of linear stochastic Galerkin systems was also examined in several previous works [16,28,23,24,35].
However, the asymptotic stability is often lost in the reduced Galerkin system for some MOR techniques.We investigate stability-preserving strategies in the case of Galerkin-type projection-based MOR like the Arnoldi method or proper orthogonal decomposition, for example.A dissipativity property guarantees the preservation of stability.If a general linear dynamical system does not satisfy the dissipativity property, then it can be transformed into a dissipative structure, see [6,21,26].The crucial part to identify a transformation consists in the solution of a Lyapunov equation.Direct methods, see [9], or approximate methods, see [19,20,30,33], yield the numerical solutions of Lyapunov equations.
We examine the stability-preserving approach in the case of linear stochastic Galerkin systems consisting of ordinary differential equations.Several variants are feasible.The high-dimensional Galerkin-projected system is transformed or, vice versa, the original systems are transformed followed by a Galerkin projection.We analyze the two strategies and another variant.
In addition, network approaches produce models consisting of differential-algebraic equations in industrial applications.Thus we extend the stability-preserving techniques to this class of problems.The Lyapunov equations have no solution now.Therefore, we use a regularization technique, which was also employed in [17].
We apply the analyzed techniques to mathematical models of two test examples: a mass-spring-damper system and an electric circuit of a band-pass filter.

Stability preservation in reduction
We review a concept for stability preservation in Galerkin-type projection-based MOR for general linear dynamical systems.

Linear dynamical systems
We consider linear dynamical systems in the form E ẋ(t) = Ax(t) + Bu(t) y(t) = Cx(t). ( with constant matrices A, E ∈ Ê n×n , B ∈ Ê n×n in , and C ∈ Ê nout×n The state variables or inner variables are x : [0, t end ] → Ê n .Inputs u : [0, t end ] → Ê n in are supplied to the system.The outputs y : [0, t end ] → Ê nout are defined as quantities of interest (QoI).Initial value problems are given by x(0) = x 0 .
If the mass matrix E is non-singular, then the system (1) consists of ordinary differential equations (ODEs).If the mass matrix E is singular, then the system (1) represents differential-algebraic equations (DAEs).Furthermore, we assume that the system satisfies the following stability condition.
The maximum real part of these eigenvalues is called the spectral abscissa.In the case of DAEs, the asymptotic stability implies that the matrix pencil (E, A) is regular.

Projection-based model order reduction
Projection matrices V, W ∈ Ê n×r of full rank are specified with r ≪ n.Concern- ing the full-order model (FOM) in (1), the reduced-order model (ROM) reads as Ē ẋ(t) = Āx(t) + Bu(t) with state variables or inner variables x : [0, t end ] → Ê r .Initial values x(0) = x0 are supposed.We obtain the matrices via which is also called a Petrov-Galerkin-type MOR.A Galerkin-type projectionbased MOR is characterized by W = V , where just one projection matrix has to be determined.Important examples are the one-sided Arnoldi method and the proper orthogonal decomposition (POD), see [1].
Each linear dynamical system is described by a transfer function in the frequency domain.The difference between the transfer functions of FOM and ROM quantifies the error of the MOR.Hardy norms like the H 2 -norm, for example, can be applied.However, a small error in the frequency domain implies a small error in the time domain only if both systems are asymptotically stable.

Dissipative systems
In balanced truncation, see [1], the ROM ( 2) is always asymptotically stable provided that the system (1) is asymptotically stable.Yet the asymptotic stability may be lost in the ROM (2) within other MOR methods like Krylov subspace techniques, see [8], and POD, for example.
Using Galerkin-type MOR, the stability is guaranteed for some classes of linear dynamical systems.In the case of ODEs, we define the following type of system.
Definition 2. A linear dynamical system (1) is called dissipative, if 1. E is symmetric as well as positive definite, and 2. A + A ⊤ is negative definite.
The above condition represents a dissipativity of the matrix A as shown in [18].Other definitions of dissipative systems are used in the literature.If a linear dynamical system (1) is dissipative with respect to Definition 2, then it is asymptotically stable as in Definition 1.In contrast, the asymptotic stability does not imply the dissipativity of Definition 2, even if the mass matrix is symmetric and positive definite.Now we consider Galerkin-type MOR.
Theorem 1.If the linear dynamical system (1) is dissipative, then a Galerkintype MOR yields a dissipative reduced system (2).Hence the reduced system is asymptotically stable.
The proof can be found in [26], for example.

Transformations
If a system of ODEs is not dissipative, then it can be converted to an equivalent dissipative form by a basis transformation in the state space, see [21].Alternatively, a basis transformation is feasible in the image space only, see [6].We require a symmetric positive definite solution M ∈ Ê n×n of the Lyapunov in- which means that the matrix on the left-hand side of ( 4) is negative definite.
If the mass matrix E is non-singular and the matrix pencil (E, A) satisfies the condition of asymptotic stability, then an infinite set of solutions exists.
We change the system of ODEs (1) into the equivalent system The transformed system exhibits the desired property.
Theorem 2. If the asymptotically stable linear dynamical system (1) has a nonsingular mass matrix, then the transformed system (5) is dissipative in view of Definition 2.
We use a Galerkin-type MOR with a projection matrix V to the transformed system (5).This approach can be written as a Petrov-Galerkin-type MOR applied to the original system (1) with matrices (3) and the projection matrix Thus we do not require to calculate the transformed system (5) explicitly.Alternatively, we compute the projection matrix (6).However, the original Galerkintype MOR of the system (1) is not equivalent to the MOR of the system (5).

Numerical solution of Lyapunov inequality
We solve the Lyapunov inequality (4) using a Lyapunov equation including a predetermined symmetric positive definite matrix F ∈ Ê n×n .This matrix represents a degree of freedom, because any choice yields a symmetric positive definite solution of the Lyapunov inequality.A simple admissible choice is the identity matrix I ∈ Ê n×n .Moreover, we do not need to solve the Lyapunov equation ( 7) with a high accuracy, because a rough approximation M often still satisfies the Lyapunov inequality (4).
There are direct methods to compute a solution M of (7) or a symmetric decomposition M = LL ⊤ , see [9,19].Their computational effort is typically O(n 3 ).In the high-dimensional case, we have to use approximate methods to decrease the computation work.The following techniques are available: i) projection methods (Krylov subspace techniques, POD, etc.), see [11,33], ii) alternating direction implicit (ADI) iteration, see [13,20], iii) frequency domain integrals, see [4,25], and others.In the cases (i) and (ii), the methods yield an approximation M = ZZ ⊤ with a low-rank factor Z ∈ Ê n×k (k ≪ n).Thus the transformation is given by a singular matrix M .It follows that the mass matrix Ē of the reduced system (2) may become singular or ill-conditioned, as shown in [26].In contrast, the method (iii) from [25] computes the projection matrix (6), where the underlying approximation M is always non-singular.

Stochastic Galerkin systems
We illustrate the concept of the stochastic Galerkin method and define our problem under investigation.

Random linear dynamical systems
We include parameters in the linear dynamical systems and obtain The matrices A, E ∈ Ê n×n , B ∈ Ê n×n in , and C ∈ Ê nout×n depend on parameters µ ∈ M ⊂ Ê q .Thus the state variables or inner variables x : [0, t end ] × M → Ê n depend on time as well as the parameters.The inputs u : [0, t end ] → Ê n in are independent of the parameters, whereas the outputs y : [0, t end ] × M → Ê nout vary with respect to the parameters.We consider a single output (n out = 1) without loss of generality.
We assume that the system ( 8) is either an ODE for all µ ∈ M or a DAE for all µ ∈ M. Let the system be asymptotically stable for each parameter with respect to Definition 1. Furthermore, initial value problems are considered including a function x 0 : M → Ê n .The initial values may be independent of the parameters.The initial values have to be consistent in the case of DAEs.
We suppose that the parameters in the linear dynamical system (1) are affected by uncertainties.A common approach is to substitute the parameters by independent random variables µ : Ω → M on a probability space (Ω, F , P ) with event space Ω, sigma-algebra F and probability measure P .We apply traditional probability distributions like uniform, Gaussian, beta, etc. Hence a joint probability density function ρ : M → Ê is available.This approach yields a stochastic model.
A measurable function f : M → Ê depending on the random variables exhibits the expected value (10) provided that the integral is finite.The expected value (10) implies the inner product for two functions in the Hilbert space The associated norm is f L 2 (Π,ρ) = f, f as usual.

Polynomial chaos expansions
In most cases, a complete orthogonal basis (Φ i ) i∈AE of polynomials Φ i : M → Ê exists, see [34].Let this basis also be normalized.We assume that the QoI of the system ( 8) is in the space (12) for each time point.Consequently, the QoI can be expanded into the series with coefficient functions w i : [0, t end ] → Ê, which is called a (generalized) polynomial chaos expansion (PCE).The series (13) converges in the norm of ( 12) pointwise in time.
Likewise, the state variables exhibit the PCE with coefficient functions v i : [0, t end ] → Ê n , provided that each state variable is in the Hilbert space (12).
We assume that the first basis polynomial is the unique constant polynomial Φ 1 ≡ 1.The orthonormality Φ i , Φ j = δ ij of the basis functions imply the formulas for the expected value and the variance of the QoI in each time point.

Stochastic Galerkin method
We truncate the series ( 13), (14) to the finite sums Typically, all basis polynomials up to some total degree d are included.The number of basis polynomials up to total degree d reads as If the number of random variables is large, then the number of basis polynomials becomes huge even for moderate degrees, say d = 3.
Inserting the approximations (15) into the linear dynamical system (8) yields a residual.The Galerkin approach requires that the residual is orthogonal to the space spanned by the basis polynomials Φ 1 , . . ., Φ m with respect to the inner product (11).Basic calculations produce a larger coupled linear dynamical system for v = (v ⊤ 1 , . . ., v⊤ m ) ⊤ and ŵ = ( ŵ1 , . . ., ŵm ) ⊤ .Hence we obtain a linear dynamical system of dimension n = mn with m outputs.The matrices exhibit the sizes A ∈ Ê n×n , B ∈ Ê n×n in , and C ∈ Ê m×n .To define the matrices, we introduce the auxiliary arrays It follows that using Kronecker products.Therein, the probabilistic integration (10) operates separately in each component of the matrices.Often it holds that rank( Ĉ) = m provided that m ≤ n.More details on the stochastic Galerkin method for linear dynamical systems are given in [22,28].
The matrix Ê may be singular even though all matrices E are non-singular due to properties of the spectrum, see [31].However, this loss of invertibility hardly occurs.We assume that Ê is always non-singular in the case of ODEs (8).Likewise, the stochastic Galerkin system ( 16) may be unstable even though the systems ( 8) are asymptotically stable for all parameters.Academic examples are given in [27].Again this loss of stability is hardy observed in practice.
Initial conditions v(0) = v0 are derived from the initial values (9) of the original dynamical system (1).If the initial values (9) are identical to zero, then the choice v(0) = 0 is obvious.The approximation of the QoI reads as where the outputs of the stochastic Galerkin system (16) yield the coefficients.
We prove a property, which will be used later.
Lemma 1.If the matrices E(µ) are symmetric and positive definite for almost all µ ∈ M, then the stochastic Galerkin projection Ê is also symmetric and positive definite.

Proof:
The Galerkin-projected matrix consists of the minors Êij = E[Φ i Φ j E] for i, j = 1, . . ., m, see (18).Hence the symmetry is obvious.We obtain because the integrand is almost everywhere non-negative in the probabilistic integration (10).The basis functions (Φ i ) i∈AE are linearly independent.Thus z = 0 implies that the above finite sum is non-zero on a subset U ⊂ M satisfying P ({ω : µ(ω) ∈ U}) > 0 for the probability measure P .It follows that z ⊤ Êz > 0 for z = 0.
Likewise, this relation applies to the case of negative definite matrices.

Stability preservation
We investigate three strategies to preserve the asymptotic stability in an MOR of a stochastic Galerkin system.Figure 1 illustrates different possibilities.

Transformation of stochastic Galerkin system
If the linear dynamical systems (8) are asymptotically stable for all parameters, then the stochastic Galerkin system ( 16) is usually asymptotically stable.However, if the system ( 16) is not dissipative, then stability may be lost in some MOR methods.In this case, we can transform the system to a dissipative form as demonstrated in Section 2.4 and Section 2.5.Consequently, the reduced system is stable due to Theorem 1. Remark that we do not have to calculate the transformed matrices of the high-dimensional system explicitly.Just an appropriate projection matrix has to be determined via (6).
In this approach, the critical part is the solution of the high-dimensional Lyapunov equation (7).A direct method of linear algebra would require O(m 3 n 3 ) operations.Thus we are restricted to approximate methods or iteration schemes.
The possibilities are listed in Section 2.5.

Transformation of parameter-dependent system
Now we transform the linear dynamical systems (8) first.

Transformation
The stochastic Galerkin system is not always asymptotically stable, even if all parameter-dependent systems (8) are asymptotically stable.However, we obtain a positive result in the case of dissipativity.
Theorem 3. If the linear dynamical systems (8) are dissipative for almost all µ ∈ M, then the stochastic Galerkin system ( 16) is also dissipative.Consequently, a Galerkin-type reduction of ( 16) yields an asymptotically stable system. Proof: Lemma 1 implies that the mass matrix Ê is symmetric and positive definite.The matrix Â + Â⊤ is the Galerkin projection of see (18).Since A(µ) + A(µ) ⊤ is negative definite for almost all µ ∈ M, Lemma 1 shows that Â + Â⊤ is negative definite.Hence both conditions of Definition 2 are satisfied.Theorem 1 guarantees the stability of a reduced system.
If the original systems (1) are not dissipative for almost all µ ∈ M, then we transform the systems appropriately.The following transformation has to be applied for almost all µ simultaneously (even if some system is already dissipative) to preserve continuity and smoothness of the functions.
We obtain the parameter-dependent Lyapunov equation for µ ∈ M.Although the matrix F may depend on the parameters, we choose a constant matrix, because no parameter-aware strategy with benefits is known yet.The Lyapunov equations ( 20) yield a unique family M(µ) of symmetric positive definite matrices.

Polynomial system matrices
Now we assume that the matrices A(µ), B(µ), E(µ) of the system (8) involve only polynomials in the variable µ, which is often given in practice.Thus the matrices Â, B, Ê of the stochastic Galerkin system ( 16) can be calculated analytically in the case of traditional probability distributions.We do not consider the output matrix C(µ), because it is not transformed.However, the matrix-valued function M(µ) satisfying (20) consists of rational functions in the variable µ.In the case of low dimensions n, the function M(µ) can be calculated explicitly by a computer algebra software.Alternatively, the solution M(µ) can be evaluated for a finite set of parameters µ.
We require the Galerkin projection of the transformed matrices in the dissipative system (5).The entries of these transformed matrices are rational functions of µ.A quadrature rule yields numerical approximations Ã, B, Ẽ of the matrices Â′ , B′ , Ĉ′ .Let {µ 1 , . . ., µ k } ⊂ M be the nodes and {γ 1 , . . ., γ k } ⊂ Ê be the weights.The approximation of the minors in Â′ reads as for i, j = 1, . . ., m.Using the auxiliary array (17), the approximation becomes Likewise, the quadrature yields B and Ẽ.The computational effort is dominated by k numerical solutions of the Lyapunov equations (20).

Proof:
The symmetry of Ẽ is obvious.Let z ∈ Ê mn \{0} with z = (z ⊤ 1 , . . ., z ⊤ m ) ⊤ .We obtain using the notation ( 22) Since E ′ (µ ℓ ) is positive definite and the weights γ ℓ are positive for each ℓ, all terms of the sum are non-negative.Hence Ẽ is positive semi-definite.The negative semi-definiteness of Ã + Ã⊤ can be concluded by the same treatment.
In addition, it is very likely that one or more terms are positive in the above sum.Thus we assume that these matrices are definite.It follows that the approximate Galerkin system has the dissipativity condition of Definition 2.

Properties
The strategy of this section looks appealing in the case of low dimensions n, because the Lyapunov equations can be solved cheap.However, the number of random parameters is typically large in this case, because otherwise the stochastic Galerkin system is not high-dimensional and an MOR is obsolete.A large number of random variables implies a quadrature in a high-dimensional space.
Using a reasonable number of nodes, the quadrature error may still be too large such that the exact solution of the approximate Galerkin system yields poor approximations ( 19) of the random-dependent QoI.Therefore, the above approach is critical.
Furthermore, there is a major loss of sparsity in this technique.Let the entries in the matrices be polynomials of low degrees depending on µ in the system (8).
Even if the matrices are dense, then the stochastic Galerkin system ( 16) exhibits sparse matrices due to the orthogonality of the basis polynomials.However, the dense matrices of the transformed systems ( 5) are rational functions depending on µ.Inner products (11) of their entries and orthogonal polynomials are nonzero and thus the matrices of ( 16) become dense in the Galerkin projection.In contrast to the approach of Section 4.1, we have to calculate the matrices of the alternative Galerkin system explicitly to determine a projection matrix V of the MOR.

Transformation using reference parameter
We derive an additional technique to decrease the computational effort and to omit quadrature errors.A single reference parameter µ * ∈ M is selected like the mean value µ * = E[µ], for example.We directly solve the Lyapunov equation (20) only for µ * using some matrix F .The solution M * = M(µ * ) ∈ Ê n×n is symmetric and positive definite.We define the larger transformation matrix using the identity matrix I m ∈ Ê m×m and the Kronecker product.Obviously, the matrix ( 23) is symmetric and positive definite again.We employ this matrix to transform the stochastic Galerkin system ( 16) into the form (5). Since the matrix ( 23) is block-diagonal, matrix-matrix multiplications with M are cheap.We do not need to compute the transformed system matrices.Alternatively, we directly compute the projection matrix (6), where a projection matrix V is determined by the original Galerkin system.
Theorem 4. Let the random variables in a linear system of ODEs (8) be of the form (24).There is a positive constant θ 0 ∈ (0, 1] such that the transformation of a stochastic Galerkin system (16) using the matrix ( 23) yields a dissipative system of type ( 5) for all θ ∈ [0, θ 0 ].

Proof:
The transformed mass matrix Ê⊤ M Ê is symmetric and positive definite for all θ, since just the non-singularity of the original mass matrix Ê is required.In the limit, we obtain lim due to the orthogonality of the basis polynomials.Hence the matrix becomes block-diagonal with identical blocks.The factor M * satisfies the Lyapunov equation (20) for the chosen positive definite matrix F .The symmetric part of the matrix ( 25) is negative definite, because it holds that The two conditions of Definition 2 are satisfied and thus the system is dissipative in the limit.Since the limit ( 25) is reached continuously with respect to the parameter θ, it follows that a sufficiently small perturbation of θ = 0 still yields matrices with the required properties.
Theorem 4 shows that the stochastic Galerkin system becomes dissipative for all sufficiently small θ ≥ 0. However, this implication does not guarantee that the system is dissipative for θ = 1, which reproduces our desired choice of random variables in (24).Nevertheless, the computational effort is low such that it is worth to try this approach.Even if this transformed stochastic Galerkin system is not dissipative, a loss of stability may happen less often in an MOR.

Differential-algebraic equations
If the linear dynamical system (1) features a singular mass matrix E, then a system of DAEs is given.MOR methods are also available for DAEs, see [5].
Yet the Lyapunov equations ( 7) do not have a solution, which is crucial in our stability-preserving technique.Let x ∈ Ê n .It follows that Choosing an x = 0 in the kernel of the matrix E implies Ex = 0 and thus a contradiction to the definiteness of the matrix F appears.We require an alternative strategy now.

Regularization
We apply a regularization of an asymptotically stable DAE system (1), which was also used in [17].The regularized system matrices read as introducing parameters α, β > 0. The matrix E reg is non-singular for all α > 0. Choosing α = β 2 , it follows that the linear dynamical system (1) with A reg , E reg is asymptotically stable for all sufficiently small parameters β.An advantage of this regularization technique is that the sparsity pattern of sE − A coincides with the sparsity pattern of sE reg − A reg for each s ∈ .Hence no loss of sparsity happens in the relevant operations.
We also choose a small parameter β (and α = β 2 ) to ensure that the difference between the original DAE system and the regularized ODE system is small.Assuming that the DAE system together with the defined outputs is strictly proper, error bounds were derived with respect to the H 2 -norm in [25].

Stochastic Galerkin projection
The following two approaches are equivalent provided that the same parameters α, β are chosen in a regularization (26): i) Regularize the parameter-dependent system (8) and then project the systems of ODEs in the stochastic Galerkin method.
ii) Project the parameter-dependent DAE system (8) in the stochastic Galerkin method and then regularize the Galerkin system (16).
In both cases, we obtain the same matrices Â and Ê in the stochastic Galerkin system (16), where the mass matrix Ê is non-singular.
Concerning the flowchart in Figure 1, the transformation to a dissipative representation is done only for a regularized system, since a system of ODEs is required for each Lyapunov equation.Now the two approaches are not equivalent any more.The stochastic Galerkin method is invariant with respect to transformations by constant non-singular matrices.However, we use the parameter-dependent matrix E(µ) ⊤ M(µ) in the transformation to the dissipative form (5) within the strategy (i).Furthermore, the drawbacks of the succession (i), see Section 4.2.3, also apply in this case.
In the approach (ii), the approximate solution of the high-dimensional Lyapunov equations ( 7) may become critical.Small regularization parameters α nearly coincide with the singular case.This problem is less pronounced in direct methods to solve the Lyapunov equations.
We can also employ the technique from Section 4.3 in this context.The linear dynamical system ( 8) is regularized for a reference parameter µ * and some α, β > 0. Solving the Lyapunov equation ( 20) including µ * yields the highdimensional transformation matrix (23).Now the regularized Galerkin system y u from strategy (ii), where the same values α, β are used, is transformed based on (23).Again the dissipativity of the transformed representation cannot be guaranteed a priori.

Illustrative examples
We investigate a system of ODEs as well as a system of DAEs in this section.The computations were executed within the software package MATLAB [15].

Mass-spring-damper system
Figure 2 depicts a mechanical configuration, which consists of 5 masses, 7 springs and 5 dampers.The single input is the excitation at the bottom spring, whereas the single output is the position of the top mass.A mathematical modeling yields a linear dynamical system (1) of n = 10 ODEs of first order including q = 17 parameters.The system matrices are affine-linear functions (polynomials of degree one) of the parameters.The Bode plot of the system is shown for a specific choice of the parameters in Figure 3.This system represents a modification of a test example used in [12].In the stochastic modeling, we replace the parameters by independent uniformly distributed random variables, which vary 10% around their mean values.The mean values are the constant parameters from above.In the truncated orthogonal expansion (15), we include all multivariate Legendre polynomials up to degree three.We obtain m = 1140 basis functions.The stochastic Galerkin system ( 16) is calculated exactly except for round-off errors.Table 1 illustrates the properties of this example.The system is asymptotically stable.The mass matrix is symmetric and positive definite.Yet the system is not dissipative.
We employ the one-sided Arnoldi method, see [1], to perform an MOR of the stochastic Galerkin system (16).The single expansion point s = 0. on a reference parameter (Section 4.3).We choose the identity matrix (F = I) as degree of freedom in each Lyapunov equation.In all three approaches, the stability is achieved for each ROM. Figure 4 illustrates the relative errors with respect to the H 2 -norm, i.e., the difference between FOM and ROM.
In the technique (a), we use the method from [25], where an integral is discretized by a quadrature rule in the frequency domain.Therein, the solution of the highdimensional Lyapunov equation is not computed but the associated projection matrix (6).The computation work is characterized by an LU-decomposition of the system matrix iω Ê − Â in each node ω.We employ the Gauss-Legendre quadrature.Table 2 shows the number of stable ROMs for different numbers of nodes.We observe that 40 nodes are sufficient to stabilize all reduced systems.Furthermore, the error of the MOR exhibits the same magnitude as in the original stochastic Galerkin system.In the technique (b), we use a sparse grid quadrature of Smolyak-type with level 3 based on the Clenshaw-Curtis rule.The scheme exhibits 7209 nodes in the 17dimensional space and negative weights arise.Table 3 illustrates the loss of sparsity in this approach.The relative errors of the MOR are demonstrated in Figure 4 (b).The inherent errors between the ROMs and the stochastic Galerkin projection of the transformed systems decay for increasing dimensions.However, the errors between these ROMs and the stochastic Galerkin system of the untransformed systems (8) are large and stagnate.This property indicates that the quadrature error is too large, even though a high number of nodes is used, because the original Galerkin system can be considered sufficiently accurate.
In the technique (c), we define the mean value of the random variables as reference parameter.We calculate the matrices of the transformed Galerkin system to analyze the definiteness.Figure 5 depicts the maximum eigenvalue of the symmetric part of Ê⊤ M Â with M from (23) in dependence on the parameter θ from (24).We observe that the negative definiteness is lost for θ ≥ θ 0 ≈ 0.04.Although our relevant case θ = 1 is not dissipative, the stability is preserved in all ROMs.Moreover, the error of the MOR is not compromised.

Band-pass filter
We examine the electric circuit of a band-pass filter shown in Figure 6.A single input voltage is supplied and a single output voltage drops at a load conductance.Modified nodal analysis [10] yields a linear system of DAEs of dimension n = 23.This DAE system exhibits the index one.The physical parameters are 7 capacitances, 7 inductances and 9 conductances (q = 23).Figure 7 depicts the Bode plot of the system for a constant selection of the parameters.Furthermore, the system matrices are affine-linear functions of the parameters.
In the stochastic modeling, we introduce independent uniformly distributed random variables varying 20% around their mean values.The truncated series (15) includes all basis polynomials up to degree two, i.e., m = 300.The stochastic  Galerkin system (16) consists of linear DAEs, which are strictly proper for the defined outputs.Thus the H 2 -norm of the system exists.
In the regularization (26), we choose the parameters β = 10 −5 and α = β 2 .Table 4 shows the properties of the regularized system.The system is asymptotically stable.Both Ê and Êreg are not symmetric.Thus the spectral abscissa of the symmetric part of Âreg is irrelevant.Multiplication by Ê⊤ reg from the left yields a system with symmetric mass matrix, which is still not dissipative.
The one-sided Arnoldi method yields the projection matrices of the MOR for dimensions r = 1, . . ., 100.Therein, the expansion point is s = 10 6 .A loss of stability occurs in the reduction of both the DAE system and the regularized system, as shown in Table 5.We apply the stabilization techniques of Section 4.1 and Section 4.3 to the regularized system.The Lyapunov equations are solved using the identity matrix as input matrix.We solve the Lyapunov equations by a direct method of linear algebra.Thus a critical behavior for small regularization parameters, as mentioned in Section 5.2, is avoided.The matrices of a transformed Galerkin system are never computed explicitly.All ROMs become Finally, we compute approximations of the H 2 -norms for the difference between the DAE system and the reduced systems.Figure 8 illustrates the relative errors of the reductions.We recognize that the errors are nearly identical in the two stabilization techniques.The errors stagnate for reduced dimensions r > 80 in (b)-(d), because the total error is dominated by the error of the regularization in this part.Most important, the stabilization approaches do not compromise the error of the MOR.

Conclusions
We examined stability-preserving model order reduction of linear stochastic Galerkin systems using transformations to dissipative forms.Three approaches were analyzed.The transformation of the conventional Galerkin system represents an adequate method.The Galerkin projection of the transformed original systems features severe drawbacks in the case of large numbers of random parameters.The approach using a cheap transformation matrix based on a reference parameter is promising.Although the dissipativity property cannot be guaranteed in the transformed system, the numerical results of test examples demonstrate that preservation of stability is achieved.Moreover, the error of the model order reduction does not increase in this stabilization.5, in band-pass filter example.

Figure 1 :
Figure 1: Flowchart for transformations, stochastic Galerkin (SG) projections and model order reduction (MOR).Inputs and outputs are not shown.

Figure 3 :
Figure 3: Bode plot of mass-spring-damper system for a constant choice of the parameters.

Figure 4 :
Figure 4: Relative errors in H 2 -norm for MOR of stochastic Galerkin system and stabilization techniques in mass-spring-damper example.

Figure 6 :
Figure 6: Electric circuit of band-pass filter.

Figure 7 :
Figure 7: Bode plot of band-pass filter for a constant choice of the parameters.

Figure 8 :
Figure 8: Relative errors in H 2 -norm for MOR of different FOMs and stabilization techniques, see Table5, in band-pass filter example.

Table 1 :
Properties of stochastic Galerkin system in mass-spring-damper example.

Table 2 :
Number of asymptotically stable systems out of 100 in stabilization using frequency domain integrals for mass-spring-damper example.

Table 3 :
Number of non-zero entries in matrices of stochastic Galerkin systems for mass-spring-damper example.

Table 4 :
Properties of regularized stochastic Galerkin system in band-pass filter example.

Table 5 :
Number of asymptotically stable systems out of 100 in MOR of stochastic Galerkin system for band-pass filter example.