Support Vector Spectrum Approximations: Efficient Calculations of Light Scattering Lineshapes
Gregory J. Hunt Cody R. Ground Robin L. Hunt


[


Support Vector Spectrum Approximations: Efficient Calculations of Light Scattering Lineshapes
    Gregory J. Hunt^1,*, 
  Cody R. Ground^2, 
  and Robin L. Hunt^2

  ^1William & Mary, PO Box 8763 Williamsburg, VA 23187

^2Hypersonic Airbreathing Propulsion Branch, NASA Langley Research Center, 12 Langley Blvd. Hampton, VA 23681.

^*ghunt@wm.edu

    December 30, 2023
=============================================================================================================================================================================================================================================================













Measuring scattered light is central to many laser-based gas diagnostic techniques, e.g., coherent anti-Stokes Raman spectroscopy (CARS) and filtered Rayleigh scattering (FRS). To produce quantitative measurements with such techniques, a computational model of the scattered spectral lineshape is necessary. While accurate, these models are often quite computationally demanding and thus cannot be used in situations where computational speed matters. To overcome this, approximations of these spectral lineshape models can be used instead. In this paper, we develop a method called support vector spectrum approximation (SVSA). This method uses machine learning to create efficient and accurate approximations of any existing spectral lineshape model. The SVSA framework improves upon existing methods by allowing efficient approximations of spectral lineshapes to be calculated in arbitrary flow regimes. We demonstrate the efficacy of SVSA in approximating coherent and spontaneous Rayleigh-Brillioun spectra. We also show that SVSA reduces the computational cost of a simulated filtered Rayleigh scattering experiment by a factor of 300.
]



§ INTRODUCTION


Measuring scattered light is central to many laser-based gas diagnostic techniques. For example, coherent anti-Stokes Raman spectroscopy (CARS) and filtered Rayleigh scattering (FRS) are techniques using Raman and Rayleigh-Brillouin  (RB) scattering, respectively, to quantify properties of a gas such as pressure, velocity, temperature, and number density <cit.>. Such laser-based gas diagnostic techniques are useful in a diverse set of fields. In atmospheric science, measurements of RB scattered light are used for remote sensing of clouds, aerosols, winds, ice crystals, etc. <cit.>. In the physics and engineering fields, Raman and RB scattering are used to nonintrusively investigate flames and other complex flows <cit.>.

To produce quantitative measurements using these techniques, one often needs to know the spectral lineshape of the scattered light. This lineshape conveys important information about the gas' properties. Consequently, accurate computational models linking gas properties to the spectral lineshape have been developed in the literature. For example, CARSFT can be used to calculate coherent anti-Stokes Raman spectra and Tenti's S6 model or Pan's S7 model allow calculation of RB spectra <cit.>. While typically quite accurate, such lineshape models are often computationally demanding. Thus, they are not ideal for situations where computational speed matters, for example, in optimization studies or real-time diagnostic techniques <cit.>.

For RB scattering, this problem has motivated methods to efficiently approximate the RB spectral lineshape. In particular, approximating Tenti's S6 model has been previously investigated <cit.>. Such methods have had good success in the regime on which they are designed. However, there does not presently exist a general approach for rapidly approximating general spectral lineshapes for any given scattering process over an arbitrary regime. Existing approximation frameworks are specific to particular types of scattering, like Rayleigh-Brillouin, and cannot be used to produce approximations of, for example, Raman spectra. Similarly, methods may only be valid for a small range of parameters or they may simplify the lineshape model by treating parameters as constant even though they are not. Consequently, while these methods may be accurate for, say, atmospheric regimes, they may not transfer well to highly dynamic conditions like those in supersonic flows.

This paper proposes a method called support vector spectrum approximation (SVSA), that creates efficient and accurate approximations of any spectral lineshape given a preexisting computational model. The method applies machine learning to produce a general lineshape approximator. The framework is applicable to preexisting lineshape models of any scattering type in an arbitrary flow regime. We show that SVSA is both fast and accurate by demonstrating its efficacy approximating coherent and spontaneous RB spectra. Furthermore, in application on a real problem, we demonstrate how SVSA might be used in an FRS study of hypersonic vehicle engines to make the computations feasible. 



§ METHODS


Given any preexisting lineshape model, we generically denote the set of input parameters for the model as φ∈ℝ^P where P denotes the number of parameters. For example, P=5 for the S6 and S7 models of RB scattering, which have as parameters y, Euken factor, internal relaxation number, internal specific heat, and translational specific heat. Given M values of the (normalized) frequency x_1,…,x_M, a lineshape model r takes φ and produces an estimate r(φ) ∈ℝ^M of the spectrum at x_1,…,x_M. Figure <ref> illustrates several examples of spontaneous Rayleigh-Brillouin scattering spectra estimated by Tenti's S6 model.



 The goal of SVSA is to create fast approximations of r(φ). SVSA does this by leveraging machine learning. Specifically, the method combines support vector regression (SVR) and the singular value decomposition (SVD) to build a flexible and general approximator.



 §.§ The SVSA Model


At the heart of SVSA is an assumption that most of the variability among spectral lineshapes can be explained by a small number of dominant lineshape modes denoted f_1,…,f_K. That is, for any , we assume that 

    ≈∑_k=1^Kw_k()f_k

so that a lineshape is approximately the linear combination of some dominant global modes f_1,…,f_K with parameter-specific weights w_k(). Using this approximation, SVSA estimates  by estimating the {f_k} and the {w_k()} and combining them according to Equation (<ref>).

Given training data generated by a preexisting computational lineshape model, SVSA does this estimation through a two-step approach:

    
  * Decompose the training data with the SVD to estimate K functional modes f̂_1,…,f̂_K∈ℝ^M. 

  * Learn the relationship between  and  using SVR. This produces, for each k, a model ŵ_k where ≈. 


Combining the estimates from these two steps, SVSA predicts  as

    = ∑_k=1^Kf̂_k.

An important feature of SVSA is that its calculation is fast. We show in the next section that  can be succinctly written as 

    = β_0 + β z()

with

    z()    = exp(-||Σ-||^2)

where β_0∈ℝ^M and β∈ℝ^M× S are the model coefficients from the SVR models expressed in the {f̂_k} basis and Σ∈ℝ^S× P is a matrix with rows {Σ_s}_s=1^S that are “support vectors” Σ_s∈ℝ^P chosen by SVR. The notation “Σ-” in Equation (<ref>) denotes the S × 1 vector where

    (Σ-)_s = Σ_s-.

Notice that we can expand -||Σ-||^2 as

    -(Σ-)'(Σ-) = 2Σ - diag^-1(Σ'Σ) - '

where the diag^-1 operator forms the vector of diagonal elements of a matrix. Thus, calculating  boils down to matrix multiplication and exponentiation. This means SVSA can be calculated very efficiently using standard scientific computing software. 



 §.§ Fitting SVSA


The SVSA method learns its model parameters β_0,β and Σ using a set of N training examples. Let φ_1,…,φ_N∈ℝ^P be a collection of training parameters that cover the range of the input space over which we want to learn the approximator and let r(φ_1),…,r(φ_N)∈ℝ^M be the corresponding training lineshapes calculated by our existing lineshape model. The first task when fitting SVSA is to learn the dominant functional modes f_1,…,f_K. We do this with the SVD. Let R∈ℝ^M × N be the lineshape matrix containing the N training spectra so that R_mn is the value of  r(_n) at x_m. If R has an SVD of R=UDV', then define W=DV' so that R=UW. This factorization allows us to decompose the n^th column of R as 

    r(_n) = ∑_k=1^M W_knU_k

where U_k is the k^th column of U. These columns of U (the noncentered principal components of the training lineshapes) capture the dominant functional modes of the lineshapes in our training data. Figure <ref> displays the top four functional modes of some example training data. The first mode approximately captures the average lineshape of the training data and subsequent modes capture typical deviations from this average.



Since the {U_k} are ordered decreasing according to how much variation they explain, if we truncate the linear combination in Equation (<ref>) after K<M terms, then ∑_k=1^K W_knU_k is the best approximation of r(_n) using only K modes. We thus let f̂_k=U_k for k=1,…,K and choose K so that these {f̂_k} capture most of the variation present in the training lineshapes. We can often achieve a very good approximation using K ≪ M. We explore the accuracy of this approximation in Section <ref>.

The second step to training SVSA is to learn the relationship between  and w_k(). To do this, we need training examples of the  and corresponding weights w_k() from which to learn. We can reuse, the training data used to estimate the {f_k} by noting that for each training parameter _n the correpsonding weight w_k(_n) is precisely W_kn. To see this, consider the similar form of Equations (<ref>), (<ref>) and (<ref>). The {W_kn}_k=1^K in Equation <ref> are essentially the weights {w_k(_n)}_k=1^K used to approximately reconstruct r(_n) from f_1,…,f_K (Equation (<ref>)). Thus, to learn the w_k, we use the training parameters _n and corresponding “training weights” W_kn=w_k(_n) to train K different SVR models. That is, for each k=1,…,K, we find a fit to predict W_kn from _n. We call this fit ŵ_k so that ŵ_k(_n)≈ w_k(_n)=W_kn.

The SVR approach used for this fit finds a flexible and regression-like function to predict w_k() from . Fitting these models amounts to using the training data to identify an important subset of training points Σ_1,…,Σ_S∈ℝ^P, called support vectors, and an associated set of coefficients α_k0,…,α_kS∈ℝ. The SVR models use these parameters to predict w_k() from  as a linear combination of the {α} weighted by the similarity between  and Σ_s. Specifically,

    ŵ_k()=α_k0+∑_s=1^Sα_ksH(Σ_s,)

where H is a kernel function measuring the similarity between Σ_s and . Essentially, the SVR models use the support vectors as a set of important landmarks in the parameter space of all  and estimates w_k()≈ w_k(Σ_s) to the extent that ≈Σ_s.

Using SVR-based estimators over more traditional approaches like multiple regression has substantial advantages. First, SVSA uses ν-SVR, a version of SVR parameterized by ν∈ [0,1]. This type of SVR finds an optimal fit for the data that uses (approximately) ν percentage of the training points as support vectors. Thus, ν directly controls the model complexity and allows for an intuitive trade-off between accuracy and speed. A low value of ν forces the fit to only retain the few most important support vectors. In this case, predictions will be less computationally demanding by sacrificing accuracy. Conversely, a high value of ν means the model can keep track of many support vectors and build a nuanced and highly accurate model. However, this prediction accuracy comes at the cost of computational speed.

A second advantage to SVR is our use of a radial basis function (RBF) kernel H in Equation (<ref>) defined as H(x,y) = exp(-||x-y||^2) for x,y∈ℝ^P. This allows SVR to fit highly nonlinear relationships between  and w_k(). While multiple regression would find a fit that is linear in , RBF kernelized SVR amounts to fitting a regression-like estimator in an expanded infinite dimensional space of all possible polynomials of . Thus, SVSA's approach allows an increased amount of flexibility in the model fit. It automatically considers all possible polynomial relationships between  and w_k() and finds the best fit among those while avoiding overfitting.  


Finally, once fit, SVSA is very fast to compute. If ŵ()=(ŵ_1(),…,ŵ_k()) then we can succinctly write ŵ() = α_0 + α z() with z as in Equation (<ref>), where α_0∈ℝ^K and α∈ℝ^K× S contain the coefficient estimates across all K of our SVR models. The parameter Σ in z (see Equation <ref>) is the combined matrix of all support vectors showing up in any of our models. Here Σ∈ℝ^S× P where S is the total number of unique support vectors across the models. Letting F̂=[f̂_1|⋯|f̂_K]∈ℝ^M× K, we have

    = F̂ = F̂(α_0 + α z()) = β_0 +β z()
  where we define β_0 = F̂α_0 and β = F̂α.
Thus, the major components of calculating  are (1) calculating Z(), which mainly involves the exponential and the product Σ, and (2) calculating a product between β and Z(). More generally, if we have a matrix Φ=[_1 | _2 | ⋯ | _L ]∈ℝ^P× L consisting of parameter values _1,…,_L∈ℝ^P for L lineshapes, then we can efficiently approximate all of the lines as R̂(Φ)∈ℝ^M × L letting

    R̂(Φ) = β_0 + β Z(Φ)

with

    Z(Φ)   =exp(-||Σ-Φ||^2).

Here, Σ-Φ denotes an S × L matrix where

    (Σ-Φ)_sℓ = Σ_s - _ℓ,

the the norm ||·|| is applied columnwise and the exponential elementwise. As Z is a Gaussian kernel, it must be bounded for all Φ and cannot be ill-conditioned.  Again, notice that we can write the exponent as

    -||Σ-Φ||^2 = 2ΣΦ - diag^-1(Σ'Σ) - diag^-1(ΦΦ')
 where the first subtraction is rowwise and the second is columnwise.  

All together, after estimating β_0, β and Σ using the training data, SVSA boils down to the following three steps codified in Algorithm <ref>:



The computationally expensive parts of each step are (1) calculating E in 𝒪(SPL) operations, (2) exponentiating the elements of E using SL exponentiations, and (3) calculating β Z(Φ) with 𝒪(MSL) operations. As each of these steps are fast to compute using modern scientific computing packages, SVSA allows us to very efficiently compute a highly flexible, nonlinear approximation of  with basically just two matrix products and exponentiation. Code for fitting and applying our method can be found on github at <gjhunt.github.io/svsa>. Currently, the fitting and lineshape approximation algorithms are implemented in Python. However, once the coefficients have been estimated in Python, they may also be exported and used to approximate the lineshape using routines we have written in Matlab and R. 



§ EVALUATION ON DATA


To assess the accuracy and speed of SVSA, we generate example lineshape data on which to evaluate our method. We generate spontaneous Rayleigh-Brillouin scattering (RBS) spectral lineshapes according to Tenti's S6 model <cit.>. Tenti's model estimates RBS using the parameters: y, Euken factor, internal relaxation number, internal specific heat and translational specific heat. Our primary interest in SVSA is application to a supersonic wind tunnel where a shock wave system is generated to simulate the flow conditions in a high speed air-breathing engine. Thus, we generated our training lineshapes using a range of parameters typically seen in these wind tunnel studies. We set the internal specific heat (c_int) to 1, the translational specific heat (c_tr) to 3/2, and vary P=3 parameters: the Euken factor (Eukenf) over a grid of 5 points from 1.8 to 2, the internal relaxation number (rlx_int) over a grid of 5 points from 1.5 to 3, and the y factor (y) over a grid of 20 points from 0 to 8. In total, this generates a set of 500 training lineshapes. Each lineshape is calculated on a grid of 500 values of the nondimensionalized frequency (x) ranging from -3 to 3.

In Figure <ref>, we plot the cumulative percentage of variance captured in the training data by the first several estimated modes f̂_k. 

We can see from this figure that almost all of the variability among spectral lineshapes can be captured using the first several dominant modes. By approximating the spectral lineshapes using K≈ 6 modes, we are able to recover >99% of the important differences among lineshapes. The dimensionality reduction achieved by these modes makes SVSA efficient to compute since we need only train 6 SVR models.

To further explore the relationship between model complexity and accuracy, we generate a collection of test lineshapes covering the same parameter ranges as the training data but over a finer grid. For the test data, we generate 32,000 example lineshapes. In Figure <ref>, we plot the maximum error of SVSA (across all 32,000 test lineshapes) against the average time needed to calculate a lineshape using the method. We do this for three values of ν, ten values of K, and two ways of calculating lineshapes.



 
 The y-axis in Figure <ref> measures the maximum difference between the SVSA approximated lineshape and the true lineshape as a percentage of the maximum height of the true line. The x-axis in Figure <ref> measures the time it takes SVSA to make the predictions. We display two sets of curves for two ways of calculating the 32,000 lineshapes. First, the lineshapes are calculated by predicting all 32,000 lineshapes simultaneously. This is done by feeding SVSA a matrix Φ∈ℝ^3 × 32000 of all parameters as per Algorithm <ref>. We divide the total time it takes to make all predictions by 32,000. Secondly, we calculate the test lineshapes in a serial fashion feeding one parameter vector ∈ℝ^3× 1 to SVSA at a time. We then loop over all 32,000 parameters and report the average computation time. These two ways of calculating the lineshapes cover two common use-cases for SVSA. The first is the batch setting where all of the parameters are known in advance of the calculation. For example, this occurs when analyzing data collected from an experiment or a computational fluid dynamics (CFD) simulation. The second situation covers doing calculations in a streaming context where new data arrives over time. In this case, we do not know the parameter values in advance and thus, can't calculate them in bulk. We must approximate each lineshape separately as the data arrives. This is what our loop over the 32,000 test cases demonstrates. An example of this situation is real-time laser diagnostics. For both approaches, the displayed time in Figure <ref> is the average number of (wall-time) seconds it takes our Python implementation to make a prediction on a 4.2 GHz processor.

Figure <ref> shows a large difference between estimating the lineshapes simultaneously and estimating them serially. It is almost an order of magnitude faster to simultaneously calculate the lineshapes. This happens because approximating the lineshapes simultaneously allows the entire calculation to be handled extremely efficiently by low-level linear algebra packages. Alternatively, the overhead when calculating the lineshapes serially slows down the average approximation speed. Nonetheless, in either case, SVSA is still extremely fast. Tenti's model takes about one-tenth of a second to calculate a single lineshape on our 4GHz processor. Thus, depending on the parameter setttings, SVSA is between one-thousand and one-hundred-thousand times faster.

We also notice in Figure <ref> that both K and ν have substantial impacts on both speed and accuracy. Increasing either K or ν increase the model complexity, and thus decrease error at the cost of speed. For ν, this effect is straight-forward. Increasing ν increases the number of support vectors. This leads to more nuanced and accurate models that are slower to compute. For K, larger values indirectly affect model complexity by increasing the complexity of the SVR models we build. This is because the weight functions w_k, corresponding to a mode f_k, are increasingly complex as we increase k since higher-order f_k pick up more subtle and complicated modes. Modeling these complicated modes requires more flexible models and hence more support vectors. Thus, increasing K correspondingly increases model complexity which decreases error at the expense of speed. Notice, however, increasing K quickly reaches diminishing returns. Indeed using more than about 6 modes seems to only increase computation time. This is because the higher order modes explain an ever-smaller proportion of the variability among the lineshapes. For example, from Figure <ref>, we see that the 7^th mode accounts for less than 0.003% of each line's shape. Thus, including these higher order terms increases model complexity for relatively little gain in approximation accuracy. For our application, it seems that we do well using K=6 and ν=0.1. This yields a maximum error in our test set of about 4% and an average computation time between about 1× 10^-5 and 8× 10^-5 seconds per line.

To test our method further, we build SVSA models to approximate coherent and spontaneous RB spectra using Tenti's S6 model and Pan's S7 model. We build these for models by training them on 500 example lineshapes over the same flow parameters previously used and setting K=6 and ν=0.1. We then evaluate these approximators on 32,000 test lineshapes generated over the same parameter space. In Figure <ref>, we plot several example predicted lineshapes for our four models. The black lines are the true spectral lineshapes and the red lines are the estimated lineshapes by SVSA. In all cases, we can see that SVSA produces highly accurate approximations of the lineshapes. 

In Figure <ref>, we plot boxplots of the percentage error for each of the four applications of SVSA across the 32,000 test examples. This figure shows that we have a median error of about 10^-2 percent for the S6 and S7 models for both spontaneous and coherent scattering. These errors are well below the accuracy of the S6 and S7 models as compared to the true scattering lineshape. This demonstrates that SVSA is a versatile estimator that is applicable across models, like S6 and S7, and across scattering types, like coherent and spontaneous RBS. SVSA can approximate these different models with relatively little tuning and produce low error predictions across a large input parameter space. 




§ APPLICATION


To evaluate the efficacy of SVSA in a real application, we apply SVSA to a computational fluid dynamics (CFD) study of filtered Rayleigh scattering (FRS). FRS is an  experimental laser-diagnostic technique whereby the total energy E of filtered Rayleigh-scattered light is used to quantify important properties of a flow like temperature, pressure, velocity, number density, etc. The data we use are Reynolds-Averaged Navier-Stokes CFD simulations of the supersonic shock wave dominated flow field in the experiment mentioned previously <cit.>. The output from these CFD simulations contains the aforementioned flow properties on a rectangular 3D grid of approximately 33 million points. Given these flow properties and the experimental setup, we can calculate the value of E at each point. Doing this allows us to study the relationship between E, the flow properties, and the experimental setup. A better understanding of these relationships will aid in experiment optimization, interpretation, and analysis of the measured energy E in future real FRS experiments. 

To calculate E at each of the 33 million points in our simulation, we need to know the flow parameters at each point (obtained from the CFD), the experimental setup (fixed in advance), and the RB scattering lineshape of the light at each point. To obtain the lineshape, one needs a lineshape model. Typically, Tenti's S6 model is used. However, because Tenti's model is comparatively slow, using Tenti's model to compute E would take about 38 days to compute across all 33 million points. Alternatively, as our method is fast, using SVSA, it takes about 3× 10^-4 seconds to compute E at each point and thus, less than three hours to compute across all 33 million points. In Figure <ref>, we plot a single cross-sectional plane of E calculated using the CFD data and either Tenti or SVSA. By eye, there is no discernible difference between Tenti and SVSA. Additionally, we plot the difference between the Tenti and SVSA calculated values of E. We can see from this plot that in absolute magnitude the difference between the two methods is never more than 2× 10^-3 (keeping in mind that the total scale for E has been normalized to range between 0 and 1). Thus, SVSA produces a value of E that is highly concurrent with Tenti's model but does so about 300 times faster. (While calculating the lineshape is 1-100 thousand times faster with SVSA, other parts of the calculation of E are identical using Tenti and SVSA. Thus, the full calculation of E doesn't see as large of an improvement as the calculation of the lineshape alone.)






§ CONCLUSION


In this paper, we developed a method called support vector spectrum approximation (SVSA) that uses machine learning to create efficient and accurate approximations of any existing spectral lineshape model in arbitrary flow regimes.The ability to produce approximations of existing lineshape models very efficiently is paramount for time-sensitive applications. In addition to our example of greatly speeding up a computational FRS study, we expect that SVSA will have application to real-time diagnostic techniques where efficiently calculating the lineshape from a continuous stream of data is necessary <cit.>. Additionally, we expect that SVSA will have application to problems outside of FRS and the RB spectra, for example, in approximating Raman spectra as part of CARS. We expect to explore such extensions in future work. 




§ FUNDING


National Aeronautics and Space Administration No.
 NNX15AI20H subaward by the Virginia Space Grant Consortium No. 19-264-100527-010. 



§ DISCLOSURES
 The authors declare that there are no conflicts of interest related to this article.

natbib


