Data-Driven Sparse Sensor Placement for Reconstruction
Krithika Manohar Bingni W. Brunton J. Nathan Kutz Steven L. Brunton




















 
 

frameNum









seqn
counteqn

countfig

sidefig
Sidebar[2] 
	@edefcurrentlabelname#2
	@edefcurrentlabel#2
	
	
	
	
	
	
	
	
	
	[backgroundcolor=lightgray!20,shadow=true,roundcorner=8pt,frametitle=#1]
		 





frameNum1













customc
  belowcaptionskip=.25,
  breaklines=true,
  frame=L,
  xleftmargin=,
  language=Matlab,
  showstringspaces=false,
  basicstyle=,
  keywordstyle=,
  identifierstyle=,  
  commentstyle=,
  stringstyle=,
  backgroundcolor=
 
 
style=customc 



























































figures/






	
> 0




















[
\begin@twocolumnfalse

	Data-Driven Sparse Sensor Placement for Reconstruction
    Krithika Manohar^*, Bingni W. Brunton, J. Nathan Kutz, and Steven L. Brunton

	^*Corresponding author: kmanohar@uw.edu
    December 30, 2023
==========================================================================================================================



	
	















Optimal sensor placement is a central challenge in the design, prediction, estimation, and control of high-dimensional systems. High-dimensional states can often leverage a latent low-dimensional representation, and this inherent compressibility enables sparse sensing. This article explores optimized sensor placement for signal reconstruction based on a tailored library of features extracted from training data.  Sparse point sensors are discovered using the singular value decomposition and QR pivoting, which are two ubiquitous matrix computations that underpin modern linear dimensionality reduction. Sparse sensing in a tailored basis is contrasted with compressed sensing, a universal signal recovery method in which an unknown signal is reconstructed via a sparse representation in a universal basis.  Although compressed sensing can recover a wider class of signals, we demonstrate the benefits of exploiting known patterns in data with optimized sensing. In particular, drastic reductions in the required number of sensors and improved reconstruction are observed in examples ranging from facial images to fluid vorticity fields. Principled sensor placement may be critically enabling when sensors are costly and provides faster state estimation for low-latency, high-bandwidth control.
	
	
\end@twocolumnfalse]






  XXXXXXXX 
   Scalars  Notation 

  n  State dimension 

  m  Number of snapshots 

    p  Number of sensors (measurements) 
  
  r  Intrinsic rank of tailored basis _r
 
    K  Sparsity of state in universal basis 

    η  Variance of zero-mean sensor noise 

  Vectors  

  ∈^n  High-dimensional state 
  
  ∈^p  Measurements of state 

  ∈^r  Tailored basis coefficients 

  _j∈^n  Canonical basis vectors for ^n 

  ∈^n  K-sparse basis coefficients

  ∈ℕ^p  Sensor placement indices 

  ∈^n  POD modes (columns of _r) 
 
  ∈^1× r  Rows of  
 [5pt]
   Matrices  

  ∈^p× n  Measurement matrix 

    Unitary QR factor matrix 

    Upper triangular QR factor matrix

  ∈^n× n  Universal basis 

  _r∈^n× r  Tailored basis of rank r 

  =  Product of measurement and basis  

  ∈ℝ^n× m  Data matrix with m snapshots

 

















Optimal sensor and actuator placement is an important unsolved problem in control theory.  
Nearly every downstream control decision is affected by these sensor/actuator locations, but determining optimal locations amounts to an intractable brute force search among the combinatorial possibilities.
Indeed, there are n p=n!/(n-p)!p! possible choices of p point sensors out of an n dimensional state .  

Determining optimal sensor and actuator placement in general, even for linear feedback control, is an open challenge.  
Instead, sensor and actuator locations are routinely chosen according to heuristics and intuition.  
For moderate sized search spaces, the sensor placement problem has well-known model-based solutions using optimal experiment design <cit.>, information theoretic and Bayesian criteria <cit.> .


We explore how to design optimal sensor locations for signal reconstruction in a framework that scales to arbitrarily large problems, leveraging modern techniques in machine learning and sparse sampling. 

Reducing the number of sensors through principled selection may be critically enabling when sensors are costly, and may also enable faster state estimation for low latency, high bandwidth control. 

[]This article explores optimized sensor placement for signal reconstruction based on a tailored library of features extracted from training data. 
In this paradigm, optimized sparse sensors are computed using a powerful sampling scheme based on the matrix QR factorization and singular value decomposition. 
Both procedures are natively implemented in modern scientific computing software, and Matlab code supplements are provided for all examples in this paper <cit.>.
These data-driven computations are more efficient and easier to implement than the convex optimization methods used for sensor placement in classical design of experiments. 
In addition, data-driven sensing in a tailored basis is contrasted with compressed sensing, a universal signal recovery method in which an unknown signal is reconstructed using a sparse representation in a universal basis. 
Although compressed sensing can recover a wider class of signals, we demonstrate the benefits of exploiting known patterns in data with optimized sensing. 
In particular, drastic reductions in the required number of sensors and improved reconstruction are observed in examples ranging from facial images to fluid vorticity fields.
The overarching signal reconstruction problem is formulated in “Sidebar1".

This paper provides a tutorial overview of current sparse sampling methods for sensor placement and reconstruction of structured signals. 
We also connect and equate certain sampling strategies with analogues in the design of experiments literature.  

Near-optimal sensor locations are obtained using fast greedy procedures that scale well with large signal dimension. This work also generalizes and extends a powerful sampling scheme based on the matrix QR factorization and demonstrates its broad applicability to image and fluid flow reconstruction as well as polynomial interpolation.
The overarching sensor placement problem is summarized in “Sidebar1", and Matlab code is provided for all examples.

There are myriad complex systems that would benefit from principled, scaleable sensor and actuator placement, including fluid flow control <cit.>, power grid optimization <cit.>, epidemiological modeling and suppression <cit.>, bio-regulatory network monitoring and control <cit.>, and high-performance computing <cit.>, to name only a few.  
In applications where individual sensors are expensive, reducing the number of sensors through principled design may be critically enabling.  
In applications where fast decisions are required, such as in high performance computing or feedback control, computations may be accelerated by minimizing the number of sensors required.  
In other words, low-dimensional computations may be performed directly in the sensor space.

Scaleable optimization of sensor and actuator placement is a grand challenge problem, with tremendous potential impact and considerable mathematical depth.  

With existing mathematical machinery, optimal placement can only be determined in general using a brute-force combinatorial search.  
Although this approach has been successful in small-scale problems <cit.>, a combinatorial search does not scale well to larger problems.  Moore's law of exponentially increasing computer power cannot keep pace with this combinatorial growth in complexity.   

Despite the challenges of sensing and actuation in a high-dimensional, possibly nonlinear dynamical system, there are promising indicators that this problem may be tractable with modern techniques.  
High-dimensional systems, such as are found in fluids, epidemiology, neuroscience, and the power grid, typically exhibit dominant coherent structures that evolve on a low-dimensional attractor.  
Indeed, much of the success of modern machine learning rests on the ability to identify and take advantage of patterns and features in high-dimensional data.  
These low-dimensional patterns are often identified using dimensionality reduction techniques <cit.> such as the proper orthogonal decomposition (POD) <cit.>, which is a variant of principal component analysis (PCA), or more recently via dynamic mode decomposition (DMD) <cit.>, diffusion maps <cit.>, etc.   
In control theory, balanced truncation <cit.>, balanced proper orthogonal decomposition (BPOD) <cit.>, and the eigensystem realization algorithm (ERA) <cit.>, have been successfully applied to obtain control-oriented reduced-order models for many high-dimensional systems.  

In addition to advances in dimensionality reduction, key developments in optimization, compression, and the geometry of sparse vectors in high-dimensional spaces are providing powerful new techniques to obtain approximate solutions to NP-hard, combinatorially difficult problems in scaleable convex optimization architectures.  
For example, compressed sensing <cit.> provides convex algorithms to solve the combinatorial sparse signal reconstruction problem with high probability.  
Ideas from compressed sensing have been used to determine the optimal sensor locations for categorical decisions based on high-dimensional data <cit.>. 
Recently, compressed sensing, sparsity-promoting algorithms such as the lasso regression <cit.>, and machine learning have been increasingly applied to characterize and control dynamical systems <cit.>.  
These techniques have been effective in modeling high-dimensional fluid systems using POD <cit.> and DMD <cit.>.  
Information criteria <cit.> has also been leveraged for the sparse identification of nonlinear dynamics <cit.>, as in  <cit.>, and may also be useful for sensor placement. 

Thus, key advances in two fields are fundamentally changing our approach to the acquisition and analysis of data from complex dynamical systems:  1) machine learning, which exploits patterns in data for low-dimensional representations, and 2) sparse sampling, where a full signal can be reconstructed from a small subset of measurements.  The combination of machine learning and sparse sampling is synergistic, in that underlying low-rank representations facilitate sparse measurements.  
Exploiting coherent structures underlying a large state space allows us to estimate and control systems with few measurements and sparse actuation.  Low-dimensional, data-driven sensing and control is inspired in part by the high performance exhibited by biological organisms, such as an insect that performs robust, high-performance flight control in a turbulent fluid environment with minimal sensing and low-latency control <cit.>.  They provide proof-by-existence that it is possible to assimilate sparse measurements and perform low-dimensional computations to interact with coherent structures in a high-dimensional system (i.e., a turbulent fluid).  


Here we explore two competing perspectives on high-dimensional signal reconstruction:  1)  the use of compressed sensing based on random measurements in a universal encoding basis, and 2) the use of highly specialized sensors for reconstruction in a tailored basis, such as POD or DMD.  
These choices are also discussed in the context of feedback control.  
Many competing factors impact control design, and a chief consideration is the latency in making a control decision, with large latency imposing limitations on robust performance <cit.>.  
Thus, for systems with fast dynamics and complex coherent structures, it is important to make control decisions quickly based on efficient low-order models, with sensors and actuators placed strategically to gather information and exploit sensitivities in the dynamics.

 

 §.§ Extensions to dynamics, control, and multiscale physics

Data-driven sensor selection is generally used for instantaneous full-state reconstruction, despite the fact that many signals are generated by a dynamical system <cit.>.  
Even in reduced-order models, sensors are typically used to estimate nonlinear terms instantaneously without taking advantage of the underlying dynamics.  
However, it is well known that for linear control systems <cit.>, the high-dimensional state may be reconstructed with few sensors, if not a single sensor, by leveraging the time history in conjunction with a model of the dynamics, as exemplified by the Kalman filter <cit.>.  
In dynamic estimation and control, prior placement of sensors and actuators is generally assumed. 
Extending the sensor placement optimization to the model reduction <cit.> and system identification <cit.> of linear control systems is an important avenue of ongoing work.  
In particular, sensors and actuators may be chosen to increase the volume of the controllability and observability Gramians, related to the original balanced truncation literature <cit.>.  
More generally, sensor and actuator placement may be optimized for robustness <cit.>, or for network control and consensus problems <cit.>.

The sensor placement algorithms discussed above are rooted firmly in linear algebra, making them readily extensible to linear control systems.  
Recent advances in dynamical systems are providing techniques to embed nonlinear systems in a linear framework through a suitable choice of measurement functions of the state, opening up the possibility of optimized sensing for nonlinear systems.  
As early as the 1930s, Koopman demonstrated that a nonlinear system can be rewritten as an infinite-dimensional linear operator on the Hilbert space of measurement functions <cit.>.  
This perspective did not gain traction until modern computation and data collection capabilities enabled the analysis of large volumes of measurement data.  



Modern Koopman theory may drive sensor placement and the selection of nonlinear measurement functions on the sensors to embed nonlinear dynamics in a linear framework for optimal nonlinear estimation and control.  
This approach is consistent with neural control systems, where biological sensor networks (e.g., strain sensors on an insect wing) are processed through nonlinear neural filters before being used for feedback control.  
Much of the modern Koopman operator theory has been recently developed <cit.>, and it has been shown that under certain conditions DMD approximates the Koopman operator <cit.>; sensor fusion is also possible in the Koopman framework <cit.>.  
Recently, Koopman analysis has been used to develop nonlinear estimators <cit.> and controllers <cit.>, although establishing rigorous connections to control theory is an ongoing effort <cit.>.  
Koopman theory has also been used to analyze chaotic dynamical systems from time-series data <cit.>, relying on the Takens embedding <cit.>, which is  related to sensor selection.  








Beyond extending sensor selection to nonlinear systems and control, there is a significant opportunity to apply principled sensor selection to multiscale systems.  
Turbulence is an important high-dimensional system that exhibits multiscale phenomena <cit.>.  
Data-driven approaches have been used to characterize turbulent systems <cit.>, including clustering <cit.>, network theory <cit.>, DMD-based model reduction <cit.>, and local POD subspaces <cit.>, to name a few.  
Recently, a multiresolution DMD has been proposed <cit.>, where a low-dimensional subspace may locally characterize the attractor, despite a high-dimensional global attractor.  
This approach may significantly reduce the number of sensors needed for multiscale problems.   







§ COMPRESSED SENSING: RANDOM MEASUREMENTS IN A UNIVERSAL BASIS


The majority of natural signals, such as images and audio, are highly compressible, meaning that when the signal is written in an appropriate coordinate system, only a few basis modes are active.  
These few values corresponding to the large mode amplitudes must be stored for accurate reconstruction, providing a significant reduction compared to the original signal size.  
In other words, in the universal transform basis, the signal may be approximated by a sparse vector containing mostly zeros.  
This inherent sparsity of natural signals is central to the mathematical framework of compressed sensing.
Signal compression in the Fourier domain is illustrated on an image example in “Sidebar4". 
Further, sparse signal recovery using compressed sensing is demonstrated on a sinusoidal example in “sb:cs_example".




The theory of compressed sensing <cit.> inverts this compression paradigm.  
Instead of collecting high-dimensional measurements just to compress and discard most of the information, it may be possible to collect  a low-dimensional subsample or compression of the data and then infer the sparse vector of coefficients in the transformed coordinate system. 



 §.§ Theory of compressed sensing


Mathematically, a compressible signal 𝐱∈ℝ^n may be written as a sparse vector 𝐬∈ℝ^n in a new basis Ψ∈ℝ^n× n such that


    = .


The vector  is called K-sparse if there are exactly K nonzero elements.  
To be able to represent any natural signal, rather than just those from a tailored category, the basis  must be complete.  

Consider a set of measurements ∈ℝ^p, obtained via a measurement matrix ∈ℝ^p× n, which satisfies


    =  =  = .


In general, for p<n (<ref>) is underdetermined, and there are infinitely many solutions.  The least least squaressquares (minimum _2) solution is not sparse, and typically yields poor reconstruction.  
Instead, knowing that natural signals are sparse, we seek the sparsest  consistent with the measurements ,


    = _''_0 , such that  = ',


where _0 is the ℓ_0 pseudo-norm corresponding to the number of non-zero entries of .  
Unfortunately, this optimization problem is intractable, requiring a combinatorial brute-force search across all sparse vectors .  
A major innovation of compressed sensing is a set of conditions on the measurement matrix  that allow the nonconvex ℓ_0-minimization in (<ref>) to be relaxed to the convex ℓ_1-minimization <cit.>


    = _''_1 , such that  = ',


where _1 = ∑_k=1^n|s_k|.  This formulation is shown schematically in Fig. <ref>.
	
For the ℓ_1-minimization in (<ref>) to yield the sparsest solution in (<ref>) with high probability, the measurements  must be chosen so that = satisfies a restricted isometry property (RIP)


    (1-δ_K)_2^2 ≤_2^2≤ (1+δ_K)_2^2,


where δ_K is a small positive restricted isometry constant <cit.>.  
In particular, there are two conditions on  for a RIP to be satisfied for all K-sparse vectors :

	
  * The measurements  must be incoherent with respect to the basis .  
	This incoherence means that the rows of  are sufficiently uncorrelated with the columns of , as quantified by μ
	
    μ(,) = √(n)max_j,k|⟨𝐜_k,_j⟩|.

	Small μ indicates better incoherent measurements, with an optimal value of μ=1.  
	Here, 𝐜_k denotes the k-th row of  and _j the j-th column of , both of which are assumed to be normalized.  A more detailed discussion about incoherence and the RIP may be found in <cit.>.  
	
  * The number of measurements p must satisfy <cit.>
	
    p∼𝒪(Klog(n/K)).

	The Klog(n/K) term above is generally multiplied by a small constant multiple of the incoherence.  
	Thus, fewer measurements are required if they are less coherent. 


Intuitively, the existence of a RIP implies that the geometry of sparse vectors is preserved through the measurement matrix .  
Determining the exact constant δ_K may be extremely challenging in practice, and it tends to be more desirable to characterize the statistical properties of δ_K, as the measurement matrix  may be randomly chosen.  
“Sidebar_6" describes why it is not possible to use QR pivot locations as optimized sensors for compressed sensing, since they fail to identify the sparse structure of an unknown signal.  


Often, a generic basis such as Fourier or wavelets may be used to represent the signal sparsely.  
Spatially localized measurements (i.e., single pixels in the case of an image) are optimally incoherent with respect to the Fourier basis, so that μ(,)=1.  
Thus, single pixel measurements are ideal because they excite a broadband frequency response.  
In contrast, a measurement corresponding to a fixed Fourier mode would be uninformative; if the signal is not sparse in this particular frequency, this measurement provides no information about the other Fourier modes.  
For many engineering applications, spatially localized measurements are desirable, as they correspond to physically realizable sensors, such as buoys in the ocean.  

One of the major results of compressed sensing is that random projection measurements of the state (i.e., entries of  that are Bernoulli or Gaussian random variables) are incoherent with respect to nearly any generic basis  <cit.>.  
This result is truly remarkable; however, the incoherence of random projections is not optimal, and typically scales as μ∼√(2log(n)).  
Moreover, it may be difficult to obtain random projections of the full state  in physical applications.  

There are many alternative strategies to solve for the sparsest solution to (<ref>).  
Greedy algorithms are often used <cit.>, including the compressed sampling matching pursuit (CoSaMP) algorithm <cit.>.  
In addition, there is additional theory about how sparse the random projections may be for compressed sensing <cit.>.  




 §.§ Compressed sensing example


As a simple example, we consider a sparse signal that is constructed as the sum of three distinct cosine waves,


    x(t) = cos(2π× 37 t) +  cos(2π× 420 t) +  cos(2π× 711 t).


The Shannon-Nyquist sampling theorem <cit.> states that for full signal reconstruction, we must sample at twice the highest frequency present, indicating a theoretical minimum sampling rate of 1422Hz.  
However, since the signal is sparse, we may sample at considerably lower than the Nyquist rate, in this case at an average of 256Hz, shown in Figure <ref>.

Note that for accurate sparse signal reconstruction, these measurements must be randomly spaced in time, so that the relative spacing of consecutive points may be quite close or quite far apart.  
Spacing points evenly with a sampling rate of 256Hz would alias the signal, resulting in poor reconstruction. Matlab code for reproducing Figure <ref> is provided below.
[firstline=3,lastline=23]MATLAB/FIG_X_CS.m








§ OPTIMAL SPARSE SENSING IN A TAILORED BASIS





The compressed sensing strategy above is ideal for the recovery of a high-dimensional signal of unknown content using random measurements in a universal basis.  
However, if information is available about the type of signal (e.g., the signal is a turbulent velocity field or an image of a human face), it is possible to design optimized sensors that are tailored for the particular signals of interest.  

Dominant features are extracted from a training dataset consisting of representative exaples, for example using the proper orthogonal decomposition (POD).  
These low-rank features, mined from patterns in the data, facilitate the design of specialized sensors that are tailored to a specific problem.  

Low-rank embeddings, such as POD, have already been used in the ROM community to select measurements in the state space that are informative for feature space reconstruction.  
The so-called empirical interpolation methods seek the best interpolation points for a given basis of POD features.  

These methods have primarily been used to speed up the evaluation of nonlinear terms in high-dimensional, parameterized systems. However, the resulting interpolation points correspond to measurements in state space, and their use for data-driven sensor selection has largely been overlooked. We will focus on this formulation of sensor selection and explore sparse, convex, and greedy optimization methods for solving it. 

We begin with brief expositions on POD and our mathematical formulation of sensor placement, followed by an overview of related work in design of experiments and sparse sampling. We conclude this section with our generalized sensor selection method that connects empirical interpolation methods, such as QR pivoting to optimize condition number, with D-optimal experimental design <cit.>. 
The QR pivoting method described in “Sec:QR" is particularly favorable, as it is fast, simple to implement, and provides nearly optimal sensors tailored to a data-driven POD basis. []Finally, the distinctions between compressed sensing and our data-driven sensing are summarized in “sb:sensing_compare".



 §.§ Proper orthogonal decomposition


POD is a widespread data-driven dimensionality reduction technique <cit.> used in many domains; it is also commonly known as the Karhunen-Loève expansion, principal component analysis (PCA) <cit.>, and empirical orthogonal functions <cit.>. 
POD expresses high-dimensional states ∈^n as linear combinations of a small number of orthonormal eigenmodes  (i.e., POD modes) that define a low-dimensional embedding space. States are projected into this POD subspace, yielding a reduced representation that can be used to streamline tasks that would normally be expensive in the high-dimensional state space. 
This low-rank embedding does not come for free, but instead requires training data to tailor the POD basis to a specific problem.  
POD is illustrated on a simple example of extracting coherent features in images of human faces in “Sidebar2".

A low-dimensional representation of  in terms of POD coefficients 𝐚 can be lifted back to the full state with a linear combination of POD modes,

    _i ≈∑_k=1^r a_k(t_i)_k(x).

For time-series data _i, the coefficients a_k(t_i) vary in time and _k(x) are purely spatial modes without time dependence, resulting in a space-time separation of variables.  
Thus, care should be taken applying POD to data from a traveling wave problem.   

The eigenmodes _k and POD coefficients a_k are easily obtained from the singular value decomposition (SVD). Given a data matrix of state space observations =[_1  _2 … _m], the resulting eigenmodes are the orthonormal left singular vectors  of  obtained via the SVD,

    = ΨΣ 𝐕^T≈_r_r_r^T.

The matrices _r and _r contain the first r columns of  and  (left and right singular vectors, respectively), and the diagonal matrix _r contains the first r× r block of  (singular values).  
 The SVD is the optimal least squares approximation to the data for a given rank r, as demonstrated by the Eckart-Young theorem <cit.>

    _⋆ = _ - _F   s. t.  rank()=r,

where _⋆ = _r_r_r^T, and ·_F is the Frobenius norm. 
The low-dimensional vector of POD coefficients for a state  is given by the orthogonal projection = _r^T. 
Thus, the POD is a widely used dimensionality reduction technique for high-dimensional systems. This reduction allows computational speedup of numerical time-stepping, parameter estimation, and control. 

Choosing the intrinsic target rank without magnifying noise in the data is a difficult task. In practice, r is often chosen by thresholding the singular values to capture some percentage of the variance in the data. An optimal hard threshold is derived in <cit.> based on the singular value distribution and aspect ratio of the data matrix, assuming additive Gaussian white noise of unknown variance.  This threshold criterion has been effective in practice, even in cases where the noise is likely not Gaussian.






 §.§ Sensor placement for reconstruction


We optimize sensor placement specifically to reconstruct high-dimensional states from point measurements, given data-driven or tailored bases. Recall that full states may be expressed as an unknown linear combination of basis vectors

    x_j = ∑_k=1^r Ψ_jk a_k,

where Ψ_jk is the coordinate form of _r from (<ref>).  
Effective sensor placement results in a point measurement matrix  that is optimized to recover the modal mixture  from sensor outputs .  Point measurements require that the sampling matrix ∈^p× n be structured in the following way

    = [ __1 __2   … __p ]^T,

where _j are the canonical basis vectors for ^n with a unit entry at index j and zeros elsewhere.  Note that point measurements are fundamentally different than the suggested random projections of compressive sensing.  

The measurement matrix results in the linear system

    y_i = ∑_j=1^n C_ij x_j = ∑_j=1^n C_ij∑_k=1^r Ψ_jk a_k,

where C_ij is the coordinate form of  from (<ref>).
The observations in  consist of p elements selected from 

    =  = [x__1 x__2 …  x__p]^T,

where = {_1,…,_p}⊂{1,…,n} denotes the index set of sensor locations with cardinality ||=p. 

When  is unknown, it can be reconstructed by approximating the unknown basis coefficients  with the Moore-Penrose pseudoinverse, = ^†=(_r)^†.  Equivalently, the reconstruction is obtained using

    = _r,  = ^-1 =  (_r)^-1,     p=r, 
    ^† = (_r)^†,     p>r
    	.

A schematic of sparse sampling in a tailored basis _r is shown in Fig. <ref>.  
The optimal sensor locations are those that permit the best possible reconstruction . 
Thus, the sensor placement problem seeks rows of _r, corresponding to point sensor locations in state space, that optimally condition inversion of the matrix . 

For brevity in the following discussion we denote the matrix to be inverted by _=^T (_= if p=r). 
Recall that  determines the structure of , i.e. the sensor locations, and hence affects the condition numbers of  and _.  
The condition number of the system may be indirectly bounded by optimizing the spectral content of _ using its determinant, trace, or spectral radius. 
For example, the spectral radius criterion for _^-1 maximizes the smallest singular value of _

    _⋆ = _,||=p_^-1_2 = _,||=pσ_min(_) .

Likewise, the sum (trace) or product of magnitudes (determinant) of its eigenvalue or singular value spectrum may be optimized

    _⋆ = _,||=p(_,||=p) = _∑_iλ_i(_),


    _⋆ = _,||=p |_|    = _,||=p∏_i|λ_i(_)| 
       =_,||=p∏_i σ_i(_).


Direct optimization of the above criteria requires a combinatorial search over n p possible sensor configurations and is hence computationally intractable even for moderate n. Several heuristic greedy sampling methods have emerged for state reconstruction specifically with POD bases. These gappy POD <cit.> methods originally relied on random sub-sampling.  However, significant performance advances where demonstrated by using principled sampling strategies for reduced-order models (ROMs) in fluid dynamics <cit.>, ocean modeling <cit.> and aerodynamics <cit.>.  More recently, variants of the so-called empirical interpolation method (EIM, DEIM and Q-DEIM) <cit.> have provided near optimal sampling for interpolative reconstruction of nonlinear terms in ROMs. This work examines an approximate greedy solution given by the matrix QR factorization with column pivoting of _r^T, which builds upon the Q-DEIM method <cit.>.


		
[] 

	
 
The goal of sparse measurement selection is to choose a measurement matrix  representing p distinct point measurements in state space, so that  consists of p rows of the identity matrix.  
Although more general linear measurements of  may be admissible in some problems, point measurements are physically appealing as spatially localized sensors.  
The row-wise sum 𝐜^+=∑_j 𝐜_k satisfies 𝐜^+ _0 = p,
where 𝐜_k is the k-th row of . 
The sensor selection problem can be formulated as choosing  to make the pseudoinverse of =_r as well-conditioned as possible, thus making the estimation of the coefficients in  robust to measurement noise on .  The condition number is discussed in more detail in “Sidebar3".  

Condition number minimization is the viewpoint taken by the reduced-order modeling community for selecting point measurements within a low-rank POD basis. These empirical interpolation methods (EIM) <cit.> assume p=r sensors so  is square and invertible. EIM and subsequent discrete variants such as DEIM <cit.> and Q-DEIM <cit.> are greedy procedures that control the condition number of (_r)^-1 by minimizing the spectral norm (i.e. largest singular value)

    _⋆ = _(_r)^-1_2,

where  is subject to the same structural constraints as above. 
Related strategies include selecting measurements at maxima of successive POD modes or iteratively seeking measurements that decrease the condition number of _r <cit.>.

This strategy is particularly advantageous when only a few modes are required to characterize the data. When singular values decay slowly, we may require p>r measurements for well-conditioned reconstruction. In “Sec:QR", we generalize one particular empirical interpolation method, Q-DEIM, which uses QR pivoting to determine sensor locations.  In particular, we extend this method to the case of p> r sensors, making it viable for more general sensor selection and significantly speeding up its computation.



  §.§.§ General formulation

Sensor selection is more generally framed as follows: given a set of n possible measurement indices V, select a subset of sensor indices S to optimize a carefully chosen function evaluating their quality. The problem is mathematically realized as 

    max_S⊆ V f(S) S.

This problem belongs to an area of research called submodular function optimization, surveyed in <cit.>. The so-called submodular function, f(S):2^V↦, evaluates these subsets to some scalar value, thus rendering a brute force search over its entire domain (the power set of all possible subsets 2^V) computationally intractable even for moderate values of n. 

Possible choices of f include mutual information, entropy, estimation error covariance, and spatial coverage of sensors. Often, even a single evaluation of f is expensive when matrix determinants or factorizations are involved. The present work aims to reconstruct states from measurements with minimal variance of reconstruction error, which requires maximizing the determinant of error covariance matrices. The minimization of various error metrics is an active topic in optimal experiment design, for instance, the determinant or volume (D-optimal), the trace or mean-squared error (A-optimal), and maximum eigenvalue or worst-case variance (e-optimal) of error covariance matrices. 

Joshi and Boyd <cit.> frame sensor selection as selecting p rows of _r that minimize the volume of the estimation error covariance matrix, (_r)^T_r=^T. Since matrix volume is the absolute value of its determinant, the naive approach requires evaluating n p determinants over a combinatorially large set of row-selected submatrices. 
Instead, the authors approach this search with a computationally tractable convex optimization,

    _1 =    _ (_r)^T_r 
       ^+_0=p,   c_i^+ ∈ [0,1].

Each intermediate iteration of this method operates on an n× n matrix, so that storage requirements scale as 𝒪(n^2). Thus this approach is costly for high-dimensional states. The authors frame _r as a general matrix whose rows are possible measurements from which they select a subset of measurements. However, this matrix is not explicitly restricted to a POD basis, as proposed here.
 




 §.§ Sparse sensor placement with QR pivoting


An original contribution of this work is extending Q-DEIM to the case of oversampled sensor placement, where the number of sensors exceeds the number of modes used in reconstruction (p>r). The key computational idea enabling oversampling is the QR factorization with column pivoting applied to the POD basis. QR pivoting itself dates back to the 1960s by Businger and Golub to solve least squares problems <cit.>, and it has found utility in various measurement selection applications <cit.>. Similar to empirical interpolation methods such as DEIM, pivots from the QR factorization optimally condition the measurement or row selected POD basis, as described below.  

	
	The reduced matrix QR factorization with column pivoting decomposes a matrix ∈^m× n into a unitary matrix , an upper-triangular matrix  and a column permutation matrix  such that ^T =. The pivoting procedure provides an approximate greedy solution method for the optimization in (<ref>), which is also known as submatrix volume maximization because matrix volume is the absolute value of the determinant. QR column pivoting increments the volume of the submatrix constructed from the pivoted columns by selecting a new pivot column with maximal 2-norm, then subtracting from every other column its orthogonal projection onto the pivot column (see Algorithm <ref>). Pivoting expands the submatrix volume by enforcing a diagonal dominance structure <cit.>
	
    σ_i^2 = |r_ii|^2 ≥∑_j=i^k |r_jk|^2;    1≤ i ≤ k ≤ m.

	This works because matrix volume is also the product of diagonal entries r_ii
	
    | | = ∏_i σ_i = ∏_i |r_ii|.

	Furthermore, the oversampled case p>r may be solved using the pivoted QR factorization of _r_r^T, where the column pivots are selected from n candidate state space locations based on the observation that
	
    ^T  = ∏_i=1^r σ_i(^T),

	where we drop the absolute value since the determinant of ^T is nonnegative.


	
	




	
		QR factorization with column pivoting of ∈^n× m
		
	
	

Thus the QR factorization with column pivoting yields r point sensors (pivots) that best sample the r basis modes _r

    _r^T^T = .

Based on the same principle of pivoted QR, which controls the condition number by minimizing the matrix volume, the oversampled case is handled by the pivoted QR factorization of _r_r^T,

    (_r_r^T)^T = .

[]Algorithm <ref>, the QR factorization, is natively implemented in most scientific computing software packages. A Matlab code implementation of Algorithm <ref> is provided in “sb:qr_code".
The oversampled case requires an expensive QR factorization of an n× n matrix, whose storage requirements scale quadratically with state dimension. 
However, this operation may be advantageous for several reasons. 
The row selection given by the first p QR pivots increase the leading r singular values of ^T, hence increasing ^T.  
This is the same maximization objective used in D-optimal experiment design <cit.>, which is typically solved with Newton iterations using a convex relaxation of the subset selection objective. 
These methods require one matrix factorization of an n× n matrix per iteration, leading to a runtime cost per iteration of at least O(n^3). The entire procedure must be recomputed for each new choice of p. 
Our proposed method only requires a single O(n^3) QR factorization and results in a hierarchical list of all n total pivots, with the first p pivots optimized for reconstruction in _r for any p>r.  
Thus, additional sensors may be leveraged if available.  

The QR factorization is implemented and optimized in most standard scientific computing packages and libraries, including Matlab, LAPACK, NumPy, among many others. In addition to software-enabled acceleration, QR runtime can be significantly reduced by terminating the procedure after the first p pivots are obtained. The operation can be accelerated further using randomized techniques, for instance, by the random selection of the next pivot <cit.> or by using random projections to select blocks of pivots <cit.>.  
“Sidebar_6" shows why QR pivoting does not find the sparsest vector in an universal basis for compressed sensing.  


The sparse sensor placement problem we have posed here is related to machine learning
concepts of variable and feature selection <cit.>. 
Such sensor (feature) selection concepts generalize to data-driven classification. For image classification using linear discriminant analysis (LDA), sparse sensors may be selected that map via POD modes into the discriminating subspace <cit.>.  Moreover, sparse classification within libraries of POD modes  <cit.> can be improved by augmenting DEIM samples with a genetic algorithm <cit.> or adapting QR pivots for classification <cit.>. Sparse sensing has additionally been explored in signal processing for sampling and estimating signals over graphs <cit.>.






 §.§ Relation to optimal experimental design

The matrix volume objective described above is closely related to D-optimal experiment design <cit.>; in fact, the two problems are identical when regarding the tailored basis _r as a set of n candidate experiments of a low-dimensional subspace. Classical experimental design selects the best p out of n candidate experiments for estimating r unknown parameters ∈^r. Each experiment, denoted _i, produces one output y_i that may include zero-mean i.i.d. Gaussian noise ξ∼𝒩(0,η^2). Again, we wish to estimate the parameters from p experiment outputs ∈^p in the following linear system,

    = [ 0.4cm0.4pt_10.4cm0.4pt; 0.4cm0.4pt_20.4cm0.4pt;                      ⋮; 0.4cm0.4pt_p0.4cm0.4pt ]· + ξ = 
    ∑_k=1^r[ Ψ__1,k; Ψ__2,k;      ⋮; Ψ__p,k ]a_k + ξ
    
    
    
    
    
     =
    _r + ξ,

which is equivalent to the state reconstruction formulation of gappy POD <cit.>. In Matlab notation we sometimes refer to _r as _r(,:).
Each possible experiment _i may be regarded as a row of  or of the tailored basis _r such that _i = _r(_i,:).
Equivalently each _i is a weighted “measurement” of the lower dimensional POD parameter space (not to be confused with the point measurement operation ). Note that when all experiments are selected the output is simply the state vector  since = _r+ξ.

Given experiment selections indexed by , the estimation error covariance is given by 

    (-) = η^2(^T)^-1 =  η^2((_r)^T_r)^-1.

D-optimal subset selection minimizes the error covariance by maximizing the matrix volume of _=^T: 

    _⋆   =
    	_,||=plog∑_i=1^p _i^T_i
       = _,||=p (_r)^T_r,

which is equivalent to (<ref>). 
Similarly, A-optimal and E-optimal design criteria optimize the trace and spectral radius of ^T, and are equivalent to (<ref>) and (<ref>), respectively.  
The exact solutions of these optimization problems are intractable, and they are generally solved using heuristics. This is most commonly accomplished by solving the convex relaxation with a linear constraint on sensor weights β,

    β_⋆ = _β∈^n   log∑_i=1^n β_i_i^T_i,  
       ∑_i=1
     0≤β_i≤ 1^nβ_i = p.




The optimized sensors are obtained by selecting the largest sensor weights from β. 
The iterative methods employed to solve this problem,  convex optimization and semidefinite programs <cit.>, require matrix factorizations of n× n matrices in each iteration. Therefore they are computationally more expensive than the QR pivoting methods, which cost one matrix factorization in total. Greedy sampling methods such as EIM and QR are practical for sensor placement within a large number of candidate locations in fine spatial grids; hence, they are the methods of choice in reduced-order modeling <cit.>. The various optimization methods for data-driven sensor selection are summarized in Table <ref>.








§ COMPARISON OF METHODS

Sensor selection and signal reconstruction algorithms are implemented and compared on data from fluid dynamics, facial images, and ocean surface temperatures.  
The examples span a wide range of complexity, exhibit both rapid and slow singular value decay, and come from both static and dynamic systems. 

In each example, optimized sensors obtained in a tailored basis with QR pivots (See “Sec:QR") outperform random measurements in a universal basis using compressed sensing (See “Sec:CS") for signal reconstruction.  
Moreover, for the same reconstruction performance, many fewer QR sensors are required, decreasing the cost associated with purchasing, placing, and maintaining sensors, as well as reducing the latency in computations.  
Thus, for a well-scoped reconstruction task with sufficient training data, we advocate principled sensor selection rather than compressed sensing. 
For example, the QR-based sampling method is demonstrated with yet another tailored basis commonly encountered in scientific computing – the Vandermonde matrix of polynomials. “Sidebar_rev" compares polynomial interpolation with QR pivots for the ill-conditioned set of equispaced points on an interval.
When the structure of the underlying signal is unknown, then compressed sensing provides more flexibility with an associated increase in the number of sensors.  



 §.§ Flow past a cylinder


Fluid flow past a stationary cylinder is a canonical example in fluid dynamics that is high-dimensional yet reveals strongly periodic, low-rank phenomena. 
It is included here as an ideal system for reduction via POD and hence, minimal sensor placement. 
The data is generated by numerical simulation of the linearized Navier-Stokes equations using the immersed boundary projection method (IBPM) based on a fast multi-domain method <cit.>. 

The computational domain consists of four nested grids so that the finest grid covers a domain of 9× 4 cylinder diameters and the largest grid covers a domain of 72× 32.  Each grid has resolution 450× 200, and the simulation consists of 151 timesteps with δ t=0.02.  The Reynolds number is 100, and the flow is characterized by laminar periodic vortex shedding <cit.>. 


Vorticity field snapshots are shown in Fig. <ref>.  
In the cylinder flow and sea surface temperature examples, each snapshot _i=(t_i) is a spatial measurement of the system at a given time t_i. Thus POD coefficients a_k(t_i) are time dependent, and _k(x) are spatial eigenmodes. The first 100 cylinder flow snapshots are used to train POD modes and QR sensors, and reconstruction error bars are plotted over 51 remaining validation snapshots in Figures <ref> and <ref>.









The POD modes of this data reflect oscillatory dynamics characterized by periodic vortex-shedding.  
The data is low-rank, and the singular values decay rapidly, as shown in Fig. <ref>.  
The singular values occur in pairs corresponding to harmonics of the dominant vortex shedding frequency.  

Most of the spectral energy in the dataset is captured by the first 42 POD modes. Thus the intrinsic rank of the dataset in POD feature space is r=42, and the minimal number of QR pivots is p=42. This yields a dramatic reduction from the initial state dimension of n≈ 90000 spatial gridpoints. Here, QR pivoting of _r^T (with O(nr^2) operations) bypasses expensive O(n^3) factorizations of large n× n matrices with alternate sampling or convex optimization methods. 
	
Reconstruction from QR sensors (Fig. <ref>) successfully captures modal content with only p=r sensors when fitting to the first 42 POD modes. The first 42 POD modes characterize nearly 100% of the system's energy, the normalized sum of the singular values. 
Using modes beyond r>42 results in overfitting, and QR pivoting selects sensors based on uninformative modes.  
Thus, accuracy stops improving beyond r=42 target modes, whereupon sensor data amplifies noise. However, these tailored sensors perform significantly better than random sensors due to the favorable conditioning properties of QR interpolation points. 
“Sidebar_7" illustrates the ability of optimized sensing to significantly reduce the number of sensors required for a given performance. 
	


 §.§ Noise comparison study

	


	Measurements of real-world data are often corrupted by sensor noise.
	The POD-based sensor selection criteria, as well as A,D and E-optimal experimental design criteria, are optimal for estimation with measurements corrupted by zero-mean Gaussian white noise. 
	We empirically demonstrate this on the cylinder flow data with increasing additive white noise. 
	Here we assume sensor noise only in the test measurements and not in the training data or features, see Eqn. (<ref>).
	The POD modes and the different sensor sets are trained on the first 100 snapshots, and these different sensor sets are used to reconstruct the remaining 50 validation snapshots, which were not used for training features. 
	
	The reconstruction accuracy of the various sampling methods are compared for increasing sensor noise in Fig. <ref>,  alongside the full-state POD approximation for illustration. 
	Here we truncate the POD expansion to r=40 eigenmodes and compare the p=r reconstruction computed with the discrete empirical interpolation method (DEIM) <cit.> against the QR pivoting reconstruction (Q-DEIM, p=r). 
	The DEIM greedy strategy places sensors at extrema of the residual computed from approximating the k-th mode with the previous k-1 sensors and eigenmodes. 
	It can be seen that QR reconstruction is slightly more accurate than that of DEIM, which is the leading sampling method currently in use for reduced-order modeling <cit.>. 
	
	QR pivoting is competitive in both speed and accuracy. 
	The speed of QR pivoting is enabled by several implementation accelerations; for example, the column norms in line 4 of Algorithm 2 are only computed once and are then reused. 
	Unlike QR pivoting, DEIM and related methods add successive sensors per iteration by similarly optimizing some metric over all locations. 
	However, this metric (e.g., the approximation residual or the largest singular value) is recomputed at every iteration. 
	The QR factorization is significantly faster than convex optimization methods used in optimal design of experiments, which typically require one matrix factorization per iteration.

	In fact, convex optimization methods that relax the subset selection to weighted sensor placement provide no bounds for deviation from the global optimum, partly because rounding procedures are employed on the weights to decide the final selection. However, reconstruction error bounds for the globally optimal selection are known for DEIM <cit.>, Q-DEIM <cit.> and related POD sampling methods <cit.>. 
	Furthermore, QR pivoting can achieve significant accuracy gains over DEIM by oversampling – when p=2r QR reconstruction error is 4x smaller than that of DEIM. It should be noted that while DEIM has not yet been extended to the p>r case, oversampling methods exist for other POD-sampling methods <cit.>. 
	However, the iterative procedures involved in the latter are typically more expensive. 
	Recent accelerated variants of greedy principled sampling <cit.> may permit oversampling for large n, when oversampled QR storage requirements would be excessive. 
	In the cylinder flow case, we bypass this storage requirement by uniformly downsampling the fine grid by a factor of 5 in each spatial direction, thus reducing the number of candidate sensor locations to n=3600 instead of n=89351.
	


 §.§ Extended Yale B eigenfaces





Image processing and computer vision commonly involve high-resolution data with dimension determined by the number of pixels. Cameras and recording devices capture massive images with rapidly increasing pixel and temporal resolution. However, most pixel information in an image can be discarded for subject identification and automated decision-making tasks. 


The extended Yale B face database <cit.> is a canonical dataset used for facial recognition, and it is an ideal test bed for recovering low-rank structure from high-dimensional pixel space. 

The data consists of 64 aligned facial images each of 38 stationary individuals in different lighting conditions. 
We validate our sensor (pixel) selection strategy by recovering missing pixel data in a validation image using POD modes or eigenfaces trained on 32 randomly chosen images of each individual.

Normalized singular values are shown in Fig. <ref>, and the optimal singular value truncation threshold <cit.> occurs at r=166, indicating the intrinsic rank of the training dataset. 
Indeed, selected eigenfaces are also shown to reveal no meaningful facial structure beyond eigenface 166. 
QR pixel selection is performed on the first 50 and first 166 eigenfaces, and selected pixels shown in Fig. <ref> cluster around important facial features – eyes, nose and mouth.

















Image reconstructions in Fig. <ref> are estimated from the same number of selected pixels as the number of modes used for reconstruction. For instance, the 50 eigenface reconstruction is uniquely constructed from 50 selected pixels out of 1024 total – 5% of available pixels. Even at lower pixel selection rates, least squares reconstruction from QR selected pixels is more successful at filling in missing data and recovering the subject's face.


For comparison, reconstruction of the same face from random pixels using compressed sensing is shown in Fig. <ref>. Compressed sensing in a universal Fourier basis demonstrates progressively improved global feature recovery. However, more than triple the pixels are required for the same quality of reconstruction as in QR selection. Moreover, the convex ℓ_1 optimization procedure is extremely expensive compared to the single ℓ_2 regression on subsampled eigenfaces. Therefore data-driven feature selection and structured measurement selection are of significant computational and predictive benefit, and occur at the small training expense of one SVD and pivoted QR operation.



The convergence of reconstruction with sensors using QR pivoting is shown in Fig. <ref>. More sensors than modes are used in reconstruction for this example. The expected error dropoff is observed with increasing number of modes and sensors, although the dropoff is slower than for the cylinder flow (Fig. <ref>) due to slower decay of singular values.







 §.§ Sea surface temperature (SST)




Next we consider the NOAA_OISST_V2 global ocean surface temperature dataset spanning the duration 1990–2016. The data is publicly available online <cit.>.  

Unlike eigenfaces, this dataset is a time series, for which a snapshot is recorded every week. Sensor selection must then track energetic temporal signatures. 

Sensors and features are trained on the first 16 years (832 snapshots), and a test snapshot is selected from the excluded validation set. The singular values are shown in Fig. <ref>.

Like the eigenfaces, localized convective phenomena have energetic contributions to otherwise globally uninformative eigenssts. This is best seen in the POD snapshot projections, in which the 100 eigensst projection already sufficiently recovers dynamics, while increasing the number of eigenssts in the projection further refines convective phenomena. These lower-energy modes containing convective effects contribute to some degree of overfitting in ℓ_2 reconstruction (Fig. <ref>). The most interesting of these is the El Niño southern oscillation (ENSO) feature that is clearly identified from QR selected sensors. El Niño is defined as any temperature increase of a specified threshold lasting greater than six months in this highlighted region of the South Pacific. It has been implicated in global weather patterns and climate change. 




Remark: Modal separation of intermittent phenomena such as the El Niño is difficult from a time-invariant POD analysis. Separation of isolated, low-energy temporal events cannot be done from a variance-characterizing decomposition such as the POD – reordering the snapshots will yield the same dominant modes. On the other hand, tensor decompositions and temporal-frequency analyses such as multiresolution dynamic mode decomposition have succeeded at identifying El Niño where POD has failed. Sensor selection using non-normal modes arising from such decompositions remains an open problem and the focus of ongoing work.






§ DISCUSSION

The efficient sensing of complex systems is an important challenge across the physical, biological, social, and engineering sciences, with significant implications for nearly all downstream tasks.  
In this work, we have demonstrated the practical implementation of several sparse sensing algorithms on a number of relevant real-world examples.  
As discussed throughout, there is no all-purpose strategy for the sparse sensing of a high-dimensional system.  
Instead, the choice depends on key factors such as the amount of training data available, the scope and focus of the desired estimation task, cost constraints on the sensors themselves, and the required latency of computations on sensor data.  
Thus, we partition the sparse sensing algorithms into two fundamental categories:  1) optimized sensing in a data-driven tailored basis, and 2) random sensing in a universal basis.  


A critical comparison of the two approaches highlights a number of relative strengths and weaknesses.  
The first strategy results in a highly optimized set of sensors that are suitable for tightly scoped reconstruction problems where sufficient training data is available.  
The second strategy requires more sensors for accurate reconstruction but also makes fewer assumptions about the underlying signal, making it more general.  
We emphasize that optimized sensing in a tailored basis typically provides more accurate signal reconstruction than random measurements, facilitating a reduction in the number of sensors by about a factor of two.  
Further, sensor selection and signal reconstruction in the tailored basis is computationally efficient and simple to implement, while compressed sensing generally requires a costly iterative algorithm.  
In problems where sensors are expensive, or when low-latency decisions are required, the reduction in the number of sensors and the associated speed-up of optimized sensing can be significant.  
Thus, when the reconstruction task is well-scoped and a sufficient volume of training data is available, we advocate principled sensor selection rather than compressed sensing.  
In addition, pivoted QR sensors may be used in conjunction with other tailored bases (polynomials, radial basis functions) when signal structure is known. Since these are not data-driven basis functions, QR optimized samples can generalize to different dynamical regimes or flow geometries.





 §.§ Potential applied impact

Many fields in science and engineering rely on sensing and imaging.  
Moreover, any application involving feedback control for stabilization, performance enhancement, or disturbance rejection relies critically on the choice of sensors.  
We may roughly categorize these sensor-critical problems into two broad categories:  1) problems where sensors are expensive and few (ocean sampling, disease monitoring, espionage, etc.), and 2) problems where sensors are cheap and abundant (cameras, high-performance computation, etc.). 
 
In the first category, where sensors come at a high cost, the need for optimized sparse sensors is clear.  However, it is not always obvious how to collect the training data required to optimize these sensors.  In some applications, high-fidelity simulations may provide insight into coherent structures, whereas in other cases a large-scale survey may be required.  It has recently been shown that it may be possible to optimize sensors based on heavily subsampled data, as long as coherent structures are non-localized <cit.>.  

In the second category, where sensors are readily available, it may still be advantageous to identify key sensors for fast control decisions.  
For example, in mobile applications, such as vision-based control of a quad-rotor or underwater monitoring of an energy harvesting site using an autonomous underwater vehicle, computational and battery resources may be limited.  
Restricting high-dimensional measurements to a small subset of key pixels speeds up computation and reduces power consumption.  
Similar performance enhancements are already exploited in high-performance computing, where expensive function evaluations are avoided by sampling at key interpolation points <cit.>.  
Finally, it may also be the case that if measurements are corrupted by noise, reconstruction may improve if uninformative sensors are selectively ignored.  



§ REPRODUCIBLE RESEARCH

A Matlab code supplement is available <cit.> for reproducing results in this manuscript, including:

	
  * Datasets in Matlab file formats, or links to data that are publicly available online;
	
  * Matlab scripts to recreate figures of results.







ieeetr





§ ACKNOWLEDGMENTS

The authors thank Eurika Kaiser, Joshua Proctor, Serkan Gugercin, Bernd Noack, Tom Hogan, Joel Tropp, and Aleksandr Aravkin for valuable discussions.  
SLB and JNK acknowledge support from the Defense Advanced Research Projects Agency (DARPA HR0011-16-C-0016).  
SLB and KM acknowledge support from the Boeing Corporation (SSOW-BRT-W0714-0004).  
BWB and SLB acknowledge support from the Air Force Research Labs (FA8651-16-1-0003).
[]SLB acknowledges support from the Air Force Office of Scientific Research (AFOSR FA9550-16-1-0650). JNK acknowledges support from the Air Force Office of Scientific Research (AFOSR FA9550-15-1-0385).






































































































































































































































































































































































































  