# Classification of Event-Related Potentials with Regularized Spatiotemporal LCMV Beamforming
## 2 Materials and Methods
### Covariance Matrix Regularization
#### 2.3.4 Kronecker-Toeplitz-Structured Covariance Estimation  
The question remains how to estimate \(\hat{S}\) and \(\hat{T}\). Although the flip-flop and non-iterative flip-flop algorithms [44; 45; 46] can estimate Kronecker or Kronecker-Toeplitz-structured covariances, new results show that a fixed point iteration is more efficient [47; 48]. After each iteration, the spatial and temporal covariance matrices are scaled to unit variance to ensure that the fixed-point iteration converges. Finally, shrinkage can also be introduced in the fixed-point iteration to improve stability and achieve more robust regularization [42; 48; 49; 50].  
The spatial and temporal covariance matrices are shrunk at every fixed-point iteration with shrinkage factors \(\beta_{k}\) and \(\gamma_{k}\) before matrix inversion in the next iteration. Combined, this leads to the iterative estimation algorithm described by the following equations:  
\[\hat{S}_{k+1}=\frac{1}{n}\sum_{i=1}^{n}X_{i}^{\intercal}\hat{T}_{k}^{+}X_{i} \tag{15a}\]\[\hat{T}_{k+1}=\frac{1}{n}\sum_{i=1}^{n}X_{i}\hat{S}_{k}^{+}X_{i}^{\intercal} \tag{15b}\] \[\hat{S}_{k+1}^{(\beta)}=(1-\beta_{k+1})\hat{S}_{k+1}+\beta_{k+1}\frac{ \operatorname{Tr}(\hat{S}_{k+1})}{c}\mathbb{I}\] (16a) \[\hat{T}_{k+1}^{(\gamma)}=(1-\gamma_{k+1})\hat{T}_{k+1}+\gamma_{k+1}\frac{ \operatorname{Tr}(\hat{T}_{k+1})}{s}\mathbb{I}\] (16b) \[\hat{S}_{k+1}=\frac{c}{\operatorname{Tr}\left[\hat{S}_{k+1}^{( \beta)}\right]}\hat{S}_{k+1}^{(\beta)}\] (17a) \[\hat{T}_{k+1}=\frac{s}{\operatorname{Tr}\left[\hat{T}_{k+1}^{( \gamma)}\right]}\hat{T}_{k+1}^{(\gamma)} \tag{17b}\]  
\(\hat{S}_{0}\) and \(\hat{T}_{0}\) can be initialized to any positive definite matrix. We choose to use the identity matrices \(\mathbb{I}^{\times<c}\) and \(\mathbb{I}^{\times s}\). After each iteration, all diagonals of \(\hat{R}_{k+1}\) are set to their mean values to ensure that \(\hat{R}_{k+1}\) and \(\hat{T}_{k+1}\) are Toeplitz-structured.  
[PERSON] et al. [51] show that the LOOCV estimates for the optimal values of \(\beta_{k+1}\) and \(\gamma_{k+1}\) also yield a closed-form solution for the Kronecker fixed-point-iteration algorithm:  
\[\beta_{k+1}=1-\frac{\frac{n}{n-1}\operatorname{[PERSON]}(\hat{S}_{k+1}^{2})-\frac{2 }{c}\left[\operatorname{[PERSON]}(\hat{S}_{k+1})\right]^{2}+\frac{1}{c} \operatorname{[PERSON]}(\hat{S}_{k+1}^{2})-\frac{1}{n(n-1)}\sum_{i=1}^{n}\left[ \operatorname{[PERSON]}(X_{i}\hat{T}_{k}^{+}X_{i}^{\intercal})^{2}\right]}{\frac{n^ {2}-2n}{(n-1)^{2}}\operatorname{[PERSON]}(\hat{S}_{k+1}^{2})-\frac{2}{c}\left[ \operatorname{[PERSON]}(\hat{S}_{k+1})\right]^{2}+\frac{1}{c}\operatorname{[PERSON]}( \hat{S}_{k+1}^{2})+\frac{1}{n(n-1)^{2}}\sum_{i=1}^{n}\left[\operatorname{[PERSON]}( X_{i}\hat{T}_{k}^{+}X_{i}^{\intercal})^{2}\right]} \tag{18a}\] \[\gamma_{k+1}=1-\frac{\frac{n}{n-1}\operatorname{[PERSON]}(\hat{T}_{k+1}^{2})-\frac{2 }{s}\left[\operatorname{[PERSON]}(\hat{T}_{k+1})\right]^{2}+\frac{1}{s} \operatorname{[PERSON]}(\hat{T}_{k+1}^{2})-\frac{1}{n(n-1)}\sum_{i=1}^{n}\left[ \operatorname{[PERSON]}(X_{i}^{\intercal}\hat{S}_{k}^{+}X_{i})^{2}\right]}{\frac{n^ {2}-2n}{(n-1)^{2}}\operatorname{[PERSON]}(\hat{T}_{k+1}^{2})-\frac{2}{s}\left[ \operatorname{[PERSON]}(\hat{T}_{k+1})\right]^{2}+\frac{1}{s}\operatorname{[PERSON]}( \hat{T}_{k+1}^{2})+\frac{1}{n(n-1)^{2}}\sum_{i=1}^{n}\left[\operatorname{Tr}( X_{i}^{\intercal}\hat{S}_{k}^{+}X_{i})^{2}\right]} \tag{18b}\]  
The shrinkage parameters \(0<\beta_{k+1}<1\) and \(0<\gamma_{k+1}<1\) should be re-determined after each iteration. The oracle approximation shrinkage method can also be used to determine \(\beta_{k+1}\) and \(\gamma_{k+1}\)[51, 52] but performs worse for spatiotemporal EEG data since not all assumptions are met.