## 2 System Model and Problem Formulation  
Let us consider a linear single-input/single-output (SISO) system, whose output at the discrete-time index \(t\) is given by  
\[d(t) =\mathbf{h}^{T}\mathbf{x}(t)+w(t) \tag{1}\] \[=y(t)+w(t),\]  
where \(d(t)\) is also known as the desired signal,  
\[\mathbf{x}(t)=\left[ \begin{array}{cccc}x(t)&x(t-1)&\cdots&x(t-L+1)\end{array}  \right]^{T}\]  
is a vector containing the most recent \(L\) time samples of the zero-mean input signal \(x(t)\), the superscript \({}^{T}\) denotes the transpose operator, \(\mathbf{h}\) is the unknown system impulse response (a column vector of length \(L\)), and \(w(t)\) is the zero-mean additive noise. It is assumed that we deal with real-valued signals, and that \(x(t)\) and \(w(t)\) are uncorrelated. Then, the objective of the so-called SISO system identification problem is to identify or estimate \(\mathbf{h}\), in the best possible way, given the input signal vector, \(\mathbf{x}(t)\), and the output signal, \(d(t)\). From (1), we see that the signal-to-noise ratio (SNR) of the SISO system is  
\[\mathrm{SNR}=\frac{\mathbf{h}^{T}\mathbf{R}_{\mathbf{x}}\mathbf{h}}{\sigma_{w }^{2}}, \tag{2}\]where \(\mathbf{R}_{\mathbf{x}}=E\big{[}\mathbf{x}(t)\mathbf{x}^{T}(t)\big{]}\) and \(\sigma_{w}^{2}=E\big{[}w^{2}(t)\big{]}\) are the covariance matrix of \(\mathbf{x}(t)\) and variance of \(w(t)\), respectively, with \(E[\cdot]\) denoting mathematical expectation. Let us assume that \(L=L_{1}L_{2}\). Without loss of generality, we may also assume that \(L_{1}\geq L_{2}\). Then, the system impulse response can be decomposed as [16]  
\[\mathbf{h}=\sum_{l=1}^{L_{2}}\mathbf{h}_{2,l}\otimes\mathbf{h}_{1,l}, \tag{3}\]  
where \(\otimes\) is the Kronecker product [31], and \(\mathbf{h}_{1,l}\) and \(\mathbf{h}_{2,l}\) are short impulse responses (column vectors) of lengths \(L_{1}\) and \(L_{2}\), respectively. Another way to write (3) is  
\[\mathbf{h} =\text{vec}\left(\sum_{l=1}^{L_{2}}\mathbf{h}_{1,l}\mathbf{h}_{2,l}^{T}\right) \tag{4}\] \[=\text{vec}\left(\mathbf{H}_{1}\mathbf{H}_{2}^{T}\right)=\text{ vec}(\mathbf{H}),\]  
where \(\text{vec}(\cdot)\) is the vectorization operation, which consists of converting a matrix into a vector,  
\[\mathbf{H}_{1} =\big{[} \begin{array}{cccc}\mathbf{h}_{1,1}&\mathbf{h}_{1,2}& \cdots&\mathbf{h}_{1,L_{2}}\end{array} \big{]},\] \[\mathbf{H}_{2} =\big{[} \begin{array}{cccc}\mathbf{h}_{2,1}&\mathbf{h}_{2,2}& \cdots&\mathbf{h}_{2,L_{2}}\end{array} \big{]}\]  
are matrices of sizes \(L_{1}\times L_{2}\) and \(L_{2}\times L_{2}\), respectively, and \(\mathbf{H}\) is the equivalent matrix (of size \(L_{1}\times L_{2}\)) representation of \(\mathbf{h}\). However, very often in practice, the column rank of \(\mathbf{H}\) is far from its maximum value, i.e., \(\text{rank}(\mathbf{H})=P\ll L_{2}\). As a result, (3) can be expressed as  
\[\mathbf{h}=\sum_{p=1}^{p}\mathbf{h}_{2,p}\otimes\mathbf{h}_{1,p} \tag{5}\]  
and \(y(t)\) in (1) becomes  
\[y(t)=\left(\sum_{p=1}^{p}\mathbf{h}_{2,p}\otimes\mathbf{h}_{1,p}\right)^{T} \mathbf{x}(t). \tag{6}\]  
Let us define the symmetric impulse response [29]:  
\[\mathbf{h}^{+}=\frac{\sum_{p=1}^{p}\big{(}\mathbf{h}_{2,p}\otimes\mathbf{h}_{ 1,p}+\mathbf{h}_{1,p}\otimes\mathbf{h}_{2,p}\big{)}}{2}. \tag{7}\]  
We should note that this particular symmetric characteristic is not related to the samples of the impulse response (e.g., like in linear-phase finite-impulse-response filters), but to a specific combination of the two sub-impulse responses \(\mathbf{h}_{1,p}\), and \(\mathbf{h}_{2,p}\). Thus, denoting by \(\mathbf{h}_{p}(\mathbf{h}_{1,p},\mathbf{h}_{2,p})\) the \(p\)th term of the sum from the numerator of (7), we can write  
\[\mathbf{h}_{p}(\mathbf{h}_{1,p},\mathbf{h}_{2,p}) =\mathbf{h}_{p}(\mathbf{h}_{2,p},\mathbf{h}_{1,p}), \tag{8}\] \[\mathbf{h}_{p}(\mathbf{h}_{1,p},\mathbf{h}_{1,p}) =2\mathbf{h}_{1,p}\otimes\mathbf{h}_{1,p}, \tag{9}\]  
which define the particular symmetry mentioned before.