Efficient Algorithms for Eigensystem Realization using Randomized SVD
Rachel Minster Arvind K. Saibaba Jishnudeep Kar Aranya Chakrabortty
math.NA cs.NA

Pressure and electric field dependence of quasicrystalline electronic states in 30 twisted bilayer graphene
    Shengjun Yuan
    December 30, 2023
===========================================================================================================

 Eigensystem Realization Algorithm (ERA) is a data-driven approach for subspace system identification and is widely used in many areas of engineering. However, the computational cost of the ERA is dominated by a step that involves the singular value decomposition (SVD) of a large, dense matrix with block Hankel structure. This paper develops computationally efficient algorithms for reducing the computational cost of the SVD step by using randomized subspace iteration and exploiting the block Hankel structure of the matrix. We provide a detailed analysis of the error in the identified system matrices and the computational cost of the proposed algorithms. We demonstrate the accuracy and computational benefits of our algorithms on two test problems: the first involves a partial differential equation that models the cooling of steel rails, and the second is an application from power systems engineering.  


  Eigensystem Realization Algorithm, System Identification, Singular Value Decomposition, Randomized Algorithms.



  93B30, 93B15, 65F15, 37M10, 15A18.




§ INTRODUCTION

Linear Time Invariant (LTI) systems form one of the most important classes of dynamic systems existing in the physical world. LTI systems closely approximate the linear behavior of any nonlinear physical process around a desired operating point, and provide insights on how the system will behave in response to any small-signal change in its equilibrium state. The state-space representation of a LTI discrete-time system can be written as 

    x_k+1   =Ax_k+Bu_k,     x_0 = x(0) 
    y_k    = Cx_k+Du_k,

where n is the order of the system, k=0, 1, 2,… is the sampling index, x_k ∈R^n is the state vector, u_k ∈R^m is the input vector, and y_k ∈R^ℓ is the output vector at sampling instant k, with m and ℓ being the number of inputs and outputs, respectively. The matrix A∈R^n × n is referred to as the state matrix, B∈R^n × m is the input matrix, C∈R^ℓ× n is the output matrix, and D∈R^ℓ× m is the input-output feedthrough matrix. In the majority of practical applications, however, these four matrices may not be exactly known to the system modeler because of different kinds of model uncertainties and operational uncertainties. In that case, one must estimate these matrices from sampled measurements of the input sequence {u_k} and the output sequence {y_k}. This process is referred to as system identification. Since the state may be represented in different coordinate frames without changing input-output characteristics, the more formal definition of system identification is given as the process of estimating (A,B,C,D) up to a similarity transformation (TAT^-1,TB,CT^-1,D), where T∈R^n× n is an invertible matrix. As mentioned, such a transformation does not change the input-output behavior or the transfer function of the system. We assume that the system to be identified is observable and reachable <cit.>.  There exists a vast literature on system identification with a variety of numerical algorithms, developed for various applications such as electric power systems, process control, mechanical systems, aerospace applications, civil and architectural applications, and chemical and biological processes to name a few. For a survey of these methods, please see  <cit.>. In this paper, we focus on identifying the system matrices when the inputs {u_k} are in the form of impulse functions. 



One of the most common identification methods used for identifying LTI models in practice is known as the Eigensystem Realization Algorithm (ERA), initially proposed in <cit.>, but can be found in many other sources such as <cit.>. ERA is a subset of a broader class of identification methods called subspace system identification, and is a purely data-driven approach consisting of two main steps. In the first step, this algorithm computes the singular value decomposition (SVD) of a block Hankel matrix, constructed using the impulse response data. In the second step, the truncated SVD is used to solve a least-squares optimization problem to identify the system matrices. The first step is computationally expensive since it has cubic complexity in the dimensions of the matrix. To reduce the computational cost, the tangential interpolation-based ERA (TERA) <cit.> was proposed. This method reduces the number of inputs and outputs, thus reducing the dimensions of the matrix whose SVD needs to be computed. Another algorithm, known as CUR-ERA <cit.>, uses the CUR decomposition to efficiently compute the low-rank decomposition of the block Hankel matrix. While these algorithms have successfully lowered the costs, important challenges remain. 





  
Contributions and Content We propose new randomized algorithms for tackling the computational costs of ERA. The first algorithm (Section <ref>) accelerates the standard randomized SVD by exploiting the block Hankel structure to efficiently perform matrix-vector products. The second algorithm (Section <ref>) employs similar ideas as the first algorithm in combination with the TERA. The resulting algorithms are efficient both in terms of storage costs and computational costs. In Section <ref>, we derive new error bounds that provide insight into the accuracy and stability of the system matrices identified using the approximation algorithms. The error analysis is not tied to any particular algorithm and in this sense is fairly general. Finally, we demonstrate the benefits of the proposed algorithms on two different applications: first, a heat transfer problem for cooling of steel rails, and second, a dynamic model of an electric power system with multiple generators and loads.





  
Comparison to related work While TERA works with smaller matrices, it still has cubic complexity in the number of time measurements. This is computationally demanding when the dynamics of the system are slow, and many measurements in time are needed to fully resolve the dynamics. The CUR-ERA algorithm has similar storage and computational complexity as the algorithms we propose but numerical experiments (see Section <ref>) suggest that our algorithm is more accurate.  Furthermore, the randomized algorithms developed here can exploit parallelism in multiple ways; the matrix-vector products can be parallelized across multiple random vectors as well as over the input/output pairs. This makes it attractive for implementations on high performance computing platforms. If the target rank for the low-rank decomposition is not known in advance, one can use randomized range finding algorithms <cit.> to estimate the target rank.  Two other features make our contributions attractive: first, the block Hankel structure that we exploit here can be adapted to any other matrix-free low-rank approximation algorithm, for example, based on the Lanczos bidiagonalization <cit.>, and second, the error analysis developed here is informative for any approximation algorithm. There are other randomized algorithms for system identification <cit.> but they tackle slightly different settings.




§ BACKGROUND

In this section, we first review the Eigensystem Realization Algorithm (Section <ref>), the associated computational costs, and motivate the need for efficient algorithms.  We also review the key ingredients needed to construct our algorithms: algorithms for storing and computing with Hankel matrices (Section <ref>), and randomized SVD (Section <ref>) for computing the low-rank factorizations.


 §.§ Eigensystem Realization Algorithm

ERA was first proposed in <cit.>, but we follow the formulation in <cit.>. Other references for the ERA include <cit.>. Now, assuming the initial condition x_0 = 0∈R^n, and given a sequence of inputs u_i for i= 1,…, the outputs observed are 

    y_k = ∑_j=0^k-1h_k-j-1u_j   k = 0,1,…,
 where the matrices {h_k} are called the Markov parameters and are given by

    h_k = {[             D         k = 0;       CA^k-1B k = 1,…,2s-1. ].


 The system is excited m times using impulse input excitations. That is, we take 

    u_k^(j) = {[   e_j  k =0;     0 k > 0 ].      j =1,…,m.
 
Here e_j is the jth column of the m× m identity matrix. The outputs from each impulse excitation can be used to construct the Markov parameters {h_k }. The ERA uses the Markov parameters to recover the system matrices. If data using impulse inputs is not available, then one can estimate the Markov parameters from the general input data <cit.>. 


These Markov parameters are first arranged to form the block Hankel matrix H_s∈R^(sℓ) × (sm) defined as 

    H_s   = h_1    h_2    …   h_s
    h_2    h_3    …   h_s+1
    ⋮   ⋮   ⋰   ⋮
    h_s    h_s+1   …   h_2s-1
       = CB   CAB   …   CA^s-1B
    CAB   CA^2B   …   CA^sB
    ⋮   ⋮   ⋰   ⋮
    CA^s-1B   CA^sB   …   CA^2s-2B.

Here, s is a chosen parameter that determines the number of block rows and columns of H_s. It is known that this matrix can be factorized into the observability matrix O_s and the controllability matrix C_s, as H_s = O_sC_s, where 

    O_s = C
    CA
    ⋮
    CA^s-1,     C_s =B   AB   …   A^s-1B.

If the system is assumed to be reachable and observable, then (H_s) = (O_s), where (·) denotes the range (or column space) of a matrix. We can obtain a basis for O_s using  the reduced SVD of H_s = U_nΣ_nV_n^⊤. Then, partition the left singular vectors as 

    U_n = Υ_f 
    * = *
    Υ_l,

where Υ_f ,Υ_l∈R^(s-1)ℓ× n. We can obtain A using the formula A = Υ_f^†Υ_l, where ^† denotes the Moore-Penrose pseudoinverse. Furthermore, we can obtain the output matrix C and the input matrix B using the formulas

    C = I_ℓ   0Υ_f,     B = Σ_nV_n^⊤I_m
    0.

The matrix D can be identified as the Markov matrix h_0 and, therefore, we will not discuss its estimation in future sections.




  
Reduced order modelIn some applications, the goal is not only to identify the system but also obtain a reduced order model of the system. That is, we seek the reduced system matrices A_r ∈R^r× r, B_r ∈R^r × m, C_r∈R^ℓ× r and D_r ∈R^ℓ× m which approximate the dynamics of the original system (<ref>).
To accomplish this, as before, we first compute a rank r approximation to the matrix H_s as 

    H_s ≈U_rΣ_rV_r^⊤,

where r ≤ n. We partition the left singular vectors U_r as 

    U_r = Υ_f^(r)
    * = *
    Υ_l^(r),

such that Υ_f^(r)∈R^(s-1)ℓ× r and Υ_l^(r)∈R^(s-1)ℓ× r. Then, we compute the reduced order model A_r = [Υ_f^(r)]^†Υ_l^(r) such it that minimizes the least squares problem 

    min_A∈R^r× rΥ_f^(r)A - Υ_l^(r).
 The reduced order output matrix C_r and the input matrix B_r are computed using the formulas

    C_r = I_ℓ   0U_r,     B_r = Σ_rV_r^⊤I_m
    0.



Other papers, such as <cit.>, use a slightly different representation for the system matrix A_r than the one used in <cit.>. In this alternate representation, A_r' is obtained using A_r' = Σ_r^-1/2Υ_f^†Υ_l Σ_r^1/2. Then, B_r' and C_r' are obtained using 
    C'_r= I_ℓ   0Υ_f Σ_r^1/2,   B_r' = Σ_r^1/2V_r^⊤I_m 
    0.
 

Note that the two formulations are equivalent up to a similarity transformation, and this results in the same Markov parameters and input-output behavior of the system. If r = n, then the ERA determines the system matrices (A,B,C,D) and there is no model reduction step. On the other hand, if the target rank r < n, then the reduced order model is guaranteed stability under the conditions of <cit.>.



  
Computational Cost We briefly review the range of possible parameters. 
    
  * State space: In the power system applications, the dimension of the state space n is around 10^2, whereas in applications with partial differential equations this dimension can be large, i.e., 10^5-10^7 (see for example, <cit.>).
    
  * Input dimension: The number of inputs m is in the range 1-10^2. 
    
  * Output dimension: The number of outputs ℓ is also in the range 1-10^2. In certain applications, the number of outputs is also the dimension of the state space. 
    
  * Sample size: When the dynamics of the system are slow, the number of samples 2s-1 can be large, i.e. 10^2-10^5.

The dominant cost of ERA is the cost of storing and factorizing the matrix H_s, which is of size sℓ× sm. The cost of storing the matrix H_s is O(s^2ℓ m), whereas the cost of computing the SVD is O(s^3ℓ m min{ℓ,m}). In <cit.>, the authors discuss applications in which the size of H_s is O(10^4). However, for the higher end of the parameters listed above, the matrix H_s is of size O(10^7). Forming and factorizing this matrix is completely infeasible. The algorithms we propose are both efficient in storage and computational cost and make ERA applicable to larger problem sizes. 



 §.§ Hankel matrices

Structured matrices such as circulant and Hankel matrices are computationally efficient to work with since matrix-vector products (matvecs) can be accelerated using the Fast Fourier Transform (FFT). We first review circulant matrices. Circulant matrices are completely determined by their first column 
    x =  x_1    …    x_N^⊤,
 and are diagonalized by the Fourier matrix. Let the circulant matrix X be defined as

    X = x_1     x_N    …    x_2   
     x_2     x_1    …    x_3 
    ⋮   ⋮   ⋱   ⋮
     x_N     x_N-1   …    x_1 ,

and let F_N be the N× N Fourier matrix with entries (F_N)_jk = e^2π ijk/n. Then the eigenvalue decomposition of X is X = F_N^*diag(𝐅_Nx) 𝐅_N, where x = X(:,1) is the first column of X. This means that the circulant matrix has eigenvalues 𝐅_Nx. This result implies that matvecs Xv can be performed efficiently using FFTs as y = IFFT( FFT(v) ⊙FFT(x)), where ⊙ is the elementwise product, FFT(·) and IFFT(·) denote the fast Fourier and inverse fast Fourier transforms respectively. The computational cost involves 2 FFTs and one inverse FFT, and can be implemented efficiently in O(Nlog_2 N) floating point operations (flops)

Hankel matrices have constant entries along every anti-diagonal; that is, given the parameters h_1,h_2,…, h_2s-1,
the Hankel matrix is 

    H_s = h_1     h_2    …    h_s 
     h_2     h_3    …    h_s+1
    ⋮   ⋮   ⋰   ⋮
     h_s     h_s+1   …    h_2s-1.

This corresponds to the single input/single output (SISO) case. Note that permuting a Hankel matrix with the reverse identity permutation matrix J_s results in a Toeplitz matrix (constant entries along every diagonal). That is, 

    J_s =  0     0    …    0      1 
    0     0    …    1     0 
    ⋮   ⋮   ⋰   ⋮   ⋮
     0      1    …    0     0 
      1     0    …    0     0  ,     H_sJ_s = h_s     h_s-1   …    h_2      h_1 
    h_s+1    h_s     …    h_3     h_2 
    ⋮   ⋮   ⋱   ⋮   ⋮
     h_2s-2   ⋮   …    h_s     h_s-1
      h_2s-1    h_2s-2   …    h_s+1    h_s .

We can use this to our advantage while computing matrix-vector products with Hankel matrices. To compute the matrix-vector product y = H_sv, we first write y = (H_sJ_s)(J_sv) so that 

    X_2sJ_sv
    0_s = H_sJ_s    B_s 
    B_s    H_sJ_sJ_sv
    0_s = y
    *,

where * denotes the part of a computation which we can ignore and 

    B_s =  0     h_2s-1   …    h_s+2    h_s+1
     h_1     0     h_2s-1   …     h_s+2
    ⋮   ⋮   ⋱   ⋮   ⋮
     h_s-2   ⋮   …    0     h_2s-1
      h_s-1    h_s-2   …    h_1    0  .



This means that we can also efficiently compute matvecs with Hankel matrices by embedding it within a 2s × 2s circulant matrix X_2s defined by the vector 

    x_2s =  h_s     h_s+1   …    h_2s-1    0      h_1    h_2   …    h_s-1^⊤.

Thus, only the first row and the last column need to be stored to compute matrix-vector products with the Hankel matrix H_s. A summary of the algorithm to compute matvecs with H_s is given in Algorithm <ref>.





 §.§ Randomized SVD

We can efficiently compute a rank r approximation of an M × N matrix X using a randomized version of the SVD <cit.> (henceforth called RandSVD). The idea is to find a matrix Q that approximates the range of X.  This is done by first drawing a Gaussian random matrix Ω∈R^N × (r+ρ), where r is the desired target rank, and ρ≥ 0 is an oversampling parameter (typically ρ≤ 20).  Then, the matrix Y = XΩ consists of random linear combinations of the columns of X.  This means that we can get a matrix Q such that (X) ≈(Q) by taking a thin QR factorization Y = QR.  If X has singular values that decay rapidly, or (X) is exactly r, then (Q) is a good approximation for (X).  We can then approximate X by the low-rank representation X≈QQ^⊤X; this can then be converted into the appropriate SVD format. The procedure is summarized in Algorithm <ref>.



RandSVD is computationally beneficial compared to the full SVD. In this paper, we use a variation of the RandSVD that uses q≥ 0 steps of the subspace iteration <cit.>. Note that for numerical stability,  we perform orthogonalization during and in between the subspace iterations. Assuming M≥ N, the cost of the full SVD is O(MN^2). On the other hand, if the cost of a matrix-vector product with X (or its transpose) is T_X, then the cost of the randomized SVD can be expressed as 

    Cost = (2q+1)(r+ρ)T_X + O(r^2(M+N)).

We will use RandSVD in different ways in the system identification algorithms that we derive. We have found RandSVD with q=0-2 subspace iterations to be computationally efficient and sufficiently accurate for our purpose; however, there are several new randomized algorithms developed that have been reviewed in the recent paper <cit.>.



§ RANDOMIZED ALGORITHMS FOR EIGENSYSTEM REALIZATION

In this section, we derive two randomized algorithms for efficient computation of the system matrices. The first algorithm accelerates the standard randomized SVD using the block Hankel structure of H_s (Section <ref>); the second algorithm is a randomized variant of the Tangential Interpolation-based ERA (TERA) and is applicable when the number of inputs and outputs are large. 


 §.§ Randomized Eigensystem Realization Algorithm

Our first approach accelerates the computation of the system matrices by combining two ingredients: first, we replace a reduced SVD of H_s by a RandSVD to obtain an approximate basis for O_s; second, we additionally exploit the block Hankel structure of H_s to accelerate the matvecs involving H_s and H_s^⊤ in the RandSVD. As we will show, each of these steps decreases the computational complexity yielding an efficient algorithm overall. 



  §.§.§ Block Hankel Matrices
 
We first explain how we exploit the block Hankel structure of the matrix H_s defined in (<ref>). The multiplication process for Hankel matrices can be extended to block Hankel matrices. Suppose we have to compute y = H_sx for a given vector x. Let us define the index sets 

    I_i =    {i,i+ℓ, …, i+ (s-1)ℓ}        i = 1,…,ℓ
    J_j =    {j,j+m, …, j+ (s-1)m}        j = 1,…,m.

If we denote H_s(I_i,J_j) as the s× s submatrix obtained by extracting the rows and the columns defined by the appropriate index sets, then it is clear that  H_s(I_i,J_j) is a Hankel matrix, and that 

    y(I_i) = ∑_j=1^m H_s(I_i,J_j)x(J_j),      i = 1,…,ℓ,

where y(I_i) and x(J_j) are the s× 1 vectors obtained from y and x respectively.



Algorithm <ref> gives the details of the procedure described here in MATLAB-like notation. It is important to note the following points. First, we need not actually form either the full block Hankel matrix or the intermediate Hankel matrices for each block element. Since the Hankel matrices are defined by the first row and the last column, we extract these quantities from the blocks {h_k }_k=1^2s-1 as and when required. Second, since each Hankel matvec costs O(slog_2 s), the overall cost of one matvec is O(ℓ m s log_2 s), compared to O(ℓ m s^2) using the naïve approach. Finally, we can easily adapt this algorithm to compute H_s^⊤x as well; the main difference involves taking as inputs the transpose of the Markov parameters {h_k^⊤}_j=1^2s-1 instead. 






  §.§.§ Randomized ERA

We now incorporate the block Hankel multiplication algorithm, Algorithm <ref>, with RandSVD in order to accelerate the system identification process. This is simple to do; every time we need to multiply H_s or H^⊤_s, we use block Hankel multiplication instead. This is beneficial in two ways: to reduce the computational cost by two orders of magnitude (one from RandSVD and one from using block Hankel structure), and to reduce storage. We need not store H_s explicitly or even form the full matrix to begin with. All we need are the blocks that make up H_s. This is a major benefit over the standard algorithms. 

 First, we use the RandSVD algorithm to compute a low-rank approximation of H_s with target rank r ≤ n to obtain U_rΣ_rV_r^⊤. As mentioned earlier, the matvecs involving  H_s or H_s^⊤ are handled using Algorithm <ref>. Then, in the system identification phase, we partition the left singular vectors as

    U_r = Υ_f 
     = * 
    Υ_l,

where Υ_f and Υ_l are both (s-1)ℓ× r. The system matrices can then be recovered as A_r = Υ_f^†Υ_l, 

    C_r = I_ℓ   0U_r,     B_r = Σ_rV_r^⊤I_m 
    0.

As before, the matrix D_r is simply the first Markov parameter h_0. We will refer to this randomized ERA that uses block Hankel multiplication as RandSVD-H. Now, we review the computational cost of this algorithm and in Section <ref>, we derive error bounds for the recovered system matrices.


  
Computational cost
We now examine the computational cost of the three system identification algorithms described so far, namely the ERA using a full SVD, ERA using RandSVD, and RandSVD-H. The dominant cost of each algorithm is the SVD step, so we focus our attention there. The complexity of the SVD step for these algorithms is shown in Table <ref>. Recall that the size of the matrix we are computing with is sℓ× sm, the target rank for each RandSVD is r, and the size of the system is n. The cost of the full SVD is then O(s^3mℓmin{ℓ,m}). To analyze the RandSVD based algorithms, we follow the analysis of cost in (<ref>). For a standard RandSVD, the cost of a matvec is O(s^2ℓ m). If the matrix H_s is stored explicitly as we do in our implementation, then the storage cost is O(ms^2ℓ); however, RandSVD can be implemented without storing H_s explicitly, in which case the cost of storage is O(msℓ). When we use block Hankel structure to accelerate the RandSVD, the cost of a matvec is reduced to O(ℓ m s log_2 s) and the storage cost is also O(msℓ). This shows that, as anticipated, replacing the SVD with a RandSVD reduces the cost, and then exploiting the block Hankel structure reduces the cost even further. We also include, for comparison purposes, the computational and storage costs of CUR-ERA <cit.>.

 



  
CUR-ERA We now include a brief description of CUR-ERA. The CUR-ERA algorithm relies on the principle of finding a maximum volume sub-matrix for computing a low-rank approximation; that is, a sub-matrix of specified dimensions with the largest determinant in absolute magnitude. Finding the maximum volume sub-matrix is a combinatorial optimization problem and, hence, one has to settle for heuristics to compute a nearly maximum volume sub-matrix in a reasonable computational time. The two approximations CUR-ERA uses are finding the cross approximation of a matrix, and finding a dominant volume submatrix. Note that in the computational cost in Table <ref>, κ and c are the number of iterations used in the cross approximation and the dominant volume submatrix, respectively. These values can be large in practice, meaning that our algorithms have a comparable computational cost. Also, the randomized SVD algorithm does not involve a combinatorial optimization problem, is known to be computationally efficient and accurate for a range of problems, and has well-developed error analysis. This makes the RandSVD-H beneficial in practical applications.  In addition, as we will show in numerical experiments in Section <ref>, our algorithms are more accurate.



 §.§ Randomized TERA

This approach is inspired by the tangential interpolation approach for model reduction. The goal of TERA is to reduce the dimension of the Markov parameters by projecting the parameters into a lower dimensional space. This reduces the size of the block Hankel matrix but preserves the block Hankel structure, therefore, reducing the computational cost of the SVD step <cit.>. We briefly review the TERA approach and describe our acceleration using RandSVD. 




  
TERA In this approach, we seek two orthogonal projection matrices

    P_1 = W_1W_1^⊤       (W_1) = ℓ' 
    P_2 = W_2W_2^⊤       (W_1) = m',

and the matrices W_1 and W_2 have orthonormal columns. To compute P_1, we first arrange the Markov parameters in the matrix 

    H_w = h_1    …   h_2s-1∈R^ℓ× m(2s-1).

We then solve the optimization problem 
    min_(P) = ℓ'PH_w - H_w _F^2 .
 The optimal solution can be constructed using the SVD of H_w = U_wΣ_wV_w^⊤. We take W_1 = U_w(:,1:ℓ'); that is, we take W_1 to be the first ℓ' left singular vectors of H_w. Similarly, to construct the matrix W_2, we first arrange the Markov parameters into the matrix 

    H_e = h_1 
    ⋮
    h_2s-1∈R^ℓ (2s-1) × m.

We then solve the optimization problem 
    min_(P) = m'H_e - H_eP_F^2 .
 The optimal solution can be constructed using the SVD of H_e = U_eΣ_eV_e^⊤. We take W_2 = V_e(:,1:m'); that is, we take W_2 to be the first m' right singular vectors of H_e. Having obtained the matrices W_1 and W_2, we construct the projected Markov parameters as 

    h_i = W_1^⊤h_i W_2  ∈R^ℓ' × m',      i = 1,…,2s-1.

The block Hankel matrix H_s is “projected” using the matrices W_1 and W_2 arranged in diagonal blocks to obtain the reduced-size block Hankel matrix 

    H_s =  W_1^⊤
       W_1^⊤
          ⋱
             W_1^⊤h_1    h_2    …   h_s
    h_2    h_3    …   h_s+1
    ⋮   ⋮   ⋰   ⋮
    h_s    h_s+1   …   h_2s-1W_2 
       W_2
          ⋱
             W_2   .

It is important to observe that due to this projection, the block Hankel structure is preserved, and we can express H_s in terms of the “projected” Markov parameters as 
 
    H_s  = h_1    h_2    …   h_s   
    h_2    h_3    …   h_s+1
    ⋮   ⋮   ⋰   ⋮
    h_s    h_s+1   …   h_2s-1∈R^(ℓ's)× (m't).

If ℓ' ≪ℓ and m' ≪ m, then the size of the matrix H_s is much less than H_s. The dimensions ℓ' and m' are determined by the singular value decay of the matrices H_w and H_e. Retaining a larger number of singular vectors (that is, large ℓ' and m') results in a more accurate approximation to H_s but results in a larger matrix H_s and in a higher computational cost.




  
Recovering system matrices The next steps mimic the standard ERA approach to reconstruct the system matrices but we must carefully account for the dimensions of the projected system. We compute the SVD of H_s = U_rΣ_rV_r^⊤, and partition the left singular vectors as 

    U_r=Υ_f 
    * = *
    Υ_l.

The system matrices can be recovered as 

    A_r = Υ_f^†Υ_l ∈R^r × r,     C_r = W_1    0Υ_f,     B_r = Σ_rV_r^⊤W_2^⊤
    0 .



  
Computational cost To motivate the need for accelerating TERA using RandSVD, we first review the computational cost which has three components. Computing the projection matrices W_1 and W_2 involves a computational cost of O(s(ℓ m^2 + mℓ^2) ). The SVD of H_s now costs 

    O( s^3ℓ'm' min{ℓ',m'}),

and, similar to earlier algorithms, the cost recovery of the system matrices is negligible compared to the cost of the SVD. Although the computational cost is significantly reduced (assuming ℓ' ≪ℓ and m' ≪ m), the cost of the SVD still dominates the computational cost and has cubic scaling with s. We now propose a new algorithm to lower this cost. 





  
RandTERA It is worth pointing out that the “projected” Hankel matrix H_s still retains its block Hankel structure. To reduce the computational cost, we combine the following ingredients that were used previously: we use the RandSVD to compute the rank-r approximation to H_s, the matvecs involving H_s can be accelerated using the block Hankel structure (as described in Section <ref>). The main difference is that we use the “projected” Markov parameters {h_i}_i=1^2s-1 rather than the Markov parameters {h_i}_i=1^2s-1. As a result, the computational cost of the SVD step is reduced to 

    O( rℓ'm' slog_2s  + r^2s(m' + ℓ') ).

We call this algorithm RandTERA, and a complete description of this algorithm is given in Algorithm <ref>. A breakdown of the computational cost of the RandTERA is given below in Table <ref>.



Regarding the storage cost, since only the projected Markov parameters need to be stored, the storage cost is O(sm'ℓ'), which is potentially lower than RandSVD-H which requires O(smℓ).



§ ERROR ANALYSIS


In this section, we analyze the accuracy of RandERA in Section <ref>, and derive bounds on the accuracy of the recovered system matrix A_r. Our derivation is not tied to the randomized algorithms used to compute the approximations; this makes the analysis applicable to more general settings. 



 §.§ Background and assumptions
 We first review the necessary background information and clearly state our assumptions needed for the analysis. 


  
Canonical Angles Given H_s = UΣ V^⊤, partition U_r = U(:,1:r), the left singular vectors corresponding to the top r singular values of H_s, as

    U_r = Υ_f
    * = *
    Υ_l∈R^sℓ× r.

Note that we have dropped the superscripts (r), compared to Section <ref>, to make the notation manageable. We can recover the matrix A_r as  A_r = Υ_f^†Υ_l. Similarly, let _s = UΣV^⊤ be an approximation to H_s with left singular vectors U_r partitioned as

    _r = Υ_f
    * = *
    Υ_l∈R^sℓ× r,

and we can compute the approximate system matrix A_r = Υ_f^†Υ_l ∈R^r× r. Let the canonical angles between the subspaces (U_r) and (U_r) be denoted by

    0 ≤θ_1 ≤…≤θ_max≤π/2.

If we collect the angles into a diagonal matrix Θ, then  U_r^⊤_r has the SVD 

    U^⊤_r_r = P(cosΘ)Q^⊤.

Let U_r denote the orthogonal projector onto (U_r); similarly, let _r denote the orthogonal projector onto (_r). Then, the distance between the two subspace (U_r) and (_r) is given by 

    U_r - _r_2 = (I-U_r)_r_2 = sinΘ_2 =  sinθ_max.

See <cit.> for more details on canonical angles between subspaces.



  
Pseudoinverses We recall some facts about the perturbation of pseudoinverses <cit.>. If M,E∈R^n× r such that (M) = (M+E) = r then

    (M + E)^† - M^†_2 ≤√(2)M^†_2 (M+E)^†_2 E_2

Furthermore, if E_2 M^†_2 < 1 then

    (M+E)^†_2 ≤M^†_2/1-E_2M^†_2.

If A∈R^m× r and B∈R^r× n such that (A) = (B) = r, then by <cit.>

    (AB)^† = B^†A^†.





  
Accuracy of eigenvalues A norm based approach, i.e., A-_r, is not meaningful since  A_r (and _r) can only be determined up to a similarity transformation. To compare the approximate system matrix _r with A_r, we compare the eigenvalues of these two matrices since the eigenvalues remain unchanged by a similarity transformation.  We measure the accuracy of the eigenvalues of  using the spectral variation, which we now define. Let A,B∈R^n× n; the spectral variation between A and B is defined as <cit.>

    (ψ(B),ψ(A)) = max_1 ≤ j ≤ nmin_ 1≤ i ≤ n |λ_i(A) -λ_j(B)|,

where ψ(·) denotes the spectrum of the matrix. 

We briefly list the various assumptions that are required for our main result.
  We assume:


  A1. System matrix:  Assume that A_r = Υ_f^†Υ_l ∈R^r× r is diagonalizable and let W∈R^r× r be the matrix of eigenvectors.

  A2. Singular vectors: Assume that the subblocks Υ_f of the singular vectors U_r are full rank; that is  
(Υ_f) =  r. 

  A3.
Markov parameters: Assume the Markov parameters converge to h_i →0 for i > s. 

  A4. Canonical angles: Assume that the canonical angles are sufficiently small and satisfy

    η≡ 2sinθ_maxΥ_f^†_2 < 1.
 






Assumption A1 is rather strong but can be weakened, if different perturbation results are used; see <cit.>. Assumption A2 ensures that the least squares problem min_AΥ_fA - Υ_l_F  
has a unique solution. Assumption A3 is also made in <cit.>, and is necessary to ensure the stability of the system. Finally, Assumption A4 ensures that the approximate singular vectors are sufficiently accurate. 



 §.§ Main result
 We are ready to state our main theorem. 
 With the notation in Section <ref> and <ref> 

    (ψ(_r), ψ(A_r)) ≤κ_2(W) η( 1 +√(2)Υ_f^†_2/1-η),

where κ_2(W) = W_2W^-1_2 is the condition number of the eigenvectors W of A_r.

The theorem identifies several factors that can cause large errors during the identification step. First,  the eigenvectors of W may be ill-conditioned, i.e., κ_2(W) can be large. The best case scenario is when A is normal, so that the condition number is 1.  Second, the norm of the pseudoinverse Υ_f^†_2 ≥ 1 can be large. However, Assumption A3 ensures that lim_s→∞Υ_f_2 =1; see <cit.> for a detailed argument.  Finally, if the singular vectors of the low-rank approximation are not computed accurately, then the canonical angles can be large, i.e., sinθ_max can be close to 1. This can be mitigated, for example, by taking several subspace iterations. The first three assumptions are intrinsic to the system (A,B,C,D) and only Assumption A4 depends on the choice of the numerical method. The theorem is applicable to the randomized algorithms developed here, namely, RandERA and RandTERA. However, it is important in a larger context, since it is applicable to any approximation algorithm for ERA, including TERA.


There are several steps involved in this proof. 

1. Introducing a similarity transformation.  Let Z∈R^r× r be an orthogonal matrix; we will leave the specific choice of this matrix to step 2. Since A_r is diagonalizable, using a perturbation theorem for diagonalizable matrices <cit.>, we have

    (ψ(_r), ψ (A_r)) = (ψ(Z_rZ^⊤), ψ(A)) ≤κ_2(W) A_r - Z_rZ^⊤_2.

 The equality in the above equation is because the eigenvalues of _r and Z_rZ^⊤ are the same. We first discuss the choice of Z before analyzing the error in the term A_r - Z_rZ^⊤_2.

2. Choosing Z. Recall that the SVD of U_r^⊤_r is P(cosΘ)Q^⊤. Using the unitary invariance of the 2-norm, U_r-_rZ^⊤_2 = U_rZ - _r _2. Now, we choose Z = PQ^⊤ and verify that Z is orthogonal. By the triangle inequality and the unitary invariance of the 2-norm,

    U_rZ - _r _2 ≤   U_rZ - U_rU_r^⊤_r _2 + (I-U_rU_r^⊤)_r _2 
    ≤   U_rZ - U_rP(cosΘ)Q^⊤_2 + (I-U_rU_r^⊤)_r_r_2
    ≤   PQ^⊤ - P(cosΘ)Q^⊤_2 + (I-U_r)_r_2 _r_2
    ≤   I - cosΘ_2 + sinΘ_2.

Since the canonical angles satisfy 0 ≤θ_i ≤π/2, we have cos^2θ_i ≤cosθ_i and, therefore, for i=1,…,r

    (1-cosθ_i)^2 = 1 -2cosθ_i + cos^2θ_i ≤ 1-cos^2θ_i = sin^2θ_i.

Therefore, I - cosΘ_2 ≤sinΘ_2. Together with (<ref>), we have 

    U_r-_rZ^⊤_2=U_rZ - _r _2 ≤ 2 sinΘ_2.


3. Ensuring  (Υ_f)=r. 
Choose M = Υ_f, which has rank r, and M+E = Υ_fZ^⊤. We now show that M+E also has rank r. 
Using Weyl's theorem on perturbation of singular values <cit.>,

    |σ_r(Υ_fZ^⊤) - σ_r(Υ_f)| ≤   E_2 = Υ_f-Υ_fZ^⊤_2 ≤U_r-_rZ^⊤_2 ≤ 2 sinΘ_2 .

We have used the fact that Υ_f and Υ_fZ^⊤ are submatrices of U_r and _rZ^⊤, respectively. Since σ_r(Υ_f) = 1/Υ_f^†_2, we can rearrange to get

    σ_r(Υ_fZ^⊤) ≥1/Υ_f^†_2 - 2sinΘ_2  > 0
 
since by assumption, η = 2 sinΘ_2Υ_f^†_2 < 1. This shows that both Υ_f and Υ_fZ^⊤ have rank r. Therefore, by (<ref>), (Υ_fZ^⊤)^† = ZΥ_f^†. Furthermore, (<ref>) applies and

    ZΥ_f^†_2 ≤Υ_f^†_2/1-2 sinΘ_2Υ_f^†_2.
 



4. Error in A_r - Z_rZ^⊤_2. Write both matrices in terms of their factors, and add and subtract the term Υ_f^†Υ_lZ^⊤ to get

    A_r -  Z_rZ^⊤_2    = Υ^†_fΥ_l - Υ_f^†Υ_lZ^⊤ + Υ_f^†Υ_lZ^⊤ -ZΥ_f^†Υ_lZ^⊤_2 
       ≤Υ^†_f _2Υ_l-Υ_lZ^⊤_2 + Υ^† - ZΥ_f^†_2Υ_lZ^⊤_2.

Note that, by Assumption A2, (Υ_f) = r. In step 3, we showed that (Υ_fZ^⊤) = r.  
As before, with M = Υ_f and M+E = Υ_fZ^⊤, (<ref>) applies and

    Υ_f^† - ZΥ_f^†_2 ≤√(2)Υ_f^†_2 ZΥ_f^†_2 Υ_f - Υ_fZ^⊤_2.

Plugging this into (<ref>), we have

    A_r -  Z_rZ^⊤_2 ≤Υ_f^†_2 (Υ_l-Υ_lZ^⊤_2  +  √(2)ZΥ_f^†_2 Υ_f - Υ_fZ^⊤_2).

We have used the fact that Υ_lZ^⊤_2 ≤_rZ^⊤_2 = 1; the inequality is because Υ_lZ^⊤ is a submatrix of _rZ^⊤ and the equality follows since _rZ^⊤ has orthonormal columns. Similarly,

    Υ_l-Υ_lZ^⊤_2 ≤U_r-_rZ^⊤_2,     Υ_f-Υ_fZ^⊤_2 ≤U_r-_rZ^⊤_2,

so that 

    A_r -  Z_rZ^⊤_2 ≤Υ_f^†_2 U_r-_rZ^⊤_2( 1  +  √(2)ZΥ_f^†_2 ).



5. Finishing the proof. Plug (<ref>) and (<ref>) into (<ref>) to obtain

    A_r -  Z_rZ^⊤_2  ≤  (2sinΘ_2 Υ_f^†_2)(1 +  √(2)Υ_f^†_2/1-2 sinΘ_2Υ_f^†_2).

Plug this bound into (<ref>) to finish the proof. 



An important aspect of the proof of Theorem <ref> is the choice of the orthogonal matrix Z that determines the similarity transformation. Our proof used this idea from <cit.>. It is worth pointing out that the matrix Z = PQ^⊤ is also the solution to the orthogonal Procrustes problem 

    min_Z∈R^r× r
    Z^⊤Z = I_r U_rZ - _r_F.

Therefore, the matrix Z = PQ^⊤ is the matrix that “best rotates” the columns U_r to align with the columns of _r. Note, however, that in the proof we use the 2-norm instead of the Frobenius norm. 




 §.§ Accuracy of the singular vectors
 Theorem <ref> shows that the accuracy of the approximation algorithms to the ERA depend on the canonical angles between the exact and approximate subspaces (U_r) and (_r). Insight into the accuracy between these subspaces can be obtained by standard perturbation theory. Suppose there are numbers α≥ 0 and δ >0 such that 

    σ_r+1(_s) ≥α + δ,     σ_r(H_s) ≤α.

Let _r,_r be the matrices containing the left and right singular vectors corresponding to the first r singular values of _s, which are taken to be the diagonals of the matrix _r. Then, by <cit.>,

    sinΘ_2 ≤max{H_s_r- _r_r _2, H_s^⊤_r- _r _r_2}/δ.

In Section <ref>, the randomized subspace iteration is used to construct the low-rank approximation _s. We can use more special purpose error bounds to quantify the accuracy of the singular vectors. Specifically, one may use techniques developed in Section 3.2 of <cit.>. We omit a detailed statement of the results.




 §.§ Stability

The discrete dynamical system with the system matrix A_r is stable if the eigenvalues of A_r lie inside the unit circle; that is, the spectral radius ρ(A_r) < 1. Assumption A3 in Section <ref> means that the Markov parameters decay after some finite time. Then, by <cit.>, there exists a positive integer s such that the identified reduced order model is a stable discrete dynamical system. Theorem <ref> can also be used to determine when the approximate discrete dynamical system with the system matrix _r is stable. 

We consider the spectral radius of _r, and write

    ρ(_r) = max_1 ≤ j ≤ n|λ_j(_r)| = |λ_j^*(_r)|,

for some index j^*. Similarly, let i^* be the index such that λ_i^*(A_r) is the closest eigenvalue to λ_j^*(_r). If i^* and j^* are not unique, break the tie in some fashion. Therefore, we have the series of inequalities

    ρ(_r) ≤ |λ_j^*(_r) - λ_i^*(A_r) | + | λ_i^*(A_r)| ≤(Ψ(_r),Ψ(A_r)) + ρ(A_r).

By Theorem <ref>, the approximate discrete dynamical system is stable if 

    κ_2(W) η( 1 +√(2)Υ_f^†_2/1-η) < 1 - ρ(A_r).

By assumption, η < 1; the above equation gives an additional condition on η for stability. This condition can be informative for the numerical method used to compute the singular vectors _r. 



§ NUMERICAL RESULTS


To test our algorithms, we consider two different test problems which pose distinct challenges. The first test problem is generated from a LTI model of a controlled heat transfer process for optimal cooling of a steel profile <cit.>. The second test problem is generated from a LTI state-space model of an electrical power generation system with 50 generators <cit.>. The dimensions of the variables related to both test problems are shown in Table <ref>. The heat transfer test problem has a very large system size n but a small number of inputs and outputs, while the power system test problem has a relatively small system size but a much larger number of inputs and outputs.  
All our experiments were run in MATLAB, and our timing experiments were run on the NCSU HPC Cluster with 72GB memory.





 §.§ Heat Transfer

We first test our algorithms on the heat transfer test problem, which has a large system size but a relatively small number of inputs and outputs. This system also has slow dynamics, so a large value of s is needed to capture the transient behavior. The model for this test problem is not in standard form, so we will need to convert it to standard form before we are able to use our algorithms. The original continuous form after spatial discretization is

    Edx/d t   = A_ox + B_ou
    y   = C_ox + D_o u.


To convert it to standard form, we use the Cholesky factorization of E = LL^⊤. Then, the continuous standard form system matrices are found by computing A_c = L^-1A_oL^-⊤, B_c = L^-1B_o, and C_c = C_oL^-⊤. Note that D_c = D_o remains unchanged. The continuous-time matrices are converted into appropriate discrete-time matrices using MATLAB's  command. It is worth noting that the continuous-time matrices are only used in the construction of the Markov parameters but are not used explicitly in the system identification algorithms.

For this heat transfer problem, we will compare accuracy and speed of our algorithms that take advantage of the block Hankel structure of H_s, i.e., RandSVD-H (see Section <ref>), to the standard ERA that uses a full SVD. We also contrast these two algorithms with another algorithm we call SVDS-H, which uses the MATLAB command , but we instead use the block Hankel multiplication described in Algorithm <ref>.



  §.§.§ Accuracy
 We consider several numerical experiments to test the accuracy of our algorithms. For each experiment to test accuracy, we used s = 1000 and a target rank of r=20. That is, we are performing model reduction in addition to the system identification. For RandSVD-H, we used q=1 for the number of subspace iterations and ρ = 20 for the oversampling parameter.  In our first experiment, we compare the first 20 singular values of the block Hankel matrix H_s ∈R^6000 × 7000 as computed by a full SVD to the first 20 singular values computed by the SVDS-H and RandSVD-H algorithms.  The singular values are plotted in Figure <ref>; we see that they are in good agreement with each other.



Next we compare the eigenvalues of the identified A_r matrices using the SVDS-H and RandSVD-H algorithms to the eigenvalues of A_r identified using the full SVD ERA algorithm.  The results are shown in Figure <ref>, where we can see that the eigenvalues of A_r computed using SVDS-H and RandSVD-H match those of A_r computed using the full SVD  accurately. 



In order to quantitatively evaluate the accuracy of the identified eigenvalues, we consider the Hausdorff distance metric, defined between the spectra of the two matrices A_r and A_r, 


    h_d(ψ(A_r),ψ(A_r)) = max{(ψ(A_r),ψ(A_r)), (ψ(A_r),ψ(A_r)) }.

Note that the Hausdorff distance is related to the spectral variation, defined in (<ref>), as it measures how close two sets are to each other. We use h_d instead of spectral variation as the spectral variation of two sets depends on the order of the sets in question, while the Hausdorff distance does not.
For both algorithms SVDS-H and RandSVD-H, we computed the Hausdorff distance between the spectra of A_r and A_r and A_r and we saw that h_d ≈  10^-14. 

The final measure of accuracy we will consider for this test problem is the relative error in the Markov parameters. Let A_r,B_r,C_r be the system matrices identified by the full SVD, and let A_r,B_r,C_r be the system matrices identified using another method. The relative error in the Markov parameters is then defined as 

    M_k = C_rA_r^kB_r - C_rA_r^k B_r_2/C_rA_r^kB_r _2,    k = 1,…,2s-1.

 We compute this error for the identified matrices using the SVDS-H and RandSVD-H algorithms, and plot the results in Figure <ref>. Both algorithms produce comparable error in the Markov parameters, and this error is very low in both cases. These results combined with those from Figure <ref> show that our algorithms are very accurate compared to the standard algorithms.






  §.§.§ Timing
 We now consider the average run time of the three algorithms on the heat transfer problem. The results, shown in Figure <ref>, clearly indicate that using the block Hankel structure as in SVDS-H and RandSVD-H significantly reduces the computational cost, and using a randomized algorithm reduces that cost further, as RandSVD-H is computationally the least expensive of the three. We also see from the plots that the SVDS and RandSVD-H algorithms, which both exploit the block Hankel structure, show approximately linear scaling with s confirming the analysis of the computational costs. 






  §.§.§ Comparison to CUR-ERA
 The final test we consider with this application is how CUR-ERA compares to the randomized algorithms we analyzed in the previous experiments. For CUR-ERA, we used the default parameters suggested in the code provided by the authors of <cit.>. The maxvol tolerance is 0.002, and the iterations stop when the relative error of the subsequent iterations is less than 10^-4. For our experiments, we start with the accuracy of the eigenvalues, plotted in Figure <ref>. 

In this figure, we see that the eigenvalues of A identified by CUR-ERA are not as accurate as those identified by RandSVD-H. To quantify this, we compute the Hausdorff distance between the spectrum of the A_r matrix identified by CUR-ERA and the spectrum of the A_r matrix identified by SVD-ERA when s = 1000. We find that this distance is approximately 8.9 × 10^-3, supporting the claim that CUR-ERA is less accurate than RandSVD-H. At the same time, the CUR-ERA algorithm took 0.83 seconds compared to 0.66 seconds for RandSVD-H (averaged over three runs). Finally, we show the relative error in the Markov parameters, computed as defined in (<ref>), for CUR-ERA, RandSVD-H, and SVDS-H. The results are shown in Figure <ref>, and we see that the error shown here confirms the accuracy comparison, thus far. The accuracy of CUR-ERA  is lower than that obtained from either RandSVD-H or SVDS-H.








 §.§ Power system

For our next application, we consider the system identification of an electric power system model. Identification is an important component of power system modeling as the operating points in a power grid keep changing due to changes in loads, renewable generation profiles, and the network topology, all of which motivate the need for fast algorithms that can quickly update the small-signal model of the grid about these changing operating points from time to time <cit.>. The updated models, in turn, enable operators to make real-time control decisions for ensuring system stability.  For the test problem, we consider a prototype model of reasonably large dimension, namely the IEEE 50-generator model with 145 buses, n=155 state variables, and m = 50 input variables, as listed in Table <ref>. Details about the physical meaning of the states and model parameters of this test system can be found in <cit.>. The original IEEE 50-generator model is nonlinear. The model is, therefore, first linearized about a chosen power flow solution. The resulting LTI model is excited by 50 impulse functions injected through the excitation input terminals of the generators. All ℓ = 155 states are measured, and are considered as outputs. The output matrix for this example is, therefore, the identity matrix.

An important feature of this problem is that it has a significantly larger number of inputs and outputs with a relatively small system size n=155. Since the computational cost of RandSVD-H scales linearly with the product of the number of inputs and outputs, this suggests that RandTERA may have be computationally advantageous if the number of inputs and outputs can be reduced. Therefore, our experiments will compare RandTERA to RandSVD-H both in terms of accuracy and computational costs. For this problem, we again perform system identification as well as model reduction, taking the target rank to be r = 75.



  §.§.§ Accuracy
 First, we will examine the accuracy of both RandTERA and RandSVD-H on the power system dataset. For all the experiments for testing accuracy, we will use s=500, and target rank r=75. For the randomized algorithms, we use an oversampling parameter ρ = 20, and q = 1 subspace iterations. Note that the full Hankel matrix for the accuracy experiments is of size 77500 × 25000.





  
Choice of ℓ' and m' In our first experiment,  we will investigate the choice of the reduced number of inputs and outputs ℓ' and m' for the tangential algorithms (TERA and RandTERA). To choose these parameters, we compute the singular values of H_w and H_e defined in (<ref>),(<ref>). The indices ℓ' and m' are chosen to be the largest integers such that


    ℓ' =    min_k {1 ≤  k ≤ℓ|σ_k(H_w) ≥ϵσ_1(H_w) . }
     m' =    min_k {1 ≤ k ≤ m |σ_k(H_e)≤ϵσ_1(H_e). }.

That is, they are chosen to be the smallest indices for which the singular values are within a threshold ϵ of the largest singular value. 
The values of ℓ' and m' for different ϵ thresholds are shown in Table <ref>.





  
Accuracy of eigenvaluesWe first compare the accuracy of the singular values of the computed block Hankel matrix H_s. First, we show in Figure <ref> the effect of different ϵ values on the accuracy of singular values computed by TERA and RandTERA. We can see that the smaller the tolerance ϵ, the more accurate the singular values. 
Also, the singular values computed using TERA and RandTERA are  less accurate compared to RandSVD-H, but the singular values computed by RandTERA are comparable with those of TERA. This shows that the randomized algorithms are accurate compared to their deterministic counterparts.






Next, we consider the eigenvalues of the identified A matrix using the RandSVD-H and RandTERA algorithms compared with the eigenvalues of the original A matrix. We use the parameters listed above, but we fix ϵ = 10^-2. These results are shown in Figure <ref>. We can see that our algorithms identify the state matrix with eigenvalues very close to the true eigenvalues. To quantify how close they are, we consider the Hausdorff distance again as defined in (<ref>). We computed the Hausdorff distance between the spectrum of the A matrix identified by an SVD-ERA and the spectra of the A matrices identified by RandSVD-H and RandTERA. We used s=500 for all algorithms, and saw that for RandSVD-H, h_d ≈ 2.7 × 10^-4. For RandTERA, h_d ≈ 3.0 × 10^-3. 


Our final measure of accuracy is the relative error in the Markov parameters as defined in (<ref>). We plot this relative error  in Figure <ref>, and make two observations. The first is that RandSVD-H has a much lower relative error than either TERA or RandTERA. The second is that the error from RandTERA matches the error from TERA very accurately, so the randomization is not causing a significant loss in the accuracy. We can see in Figure <ref>, that the relative error increases with increasing k. This is because the identification techniques using impulse response fail to identify the well damped poles (i.e., the poles whose real part is significantly lesser than 1) accurately, since information corresponding to these modes have to be obtained from the first few Markov parameters only.  In Figure <ref>, we compare the impulse response of the original full, discrete-time model with that of the models identified using RandSVD-H and RandTERA (an impulse input was fed into Generator 1). The top plot shows the speed response of Generator 1 while the bottom plot shows that for Generator 2. We see that both the identified models have similar performance for Generator 1 and match the original trajectory well, whereas in Generator 2, the trajectories from all  algorithms match each other closely, with only minor deviations from the “true” trajectory in the the first 2 seconds.
 The results for the other generators were observed to follow similar matching trends.





  §.§.§ Timing
 
We next examine the computational cost of RandSVD-H and RandTERA on the power system test problem. We show in Table <ref> the average CPU time of both of these algorithms compared to a full SVD as s increases. As before, we averaged over three different runs, and used the same parameters for the randomized algorithms, r=75, ρ = 20 and q=1. We see that RandTERA is much cheaper compared to RandSVD-H, and both are much cheaper than a full SVD as expected. Note that for s = 1000, the matrix size was 155,000 × 50,000 and since this matrix was too large to store explicitly, we do not report the computational cost using the SVD.






§ CONCLUSIONS AND FUTURE WORK

We have presented two different randomized algorithms for accelerating the computational cost of system identification using ERA. The first algorithm accelerates the randomized subspace iteration by efficiently computing matrix vector products using the block Hankel structure. The second algorithm uses the previous algorithm, but on a matrix with a reduced number of inputs and outputs using tangential interpolation. The algorithms no longer have the cubic scaling with s, which controls the number of time steps. The error analysis relates the accuracy of the eigenvalues of the  identified system matrix to the accuracy of the singular vectors of the block Hankel matrix H_s. Numerical experiments from different applications show that our algorithms are both accurate and computationally efficient. 


There are several avenues for future work. First, the analysis in Section <ref> is general and does not make assumptions on the specific algorithm used to compute the low-rank approximation to H_s. The analysis can be specialized by using specific results for the accuracy of the singular vectors, for example, following the approaches in <cit.>. Second, the paper assumes that we have access to the Markov parameters. This is not the case when it is possible to excite the system using general inputs. One possibility, as mentioned earlier, is to use the approach in <cit.> which may not be feasible when a large number of Markov parameters need to be estimated. Algorithms such as MOESP and N4SID (see <cit.>) can be used for general inputs but also suffer from high computational costs. However, we can combine the block Hankel structure with randomized algorithms for the general input case as well. Such techniques which involve identification using general input can be extremely accurate in identifying the well-damped poles as they constant excite the system with a persistently exciting signal. This is an ongoing research direction. 



§ ACKNOWLEDGEMENTS

We would also like to thank Eric Hallman for useful comments and for generously sharing his code with us.
siamplain



