Second-Order Arnoldi Reduction Using Weighted Gaussian Kernel

Modeling and design of on-chip interconnect continue to be a fundamental roadblock for high-speed electronics. The continuous scaling of devices and on-chip interconnects generates self and mutual inductances, resulting in generating second-order dynamical systems. The model order reduction is an essential part of any modern computer-aided design tool for prefabrication verification in the design of on-chip components and interconnects. The existing second-order reduction methods use expensive matrix inversion to generate orthogonal projection matrices and often do not preserve the stability and passivity of the original system. In this work, a second-order Arnoldi reduction method is proposed, which selectively picks the interpolation points weighted with a Gaussian kernel in the given range of frequencies of interest to generate the projection matrix. The proposed method ensures stability and passivity of the reduced-order model over the desired frequency range. The simulation results show that the combination of multi-shift points weighted with Gaussian kernel and frequency selective projection dynamically generates optimal results with better accuracy and numerical stability compared to existing reduction techniques.


I. INTRODUCTION
In the past several decades, the complexity of Integrated Circuits (IC) and on-chip interconnects presented great challenges for the prefabrication validation of design [1]- [3]. While the technology advances for the manufacturing of ICs, the accurate modeling of interconnect becomes an important factor in measuring the performance and validation of design [1], [2], [4]. This forms the basis for the development of accurate representation and precise design tools as an essential part of the current VLSI (Very Large Scale Integration) design softwares [5]. The current CAD (Computer-Aided Design) design tools, using 3-D extraction, rely on the partial inductance along with related other components which usually results in a very large and dense extracted inductance matrix [5]- [7]. This often results in numerical inaccuracies The associate editor coordinating the review of this manuscript and approving it for publication was Wen-Sheng Zhao .
in the simulation process and subsequent difficulties in the approximation of large-scale dense systems. Researchers are now using the alternative methods by extracting susceptance, which is a newly introduced circuit element and also accounts for the coupling effect essential for the accurate representation of on-chip interconnects and structures [4], [8], [9]. It is important to note that the susceptance-based extraction results in a diagonally dominant state-space as its effects decrease with an increase in the distance from the source. This newly formulated extraction method including the susceptance effect gets transformed into state-space often known as RCS (Resistor Capacitor and Susceptance) formulation and is usually solved using the Second-Order MOR method, as it retains many of the inherent properties in its second-order formulation [10], [11]. In the classical design flow, the extraction is carried out using various field solvers, which later combine together to formulate a largescale state-space representation [12]- [15]. The MOR (Model Order Reduction) is part of the design tools and provides the much-needed solutions for the reduction of the complexity of state formulation in the simulation of ICs along with its associated parasitics [3], [16]- [18]. It is important to observe here that the state formulation as part of ICs design can be expressed either as a first or second-order system.
In second-order MOR for on-chip interconnects, the ENOR (Efficient Nodal Order Reduction) method was first proposed, which uses positive definiteness of the original system to give passive reduced-order model [19]. However, the recursive nature of the ENOR to compute moments leads to numerical instability. In order to solve this problem, SMOR (Second Model Order Reduction) improves the numerical stability and has been generally used for approximation of second-order dynamical systems [2], [20], [21]. It is important to observe here that SMOR actually approximates the moments of the original system and quite often does not give the required accuracy of the reduced system [2], [22]. In considering second-order systems, another notable contribution is SAPOR (Second-order Arnoldi method for Passive Order Reduction) [2]. The SAPOR works on RCS (Resistance Capacitance and Susceptance) extracted staterepresentation, however, requires expensive matrix inversion to create an orthonormal projection matrix for MOR. In addition, moment matching at interpolation points needs matrix inversion for each Arnoldi iterations [23], [24]. Therefore, the use of SAPOR may be numerically very expensive for order-reduction of large-scale second-order systems [2], [3], [25], [26]. Although these considerations of passivity, stability of reduced-order model, and numerical efficiency are of theoretical as well as of practical interest.
The study of various second-order MOR methods in some details to determine their limitations, is essential to have a clear understanding of second-order MOR of largescale dynamical systems. Another passivity preserving MOR method based on Arnoldi iteration is PRIMA [27]. It is important to note that PRIMA is not a numerically efficient technique and also works on first order state representation and like many other methods will not guarantee passivity [28]. Therefore, MOR of RCS circuit using first-order formulation may not ensure stability and passivity [2], [4], [29]. On contrary, second-order formulations of systems matrices accounting for susceptance not only preserve the theoretical properties but also ensure stability. In this work, the main contributions are 1) Formulation of a new Second-order MOR method SOAR-GK (Second-Order Arnoldi Reduction method using Gaussian Kernel), which adaptively selects expansion points weighted with Gaussian kernel. The proposed method is numerically efficient and also ensures passivity and preserves the stability of the system. 2) The simulation results of the SOAR-GK compared to existing state of the art first and second-order MOR methods gives better moment matching with less absolute error for various examples over a specific range of frequencies. In the next section, the problem of second-order model order reduction is formulated. The proposed Second-order Arnoldi Reduction using Gaussian Kernel (SOAR-GK) method is explained in Section 3. Simulations and results are discussed in Section 4 and Section 5 conclude the paper.

II. SECOND-ORDER MOR: PROBLEM FORMULATION
Model order reduction of on-chip interconnects falls into the main category of approximation of large-scale dynamical systems. In large-scale systems, the order of the problem representing on-chip interconnects and structures are quite often as large as a few million [30]. In addition, the operating frequency of the integrated circuit design is continuously increasing, which also affects the magnetic coupling between various interconnects. In order to correctly extract the magnetic coupling, the effect of susceptance modifies the matrices and generates second-order dynamical systems. The resulting modified nodal equations for an RCS circuit with susceptance can be written as follow [2], where matrices G, C ∈ R n×n , and S ∈ R m×m represent conductance, capacitance, and susceptance, respectively; I s (t) ∈ R m and V (t) ∈ R n represent the unknown vectors of susceptance currents and nodal voltages; I ∈ R m×m is the identity matrix; E s ∈ R n×m is incidence matrix for susceptance, and b, c ∈ R n are configuration vectors for current and voltage sources; u ∈ R and y ∈ R represent the input source current and output voltage, respectively. The modified nodal equations in the frequency domain can be written as follows where the Laplace transforms of V (t), I s (t), u(t), and y(t) are V (s), I s (s), U (s) and Y (s), respectively. This is a first-order system in the s domain. It is observed here that the susceptance currents I s (s) are usually the intermediate currents, therefore, the state representation can be modified by eliminating I s (s) from equation (2). Thus we get where T = E s SE T s . This expression obtained by eliminating I s (s) is equivalent to the state-space representation given by (1) and forms the basis of our second-order problem formulation [2].

A. SECOND-ORDER KRYLOV SUBSPACE REDUCTION
The general second-order Krylov subspace K(A, B, q) consists of a pair of square matrices (A, B) and a nonzero vector q. This section explains the general second-order Arnoldi reduction method, which generates the orthonormal basis for Krylov subspace reduction K(A, B, q). In the following, we give the definition of second-order Krylov subspace.

Definition (Second-Order Krylov Subspace):
Let A and B be square matrices of order n, and suppose q = 0 be a sequential order of n vectors given by where · · · We can write generalized formula of u i as which forms the basis of constructing second-order Krylov subspace based on A and B. The final representation of the nth second-order Krylov subspace using q vector is given as [21], [31] K(A, B, q) = span{u 0 , u 1 , u 2 , . . . , u n−1 }.
In order to obtain the moment vectors, a new variable ζ is introduced by substitution s = s 0 + ζ in (3). The substitution yields the following equation where D = 2s 0 C + G, K = s 2 0 C + s 0 G + T , b 0 = s 0 bu and b 1 = bu. To calculate the moments, we expand V (ζ ) using Taylor series expansion around s 0 . The resultant expression can be written as . By matching coefficients of ζ on both sides of equation (10), we get · · · · · · The above discussion leads us to calculate the individual moments and we can write the generalized recurrence formula of V j as follows The above recurrence relationship clearly expresses the calculation of moments of V . However, the moment calculation using (13) may lead to instability and a Krylov subspace-based method is suitable to yield stable results [2]. Note the introduction of the new variable ζ in (9), which shifts the second-order Krylov subspace from its first moment s 0 . It is observed here that a shifted Krylov subspace is better dealt by the modified second-order MOR, which can account for the moment matching at selected shifts. Therefore, the quadratic form of (9) is suitable for a modified second-order MOR method, which reduces the original RCS system from n to r. In this connection, the first step is to construct the orthonormal basis using moment vectors of nodal voltage up to order r. To illustrate this, let us introduce a new variable J (σ ) satisfying Using (14) into (9), we have a form appropriate for linearization as Combining (14) and (15), we get where By moving (I − σ A) to the RHS of (16) and performing a Maclaurin series expansion, we have Clearly, A i−1 q 0 p 0 is the ith moment of V J , with q 0 and p 0 be the first moments of V and J respectively. It is important to note that (17) is a linearized form of (9). It is obvious that if V J is the solution of (17), then V must be the solution of (9), which shows that the upper part of the ith moment of The key to the solution of this problem now is to find out the interpolation points for the projection matrix V , which should be selected in some sense to reduce the overall error. The selected interpolation points s k (n) are then used to obtain the orthogonal projection matrix V ∈ R n×r . To obtain the reduced-order model, V can be applied to (3), and the desired system can be obtained. Therefore, by performing an orthogonal projection using V , we obtain the following reduced-order system as are the reduced system matrices.

III. SECOND ORDER ARNOLDI REDUCTION USING WEIGHTED GAUSSIAN KERNEL (SOAR-GK)
The key objective of the second-order model order reduction is to create a ROM, which preserves the significant moments of the original system. The classical moment matching technique is central to many MOR methods used to approximate large systems [31]. In this section, we proposed a Second-Order Arnoldi Reduction using the weighted Gaussian Kernel (SOAR-GK) method, which is used to adaptively select multi-shift points weighted with normal distribution kernel in the selected range of frequencies to generate a projection matrix. The SOAR-GK method is stable and ensures passivity over the selected range of frequencies. It is important to note that the combination of multi-shift points weighted with Gaussian kernel and frequency selective projection dynamically generates optimal results with better accuracy and numerical stability compared to existing second model order reduction techniques.
First of all, in the proposed SOAR-GK method, the multi-shift points weighted with Gaussian kernel are selected in the given range as interpolation points. We begin by considering the original system transfer function P(s) : C → C of second-order dynamical system. This section explains the complete method of MOR using Second Order Arnoldi Reduction using the weighted Gaussian Kernel (SOAR-GK) method. The main purpose is to construct a projection matrix V by picking expansion points for Arnoldi iteration to correctly project the original system into reduced subspace. Let k max be the maximum number of expansion points for a given system, then the first step is to select shifts in some sense to reduce the absolute error between the original and reduced systems. Hence, given the original transfer function P(s), the selection of expansion points in SOAR-GK are weighted with Gaussian kernel for a given frequency range of interest to give optimal matching of the projected reduced-order model with the original system. The original system function P(s) can be written as where P(s) can be expressed using moments of second-order Krylov subspace. Expanding P(s) at zero-shift, i.e. the moments at s 0 = 0 are and for any shift at s k , we can write Usually s k ∈ C are the expansion points and in this case are chosen using Gaussian kernel. Similarly, expanding P(s) around infinity using Laurent series, then P(s) can be expressed as [32] P(s) = D + cBs −1 + cTBs −2 + · · · + cT n Bs −(n+1) + · · · .
It is important to remind here of the important concept of reachability, which is central to developing the Krylov subspace. It determines the extent to which the internal states can be manipulated with the excitation (input). For example, let a linear system be represented by := A B (neglecting C and D) and can be expressed as follows Expressing P(s) with moments around s k as The expression (24) gives moments of P(s) around selected expansion points s k [32]. The development of reduce-order model includes expanding P(s) at the expansion point s k as P(s) = P 0 (s k ) + P 1 (s k )(s − s k ) + P 2 (s k ) (s−s k ) 2 2! + · · · + P n (s k ) i.e. for a suitable r P k (s k ) ≈ P(s k ), k = 1, 2, 3, · · · , r.
The problem now comes down to explaining the criteria of selected expansion points. It is important to observe here that the inspiration comes from the fact that many natural phenomena follow a normal distribution. Normal distribution, commonly known as Gaussian distributions can differ in their means (µ) and standard deviations (σ ). The density of the Gaussian distribution function is plotted along x−axis and is governed by µ and σ for a given function. The probability density function of Gaussian distribution is given by In this work, we selectively pick the frequency points weighted with a Gaussian kernel, which serve as interpolation points for approximation using SOAR-GK. The mean and standard deviation of the Gaussian distribution is driven by the application-specific solution for a given problem, which also affects the range of the data. In our case, it is the interconnect of the modern processor with the difference between low-power and turbo-mode clock frequencies determining VOLUME 10, 2022 the range of the density function. The range z of the normal distribution act as input to the function to generate frequency points t p for a given µ and σ and is given by The interpolation points s k can now be defined as s k (n) = sort(t p [1 : k max ]) (slected shifts at t p ), for p = 1, 2, 3, · · · , r.
Note that t p with p = 1, 2 · · · r are the chosen expansion points using Gaussian kernel, where r is the order of the reduced system. This explains the procedure of selecting shifts, and it is now convenient to explain the complete proposed SOAR-GK method with the help of pseudocode in the next section to allow easy implementation in any programming language.

A. PSEUDOCODE FOR SOAR-GK
In this work, a second-order Arnoldi reduction using interpolation points selection weighted with Gaussian kernel is proposed. In terms of the basic procedure, the SOAR-GK is similar to the second-order Arnoldi reduction method. However, the projection using interpolation points selection weighted with Gaussian kernel leads to the reduced model with better numerical accuracy as well as retention of the stability and passivity inherent to the original system. In this section, the pseudocode of the basic algorithm used in the SOAR-GK method with all the necessary steps needed to implement in any programming language is explained. The input to the main function is the original system matrices i.e. C, G, T , and B, whereas r and k max represent the order of the reduced system and a maximum number of selected expansion points, respectively. The main function in the proposed method acts on the input and initially finds the interpolation points for the selected range, mean and standard deviation using Gaussian distribution. Note that the SOAR-GK method is a frequency selective estimation of the large scale-system. In application specific solutions, it is the range of frequency of interest (s 1 to s 2 ) i.e. the difference between low-power and turbo-mode clock frequencies. In our proposed model, the frequency selection is performed by normally distributing the frequency of interest and selecting the interpolation points, which follow a normal distribution in the reduced norm sense. The selected frequency becomes the interpolation point for constructing the projection matrix V [33], [34]. The normal Gaussian kernel can be defined as Considering the selected shift t p , we generate s k as the interpolation points to construct the orthonormal matrix. The matrix is constructed by following Gram-Schmidt orthogonalization process as q k = (s 2 k (n)C +s k (n)G+T ) −1 b k , using s k . Finally, we construct the projection matrix V using q k and Krylov matrix H . Following the procedure, 2,3, . . . , n. To construct the orthonormal basis matrix V , each column corresponds to the selected expansion point i.e. s 1 , s 2 , s 3 , . . . , s n . Using projection matrix V , this generates the reduced order modelG,B,C,D andT . At this point, let us take an example of order n = 5 to show the implementation of SOAR-GK. Note that the selected range of frequency for a given application is a subset of frequencies of the original large scale dynamical system. Example 3.1: Consider a second-order dynamical system of an inductor with n = 5. The state representation can be written as [33], [34], Using SOAR-GK method and the procedure explained in Algorithm 1 for every iteration with r = 1, 2, 3, the selected multi-shift interpolation points generates the following projection matrix V as

IV. SIMULATION AND RESULTS
This section discusses the performance as well as numerical efficiency of the proposed SOAR-GK (Second Order Arnoldi Reduction Gaussian Kernel) method. The problem of showing simulation and analysis is solved by selecting two examples of on-chip interconnect structures i.e. an interconnect inductor and an RLC (Resistance, Inductance, Capacitance) network with their second-order state representation. The models are generated by applying the field solvers for parasitic extraction [12], [13]. The nodal analysis equations are used to create state-space representation as explained in (19).
The generated examples are of order 11 and 53 for the interconnect inductor and RLC network, respectively. All simulations are carried out using MATLAB. In addition, the order of the reduced model is kept the same for the two examples i.e. 7 and 24 for the spiral and RLC network  respectively. This enables to keep the same time for carrying out the simulation of reduced-order models, thus, allowing for comparison of minimum, maximum, and mean error of the proposed SOAR-GK with all techniques. In the first example of a lower order interconnect inductor, the reduced-order modeling results of the proposed model are compared with the PRIMA [27], Single Shift Arnoldi [35] and Uniform Multi-shift [35]. The original system response of interconnect inductor is shown in Figure 1. The comparison results in Figure 2 with projection showing the selected frequency of interest from 1.8 GHz to 4 GHz, gives a close match of SOAR-GK with an overlapping peak and response compared to the original system. The absolute error in dB showing comparison with state of the art MOR methods is plotted in Figure 3. The error plot shows that our proposed SOAR-GK method has considerably reduced error relative to the PRIMA, single shift Arnoldi and uniform multi-shift for the selected range of frequency.
The second example is a state-space model of the order 53 RLC network. The original response of the RLC network is plotted in Figure 4. The example is taken keeping in mind state of the art Intel 7 (i7), 8th generation modern processor [36]. The processor works on dual frequencies, which is defined in Figure 5 with the selected range of frequency. Figure 5 plots the magnitude response of original order n = 53 interconnect RLC network, with   reduced order r = 24 response using PRIMA [27], singleshift Arnoldi [35], uniform multi-shift [35] and proposed SOAR-GK. The result of the comparison in term of absolute error plotted in Figure 6 shows that a reduced order 24 using the proposed SOAR-GK the method correctly captures the behavior of the original system. However, the PRIMA implementation misses most of the resonance peaks and does not accurately characterize the original system response. It is observed that the Arnoldi-based PRIMA matches the low-frequency response first and a reduced order of 24 is not sufficient to accurately cover the selected frequency range. The absolute error plot in Figure 6 also highlights the accuracy and numerical efficiency of our proposed SOAR-GK method, showing the minimum error compared to other techniques over the entire range of frequencies i.e. from 1.8 GHz to 4 GHz.
The error comparison of the two examples for various methods is summarized in Table 1 and Table 2. In Table 1, we compare the minimum, maximum, and mean error for the   Table 2 shows RLC network errors for various techniques, with SOAR-GK again showing lowest mean and maximum error of −161.55dB and −7.32dB. However, we note that the minimum error on PRIMA is −294.16dB compared to −272.05dB for the SOAR-GK, which is due to the Arnoldi-based PRIMA matching the low-frequency response first. Nevertheless, the overall response of SOAR-GK is far better as shown by the maximum and mean error in Table 2 for the RLC example.

V. CONCULSION
Second-order MOR plays an important role in the modeling and design of the integrated circuits to account for the parasitic effects created in the design of on-chip components and interconnects. In this work, we developed an accurate SOAR-GK method for the MOR of second-order systems, which selectively picks the expansion points weighted with Gaussian kernel in the given range of frequencies of interest. The proposed method ensures stability and passivity of the ROM over the desired frequency range, with simulation results showing better accuracy and numerical stability compared to the existing reduction methods.