Implementation of algorithms for tuning parameters in regularized least squares problems in system identification

TLDR: This work investigates implementation of algorithms for solving the hyper-parameter estimation problem that can deal with both large data sets and possibly ill-conditioned computations and proposes a QR factorization based matrix-inversion-free algorithm to evaluate the cost function in an efficient and accurate way.

About: This article is published in Automatica.The article was published on 2013-07-01 and is currently open access. It has received 126 citations till now. The article focuses on the topics: QR decomposition & System identification.

Figures

Fig. 1. Histogram-plot of log10 cond(σ2In +LT ΦNΦTNL) over 1000 data sets in the second experiment of (Pillonetto & Nicolao, 2011; Chen et al., 2011). Here, we consider the stable spline kernel (Pillonetto & Nicolao, 2010) with the optimal hyper-parameter estimate α̂ from solving (6).

Table 2 Average total time required to compute the cost function in (6) while solving (6) with the DC kernel (31c) over the 2500 data sets in S1D1. The unit of all figures in the table is “second”. Items for Algorithms 1 and 2 takes the form of “a+b”, where a is the preprocessing time and b is the total time to compute the reformulated cost function while solving (6).

Frequently asked questions

Q1. What have the authors contributed in "Implementation of algorithms for tuning parameters in regularized least squares problems in system identification" ?

There is recently a trend to study linear system identification with high order finite impulse response ( FIR ) models using the regularized least-squares approach. 

Q2. What is the numerical accuracy of the cost function in (6)?

When seeking efficient algorithms to compute the costfunction in (6) for very large N and N À n, the numerical accuracy depends on the conditioning and the magnitude of the matrices P(α) and ΦNΦTN . 

Q3. What is the problem of the matrixAT A?

For the impulse response estimation problem, the matrixAT A = σ2In +LT ΦNΦTNL (17)can be ill-conditioned due to the following two problems:• 

Q4. What is the command used to solve the nonconvex optimization problem (6)?

The command fmincon is used here to solve the nonconvex optimization problem (6) with the trust region reflective algorithm selected. 

Q5. How many data sets are used to estimate the FIR model?

For each data set, the authors aim to estimate FIR model (3) with n = 125 using the regularized least squares (5) including the empirical Bayes method (6). 

Q6. What is the condition number of the “DC” kernel?

For the “DC” kernel (31c), further assume 0.72≤ λ < 1 and−0.99≤ ρ ≤ 0.99 so that the condition number of the DC kernel is smaller than 2.0×1020. • 

Q7. What is the way to solve the least squares problem?

It is well-known, e.g. (Golub & Van Loan, 1996, Section 5), that the least squares problem (15) can be solved more accurately with the QR factorization method than the Cholesky factorization method, when the condition number of AT A defined in (17) is ill-conditioned. 

Q8. What is the basic problem to estimate linear regressions?

The linear least squares problem to estimate linear regressions is one of the most basic estimation problems, and there is an extensive literature around it, e.g. (Rao, 1973; Daniel & Wood, 1980; Draper & Smith, 1981). 

Q9. What is the reason why the extra constraints are used in Monte Carlo?

So one may question if the extra constraints can eventually cause performance loss in the regularized least squares estimate (5b). 

Q10. What is the magnitude of LT NTNL?

The magnitude of LT ΦNΦTNL can be very large if the element of α that controls the magnitude of P(α) is large, which is often the case for the stable spline kernel (Pillonetto & Nicolao, 2010; Pillonetto et al., 2011). 

Q11. what is the maximum likelihood method to estimate from (4)?

It is the maximum likelihood method to estimate α from (4) under the (Bayesian) assumptions that θ is Gaussian with zero mean and covariance matrix P(α) and VN is Gaussian with zero mean and covariance matrix σ2IN . 

Q12. What is the algorithm for estimating the impulse response?

In this case, as discussed in (Chen, Ohlsson & Ljung, 2012), the authors can first estimate, with the Maximum Likelihood/Prediction Error Method e.g. (Ljung, 1999), a low-order “base-line model” that can take care of the dominating part of the impulse response. 

Q13. What is the algorithm for estimating the FIR model?

The authors then use regularized least squares (based on Algorithm 2) to estimate an FIR model with reasonably large n, which should capture the residual (fast decaying) dynamics (Chen, Ohlsson & Ljung, 2012).