Data-Driven Impulse Response Regularization via Deep Learning
Carl Andersson Niklas Wahlström Thomas B. Schön
cs.SY cs.LG stat.ML




[breaklines,basicstyle=]
@inproceedingsCarlAndersson2018,
author    = Carl Andersson and Niklas Wahlström and Thomas Bo Schön,
title     = Data-Driven Impulse Response Regularization via Deep Learning,
booktitle = Proceedings of 18th IFAC Symposium on System Identification (SYSID),
address   = Stockholm, Sweden,
pages     = 1-6,
year      = 2018





[
\begin@twocolumnfalse
	C2MSNet: A Novel approach for single image haze removal
    Akshay Dudhane1  Subrahmanyam Murala2

Computer Vision and Pattern Recognition Lab

Indian Institute of Technology Ropar

India

12017eez0001@iitrpr.ac.in

2subbumurala@iitrpr.ac.in

    Accepted January 20, 2018. Received January 3, 2018; in original form October 30, 2017
==========================================================================================================================================================================================



	Please cite this version:

Carl Andersson, Niklas Wahlström, Thomas B. Schön.   Data-Driven Impulse Response Regularization via Deep Learning. In Proceedings of 18th IFAC Symposium on System Identification (SYSID), Stockholm, Sweden,  2018.



.75
	



	
	

	
	  
	We consider the problem of impulse response estimation of stable linear single-input single-output systems. It is a well-studied problem where flexible non-parametric models recently offered a leap in performance compared to the classical finite-dimensional model structures. Inspired by this development and the success of deep learning we propose a new flexible data-driven model. Our experiments indicate that the new model is capable of exploiting even more of the hidden patterns that are present in the input-output data as compared to the non-parametric models.





\end@twocolumnfalse
]	

 
C2MSNet: A Novel approach for single image haze removal
    Akshay Dudhane1  Subrahmanyam Murala2

Computer Vision and Pattern Recognition Lab

Indian Institute of Technology Ropar

India

12017eez0001@iitrpr.ac.in

2subbumurala@iitrpr.ac.in

    Accepted January 20, 2018. Received January 3, 2018; in original form October 30, 2017
==========================================================================================================================================================================================

  
	We consider the problem of impulse response estimation of stable linear single-input single-output systems. It is a well-studied problem where flexible non-parametric models recently offered a leap in performance compared to the classical finite-dimensional model structures. Inspired by this development and the success of deep learning we propose a new flexible data-driven model. Our experiments indicate that the new model is capable of exploiting even more of the hidden patterns that are present in the input-output data as compared to the non-parametric models.




§ INTRODUCTION


Impulse response estimation has for a long time been at the core of system identification. Up until some five to seven years ago, the generally held belief in the field was indeed that we knew all there was to know about this topic. However, the enlightening work by <cit.> changed this by showing that the estimate can in fact be improved significantly by placing a Gaussian Process (GP) prior on the impulse response, which acts as a regularizer. This model-driven approach has since then been further refined <cit.>, where the prior in this case could be interpreted to encode not only smoothness information, but also information about the exponential decay of the impulse response. In this paper we employ deep leaning (DL) to find a suitable regularizer via a method that is driven by data. Fig. <ref> depicts the general idea and the similarity of our method compared to the method based on Gaussian processes.

Deep learning is a  fairly new area of research that continues the work on neural networks from the 1990's. To get a brief, but informative, overview of the field of deep learning we recommend the paper by <cit.> and for a more complete snapshot of the field we refer to the monograph by <cit.>. Deep learning has recently revolutionized several fields, including image recognition (e.g. <cit.>) and speech recognition (e.g. <cit.>). In both fields, deep learning has surpassed domain specific methods and hand-crafted feature design, by making use of large quantities of data in order to learn data-driven neural network models as general function approximators. 

The idea of using neural networks within system identification is certainly not new and they have been a standard tool for a long time, see e.g. <cit.>. However, the current development in deep learning is  different from the past due to advances in computational software and hardware. As a consequence, contemporary neural networks are more reproducible than before which has increased the credibility of the area. Finally, the amount of available data has sky-rocketed which has made it possible to train larger models with better results. We believe that system identification still has lots to gain from using deep learning and this paper is just one concrete example of what can be done.


 



§ BACKGROUND & PROBLEM FORMULATION

Consider a stable single-input single-output time-invariant linear system G^0 relating an input sequence u(t) to an output sequence y(t) according to

    y(t) = G^0(q) u(t) + v(t),

where q denotes the shift operator qu(t) = u(t+1) and v(t) denotes additive noise. The noise is assumed to be Gaussian white noise with zero mean and variance σ^2. The system G^0(q) is represented using a  transfer function

    G^0(q) = ∑_k=1^∞ g_k^0 q^-k,

where g_1^0, g_2^0, … denote the impulse response coefficients of the system. We make use of superscript 0 to denote the true system.

Based on data, i.e. an input sequence u = (u(1),…,u(N))^T and an output sequence y = (y(1),…,y(N))^T both of length N, the task is to compute an estimate that represents the true system G^0(q) as well as possible. 

Traditionally, the transfer function is encoded using a finite number of parameters  (dim() = n ≪ N), for example, via a finite impulse response (FIR) model,

    G(q; ) = B(q) = b_1 q^-1 + … + b_n_b q^-n_b,
	 
or an output-error (OE) model,

    G(q; ) = B(q)/F(q) = b_1 q^-1 + … + b_n_b q^-n_b/f_1 q^-1 + … + f_n_f q^-n_f,

where =  (b_1,…,b_n_b  )^ or =  ( b_1,…,b_n_b,f_1,…,f_n_f )^, respectively.
See <cit.> for a comprehensive list of model structures. With these model structures, the prediction error method can be used to compute a point estimate  of the unknown parameters by solving the following optimization problem,

    = _∑_t=1^N (y(t) - G(q, )u(t))^2 + ^ D ,

where ^ D describes the added regularization term, governed by the regularization matrix D. It has been shown by <cit.> that the use of an effective regularization is more important than the specific choice of model-order, n. Intuitively the use of a GP model and the regularization term ^ D opens up for model selection at a "finer" scale compared to what is possible with the classical finitely parameterized model structures. Adding this kind of regularization is especially important when the number of samples, N, is low compared to the parameters we wish to estimate, n. 

In the case of the FIR model, the optimization problem reduces to a linear least squares problem to find the optimal parameters . For numerical reasons the inverse regularization matrix P=D^-1 is often used in place of D. This gives rise to the following analytic solution,

	
	
    =   (PR+I_n ) ^-1 P F_N,

	where 
	
    R    = ∑_t = n+1^N φ(t) φ(t)^, 
    φ(t)    =  { u(t-1) … u(t-n)  }^, 
    
    		F_N    = ∑_t = n+1^N φ(t) y(t).


We will throughout this paper denote then estimate (<ref>) by (P) to stress its dependence on the inverse regularization matrix P. As a special case, with P = 0, we have the least squares solution, which we denote .

One question still remains though; how do we find the inverse regularization matrix P? One of the most general ideas is to let it  depend on u and y. This was done in <cit.> by modelling the impulse response as a Gaussian process. A Gaussian process is known to be a very flexible prior, even so, since the model only depends on a low number of hyper-parameters, typically one for a lengthscale and one for some noise, it heavily depends on the specific model we choose for the covariance function. These hyper-parameters are replaced by a point estimate obtained by maximizing the marginal likelihood of the observed data, a procedure known as Empirical Bayes <cit.>. The regularization matrix will thus implicitly depend on u and y via the hyper-parameters. We explicitly denote this dependence by P=f_GP(u,y). This method is also explained in more detail in Section <ref>.

We instead propose an arguably even more flexible model by parametrizing the inverse regularization matrix P with a neural network P = f_DL(u,y). In contrast to the Gaussian process model these parameters are computed by training the model with lots of training data consisting of an input sequence, an output sequence and the true impulse response for either real or simulated systems. Compared to the GP model this is a more data-driven approach to the problem which also makes it possible to use existing techniques from deep learning when building and training the model.







§ REGULARIZATION USING 
 GAUSSIAN PROCESS


The Gaussian prior offers a natural way of encoding the smoothness and decay characteristics that we find in the impulse response from a stable linear system. The specific details of these characteristics are tuned via the hyper-parameters λ. The resulting GP prior can be written as

    ∼ p_λ() =  𝒩( 0, P^λ),

where P^λ in this Bayesian setting is equal to the inverse regularization matrix in (<ref>). The matrix P^λ is related to the covariance function of the Gaussian Process as k_λ (i, j)= P_ij^λ,

    P_ij^λ = k_λ (i, j),
 
where P_ij^λ denotes then entry on row i and column j of P^λ. With the aid of measured data we can select a point estimate of the hyper-parameters λ by maximizing the marginal log-likelihood,

    λ = _λlog∫  p(y;u) p_λ() d,

where 

    p(y;u)  = ∏_t=1^N p(y(t) ;u)  =   ∏_t=1^N 𝒩(y(t) ^φ(t),σ^2).

As a direct consequence of this, the optimal inverse regularization matrix P^λ implicitly depends on the input sequence u and the output sequence y. The equation (<ref>) then describe the function P = f_GP(u,y). Using this inverse regularization matrix together  with (<ref>), we obtain an estimate of the FIR parameters that has better accuracy and is more robust than the non-regularized approach <cit.>.





§ REGULARIZING USING DEEP LEARNING


In contrast to previous work where the regularization matrix only depends on a few hyper-parameters we instead model the regularization matrix directly with a neural network that depends on a large number of parameters η  according to

    P = f_DL(u,y;η).

To select specific values for all these parameters we start by formulating the mean squared error (MSE) of the estimate (P) in (<ref>), the so-called estimation error,

    MSE ( ) = ‖( f_DL(u,y;η) ) - _0  ‖^2,

where ‖·‖^2, denotes the 2-norm and _0 denotes the true FIR parameters which corresponds to the truncated true impulse response. Here we make use of the neural network model (<ref>) of P when forming the MSE. The parameters η to use in (<ref>) are found by simply minimizing  this estimation error,

    η = _η1/M∑_i=1^M  ‖ (f_DL(u^(i),y^(i);η)  ) - _0^(i)‖^2,
 
where M is the number of training examples. To set the terminology straight we use the term training for the minimization of the estimation error w.r.t. η whereas estimation refers to the computation of a point estimate of the FIR parameters by minimizing the one step prediction error w.r.t. .

The following sections describe our method and an overview is provided in Algorithm <ref>.



 §.§ Regularization model


We know that P must be positive  semi-definite. Inspired by this fact we choose P as a weighted sum of rank-one positive definite matrices. The idea behind this choice is to have a representation that is flexible enough to represent all possible regularization matrices and at the same time encode the knowledge that the optimal regularization matrix is a rank-one matrix (see Appendix <ref>).

Our rank-one matrices are constructed as outer products of a vector s_i with itself where the elements of the vectors are free parameters. The weights, w_i, that weight our rank-one matrices, are modelled as the output of a softmax layer from a neural network. The input to this neural network is u and y as well as the nonregularized least squares solution, which only depends on u and y. Hence, we have,

    P  = ∑_i=1^n_m w_i(u,y,;η') s_i s_i^,

where n_m is the number of rank-one matrices used to represent the regularization matrix and η' is the parameters of the neural network. We use η = {η',s_1,…,s_n_m} to collect all the parameters in the regularization matrix model.



 §.§ Neural network model


We have made use of four fully connected layers in our neural network, where the final layer is a softmax layer producing  weights between 0 and 1. 
The other activation functions for the fully connected layers are Rectified Linear Units (ReLU) which are defined as ReLU(x) = max(0,x). A typical layer of the network is thus written as 

h^(i+1) = ReLU(W h^(i) + b), 
where h^(i) denotes the input to the layer, h^(i+1) denotes the output from the layer, W denotes a so-called weight matrix of dimensions dim(h^(i+1))×dim(h^(i)) 
and b denotes a so-called bias term of dimension dim(h^(i)). Both the weight matrices and the bias terms are part of the parameters η' describing the network, i.e. W^(i)∈η' and b^(i)∈η'. To regularize the training procedure we add a dropout layer after the softmax layer which with 30% probability sets a  weight to zero during training. This is a standard technique to avoid overfitting in NN <cit.>. The input to the network is the input u and output y sequence of the system we want to estimate and the corresponding non-regularized least squares solution . 

The resulting network is schematically illustrated in Fig. <ref>.






 §.§ Normalizing the data


A key aspect to successfully train neural networks is the normalization of the input and output of the network by subtracting the mean and dividing with the standard deviation. We notice that we can, without loss of generality, normalize each data example of y and u if we at the same time do the corresponding scaling of the impulse response  to keep the analytic relationship intact.  

The non-regularized least square solution, , that we use as an input to the network is also straightforward normalize with the statistics form the true impulse response calculated from the training data. Finally, we want to normalize the networks dependence of the vectors s_i. The outer product of these vectors should correspond to the optimal regularization (see Section <ref> and Appendix <ref>). To enforce that they follow the same statistics, we initialize the vectors s_i with unit Gaussian and then multiply we the standard deviation and add the mean of the true impulse response. 






§ EXPERIMENT


The model explained in Section <ref> is implemented using Tensorflow <cit.> and our implementation of the model is available on GitHub[https://github.com/carl-andersson/impest]. 



 §.§ Simulate data using 

To train the model, we use an artificial distribution over real systems to produce input and output sequences along with their true impulse responses. The artificial distribution over systems we use is MATLAB's  function. This function is not ideal as was recently pointed out by for example <cit.>. Hence, there is potential to further improve the results by making use of better data.

To generate data we use the same method as <cit.> with some minor alterations concerning the signal to noise ratio (SNR). The full procedure follows as:

	
  * Sample a system of order 30 using a slightly modified version of MATLAB's  where the probability of having an integrating pole is zero. 
	
  * Sample N = 125 timesteps from the system at a sampling rate of 3 times the bandwidth, i.e.,
		
   	
   	
  * Calculate the true FIR parameters as the truncated impulse response to length n=50.
   	
  * Sample white noise as an input sequence u and get the corresponding output sequence ^*. Add noise to the output with the SNR drawn randomly from a uniform distribution between 1 and 10, i.e., the noise added has variance that is between 1 and 1/10 of the variance in the output sequence. 
   	
  * Repeat all steps until we have M examples from M different systems.


These M training examples are then used in Equation (<ref>) and are called the training data. Using the same method we generate a validation set which we also split in two sets depending on the SNR, one with SNR larger than 5.5 and one set with SNR less than 5.5 with roughly M_v≈ 5000 examples in each.





 §.§ Evaluation metrics

To evaluate the performance of the model we use a metric that is calculated as the mean squared error of the estimate  normalized with the mean squared error of the least squares solution without any regularization. We denote this metric by S, i.e.,

    S = 1/M_v∑_i=1^M_v [ ‖ (f_DL(u^(i),y^(i);η)  ) - ^(i)_0 ‖^2/‖^(i)-^(i)_0  ‖^2 ],

where M_v denotes the number of examples in the validation set. This metric makes sure that each impulse response gets equal weighting when computing the performance of the algorithm and measures the average effect of the regularization. A perfect match for this measure corresponds to a measure of 0. 

<cit.> make use of a slightly different metric defined as

    S = 1/M_v∑_i=1^M_v 100 ( 1 -  [ ‖ (f_DL(u^(i),y^(i);η)  ) - ^(i)_0  ‖^2 /‖_0^(i) - _0^(i)‖^2 ]  )

where _0^(i) is  the mean of _0^(i). This metric averages over a so-called 'model fit', i.e. how well the estimated parameters fit the true impulse response. Besides the shifting and scaling in S, the only difference between the two metrics is the normalization factor used, where S is normalized with the least squares estimation error and S is normalized with the variance of the true impulse response. We have empirically observed that the terms in S might vary a lot between different examples and the average might thus be dominated by a few examples leading to measures that are hard to interpret. Our slightly modified metric S will on the other hand measure the average effect of using a regularization method compared to not using a regularization method, which seem to result in a more stable performance indicator.





 §.§ Simulation results


The model is trained using M = 1000000 training examples for roughly 2.5 hours using a desktop computer with an Nvidia Titan Xp GPU. The chosen hyperparameters of the model is described in Figure <ref>. Note that while training require a GPU with large memory, the evaluation can easily be done in CPU on an ordinary computer. We are using early stopping as a stopping criteria even though the model essentially does not seem to overfit with this amount of training data. The model is not very dependent on the number of examples in the training data either. Even with M=10000 it managed to achieve comparable results to previous methods.

Fig. <ref> shows a subset of the matrices s_is_i^T from (<ref>) after training. Note that the matrices have an oscillating pattern with different periods and a decay towards zero for parameters with high index (lower right corner). Fig. <ref> shows three different regularization matrices for an example in the validation dataset. We can see that the deep learning regularization seems to capture the behaviour of the optimal regularization matrix fairly well. In Fig. <ref> we can compare the estimates from the different regularization approaches.  



The trained model does not produce a useful regularization matrix for all examples. In cases where it fails the neural network seems to fail in the same way for all examples by producing a similar regularization matrix for each example with a bad performance for S as consequence. Despite this problem, the model manages to produce average results which are comparable or better than previous methods which is reflected in Table <ref>. 
We can see that the performance decreases when the SNR gets larger. This is due to the improved effectiveness of the least squares estimate and we are thus less dependent on the regularization. For comparison we also present the result for the optimal regularization which of course is unachievable since it depends on the true impulse response, but it is still an useful lower bound.




 §.§ Real data results

To show that our method at least does not give unreasonable estimates for real systems and data we test our method on data measured at a processing plant. See e.g. <cit.> for an introduction to this problem area.  We do not know the true parameter values for these systems, implying that we cannot evaluate the performance of the estimates. We use the same trained network as we evaluated in the previous section. In Fig. <ref> we simply present an input sequence and the corresponding output sequence from a real system together with the estimates produced by our method compared to using least squares without regularization. We note also that the results seems reasonable and that the regularization removes a lot of noise in the estimation. 














§ CONCLUSION

In this paper we present a method to regularize an impulse response estimation problem. We train a model with simulated data in a data-driven fashion to encode this regularization with parameters in a neural network. A trained model can then be used to improve the mean squared error of the new estimations. 
The results of our method seems promising and there is plenty of scope for future work along this line of research, both when it comes to impulse response estimation, but also for other problems. We find it especially interesting that this model can mimic the optimal regularization matrix to higher degree than previous methods which we believe is the reason for why it sometimes produce better estimates.

Although training our model is quite time consuming, estimating the impulse response using our method is very fast since it only involves a couple of matrix  multiplications to compute the regularization matrix, whereas the method of <cit.> needs to solve an optimization problem for each example. 



§ FUTURE WORK

We are planning to further investigate how one can make use of real data from e.g. the process industry. This would make it possible to use large amounts of collected data to improve the estimated parameters in a data-driven manner. The process industry has a lot of data available and makes extensive use of linear models.

The idea of learning a prior by representing it with a regularization matrix in the form of a neural network is not unique to the problem of estimating the impulse response. It could easily be generalized to other situations where the least squares solution is available but the prior of the solution is either unknown or intractable. If one can simulate many such systems at low cost, or have data from true systems available, formulating a regularization matrix as a neural network might be a tractable way of regularizing the estimate.

The presented approach can easily be extended to multi-input multi-output systems where the only difference is that the dimension of the input and output sequences and the parameters  increases. The deep structure of the model automatically induces any relevant connection between the different system input and output components present in the training data. 

Finally we want to stress that this is just one example of what one might do with deep learning in system identification. There might and should be other areas where it is possible to make use of either simulated or real data to improve standard methods, or invent new methods for system identification. For example it might be worth looking into Recurrent Neural Networks (RNN) such as Long Short Term Memory (LSTM)<cit.>, Gated Recurrent Units (GRU) <cit.> or Stocastic Recurrent Neural Network (SRNN)<cit.> and apply it in a system identification setting or even to bring some of the system identification knowledge of dynamical systems to the field of deep learning. 



§ ACKNOWLEDGEMENTS

We would like to thank Professor Krister Forsman at Perstorp and the Norwegian University of Science and Technology for providing us with the real world data. We would also like to thank Dr Jack Umenberger for useful feedback on an early draft of this work.


	abbrvnat
	             





§ OPTIMAL REGULARIZATION


The optimal regularization matrix is a term coined by <cit.>. It corresponds to the regularization matrix that is optimal in the sense that it minimizes the mean squared error. The optimal regularization matrix can be written as, P = σ^-2_0 _0^ <cit.>, where σ^2 is the variance of the additive noise v(t) in (<ref>), _0 is the true impulse response of the system without noise. 



