Interplay of Sensor Quantity , Placement and System Dimension in POD-based Sparse Reconstruction of Fluid Flows

Sparse recovery of fluid flows using data-driven proper orthogonal decomposition (POD) basis is systematically explored in this work. Fluid flows are manifestations of nonlinear multiscale PDE dynamical systems with inherent scale separation that impact the system dimensionality. Given that sparse reconstruction is inherently an ill-posed problem, the most successful approaches require the knowledge of the underlying basis space spanning the manifold in which the system resides. In this study, we employ an approach that learns basis from singular value decomposition (SVD) of training data to reconstruct sparsely sensed information. This results in a set of four control parameters for sparse recovery, namely, the choice of basis, system dimension required for sufficiently accurate reconstruction, sensor budget and their placement. The choice of control parameters implicitly determines the choice of algorithm as either l2 minimization reconstruction or sparsity promoting l1 norm minimization reconstruction. In this work, we systematically explore the implications of these control parameters on reconstruction accuracy so that practical recommendations can be identified. We observe that greedy-smart sensor placement provides the best balance of computational complexity and robust reconstruction for marginally oversampled cases which happens to be the most challenging regime in the explored parameter design space.


Introduction
Multiscale fluid flow phenomena in engineering and geophysical settings are invariably data-sparse, i.e. there are more scales to resolve than there are sensors.A major goal is to recover more information Sparse reconstruction (SR) is inherently ill-posed, under-determined inverse problem where the number of constraints (i.e., sensor quantity) are much less than the number of unknowns (i.e., high resolution field).However, if the underlying system is sparse in a feature (space of basis coefficients) space then the probability of recovering a unique solution increases by solving the reconstruction problem in this lower-dimensional space.The core theoretical developments of such ideas and their first practical applications happened in the realm of image compression and restoration [10,21].
Sparse Reconstruction with Data-driven Basis: Data recovery techniques based on Karhunen-Loeve (K-L) procedure with least squares (l 2 ) error minimization , also known as Gappy POD or GPOD [12,13,17], was originally developed in the nineties to recover marred faces [17] in images.The fundamental idea is to utilize the POD basis computed offline from the data ensemble to encode the reconstruction problem into a low-dimensional feature space.This way, the sparse data can be used to recover the sparse unknowns in the space of POD coefficients by minimizing the l 2 errors.If the data-driven POD basis are not known a priori, an iterative formulation [12,17] to successively approximate the POD basis and the coefficients was proposed with limited success.While this approach has been shown to work in principle [12,14,22], it is prone to numerical instabilities and inefficiency.Advancements in the form of a progressive iterative reconstruction framework [14] are effective, but impractical for real-time application.In fact, all the aforementioned issues are related to the POD-basis being data-driven and therefore, can represent the data optimally but not easily generalizable.This requires that they be known a priori -a stringent requirement for many practical systems as training data is rarely available and even when it is, it may not effectively (i.e.relevant to) the prediction regime.Such limitations make data-driven basis hard to use for practical sparse sensing-based reconstruction.Nevertheless, they find tremendous value in data-driven modeling (machine learning, Koopman operator models [23,24]) applications and nonlinear model order reduction [4] of systems that are statistically stationary where training data is available and tends to remain relevant to the system as it evolves.
Sparse Reconstruction with Generic Basis: A way to bypass the above limitations is to use generic basis such as wavelets [25] or Fourier functions.Such choices assume that most practical systems are relatively sparse with respect to these basis choices although not optimal.This is particularly true for image processing applications but may not be optimal for fluid flows whose dynamics obey PDEs (and may include sharp gradients).While avoiding the cost of computing the basis offline, such approaches run into sparsity issues as the basis do not optimally encode the underlying dynamical system.In other words, larger the basis, more sensors are needed for complete and accurate reconstruction.Thus, once again the reconstruction problem is ill-posed because the basis space is not sufficiently low-dimensional resulting in the flow being under sampled.To bypass this, one needs to look for a sparse solution which is usually not realized through l 2 error minimization approaches.The magic of Compressive Sensing (CS) [8][9][10][11] is in its ability to overcome this constraint by seeking a solution that can be less sparse than the dimensionality of the chosen feature space using l 1 -norm regularized least-squares reconstruction.Such methods have been successfully applied in image processing using Fourier or wavelet basis and also to fundamental fluid flows [6,7,[26][27][28][29].Compressive sensing essentially looks for a regularized sparse solution using l 1 norm minimization of the sparse coefficients by solving a computationally manageable convex optimization problem.Such sparsity promoting l 1 regularized reconstruction can also be combined with data-driven POD basis to overcome the aforementioned limitations.Successful attempts at this include reconstruction of sparse PIV data [6] and pressure measurements around a cylinder surface [7].To bypass the limitations from using data-driven POD basis for SR, one often tends to build a library of POD modes from precursor training simulation or measurement data over a range of flow parameters (e.g.Reynolds numbers) and leverage l 1 regularized reconstruction to identify the most relevant K-sparse solution.Building such a library of POD modes is a one-time investment for a chosen class of flows.Such a framework has been attempted in [7] where POD modes from simulations over a range of Reynolds (Re) numbers of a cylinder wake flow were used to populate a library of bases and then used to classify the flow regime based on sparse measurements.
In order to carry out a systematic analysis of interplay between the various SR control parameters, namely, desired reconstruction dimension, sensor budget and their placement, we limit ourselves to primarily l 2 minimization approaches (except for comparative analysis where we use l 1 ) and four different sensor placement methods including random sampling (as used in [30]), QR with column pivoting [31,32], Discrete Empirical Interpolation Method (DEIM) [18] and explicit condition number minimization [13].The suitability of these choices are numerically explored in a discretized parameter space of sensor budget and system dimension.The rest of manuscript is organized as follows.In section 2 we review the basics of sparse reconstruction theory (sec.2.1), computation of data-driven POD basis (sec.2.2), role of measurement locations (sec.2.3), algorithms for sensor placement (sec.2.4) and sparse recovery (sec.2.5).In section 3 we discuss how the training data is generated and also summarize the outcomes form sensor placement using the different approaches.In section 4 we discuss the results from our analysis of the SR of the cylinder wake flow using both POD and ELM basis.This is followed by a summary of major conclusions from this study in section 5.

Formulating the Sparse Reconstruction Problem
Given a high resolution data representing the state of the flow system at any given instant denoted by x ∈ R N , its corresponding sparse representation given by x ∈ R P with P N .Then, the sparse reconstruction problem is to recover x, when given x along with information of the sensor locations in the form the measurement matrix C ∈ R P ×N as shown in eqn.(1).The measurement matrix C determines how the sparse data x is collected from x. Variables P and N are the number of sparse measurements and the dimension of the high resolution field, respectively.
In this article, we focus on vectors x that have a sparse representation in a basis space Φ ∈ R N ×K such that K N and yielding x = Φa.Naturally, when one loses the information about the system, the recovery of said information is not absolute as the reconstruction problem is ill-posed, i.e., more unknowns than equations in eqn.(1).Thus the most straightforward approach to recover x is by computing the inverse of C using a least-squares error minimization procedure as shown in eqn.(2).However, this does not result in stable solutions as the system is ill-posed (under-determined). (2)

Sparse Reconstruction Theory
The theory of sparse reconstruction has strong foundations in the field of inverse problems [33] with applications in diverse fields of study such as a geophysics [34,35] and image processing [36,37].In this section, we formulate the reconstruction problem which has been presented in CS literature [8,10,[38][39][40].Many signals tend to be "compressible", i.e., they are sparse in some K-sparse basis Φ as shown below: where Φ ∈ R N ×N b and a ∈ R N b with K non-zero elements.In the sparse reconstruction formulation above, Φ ∈ R N ×N b is used instead of Φ ∈ R N ×K as the K most relevant basis vectors for a given data are not identified a priori.Consequently, a more exhaustive basis set of dimension N b ≈ P > K is typically employed.To recover N -dimensional data, one can atmost use N basis vectors, i.e., N b ≤ N .
In practice, the dimension of the candidate basis space need not be N and can be represented by N b N as only K(≤ N b ) of them are needed to represent the acquired signal up to a desired quality.This is typically the case when Φ is composed of optimal data-driven basis vectors such as POD modes.The reconstruction problem is then recast as identification of these K sparse coefficients.
In many practical situations, Φ and K along with N b , N are input by the user.Standard transform coding [25] practice in image compression involves collecting a high resolution sample, transforming it to a Fourier or wavelet basis space where the data is sparse, retain the K most relevant coefficients while discarding the rest of the information.This is the basis of JPEG and JPEG-2000 compression standards [25].The sample and then compress mechanism still requires acquisition of high resolution samples and processing them before reducing the dimensionality.This is highly challenging as handling large amounts of data is difficult in practice due to demands on processing power, storage, and time.
Compressive sensing [8,10,[38][39][40] focuses on direct sparse sensing based inference of the K-sparse coefficients by essentially combining the steps in equations 1 and 3 as below: where Θ ∈ R P ×N b is the map between the basis coefficients a that represent the data in a feature space and the sparse measurements, x in physical space.The challenge in solving for x using the under determined eqn.( 1) is that C is ill-conditioned and x in itself is not sparse.However, when x is sparse in Φ, the reconstruction using Θ in eqn.(4) becomes practically feasible (for P K) by solving for a.Thus, one effectively seeks a K-sparse a with P constraints (given by x) using established methods from linear algebra and constrained optimization.

Case 1: For
For the over determined system with P > K = N b , a is estimated using a regularized least squares solution based on the normal equation as a = (Θ) L+ x = Θ T Θ + λI −1 Θ T x for a chosen λ.This is obtained by minimizing the appropriate cost function given by J cost = x − Θa 2 2 + λ a 2 2 .This regularized least-squares solution procedure for the overdetermined case is nearly identical to the original GPOD algorithm developed by Everson and Sirovich [17] if Φ is chosen as the POD basis.However, x in GPOD contains zeros as placeholders for all the missing elements whereas the above formulation retains only the measured data points.
When P ≤ K = N b , and the system is under-determined with non-unique solutions, one looks for a minimum norm reconstruction of a.This is achieved by minimizing the corresponding s-norm of a (s chosen appropriately) and x is then recovered from eqn. (3).s chosen as 2 yields the minimum l 2 norm reconstruction of x by penalizing the larger elements of a.The l 2 -regularized method finds an a that minimizes the expression shown in eqn.(5).
One can solve for T and a in eqn.( 5) using Lagrange multipliers approach to yield a solution that is a right pseudo-inverse of Θ as in eqn.( 6) below as long as Θ has a full rank of P : 2.1.2.Case 2: For K < N b When K N b , one typically looks for a sparse solution of a.The l 2 norm minimization and l 2 regularization approaches discussed above provide numerical stability, but do not produce sparse solutions.A natural way to enhance sparsity of a is to minimize a 0 , i.e., minimize the number of non-zero elements such that Θa = x is satisfied.It has been shown [41] that with P = K +1 (P > K in general) independent measurements, one can recover the sparse coefficients with high probability using minimum l 0 norm reconstruction.This condition can be heuristically interpreted as each measurement needing to excite a different basis vector φ i so that its coefficient a i can be optimally identified.If two or more measurements excite the same basis φ j then additional measurements may be needed to produce acceptable reconstruction.On the other hand, for P ≤ K independent measurements, the probability of recovering the sparse solution is highly diminished.Nevertheless, l 0 -reconstruction is a computationally complex, np-complete and poorly conditioned problem with no stability guarantees.
The popularity of compressed sensing arises on account of the theoretical advances [42][43][44][45] guaranteeing near-exact reconstruction of the uncompressed information by solving for the K sparsest coefficients using l 1 norm minimization methods.The l 1 reconstruction is a relatively simpler convex optimization problem (as compared to l 0 ) and solvable using linear programming techniques for basis pursuit [8,38,46] and shrinkage methods [47].Theoretically, one can perform the simplistic brute force search to locate the largest K coefficients of a that satisfy eqn.(7), but the computational effort increases exponentially with dimension.To overcome this burden, a host of greedy algorithms [9,11,48] have been developed to solve the l 1 norm minimization problem in eqn.(7) with complexity O(N 3 ) for N b ≈ N .However, this approach requires P > O(Klog(N b /K)) measurements [8,38,42] to reconstruct the K-sparse vectors with high probability.The schematic illustration of both l 2 and l 1 -based formulations are presented in Figure 1.
Solving the l 1 minimization problem shown in eqn.( 7) is complicated relative to the l 2 minimization solution described in eqn.(5).This is because, unlike the cost function to be minimized in eqn.(5), the cost function in eqn.(7) is not differentiable when a i = 0 which necessitates an iterative solution instead of a closed form solution.Further, the minimization of the l 1 cost function is also an unconstrained optimization problem that is commonly converted into a constrained optimization as shown in eqn.(8).Here t is a user defined sparsity knob to 'shrink' the coefficients.
This constrained optimization in eqn.( 8) is quadratic in a and therefore, a quadratic programming problem with the feasible region bounded by polyhedron (in the space of a).There exists two classes of l 1 solution methodologies: (i) least absolute selection and shrinkage operator or LASSO [47] and (ii) basis pursuit denoising [46].LASSO and its variant essentially convert the constrained formulation into a set of linear constraints.Recently popular approaches include greedy methods such as optimal matching pursuit(OMP) [11,27] and interior point methods [49].An intuitive iterative sequential least-squares thresholding framework is used by Brunton et al. [50].The idea here is to achieve 'shrinkage' by repeatedly zeroing out the coefficients smaller than a given choice of hyperparameter.In summary, the reconstruction framework is characterized by three parameters, N b , K, P .N b is the candidate basis dimension employed for this reconstruction which is bounded by N , i.e.N b N .K represents the desired reconstruction sparsity and chosen such that if these features are predicted accurately, then the achieved reconstruction meets a desired quality.The more sparse a system, the smaller K is for a desired reconstruction quality and is specified by the user based on knowledge of the system.P represents the available sensor quantity as input to the problem.N is the preferred dimension of the reconstructed state.The interplay of N b , K, and P determine the choice of algorithm employed, i.e., whether the reconstruction is based on least squares minimization, l 2 norm minimization or sparsity enabling l 1 approaches as summarized in Table 1.
Table 1.The choice of sparse reconstruction algorithm based on problem design using parameters P (sensor sparsity), K (targeted reconstruction sparsity) and N b (candidate basis dimension).
All of the above sparse recovery estimations are conditional upon the measurements (rows of C) being incoherent with respect to the sparse basis Φ.This is usually accomplished by using a random choice of sensor placement, especially when Φ is made up of Fourier functions or wavelets.If the basis functions Φ are orthonormal, such as wavelet or POD basis (with inherent hierarchy), one can discard the majority of the small coefficients in a (setting them as zeros) and still retain reasonably accurate reconstruction.The mathematical explanation for this has been previously shown in [10].However, it should be noted that incoherency is a necessary, but not sufficient condition for exact reconstruction.Exact reconstruction requires optimal sensor placement to capture the most information for a given flow field.In this study, we assess the interplay between sensor quantity, quality and desired system sparsity for POD-based sparse reconstruction.

Data-driven Basis Computation using POD
In the SR framework, basis such as POD modes, Fourier functions, and wavelets [8,10] can be used to generate low-dimensional representations for both l 2 and l 1 -based methods.While an exhaustive study on the effect of different choices on reconstruction performance is potentially useful, in this study we focus on POD-based SR.A similar effort has been reported in [6] where a comparison between discrete cosine transform and POD bases was performed.
Proper orthogonal decomposition (POD), also known as Principal Components Analysis (PCA) or Singular Value Decomposition (SVD), is a dimensionality reduction technique that computes a linear combination of low-dimensional basis functions (POD modes) and weights (POD coefficients) from snapshots of experimental or numerical data [2,4] through eigen-decomposition of the spatial correlation tensor of the data.It was introduced in the turbulence community by Lumley [51] to extract coherent structures in turbulent flows.The resulting singular vectors or POD modes represent an orthogonal basis that maximizes the energy capture from the flow field.For this reason, such eigenfunctions are considered optimal in terms of energy capture and other optimality constraints are theoretically possible.Taking advantage of the orthogonality, one can project these POD basis onto each snapshot of data in a Galerkin sense to deduce coefficients that represent evolution over time in the POD feature space.The optimality of the POD basis also allows one to effectively reconstruct the full field information with knowledge of very few coefficients, a feature that is attractive for solving sparse reconstruction problems such as in eqn.(4).However, this is contingent on the spectrum of the spatial correlation tensor of the data having sufficiently rapid decay of the eigenvalues, i.e. it supports a low-dimensional representation.This is typically not true in the case of turbulent flows with very gradual decay of energy across the singular values.Further, in such dynamical systems, the small scales with low-energy can still be dynamically important and will need to be reconstructed, thus requiring significant number of sensors.
Since the eigen-decomposition of the spatial correlation tensor of the flow field requires handling a system of dimension N , it requires significant computational expense.An alternative method is to compute the POD modes using the method of snapshots [52] where the eigen-decomposition problem is reformulated in a reduced dimension (assuming the number of snapshots in time is smaller than the spatial dimension) framework as summarized below.Consider that X ∈ R N ×M (different from x ∈ R N ) is the full field representation with only the fluctuating part, i.e., the temporal mean is taken out of the data.N is the dimension of the full field representation and M is the number of snapshots.The procedure involves computation of the temporal correlation matrix CM as: The resulting correlation matrix CM ∈ R M ×M is symmetric and an eigendecomposition problem can be formulated as: where the eigenvectors are given in Typically, both the eigenvalues and corresponding eigenvectors are sorted in descending order such as λ 1 > λ 2 > ... > λ M .The POD modes Φ and coefficients a can then be computed as One can represent the field X as a linear combination of the POD modes Φ as shown in eqn.(3) and leverage orthogonality to compute the Moore-Penrose pseudo-inverse, i.e., Φ † = Φ T , and compute the POD coefficients, a ∈ R M ×M as shown in eqn.(12), It is worth mentioning that subtracting the temporal mean from the input data is not critical to the success of this procedure.In fact, retaining the mean of the data during SVD computation generates an extra mean mode which modifies the energy spectrum.For the low-dimensional cylinder wake used in this study,the dominant mode when performing SVD with mean captures 98% energy whereas for the SVD without mean, the dominant mode captures 49% energy.Based on our experience over the course of this research study, these choices have minimal impact on the reconstruction performance.In fact, during practice, one does not use POD basis with the mean removed as the sparse data includes the mean that is not known a priori.
Using the snapshot procedure for the POD/SVD computation fixes the maximum number of POD basis vectors to M which is typically much smaller than the dimension of full state vector, N .If one wants to reduce the dimension further, then a criterion based on energy capture is devised so that the modes carrying the least amount of energy truncated to dimension K < M .For many common fluid flows, using the first few POD modes and coefficients are sufficient to capture almost all the relevant dynamics.However, for turbulent flows with large-scale separation, a significant number of POD modes will need to be retained.

Measurement Locations, Data Basis and Incoherence
Recalling from subsection 2.1, the reconstruction performance is strongly tied to the measurement matrix, C being necessarily incoherent with respect to the low-dimensional basis Φ [10] and this is usually accomplished by using random elements to populate the matrix.In practice, one can adopt two types of random sampling of the data, namely, single-pixel measurement [7,53,54] or random projections used in compressive sensing [8,10,39].Typically, single-pixel measurement refers to measuring information at chosen spatial locations such as measurements by unmanned aerial systems (UAS) in the atmospheric fields or point probes.The resulting matrix C is structured as C ← [e 1 , e 2 , ..., e p ] T , where e p is column vector with zeros except at the index p where it assumes a value of 1.
Another popular choice of measurements employed in the compressive sensing or image processing community is random projections where the measurement matrix is populated using random numbers based on a chosen distribution (Gaussian, Bernoulli) on to which the full state data is projected.In theory, the random matrix is highly likely to be incoherent with any generic basis [10], and hence suitable for sparse recovery purposes.However, for most fluid flow and other practical sensing applications, the sparse data is usually sourced from point measurements.Irrespective of the approach adopted, the measurement matrix C (eqns.( 1)-( 2)) and basis Φ should be incoherent to ensure optimal sparse reconstruction.This essentially implies that one should have sufficient measurements distributed in space to excite the different modes relevant to the data being reconstructed.Mathematically, this is related to CΦ being full rank and invertible.There exist metrics to estimate the extent of coherency between C and Φ in the form of an coherency number, µ as shown in eqn.(13) [55], where c i is a row vector in C (i.e.c i = e j ) and φ j is a column vector of Φ. µ typically ranges from 1 (incoherent) to √ N (coherent).The smaller the µ, the less measurements one needs to reconstruct the data using minimum norm approaches.This is because the coherency parameter enters as the pre-factor in the lower-bound for the sensor quantity in l 1 -based CS for accurate recovery.There exist optimal sensor placement algorithms such as K-means clustering, the data-driven Online sparse Gaussian Processes [56], physics-based approaches based on inflection in the modal shapes [57] and mathematical approaches [13] that minimize condition number of Θ and maximize determinant of the Fisher information matrix [58].A thorough study on the role of sensor placement on reconstruction quality is much needed and is an active topic of research.In this study, we assess three different greedy approaches for nearly optimal sensor placement in sparse recovery applications, namely, the Discrete Empirical Interpolation Method (DEIM) [18,59], iterative condition number minimization [60] and the reduced matrix QR-factorization with column pivoting as reviewed in [61].We compare these approaches with outcomes from a random sensor placement that can be easily implemented using random number generators in Matlab.To ensure the resulting measurement matrices are incoherent with respect to the POD basis, we also estimate the coherency numbers in the following discussion.In this section, we summarize the different sensor placement algorithms employed in this research study.Optimal sensor placement for a given data, especially for a flow field that evolves over time is highly challenging and is an ongoing topic of active research.The goal of optimal point sensor placement for reconstruction is to identify and activate only a few rows of the basis matrix Φ that effectively conditions the matrix Θ (for P = K = N b ) or its variants M = Θ T Θ or M = ΘΘ T )depending on if P > K = N b or P < K = N b respectively).This is schematically illustrated in fig. 2. To design smart sensor placement, one needs an optimization criteria which in this case is to minimize reconstruction error when using a small number of sensors which, of course depends on the choice of basis Φ.Since reconstruction from sparse data in general requires inversion of Θ or M, most smart sensing strategies are designed to improve the condition of Θ, M for inversion purposes by optimizing their spectral content in the form of its determinant, trace, spectral radius or condition number.A direct method of optimizing such metrics requires searching over the different possible sensor selections resulting in combinatorial complexity.Thankfully, there exist a variety of greedy algorithms [13,18,62] that have been shown to be successful in the context of fluid flow data.

Random Sensor Placement
The most simple and efficient sensor placement strategy is to sample at random locations.This is commonly accomplished using a random number generator for the sensor locations, accomplished easily using packages such as Matlab or Python with functions such as randperm().In this study, we sample from simulation data by choosing the first P values from a random permutation of the entire data of dimension N .We note that it is equally effective to adopt ideas such as K-means clustering as was employed in [63].Since the experiments performed in this study use highly sparse measurements relative to the total number of grid points, i.e., P N to mimic practical data acquisition, we focus on random single-pixel measurements.To better assess the effectiveness of such sensor placement methods, we generate multiple realizations by changing the seed of the random number generator.The outcomes are then quantified in terms of the mean as well as the outliers.This particular choice of sensor placement is designed to serve as an inexpensive benchmark to compare against other greedy sampling methods.

Minimization of Matrix Condition Number (MCN)
As shown in sec.2.1.1,the success of the reconstruction effort for K = N b is tied to the accuracy of the inverse computation of M = Θ T Θ or M = ΘΘ T (depending on whether the system is over determined (P > K = N b ) or under determined (P < K = N b )).Therefore, if M = Θ T Θ or M = ΘΘ T has full column or row rank respectively along with a reasonable condition number, then the estimated inverse has a chance to be accurate.This approach focuses on sensor placement (through the construction of C) that minimizes the condition number of M or κ(M).The condition number is directly related to the orthogonality of Θ and the presence of significant diagonal entries in M .Therefore, this algorithm can be viewed as placing sensors at locations that display significant flow dynamics and also preserve the orthogonality of the POD modes.From a mathematical perspective, the condition number represents the ratio of maximum and minimum singular values of Θ or M. Therefore, for large κ(M), the errors tend to be amplified with respect to the signal.As shown by Willcox [13], and Yildrim et al. [62] such a method compares favorably to more physics-based approaches [57] where sensors are placed at the extrema of dominant POD modes although the latter method is computationally efficient.Such a heuristic approach has been used by Bayon de Noyer [64] to locate effective sensor locations for feedback control to alleviate tail buffeting of a high performance twin tail craft.However, the formulation based on condition number minimization of M (see eqn. ( 26)) is considered more reliable.The key steps of this MCN algorithm is listed below for completeness.
(i) Consider all possible placement points, evaluate M for each point, choose the point that minimizes condition number of M i.e. κ(M).
(ii) With the previous sensor locations set, loop over all possible remaining placement points.For each point, update the mask vector, evaluate M, and choose the point that minimizes κ(M).
(iii) In same way as (ii) find all remaining sensor locations.
A slightly more efficient version of this algorithm is presented by Willcox [60] where the first sensor point is chosen as the location that maximizes the sum of the difference between the diagonal and off-diagonal entries of M. The rest of the algorithm is similar to that presented above.

QR Factorization with Column Pivoting
The reduced matrix QR factorization [31] decomposes any given real matrix A ∈ R S×T with full column rank into a unitary matrix Q ∈ R S×T and an upper triangular matrix R ∈ R T ×T .Therefore, where r ii are the diagonal entries of R and λ i are the eigenvalues.It is easy to show that minimizing the condition number of A is related to optimizing the spectral characteristics of the matrix such as the determinant or spectral radius, i.e., maximize i r ii .In general, the R from a reduced matrix QR factorization has diagonal values, r ii in no particular sequence.However, when combined with column pivoting, we have AD = QR, where D ∈ R T ×T is a square column permutation matrix containing ones and zeros.
The resulting QR decomposition outcome can be controlled through the pivoting procedure such that the diagonal values of R, r ii form a decreasing sequence.Therefore, pivoting provides a smart approach for 'submatrix volume maximization' and in turn maximize the absolute value of the determinant [32] by reordering the columns of A. This approach can easily be extended for sensor placement by leveraging the connections between the permutation matrix D and the point sensor measurement matrix C in fig. 2 and eqn.(4).
For the case with P = K, the reconstruction problem in eqn.( 4) requires inversion of the square matrix Θ = CΦ k .For improved reconstruction, the determinant of Θ needs to be maximized through sensor placement which in turn is expected to reduce (and maybe minimize) the condition number.One can see that for a square matrix the following relationship is true.
Therefore, reduced matrix QR factorization of Φ T ∈ R K×N with column pivoting will yield where D ∈ R N ×N is a square permutation matrix.The right hand side of eqn.( 14) will be maximized for a given sensor quantity P if C is chosen as the first P rows of D T .The index locations of ones in each row of C are denoted by [ 1 , 2 , 3 , . . ., P ].
For the case of oversampled case with P > K, Θ is a tall and slender matrix whose inversion (section 2.1.1)is typically handled using a least squares minimization.The left (Moore-Penrose) pseudoinverse requires computing the inverse of the square matrix.M = Θ T Θ ∈ R K×K .Therefore, the sensor placement that increases the probability of accurate reconstruction should maximize det Θ T Θ so that condition number of M is bounded.Specifically, we will have the following: Leveraging the relationships in eqn.( 16), we see that maximizing the determinant of M = Θ T Θ can be realized by maximizing det Φ K Φ T K C T by using a reduced matrix QR factorization as shown below and once again choosing C as the first P rows of the N × N square matrix given by D T .The index locations of ones in each row of C are denoted by [ 1 , 2 , 3 , . . ., P ].
The algorithm of greedy sensor selection for an oversampled case using a given tailored basis Φ K and number of sensors P is summarized below:

Discrete Empirical Interpolation Method (DEIM)
All the above smart sensor placement methods focus on minimizing the condition number directly or indirectly controlling it through the determinant of the matrix to the inverted.The discrete empirical interpolation method or DEIM is a discrete variant of the Empirical Interpolation Method (EIM) originally presented by Barrault et al. [65] and subsequently extended to nonlinear model order reduction applications by Chaturantabut and Sorensen [18,59].Here, one recursively learns the interpolation points (with indices j ) at locations carrying maximum linear dependence error estimated at previously estimated interpolation points.The primary idea behind DEIM is to estimate a high-dimensional state using information at sparsely sampled interpolation points.Such techniques (other examples being Gappy POD [17] and missing point estimation or MPE [20]) are popular as hyper-reduction tools that bypass expensive nonlinear term computations in model order reduction.Naturally, one can adopt these interpolation points for sensor placement in sparse reconstruction applications.To illustrate this, the POD-based DEIM approximation of order M (number of interpolation points) for u(t) in the space spanned by the basis Ψ ∈ R N ×M is given by u = Ψb(t) where b(t) is coefficient vector, b(t) ∈ R M .When using POD bases, Ψ, one can simply estimate b(t) = Ψ T u(t), but this requires dealing with the higher dimensional state vectors ∈ R N that are computationally comparable to high-fidelity models even when using projection-based reduced order models.Hyper-reduction strategies [19] bypass this issue by estimating the approximate coefficients b(t) using carefully chosen set of M interpolation points instead of the full (N -dimensional) state so that computational cost scales with M .Specifically, one chooses M interpolation points corresponding to indices , where the interpolation or measurement matrix is given by D = [e 1 , ...., e M ] ∈ R N ×M with columns where Ψ(D T Ψ) −1 can be precomputed once to yield a N × M matrix while D T u(t) represents M -dimensional representation of the state denoting the interpolation points.This way, one avoids repeated computation of the high-dimensional u(t).One can easily see the connections between DEIM approximation and the sparse reconstructed state x SR = Φa with a obtained using eqns.( 5)-( 8) using C = D T .Therefore, estimating the interpolation points (D T ) is similar to estimating the sparse measurement locations in C. The indices 1 ....., M are estimated recursively from the input basis {Φ j } M j=1 using the algorithm 2 from [18].The process starts from selecting the first interpolation index 1 ∈ {1, 2, ..., N } using the first input POD basis |Φ 1 |.The remaining interpolation indices { j , j = 2, 3, ..., M } are selected such that they correspond to the largest magnitude of |r| where r = Φ j − Φa is the residual error between current input basis and its interpolation Φa obtained using {Φ 1 , Φ 2 , Φ 3 . . .Φ j−1 } at the indices { 1 , 2 , 3 . . .j−1 }.This step is given by line 5 of the algorithm 2. In lines 1 & 6, the 'max' function is similar to that available in MATLAB and |ρ| = |r j |.The residual represents a measure of the linear independence of Φ j with respect to the basis vectors ahead of this in the sequence and the interpolation point is at the maximum absolute value of this measure.Naturally, the realized j depends on the choice of basis Φ j and their sequence.The ordering of the input basis is not critical for the QR-pivoting based approach.The linear independence of the input basis ensures the above procedure is well-defined, i.e.D T Φ is non singular and ρ = 0 for all iterations.By using POD basis as the input, the linear independence and hierarchy (i.e.basis ordered in terms of decreasing singular values) characteristics are guaranteed which in turn ensures that the sparse interpolation indices are hierarchical and non repeating.

Sparse Recovery (SR) Framework
The reconstruction algorithm used in this work based on the Gappy POD or GPOD framework [17] and can be viewed as an l 2 minimization solution of the sparse recovery problem summarized through eqns.( 4),( 5) and ( 6) with Φ composed of K ≤ M basis vectors, i.e. dimension of a is K ≤ M .At this point, we remind the reader of naming conventions adopted in this paper: the instantaneous j th full flow state is denoted by x j ∈ R N , whereas the entire set consisting of M snapshots is assembled into a matrix form denoted by X ∈ R N ×M .This discussion focuses on single snapshot reconstruction as the extension to multiple snapshots is trivial, i.e. each snapshot can be reconstructed sequentially or in some cases be grouped together as a batch.This allows for such algorithms to be parallelized easily.
The primary difference between the SR framework in eqn.(4) as used in compressive sensing or DEIM-based approaches and GPOD [12,13,17,66] as shown in eqn.(21) is the construction of the measurement matrix C and the sparse measurement vector xj .In SR (eqn.( 4)), the down sampled state xj ∈ R P is a compressed version containing only the measured data, whereas in GPOD, xj ∈ R N is a masked version of the full state vector, i.e. values outside of the P measured locations are zeroed out to generate a filtered version of x j .For high resolution data x j ∈ R N with chosen basis Φ j ∈ R N , the low-dimensional features, a j ∈ R K are obtained from the relationship shown below: The masked (incomplete) data xj ∈ R N , corresponding measurement matrix C and mask vector m ∈ R N are related by: where C ∈ R N ×N .Therefore, the GPOD algorithm results in a larger measurement matrix with numerous rows full of zeros as shown in fig. 3 (compare with fig.1).To bypass the complexity of handling this N × N matrix, a mask vector, m ∈ N × 1 with 1s and 0s operates on x j through a point-wise multiplication operator < • >.To illustrate, the point-wise multiplication is represented as xj =< m j • x j > for each snapshot j = 1..M where each element of x j multiplies with the corresponding element of m j .This is applicable when each data snapshot, x j has its own measurement mask m j which is a useful way to represent the evolution of sparse sensor locations over time.The SR formulation in eq. ( 4) can also support time varying sensor placement, but would require a compression matrix, C i ∈ R P ×N that changes with each snapshot.The goal of the SR procedure is to recover the full data from the masked data given in eqn.( 22) by approximating the coefficients āj (in the l 2 sense) with basis, Φ K , learned offline using training data (snapshots of the full field data).
The coefficient vector āj cannot be computed by direct projection of xj onto Φ as these are not designed to optimally represent the sparse data.Instead, one needs to obtain the "best" approximation of the coefficient āj , by minimizing the error E j in the l 2 sense as show in eqn.(23).
In eqn.(23) we see that m j acts on each column of Φ through a point-wise multiplication operation which is equivalent to masking each basis vector Φ i .The above formulation is valid for a single snapshot reconstruction where the mask vector, m j could change with every snapshot xj for j = 1..M and the error E j represents the single snapshot reconstruction error to be minimized in order to estimate the approximate coefficients āj .It can easily be seen from below that one will have to minimize the different E j 's sequentially to learn the entire coefficient matrix, ā ∈ R K×M for all the M snapshots.Denoting the masked basis functions as Φi =< m j • Φ i >, eq. ( 23) is rewritten as in eqn.(24).
In the above, Φ is analogous to Θ = CΦ in eq. ( 4).To minimize E j , one computes the derivative with respect to āj and equated to zero as below: The result is the linear normal equation given by eqn.( 26) where M k1,k2 = Φk1 , Φk2 or M = ΦT Φ and f j i = xj , Φi or f j = ΦT xj .The reconstructed solution is given by eqn.( 27) below: Algorithm 3 summarizes the above SR framework assuming the basis functions (Φ i ) are known.The above solution procedure for sparse recovery is the same as that described in section 2.1 except for the dimensions of the x and C.

Algorithmic Complexity
The cost of computing the POD basis is O(N × M 2 ) where N is the full state dimension and M , the number of snapshots.The cost of sparse reconstruction turns out to be O(N × K × M ) for both the methods, where the K ≤ M is the system dimension chosen for reconstruction.Naturally, for many practical problems with an underlying low-dimensional structure, POD is expected to result in a smaller K than any other choice of basis.This helps reduce the sensor quantity requirement and reconstruction cost.Further, since the number of snapshots (M ) is also tied to the desired basis dimension (K), a larger K requires a larger M and in turn a higher cost to generate the basis.
The complexity of the sensor placement depends on the choice of method adopted.Using QR factorization with column pivoting, the cost of sensor placement for a N × N matrix is O(N 3 ) and O(N M 2 ) for a N × M matrix.The DEIM method carries a complexity of O(N M 3 ) when retaining M POD modes and identifying M sensors with a state dimension of N .The MCN algorithm requires an expensive combinatorial search to find the sensor location that minimizes the condition number of M ∈ R M ×M in each possible sweep.This results in a computational cost of O(N 2 M 3 ) for identifying M sensors.These estimates are consistent with our observation that the DEIM and QR-pivoting approaches yield results in very short time.

Low-dimensional Cylinder Wake Flows
Studies of cylinder wakes [24,[67][68][69] have attracted considerable interest from the model reduction and dynamical systems communities for its particularly rich flow physics content, encompassing many of the complexities of nonlinear flow systems, while easy to simulate accurately on the computer using established CFD tools.In this study, we explore data-driven sparse reconstruction for the cylinder wake flow at Reynolds numbers with unsteady wakes in the range of Re = 100 − 1000.To generate two-dimensional cylinder flow data, we adopt the spectral Galerkin method [70] to solve incompressible Naiver-Stokes equations, as shown in Eq. ( 28) below: where u and v are horizontal and vertical velocity components.P is the pressure field, and ν is the fluid viscosity.The rectangular domain used for this flow problem is −25D < x < 45D and −20D < y < 20D, where D is the diameter of the cylinder.For the purposes of this study, data from a reduced domain, i.e., −2D < x < 10D and −3D < y < 3D, is used.The mesh was designed to sufficiently resolve the thin shear layers near the surface of the cylinder and transit wake physics downstream.For the case of Re = 100 the grid includes 24, 000 points whereas for Re = 1000 the grid is refined to include approximately 95, 000 points for the sampled flow region.The computational method employed is a fourth order spectral expansion within each element in each direction.The sampling rate for each snapshot output is chosen as ∆t = 0.2 seconds.

Cylinder Wake Limit-cycle Dynamics
In this section, we explore sparse reconstruction of unsteady wake flows with well-developed periodic vortex shedding behavior.The GPOD type algorithm is chosen over the traditional Compressive sensing-based SR formulations to bypass the need for maintaining a separate measurement matrix, especially when employing point sensors to mimic realistic data acquisition.The first three POD modes and coefficients are computed and shown in figure 5 and 6, respectively.Similar to the velocity field visualized in figure 4, the dominant POD modes (mode 1 and mode 2) capture the symmetric vortex shedding patterns at various length scales for all the cases.The vestiges of onset of instability at higher Re numbers (Re = 1000) is observed from the asymmetry in mode 3.This is consistent with the observations in [71], where the laminar vortex shedding happens at around Re = 47, and becomes unstable at ∼ 190 which deforms the small-scale structures at higher POD modes.On the other hand, the temporal evolution of POD coefficients show periodic limit-cycle behavior for all the Re numbers with the characteristic frequencies increasing with Re as shown in fig.6.The dependence of sparse reconstruction performance on Re is beyond the scope this article although we expect the system dimensionality to increase with Re as the shear layers become thinner and consequently, require increased spatial sampling to capture all the relevant scales of motion.Instead, in this work we operate with narrower focus, i.e., explore and quantify how sparse reconstruction performance (in terms of error metrics) depend on sensor quantity, their placement and model dimensionality.With this in mind, we focus the analysis in the rest of the paper to a single Reynolds number, Re = 100.In this study, we choose 300 snapshots of data corresponding to a Reynolds number, Re = 100, corresponding to a non-dimensional time (T = U t D ) of T = 60 with uniform temporal spacing of dT = 0.2s.T = 60 corresponds to multiple (≈ 10 − 15) cycles of periodic vortex shedding behavior.

Smart Sensor Placement in Cylinder Wake
Figure 7 summarizes the realized sensor placement with a budget of 200 sensors using four different algorithms: (i) random algorithm in sec.2.4.1;(ii) QR factorization with column pivoting that maximizes the determinant capture for a given P as discussed in sec.2.4.3;(iii) discrete empirical interpolation method (DEIM) in sec.2.4.4 and (iv) explicit condition number minimization through an iterative combinatorial search (sec.2.4.2).All these sensor placement algorithms use POD bases computed using the 300 snapshots of data containing both u and v velocity components.As expected, the random sensor placement in fig.7(a) samples the entire domain without bias although the regions with higher grid density tend to be sampled more frequently.This includes the region close to the cylinder surface where the mesh resolution is high to capture the thin shear layers accurately.On the other hand, the three 'smart' sensor placement algorithms and shown fig.7(b)-7(d) identify locations either in the wake and also on the cylinder surface.The QR-Pivot and the DEIM methods identify a cluster of centers about two diameters downstream of the cylinder.Further downstream, the computed centers split into two narrow 'tail'-like regions on either side of the wake centerline.On the other hand, the MCN framework (fig.7(d)) identifies only five sensors in the wake, mostly downstream while a overwhelming number are predicted on the cylinder surface.This is potentially indicative of the low-dimensional characteristic of the wake flow where only a few POD modes are needed to capture a significant portion of the energy.

Sparse Reconstruction of Canonical Fluid Flows
In this section, we explore sparse reconstruction of cylinder wake flow at Re = 100 using the above SR infrastructure based on the GPOD formulation as against the traditional SR formulation.This choice is purely a matter of convenience and helps bypass the need for maintaining a separate measurement matrix (resulting in memory savings).In all the oversampled cases reported in this section, Tikhonov regularization is employed to reduce overfitting.

Sparse Reconstruction (SR) Experiments and Analysis
For this a priori assessment of SR performance we reconstruct sparse data from simulations where the full field representation is available.The sparse sensor locations are chosen as single point measurements using a variety of sensor placement methods as discussed in section 3.1.2and these locations are fixed for the ensemble of snapshots used for the reconstruction (i.e.we do not consider dynamic sensor placement).Reconstruction performance is evaluated by comparing the original simulated field with that from SR using POD as the choice of basis across the entire ensemble of numerical experiments.Of course, using such a data-driven basis requires availability of data so that one can compute the basis vectors a priori.In practice, one would have to build a library of basis vectors from data that can in turn be used for sparse recovery and is currently being developed within our team.In this we undertake this approach in order to focus on assessing the relative roles of system dimensionality (K), sensor sparsity (P ) and sensor placement (C) for the POD-based SR.In particular, we aim to accomplish the following through this study: (i) check if oversampling relative to desired system dimension (P > K) is sufficient condition for accurate reconstruction of the fluid flow for this POD-based l 2 reconstruction; (ii) understand how sensor placement impacts reconstruction quality for different choice of basis.
To learn the data-driven POD basis we employ the method of snapshots [52] as shown in eqns.( 9)- (12).For the numerical experiments described here, the data-driven basis and coefficients are obtained from the full data ensemble, i.e., M = 300 snapshots corresponding to T = 60 non-dimensional times.This gives rise to at most M basis for use in the reconstruction process in eqn.(3), i.e. a candidate basis dimension of N b = M .As shown in table 1, the choice of algorithms depend on the choice of system dimensionality (K), data sparsity (P ) and dimension of the candidate basis space, N b .Recalling from the earlier discussion in section 2, we see that P ≥ K would require an l 2 method for a desired reconstruction sparsity K as long as P ≥ N b .In case of POD-based SR, the basis are energy optimal for the training data and hence, contain built-in sparsity.That is, as long as the basis is relevant for the flow to be reconstructed, retaining only the most energetic modes (basis) should generate the best possible reconstruction for the given sensor placement and quantity.Therefore, the POD basis need to be generated once and the reconstructed system dimension is determined by just choosing to retain the first few modes in the sequence.For generic basis with no built-in ordering of when using a library of basis with no clear classification available a priori, one requires to search for the K most significant basis amongst the maximum possible dimension of N b = M using sparsity promoting l 1 methods.l 1 methods such as the iterative thresholding algorithms [27,72] require multiple solutions of the least-squares problem before a sparse solution is realized.This results in necessary, but increased computational cost.

Comparison of l 2 and l 1 Sparse Reconstruction using Energy-ordered POD Basis
In this subsection, we verify the equivalence between l 2 reconstruction using energy ordered POD basis and sparsity promoting l 1 − minimization reconstruction for the cylinder wake data when presented the same sparse sensor data.The success of this verification will imply that the chosen POD basis ordered in terms of energy of the training data is also relevant to the data being reconstructed and also validate the accuracy of our l 1 implementation.In such cases, choosing the K-most relevant basis for reconstruction is trivial (and efficient) as one chooses the first or last K modes (as against searching for the sparsest K modes using l 1 ) as the candidate basis for l 2 reconstruction provided sufficient sparse data, P > K, is available.To verify this, we consider two cases of reconstruction for the limit-cycle cylinder wake at Re = 100.Case (a): P > K, K = N b and using least squares reconstruction; Case (b): P > K, K < N b ≤ M and employs l 1 reconstruction using optimal matching pursuit.Figure 8 and table 2 compares the reconstructed fields and the K predicted weights for the two cases for a single snapshot corresponding to T = 0.2.In this study we use a greedy optimal matching pursuit algorithm for the l 1 norm minimization, CoSAMP as described by Needell and Tropp in [11].
The above results confirm the equivalence of these two approaches for POD-based reconstruction of cylinder wake flows considered in this study.It is worth noting that l 1 should be the preferred choice when the relevance of sparse data to the candidate basis space is unclear.For the rest of this paper, we will leverage energy hierarchy of the candidate POD basis and adopt the SR problem formulation using a l 2 minimization algorithm with K = N b .In other words, the candidate basis space is chosen based on the desired reconstructed field dimension and not fixed.

Sparsity and Energy Metrics
For this SR study, we explore the conditions for accurate recovery of information in terms of data sparsity (P ) and system sparsity (K) which also represents the dimensionality of the system in a given basis space.In other words, sparsity represents the size of a given basis space needed to capture a desired amount of energy.As long as the measurements are incoherent with respect to the basis Φ and the system is overdetermined, i.e., P > K, one should be able to invert Θ to recover the higher dimensional state, X.From earlier discussions in section 2, we know P > K is a sufficient condition for accurate reconstruction using l 0 minimization.Thus, both interpretations require a minimum quantity of sensor data for accurate reconstruction and is verified through numerical experiments in section 4.4.In this section, we describe how the different system sparsity metrics, K = N b , are chosen for the numerical experiments with both POD and ELM basis.Since the basis spaces are different, a good way to compare system dimensions is through the energy captured by the respective basis representations.For POD one easily defines a cumulative energy fraction captured by the K most energetic modes, E K , using the corresponding singular values of the data as shown in eqn.(29).
In the above, the singular values, λ are computed from eqn. (10), and M is the total number of possible eigenvalues for M snapshots.As a result, the energy fraction E K corresponding to sparsity K for the different flow Reynolds numbers, Re, is shown in Figure 9.For the cylinder flow case with Re = 100, one requires two and five POD modes to capture 95% and 99% of the energy content respectively, indicative of the low dimensional dynamics.We denote the mode numbers corresponding to 95% and 99% of the energy capture as K 95 and K 99 respectively.In this case, we compute err P OD K = X − Φ P OD 1..K a P OD 1..K 2 , where Φ P OD 1..K , and a P OD 1..K represent the matrix comprising of K POD bases and the corresponding coefficients for the different snapshots respectively.To assess SR performance across different flow regimes (that have different K 95 ) with different values of K we define a normalized system dimension, K * = K/K 95 and a normalized sensor sparsity, P * = P/K 95 .This allows us to design an ensemble of numerical experiments in the discretized P * − K * space and the outcomes can be generalized beyond the current study.In this work, the design space is populated over the range 1 < K * < 6 and 1 < P * < 12 for POD-based SR with K bounded by the total number of snapshots, M = 300.The lower bound of one is chosen such that the minimally accurate reconstruction captures 95% of the energy.If one desires a different reconstruction norm, then K 95 can be changed to K xx without loss of generality and the corresponding K-space modified accordingly.Alternately, one can choose E K , the normalized energy fraction metric to represent the desired energy capture as a fraction of E K 95 , but is not used in this study.To quantify the l 2 reconstruction performance, we define the mean squared error as shown in eqn.(30) below, where X is the true data, and XSR is the reconstructed field from sparse measurements as per algorithm 3.
N and M represent the state and snapshot dimensions and related to indices i and j, respectively.Similarly, the mean squared error F R K * 95 and F R K * for the full reconstruction from the POD basis can be computed as: where XFR is the full field reconstruction using exactly computed POD coefficients a.For POD this is simply a = Φ T X as per eqn.(12).K * 95 is the normalized system dimension (i.e.number of POD modes normalized by K 95 ) corresponding to 95% energy capture and K * = K/K 95 represents the desired reconstructed system dimension.Note that K * 95 = K 95 /K 95 = 1 is trivially seen to be one for this case.Using the above definitions, we can now generate normalized versions of the absolute ( 1 ) and relative ( 2 ) errors as shown in eqn.(33). 1 represents the SR error normalized by the corresponding full reconstruction error for 95% energy capture. 2 represents the normalized error relative to the desired reconstruction accuracy for the chosen system reconstruction dimension, K.These two error metrics are chosen so as to achieve the twin goals of assessing the overall quality of the SR in a normalized sense ( 1 ) and the best possible reconstruction accuracy for the chosen problem set-up (i.e P, K).Thus, if the best possible reconstruction for a given K is realized then 2 will take the same value across different K * .This error metric is used to assess relative dependence of P * on K * for the chosen flow field.On the other hand, 1 provides an absolute estimate of the reconstruction accuracy for a given flow system so that minimal values of P * , K * needed to achieve a desired accuracy can be identified.

Sparse Reconstruction of Cylinder Wakes
As one of the goals of this study is to establish the baseline SR quality using POD basis we performed a series of nearly one thousand SR experiments corresponding to different points in the P * − K * design space and spread over the different sensor placements.The underlying goal through this work is to bare the aspects of interplay between the various SR control variables, namely, reconstruction dimension, K, sensor quantity, P and the choice of sensor placement method (which affects C).In all these numerical experiments, the sensor placement is fixed across the range of snapshots as is the case in stationary sensing.

Sparse Reconstruction Accuracy
By computing the errors as described in section 4.3 across the K * − P * space, the contours of 1 and 2 are generated for both the random as well as greedy sensor placements in figures 10 and 11.As the random sensor placement can lead to high variability between realizations, we estimate multiple sets of sensor locations corresponding to different seeds (as denoted by β in this article).Specifically, we compute the SR errors from ten different random seeds and the corresponding reconstruction errors are presented in terms of the mean, maximum and minimum (based on the average over the K * − P * space) in fig.10.For the greedy 'smart' sensor placements a single realization is representative of the method and corresponding error contours for 1 and 2 are shown (fig.11).For ease of interpretation, the contour levels in both these figures are made consistent along with appropriately normalizing the error metrics.POD weights, a comparison with that for the full data reconstruction in fig.12.In particular, we show results for an over-sampled case P * = 10 for different reconstruction dimension, K * = 1, 2, 3 where all the different error contours in figs.10 and 11 show low error metrics.As expected, the quality of the isocontour reconstruction improves with increasing K * .Further, amongst the different experiments for a given over-sampled P * , K * , the DEIM and QR with column pivoting provide the most accurate sparse estimation of the POD coefficients (a), especially for the higher K * .The relative inaccuracy of the MCN framework in the estimation of a is observed even for these carefully chosen design points with low average error metrics although the discrepancies are barely noticeable.For such low-dimensional wake flows, these small errors in a amplify the discrepancy in the full field reconstruction.The relevant quantifications including the sensor quantity, placement method, reconstruction dimension and error metrics for these select dissection cases are summarized in table 3.In addition, we also estimate the coherency parameter for each of these cases.The computed µ values are O(1) indicating that the rows of the measurement matrix are sufficiently incoherent with respect to the POD basis on average.On careful, examination, we observe that the coherency parameters for the QR-pivoting and DEIM placements are higher than the random sensor placement.This again is not surprising as these approaches take advantage of the underlying physical structure as contained in the POD modes to determine sensor placement.

Sparse Reconstruction with Marginally Oversampled Sensors
We observed earlier that the importance of data-driven sensor placement will be felt most when the sensor quantity is comparable to the targeted reconstruction dimension, i.e.P * ≈ K * .To verify this, we dissect the instantaneous snapshot reconstruction corresponding to K * = 5 and P * = 6 as shown in fig.13.The plots in the left column of this figure denote the sensor locations, the middle column represents the instantaneous reconstructed fields (compared with the exact contours) and the right column shows the predicted and exact POD coefficients.As expected, the greedy smart sensor placements perform better at reconstruction when compared to random sensor placement which is evident from inaccuracy in the reconstructed coefficients (see figs.13(c,f,i,l) .This is because, the smart sensors are placed at locations that are physically relevant while the random placement as shown in fig.13 results in under-sampling of the wake region.While not all random sensor placement will generate bad reconstructions, we see strong variability in the error metrics between realizations as a consequence.This is clearly evident from fig. 10 where the variability from different choices of random sensor placement is captured primarily in the under-smapled and marginally oversampled regions of the P * − K * design space.

Sparse Reconstruction with Highly Oversampled Sensors
In the limit of oversampling, i.e.P * > K * , one expects all these different methods to perform accurately as it is not hard even using purely random placement to populate the physically relevant regions.However, we observe surprising inaccuracies with the MCN method in the limit of large P .We dissect this observation using a single design point, K * = 6, P * = 10 (as shown in fig.14) for which the sparse recovery should be close to exact when applied to this low-dimensional wake flow.Figure 14 shows the sensor locations for each of the different algorithms in the left column, instantaneous reconstructed fields (compared with the exact contours) in the middle column and the predicted POD coefficients compared with the exact values in the right column.We easily see that the sensor combinations corresponding to random, DEIM and QR with column pivoting produce identical results while the sensors based on condition number minimization generates erroneous outcomes.The reasons for this is not hard to decipher as the MCN estimates only six sensors in the wake of the cylinder and any additional sensors that are generated happen to be located on the cylinder surface.This in turn effectively under samples the system although P = 20, i.e. 20 sensors were used.This is clearly illustrated in fig.15 where the sensor locations identified for different budgets with P * = 4, 7, 10 are shown for QR with column pivoting, DEIM and MCN.These plots show that MCN generates all the sensors in the wake region for P * = 4 and subsequent increase in the budget does not produce meaningful additions.For this low-dimensional system, K = 5 modes captures 99% of the energy, the least condition number of M appears to be estimated when sensors are placed at locations where the field values are zero.More investigation is needed to identify the causes that underlie this observation.

Conclusion
In this work, we have developed a framework for systematic assessment of sparse reconstruction (SR) performance using flow normalized error metrics and basis/sensor dimension.Using this, we explore the interplay of sensor quantity, their placement and system dimension using optimal data-driven basis from proper orthogonal decomposition for reconstruction from sparse data of canonical low-dimensional wake flows.It is well known that sparse solutions are sought in such SR problems as the candidate basis is typically generic and exhaustive.While such generic candidate basis have a good probability of representing the data, the hierarchy of its relevance to the flow system of interest is not known a priori.In this study, we used a priori generated POD-basis along with knowledge of its energy hierarchy and relevance to the given (full state) flow system for sparse recovery.This allows us to bypass the need for more expensive l 1 norm minimization with l 2 approaches of complexity O(N KM ) as long as the same relevance exists for the sparse measurements.To verify this, we explicitly show that given basis relevance and ordering for the full flow state, l 1 methods using optimal matching pursuit and l 2 minimization (with the first K basis) using least squares yield the same K-sparse reconstruction from very few random measurements.In addition, we also observe from systematic analysis of error metrics over a carefully designed parameter (P * − K * ) space that even marginal oversampling, i.e.P K is sufficient for accurate full state recovery using l 2 SR with high probability which is akin to the requirement for l 0 minimization.We reiterate here that this study is not intended to advocate the use l 2 over sparsity promoting l 1 methods and acknowledge that the latter is the sensible choice in practical situations where knowledge of the underlying system is scarce.We further expand the P − K design space to include the effect of data-driven sensor placement with the following candidates: random sensing and greedy-smart sensing algorithms such as DEIM, QR with column pivoting and explicit condition number minimization or MCN.We observe that while the random sampling methods show highly variable errors for marginal oversampling, greedy-smart sensor placement show consistently accurate recovery under these conditions.In the limit of heavy oversampling, the MCN method produces diminishing returns as seen from the wake flow example reported in this study.This can be ascribed to the inability of MCN to generate meaningful sensors in dynamically relevant regions of the flow beyond a small threshold.More research is necessary to delineate the causes for this behavior.Considering the computational complexity of learning sensors from data and the propensity for reliable sparse reconstruction, QR factorization with column pivoting and DEIM (in that order) turn out to be the best alternatives to purely random sampling.

Figure 1 .
Figure 1.Schematic illustration of l 2 (left) and l 1 (right) minimization reconstruction for sparse recovery using a single-pixel measurement matrix.The numerical values in C are represented by colors: black (1), white (0).The other colors represent numbers that are neither 0 nor 1.In the above schematics x ∈ R P , C ∈ R P ×N , Φ ∈ R N ×N b and a ∈ R N b , where N b ≤ N .The number of colored cells in a represents the system sparsity K. K = N b for l 2 and K < N b for l 1 .

Figure 2 .
Figure 2. Schematic illustration of sparse sensor placement.The pastel colored rectangles represent rows activated by the sensors denoted in the measurement matrix through dark squares.

Figure 4 .Figure 5 . 3 Figure 6 .
Figure 4. Color contour snapshots of the instantaneous stream-wise velocity component for the cylinder flow at Re = 100 and Re = 1000.

( a )Figure 7 .
Figure 7. Sensor locations learned form different choices of methods considered in this work including both random as well as smart algorithms for the cylinder wake flow (Re = 100).These image show the most relevant 200 sensors, i.e. (P = 200).

Figure 8 .
Figure 8.Comparison of the sparse reconstruction using both l 1 and l 2 minimization methods for basis that is ordered in terms of energy content.The reconstructed and actual flowfields at T = 0.2 are compared in (a) for 2 and (b) l 1 .The corresponding POD features from both methods are shown in (c).

(a) 1 (Figure 10 .
Figure 10.Isocontours of the normalized mean squared POD-based sparse reconstruction errors (l 2 norm) corresponding to the sensor placement with maximum and minimum errors from the chosen ensemble of random sensor arrangements.The average error across the entire ensemble of ten random sensor placements is also shown.Left: normalized absolute error metric, 1 .Right: normalized relative error metric, 2 .

Figure 15 .
Figure 15.Incremental evolution of sensor placement using the different greedy sensor placement methods (QR-pivoting (left), DEIM (middle) and MCN (right)) as P * = P/K 95 is increased.We show plots for three different cases corresponding to P * = 4, 7, 10 in each of the rows.

Table 3 .
Sparse reconstruction performance quantification for different sensor location selection method for periodic cylinder flows at Re = 100. 1 is the SR error normalized by the exact reconstruction error corresponding to a dimension of K 95 . 2 is the SR error normalized by the exact reconstruction error corresponding to a dimension of K.