# A Robust Approach to ARMA Factor Modeling

Lucia Falconi, Augusto Ferrante, Mattia Zorzi L. Falconi, A. Ferrante, and M.
Zorzi are with the Department of Information Engineering, University of
Padova, Padova, Italy; e-mail<EMAIL_ADDRESS>(L. Falconi);
<EMAIL_ADDRESS>(A. Ferrante<EMAIL_ADDRESS>(M. Zorzi).

###### Abstract

This paper deals with the dynamic factor analysis problem for an ARMA process.
To robustly estimate the number of factors, we construct a confidence region
centered in a finite sample estimate of the underlying model which contains
the true model with a prescribed probability. In this confidence region, the
problem, formulated as a rank minimization of a suitable spectral density, is
efficiently approximated via a trace norm convex relaxation. The latter is
addressed by resorting to the Lagrange duality theory, which allows to prove
the existence of solutions. Finally, a numerical algorithm to solve the dual
problem is presented. The effectiveness of the proposed estimator is assessed
through simulation studies both with synthetic and real data.

###### Index Terms:

Convex optimization, duality theory, dynamic factor analysis, nuclear norm.

##  I Introduction

We deal with the problem of constructing a dynamical model from a high
dimensional stream of data that are assumed to be noisy observations of a
process depending on a small number of hidden variables. In the static case,
this problem is known as factor analysis. Its origins can be traced back to
the beginning of the last century and the amount of literature produced on
this topic is impressive: we refer the readers to the recent papers [1, 2, 3,
4] for an overview of the literature and a rich list of references. The
solution of factor analysis problems may be obtained by decomposing the
covariance matrix of the observed data as the sum of a diagonal positive
definite matrix (accounting for the noise covariance) and a positive
semidefinite matrix whose rank must be as small as possible since it equals
the number of hidden variables in the model. The main problem of this solution
is that it is inherently fragile; in fact, even a minuscule variation in the
covariance matrix of the observed data usually leads to a substantial
variation of the number of hidden variables which is the key feature of the
modelling procedure. On the other hand such a matrix must be estimated and is
therefore subject to errors. To address this fragility issue, a robust method
has been recently proposed, [5], that has been generalized with good results
also to the dynamic framework of learning latent variable dynamic graphical
models, [6].

Dynamic Factor Analysis (DFA) has been addressed much more recently, the first
contribution in this field being apparently [7]. We refer to the surveys [8,
9] and to the recent paper [10] for an overview of the literature on this
subject. In [11] an interesting generalization is applied to modelization of
dynamical systems.

In this paper, we address the dynamic autoregressive moving average (ARMA)
case with the aim of extracting, from the observed data, a model featuring a
small number of hidden variables. This is important both from the point of
view of the model simplicity and to uncover the structure of the mechanism
generating the data. The problem may be mathematically formulated as that of
decomposing the spectral density of the process generating the data as the sum
of a diagonal spectral density and a low-rank one. The fragility issue in this
case is even more severe. We address the problem as follows:

  * •

Given the observed data, we compute by standard methods (e.g. truncated
periodogram) a raw estimate $\hat{\Phi}$ of the spectral density $\Phi$
generating the data.

  * •

We compute a neighbourhood $\cal N$ of $\hat{\Phi}$ that contains $\Phi$ with
prescribed probability; clearly the size of $\cal N$ depends on the sample
size.

  * •

We compute a refined estimate $\Phi^{\circ}\in{\cal N}$ by imposing that it
admits an additive decomposition as a diagonal spectral density and a spectral
density with the lowest possible rank. To this end we set up an optimization
problem that we address by resorting to duality. In particular, we prove
existence of solutions and provide a numerical algorithm to compute a
solution.

Our work may be cast in the rich stream of literature devoted to learning
dynamic models having a topological structure describing the presence or the
absence of interactions among the variables of the systems; see the former
works [12, 13, 14] as well as their extensions to reciprocal processes [15,
16], sparse plus low rank graphical models [17, 18, 6], the Bayesian viewpoint
proposed in [19, 20] and the case of oriented graphical models [21, 22].

The contribution of this paper is twofold; first, we propose a procedure to
estimate the number of latent factors in dynamic ARMA factor models: this is
the most delicate aspect of factor analysis problems; second, we derive an
identification method to estimate the parameters of a factor model describing
the observed data.

The outline of the paper is as follows. In Section II we introduce the DFA
problem for moving average (MA) models. Section III shows that such a problem
admits solution by means of duality theory, while Section IV shows how to
reconstruct the solution of the primal problem from the dual one. In Section V
we propose an algorithm to compute the solution of the dual problem. In
Section VI we extend the previous ideas to ARMA models. Section VII presents
some numerical results. Finally, in Section VIII we draw the conclusions.

###  I-A Notation

Given a matrix $M$, we denote its transpose by $M^{\top}$ and by $M_{(i,j)}$
the element of $M$ in the $i$-th row and $j$-th column. If $M$ is a square
matrix, $\operatorname{tr}(M)$, $|M|$ and $\sigma(M)$ denote its trace, its
determinant and its spectrum, respectively. The symbol $\|\cdot\|$ stands for
the Frobenius norm. For $A,B\in\mathbb{R}^{m\times m}$, we define their inner
product as $\langle A,B\rangle:=\operatorname{tr}(A^{\top}B)$. Let
$\mathbf{Q}_{m}$ be the space of real symmetric matrices of size $m$; if
$M\in\mathbf{Q}_{m}$ is positive definite or positive semidefinite, then we
write $M\succ 0$ or $M\succeq 0$, respectively. We denote by $(\cdot)^{*}$ the
complex conjugate transpose. $\Phi(e^{i\vartheta})$ for
$\vartheta\in[-\pi,\pi]\\}$ denotes a function defined on the unit circle
$\\{e^{i\vartheta}:\vartheta\in[-\pi,\pi]\\}$, and the dependence on
$\vartheta$ is dropped if necessary. If $\Phi(e^{i\vartheta})$ is positive
(semi-)definite $\forall\vartheta\in[-\pi,\pi]$ we write
$\Phi(e^{i\vartheta})\succ 0$ ( $\succeq 0$). Integrals are always defined
from $-\pi$ to $\pi$ with respect to the normalized Lebesgue measure
$d\vartheta/2\pi$.

##  II Identification of MA factor models

Consider the MA factor model whose order is $n$:

$\displaystyle y(t)=W_{L}u(t)+W_{D}w(t)$ (1)

where

$W_{L}(e^{i\vartheta})=\sum_{k=0}^{n}W_{L,k}e^{-i\vartheta k},\quad
W_{D}(e^{i\vartheta})=\sum_{k=0}^{n}W_{D,k}e^{-i\vartheta k},$

$W_{L,k}\in\mathbb{R}^{m\times r}$, $W_{D,k}\in\mathbb{R}^{m\times m}$
diagonal; $u=\\{u(t),\;t\in\mathbb{Z}\\}$ and $w=\\{w(t),\;t\in\mathbb{Z}\\}$
are normalized white Gaussian noises of dimension $r$ and $m$, respectively,
such that $\mathbb{E}[u(t)w(s)^{\top}]=0$ $\forall t,s.$ The aforementioned
model has the following interpretation: $u$ is the process which describes the
$r$ factors, with $r\ll m$, not accessible to observation; $W_{L}$ is the
factor loading transfer matrix; $W_{L}u(t)$ is the latent variable;
$W_{D}w(t)$ is idiosyncratic noise. Accordingly,
$y=\\{y(t),\;t\in\mathbb{Z}\\}$ is a $m$-dimensional Gaussian stationary
stochastic process with power spectral density

$\displaystyle\Phi=\Phi_{L}+\Phi_{D}$ (2)

where $\Phi_{L}=W_{L}W_{L}^{*}\succeq 0$ and $\Phi_{D}=W_{D}W_{D}^{*}\succeq
0$ belong to the finite dimensional space:

$\mathcal{Q}_{m,n}=\left\\{\sum_{k=-n}^{n}R_{k}e^{-i\vartheta
k},\;\;R_{k}=R_{-k}^{T}\in\mathbb{R}^{m\times m}\right\\}.$

By construction, $\mathrm{rank}(\Phi_{L})=r$, where $\mathrm{rank}$ denotes
the normal rank (i.e. the rank almost everywhere), and $\Phi_{D}$ is diagonal.
Therefore, $y$ represents a factor model if its spectral density can be
decomposed as “low rank plus diagonal” as in (2).

Assume to collect a finite length realization of $y$ defined in (1), say
$\mathrm{y}^{N}=\\{\,\mathrm{y}(1)\ldots\mathrm{y}(N)\,\\}$ where the order
$n$ is known. We want to estimate the corresponding factor model, that is the
decomposition in (2) as well as the number of factors $r$. To this aim, given
our data $\mathrm{y}^{N}$, we first compute the sample covariance lags
$\hat{R}_{j}$ as

$\hat{R}_{j}=\dfrac{1}{N}\sum_{t=0}^{N-j}\mathrm{y}(t+j)\mathrm{y}(t)^{\top},\;\;\;j=0\ldots
n.$

Then, an estimate $\hat{\Phi}$ of $\Phi$ is obtained by the truncated
periodogram:

$\displaystyle\hat{\Phi}=\sum_{k=-n}^{n}\hat{R}_{k}e^{i\vartheta k}.$ (3)

Notice that $\hat{\Phi}$ could be not positive definite for all $\vartheta$;
in that case, we can add $\varepsilon I_{m}$ to the right side of Equation
(3), with the constant $\varepsilon>0$ chosen in such a way as to ensure the
positivity of $\hat{\Phi}$. On the other hand, $\hat{\Phi}$ may not admit a
low rank plus diagonal decomposition. Thus, we estimate directly the two terms
$\Phi_{L}$ and $\Phi_{D}$ of the decomposition (2) by solving the following
optimization problem:

$\displaystyle\min_{\Phi,\Phi_{L},\Phi_{D}\in\mathcal{Q}_{m,n}}$
$\displaystyle\operatorname{tr}{\int}\Phi_{L}$ (4) subject to
$\displaystyle\Phi_{L}+\Phi_{D}=\Phi,$ $\displaystyle\Phi\succ 0\text{
a.e.},\;\;\Phi_{L},\Phi_{D}\succeq 0,$ $\displaystyle\Phi_{D}\hbox{
diagonal},$ $\displaystyle\mathcal{S}_{IS}(\Phi||\hat{\Phi})\leq\delta.$

Here, the objective function promotes a solution for $\Phi_{L}$ having low
rank, see [17]. The first three constraints impose that $\Phi_{L}$ and
$\Phi_{D}$ provide a genuine spectral density decomposition of type (2). The
last constraint, in which $\mathcal{S}_{IS}(\Phi||\hat{\Phi})$ is the Itakura-
Saito divergence defined by

$\mathcal{S}_{IS}(\Phi||\hat{\Phi})=\int\log|\hat{\Phi}{\Phi}^{-1}|+\operatorname{tr}[\hat{\Phi}^{-1}\Phi-
I_{m}],$

imposes that $\Phi$ belongs to a set “centered” in the nominal spectral
density $\hat{\Phi}$ and with prescribed tolerance $\delta$. Notice that
$\Phi_{D}$ is uniquely determined by $\Phi$ and $\Phi_{L}$. Thus, Problem (4)
can be rewritten by removing $\Phi_{D}$:

$\displaystyle(\Phi^{\circ},\Phi_{L}^{\circ})=\operatorname*{arg\,min}_{\Phi,\Phi_{L}\in\mathcal{Q}_{m,n}}$
$\displaystyle\operatorname{tr}{\int}\Phi_{L}$ (5) subject to
$\displaystyle\Phi\succ 0\text{ a.e.},\;\;\Phi_{L},\Phi-\Phi_{L}\succeq 0,$
$\displaystyle\Phi-\Phi_{L}\hbox{ diagonal},$
$\displaystyle\mathcal{S}_{IS}(\Phi||\hat{\Phi})\leq\delta.$

###  II-A The Choice of $\delta$

Before solving our problem, we deal with the choice of the tolerance parameter
$\delta$ appearing in the constraint of (5). This choice should reflect the
accuracy of the estimate $\hat{\Phi}$ of $\Phi$. This can be accomplished by
choosing a desired probability $\alpha\in(0,1)$ and considering a ball of
radius $\delta_{\alpha}$ (in the Itakura-Saito topology) centered in
$\hat{\Phi}$ and containing the true spectrum $\Phi$ with probability
$\alpha$. The estimation of $\delta_{\alpha}$ is not an easy task because we
do not know the true power spectral density $\Phi$. Next, we propose a
resampling-based method to estimate it.
The idea is to approximate $\Phi$ with $\hat{\Phi}$, and use this model to
perform a resampling operation. Let

$W(e^{i\vartheta})=\sum_{k=0}^{n}W_{k}e^{-i\vartheta
k},\;W_{k}\in\mathbb{R}^{m\times m}$

be the minimum phase spectral factor of $\hat{\Phi}$ and define the process
$\hat{y}=\\{\hat{y}(t),\;t\in\mathbb{Z}\\}$ as
$\hat{y}(t):=W(e^{i\vartheta})e(t),$ where $e(t)$ is an $m$-dimensional
normalized white noise. The truncated periodogram (understood as estimator)
based on a sample of the process $\hat{y}$ of length $N$ is

$\hat{\mathbf{\Phi}}_{r}(e^{i\vartheta})=\sum_{k=-n}^{n}e^{-i\vartheta
k}\frac{1}{N}\sum_{t=0}^{N-k}\hat{y}(t+k)\hat{y}(t)^{T},$

where the subscript “$r$” stands for resampling, as it is the means by which
we perform the resampling operation. By generating a realization
$\hat{\mathrm{y}}^{N}=\\{\,\hat{\mathrm{y}}(1)\ldots\hat{\mathrm{y}}(N)\,\\}$
from $\hat{\Phi}$ (i.e. by resampling the data), we can easily obtain a
realization of the random variable
$\mathcal{S}_{IS}(\hat{\Phi}||\hat{\mathbf{\Phi}}_{r})$. Accordingly, it is
possible to compute numerically $\delta_{\alpha}$ such that
$\text{Pr}(\mathcal{S}_{IS}(\hat{\Phi}||\hat{\mathbf{\Phi}}_{r})\leq\delta_{\alpha})=\alpha$
by a standard Monte Carlo procedure. Numerical simulations show that this
technique indeed provides a good estimate of $\delta$.

It is worth noting that if the chosen $\alpha$ is too large with respect to
the data length $N$, the resulting $\delta_{\alpha}$ may be too generous
yielding to a diagonal $\Phi$ obeying
$\mathcal{S}_{IS}(\Phi||\hat{\Phi})\leq\delta_{\alpha}$. In this case Problem
(5) admits the trivial solution $\Phi_{L}=0$ and $\Phi_{D}=\Phi$. To rule out
this trivial case, $\delta$ in (5) must be be strictly smaller than the upper
bound

$\delta_{\max}:=\min_{\begin{subarray}{c}\Phi\in\mathcal{S}_{m}^{+}\\\
\Phi\hbox{ diagonal}\end{subarray}}\mathcal{S}_{IS}(\Phi||\hat{\Phi})$

where $\mathcal{S}_{m}^{+}$ denotes the family of bounded and coercive
functions defined on the unit circle and taking values in the cone of positive
definite $m\times m$ Hermitian matrices. Since $\Phi$ must be diagonal, by
denoting with $\phi_{i}$ and by $\hat{\gamma}_{i}$ the $i$-th element in the
diagonal of $\Phi$ and of $\hat{\Phi}^{-1}$, respectively, we have

$\displaystyle\delta_{\max}$
$\displaystyle=\left[\sum_{i=1}^{m}\min_{\phi_{i}\in\mathcal{S}^{+}_{1}}\mathcal{S}_{IS}({\phi}_{i}||\hat{\gamma}_{i}^{-1})\right]+\int\log|\hat{\Phi}\operatorname{diag}^{2}(\hat{\Phi}^{-1})|$

where $\operatorname{diag}^{2}(\cdot)$ is the (orthogonal projection) operator
mapping a square matrix $M$ into a diagonal matrix of the same size having the
same main diagonal of $M$. Therefore, since the Itakura-Saito divergence is
nonnegative, the solution corresponds to
$\phi^{opt}_{i}(e^{i\vartheta})=(\hat{\gamma}_{i}(e^{i\vartheta}))^{-1}$,
$i=1,...,m$ for which
$\mathcal{S}_{IS}(\phi_{i}^{opt}||\hat{\gamma}_{i}^{-1})=0$. Accordingly,

$\delta_{max}=\int\log|\hat{\Phi}\operatorname{diag}^{2}(\hat{\Phi}^{-1})|.$
(6)

The derivation of the aforementioned result is based on reasonings similar to
[6, Section IV].

A more generous upper bound can be derived by assuming that $\Phi$ is the
spectrum of an MA process of order $n$. However, numerical experiments showed
that $\delta_{max}\gg\delta_{\alpha}$ even in the case that $N$ is relatively
small.

##  III Problem solution

In this section we first provide a finite dimensional matrix parametrization
of Problem (5). The latter is then analyzed by resorting to the Lagrange
duality theory, which allows us to prove the existence of a solution.

###  III-A Matricial Reparametrization of the Problem

To study Problem (5) it is convenient to introduce the following matrix
parametrization for $\Phi,\Phi_{L}$ and $\Phi-\Phi_{L}$:

$\displaystyle\Phi=\Delta X\Delta^{*}$ $\displaystyle\in\mathcal{Q}_{m,n}$ (7)
$\displaystyle\Phi_{L}=\Delta L\Delta^{*}$ $\displaystyle\in\mathcal{Q}_{m,n}$
$\displaystyle\Phi-\Phi_{L}=\Delta(X-L)\Delta^{*}$
$\displaystyle\in\mathcal{Q}_{m,n}$

where $\Delta(e^{i\theta})$ is the so-called shift operator:

$\Delta(e^{i\vartheta}):=[I_{m}\quad e^{i\vartheta}I_{m}\quad\dots\quad
e^{in\vartheta}I_{m}];$ (8)

$X$ and $L$ are matrices in $\mathbf{Q}_{m(n+1)}$ and $X_{ij}$ denotes the
block of $X$ in position $i,j$ with $i,j=0,\dots,n$, so that

$X=\begin{bmatrix}X_{00}&X_{01}&\dots&X_{0n}\\\
X_{01}^{\top}&X_{11}&\dots&\vdots\\\ \vdots&\vdots&\vdots&\vdots\\\
X_{0n}^{\top}&X_{1n}^{\top}&\dots&X_{nn}\end{bmatrix}.$

Moreover, $\textbf{M}_{m,n}$ denotes the vector space of matrices of the form

$Y:=[Y_{0}\quad Y_{1}\quad...\quad Y_{n}],\quad Y_{0}\in\mathbf{Q}_{m},\ \
Y_{1},...,Y_{n}\in\mathbb{R}^{m\times m}.$ (9)

The linear mapping $T:\mathbf{M}_{m,n}\rightarrow\mathbf{Q}_{m(n+1)}$
constructs a symmetric block-Toeplitz matrix from its first block row so that
if $Y$ is given by (9),

$T(Y)=\begin{bmatrix}Y_{0}&Y_{1}&\dots&Y_{n}\\\
Y_{1}^{\top}&Y_{0}&\ddots&\vdots\\\ \vdots&\ddots&\ddots&Y_{1}\\\
Y_{n}^{\top}&\dots&Y_{1}^{\top}&Y_{0}\end{bmatrix}.$

The adjoint of $T$ is the mapping
$D:\mathbf{Q}_{m(n+1)}\rightarrow\mathbf{M}_{m,n}$ defined by
$D(X)=[[D(X)]_{0}\quad\dots\quad[D(X)]_{n}]$ with

$[D(X)]_{0}=\sum_{h=0}^{n}X_{hh},\quad[D(X)]_{j}=2\sum_{h=0}^{n-j}X_{hh+j},\
j=1,...,n.$

Next, the objective is to provide a more convenient formulation of Problem (5)
in terms of $X$ and $L$. To this end, we have to take into account the
following points.

1) Positivity Constraints $\Phi\succ 0\text{ a.e.}$ and
$\Phi_{L},\Phi-\Phi_{D}\succeq 0:$ It can been shown (see, for example, [17,
Appendix A]) that, for any $\Psi\in\mathcal{Q}_{m,n}$, $\Psi\succeq 0$ if and
only if there exists a matrix $P\in\mathbf{Q}_{m(n+1)}$ such that $\Delta
P\Delta^{*}$ and $P\succeq 0$. Therefore, we replace the conditions
$\Phi_{L}\succeq 0$ with $L\succeq 0$, the condition $\Phi-\Phi_{L}\succeq 0$
with $X-L\succeq 0$. Note that these conditions only guarantees $X\succeq 0$
and thus $\Phi$ to be positive semidefinite, however we will show that this is
sufficient to guarantee that $\Phi\succ 0\text{ a.e.}$ at the optimum.

2) Constraint $\Phi-\Phi_{L}$ diagonal: Let
$\operatorname{ofd}\,:\,\mathbb{R}^{m\times m}\rightarrow\mathbb{R}^{m\times
m}$ denote the linear operator such that, given $A\in\mathbb{R}^{m\times m}$,
$\operatorname{ofd}(A)$ is the matrix in which each off-diagonal element is
equal to the corresponding element of $A$ and each diagonal element is zero.
We define the “block ofd” linear operator
$\operatorname{ofd}_{B}\,:\,\mathbf{M}_{m,n}\rightarrow\mathbf{M}_{m,n}$ as
follows. Given $Z=[\,Z_{0}\;Z_{1}\ldots Z_{n}\,]\in\mathbf{M}_{m,n}$, then

$\operatorname{ofd}_{B}(Z)=[\,\operatorname{ofd}(Z_{0})\;\operatorname{ofd}(Z_{1})\ldots\operatorname{ofd}(Z_{n})\,].$

It is not difficult that $\operatorname{ofd}_{B}$ is a self-adjoint operator,
since $\operatorname{ofd}$ is self-adjoint as well. Then, it is easy to see
that the condition $\Phi-\Phi_{L}$ diagonal is equivalent to the condition
$[D(X-L)]_{j}$ diagonal for $j=0,\dots,n$, that is
$\operatorname{ofd}_{B}(D(X-L))=0$.

3) The Low Rank Regularizer: We have

$\displaystyle\operatorname{tr}\int\Phi_{L}$
$\displaystyle=\operatorname{tr}\int\Delta
L\Delta^{*}=\operatorname{tr}\left(L\int\Delta^{*}\Delta\right)=\operatorname{tr}(L)$

where we exploited the fact that $\int e^{ij\vartheta}=1$ if $j=0$, and $\int
e^{ij\vartheta}=0$ otherwise.

4) The Divergence Constraint: A convenient matrix parameterization of the
Itakura-Saito divergence $\mathcal{S}_{IS}(\Phi||\hat{\Phi})$ can be obtained
by making use of the following facts.

First, since $\Phi=\Delta X\Delta^{*}$ with $X\succeq 0$, there exists
$A\in\mathbb{R}^{m\times m(n+1)}$ such that $X=A^{\top}A$. Then, by using the
Jensen-Kolmogorov formula we obtain

$\int\log|\Phi|=\\!\\!\int\log|\Delta
A^{\top}A\Delta^{*}|=\log|A_{0}^{\top}A_{0}|=\log|X_{00}|$ (10)

which holds provided that $X_{00}\succ 0$ and $\Phi$ is coercive (i.e.
$|\Phi|$ is bounded away from zero on the unit circle). We need to generalize
this result to spectral densities that may be singular on the unit circle.
This is possible because the zeros of a rational spectral density, if any,
have finite multiplicity so that the logarithm of the determinant of a
rational spectral $\Phi$ is integrable as long as the normal rank of $\Phi$ is
full.

###### Lemma III.1

Consider a power spectral density $\Phi\in\mathcal{Q}_{m,n}$ having full
normal rank. Let $X\in\mathbf{Q}_{m(n+1)}$ be such that $X\succeq 0$,
$X_{00}\succ 0$, and $\Phi=\Delta X\Delta^{*}$. Then

$\int\log|\Phi|=\log|X_{00}|.$

The proof is deferred to the appendix.

A second observation in order to conveniently parameterize the Itakura-Saito
divergence constraint is that, by exploiting the cyclic property of the trace,

$\displaystyle\int\operatorname{tr}(\hat{\Phi}^{-1}\Phi)$
$\displaystyle=\int\operatorname{tr}(\hat{\Phi}^{-1}\Delta X\Delta^{*})$
$\displaystyle=\operatorname{tr}\left(X\int\Delta^{*}\hat{\Phi}^{-1}\Delta\right)=\langle
X,T(\hat{P})\rangle,$

where $\hat{P}$ is defined from the expansion

$\hat{\Phi}^{-1}=\sum_{k=-\infty}^{\infty}\hat{P}_{k}e^{-i\vartheta k}$

as $\hat{P}:=[\hat{P}_{0}\dots\hat{P}_{n}].$

Summing up, we get the following matrix re-parametrization of Problem (5):

$\displaystyle(X^{\circ},L^{\circ})=\operatorname*{arg\,min}_{X,L\in\mathbf{Q}_{m(n+1)}}$
$\displaystyle\operatorname{tr}(L)$ (11) subject to $\displaystyle X_{00}\succ
0,\,L\succeq 0,\,X-L\succeq 0,$
$\displaystyle\operatorname{ofd}_{B}(D[X-L])=0,$
$\displaystyle-\log|X_{00}|+\int\log|\hat{\Phi}|$ $\displaystyle+\langle
X,T(\hat{P})\rangle-m\leq\delta.$

We remark once again that to prove the equivalence between (5) and (11) we
still need to show that $\Phi\succ 0\text{ a.e.}$ at the optimum: this fact
will be established after the variational analysis.

###  III-B The Dual Problem

We reformulate the constrained minimization problem in (11) as an
unconstrained problem by means of Duality Theory.
If we use $V,U\in\mathbf{Q}_{m(n+1)},\;V,U\succeq 0$ as the multipliers
associated with the constraints on the positive semi-definiteness of $X-L$ and
$L$, respectively; $Z\in\mathbf{M}_{m,n}$ as the multiplier associated with
the constraint $\operatorname{ofd}_{B}(D(X-L))=0$ and
$\lambda\in\mathbb{R},\lambda\geq 0$, as the multiplier associated with the
Itakura-Saito divergence, then the Lagrangian of Problem (11) is

$\displaystyle\mathcal{L}(X,L,\lambda,$ $\displaystyle
U,V,Z)=\operatorname{tr}(L)-\langle V,X-L\rangle-\langle U,L\rangle+$ (12)
$\displaystyle\langle
Z,\operatorname{ofd}_{B}(D(X-L))\rangle+\lambda\big{(}-\log|X_{00}|+$
$\displaystyle\int\log|\hat{\Phi}|+\langle
X,T(\hat{P})\rangle-m-\delta\big{)}\ $ $\displaystyle=\langle
L,I-U+V-T(\operatorname{ofd}_{B}(Z))\rangle+$ $\displaystyle\langle
X,T(\operatorname{ofd}_{B}(Z))-V+\lambda T(\hat{P})\rangle-$
$\displaystyle\lambda\big{(}\log|X_{00}|-\int\log|\hat{\Phi}|+m+\delta\big{)}.$

Note that we have not included the constraint $X_{00}\succ 0$ because, as we
will show later on, this condition is automatically met by the solution of the
dual problem.

The dual function is defined as the infimum of $\mathcal{L}$ over $X$ and $L$.
Thanks to the convexity of the Lagrangian, we rely on standard variational
methods to characterize the minimum.

  * •

Partial minimization with respect to $L$: $\mathcal{L}$ depends on $L$ only
through $\langle L,I-U+V-T(\operatorname{ofd}_{B}(Z))\rangle$ which is bounded
below only if

$I-U+V-T(\operatorname{ofd}_{B}(Z))=0.$ (13)

Thus, we get that

$\displaystyle\inf_{L}\mathcal{L}=\begin{cases}\langle
X,T(\operatorname{ofd}_{B}(Z))-V+\lambda T(\hat{P})\rangle-\\\
\lambda\big{(}\log|X_{00}|-\int\log|\hat{\Phi}|+m+\delta\big{)}&\text{if
\eqref{over_S}}\\\ -\infty&\text{otherwise.}\end{cases}$

  * •

Partial minimization with respect to $X$: The terms in $X_{00}$ are bounded
below only if

$\left[T(\operatorname{ofd}_{B}(Z))-V+\lambda T(\hat{P})\right]_{00}\succ 0$
(14)

and are minimized if $\lambda>0$ and

$X_{00}=\left(\big{[}T(\hat{P})+\lambda^{-1}(T(\operatorname{ofd}_{B}(Z))-V)\big{]}_{00}\right)^{-1}.$
(15)

The Lagrangian is linear in the remaining variables $X_{lh}$, for
$(l,h)\neq(0,0)$, and therefore bounded below only if

$\left[T(\operatorname{ofd}_{B}(Z))-V+\lambda
T(\hat{P})\right]_{lh}=0\quad\forall(l,h)\neq(0,0).$ (16)

Therefore, the minimization of the Lagrangian with respect to $X$ and $L$ is
finite if and only if (13), (14), and (16) hold in which case

$\displaystyle\min_{X,L}\mathcal{L}=-\lambda\big{(}-\log\big{|}\big{[}T(\hat{P})+\lambda^{-1}(T(\operatorname{ofd}_{B}(Z))$
$\displaystyle-V)\big{]}_{00}\big{|}-\int\log\big{|}\hat{\Phi}\big{|}+\delta\big{)}.$

Otherwise the Lagrangian has no minimum and its infimum is $-\infty$.

To simplify the notation, let us define the vector space $\mathcal{O}$ as:

$\displaystyle\mathcal{O}:=\\{Z\in\mathbf{M}_{m,n}:\operatorname{ofd}_{B}(Z)=Z,\;j=0,...,n\\};$

since $Z$ always appears in the form $\operatorname{ofd}_{B}(Z)$, we can
replace it with $Z\in\mathcal{O}$. Then, we can formulate the dual problem for
the Lagrangian (12) as

$\underset{(\lambda,U,V,Z)\in\tilde{\mathcal{C}}}{\max}\tilde{J}$ (17)

where

$\tilde{J}:=\lambda\Big{(}\log\big{|}\big{[}T(\hat{P})+\lambda^{-1}(T(Z)-V)\big{]}_{00}\big{|}+\int\log|\hat{\Phi}|-\delta\Big{)}$

and the feasible set $\tilde{\mathcal{C}}$ is given by:

$\displaystyle\tilde{\mathcal{C}}$
$\displaystyle:=\\{(\lambda,U,V,Z):U,V\in\mathbf{Q}_{m(n+1)},U,V\succeq
0,Z\in\mathcal{O},$
$\displaystyle\lambda\in\mathbb{R},\lambda>0,I-U+V-T(Z)=0,[\lambda
T(\hat{P})+T(Z)-$ $\displaystyle V]_{00}\succ 0,[\lambda
T(\hat{P})+T(Z)-V]_{lh}=0\;\;\forall(l,h)\neq(0,0)\\}.$

Note that the constraints $I-U+V-T(Z)=0$ and $U\succeq 0$ are equivalent to
the constraint $I+V-T(Z)\succeq 0$. Thus, we can eliminate the redundant
variable $U$; moreover, by changing the sign to the objective function
$\tilde{J}$ and observing that
$\big{[}T(\hat{P})+\lambda^{-1}(T(Z)-V)]_{00}=\hat{P}_{0}+\lambda^{-1}(Z_{0}-V_{00})$,
we can rewrite (17) as a minimization problem:

$\min_{(\lambda,V,Z)\in\mathcal{C}}J$ (18)

where

$J:=\lambda\Big{(}-\log\big{|}\hat{P}_{0}+\lambda^{-1}(Z_{0}-V_{00})\big{|}-\int\log|\hat{\Phi}|+\delta\Big{)}.$

and the corresponding feasible set $\mathcal{C}$ is:

$\displaystyle\mathcal{C}$
$\displaystyle:=\\{(\lambda,V,Z):V\in\mathbf{Q}_{m(n+1)},V\succeq
0,Z\in\mathcal{O},$ $\displaystyle I+V-T(Z)\succeq
0,\lambda\in\mathbb{R},\lambda>0,[\lambda\hat{P}_{0}+Z_{0}-V_{00}]\succ 0,$
$\displaystyle[\lambda(T(\hat{P}))+T(Z)-V]_{lh}=0\;\;\forall(l,h)\neq(0,0)\\}.$

###  III-C Existence of solutions

The aim of this section is to show that (18) admits solution. The set
$\mathcal{C}$ is not compact, as it is neither closed nor bounded. We show
that we can restrict the search of the minimum of $J$ over a compact set.
Then, since the objective function is continuous over $\mathcal{C}$ (and hence
over the restricted compact set), we can use Weierstrass’s Theorem to conclude
that the problem does admit a minimum.

The first step consists in showing that we can restrict $\mathcal{C}$ to a
subset where $\lambda\geq\varepsilon$ with $\varepsilon>0$ a positive
constant.

###### Proposition III.1

Let $(\lambda^{(k)},V^{(k)},Z^{(k)})_{k\in\mathbb{N}}$ be a sequence of
elements in $\mathcal{C}$ such that

$\lim_{k\to\infty}\lambda^{(k)}=0.$

Then, such a sequence cannot be an infimizing sequence.

The proof is essentially the same as the proof of Proposition 6.1 in [6] and
it is therefore omitted.
As a consequence, minimizing the dual functional over the set $\mathcal{C}$ is
equivalent to minimize it over the set:

$\displaystyle\mathcal{C}_{1}$
$\displaystyle:=\\{(\lambda,V,Z):V\in\mathbf{Q}_{m(n+1)},V\succeq
0,Z\in\mathcal{O},$ $\displaystyle I+V-T(Z)\succeq
0,\lambda\in\mathbb{R},\lambda\geq\varepsilon,[\lambda\hat{P}_{0}+Z_{0}-V_{00}]\succ
0,$
$\displaystyle[\lambda(T(\hat{P}))+T(Z)-V]_{lh}=0\;\forall(l,h)\neq(0,0)\\}.$

Next we show that we can restrict $\mathcal{C}_{1}$ to a subset in which both
$(T(Z)-V)$ and $\lambda$ cannot diverge.

###### Proposition III.2

Let $(\lambda^{(k)},V^{(k)},Z^{(k)})_{k\in\mathbb{N}}$ be a sequence of
elements in $\mathcal{C}_{1}$ such that either

$\lim_{k\to\infty}\parallel T(Z^{(k)})-V^{(k)}\parallel=+\infty$

or

$\lim_{k\to\infty}\lambda^{(k)}=+\infty$

or both. Then, such a sequence cannot be an infimizing sequence.

The above result is obtained by following arguments similar to the proof of
Proposition 6.2 in [6] with a few small differences; we refer the interested
reader to [23, Appendix C] for the detailed proof.
It follows from the previous proposition that there exists
$\beta\in\mathbb{R}$ with $\mid\beta\mid<\infty$ such that $T(Z)-V\succeq\beta
I,$ and $0<\gamma<\infty$ such that $\lambda\leq\gamma$. Therefore, the set
$\mathcal{C}_{1}$ can be further restricted to the set:

$\displaystyle\mathcal{C}_{2}$
$\displaystyle:=\\{(\lambda,V,Z):V\in\mathbf{Q}_{m(n+1)},V\succeq
0,Z\in\mathcal{O},\lambda\in\mathbb{R},$ $\displaystyle\beta I\preceq
T(Z)-V\preceq
I,\gamma\geq\lambda\geq\varepsilon,[\lambda\hat{P}_{0}+Z_{0}-V_{00}]\succ 0,$
$\displaystyle[\lambda(T(\hat{P}))+T(Z)-V]_{lh}=0\;\forall(l,h)\neq(0,0)\\}.$

In addition, it is not possible for $V$ and $Z$ to diverge while keeping the
difference $T(Z)-V$ finite. Accordingly, we can further restrict the search
for the optimal solution to a subset $\mathcal{C}_{3}$ in which neither $V$
nor $Z$ can diverge:

###### Proposition III.3

Let $(\lambda^{(k)},V^{(k)},Z^{(k)})_{k\in\mathbb{N}}$ be a sequence of
elements in $\mathcal{C}_{2}$ such that

$\lim_{k\to\infty}\parallel V^{(k)}\parallel=+\infty$ (19)

or

$\lim_{k\to\infty}\parallel Z^{(k)}\parallel=+\infty$ (20)

or both. Then, such a sequence cannot be an infimizing sequence.

The proof can be found in the appendix.
Thus, the minimization over $\mathcal{C}_{2}$ is equivalent to the
minimization over the subset:

$\displaystyle\mathcal{C}_{3}$
$\displaystyle:=\\{(\lambda,V,Z):V\in\mathbf{Q}_{m(n+1)},\alpha I\succeq
V\succeq 0,Z\in\mathcal{O},\lambda\in\mathbb{R},$ $\displaystyle\beta I\preceq
T(Z)-V\preceq
I,\gamma\geq\lambda\geq\varepsilon,[\lambda\hat{P}_{0}+Z_{0}-V_{00}]\succ 0,$
$\displaystyle[\lambda(T(\hat{P}))+T(Z)-V]_{lh}=0\;\forall(l,h)\neq(0,0)\\}$

for a certain $\alpha>0$ positive constant.

Finally, we consider a sequence
$(\lambda^{(k)},V^{(k)},Z^{(k)})_{k\in\mathbb{Z}}\in\mathcal{C}_{3}$ such that
$[(\lambda^{(k)})^{-1}\big{(}[Z^{(k)}]_{0}-[V^{(k)}]_{00}\big{)}+\hat{P_{0}}]$
tends to be singular as $k\to\infty$. This implies that
$|(\lambda^{(k)})^{-1}\big{(}[Z^{(k)}]_{0}-[V^{(k)}]_{00}\big{)}+\hat{P_{0}}|$
tends to zero and hence $J\to+\infty$. Thus, such a sequence cannot be an
infimizing sequence. Therefore, the final set $\mathcal{C}_{C}$ is:

$\displaystyle\mathcal{C}_{C}$
$\displaystyle:=\\{(\lambda,V,Z):V\in\mathbf{Q}_{m(n+1)},\alpha I\succeq
V\succeq 0,Z\in\mathcal{O},\beta I\preceq$ $\displaystyle T(Z)-V\preceq
I,\lambda\in\mathbb{R},\gamma\geq\lambda\geq\varepsilon,[\lambda\hat{P}_{0}+Z_{0}-V_{00}]\succeq\mu
I,$
$\displaystyle[\lambda(T(\hat{P}))+T(Z)-V]_{lh}=0\;\forall(l,h)\neq(0,0)\\}$

where $\alpha,\beta,\gamma,\varepsilon$ and $\mu$ such that
$|\alpha|,|\beta|,|\gamma|,|\varepsilon|$ and $|\mu|<+\infty.$

###### Theorem III.1

Problem (18) is equivalent to

$\min_{(\lambda,V,Z)\in\mathcal{C}_{C}}J(\lambda,V,Z)$

and it admits solution.

###### Proof:

Equivalence of the two problems has already been proven by the previous
arguments. Since $\mathcal{C}_{C}$ is closed and bounded, hence compact, and
$J$ is continuous over $\mathcal{C}_{C}$, by the Weierstrass’s Theorem the
minimum exists. $\;\blacksquare$ ∎

##  IV Solution of the primal problem

In this section, after proving that the primal problem (5) and its matrix
reformulation (11) are equivalent, we show how to recover the solution of the
primal problem.

Let $(\lambda^{\circ},V^{\circ},Z^{\circ})$ be a solution of (18) and
$(X^{\circ},L^{\circ})$ be the corresponding solution of (11). Since
$X_{00}^{\circ}$ is positive definite, $\log|X_{00}^{\circ}|$ is finite. By
Lemma III.1, at the optimum $\int\log|\Phi|$ must be finite as well; this
implies that $\Phi(e^{i\vartheta})$, $\vartheta\in\left[-\pi,+\pi\right],$ may
be singular at most on a set of zero measure, or, in other terms, $\Delta
X^{\circ}\Delta^{*}\succ 0\text{ a.e.}$. This observation leads to the
following proposition:

###### Proposition IV.1

Let $(X^{\circ},L^{\circ})$ be a solution of (11). Then $\Delta
X^{\circ}\Delta^{*}\succ 0\text{ a.e.}$. Accordingly, (5) and (11) are
equivalent.

Now we are ready to show how to recover the solution of the primal problem; to
this aim we need the following result, see [24].

###### Lemma IV.1

Let $Z\in\mathbf{M}_{m,n}$ and $W\in\mathbf{Q}_{m}$. If $W\succ 0$ is such
that

$\displaystyle T(Z)\succeq\left[\begin{array}[]{cc}W&0\\\
0&0\end{array}\right]$ (23)

then $T(Z)\succ 0$.

Exploiting the constraints
$[\lambda(T(\hat{P}))+T(Z)-V]_{lh}=0,\forall(l,h)\neq(0,0)$ and
$[\lambda\hat{P}_{0}+Z_{0}-V_{00}]\succ 0$, it is not difficult to see that

$\displaystyle
V^{\circ}=\lambda^{\circ}T(\hat{P})+T(Z^{\circ})-\left[\begin{array}[]{cc}W^{\circ}&0\\\
0&0\end{array}\right]$ (26)

where

$\displaystyle
W^{\circ}:=Z^{\circ}_{0}-V^{\circ}_{00}+\lambda^{\circ}\hat{P}_{0}\succ 0.$
(27)

Since $V^{\circ}\succeq 0$ and in view of Lemma IV.1,
$\lambda^{\circ}T(\hat{P})+T(Z^{\circ})\succ 0$. Hence, $V^{\circ}$ has rank
at least equal to $mn$.

Since the duality gap between (11) and (18) is equal to zero, we have that
$\langle V^{\circ},X^{\circ}-L^{\circ}\rangle=0$, which in turn implies

$\displaystyle V^{\circ}(X^{\circ}-L^{\circ})=0$ (28)

because $V^{\circ},X^{\circ}-L^{\circ}\succeq 0$. Recalling that
$\mathrm{rank}(V^{\circ})\geq mn,$ in view of (28) the matrix
$X^{\circ}-L^{\circ}$ has rank at most equal to $m$. Let
$\mathrm{rank}(X^{\circ}-L^{\circ})=\tilde{m}\leq m$. Then, there exists a
full-row rank matrix $A\in\mathbb{R}^{\tilde{m}\times m(n+1)}$ such that

$X^{\circ}-L^{\circ}=A^{\top}A.$ (29)

By (28), it follows that $V^{\circ}A^{\top}=0.$ Let $Y_{D}:=[v_{o}\quad
v_{1}\quad...\quad v_{l}]\in\mathbb{R}^{m(n+1)\times l}$ denote the matrix
whose columns form a basis of $\ker(V^{\circ})$. Note that the dimension $l$
of the null space of $V^{\circ}$ is at least $\tilde{m}$ because
$\mathrm{Im}(A^{\top})\subseteq\ker(V^{\circ})$ and
$\mathrm{rank}(A^{\top})=\tilde{m}$; also $l\leq m$ because
$\mathrm{rank}(V^{\circ})\geq mn$. Rewriting the matrix $A^{\top}$ as
$A^{\top}=Y_{D}S$ with $S\in\mathbb{R}^{l\times\tilde{m}}$, from (29) we
obtain

$X^{\circ}-L^{\circ}=Y_{D}Q_{D}Y_{D}^{\top},$ (30)

with $Q_{D}:=SS^{\top}\in\mathbf{Q}_{l}$ unknown.

In a similar fashion, by the zero duality gap between (11) and (18), the
complementary slackness condition for the multiplier associated to the
positive semi-definiteness of $L$ reads as $\langle
U^{\circ},L^{\circ}\rangle=0,$ which in turn implies $U^{\circ}L^{\circ}=0.$
Repeating the same reasoning as before, it can be seen that, if the dimension
of the null space of ${U^{\circ}}$ is $\tilde{r}$ with $\tilde{r}\geq r$ and
$Y_{L}:=[u_{o}\quad u_{1}\quad...\quad
u_{\tilde{r}}]\in\mathbb{R}^{m(n+1)\times\tilde{r}}$ is a matrix whose columns
form a basis of $\ker(U^{\circ})$, then $L^{\circ}$ can be written as

$\displaystyle L^{\circ}=Y_{L}Q_{L}Y_{L}^{\top}$ (31)

with $Q_{L}\in\mathbf{Q}_{\tilde{r}}$ unknown. Plugging (31) into (30), we
then obtain

$\displaystyle X^{\circ}-Y_{L}Q_{L}Y_{L}^{\top}=Y_{D}Q_{D}Y_{D}^{\top}.$ (32)

Assume now that each block of $X^{\circ}-L^{\circ}$ is diagonal, namely

$\displaystyle\operatorname{ofd}(\left[Y_{D}Q_{D}Y_{D}^{\top}\right]_{hk})=0\quad
h,k=0,...,n.$ (33)

###### Remark 1

We can make the previous assumption without loss of generality. Indeed, let
$(\Phi^{\circ}$, $\Phi_{L}^{\circ})$ be the solution of Problem (5) and
$\Phi_{D}^{\circ}=\Phi^{\circ}-\Phi_{L}^{\circ}$; $X$, $L$ and $D=X-L$ are any
matrices in $\mathbf{Q}_{m(n+1)}$ such that $\Phi^{\circ}=\Delta X\Delta^{*}$,
$\Phi^{\circ}_{L}=\Delta L\Delta^{*}$ and $\Phi^{\circ}_{D}=\Delta
D\Delta^{*}$. We can always consider a different matrix parametrization
$(\tilde{X},\tilde{L},\tilde{D})$ for $\Phi^{\circ}$, $\Phi^{\circ}_{L}$ and
$\Phi^{\circ}_{D}$ as follows. First notice that there always exists a matrix
$\tilde{D}$ with all diagonal blocks such that
$\Phi^{\circ}_{D}=\Delta\tilde{D}\Delta^{*}$; in other words, we can always
find $\delta D\in\mathbf{Q}_{m(n+1)}$ such that $\Delta\delta D\Delta^{*}=0$
and $\tilde{D}:=D+\delta D$ satisfies
$\operatorname{ofd}(\big{[}\tilde{D}\big{]}_{hk})=0$ for $h,k=0,...,n.$ Now,
let $\delta X\in\mathbf{Q}_{m(n+1)}$ such that $\Delta\delta X\Delta^{*}=0$
and $\tilde{X}:=X+\delta X$ satisfies (15). Define
$\tilde{L}=\tilde{X}-\tilde{D}=X-D+\delta L$ with $\delta L:=\delta
X-\delta{D}$. It is easy to see that $\Phi^{\circ}=\Delta\tilde{X}\Delta^{*}$
and $\hat{\Phi}_{L}=\Delta\tilde{L}\Delta^{*}$. This means that
$(\tilde{X},\tilde{L})$ is still a solution of Problem (11) and it allows us
to restrict to solutions of (11) for which (33) holds.

By applying the $\operatorname{ofd}$ operator to both sides of (32) and
exploiting the assumption (33), it is not difficult to obtain:

$\operatorname{ofd}(\left[Y_{L}Q_{L}Y_{L}^{\top}\right]_{00})=\operatorname{ofd}(X^{\circ}_{00})$
(34)

which is a system of $m(m-1)/2$ linear equations in the
$\tilde{r}(\tilde{r}+1)/2$ unknowns $Q_{L}$. Notice that $X_{00}$ is given by
(15). Finally, once $L^{\circ}$ is computed, in order to retrieve $Q_{D}$ we
exploit (33) and the following system of $m(m+1)/2$ linear equations:

$\displaystyle\left[Y_{D}Q_{D}Y_{D}^{\top}\right]_{00}=X^{\circ}_{00}-L^{\circ}_{00}.$
(35)

Since both the dual and the primal problem admit solution, the resulting
systems of equations (33), (34) and (35) do admit solutions.

##  V The proposed algorithm

In this section we propose an algorithm to solve numerically the dual problem.
To start with, as observed in Section IV, we rewrite (18) in a different
fashion by getting rid of the slack variable $V\in\mathbf{Q}_{m(n+1)}$. This
is done by introducing a new variable $W\in\mathbf{Q}_{m}$ defined, similarly
to (27), as

$\displaystyle W:=Z_{0}-V_{00}+\lambda\hat{P}_{0}\succ 0$ (36)

such that, as in (26), the variable $V$ can be expressed as

$\displaystyle V=\lambda T(\hat{P})+T(Z)-\left[\begin{array}[]{cc}W&0\\\
0&0\end{array}\right].$ (39)

Accordingly, the dual problem (18) can be expressed in terms of the variables
$\lambda$, $W$ and $Z$ as follows:

$\min_{(\lambda,W,Z)\in\mathcal{C}}J$ (40)

where

$J:=\lambda\Big{(}-\log\big{|}\lambda^{-1}W\big{|}-\int\log|\hat{\Phi}|+\delta\Big{)}$

and the corresponding feasible set $\mathcal{C}$ is:

$\displaystyle\mathcal{C}:=\\{$
$\displaystyle(\lambda,W,Z):W\in\mathbf{Q}_{m},W\succ
0,Z\in\mathcal{O},\lambda\in\mathbb{R},$ $\displaystyle\lambda>0,\lambda
T(\hat{P})+T(Z)-\left[\begin{array}[]{cc}W&0\\\ 0&0\end{array}\right]\succeq
0,$ $\displaystyle I+\lambda T(\hat{P})-\left[\begin{array}[]{cc}W&0\\\
0&0\end{array}\right]\succeq 0\\}.$

We can further simplify our problem as follows. First, we observe that the
constraint

$V=\lambda T(\hat{P})+T(Z)-\left[\begin{array}[]{cc}W&0\\\
0&0\end{array}\right]\succeq 0$ (41)

implies

$\lambda T(\hat{P})+T(Z)\succeq\left[\begin{array}[]{cc}W&0\\\
0&0\end{array}\right]$

and then, by Lemma IV.1, $\lambda T(\hat{P})+T(Z)\succ 0$. Now, we can easily
rewrite (41) recalling the characterization of a symmetric positive
semidefinite matrix using the Schur complement. To this aim, it is convenient
to introduce the linear operators $T_{0,0}:\mathbf{M}_{m,n}\to\mathbf{Q}_{m},$
$T_{0,1:n}:\mathbf{M}_{m,n}\to\mathbf{M}_{m,n-1}$ and
$T_{1:n,1:n}:\mathbf{M}_{m,n}\to\mathbf{Q}_{mn}$ that, for a given matrix
$H\in\mathbf{M}_{m,n}$ construct a symmetric block-Toeplitz matrix and extract
the blocks in position $(0,0)$, $(0,1:n)$ and $(1:n,1:n)$, respectively. With
this notation, we have

$T(Z+\lambda\hat{P})=\left[\begin{array}[]{cc}T_{0,0}(Z+\lambda\hat{P})&T_{0,1:n}(Z+\lambda\hat{P})\\\
T_{0,1:n}^{\top}(Z+\lambda\hat{P})&T_{1:n,1:n}(Z+\lambda\hat{P})\end{array}\right]$

and the constraint (41) is equivalent to require
$T_{1:n,1:n}(Z+\lambda\hat{P})\succ 0$ and $W\preceq Q(\lambda,Z)$ with

$\displaystyle
Q(\lambda,Z):=T_{0,0}(Z+\lambda\hat{P})-T_{0,1:n}(Z+\lambda\hat{P})\times
T_{1:n,1:n}^{-1}(Z+$
$\displaystyle\lambda\hat{P})T_{0,1:n}^{\top}(Z+\lambda\hat{P}).$

In a similar fashion, the last matricial inequality constraint in
$\mathcal{C}$ can be equivalently expressed as $W\preceq R(\lambda)$ where

$\displaystyle
R(\lambda):=I+T_{0,0}(\lambda\hat{P})-T_{0,1:n}(\lambda\hat{P})\big{(}I+T_{1:n,1:n}(\lambda\hat{P})\big{)}^{-1}\times$
$\displaystyle T_{0,1:n}^{\top}(\lambda\hat{P}).$

Therefore, Problem (18) can be formulated as

$\min_{(\lambda,W,Z)\in\mathcal{C}}J=\lambda\Big{(}-\log\big{|}\lambda^{-1}W\big{|}-\int\log|\hat{\Phi}|+\delta\Big{)}$
(42)

where

$\displaystyle\mathcal{C}:=$
$\displaystyle\\{(\lambda,W,Z):Z\in\mathcal{O},\lambda\in\mathbb{R},\lambda>0,T_{1:n,1:n}(Z+\lambda\hat{P})\succ
0,$ $\displaystyle W\in\mathbf{Q}_{m},W\succ 0,W\preceq
Q(\lambda,Z),\;W\preceq R(\lambda)\\}.$

Solving Problem (42) simultaneously for $\lambda$, $W,$ and $Z$ is not trivial
because the inequality constraints $W\preceq Q(\lambda,Z)$ and $W\preceq
R(\lambda)$ both depend on $\lambda$. On the other hand, once we fix the dual
variable $\lambda$ to a positive constant $\bar{\lambda}>0$, the problem:

$\min_{(W,Z)\in\mathcal{C}_{\bar{\lambda}}}J(\bar{\lambda},W,Z)$ (43)

with

$\displaystyle\mathcal{C}_{\bar{\lambda}}:=\\{$
$\displaystyle(W,Z):Z\in\mathcal{O},\;W\in\mathbf{Q}_{m},\;T_{1:n,1:n}(Z+\bar{\lambda}\hat{P})\succ
0,$ $\displaystyle W\succ 0,\;W\preceq Q(\bar{\lambda},Z),\;W\preceq
R(\bar{\lambda})\\}.$

can be efficiently solved by resorting to the ADMM algorithm [25]. To this
aim, we rewrite Problem (43) by introducing a new variable
$Y\in\mathbf{Q}_{m}$ defined as $Y=Q(\bar{\lambda},Z)-W:$

$\displaystyle\min_{\begin{subarray}{c}(W,Z)\in\mathcal{C}_{W,Z},\\\
{Y\in\mathbf{Q}_{m}^{+}}\end{subarray}}$ $\displaystyle
J=\bar{\lambda}\big{(}-\log\big{|}\bar{\lambda}^{-1}W\big{|}-\int\log|\hat{\Phi}|+\delta\big{)}$
(44) subject to $\displaystyle Y=Q(\bar{\lambda},Z)-W$

where

$\displaystyle\mathcal{C}_{W,Z}:=\\{$
$\displaystyle(W,Z):Z\in\mathcal{O},W\in\mathbf{Q}_{m},\;W\succ 0,$
$\displaystyle W\preceq
R(\bar{\lambda}),\;T_{1:n,1:n}(Z+\bar{\lambda}\hat{P})\succ 0\\}$

and $\mathbf{Q}_{m}^{+}$ denotes the cone of symmetric positive semidefinite
matrices of size $m\times m$. The augmented Lagrangian for (44) is:

$\displaystyle\mathcal{L}_{\rho}(W,Z,Y,M):=\bar{\lambda}\Big{(}-\log\big{|}\bar{\lambda}^{-1}W\big{|}-\int\log|\hat{\Phi}|+\delta\Big{)}+$
$\displaystyle\langle M,Y-Q(\bar{\lambda},Z)+W\rangle+\frac{\rho}{2}\parallel
Y-Q(\bar{\lambda},Z)+W\parallel^{2}$

where $M\in\mathbf{Q}_{m}$ is the Lagrange multiplier, and $\rho>0$ is the
penalty parameter. Accordingly, given the initial guesses $W^{(0)},$
$Z^{(0)},$ $Y^{(0)}$ and $M^{(0)}$, the ADMM updates are:

$\displaystyle(W^{(k+1)},Z^{(k+1)})=\operatorname*{arg\,min}_{(W,Z)\in\mathcal{C}_{W,Z}}\mathcal{L}_{\rho}(W,Z,Y^{(k)},M^{(k)})$
(45) $\displaystyle
Y^{(k+1)}=\operatorname*{arg\,min}_{Y\in\mathbf{Q}_{m}^{+}}\mathcal{L}_{\rho}(W^{(k+1)},Z^{(k+1)},Y,M^{(k)})$
(46) $\displaystyle
M^{(k+1)}=M^{(k)}+\rho\big{(}Y^{(k+1)}-Q(\bar{\lambda},Z^{(k+1)})+W^{(k+1)}\big{)}.$

Problem (45) does not admit a closed form solution, therefore we approximate
the optimal solution by a gradient projection step:

$\displaystyle W^{(k+1)}$
$\displaystyle=\Pi\big{(}W^{(k)}-t_{k}\nabla_{W}\mathcal{L}_{\rho}(W^{(k)},Z^{(k)},Y^{(k)},M^{(k)})\big{)}$
$\displaystyle Z^{(k+1)}$
$\displaystyle=\Pi_{\mathcal{O}}\;\big{(}Z^{(k)}-t_{k}\nabla_{Z}\mathcal{L}_{\rho}(W^{(k)},Z^{(k)},Y^{(k)},M^{(k)})\big{)}$

where:

  * •

$\nabla_{W}\mathcal{L}_{\rho}(W,Z,Y,M)$ denotes the gradient of the augmented
Lagrangian with respect to $W$:

$\nabla_{W}\mathcal{L}_{\rho}=-\bar{\lambda}W^{-1}+M+\rho(Y-Q+W).$

  * •

$\nabla_{Z}\mathcal{L}_{\rho}(W,Z,Y,M)$ denotes the gradient of the augmented
Lagrangian with respect to $Z$:

$\displaystyle\nabla_{Z}\mathcal{L}_{\rho}=D\big{(}\left[\begin{array}[]{c}I_{m}\\\
-T^{-1}_{1:n,1:n}T^{\top}_{0,1:n}\end{array}\right]\big{(}-M-\rho(Y-$
$\displaystyle
Q+W)\big{)}\left[\begin{array}[]{cc}I_{m}&-T_{0,1:n}T^{-1}_{1:n,1:n}\end{array}\right]\big{)}$

where the omitted argument of the operators $T_{0,1:n}$ and $T_{1:n,1:n}$ is
intended to be equal to $(Z+\bar{\lambda}\hat{P}).$

  * •

$\Pi_{\mathcal{O}}$ denotes the projection operator onto $\mathcal{O}$:

$\Pi_{\mathcal{O}}(A)=\operatorname{ofd}_{B}(A).$

  * •

$\Pi$ denotes the projection operator onto the convex cone
$\\{S\in\mathbf{Q}_{m}:S\preceq R(\bar{\lambda})\\}.$ It is not difficult to
see that

$\Pi(A)=R(\bar{\lambda})-\Pi_{+}(R(\bar{\lambda})-A),$

where $\Pi_{+}$ is the projection operator onto the cone $\mathbf{Q}_{m}^{+}$.

  * •

the step-size $t_{k}$ is determined at each step $k$ in an iterative fashion:
we start by setting $t_{k}=1$ and we decrease it progressively of a factor
$\beta,$ with $0<\beta<1,$ until the conditions $W^{(k+1)}\succ 0$ and
$T_{1:n,1:n}(Z^{(k+1)}+\bar{\lambda}\hat{P})\succ 0$ are met and the Armijo’s
condition [26] is satisfied.

Problem (46) admits a closed form solution, which can be easily computed as:

$Y^{(k+1)}=\Pi_{+}\Big{(}Q(\bar{\lambda},Z^{(k+1)})-W^{(k+1)}-\frac{1}{\rho}M^{(k)}\Big{)}.$

To define the stopping criterion, we need to introduce the following
quantities

$\displaystyle R^{P}$ $\displaystyle=Y-Q(\bar{\lambda},Z^{(k+1)})+W^{(k+1)}$
$\displaystyle R^{D}$ $\displaystyle=\begin{aligned}
D\big{(}\left[\begin{array}[]{c}I_{m}\\\
-T^{-1}_{1:n,1:n}T^{\top}_{0,1:n}\end{array}\right](\rho(Y^{(k+1)}-Y^{(k)}))\times\\\
\left[\begin{array}[]{cc}I_{m}&-T_{0,1:n}T^{-1}_{1:n,1:n}\end{array}\right]\big{)}\end{aligned}$

which are referred to as the primal and dual residual, respectively. Notice
that the omitted argument of the operators $T_{0,1:n}$ and $T_{1:n,1:n}$ is
intended to be equal to $(Z^{(k+1)}+\bar{\lambda}\hat{P})$.
Then, the algorithm stops when the following conditions are met:

$\displaystyle\|R^{P}\|$ $\displaystyle\leq\begin{aligned}
m\varepsilon^{\text{ABS}}+\varepsilon^{\text{REL}}\max&\\{\|W^{(k)}\|,\|Q(\bar{\lambda},Z^{(k)})\|,\|Y^{(k)}\|\\}\end{aligned}$
$\displaystyle\|R^{D}\|$ $\displaystyle\leq\begin{aligned}
m\sqrt{(n+1)}\varepsilon^{\text{ABS}}+\varepsilon^{\text{REL}}\|D\big{(}\left[\begin{array}[]{c}I_{m}\\\
-T^{-1}_{1:n,1:n}T^{\top}_{0,1:n}\end{array}\right]\times\\\
M^{(k)}\left[\begin{array}[]{cc}I_{m}&-T_{0,1:n}T^{-1}_{1:n,1:n}\end{array}\right]\big{)}\|\end{aligned}$

where $\varepsilon^{\text{ABS}}$ and $\varepsilon^{\text{REL}}$ are the
desired absolute and relative tolerances.

It remains to determine the optimal value $\lambda^{\circ}$ for $\lambda$
which solves Problem (42). To this aim, we exploit the following result (see
[26, pp.87-88]):

###### Proposition V.1

If $f$ is convex in $(x,y)$ and $\mathcal{C}$ is a convex non-empty set, then
the function

$g(x)=\inf_{y\in\mathcal{C}}f(x,y)$ (47)

is convex in $x$, provided that $g(x)>-\infty$ for some $x$. The domain of $g$
is the projection of $\operatorname{dom}(f)$ on its $x$-coordinates.

This result guarantees that the function

$g(\lambda)=\min_{(W,Z)\in\mathcal{C}_{\lambda}}J(\lambda,W,Z)$

is convex in $\lambda$. Hence, in order to determine
$\lambda^{\circ}=\operatorname*{arg\,min}_{\lambda>0}g(\lambda)$ we can choose
an initial interval of uncertainty $[a,b]$ containing $\lambda^{\circ}$, and
we progressively reduce it by evaluating $g(\lambda)$ at two points within the
interval placed symmetrically, each at distance $h>0$ from the midpoint. This
is repeated until the width of the uncertainty interval is smaller than a
certain tolerance $l>0$.

The overall procedure to solve the dual problem (42) is summarized in
Algorithm 1.

Algorithm 1

Input: $b>a>0$, $l>0$, $h>0$
Output: $(\lambda^{\circ},W^{\circ},Z^{\circ})$

1:repeat

2: $\tilde{a}=(a+b)/2-h;\;\tilde{b}=(a+b)/2+h.$

3: Compute $g(\tilde{a})$ by applying the ADMM with $\lambda=\tilde{a}.$

4: Compute $g(\tilde{b})$ by applying the ADMM with $\lambda=\tilde{b}.$

5: if $g(\tilde{a})<g(\tilde{b})$ then

6: $b=\tilde{b}$

7: else

8: $a=\tilde{a}$

9: end if

10:until $b-a<l$

11:$\lambda^{\circ}=(a+b)/2$.

12:Compute $(W^{\circ},Z^{\circ})$ by applying the ADMM with
$\lambda=\lambda^{\circ}.$

##  VI Identification of ARMA factor models

In this section we extend the proposed approach to ARMA processes. Consider
the ARMA factor model:

$y(t)=a^{-1}(W_{L}u(t)+W_{D}w(t))$ (48)

where

$a(e^{i\vartheta})=\sum_{k=0}^{p}a_{k}e^{-i\vartheta k},\quad
a_{k}\in\mathbb{R}$

and $W_{L},W_{D}$, $u$ and $w$ are defined analogously to (1). Notice that
$y_{MA}(t):=ay(t)=W_{L}u(t)+W_{D}w(t)$ is a MA process of order $n$ whose
spectral density $\Phi=W_{L}W_{L}^{*}+W_{D}W_{D}^{*}\in\mathcal{Q}_{m,n}$
admits a low rank plus diagonal decomposition. Finally, it is worth noting
that it is not restrictive to assume that the autoregressive part in (48) is
characterized by a scalar filter $a$; Indeed, any ARMA factor model can be
written in the form of (48).

Assume now to collect a realization
$\mathrm{y}^{N}=\\{\,\mathrm{y}(1)\ldots\mathrm{y}(N)\,\\}$ of numerosity $N$
of the process $y$. Our aim is to estimate the factor model (48) and the
number of factors $r$. Before proceeding, the following observation needs to
be made: there is an identifiability issue in the problem. Indeed, if we
multiply $a(z)$, $W_{L}$ and $W_{D}$ by an arbitrary non-zero real number $c$,
the model remains the same. We can easily eliminate this uninteresting degree
of freedom by normalizing the polynomial $a(z)$, so that from now on we assume
$a_{0}=1$.

The idea is to estimate first $a$, and then $\Phi_{L}$ and $\Phi_{D}$ by
preprocessing $\mathrm{y}^{N}$ through $a$. In more detail, the proposed
solution consists of the following two steps:

  1. 1.

The _AR dynamic estimation._ Given the realization $\mathrm{y}^{N}$, we
estimate the $p$ parameters of the filter $a$ by applying the maximum
likelihood estimator proposed in [27, Section II.b]. In doing so, we are
estimating an AR process whose spectral density is $a^{-1}(a^{-1})^{*}I_{m}$.

  2. 2.

The _MA dynamic factor analysis._ Let $\mathrm{y}_{MA}^{N}$ be the finite
length trajectory obtained by passing through the filter
$a^{\circ}(e^{i\vartheta})$ the trajectory $\mathrm{y}^{N}$ with zero initial
conditions. After computing the truncated periodogram
$\hat{\Phi}\in\mathcal{Q}_{m,n}$ from $\mathrm{y}_{MA}^{N}$, we solve Problem
(5) with $\hat{\Phi}$ in order to recover the number of latent factors.

Although the above procedure is suboptimal, the numerical simulations showed
that the resulting estimator of the number of factors performs well, see
Section VII-B.

##  VII Numerical simulations

In this section, we test the performance of the proposed approach both for MA
and ARMA factor models. In all the simulations, the parameter $\delta$ is
computed according to the empirical procedure of Section II-A for
$\alpha=0.5$. Then, Problem (42) is solved by applying Algorithm 1 with $l=7$
and $h=3.$ In regard to the ADMM algorithm, we set
$\varepsilon^{\text{ABS}}=10^{-4}$, $\varepsilon^{\text{REL}}=10^{-4}$ and the
penalty term $\rho=0.05.$

###  VII-A Synthetic Example - MA factor models

Figure 1: Estimated MA factor model with $n=2$, $m=40$, and $r=2$. Integral
over the unit circle of the first 20 normalized singular values of
$\Phi_{L}^{\circ}$ with $N=5000$. Figure 2: Estimated MA factor model with
$n=2$, $m=40$, and $r=4$. Integral over the unit circle of the first 30
normalized singular values of $\Phi_{L}^{\circ}$ with $N=5000$

We consider an MA factor model (1) of order $n=2$, with $m=40$ manifest
variable and $r=2$ latent factors, computed by randomly generating the zeros
of the transfer functions $[W_{L}]_{(i,j)}$’s and $[W_{D}]_{(i,i)}$’s for
$i=1,\dots,m,$ $j=1,\dots,r$ within the circle with center at the origin and
radius $0.95$ on the complex plane. It is worth noting that
$\int\|\Phi_{L}(e^{i\theta}))\|=141.83$ and
$\int\|\Phi_{D}(e^{i\theta}))\|=31.29$, that is the idiosyncratic component is
not negligible with respect to the latent variable. We generate from the model
a sample $\mathrm{y}^{N}$ of length $N=5000$ and we apply the proposed
identification procedure to estimate the number of common factors. We define

$s_{j}:=\int\frac{\sigma_{j}(\Phi_{L}^{\circ}(e^{i\theta}))}{\sigma_{1}(\Phi_{L}^{\circ}(e^{i\theta}))}$

where $\sigma_{j}(\Phi_{L}^{\circ}(e^{i\theta})$ denotes the $j-th$ largest
eigenvalue of $\Phi_{L}^{\circ}$ at frequency $\theta$. It is clear that
$s_{j}$ represents the integral of the $j-th$ largest normalized singular
value of $\Phi_{L}^{\circ}$ over the unit circle. The quantities $s_{j}$ are
plotted in Figure 1; we can notice that there is a knee point at $j=2$, so
that the numerical rank of $\Phi_{L}^{\circ}$ is equal to $2$ and, in doing
so, we can recover the exact number of common factors.

The effectiveness of the proposed estimator is also tested by considering a
data sample of numerosity $N=5000$ generated by a MA factor model with $n=2$,
$m=40$ and $r=4$ common factors; the integrals of the Frobenius norm of
$\Phi_{L}$ and $\Phi_{D}$ are $80.76$ and $8.89,$ respectively. As showed in
Figure 2, even in this case we are able to estimate the correct number of
latent variables.

Finally, we obtained similar results with different samples and by changing
the “true” factor model from which we generated the data.

###  VII-B Synthetic Example - ARMA factor models

Figure 3:  Estimated ARMA factor model with $m=40,$ $r=2,$ $n=2$ and $p=2.$
Box-plot of the integral over the unit circle of the first 15 normalized
singular values of $\Phi_{L}^{\circ}$ with $N=5000$.

To provide empirical evidence of the estimation performance of the algorithm
of Section VI, a Monte Carlo simulation study composed of 50 experiments is
performed. We randomly build an ARMA factor model (48) with $m=40$, $r=2$,
$n=2$ and $p=2$; without loss of generality we fix $a_{0}=1$. Then, for each
Monte Carlo experiment a data sequence of length $N=5000$ is randomly
generated from the model and the ARMA factor model identification procedure is
performed. The boxplot of the quantities $s_{j}$ for the estimated
$\Phi_{L}^{\circ}$’s are shown in Figure 3 and it reveals that the proposed
identification procedure is able to successfully recover the number of latent
factors.

###  VII-C Smart Building Dataset

The SMLsystem is a house built in Valencia at the Universidad CEU Cardenal
Herrera (CEU-UCH). It is a modular house that integrates a whole range of
different technologies to improve energy efficiency, with the objective to
construct a near zero-energy house. A complex monitoring system has been used
in the SMLsystem: it has indoor sensors for temperature, humidity and carbon
dioxide; outdoor sensors are also available for lighting measurements, wind
speed, rain, sun irradiance and temperature. We refer the reader to [28] for a
detailed description of the building and its monitoring system. Two datasets
from the SMLsystem are available for download at the UCI Machine Learning
repository http://archive.ics.uci.edu/ml. We take into account $m=17$ sensor
signals extracted from these datasets: the indoor temperature (in
${}^{\circ}\text{C}$) of the dinning-room and of the room, the weather
forecast temperature (in ${}^{\circ}\text{C}$), the carbon dioxide (in ppm) in
the dinning room and in the room, the relative humidity (in %) in the dinning
room and the room, the lighting in the dinning room and the room (in lx), the
sun dusk, the wind (in cm/sec), the sun light (in klx) in the west, east and
south facade, the sun irradiance (in dW), the outdoor temperature (in
${}^{\circ}\text{C}$) and finally the outdoor relative humidity (in %). The
data are sampled with a period of $T=15\text{min}$ and each sample is the mean
of the last quarter, reducing in this way the signal noise. The first dataset
$\mathrm{y}^{N_{1}}=\\{\,\mathrm{y}(1),\ldots,\mathrm{y}(N_{1})\,\\}$ was
captured during March 2011 and has $N_{1}=2764$ points ($\approx 28$ days),
while the second dataset
$\mathrm{y}^{N_{2}}=\\{\,\mathrm{y}(N_{1}+1),\ldots,\mathrm{y}(N_{1}+N_{2})\,\\}$
has $N_{2}=1373$ points ($\approx 14$ days) collected in June 2011.

It is reasonable to expect that the variability of the considered signals may
be successfully explained by a smaller number of factors. Motivated by this
reason, we apply the ARMA factor model identification procedure with
parameters $n=2$ and $p=2$ using the realization $\mathrm{y}^{N_{1}}$. As
shown in Figure 4, we obtain an estimate of $r=4$ latent factors.

For the sake of comparison, we also use the Matlab function armax() of the
System Identification Toolbox to compute the prediction-error method (PEM)
estimate for an ARMA model with polynomials $A(z)$ and $B(z)$ of order 2 and
$A(z)$ diagonal from the realization $\mathrm{y}^{N_{1}}$.

Finally, the second dataset $\mathrm{y}^{N_{2}}$ is used in the validation
step to test the prediction capability of the two estimated ARMA models. The
results are summarized in Figure 5 which displays for each output channel
$j=1,\dots,m$ the fit (percentage) term:

$J_{FIT,j}:=100\left(1-\frac{\sqrt{\sum_{t=N_{1}+1}^{N_{1}+N_{2}}(\mathrm{y}_{j}(t)-\hat{\mathrm{y}}_{j}(t|t-1))^{2}}}{\sqrt{\sum_{t=N_{1}+1}^{N_{1}+N_{2}}(\mathrm{y}_{j}(t)-\bar{\mathrm{y}}_{j})^{2}}}\right)$

where
$\bar{\mathrm{y}}_{j}:=\frac{1}{N_{2}}\sum_{t=N_{1}+1}^{N_{1}+N_{2}}\mathrm{y}_{j}(t)$
and $\hat{\mathrm{y}}_{j}(t|t-1)$ is the one-step ahead prediction at time $t$
computed with zero initial conditions. The figure shows that the ARMA factor
model matches quite well the measurement data $\mathrm{y}^{N_{2}}$, reaching
fit values similar to the PEM estimate. This is a remarkable result as the
performances of the two approaches are essentially the same, but the factor
model is parameterized by 257 coefficients, much less than the 869
coefficients of the PEM estimate.

Figure 4: Application of the ARMA factor models identification procedure by
using the measurements $\mathrm{y}^{N_{1}}$ from the SMLsystem as training
data. The figure shows the integral over the unit circle of the normalized
singular values of $\Phi_{L}^{o}$.  Figure 5: Fit (in percentage) term
$J_{FIT,j}$ for each output channel for the model estimated via factor
analysis and via PEM. The fit values are computed by using the measurements
$\mathrm{y}^{N_{2}}$ from the SMLsystem as validation data.

##  VIII Conclusion

A procedure to estimate the number of factors and to learn ARMA factor models
has been proposed. This method is based on the solution of an optimization
problem whose solution has been proven to exist via dual analysis. The
simulations results applying the procedure both to synthetic and real data
provide evidence of a good performance.

### Proof of Lemma III.1

Since $\Phi=\Delta X\Delta^{*}$ with $X\succeq 0$, there exists
$A\in\mathbb{R}^{m\times m(n+1)}$ such that $X=A^{\top}A$. The matrix $A$ is
such that $\Phi\succeq 0$ admits the spectral factorization $\Phi=WW^{*}$
where $W:=\Delta A^{\top}$. Now, define $\Phi_{n}:=\Phi+\frac{1}{n}I$ with
$n\in\mathbb{N}$ and let $W_{n}:=\Delta A_{n}$ be a spectral factor of
$\Phi_{n}$ with $A_{n}\in\mathbb{R}^{m\times m(n+1)}.$ Clearly,
$\lim_{n\to+\infty}\Phi_{n}=\Phi$; accordingly, $\lim_{n\to+\infty}W_{n}=W$
and $\lim_{n\to+\infty}A_{n}=A$. Since $\Phi_{n}\succ 0$ $\forall\vartheta$ we
can exploit (10) to obtain

$\int\log|\Phi_{n}|=\log|A_{n_{0}}^{\top}A_{n_{0}}|.$

Then, applying the limit operator to both sides, we have

$\displaystyle\lim_{n\to+\infty}\int\log|\Phi_{n}|=\log|A_{0}^{\top}A_{0}|=\log|X_{00}|.$

To conclude the proof, it remains to show that in the left side of the
previous equation it is possible to interchange the limit and the integral
operators. To this aim, we introduce the sequence
$\\{f_{n}\\}_{n=1}^{+\infty}$ where $f_{n}(t):=\log|\Phi_{n}(\vartheta)|$ and
the function
$f(\vartheta):=\lim_{n\rightarrow+\infty}f_{n}(t)=\log|\Phi(\vartheta)|$.
Observe that, since the interval of integration $[-\pi,\pi]$ is bounded and
$f_{1}(\vartheta)<+\infty$ for any $\vartheta\in[-\pi,\pi]$, then $\int
f_{1}(\vartheta)d\vartheta<+\infty.$ We also define the sequence
$\\{g_{n}\\}_{n=1}^{+\infty}$ as
$g_{n}(\vartheta):=f_{n}(\vartheta)-f_{1}(\vartheta)$ and
$g(\vartheta):=\lim_{n\rightarrow+\infty}g_{n}(\vartheta).$ $\\{g_{n}\\}$ is a
pointwise non-increasing sequence of measurable non-positive functions,

$\dots\leq g_{2}(\vartheta)\leq g_{1}(\vartheta)\leq
0,\quad\forall\vartheta\in\left[-\pi,+\pi\right]$

converging to $g(\vartheta)$ from above. Hence, it satisfies all the
hypotheses of Beppo-Levi’s monotone convergence theorem (applied with opposite
signs), from which it immediately follows that

$\lim_{n\rightarrow+\infty}\int g_{n}(\vartheta)=\int g(\vartheta),$

and consequently

$\lim_{n\rightarrow+\infty}\int f_{n}(\vartheta)=\int g(\vartheta)+\int
f_{1}(\vartheta).$ (49)

Now, since $f_{1}(\vartheta)<+\infty$ for all $\vartheta$,

$g(\vartheta)=f(\vartheta)-f_{1}(\vartheta),$ (50)

and, by plugging (50) into (49), we finally obtain

$\lim_{n\rightarrow+\infty}\int f_{n}(\vartheta)=\int
f(\vartheta).\quad\blacksquare$

### Proof of Proposition III.3

Consider a sequence $(\lambda^{(k)},V^{(k)},Z^{(k)})_{k\in\mathbb{N}}$ in
$\mathcal{C}_{2}$.

We first show that $[Z^{(k)}]_{0}$ cannot diverge. Indeed, assume by
contradiction that $\lim_{k\to\infty}\parallel[Z^{(k)}]_{0}\parallel=+\infty$.
Since it is a symmetric and traceless matrix, this implies

$\lim_{k\to\infty}\quad\min_{\alpha^{(k)}\in\sigma\big{(}[Z^{(k)}]_{0}\big{)}}\alpha^{(k)}=-\infty.$
(51)

In view of (51), since $\lambda^{(k)}\hat{P}_{0}$ is bounded and $V^{(k)}$
positive semidefinite $\forall k$, then
$(\lambda^{(k)}\hat{P}_{0}+[Z^{(k)}]_{0}-[V^{(k)}]_{00})$ has at least a
negative eigenvalue for $k$ sufficiently large, so that the sequence
$(\lambda^{(k)},V^{(k)},Z^{(k)})$ is not in $\mathcal{C}_{2}$. We conclude
that

$\lim_{k\to\infty}\parallel[Z^{(k)}]_{0}\parallel<\infty.$

As a consequence, since $\beta I\preceq T(Z^{(k)})-V^{(k)}\preceq I$ (which is
one of the condition for the sequence to be in $\mathcal{C}_{2}$ ), and
$[T(Z^{(k)})]_{hh}=[Z^{(k)}]_{0}$ by construction, it holds that $\forall k$

$\parallel[V^{(k)}]_{hh}\parallel<\infty,\quad h=0,\dots,n.$

Then, from $V^{(k)}\succeq 0$ it follows that also the off-diagonal blocks of
$V^{(k)}$ must be bounded $\forall k$, i.e.

$\parallel[V^{(k)}]_{hl}\parallel<\infty,\quad l\neq h,\;\;l,h=0,\dotsc,n.$
(52)

Finally, by the boundedness of $(T(Z^{(k)})-V^{(k)})$ and by (52) we obtain
that $\forall k$

$\parallel[Z^{(k)}]_{h}\parallel<\infty\quad h=1,\dots,n,$ (53)

which concludes the proof. $\;\blacksquare$

## References

  * [1] L. Ning, T. T. Georgiou, A. Tannenbaum, and S. P. Boyd, “Linear models based on noisy data and the Frisch scheme,” _SIAM Review_ , vol. 57, no. 2, pp. 167–197, 2015.
  * [2] D. Bertsimas, M. S. Copenhaver, and R. Mazumder, “Certifiably optimal low rank factor analysis,” _Journal of Machine Learning Research_ , vol. 18, no. 29, pp. 1–53, 2017.
  * [3] V. Ciccone, A. Ferrante, and M. Zorzi, “Learning latent variable dynamic graphical models by confidence sets selection,” _Kybernetika_ , vol. 55, no. 4, pp. 74–754, 2019.
  * [4] M. Zorzi and R. Sepulchre, “Factor analysis of moving average processes,” in _European Control Conference (ECC)_ , Linz, 2015, pp. 3579–3584.
  * [5] V. Ciccone, A. Ferrante, and M. Zorzi, “Factor models with real data: A robust estimation of the number of factors,” _IEEE Transactions on Automatic Control_ , vol. 64, no. 6, pp. 2412–2425, June 2019.
  * [6] ——, “Learning latent variable dynamic graphical models by confidence sets selection,” _IEEE Transactions on Automatic Control_ , vol. 65, no. 12, pp. 5130–5143, 2020.
  * [7] J. Geweke, “The dynamic factor analysis of economic time series,” in _Latent variables in socio-economic models_ , D. Aigner and A. Goldberger, Eds. Amsterdam: North-Holland, 1977.
  * [8] M. Deistler and C. Zinner, “Modelling high-dimensional time series by generalized linear dynamic factor models: An introductory survey,” _Communications in Information & Systems_, vol. 7, no. 2, pp. 153–166, 2007\.
  * [9] J. Stock and M. Watson, “Dynamic factor models,” 2010, internal report. [Online]. Available: https://www.princeton.edu/~mwatson/papers/dfm_oup_4.pdf
  * [10] G. Figá-Talamanca, S. Focardi, and M. Patacca, “Common dynamic factors for cryptocurrencies and multiple pair-trading statistical arbitrages,” _Decisions in Economics and Finance_ , 2021.
  * [11] G. Bottegal and G. Picci, “Modeling complex systems by generalized factor analysis,” _IEEE Transactions on Automatic Control_ , vol. 60, no. 3, pp. 759–774, 2015.
  * [12] J. Songsiri and L. Vandenberghe, “Topology selection in graphical models of autoregressive processes,” _Journal of Machine Learning Research_ , vol. 11, no. Oct, pp. 2671–2705, 2010.
  * [13] E. Avventi, A. Lindquist, and B. Wahlberg, “ARMA identification of graphical models,” _IEEE Trans. Autom. Control_ , vol. 58, no. 5, pp. 1167–1178, May 2013.
  * [14] S. Maanan, B. Dumitrescu, and C. Giurcăneanu, “Conditional independence graphs for multivariate autoregressive models by convex optimization: Efficient algorithms,” _Signal Processing_ , vol. 133, pp. 122–134, 2017\.
  * [15] D. Alpago, M. Zorzi, and A. Ferrante, “Identification of sparse reciprocal graphical models,” _IEEE Control Systems Letters_ , vol. 2, no. 4, pp. 659–664, Oct 2018.
  * [16] ——, “A scalable strategy for the identification of latent-variable graphical models,” _Submitted_ , 2018.
  * [17] M. Zorzi and R. Sepulchre, “AR identification of latent-variable graphical models,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 9, pp. 2327–2340, Sept 2016.
  * [18] S. Maanan, B. Dumitrescu, and C. Giurcăneanu, “Maximum entropy expectation-maximization algorithm for fitting latent-variable graphical models to multivariate time series,” _Entropy_ , vol. 20, p. 76, 01 2018\.
  * [19] M. Zorzi, “Empirical Bayesian learning in AR graphical models,” _Automatica_ , vol. 109, 2019.
  * [20] M. Zorzi, “Autoregressive identification of Kronecker graphical models,” _Automatica_ , vol. 119, p. 109053, 2020.
  * [21] M. S. Veedu, H. Doddi, and M. V. Salapaka, “Topology learning of linear dynamical systems with latent nodes using matrix decomposition,” 2020.
  * [22] M. S. Veedu and M. V. Salapaka, “Topology identification under spatially correlated noise,” 2020.
  * [23] L. Falconi, “Robust factor analysis of moving average processes,” 2021, master thesis. [Online]. Available: http://tesi.cab.unipd.it/65215/
  * [24] J. Songsiri, J. Dahl, and L. Vandenberghe, “Graphical models of autoregressive processes,” _Convex optimization in signal processing and communications_ , pp. 89–116, 2010.
  * [25] S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” _Found. Trends Mach. Learn._ , vol. 3, no. 1, pp. 1–122, Jan. 2011. [Online]. Available: http://dx.doi.org/10.1561/2200000016
  * [26] S. Boyd and L. Vandenberghe, _Convex optimization_. Cambridge, United Kingdom: Cambridge University Press, 2004\.
  * [27] F. Crescente, L. Falconi, F. Rozzi, A. Ferrante, and M. Zorzi, “Learning ar factor models,” in _2020 59th IEEE Conference on Decision and Control (CDC)_ , 2020, pp. 274–279.
  * [28] F. Zamora-Martínez, P. Romeu, P. Botella-Rocamora, and J. Pardo, “On-line learning of indoor temperature forecasting models towards energy efficiency,” _Energy and Buildings_ , vol. 83, pp. 162 – 172, 2014. [Online]. Available: http://www.sciencedirect.com/science/article/pii/S0378778814003569


# Tailoring Nanowire Lasing Modes via Coupling to Metal Gratings

F. Vitale  Institute of Solid State Phyisics, Friedrich Schiller University
Jena, Max-Wien-Platz 1, 07743 Jena, Germany<EMAIL_ADDRESS>D.
Repp  Institute of Applied Physics, Abbe Center of Photonics,Friedrich
Schiller University Jena, Albert-Einstein-Straße 15, 07745 Jena, Germany  T.
Siefke  Institute of Applied Physics, Abbe Center of Photonics,Friedrich
Schiller University Jena, Albert-Einstein-Straße 15, 07745 Jena, Germany  U.
Zeitner  Institute of Applied Physics, Abbe Center of Photonics,Friedrich
Schiller University Jena, Albert-Einstein-Straße 15, 07745 Jena, Germany
Fraunhofer Institute for Applied Optics and Precision Engineering, Albert-
Einstein-Straße 7, 07745 Jena, Germany  U. Peschel  Institute of Condensed
Matter Theory and Solid State Optics, Friedrich Schiller University Jena, Max-
Wien-Platz 1, Jena 07743, Germany  T. Pertsch  Institute of Applied Physics,
Abbe Center of Photonics,Friedrich Schiller University Jena, Albert-Einstein-
Straße 15, 07745 Jena, Germany  C. Ronning  Institute of Solid State Phyisics,
Friedrich Schiller University Jena, Max-Wien-Platz 1, 07743 Jena, Germany

###### Abstract

Tailoring the emission of plasmonic nanowire-based lasers represents one of
the major challenges in the field of nanoplasmonics, given the envisaged
integration of such devices into on-chip all-optical circuits. In this study,
we proposed a mode selection scheme based on distributed feedback, achieved
via the external coupling of single zinc oxide nanowires to an aluminum
grating, which enabled a quasi-single mode lasing action. The nano-
manipulation of a single nanowire allowed for a reliable comparison of the
lasing emission characteristics in both planar (i.e. nanowire on the metallic
substrate) and on-grating configurations. We found that, by varying the
orientation of the nanowire on the grating, only when the nano-cavity was
perpendicular to the ridge direction, an additional peak emerged in the
emission spectrum on the low-energy side of the gain envelope. As a
consequence of the fulfillment of the Bragg condition, such a peak was
attributed to a hybrid mode dominating the mode competition. Simulation
results showed that the hybrid mode could be efficiently waveguided along the
nanowire cavity and supported by localized plasmon polaritons building up at
the raised features ("fences") on top of the metal grating ridges. Moreover,
the hybrid mode was found to experience an extra reflectance of nearly 50%
across the grating periods, in addition to that provided by the nanowire end
facets.

††preprint: AIP/123-QED

Since the theoretical formulation of the SPASER Bergman and Stockman (2003)
and following the first experimental realization Oulton  _et al._ (2008,
2009), nanowire-based plasmonic lasers - mostly designed as a semiconductor-
insulator-metal (SIM) scheme - have gained an ever-increasing scientific
appeal Lu  _et al._ (2012); Chou  _et al._ (2016a); Chu  _et al._ (2022). In
particular, due to the buildup of surface plasmon polaritons (SPPs) at the
interface between the nanowire and the metal film, plasmonic nanolasers can
overcome the diffraction limit via the sub-wavelength confinement of the modal
field Ma  _et al._ (2011). As a consequence, such devices are characterized by
smaller modal volumes Bergman and Stockman (2003); Oulton (2012), a larger
Purcell-enhanced stimulated emission Zhang  _et al._ (2014); Chou  _et al._
(2015) and, thus, a faster lasing dynamics Stockman (2010); Sidiropoulos  _et
al._ (2014), when compared to photonic nanolasers Lu  _et al._ (2017); Wille
_et al._ (2016). However, within the framework of next-generation all-optical
circuitry De Leon  _et al._ (2012); Davis  _et al._ (2017); Fukuda  _et al._
(2020), the tailoring of the emission characteristics and the lowering of the
lasing threshold still represent the major challenges in the field for the
integration of plasmonic nanolasers Xu  _et al._ (2019); Liang  _et al._
(2020). Along with improving the metal quality and smoothening the surface of
the metallic substrate Yu  _et al._ (2017); Chou  _et al._ (2016b); Liao  _et
al._ (2019), researchers have started focusing on the modification of the
plasmonic cavities to boost their lasing performances, either by structuring
the dielectric/metal substrate Bermúdez-Ureña  _et al._ (2017); Chou  _et al._
(2018); Cheng  _et al._ (2018) or by hybridizing the dielectric spacer between
the nanowire and the metal film Li  _et al._ (2019, 2020).
In this work, we made use of a mechanical nano-manipulation approach, to study
the effects of a metal grating on the lasing properties of semiconductor
nanowire-based plasmonic lasers. Namely, the same nanowire could be
manipulated and “switched” among different configurations: from placed on the
dielectric/metal substrate (planar plasmonic configuration) to overlaid on a
metal grating (nanowire-metal grating configuration) at different orientations
with respect to the ridge (or trench) direction. To ensure reliable
comparability among the different configurations, the single-nanowire
experiments compensated for the statistical fluctuations arising from the
wire-to-wire variations in crystal quality, cavity size (length and diameter),
and morphology (if tapered, strained, roughened or cleaved at the end facets):
all parameters that significantly affect the cavity Q factor and the overall
emission characteristics Zimmler  _et al._ (2010); Couteau  _et al._ (2015);
Zapf  _et al._ (2019). To this end, we used Al gratings characterized by
geometrical parameters tailored to the lasing emission wavelength, in order to
externally induce additional distributed feedback to the gain spectrum.
Zinc oxide (ZnO) nanowires were produced by an Au-catalyzed vapor-liquid-solid
(VLS) growth inside a horizontal three-zone tube furnace: further details can
be found in reference Borchers  _et al._ (2006). Afterward, the nanowires were
dry-imprinted from the growth sample onto the dielectric/metal substrates.
Such substrates consisted of a native Al2O3 layer (thickness $t\approx$ 3 nm),
serving as a spacer for enhancing the modal field confinement and for
mitigating the metal losses, on top of a polycrystalline Al layer ($t\approx$
70 nm), which featured the co-deposition of 10 at% Si for decreasing the
surface roughness. Subsequently, these substrates were subdivided into etching
spots that were patterned via electron beam lithography (EBL) for the
fabrication of the metal gratings Kley  _et al._ (2012): for further details
see the Supporting Material, section S1.
The optical characterization of single nanowires in a given configuration was
accomplished by focusing the emission from a frequency-tripled Nd:YAG laser
($\lambda_{\textup{exc}}$ = 355 nm, $f_{\textup{rep}}$ = 100 Hz,
$t_{\textup{pulse}}$ = 7 ns) through a 50x NUV objective (NA = 0.43), down to
a spot size of about 20 $\upmu$m, to encompass the whole nanowire. The
spectral acquisition was carried out by using a spectrometer (Princeton
Instruments SP-2500i) equipped with a 1200 lines/mm grating (blazed at 300 nm)
and connected to a liquid-nitrogen-cooled, front-illuminated CCD camera. The
photoluminescence measurements were performed at a temperature of 150 K inside
a liquid He flow cryostat (Janis ST-500).

Figure 1: (a) SEM image of a “fenced” Al grating, featuring a pitch $p$ = 300
nm and fill factor (i.e. the trench-width/pitch ratio) $FF$ = 0.3, used for
the coupling with a single nanowire. (b) Sketch of the nano-manipulation
performed in air: after the characterization on the bare dielectric/metallic
substrate, the nanowire was moved onto the metal grating, and oriented at a
given angle $\theta$ with respect to the ridge (trench) direction. (c)
Exemplary SEM image of a nanowire used in one of the single-nanowire
experiments, which was eventually oriented at $\theta\approx 60\degree$ via
the nano-manipulation.

We started our characterization by screening the randomly distributed
nanowires with respect to their morphology (heavily tapered, bundled, and
inhomogeneous nanowires were discarded) and proximity to the metal gratings:
suitable nanowires were pre-characterized via low-power PL measurements. The
nanowire manipulation had to be performed in air under the PL microscope
rather than in a scanning electron microscope (SEM), in order to avoid the
modification of the optical properties induced by the electron beam Zapf
(2021); Lähnemann  _et al._ (2016). Therefore, SEM investigations were
performed after all the optical characterizations of each cavity
configuration.
Figure 1(a) exemplarily shows a top-view SEM image of an Al grating with pitch
$p$ = 300 nm, along with the additional “fence-like” structures on top of the
metallic ridges. Such structures originated from the pattern transfer
processes using ion beam etching (for more details see Supporting Material,
section S1). The nanowires were overlaid onto such structures via nano-
manipulation, as schematically shown in Figure 1(b), which was accomplished by
exploiting a single piezo-controlled nano-manipulator (Kleindiek GmbH)
equipped with a tungsten tip (further details can be found in the Supporting
Material, section S2). Single nanowires were accurately placed on the metal
grating and oriented at a certain angle $\theta$ with respect to the ridge
direction with an uncertainty of $\approx\pm 5\degree$: the result can be seen
in the SEM image of Figure 1(c), which exemplarily shows a nanowire used in
the PL experiments.

Figure 2: (a) PL spectra of the planar (top panel) and nanowire-grating
(bottom panel) configurations acquired for different pump power densities at
150 K; the NW was oriented perpendicular to the ridge direction $\theta\approx
90\degree$. (b) Corresponding light-in-light-out curves (circles) of the
emission intensity and spectral width trends (triangles) of the emission for
the planar (top panel) and nanowire-grating (bottom panel) configurations, as
a function of the excitation power density. The dashed lines indicate the
respective lasing thresholds estimated by fitting the experimental data points
with an adapted multimode lasing model. (c) SEM cross-section of an Al grating
with a nominal pitch $p$ = 300 nm and fill factor $FF$ = 0.3. The extracted
geometrical parameters for the numerical simulations, within the fabrication
tolerances and measurement errors, were: fence separation $s$ = 170 nm, fence
width $b$ = 20 nm, and trench width $w$ = 90 nm. (d) Simulated reflectance of
the “fenced” Al grating as a function of the wavelength and nanowire index.
The intersection point of the dashed white lines indicates the reflectance at
$\lambda$ = 378 nm for the corresponding material index value
($n_{\textup{ZnO}}$ = 2.26).

Figure 2(a) shows the photoluminescence spectra acquired as a function of the
excitation power density, with the pump polarization kept perpendicular to the
nanowire axis, at a temperature of 150 K, which turned out to be the upper-
temperature limit for ensuring good device stability throughout the whole
experiment. The spectra were recorded for a ZnO nanowire with length
$L\approx$ 12 $\upmu$m and average diameter $\langle d\rangle\approx$ 200 nm,
first lying planar on the plasmonic substrate (top panel), and subsequently
moved onto the Al grating with an orientation of $\theta\approx 90\degree$
with respect to the ridge direction (bottom panel). The diameter was averaged
over the width variations ($\approx\pm 20-30$ nm) along the nanowire, while
the widening, usually occurring at one of the end facets as a consequence of
tapering, was not considered in the averaged sum. The PL spectrum of the
nanowire on the grating clearly shows the emergence of an additional peak on
the low-energy side of the gain envelope at $\lambda$ = 378 nm, when compared
to the case of the planar configuration.
The geometrical parameters were extracted from the SEM cross-section, shown in
Figure 2(c), and the structure was reproduced in the commercial FDTD software
suite Lumerical Ansys Optics (2019) to calculate the corresponding reflectance
spectrum, shown in Figure 2(d), as a function of wavelength and refractive
index. We found that this structure supports reflectance resonances within the
ZnO gain spectrum around 380 nm, where we assumed the refractive index of the
nanowire to be $n_{\textup{ZnO}}$ = 2.26, as marked by the white dashed lines
in Figure 2(d). Thus, the emergence of the additional peak can be explained
based on the reflectivity characterizing different waveguided modes: indeed,
we observed a sharp reflectance peak for a mode polarized perpendicularly to
the grating plane. Concerning this mode, we varied the material index until
the reflectance peak coincided spectrally with that observed in the
experiment, and verified that the Bragg condition:

$2\times\frac{2\pi}{\lambda}n_{eff}=m\frac{2\pi}{p},$ (1)

\- where $\lambda$ is the peak wavelength, $n_{\textup{eff}}$ the effective
refractive index of the mode, and $p$ the grating pitch -, is fulfilled for a
grating order of $m$ = 3.
Figure 2(b) shows the corresponding dependence of the photoluminescence
intensity on the pump intensity together with that of the full-width-half-
maximum (FWHM) of the dominant lasing peak of the nanowire, on the planar
substrate (top panel) and on the metal grating (bottom panel). The light-in-
light-out curves demonstrate solid device robustness over a large pumping
range (about 2 orders of magnitude), with no sign of optical degradation,
which could affect the comparison between the two different configurations.
These curves were fitted using an adapted multimode lasing model Casperson
(1975), which returned as fitting parameters the lasing threshold
$I_{\textup{th}}$ and the spontaneous emission factor $\beta$. We found
$I_{\textup{th}}=27\pm 5$ kW/cm2 and $\beta=0.06\pm 0.01$ for the planar
configuration and $I_{\textup{th}}=39\pm 8$ kW/cm2 and $\beta=0.10\pm 0.02$
for the on-grating configuration, respectively. These fitted values for the
lasing threshold also coincide with the characteristic kink observed in the
respective “S-shaped” curves, marking the transition from the amplified
spontaneous emission to the lasing regime. It is apparent from Figure 2(b)
that the values for the two configurations are comparable within the
experimental errors. From the previously-reported calculated reflectance
values, which turned out to be higher than the ones characteristic of the
pristine nanowire end-facet cavity Zimmler  _et al._ (2010), we expected a
lower threshold for the lasing of the nanowire on the grating. However, one
should take into account that the peak associated with the mode dominating the
lasing action, appears at the low-energy side at $\lambda\approx$ 378 nm,
where the gain of ZnO is significantly lower. Furthermore, we attribute the
comparable threshold values also to the fabrication imperfections, also
visible in Figure 1(b), resulting in additional losses at the interfaces
between the nanowire and the grating Wright  _et al._ (2014).

Figure 3: (a) PL spectra of a single ZnO nanowire placed on an Al grating with
$p$ = 300 nm and $FF$ = 0.3 and switched between two different orientations
($\theta\approx 90\degree$ and $\theta\approx 60\degree$), then moved onto the
metallic substrate, as illustrated by the sketches in the insets. Simulated
field distributions of the lasing mode for the experimentally investigated
geometry across the nanowire at the (b) fence and (c) free-standing positions.
Field intensity distributions of the lasing cavity mode - cut perpendicular to
the grating - along the nanowire axis in (d) zoomed in and (e) full view. The
scale bars on the x and y axes are chosen differently for display purposes,
while the white solid lines delimit the nanowire diameter. The intensity
increase in the mode distribution at the center of the nanowire, when compared
to the two end facets, visible in (e), can be likely attributed to the field
enhancement given by the metallic fences acting as nano-resonators. As
evidenced in (d), 3 maxima appear across one grating period, further
corroborating our hypothesis of a third-order distributed feedback mechanism.

To further verify the effects of the extra feedback provided by the metal
grating, we conducted other single-nanowire experiments (see also Supporting
Material, section S3), featuring an additional re-orientation of the ZnO
nanowire on top of the Al grating. Figure 3(a) shows the lasing spectra of a
nanowire, with cavity dimensions $L\approx 13$ $\upmu$m and $\langle
d\rangle\approx$ 250 nm, which was first measured on the Al grating with a
ridge orientation of $\theta\approx 90\degree$, then re-oriented at
$\theta\approx 60\degree$ and, ultimately, moved off the grating onto the
planar substrate. This experiment further corroborates the fact that the
appearing emission peak in the low-energy side of the gain spectrum, is caused
by the external resonator, which provides additional distributed feedback to
the nanowire cavity. Indeed, compared to the other single-nanowire experiment,
not only the emission spectrum changes between the planar and grating
configurations, but also the further re-orientation at $\theta\approx
60\degree$ leads to the disappearance of the additional peak. This proves that
only for an efficient coupling of the whole nanowire cavity with the fence-
like structures, i.e. for an effective pitch matching the coupling length
between the nanowire and the metal grating, a “grating-supported” mode starts
to dominate the lasing action in the high power regime.

To quantify and obtain more detailed insights into the feedback provided by
the coupling of the ZnO nanowires with the metal grating, we performed active
lasing simulations in Lumerical (for details see Supporting Material, section
S4) by making use of the same above-mentioned parameters. The results show
that the excited modes, polarized perpendicular to the substrate plane, attain
a field enhancement at the edges of the metallic fences, thanks to the buildup
of localized plasmon polaritons Chou  _et al._ (2018). This, in turn, leads to
a large modal mismatch between the waveguided modes for the free-standing
parts of the nanowire, shown in Figure 3(b), and those in contact with the
metal, shown in Figure 3(c). Hence, this results in reflectance values of
about 50%, which add up to the overall end facet reflectivity even for a small
number of periods ($N$ = 40), as already shown in Figure 2(d). The resulting
increase of the quality factor of the cavity mode, thus, causes a considerable
line narrowing compared with the nanowire cavity when lying on the planar
dielectric/metal substrate. Indeed, the mode dominating the lasing action is
given by the fenced grating rather than by the Fabry-Pérot modes building up
between the end facets of the nanowire. Furthermore, one can see that the
field intensity distribution along the nanowire is highest in the center of
the nanowire, due to the larger and symmetric feedback of the external grating
cavity, conversely to the poorer confinement toward the end facets, where one
would expect most of the radiation to be coupled out in the absence of the
additional grating-supported resonator. This also suggests that an optimal
coupling between the nanowire and the grating could be attained only for an
orientation of the nanowire perpendicular to the metal grating. Ultimately,
due to the action of the grating, the rather broad-band plasmonic lasing
emission was converted into narrow-line quasi-single mode lasing in the high-
power regime, where the device proved to be robust and stable.
In summary, the desired modification of nanowire-based lasers represents a
promising tool for the realization of next-generation plasmonic circuits. In
this work, we demonstrated that the coupling of ZnO nanowires with “fenced” Al
gratings, enables the tailoring of their lasing emission. Namely, for an
orientation of the nanowire perpendicular to the grating ridges (trenches),
the experimental and theoretical results showed that the lasing action is
dominated by a hybrid mode which is supported by the effective coupling with
the localized plasmon polaritons at the fence edges. This enables a
distributed-feedback-like behavior for the grating-supported resonator, which
provides additional feedback to the nanowire cavity, with reflectance values
as high as 50%. With this work, we intend to foster efforts toward the
optimization of tailoring schemes for such a class of nanolasers, given the
ultimate goal of attaining a stable room-temperature operation, which, as of
now, seems to be still highly hindered by the strong plasmonic losses.
Nonetheless, the compactness of the grating structures paves the way to the
integration of nanowire–metal grating arrays, for which, further mode
selection schemes can be accordingly tailored and realized.

See the supplementary material for more information about the fabrication
process, the nanowire manipulation procedure, additional measurements and
simulation details.

## AUTHOR DECLARATIONS

### Acknowledgments

The authors acknowledge the financial support via the projects C5, A2, and Z3
in the CRC 1375 “NOA–Nonlinear optics down to atomic scales” funded by the
Deutsche Forschungsgemeinschaft (DFG) under project-ID 398816777.

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

Francesco Vitale: conceptualization (equal); investigation – experimental
(lead); writing – original draft (lead). Daniel Repp: formal analysis,
software – simulations (lead); writing – review & editing (equal). Thomas
Siefke: resources – sample fabrication (lead); writing – review & editing
(equal). Uwe Zeitner: writing – review & editing (equal). Ulf Peschel: writing
– review & editing (equal). Thomas Pertsch: funding acquisition (equal);
writing – review & editing (equal). Carsten Ronning: funding acquisition
(equal); conceptualization (equal); project administration (lead); supervision
(lead); writing – review & editing (lead).

## References

  * Bergman and Stockman (2003) D. J. Bergman and M. I. Stockman, Phys. Rev. Lett. 90, 027402 (2003).
  * Oulton  _et al._ (2008) R. F. Oulton, V. J. Sorger, D. Genov, D. Pile, and X. Zhang, Nat. Photon. 2, 496 (2008).
  * Oulton  _et al._ (2009) R. F. Oulton, V. J. Sorger, T. Zentgraf, R.-M. Ma, C. Gladden, L. Dai, G. Bartal, and X. Zhang, Nature 461, 629 (2009).
  * Lu  _et al._ (2012) Y.-J. Lu, J. Kim, H.-Y. Chen, C. Wu, N. Dabidian, C. E. Sanders, C.-Y. Wang, M.-Y. Lu, B.-H. Li, X. Qiu, _et al._ , Science 337, 450 (2012).
  * Chou  _et al._ (2016a) Y.-H. Chou, Y.-M. Wu, K.-B. Hong, B.-T. Chou, J.-H. Shih, Y.-C. Chung, P.-Y. Chen, T.-R. Lin, C.-C. Lin, S.-D. Lin, _et al._ , Nano Lett. 16, 3179 (2016a).
  * Chu  _et al._ (2022) K. Chu, F. Ma, X. Zhu, X. Jia, Z. Huang, K. Dong, J. Sun, K. Liu, P. Jin, Z. Wang, _et al._ , J. Phys. D 55, 215104 (2022).
  * Ma  _et al._ (2011) R.-M. Ma, R. F. Oulton, V. J. Sorger, G. Bartal, and X. Zhang, Nat. Mater 10, 110 (2011).
  * Oulton (2012) R. F. Oulton, Mater. Today Commun. 15, 26 (2012).
  * Zhang  _et al._ (2014) Q. Zhang, G. Li, X. Liu, F. Qian, Y. Li, T. C. Sum, C. M. Lieber, and Q. Xiong, Nat. Comm. 5, 1 (2014).
  * Chou  _et al._ (2015) Y.-H. Chou, B.-T. Chou, C.-K. Chiang, Y.-Y. Lai, C.-T. Yang, H. Li, T.-R. Lin, C.-C. Lin, H.-C. Kuo, S.-C. Wang, _et al._ , ACS Nano 9, 3978 (2015).
  * Stockman (2010) M. I. Stockman, J. Opt. 12, 024004 (2010).
  * Sidiropoulos  _et al._ (2014) T. P. Sidiropoulos, R. Röder, S. Geburt, O. Hess, S. A. Maier, C. Ronning, and R. F. Oulton, Nat. Phys. 10, 870 (2014).
  * Lu  _et al._ (2017) J. Lu, M. Jiang, M. Wei, C. Xu, S. Wang, Z. Zhu, F. Qin, Z. Shi, and C. Pan, ACS Photonics 4, 2419 (2017).
  * Wille  _et al._ (2016) M. Wille, C. Sturm, T. Michalsky, R. Röder, C. Ronning, R. Schmidt-Grund, and M. Grundmann, Nanotechnology 27, 225702 (2016).
  * De Leon  _et al._ (2012) N. P. De Leon, M. D. Lukin, and H. Park, IEEE J. Sel. Top. Quantum Electron. 18, 1781 (2012).
  * Davis  _et al._ (2017) T. J. Davis, D. E. Gómez, and A. Roberts, Nanophotonics 6, 543 (2017).
  * Fukuda  _et al._ (2020) M. Fukuda, S. Okahisa, Y. Tonooka, M. Ota, T. Aihara, and Y. Ishikawa, IEEE Access 8, 142495 (2020).
  * Xu  _et al._ (2019) L. Xu, F. Li, Y. Liu, F. Yao, and S. Liu, Appl. Sci. 9, 861 (2019).
  * Liang  _et al._ (2020) Y. Liang, C. Li, Y.-Z. Huang, and Q. Zhang, ACS Nano 14, 14375 (2020).
  * Yu  _et al._ (2017) H. Yu, T. P. Sidiropoulos, W. Liu, C. Ronning, P. K. Petrov, S.-H. Oh, S. A. Maier, P. Jin, and R. F. Oulton, Adv. Opt. Mat. 5, 1600856 (2017).
  * Chou  _et al._ (2016b) B.-T. Chou, Y.-H. Chou, Y.-M. Wu, Y.-C. Chung, W.-J. Hsueh, S.-W. Lin, T.-C. Lu, T.-R. Lin, and S.-D. Lin, Sci. Rep. 6, 1 (2016b).
  * Liao  _et al._ (2019) Y.-J. Liao, C.-W. Cheng, B.-H. Wu, C.-Y. Wang, C.-Y. Chen, S. Gwo, and L.-J. Chen, RSC Adv. 9, 13600 (2019).
  * Bermúdez-Ureña  _et al._ (2017) E. Bermúdez-Ureña, G. Tutuncuoglu, J. Cuerda, C. L. Smith, J. Bravo-Abad, S. I. Bozhevolnyi, A. Fontcuberta i Morral, F. J. García-Vidal, and R. Quidant, Nano Lett. 17, 747 (2017).
  * Chou  _et al._ (2018) Y.-H. Chou, K.-B. Hong, C.-T. Chang, T.-C. Chang, Z.-T. Huang, P.-J. Cheng, J.-H. Yang, M.-H. Lin, T.-R. Lin, K.-P. Chen, _et al._ , Nano Lett. 18, 747 (2018).
  * Cheng  _et al._ (2018) C.-W. Cheng, Y.-J. Liao, C.-Y. Liu, B.-H. Wu, S. S. Raja, C.-Y. Wang, X. Li, C.-K. Shih, L.-J. Chen, and S. Gwo, ACS Photonics 5, 2624 (2018).
  * Li  _et al._ (2019) H. Li, J.-H. Li, K.-B. Hong, M.-W. Yu, Y.-C. Chung, C.-Y. Hsu, J.-H. Yang, C.-W. Cheng, Z.-T. Huang, K.-P. Chen, _et al._ , Nano Lett. 19, 5017 (2019).
  * Li  _et al._ (2020) H. Li, Z.-T. Huang, K.-B. Hong, C.-Y. Hsu, J.-W. Chen, C.-W. Cheng, K.-P. Chen, T.-R. Lin, S.-J. Gwo, and T.-C. Lu, Adv. Sci. 7, 2001823 (2020).
  * Zimmler  _et al._ (2010) M. A. Zimmler, F. Capasso, S. Müller, and C. Ronning, Semicond. Sci. Technol. 25, 024001 (2010).
  * Couteau  _et al._ (2015) C. Couteau, A. Larrue, C. Wilhelm, and C. Soci, Nanophotonics 4, 90 (2015).
  * Zapf  _et al._ (2019) M. Zapf, T. Sidiropoulos, and R. Röder, Adv. Opt. Mat. 7, 1900504 (2019).
  * Borchers  _et al._ (2006) C. Borchers, S. Müller, D. Stichtenoth, D. Schwen, and C. Ronning, J. of Phys. Chem. B 110, 1656 (2006).
  * Kley  _et al._ (2012) E.-B. Kley, H. Schmidt, U. Zeitner, M. Banasch, and B. Schnabel, in  _28th European Mask and Lithography Conference_ , Vol. 8352 (SPIE, 2012) pp. 187–194.
  * Zapf (2021) M. Zapf,  _Performance and stability of semiconductor nanowire devices_ , Ph.D. thesis, Friedrich-Schiller-Universität Jena (2021).
  * Lähnemann  _et al._ (2016) J. Lähnemann, T. Flissikowski, M. Wölz, L. Geelhaar, H. T. Grahn, O. Brandt, and U. Jahn, Nanotechnology 27, 455706 (2016).
  * Ansys Optics (2019) Ansys Optics, “Lumerical Inc. FDTD Solutions v.2019aR2,” https://www.lumerical.com (2019).
  * Casperson (1975) L. W. Casperson, J. of Appl. Phys. 46, 5194 (1975).
  * Wright  _et al._ (2014) J. B. Wright, S. Campione, S. Liu, J. A. Martinez, H. Xu, T. S. Luk, Q. Li, G. T. Wang, B. S. Swartzentruber, L. F. Lester, _et al._ , Appl. Phys. Lett. 104, 041107 (2014).