TITLE: Language around Pointwise convergence v. Uniform convergence
QUESTION [3 upvotes]: I've heard uniform convergence described loosely as convergence "occurring at the same rate at every point." Is that really accurate? If the absolute difference between a member of a sequence of functions and the function to which that sequence converges pointwise can be represented by a function of the point with a finite maximum, then isn't it true that the sequence of functions uniformly convergent AND the rate of convergence varies with the point?

REPLY [9 votes]: Given an $\epsilon>0$ and $x$, we can define $N_x(\epsilon)$ as the smallest $N$ such that $|f_n(x)-f(x)|<\epsilon$ for all $n\geq N.$
In some sense, $N_x(\epsilon)$ is a measure of the rate convergence at $x.$
Uniform convergence means that for any $\epsilon>0,$ $N(\epsilon)=\sup_{x} N_x(\epsilon)<+\infty.$
So this doesn't strictly mean they converge at the same speed, but that there is some speed of convergence which is slower than (or equal to) all of them.


TITLE: Generalized Josephus problem and solutions for higher cases
QUESTION [0 upvotes]: So from the Wikipedia page about this problem, there is a variant of this problem where for $n$ people in a circle, instead of killing the next person, kills the $k$th person instead. (Thus the original problem will have $k=2$.)
Apparently there are solutions for the original $k=2$ which is $f(n)=2(n-2^{\lfloor\log_2n\rfloor})+1$ and for $k=3$ which is
$f(n)=3(n-round(\alpha\cdot(3/2)^m))+2$ with $m$ being the biggest number such that $round(\alpha\cdot(3/2)^m)\le n$ and this constant $\alpha\approx0.8111$.
So I was wondering how is the solution for $k=3$ case found (lack of citation) and is there any way to search for solutions for higher $k$?
Any research/paper links will be appreciated.
Thanks!

REPLY [1 votes]: For the k=3 case, the paper is:
http://citeseerx.ist.psu.edu/viewdoc/download;jsessionid=0EF1E16177A3F71FB825A92AE8245B5A?doi=10.1.1.46.1946&rep=rep1&type=pdf
general problem:
https://www.jstor.org/stable/43686318?casa_token=flGruEZzecEAAAAA%3A-H4adfhYl4mJ7kTdv61n0VclU2xdXE_OcBOk7qK2ewmko-HbXSdCMfJlZpQVGbNZwglHDQo-iEzoSwfdIPcnfZNmG9FfAyNWaZQsK-BrZoU0t6QwZg&seq=1#metadata_info_tab_contents
Another variant where each person has $l$ lives:
https://link.springer.com/article/10.1007/s00224-011-9343-6


TITLE: Why is $T_{\sigma}(\sigma H)$ the kernel of $\pi_{*\sigma} : T_\sigma G \rightarrow T_{\sigma H} G/H$?
QUESTION [2 upvotes]: I'm trying to understand the proof of theorem $3.62$ page 123 from "Foundations of differentiable manifolds and Lie groups" by Frank F. Warner

Let $\eta : G \times M \rightarrow M$ be a transitive action of a Lie group on the manifold $M$ on the left. Let $m_0 \in M$ and $H = \{g \in G : \eta(g,m_0) = m_0 \}$ be the isotropy group at $m_0$. Define a mapping
$$ \tilde \beta : G/H \rightarrow M : gH \mapsto \eta(g,m_0).$$
Then $\tilde \beta $ is a diffeomorphism.

In the proof the author states that the kernel of the differential map of $\pi$ at $\sigma \in G $ is the tangent space to $\sigma H$ at $g$. That is
$$ \ker \pi_{*\sigma} = T_{\sigma}(\sigma H) \subset T_\sigma G.$$
where
$$ \pi_{*\sigma} : T_\sigma G \rightarrow T_{\sigma H} G/H $$
I don't understand why this is the case. I presume that it is obvious since the author doesn't give a detailed explanation.

My (updated) attempt (Using Thorgott's hint.)
From previous results I know that $H$ is closed which means $\sigma H$ is closed. From this I deduce that $\sigma H$ is a Lie subgroup of $G$ which is embedded in $G$.
I can therefore identify the tangent space $T_{\sigma}{\sigma H}$ with a subspace of $T_\sigma G$ by using the differential of the inclusion map.
Let $X_\sigma \in T_\sigma \sigma H$ and $X_{i(\sigma)} = i_{*\sigma} X_{\sigma} \in T_\sigma G.$ Let $f \in C^\infty (G/H)$. Then
$$ \pi_* X_{i(\sigma)} f = X_{\sigma} (f \circ \pi \circ i).$$
But $ f \circ \pi \circ i$ is constant on $\sigma H$ since $f \circ \pi \circ i(\sigma h) = f(\sigma h H) = f (\sigma H)$. Therefore
$$\pi_* X_{i(\sigma)} = 0$$ and $T_{\sigma} \sigma H \subset \ker \pi_{*\sigma}.$
Using ThorGott's hint again I want to show that $\ker \pi_{*\sigma}$ has the same dimension as $T_\sigma \sigma H = \dim H$ to show the equality.
Since $\dim T_{gH} G/H = \dim G - \dim H$ I can conclude from the rank nullity theorem that $\dim \ker \pi_* = \dim H$ only when $\pi_*$ is surjective.
How could I show surjectivity ?

REPLY [1 votes]: Hint: The map $\pi$ is constant on $\sigma H$, which gives one inclusion. Equality then follows by counting dimensions.


TITLE: How to show that smallest value of $x(x-a)$ (if $0\leq x\leq a$) is when $x=\frac{a}{2}$?
QUESTION [0 upvotes]: How to show that smallest value of $x(x-a)$ is when $x=\frac{a}{2}$?
$0\leq x\leq a$

REPLY [1 votes]: We have the expression $x(x-a)$ with $a \gt 0$.
To answer the OP's question using only elementary mathematics, we first calculate
$\tag 1 x(x-a) = \frac{-a^2}{4} \; \; \text{ when } \; x =\frac{a}{2}$
For any $x \in \mathbb R$, we can write $x = \lambda a$ for some $\lambda\in \mathbb R$. So,
$\quad x(x-a) =  \lambda a (\lambda a - a) \gt \frac{-a^2}{4} \; \; \text{ iff}$
$\quad \lambda^2 a^2 -\lambda a^2 \gt \frac{-a^2}{4} \; \; \text{ iff}$
$\quad 4 \lambda^2 - 4 \lambda + 1 \gt 0  \; \; \text{ iff}$
$\quad (2 \lambda - 1)^2 \gt 0  \; \; \text{ iff}$
$\quad \lambda \ne \frac{1}{2}$


TITLE: Solve the complex system :
QUESTION [1 upvotes]: Find $x$ and $y$ where :
\begin{cases}
ix-(1+i)y=3 \\
(2+i)x+iy=4 & \text{}
\end{cases}
I have used the rectangular form of complex number by expressing $x=(a,b)$ and $y=(c,d)$ and reaching in a system with 4 equations and 4 unknowns in this form :
\begin{cases}
-b=3+c+-d \\
a=d+c \\
2a=4+d+c \\
2b=-c-a \\
\end{cases}
where we get the right answer $x=(\frac{7}{3}, \frac{-4}{3})$ and $y=(\frac{1}{3},2)$
I would like to know if there is any simpler approach with some more clever manipulations cause my solving method is a bit long.

REPLY [1 votes]: Just express $x$ as function of $y$ from both equations.
Then set the two expressions to be equal. Then find $y$.
By using this approach I get:
$y=1/3 + 2i$
$x=(7-4i)/3$
So it seems your result is correct.

REPLY [1 votes]: Use the Cramer rule, evaluating the $2\times2$ determinants in complex numbers.
$$\Delta=\left|\begin{array}&i&-(1+i)\\2+i&i\end{array}\right|=3i$$
and so on.


text:
\begin{document}
\title[Mass transport equations]{A McLean Theorem for the moduli space of Lie solutions to mass transport
equations. }
\author{Micah Warren}
\address{Department of Mathematics, Princeton University, Princeton NJ, USA}
\email{mww@math.princeton.edu}
\thanks{The author is supported in part by NSF Grant DMS-0901644. }
\date{\today}
\maketitle

\begin{abstract}
On compact manifolds which are not simply connected, we prove the existence of
\textquotedblleft fake\textquotedblright\ solutions to the optimal
transportion problem. These maps preserve volume and arise as the exponential
of a closed 1-form, hence appear geometrically like optimal transport maps.
The set of such solutions forms a manifold with dimension given by the first
Betti number of the manifold. In the process, we prove a Hodge-Helmholtz
decomposition for vector fields. The ideas are motivated by the analogies
between special Lagrangian submanifolds and solutions to optimal transport problems.

\end{abstract}

\section{\bigskip Introduction}

In \cite{KMW}, graphs of solutions to the optimal transportation were shown to
solve a volume maximization problem, using a calibration argument. With the
appropriate metric on the product space $M\times\bar{M}$, maximality is
equivalent to the vanishing of certain differential forms along the graph of
the optimal map. \ In this note, we discuss the converse and see that, at
least in the smooth case, topology allows for maximizers of the volume problem
which do not arise as solutions to the optimal transportation problem, despite
locally having the same geometric properties. \ These maximizers are special
Lagrangian in the sense of Hitchin \cite{H} and Mealy \cite{Me}, a pseudo
Riemannian analogue of the special Lagrangian geometry of Harvey and Lawson
\cite{HL}. \

In the case when the cost is given by Riemannian distance squared, Delano\"{e}
\cite{D} introduced the notion of ``Lie solutions of Riemannian transport
equations" (see also [L].) \ The graphs of Delano\"{e}'s solutions are
maximizers of the volume maximization problem discussed in \cite{KMW}.

We recall the McLean Theorem \cite[Theorem 3.6]{McL} \cite[Theorem 3.21]{Ma}

\begin{theorem}
[McClean]\ Suppose $L$ is an smooth embedded special Lagrangian submanifold of
a Calabi-Yau\ manifold. \ The moduli space $\mathcal{M}$ of special Lagrangian
submanifolds near $L$ is a manifold of dimension $b_{1}(L).$ \ The tangent
space to $\mathcal{M}$ is identified with the harmonic $1$-forms on $L,$ which
has a naturally induced $L^{2}~$metric.
\end{theorem}

\bigskip The seminal paper of Harvey and Lawson shows that special Lagrangian
submanifolds are minimal submanifolds in Calabi-Yau manifolds, and that a
manifold is special Lagrangian if and only if the K\"{a}hler form and a
certain $n$-form vanishes along the submanifold. \ Hitchin \cite{H} analyzed
the metric on the moduli space in McLean's theorem and showed that it arises
via a Lagrangian embedding into a pseudo-Riemannian space. Hitchin's notion of
special Lagrangian, that the K\"{a}hler form and certain combinations of
$n$-forms vanish, describes Lie solutions of the mass transport problems. \

The optimal transportation problem is the following. \ Given probability
volume forms $\rho$ and $\bar{\rho}$ on manifolds $M$ and $\bar{M}$, and a
cost function $c:$ $M\times\bar{M}\rightarrow

\mathbb{R}

,$ find a map $T:M\rightarrow\bar{M}$ which minimizes a cost integral,
\[
\int_{M}c(x,T(x))d\rho
\]
among all maps $T$ which preserve the volume, i.e.
\begin{equation}
T_{\ast}\rho=\bar{\rho}.\ \label{meas}
\end{equation}
The work of Brenier \cite{B} and McCann \cite{McC} shows that given standard
conditions on the cost function, the unique solution will be the map
satisfying (\ref{meas}) and arising as the cost exponential of the gradient of
a potential function $u:$
\begin{equation}
T(x)=c\text{-}\exp_{x}du(x). \label{cexpdu}
\end{equation}
\ A local version of (\ref{meas}), namely,
\begin{equation}
T^{\ast}\bar{\rho}=\rho\ \label{pullback}
\end{equation}
is equivalent to the form
\[
\bar{\rho}d\bar{x}-\rho dx
\]
vanishing on the graph $(x,T(x))\subset M\times\bar{M},$ and (\ref{cexpdu}
)\ implies the vanishing of a K\"{a}hler form defined in \cite{KM}. \ The more
general problem we attack here is to find maps which locally solve the optimal
transport problem, that is, satisfy (\ref{pullback})\ and arise as the
exponential of a closed form
\begin{equation}
T(x)=c\text{-}\exp_{x}\left(  \eta(x)\right)  ,\text{ \ }d\eta=0. \label{Lie}
\end{equation}
\ \

\ The result is the following. \

\begin{theorem}
\label{main thm} Suppose $M,\bar{M}$ are compact manifolds with nowhere
vanishing smooth densities $\rho$ and $\bar{\rho}$. \ Let $c$ be a continuous
cost function on $M\times\bar{M},$ which away from a cut locus $\mathcal{C}$,
is smooth and satisfies local and global twist assumptions (A2) and (A1). \ If
a diffeomorphism $T$ is a smooth Lie solution of the mass transport equation
(i.e. satisfying (\ref{pullback}) and (\ref{Lie}))\ avoiding the cut locus,
then there is a moduli space of smooth Lie solutions near $T$ which is a
smooth manifold of dimension $b_{1}(M).$ \ \
\end{theorem}

Remark:\ \ There are cost functions such that even for some smooth densities,
the optimal solutions will not be smooth, as was demonstrated by Loeper [Lo],
when the (A3) assumption on the cost is not satisfied. \ \ Also, in the
absence of the (A3) assumption, local $c$-convexity does not imply global, so
even on simply connected domains, Lie solutions are not necessarily solutions
to the optimal transport problem.

The core of our result is a Hodge-Helmholtz type decomposition holding at Lie
solutions, which decomposes deformations of maps into those which preserve
property (\ref{pullback}), those which preserve property (\ref{Lie}) and the
harmonic deformations, which preserve both.  The metric used is the one
induced by the linearized operator to the optimal transport equation,
multiplied by an explicit conformal factor. This metric and concomitant
Laplacian are defined in section 2, and their Hodge-Helmholtz properties are
shown in section 3. \ In section 4 we give a possible politico-economic
application of solutions.

\begin{acknowledgement}
The author is very grateful for advice from Phillipe Delano\"{e}.
\end{acknowledgement}

\section{Preliminaries}

For $n$-dimensional manifolds $M$ and $\bar{M}$, let $c:N\subset M\times
\bar{M}\rightarrow\mathbf{R}$ be a continuous cost function which is smooth
almost everywhere, except on a set $\mathcal{C}\mathfrak{=}M\times\bar{M}-N$
which we call the \textquotedblleft cut locus." \ (The reason for this
terminology is clear if we use the distance squared function on a manifold as
the cost function.) \ \ Let $D\bar{D}c$ be the $n\times n$ matrix\ given by
(in local coordinates)
\[
\left(  D\bar{D}c\right)  _{i\bar{j}}(x,\bar{x})=\frac{\partial^{2}}{\partial
x^{i}\partial\bar{x}^{j}}c(x,\bar{x}).
\]
On $N$ we will require that
\begin{equation}
\det\left(  D\bar{D}c\right)  \neq0 \tag{A2}
\end{equation}
which is a local version of the standard \emph{twist }condition: \ For each
$x$ and $\bar{x}$ respectively, \
\begin{equation}
\bar{x}\rightarrow Dc(x,\bar{x}) \tag{A1}
\end{equation}
is invertible, with an inverse depending continuously on $x$. \

To be clear with our conventions, we recall the following Kantorovich problem:
If
\[
J(u,v)=\int_{M}(-u)d\rho+\int_{\bar{M}}vd\bar{\rho},
\]
the problem is to maximize $J$ over all $-u(x)+v(\bar{x})\leq c(x,\bar{x}).$
\ One also considers a dual problem:\ If
\[
I(\pi)=\int_{M\times\bar{M}}c(x,y)d\pi,
\]
find the minimum of $I$ over all measures $\pi$ on the product space
$M\times\bar{M}$ which have marginals $\rho$ and $\bar{\rho}.$ It is well
known (cf \cite{V}) that
\[
\sup_{-u(x)+v(\bar{x})\leq c(x,\bar{x})}J(u,v)=\inf_{\pi\epsilon\Pi(\rho
,\bar{\rho})}I(\pi).
\]
With this setup in mind we can derive the optimal map $T$ from $u$ as
follows:\ \ Suppose $\left(  x_{0},\bar{x}_{0}\right)  $ is a point where the
equality $-u(x_{0})+v(\bar{x}_{0})=c(x_{0},\bar{x}_{0})$ occurs. The function
\[
z_{\bar{x}_{0}}(x)=c(x,\bar{x}_{0})+u(x)\text{ }
\]
must have a minimum at $x_{0}.$\ \ Then define the cost exponential
$T(x_{0},du)=\bar{x}_{0}$. \ If differentiable, from the fact that $z_{\bar
{x}_{0}}(x)$ is at a minimum we have
\begin{equation}
u_{i}+c_{i}(x,T(x,du))=0\label{eq:contact}
\end{equation}
where $c_{i}$ refers to differentiation in the first variable. \ (One can
check that $T(x,u)=-\exp_{x}\nabla u$ when $c(x,\bar{x})=d^{2}(x,\bar{x})/2$).
\ This only depends locally on the function $u,$ and clearly requires that
$Du$ stay inside the range of $Dc(x,\cdot).$ \ \ The elliptic optimal
transportation equation can be derived by taking another derivative and then a
determinant:
\begin{align}
u_{ij}+c_{ij}+c_{i\bar{s}}T_{j}^{\bar{s}} &  =0\label{id0}\\
\det(u_{ij}+c_{ij}) &  =\det(-c_{i\bar{s}}T_{j}^{\bar{s}})=\det(-c_{i\bar{s}
})\frac{\rho(x)}{\bar{\rho}(T(x))}\label{eq:elliptic}
\end{align}
using the fact that $T$ locally pushes $\rho$ forward to $\bar{\rho},$ so
satisfies
\begin{equation}
\bar{\rho}(T(x))\det DT=\rho(x).\label{mp}
\end{equation}
Here and in the sequel we use the following conventions:\ Indices $i,j,k$ etc
will be coordinates on $M$ while $\bar{p},\bar{s},\bar{v},$ etc will be
coordinates on $\bar{M}.$The variable $t$ will be reserved for variations.$\ $
\ \ We assume on any local product chart that $c_{i\bar{s}}$ is negative
definite, and let $b_{i\bar{s}}=-c_{i\bar{s},}$ with $b^{\bar{s}i\text{ }}
$\ its inverse. \ The tangent space is denoted $\mathfrak{T}M.$ \ We use
$\ w_{ij}(x)=u_{ij}(x)+c_{ij}(x,T(x)).$ \ \ To be clear, we will use
$\frac{\partial}{\partial x^{k}}w_{ij}$ to denote derivatives of $w$,
otherwise a subscript will denote differentiation. \ We note that most of the
quantities we deal with are coordinate invariant. \ In particular $w_{ij}$
actually is an honest tensor. \ For future use, we make note of the following
restatement of the above identity (\ref{id0})
\begin{align}
w_{ij} &  =b_{i\bar{s}}T_{j}^{\bar{s}}\label{id1}\\
w^{kj}T_{j}^{\bar{s}} &  =b^{\bar{s}k}.\nonumber
\end{align}

\subsection{ Linearizing the elliptic equation.}

In order to linearize (\ref{eq:elliptic}) take a variation,
\begin{equation}
u(x)+tv(x).\label{variation}
\end{equation}
First, we insert (\ref{variation}) into (\ref{eq:contact}) and differentiate
to obtain
\[
v_{i}-b_{i\bar{s}}(x,T(x,du))D_{t}T^{\bar{s}}=0
\]
hence
\[
D_{t}T^{\bar{s}}=T_{t}^{\bar{s}}=v_{i}b^{\bar{s}i}.
\]
Now, taking a logarithm of (\ref{eq:elliptic})
\begin{align*}
F(x,Du,D^{2}u) &  =\\
&  \ln\det(u_{ij}+c_{ij}(x,T(x))-\ln\det(b_{i\bar{s}}(x,T(x))-\ln\rho
(x)+\ln\bar{\rho}(T(x))
\end{align*}
which is linearized
\begin{align}
Lv &  =\frac{d}{dt}F(u+tv)=w^{ij}(v\,_{ij}+c_{ij\bar{s}}T_{t}^{\bar{s}
})-b^{\bar{s}i}b_{\bar{s}i\bar{p}}T_{t}^{\bar{p}}+\left(  \ln\bar{\rho
}\right)  _{\bar{s}}T_{t}^{\bar{s}}\label{eq:linearized}\\
&  =w^{ij}v\,_{ij}-w^{ij}b_{ij\bar{s}}b^{\bar{s}k}v_{k}-b^{\bar{s}i}b_{\bar
{s}i\bar{p}}b^{pk}v_{k}+\left(  \ln\bar{\rho}\right)  _{\bar{s}}b^{\bar{s}
k}v_{k}\nonumber
\end{align}
A version of this linearized operator was introduced by Trudinger and Wang
\cite[2.18]{TW}.

\subsection{The KM and modified KM\ metrics and a related Laplace-Beltrami}

In \cite{KM} Kim and McCann considered the following pseudo metric on the
product space $M\times\bar{M}:$
\[
h=\left(
\begin{array}
[c]{cc}
0 & b_{i\bar{s}}\\
b_{i\bar{s}}^{T} & 0
\end{array}
\right)
\]
and symplectic form
\[
\omega=b_{i\bar{s}}dx^{i}\wedge d\bar{x}^{\bar{s}}.
\]
In this metric, the graph $(x,T(x))$ of the cost exponential of any locally
$c$-convex potential $u$ is space-like and Lagrangian. \ The induced metric,
in terms of coordinates on $M$, is given by
\[
(x,T(x))^{\ast}h=u_{ij}+c_{ij}=w_{ij}.
\]
For given mass densities $\rho$ and $\bar{\rho}$, consider the following
metric in \cite{KMW}

\[
h=\frac{1}{2}\left(  \frac{\rho\bar{\rho}}{\det b_{is}}\right)  ^{1/n}\left(
\begin{array}
[c]{cc}
0 & b_{i\bar{s}}\\
b_{i\bar{s}}^{T} & 0
\end{array}
\right)
\]
and the calibrating form
\[
\Omega=\frac{1}{2}(\rho+\bar{\rho}).
\]
The graph of any solution to the optimal transportation problem will be a
calibrated maximal Lagrangian surface with respect to this metric, but the
converse is not true in general. \ Calibrated maximal Lagrangian surfaces are
what we are studying in this paper. \

To begin, we define another metric. \ \ We will show that with respect to this
metric, the tangent space of deformations of calibrated submanifolds coincide
with harmonic $1$-forms, which is the same situation that occurs in McLean's
theorem. \ \ The metric we use is yet another conformal factor of the metric
used by Kim and McCann, differing by a power of the conformal factor from the
metric associated to the calibration in \cite{KMW}. \ \ For a given frame on
$M$, and a solution of the problem $T$ define
\begin{equation}
g_{ij}(x)=w_{ij}(x)\left(  \frac{\rho(x)\bar{\rho}(T(x))}{\det b_{i\bar{s}
}(x,T(x))}\right)  ^{1/(n-2)}.\label{themetric}
\end{equation}
This metric appears naturalyl suited for our problem, as we will use the Hodge
decomposition with respect to this metric in the sequel.

Before proving the next claim, note that at any point $x\in M,$ one has by the
twist condition that for any covector $\eta\in\mathfrak{T}_{x}^{\ast}M$ there
is at most one point $\bar{x}=T(x,\eta)$ in $\bar{M}$ so that $\eta
=-D(x,T(x,\eta))$ \ This defines $T$ whenever $\eta$ is in an appropriate
domain. \

Our first claim is the following.\

\begin{proposition}
{\label{LB}} Let $n\geq3.$ \ Let $T=T(x,\eta)$ be a cost exponential of a
$1$-form locally given by $du.$ \ \ Define
\[
\theta=\ln\det w_{ij}-\ln\rho(x)+\ln\bar{\rho}(x,T(x))-\ln\det b_{i\bar{s}},
\]
and
\[
\lambda=\left(  \frac{\rho(x)\bar{\rho}(T(x))}{\det b_{i\bar{s}}
(x,T(x))}\right)  ^{1/(n-2)},
\]
then with the above metric (\ref{themetric}) on $M,$
\[
Lz=\lambda\left(  \triangle_{g}z+\frac{1}{2}\langle\nabla\theta,\nabla
z\rangle_{g}\right)
\]
where $\triangle_{g}z$ is Laplace-Beltrami operator with respect to $g$ and
$L$ is a the linearized operator defined by (\ref{eq:linearized}) . \
\end{proposition}

\begin{proof}
\ Defining
\[
g_{ij}=\lambda w_{ij}
\]
with
\[
\lambda=\left(  \frac{\rho(x)\bar{\rho}(T(x))}{\det b_{is}(x,T(x))}\right)
^{1/(n-2)}
\]
we compute in local coordinates:
\begin{align*}
\triangle_{g}z &  =\frac{1}{\lambda^{n/2}\sqrt{\det w_{ij}}}(\lambda
^{n/2-1}\sqrt{\det w_{ij}}w^{ij}z_{i})_{j}\\
&  =\frac{1}{\lambda}w^{ij}z_{ij}+\frac{1}{\lambda}(-w^{ik}w^{lj}
\frac{\partial}{\partial x^{j}}w_{kl}+\frac{n-2}{2}w^{ij}\left(  \ln
\lambda\right)  _{j}+\frac{1}{2}w^{kl}\frac{\partial}{\partial x^{j}}
w_{kl}w^{ij})z_{i}.
\end{align*}
Locally, because $\eta=du,$ we can write $w_{ab}=u_{ab}+c_{ab}$, thus
\[
\lambda\triangle_{g}z=w^{ij}z_{ij}+(-w^{ik}w^{lj}(\frac{\partial}{\partial
x^{j}}w_{kl}-\frac{\partial}{\partial x^{k}}w_{lj})-w^{ik}w^{lj}\frac
{\partial}{\partial x^{k}}w_{lj}+\frac{n-2}{2}w^{ij}\left(  \ln\lambda\right)
_{j}+\frac{1}{2}w^{kl}\frac{\partial}{\partial x^{j}}w_{kl}w^{ij})z_{i}
\]
\begin{align}
&  =w^{ij}z_{ij}+\left\{  w^{ik}b_{kl\bar{s}}b^{\bar{s}l}-b_{lj\bar{s}}
b^{\bar{s}i}w^{lj}-\frac{1}{2}w^{ij}\left(  \ln\det w\right)  _{j}+\frac
{n-2}{2}w^{ij}\left(  \ln\lambda\right)  _{j}\right\}  z_{i}\nonumber\\
&  =w^{ij}z_{ij}+\left[
\begin{array}
[c]{c}
w^{ik}b_{kl\bar{s}}b^{\bar{s}l}-b_{lj\bar{s}}b^{\bar{s}i}w^{lj}\\
-\frac{1}{2}w^{ij}\left(  \theta_{j}-2\frac{\bar{\rho}_{\bar{s}}T_{j}^{\bar
{s}}}{\bar{\rho}}+2b^{\bar{s}k}(b_{k\bar{s}j}+b_{k\bar{s}\bar{p}}T_{j}
^{\bar{p}})\right)
\end{array}
\right]  z_{i}\label{line3}\\
&  =w^{ij}z_{ij}+\left\{  -b_{lj\bar{s}}b^{\bar{s}i}w^{lj}+\frac{\bar{\rho
}_{\bar{s}}}{\bar{\rho}}b^{\bar{s}i}-\frac{1}{2}w^{ij}\theta_{j}-b^{\bar{s}
k}b_{k\bar{s}\bar{p}}b^{\bar{p}i})\right\}  z_{i}\nonumber\\
&  =Lz-\frac{1}{2}\lambda\langle\nabla\theta,\nabla z\rangle_{g}\nonumber
\end{align}
where we have used (\ref{id1}) repeatedly. \ In (\ref{line3}) we combined the
following two identities \
\[
\ln\lambda=\frac{1}{(n-2)}\left\{  \ln\rho(x)+\ln\bar{\rho}(x,T(x))-\ln\det
b_{is}\right\}
\]
\[
\ln\det w=\theta+\ln\rho(x)-\ln\bar{\rho}(x,T(x))+\ln\det b_{is}
\]
before taking the derivative. \
\end{proof}

\section{\bigskip Deformations}

Let $T:M\rightarrow\bar{M}$ be a Lie solution of the transport equation. \ As
was used in the previous calculation, this means that $T$ is locally a
solution of the optimal transport equation, or equivalently, $T$ is a cost
exponential of a closed $1$ form, satisfying

\[
\det DT=\frac{\rho(x)}{\bar{\rho}(T(x))}.
\]
Throughout this section we will use the metric $g$ on $M$ defined by
(\ref{themetric}).

Given a $1$-form $\eta$ on $M$ we can define a (vertical)\ deformation vector
field
\[
V=-\eta_{i}c^{\bar{s}i}\partial_{\bar{s}}
\]
where $\partial_{\bar{s}}$ is the coordinate tangent frame for\ $\mathfrak{T}
\bar{M}$ \ \ (to be precise, $V$ is a section of the pullback bundle $T^{\ast
}\mathfrak{T}\bar{M}).$ \ For a vector field in this bundle, define a
deformation of the map $T$ via
\[
T_{V}(x)=\exp_{T(x)}^{\bar{M}}(V(x)).
\]
\ (We have fixed arbitrarily a metric on $\bar{M}$ in order to define
$\exp^{\bar{M}}$, which is no more than a convenience.) \

Along $N=$ $M\times\bar{M}-\mathcal{C}$ we define the $1$-form $\sigma
=c_{i}(x,\bar{x})dx^{i}$ (which differs from the total differential of $c$ by
$c_{\bar{s}}(x,\bar{x})d\bar{x}^{\bar{s}}$ ). \ Then we get the following
exact K\"{a}hler form
\[
\omega=d\sigma=-c_{i\bar{s}}(x,\bar{x})dx^{i}\wedge d\bar{x}^{\bar{s}}
\]
on $N.$

At a solution $T$ we define the map
\[
\Phi:\Lambda^{1}(M)\rightarrow\Lambda^{2}(M)\oplus\Lambda^{0}(M)
\]
via
\[
\Phi(\eta)=((Id\times T_{V})^{\ast}\omega,\ast\left(  T_{V}(x)^{\ast}\bar
{\rho}d\bar{x}-\rho dx\right)  )
\]

Note that the level set $\Phi^{-1}(0,0)$ consists of forms whose corresponding
maps $T_{V}$ are Lie solutions. \ \ Equivalently, the map $(Id\times T_{V})$
is a calibrated Lagrangian submanifold of $M\times\bar{M}$ with respect to the
metric in \cite{KMW}.

\begin{lemma}
\bigskip The image of $\Phi$ lies in exact $2$-forms and coexact $0$-forms. \
\end{lemma}

\begin{proof}
The first factor is a pullback of an exact form. For the second factor,
\[
\int T_{V}(x)^{\ast}\bar{\rho}d\bar{x}-\rho dx=0,
\]
which follows easily from the fact that the diffeomorphism $T$ is simply a
change of integration variables. We are assuming the total mass of both
densities is $1.$
\end{proof}

\begin{lemma}
\label{main lemma} At $0\in C^{k+1,\alpha}(\Lambda^{1}(M))$, the differential
$D\Phi$ satisfies
\[
D\Phi(\eta)=(d\eta,d^{\ast}\eta).
\]
In particular, $D\Phi$ is a topological linear isomorphism
\[
d^{\ast}C^{k+2,\alpha}(\Lambda^{2}(M))\oplus dC^{k+2,\alpha}(\Lambda
^{0}(M))\longrightarrow dC^{k+1,\alpha}(\Lambda^{1}(M))\oplus d^{\ast
}C^{k+1,\alpha}(\Lambda^{1}(M)).
\]

\end{lemma}

\begin{proof}
To compute the derivative, for the first factor of the map $\Phi$ we use the
Lie derivative on the manifold $M\times\bar{M}$ \ (cf \cite[section
2.2.2.]{Ma})
\begin{align*}
\frac{d}{dt}(Id\times T_{V})^{\ast}\omega|_{t=0} &  =

\mathcal{L}

_{\left(  -c^{\bar{s}m}\eta_{m}\partial_{\bar{s}}\right)  }(Id\times
T_{V})^{\ast}\omega\\
&  =(Id\times T_{V})^{\ast}\left[  -c^{\bar{s}m}\eta_{m}\partial_{\bar{s}
}\lrcorner d\omega-d(c^{\bar{s}m}\eta_{m}\partial_{\bar{s}}\lrcorner
\omega)\right]  \\
&  =(Id\times T_{V})^{\ast}\left[  d(-c^{\bar{s}m}\eta_{m}\partial_{\bar{s}
}\lrcorner\omega)\right]  \\
&  =(Id\times T_{V})^{\ast}\left[  d\left(  \eta_{m}dx^{m}\right)  \right]  \\
&  =d\eta.
\end{align*}
For the second factor, to begin we note that because the volume element is
given by
\[
Vol_{g}=\sqrt{\det w_{ij}\left(  \frac{\rho\bar{\rho}}{\det b_{i\bar{s}}
}\right)  ^{n/n-2}}=\sqrt{\frac{\rho}{\bar{\rho}}\det b_{i\bar{s}}\left(
\frac{\rho\bar{\rho}}{\det b_{i\bar{s}}}\right)  ^{n/n-2}}=\rho\lambda,
\]
we have
\[
\ast(T_{v}^{\ast}\bar{\rho}d\bar{x}-\rho dx)=\frac{1}{\rho\lambda}\left(  \det
DT_{V}\bar{\rho}(T)-\rho\right)  .
\]
Differentiating
\begin{align}
&  \left(  \frac{1}{\rho\lambda}\det DT_{V}\bar{\rho}(T)-\frac{1}{\lambda
}\right)  ^{\prime}\nonumber\\
&  =\det DT_{V}\bar{\rho}(T)\frac{1}{\rho\lambda}\left\{  (\ln\det
DT)^{\prime}+\left(  \ln\bar{\rho}\right)  ^{\prime}-\left(  \ln\rho\right)
^{\prime}-\ln\lambda^{\prime}\right\}  +\frac{1}{\lambda}\left(  \ln
\lambda\right)  ^{\prime}.\label{prime}
\end{align}
Noting that at $0$
\[
\det DT=\frac{\rho}{\bar{\rho}(T)}
\]
the expression (\ref{prime}) becomes
\[
\frac{1}{\lambda}\left\{  (\ln\det DT)^{\prime}+\left(  \ln\bar{\rho}\right)
^{\prime}-\left(  \ln\rho\right)  ^{\prime}\right\}  .
\]
Now in particular, differentiating with respect to $t$ we have
\begin{align}
&  \frac{d}{dt}\ast(T_{v}^{\ast}\bar{\rho}d\bar{x}-\rho dx)|_{t=0}=\frac
{1}{\lambda}\left[  \left(  DT^{-1}\right)  _{\bar{s}}^{j}T_{j,t}^{\bar{s}
}+\left(  \ln\bar{\rho}\right)  _{\bar{s}}T_{t}^{\bar{s}}\right]  \nonumber\\
&  =\frac{1}{\lambda}\left[  \left(  DT^{-1}\right)  _{\bar{s}}^{j}\left(
b^{\bar{s}k}\eta_{k}\right)  _{j}+\left(  \ln\bar{\rho}\right)  _{\bar{s}
}b^{\bar{s}k}\eta_{k}\right]  \nonumber\\
&  =\frac{1}{\lambda}\left[  b_{i\bar{s}}w^{ij}b^{\bar{s}k}\eta_{k,j}
-b_{i\bar{s}}w^{ij}b^{\bar{p}k}b^{\bar{s}m}(b_{m\bar{p}j}+b_{m\bar{p}\bar{r}
}b^{\bar{r}l}w_{lj})\eta_{k}+\left(  \ln\bar{\rho}\right)  _{\bar{s}}
b^{\bar{s}k}\eta_{k}\right]  \nonumber\\
&  =\frac{1}{\lambda}\left[  w^{ij}\eta_{i,j}-w^{ij}b^{\bar{p}k}b_{i\bar{p}
j}\eta_{k}-b^{\bar{p}k}b_{i\bar{p}\bar{r}}b^{\bar{r}i}\eta_{k}+\left(  \ln
\bar{\rho}\right)  _{\bar{s}}b^{\bar{s}k}\eta_{k}\right]  \label{delta}\\
&  =d^{\ast}\eta
\end{align}
using (\ref{id1}) repeatedly. \ The\ last equality follows from the proof of
Proposition \ref{LB}
\[
\lambda d^{\ast}\eta=w^{ij}\eta_{i,j}+\left\{  -b_{lj\bar{s}}b^{\bar{s}
i}w^{lj}+\left(  \ln\bar{\rho}\right)  _{\bar{s}}b^{is}-b^{\bar{s}k}
b_{k\bar{s}\bar{p}}b^{\bar{p}i}\right\}  \eta_{i}.
\]
The latter conclusion of the Lemma follows from standard algebra using the
Hodge decomposition
\[
C^{k+1,\alpha}(\Lambda^{1}(M))=\mathcal{H}^{1}\oplus dC^{k+2,\alpha}
(\Lambda^{0}(M))\oplus d^{\ast}C^{k+2,\alpha}(\Lambda^{2}(M)).
\]
\
\end{proof}

\textbf{Proof of Main Theorem when $n \geq3$} . \ \ Lemma\ \ref{main lemma}
shows that at $0$, the smooth map $\Phi$ has surjective differential
\[
D\Phi:C^{k+1,\alpha}(\Lambda^{1}(M))\longrightarrow d^{\ast}C^{k+1,\alpha
}(\Lambda^{1}(M))\oplus dC^{k+1,\alpha}(\Lambda^{1}(M))
\]
onto the product of exact and coexact forms. \ The kernel at $0$ is the
harmonic forms, which splits by the Hodge decomposition. \ It follows by the
Implicit Function Theorem (see [Ma, Thm 2.11]) that for each harmonic form
$\eta$ close enough $0$ there is a unique form $\chi(\eta)$ lying in the
orthogonal complement of the harmonic forms so that $\Phi(\eta+\chi(\eta))=0.$
\ Thus a neighborhood of $0$ in the harmonic forms on $M$ parametrizes the
moduli space near $T.$

\bigskip

\subsection{$n=2.$}

When $n=2$ we are missing the power of the Hodge theorem, which asserts that
the kernel of the operator has the same dimension as the first cohomology
class. \ We must somehow show that the set of solutions of \{(\ref{delta}
)$=0\}$ has this same dimension. \ Forming the elliptic operator
\[
L=d^{\ast}d+\delta^{\ast}\delta
\]
where $\delta:\Lambda^{1}\rightarrow\Lambda^{0}$ is the operator defined by
(\ref{delta}), the Implicit Function Theorem arguments as used above provide
that the dimension of the moduli space near a solution is the dimension of the
kernel of $L.$

\bigskip For a given $M^{2},\bar{M}^{2} $ consider the transportation problem
of finding
\begin{equation}
F\left(  x,\theta\right)  =(T\left(  x,\theta\right)  ,\Theta\left(
x,\theta\right)  ):M^{2}\times S^{1}\rightarrow\bar{M}^{2}\times S^{1}
\label{lieprod}
\end{equation}
minimizing the cost
\[
\tilde{c}[\left(  x,\theta\right)  ,(\bar{x},\bar{\theta})]=c(x,\bar
{x})+dist_{S^{1}}^{2}(\theta,\bar{\theta})
\]
among all maps $(F,\Theta)$ pushing $\rho\wedge d\theta$ forward to $\bar
{\rho}\wedge d\theta$. \

For any Lie solution $T:M^{2}$ $\rightarrow\bar{M}^{2}$ \ it is clear that
\begin{equation}
F\left(  x,\theta\right)  =(T(x),\theta+\theta_{0}) \label{lieF}
\end{equation}
is a Lie solution to the problem (\ref{lieprod}). \ \ By Theorem
\ref{main thm} for $n=3$ there is a space of deformations of solutions, and it
has dimension equal to $1+b_{1}(M).$ \ Suffice then to show that these
deformations decompose into deformations in each factor. \

Let $T$ be a Lie solution on $M$ and let (\ref{lieF}) be a Lie solution on
$M\times S^{1}.$ Let $\eta$ be a $1$-form which defines a tangent vector to
the space of deformations on $M^{2}\times S^{1},$ at the solution $F,$ \ and
write
\[
\eta=\eta_{1}\omega^{1}+\eta_{2}\omega^{2}+\eta_{\theta}d\theta
\]
for some cotangent frame $\omega^{1},\omega^{2}$ for $M^{2}.$
\ Differentiating the equation
\[
d^{\ast}\eta=0
\]
in the $\theta$ direction, \ we have that

\[
w^{ij}\eta_{i,j\theta}+\left\{  -b_{bj\bar{s}}b^{\bar{s}i}w^{bj}+\left(
\ln\bar{\rho}\right)  _{\bar{s}}b^{\bar{s}i}-b^{\bar{s}k}b_{\bar{s}k\bar{p}
}b^{\bar{p}i}\right\}  \eta_{i},_{\theta}=0,
\]
using the fact that the warped product metric does not depend on $\theta.$
Also, because $\eta$ is a closed form, locally we have $\eta_{\theta,ij}
=\eta_{i,j\theta}$ and $\eta_{i},_{\theta}=\eta_{\theta,i}$. \ Thus (the
honest function) $z=\eta_{\theta}$ locally satisfies an elliptic equation of
the form

\[
w^{ij}z_{ij}+A^{i}z_{i}=0
\]
so enjoys a maximum principle. \ We conclude that on the compact manifold
$M^{2}\times S^{1},$ $\ $the $S^{1}$ $\ $component $\eta_{\theta}$ of the
deformation must be constant. \ In particular, the other two components
satisfy
\[
\delta(\eta_{1}\omega^{1}+\eta_{2}\omega^{2})=0
\]
so we conclude that the kernel of $L$ must have dimension $b_{1}(M).$

\section{A politician's optimal transportation problem}

We give the following application of solutions to the above problem. Suppose
you are the political leader of a compact manifold with non simply connected
topology, and you are in charge of solving an optimal transportation problem.
Due to forces in effect before you came into power, you discover that the
distributions which need to be paired are in fact the same (this fact is known
only to your office.) You are aware of the unspeakably high political
\textquotedblleft cost of doing nothing,\textquotedblright\ as well as the
clout of the Transportation Lobby, so proposing the trivial
identity\ transference plan is not an option. Instead, you must come up with a
plan which is nontrivial, but does not appear arbitrary. Deforming according
to a harmonic $1$-form as described in section 3, you can find a
transportation plan which has positive cost, and locally, to each supplier in
simply connected districts, is given by the cost exponential of a potential
function, so appears optimal. \ As long as no one exhibits a set of points for
which the plan is not cyclically monotone (cf \cite[Chapter 5]{V}), it will be
believed that the plan is a good one . \ If the deformation is small, such a
cycle will involve a large set of carefully selected points and will be quite
difficult to exhibit.

config: arxiv


text:
\begin{document}

\maketitle
\thispagestyle{empty}
\pagestyle{empty}

\begin{abstract}
In this paper we propose an identification procedure of a sparse graphical model associated to a Gaussian stationary stochastic process. The identification paradigm exploits the approximation of autoregressive processes through reciprocal processes in order to improve the robustness of the identification algorithm, especially when the order of the autoregressive process becomes large. We  show that the proposed paradigm leads to a regularized, circulant matrix completion problem whose solution only requires computations of the eigenvalues of matrices of dimension equal to the dimension of the process.
\end{abstract}

\begin{IEEEkeywords}
	 Stochastic systems, Identification, Optimization.
\end{IEEEkeywords}

\section{INTRODUCTION}\label{sec:intro}

\IEEEPARstart{I}{n} the last decade miniaturisation led to an ubiquitous pervasiveness of technology. As a consequence, the number of available high-dimensional data is skyrocketing in every scientific and applicative domain. Several methods to deal with problems involving high-dimensional data have been recently proposed in the literature \cite{BDS}. In this paper we are focusing on graphical models, that represent a possible tool to deal with high-dimensionality of the data \cite{Lauritzen}. Graphical representations provide an immediate visual intuition on the data interdependence. The simplest  graphical model is an undirected graph that can be associated with a Gaussian random vector
\cite{candes2012exact,candes2010matrix,ZorzCh}: nodes correspond to the components of the random vector, and there is an edge between two nodes if the corresponding components are conditionally dependent given all the others.
Very often data are given as time-series and can thus be modelled as stochastic processes. Also such processes  can  be represented as graphical models: in \cite{songsiri2010graphical} a maximum likelihood approach has been proposed for graphical model estimation of autoregressive (AR) Gaussian processes, exploiting the fact that conditional independence relations translates in zero entries in the inverse of the power spectral density. Since the conditional independent pairs are not known a priori, a sparsity-inducing regularizer can be introduced in the problem \cite{SongVan}, leading to a graphical model with a \emph{sparse} structure, i.e. with few edges. The reduced number of edges gives a double advantage: the graph gives a clearer representation of the phenomenon we are observing, and the number of parameters that have to be estimated is reduced. Moreover, sparsity makes the identification procedure less subject to overfitting, thus leading to a parsimonious estimated model. In \cite{LindqARMA} the aforementioned paradigm has been extended to ARMA Gaussian processes wherein the moving average (MA) part is introduced by a scalar ``prior'' power spectral density. The MA part can be determined from logarithmic moments (cepstral coefficients) of the spectrum \cite{PICCI_LINDQUIST, BEL01, BEL02, MK85} but it is not clear how to impose such moments together with the constraints on the covariance lags, so that the estimated spectrum reflects the underlying graphical structure.
The paradigm in \cite{SongVan} can be extended to latent-variable graphical models \cite{ChaParWill,ZorzSep}.

The identification of Gaussian ARMA graphical models can be
performed, in principle, by using the method in \cite{SongVan}: indeed, an ARMA process is well approximated by an high-order AR process. The problem, however, is that the optimization procedure involves the inversion and the eigenvalue decomposition of matrices whose dimensions are proportional to the product of the order of the
AR process by the dimension of the data. As a consequence, the procedure becomes numerically less robust when the process is high-dimensional and the AR approximation is sufficiently accurate and hence ``long''.

In this paper we consider the problem of identifying sparse graphical models for Gaussian reciprocal processes defined in the ``discrete circle'' ${\ZZ}_N$ (the group of the integers modulo $N$). The latter constitute a particular class of periodic processes, \cite{CarFerrPav}, \cite{Rec1}, \cite{Rec2}, \cite{Rec3}, \cite{ringh2016multidimensional}, \cite{PICCI_LINDQUIST}.  It is possible to show that a Gaussian AR process of order $n$ can be well approximated by a reciprocal process of the same order $n$ as long as the  period $N$ is sufficiently large \cite{FerrLinAlg}, \cite{PICCI_LINDQUIST}.
We will show that the identification problem involves block-circulant matrices. This is a big numerical advantage because the inversion and the eigenvalue decomposition of such matrices can be performed robustly even if $n$ is large \cite{ringh2015fast}. Accordingly, the proposed  paradigm can represent potentially a robust method for estimating ARMA graphical models.

The paper is organized as follows: In Section \ref{sec:recpr} we recall the fundamental results regarding the identification of reciprocal  processes. In Section \ref{sec:gms} we introduce graphical models for the reciprocal processes. In Section \ref{sec:id} we present the problem of estimating a sparse graphical model for a reciprocal process and in Section \ref{sec:ne} we show how the proposed method behaves in a numerical example. Finally, in Section \ref{sec:concl} we draw the conclusions.

\emph{Notation and background.} We denote by $\NN,\,\ZZ,\,\RR,\,\RR^{p\times q}$ the set of natural, integers, real numbers and $p\times q$ real matrices, respectively. For any sub-interval $(x_1,x_2):=\{x:\,x_1<x<x_2\}$ of an interval $(a,b)\subset\RR$, we denote with $(x_1,x_2)^c$ the complement set of $(x_1,x_2)$ in $(a,b)$. Given a matrix $A\in\RR^{p\times q}$, we denote by $A^\top$ its transpose and by $\ker(A)$ its kernel, while $I_p$ denotes the identity matrix of order $p$. If $A\in\RR^{p\times p}$ is a square matrix, $\tr(A),\,\det(A)$ and $A^{-1}$ denote, respectively, the trace of $A$, the determinant of $A$ and its inverse, while $\text{diag}(A)\in\RR^p$ is the vector whose entries are the diagonal elements of $A$. If $A$ is symmetric, $A>0$ and $A\ge 0$ indicate that it is positive definite or positive semi-definite, respectively. Moreover, we will use $\EE[\cdot]$ to denote the expectation operator.

In this paper $n$ will always represent the order of the original AR process which will be of dimension $m$. We will denote by $N$  the length the period of the reciprocal
process. Thus $n$, $m$ and $N$ are fixed numbers.  For simplicity we  assume that $N\in\NN$ is an even number, however the results in this paper can be easily adapted to the case when $N$ is odd.
Also, we assume $N>2n$.
A central role in this paper will be played by the vector space $\C\subset\RR^{mN\times mN}$ containing all the (real) symmetric, block-circulant matrices
\[
   \vect{C} = \text{circ}\{C_0,C_1,\dots,C_{\frac{N}{2}-1},C_{\frac{N}{2}},C_{\frac{N}{2}-1}^\top,\dots,C_1^\top\},
\]
whose first block-column is composed by the $m\times m$ blocks $C_0,C_1,\dots,C_{\frac{N}{2}-1},C_\frac{N}{2},C_{\frac{N}{2}-1}^\top,\dots,C_1^\top$. For any $\vect{C},\,\vect{D}\in\C$, the inner product on $\C$ is defined by $\left<\vect{C},\vect{D}\right>_\C:= \tr(\vect{C}^\top\vect{D})$. We define the \emph{symbol} of the circulant matrix $\vect{C}\in\C$ as the $m\times m$ pseudo-polynomial
\begin{equation}\label{eq:symb}
\Phi(\zeta) := \sum_{k=0}^{N-1}\,C_k\,\zeta^{-k},\quad \text{ with } \quad C_k = C_{N-k}^\top\text{ for } k>\frac{N}{2},
\end{equation}\noindent
where $\zeta:=e^{i\frac{2\pi}{N}}$ is the $N$-th root of unity. It is useful to recall the following result on block-circulant matrices.
\begin{proposition}\label{prop:diagsymb}
	Let $\vect{C}$ be a block-circulant matrix with symbol $C(\zeta)$ defined by \eqref{eq:symb}. Then
	\begin{equation}
	\vect{C} =
	\vect{F}^*
	\text{diag}\left\{\Phi(\zeta^0),\,\Phi(\zeta^1),\,\dots,\,\Phi(\zeta^{N-1})\right\}
	\vect{F},
	\end{equation}
	where $\vect{F}$ is the $mN\times mN$ (Fourier) unitary matrix
	{\small\begin{equation*}
	\vect{F}=
	\frac{1}{\sqrt{N}}
	\begin{bmatrix}
	\zeta^{-0 \cdot 0}I_m    & \zeta^{-0 \cdot 1}I_m     & \ldots & \zeta^{-0 \cdot (N-1)}I_m    \\
	\zeta^{-1 \cdot 0}I_m     & \zeta^{-1 \cdot 1}I_m     & \ldots & \zeta^{-1 \cdot (N-1)}I_m     \\
	\vdots                   & \vdots                   & \ddots & \vdots                       \\
	\zeta^{-(N-1) \cdot 0}I_m & \zeta^{-(N-1) \cdot 1}I_m & \ldots & \zeta^{-(N-1) \cdot (N-1)}I_m \\
	\end{bmatrix}.
	\end{equation*}}
\end{proposition}
This is a classical result in the scalar case; the general analysis for block-circulant matrices can be found, for instance, in \cite[page 6]{tesi}. We introduce the subspace $\B\subseteq\C$ of symmetric, banded block-circulant $mN\times mN$ matrices of bandwidth $n$, with $N>2n$, containing matrices of the form
\begin{equation}\label{eq:bandmtx}
\vect{B} = \text{circ}\{B_0,B_1,\dots,B_n,0,\dots,0,B_n^\top,\dots,B_1^\top\},
\end{equation}
that inherits the inner product defined on $\C$. Note that, according to definition \eqref{eq:symb}, the symbol of a banded matrix $\vect{B}\in\B$ is
\begin{equation*}
\Psi(\zeta) = \sum_{k=-n}^n\, B_k\,\zeta^{-k}, \qquad B_{-k} = B_k^\top.
\end{equation*}
If $\mathscr{B} := [0,n]$ denotes the set of indexes of the blocks in the banded structure, then the projection operator $\mathsf{P}_\mathscr{B}:\C\to\B$ is defined as
\[
\mathsf{P}_\mathscr{B}(\vect{C}) := \text{circ}\{C_0,C_1,\dots,C_n,0,\dots,0,C_n^\top,\dots,C_1^\top\}.
\]

\section{RECIPROCAL PROCESS IDENTIFICATION}\label{sec:recpr}
Let $\{\vect{y}(k),\,k=1,2,\dots,N\}$, be a zero-mean, $m$-dimensional Gaussian stationary stochastic process defined on a finite interval $[1,N]$. More explicitly we have $\vect{y}(k) := [\vec{y}_1(k)\,\dots\,\vec{y}_m(k)]^\top\in\RR^m$, $k=1,\dots,N$, therefore the process is completely characterized by the random vector $\vect{y} := [\vec{y}_1(1)\,\dots\,\vec{y}_m(1)\,\dots\,\dots\,\vec{y}_1(N)\,\dots\,\vec{y}_m(N)]^\top\in\RR^{mN}$. In \cite{CarFerrPav} it has been shown that $\vect{y}$ is a restriction of a wide-sense stationary periodic process of period $N$ defined on the whole integer line $\ZZ$ if and only if the $mN\times mN$ covariance matrix $\vect{\Sigma}$ of $\vect{y}$ is symmetric block-circulant:
\begin{equation}\label{eq:circmtx}
   \vect{\Sigma} =  \text{circ}\{\Sigma_0,\Sigma_1,\dots,\Sigma_\frac{N}{2},\dots,\Sigma_1^\top\},
\end{equation}
where $\EE[\vect{y}(i)\vect{y}(j)^\top] = \Sigma_{i-j}$, $i,j = 1,\dots,N$, are the covariance lags of the process such that $\Sigma_k = \Sigma_{N-k}^\top$ for $k>N/2$. In view of the above equivalence, we will denote with $\vect{y}$ both the wide-sense stationary periodic process defined in the whole line $\ZZ$ and its restriction, depending on the context.
A particular class of stationary periodic processes is represented by reciprocal processes.

\begin{definition}
	We say that $\vect{y}$ is a (periodic) reciprocal process of order $n$ on $[1,N]$ if the random variables of the process in the interval $(t_1,t_2)\subset[1,N]$ are conditionally orthogonal to the random variables in $(t_1,t_2)^c$, given the $2n$ boundary values $\vect{y}(t_1-n+1),\dots,\vect{y}(t_1),\vect{y}(t_2),\dots,\vect{y}(t_2+n-1)$, where the sums $t-k$ and $t+k$ are to be understood modulo $N$.
\end{definition}

The following result has been proved in \cite[Theorem 3.3]{CarFerrPav}, and it simply states that a reciprocal model is completely specified by a block-circulant matrix of type \eqref{eq:circmtx} whose inverse is banded, block-circulant as in \eqref{eq:bandmtx}.

\begin{theorem}\label{thm:conc}
	A non-singular $mN \times mN$-dimensional matrix $\vect{\Sigma}$ is the covariance matrix of a periodic reciprocal process of order $n$ if and only if its inverse is a positive definite symmetric block-circulant matrix which is banded of bandwidth $n$, namely $\vect{\Sigma}^{-1}\in\B$.
\end{theorem}
\mbox{\\}

We are now ready to deal with the identification problem of a reciprocal process. Let $\hat{\Sigma}_0,\dots,\hat{\Sigma}_n$ be given estimates of the first $n+1$ covariance lags $\Sigma_0,\dots,\Sigma_n$ of the underlying reciprocal process. According to Theorem \ref{thm:conc}, the identification of a reciprocal model can be formulated as a matrix completion problem.

\begin{problem}\label{pb:covext}
	Given the $n+1$ estimates $\hat{\Sigma}_0,\dots,\hat{\Sigma}_n$, compute a sequence $\Sigma_{n+1},\dots,\Sigma_\frac{N}{2}$, in such a way to form a symmetric, positive definite block-circulant matrix
	\[
	   \vect{\Sigma} = \text{circ}\{\hat{\Sigma}_0,\dots,\hat{\Sigma}_n,\Sigma_{n+1},\dots,\Sigma_\frac{N}{2},\dots,\Sigma_{n+1}^\top,\hat{\Sigma}_n^\top,\dots,\hat{\Sigma}_1^\top\},
	\]
	with $\Sigma^{-1}\in\B$.
\end{problem}
It has been shown in \cite{CarFerrPav} that the condition $\Sigma^{-1}\in\B$ is equivalent to maximizing the entropy of the process so that the previous problem is equivalent to the following optimization program:
\begin{equation}\label{op:mep}
\begin{aligned}
\argmax_{\vect{\Sigma}\in\C} &\quad \log\det\vect{\Sigma}\\
\text{subject to } &\quad \vect{\Sigma}>0\\
                   &\quad \mathsf{P}_{\mathscr{B}}(\vect{\Sigma}-\hat{\vect{\Sigma}}) = 0.
\end{aligned}
\end{equation}
whose dual problem has been proven to be
\begin{equation}\label{op:mepdu}
\begin{aligned}
\argmin_{\vect{X}\in\B} &\quad -\log\det\vect{X} + \left<\vect{X},\,\hat{\vect{\Sigma}}\right>_\C\\
\text{subject to } &\quad \vect{X}>0
\end{aligned}
\end{equation}
where $\hat{\vect{\Sigma}}\in\B$ is the symmetric, banded block-circulant matrix of bandwidth $n$,
$
   \hat \Sigma= \text{circ}\{\hat \Sigma_0,\hat \Sigma_1,\dots,\hat \Sigma_n,0,\dots,0,\hat \Sigma_n^\top,\dots,\hat \Sigma_1^\top\},
$
containing the covariance lags estimated from the data and
the optimal value of dual variable $\vect{X}$ is indeed equal to $\Sigma^{-1}$, i.e. the inverse of the solution of problem \ref{pb:covext}.
It can be shown that strong duality holds between \eqref{op:mep} and \eqref{op:mepdu}, so that \eqref{op:mep} and \eqref{op:mepdu} are equivalent. In what follows we assume that $\hat{\vect{\Sigma}}>0$ as it is a necessary condition for problem \eqref{op:mep} to be feasible. In the case that $\hat{\vect{\Sigma}}$ is not positive definite, we can consider a positive definite banded block-circulant matrix sufficiently close to $\hat{\vect{\Sigma}}$ which can be obtained by solving a structured covariance estimation problem, see \cite{structcov1}, \cite{structcov2}.

\begin{remark}\label{rmk:approx}
	Recall that, for $N\to\infty$, Toeplitz matrices can be approximated arbitrarily well by circulant matrices, see \cite[Lemma 4.2]{Gray}. Accordingly, for $N\to\infty$, problem \ref{pb:covext} consists in searching a completion that leads to an infinite positive definite block-Toeplitz covariance matrix, i.e. such that the Fourier transform of the resulting extended sequence is a power spectral density. By Theorem 3.1 of \cite{FerrLinAlg}, for $N\to\infty$, problem \eqref{op:mep} is the classical Burg's maximum entropy problem \cite{BurgPhd,family_zorzi,ZORZI201587,PICCI_LINDQUIST}, whose solution is an autoregressive process of order $n$. In light of this observation, we can understand the reciprocal process solution of \eqref{op:mepdu} as an approximation of the AR process solution of the Burg's maximum entropy problem. In the following sections we will exploit this approximation for the identification of sparse AR graphical models.
\end{remark}

The reciprocal approximation described in Remark \ref{rmk:approx} has also an interesting interpretation in the frequency domain. Let $\Phi$ denote the power spectrum of the autoregressive, wide-sense stationary process $\vect{x}$,
\begin{equation*}
\Phi(e^{i\theta}) = \sum_{k=-\infty}^{\infty}\,R_k\,e^{-i\theta k},\qquad R_{-k}=R_k^\top,\ \ \theta\in[-\pi,\pi].
\end{equation*}
The reciprocal approximation translates in sampling the power spectrum over the interval $[-\pi,\pi]$, with sample frequency $2\pi /N$, obtaining the \emph{symbol} of the covariance matrix of the corresponding reciprocal process:
\[
\Phi(\zeta) = \sum_{k=0}^{N-1}\, \Sigma_k \,\zeta^{-k},\qquad \Sigma_{k}=\Sigma_{N-k}^\top\text{ for } k>\frac{N}{2}.
\]
Figure \ref{fig:recapx} gives an intuitive idea of this approximation. According to Proposition \ref{prop:diagsymb}, the covariance matrix $\vect{\Sigma}$ of the reciprocal process $\vect{y}$ that approximates $\vect{x}$ writes as
\begin{equation}\label{eq:covrec}
\vect{\Sigma} = \vect{F}^*\text{circ}\{ \Phi(\zeta^0),\,\Phi(\zeta^1),\dots,\Phi(\zeta^{N-1}) \}\vect{F},
\end{equation}
hence, its inverse
\begin{equation}\label{eq:invcovrec}
\vect{\Sigma}^{-1} = \vect{F}^*\text{circ}\{ \Phi(\zeta^0)^{-1},\,\Phi(\zeta^1)^{-1},\dots,\Phi(\zeta^{N-1})^{-1} \}\vect{F},
\end{equation}
can be robustly computed by inverting the $N$ blocks $\Phi(\zeta^0),\,\Phi(\zeta^1),\dots,\Phi(\zeta^{N-1})$, all of size  of size $m\times m$. As a final remark, we recall that eigevalues and eigenvectors of circulant matrices can be robustly computed as well, thanks to the availability of closed-form formulas, see for instance \cite{Gray}.
\vspace{2mm}
\begin{figure}[h!]\centering
	\includegraphics[scale=0.8]{recapx.pdf}
	\caption{Power spectrum $\Phi(e^{i\theta})$ and its sampled version $\Phi(\zeta)$ with $N=12$ samples.}\label{fig:recapx}
\end{figure}\\
As highlighted by the frequency-domain interpretation, the goodness of the approximation strictly depends on the regularity of the spectrum: the less the spectrum is regular, the larger $N$ has to be chosen in order to get a good approximation of the AR process.

\section{GRAPHICAL MODELS}\label{sec:gms}
Consider a Gaussian random vector $\vect{x}=[\vec{x}_1\,\dots\,\vec{x}_m]^\top$ with covariance matrix $\vect{\Sigma}$ and let $K := \Sigma^{-1}$ be its \emph{concentration matrix}. The notation
\[
   \vec{x}_i \indep \vec{x}_j \mid \{\vec{x}_k\}_{k\ne i,j}
\]
means that the random variable $\vec{x}_i$ is \emph{conditionally independent} from the random variable $\vec{x}_j$ given the remaining random variables $\vec{x}_k$, $k\ne i,j$. It can be proven that, \cite{Lauritzen}:
\begin{equation}\label{eq:cic}
	\vec{x}_i \indep \vec{x}_j \mid \{\vec{x}_k\}_{k\ne i,j}
	\iff
	k_{ij}=0,
\end{equation}
where $k_{ij}:=(K)_{ij}$ is the element in position $(i,j)$ in the concentration matrix $K$, $i,\,j=1,\dots,m$. Relation \eqref{eq:cic} defines an undirected graph $\G = (V,E)$, $E\subset V\times V$, associated to the random vector $\vect{x}$, whose nodes are the components $\vec{x}_1,\dots,\vec{x}_m$ of $\vect{x}$, and the absence of edges describes conditional independence  between the components, namely for $i\ne j$,
\[
    (i,j)\notin E \iff \vec{x}_i \indep \vec{x}_j \mid \{\vec{x}_k\}_{k\ne i,j}.
\]
The graph $\G$ is called the \emph{graphical model} associated to $\vect{x}$. Property \eqref{eq:cic} provides a complete characterization of the graphical model associated to a certain Gaussian random vector in terms of its concentration matrix.
In practice, there is a large interest in \emph{sparse graphical models}, i.e. graphs that describe the interactions between a large number of components $x_i$s with few edges (equivalently with $K$ being a sparse matrix), and thus give an easily understandable description of the underlying system we are modeling.

Although there is a large literature that deals with sparse graphical models \cite{SongVan}, \cite{ZorzSep,ZorzCh,songsiri2010graphical,LindqARMA}, \cite{ChaParWill,candes2012exact,candes2010matrix}, the problem of deriving such models for the case in which the underlying process is a reciprocal process seems to not have been considered till now. In what follows we present sparse graphical models associated to reciprocal processes introduced in Section \ref{sec:recpr}. We will describe how the combination of the underlying reciprocal structure and the sparsity constraint on the concentration matrix of the process impact the properties of the resulting sparse graphical model. The proofs of the following results can be found in \cite{tesi}.\\
Let $\vect{y}$ be a Gaussian, periodic, reciprocal process of order $n$ defined on $[1,N]$ with covariance matrix $\vect{\Sigma}\in\C$ and let $\vect{S}:=\vect{\Sigma}^{-1}$ denote its \emph{concentration matrix} so that, according to Theorem \ref{thm:conc},
\begin{equation}\label{eq:Smtx}
	\vect{S} =\text{circ}\{S_0,S_1,\dots,S_n,0,\dots,0,S_n^\top,\dots,S_1^\top\}.

\end{equation}
In the following we will generalize the characterization of conditional independence we have given in the classical setting of Gaussian random vectors to the case of graphical models associated to Gaussian reciprocal processes. For this purpose, it is useful to define the $j$-th component of the process $\vect{y}$ as the $\RR^N$-valued vector $\vect{y}_j:=[\vec{y}_j(1) \dots \vec{y}_j(N)]^\top$, obtained by stacking all the $j$-th components of the process for each $k=1,\dots,N$. The components of the reciprocal process are defined \emph{for any} $k\in\ZZ$. The process, however, is periodic of period $N$ so that we can impose conditional independence only for $k\in[1,N]$. By Property \eqref{eq:cic}, this implies that $\vect{S}$ is a \emph{sparse matrix} and that the blocks $S_0, S_1,\dots,S_n$ have common support $\Omega\subseteq\{(i,j):\,i,j=1,\,\dots,\,m\}$ namely,
\begin{equation}\label{eq:cs}
(S_k)_{ij} = (S_k)_{ji} = 0,\qquad k=0,\,\dots\,,n,\quad\forall\,(i,j)\in\Omega^c,
\end{equation}
where $\Omega$ is the set of conditionally dependent pairs that necessarily contains all the pairs $(i,\,i),\,i=1,\dots,m$, since conditional independence is not defined between one variable and itself. The above relation is equivalent to
\begin{equation}\label{eq:cs1}
\begin{aligned}
\EE\big[\vec{y}_i(t_1)\,\vec{y}_j(t_2)\mid\,&\vec{y}_h(s),\,h\ne i,j,\,\,s=1,\dots,N,\\
& \vec{y}_i(s_1),\,s_1\ne t_1,\,\vec{y}_j(s_2),\,s_2\ne t_2\big]=0,
\end{aligned}
\end{equation}
for any $t_1,\,t_2\in[1,N]$ and for any pair $(i,j)\in\Omega^c$. The following result has been proved in \cite{tesi}.

\begin{proposition}\label{thm:tesi}
	Condition \eqref{eq:cs1} is equivalent to
	\begin{equation}\label{eq:condindp}
	\EE\left[ \vec{y}_i(t_1)\,\vec{y}_j(t_2)\mid
	\vec{y}_h(s),\,h\ne i,j,\,\,s=1,\dots,N\right] = 0,
	\end{equation}
	for any $t_1,\,t_2\in[1,N]$ and for any $(i,j)\in\Omega^c$.
\end{proposition}

The above result reflects the fact that the random variables $\{\vec{y}_i(s_1),\vec{y}_j(s_2),\,s_1\ne t_1,\,s_2\ne t_2\}$ do not play any role in the conditioning \eqref{eq:cs1}.

According to \eqref{eq:condindp} we can associate a graphical model to the process $\vect{y}$, whose nodes are its $m$, $N$-dimensional components $\vect{y}_1,\dots,\vect{y}_m$ and the presence of an edge between two vectors $\vect{y}_i$ and $\vect{y}_j$, $i\ne j$, means that $\vect{y}_i$ and $\vect{y}_j$ are conditionally dependent. According to the characterization of conditional dependence, such an edge is described by the quantities
\[
\EE\left[ \vec{y}_i(t_1)\,\vec{y}_j(t_2)\mid
\vec{y}_h(s),\,h\ne i,j,\,\,s=1,\dots,N\right]
\]
for any $t_1,\,t_2\in[1,N]$ and for $i,j=1,\dots,m$, i.e. the edge $(i,j)$ is completely characterized by the vector
\[
\left[ (S_0)_{ij}\,\,(S_1)_{ij}\,\,\dots\,\,(S_n)_{ij}\,\, 0 \,\, \dots \,\, 0 \,\, (S_n)_{ji}\,\,\dots\,\, (S_1)_{ji}\right].
\]
By Property \eqref{eq:cic}, the identification of such graphical models consists in the estimation of the sequence $S_0,S_1,\dots,S_n$ obeying the structural constraint (\ref{eq:cs}).

\begin{example}
	Consider the case in which $N=6$, $n=2$ and $m=4$, and suppose that the graphical model associated to process $\vect{y}$ is the one depicted in
	Figure \ref{fig:graphprop}.
	\begin{figure}[h!]\centering
		\includegraphics[scale=1.2]{graphprop.pdf}
		\caption{Example of a sparse reciprocal graphical model.}\label{fig:graphprop}
	\end{figure}\\
    In this case, the concentration matrix of vector $\vect{y}$ is
    \[
       \vect{S} =\text{circ}\{S_0,\,S_1,\,S_2,\,0,\,S_2^\top,\,S_1^\top\},
    \]
    where $(S_j)_{12},\,(S_j)_{24},\,(S_j)_{34}$, $j=0,1,2$, are the unique entries different from zero. In this case, $\Omega=\{ (i,i):\,i=1,\dots,4\}\cup\{(1,2),\,(2,4),\,(3,4)\}$.
\end{example}

\section{IDENTIFICATION OF SPARSE RECIPROCAL GRAPHICAL MODELS}\label{sec:id}
Let $\vect{x}:=\{\vect{x}(t):\,t\in\ZZ\}$ be an $m$-dimensional, AR Gaussian stationary process of order $n$,
\begin{equation}\label{eq:AR}
\sum_{k=0}^{n}\,B_k\,\vect{x}(t-k) = \vect{e}(t),\qquad\vect{e}(t)\sim\N(0,I_m),
\end{equation}
defined in the whole integer line $\ZZ$, and let $R_k := \EE[\vect{x}(t)\vect{x}(t-k)^\top]$, $k\in\ZZ$, be its $k$-th covariance lag. Suppose now that $T$ observations $x(1),\dots,x(T)$ are available, and let
\begin{equation}\label{eq:estcov}
\hat{R}_k = \frac{1}{T}\sum_{t=k}^T\,x(t)x(t-k)^\top,\qquad k = 0,1,\dots,n,
\end{equation}
be estimates of the first $n+1$ covariance lags $R_0,\dots,R_n$.
In view of Remark \ref{rmk:approx}, the idea is to approximate process $\vect{x}$ with a Gaussian reciprocal process $\vect{y}$ of order $n$ defined over the interval $[1,N]$, with $N$ sufficiently large, having a sparse graphical model.

By what we have explained in Section \ref{sec:recpr}, the introduction of the reciprocal approximation allows to obtain a robust procedure even in the case that $n$ is large. In fact, since the matrices that are involved in the optimization are symmetric and block-circulant, according to relations \eqref{eq:covrec} and \eqref{eq:invcovrec}, we can compute the corresponding inverse matrices and eigenvalues in a  robust way. Moreover, it is apparent from \eqref{eq:invcovrec} that the identification algorithm we are proposing scales with respect to $n$ because the dimensions of the matrices, whose eigenvalues must be computed in the optimization procedure, depend only on $m$. This ensures robustness in the results even if the order of the AR process we are considering is large.\\
Now we can formally state the identification problem.

\begin{problem}\label{prob2}
Consider an $m$-dimensional process $\vect{x}$ and let $\hat{R}_0,\dots,\hat{R}_n$ be the estimates of the first $n+1$ covariance lags of $\vect{x}$ given by \eqref{eq:estcov}. Set $\Sigma_0:=\hat{R}_0,\dots,\Sigma_n:=\hat{R}_n$. Compute the blocks $\Sigma_{n+1},\dots,\Sigma_\frac{N}{2}$ of the  block-circulant covariance matrix
$
\Sigma = \text{circ}\{\Sigma_0,\Sigma_1,\dots,\Sigma_{\frac{N}{2}-1},\Sigma_{\frac{N}{2}},\Sigma_{\frac{N}{2}-1}^\top,\dots,\Sigma_1^\top\}
$
such that its inverse $\vect{S}\in\B$, i.e.
\[
   \vect{S} =\text{circ}\{S_0,S_1,\dots,S_n,0,\dots,0,S_n^\top,\dots,S_1^\top\},
\]
and the blocks $S_0,\dots,S_n$ have common support $\Omega$ as small as possible.
\end{problem}

Clearly, the matrix $\Sigma$ solving problem \ref{prob2} is the covariance of the
reciprocal process $\vect{y}$ approximating $\vect{x}$ and featuring a sparse graphical model.

Since we are going to identify a model for a reciprocal process, we can exploit the maximum entropy dual problem \eqref{op:mepdu} recalled before. It is worth noting that the support $\Omega$ is not known in advance, thus it has to be estimated from the data. In order to do that, inspired by \cite{SongVan}, we consider the following regularizer
\[
   h_\infty(\vect{S}) = \sum_{k>h}\,\max\left\{ |(S_0)_{hk}|,2\max_{j=1,\dots,n}|(S_j)_{hk}|,2\max_{j=1,\dots,n}|(S_j)_{kh}|\right\},
\]
which is basically a generalization of the $\ell^\infty$-norm used to induce sparsity on vectors. The optimization problem for the estimation of a sparse reciprocal model for the process is a regularized version of problem \eqref{op:mepdu}:
\begin{equation}\label{op:regop}
\begin{aligned}
\argmin_{\vect{S}\in\B} &\quad-\log\det\vect{S} + \left<\hat{\vect{\Sigma}},\vect{S}\right>_\C + \lambda_S\,h_\infty(\vect{S})\\
\text{subject to } &\quad \vect{S}>0,
\end{aligned}
\end{equation}
where $\lambda_S>0$ is the regularization parameter. Further research is needed to understand weather \eqref{op:regop} can be seen as the dual of some kind of entropy-related optimization problem.

Notice that, although the objective function in \eqref{op:regop} is strictly convex in $\vect{S}$, it is non-differentiable due to the presence of the regularizer $h_\infty$. For this reason, we consider the dual of problem \eqref{op:regop}: as we will see below, the dual objective function is smooth and therefore it is suitable to be minimized by a projected gradient approach making the implementation of the optimization algorithm easy. Introducing the auxiliary variable $\vect{Y}\in\B$, problem \eqref{op:regop} can be rewritten as
\begin{equation}\label{op:regop1}
\begin{aligned}
\argmin_{\vect{S}\in\C,\vect{Y}\in\B} &\quad-\log\det\vect{S} + \left<\hat{\vect{\Sigma}},\vect{S}\right>_\C + \lambda_S\,h_\infty(\vect{Y})\\
\text{subject to } &\quad \vect{S}>0,\quad \vect{Y} = \vect{S}.
\end{aligned}
\end{equation}
Exploiting strong-duality between \eqref{op:regop1} and its dual, we address problem \eqref{op:regop1} using Lagrange multipliers theory. The Lagrangian for this problem is
\begin{align*}
\mathcal{L}(\vect{S},\vect{Y},\vect{Z}) = &-\log\det\vect{S} + \left<\hat{\vect{\Sigma}},\,\vect{S}\right>_\C + \lambda_S\,h_\infty(\vect{Y})+\left< \vect{Z},\vect{S}-\vect{Y}\right>_\C\\
= &-\log\det\vect{S} + \left<\hat{\vect{\Sigma}}+\vect{Z},\,\vect{S}\right>_\C + \lambda_S\,h_\infty(\vect{Y}) - \left<\vect{Z},\,\vect{Y}\right>_\C
\end{align*}
where $\vect{Z}\in\C$ is the Lagrange multiplier. The dual objective function is the infimum over $\vect{S}$ and $\vect{Y}$ of the Lagrangian. The unique term on $\mathcal{L}$ that depends on $\vect{Y}$ is
$ \lambda_S\,h_\infty(\vect{Y}) - \left<\vect{Z},\,\vect{Y}\right>_\C$. The latter
is bounded below if and only if
\begin{align}
	&\text{diag}(Z_j)=0,\quad  j=0,\,\dots,\,n, \label{eq:minY1}\\
	&2|(Z_0)_{kh}|+\sum_{j=1}^n\,|(Z_j)_{kh}|+|(Z_j)_{hk}|\le\frac{\lambda_S}{N},\quad k>h,\label{eq:minY2}
\end{align}
in which case the infimum is zero. Accordingly,
\[
\inf_{\vect{Y}\in\B}\,\mathcal{L}=
\left\{
\begin{split}
&-\log\det\vect{S} + \left<\hat{\vect{\Sigma}}+\vect{Z},\,\vect{S}\right>_\C \quad\text{ if }\quad\eqref{eq:minY1},\,\eqref{eq:minY2}\text{ hold,}\\\\
&-\infty\quad\text{ otherwise.}
\end{split}
\right.
\]
If \eqref{eq:minY1} and \eqref{eq:minY2} hold, it remains to minimize the strictly convex function (of $\vect{S}$)
$
  \bar{\mathcal{L}}(\vect{S}) := -\log\det\vect{S} + \left<\hat{\vect{\Sigma}}+\vect{Z},\,\vect{S}\right>_\C
$
over the symmetric, positive definite, banded block-circulant matrices. Observe that, $\forall\,\vect{Z}\in\C$, and for any sequence $\vect{S}_k>0$ converging to a singular matrix,
\[
\lim_{k\rightarrow \infty}\bar{\mathcal{L}}(\vect{S}_k)=\infty.
\]
Accordingly, we can assume that the solution lies in the interior of the cone so that a necessary and sufficient condition for $\vect{S}_o$ to be a minimum point for $\bar{\mathcal{L}}$ is that its first Gateaux derivative computed at $\vect{S}=\vect{S}_o$ is equal to zero in every direction $\vect{\delta}\vect{S}$, namely
\begin{equation}\label{eq:ddcond}
   \vect{\delta}\bar{\mathcal{L}}(\vect{S}_o;\vect{\delta}\vect{S}) = \tr\left[ \left(-\vect{S}_o^{-1} + \hat{\vect{\Sigma}} + \vect{Z}\right)\vect{\delta}\vect{S}\right] = 0,
   \qquad \forall\,\vect{\delta}\vect{S}\in\C.
\end{equation}
Assuming that $\vect{Z}\in\C$ is such that
\begin{equation}\label{eq:pdcstr}
	\hat{\vect{\Sigma}}+\vect{Z} > 0,
\end{equation}
condition \eqref{eq:ddcond} is satisfied if and only if $\vect{S}_o = (\hat{\vect{\Sigma}}+\vect{Z})^{-1}$. Finally, we have that
\[
\inf_{\vect{Y}\in\B,\vect{S}\in\C}\,\mathcal{L}=
\left\{
\begin{split}
&\log\det(\hat{\vect{\Sigma}}+\vect{Z}) + mN,  && \text{ if } \eqref{eq:minY1},\,\eqref{eq:minY2},\,\eqref{eq:pdcstr}\text{ hold,}\\\\
&-\infty&&\text{ otherwise.}
\end{split}
\right.
\]
The dual problem of problem \eqref{op:regop} follows straightforward
\begin{equation}\label{op:dual}
\begin{aligned}
\argmin_{\vect{Z}\in\C} &\quad -\log\det(\hat{\vect{\Sigma}}+\vect{Z}) - mN\\
\text{subject to } &\quad \eqref{eq:minY1},\,\eqref{eq:minY2}\\
&\quad \hat{\vect{\Sigma}}+\vect{Z}>0.
\end{aligned}
\end{equation}

\begin{proposition}
	Under the assumption that $\hat{\vect{\Sigma}}\in\B$ and $\hat{\vect{\Sigma}}>0$, problem \eqref{op:dual} admits a unique solution.\\
	\begin{IEEEproof}

		Define $f(\vect{Z}):=\log\det(\hat{\vect{\Sigma}}+\vect{Z})$. Let
		\[
			\Q :=\left\{ \vect{Z}\in\C\,|\,\hat{\vect{\Sigma}}+\vect{Z}>0,\text{ and }\eqref{eq:minY1},\,\eqref{eq:minY2}\text{ hold}\right\}
		\]
		be the set of constraints of problem \eqref{op:dual}. First of all, notice that constraints \eqref{eq:minY1} and \eqref{eq:minY2} ensure that $\Q$ is a bounded subset of $\C$. Indeed, the entries of any $\vect{Z}\in\Q$ are bounded by $\lambda_S/N$ in the element-wise max-norm of the matrix. By the equivalence of norms in finite-dimensional spaces, this implies in particular that $\|\vect{Z}\|_\C<\infty$ for any $\vect{Z}\in\Q$. Let now $(\vect{Z}^{(k)})_{k\in\NN}$ be a generic sequence of elements of $\Q$ converging to some $\bar{\vect{Z}}\in\C$, such that $\hat{\vect{\Sigma}}+\bar{\vect{Z}}\ge0$ singular. Then
		\[
			\lim_{k\to\infty}\,-\log\det(\hat{\vect{\Sigma}}+\vect{Z}^{(k)})=+\infty,
		\]
		and therefore $\vect{Z}^{(k)}$ is not an infimizing sequence. Hence, we can restrict our attention to the closed subset of $\Q$ defined by
		\[
		\bar{\Q}:=\left\{ \vect{Z}\in\C\,|\,\hat{\vect{\Sigma}}+\vect{Z}\ge\epsilon I_{mN},\text{ and }\eqref{eq:minY1},\,\eqref{eq:minY2}\text{ hold}\right\}
		\]
		with $\epsilon>0$ small enough. By what we have shown till now, the function $f$ is continuous on the compact set $\bar{\Q}$ and therefore admits at least one minimum point. Since $f$ is strictly convex, the minimum is unique.
	\end{IEEEproof}
\end{proposition}

\begin{proposition}
	Under the assumption that $\hat{\vect{\Sigma}}\in\B$ and $\hat{\vect{\Sigma}}>0$, problem \eqref{op:regop} admits a unique solution $\vect{S}_o$.
\end{proposition}
\begin{IEEEproof}
	Notice that problem \eqref{op:regop} is a strictly feasible (for instance, pick $X=I_{mN}$) convex optimization problem. Accordingly, Slater's condition holds, hence strong duality holds between \eqref{op:regop} and its dual. The strong duality between problems \eqref{op:regop} and \eqref{op:dual} and the existence of a unique optimum $\vect{Z}_o$ for the dual problem \eqref{op:dual}, imply that there exists a unique $\vect{S}_o\in\C$ so that
	$
	\vect{S}_o = \left(\hat{\vect{\Sigma}} + \vect{Z}_o\right)^{-1}
	$
	which solves the primal problem \eqref{op:regop}.
\end{IEEEproof}

\section{NUMERICAL EXAMPLE}\label{sec:ne}
Here we present a numerical example that shows how the algorithm works in practice. We consider the AR model $\vect{x}(t)=A(z)\vect{x}(t)+\vect{w}(t)$ with $A(z)=\sum_{k=1}^n A_k z^{-k}$, $m=15$, $n=8$, $\vect{w}(t)$ is white Gaussian noise with variance equal to $15.9$ and the poles $p_j$ of the shaping filter $[I-A(z)]^{-1}$ are depicted in Figure \ref{fig:omega} (right) and they are such that $|p_j|\le 0.9$. Figure \ref{fig:omega} also shows the sparsity pattern of the true inverse spectrum with the sparsity pattern estimated by the proposed algorithm when the regularization parameter is $\lambda_S=125$ and $T=1000$ samples are used to estimate the covariance lags $\hat{R}_k$.
\begin{figure}[h!]
		\begin{minipage}[.c]{0.15\textwidth}
			\centering
			\includegraphics[scale=0.24]{OmegaV.pdf}
		\end{minipage}
		\begin{minipage}[.c]{0.15\textwidth}
			\centering
			\includegraphics[scale=0.24]{OmegaEst.pdf}
		\end{minipage}
		\begin{minipage}[.c]{0.15\textwidth}
			\centering
			\includegraphics[width=3.0cm,height=2.9cm]{pmap.pdf}
		\end{minipage}
		\caption{Sparsity pattern of the inverse spectrum: true (left) and estimated (center). The yellow squares represent the conditional dependent pairs while the blue ones stand for the conditional independent pairs. Poles of the model used to generate the data (right).}
		\label{fig:omega}
\end{figure}\\
From Figure \ref{fig:omega} we can see that only one non-zero element has been estimated in a wrong position. The average relative error between the estimated and the true spectra is about $4\%$.

\section{CONCLUSIONS AND FUTURE DEVELOPMENTS}\label{sec:concl}
In this paper we presented an identification procedure  for a sparse graphical model associated with a reciprocal process.
As discussed in detail in \cite{PICCI_LINDQUIST} and recalled in the introduction, the reciprocal approximation must be understood as an attempt in the direction of the development of an effective procedure for the identification of ARMA graphical models that scales robustly with the product of the process dimension by the length of an accurate AR approximation of the original process. This is a promising theoretical idea that will be tested in simulation and in real examples. Moreover, this approach can be pushed forward in many directions: for example there is the possibility of adding a (small) number of latent variables
to the picture in order  to provide a better approximation of the dynamics of the original process.