$\underline{F}_{\mathcal{A}_{j^{\prime}}}(w_{i})=0$ for all $j<j^{\prime}\leq
k$, we have $\hat{x}_{i^{\prime}j^{\prime}}=\hat{x}_{i^{\prime}j}$ for all
$i^{\prime}\leq i$ and $j^{\prime}\in[j+1,k]$.

We will prove the above theorem by showing a process that gradually converts
any optimal solution $\\{x^{\ast}_{ij}\\}$ to such a solution
$\\{\hat{x}_{ij}\\}$ so that the feasibility is preserved and the objective is
weakly improved. Consider the following process, which we refer to as
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\},\hat{i},\hat{j})$, that takes the
prior $\mathcal{A}$, the allocation variables $\\{x_{ij}\\}$, a valuation type
$\hat{i}\in[0,m]$, and a deadline $\hat{j}\in[1,k-1]$ as inputs, and modifies
$\\{x_{ij}\\}$ into another solution $\\{x^{\prime}_{ij}\\}$ as follows:

  1. (1)

Let $x^{\prime}_{ij}=x_{ij}$ for all $i,j$, except for $i\geq\hat{i}$ and
$j=\hat{j}+1$.

  2. (2)

Initialize $x^{\prime}_{\hat{i}(\hat{j}+1)}=x_{\hat{i}\hat{j}}$.

  3. (3)

If $\hat{i}\neq m$, let $y_{\hat{i}\hat{j}}$ be the smallest $y\geq 0$ such
that

(23)
$\displaystyle\textsf{Area}_{x}(w_{\hat{i}+1},\hat{j})+(w_{i}-w_{\hat{i}+1})\cdot
y=\textsf{Area}_{x}(w_{i^{\prime}},\hat{j}+1)$

is satisfied for some $i^{\prime}>\hat{i}$; equivalently, $y_{\hat{i}\hat{j}}$
is the smallest $y$ such that the straight line $\ell_{y}(\cdot)$ passing
$\big{(}w_{\hat{i}+1},\textsf{Area}_{x}(w_{\hat{i}+1},\hat{j})\big{)}$ with
slope $y$ is tangent to $\textsf{Area}_{x}(v,\hat{j}+1)$ at some
$v=v^{\prime}=w_{i^{\prime}}$ (see Figure 3).

  4. (4)

Let $y^{\prime}_{\hat{i}\hat{j}}=\min\\{y_{\hat{i}\hat{j}},1\\}$, and let
$x^{\prime}_{i(\hat{j}+1)}=\max\\{y^{\prime}_{\hat{i}\hat{j}},x_{i(\hat{j}+1)}\\}\geq
x_{i(\hat{j}+1)}$ for all $i\in[\hat{i}+1,m]$.

  5. (5)

Return $\\{x^{\prime}_{ij}\\}$ as the modified solution.

$v$$w_{2}$$w_{3}$$w_{4}$$w_{1}=w_{\hat{i}}$$\textsf{Area}_{x}(v,\hat{j}+1)$$\textsf{Area}_{x}(v,\hat{j})=\textsf{Area}_{x^{\prime}}(v,\hat{j})$$\ell_{y_{\hat{i}\hat{j}}}$$\textsf{Area}_{x^{\prime}}(v,\hat{j}+1)$slope
$=y_{\hat{i}\hat{j}}$ Figure 3. Illustration of
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\},\hat{i},\hat{j})$ with $\hat{i}=1$.
In step (3), $y_{\hat{i}\hat{j}}$ is the smallest $y$ such that the straight
line $\ell_{y}(\cdot)$ passing
$\big{(}w_{2},\textsf{Area}_{x}(w_{2},\hat{j})\big{)}$ with slope $y$ is
tangent to $\textsf{Area}_{x}(v,\hat{j}+1)$. The tangent point is
$v^{\prime}=w_{4}$, and thus
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\},\hat{i},\hat{j})$ sets
$x^{\prime}_{2(\hat{j}+1)}=x^{\prime}_{3(\hat{j}+1)}=y_{\hat{i}\hat{j}}$, and
$x^{\prime}_{i(\hat{j}+1)}=x_{i(\hat{j}+1)}$ for $i\geq 4$. The resulting new
area curve $\textsf{Area}_{x^{\prime}}(v,\hat{j}+1)$ shown in red remains
convex, and lies between $\textsf{Area}_{x}(v,\hat{j})$ and
$\textsf{Area}_{x}(v,\hat{j}+1)$ at all $v$.

Recall from Definition 4.2 that
$\hat{i}_{j}\coloneqq\max\\{i:\underline{F}_{\mathcal{A}_{j^{\prime}}}(v_{i})=0\quad\forall
j^{\prime}\in[j,k]\\}$ is the highest type $i$ such that no buyer with
valuation at most $v_{i}$ and deadline at least $j$ exists in $\mathcal{A}$.

Define $\mathsf{Align}(\mathcal{A},\\{x_{ij}\\})$ as follows: Apply
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\},\hat{i},\hat{j})$ in increasing order
of $\hat{j}$ from $\hat{j}=1$ to $k-1$, and for each $\hat{j}$, in increasing
order of $\hat{i}$ from $\hat{i}=0$ to $\hat{i}_{\hat{j}}$. We next show that
(i) If $\\{x_{ij}\\}$ input to
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\},\hat{i},\hat{j})$ is feasible, so is
the constructed $\\{x^{\prime}_{ij}\\}$ (via Lemmas B.3\- B.4); and (ii) If
the input $\\{x_{ij}\\}$ is optimal for D4$(\mathcal{A})$, then so is
$\\{\hat{x}_{ij}\\}$ (via Lemma B.5).

Consider the execution of $\mathsf{Align}(\mathcal{A},\\{x_{ij}\\})$, and the
instant just before the execution of
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\},\hat{i},\hat{j})$ for some
$\hat{i}\in[0,m]$ and $\hat{j}\in[1,k-1]$. The next lemma assumes the
following inductive scenario:

Inductive Scenario $P(\hat{i},\hat{j})$: Assume the input $\\{x_{ij}\\}$ to
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\},\hat{i},\hat{j})$ is feasible for
D4$(\mathcal{A})$. Further, by the execution of
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\})$, we have
$x_{i(\hat{j}+1)}=x_{i\hat{j}}$ holds for all $i\in[0,\hat{i}-1]$.

###### Lemma B.3.

Assuming $P(\hat{i},\hat{j})$, the procedure
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\},\hat{i},\hat{j})$ is well-defined,
and guarantees that the output $\\{x^{\prime}_{ij}\\}$ of
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\},\hat{i},\hat{j})$ is feasible.

###### Proof.

Assume $P(\hat{i},\hat{j})$ holds. Then
$(w_{\hat{i}},\textsf{Area}_{x}(w_{\hat{i}},\hat{j}))$ lies weakly below the
curve $\textsf{Area}_{x}(v,\hat{j}+1)$. By the convexity of
$\textsf{Area}_{x}(v,\hat{j}+1)$, $y_{\hat{i}\hat{j}}$ uniquely exists, i.e.,
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\},\hat{i},\hat{j})$ is well-defined.
Furthermore, we have $y_{\hat{i}\hat{j}}\geq x_{(\hat{i}+1)\hat{j}}$, since
the latter is the slope of $\textsf{Area}_{x}(v,\hat{j})$ at
$v=(w_{\hat{i}+1})^{+}$. Therefore we have
$x^{\prime}_{(\hat{i}+1)(\hat{j}+1)}=\min\\{y_{\hat{i}\hat{j}},1\\}\geq
x_{(\hat{i}+1)\hat{j}}\geq
x_{\hat{i}\hat{j}}=x^{\prime}_{\hat{i}(\hat{j}+1)}.$ By $P(\hat{i},\hat{j})$,
we also have $x^{\prime}_{i(\hat{j}+1)}=x_{i\hat{j}}$ for $i\in[0,\hat{i}]$,
which implies $\\{x^{\prime}_{i(\hat{j}+1)}\\}$ is non-decreasing for
$i\in[0,\hat{i}]$. Finally, notice that $\\{x^{\prime}_{i(\hat{j}+1)}\\}$ is
the max of two non-decreasing functions for $i\in[\hat{i}+1,m]$. Combining the
above, $\\{x^{\prime}_{i(\hat{j}+1)}\\}$ is non-decreasing for $i\in[0,m]$.

Suppose $y_{\hat{i}\hat{j}}\leq 1$, and thus
$y^{\prime}_{\hat{i}\hat{j}}=y_{\hat{i}\hat{j}}$. Recall that
$\ell_{y_{\hat{i}\hat{j}}}$ is tangent to $\textsf{Area}_{x}(v,\hat{j}+1)$ at
$v=w_{i^{\prime}}$. By convexity of $\textsf{Area}_{x}(v,\hat{j}+1)$, if
$i^{\prime}<m$, we must have $y_{\hat{i}\hat{j}}\leq
x_{i^{\prime}(\hat{j}+1)}$, so we have

$x^{\prime}_{i(\hat{j}+1)}=\begin{cases}x_{i\hat{j}},&i\in[0,\hat{i}],\\\
y_{\hat{i}\hat{j}},&i\in[\hat{i}+1,i^{\prime}-1],\\\
x_{i(\hat{j}+1)},&i\in[i^{\prime},m-1],\end{cases}\quad\textsf{Area}_{x^{\prime}}(v,\hat{j}+1)=\begin{cases}\textsf{Area}_{x}(v,\hat{j}),&v\in[0,w_{\hat{i}+1}),\\\
\ell_{y_{\hat{i}\hat{j}}}(v),&v\in[w_{\hat{i}+1},w_{i^{\prime}}),\\\
\textsf{Area}_{x}(v,\hat{j}+1),&v\in[w_{i^{\prime}},w_{m}],\end{cases}$

where the interval $[i^{\prime},m-1]$ may be empty when $i^{\prime}=m$.
Otherwise, if $y_{\hat{i}\hat{j}}>1$, then we have

$x^{\prime}_{i(\hat{j}+1)}=\begin{cases}x_{i\hat{j}},&i\in[0,\hat{i}],\\\
1,&i\in[\hat{i}+1,m-1],\end{cases}\quad\textsf{Area}_{x^{\prime}}(v,\hat{j}+1)=\begin{cases}\textsf{Area}_{x}(v,\hat{j}),&v\in[0,w_{\hat{i}+1}),\\\
\ell_{1}(v)<\textsf{Area}_{x}(v,\hat{j}+1),&v\in[w_{\hat{i}+1},w_{m}].\end{cases}$

Recall that $y_{\hat{i}\hat{j}}\geq x_{(\hat{i}+1)\hat{j}}$ and both
$\ell_{y_{\hat{i}\hat{j}}}(v)$ and $\ell_{1}(v)$ goes through
$\big{(}w_{\hat{i}+1)},\textsf{Area}_{x}(w_{\hat{i}+1)},\hat{j})\big{)}$.
Therefore, in every interval in both cases above,
$\textsf{Area}_{x^{\prime}}(v,\hat{j}+1)$ is equal to some function upper
bounded by $\textsf{Area}_{x}(v,\hat{j}+1)$ and lower bounded by
$\textsf{Area}_{x}(v,\hat{j})=\textsf{Area}_{x^{\prime}}(v,\hat{j})$. This
implies

(24)
$\displaystyle\textsf{Area}_{x^{\prime}}(w_{i},\hat{j})\leq\textsf{Area}_{x}(w_{i},\hat{j}+1)\leq\textsf{Area}_{x}(w_{i},\hat{j}+1)\leq\textsf{Area}_{x}(w_{i},\hat{j}+2)=\textsf{Area}_{x^{\prime}}(w_{i},\hat{j}+2)$

holds for all $i\in[0,m]$, i.e., $\\{x^{\prime}_{ij}\\}$ satisfies all inter-
deadline IC constraints (as
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\},\hat{i},\hat{j})$ only modifies the
$(\hat{j}+1)$-th curve). Finally, since $y^{\prime}_{\hat{i}\hat{j}}\leq 1$,
we have $\\{x^{\prime}_{i(\hat{j}+1)}\\}\in[0,1]$, and the lemma follows. ∎

The next lemma now follows by induction on the execution of
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\})$ using the observation from the
proof of Lemma B.3 that any Area curve cannot increase in an iteration.

###### Lemma B.4.

Suppose the input $\\{x_{ij}\\}$ to $\mathsf{Align}(\mathcal{A},\\{x_{ij}\\})$
is feasible to D4$(\mathcal{A})$. Then $\\{\hat{x}_{ij}\\}$ (the output of
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\})$ is also feasible in
D4$(\mathcal{A})$. Furthermore, for all $i\in[0,m],j\in[1,k]$ it holds that

(25)
$\displaystyle\textsf{Area}_{\hat{x}}(w_{i},j)\leq\textsf{Area}_{x}(w_{i},j).$

The next lemma shows the optimality of $\\{\hat{x}_{ij}\\}$ if the input to
$\mathsf{Align}(\mathcal{A},\\{x_{ij}\\})$ is optimal.

###### Lemma B.5.

If the input $\\{x_{ij}\\}$ to $\mathsf{Align}(\mathcal{A},\\{x_{ij}\\})$ is
optimal for D4$(\mathcal{A})$, then the output $\\{\hat{x}_{ij}\\}$ is also
optimal for D4$(\mathcal{A})$.

###### Proof.

Recall that $\hat{p}_{ij}\coloneqq
w_{i}\cdot\hat{x}_{ij}-\textsf{Area}_{\hat{x}}(i,j)$ and $p_{ij}\coloneqq
w_{i}\cdot x_{ij}-\textsf{Area}_{x}(i,j)$ are the prices associated with the
buyer with valuation $w_{i}$ and deadline $j$ when the allocation variables
are $\\{\hat{x}_{ij}\\}$ and $\\{x_{ij}\\}$, respectively. Since
$\\{\hat{x}_{i1}\\}=\\{x_{i1}\\}$, we have $\hat{p}_{i1}=p_{i1}$ for all $i$.

Fix some $j\in[2,k]$. For all $i\in[0,\hat{i}_{j-1}]$ we have
$f_{\mathcal{A}_{j}}(w_{i})=0$ and thus
$f_{\mathcal{A}_{j}}(w_{i})\cdot\hat{p}_{ij}=f_{\mathcal{A}_{j}}(w_{i})\cdot
p_{ij}=0$. For $i\in[\hat{i}_{j-1}+1,m]$, we have $\hat{x}_{ij}\geq x_{ij}$ by
the execution of $\mathsf{Align}(\mathcal{A},\\{x_{ij}\\})$, and thus

$\displaystyle\hat{p}_{ij}$
$\displaystyle=w_{i}\cdot\hat{x}_{ij}-\textsf{Area}_{\hat{x}}(i,j)\geq
w_{i}\cdot x_{ij}-\textsf{Area}_{\hat{x}}(i,j)\geq w_{i}\cdot
x_{ij}-\textsf{Area}_{x}(i,j)=p^{\ast}_{ij},$

where the last inequality follows from Lemma B.4. Therefore, for all
$i\in[1,m]$ we have $f_{\mathcal{A}_{j}}(w_{i})\cdot\hat{p}_{ij}\geq
f_{\mathcal{A}_{j}}(w_{i})\cdot p_{ij}$. Since this holds for all $j\in[1,k]$,
it implies that $\\{\hat{x}_{ij}\\}$ achieves (in D4$(\mathcal{A})$) at least
the objective achieved by $\\{x_{ij}\\}$. By Lemma B.4, it is also feasible,
hence optimal. ∎

##### Proof of Theorem B.2.

Consider the output $\\{\hat{x}_{ij}\\}$ by applying
$\mathsf{Align}(\mathcal{A},\\{x^{\ast}_{ij}\\})$ to any arbitrary optimal
solution $\\{x^{\ast}_{ij}\\}$. For any $i,j$ such that
$\underline{F}_{\mathcal{A}_{j^{\prime}}}(w_{i})=0$ for all
$j^{\prime}\in[j+1,k]$, we have
$i\leq\hat{i}_{j}\leq\hat{i}_{j+1}\leq\cdots\leq\hat{i}_{k-1}$, and thus
$\hat{x}_{ij}=\hat{x}_{i(j+1)}=\cdots=\hat{x}_{ik}$. ∎

####  B.0.3. Matching the Allocation Curves to the Lower Envelope

Next, we show that there exists some optimal solution to D4$(\mathcal{A})$
that satisfies an additional set of properties required for proving Theorem
4.11. We will prove the following theorem that shows the allocation curve only
has breakpoints at values with non-zero mass in the lower envelope.
Furthermore, the allocation for the highest type is deterministic.

###### Theorem B.6.

For any prior $\mathcal{A}=\mathcal{G}(t)$, there is an optimal solution
$\\{\hat{x}_{ij}\\}$ to D4$(\mathcal{A})$ that satisfies the following:

  1. (1)

$\hat{x}_{0j}=0$ for all $j\in[1,k]$;

  2. (2)

If $(w_{a},r)\in\textsf{LE}(\mathcal{A})$, then $\hat{x}_{ij}=\hat{x}_{ir}$
for all $i\in[0,a]$ and $j\geq r$;

  3. (3)

If $(w_{a},r)$ and $(w_{b},r^{\prime})$ are consecutive points on
$\textsf{LE}(\mathcal{A})$ where $a<b$ and $r^{\prime}\geq r$, then
$\hat{x}_{ij}=\hat{x}_{aj}$ for all $i\in[a,b)$ and $j\geq r$;

  4. (4)

$\hat{x}_{mk}=1$.

###### Proof.

$v$$w_{1}$$w_{2}$$w_{3}$$w_{4}$$\textsf{Area}_{x^{\ast}}(v,k)$$\textsf{Area}_{\hat{x}}(v,k)$
Figure 4. Illustration of the modification process in proof of Theorem B.6.
The black and red curves illustrate $\textsf{Area}_{x^{\ast}}(v,k)$ and
$\textsf{Area}_{\hat{x}}(v,k)$ (corresponding to some optimal solution
$\\{x^{\ast}_{ij}\\}$ and the final solution $\\{\hat{x}_{ij}\\}$,
respectively) for the last deadline $k$. $(w_{4},k)$ is the only value-
deadline pair in the support of $\mathcal{A}$, and no value-deadline pair
$(w_{3},j)$ is on the lower envelope of $\mathcal{A}$. This implies (i)
$\hat{x}_{0j}=0$, and thus $\textsf{Area}_{\hat{x}}(w_{1},j)=0$; (ii)
$\textsf{Area}_{\hat{x}}(v,k)$ is linear in $v\in[w_{2},w_{4}]$ without a
breakpoint at $v=w_{3}$; and (iii) $\hat{x}_{ik}=1$ for all $i\geq 4$, i.e.,
$\textsf{Area}_{\hat{x}}(v,k)$ has a slope of $+1$ in $[w_{4},w_{m}]$.

Similar to the proof of Theorem B.2, we show a process that gradually modifies
an optimal solution until it satisfies all the properties above.

Consider any optimal solution $\\{x^{\ast}_{ij}\\}$ to D4$(\mathcal{A})$ in
which $x^{\ast}_{01}>0$. Define another solution $\\{x^{\prime}_{ij}\\}$ such
that $x^{\prime}_{01}=0$ and $x^{\prime}_{ij}=x^{\ast}_{ij}$ for all
$(i,j)\neq(0,1)$. Since the curve $\textsf{Area}_{x^{\prime}}(v,1)$ is upper
bounded by $\textsf{Area}_{x^{\ast}}(v,1)$, it is also upper bounded by
$\textsf{Area}_{x^{\ast}}(v,2)=\textsf{Area}_{x^{\prime}}(v,2)$. Hence all
inter-deadline IC constraints are preserved by $\\{x^{\prime}_{ij}\\}$, so it
is feasible. Since the objective of D4$(\mathcal{A})$ does not depend on
$x_{01}$, $\\{x^{\prime}_{ij}\\}$ achieves (in D4$(\mathcal{A})$) exactly the
objective achieved by $\\{x^{\ast}_{ij}\\}$, and thus is optimal. Take
$\\{\bar{x}_{ij}\\}$ to be the output of
$\mathsf{Align}(\mathcal{A},\\{x^{\prime}_{ij}\\})$. Since
$\\{\bar{x}_{ij}\\}$ satisfies Theorem B.2, it satisfies properties (1) and
(2).

Now consider another solution $\\{\tilde{x}_{ij}\\}$ obtained by the
following:

  * •

Initialize $\tilde{x}_{ij}=\bar{x}_{ij}$ for all $i,j$.

  * •

For all consecutive $(w_{a},r),(w_{b},r^{\prime})\in\textsf{LE}(\mathcal{A})$
such that $b>a+1$ and $r^{\prime}\geq r$, set

$\tilde{x}_{ij}=\frac{\textsf{Area}_{\bar{x}}(w_{a},w_{b},r)}{w_{b}-w_{a}}$

for all $i\in[a,b)$ and $j\in[r,k]$.

Clearly, $\\{\tilde{x}_{ij}\\}$ satisfies property (3) while preserving
properties (1) and (2) from $\\{\bar{x}_{ij}\\}$; for property (1), observe
that $w_{1}\in\text{values}(\mathcal{A})$, and thus
$\tilde{x}_{0j}=\bar{x}_{0j}$ for all $j$. Hence it suffices to prove that
$\\{\tilde{x}_{ij}\\}$ is feasible and optimal to $\textsf{D4}(\mathcal{A})$.

Consider any consecutive
$(w_{a},r),(w_{b},r^{\prime})\in\textsf{LE}(\mathcal{A})$ such that $b>a+1$
and $r^{\prime}\geq r$. For all $j\geq r$ we have

(26)
$\displaystyle\textsf{Area}_{\tilde{x}}(v,j)=\begin{cases}\textsf{Area}_{\bar{x}}(v,j),&v\notin(w_{a},w_{b}),\\\
\frac{w_{b}-w_{i}}{w_{b}-w_{a}}\cdot\textsf{Area}_{\bar{x}}(w_{a},j)+\frac{w_{i}-w_{a}}{w_{b}-w_{a}}\cdot\textsf{Area}_{\bar{x}}(w_{b},j),&v\in(w_{a},w_{b}).\end{cases}$

By convexity of $\textsf{Area}_{\bar{x}}(v,j)$, this implies

(27)
$\displaystyle\textsf{Area}_{\tilde{x}}(v,j)\geq\textsf{Area}_{\bar{x}}(v,j)$

for all $v\in[0,v_{n}]$ and $j\geq r$. By taking $j=r$ in Eq. (27) above, all
inter-deadline IC constraints between deadlines $(r-1)$ and $r$ are preserved
in $\\{\tilde{x}_{ij}\\}$. Furthermore, since
$(w_{a},r),(w_{b},r^{\prime})\in\textsf{LE}(\mathcal{A})$, we have
$\textsf{Area}_{\bar{x}}(v,j)=\textsf{Area}_{\bar{x}}(v,r)$ for all
$v\in[0,b]$, which implies all inter-deadline IC constraints for $i\in[0,b]$
and $j\geq r$ are satisfied by $\\{\tilde{x}_{ij}\\}$ as well. Repeating this
argument for all $a,b$, and $r$ shows all inter-deadline IC constraints in
$\textsf{D4}(\mathcal{A})$ are preserved in $\\{\tilde{x}_{ij}\\}$. Since
$\\{\tilde{x}_{ij}\\}\in[0,1]$ and is monotone non-decreasing for all $j$, it
is feasible.

We next show $\\{\tilde{x}_{ij}\\}$ is optimal. Fix any deadline $r$, and
recall that
$\tilde{p}_{ir}=w_{i}\cdot\tilde{x}_{ir}-\textsf{Area}_{\tilde{x}}(w_{i},r)$
and $\bar{p}_{ir}=w_{i}\cdot\bar{x}_{ir}-\textsf{Area}_{\bar{x}}(w_{i},r)$ are
the prices associated with the buyer $(v,d)=(w_{i},r)$ given the allocation
variables $\\{\tilde{x}_{ij}\\}$ and $\\{\bar{x}_{ij}\\}$, respectively. By
Eq. (26) and the construction of $\\{\tilde{x}_{ij}\\}$, for all $i$ such that
$f_{\mathcal{A}_{r}}(w_{i})>0$, it holds that (i)
$\tilde{x}_{ir}=\bar{x}_{ir}$, (ii)
$\textsf{Area}_{\tilde{x}}(w_{i},r)=\textsf{Area}_{\bar{x}}(w_{i},r)$, and
thus (iii) $\tilde{p}_{ir}=\bar{p}_{ir}$. Thus the revenue raised by
$\\{\tilde{x}_{ij}\\}$ for deadline $r$ equals that of $\\{\bar{x}_{ij}\\}$.
Since this is true for all $r$, $\\{\tilde{x}_{ij}\\}$ raises the same revenue
as $\\{\bar{x}_{ij}\\}$, and thus is optimal.

As the last step, we modify $\\{\tilde{x}_{ij}\\}$ such that it satisfies
property (4) and remains optimal. Let $k^{\prime}\leq k$ be the largest $j$
such that some $(w_{i},j)\in\textsf{LE}(\mathcal{A})$, and let $m^{\prime}$ be
the largest such $i$. Consider the following solution $\\{\hat{x}_{ij}\\}$:

$\displaystyle\hat{x}_{ij}\coloneqq\begin{cases}1,\quad&i\in[m^{\prime},m],j\in[k^{\prime},k];\\\
\tilde{x}_{ij},\quad&\text{otherwise.}\end{cases}$

Therefore, $\hat{x}_{mk}=1$, i.e., property (4) is satisfied. It is easy to
see properties (1) - (3) are preserved by $\\{\hat{x}_{ij}\\}$. Hence it
suffices to show $\\{\hat{x}_{ij}\\}$ is feasible and optimal to
$\textsf{D4}(\mathcal{A})$.

The curve $\textsf{Area}_{\hat{x}}(v,k^{\prime})$ is weakly above
$\textsf{Area}_{\tilde{x}}(v,k^{\prime})$, which by the feasibility of
$\\{\tilde{x}_{ij}\\}$, is weakly above
$\textsf{Area}_{\tilde{x}}(v,k^{\prime}-1)=\textsf{Area}_{\hat{x}}(v,k^{\prime}-1)$.
Hence $\\{\hat{x}_{ij}\\}$ satisfies all inter-deadline IC constraints (note
that for all $j>k^{\prime}$ we have
$\textsf{Area}_{\hat{x}}(v,j)=\textsf{Area}_{\hat{x}}(v,k^{\prime})$ for all
$v\in[0,v_{n}]$). Since $\\{\hat{x}_{ij}\\}\in[0,1]$ and remains monotone non-
decreasing for all $j$, it is feasible.

Again let
$\tilde{p}_{ij}=w_{i}\cdot\tilde{x}_{ij}-\textsf{Area}_{\tilde{x}}(w_{i},j)$
and $\hat{p}_{ij}=w_{i}\cdot\hat{x}_{ij}-\textsf{Area}_{\hat{x}}(w_{i},j)$
denote the price associated with the buyer $(v,d)=(w_{i},j)$ given the
allocation variables $\\{\tilde{x}_{ij}\\}$ and $\\{\hat{x}_{ij}\\}$,
respectively. We have $\hat{p}_{ij}=\tilde{p}_{ij}$ for all $j<k^{\prime}$ or
$i<m^{\prime}$, and the only other $(w_{i},j)$ such that
$f_{\mathcal{A}_{j}}(w_{i})>0$ is $(w_{m^{\prime}},k^{\prime})$. Note that the
two curves $\textsf{Area}_{\hat{x}}(v,k^{\prime})$ and
$\textsf{Area}_{\tilde{x}}(v,k^{\prime})$ are identical for
$v\in[0,w_{m^{\prime}}]$. Thus, we have

(since $\hat{x}_{m^{\prime}k^{\prime}}=1$)
$\displaystyle\hat{p}_{m^{\prime}k^{\prime}}$
$\displaystyle=w_{m^{\prime}}\cdot\hat{x}_{m^{\prime}k^{\prime}}-\textsf{Area}_{\hat{x}}(w_{m^{\prime}},k^{\prime})=w_{m^{\prime}}-\textsf{Area}_{\hat{x}}(w_{m^{\prime}},k^{\prime})$
$\displaystyle=w_{m^{\prime}}-\textsf{Area}_{\tilde{x}}(w_{m^{\prime}},k^{\prime})\geq
w_{m^{\prime}}\cdot\tilde{x}_{m^{\prime}k^{\prime}}-\textsf{Area}_{\tilde{x}}(w_{m^{\prime}},k^{\prime})=\tilde{p}_{m^{\prime}k^{\prime}},$

where the inequality follows from $\tilde{x}_{m^{\prime}k^{\prime}}\leq 1$.
Therefore, we have $\hat{p}_{ij}\geq\tilde{p}_{ij}$ for all $(w_{i},j)$ such
that $f_{\mathcal{A}_{j}}(w_{i})>0$, which implies $\\{\hat{x}_{ij}\\}$ raises
at least the revenue achieved by $\\{\tilde{x}_{ij}\\}$. Hence
$\\{\hat{x}_{ij}\\}$ is optimal, and the theorem follows. ∎

####  B.0.4. Completing proof of Theorem 4.11

Consider an optimal solution $\\{\hat{x}_{ij}\\}$ to D4$(\mathcal{A})$
satisfying Theorem B.6. For all $j\in[1,k]$, let $\hat{x}_{j}(w)=\hat{x}_{ij}$
for all $w\in[w_{i},w_{i+1})$. Let $w_{0}=0$. Then for
$i^{\prime}\in\\{1,\ldots,m\\}$ denote allocation function
$x^{j}_{i^{\prime}}(\cdot)$ such that

$x^{j}_{i^{\prime}}(w_{i})=\begin{cases}1,\quad&w_{i}\geq w_{i^{\prime}};\\\
0,\quad&w_{i}<w_{i^{\prime}}.\end{cases}$

This corresponds to posting a price $w_{i^{\prime}}$ and selling to only the
buyers with valuation $w_{i}\geq w_{i^{\prime}}$ with probability $1$.
Therefore, we have

(28) $\displaystyle\sum_{i=1}^{m}\big{(}f_{\mathcal{A}_{j}}(w_{i})\cdot
p^{j}_{i^{\prime}}(w_{i})\big{)}$
$\displaystyle=\overline{F}_{\mathcal{A}_{j}}(w_{i^{\prime}})\cdot
w_{i^{\prime}},$

where $p^{j}_{i^{\prime}}\coloneqq w_{i}\cdot
x^{j}_{i^{\prime}}(w_{i})-\int_{w=0}^{w_{i}}x^{j}_{i^{\prime}}(w)\differential{w}$
is the price associated with the buyer with valuation $w_{i}$ and deadline $j$
given the allocation curve $x^{j}_{i^{\prime}}(\cdot)$.

Let
$\delta^{j}_{i^{\prime}}=\hat{x}_{j}(w_{i^{\prime}})-\hat{x}_{j}(w_{i^{\prime}-1})=\hat{x}_{i^{\prime}j}-\hat{x}_{(i^{\prime}-1)j}$
for all $j\in[1,k]$ and $i^{\prime}\in[1,m]$. Then for all $i\in[1,m]$ we can
write

$\hat{x}_{j}(w_{i})=\sum\limits_{i^{\prime}=1}^{m}\big{(}\delta^{j}_{i^{\prime}}\cdot
x^{j}_{i^{\prime}}(w_{i})\big{)},\quad\hat{p}_{j}(w_{i})=\sum\limits_{i^{\prime}=1}^{m}\big{(}\delta^{j}_{i^{\prime}}\cdot
p^{j}_{i^{\prime}}(w_{i})\big{)},$

where $\hat{p}_{j}(w_{i})\coloneqq
w_{i}\cdot\hat{x}_{j}(w_{i})-\int_{w=0}^{w_{i}}\hat{x}(w)\differential{w}$.
Therefore, any optimal auction for D4$(\mathcal{A})$ raises a revenue of

$\displaystyle\sum\limits_{i=1}^{m}\big{(}f_{\mathcal{A}_{j}}(w_{i})\cdot\hat{p}_{j}(w_{i})\big{)}$
$\displaystyle=\sum\limits_{i=1}^{m}\Big{(}f_{\mathcal{A}_{j}}(w_{i})\cdot\sum\limits_{i^{\prime}=1}^{m}\big{(}\delta^{j}_{i^{\prime}}\cdot
p^{j}_{i^{\prime}}(w_{i})\big{)}\Big{)}$
$\displaystyle=\sum\limits_{i^{\prime}=1}^{m}\Big{(}\delta^{j}_{i^{\prime}}\cdot\sum\limits_{i=1}^{m}\big{(}f_{\mathcal{A}_{j}}(w_{i})\cdot
p^{j}_{i^{\prime}}(w_{i})\big{)}\Big{)}=\sum\limits_{i^{\prime}=1}^{m}\Big{(}\delta^{j}_{i^{\prime}}\cdot\overline{F}_{\mathcal{A}_{j}}(w_{i^{\prime}})\cdot
w_{i^{\prime}}\Big{)}$

for the buyer with deadline $j$, where the last equality follows from Eq.
(28). Summing this up for all $j\in[1,k]$ (and weighting each term by
$\Pr_{(v,d)\sim\mathcal{A}}[d=j]$) gives the expression in the theorem.

Consider any arbitrary $(w_{a},r)\in\textsf{LE}(\mathcal{A})$. By property 2
in Theorem B.6, we have $\hat{x}_{ij}=\hat{x}_{ir}$ for all $i\in[0,a]$ and
$j\in[r,k]$. Therefore for all $j\in[r,k]$ we have
$\delta^{j}_{a}=\hat{x}_{aj}-\hat{x}_{(a-1)j}=\hat{x}_{ar}-\hat{x}_{(a-1)r}=\delta^{r}_{a}$.

Next, consider all consecutive
$(w_{a},r),(w_{b},r^{\prime})\in\textsf{LE}(\mathcal{A})$ such that $a<b$,
$r^{\prime}\geq r$. Since $f_{\mathcal{A}_{j}}(w_{i})=0$ for all $i\in(a,b)$
and $j\geq r$, by property 3 in Theorem B.6, we have
$\hat{x}_{ij}=\hat{x}_{aj}$ for all $i\in[a,b)$ and $j\in[r,k]$. Hence for all
$i\in(a,b)$ and $j\in[r,k]$ we have
$\delta^{j}_{i}=\hat{x}_{ij}-\hat{x}_{(i-1)j}=\hat{x}_{aj}-\hat{x}_{aj}=0$.

Finally, notice that $\hat{x}_{0k}=0$ and $\hat{x}_{mk}=1$. Therefore we have

$\displaystyle\sum\limits_{i^{\prime}=1}^{m}\delta^{k}_{i^{\prime}}$
$\displaystyle=\sum\limits_{i^{\prime}=1}^{m}\big{(}\hat{x}_{i^{\prime}k}-\hat{x}_{(i^{\prime}-1)k}\big{)}=\hat{x}_{mk}-\hat{x}_{0k}=1.$

The last property then follows from observing that $\hat{x}_{ik}=0$ for every
$i$ such that $\mathbbm{1}^{\textsf{LE}(\mathcal{A})}_{i}=0$. ∎


# Graph Learning from Gaussian and Stationary Graph Signals

###### Abstract

Graphs have become pervasive tools to represent information and datasets with
irregular support. However, in many cases, the underlying graph is either
unavailable or naively obtained, calling for more advanced methods to its
estimation. Indeed, graph topology inference methods that estimate the network
structure from a set of signal observations have a long and well established
history. By assuming that the observations are both Gaussian and stationary in
the sought graph, this paper proposes a new scheme to learn the network from
nodal observations. Consideration of graph stationarity overcomes some of the
limitations of the classical Graphical Lasso algorithm, which is constrained
to a more specific class of graphical models. On the other hand, Gaussianity
allows us to regularize the estimation, requiring less samples than in
existing graph stationarity-based approaches. While the resultant estimation
(optimization) problem is more complex and non-convex, we design an
alternating convex approach able to find a stationary solution. Numerical
tests with synthetic and real data are presented, and the performance of our
approach is compared with existing alternatives.

Index Terms—  Network topology inference, Gaussian Markov random fields, graph
stationarity, graphical lasso, graphical models.

##  1 Introduction

Data defined over non-Euclidean supports is nowadays ubiquitous, with graphs,
due to their versatility and tractability, having emerged as a suitable
alternative to represent those irregular domains. This has prompted the
emergence of works investigating how to generalize classical information
processing architectures to graph-based signals. The effort has permeated also
to applications, with a number of graph-based data-processing schemes being
successfully applied to electrical, communication, social, geographic,
financial, genetic, and brain networks [1, 2, 3, 4, 5], to name a few.

When dealing with network data, the workhorse assumption is that the
underlying graph is known. Unfortunately, this is not always the case. In many
scenarios, the graph structure is not well defined, either because there is no
underlying physical network (e.g., when the graph models logical relations
between variables) or because there is no unique metric to measure the
relationship between the different nodes of the network. Since in most cases
such relationships cannot be known, the graph topology is usually inferred
from the set of available observations. To estimate the graph topology, a
first step is to formalize the relation between the network edges and the
properties of the signals that live in the network. The early graph topology
inference methods [6, 7] adopted a statistical approach, such as the
correlation network [1, Ch. 7.3.1], partial correlations or Gaussian Markov
Random Fields (MRFs) [8, 1], with graphical Lasso (GL) [9] being the workhorse
scheme for their estimation. In recent years, graph signal processing (GSP)
based models exploiting properties such as smoothness or graph stationarity
[10, 11, 12, 13] have become increasingly important, along with their
nonlinear generalizations[14]. The main advantage of heuristic pairwise
distance-based and correlation-based approaches, which infer each link
separately, is that the graph can be learned even in setups where only a few
observations are available. In contrast, theoretically sounded methods that
postulate a more advanced relation (mapping) between the observations and the
totality of the graph require a significantly larger number of observations.
Remarkably, GL achieves a sweet spot in this regard, postulating a relatively
sophisticated model (a Gaussian MRF) and being able to learn the graph using
only a moderate number of observations.

In this context, this paper proposes a new graph-learning method to estimate a
network from a set of signal observations. Our main contributions are: i)
building a new graph-learning approach around the assumption that the
observations at hand are both Gaussian and stationary in the graph; ii)
formulating a joint optimization problem that yields as a result the sough
graph along with an enhanced estimate of the precision matrix of the Gaussian
process; and iii) designing an efficient provably convergent algorithm to
address the (non-convex) optimization in ii). Relative to GL, which assumes
Gaussianity but not stationarity, we end up with a more general graph learning
algorithm able to recover the graph in a wider range of scenarios and that
includes GL as a particular case. For our problem formulation we made
assumptions related to 1) the graph (to be sparse), 2) the signals (to be
Gaussian) and 3) the relationship between the graph and the signals (to be
graph stationary). Relative to graph-stationary-based approaches, the scheme
in this paper is able to work with less samples, is more robust to noise and
perturbations, and can be used to draw links with classical statistical
approaches. The price to pay is that the resultant algorithm incurs higher
computational complexity and that, due to non-convexity, only convergence to a
stationary point is guaranteed. Interestingly, the numerical results
presented, while limited, show that our algorithm either outperforms or
matches the results achieved by its convex counterparts.

Outline. The remaining of the paper is organized as follows. Sec. 2 surveys
the required graph background, making special emphasis on the concepts of
Gaussian and graph stationary signals. In Sec. 3 we formalize the problem of
network topology inference from signals that are Gaussian and graph
stationary. Sec. 4 presents a bi-convex optimization formulation of the
problem in Sec. 3 and develops a provably convergent alternating iterative
algorithm to solve it. Sec. 5 illustrates the recovery performance of the
proposed algorithm using both synthetic and real data.

##  2 Fundamentals of graph signal processing

This section introduces GSP concepts that help to formalize the learning
problem and explain the relation between the available signals and the
topology of the underlying graph.

GSO and graph signals. Let ${\mathcal{G}}=\\{{\mathcal{V}},{\mathcal{E}}\\}$
denote an _undirected_ and weighted graph with $N$ nodes where ${\mathcal{V}}$
and ${\mathcal{E}}$ represent the vertex and edge set, respectively. For any
${\mathcal{G}}$, its weighted adjacency matrix
${\mathbf{A}}\in{\mathbb{R}}^{N\times N}$ is a sparse matrix with $A_{ij}$
capturing the weight of the edge between the nodes $i$ and $j$. A general
representation of the graph is the graph-shift operator GSO
${\mathbf{S}}\in\mathbb{R}^{N\times N}$, where $S_{ij}\neq 0$ if and only if
$i=j$ or $(i,j)\in{\mathcal{E}}$ [4]. Typical choices for the GSO are the
adjacency matrix ${\mathbf{A}}$, the combinatorial graph Laplacian
${\mathbf{L}}:=\text{diag}({\mathbf{A}}\textbf{1})-{\mathbf{A}}$, and their
respective generalizations. The signals defined on the nodes of the graph
${\mathcal{G}}$ are called graph signals. A graph signal can be denoted as a
vector ${\mathbf{x}}\in\mathbb{R}^{N}$ with $x_{i}$ being the signal value
observed at node $i$.

Gaussian graph signals. Since graph signals can be represented by
$N$-dimensional vectors, the random (zero-mean) graph signal ${\mathbf{x}}$ is
said to be Gaussian, if its joint distribution has the form

$f_{\bm{\Theta}}({\mathbf{x}})={(2\pi)^{-N/2}}\cdot\det(\bm{\Theta})^{\frac{1}{2}}\cdot\exp\big{(}{-\frac{1}{2}{\mathbf{x}}^{T}\bm{\Theta}{\mathbf{x}}}\big{)}.\vspace{-0.2cm}$
(1)

with $\bm{\Theta}$ being the precision matrix (the inverse of the covariance
matrix $\bm{\Sigma}$). Suppose now that we have a collection of $R$ Gaussian
signals
${\mathbf{X}}=[{\mathbf{x}}_{1},...,{\mathbf{x}}_{R}]\in\mathbb{R}^{N\times
R}$ each of them independently drawn from the distribution in (1). The log-
likelihood associated with ${\mathbf{X}}$ is

$L({\mathbf{X}}|\bm{\Theta})=\prod_{r=1}^{R}f_{\bm{\Theta}}({\mathbf{x}}_{r}),\
\ \
{\mathcal{L}}({\mathbf{X}}|\bm{\Theta})=\sum_{r=1}^{R}\log(f_{\bm{\Theta}}({\mathbf{x}}_{r})).\vspace{-0.1cm}$
(2)

This expression will be exploited when formulating our proposed inference
approach and establishing links with classical methods.

Graph filters. A versatile tool to model the relationship between the signal
${\mathbf{x}}$ and its underlying graph are graph filters. A graph filter
${\mathbf{H}}\in{\mathbb{R}}^{N\times N}$ is a linear operator that accounts
for the topology of the graph and is defined as a polynomial of the GSO
${\mathbf{S}}$ of the form

${\mathbf{H}}=\sum_{l=0}^{L-1}h_{l}{\mathbf{S}}^{l}=\sum_{l=0}^{L-1}h_{l}{\mathbf{V}}\bm{\Lambda}^{l}{\mathbf{V}}^{\top}={\mathbf{V}}\Big{(}\sum_{l=0}^{L-1}h_{l}\bm{\Lambda}^{l}\Big{)}{\mathbf{V}}^{\top},\vspace{-0.1cm}$
(3)

where $L-1$ is the filter degree, $\\{h_{l}\\}_{l=0}^{L-1}$ represent the
filter coefficients, and ${\mathbf{V}}$ and $\bm{\Lambda}$ are the
eigenvectors and eigenvalues of the GSO respectively. Since ${\mathbf{H}}$ is
a polynomial of ${\mathbf{S}}$, it readily follows that both matrices have the
same eigenvectors ${\mathbf{V}}$.

Stationary graph signals. A random graph signal ${\mathbf{x}}$ is stationary
in the GSO ${\mathbf{S}}$ if it can be represented as the output generated by
a graph filter ${\mathbf{H}}$ [cf.(3)] whose input is a zero mean white signal
${\mathbf{w}}\in\mathbb{R}^{N}$ with
$\mathbb{E}[{\mathbf{w}}{\mathbf{w}}^{\top}]={\mathbf{I}}$. Under this
definition, it follows that if ${\mathbf{x}}$ is stationary and
${\mathcal{G}}$ is undirected, the covariance of ${\mathbf{x}}$ is given by

$\bm{\Sigma}\\!=\\!\mathbb{E}[{\mathbf{x}}{\mathbf{x}}^{\\!\top}\\!]\\!=\\!\mathbb{E}[{\mathbf{H}}{\mathbf{w}}{\mathbf{w}}^{\\!\top}{\mathbf{H}}^{\\!\top}]\\!=\\!{\mathbf{H}}\mathbb{E}[{\mathbf{w}}{\mathbf{w}}^{\\!\top}\\!]{\mathbf{H}}^{\\!\top}\\!=\\!{\mathbf{H}}{\mathbf{H}}^{\\!\top}\\!=\\!{\mathbf{H}}^{2},\vspace{-0.1cm}$
(4)

where the last equality follows from the symmetry of the GSO. Since
${\mathbf{H}}$ is a polynomial in ${\mathbf{S}}$, so is $\bm{\Sigma}$. This
implies that: i) $\bm{\Theta}$ is also a polynomial of ${\mathbf{S}}$, and ii)
${\mathbf{S}}$, $\bm{\Sigma}$ and $\bm{\Theta}$ share the same eigenvectors
${\mathbf{V}}$ [15, 16, 17]. As a byproduct, we also have that matrices
${\mathbf{S}}$ and $\bm{\Sigma}$ commute, which are properties to be exploited
later on.

$\displaystyle\vspace{-0.1cm}\bm{\Theta}_{1}^{(t+1)}$
$\displaystyle=\operatornamewithlimits{argmin}_{\bm{\Theta}_{1}\succeq
0}\text{tr}(\hat{\bm{\Sigma}}\bm{\Theta}_{1})+\frac{\lambda_{1}}{2}\|\bm{\Theta}_{1}{\mathbf{S}}^{(t)}\\!-\\!{\mathbf{S}}^{(t)}\bm{\Theta}_{1}\|_{F}^{2}+\frac{\lambda_{2}}{2}\|\bm{\Theta}_{1}\\!-\\!\bm{\Theta}_{2}^{(t)}\|_{F}^{2}$
(10) $\displaystyle\vspace{-0.1cm}\bm{\Theta}_{2}^{(t+1)}$
$\displaystyle=\operatornamewithlimits{argmin}_{\bm{\Theta}_{2}\succeq 0}\
-\log\det(\bm{\Theta}_{2})+\frac{\lambda_{2}}{2}\|\bm{\Theta}_{1}^{(t+1)}-\bm{\Theta}_{2}\|_{F}^{2}$
(11) $\displaystyle\vspace{-0.1cm}{\mathbf{S}}^{(t+1)}$
$\displaystyle=\operatornamewithlimits{argmin}_{{\mathbf{S}}\in{\mathcal{S}}}\
\rho\|{\mathbf{S}}\|_{1}+\frac{\lambda_{1}}{2}\|\bm{\Theta}_{1}^{(t+1)}{\mathbf{S}}-{\mathbf{S}}\bm{\Theta}_{1}^{(t+1)}\|_{F}^{2}\vspace{-0.4cm}$
(12)

##  3 Graph learning problem formulation

To formally state the graph learning problem, let
${\mathbf{X}}:=[{\mathbf{x}}_{1},...,{\mathbf{x}}_{R}]\in{\mathbb{R}}^{N\times
R}$ denote a collection of $R$ independent realizations of a random process
defined on the $N$ nodes of the unknown graph ${\mathcal{G}}$. The goal of
graph learning problems is to find the optimal graph descriptor (i.e., the GSO
${\mathbf{S}}$ associated with ${\mathcal{G}}$) given the set of signals
${\mathbf{X}}$. There exist a range of approaches in the literature to address
this problem [9, 18, 15]. In this paper, we approach it by assuming that the
signals in ${\mathbf{X}}$ are _simultaneously_ Gaussian and stationary in the
sought graph. Based on these two assumptions (along with the standard
consideration of the graph being sparse), we propose a new approach more
robust and general than those existing in the literature, able to recover the
graph in a wider range of scenarios. This section starts by formalizing the
problem we want to address and then, formulates it as an optimization problem.
After that, we discuss the main features of our approach and elaborate in
detail on the differences between GL and graph-stationary approaches.

Problem 1. Given the set of signals ${\mathbf{X}}\in{\mathbb{R}}^{N\times R}$,
find the underlying sparse graph structure encoded in ${\mathbf{S}}$ under the
assumptions:

(AS1): The graph signals in ${\mathbf{X}}$ are i.i.d. realizations of
${\mathcal{N}}({\mathbf{0}},\bm{\Sigma})$

(AS2): The graph signals in ${\mathbf{X}}$ are stationary in ${\mathbf{S}}$.

Our approach is to recast Problem 1 as the following optimization

$\displaystyle\\!\hat{\bm{\Theta}},{\hat{\mathbf{S}}}=\operatornamewithlimits{argmin}_{\bm{\Theta}\succeq
0,{\mathbf{S}}\in{\mathcal{S}}}\ $
$\displaystyle-\log(\det(\bm{\Theta}))+\text{tr}(\hat{\bm{\Sigma}}\bm{\Theta})+\rho\|{\mathbf{S}}\|_{0}.$
$\displaystyle\\!\mathrm{\;\;s.\;to}$
$\displaystyle\bm{\Theta}{\mathbf{S}}={\mathbf{S}}\bm{\Theta},\vspace{-0.6cm}$
(5)

where ${\mathcal{S}}$ is a generic set representing additional constraints
that ${\mathbf{S}}$ is known to satisfy (e.g., the GSO being symmetric on its
entries being between zero and one). The minimization takes as input the
sample covariance matrix
$\hat{\bm{\Sigma}}=\frac{1}{R}{\mathbf{X}}{\mathbf{X}}^{\mathsf{T}}$ and
generates as output the estimate for ${\mathbf{S}}$ and, as a byproduct, the
estimate for $\bm{\Theta}$. Next, we explain the motivation for each of the
terms in (3) and, especially, for the constraint
$\bm{\Theta}{\mathbf{S}}={\mathbf{S}}\bm{\Theta}$, which is critical in our
approach. The first two terms in the objective function come from leveraging
(AS1) and are the result of minimizing the negative log-likelihood expression
in (2). In contrast, the term $\rho\|{\mathbf{S}}\|_{0}$ accounts for the fact
of ${\mathbf{S}}$ being a GSO (hence, sparse) with $\rho$ being a
regularization constant / multiplier/ prior related to the sparsity level of
the graph. Finally, the equality constraint codifies (AS2). Note that the
polynomial relation between $\bm{\Sigma}$ and ${\mathbf{S}}$ implied by (AS2)
is typically codified in estimation/optimization problems by either: i)
extracting the eigenvectors of $\bm{\Sigma}$ and imposing than those should be
the eigenvectors of ${\mathbf{S}}$ [12] or ii) imposing the constraint
$\bm{\Sigma}{\mathbf{S}}={\mathbf{S}}\bm{\Sigma}$ [13]. Differently, here we
codify the polynomial relation implied by (AS2) by imposing commutativity
between $\bm{\Theta}$ and ${\mathbf{S}}$. This is not only equivalent (if the
matrices are full rank and $\bm{\Sigma}$ and ${\mathbf{S}}$ commute, the same
is true for $\bm{\Theta}$ and ${\mathbf{S}}$ –for more details check the
Cayley-Hamilton theorem), but also more convenient when dealing with Gaussian
signals, since their likelihood is convex with respect to $\bm{\Theta}$. It is
important to note that while the assumption of stationarity may appear strong,
it is in fact more lenient (more general) than an MRF, as it allows the
precision matrix to be a polynomial function of ${\mathbf{S}}$. This stands in
contrast to the GL framework, which limits this mapping to the equality
constraint, ${\mathbf{S}}=\bm{\Theta}$.

From a conceptual point of view, our formulation reaches a sweet spot between
GL and graph-stationarity approaches. The price to pay is that the resultant
algorithm (even if the $\ell_{0}$ norm is relaxed) is not convex due to the
bilinear constraint coupling the _optimization_ variables $\bm{\Theta}$ and
${\mathbf{S}}$. The remaining of this section provides additional details
about the relation between (3), GL and graph-stationary approaches, while Sec.
4 designs an efficient algorithm to deal with the non-convexities.

Suppose first that we simplify (AS2) imposing that
$\sigma{\mathbf{I}}+\delta{\mathbf{S}}=\bm{\Sigma}^{-1}$. Then, up to the
diagonal values and scaling issues, the sparse matrix ${\mathbf{S}}$ to be
estimated and $\bm{\Theta}=\bm{\Sigma}^{-1}$ are the same and, as a result, it
suffices to optimize over one of them. Upon defining
${\mathcal{S}}_{\bm{\Theta}}:=\\{\bm{\Theta}=\bm{\Theta}^{T};\
\bm{\Theta}\succeq 0\\}$, this leads to

$\displaystyle\\!\hat{\bm{\Theta}}=\operatornamewithlimits{argmin}_{\bm{\Theta}\succeq
0,\bm{\Theta}\in{\mathcal{S}}_{\bm{\Theta}}}\ $
$\displaystyle-\log(\det(\bm{\Theta}))+\text{tr}(\hat{\bm{\Sigma}}\bm{\Theta})+\rho\|\bm{\Theta}\|_{0},\vspace{-0.3cm}$
(6)

which coincides exactly with the GL formulation [9]. The main advantages of
(6) relative to (3) are that the number of variables is smaller and the
resultant problem (after relaxing the $\ell_{0}$ norm) is convex. The main
drawback is that by forcing the support of ${\mathbf{S}}$ and $\bm{\Theta}$ to
be the same, the set of feasible graphs (and their topological properties) is
more limited. Remarkably, when the model assumed in (6) holds true (i.e., data
is Gaussian and $\bm{\Theta}$ is sparse), GL is able to find reliable
estimates of ${\mathbf{S}}$ even when the number of samples $R$ is fairly low.
On the other hand, simulations will show that GL does a poor job estimating
${\mathbf{S}}$ when the relation between the precision matrix and
${\mathcal{G}}$ is more involved.

Suppose now that we drop (AS1), so that we approach the graph learning problem
assuming only that signals are stationary in ${\mathcal{G}}$. Since
$\bm{\Sigma}$ shares the same eigenvectors, if the covariance if full rank the
expression $\bm{\Theta}{\mathbf{S}}={\mathbf{S}}\bm{\Theta}$ is equivalent to
requiring $\bm{\Sigma}{\mathbf{S}}={\mathbf{S}}\bm{\Sigma}$. This allows us to
rewrite (3) as

$\displaystyle\\!{\hat{\mathbf{S}}}=\operatornamewithlimits{argmin}_{{\mathbf{S}}\in{\mathcal{S}}}\
\rho\|{\mathbf{S}}\|_{0}\;\;\;\mathrm{\;\;s.\;to}\;\;\;\hat{\bm{\Sigma}}{\mathbf{S}}={\mathbf{S}}\hat{\bm{\Sigma}},\vspace{-0.2cm}$
(7)

where the constraint is typically relaxed as
$\|\hat{\bm{\Sigma}}{\mathbf{S}}-{\mathbf{S}}\hat{\bm{\Sigma}}\|_{F}^{2}\leq\epsilon$
to account for the fact that we have $\hat{\bm{\Sigma}}\approx\bm{\Sigma}$.
While more general than either (3) or (6), the absence of Gaussianity implies
that (7) being able to identify the ground truth ${\mathbf{S}}$ requires very
reliable estimates of $\hat{\bm{\Sigma}}$, so that $\epsilon$ can be set close
to zero. This is indeed a challenge, especially in setups when the number of
nodes is large.

##  4 Graph learning algorithm

The formulation in (3), which optimizes jointly over $\bm{\Theta}$ and
${\mathbf{S}}$, is non-convex due to the presence of the $\ell_{0}$ norm and
the bilinear constraint. In this section, we design an algorithm to solve a
relaxed version of (3) in an efficient way. The most important design steps
include: i) relaxing the $\ell_{0}$ norm in the objective with the $\ell_{1}$
norm; ii) renaming variable $\bm{\Theta}$ as $\bm{\Theta}_{1}$, creating an
auxiliary variable $\bm{\Theta}_{2}$, and adding a new constraint that
enforces $\bm{\Theta}_{1}=\bm{\Theta}_{2}$; iii) augmenting the cost with
quadratic terms that account for the constraints; and iv) proposing a block
coordinate descent framework that optimizes each of the blocks of variables
($\bm{\Theta}_{1}$, $\bm{\Theta}_{2}$, and ${\mathbf{S}}$) in an iterative
manner. Specifically, after incorporating steps i) and ii), (3) is re-written
as

$\displaystyle\\!\hat{\bm{\Theta}}_{1},\hat{\bm{\Theta}}_{2},{\hat{\mathbf{S}}}=$
$\displaystyle\\!\\!\\!\\!\\!\operatornamewithlimits{argmin}_{\bm{\Theta}_{1},\bm{\Theta}_{2}\succeq
0;{\mathbf{S}}\in{\mathcal{S}}}\
-\log(\det(\bm{\Theta}_{1}))+\text{tr}(\hat{\bm{\Sigma}}\bm{\Theta}_{2})+\rho\|{\mathbf{S}}\|_{1},$
$\displaystyle\\!\mathrm{\;\;s.\;to}$
$\displaystyle\bm{\Theta}_{2}{\mathbf{S}}={\mathbf{S}}\bm{\Theta}_{2}\;\;\text{and}\;\;\bm{\Theta}_{1}=\bm{\Theta}_{2}.\vspace{-0.4cm}$
(8)

Input: $\hat{\bm{\Sigma}}$

Outputs : $\hat{\bm{\Theta}}_{1}$, $\hat{\bm{\Theta}}_{2}$, and
${\hat{\mathbf{S}}}$

1

2Initialize $\hat{\bm{\Theta}}_{2}^{(0)}$, ${\hat{\mathbf{S}}}^{(0)}$, $=0$.

3

4while  _not converged_ do

5

6 Update $\hat{\bm{\Theta}}_{1}^{(t+1)}$ by solving (10) (Step I).

7

8 Update $\hat{\bm{\Theta}}_{2}^{(t+1)}$ by solving (11) (Step II).

9

10 Update ${\hat{\mathbf{S}}}^{(t+1)}$ by solving (12) (Step III).

11

12 $t=t+1$

13

14 end while

$\hat{\bm{\Theta}}_{1}=\hat{\bm{\Theta}}_{1}^{(t)}$,$\hat{\bm{\Theta}}_{2}=\hat{\bm{\Theta}}_{2}^{(t)}$,${\hat{\mathbf{S}}}={\hat{\mathbf{S}}}^{(t)}$

Algorithm 1 BSUM graph-learning method for Gaussian and stationary signals
(BSUM-GGSR)

Next, to deal with equality constraints
$\bm{\Theta}_{2}{\mathbf{S}}={\mathbf{S}}\bm{\Theta}_{2}$ and
$\bm{\Theta}_{1}=\bm{\Theta}_{2}$, we augment the objective in (4) with
quadratic terms whose weights are controlled by regularization parameters
$\lambda_{1}$ and $\lambda_{2}$, yielding

$\displaystyle\\!\hat{\bm{\Theta}}_{1},\hat{\bm{\Theta}}_{2},{\hat{\mathbf{S}}}=$
$\displaystyle\operatornamewithlimits{argmin}_{\bm{\Theta}_{1},\bm{\Theta}_{2}\succeq
0;{\mathbf{S}}\in{\mathcal{S}}}\text{tr}(\hat{\bm{\Sigma}}\bm{\Theta}_{1})-\log(\det(\bm{\Theta}_{2}))+\rho\|{\mathbf{S}}\|_{1}$
$\displaystyle+\frac{\lambda_{1}}{2}\|\bm{\Theta}_{1}{\mathbf{S}}-{\mathbf{S}}\bm{\Theta}_{1}\|_{F}^{2}+\frac{\lambda_{2}}{2}\|\bm{\Theta}_{1}-\bm{\Theta}_{2}\|_{F}^{2}.$
(9)

Fig. 1: (Left) Ratio of recovered graphs while increasing the number of
signals $R$ used to estimate $\hat{\bm{\Sigma}}$. The values in the y-axis
correspond to averages across $100$ realizations of Erdös-Rényi graphs with
$N=20$ and link probability $p=0.1$. (Center) Setup similar to the previous
one, but setting $R=10^{6}$ and with the x-axis representing the normalized
power noise present in the observations. (Right) Visual representations of
recovered graphs using as input financial data signals of 40 companies from 4
different sectors.

The pseudocode describing how to handle the optimization of (9) is detailed in
Alg. 1. The overall idea is to solve the problem via a Block Successive Upper
Bound Minimization (BSUM) alternating minimization approach for the three
blocks of variables. The iterative algorithm, referred to as graph Gaussian
stationary recovery (GGSR), consists in successively optimizing the three
blocks of variables with the other two blocks fixed. Note that, thanks to the
auxiliary variables $\bm{\Theta}_{1}$ and $\bm{\Theta}_{2}$, the optimization
for each of the three blocks is relatively simple and easy to solve (details
are omitted due to the page limit constraints, but will be provided in the
journal version of this work).

Convergence guarantees of the iterative algorithm to a stationary point of
function in (9) are formally stated next.

Proposition 1. Denote with $f$ the objective function in (4). Let
${\mathcal{Y}}^{*}$ be the set of stationary points of (4), and let
${\mathbf{y}}^{(t)}=[\operatorname{vec}(\bm{\Theta}_{1}^{(t)})^{\top}\\!,$
$\operatorname{vec}(\bm{\Theta}_{2}^{(t)})^{\top}\\!,\operatorname{vec}({\mathbf{S}}^{(t)})^{\top}]^{\top}$
be the solution generated after running the three steps in Table 1 $t$ times.
Then, the solution generated by the iterative algorithm presented in Table 1
converges to a stationary point of $f$ as $t$ goes to infinity, i.e.,

$\lim_{t\to\infty}d({\mathbf{y}}^{(t)},{\mathcal{Y}}^{*})=0,$

with
$d({\mathbf{y}},{\mathcal{Y}}^{*}):=\min_{{\mathbf{y}}^{*}\in{\mathcal{Y}}^{*}}\|{\mathbf{y}}-{\mathbf{y}}^{*}\|_{2}$.

The proof relies on the convergence results shown in [19, Th. 1b]. It can be
verified that our iterative algorithm meets the required conditions for the
convergence to a stationary point, which are the following: 1)
$f({\mathbf{y}})$ is regular for all feasible points; 2) the first-order
behavior of each step (when the remaining steps are fixed) is the same as
$f({\mathbf{y}})$; 3) the set
${\mathcal{Y}}^{(0)}=\\{{\mathbf{y}}|f({\mathbf{y}})\leq
f({\mathbf{y}})^{(0)}\\}$ is compact; and 4) Steps I and II have a unique
solution. Note that the convergence of the BSUM-GGSR algorithm was not obvious
since the original function in (4) is non-convex and Step III has no unique
solution, we refer readers to [20] for more details.

The computational complexity required by Alg. 1 is polynomial but, in the
worst case, it can scale as $O(N^{7})$, which may hinder its application to
large graphs. Although the implementation of more efficient and scalable
algorithms is beyond the scope of this conference paper, we find it relevant
and it is subject of our current work.

##  5 Numerical experiments

This section numerically assesses the performance of our proposed algorithm,
GGSR, comparing it with the two most related alternatives in literature: GL
[9] and GSR [12]. Fig. 1 summarizes the results for synthetic data (left and
center panels) and real data (right panel). For the synthetic experiments, the
network topology is considered to be successfully recovered (score 1) if each
and every link is successfully estimated and zero otherwise. This score is
averaged over $100$ realizations of random graphs with $N=20$ nodes following
an Erdös-Rényi model with link probability $p=0.1$. Regarding the generation
of the graph signals, two different models for the covariance matrices have
been studied: 1) ${\mathbf{C}}_{poly}={\mathbf{H}}^{2}$ with
${\mathbf{H}}=\sum_{l=0}^{L}h_{l}{\mathbf{S}}^{l}$ being a random polynomials
in ${\mathbf{S}}$ and $h_{l}\sim{\mathcal{N}}(0,1)$; and 2)
${\mathbf{C}}_{MRF}=(\sigma{\mathbf{I}}+\delta{\mathbf{S}})^{-1}$ with
$\sigma$, $\delta>0$ selected to guarantee that ${\mathbf{C}}_{MRF}^{-1}$ is
positive definite. Finally, three different algorithms are compared: (i) “GL”,
which assumes that the signals are Gaussian in the graph; (ii) “GSR”, which
assumes that the signals are stationary in the graph; and (iii) “GGSR”, which
assumes both stationarity and Gaussianity for the graph signals.

Performance for limited number of samples. In this experiment, we quantify the
recovery performance of the different algorithms as the number of signals $R$
increases (left panel of Fig. 1). Moreover, we consider the two covariance
models mentioned above: ${\mathbf{C}}_{MRF}$ and ${\mathbf{C}}_{poly}$.
Looking at the results obtained when using the ${\mathbf{C}}_{MRF}$ covariance
model, we observe that GGSR provides a ratio of recovered graphs very similar
to GL (which is the maximum-likelihood estimator under ${\mathbf{C}}_{MRF}$)
and clearly outperforms GSR. Regarding the results when using the covariance
model ${\mathbf{C}}_{poly}$ it can be seen that GGSR outperforms both GSR and
GL, the latter being unable to estimate the topology of any graph. In summary,
when considering the ${\mathbf{C}}_{MRF}$ model, GL offers the best estimation
because its assumption fits the ${\mathbf{C}}_{MRF}$ covariance model
perfectly and GGSR, which incurs higher computational complexity, obtains very
similar results. Differently, a big performance gap is observed under the
${\mathbf{C}}_{poly}$ data generation model, where GGSR achieves a recovery
performance similar to that of GSR while requiring 10 times less samples. This
experiment validates our claims (laid out in Sec. 3) about the ability of GGSR
to find a “sweet spot” between GL and GSR.

Performance for noisy samples. In most setups, the signals at hand may contain
a certain amount of noise. In this test case, we evaluate the robustness of
the different algorithms for learning the graph topology from a set of noisy
samples, keeping $R$ constant and increasing the noise present in
${\mathbf{X}}$ (center panel of Fig. 1). For the ${\mathbf{C}}_{MRF}$ model,
we observe that the performance is fairly robust for all algorithms, showing
only a slight deterioration as the noise increases. The two main findings from
the curves associated with the ${\mathbf{C}}_{poly}$ model are: i) the impact
of the noise is much more noticeable than for ${\mathbf{C}}_{MRF}$ and ii)
GGSR is less sensitive to noise than GSR, with GGSR being able to correctly
recover approximately twice as many graphs as GSR.

Learning clustered graph from financial stocks. Finally, we test our algorithm
on a dataset with real data. We selected 40 companies from 4 different sectors
of the S&P 500 index (10 companies from each sector). We collected the daily
log returns of each company in the period Jan. 1st, 2010 to Jan. 1st, 2016
totaling 1510 days worth of data. Our goal for this experiment is to estimate
the graph structure from ${\mathbf{X}}\in{\mathbb{R}}^{40\times 1510}$ using
GL, GSR and GGSR. The first column of the right panel in Fig. 1 represents the
graph estimation as the thresholded version of the correlation matrix obtained
from the data. The remaining 3 columns correspond to the estimation performed
by GSR, GL and GGSR, respectively. The first row of the right panel in Fig. 1
represents the graph estimated by each of the 4 algorithms, with the different
colors representing the sector to which each company belongs. Since we do not
have a ground truth graph, to asses the quality of the learned graphs, we
leverage the fact that the companies belonged to 4 financial sectors and run a
clustering algorithm using as input each of the learned graphs. The second row
of the right panel in Fig. 1 represents the output of a spectral clustering
algorithm carried out over each of the graphs in the first row and setting the
number of clusters to 4. From the results, we observe that the clustering
based on the graph estimated by GGSR is the one leading to less clustering
errors. In addition, note that GGSR shows higher intra-cluster connections and
lower inter-cluster connections compared to the other algorithms.

##  6 Conclusions

This paper addressed the problem of estimating a graph from nodal observations
under the assumption that the observed signals were Gaussian and graph
stationary. Relative to classical Gaussian Markovian graphical models, the
stationary assumption covers a wider range of scenarios where GL fails to
obtain the true graph. The estimation problem was formulated as a constraint
optimization where both the graph and the precision matrix of the Gaussian
distribution were optimization variables, resulting in a non-convex
formulation. An alternating method with convergence guarantees was designed
and the recovery performance was assessed using synthetic and real data.

## References

  * [1] E. D. Kolaczyk,  Statistical Analysis of Network Data: Methods and Models,  Springer, New York, NY, 2009.

  * [2] O. Sporns,  Discovering the Human Connectome,  MIT Press, Boston, MA, 2012.
  * [3] D.I. Shuman, S.K. Narang, P. Frossard, A. Ortega, and P. Vandergheynst,  “The emerging field of signal processing on graphs: Extending high-dimensional data analysis to networks and other irregular domains,”  IEEE Signal Process. Mag., vol. 30, no. 3, pp. 83–98, May 2013\.
  * [4] A. Ortega, P. Frossard, J. Kovacevic, J. M. F. Moura, and P. Vandergheynst,  “Graph signal processing: Overview, challenges, and applications,”  Proc. IEEE, vol. 106, no. 5, pp. 808–828, 2018.
  * [5] J. V. D. M. Cardoso, J. Ying, and D. P Palomar,  “Algorithms for learning graphs in financial markets,”  arXiv preprint arXiv:2012.15410, 2020.
  * [6] G. Mateos, S. Segarra, A. G. Marques, and A. Ribeiro,  “Connecting the dots: Identifying network structure via graph signal processing,”  IEEE Signal Process. Mag., vol. 36, no. 3, pp. 16–43, 2019.
  * [7] S. Sardellitti, S. Barbarossa, and P. Di Lorenzo,  “Graph topology inference based on sparsifying transform learning,”  IEEE Trans. Signal Process., vol. 67, no. 7, pp. 1712–1727, 2019\.
  * [8] N. Meinshausen and P. Buhlmann,  “High-dimensional graphs and variable selection with the lasso,”  Ann. Statist., vol. 34, pp. 1436–1462, 2006.
  * [9] J. Friedman, T. Hastie, and R. Tibshirani,  “Sparse inverse covariance estimation with the graphical lasso,”  Biostatistics, vol. 9, no. 3, pp. 432–441, 2008.
  * [10] J. Mei and J.M.F. Moura,  “Signal processing on graphs: Estimating the structure of a graph,”  in IEEE Intl. Conf. Acoust., Speech and Signal Process. (ICASSP), 2015, pp. 5495–5499.
  * [11] H. E. Egilmez, E. Pavez, and A. Ortega,  “Graph learning from data under laplacian and structural constraints,”  IEEE J. Sel. Topics Signal Process., vol. 11, no. 6, pp. 825–841, 2017.
  * [12] S. Segarra, A. G. Marques, G. Mateos, and A. Ribeiro,  “Network topology inference from spectral templates,”  IEEE Trans. Signal Info. Process. Networks, vol. 3, no. 3, pp. 467–483, Sep. 2017.
  * [13] A. Buciulea, S. Rey, and A. G Marques,  “Learning graphs from smooth and graph-stationary signals with hidden variables,”  arXiv preprint arXiv:2111.05588, 2021.
  * [14] G. V. Karanikolas, G. B. Giannakis, K. Slavakis, and R. M. Leahy,  “Multi-kernel based nonlinear models for connectivity identification of brain networks,”  in IEEE Intl. Conf. Acoust., Speech and Signal Process. (ICASSP), Shanghai, China, Mar. 20-25, 2016.
  * [15] A. G. Marques, S. Segarra, G. Leus, and A. Ribeiro,  “Stationary graph processes and spectral estimation,”  IEEE Trans. Signal Process., vol. 65, no. 22, pp. 5911–5926, 2017\.
  * [16] Nathanaël Perraudin and Pierre Vandergheynst,  “Stationary signal processing on graphs,”  IEEE Trans. Signal Process., vol. 65, no. 13, pp. 3462–3477, 2017\.
  * [17] B. Girault,  “Stationary graph signals using an isometric graph translation,”  in European Signal Process. Conf. (EUSIPCO), Aug 2015, pp. 1516–1520.
  * [18] X. Dong, D. Thanou, M. Rabbat, and P. Frossard,  “Learning graphs from data: A signal representation perspective,”  IEEE Signal Process. Mag., vol. 36, no. 3, pp. 44–63, 2019.
  * [19] M. Hong, M. Razaviyayn, Z. Luo, and J. Pang,  “A unified algorithmic framework for block-structured optimization involving big data: With applications in machine learning and signal processing,”  IEEE Signal Process. Mag., vol. 33, no. 1, pp. 57–77, 2015.
  * [20] R. Shafipour and G. Mateos,  “Online topology inference from streaming stationary graph signals with partial connectivity information,”  Algorithms, vol. 13, no. 9, pp. 228, 2020.