PySensors: A Python Package for Sparse Sensor Placement

PySensors is a Python package for selecting and placing a sparse set of sensors for classification and reconstruction tasks. Specifically, PySensors implements algorithms for data-driven sparse sensor placement optimization for reconstruction (SSPOR) and sparse sensor placement optimization for classification (SSPOC). In this work we provide a brief description of the mathematical algorithms and theory for sparse sensor optimization, along with an overview and demonstration of the features implemented in PySensors (with code examples). We also include practical advice for user and a list of potential extensions to PySensors. Software is available at https://github.com/dynamicslab/pysensors.


Introduction
The success of predictive models and controllers in engineering and natural processes is largely determined by critical in situ measurements and feedback from sensors [3]. However, the deployment of sensors into complex environments, such as manufacturing [4], geophysical [5] and biological processes [6,7], is often expensive and challenging. Moreover, performance outcomes are extremely sensitive to the location and number of sensors deployed, motivating the optimal placement of sensors for diverse decision-making tasks. In general, choosing globally optimal placements within the search space of a large-scale complex system is an intractable computation, in which the number of possible placements grows combinatorially with the number of candidates [8]. While sensor placement, or sensor selection, has traditionally been guided by expert knowledge and first principles models, the continued growth in system complexity has motivated new mathematical paradigms and optimization algorithms for data collection and data-driven modeling that aim to automate the optimal, or near-optimal, sensor selection task.
PySensors is a Python package for the scalable optimization of sensor placements from data. In particular, PySensors provides tools for sparse sensor placement optimization approaches that employ data-driven dimensionality reduction [1,2]. This approach results in near-optimal placements for various decision-making tasks and can be readily customized using different optimization algorithms and objective functions. The PySensors package is aimed at researchers and practitioners alike, enabling anyone with access to measurement data to engage in scientific model discovery. The package is designed to be accessible to inexperienced users, adhering to scikit-learn standards, while also including customizable options for more advanced users. A number of popular sensor placement variants are implemented, but PySensors is also designed to enable further extensions for research and experimentation.
Maximizing the impact of sensor placement algorithms requires tools to make them accessible to scientists and engineers across various domains and at various levels of mathematical expertise and sophistication. PySensors unifies the algorithms developed in the recent papers [1,2,30] and their accompanying codes SSPOR pub and SSPOC pub into one software package. The only other packages in this domain of which we are aware are Chama [31] and Polire [32]. While these packages and PySensors all enable sparse sensor placement optimization, Chama and Polire are geared towards event detection and Gaussian processes respectively, whereas PySensors is aimed at signal reconstruction and classification tasks. As such, there are marked differences in the objective functions optimized by PySensors and its precursors. In addition to these two packages, researchers and practitioners have made available various custom scripts for sensor placement. Currently, researchers seeking to employ modern sensor placement methods must choose between implementing them from scratch or manually augmenting existing unpolished codes.

Background
PySensors was designed to solve reconstruction and classification tasks, which often arise in the modeling, prediction, and control of complex processes in geophysics, fluid dynamics, biology, and manufacturing.

Reconstruction
PySensors implements the sparse sensor placement optimization for reconstruction (SSPOR) method for recovering high-dimensional signals x from linear sensor measurements of the form Given data in the form of state measurements x k ∈ R n , k = 1, . . . , m, SSPOR identifies the optimal measurements of x, given by the operator C, which describes which components of x to observe. The SSPOR framework aims to find the best subset of the available measurements (usually, components of x) from which the full signal can be recovered in the estimation problem We assume that C is a mostly sparse subset selection operator consisting of rows of the identity, with nonzero entries designating the selected measurements. Given a p sensor budget and n candidates state components, C is constrained to have the following structure where e j is the canonical basis vector with a unit entry at the jth component and zeros elsewhere. The action of this measurement operator extracts the selected components of the signal

Classification
The SSPOC formulation seeks the placement of a small number of point sensors that classify highdimensional signals x ∈ R n as one of c classes. We start with labeled training data {(x i , y i )|i = 1, 2, . . . , m}, where training examples x i ∈ R n are concatenated into a matrix X ∈ R n×m and labels y i ∈ {0, 1, . . . , c − 1} are collected in a vector y ∈ R m .
An r-dimensional feature basis of the training data Ψ ∈ R n×r is extracted, where r m. Next, a decision space that best separates classes of training data w is computed on training data projected onto Ψ. w is obtained as the weights of a linear classifier fit to predict labels y from features Ψ † X. SSPOC seeks a measurement vector s that satisfies Ψ † s = w, where Ψ † denotes the Moore-Penrose pseudoinverse of Ψ.
In particular, we seek the sparse solution s s = argmin where is a small error tolerance. The non-zero elements of s are the sensor locations.

Features
PySensors enables the sparse placement of sensors for two classes of problems: reconstruction and classification. Each problem has a dedicated class with a number of options for tailoring the package to different applications of interest. Additionally, PySensors provides methods to enable straightforward exploration of the impacts of critical hyperparameters like the number of sensors or basis modes. Furthermore, because PySensors was built with scikit-learn compatibility in mind, it is easy to use cross-validation to select among possible choices of bases, basis modes, and other hyperparameters. The package is divided into three primary submodules: reconstruction, classification, and basis. The reconstruction and classification submodules contain classes and methods for reconstruction and classification problems, respectively. And basis houses different implementations of various bases. There are also two helper submodules: optimizers holds optimization routines used by the sensor selectors and utils has a variety of utility functions. In the subsections that follow we cover each of the main submodules in more depth.

Reconstruction
For reconstruction problems the package implements a unified SSPOR class (SSPOR is an acronym for Sparse Sensor Placement Optimization for Reconstruction), with methods for efficiently analyzing the effects that data or sensor quantity have on reconstruction performance [1]. When a SSPOR instance is fit to measurement data, internally it first fits a basis object to the data (see Section 3.3), then employs the computationally efficient and flexible QR algorithm [33][34][35], which has recently been used for hyperreduction in reduced-order modeling [29] and for sparse sensor selection [1]. The learned sensors are specific to the basis that was chosen; the method leverages structure present in the basis to decide which sensors are most important.
Often different sensor locations impose variable costs, e.g. if measuring sea-surface temperature, it may be more expensive to place buoys/sensors in the middle of the ocean than close to shore. These costs can be taken into account during sensor selection via a built-in cost-sensitive optimization routine [30]. This CCQR algorithm is found in the optimizers submodule along with the standard QR algorithm.

Classification
For classification tasks, the package implements the Sparse Sensor Placement Optimization for Classification (SSPOC) algorithm [2], allowing one to optimize sensor placement for classification accuracy. The algorithm is related to the compressed sensing optimization [36][37][38], but identifies the sparsest set of sensors that reconstructs a discriminating plane in a feature subspace. To instantiate a SSPOC object, one specifies a basis and a linear classifier. The implementation is fully general in the sense that it can be used in conjunction with any linear classifier and any basis from the basis submodule, however a linear discriminant analysis (LDA) classifier and Identity basis are used by default. When the fit method is called, a SSPOC instance (a) fits a basis object to the data, (b) fits the classifier to a set of examples dependent on the newly learned basis, (c) solves an optimization problem involving the weights of the classifier and the basis, (d) the sensors are selected based on the output of (c), and (e) the classifier is optionally refit using data sampled at the chosen sensor locations. The learned sensors depend on the combination of basis and classifier and on additional hyperparameters.
SSPOC employs different optimization methods depending on whether a binary or multi-class classification problem is being solved. The Scikit-learn orthogonal matching pursuit implementation is used to solve binary problems and the multi-task Lasso implementation is used for multi-class problems. Note that the CVX 1 package was used in the original SSPOC formulation [2].

Basis
It is well known [1] that the basis in which one represents measurement data can have a pronounced effect on the sensors that are selected and the quality of the reconstruction. Users can readily switch between different bases typically employed for sparse sensor selection: • Identity: Use the raw measurement data directly, without modification. This class also empowers the user to work with other bases than those provided by PySensors-one can map the data to a different basis before feeding it to a PySensors class instantiated with an Identity basis.
• SVD: Use the left singular vectors from a truncated singular value decomposition. Only the specified number of modes are computed to minimize the computational footprint. A randomized SVD can also be computed to cut costs even further.
• RandomProjection: Multiply measurements with random Gaussian vectors to project them to a new space. This basis is related to compressed sensing approaches [36][37][38].
Each of these classes is housed in the basis submodule.

Reconstruction examples
Consider the problem of interpolating a real-valued function f on [0, 1]. Suppose we wish to use polynomials to perform this interpolation. The sparse sensor placement problem is then equivalent to the problem of selecting points in [0, 1] at which to sample f . Before we can choose these optimal interpolation points, we must choose a polynomial basis to use. Ideally we should work with a numerically stable basis such as Chebyshev polynomials, but suppose we choose something simple such as monomials: import numpy as np x = np . linspace ( 0 , 1 , 1 0 0 1 ) r = 11 psi_r = np . vander ( x , r , increasing=True ) . T In keeping with Scikit-learn conventions, each row corresponds to an example and each column to a feature or sensor location. This choice departs from the mathematical literature wherein examples are typically represented as column vectors. We can then use the SSPOR class to find close-to-optimal sensor locations tailored to this basis. The fit method is used to learn the locations.
from pysensors . reconstruction import SSPOR selector = SSPOR ( ) selector . fit ( phi_r ) The SSPOR object now contains a list of sensor locations ranked in descending order of importance stored in its ranked sensors attribute. Note that because the model was fit on 11 examples, only the first 11 sensor locations are meaningful. The remaining indices are given in random order. We can specify the number of sensors we would like to see via the set n sensors function: which prints the following locations These approximate the Fekete (or Gauss-Lobatto) points, which are known to be optimal for interpolation via monomials [39]. Alternatively we could have initialized the object with this preference SSPOR(n sensors=10).
The fitted SSPOR object can also be used to reconstruct the signal (perform interpolation) from sparse measurements. The predict function is used to accomplish this. We will apply the method to a function not in the range of the monomial basis (see Figure S6 of [1]): For comparison, in Figure 1 we plot the interpolants generated by both our method and equispaced points.
from numpy . linalg import lstsq A common question to ask at this point is "How does the reconstruction error depend on the number of sensors I use?" We can use the reconstruction error method to get an answer. We simply specify an array of values of n sensors to try and the test data on which to measure the error and the function will compute a metric of interest for each number of sensors. The default metric or scoring function is the root-mean-square error. sensor_range = np . arange ( 2 , r + 1 ) recon_error = selector . reconstruction_error ( f , sensor_range ) A plot of the reconstruction error is given in Figure 2.

Classification examples
Next we turn to the problem of identifying a sparse set of sensor positions optimized for classification tasks. For the sake of simplicity, we will work with the digits dataset of Scikit-learn, which consists of eight-by-eight images of handwritten digits. See Figure 3a for example images. Our overall objective is to train a classifier to predict which digit is drawn in each image. We will add the restriction that the classifier is limited in the pixels it is allowed to see. The class designed to select the most salient sensor locations is named after the algorithm it employs: Sparse Sensor Placement Optimization for Classification, or SSPOC for short. SSPOC instances can be used in much the same way as standard Scikit-learn estimators. They are fit to the data with a fit method:   Note that once the estimator has been fit it expects subsequent samples to consist only of measurements taken at the sensors it has chosen. This is because it refits its internal classifier-in this case linear discriminant analysis (LDA)-on the subsampled data upon deciding on a set of sensor locations. Figure 3b visualizes the 10 pixels that were picked. Which pixels optimize classification accuracy depends on the classifier that is used. The PySensors SSPOC implementation is compatible with any linear classifier.
The number of sensors can be modified after fitting, but the training data are needed to refit the classifier for the new set of sensors: classifier . update_sensors ( n_sensors=5 , xy=( X_train , y_train ) ) There are many other parameters affecting the performance of the SSPOC class that we omit from this discussion. Please see the documentation and examples for a more comprehensive exploration of such options.

Basis examples
Both the SSPOR and SSPOC constructors accept a basis parameter, which specifies the basis in which to represent the data: It is often useful to track how the performance of the model changes as the number of basis modes is varied. One option would be to refit the model for each number of basis modes. However, for bases such as SVD, which are computed based on input data, this approach is wasteful. Each time the fit method is called, the SVD modes will be recomputed. For this reason PySensors has convenience functions for efficiently updating the number of basis modes. They allow us to avoid refitting the basis, but not re-computing the sensor locations. The methods are slightly different for SSPOR and SSPOC objects in that the training data are required for SSPOC but not for SSPOR: It is important to start with a large number of basis modes, then update to smaller numbers of basis modes.

Reconstruction
It is important to select the appropriate basis for a given problem, such as choosing polynomials for function approximation in the example above. For a general large, data-driven problem such as sea surface temperatures or photographs of faces, the identity basis (performing QR on the raw snapshots) will produce the lowest reconstruction error at a given number of sensors. This is because no information is lost to construct a low-rank approximation of the data, but by that same virtue, the identity basis can lead to impractically long run times for a large data set. Hence the built-in options to use SVD or randomized projections, both of which provide optimal or nearoptimal low-rank approximations. Other basis options that could be manually employed include the dynamic mode decomposition basis, Fourier modes, and basis modes arising from the solution of a system's equations of motion, if known.
In general, the reconstruction error will decrease as the number of sensors and basis modes increases, but the number of sensors relative to the number of modes is also important and depends on the basis. With an SVD basis, as the number of modes is increased, any additive noise in the measurements begins to dominate the reconstruction error, leading to the unintuitive result that reconstruction error increases with the number of modes, also see [40]. This can be mitigated by oversampling, i.e. using more sensors than modes. Note that when oversampling, PySensors randomly selects the sensors beyond the number of modes r, which has been shown to be fast and effective [40,41].
Conversely, with random projections, the quality of the reduced-order approximation depends on the presumed rank of the system [42,43]. In order for the randomized approximation to have a high probability of accurately representing the system, the number of modes should be at least five or ten more than the system's presumed rank. When sparsely sampling, the rank of the system can be at most equal to the number of sensors p, and so we recommend choosing at least p + 10 basis modes.
If using a cost function, determine the trade-off between cost savings and reconstruction accuracy by multiplying the cost function by a constant factor. If this factor is set to zero, unmodified QR is performed and the sensor locations with the lowest reconstruction error will be returned. If the weighting is large (what constitutes "large" depends on the system, basis, and cost function), low-cost sensors will be selected, regardless of their effectiveness for reconstruction.

Classification
Just as for reconstruction, the choice of an appropriate low-dimensional basis for the data of interest is crucial. Generally, the SVD or a random projection produce low-rank approximations that are nearly optimal for most datasets. The choice of r, the rank of the projection, determines the number of sensors chosen. For binary classification, the number of sensors chosen is approximate r; for c > 2, the number of sensors is at most r(c − 1).
For certain datasets, some reweighing of the bases improves the performance of the sparse classification. When the most discriminating directions in the data are not well captured by the most energetic modes, it is helpful to construct a biased basis by reweighing the basis vectors by largest magnitude elements in w [7]. In other words, instead of using the largest singular values Σ r to choose the Ψ r basis, we may start with a larger r, compute w, and then use the largest elements of Σ r |w| to learn a modified SVD basis better tailored to the training data.
Since the learning of sparse sensor locations and the execution of classification on those sparse sensors are separate, decoupled tasks, we typically re-train a separate classifier on the sparse sensor locations. While the theoretical derivation of SSPOC relies on linear projections and discriminates, this ultimate step can take many forms, including any nonlinear classifier.