# Coherence-based Input Design for Sparse System Identification

Javad Parsa∗, Cristian R. Rojas∗, and Håkan Hjalmarsson∗ The authors are with
the Division of Decision and Control Systems, School of Electrical Engineering
and Computer Science, KTH Royal Institute of Technology, 100 44 Stockholm,
Sweden javadp, crro<EMAIL_ADDRESS>work was supported by the Swedish
Research Council research environment NewLEADS—New Directions in Learning
Dynamical Systems, contract 2016-06079; and Wallenberg AI, Autonomous Systems
and Software Program (WASP), funded by Knut and Alice Wallenberg Foundation.

###### Abstract

The maximum absolute correlation between regressors, which is called mutual
coherence, plays an essential role in sparse estimation. A regressor matrix
whose columns are highly correlated may result from optimal input design,
since there is no constraint on the mutual coherence, so when this regressor
is used to estimate sparse parameter vectors of a system, it may yield a large
estimation error. This paper aims to tackle this issue for fixed denominator
models, which include Laguerre, Kautz, and generalized orthonormal basis
function expansion models, for example. The paper proposes an optimal input
design method where the achieved Fisher information matrix is fitted to the
desired Fisher matrix, together with a coordinate transformation designed to
make the regressors in the transformed coordinates have low mutual coherence.
The method can be used together with any sparse estimation method and in a
numerical study we show its potential for alleviating the problem of model
order selection when used in conjunction with, for example, classical methods
such as AIC and BIC.

###### Index Terms:

System identification, input design, sparse estimation, mutual coherence

## Notation

Vectors and matrices are denoted by bold lower-case and upper-case letters,
respectively. The support of a vector $\theta$, which is denoted by
supp$({\mbox{\boldmath$\theta$}})$, is the set of indices of non-zero elements
in $\theta$. $\|{\bf B}\|_{\infty}$ and $\|{\bf B}\|_{F}$ are the maximum
absolute value between all elements in matrix ${\bf B}$ and the Frobenius norm
of the matrix ${\bf B}$, respectively. The vectorization operator, i.e.,
$\text{vec}\colon{\mathbb{R}}^{a\times b}\rightarrow{\mathbb{R}}^{ab}$,
converts a matrix to a vector by stacking the columns of a matrix on top of
one another, and
$\text{vec}^{-1}\colon{\mathbb{R}}^{ab}\rightarrow{\mathbb{R}}^{a\times b}$
refers to the inverse of the vectorization operator. Furthermore,
${\mbox{\boldmath$\phi$}}_{i}$ denotes the i-th column of matrix $\Phi$, and
$|\cdot|$ represents the absolute value.

##  I Introduction

System identification is the science of building dynamic models of systems
from data [1]. A model is an approximation of the real-life behavior of the
system, and it does not provide an exact description of the system. As a
consequence, different models can be found under various assumptions. A simple
but widely used structure is the linear regression model, given by

${\bf y}={\mbox{\boldmath$\Phi$}}{\mbox{\boldmath$\theta$}}+{\bf e},$ (1)

where ${\bf y}\in{\mathbb{R}}^{N\times 1}$,
${\mbox{\boldmath$\theta$}}\in{\mathbb{R}}^{n_{\theta}\times 1}$,
${\mbox{\boldmath$\Phi$}}\in{\mathbb{R}}^{N\times n_{\theta}}$ and ${\bf
e}\in{\mathbb{R}}^{N\times 1}$ denote the observations, parameter vector,
regressor matrix and noise, respectively.  This paper considers models where
$\Phi$ has a lower Toeplitz structure. This covers finite impulse response
(FIR) models

$y(t)=\sum_{k=1}^{n_{\theta}}\theta_{k}u(t-k)+e(t),\quad t=1,\ldots,N$ (2)

in which case $\Phi$ is constructed from the input sequence $\\{u(t)\\}$.
However, this structure also covers models of the form

$y(t)=\sum_{k=1}^{n_{\theta}}\theta_{k}u^{\prime}(t-k)+e(t),\quad
t=1,\ldots,N$

where

$u^{\prime}(t)=F(q)u(t)$ (3)

where $F(q)$ is a fixed transfer function ($q$ is the discrete-time time-shift
operator $qu(t)=u(t+1)$). In this case $\Phi$ is constructed from
$\\{u^{\prime}(t)\\}$. This means that the considered model covers, among
others, Laguerre, Kautz, and generalized orthonormal basis function expansion
models in which case $F(q)$ has no zeros and the denominator contains the
poles of the Laguerre, Kautz, or generalized basis function expansion model
[2, 3]. Thus the model covers several of the most commonly used models in
system identification.

In addition, we will consider the case where $\theta$ is sparse, meaning that
only, say, $s$ elements of this vector are non-zero, and their position is
unknown. This setting appears in numerous applications. For example, in
wireless communication, a common challenge is channel fading, where the
received signal comprises multiple attenuated and delayed copies of the
transmitted signal [4]. In such scenarios, the channel can be modeled using a
linear FIR model with inherently sparse impulse response [5]. Sparse
estimation is a family of methods tailored to this problem, balancing the
trade-off between model accuracy and complexity in system identification by
finding a parameter vector with a suitable number of zero elements, which
results in a model with low order and reduced complexity. A substantial
literature has been published on sparse system identification, including works
such as [6, 7, 8]. For example, in [8], a Bayesian approach for sparse system
identification is proposed, while the parameters of a linear regression model
are estimated using a coordinate transformation in [6]. It is also worth
pointing out that there may be other prior information present than that the
parameter vector is sparse. A very common problem is that of model order
selection, in which case it is known that the first elements of the parameter
vector are non-zero but the exact number of non-zero elements is not known. We
will illustrate in Section V that also this problem is alleviated by the
proposed method.

The mutual coherence of the regressor matrix $\Phi$, defined as

$\mu_{{\mbox{\boldmath$\Phi$}}}=\max_{i\neq
j}\frac{|\langle{\mbox{\boldmath$\phi$}}_{i},{\mbox{\boldmath$\phi$}}_{j}\rangle|}{\|{\mbox{\boldmath$\phi$}}_{i}\|_{2}\|{\mbox{\boldmath$\phi$}}_{j}\|_{2}},$
(4)

where $\\{{\mbox{\boldmath$\phi$}}_{j}\\}$ are the columns of $\Phi$, plays an
important role in sparse estimation, as evidenced by extensive publications on
the subject, e.g. [9, 10, 11, 12]. These works consistently highlight the
importance of maintaining low mutual coherence to achieve accurate sparse
estimation. Now, in system identification, the columns of the regressor matrix
in (1) are typically taken as time-shifted versions of the same signal, e.g.
the input. When this signal has strong temporal correlation, this leads to a
high mutual coherence, cf. the definition (4). This problem may be exacerbated
by optimal input design which often results in an input that is tuned to
rather narrow frequency bands, implying high temporal correlation, aimed at
revealing system properties of importance for the intended model use. In
Section V this is illustrated for a control design problem where a frequency
band covering the cross-over frequency is excited. The purpose of this
contribution is to, in addition to other objectives, include the requirement
that sparse model estimation should be promoted in optimal input design.

Optimal input design typically shapes the Fisher Information Matrix (FIM)
using the input signal as a decision variable to meet some design objectives.
In model (1) with ${\bf e}\sim{\cal N}(0,\sigma^{2}{\bf I})$, the FIM is given
by111We will not need the Gaussian assumption in the following but, with some
misuse of notation, will continue to use the term FIM for ${\bf I}_{F}$, even
though without the Gaussian noise assumption, ${\bf I}_{F}$ is not the FIM but
the inverse covariance matrix of the least-squares estimate of $\theta$ in the
model (1) which can also be used in input design.:

${\bf
I}_{F}=\frac{1}{\sigma^{2}}{\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}}.$
(5)

This implies that the mutual coherence $\mu_{\mbox{\boldmath$\Phi$}}$ in (4)
is directly influenced by the Fisher matrix. This paper introduces an optimal
input design method that considers the mutual coherence by leveraging this
connection. A previous study [13] introduced a method for input design based
on a linear coordinate transformation with two terms. This method employs an
approximation of mutual coherence using the Frobenius norm. Although this
method utilized convex optimization with three sub-steps, it was found to be
slow using generic software such as CVX, especially when working with large
scale problems. Additionally, the use of the Frobenius norm approximation was
not ideal, focusing on reducing correlation between rows rather than columns
[12].

### Contribution

As a remedy to the shortcomings in [13], we propose an alternative approach
with a linear coordinate transformation that instead uses the exact
formulation of mutual coherence. According to Section V, exact formulation
results in a significant reduction in mutual coherence and subsequently lower
error compared to [13]. Also the proposed method provides a closed-form
solution in each step that significantly reduces the computation time, as
evidenced in Section V. Additionally, our proposed approach is versatile,
capable of handling commonly used linear models in system identification, such
as FIR, Laguerre, Kautz, and generalized orthonormal basis function expansion
models, in contrast to the method presented in [13], which is specifically
designed for the FIR model. It is worth noting that the proposed method can be
combined with any standard sparse estimation algorithm. In summary, the
contributions of our paper are:

  1. 1.

A new optimal input design criterion suitable for sparse estimation problems
in the linear regression model (1) is proposed. The method combines a
coordinate transformation of the model (1) with quality and low mutual
coherence constraints.

  2. 2.

A numerical algorithm is proposed to solve the optimization problem. The
method is a three step iterative procedure using proximal mappings and
alternating minimization. In comparison with [13], the algorithm can handle
fixed denominator types of models (see above), and is computationally much
more efficient.

  3. 3.

Similar to [13], the proposed method is numerically evaluated in a model
reference control problem where the actual model order is unknown. However, in
addition to its use together with standard sparse algorithms, as done in [13],
it is also combined with the classical AIC and BIC model order selection
methods. In addition, the proposed method is compared with the method in [13],
both in terms of accuracy and computation time.

### Structure of the Paper

The remaining parts of the paper can be summarized as follows: In the next
section, we give a brief recap of sparse estimation and optimal input design.
The novel optimal input design, incorporating low mutual coherence, is
formulated as an optimization problem in Section III. In Section IV, an
iterative algorithm is proposed for solving the optimization problem. In the
numerical study part, Section V, the proposed method is evaluated against
state-of-the-art algorithms.

##  II Recap of sparse estimation and optimal input design

Before proceeding with the core of the paper, a brief review of sparse
estimation and optimal input design is in place.

###  II-A Sparse estimation

To estimate a sparse parameter vector $\theta$ in the model (1), the error
$\|{\bf y}-{\mbox{\boldmath$\Phi$}}{\mbox{\boldmath$\theta$}}\|_{2}^{2}$ is
minimized with respect $\theta$ while the $\ell_{1}$ norm of $\theta$,
$\|{\mbox{\boldmath$\theta$}}\|_{1}$, is constrained to be smaller than $\tau$
for achieving the desired sparsity. This method is known as LASSO (least
absolute shrinkage and selection operator), and will be denoted by
$P_{1}^{\epsilon}$ in the following. Greedy methods, such as Orthogonal
Matching Pursuit (OMP) [14], and Generalized OMP (GOMP) [15], is another
family of sparse estimation algorithms.

###  II-B Optimal input design

Experimental data are noisy, resulting in uncertainty in the estimated model.
Input design may be utilized to improve the accuracy of an estimator so that
desired criteria are met. Many different methods exist, and we refer to [16,
1, 17, 18] and references therein for a wide range of examples. As for linear
time-invariant systems the per sample Fisher matrix is linear in the input
spectrum, many methods proceed in two steps: i) The optimal input spectrum,
together with the optimal FIM ${\bf I}_{F}^{d}$, is determined using (convex)
optimization, ii) filtering a realization of white noise through a stable
minimum phase spectral factor of the optimal input spectrum. This approach is
called the Frequency Domain Method (FDM) in the following. Least-costly
identification [17] and application-oriented input design [18] are examples of
methods of this type. They minimize experimental effort, such as input power,
subject to quality constraints related to the intended use of the model, where
these constraints are expressed in terms of the per sample FIM. We will employ
this method in the numerical example in Section V where more details will be
provided.

##  III An input design criterion incorporating mutual coherence

The following will be our standing assumption throughout the paper.

###### Assumption 1

The data obeys model (1), with
${\mbox{\boldmath$\theta$}}={\mbox{\boldmath$\theta$}}_{0}$ denoting the true
model parameter. It holds that ${\rm E}[{\bf e}]=0$, ${\rm E}[{\bf e}{\bf
e}^{T}]=\sigma^{2}{\bf I}$, $N\geq n_{\theta}$,
${\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}}>0$ with $\Phi$ being
lower Toeplitz.

As outlined in Section II-B, many optimal input design methods in the first
step determine the optimal input spectrum and the desired FIM ${\bf
I}_{F}^{d}$. Alternatively to realizing the input by way of filtering (see
Section II-B), one may obtain the optimal regressor by minimizing the
difference between the two sides in (5), with ${\bf I}_{F}={\bf I}_{F}^{d}$.
In [19] the Frobenius norm is used to this end. Consequently,

$\min_{{\mbox{\boldmath$\Phi$}}\in{\cal
P}}\|{\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}}-\sigma^{2}{\bf
I}_{F}^{d}\|_{F}^{2},$ (6)

in which ${\cal P}$ denotes the set of $N\times n_{\theta}$ Toeplitz matrices
due to the considered linear FIR model in [19]. The achieved regressor from
(6) can then be used in the identification experiment. An intriguing question
arises: at what level of mutual coherence for the derived regressor from (6)
does it lead to accurate sparse estimation of $\theta$? The following result
addresses this question.

######  Theorem 1 (Stability of $P_{1}^{\epsilon}$ [11])

Consider the model (1), where ${\bf e}\sim{\cal N}(0,\sigma^{2}{\bf I})$.
Here, ${\mbox{\boldmath$\theta$}}_{0}$ denotes a parameter vector with
sparsity characterized by the support set $\Lambda_{0}$, and
$\|{\mbox{\boldmath$\theta$}}_{0}\|_{0}=s$. Assume that the mutual coherence
of the matrix $\Phi$ satisfies the condition
$0\leq\mu_{{\mbox{\boldmath$\Phi$}}}<\frac{1}{3s}$. Thus, under these
conditions, with a probability of at least

$\Big{(}1-\frac{(n_{\theta}-s)}{n_{\theta}^{1+\nu}}\Big{)}\Big{(}1-e^{-\frac{s}{7}}\Big{)},\quad\forall\nu>0,$
(7)

the support set of the solution ${{\mbox{\boldmath$\theta$}}}^{\epsilon}_{1}$
obtained through $P_{1}^{\epsilon}$ is $\Lambda_{0}$.

According to Theorem 1, accurate sparse estimation of parameter vector is
achievable using $P_{1}^{\epsilon}$ if there is a small mutual coherence.
However, since there is no mutual coherence constraint in (7), its optimal
solution often does not meet the conditions of the above theorem. This is
evident from the optimal regressor in the numerical study of Section V.

To incorporate the need for low mutual coherence, we borrow the idea of
employing a linear coordinate transformation from [6], and write the model (1)
as

$\displaystyle\left\\{\begin{array}[]{lr}{\bf y}={\mbox{\boldmath$\Phi$}}{\bf
H}^{-1}{\bf x}+{\bf e}\\\ {\bf x}={\bf
H}{\mbox{\boldmath$\theta$}},\end{array}\right.$ (10)

in which the matrix ${\bf H}$ is utilized to transform the coordinates such
that the sparse parameter vector $\theta$ is estimated in the new coordinates
${\bf x}$. In (10) the vector ${\bf x}$ is not necessarily sparse. Hence, we
can use least squares to estimate this vector as follows:

${\mathbf{\hat{x}}}=\mathop{{\rm argmin}}_{{\bf x}}\|{\bf
y}-{\mbox{\boldmath$\Phi$}}{\bf H}^{-1}{\bf x}\|_{2}^{2}.$ (11)

Solving the above optimization problem leads to

$\displaystyle{\mathbf{\hat{x}}}={\bf
H}({\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}})^{-1}{\bf
H}^{T}({\mbox{\boldmath$\Phi$}}{\bf H}^{-1})^{T}{\bf y}={\bf
x}+{\mbox{\boldmath$\omega$}}={\bf
H}{\mbox{\boldmath$\theta$}}+{\mbox{\boldmath$\omega$}},$ (12)

where ${\mbox{\boldmath$\omega$}}={\bf
H}({\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}})^{-1}{\bf
H}^{T}({\mbox{\boldmath$\Phi$}}{\bf H}^{-1})^{T}{\bf e}$ and
$\operatorname{\mathbb{E}}[{\mbox{\boldmath$\omega$}}{\mbox{\boldmath$\omega$}}^{T}]=\sigma^{2}{\bf
H}({\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}})^{-1}{\bf H}^{T}$. In
the new linear regression model (12), the parameter vector $\theta$ is
estimated by using the regressor ${\bf H}$. Hence, this new regressor should
have a low mutual coherence for accurately estimating the parameters when the
vector $\theta$ is sparse. On the other hand, the coordinate transformation
has an undesired effect on the noise since the covariance matrix of $\omega$
in (12) is not an identity matrix. The following lemma provides an expression
for the mutual coherence of ${\bf H}$.

######  Lemma 2 ([6])

Consider the matrix ${\bf H}$ with columns normalized to unit length. The
mutual coherence of ${\bf H}$, denoted as $\mu_{{\bf H}}$, then is given by

$\mu_{{\bf H}}=\|{\bf H}^{T}{\bf H}-{\bf I}\|_{\infty}.$ (13)

In [6] where the purpose was to design ${\bf H}$ for facilitating sparse
estimation without input design, to find the optimal ${\bf H}$, the following
optimization problem was proposed:

$\displaystyle\min_{{\bf H}\in{\cal D}}\|{\bf
H}({\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}})^{-1}{\bf H}^{T}-{\bf
I}\|_{\infty}$ $\displaystyle\text{subject to}\leavevmode\nobreak\
\leavevmode\nobreak\ \|{\bf H}^{T}{\bf H}-{\bf I}\|_{\infty}\leq\alpha_{1},$
(14)

where ${\cal D}=\\{{\bf H}\in\mathbb{R}^{n_{\theta}\times n_{\theta}}:\|{\bf
H}(:,i)\|_{2}=1,\,\forall\leavevmode\nobreak\ 1\leq i\leq n_{\theta}\\}$. In
the above problem, the cost function is used to minimize the undesired effect
of the coordinate transformation on the noise. Simultaneously, the constraint,
as indicated by Lemma 13, guarantees that the newly transformed regressor
exhibits low mutual coherence provided $\alpha_{1}$ is chosen small enough.

In (14), the cost function minimizes the difference between ${\bf
H}({\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}})^{-1}{\bf H}^{T}$ and
the identity matrix. Using the Frobenius norm instead of the infinity norm
simplifies the process of solving the optimization problem. Consequently, here
we propose to instead use

$\displaystyle\min_{{\bf H}\in{\cal D}}\|{\bf
H}({\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}})^{-1}{\bf H}^{T}-{\bf
I}\|_{F}^{2}$ $\displaystyle\text{subject to}\leavevmode\nobreak\
\leavevmode\nobreak\ \|{\bf H}^{T}{\bf H}-{\bf I}\|_{\infty}\leq\alpha_{1}.$
(15)

Finally, to incorporate also the degrees of freedom offered by optimal input
design, we combine (6) with (15). This results in the following optimization
problem, where the coordinate transformation and input are jointly designed to
meet the objectives of low mutual coherence and achieving the optimal FIM:

$\displaystyle\min_{{\bf H}\in{\cal D},{\mbox{\boldmath$\Phi$}}\in{\cal P}}$
$\displaystyle\|{\bf
H}({\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}})^{-1}{\bf H}^{T}-{\bf
I}\|_{F}^{2}$ subject to $\displaystyle\|{\bf H}^{T}{\bf H}-{\bf
I}\|_{\infty}\leq\alpha_{1}$
$\displaystyle\|{\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}}-{\bf
I}_{F}^{d}\|_{F}^{2}\leq\alpha_{2}.$ (16)

In (6), the constraint ${\mbox{\boldmath$\Phi$}}\in{\cal D}$ was utilized to
enforce the structure of a linear FIR model. From $\Phi$ we can recover the
sequence $u^{\prime}(t)$, $t=1-n_{\theta},\ldots,N-1$. However, unless
$F(q)=1$, i.e. a FIR model is used, an additional step has to be employed to
obtain the input $u(t)$. From (3) it follows that

$u(t)=F^{-1}(q)u^{\prime}(t).$ (17)

Notice that in the case that $F(q)$ has zeros outside the unit circle, the
filtering needs to be non-causal to enforce stability.

##  IV Algorithm

The optimization problem (16) is non-convex and in this section, we will
outline an iterative numerical method for solving it. Each iteration consists
of three steps using alternating minimization where the first step is based on
the proximal mapping of the infinity norm. We begin by transforming (16) into
an unconstrained version by adding the constraints as penalties to the cost
function, resulting in the Lagrangian

$\displaystyle\min_{{\bf H}\in{\cal D},{\mbox{\boldmath$\Phi$}}\in{\cal P}}$
$\displaystyle\|{\bf
H}({\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}})^{-1}{\bf H}^{T}-{\bf
I}\|_{F}^{2}+\lambda\|{\bf H}^{T}{\bf H}-{\bf
I}\|_{\infty}+\lambda^{\prime}\|{\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}}-{\bf
I}_{F}^{d}\|_{F}^{2},$ (18)

where $\lambda$ and $\lambda^{\prime}$ denote positive penalty parameters. The
penalty parameters can be adjusted in an iterative procedure so that the
constraints in (16) are met. To solve the optimization problem (18) we define
the new variable ${\bf N}={\bf H}^{T}{\bf H}-{\bf I}$, which converts this
problem to the constrained problem

$\displaystyle\begin{array}[]{cl}\min\limits_{{\mbox{\boldmath$\Phi$}}\in{\cal
P},{\bf N}\in{\cal N},{\bf H}\in{\cal D}}&\|{\bf
H}({\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}})^{-1}{\bf H}^{T}-{\bf
I}\|_{F}^{2}+\lambda^{\prime}\|{\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}}-{\bf
I}_{F}^{d}\|_{F}^{2}+\lambda\|{\bf N}\|_{\infty}\\\ \text{subject to }&{\bf
N}={\bf H}^{T}{\bf H}-{\bf I},\end{array}$

in which ${\cal N}$ indicates that ${\bf N}$ is constrained to be a a
symmetric matrix for which ${\bf N}+{\bf I}>0$. We proceed and add also the
new constraint to the cost function, resulting in

$\displaystyle\min\limits_{{\mbox{\boldmath$\Phi$}}\in{\cal P},{\bf N}\in{\cal
N},{\bf H}\in{\cal D}}$ $\displaystyle\|{\bf
H}({\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}})^{-1}{\bf H}^{T}-{\bf
I}\|_{F}^{2}+\lambda^{\prime}\|{\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}}-{\bf
I}_{F}^{d}\|_{F}^{2}+\lambda\|{\bf N}\|_{\infty}+\frac{1}{2\varsigma}\|{\bf
N}+{\bf I}-{\bf H}^{T}{\bf H}\|_{F}^{2},$ (19)

where $\varsigma$ is a positive penalty parameter. In the above optimization
problem, the three matrices, ${\bf N}$, ${\bf H}$, and $\Phi$, act as decision
variables. Next, alternating minimization with three steps is used, where in
each step one of these three matrices is optimized, while the other variables
are kept fixed. This gives the following procedure:

1) Update ${\bf N}$:

${\bf N}_{k+1}^{c}=\mathop{{\rm argmin}}_{{\bf N}\in{\cal N}}\lambda\|{\bf
N}\|_{\infty}+\frac{1}{2\varsigma}\|{\bf N}+{\bf I}-{\bf H}_{k}^{T}{\bf
H}_{k}\|_{F}^{2}.$ (20)

2) Update ${\bf H}$:

$\displaystyle{\bf H}_{k+1}=\mathop{{\rm argmin}}_{{\bf H}\in{\cal D}}\,\|{\bf
H}({\mbox{\boldmath$\Phi$}}^{T}_{k}{\mbox{\boldmath$\Phi$}}_{k})^{-1}{\bf
H}^{T}-{\bf I}\|_{F}^{2}+\frac{1}{2\varsigma}\|{\bf N}_{k+1}+{\bf I}-{\bf
H}^{T}{\bf H}\|_{F}^{2}.$ (21)

3) Update $\Phi$:

$\displaystyle{\mbox{\boldmath$\Phi$}}_{k+1}=\mathop{{\rm
argmin}}_{{\mbox{\boldmath$\Phi$}}\in{\cal P}}\|{\bf
H}_{k+1}({\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}})^{-1}{\bf
H}^{T}_{k+1}-{\bf
I}\|_{F}^{2}+\lambda^{\prime}\|{\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}}-{\bf
I}_{F}^{d}\|_{F}^{2}.$ (22)

We use proximal mapping in the first step because there is a closed-form
solution for the proximal mapping of the infinity norm. We minimize the cost
function (20) without the constraint ${\cal N}$, and then the solution is
projected into ${\cal N}$.

######  Definition 1 ([20])

Let $p\colon\text{dom}_{p}\rightarrow(-\infty,+\infty]$ be a proper and lower
semicontinuous function. Then the proximal mapping of this function at ${\bf
x}\in{\mathbb{R}}^{n}$ is defined as

$\text{prox}_{p}({\bf x})=\mathop{{\rm
argmin}}_{{\mathbf{u}}\in\text{dom}_{p}}\left\\{\frac{1}{2}\|{\bf
x}-{\mathbf{u}}\|_{2}^{2}+p({\mathbf{u}})\right\\}.$ (23)

Since in the definition of proximal mapping, we deal with vectors, the
vectorization operator vec is applied to the first step of the above
alternating minimization, which gives

${\bf n}_{k+1}=\mathop{{\rm argmin}}_{{\bf n}}\lambda\|{\bf
n}\|_{\infty}+\frac{1}{2\varsigma}\|{\bf n}-{\bf f}_{k}\|_{2}^{2},$ (24)

in which ${\bf n}\triangleq\text{vec}({\bf N})$ and ${\bf
f}_{k}\triangleq\text{vec}({\bf F}_{k})$, with ${\bf F}_{k}\triangleq{\bf
H}_{k}^{T}{\bf H}_{k}-{\bf I}$. According to the definition of proximal
mapping in (23), the above problem is the proximal mapping of the function
$\lambda\leavevmode\nobreak\ \varsigma\|.\|_{\infty}$ at ${\bf f}_{k}$. Hence,
we can write

$\displaystyle{\bf n}_{k+1}=\text{prox}_{\lambda\leavevmode\nobreak\
\varsigma\|.\|_{\infty}}({\bf f}_{k}).$ (25)

Due to the presence of the infinity norm, a direct computation of this
proximal mapping may be difficult. However, the following lemma is of help.

######  Lemma 3 (Moreau decomposition [20])

For a proper and lower semicontinuous function
$p\colon{\mathbb{R}}^{m}\rightarrow(-\infty,+\infty]$ whose conjugate is
denoted by $p^{*}$, the following relation holds at ${\bf
h}\in{\mathbb{R}}^{m}$ and $\eta>0$:

${\bf h}=\text{prox}_{\eta p}({\bf
h})+\eta\,\text{prox}_{\frac{p^{*}}{\eta}}\left(\frac{{\bf h}}{\eta}\right).$
(26)

Assuming that $p=\|.\|_{\infty}$ and $\eta=\lambda\leavevmode\nobreak\
\varsigma$ in (26), and noting that the conjugate function of $\|.\|_{\infty}$
is the $\ell_{1}$-norm, the proximal mapping in (25) can be computed as

$\text{prox}_{\lambda\leavevmode\nobreak\ \varsigma\|.\|_{\infty}}({\bf
f}_{k})={\bf f}_{k}-\lambda\leavevmode\nobreak\
\varsigma\,\text{prox}_{\frac{\|.\|_{1}}{\lambda\leavevmode\nobreak\
\varsigma}}\left(\frac{{\bf f}_{k}}{\lambda\leavevmode\nobreak\
\varsigma}\right).$ (27)

Furthermore, it is straightforward to show that

$\text{prox}_{\frac{\|.\|_{1}}{\lambda\leavevmode\nobreak\
\varsigma}}\left(\frac{{\bf f}_{k}}{\lambda\leavevmode\nobreak\
\varsigma}\right)=\frac{1}{\lambda\leavevmode\nobreak\ \varsigma}{\cal
P}_{{\cal D}_{2}}({\bf f}_{k}),$ (28)

where ${\cal D}_{2}$ denotes the $\ell_{1}$-norm ball with radius
$\lambda\leavevmode\nobreak\ \varsigma$, and ${\cal P}_{{\cal D}_{2}}$ denotes
the orthogonal projection onto ${\cal D}_{2}$. Using (28) into (27) yields

$\text{prox}_{\lambda\leavevmode\nobreak\ \varsigma\|.\|_{\infty}}({\bf
f}_{k})={\bf f}_{k}-{\cal P}_{{\cal D}_{2}}({\bf f}_{k}).$ (29)

Next the inverse of the vectorization operator $\text{vec}^{-1}$ is applied to
(29), resulting in

${\bf N}_{k+1}=\text{vec}^{-1}(\text{prox}_{\lambda\leavevmode\nobreak\
\varsigma\|.\|_{\infty}}({\bf f}_{k}))={\bf F}_{k}-\text{vec}^{-1}({\cal
P}_{{\cal D}_{2}}({\bf f}_{k})).$ (30)

Next the above solution is projected into ${\cal N}$. Since in (30) the
projection ${\cal P}_{D_{2}}$ operates element-wise and the matrix ${\bf
F}_{k}$ is symmetric, ${\bf N}_{k+1}$ is symmetric. To enforce ${\bf N}+{\bf
I}>0$, we set

${\bf N}_{k+1}^{c}={\bf A}^{+}-{\bf I}$ (31)

where ${\bf N}_{k+1}+{\bf I}={\bf E}{\bf L}{\bf E}^{T}$ denotes the SVD
(Singular Value Decomposition) of ${\bf N}_{k+1}+{\bf I}$, ${\bf A}^{+}={\bf
E}{\bf L}^{+}{\bf E}^{T}$ and ${\bf L}^{+}=\max(\epsilon{\bf I},{\bf L})$ in
which the $\max$ operator takes the element-wise maximum between the diagonal
entries of two matrices and $\epsilon$ is a very small positive real value.
This gives the nearest positive definite matrix to ${\bf N}_{k+1}$ in the
Frobenius norm sense.

To summarize by using proximal mapping, a closed-form solution to update the
matrix ${\bf N}_{k+1}^{c}$ in the first step of cyclic minimization was found.

In the second and third steps of the alternating minimization, the cost
functions (21) and (22) are minimized with respect to ${\bf H}$ and $\Phi$,
respectively, while the matrix ${\bf N}_{k+1}$ is known from (30). Although
the main challenge in the optimization problem (19) with respect to ${\bf H}$
and $\Phi$ is that the problem is non-convex, the class of unconstrained
matrices when the cost function is zero satisfies (see [21])

${\bf H}({\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}})^{-1}{\bf
H}^{T}={\bf I},{\bf H}^{T}{\bf H}={\bf N}_{k+1}^{c}+{\bf
I},\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
{\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}}={\bf I}_{F}^{d}.$

The above equations are equivalent to the formulas

${\mbox{\boldmath$\Phi$}}={\mathbf{U}}{\bf H},{\bf H}={\bf Q}({\bf
N}_{k+1}^{c}+{\bf I})^{\frac{1}{2}},\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ {\mbox{\boldmath$\Phi$}}={\bf Z}({\bf
I}_{F}^{d})^{\frac{1}{2}},$ (32)

where ${\mathbf{U}}\in{\mathbb{R}}^{N\times n_{\theta}}$ and ${\bf
Z}\in{\mathbb{R}}^{N\times n_{\theta}}$ are semi-unitary matrices while ${\bf
Q}\in{\mathbb{R}}^{n_{\theta}\times n_{\theta}}$ is a unitary matrix. Thus,
the optimization problems (21) and (22) can be relaxed to

2) Update ${\bf H}$:

$\displaystyle({\bf H}_{k+1},{\bf Q}_{opt},{\mathbf{U}}_{opt})=\mathop{{\rm
argmin}}_{{\bf H}\in{\cal D},{\bf Q}\in{\cal Q},{\mathbf{U}}\in{\cal
U}}\,\|{\mathbf{U}}{\bf
H}-{\mbox{\boldmath$\Phi$}}_{k}\|_{F}^{2}+\frac{1}{\varsigma}\|{\bf H}-{\bf
Q}({\bf N}_{k+1}^{c}+{\bf I})^{\frac{1}{2}}\|_{F}^{2}.$ (33)

3) Update $\Phi$:

$\displaystyle({\mbox{\boldmath$\Phi$}}_{k+1},{\mathbf{U}}_{opt},{\bf
Z}_{opt})\quad=\mathop{{\rm argmin}}_{{\mbox{\boldmath$\Phi$}}\in{\cal
P},{\mathbf{U}}\in{\cal U},{\bf Z}\in{\cal
Z}}\|{\mbox{\boldmath$\Phi$}}-{\mathbf{U}}{\bf
H}_{k+1}\|_{F}^{2}+\lambda^{\prime}\|{\mbox{\boldmath$\Phi$}}-{\bf Z}({\bf
I}_{F}^{d})^{\frac{1}{2}}\|_{F}^{2},$ (34)

where ${\cal U}$, ${\cal Q}$ and ${\cal Z}$ are defined as ${\cal
U}=\\{{\mathbf{U}}\in{\mathbb{R}}^{N\times
n_{\theta}}\colon{\mathbf{U}}^{T}{\mathbf{U}}={\bf I}\\}$, ${\cal Q}=\\{{\bf
Q}\in{\mathbb{R}}^{n_{\theta}\times n_{\theta}}\colon{\bf Q}^{T}{\bf Q}={\bf
I}\\}$ and ${\cal Z}=\\{{\bf Z}\in{\mathbb{R}}^{N\times n_{\theta}}\colon{\bf
Z}^{T}{\bf Z}={\bf I}\\}$, respectively.

To solve the optimization problem (33), we have three variables. Hence, to
find these variables, the following cyclic minimization with two sub-steps is
proposed:

Step 2.1: The semi-unitary matrices ${\mathbf{U}}$ and ${\bf Q}$ are kept
fixed to their most recent values, and the problem (33) is solved with respect
to ${\bf H}$:

$\min_{{\bf H}\in{\cal D}}\|{\mathbf{U}}_{opt}{\bf
H}-{\mbox{\boldmath$\Phi$}}_{k}\|_{F}^{2}+\frac{1}{\varsigma}\|{\bf H}-{\bf
Q}_{opt}({\bf N}_{k+1}^{c}+{\bf I})^{\frac{1}{2}}\|_{F}^{2}.$

To solve the above minimization problem, the first derivative of this cost
function should be zero, which gives:

$\displaystyle{\bf
H}_{k+1}=\big{(}1+\frac{1}{\varsigma}\big{)}^{-1}\bigg{(}{\mathbf{U}}_{opt}^{T}{\mbox{\boldmath$\Phi$}}_{k}+\frac{1}{\varsigma}{\bf
Q}_{opt}({\bf N}_{k+1}^{c}+{\bf I})^{\frac{1}{2}}\bigg{)}.$ (35)

Due to the restriction to ${\cal D}$, each column of ${\bf H}_{k+1}$ should be
normalized, i.e.,

${\bf H}_{k+1}(:,i)=\frac{{\bf H}_{k+1}(:,i)}{\|{\bf
H}_{k+1}(:,i)\|_{2}},\;\forall\leavevmode\nobreak\ 1\leq i\leq n_{\theta}.$
(36)

Step 2.2: The optimal ${\mathbf{U}}$ and ${\bf Q}$ in (33) can be found
through SVD when ${\bf H}$ is known from (35), i.e.,

$\displaystyle{{\bf
H}_{k+1}{\mbox{\boldmath$\Phi$}}_{k}^{T}}=\bar{{\mathbf{U}}}{\mbox{\boldmath$\Sigma$}}{\tilde{{\mathbf{U}}}}^{T},{\mathbf{U}}_{opt}={\tilde{{\mathbf{U}}}}\bar{{\mathbf{U}}}^{T},$
$\displaystyle{{\bf H}_{k+1}}^{T}({\bf N}_{k+1}^{c}+{\bf
I})^{\frac{1}{2}}=\bar{{\mathbf{V}}}\tilde{{\mbox{\boldmath$\Sigma$}}}{\tilde{{\mathbf{V}}}}^{T},{\bf
Q}_{opt}=\bar{{\mathbf{V}}}{\tilde{{\mathbf{V}}}}^{T}.$ (37)

Finally, in the third step of the cyclic algorithm, i.e. the solution of (22),
the cost function (34) is minimized with respect to $\Phi$, ${\bf Z}$ and
${\mathbf{U}}$. Now, each Toeplitz matrix, e.g., $\Phi$, can be written as
${\mbox{\boldmath$\Phi$}}=\sum\limits_{l=-N+1}^{n_{\theta}-1}c_{l}{\bf
R}_{l}$, where the matrix ${\bf R}_{l}$ is defined by:

$\displaystyle{\bf
R}_{l}(i,j)=\left\\{\begin{array}[]{lr}1,\nobreak\leavevmode\nobreak\leavevmode&\text{if}\leavevmode\nobreak\
j=i+l,\leavevmode\nobreak\ \leavevmode\nobreak\ 1\leq i\leq N,\\\
&\max(1,1+l)\leq j\leq\min(N+l,n_{\theta})\\\
0,&\text{otherwise.}\end{array}\right.$

Also, it is straightforward to show that

$\text{vec}(\sum\limits_{l=-N+1}^{n_{\theta}-1}c_{l}{\bf R}_{l})={\bf W}{\bf
c},$ (38)

in which

${\bf W}=\begin{bmatrix}\text{vec}({\bf R}_{-N+1})&\text{vec}({\bf
R}_{-N+2})&\cdots&\text{vec}({\bf R}_{n_{\theta}-1})\end{bmatrix},$ ${\bf
c}=\begin{bmatrix}c_{-N+1}&c_{-N+2}&\cdots&c_{n_{\theta}-1}\end{bmatrix}^{T}.$

Substituting (38) in (34) results in

$\displaystyle({\bf c}_{k+1},{\mathbf{U}}_{opt},{\bf Z}_{opt})=$
$\displaystyle\mathop{{\rm argmin}}_{{\bf c},{\mathbf{U}}\in{\cal U},{\bf
Z}\in{\cal Z}}\|{\bf W}{\bf c}-\text{vec}({\bf Z}{\bf
T}^{\frac{1}{2}})\|_{2}^{2}+\lambda^{\prime}\|{\bf W}{\bf
c}-\text{vec}({\mathbf{U}}{\bf H}_{k+1})\|_{2}^{2}.$ (39)

To solve (39), we use cyclic minimization as follows:

3.1) Update ${\bf c}$:

${\bf c}_{k+1}=\mathop{{\rm argmin}}_{{\bf c}}\|{\bf W}{\bf c}-{\bf
z}_{opt}^{\prime}\|_{F}^{2}+\lambda^{\prime}\|{\bf W}{\bf
c}-{\mathbf{u}}_{opt}\|_{2}^{2},$ (40)

in which ${\bf z}^{\prime}_{opt}=\text{vec}({\bf Z}_{opt}{\bf
T}^{\frac{1}{2}})$ and ${\mathbf{u}}_{opt}=\text{vec}({\mathbf{U}}_{opt}{\bf
H}_{k+1})$. The first derivative of the above cost function should be zero.
Hence, we have:

${\bf c}_{k+1}=\frac{({\bf W}^{T}{\bf W})^{-1}}{1+\lambda^{\prime}}{\bf
W}^{T}({\bf z}^{\prime}_{opt}+\lambda^{\prime}{\mathbf{u}}_{opt}).$ (41)

Using the above equation, the regressor matrix $\Phi$ can be updated as:

$[c_{-N+1}^{k+1}\leavevmode\nobreak\ c_{-N+2}^{k+1}\leavevmode\nobreak\
\cdots\leavevmode\nobreak\ c_{n_{\theta}-1}^{k+1}]={\bf
c}_{k+1},{\mbox{\boldmath$\Phi$}}_{k+1}=\sum_{l=-N+1}^{n_{\theta}-1}c_{l}^{k+1}{\bf
R}_{l}.$ (42)

3.2) Update ${\bf Z}$ and ${\mathbf{U}}$:

In this step, we update the matrices ${\bf Z}_{opt}$ and ${\mathbf{U}}_{opt}$
through the SVD, i.e.:

$\displaystyle({\bf
I}_{F}^{d})^{\frac{1}{2}}{\mbox{\boldmath$\Phi$}}_{k+1}^{T}={\mathbf{S}}^{\prime}{\mbox{\boldmath$\Sigma$}}^{\prime}{{\mathbf{U}}^{\prime}}^{T}\rightarrow{\bf
Z}_{opt}={{\mathbf{U}}^{\prime}}{{\mathbf{S}}^{\prime}}^{T}$
$\displaystyle{{\bf
H}_{k+1}{\mbox{\boldmath$\Phi$}}_{k+1}^{T}}=\bar{{\mathbf{U}}}{\mbox{\boldmath$\Sigma$}}{\tilde{{\mathbf{U}}}}^{T}\rightarrow{\mathbf{U}}_{opt}={\tilde{{\mathbf{U}}}}\bar{{\mathbf{U}}}^{T}.$
(43)

The pseudo-code of the final algorithm is summarized in Algorithm 1.

Algorithm 1 Low Coherence Input Design (LCID)-OMP

Input: ${\bf I}_{F}^{d}$, ${\mbox{\boldmath$\Phi$}}_{1}$, ${\bf H}_{1}$,
$\varsigma$, $\lambda$, $\lambda^{\prime}$ and $0<c<1$

for $k=1$ to MaxIteration do

$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\varsigma=c\varsigma$

$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ {\bf
F}_{k}\triangleq{\bf H}_{k}^{T}{\bf H}_{k}-{\bf I},{\bf
f}_{k}\triangleq\text{vec}({\bf F}_{k})$

$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
1)\leavevmode\nobreak\ {\bf N}_{k+1}={\bf F}_{k}-\text{vec}^{-1}({\cal
P}_{{\cal D}_{2}}({\bf f}_{k}))$ and update ${\bf N}_{k+1}^{c}$ by

(31)

$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \textbf{for }i=1\text{ to }3\textbf{
do}$

$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
2.1)\leavevmode\nobreak\ \text{Update}\leavevmode\nobreak\ {\bf
H}_{k+1}\leavevmode\nobreak\ \text{by \eqref{eq:udateh} and
\eqref{eq:normali}}$

$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
2.2)\leavevmode\nobreak\ \text{Update}\leavevmode\nobreak\
{\mathbf{U}}_{opt}\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ {\bf
Q}_{opt}\leavevmode\nobreak\ \text{by \eqref{eq:updateuq}}$

end for

$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \textbf{for }i=1\text{ to }3\textbf{
do}$

$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
3.1)\leavevmode\nobreak\ \text{Update}\leavevmode\nobreak\
\leavevmode\nobreak\ {\mbox{\boldmath$\Phi$}}_{k+1}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{by}\leavevmode\nobreak\ \leavevmode\nobreak\
\eqref{eq:finphi}.$

$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
3.2)\leavevmode\nobreak\ \text{Update}\leavevmode\nobreak\ {\bf
Z}_{opt}\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
{\mathbf{U}}_{opt}\leavevmode\nobreak\ \text{by \eqref{eq:svdupdate}}$

end for

end for

Output: ${\mbox{\boldmath$\Phi$}}_{k+1}$ and ${\bf H}_{k+1}$.

$\leavevmode\nobreak\ \leavevmode\nobreak\
{\mathbf{\hat{x}}}=(({\mbox{\boldmath$\Phi$}}_{k+1}{\bf
H}_{k+1}^{-1})^{T}({\mbox{\boldmath$\Phi$}}_{k+1}{\bf
H}_{k+1}^{-1}))^{-1}({\mbox{\boldmath$\Phi$}}_{k+1}{\bf H}_{k+1}^{-1})^{T}{\bf
y}$

Support estimation:
${{\mbox{\boldmath$\theta$}}^{\prime}}=\text{OMP}({\mathbf{\hat{x}}},{\bf
H}_{k+1},s)$, ${\Lambda}=\mbox{{supp}}({\mbox{\boldmath$\theta$}}^{\prime})$.

Sparse parameter estimation:
${{\mbox{\boldmath$\theta$}}}=({{\mbox{\boldmath$\Phi$}}_{k+1}}_{\Lambda_{0}})^{\dagger}{\bf
y}$

###### Remark 4

We currently lack proof of convergence for the suggested approach,
particularly concerning Step 3. Nevertheless, the subsequent section’s
numerical simulations demonstrate favorable convergence outcomes; see also
[21] and the references therein for more examples. For further insights and
characteristics of alternating minimization, we also refer to [22, 23].
Additionally, it is worth noting that the proposed algorithm may also serve as
a means to determine initial values when addressing (16) with a general
purpose constrained nonlinear optimization algorithm.

##  V NUMERICAL STUDY

This section presents experimental evaluations of the proposed method in a
model reference control problem studied in [18]. The proposed method can be
utilized in conjunction with any sparse estimation algorithm, so that first
the regressors $\Phi$ and ${\bf H}$ are designed through Algorithm 1, and then
the sparse estimation algorithm is applied to the model (12). In this study,
we use the proposed method in conjunction with the OMP, LADMM [24] (Lasso via
Alternating Direction Method of Multipliers) and BIC [1], we denote these
combinations by LCID-OMP, LCID-LADMM, and LCID-BIC, respectively.

The simulations were performed on a laptop with a 4.00 GHz I7 CPU and 16GB
RAM, using MATLAB R2022a and Microsoft Windows 10.

The objective of model reference control problem is to design a feedback
controller to achieve a sensitivity function that is as close as possible to
the desired sensitivity $S$. The true system is assumed to be linear and time-
invariant and given by

$y(t)=G_{0}(q)u(t)+e(t),$ (44)

The true system has the structure
$G(q;{\mbox{\boldmath$\theta$}})=\sum\limits_{k=1}^{n_{\theta}}{\theta}_{k}q^{-k}$
such that $G(q,{\mbox{\boldmath$\theta$}}_{0})=G_{o}(q)$ and the exact order
$n_{\theta}$ is unknown. While this imposes more information than that the
parameter vector is sparse we will test if the proposed method is able to help
also with the model order selection problem. Two different values for
$n_{\theta}$ are considered, namely $n_{\theta}=40$ and $n_{\theta}=30$.
However, in both cases, only the first 10 coefficients of $\theta$ are non-
zero and the number of observations is $N=100$. Also, the noise ${e(t)}$ is
modeled as a realization of an independent zero-mean Gaussian random variable
with variance $\sigma^{2}$, which is chosen to obtain Signal-to-Noise Ratios
(SNRs) ranging from 5 to 29. The SNR is defined as
$\mbox{{SNR}}=10\log\frac{\|{\mbox{\boldmath$\Phi$}}{\mbox{\boldmath$\theta$}}_{0}\|_{2}^{2}}{\operatorname{\mathbb{E}}({{\bf
e}^{T}{\bf e}})}$.

A model reference controller when the given model $G$ is minimum phase, is as
follows:

$C(G)=\frac{1}{G}\frac{1-S}{S}.$ (45)

Here $S$ is the desired sensitivity function which we assume is given by

$S=\frac{1-q^{-1}}{1-\frac{1-\eta}{1+\eta}q^{-1}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \eta\in[0,1],$ (46)

where $\eta$ controls the desired closed-loop bandwidth. Using the reference
model controller (45) results in the following achieved sensitivity function

$S(G)=\frac{1}{1+G_{0}C(G)}.$ (47)

The study uses an application-oriented input design method to design the
desired FIM and input spectrum. Application-oriented input design requires the
performance degradation cost of the parameters estimate, which directly
represents the performance of the application, to be less than a pre-specified
bound $\frac{1}{2\gamma}$. One way to ensure this constraint is expressed in
[18] as

${\bf
I}_{F}\geq\gamma\chi_{\alpha}^{2}(n_{\theta})V^{{}^{\prime\prime}}_{app}({\mbox{\boldmath$\theta$}}_{0}),$
(48)

in which $V^{{}^{\prime\prime}}_{app}({\mbox{\boldmath$\theta$}}_{0})$ and
$\chi_{\alpha}^{2}(n_{\theta})$ denote the Hessian of the performance
degradation cost and $\alpha$-percentile of the chi-square distribution with
$n_{\theta}$ degrees of freedom, respectively.

In this study, the following performance degradation cost based on the $H_{2}$
error is considered

$V_{app}(G)=\frac{1}{2}\Big{\|}\frac{S(G)-S}{S}\Big{\|}_{2}^{2},$ (49)

so that the set of acceptable models is given by

$\varepsilon_{app}=\\{{\mbox{\boldmath$\theta$}}:V_{app}({\mbox{\boldmath$\theta$}})\leq\frac{1}{2\gamma}\\}.$
(47.b)

The proposed method is evaluated using two performance measures: the
performance degradation cost (49) and the normalized mean square error (NRMSE)
defined as
NRMSE$=\frac{\|{\mbox{\boldmath$\theta$}}_{0}-{{\mbox{\boldmath$\theta$}}}\|_{2}}{\|{\mbox{\boldmath$\theta$}}_{0}\|_{2}}$,
where $\theta$ denotes the estimated parameters.

Additionally, the input spectrum is taken to have a finite representation with
an order of $J$ given by

$\Phi_{u}(\omega)=\sum_{k=-J+1}^{J-1}r_{k}e^{j\omega k},$ (50)

where $\\{r_{k}\\}$ denotes the autocorrelation coefficients of the input
signal, i.e.,
$r_{k}=\operatorname{\mathbb{E}}\\{u(t)u(t-k)\\}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ \leavevmode\nobreak\
r_{k}=r_{-k}.$

The order of the input spectrum $J$ is set to 50 throughout this study.

With the experimental measured by the used input power, the resulting
application-oriented input design problem is given by

$\displaystyle\min_{{\Phi_{u}(\omega)}}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $
$\displaystyle\frac{1}{2\pi}\int_{-\pi}^{\pi}\Phi_{u}(\omega)\mathrm{d}\omega$
subject to $\displaystyle\Phi_{u}(\omega)=\sum_{k=-J+1}^{J-1}r_{k}e^{j\omega
k}\geq 0$ $\displaystyle{\bf
I}_{F}\geq\gamma\chi_{\alpha}^{2}(n_{\theta})V^{{}^{\prime\prime}}_{app}({\mbox{\boldmath$\theta$}}_{0}),$
(51)

where

$\displaystyle
V_{app}^{{}^{\prime\prime}}({\mbox{\boldmath$\theta$}})=\frac{1}{2\pi}\int_{-\pi}^{\pi}\Big{|}\frac{1-S(e^{j\omega})}{G_{0}(e^{j\omega})}\Big{|}^{2}G^{{}^{\prime}}(\theta_{0},e^{j\omega})(G^{{}^{\prime}}(\theta_{0},e^{j\omega}))^{*}d\omega.$

The above Hessian matrix is computed by numerical differentiation in MATLAB.

To solve the optimization problem (51), the MOOSE toolbox222For MOOSE2, we
used the MATLAB code available at
‘https://www.kth.se/polopoly_fs/1.195834.1600688635!/Menu/general/column-
content/attachment/moose_ver01f.zip’. [25] is utilized. The toolbox provides
the corresponding optimal FIM, denoted as ${\bf I}_{F}^{d}$. In the Frequency
Domain Method (FDM), the input is generated by filtering white noise through
the stable minimum phase spectral factor of the obtained input spectrum.

In this study, hyperparameters are used in all algorithms and are tuned using
Cross Validation (CV). The data is split into a training set (80$\%$ of the
data) and a validation set (20$\%$ of the data). Initially, the parameter
vector $\theta$ is estimated using the training set based on 200 different
hyperparameters that are logarithmically spaced between $10^{-5}$ and $10$.
Subsequently, the hyperparameters that result in the minimum error $\|{\bf
y}_{v}-{\mbox{\boldmath$\Phi$}}_{v}\hat{{\mbox{\boldmath$\theta$}}}\|_{2}$ are
selected using the validation set, where ${\bf y}_{v}$ and
${\mbox{\boldmath$\Phi$}}_{v}$ indicate the observations and regressor matrix
for the validation set, respectively. To obtain the minimum error, a
sufficiently small and large range of hyperparameters, i.e., $10^{-5}$ and
$10$, respectively, are chosen.  The hyperparameter of the algorithms are:
$\lambda$ (LADMM), $\lambda$ and $\lambda^{\prime}$ (proposed method), and
$\lambda$ (CORED-OMP as in [6]).

Other parameters of the proposed algorithms (see Algorithm 1) are set as
$\varsigma=1$, and $c_{1}=0.97$.

In the simulations, we compare the performance of the proposed methods, i.e.,
LCID-OMP, LCID-LADMM, and LCID-BIC with standard least-squares (LS) used in
conjunction with AICc (small sample-size corrected AIC) and BIC to select the
model order, resulting in LS-AICc and LS-BIC, respectively. Additionally, the
comparison includes OMP, LADMM, the method proposed in [13] referred to as
ARCID-OMP (Approximating the Reduction of Coherence in Input Design via OMP),
and CORED-OMP [6]. In LS-AICc, LS-BIC, OMP, and LADMM, the input is designed
using the FDM explained above. In CORED-OMP, the parameters of the system are
estimated using the coordinate transformation (14) in which the regressor
$\Phi$ is known from (6).

Additionally, we benchmark the performance with an oracle knowing the sparsity
pattern of ${\mbox{\boldmath$\theta$}}_{0}$. This method denoted FDM-Known
Sparsity, estimates the non-zero parameters using standard least-squares with
the input obtained from (51).

In order to assess the performance of the different methods, 400 experiments
are conducted per $\text{SNR}=[5\leavevmode\nobreak\ 8\leavevmode\nobreak\
11\leavevmode\nobreak\ 14\leavevmode\nobreak\ 17\leavevmode\nobreak\
20\leavevmode\nobreak\ 23\leavevmode\nobreak\ 26\leavevmode\nobreak\ 29]$. The
mean values of NRMSE and $V_{app}({\mbox{\boldmath$\theta$}})$ are then
computed for each method. It is worth noting that although the different
methods may result in varying input powers, the same SNRs are used for each
method for a fair comparison.

The mutual coherence of the regressor $\Phi$ in this study, denoted by
$\mu_{{\mbox{\boldmath$\Phi$}}}$, is 0.98 which is very high and fails to meet
the required upper bound in Theorem 1, whereas the mutual coherence of the
regressor in the new coordinates, represented by $\mu_{{\bf H}}$ for CORED-OMP
and LCID, is 0.05. Additionally, the obtained $\mu_{{\bf H}}$ for ARCID-OMP is
0.73. Therefore, this numerical study demonstrates that while the optimal
regressor $\Phi$ exhibits a high mutual coherence, we can obtain a small
mutual coherence through the proposed method by using the coordinate
transformation matrix ${\bf H}$.

Figure 2 displays the ideal input spectrum and sensitivity function for
${\mbox{\boldmath$\Phi$}}\in{\mathbb{R}}^{100\times 40}$, where $s=10$ and
$\eta=0.1$ obtained by solving (51). This figure indicates that the crucial
frequency interval of the optimal input spectrum, i.e., $[0.1,0.3]$rad/s,
coincides with the cross-over frequency interval of the desired sensitivity
function. Furthermore, the optimal input spectrum is negligible at other
frequencies. The FDM thus generates a highly correlated input sequence.

The results of performance degradation cost and NRMSE for
${\mbox{\boldmath$\Phi$}}\in{\mathbb{R}}^{100\times 40}$ and $s=10$ are
reported in Figures 2 and 4, respectively. In Figure 2, we see that FDM with
known sparsity performs the best. This is of course expected since the optimal
input design is used together with known sparsity pattern. However,
significant performance improvements are observed when the proposed method is
utilized with OMP, LADMM and BIC. As a result, LCID-OMP and LCID-BIC
outperform LS-AICc and LS-BIC significantly to reduce the model order, for low
SNRs. Turning to Figure 4, we see, perhaps surprisingly, that the three new
methods (LCID-OMP, LCID-LADMM, and LCID-BIC) have lower NRMSE than the oracle
method FDM-Known Sparsity. This is because in FDM-Known Sparsity, although the
sparsity pattern is known, the obtained optimal regressors are highly
correlated, making it difficult for an estimator to distinguish between
different regressors. On the other hand, the mutual coherence in LCID-OMP is
low, resulting in better estimation performance. We also see that the three
new methods outperform all other methods, including LS-AICc and LS-BIC.
Therefore, this study highlights the efficacy of the proposed method in
improving estimation accuracy and selecting model order.

Figures 4 and 5 depict the results of performance degradation cost and NRMSE,
respectively, for ${\mbox{\boldmath$\Phi$}}\in{\mathbb{R}}^{100\times 30}$ and
$s=10$. In this study, the LCID-OMP outperforms CORED-OMP significantly. A
noteworthy observation from these figures is that the proposed method yields
lower performance degradation cost and NRMSE compared to ARCID-OMP. This
difference can be attributed to the exact formulation of mutual coherence
within the proposed method, resulting in a calculated $\mu_{{\bf H}}$ of 0.05
for LCID, in contrast to the value of 0.73 obtained for ARCID. This numerical
analysis thus provides compelling evidence for the substantial improvement
achieved by LCID over ARCID.

According to the all figures, the proposed method has a significant effect on
both OMP and LADMM, so that using LCID greatly reduces the NRMSE and
$V_{app}({\mbox{\boldmath$\theta$}})$.

To measure the complexity of the methods, we report computational time which
is required to design the optimal regressor and estimate the parameters of
model for ${\mbox{\boldmath$\Phi$}}\in{\mathbb{R}}^{100\times 30}$, $s=10$ and
$\mbox{{SNR}}=15$ in Table I. As shown in this table, the computational
requirement is higher using LCID than the other methods, except for ARCID,
still the computational overhead is in no way dramatic. An intriguing
observation from this table is that utilizing the closed-form solution in
LCID, as opposed to employing the CVX toolbox in ARCID, leads to a significant
reduction in computational time.

The proposed method involves designing the regressor and coordinate
transformation jointly, whereas in CORED-OMP, the optimal regressor design is
independent of the coordinate transformation. Therefore, in LCID-OMP and
CORED-OMP, a regressor with low mutual coherence is obtained in the new
coordinates. However, compared to CORED-OMP, the coordinate transformation in
LCID-OMP has a lower undesired effect on the noise. This is why LCID-OMP
exhibits a better ability to estimate the sparse parameter vector than other
methods. Additionally, the error
$\frac{\|{\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}}-{\bf
I}_{F}^{d}\|_{F}}{\sqrt{Nn_{\theta}}}$ for LCID and other methods when
$n_{\theta}=40$ is 0.018 and 0.0001, respectively, which means the proposed
method sacrifices the error
$\frac{\|{\mbox{\boldmath$\Phi$}}^{T}{\mbox{\boldmath$\Phi$}}-{\bf
I}_{F}^{d}\|_{F}}{\sqrt{Nn_{\theta}}}$ slightly to achieve higher accuracy in
estimating the parameters.

Figure 1: Evaluation with assumption:
${\mbox{\boldmath$\Phi$}}\in{\mathbb{R}}^{100\times 40}$, $s=10$ and
$\zeta=0.1$.

Figure 2: Evaluation of performance degradation with assumptions:
${\mbox{\boldmath$\Phi$}}\in{\mathbb{R}}^{100\times 40}$ and $s=10$.

Figure 3: Evaluation of NRMSE with assumptions:
${\mbox{\boldmath$\Phi$}}\in{\mathbb{R}}^{100\times 40}$ and $s=10$.

Figure 4: Evaluation of performance degradation with assumptions:
${\mbox{\boldmath$\Phi$}}\in{\mathbb{R}}^{100\times 30}$ and $s=10$.

Figure 5: Evaluation of NRMSE with assumptions:
${\mbox{\boldmath$\Phi$}}\in{\mathbb{R}}^{100\times 30}$ and $s=10$.

Algorithm $\mbox{{SNR}}=15$ LCID-OMP 18.9s CORED-OMP 12.1s LS 9.1s FDM-Known
Sparsity 7.53s OMP 9.15s LADMM 9.12s LCID-LADMM 19s ACRID-OMP 548s

Table I: Average running times (in seconds) to design optimal input with
assumptions:${\mbox{\boldmath$\Phi$}}\in{\mathbb{R}}^{100\times 30}$ and
$s=10$.

##  VI Conclusion

The paper presents an algorithm for designing optimal inputs in models with a
sparse parameter vector. The method is applicable to various commonly used
models in system identification, such as FIR, Laugerre, Kautz, and generalized
basis function expansion models. The proposed method is based on minimizing
the difference between the Fisher information matrix and the desired Fisher
matrix, along with a coordinate transformation.

To compute the optimal input, we utilize cyclic minimization involving three
sub-steps and proximal mapping.

We have also applied the method to a model reference control problem where the
actual model order is unknown. Here, the proposed method, LCID, combined with
OMP, LADMM, and BIC, significantly improves the accuracy of estimating the
sparse parameter vector compared to other algorithms. When there is no prior
knowledge about the location of non-zero entries in the parameter vector,
LCID-OMP, and LCID-ADMM have lower performance degradation costs than other
methods. The numerical study also showed that model order selection criteria
such as AICc and BIC benefit significantly from the method for low SNRs.

The proposed method can be extended to incorporate input design for nonlinear
systems, where the regressor $\Phi$ includes the nonlinear components of input
sequences. It is noteworthy that many nonlinear models, including the Volterra
series and non-linear Auto-Regressive models with eXogenous inputs (NARX), are
inherently sparse. This important topic is under current research.

## References

  * [1] L. Ljung, _System Identification: Theory for the User_ , 2nd, Ed. Prentice-Hal, 1999.
  * [2] B. Wahlberg, “System identification using kautz models,” _IEEE Transactions on Automatic Control_ , vol. 39, no. 6, pp. 1276–1282, 1994.
  * [3] P. Heuberger, P.M.J. Van den Hof, and O. Bosgra, “A generalized orthonormal basis for linear dynamical systems,” _IEEE Transactions on Automatic Control_ , vol. 40, no. 3, pp. 451–465, March 1995.
  * [4] J. Kennington, E. Olinick, and D. Rajan, Eds., _Wireless Network Design_. Springer, June 2011, no. 978-1-4419-6111-2.
  * [5] A. P. Kannu and P. Schniter, “On communication over unknown sparse frequency-selective block-fading channels,” _IEEE Transactions on Information Theory_ , vol. 57, no. 10, pp. 6619–6632, 2011.
  * [6] J. Parsa, C. R. Rojas, and H. Hjalmarsson, “Transformation of regressors for low coherent sparse system identification,” _IEEE Transactions on Automatic Control_ , pp. 1–16, 2023.
  * [7] W. Pan, Y. Yuan, J. M. Gonçalves, and G.-B. V. Stan, “A sparse bayesian approach to the identification of nonlinear state-space systems,” _IEEE Transactions on Automatic Control_ , vol. 61, pp. 182–187, 2014.
  * [8] A. Chiuso and G. Pillonetto, “A bayesian approach to sparse dynamic network identification,” _Automatica_ , vol. 48, no. 8, pp. 1553–1565, 2012.
  * [9] M. Elad, _Sparse and Redundant Representations_. Springer, 2010.
  * [10] D. L. Donoho, M. Elad, and V. Temlyakov, “Stable recovery of sparse overcomplete representations in the presence of noise,” _IEEE Trans. Info. Theory_ , vol. 52, no. 1, pp. 6–18, 2006.
  * [11] Z. Ben-Haim, Y. C. Eldar, and M. Elad, “Coherence-based performance guarantees for estimating a sparse vector under random noise,” _IEEE Transactions on Signal Processing_ , vol. 58, no. 10, pp. 5030–5043, oct 2010\.
  * [12] J. Parsa, M. Sadeghi, M. Babaie-Zadeh, and C. Jutten, “Joint low mutual and average coherence dictionary learning,” in _2018 26th European Signal Processing Conference (EUSIPCO)_ , 2018, pp. 1725–1729.
  * [13] J. Parsa, C. R. Rojas, and H. Hjalmarsson, “Optimal input design for sparse system identification,” in _European Control Conference (ECC)_ , 2022, pp. 1999–2004.
  * [14] Y. C. Pati, R. Rezaiifar, and P. S. Krishnaprasad, “Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition,” in _In Proc. Asilomar Conf. Signal Syst. Comput._ , 1993\.
  * [15] J. Wang, S. Kwon, and B. Shim, “Generalized orthogonal matching pursuit,” _IEEE Trans. on Signal Proc._ , vol. 60, no. 12, pp. 6202–6216, 2012.
  * [16] G. Goodwin and R. Payne, _Dynamic System Identification: Experiment Design and Data Analysis_. Academic Press, New York, 1977.
  * [17] X. Bombois, G. Scorletti, M. Gevers, P. Van den Hof, and R. Hildebrand, “Least costly identification experiment for control,” _Automatica_ , pp. 1651–1662, 2006.
  * [18] H. Hjalmarsson, “System identification of complex and structured systems,” _European Journal of Control_ , pp. 275–310, 2009.
  * [19] A. Ebadat, B. Wahlberg, H. Hjalmarsson, C. Rojas, P. Hägg, and C. Larsson, “Applications oriented input design in time-domain through cyclic methods,” _IFAC Proceedings Volumes_ , vol. 47, no. 3, pp. 1422–1427, 2014, 19th IFAC World Congress.
  * [20] N. Parikh and S. Boyd, “Proximal algorithms,” _Foundations and Trends in Optimization_ , vol. 1, no. 3, pp. 123–231, 2014.
  * [21] P. Stoica, J. Li, X. Zhu, and B. Guo, “Waveform synthesis for diversity-based transmit beampattern design,” _IEEE Transactions on Signal Processing_ , vol. 56, pp. 2593–2598, 2008.
  * [22] J. A. Tropp, I. S. Dhillon, R. W. Heath, and T. Strohmer, “Designing structured tight frames via an alternating projection method,” _IEEE Trans. Info. Theory_ , vol. 51, no. 1, pp. 188–209, 2005.
  * [23] M. J. D. Powell, “On search directions for minimization algorithms,” _Mathematical Programming_ , vol. 4, no. 1, pp. 193–201, 1973.
  * [24] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Foundations and Trends in Machine Learning_ , vol. 3, no. 1, pp. 1–122, 2011.
  * [25] M. Annergren and C. Larsson, “MOOSE2—a toolbox for least-costly application-oriented input design,” _SoftwareX_ , vol. 5, pp. 96–100, 2016\.