\section{Introduction}
\label{sec:introduction}

The objective of stellar interferometric imaging is to recover an
approximation of the specific brightness distribution $I_\lambda(\V{\theta})$
of observed astronomical objects given measurements providing incomplete
samples of the spatial Fourier transform of $I_\lambda(\V{\theta})$.
Reconstruction of a monochromatic image from optical interferometry data is a
challenging task which has been the subject of fruitful research and resulted
in various algorithms (\emph{e.g.}\xspace, Mi\textsc{ra}\xspace \cite{Thiebaut-2008-Marseille}, \textsc{Bsmem}\xspace
\cite{Buscher-1994-BSMEM,Baron_Young-2008-Marseille}, \textsc{Wisard}\xspace
\cite{Meimon_et_al-2005-weak_phase_imaging}, the \emph{building-block method}
\cite{Hofmann_Weigelt-1993-building_blocks}).  When dealing with
multi-spectral data, a first possibility is to process each wavelength
independently and reconstruct a monochromatic image for each subset of
measurements from a given spectral channel.  For instance, this is what have
been done by le~Bouquin et al. \cite{leBouquin_et_al-2009-TLep} for the
multi-spectral images of the Mira star T~Lep.  Another possibility is to
exploit some assumed spectral continuity of $I_\lambda(\theta)$ and process
the multi-spectral data globally to reconstruct an approximation of the 3-D
distribution $I_\lambda(\theta)$.  This computationally more challenging
approach can potentially lead to better reconstructions. Significant
improvements have been shown when following such spatio-spectral processing in
the context of integral field spectral spectroscopy
\cite{Soulez_et_al-2011-Lisbon,
  Bongard_et_al-2011-hyperspectral_deconvolution, Bourguignon2011}.  This
paper describes a method to jointly reconstruct multi-spectral optical
interferometric data.

In order to simplify the problem, we restricted our study to the cases where
the complex visibilities are observed and where the observed scene is a
collection of point-like sources. This correspond, for instance, to the
science case of the instrument \textsc{Gravity}\xspace which will be installed at the Very
Large Telescope Interferometer (VLTI) to carry out astrometry with absolute
phase reference of stars in the galactic center or in globular clusters
\cite{Gillessen_et_al-2010-GRAVITY}. In some sense, this latter assumption
makes our algorithm a successor of the \textsc{Clean}\xspace algorithm
\cite{Hogbom-1974-CLEAN, Schwarz-1978-Clean} and the \emph{building-block
  method} \cite{Hofmann_Weigelt-1993-building_blocks} developed for recovering
monochromatic images from radio and optical interferometric data respectively.
The \textsc{Clean}\xspace algorithm have been proposed for the processing of \textsc{Gravity}\xspace
interferometric data \cite{Vincent_el_al-2011-hotspot} but only considering
``gray'' data and not more than three stars in the field of view.  In addition
to processing multi-variate data, we also introduce the explicit minimization
of a non-differentiable regularization term so as to favor spatial sparsity of
the reconstructed brightness distribution in a way which is known to be more
efficient \cite{Donoho-2006-ell1_approx, Fornasier_Rauhut-2008-joint_sparsity}
than greedy algorithms like \textsc{Clean}\xspace
\cite{Marsh_Richardson-1987-CLEAN_objective_function} or the
\emph{building-block method} \cite{Hofmann_Weigelt-1993-building_blocks}.  The
method presented in this paper improves over early developments presented as
an invited paper at the 2012 SPIE Conf. on Astronomical Telescopes \&
Instrumentation in Amsterdam\cite{Thiebaut_Soulez-2012-Amsterdam}.

Our paper is organized as follows: we first summarize the inverse approach for
image reconstruction from interferometric data and discuss various
possibilities to impose spatial sparsity, we then detail our algorithm for
minimizing the objective function; finally we present some results on
simulated data and discuss the advantages of our approach.

\section{Method}
\label{sec:method}

\subsection{General principle of image reconstruction}

Following an inverse approach, we state image reconstruction as a constrained
optimization problem \cite{Thiebaut-2009-interferometric_imaging}:
\begin{equation}
  \V{x}^{+} = \MathFuncWithLimits{arg\,min}_{\V{x} \in \Set{X}}
  \left\{ \Fcost_\DataTag(\V{x}) + \mu \, \Fcost_\PriorTag(\V{x}) \right\}
  \label{eq:image-problem}
\end{equation}
where $\V{x} \in \Set{R}^{\abs{\V{x}}}$ are the sought image parameters,
$\abs{\V{x}} = \Card(\V{x})$ is the number of parameters, $\Set{X} \subset
\Set{R}^{\abs{\V{x}}}$ is the subset of feasible parameters, $\Fcost_\DataTag(\V{x})$
is a data fitting term which enforces agreement of the model with the
measurements $\V{y} \in \Set{R}^{\abs{\V{y}}}$, $\Fcost_\PriorTag(\V{x})$ is a
regularization term and $\mu > 0$ is a so-called hyper-parameter used to tune
the relative weight of the regularization.  Following a Bayesian
interpretation, $\Fcost_\DataTag$ is the opposite of the log-likelihood of the
measurements, $\Fcost_\PriorTag$ is the opposite of the log of the prior distribution
of the parameters, and $\V{x}^{+}$ defines the maximum a posteriori (MAP)
estimate.

Constraining the solution to belong to the feasible set $\Set{X}$ is a mean
to impose \emph{strict constraints} such as the non-negativity:
\begin{equation}
  \Set{X} = \{\V{x} \in \Set{R}^{n} : \V{x} \ge 0\}
  \label{eq:feasible-set}
\end{equation}
where the inequality $\V{x} \ge 0$ is to be taken element-wise.

\subsection{Direct model and likelihood}

In this study, we assume that the optical interferometric data consist in
Fourier transform of the brightness distribution $\FT{I}_{\lambda}(\V{\nu})$
measured for a finite set of spatial frequencies $\V{\nu} = \V{B}/\lambda$
with $\V{B}$ the interferometric baseline (projected in a plane perpendicular
to the line of sight) and $\lambda$ the wavelength
\cite{Thiebaut_Giovannelli-2010-interferometry}.

The most practical representation of a multi-variate distribution such as
$I_\lambda(\theta)$ by a finite number of parameters $\V{x}$ consists in
sampling $I_\lambda(\theta)$ separately along its spatial and spectral
dimensions.  The image parameters are then:
\begin{equation}
  x_{n,\ell} \approx I_{\lambda_\ell}(\V{\theta}_n)
  \label{eq:multi-spectral-parameters}
\end{equation}
for $\lambda_\ell \in \Set{W}$ the list of sampled wavelengths and
$\V{\theta}_n \in \Set{A}$ the list of angular directions, the so-called
\emph{pixels}.

For the sake of notational simplicity, we use the same wavelengths in $\Set{W}$ as the
ones of the data spectral channels and we denote by $y_{p,m,\ell}$ the real
($p=1$) or imaginary ($p=2$) part of the complex visibility obtained with
$m$th baseline in $\ell$th spectral channel.  This notation is intended to
clarify the equations and does not impose or assume that all baselines have
been observed in all spectral channels.  By considering that complex numbers
are just pairs of real values, our notation also avoids dealing with complex
arithmetic.  In our framework, the model of the data is affine:
\begin{align}
  y_{p,m,\ell}
  &= (\M{H}\cdot\V{x})_{p,m,\ell} + e_{p,m,\ell} \notag \\
  &= \sum\nolimits_{n} H_{p,m,n,\ell} \, x_{n,\ell} + e_{p,m,\ell} \, ,
  \label{eq:direct-model}
\end{align}
where the term $\V{e}$ accounts for noise and modeling approximations.
Formally, the coefficients of the operator $\M{H}$ are given by
\cite{Thiebaut_Giovannelli-2010-interferometry}:
\begin{equation}
  H_{p,m,n,\ell} =
  \left\{\begin{array}{ll}
  {\displaystyle + \cos(\V{\theta}_n\T\!\!\cdot\V{B}_{m}/\lambda_\ell)}
  & \text{for $p=1$}\\[1ex]
  {\displaystyle - \sin(\V{\theta}_n\T\!\!\cdot\V{B}_{m}/\lambda_\ell)}
  & \text{for $p=2$}\\
  \end{array}\right.
  \label{eq:multi-spectral-model-matrix}
\end{equation}
with $\V{B}_m$ the $m$th observed baseline and $\V{\theta}_n\T\!\!\cdot\V{B}_{m}$
the usual scalar product between $\V{B}_{m}$ and $\V{\theta}_n$.

At least because of the strict constraints imposed by the feasible set
$\Set{X}$, solving the image reconstruction problem in \Eq{eq:image-problem}
must be carried out by an iterative algorithm.  Owing to the size of the
problem, a fast version of $\M{H}$ has to be implemented.  First, we note that
the model is separable along the spectral dimension (using a conventional
matrix representation, $\M{H}$ would have a block diagonal structure):
\begin{equation}
  \V{y}_{\ell} = \M{H}_{\ell} \cdot \V{x}_{\ell} + \V{e}_{\ell} \, ,
  \label{eq:direct-model-in-a-channel}
\end{equation}
where the index $\ell$ denotes the sub-vector or the sub-operator restricted
to the coefficients corresponding to the $\ell$th spectral channel.  With the
generalization of multi-processor computers or multi-core processors, this
property of the operator $\M{H}$ may be easily exploited to parallelize the
code to apply $\M{H}$ (or its adjoint $\M{H}\T$) to a given argument.  Second,
an algorithm such as the \emph{nonuniform fast Fourier transform} (NU-FFT)
\cite{Fessler_Sutton-2003-NUFFT} can be implemented to speed up the
computations by approximating the operator $\M{H}_{\ell}$ by:
\begin{equation}
  \M{H}_{\ell} \approx \M{R}_{\ell} \cdot \M{F} \cdot \M{S}
  \label{eq:multi-spectral-nufft}
\end{equation}
where $\M{F}$ is the discrete Fourier transform (DFT), $\M{R}_{\ell}$
interpolates the discrete spatial frequencies resulting from the DFT at the
frequencies observed in $\ell$th channel and $\M{S}$ is a zero-padding and
apodizing operator.  Zero-padding improves the accuracy of the approximation,
while apodization pre-compensates for the convolution by the interpolation
kernel used in $\M{R}_{\ell}$ \cite{Fessler_Sutton-2003-NUFFT}.  Note that
only the interpolation in the Fourier domain $\M{R}_{\ell}$ depends on the
spectral channel.  In NU-FFT, $\M{S}$ is diagonal, $\M{R}_{\ell}$ is very
sparse and $\M{F}$ is implemented by a fast Fourier transform (FFT) algorithm,
thus the approximation in \Eq{eq:multi-spectral-nufft} is very fast to
compute.

Assuming Gaussian noise distribution, the likelihood term writes:
\begin{equation}
  \Fcost_\DataTag(\V{x}) = \frac{1}{2} \, (\M{H} \cdot \V{x} - \V{y})\T
  \cdot \M{W} \cdot (\M{H} \cdot \V{x} - \V{y})
  \label{eq:fdata}
\end{equation}
where $\M{W} \in \Set{R}^{\abs{\V{y}} \times \abs{\V{y}}}$ is a statistical
weighting matrix; in principle, $\M{W}$ is the inverse of the covariance
matrix of the measurements: $\M{W} = \Cov(\V{y})^{-1}$.

\subsection{Regularization based on spatial sparsity}

Due to the voids in the spatial frequencies covered by the observations, the
constraints provided by the data alone do not suffice to define a unique
image.  The prior constraints imposed by $\Fcost_\PriorTag(\V{x})$ are then required to
help choosing a unique solution among all the images that are compatible with
the measurements.

In this paper, we focus on a particular type of astronomical targets which
consist in a number of point-like sources with different spectral energy
distributions.  This includes the case of multiple stars, globular clusters,
or groups of stars as observed in the center of our galaxy.  For such objects,
the most effective means to regularize the problem is to favor spatially
sparse distributions, \emph{i.e.}\xspace images with as few sources as possible on a dark
background.  In this section, we derive expressions of the regularization term
$\Fcost_\PriorTag(\V{x})$ suitable to favor spatially sparse distributions.

\subsubsection{Fully separable sparsity}

It is now well known that using the $\ell_1$ norm as the regularization term
is an effective mean to impose the sparsity of the solution while
approximating the data \cite{Donoho-2006-ell1_approx}.  This leads to take
$\Fcost_\PriorTag(\V{x}) = \Fcost_\SparseTag(\V{x})$ with:
\begin{equation}
  \Fcost_\SparseTag(\V{x}) = \norm{\V{x}}_1 \stackrel{\text{\tiny def}}{=} \sum_{k,\ell} \Abs{x_{k,\ell}}
  = \sgn(\V{x})\T\!\!\cdot\V{x} \, ,
  \label{eq:sparse-prior}
\end{equation}
where $\sgn(\V{x})$ is the sign function applied element-wise to the
parameters $\V{x}$.  When the parameters are non-negative, $\sgn(\V{x}) =
\V{1}$ with $\V{1} \stackrel{\text{\tiny def}}{=} (1,\ldots,1)\T$.

The regularization term $\Fcost_\SparseTag(\V{x})$ in \Eq{eq:sparse-prior} is
completely separable.  In our framework where the model is spectrally
separable, the global criterion defined in \Eq{eq:image-problem} is therefore
separable along the spectral dimension.  Provided data from different spectral
channels are statistically independent, the image reconstruction can be solved
independently for each spectral channel.

If the wavelength samples $\lambda_\ell$ of the discrete model $\V{x}$ of
$I_\lambda(\theta)$ do not coincide with the effective wavelengths of the
data, spectral interpolation of the model is required to match the observed
wavelengths.  In this case, a certain spectral correlation is intrinsic to the
model and the 3-D image reconstruction has to be performed globally even if
the regularization does not impose any kind of spectral continuity.

\subsubsection{Non-separable spatial-only sparsity}

Physically, sources emit light at all wavelengths and we expect better image
reconstruction if we can favor restored sources having the same position
whatever the wavelength.  Clearly, this is not achieved by the regularization
$\Fcost_\SparseTag(\V{x})$ in \Eq{eq:sparse-prior} which is fully separable.  In
order to impose some spectral continuity while favoring spatial sparsity, we
consider the following regularization instead
\cite{Fornasier_Rauhut-2008-joint_sparsity, Soulez_et_al-2011-Lisbon}:
\begin{equation}
  \Fcost_\JointTag(\V{x}) = \sum\nolimits_{n}
  \left(\sum\nolimits_{\ell} x_{n,\ell}^2\right)^{\!1/2}
  \label{eq:joint-sparse-prior}
\end{equation}
with $n$ the spatial index (pixel) and $\ell$ the spectral channel.  The fact
that such a regularization favors spatial sparsity and spectral grouping is a
consequence of the triangular inequality
\cite{Fornasier_Rauhut-2008-joint_sparsity}.  This sparsity prior is a special
case of several recent generalizations such as group \textsc{Lasso}\xspace
\cite{yuan2006model}, mixed norms \cite{Kowalski2009mixednorms} or structured
sparsity \cite{jenatton2011structured}.

\subsubsection{Explicit spectral continuity and gray object}

The penalization defined in \Eq{eq:joint-sparse-prior} can be seen as a
spatial regularization term which favors spectral grouping but no spectral
continuity nor spectral smoothness.  Some authors
\cite{Thiebaut_Mugnier-2006-nulling, Soulez_et_al-2011-Lisbon,
  Bongard_et_al-2011-hyperspectral_deconvolution} have shown the efficiency of
exploiting the spectral continuity of the sought distribution by using, in
addition to a spatial regularization term, an additional spectral
regularization like:
\begin{equation}
  \Fcost_\SpectralTag(\V{x}) = \sum_n \mu_n
  \sum_\ell (x_{n,\ell+1} - x_{n,\ell})^2
  \label{eq:spectral-regul}
\end{equation}
with $\mu_n > 0$ suitable regularization weights.  In the limit $\mu_n
\rightarrow \infty, \forall n$, the regularization in \Eq{eq:spectral-regul}
amounts to assuming that the spectral energy distributions of all sources are
flat, that is:
\begin{equation}
  x_{n,\ell} = g_{n} \,, \quad \forall (n, \ell)
  \label{eq:gray-image}
\end{equation}
where $\V{g}$ is a \emph{gray image} of the object which does not depend on
the spectral index $\ell$.  To speed up the reconstruction, only the gray
image has to be reconstructed, using the model:
\begin{equation}
  y_{p,m,\ell} = \sum_{n} H_{p,m,n,\ell} \, g_{n} + e_{p,m,\ell} \, ,
  \label{eq:gray-model}
\end{equation}
and, to impose the spatial sparsity, the regularization on $\V{g}$ is:
\begin{equation}
  \label{eq:gray-regularization}
  \Fcost_\SparseTag(\V{g}) = \norm{\V{g}}_1 = \sum\nolimits_n \abs{g_n}
  = \sgn(\V{g})\T\!\!\cdot\V{g} \, .
\end{equation}

\subsection{Optimization algorithm}

Most existing image reconstruction algorithms for optical interferometry (\emph{e.g.}\xspace,
Mi\textsc{ra}\xspace \cite{Thiebaut-2008-Marseille}, \textsc{Bsmem}\xspace \cite{Buscher-1994-BSMEM,
  Baron_Young-2008-Marseille} and \textsc{Wisard}\xspace
\cite{Meimon_et_al-2005-weak_phase_imaging}) were designed for minimizing a
smooth cost function.  For that purpose, non-linear conjugate gradient method
\cite{Nocedal_Wright-2006-numerical_optimization} or limited memory
quasi-Newton methods such as VMLM-B \cite{Thiebaut-2002-optim_bdec} are quite
efficient and easy to use as they only require computing the cost function and
its gradient.  A notable exception is
\textsc{Macim}\cite{Ireland_et_al-2006-MACIM} which is based on a
Markov-Chain-Monte-Carlo (MCMC) optimization strategy suitable, in theory, for
any type of criteria, in particular the non-smooth and non-convex ones; in
practice, this is however too computationally intensive for estimating a large
number of parameters as it is the case for image reconstruction. When using
non smooth regularizations as the ones in \Eq{eq:sparse-prior} and
\Eq{eq:joint-sparse-prior} to impose spatial sparsity, optimization algorithms
based on Newton method (that is, on a quadratic approximation of the cost
function) are inefficient and completely different optimization strategies
must be followed to solve the problem in \Eq{eq:image-problem} with
non-differentiable cost functions.  In our algorithm, we introduce variables
splitting \cite{Combettes_Pesquet-2011-proximal_methods} to handle the two
terms of the cost function as independently as possible and we implement an
\emph{alternating direction method of multipliers}
\cite{Boyd_et_al-2010-method_of_multipliers} (ADMM) to solve the resulting
constrained problem.  The augmented Lagrangian with ADMM emerges as the most
effective in the family of decomposition methods that includes proximal
methods \cite{Combettes_Pesquet-2011-proximal_methods}, variable splitting
with quadratic penalty \cite{wang2008new}, iterative Bregman
\cite{yin2008bregman}.  See \cite{Afonso_et_al-2011-C_SALSA} for detailed
comparisons.

\subsubsection{Variable Splitting and ADMM}
\label{sec:ADMM}

Introducing auxiliary variables $\V{z}$, minimization of the two-term cost
function in \Eq{eq:image-problem} can be recast in the equivalent constrained
problem:
\begin{equation}
  \min_{\V{x} \in \Set{X}, \V{z}} \left\{ \Fcost_\DataTag(\V{z})
  + \mu\,\Fcost_\PriorTag(\V{x}) \right\}
  \quad \text{s.t.}\quad \V{x} = \V{z} \, .
  \label{eq:split-problem}
\end{equation}
Imposing that $\V{x} \in \Set{X}$ (\emph{i.e.}\xspace $\V{x} \ge \V{0}$) rather than $\V{z}
\in \Set{X}$ is not arbitrary and our motivation for that choice is explained
in what follows.  Another possible splitting would have been to choose the
auxiliary variables as $\V{z} = \M{H}\cdot\V{x}$ but this would have prevented
us to exploit the separability of the resulting penalty with respect to
$\V{x}$.

The augmented Lagrangian \cite{Nocedal_Wright-2006-numerical_optimization} is
a very useful method to deal with constrained problems such as the one in
\Eq{eq:split-problem}.  In our case, the augmented Lagrangian writes:
\begin{align}
  \mathcal{L}_{\rho}(\V{x},\V{z},\V{u}) &= \Fcost_\DataTag(\V{z})
  + \mu\,\Fcost_\PriorTag(\V{x}) \notag \\
  &\quad + \V{u}\T\!\!\cdot(\V{x} - \V{z})
  + \frac{\rho}{2} \, \norm{\V{x} - \V{z}}_2^2 \, ,
  \label{e:augmented-Lagrangian}
\end{align}
with $\V{u}$ the Lagrange multipliers associated with the constraint $\V{x} =
\V{z}$, $\rho > 0$ the quadratic weight of the constraints, and
$\norm{\V{v}}_2$ the Euclidean ($\ell_2$) norm of $\V{v}$.  Note that taking
$\rho=0$ yields the classical Lagrangian of the constrained problem.

The alternating direction method of multipliers
\cite{Boyd_et_al-2010-method_of_multipliers} (ADMM) consists in alternatively
minimizing the augmented Lagrangian for $\V{x}$ given $\V{z}$ and $\V{u}$,
then for $\V{z}$ given $\V{x}$ and $\V{u}$, and finally updating the
multipliers $\V{u}$.  This scheme, adapted to our specific problem in
\Eq{eq:split-problem}, is detailed by the following algorithm with the
convention that $\V{v}^{(t)}$ is the value of $\V{v}$ at iteration number $t$:

\smallskip

\textbf{Algorithm 1.}  \emph{Resolution of problem~(\ref{eq:split-problem}) by
  alternating direction method of multipliers}.  Choose initial variables
$\V{z}^{(0)}$ and Lagrange multipliers $\V{u}^{(0)}$.  Then repeat, for $t =
1, 2, \ldots$ until convergence, the following steps:
\begin{align}
  \intertext{1.~choose $\rho^{(t)} > 0$ and update variables $\V{x}$:}
  \V{x}^{(t)}
  &= \MathFuncWithLimits{arg\,min}_{\V{x} \in \Set{X}}
  \mathcal{L}_{\rho^{(t)}}\bigl(\V{x}, \V{z}^{(t-1)}, \V{u}^{(t-1)}\bigr) \notag \\
  &= \MathFuncWithLimits{arg\,min}_{\V{x} \in \Set{X}}
  \Bigl\{ \Fcost_\PriorTag(\V{x})
  + \frac{\rho^{(t)}}{2\,\mu} \, \bigNorm{ \V{x} - \proxy{\V{x}}^{(t)} }_2^2 \Bigr\}
  \label{eq:x-update}
  \intertext{\hspace*{2.5ex}with:}
  \proxy{\V{x}}^{(t)} &= \V{z}^{(t-1)} - \V{u}^{(t-1)}/\rho^{(t)} \, ;
  \label{eq:x-proxy-def}
  \intertext{2.~update variables $\V{z}$:}
  \V{z}^{(t)}
  &= \MathFuncWithLimits{arg\,min}_{\V{z}}
  \mathcal{L}_{\rho^{(t)}}\bigl(\V{x}^{(t)}, \V{z}, \V{u}^{(t-1)}\bigr) \notag \\
  &= \MathFuncWithLimits{arg\,min}_{\V{z}}
  \Bigl\{ \Fcost_\DataTag(\V{z})
  + \frac{\rho^{(t)}}{2} \, \bigNorm{\V{z} - \proxy{\V{z}}^{(t)}}_2^2 \Bigr\}
  \label{eq:z-update}
  \intertext{\hspace*{2.5ex}with:}
  \proxy{\V{z}}^{(t)} &= \V{x}^{(t)} + \V{u}^{(t-1)}/\rho^{(t)} \, ;
  \label{eq:z-proxy-def}
  \intertext{3.~update multipliers $\V{u}$:}
  \V{u}^{(t)} &= \V{u}^{(t-1)} + \rho^{(t)} \,
  \left(\V{x}^{(t)} - \V{z}^{(t)}\right) \, .
  \quad\ensuremath{Box{}}
  \label{eq:u-update}
\end{align}

Our algorithm can be seen as an instance of SALSA
\cite{Afonso_et_al-2010-SALSA} with however some improvements.  First, we deal
with the additional constraints that the variables are non-negative. Second,
we allow for changing the weight of the augmented penalty at every iteration
which can considerably speed-up convergence. Third, for real observations, the
operator $\M{H}$ cannot be easily diagonalized (\emph{e.g.}\xspace by using FFT) thus the
updating of variables $\V{z}$ cannot be exactly carried out. Finally, we
consider the possibility of warm starting the algorithm with a solution
previously computed.  This latter feature is of interest to improve a solution
if too few iterations have been performed or to find the solution of the
problem with a slightly different value of the hyper-parameter $\mu$.

In Appendix~\ref{sec:proximity-operators}, we show how the updating of
the variables $\V{x}$ (step~1 of Algorithm~1) can be implemented taking into
account the constraints that the parameters are non-negative.  This is our
motivation for imposing $\V{x}\in\Set{X}$ on the variables $\V{x}$ and not on
the variables $\V{z}$.  This avoids introducing additional auxiliary variables
for the sole purpose of accounting for the feasible set.  The formulae to
update the variables $\V{x}\in\Set{X}$ for $\Fcost_\SparseTag(\V{x})$ and
$\Fcost_\JointTag(\V{x})$ are respectively given by \Eq{eq:sparse-prox-plus} and
\Eq{eq:joint-sparse-prox-plus} in Appendix~\ref{sec:proximity-operators}.

\subsubsection{Solving for the auxiliary variables}
\label{sec:solving-auxiliary-variables}

Since $\mathcal{L}_{\rho}(\V{x},\V{z},\V{u})$ is quadratic with respect to
$\V{z}$, updating these variables (step~2 of Algorithm~1) amounts to solving
the linear problem:
\begin{equation}
  \M{A}^{(t)}\cdot\V{z}^{(t)} = \V{b}^{(t)}
  \label{eq:z-update-linear-problem}
\end{equation}
with:
\begin{align}
  \M{A}^{(t)} &= \M{H}\T\!\!\cdot\M{W}\cdot\M{H} + \rho^{(t)}\,\M{I} \\
  \V{b}^{(t)} &= \M{H}\T\!\!\cdot\M{W}\cdot\V{y} + \rho^{(t)}\,\V{x}^{(t)} + \V{u}^{(t-1)}
\end{align}
with $\M{I}$ the identity matrix (of suitable size).  Since the augmented term
is diagonal and provided data from different spectral channels are
statistically independent, this problem can be solved separately for each
spectral channel:
\begin{equation}
  \M{A}^{(t)}_\ell\cdot\V{z}^{(t)}_\ell = \V{b}_\ell^{(t)} \, ,\quad \forall
  \ell\, ,
\end{equation}
with:
\begin{align}
  \M{A}^{(t)}_\ell &= \M{H}_\ell\T\!\!\cdot\M{W}_\ell\cdot\M{H}_\ell + \rho^{(t)}\,\M{I} \\
  \V{b}^{(t)}_\ell &= \M{H}_\ell\T\!\!\cdot\M{W}_\ell\cdot\V{y}_\ell
  + \rho^{(t)}\,\V{x}_\ell^{(t)} + \V{u}_\ell^{(t-1)}
\end{align}
using the same conventions as in \Eq{eq:direct-model-in-a-channel}.

In practice, we (approximately) solve these problems by means of the conjugate
gradients algorithm \cite{Hestenes_Stiefel-1952-conjugate_gradients,
  Nocedal_Wright-2006-numerical_optimization} and starting with the previous
solution $\V{z}^{(t-1)}$.  Since, in $\M{A}^{(t)}$, the Hessian matrix
$\M{H}\T\!\!\cdot\M{W}\cdot\M{H}$ is regularized by the term
$\rho^{(t)}\,\M{I}$, its condition number is better than that of
$\M{H}\T\!\!\cdot\M{W}\cdot\M{H}$.  We therefore expect that the conjugate
gradients algorithm has a better convergence rate with $\M{A}^{(t)}$ than with
$\M{H}\T\!\!\cdot\M{W}\cdot\M{H}$.  Moreover, according to Eckstein-Bertsekas
theorem \cite{Eckstein_Bertsekas-1992-Douglas-Rachford_splitting}, the ADMM
algorithm is proved to converge provided that approximations in the update of
auxiliary variables $\V{z}$ be absolutely summable.  That is:
\begin{equation}
  \sum_{t = 1}^{\infty} \bigNorm{\V{z}^{(t)} - \V{z}^{(t)}_\Tag{exact}}_2 < \infty
  \label{eq:Eckstein-Bertsekas-constraint}
\end{equation}
must hold with $\V{z}^{(t)}_\Tag{exact}$ the exact solution of
\Eq{eq:z-update} and $\V{z}^{(t)}$ the approximate solution of
\Eq{eq:z-update} returned by the conjugate gradient iterations.  The
demonstration in \cite{Eckstein_Bertsekas-1992-Douglas-Rachford_splitting}
consider ADMM iterations with a fixed quadratic penalty parameter $\rho$, so
it may not strictly apply to our method where $\rho$ is allowed to vary (see
Section~\ref{sec:tuning-rho}).  An easy solution to warrant convergence is to
fix the value of $\rho$ after a certain number of ADMM iterations.
Nevertheless, we observed in our tests that the value of $\rho$ stabilizes to
a fixed value without imposing this.  In
Appendix~\ref{sec:Eckstein-Bertsekas-constraint}, we show how to set the
stopping criterion of the conjugate gradient method so that the constraint in
\Eq{eq:Eckstein-Bertsekas-constraint} holds.  In our tests, we however simply
stop the conjugate gradient iterations when the Euclidean norm of the
residuals of \Eq{eq:z-update-linear-problem} becomes significantly smaller
than its initial value:
\begin{equation}
  \norm{\M{A}^{(t)}\cdot\V{z}^{(t)} - \V{b}^{(t)}}_2 \le \epsilon_\Tag{CG} \,
  \norm{\M{A}^{(t)}\cdot\V{z}^{(t-1)} - \V{b}^{(t)}}_2
\end{equation}
with $\epsilon_\Tag{CG} \in (0,1)$.  For our tests, we took
$\epsilon_\Tag{CG}=10^{-2}$ and allowed a maximum of 5 conjugate gradient
iterations.  With this simple prescription, we did not experiment any
divergence of the global algorithm although it may depend on the problem at
hands and could require to be adapted.

\subsubsection{Stopping criteria}
\label{sec:convergence}

At the solution $\{\V{x}^{*}, \V{z}^{*}, \V{u}^{*}\}$ of problem
(\ref{eq:split-problem}), Karush-Kuhn-Tucker (KKT) conditions of optimality
\cite{Nocedal_Wright-2006-numerical_optimization} stipulate that, the
constraints must be satisfied and that the solution must be a stationary point
of the Lagrangian:
\begin{align}
  &\V{x}^{*} = \V{z}^{*} \label{eq:KKT-1} \\
  &\V{0} \in \partial_{\V{x}} \mathcal{L}_0(\V{x}^{*}, \V{z}^{*}, \V{u}^{*})
  = \mu\,\partial \Fcost_\PriorTag(\V{x}^{*}) + \V{u}^{*}
  \label{eq:KKT-2} \\
  &\V{0} \in \partial_{\V{z}} \mathcal{L}_0(\V{x}^{*}, \V{z}^{*}, \V{u}^{*})
  = \partial \Fcost_\DataTag(\V{z}^{*}) - \V{u}^{*}
   \label{eq:KKT-3}
\end{align}
where $\partial$ denotes the subdifferential operator
\cite{Boyd_et_al-2010-method_of_multipliers}. Since $\Fcost_\DataTag$ is
differentiable, $\in$ and $\partial \Fcost_\DataTag$ can be replaced by $=$ and
by $\nabla \Fcost_\DataTag$, the gradient of $\Fcost_\DataTag$ in the third KKT
condition (\ref{eq:KKT-3}) which becomes:
\begin{equation}
  \V{u}^{*} = \nabla \Fcost_\DataTag(\V{z}^{*}) \, .

  \label{eq:KKT-3-alt}
\end{equation}
If $\V{z}^{(t)}$ exactly minimizes $\mathcal{L}_{\rho^{(t)}}(\V{x}^{(t)},
\V{z}, \V{u}^{(t-1)})$, we have:
\begin{align}
  &\nabla \Fcost_\DataTag(\V{z}^{(t)}) - \V{u}^{(t-1)}
  + \rho^{(t)} \, (\V{z}^{(t)} - \V{x}^{(t)}) = \V{0} \notag \\
  \Longrightarrow\
  &\nabla \Fcost_\DataTag(\V{z}^{(t)}) = \V{u}^{(t)} \, ,
  \label{eq:KKT-3-admm}
\end{align}
thus the 3rd KKT condition in \Eq{eq:KKT-3-alt} is automatically satisfied at
the end of an \emph{exact} ADMM iteration.

As we use the method of conjugate gradients to solve
\Eq{eq:z-update-linear-problem}, \Eq{eq:KKT-3-admm} is only approximately
satisfied.  Moreover, updating the multipliers $\V{u}$ according to step~3 of
the algorithm may be subject to accumulation of rounding errors.  The
stability of the algorithm or its convergence rate may be improved by taking
$\V{u}^{(t)} = \nabla \Fcost_\DataTag(\V{z}^{(t)})$.  In our tests, tough, we have not seen
significant differences between updating the Lagrange multipliers according to
\Eq{eq:u-update} or according to \Eq{eq:KKT-3-admm}.

Since $\V{x}^{(t)}$ minimizes $\mathcal{L}_{\rho^{(t)}}(\V{x},
\V{z}^{(t-1)}, \V{u}^{(t-1)})$, we have:
\begin{align}
  \V{0}
  &\in \mu \, \partial \Fcost_\PriorTag(\V{x}^{(t)}) + \V{u}^{(t-1)}
  + \rho^{(t)} \, (\V{x}^{(t)} - \V{z}^{(t-1)}) \notag \\
  &\in \mu \, \partial \Fcost_\PriorTag(\V{x}^{(t)}) + \V{u}^{(t)}
  + \rho^{(t)} \, (\V{z}^{(t)} - \V{z}^{(t-1)}) \notag \\
  \Longrightarrow&\quad
  - \rho^{(t)} \, (\V{z}^{(t)} - \V{z}^{(t-1)})
  \in \mu \, \partial \Fcost_\PriorTag(\V{x}^{(t)}) + \V{u}^{(t)}
  \notag
\end{align}
thus:
\begin{equation}
  \V{s}^{(t)} = \rho^{(t)} \, \bigl(\V{z}^{(t)} - \V{z}^{(t - 1)}\bigr)
  \label{eq:dual-residual}
\end{equation}
can be seen as the residuals for the 2nd KKT condition in \Eq{eq:KKT-2},
while:
\begin{equation}
  \V{r}^{(t)} = \V{x}^{(t)} - \V{z}^{(t)}
  \label{eq:primal-residual}
\end{equation}
are the residuals for the primary constraint in \Eq{eq:KKT-1}.

Finally, the KKT conditions imply that the so-called \emph{primal} and
\emph{dual residuals} \cite{Boyd_et_al-2010-method_of_multipliers} defined in
\Eq{eq:primal-residual} and \Eq{eq:dual-residual} must converge to zero.
Following \cite{Boyd_et_al-2010-method_of_multipliers}, we therefore stop the
algorithm when:
\begin{equation}
  \bigl\Vert \V{r}^{(t)} \bigr\Vert_2 \le \tau_\Tag{prim}^{(t)}
  \quad\text{and}\quad
  \bigl\Vert \V{s}^{(t)} \bigr\Vert_2 \le \tau_\Tag{dual}^{(t)} \, ,
  \label{eq:convergence-criterion}
\end{equation}
where the convergence thresholds are given by:
\begin{align}
  \tau_\Tag{prim}^{(t)} &\stackrel{\text{\tiny def}}{=} \sqrt{N}\,\epsilon_\Tag{abs} + \epsilon_\Tag{rel}
  \, \max\bigl(\bigl\Vert \V{x}^{(t)} \bigr\Vert_2,
               \bigl\Vert \V{z}^{(t)} \bigr\Vert_2 \bigr) \, ,
  \label{eq:primal-threshold} \\
  \tau_\Tag{dual}^{(t)} &\stackrel{\text{\tiny def}}{=} \sqrt{N}\,\epsilon_\Tag{abs} +
  \epsilon_\Tag{rel} \, \bigl\Vert \V{u}^{(t)} \bigr\Vert_2 \, ,
  \label{eq:dual-threshold}
\end{align}
where $N = \Card(\V{x})$ is the number of sought parameters,
$\epsilon_\Tag{abs} \ge 0$ and $\epsilon_\Tag{rel} \in (0,1)$ are absolute and
relative convergence tolerances.  For our tests, we found that
$\epsilon_\Tag{abs} = 0$ and $\epsilon_\Tag{rel} = 10^{-3}$ yield sufficient
precision for the solution.

\subsubsection{Initialization and warm start}
\label{sec:warm-start}

Given an initial estimate $\V{z}^{(0)}$ for the auxiliary variables, the 3rd
KKT condition in \Eq{eq:KKT-3-alt} suggests to start the iterative algorithm
with initial Lagrange multipliers given by:
\begin{equation}
  \V{u}^{(0)} = \nabla \Fcost_\DataTag\bigl(\V{z}^{(0)}\bigr)
  = \M{H}\T\!\!\cdot\M{W}\cdot\bigl(\M{H}\T\!\!\cdot\V{z}^{(0)} - \V{y}\bigr) \, .
\end{equation}
Starting the algorithm as suggested has several advantages.  First, there are
no needs for the initial variables $\V{z}^{(0)}$ to belong to the feasible
set.  Second, this readily provides initial Lagrange multipliers.  Note that
starting with an initial estimate $\V{x}^{(0)}$ for the variables $\V{x}$
would have the double drawback that $\V{x}^{(0)}$ must be feasible and that
since $\Fcost_\PriorTag(\V{x})$ may be non-differentiable it does not yields an
explicit expression for the initial Lagrange multipliers.

The other required initial setting is the value of the augmented penalty
parameter $\rho^{(1)}$ used to compute the first estimate $\V{x}^{(1)}$ of the
variables $\V{x}$ given $\V{z}^{(0)}$ and $\V{u}^{(0)}$.  For the subsequent
iterations, $\rho$ can be kept constant or updated according to the
prescription described in Sec.~\ref{sec:tuning-rho}.

To continue the iterations or compute a solution with slightly different
parameters (\emph{e.g.}\xspace the regularization parameter $\mu$), the possibility to
restart the algorithm with the output of a previous run with no loss of
performances regarding the rate of convergence is a needed feature.  This is
called \emph{warm restart} and is simply achieved by saving a minimal set of
variables upon return of the algorithm.  Since each iteration of our algorithm
starts by computing the variables $\V{x}$ given the auxiliary variables
$\V{z}$, the Lagrange multipliers $\V{u}$ and the augmented penalty parameter
$\rho$, it is sufficient to save $\{\V{z},\V{u},\rho\}$ for being able to warm
restart the method.

\subsubsection{Tuning the augmented penalty parameter $\rho$}
\label{sec:tuning-rho}

One of the important settings of the ADMM method is the value of the augmented
penalty parameter $\rho$: if it is too small, the primal constraints $\V{x} =
\V{z}$ will converge too slowly; while the cost functions will decrease too
slowly if $\rho$ is too large.  Some authors, \emph{e.g.}\xspace
\cite{Afonso_et_al-2010-SALSA}, use a constant augmented penalty parameter for
all the iterations which requires trials and errors to find an efficient value
for $\rho$.  In fact, it is worth using a good value of $\rho$ at every
iteration of ADMM to accelerate the convergence
\cite{Boyd_et_al-2010-method_of_multipliers}.  In this section, we describe
means to automatically derive a suitable value for the augmented penalty
parameter following a simple reasoning.

The convergence criterion defined in \Eq{eq:convergence-criterion} is
equivalent to have:
\begin{equation}
  \phi^{(t)} \le 1
  \quad\text{with}\quad
  \phi^{(t)} \stackrel{\text{\tiny def}}{=} \max\left(
  \frac{\bigl\Vert \V{r}^{(t)} \bigr\Vert_2}{\tau_\Tag{prim}^{(t)}},
  \frac{\bigl\Vert \V{s}^{(t)} \bigr\Vert_2}{\tau_\Tag{dual}^{(t)}}
  \right) .
  \label{eq:convergence-criterion-alt}
\end{equation}
According to the updating rules in one ADMM iteration (see Algorithm~1),
$\phi^{(t)}$ does only depend on $\V{z}^{(t-1)}$, $\V{u}^{(t-1)}$ and
$\rho^{(t)}$.  The augmented penalty parameter $\rho^{(t)}$ is therefore the
only tunable parameter that has an incidence on the value of $\phi^{(t)}$ for
the $t$\textsup{th} ADMM iteration. The idea is then to chose the value of
$\rho^{(t)}$ so as to approximately minimize $\phi^{(t)}$.  In terms of number
of ADMM iterations, we expect to achieve the faster convergence of the
algorithm in that way.  However, tuning $\rho$ at every ADMM iteration
requires to repeat each iteration for different values of $\rho$ and has
therefore the same computational cost as several ADMM iterations.  A
compromise has to be found between the accuracy on $\rho$ and the number of
trials.

Our objective is to derive an economical way to find:
\begin{equation}
  \rho^{(t)} \approx \rho^{(t)}_{*} \stackrel{\text{\tiny def}}{=} \MathFuncWithLimits{arg\,min}_{\rho} \phi_{t}(\rho) \, .
  \label{eq:best-rho}
\end{equation}
where $\phi_{t}(\rho)$ is the value taken by $\phi^{(t)}$ when $\rho^{(t)} =
\rho$.  All the quantities, $\norm{\V{r}^{(t)}}_2$, $\norm{\V{s}^{(t)}}_2$,
$\tau_\Tag{prim}^{(t)}$, and $\tau_\Tag{dual}^{(t)}$, involved in $\phi^{(t)}$
vary continuously (though not necessarily smoothly) with respect to
$\rho^{(t)}$; hence, considering the definition of $\rho^{(t)}_{*}$ and
$\phi^{(t)}$, we obtain the following implication:
\begin{equation}
  \rho^{(t)} = \rho^{(t)}_{*}
  \Longrightarrow
  \frac{\bigl\Vert \V{r}^{(t)} \bigr\Vert_2}{\tau_\Tag{prim}^{(t)}} =
  \frac{\bigl\Vert \V{s}^{(t)} \bigr\Vert_2}{\tau_\Tag{dual}^{(t)}} \, .
  \label{eq:best-rho-necessary-condition}
\end{equation}
Besides, the norm of the primal residuals $\norm{\V{r}^{(t)}}_2$ is a
decreasing function of $\rho^{(t)}$; while the norm of the dual residuals
$\norm{\V{s}^{(t)}}_2$ is an increasing function of $\rho^{(t)}$
\cite{Boyd_et_al-2010-method_of_multipliers} and, close enough to the
solution, the values of $\tau_\Tag{prim}$ and $\tau_\Tag{dual}$ should
converge to their final values and thus not depend too much on $\rho$.  Under
these assumptions, the ratios $\norm{\V{r}^{(t)}}_2/\tau_\Tag{prim}^{(t)}$ and
$\norm{\V{s}^{(t)}}_2/\tau_\Tag{dual}^{(t)}$ should also be decreasing and
increasing functions of $\rho$ respectively.  Close to the solution
$\{\V{x}^{*}, \V{z}^{*}, \V{u}^{*}\}$ of the problem, the necessary condition
in \Eq{eq:best-rho-necessary-condition} is therefore also a sufficient
condition to define the optimal value $\rho^{(t)}_{*}$.  These considerations
lead us to choose $\rho^{(t)}$ such that:
\begin{equation}
  \eta^{(t)} \approx 1
  \quad\text{with}\quad
  \eta^{(t)} \stackrel{\text{\tiny def}}{=}
  \frac{\bigl\Vert \V{r}^{(t)} \bigr\Vert_2 \, \tau_\Tag{dual}^{(t)}}
       {\bigl\Vert \V{s}^{(t)} \bigr\Vert_2 \, \tau_\Tag{prim}^{(t)}} \, ,
  \label{eq:eta}
\end{equation}
which is expected to be a decreasing function of $\rho^{(t)}$ close to the
solution.  A better alternative may be to choose $\rho^{(t)}$ such that:
\begin{equation}
  \eta^{(t)}_\Tag{alt} \approx 1
  \quad\text{with}\quad
  \eta^{(t)}_\Tag{alt} \stackrel{\text{\tiny def}}{=}
  \frac{\bigl\Vert \V{r}^{(t)} \bigr\Vert_2 \, \tau_\Tag{dual}^{(t - 1)}}
       {\bigl\Vert \V{s}^{(t)} \bigr\Vert_2 \, \tau_\Tag{prim}^{(t - 1)}} \, .
  \label{eq:eta-alt}
\end{equation}
Indeed, as $\tau_\Tag{dual}^{(t - 1)}$ and $\tau_\Tag{prim}^{(t - 1)}$ do not
depend on $\rho^{(t)}$, $\eta^{(t)}_\Tag{alt}$ is always decreasing function
of $\rho^{(t)}$ while approaching $\eta^{(t)}$ when the algorithm is close to
the solution.  The following algorithm implements our safeguarded strategy to
find $\rho^{(t)}>0$ such that $\eta^{(t)} \approx 1$ or $\eta^{(t)}_\Tag{alt}
\approx 1$.

\smallskip

\textbf{Algorithm 2.}  \emph{Tuning of the augmented penalty parameter
  $\rho$ so that $\eta \approx 1$.}  Choose $\sigma \ge 0$, $\tau > 1$, $\gamma
> 1$ and an initial value for $\rho$, and set $\rho_\Tag{min} = 0$,
$\rho_\Tag{max} = +\infty$.  Then, until convergence, repeat the following
steps:
\begin{enumerate}
\item Update $\V{x}$, $\V{z}$, and $\V{u}$ according to ADMM updating rules.
  Compute $\eta$, defined in \Eq{eq:eta} or in \Eq{eq:eta-alt}, and $\phi$,
  defined in \Eq{eq:convergence-criterion-alt}.
\item If $1/\tau \le \eta \le \tau$ or $\phi < \sigma\,\phi^{(t-1)}$, accept
  the current solution and stop.
\item If $\eta < 1/\tau$, then $\rho$ is too large; let $\rho_\Tag{max} :=
  \rho$ and
  \begin{displaymath}
    \rho := \left\{
    \begin{array}{ll}
      \sqrt{\rho_\Tag{min}\,\rho_\Tag{max}} & \text{if } \rho_\Tag{min} > 0\\
      \rho_\Tag{max}/\gamma & \text{otherwise}\\
    \end{array}\right.
  \end{displaymath}
  then go to step 1.
\item If $\eta > \tau$, then $\rho$ is too small; let
  $\rho_\Tag{min} := \rho$ and
  \begin{displaymath}
    \rho := \left\{
    \begin{array}{ll}
      \sqrt{\rho_\Tag{min}\,\rho_\Tag{max}} & \text{if } \rho_\Tag{max} < \infty\\
      \gamma\,\rho_\Tag{min} & \text{otherwise}\\
    \end{array}\right.
  \end{displaymath}
  then go to step 1. \ensuremath{Box{}}
\end{enumerate}

The following remarks clarify some aspects of this algorithm:
\begin{itemize}
\item To simplify the notations, we dropped the index $t$ of the ADMM
  iteration in the equations of Algorithm~2.  The updating of variables
  (step~1) must be understood as computing $\V{x}^{(t)}$, $\V{z}^{(t)}$, \emph{etc.}\xspace
  given $\V{z}^{(t-1)}$ and $\V{u}^{(t-1)}$, and assuming $\rho^{(t)} = \rho$.
\item Except for the very first ADMM iteration ($t=1$), the initial value for
  $\rho$ is the previous selected value $\rho^{(t-1)}$.  For the first iteration,
  we derive an initial value of $\rho$ such that:
  \begin{align}
    \rho^{(1)}
    & = \MathFuncWithLimits{arg\,min}_{\rho} \Fcost_\DataTag\!\left(\V{z}^{(1)} - \V{u}^{(1)}/\rho\right)
    \notag\\
    & = \frac{
      \V{u}^{(1)\top}\!\!\cdot
      \M{H}\T\!\!\cdot\M{W}\cdot\M{H}\cdot
      \V{u}^{(1)}
    }{
      \V{u}^{(1)\top}\!\!\cdot
      \V{u}^{(1)}
    }
  \end{align}
  with $\V{u}^{(1)} = \nabla\Fcost_\DataTag(\V{z}^{(1)}) =
  \M{H}\T\!\!\cdot\M{W}\cdot(\M{H}\cdot\V{z}^{(1)} - \V{y})$ and which amounts
  to have the best $\proxy{\V{x}}^{(1)}$, defined in \Eq{eq:x-proxy-def}, with
  respect to $\Fcost_\DataTag$.  This choice has the advantage of avoiding an
  initialization with an arbitrary value for $\rho$.  The rule does however
  not yield an efficient strategy for tuning $\rho$ at every iteration.
\item Algorithm~2 is safeguarded in the sense that it maintains a strict
  bracketing $\rho_\Tag{min} < \rho^{(t)}_{*} < \rho_\Tag{max}$ of the
  solution.
\item In step~2 of Algorithm~2: The value of $\rho$ is accepted when $1/\tau
  \le \eta \le \tau$ which, with $\tau > 1$, is how we express that $\eta
  \approx 1$.  We achieved good results with $\tau=1.2$ in our tests.  The
  current value of $\rho$ is also accepted, if the relative reduction in the
  convergence criterion $\phi^{(t)}$, defined in
  \Eq{eq:convergence-criterion-alt}, is better than $\sigma$ with respect to
  the previous iteration.  This shortcut helps to reduce the number of inner
  iterations.  To avoid this shortcut, it is sufficient to take: $\sigma=0$.
  We took $\sigma=0.9$ in our tests.
\item In step~3 and step~4 of Algorithm~2: When $\rho$ has been bracketed by
  $(\rho_\Tag{min},\rho_\Tag{max})$, taking $\rho =
  \sqrt{\rho_\Tag{min}\,\rho_\Tag{max}}$, that is the geometrical means of the
  end points, is similar to a bisection step in a zero finding algorithm.
\item For the first ADMM iteration, the magnitude of $\rho$ is not yet known
  so to avoid too many iterations, we use a larger value of the loop gain
  $\gamma$, say $\gamma = 10$ when $t=1$ and $\gamma = 1.5$ for $t>1$.
\end{itemize}

Another possibility is to always accept an ADMM iteration and simply use the
value of $\eta^{(t)}$ to determine whether $\rho$ should be reduced, kept the
same, or augmented for the next iteration.  For instance:
\begin{equation}
  \rho^{(t+1)} = \left\{
  \begin{array}{ll}
    \gamma \, \rho^{(t)} & \text{if $\eta_t > \tau$,} \\
    \rho^{(t)} / \gamma  & \text{if $\eta_t < 1/\tau$,} \\
    \rho^{(t)}           & \text{else;} \\
  \end{array}
  \right.
  \label{eq:rho-forward-update}
\end{equation}
with $\tau \ge 1$ and $\gamma > 1$.  In words, $\rho$ is augmented (by
multiplying it by a factor $\gamma$) whenever the relative size of the primal
residuals is significantly larger than that of the dual residuals; while
$\rho$ is reduced (by dividing it by a factor $\gamma$) whenever the relative
size of the primal residuals is significantly smaller than that of the dual
residuals.  This strategy is similar to the one described by Boyd \emph{et
  al.}  \cite{Boyd_et_al-2010-method_of_multipliers} except that our
prescription properly scales with the magnitudes of the residuals and of the
objective function of the problem so we expect a better behavior.

\subsection{Debiasing the solution}
\label{sec:debiasing}

One of the drawback of sparsity imposed by means of the $\ell_1$ norm is that
it yields a result which is biased toward zero
\cite{Figueiredo_et_al-2007-gradient_projection}.  For the simulations
presented in Fig.~\ref{fig:integrated-flux}, this bias can be seen on the
recovered spectra in Fig.~\ref{fig:spectra} and in the brightness
distributions of the pixels in Fig.~\ref{fig:mean-flux-histo}.  Since the
sparsity constraint really improves the detection of the sources, the
resulting image can be used to decide where the sources are.  By thresholding
the gray image or the wavelength integrated multi-spectral image resulting
from a reconstruction with a spatial sparsity constraint, we define a sparse
spatial support $\mathcal{S}$ containing all detected sources.  Then, as
proposed in \cite{Wright_et_al-2009-SPARSA}, for the debiasing step, we
minimize the likelihood function $\Fcost_\DataTag(\V{x})$ with non-negative constraints
only over $\V{x}_\mathcal{S}$ defined as the parameters $\V{x}$ restricted to
the support $\mathcal{S}$ while keeping all other parameters equal to zero.
As the sub-matrix $\M{H}_\mathcal{S}$ containing the columns of $\M{H}$
restricted by $\mathcal{S}$ is well conditionned, no additionnal prior is
needed to define the debiased solution.

\begin{figure}
  \centering
  \includegraphics[width=42mm]{star-cluster-true-object}
  \includegraphics[width=42mm]{star-cluster-reconst-sparse}\\
  \includegraphics[width=42mm]{star-cluster-reconst-gray-sparse}
  \includegraphics[width=42mm]{star-cluster-reconst-joint-sparse}
  \caption{Integrated flux for the star cluster. From top to bottom, left to
    right: true object; reconstruction with fully separable sparsity prior
    $\Fcost_\SparseTag(\V{x})$ defined in Eq.~(\ref{eq:sparse-prior}); reconstruction
    assuming a gray object and with spatial sparsity $\Fcost_\SparseTag(\V{g})$ defined
    in Eq.~(\ref{eq:gray-regularization}); reconstruction with joint-sparsity
    prior $\Fcost_\JointTag(\V{x})$ defined in Eq.~(\ref{eq:joint-sparse-prior}).  The
    spectra of the sources encircled by the boxes are shown in
    Fig.~\ref{fig:spectra}.  Axes units are in
    milliarcseconds. \label{fig:integrated-flux}}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=80mm]{spectra-001}
  \includegraphics[width=80mm]{spectra-003}
  \caption{Spectra of two selected sources.  Each panel show the spectra of
    one of the sources encircled by the boxes in
    Fig.~\ref{fig:integrated-flux}. Thick lines are for the true spectra and
    thin lines with markers for the restored spectra.  The open squares
    indicate the reconstruction with fully separable sparsity prior
    $\Fcost_\SparseTag(\V{x})$; the open triangles indicate the reconstruction with
    joint sparsity prior $\Fcost_\JointTag(\V{x})$; the filled triangles indicate the
    restored spectra after debiasing (which are virtually indistinguishable
    from the true ones).  \label{fig:spectra}}
\end{figure}

\section{Results}

To check the proposed algorithm, we simulated a cluster of 50 stars with
random positions and luminosities and with spectra randomly taken from the
library compiled by Jacoby et al. \cite{Jacoby_et_al-1984-stellar_spectra}.
The field of view is $128\times128$ pixels with
$0.5\,\mathrm{milliarcseconds/pixel}$ and we took $100$ spectral channels from
$\lambda = 493\,\mathrm{nm}$ to $\lambda = 507 \,\mathrm{nm}$ by steps of
$\Delta\lambda=0.14\,\mathrm{nm}$. To simulate the observations, we took $100$
random interferometric baselines with a maximum baseline of $180\,\mathrm{m}$.
We added Gaussian white noise to the complex visibilities with a level such
that the maximum signal-to-noise ratio (SNR) is equal to 100.

For the image reconstructions, we considered three different cases: the
reconstruction of a multi-spectral distribution with the regularization
$\Fcost_\SparseTag(\V{x})$ in \Eq{eq:sparse-prior}, or the regularization
$\Fcost_\JointTag(\V{x})$ in \Eq{eq:joint-sparse-prior}, and the reconstruction
of a gray object $\V{g}$ with the regularization $\Fcost_\SparseTag(\V{g})$.  In
order to set the relative weight of the priors, we choose the value of the
hyper-parameter $\mu$ which yields an image which has the least mean square
error with the true distribution.  Once the values of $\mu$ and $\rho$ are
chosen, the reconstruction of a $128\times128\times100$ distribution from
$\sim2\times10^4$ measurements takes about 4 minutes on a GNU/Linux
workstation with a quad-core processor at $3\,\mathrm{GHz}$ and using a
multi-threaded version of FFTW \cite{FFTW3} to compute the discrete Fourier
transforms.

Figure~\ref{fig:integrated-flux} shows the integrated flux, \emph{i.e.}\xspace $\sum_\ell
x_{n,\ell}$, for the true distribution and for the reconstructed ones.  In all
cases, sparsity priors effectively yield a solution with point-like
structures.  However, when there is no trans-spectral constraints, only a few
sources are correctly found and there are many more spurious sources.  When
using $\Fcost_\JointTag(\V{x})$ or assuming a gray object, the estimated integrated
luminosity is much more consistent with that of the true object: all existing
sources are found and the spurious sources are not only less numerous but also
much fainter than the true ones.  This is shown by the brightness
distributions depicted by Fig.~\ref{fig:mean-flux-histo} and
Fig.~\ref{fig:detect}.

Figure~\ref{fig:spectra} shows the spectra of the two stars encircled by boxes
in Fig.~\ref{fig:integrated-flux}, clearly the spectra recovered with
$\Fcost_\JointTag(\V{x})$ (thin curves marked with open triangles) are of much
higher quality than the spectra estimated when treating the spectral channels
independently (thin curves marked with open boxes).  Compared to the true
spectra (thick lines) there is however a small but significant bias in the
spectra obtained with $\Fcost_\JointTag(\V{x})$.  This is not unexpected as the
mixed norm implemented by $\Fcost_\JointTag(\V{x})$ results in an attenuation as
shown by \Eq{eq:joint-sparse-prox}.

As shown by Fig.~\ref{fig:detect}, in the reconstructed gray image or in the
image reconstructed with $\Fcost_\JointTag(\V{x})$, all sources whose mean flux is
greater than 1 are true positive detections while all false positives have a
smaller mean flux.  We therefore select the sources with mean fluxes greater
than this level to apply the debiasing method described in
Sec.~\ref{sec:debiasing} to effectively remove this bias as shown by the thin
curves with filled triangles in Fig.~\ref{fig:spectra}.

In our reconstructions with the regularization $\Fcost_\JointTag(\V{x})$, we found that
$\phi^{(t)} < 10^{-3}$ was a good threshold for the global convergence of the
algorithm and we compared the different strategies proposed in
Sec.~\ref{sec:tuning-rho} to set the augmented penalty parameter $\rho$.  The
evolution of the convergence criterion for some of these strategies is plotted
in Figure~\ref{fig:convergence-criterion}.  With a constant value for $\rho$,
we observed that the rate of convergence is quite sensitive to the value of
the augmented penalty parameter.  Indeed with $\rho=3\times10^{2}$ (which is
the best value we found), the algorithm converged in 455\,s, while it took
1\,071\,s and 847\,s with $\rho=10^{2}$ and $\rho=10^{3}$ respectively.
Although we did not try many different values for the parameters $\tau$ and
$\gamma$, we found that the automatic strategies for setting $\rho$ with
Algorithm~2 and $\eta^{(t)}$ or according to \Eq{eq:rho-forward-update} failed
with their convergence criterion oscillating with $\phi^{(t)} \approx
10^{-2}$.  In fact, in spite of the loss of time due to the number of retries
needed to find a correct value for $\rho$ at each ADMM iteration, we found
that the best strategy was to use $\eta^{(t)}_\Tag{alt}$ in Algorithm~2.  In
this case, it seemed to be better to use a tighter tolerance for
$\eta^{(t)}_\Tag{alt} \approx 1$ as with $\tau=3$ the algorithm converged in
402\,s, while it took only 250\,s with $\tau=1.2$ (see
Fig.~\ref{fig:convergence-criterion}).

\begin{figure}
  \centering
  \includegraphics[width=80mm]{mean-flux-histo}
  \caption{Histograms of the mean fluxes of the sources for the true object
    (in \emph{black}), for the 3-D images restored with fully separable
    sparsity (in \emph{white}) and with joint-sparsity (in \emph{dark gray})
    priors, and for the 2-D gray image restored with sparsity prior (in
    \emph{light gray}).  The vertical scale has been truncated to focus on the
    distributions of the brightest sources. \label{fig:mean-flux-histo}}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=80mm]{true-false-histo}
  \caption{Histograms of the mean flux of the true and false positive
    detection in the reconstructions under joint sparsity and gray sparsity
    priors.  A positive detection is defined as a pixel with non-zero mean
    flux in the reconstruction.  The vertical scale has been truncated to
    focus on the distributions of the true positive
    sources. \label{fig:detect}}
\end{figure}

\begin{figure}
  \centering
  \includegraphics[width=84mm]{convergence-criterion}
  \caption{Evolution of the convergence criterion $\phi^{(t)}$, defined in
  \Eq{eq:convergence-criterion-alt}, for different strategies to
    choose the augmented penalty parameter. Dashed curves are for a constant
    $\rho$.  Solid curves are for $\rho$ automatically set to have
    $\eta^{(t)}_\Tag{alt} \approx 1$. \label{fig:convergence-criterion}}
\end{figure}

\section{Discussion and Perspectives}

We have shown the importance of using trans-spectral constraints to improve
the quality of the restoration of the multi-spectral brightness distribution
$I_\lambda(\V{\theta})$ of an astronomical target from optical interferometric
data.  These results confirm what has been observed for other types of data
(like integral field spectroscopy).

For the moment, our demonstration is restricted to specific objects which are
spatially sparse (\emph{e.g.}\xspace point-like sources) and must be generalized to other
types of spatial distributions.  Being implemented by non-differentiable cost
functions, spatial sparsity requires specific optimization algorithms.  We
have shown that variable splitting by the alternate direction method of
multipliers (ADMM) is suitable to solve the optimization problem in a short
amount of time.  In addition to being able to deal with non-differentiable
criteria, the ADMM method leads to splitting the full problem in sub-problems
that are easier to solve and that may be independent.  This straightforwardly
gives the opportunity of speeding up the code, \emph{e.g.}\xspace by means of
parallelization.  This possibility remains if other priors are used, \emph{e.g.}\xspace to
account for a smooth spatial distribution.

To simplify the problem at hand, we considered that complex visibilities have
been measured.  At optical wavelengths, this is only possible with phase
referencing \cite{Delplancke_at_al-2003-Prima}.  In order to process most
existing interferometric data, we will have to modify the likelihood term
$\Fcost_\DataTag(\V{x})$ and use a non-linear method (\emph{i.e.}\xspace not the linear conjugate
gradients) to update the auxiliary variables $\V{z}$.  The new algorithm that
we proposed, because it splits the two cost functions, $\Fcost_\DataTag(\V{z})$ and
$\Fcost_\PriorTag(\V{x})$, may however be an efficient alternative to the variable
metric method used in Mi\textsc{ra}\xspace \cite{Thiebaut-2008-Marseille,
  Thiebaut-2002-optim_bdec} or the non-linear conjugate gradients method in
\textsc{Bsmem}\xspace \cite{Buscher-1994-BSMEM, Baron_Young-2008-Marseille} which considers
directly the sum of the cost functions.

\section*{Acknowledgments}

This work is supported by the French ANR (\emph{Agence Nationale de la
  Recherche}), Éric Thiébaut and Lo\"ic Denis work for the MiTiV project
(\emph{Méthodes Inverses de Traitement en Imagerie du Vivant},
ANR-09-EMER-008) and Ferréol Soulez is funded by the POLCA project
(\emph{Percées astrophysiques grâce au traitement de données
  interférométriques polychromatiques}, ANR-10-BLAN-0511).

Our algorithm has been implemented and tested with \textsc{Yorick}\xspace
(\myURL{http://yorick.sourceforge.net/}) which is freely available.


\section{Introduction}
\label{sec:Introduction}

In this year we celebrate the 100th anniversary of the Weyl law, which in its earliest version gives the leading term in the asymptotic description of the counting function of the Dirichlet or Neumann Laplacien on a bounded domain in the Euclidean space \cite{wey12a}. Weyl's original intention was to justify the derivation of the Rayleigh-Jeans law \cite{jea1905} for an arbitrary domain, not just boxes. Following works by many mathematicians, among them Courant, Hilbert, Avakumovic, and Levitan, in 1968 H\"{o}rmander obtained a very general Weyl law for elliptic operators \cite{hoer68}. The improved remainders under dynamical assumptions were obtained by Duistermaat-Guillemin on manifolds without boundary \cite{dui75} and by Ivrii, who proved the Weyl conjecture on the second term in asymptotics [see (\ref{eq:WeylClosed}) below] \cite{ivr80a,ivr82a}. In physics literature, higher order terms for smoothed out counting function were introduced already in the 1970s by Balian and Bloch in the context of electromagnetic eigenmodes in cavities with perfectly conducting smooth walls \cite{bal70}. These results were experimentally verified in three and two dimensional microwave cavities \cite{deu95,grae92}. It is intriguing that for typical ``closed'' cavities the Weyl law agrees well even down to the ground state. In more recent developments the Weyl formula has been extended to take into account bouncing ball fluctuations \cite{grae92,sie93}, fractal boundaries \cite{ber79b,bro86}, and ray-splitting \cite{pra96b}.

As soon as a wave-mechanical system is open the eigenvalues turn into resonances. The study of eigenvalues, resonances, and quasibound states has a long tradition in theoretical, numerical, and experimental chaotic scattering (see, for instance, \cite{gas89a} and references therein). One argument for the importance of open systems is the fact that for a measurement it is inevitable to allow interaction with the outside of the system which effectively makes it open. If the system is opened only weakly, the Weyl formula is a very good approximation, as one can already see from the microwave experiments mentioned before. However, if the system is strongly coupled to the environment the resonances cannot be related to individual eigenvalues of the closed systems. If the number of attached channels is finite the framework of random matrix theory (RMT) using non-Hermitian Hamiltonians predicts a separation of nearby resonances \cite{joh91,haa92,pol12}. By varying the coupling resonance trapping effects show up, that is, by
increasing coupling many resonances may become sharper and only a few become much broader \cite{oko03,rot09,per00}. All these effects increase the difficulty of counting resonances compared to the counting of eigenvalues.

If one is interested in the description of the resonances in the semiclassical limit the number of channels diverges and random matrix description breaks down. In this case many quantum mechanical quantities are determined by quantities of the underlying classical dynamics.  For example, Gaspard and Rice \cite{gas89a} derived a semiclassical lower bound on the resonance width and Sridhar observed experimentally that the quantum escape rate for $n$-disk systems coincides with half of the classical escape rate \cite{lu99}. For open dynamical systems with a fractal repeller, Sj\"{o}strand showed that the counting function is polynomially bounded with the power given by the box dimension of the classical repeller \cite{sjo90}.
In systems with only two degrees of freedom, such as considered here, this corresponds to the Hausdorff dimension. The counting function will be introduced in Sec.~\ref{sec:CF}. Numerical investigations in \cite{lin02b} and \cite{lu03} lead to a conjecture that this bound is, in fact, optimal. Numerous numerical and theoretical studies in physics and mathematics followed \cite{zwo99b,gui04,sch04e,non05,kea06b,sjo07}. In particular, the Weyl law has been investigated numerically on maps, the kicked rotator, or the three-disk billiard \cite{lu03,sch04e,non05,kea06b,kop10}. The fractal dimension of the chaotic set also occurs in systems with mixed phase space \cite{kop10,erm10,erm11,ebe10,ped09,ped12,ram09,str04,wie08a,sch09,chr07}. Until now there exists no experimental verification for this law as it is very challenging to extract resonances of strongly open systems in the semiclassical regime. The most recent mathematical advances include fractal upper bounds for several convex bodies \cite{arXnon1} and for
arbitrary manifolds with hyperbolic ends \cite{arXdat1}.

The paradigmatic physical example of an open system with fractal repeller is the $n$-disk system which has been introduced in the context of resonances by Ikawa in mathematics \cite{ika88} and by Gaspard and Rice \cite{gas89b,gas89a} and Cvitanovi{\'c} and Eckhardt \cite{cvi89} in physics. Theoretically it has been studied in the classical, semiclassical, and quantum mechanical regime and there have also been an experimental study on the spectral autocorrelation \cite{pan00a,lu00}. In this article we want to investigate the counting function of a symmetry reduced $n$-disk system experimentally. We realize the ``quantum'' $n$-disk system by a microwave cavity and use the well established equivalence between Helmholtz and Schr\"odinger equation \cite{stoe99} which makes it possible to measure wave characteristics of quantum single-particle systems with a table-top microwave experiment. As we are interested in the transition from a closed to an open system, we chose the five-disk system, as for this system the completely closed system is already sufficiently large for performing experiments. Experimentally, we are restricted to a finite frequency range leading to about 150 resonances. We start from a closed system and by changing a parameter the system opening increases (see Sec.~\ref{sec:ExpSetup}). For small openings the system only couples via few channels, whereas for large opening a ``semiclassical'' coupling will be realized. The resonances are extracted by the method of harmonic inversion (see Sec.~\ref{sec:HI}) and in Sec.~\ref{sec:ExpResults} we present the experimental findings. Concluding remarks are given in Sec.~\ref{sec:Conclusion}.

\section{Eigenvalue and resonance counting functions}
\label{sec:CF}
As we are dealing with a two dimensional system with hard wall potentials we will restrict the discussion of the counting function to this case. If we consider a two dimensional quantum mechanical closed system with hard wall potential, then it is described by the Laplace operator with Dirichlet boundary conditions which has a discrete real spectrum $\Delta \psi_n = -k_n^2\psi_n$. The system's counting function $N(k)$ is then defined as
\begin{equation} \label{eq:DefNkReal}
 N(k):=\#\{k_n: k_n\leq k\}
\end{equation}
If the area of the system is denoted by $A$ and its circumference by $U$ the Weyl law predicts, that the counting function is given by
\begin{equation} \label{eq:WeylClosed}
 N_\mathrm{Weyl}(k) = \frac{ 1 }{ 4 \pi } A k^2\pm \frac{ 1} { 4 \pi } U k+c
\end{equation}
where $\pm$ depends on the boundary condition: $+$ for Neumann and $-$ for Dirichlet boundary conditions. The constant $c$ is defined by the curvature and the corners in the system \cite{bal76}. In general the leading term is proportional to $k^d$, where $d$ is the dimensionality of the system. The Weyl law (\ref{eq:WeylClosed}) corresponds to the counting function (\ref{eq:DefNkReal}) only after appropriate smoothing of the counting function.

By coupling the system to the environment, eigenvalues and eigenfunctions turn into resonances and scattering states. One possible description is via a non-Hermitian effective Hamiltonian \cite{oko03,rot09}. If the coupling is performed only via a few weakly coupled channels all resonances will only acquire a small imaginary part. Thus, counting can be performed for the real part of the resonances and the number of resonances will satisfy the Weyl law (\ref{eq:WeylClosed}). To keep the coupling small and the number of channels constant in the semiclassical limit is not realistic.

Let us first assume that the number of channels $M$ stays constant but the coupling of all channels increases. This corresponds to a RMT where the restriction of Hermiticity is dropped. Its properties have been investigated and for strong coupling a separation of resonance by their imaginary part is given \cite{haa92}. If the system contains $N>M$ resonances then one finds $M$ broad resonances, that is, resonances with large imaginary parts, and $N-M$ resonances with small imaginary parts. Now one has to define whether all resonances are to be counted or only resonances up to a certain imaginary part $C$. The counting statistics of resonances is typically defined by
\begin{equation}\label{eq:DefCounting}
 N(k) := \#\{ \tilde{k}_n: \mathrm{Im}\,\tilde{k}_n > -C, \mathrm{Re}\,\tilde{k}_n \leq k \}
\end{equation}
where $C$ is a fixed finite positive constant, which should not be too small. Thus depending on $C$ one would either count $N$ or $N-M$ resonances. This effect is not restricted to RMT models but holds for any wave system strongly coupled to the environment. It has been phrased as resonance trapping in the framework of the effective Hamiltonian theory (for details, see \cite{rot09} and references therein). It has also been observed experimentally in microwave cavities with variable coupling \cite{per00,stoe02c}. If the maximal wave number $k_\mathrm{max}$ is not too large, this might be observable in the counting function as well (see, e.g., Fig.~5 of \cite{stoe02c}). If waveguides with width $d$ are attached, each of them supports $M_w=d/(\lambda/2)=dk/\pi$ modes, where $\lambda$ is the wavelength. Thus, the number of waveguide modes and therefore also the number of channels increases linearly in the semiclassical limit, whereas the number of eigenvalues of the closed cavity increases quadratically. This would
still lead to a dominating Weyl term of $k^2$. Another possibility is to couple as many channels as eigenvalues (or more) to the system. In the framework of RMT this is related to the Ginibre ensemble \cite{gin65,haa01b}.

Until now we have neglected any additional variation of the real part of the resonances that might be induced by changing the coupling as well as the internal classical dynamics. The structure of this internal dynamics will induce special phase space structures leading to deviations from RMT predictions. The classical phase space of open chaotic systems is characterized by the forward and backward trapped sets. The repeller is the intersection of the fractal sets of trajectories that stay trapped forever in the future and in the past. In the semiclassical limit the wave functions of the long-lived resonances, that is, resonances with small imaginary parts, will localize on the trapped set \cite{erm09,ped12}. Thus, they will avoid the coupling to the channels which is in correspondence with the resonance trapping effect.

As indicated in the Introduction the mathematical and numerical works suggest that the counting statistics in case of open systems with a classical fractal repeller correspond to a fractal Weyl law:
\begin{equation}\label{eq:FractalWeyl}
 N(k) \sim k^{1+d_H},
\end{equation}
where $d_H$ is the reduced fractal dimension of the repeller. We stress that the rigorous results so far prove only an upper bound of this form, or finer bounds in smaller intervals \cite{sjo07,arXnon1} and for arbitrary manifolds with hyperbolic ends \cite{arXdat1}. The numerical papers on maps \cite{non05}, open cavities \cite{lu03}, and microdisk lasers \cite{wie08a} showed the asymptotic behavior.

In this article we investigate experimentally a transition from a closed system to an open five-disk system. Experimentally, one is typically restricted to the wave number $k$ range to $0 \le k \le k_\mathrm{max}$, corresponding to approximate 150 resonances in our case. The billiard of interest consists of two straight walls inducing Dirichlet boundary conditions with an angle of $\theta=36^\circ$. A half circle is attached to only one wall with radius $a$ and distance $d$ from the corner. The position $d$ can be varied but angle and radius of the disk are fixed. In the following section we describe the experimental setup in detail.

\section{Experimental Setup}
\label{sec:ExpSetup}

\begin{figure}
\includegraphics[width=.8\columnwidth]{fig1-Setup.eps}
\caption{\label{fig:setup}
(a) A sketch of the five-disk scattering system is given. The gray shaded part shows the fundamental domain of the symmetry reduction. The radius of the disks is denoted by $a$ and their distance by $R$. (b) A photograph of the experimental setup is presented.
}
\end{figure}

In the experiments we use a flat microwave resonator with a single wire antenna inserted through a hole in the top plate [see Fig.~\ref{fig:setup}(b)]. The baseplate is an aluminum triangle whose largest side is $100$\,cm long. Two side walls with a height of $6$\,mm are set atop; the angle between them is $36^\circ$. Thus the resonator forms the fundamental domain of a five disk system with metallic sidewalls, which act as mirrors reducing the symmetry of the system. The metallic walls induce Dirichlet boundary conditions for the electric field of the TM$_0$ mode, thereby restricting the measurements to a single representation $A_{1u}$ of the underlying symmetry group $C_{5v}$ [see sketch in Fig.~\ref{fig:setup}(a)]. The third side is left open, but additionally covered with a wide strip of microwave absorber to avoid reflections at the open end. Along the longest triangle side we move a half-disk inset of radius $a=19.5\,$cm, with the same height as the side walls, in steps of 10\,
mm. The cover plate is not shown in Fig.~\ref{fig:setup}. An antenna is inserted through a hole in the upper corner of the cover. The radius of the antenna ($r$=0.7\,mm) is much smaller than the wavelength and the antenna is sufficiently short not to touch the bottom plate. The reflection coefficient is measured with a vector-network-analyzer (VNA), revealing the complex $S$ matrix. The height of the cavity $h=6$\,mm (along the vertical $z$ direction) leads to a cutoff frequency of 25\,GHz, below which the $x$- and $y$ component of the electro-magnetic field must vanish and only $E_z$ remains. These electromagnetic modes with vanishing $E_x$ and $E_y$ components are also called TM$_0$ modes. By limiting the frequency range to be analyzed to 24\,GHz, we thus make sure that only the TM$_0$ mode can propagate and the cavity may be considered as two dimensional. The Helmholtz equation describing the electromagnetic wave can then be written as
\begin{equation}\label{eq:Helmholtz}
 \Delta_{xy} \psi_n = -k_n^2\psi_n
\end{equation}
where the wave function $\psi_n = E_z$ and the wavenumber $k_n = 2\pi\nu_n/c$ with $\nu_n$ being the eigenfrequency of the electromagnetic wave. For Dirichlet boundary conditions ($E_z$ must vanish at the side walls) there exists a one-to-one correspondence between wave mechanics and quantum mechanics, that is, between the time independent Helmholtz and Schr\"{o}dinger equation \cite{stoe99}.

The disk positions are characterized by a distance-to-radius parameter $R/a$, the ratio between the distance $R$ of two disks in the full system and the disk radius $a$ (see Fig.~\ref{fig:setup}). The antenna is placed in between the acute angle of the aluminum walls and the movable disk inset. The position is chosen such, that the antenna is not to close to any wall even in the closed case for $R/a =2$. Thus, always a coupling to the interior of the scattering system is guaranteed, where the long-living states live. While increasing the ratio $R/a$ from 2 to 3.9 the fractal dimension of the underlying classical repeller changes and accordingly should the exponent of the counting function. The whole setup is covered by a second aluminum plate that is firmly pressed onto the billiard, leaving no gap between walls and cover plate.

\section{Harmonic Inversion}
\label{sec:HI}

\begin{figure}
\includegraphics[width= 0.45\textwidth]{fig2-Spectra.eps}
\caption{\label{fig:Spec} (Color online)
The plot shows measured reflection spectra for the closed system $R/a=2$ (black solid line),
a system in the transition region with $R/a=2.25$ (red dotted line), and an open system with $R/a = 3.83$ (blue dashed line).}
\end{figure}

Under the assumption that the antenna couples point-like to the system the measured reflection signal is of the form \cite{ste95}
\begin{equation}\label{eq:SumSquaredLorentzians}
  S_{11}(\nu) =1+ \sum\limits_j \frac{\tilde A_j}{\nu^2 - \nu_j^2}
\end{equation}
where the $\nu_j$ are the complex valued resonance positions. As for all the resonances that will be studied in the sequel, one has for $|\textup{Re}(\nu_j)| \gg |\textup{Im}(\nu_j)|$ the approximation
\begin{equation}\label{eq:ApproxLorentzians}
  \frac{\tilde A_j}{\nu^2 - \nu_j^2} = \frac{\tilde A_j}{(\nu - \nu_j)(\nu + \nu_j)} \approx \frac{A_j}{\nu - \nu_j},
\end{equation}
where $A_j= \frac{\tilde A_j}{2\nu_j}$ is good for $\textup{Re}(\nu) >0 $. This leads to
\begin{equation}\label{eq:SumLorentzians}
  S_{11}(\nu) =1+ \sum\limits_j \frac{A_j}{\nu - \nu_j},
\end{equation}
which is more convenient in the following. In Fig.~\ref{fig:Spec} the spectra for three different $R/a$ parameters are shown. For the closed system ($R/a$=2, black solid line) even in the high frequency regime some separate resonances are visible. For $R/a$=2.25 the opening between the half circle and the straight wall is approximately 24\,mm. Thus, the opening supports approximately four modes in the shown frequency range and the resonances can still be recognized but are sufficiently broadened (red dotted line). By further increasing the opening ($R/a$=3.83, blue dashed line) the resonances are strongly overlapping.

In order to extract the resonances from this frequency signal we use the Harmonic Inversion (HI) algorithm as it has been presented by Main et\,al.~\cite{mai00}. To illustrate the challenges of applying this algorithm to experimental data, we first give a short summary of the algorithm described in \cite{mai00} and refer to \cite{mai00,kuh08b} for further details. Afterwards we propose some new tools to circumvent the occurring problems.

The HI algorithm extracts the resonances from the signal in four steps: windowing, truncating the Fourier transformed signal, Pad\'{e} approximation, and finally filtering. In the first step the measured signal is divided into overlapping windows. This is necessary to reduce the number of resonances in each separately analyzed window, which is crucial for applying the Pad\'{e} approximation. In order to reduce boundary effects of the windowing, only resonances within a given range $\Delta\nu_r$ in the center of each of the windows will be taken into account. The windows must therefore overlap correspondingly. In the following $\Delta\nu_r$=1\,GHz with a buffer region of 1\,GHz on each side leading to a total window size of 3\,GHz.

\begin{figure}
\includegraphics[width= 0.45\textwidth]{fig3-FFT.eps}
\caption{\label{fig:FFT} (Color online)
Fast Fourier transform of the reflection signal of the frequency window from 20 to 23\,GHz for the closed system ($R/a =2$, upper black curve) and an open system ($R/a=3.83$, lower orange curve). The unit of the $x$ axis is the number of data points of this discrete time signal. The inset shows the decay for the closed system in a double logarithmic plot. The red straight line correspond to a $1/n$-decay. The light blue shaded region indicates the variation of the cutoff parameter between 80 and 120 data points, which is used in the data analysis described in Sec.~\ref{sec:ExpResults}.}
\end{figure}

For each of these windows a fast Fourier transform (FFT) is applied. Using (\ref{eq:SumLorentzians}) and neglecting the window effects, the FFT signal is given by
\begin{equation}\label{eq:Ct}
 C(t) = \delta(t)+\sum_j d_j e^{-2\pi i\nu_j t}.
\end{equation}
As the system is open, the resonance frequencies have a negative imaginary part and $C(t)$ decays exponentially with $t$. In Fig.~\ref{fig:FFT} examples for two different openings are shown. The decay at the beginning given by a sum of exponentials is due to the widths of the resonances and then changes to a $1/n$ decay which is due to the windowing, which can be seen in the inset. Applying the FFT to a measured window of 3\,GHz width and a step width of 0.1\,MHz we obtain a discrete time series of 30\,000 data points, which is still too large to apply the Pad\'{e} approximation. This problem is solved by truncating the time signal after $N_\mathrm{trunc}$ data points, which is justified by the fact that for high $t$ values the signal is governed by the $1/t$ decay which is due to the discontinuities at the window boundaries. The physical information of the decay is hence contained in the first part of the time signal. As it will be crucial in the sequel that the analyzed time series is of form (\ref{eq:Ct}) it
is important that one cuts off the signal before the $1/t$ decay takes over. On the other hand, $N_\mathrm{trunc}$ should also be chosen sufficiently large such that sufficient information is contained in the truncated signal. Figure~\ref{fig:FFT} shows that for the analyzed microwave spectra and a total window size of 3\,GHz a truncation at $N_\mathrm{trunc} = 80\dots 120$ is reasonable.

After this truncation an equidistant discrete time series of $N_\mathrm{trunc}$ data points $C(\Delta t\cdot n) = \sum_j d_j e^{-2\pi i\nu_j \Delta t\cdot n}$ with $n=0\dots N_\mathrm{trunc}-1$ remains and the central step is now to interpret this signal as a system of nonlinear equations of the form
\begin{equation}\label{eq:EqSys}
  c_n = \sum_{k=1}^{N_\mathrm{trunc}/2} d_k z_k^n \quad \textrm{ with } n=1\dots N_\mathrm{trunc}/2
\end{equation}
which is solved by the Pad\'{e} approximation \cite{mai00, pad92}.

\begin{figure}
\includegraphics[width= 0.45\textwidth]{fig4-ImVsRe.eps}\\
\caption{\label{fig:Filter} (Color online)
This figure shows the resonances of the five disk system with $R/a=3.83$ in the complex plane obtained via the HI for one set of parameters. The light blue circles are resonances that are removed by the filter, whereas the black dots correspond to resonances that are assumed to be true.}
\end{figure}

An obvious problem with this approach is that $N_\mathrm{trunc}$ determines the number of resonances which are returned. Choosing this parameter according to the criteria discussed above $\frac{N_\mathrm{trunc}}{2}$ will be larger than the true number of resonances and the Pad\'{e} approximation will generate spurious resonances which have to be filtered out. There are two filter mechanisms: (i) The Pad\'{e} approximation is performed a second time with a time signal shifted by one point, that is, $C(\Delta t\cdot n)$ with $n = 1\dots N_\mathrm{trunc}$. Resonances that are unstable with respect to this point shift are rejected. (ii) Only resonances whose heights $\frac{|A_j|}{\left|\mathrm{Im}{\nu_j}\right|}$ exceed sufficiently the noise level are accepted. For the data analyzed in this article it has proven to be convenient to set this filter to 0.05. The effect of the filtering is illustrated in Fig.~\ref{fig:Filter}. Without filtering, the HI returns too many resonances; especially there appear resonances
with very small widths which are not realistic for such an open system. After filtering, only the resonances shown as black dots survive and the other resonances (light blue circles) are sorted out. For the remaining resonances a spectral gap and a conglomeration around $\mathrm{Im}(\nu) \approx -0.03$ can be observed. Have in mind that the imaginary part includes terms coming from the absorption in the top, bottom, and side walls and from the antenna. The description of the algorithm above shows that the HI algorithm requires the choice of several parameters, such as $N_\mathrm{trunc}$ or the window overlap. For all these parameters there is a range of plausible values; their exact values, however, are arbitrary. The assumption that the time series is a superposition of decaying exponentials only is another problem. For experimental data this will never be exactly true due to inevitable experimental noise and errors.

\begin{figure}
\includegraphics[width= 0.45\textwidth]{fig5-SpectrumRecalc.eps}\\
\caption{\label{fig:BeWoRec} (Color online)
The top plot shows the measured microwave spectrum ($R/a = 3.83$) (black solid line) as well as two reconstructions obtained by the HI results for two different parameter sets (orange dashed and red dotted lines). The bottom plot shows the corresponding complex resonance positions. The red circles correspond to the resonances of the reconstructed signal shown above as red dotted line and the blue triangles to the orange dashed line. Light blue dots represent the resonance positions of the HI results for all parameter sets. The two vertical dashed lines mark the valid range of the analyzed window.}
\end{figure}

The possible problems occurring with the harmonic inversion of experimental data are visualized in Fig.~\ref{fig:BeWoRec}. The measured microwave spectrum for $R/a=3.83$ (black solid line) was analyzed by the harmonic inversion in order to extract the resonances in the valid range between 21 and 22\,GHz. Therefore, the HI was applied to the data with 160 different parameter sets by varying $N_\mathrm{trunc}$ between 80 and 120 and changing the buffer region on both sides between 1 and 0.98\,GHz, respectively. According to Fig.~\ref{fig:FFT} all these parameter sets are plausible. If the complex resonance positions of these 160 HI results are plotted in one plot (light blue dots in the lower part of Fig.~\ref{fig:BeWoRec}) one observes that the resonances of the different parameter sets form clusters in the complex plane. If one looks, however, at an individual result of a particular parameter set marked by the red circles one observes, that for this parameter set several clusters are missing. For other
parameter sets, however, nearly all the clusters within the valid interval are matched in the resulting resonances. The blue triangles mark one example.

Useful to reject bad HI results is the quality of the reconstruction which is calculated as follows: Suppose one is interested in resonances in the frequency interval from 21 to 22\,GHz as in Fig.~\ref{fig:BeWoRec} and has chosen to take the overlap on each side to be 1\,GHz. The HI algorithm thus analyzes the whole interval from 20 to 23\,GHz, Fourier transforms it, applies the Pad\'{e} approximation, and filters the resonances as described above. As a result it returns a set of resonances and amplitudes $\{(A_j, \nu_j)\}$, where the real parts of the found resonances range from 20 to 23 GHz. One now compares the superposition $\sum_j \frac{A_j}{\nu - \nu_j}$ with the measured signal $S_{11}(\nu)$ between 21 and 22\,GHz. Note that all the resonances in the large window range between 20 and 23 GHz are already contained in the superposition. However, influences of resonances outside this window and very broad resonances which cannot be extracted from the signal are not yet taken into account. We suppose that
these influences can be approximated over the size of the valid region by a complex valued linear function. This background is thus fitted to the difference between superposition and measured signal. The result of the reconstruction (which now contains the superposition of Lorentzians and the fitted background) is plotted in the top part of Fig.~\ref{fig:BeWoRec}. While the reconstructed signal (orange dashed line) of one set of resonances (blue triangles) agrees well with the measured data, the deviation of another reconstruction (red dotted line) using the other set (red circles) is enormous.

Based on these observations it is not sufficient to use only one parameter set for applying the HI on such kind of experimental data, as even the set leading to the best reconstruction might overlook single resonances. Therefore, we perform the HI for several parameter sets, reject those with unconvincing reconstruction by means of their $\chi^2$ value, and average over all the others.

\section{Counting Resonances: Experimental Results}
\label{sec:ExpResults}

In order to study the spectral asymptotics of the five-disk system, we have to determine the counting function (\ref{eq:DefCounting}) from the measured reflection spectrum. In this section we first give the detailed parameters used for the data analysis. Afterwards we present the results on the spectral asymptotics.

The resonance positions which are needed for the calculation of the counting function are extracted from the reflection spectra by the harmonic inversion. Thus according to Sec.~\ref{sec:HI} we start the analysis by decomposing the complete signal ranging from 1 to 25\,GHz into smaller windows. The width of the valid interval is chosen to be 1\,GHz and the buffer region on both sides as 1\,GHz as well. Thus, we obtain in total 22 windows allowing us to extract the resonances between 2 and 24\,GHz. For each of these windows we applied the HI with various different parameter sets for the reasons given in the preceding section. We varied $N_\mathrm{trunc}$ between 80 and 120, the shaded area in Fig.~\ref{fig:FFT}. Additionally, we slightly shifted the buffer region on both sides by 0.02\,GHz. The combination of these variations thus yields 160 different parameter sets.

For each window and each parameter set the HI returned a set of complex resonances and amplitudes. For each of these results the reconstruction had been compared to the original signal as described above and all parameter sets leading to a reconstruction with a $\chi^2$ error worse than the variance of the signal had been rejected. Now the density of states in the valid range of each window is obtained by averaging over the results of all accepted parameter sets. The \emph{averaged} counting function finally was obtained by the averaged density of states of all windows.

\begin{figure}
\includegraphics[width=.8\columnwidth]{fig6-Nk.eps}
\caption{\label{fig:density}(color online)
The counting functions for $R/a$=2, 2.25, and 3.9 (in this order from bottom to top) are plotted in black (histograms). Fits of their slope in the frequency range 15-24\,GHz, corresponding to a $k$ range of 315-500\,mm$^{-1}$ (dotted vertical line) are shown in blue (straight lines). The light orange curve over the lower histogram corresponds to the Weyl formula with 12\% loss for the closed system. Plotted in the inset is the difference between the Weyl formula with 12\% loss and the experimental counting function for the closed system ($R/a$=2).}
\end{figure}

In Fig.~\ref{fig:density} three examples of the experimentally obtained counting functions for different $R/a$-parameters are shown: for the closed ($R/a = 2$) and the most open system ($R/a = 3.9$) as well as for the transition region ($R/a = 2.25$). The closed curve is compared to the prediction of the Weyl law (\ref{eq:WeylClosed}) with area and boundary term. Measuring only with one fixed antenna we cannot expect to find all resonances as several of them have an amplitude of the order of the noise. In closed room temperature aluminum cavities of similar sizes this typically lead to a loss of 5\%-10\%. For $R/a=2$ the orange line shows the Weyl law (\ref{eq:WeylClosed}) with 12\% loss and is in good agreement with the extracted counting function. The inset shows the difference between Weyl law with 12\% loss and counting function. The experimental data fluctuates around the theory with a deviation of less than four resonances, but there is no overall tendency.

For the open systems the fractal Weyl law only predicts the asymptotic exponent of the counting function. To compare our results with this prediction the slope was fitted with a standard regression (blue lines) to our experimental counting functions in the interval 15-24\,GHz (marked by the dotted vertical line). For the fit we need, on the one hand, a sufficiently large range to extract the slope reliably; on the other hand, we need large $k$ to get into the semiclassical regime. In the chosen fit range all counting functions show an approximately linear behavior.

\begin{figure}
\includegraphics[width=.8\columnwidth]{fig7-Dimension.eps}
\caption{\label{fig:exponent} (Color online)
The data points correspond to the fitted exponent of the counting function in dependence of the $R/a$ parameter. The solid line shows the asymptotic exponent $1+d_H$ predicted by the fractal Weyl law and the dashed curve is the same shifted by 0.4 as a guide for the eye. The three squares mark the examples which have already been presented in the previous figures. The darker shaded blue region indicates the $R/a$ values without open channels whereas in the range marked by the lighter shaded blue region only a few open channels (1, \dots, 8) exist.}
\end{figure}

In Fig.~\ref{fig:exponent} the slopes of the counting functions are plotted versus their $R/a$ parameter. For the quantum-mechanically closed system the classical Weyl law predicts a value of 2. As long as the slit between disk and metallic wall is smaller than half a wavelength of the maximal frequency 24\,GHz, the system is quantum mechanically closed and couples only by tunneling to the exterior. This region is highlighted by the darker blue shading and an exponent of 2 is still expected within the wave number range investigated. For large $R/a$ parameter the fractal Weyl law predicts a fractional exponent. The black solid line shows the predicted exponent $1+d_H$ of the fractal Weyl law, where $d_H$ is the reduced fractal dimension of the repeller. This dimension has been calculated via the topological pressure and cycle expansion in dependence of $R/a$ (see \cite{gas89b,gas89a,cvi89} for details). Due to pruning effects this method fails to calculate $d_H$ for small $R/a$ parameters; thus, the black theoretical line stops at $R/a=2.4$. In between there will be a transition region. In Fig.~\ref{fig:exponent} the region between one and eight open channels is shaded in light blue. Experimentally, the expected start value of about 2 for the closed system and a smooth transition to lower noninteger values between 1 and 2 is seen. For the quantum-mechanically closed system the fitted exponent agrees well with the predicted value of 2 of the classical Weyl law. In the transition region with a few open channels the exponent decreases smoothly. For the open system the fitted exponents do not match the theoretically predicted curve, but the parametric dependence still can be seen though the theoretical curve has to be shifted down by 0.4 in order to match the fitted exponents. However, the fitted exponents are significantly larger than one and definitely non integer.

A possible explanation for the down shifted fitted exponents could be the following. Suppose that the probability that the harmonic inversion overlooks some resonances, increases when the resonances become more and more overlapping. The density of states for the five-disk system is supposed to increase like $k^{d_H}$; thus, with increasing $k$ the resonances become stronger overlapping and the supposed loss of resonances increases as well. Such a  loss which increases with $k$ would, however, lead to a systematically lower fitted exponent. Even if the procedure of averaging over many HI parameter sets and checking the reconstruction significantly increases the reliability of the results, those ambiguities will remain. This process might also explain why for the closed systems this shift does not occur: As the resonance for the more closed systems are significantly less overlapping even for high frequencies (see Fig.~\ref{fig:Spec}) the loss mechanism might not yet have set in.

Moreover, there are, of course, the experimental uncertainties influencing the findings: The exponential behavior of the counting function is predicted in the semiclassical limit, corresponding to the limit of infinite frequencies. Experimentally, we are, however, restricted to a finite frequency and width range. Additionally, the resonance structure of such open microwave systems is very sensitive to inevitable reflections on imperfect absorbers.

\section{Conclusion}
\label{sec:Conclusion}

In this article we experimentally examined the behavior of the counting function in a five-disk microwave system which was successively transformed from a closed to an open system with a classically fractal repeller. The counting function was obtained from the measured spectrum by the harmonic inversion. We pointed out the occurring problems and proposed some tools to circumvent them. In the closed regime we found good agreement between the extracted counting function and theory. In the range of many open modes the functional dependence showed an agreement with theoretical predictions, though the experimental exponents had been found to be approximately 0.4 below the prediction. This certainly is in agreement with mathematics upper bounds though falls short of demonstrating asymptotics (\ref{eq:FractalWeyl}). It is, however, clear that the growth of resonances is not linear as predicted in \cite{gas89b}. In between there seems to exist a smooth monotonic transition between the two regimes, exhibiting the existence of a small number of open channels for the investigated wave-number range.

\begin{acknowledgments}
We are thankful for St\'{e}phane Nonnenmacher and Bruno Eckhardt for intensive discussions. This work was supported by the Deutsche Forschungsgemeinschaft via an individual grant and the Forschergruppe 760 `Scattering systems with complex dynamics'. TW acknowledges financial support of the `German National Academic Foundation'.
\end{acknowledgments}