TITLE: Showing a 2-by-2 matrix is a contraction
QUESTION [2 upvotes]: Let $S\subseteq\mathbb{T}:=\{z\in\mathbb{C}:\vert z\vert=1\}$ be a compact set such that $\operatorname{conv}S\supseteq\{z\in\mathbb{C}:\vert z\vert\leq\frac{1}{\sqrt{2}}\}$ and $B\in M_2(\mathbb{C})$. Suppose $I_2+\lambda B+\overline{\lambda}B^*\geq 0$ for all $\lambda\in S$. Is $B$ a contraction?

Comments: I can see this question has an affirmative answer for some particular choice of $S$. So I am hoping this may be true. But I am unable to prove it in general.
I apologize that I forgot to mention $S\subseteq\mathbb{T}$ before.
Any comment is highly appreciated.

REPLY [1 votes]: Taking $B=\begin{pmatrix}0&\sqrt{2}\\0&0\end{pmatrix}$ and $S=\{z\in\mathbb{C}, |z|\le \frac{1}{\sqrt{2}}\}$ is a counterexample. The bound is $\sigma_1(B)\le \sqrt{2}$.
For any matrix $B,M\in \mathbb{M}_n(\mathbb{C})$ the numerical range of $B$ is the set $W(B)=\{x^*Bx, x\in \mathbb{C}^n, x^*x=1\}$. Also $W(\dfrac{M+M^*}{2})=\mathfrak{R}(W(M))$,  $W(\lambda B)=\lambda W(B)$ (rotation and homothety) and $\sigma_1(M)\le 2w(M)$ where $w(M)$ is the numerical radius of $M$, (these are known),
with $\lambda B=M$. We see that if $\dfrac{M+M^*}{2}\ge -0.5I$, for any $\lambda\in S$, then $w(M)\le 0.5$. Say $\lambda=e^{i\theta}\frac{1}{\sqrt{2}}$, if there is a point $z\in W(M)$ with $|z|>0.5$ applying a certain rotation  to $B$ this point intersect the $y=0$ line contradicting the fact that the real part of $W(M)$ is in $ [-0.5;*[$.


TITLE: Equation for $N$-body problem using Jacobi Coordinate
QUESTION [0 upvotes]: For reference on Jacobi Coordinate used for solving 2-Body problem, I referred Wikipedia Jacobi Coordinate, and on looking at those equation
I can't get the meaning of the symbol q in the equation for Vector to center of mass,
Rjk = ( mjqj + mkqk )( mj + mk )
It would be helpful if some one could help.

REPLY [1 votes]: It should be an r in the notation that they used before. These are the position vectors of the two objects. And a hint: If you are confused next time, try finding different sources. I'm sure the definition of the center of mass vector can be found elsewhere.


TITLE: Is it possible for an integral between $a$ and $a$ to have a value other than 0
QUESTION [4 upvotes]: Does there exist a function $f: \mathbb{R} \to \mathbb{R}$ with a constant $a \in \mathbb{R}$ such that
$$\int_{a}^{a} f(x) \, dx  \neq 0 \quad$$
holds?

REPLY [4 votes]: As others have explained, in general the answer is 'no'.
However in physics it is sometimes nice to have something 'function-like' that is zero everywhere but in one point (a null set) and still has a integral of 1. This is not a function but a more generalized object - a distribution, the Dirac Distribution.


TITLE: If $p : E \to B$ is a covering map and $E$ is simply connected, then each fiber of $p$ has the cardinality of that of $\pi_1(B)$.
QUESTION [1 upvotes]: Can anyone provide a source or a proof of the following fact:
If $p : E \to B$ is a covering map and $E$ is simply connected, then each fiber of $p$ has the cardinality of that of $\pi_1(B)$?

REPLY [0 votes]: I like the view that algebraic topology is about modelling topology by algebra; in these terms, a useful model of a covering map  $p: E \to B$,  given in Topology and Groupoids,  and, following earlier editions of that book,  to some extent in Peter May's "Concise ...." book, is the "covering morphism" of fundamental groupoids  $\pi_1 p: \pi_1 E \to \pi_1 B$.
A morphism $q: H \to G$ of groupoids is a covering morphism if for each object $x$ of $H$ and element $g: q(x) \to y$ in $G$ there is a unique element $h$ in $H$ starting at $x$ such that $q(h)=g$. It is helpful to draw pictures of this.
The traditional arguments show that a covering map determines a covering morphism, and a version of traditional arguments show that if the space $B$ has appropriate local conditions then a covering morphism $H \to \pi_1 B$ determines a topology on $E= Ob(H)$ so that $Ob(q) :E \to B$ is a covering map with $\pi_1 E \cong H$.
The standard argument can now be given in the algebraic model, as explained in the books and by squirrel,  to show the result required.


text:
\begin{document}

\title{Overcomplete representation in a hierarchical Bayesian framework}

\author[1]{Monica Pragliola\thanks{monica.pragliola2@unibo.it}}
\author[2]{Daniela Calvetti\thanks{dxc57@case.edu}}
\author[2]{Erkki Somersalo\thanks{ejs49@case.edu}}
\affil[1]{Department of Mathematics, University of Bologna, Italy}
\affil[2]{Case Western Reserve University, Department of Mathematics, Applied Mathematics and Statistics, 10900 Euclid Avenue, Cleveland, OH 4410, USA}

\maketitle

\begin{abstract}
	A common task in inverse problems and imaging is finding a solution that is sparse, in the sense that most of its components vanish. In the framework of compressed sensing, general results guaranteeing exact recovery have been proven. In practice, sparse solutions are often computed combining $\ell_1$-penalized least squares optimization with an appropriate numerical scheme to accomplish the task. A computationally efficient alternative for finding sparse solutions to linear inverse problems is provided by Bayesian hierarchical models, in which the sparsity is encoded by defining a conditionally Gaussian prior model with the prior parameter obeying a generalized gamma distribution. An iterative alternating sequential (IAS) algorithm has been demonstrated to lead to a computationally efficient scheme, and combined with Krylov subspace iterations with an early termination condition, the approach is particularly well suited for large scale problems. Here the Bayesian approach to sparsity is extended to problems whose solution allows a sparse coding in an overcomplete system such as composite frames. It is shown that among the multiple possible representations of the unknown, the IAS algorithm, and in particular, a hybrid version of it, is effectively identifying the most sparse solution. Computed examples show that the method is particularly well suited not only for traditional imaging applications but also for dictionary learning problems in the framework of machine learning.
\end{abstract}

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

Sparsity promoting methods and algorithms for inverse problems and imaging applications have been extensively studied in the past decades, and they continue to be a very active field of research. The interest in compressed sensing has motivated a significant part of the works on the topic.
The starting point of many sparse reconstruction problems is a  \emph{dictionary},
intended as a collection of elements in the ambient space referred to as \emph{atoms} \cite{mallat}, used to represent the unknown quantity of interest.  The dictionary may be selected according to some \emph{a priori} information available on the problem of interest or, alternatively, its formation can be data-driven - see, e.g., \cite{dict1,dict2,dict3}. Recently, approaches aimed at learning the dictionary while jointly recovering the signal have also been developed \cite{chamb}.
Typically, the cardinality of the dictionary is significantly larger than the dimension of the ambient space.
When the atoms in the dictionary do not form a basis for the ambient space, the dictionary is called \emph{redundant} or \emph{overcomplete}. The use of redundant dictionaries has proved to be a useful strategy in terms of artifact reduction, especially in the framework of signal denoising problems \cite{SED,SFM}.

Let $x\in\R^n$ be an unknown signal, and let $W=\{w_i\}_{i=1}^N$ be a dictionary with atoms $w_i\in\R^n$. We arrange the atoms as columns of the dictionary matrix, $\mW\in\R^{n \times N}$, with $n\ll N$, and refer to this matrix as the dictionary. In a synthesis perspective, a sparse reconstruction problem is the task of recovering a sparse vector $\alpha\in\R^N$, with most of its components vanishing,  that represents the original signal in terms of $\mW$, $x = \mW\alpha$, starting from a corrupted and possibly poorly sampled indirect observation  $b\in\R^{m}$ of $x$, with $m \leq n$. Assuming that the observation is linear in $x$ and the noise is additive, the sparse dictionary representation can be formulated as an optimization problem of the form
\begin{equation}
\mbox{minimize  $\|\alpha\|_0$  such that  $b = \mA x+\varepsilon$, with $x = \mW \alpha$,}
\label{eq:cs}
\end{equation}
where $\|\alpha\|_0=\text{card}(\text{supp}(\alpha))$,  $\mA\in\R^{m \times n}$ is the forward model operator, and $\varepsilon\in\R^m$ is an additive noise vector. Addressing the minimization problem \eqref{eq:cs}  directly is a challenge due to its NP-hardness, thus explaining the need for alternative approaches. A strategy which has been widely explored and goes under the name of \emph{basis pursuit}  replaces the $\ell_0$-(semi)norm with its $\ell_1$ convex relaxation \cite{donoho}.

When the signal itself is known to be compressible, i.e., $\mathrm{card}\left(\left\{x_i \mid |x_i| < \epsilon\right\} \right)\ll n$ with $\epsilon>0$ arbitrarily small, and the forward model operator $\mA$ satisfies the {\em restricted isometry property} (RIP) condition, an optimal bound for the error $\|x - \hat{x}\|_2$, with $\hat{x}$ denoting the recovered signal, has been derived. Moreover, if $x$ is sufficiently sparse, the signal can be recovered exactly \cite{candes1}. When the signal itself is not sparse, but it allows a sparse or compressible representation in a given dictionary, the exact recovery results still hold, provided that $\mA$ satisfies a {\em restricted isometry property adapted to a dictionary} (D-RIP) condition \cite{drip}. Given the theoretical motivation, a significant amount of research is devoted to identifying classes of operators to which these results could be applied.

Besides the convex approaches, $\ell_p$-norms with $p<1$ have also been considered in place of the $\ell_0$-(semi)norm in problem \eqref{eq:cs}, as they are known to promote sparsity more strongly than the case $p=1$. Nonetheless, the presence of local minima is a clear limitation to the reliability of non-convex strategies.

The sparse reconstruction problem for linear inverse problems allows a natural formulation in the Bayesian computational framework, with the notion of sparsity promoting priors.

In a number of previous works \cite{CSS,CPrSS,CPrS}, the  recovery of a sparse signal $x$ - or of a signal admitting a sparse representation in a given basis - has been addressed by modeling its entries in a hierarchical Bayesian framework as conditionally Gaussian random variables with unknown variances, with a generalized gamma hyperprior distribution.
The sparsity promotion and the convexity properties of the corresponding class of hypermodels have been  studied in \cite{CSS,CPrSS}. The derived results,
and in particular the considerations on the convexity properties of the resulting maximum a posteriori estimation problems, motivated the introduction of a hybrid hypermodel combining the strong sparsity promotion that typically characterizes non-convex settings with the convexity guarantees \cite{CPrS}.

In this article, the hierarchical Bayesian framework outlined in the previous articles, and in particular the use of the hybrid algorithm of \cite{CPrS}, is extended to address sparse recovery problems in presence of redundant dictionaries. We consider a version of the iterated alternating sequential (IAS) algorithm that combines ideas from the Bayesian inference and iterative Krylov subspace methods, suitable for large scale problems, and is therefore particularly attractive for problems with large dictionaries. Numerical examples demonstrate the computational efficiency of the approach, and most importantly, show that for composite frame dictionaries, where each subdictionary would provide a sufficient representation of the signal, the method is capable of identifying an optimally sparse representation.

\section{Hierarchical Bayesian formulation}
\label{sec:bayes}
Consider the linear inverse problem
\begin{equation}
\label{eq:model}
b = \mA x + \varepsilon\,,\; \varepsilon \sim \mathcal{N}(0,\mSigma) \quad\text{ such that }\quad x = \mW \alpha\,,
\end{equation}
where $\mA\in\R^{m \times n}$, with $m \leq n$, is the known forward model operator, $x\in\R^n$ is the unknown of interest and $\mSigma\in\R^{m \times m}$ is the symmetric positive definite covariance matrix of the additive Gaussian noise. In addition, we assume that $x$ admits a representation in the redundant dictionary $\mW\in \R^{n \times N}$, with $n \ll N$, where the unknown vector $\alpha\in\R^N$ is sparse, i.e. $\|\alpha\|_0 \ll N$, either because $x$ can be naturally described by few atoms in the dictionary or because $x$ needs to be compressed.
\\In a number of previous contributions \cite{CSS,CPrSS,CPrS}, the \emph{a priori} sparsity belief on the unknown has been exploited by modeling its entries as independent random variables following a \emph{conditionally Gaussian} distribution, i.e.,
\begin{equation}
\label{eq:cond gauss}
\alpha_j\mid \theta_j \sim \mathcal{N}(0,\theta_j)\,,\; 1 \leq j \leq N\,,
\end{equation}
or, in an equivalent compact form,
\begin{equation}
\label{eq:cond gauss 2}
\alpha\mid \theta \sim \mathcal{N}(0,\mD_{\theta})\,,\; \mD_{\theta}=\mathrm{diag}\left(\theta_1,\ldots,\theta_N\right)\in\R^{N\times N}\,.
\end{equation}
The conditional Gaussian prior on $\alpha$ given the vector $\theta$ takes the form
\begin{align}
\label{eq:cond prior}
\pi_{\alpha\mid\theta} (\alpha\mid \theta) \:{\propto}\:&\displaystyle{ \frac{1}{\prod_{j=1}^{N}\sqrt{\theta_j}}\exp\bigg({-\frac{1}{2}\| \mD_{\theta}^{-1/2} \alpha \|^{2}}\bigg)}\\
\nonumber
\:{=}\:&\displaystyle{\exp\bigg({-\frac{1}{2}\| \mD_{\theta}^{-1/2} \alpha \|^{2}} -\frac{1}{2}\sum_{j=1}^N \log\theta_j\bigg)}\,.
\end{align}
According to the Bayesian paradigm, the unknown vector of variances $\theta$ is also modeled as a random variable. The \emph{a priori} beliefs about $\theta$ are encoded in the \emph{hyperprior} $\pi_{\theta}(\theta)$, and the joint prior on the coupled vector of unknowns $(\alpha,\theta)$ reads
\begin{eqnarray}
\label{eq:joint pr}
\pi_{(\alpha,\theta)}(\alpha,\theta) = \pi_{\alpha\mid\theta}(\alpha\mid \theta)\pi_{\theta}(\theta)\,.
\end{eqnarray}
In \cite{CPrSS}, the authors propose to model the unknown variances $\theta_j$ as mutually independent random variables following a generalized gamma distribution,
\begin{equation}
\pi_{\theta}(\theta) =
\pi_{\theta}(\theta \mid r,\beta,\vartheta) = \frac{|r|^n}{\Gamma(\beta)^n} \prod_{j=1}^{N}\frac 1{\vartheta_j} \left(\frac{\theta_j}{\vartheta_{j}}\right)^{r\beta - 1} \exp\bigg(-\left(\frac{{\theta_j}}{{\vartheta_{j}}}\right)^r\;\bigg)\,,
\label{eq:hyper}
\end{equation}
where $r\in \R\setminus\{0\}$, $\beta >0$, $\vartheta_j>0$. This choice is motivated by the observation that generalized gamma distributions tend to favor values which are close to the expected value while also allowing for few outliers very far from the mean. Presumably, the outlier variances give rise to the few non-zero values of $\alpha$, or values above a tiny threshold.

The information about the observation process is encoded in the likelihood distribution, which in view of the additive Gaussian noise model, takes the form
\begin{equation}
\label{eq:likelihood}
\pi_{b\mid \alpha}(b \mid \alpha) \propto \exp\left(-\frac{1}{2}\|\mS (\mA \mW \alpha - b)\|_2^{2}\right)\,,
\end{equation}
where $\mS$ is the Cholesky factor of the precision matrix $\mSigma^{-1}$, i.e. $\mSigma^{-1} = \mS^{\mT}\mS$. If the matrix $\mSigma$ and thereby $\mS$, are known, without loss of generality we can assume the noise to be white, i.e. $\mSigma = \mI$, because it can be whitened by a linear transform on $\mA$ and $b$, namely
\begin{equation}
\mA \longrightarrow \mS \mA\,,\quad b\longrightarrow \mS b\,.
\end{equation}
Under the white normal noise assumption, the likelihood distribution is of the form
\begin{equation}
\label{eq:likelihood 2}
\pi_{b\mid \alpha}(b \mid \alpha) \propto \exp\left(-\frac{1}{2}\| \mA \mW \alpha - b\|_2^{2}\right)\,.
\end{equation}
The conditional prior and the hyperprior are coupled to the posterior distribution via Bayes' formula, yielding the following expression for the posterior distribution
\begin{equation}
\label{eq:bayes}
\pi_{(\alpha,\theta)\mid b}(\alpha,\theta \mid b) \propto \pi_{b \mid \alpha}(b \mid \alpha)\pi_{(\alpha,\theta)}(\alpha \mid \theta) \pi_{\theta}(\theta).
\end{equation}
In the Bayesian framework, the posterior distribution is the complete solution to the inverse problem, that can be used to produce representative estimates of the unknown of interest, and quantify the uncertainty.
Here, we chose to summarize the posterior with the Maximum A Posteriori (MAP) estimate,
\begin{equation}
(\alpha^{*},\theta^*) \in \arg\max_{\alpha,\theta}\left\{\pi_{(\alpha,\theta)\mid b}(\alpha,\theta \mid b)\right\}
\end{equation}
or equivalently, by taking the negative logarithm of the density and ignoring the additive constants,
\begin{equation}
\label{eq:map}
(\alpha^{*},\theta^*) \in \arg\min_{\alpha,\theta}\left\{\mathcal{F}(\alpha,\theta)\right\},
\end{equation}
where
\begin{align}
\nonumber
&\mathcal{F}(\alpha,\theta) = \mathcal{F}(\alpha,\theta \mid r,\vartheta,\beta)\\
=&\quad \lefteqn{ \phantom{\frac 12 \|b - \mA  \mW\alpha\| ^2 + \frac 12 \sum_{j=1}^N\frac{\alpha_j^2}{\theta_j} }}
{\displaystyle\frac 12} \|b - \mA  \mW\alpha\|^2+
\underbrace{\frac 12\sum_{j=1}^N\frac{\alpha_j^2}{\theta_j}
	-\eta
	\sum_{j=1}^N\log\frac{\theta_j}{\vartheta_j} +\sum_{j=1}^N\left(\frac{\theta_j}{\vartheta_j}\right)^r}_{\mathcal{P}(\alpha,\theta \mid r,\beta,\vartheta)}\,,\quad \eta = \bigg(r\beta - \frac 32\bigg). \label{eq:obj map}
\end{align}
Echoing the terminology of classical regularization schemes, we refer to $\mathcal{P}(\alpha,\theta \mid r,\beta,\vartheta)$ as the \emph{penalty term}.

\section{The IAS algorithm}
\label{sec:ias}
The search for the minimizer of the MAP objective function in \eqref{eq:obj map} is carried out with the {\em global hybrid} scheme introduced in \cite{CPrS}, based  on the iterative alternating sequential (IAS) algorithm described below.  Details of the hybrid scheme that ensues are reviewed in Section \ref{sec:ghias}.

Given a suitable initialization of the variances $\theta^{0}$,  at each iteration step the IAS algorithm updates the iterates $\alpha^t$, $\theta^t$ by solving the minimization problem in alternating directions, that is
\begin{equation}
\label{eq:ias}
\alpha^{t+1}\in\arg\min_{\alpha}\left\{\mathcal{F}(\alpha,\theta^{t})\right\}\,,\quad
\theta^{t+1}\in\arg\min_{\theta}\left\{\mathcal{F}(\alpha^{t+1},\theta)\right\}\,.
\end{equation}
Because of the particular form of the objective function, both variables can be updated efficiently as follows.

\paragraph{\textbf{Update of }$\bm{\alpha}$}
The $\alpha$-update reduces to solving a quadratic minimization problem, i.e.,
\begin{equation}
\label{eq:get alpha}
\alpha^{t+1}\in\arg\min_{\alpha}   \left\{ \|b-\mA\mW\alpha\|_2^2  + \|\mD_\theta^{-1/2}\alpha\|^2 \right\},\quad \theta = \theta^t,

\end{equation}
or, equivalently, finding the solution in the least squares sense of the linear system
\begin{equation}
\label{eq:get alpha LS}
\begin{bmatrix}
\mA\mW\\
\mD_{\theta}^{-1/2}
\end{bmatrix}\alpha = \begin{bmatrix}
b\\
0
\end{bmatrix}\,.
\end{equation}

After performing the change of variable
\begin{equation}
\mD_{\theta}^{-1/2}\alpha = \gamma\,,
\end{equation}
we can write  \eqref{eq:get alpha LS} as
\begin{equation}
\label{eq:get gamma}
\begin{bmatrix}
\mA\mW\mD_{\theta}^{1/2}\\
\mI
\end{bmatrix}\gamma = \begin{bmatrix}
b\\
0
\end{bmatrix}\,,
\end{equation}
where $\mI$ is an $n\times n$ unit matrix. The solution of this least squares problem is also Tikhonov regularized solution of
\begin{equation}
\mA\mW\mD_{\theta}^{1/2}\gamma = b,\,\; \alpha = \mD_{\theta}^{1/2}\gamma\,,
\end{equation}
with regularization parameter equal to one. An alternative to Tikh\-on\-ov regularization yielding a similar solution is to solve the underlying linear system with an iterative solver equipped with an early stopping criterion. The stopping condition is usually based on a variant of Morozov discrepancy principle, whereas the iterations terminate as soon as the discrepancy is of the order of the observation noise. In the statistical framework, under the Gaussian noise assumption, the noise level can be expressed in terms of the standard deviation of the noise.  In our case, where we assume  $m$-dimensional white noise, this quantity is equal to $\sqrt{m}$.   Following \cite{SIREV,CPrSS}, we solve the linear system using the Conjugate Gradient for Least Squares (CGLS) algorithm with the early stopping at noise level $\sqrt{m}$; see \cite{SIREV}  for more details.

\paragraph{\textbf{Update of }$\bm{\theta}$} Due to the mutual independence of the entries of $\theta$, each variance $\theta_j$ can be updated separately by imposing the
component-wise first order optimality condition on \eqref{eq:obj map}. More specifically, $\theta_j^{t+1}$ is the solution of the non-linear equation
\begin{equation}
\label{eq:get th}
\frac{\partial{\mathcal F}}{\partial\theta_j} =
-\frac{1}{2} \frac{\alpha_j^2}{\theta_j^2}-\left(r \beta -\frac{3}{2}\right)\frac{1}{\theta_j} + r \frac{\theta_j^{r-1}}{\vartheta_j^r} = 0\,,\; \alpha=\alpha^{t+1}\,.
\end{equation}
For some values of $r$, e.g., $r=\pm 1$, \eqref{eq:get th} admits an analytic solution. However, in general we need to solve it numerically. It was shown in \cite{CPrSS} that after the changes of variables  $\theta_j = \vartheta_j \xi_j$, $\alpha_j = \sqrt{\vartheta_j} z_j$, we may write $\xi_j = \varphi(|z_j|)$, and  via implicit differentiation, the function $\varphi$ satisfies the initial value problem
\begin{equation}\label{get xi 2}
\varphi'(z) = \frac{2z \varphi(z)}{2r^2 \varphi(z)^{r+1}+z^2} , \quad \varphi(0) = \left(\frac{\eta}{r}\right)^{1/r}.
\end{equation}
Therefore the updated value of $\theta_j$ can be computed by a numerical time integrator. Since the same type of differential equation is satisfied by all components, an efficient way to update $\theta$ is to sort the current values $z_j$ in an ascending order, and integrate sequentially over the gaps between the values by a suitable time integrator.

We point out that unlike in the formally similar alternating direction method for multipliers (ADMM) algorithm \cite{boyd2011distributed} that is often used to solve regularized inverse problems with sparsity promoting priors, the IAS algorithm does not require the introduction of an artificial decoupling term of the fidelity and penalty terms, as the partial decoupling in IAS is automatic and exact.

From the point of view of statistical analysis, the proposed algorithm is not aiming at exploring the posterior density, and the MAP estimate might not be the best single point estimate to characterize the posterior. Other strategies of interest include the marginalization of the posterior density with respect to the hyperparameter $\theta$, or estimating an optimal $\theta$ by first marginalizing $\alpha$. These alternative strategies have been discussed in literature, see, e.g., \cite{Vidal} for a recent and comprehensive contribution.

\section{Parameter selection strategies}
Before presenting the details of the hybrid scheme used in the numerical tests, we briefly review some of the main results related to the selection  of the hyperparameters $(r,\beta,\vartheta)$ appearing in the expression of the hyperprior in \eqref{eq:hyper}.

We start recalling a theorem, whose proof can be found in \cite{CPrSS}, summarizing how $r$ and $\beta$ affect the convexity properties of the functional $\mathcal{F}$ .

\begin{theorem}\label{th:convexity}
	Let $\beta>0$ and $r\neq 0$, and let $\mathcal{F}(\alpha,\theta)$ be the objective function for the minimization problem in (\ref{eq:map}).
	\begin{itemize}
		\item[(a)]
		If $r\geq 1$ and $\eta = r\beta  -3/2>0$, the function $\mathcal{F}(\alpha,\theta)$ is globally convex.
		\item[(b)]
		If $0<r<1$ and $\eta = r\beta  -3/2>0$, or, if $r<0$ and $\beta>0$,
		the function $\mathcal{F}(\alpha,\theta)$ is convex provided that
		\begin{equation} \label{th bar}
		\theta_j < \overline{\theta}=\vartheta_j \left(\frac{\eta}{r|r-1|}\right)^{1/r}.
		\end{equation}
	\end{itemize}
\end{theorem}

The convexity of the MAP objective function, guaranteed for $r \geq 1$, is very convenient, however some of the configurations attained for $r<1$ can be very attractive in terms of sparsity promotion and rate of convergence. To better understand the connection between the parameter of the hyperprior and sparsity, consider the updating formula  \eqref{get xi 2}, expressing $\theta_j$ as a function of $\alpha_j$ as
\begin{equation}\label{def g}
\theta_j =  g_j(\alpha_j) =  \vartheta_j \varphi\left(\frac{|\alpha_j|}{\sqrt{\vartheta_j}}\right)\,.
\end{equation}
We review some recent results \cite{CHPS,CSS,CPrSS} about the connections that can be drawn between the generalized gamma hyperpriors and classical sparsity promoting penalty terms, assuming that $(\theta_j,\alpha_j)$ satisfies the above identity.
\begin{itemize}
	\item[(i)] For the gamma hypermodel, i.e. $r=1$, as $\eta=r\beta-\frac{3}{2}\to 0^{+}$ the penalty term approaches a weighted $\ell_1$-penalty term \cite{CHPS,CSS},
	\begin{equation}
	\label{eq:lim ell1}
	\lim_{\eta \to 0^+}\mathcal{P}\left(\alpha,g(\alpha)\mid 1,\frac{3}{2}+\eta,\vartheta\right) = \sqrt{2} \sum_{j=1}^n \frac{|\alpha_j|}{\sqrt{\vartheta_j}}\,.
	\end{equation}
	\item[(ii)] If $r\beta=\frac{3}{2}$,
	the penalty term coincides with the weighted $\ell_p$-norm, with $p =2r/(r+1)$ \cite{CPrSS},
	\begin{equation}
	{\mathcal P}\left(\alpha, g(\alpha)\mid r,\frac{3}{2r},\vartheta\right)  = C_r\sum_{j=1}^n \frac{|\alpha_j|^p}{\sqrt{\vartheta_j}^p}, \quad C_r = \frac{r+1}{(2r)^{r/(r+1)}}.
	\end{equation}
	\item[(iii)]For the inverse gamma hypermodel, corresponding to $r=-1$, the penalty term approaches the Student distribution, a prominently fat tailed distribution favoring large outliers, and leading to a greedy algorithm that strongly promotes sparsity \cite{CPrSS}.
\end{itemize}

To summarize, the above results indicate that
the hyperpriors for which the global convexity of the corresponding hypermodel is not guaranteed ($r<1$) are expected to promote sparsity more effectively than the limit case $r=1$ that can be seen as a counterpart of the $\ell_1$-penalized case.

While the hyperparameters $r$ and $\beta$ determine the strength of the sparsity promotion and the convexity properties of the MAP objective function, the vector of the scale parameters $\vartheta$ can be set automatically once the operator $\mA\mW$ is given. More specifically, for each $j$, $\vartheta_{j}$ can be related to the sensitivity of the data to $x_j$, given by the quantity  $\|\mA\mW e_j\|^2$, where $e_j\in\R^N$ denotes the canonical $j$-th Cartesian unit vector. It was proven, for $r=1$ in \cite{SIREV,CSS} and in more general settings in \cite{CPrSS}, that under the assumption that the signal-to-noise ratio is given, and that the prior satisfies an exchangeability condition guaranteeing that no particular sparse combinations of components of $x$ are favored over others,
the entries of $\vartheta$ must be chosen as
\begin{equation}
\label{eq:sens}
\vartheta_j = \frac{C}{\|\mA\mW e_j\|^2}\,,
\end{equation}
where $C>0$ is a constant encoding the expected sparsity on the solution and an estimate of the signal-to-noise-ratio. For details, we refer to the cited articles.
We remark that, in general, sensitivity weights are introduced to compensate for the possible non-uniform design of the forward model operator $\mA$.  Sensitivity weights play an important role in, e.g.,  inverse source problems, in which sources near the observation points may be favored over far away sources unless the exchangeability condition is imposed. In the current setting, when the dictionary consists of sub-frames with possibly different column norms, we expect the different weights $\vartheta_j$ to prevent the representation of the signal
the frames with larger column norms to dominate.

\section{Local and global hybrid IAS}
\label{sec:ghias}

In the following discussion, we write the penalty function ${\mathcal P}(x,\theta\mid r,\beta,\vartheta)$ in terms of components,
\begin{eqnarray*}
	{\mathcal P}(\alpha,\theta\mid r,\beta,\vartheta) &=& \sum_{j=1}^N\left(\frac 12 \frac{\alpha_j^2}{\theta_j}  -\eta\log\frac{\theta_j}{\vartheta_j} +\left(\frac{\theta_j}{\vartheta_j}\right)^r\right)\\
	&=&\sum_{j=1}^N {\mathcal P}_j(\alpha_j,\theta_j\mid r,\beta,\vartheta_j).
\end{eqnarray*}

In \cite{CPrS}, two different hybrid strategies to speed up and enhance sparsity promotion in the IAS algorithm were proposed. In both versions, the IAS iterations are initiated by selecting a conservative set of hyperparameters for which the objective function is convex, thus guaranteeing global convergence to a unique minimizer. We denote this set of parameters by $(r^{(1)},\beta^{(1)},\vartheta^{(1)})$. For the second phase of the hybrid algorithm, we select another set of parameters, $(r^{(2)},\beta^{(2)},\vartheta^{(2)})$, for which the global convexity of the objective function is not valid.
To match the models so that they express coherent prior beliefs, we adjust the scale parameters $\vartheta^{(j)}$ so as to satisfy the compatibility condition
\begin{equation}\label{compatibility}
\left(\frac{\eta^{(1)}}{r^{(1)}}\right)^{1/r^{(1)}}\vartheta_j^{(1)} = \left(\frac{\eta^{(2)}}{r^{(2)}}\right)^{1/r^{(2)}}\vartheta_j^{(2)},
\end{equation}
that guarantees that the parameter $\theta_j$ computed at $\alpha_j=0$ returns the same value regardless of the model.
For further discussion, we refer to \cite{CPrS}.

In the local hybrid version,  the IAS algorithm is initially run with hyperparameters  $(r^{(1)},\beta^{(1)},\vartheta^{(1)})$, and after each iteration step, we check which $\theta_j$, if any, satisfies the condition (\ref{th bar}), where $\overline\theta$ is computed using the hyperparameter set $(r^{(2)},\beta^{(2)},\vartheta^{(2)})$. In correspondence of those which do, we modify the local objective function so that
\begin{equation}
{\mathcal P}_j(\alpha_j,\theta_j\mid r^{(1)},\beta^{(1)},\vartheta_j^{(1)}) \rightarrow {\mathcal P}_j(\alpha_j,\theta_j\mid r^{(2)},\beta^{(2)},\vartheta_j^{(2)}).
\end{equation}

The global hybrid scheme is based on the idea that after a number of IAS iteration rounds, the iterate of the globally convex objective function with hyperparameters
$(r^{(1)},\beta^{(1)},\vartheta^{(1)})$  is near the unique global minimum of that objective function. Restarting the IAS from the current point with the parameters $(r^{(2)},\beta^{(2)},\vartheta^{(2)})$ may quickly find a local minimizer near the global minimizer of the original objective function. While the two minimizers are likely not far apart, the local minimizer is typically sparser, and the convergence to it is faster.

\section{Computed Examples}

In this section, we demonstrate the viability of the hybrid IAS algorithm in the context of overcomplete representations. More specifically,  we restrict ourselves to the global hybrid strategy, switching from the first to the second hyperprior after 10 iterations, if not differently specified. The main goal of the following examples is to demonstrate that the global hybrid IAS is capable of selecting from a dictionary of sub-frames, where several representations are admissible, a set of atoms that make the representation as sparse as possible.

\medskip

\paragraph{\textbf{Signal restoration from convolution data}}
The first test case is a one-dimensional deconvolution problem. The generative model is a piecewise constant signal  $f:[0,1]\to \R$,  $f(0)=0$, and the data consist of a few discrete observations,
\begin{equation}
b_j = \int_0^1 A(s_j - t) f(t) dt + \varepsilon_j,\quad 1\leq j\leq m,  \quad A(s_j,t) = \frac{1}{\sqrt{2\pi w^{2}}}e^{-\frac{(s_j - t)^{2}}{2w^{2}}},
\end{equation}
corrupted by Gaussian blur with $w=0.02$ and additive scaled white Gaussian noise, with standard deviation $\sigma$ set to $2\%$ of the maximum of the noiseless signal.
The data has been generated using a discretization of the unit interval with $n=n_{\rm dense} = 1253$ nodes, while in the forward model used for solving the inverse problem, we set $n=500$. The number of equidistant observation points in the signal domain is $m=46$. The generative signal and the data are shown in Figure~\ref{fig:1D_Data}.

\begin{figure}[!t]
	\centering
	\centerline{		\includegraphics[scale = 0.25]{images/over/ex1_new/true.png}\quad
		\includegraphics[scale = 0.25]{images/over/ex1_new/data.png}	}
	\caption{The generative model (left) and the blurred and noisy data vector $b \in \R^{46}$ (right).}
	\label{fig:1D_Data}
\end{figure}

The generative signal admits a natural sparse representation in terms of its increments  $z_j = z_j - z_{j-1}$ over the interval of definition. Assuming  $x_0=0$, then
\begin{equation}\label{difference}
z = \mB  x\,,\quad \mB=\left[\begin{array}{cccc}
1&0&\ldots&0\\
-1&1&\ldots&0\\
&&\ddots&\\
0&\ldots&-1&1\\
\end{array}\right]
\in \R^{n\times n },
\end{equation}
hence
\begin{equation}
x = \mL z\quad\mathrm{with}\quad  \mL = \mB^{-1} = \left[\begin{array}{cccc}
1&0&\ldots&0\\
1&1&\ldots&0\\
\vdots &&\ddots&\\
1&\ldots&1&1\\
\end{array}\right]\in\R^{n\times n}.
\end{equation}
Our goal is to test the effectiveness of the outlined framework in recovering the most natural sparse representation of the given signal.
Let $\mC$ denote the discrete cosine transform matrix, providing an alternative and accurate way of representing the signal,
\begin{equation}
x =  \mC^\mT y, \quad y = \mC x,
\end{equation}
which is, however, not sparse. To test wether the algorithm is able to identify the frame that allows a sparse representation,
we consider the overcomplete dictionary,
\begin{equation}
\label{eq:dict 1d}
\mW= [\mW_1\,,\,\mW_2]\in\R^{n \times 2n} \quad\text{with}\quad \mW_1 = \mL \in \R^{n \times n }\;\text{and}\;\mW_2 = \mC^{\mT}\in\R^{n \times n}\,,
\end{equation}
and formulate the underlying linear inverse problem as
\begin{equation}
b = \mA\mW \alpha + \varepsilon = \mA [\mW_1\,,\,\mW_2] \begin{bmatrix}
\alpha_1\\
\alpha_2
\end{bmatrix}+\varepsilon, \quad \varepsilon \sim \mathcal{N}(0,\sigma^2 \mI_{m})\,,
\end{equation}
where $\mA$ is the discrete blur operator.

In this example, the global hybrid IAS is run with parameters $(r^{(1)},\eta^{(1)})=(1,10^{-4})$, $(r^{(2)},\eta^{(2)})=(1/2,10^{-3})$; we recall that the sensitivity weights $\vartheta^{(1)}$ are set automatically according to \eqref{eq:sens}, while the vector $\vartheta^{(2)}$ is fixed so that condition \eqref{compatibility} is satisfied.

The signal reconstructed by the global hybrid IAS scheme is shown in Figure \ref{fig:rec 1d}. The restored $\alpha_1$ and $\alpha_2$ and their contribution in the estimated signal are shown in Figure \ref{fig:rec 1d 2}, together with the scaled variances corresponding to $\alpha_1$ and $\alpha_2$, i.e.
\begin{equation}
\frac{\theta_j}{\vartheta_j^{(2)}}\,,\;1 \leq j \leq n \quad\text{and}\quad \frac{\theta_j}{\vartheta_j^{(2)}}\,,\;n+1 \leq j \leq 2n\,.
\end{equation}
Notice that the output variances are scaled by the sensitivities corresponding to the second hyperprior  used to design the hybrid scheme.

\noindent Despite the relatively high level of degradation (blur and noise) and down-sampling in the observed data $b$, the algorithm has no problem detecting the basis that provides a more natural and sparse representation for the original signal. In fact, the coefficients $\alpha_2$ are five to six orders of magnitude smaller than the non-vanishing components of $\alpha_1$. The degree of sparsity in the final representation is also reflected in the number of CGLS steps per outer iteration of the global hybrid IAS - see Figure \ref{fig:rec 1d} - which quickly settles around the cardinality of the support of $\alpha$.

\begin{figure}[!t]
	\centering
	\centerline{
		\includegraphics[scale = 0.23]{images/over/ex1_new/rec_leg.png}\quad
		\includegraphics[scale = 0.23]{images/over/ex1_new/cgls.png}
	}
	\caption{Reconstruction of the signal $x$ (left) and the count of CGLS steps per outer iteration of the global hybrid IAS (right).}
	\label{fig:rec 1d}
\end{figure}

\begin{figure}[!t]
	\centering
	\centerline{\includegraphics[width = 4.3cm]{images/over/ex1_new/alpha1.png}		\includegraphics[width = 4.3cm]{images/over/ex1_new/rec1.png}		\includegraphics[width = 4.3cm]{images/over/ex1_new/var1.png}
	}
	\centerline{	\includegraphics[width=4.3cm]{images/over/ex1_new/alpha2.png}	\includegraphics[width=4.3cm]{images/over/ex1_new/rec2.png}		\includegraphics[width=4.3cm]{images/over/ex1_new/var2.png}
	}
	\caption{Vector $\alpha_i$ (left panles), corresponding scaled variances (middle panels) and contribution of the signal $\mW_i \alpha_i$ (right panels) for $i=1$, i.e. representation in terms of increments, (top row) and $i=2$, i.e. representation in terms of cosine transform, (bottom row).}

	\label{fig:rec 1d 2}
\end{figure}

\medskip

\paragraph{\textbf{Image denoising on a synthetic image}}
In the second example, we consider the problem of denoising a blocky gray scale test image $x\in\R^{n \times n}$, $n=200$. The pixel values, which are between $0$ and $1$, are corrupted by scaled white Gaussian noise with standard deviation $\sigma$ set to $10\%$ of the maximum of the noiseless image, i.e. $\sigma=0.1$ - see Figure \ref{fig:Data_2d}.

\begin{figure}[!t]
	\centering
	\centerline{
		\includegraphics[height=3.75cm]{images/over/ex2_new/rects_or.png}\;
		\includegraphics[height=3.75cm]{images/over/ex2_new/rects_cor.png}\;	\includegraphics[height=3.75cm]{images/over/ex2_new/rects_rec_2.pdf}	}
	\caption{Original image (left), observed data (middle) and reconstructed image (right).}5
	\label{fig:Data_2d}
\end{figure}

The test image presents sharp edges lying along the horizontal and vertical axes. Therefore, $x$ admits a sparse representation both in the vertical and horizontal increment bases, the latter being slightly less sparse than the former. After representing the image in vector form $x\in\R^{n^2}$ by stacking the pixel values columnwise, we introduce the redundant dictionary $\mW= [\mW_1\,,\,\mW_2]\in\R^{n^2 \times 2n^2}$ with
\begin{equation}
\label{eq:dict 2d}
\mW_1 = \left(\mI_{n} \otimes \mB\right)^{-1} \in \R^{n^2 \times n^2 }\;\text{and}\;\mW_2 = \left(\mB \otimes \mI_n\right)^{-1}\in\R^{n^2 \times n^2}\,,
\end{equation}
where $\mB$ is defined as in \eqref{difference}, and $\otimes$ stands for the Kronecker product.  Homogenous Dirichlet boundary conditions are assumed on the left and top edges of the image. We want to estimate the sparse vector $\alpha = [\alpha_1\,,\,\alpha_2]^{\mT}$, with $\alpha_i\in\R^{n^2}$, $i=1,2$, from the data vector $b\in\R^{n^2}$, given the forward model
\begin{equation}
b = \mW \alpha + \varepsilon = [\mW_1\,,\,\mW_2] \begin{bmatrix}
\alpha_1\\
\alpha_2
\end{bmatrix}+\varepsilon, \quad \varepsilon \sim \mathcal{N}(0,\sigma^2 \mI_{n^2})\,.
\end{equation}
It is worth remarking here that we require $\alpha$ to be not only sparse, but as sparse as possible.

The hyperparameters of the global hybrid IAS are set as $(r^{(1)},\eta^{(1)})=(1,10^{-3})$ and $(r^{(1)},\eta^{(1)})=(1/2,10^{-2})$, while, as before, $\vartheta^{(1)},\vartheta^{(2)}$ are automatically fixed according to \eqref{eq:sens} and \eqref{compatibility}, respectively.

The restored image is shown in Figure \ref{fig:Data_2d}, while the contribution of the vertical and horizontal increment bases together with the output scaled variances corresponding to vectors $\alpha_1$ and $\alpha_2$ are shown in Figure \ref{fig:rec 2d}. We observe that the image is almost completely restored in terms of the basis vectors corresponding to increments in the vertical direction ($\alpha_1$), whereas the entries of $\alpha_2$, corresponding to increments in the horizontal direction is negligible. The representation in terms of $\mW_1$ is indeed sparser than that in terms of $\mW_2$, due to the shorter horizontal boundary of the white inclusion compared to the vertical boundary.

\begin{figure}[!t]
	\centering
	\centerline{\includegraphics[height=3.75cm]{images/over/ex2_new/alpha1.pdf}\;		\includegraphics[height=3.75cm]{images/over/ex2_new/th_alpha1.pdf}\; 		\includegraphics[height=3.75cm]{images/over/ex2_new/walpha1.pdf} }
	\centerline{\includegraphics[height=3.75cm]{images/over/ex2_new/alpha2.pdf}\;	\includegraphics[height=3.75cm]{images/over/ex2_new/th_alpha2.pdf}\; 		\includegraphics[height=3.75cm]{images/over/ex2_new/walpha2.pdf} }
	\caption{Vector $\alpha_i$ (left panels), base-10 logarithmic plot of the corresponding scaled variances (middle panels) and vectors $\mW_i\alpha_i$ contributing to the final restoration (right panels) for $i=1$, i.e. vertical increments representation, (top row), $i=2$, i.e. horizontal increments representation, (bottom row).}
	\label{fig:rec 2d}
\end{figure}

\medskip

	\paragraph{\textbf{Image denoising on a natural image}}
	We demonstrate the scalability of our approach with large-scale denoising problem. Consider the $512 \times 512$ gray-scale natural image with pixel values between $0$ and $1$ shown in the top left panel of Figure \ref{fig:bee_data}. The observed data $b$ is a version corrupter by added white Gaussian noise with standard deviation $\sigma$ set to $5\%$ of the maximum of the noiseless image, i.e. $\sigma=0.05$,  shown in the top middle panel of Figure \ref{fig:bee_data}.

	We consider an overcomplete basis $\mathsf{W}=[\mathsf{W}_1,\,\mathsf{W}_2]\in\R^{n^2\times 2n^2}$, with $\mathsf{W}_1,\mathsf{W}_2\in\R^{n^2\times n^2}$, $n=512$, defined in \eqref{eq:dict 2d}; the vector describing the 2D signal in the selected dictionary can thus be written as $\alpha=[\alpha_1,\,\alpha_2]^\mT$, with $\alpha_1,\alpha_2\in\R^{n^2}$ representing the vertical and horizontal increments, respectively.
	Note that since this test image is not piecewise constant, but rather a mixture of jumps, smooth, and textured parts, we do not expect the vertical and horizontal increments to be \emph{naturally} sparse. In fact, when representing the original image in the chosen dictionary, about 47\% of the coefficients $\alpha_1,\alpha_2$ in both horizontal and vertical directions are non-zero.

We run the global hybrid IAS with $(r^{(1)},\eta^{(1)})=(1,10^{-4})$, $(r^{(2)},\eta^{(2)})=(1/2,10^{-3})$, and letting  $\vartheta^{(1)},\vartheta^{(2)}$ be set automatically.

The restored image is shown in the top right panel of Figure \ref{fig:bee_data}, together with a close-up to facilitate the visual comparison. In Figure \ref{fig:bee_th}, we also show the output representation vectors $\alpha_1,\alpha_2$, with the corresponding scaled output variances and their contribution in the final restoration.

\begin{figure}[!t]
	\centering
	\centerline{\includegraphics[width = 3.8cm]{images/over/bee/bee_true_rect.png}\quad \includegraphics[width = 3.8cm]{images/over/bee/bee_data2_rect.png}\quad\includegraphics[width = 3.8cm]{images/over/bee/out_it.png}}
	\centerline{\includegraphics[width = 3.8cm]{images/over/bee/bee_true_zoom.png}\quad\includegraphics[width = 3.8cm]{images/over/bee/bee_data2_zoom.png}\quad\includegraphics[width = 3.8cm]{images/over/bee/bee_out2_zoom.png}}
	\caption{\label{fig:bee_data}\emph{Top row}: original $512\times512$ test image (left), observed data (middle) and denoised image (right). \emph{Bottom row}: respective close-up(s).}
\end{figure}

\begin{figure}
	\centerline{\includegraphics[height =3.7cm]{images/over/bee/alpha_v.pdf}\;\includegraphics[height = 3.7cm]{images/over/bee/th_v.pdf}\;\includegraphics[height = 3.7cm]{images/over/bee/feat_v.pdf}}
	\centerline{\includegraphics[height =3.7cm]{images/over/bee/alpha_h.pdf}\;\includegraphics[height = 3.7cm]{images/over/bee/th_h.pdf}\;\includegraphics[height = 3.7cm]{images/over/bee/feat_h.pdf}}
	\caption{Vector $\alpha_i$ (left panels), base-10 logarithmic plot of the corresponding scaled variances (middle panels), and vectors $\mathsf{W}_i\alpha_i$ (right panels) for $i=1$, i.e. vertical increments representation, (top row) and $i=2$, horizontal increments representation, (bottom row).}
	\label{fig:bee_th}
\end{figure}

	Besides the quality of the restored image, we are interested in highlighting the compressing capability of our approach. Consider first the coefficients $\alpha=[\alpha_1,\,\alpha_2]^\mT$ of the original image, and choose a threshold value $\beta>0$ to be the smallest non-zero coefficient,
	\begin{equation}\notag
	\beta = \min\{|\alpha_j| \mid |\alpha_j|>0\}.
	\end{equation}
	After computing the restored coefficients using the sparsity promoting hybrid IAS, setting to zero those whose absolute value is below a threshold value, we find that only 3.92\% of those for the vertical components, and 4.07\% of those for the horizontal components are nonvanishing. Thus, the representation of the image in this basis was compressed by a factor more than ten from the original image without a significant deterioration in the image quality. Figure~\ref{fig:bee_hist} shows the value distributions of the original and restored coefficients, plotted as histograms in logarithmic scale: clearly the coefficient values are significantly compressed towards zero in both directions.  It is worth remarking here that the Bayesian target in sparsity promoting problems has to be understood as compressibility, as the entries of $\alpha$ cannot vanish by construction of the algorithm, but they can be made arbitrary small by suitable parameter selection.

\begin{figure}
	\centerline{\includegraphics[height =5.2cm]{images/over/bee/ValueDist_vert.eps}\;\includegraphics[height = 5.2cm]{images/over/bee/ValueDist_hor.eps}}

	\caption{Value distributions of the coefficients corresponding to vertical (left) and horizontal (right) increments in logarithmic scale. The blue distributions correspond to the original image which does not allow sparse representation in the basis, which is reflected in the high percentage of non-vanishing coefficients, while the red distribution is the sparse denoised reconstruction, in which the percentage of coefficients above the negligible threshold value is reduced by an order of magnitude.}
	\label{fig:bee_hist}
\end{figure}

\paragraph{\textbf{Image restoration}} In the fourth example, we consider the restoration problem of the $n \times n$ generative image in Figure \ref{fig:Data_2d_2}, with $n=100$,  with values in $[0,1]$. The image has been corrupted by Gaussian blur of width $w=0.006$ and additive scaled white Gaussian noise with standard deviation $\sigma$ set to $1\%$ of the maximum of the noiseless signal, i.e. $\sigma=0.01$ - see Figure \ref{fig:Data_2d_2}. The test image presents three distinctive features, namely point-wise stars, the blocky moon and the smooth cloud. After re-arranging the original $x$ in a vectorized form by stacking its entries in columnwise order, we hypothesize that a suitable dictionary for the problem of interest is $\mW = [\mW_1 \,,\, \mW_2 \,,\, \mW_3 \,,\, \mW_4] \in \R^{n^2 \times 4n^2}$, with
\begin{equation}
\mW_1 {=} \mI_{n^2}\,,\;\mW_2 {=} (\mI_n \otimes \mB)^{-1} \in \R^{n^2 \times n^2}\,,\; \mW_3 {=} (\mB \otimes \mI_n)^{-1} \in \R^{n^2 \times n^2} \text{ and } \mW_4 {=} \mC^\mT\,,
\end{equation}
where $\mB$ is as in \eqref{difference} and $\mC$ is the 2D cosine transform matrix. The problem is to estimate the sparse vector $\alpha = [\alpha_1\,,\,\alpha_2\,,\,\alpha_3\,,\,\alpha_4]^{\mT}$, with $\alpha_i\in\R^{n^2}$, for $i=1,2,3,4$, from the data vector $b\in\R^{n^2}$, given the forward model
\begin{equation}
b = \mA\mW \alpha + \varepsilon =  [\mW_1 \,,\, \mW_2 \,,\, \mW_3 \,,\, \mW_4] \begin{bmatrix}
\alpha_1\\
\alpha_2\\
\alpha_3\\
\alpha_4
\end{bmatrix}+\varepsilon, \quad \varepsilon \sim \mathcal{N}(0,\sigma^2 \mI_{n^2})\,,
\end{equation}
with $\mA$ representing the discrete blur operator.

The global hybrid IAS is run with hyperparameters $(r^{(1)},\eta^{(1)})=(1,10^{-4})$, $(r^{(2)},\eta^{(2)})=(1/2,10^{-4})$ and $\vartheta^{(1)},\vartheta^{(2)}$ automatically fixed as in the previous examples.

The image restored via the global hybrid IAS algorithm is shown in Figure \ref{fig:Data_2d_2}, while Figure \ref{fig:rec 2d 2} shows the reconstructions of the representation vectors $\alpha_i$, the corresponding variances scaled by the sensitivities, and the contribution of the vectors $\mW_i\alpha_i$ in the final restoration, for $i=1,2,3,4.$
\begin{figure}[!t]
	\centering
	\centerline{
		\includegraphics[height=3.75cm]{images/over/ex3_new/sky_or.pdf}\;
		\includegraphics[height=3.75cm]{images/over/ex3_new/sky_cor.pdf}\;	\includegraphics[height=3.75cm]{images/over/ex3_new/res_sky.pdf}	}
	\caption{Original image (left), observed data (middle) and reconstructed image (right).}
	\label{fig:Data_2d_2}
\end{figure}

We point out that, as in the previous example, the representation vectors in both the vertical and horizontal increment bases are sparse. Nonetheless, the hybrid hypermodel selects the one with fewer non-zero entries.

\begin{figure}[!t]
	\centering
	\centerline{\includegraphics[height=3.75cm]{images/over/ex3_new/a_1.pdf}\;	\includegraphics[height=3.75cm]{images/over/ex3_new/th1.pdf}\; 		\includegraphics[height=3.75cm]{images/over/ex3_new/wa_1.pdf} }
	\centerline{\includegraphics[height=3.75cm]{images/over/ex3_new/a_2.pdf}\;		\includegraphics[height=3.75cm]{images/over/ex3_new/th2.pdf}\; 		\includegraphics[height=3.75cm]{images/over/ex3_new/wa_2.pdf} }
	\centerline{\includegraphics[height=3.75cm]{images/over/ex3_new/a_3.pdf}\;		\includegraphics[height=3.75cm]{images/over/ex3_new/th3.pdf}\; 		\includegraphics[height=3.75cm]{images/over/ex3_new/wa_3.pdf} }
	\centerline{\includegraphics[height=3.75cm]{images/over/ex3_new/a_4.pdf}\;		\includegraphics[height=3.75cm]{images/over/ex3_new/th4.pdf}\; 		\includegraphics[height=3.75cm]{images/over/ex3_new/wa_4.pdf} }
	\caption{Representation vectors $\alpha_i$ (left panels), base-10 logarithmic plot of the corresponding scaled variances (middle panels) and vectors $\mW_i\alpha_i$ contributing to the final restoration (right panels) for $i=1$ (first row), $i=2$ (second row), $i=3$ (third row) and $i=4$ (fourth row). }
	\label{fig:rec 2d 2}
\end{figure}

\medskip
\paragraph{\textbf{Dictionary learning}} The final example, coming from machine learning, is concerned with the sparse identification of hand-written digits based on a dictionary of annotated data.
Consider the MNIST data set of hand-written digits $0,1,\ldots,9$ digitized as $16\times 16$ black-and-white images. Denoting by $w^{(j)}\in\R^{(16)^{2}}$, $1\leq j\leq N$ the vectorized image vectors of $N = 1\,707$ handwritten digits constituting the atoms of the dictionary, and by $c_j\in\{0,1,\ldots,9\}$ the corresponding annotations, we form the dictionary matrix
\begin{equation}
\mW = \left[\begin{array}{ccc} w^{(1)} &\cdots & w^{(N)}\end{array}\right] \in \R^{n\times N}, \quad n=256,\; N= 1\,705.
\end{equation}
To identify an handwritten digit $b$ drawn from an independent set of handwritten digits, we seek to represent it in a sparse manner in terms of the given dictionary,
\begin{equation}
b = \mW \alpha +\varepsilon,
\end{equation}
where $\alpha\in\R^N$ is a sparse vector, and $\varepsilon$ represents the discrepancy between the data and its representation. The idea is represented schematically in Figure~\ref{fig:DLcartoon}.
We point out that in the dictionary consisting of all handwritten digits, the digits with same annotation can each be thought of representing a sub-dictionary, and as the proposed algorithm seeks the most economic representation, it is natural that the representation corresponds to picking the representing atoms from the sub-dictionary with greatest affinity with the  digit that represents the data.
\begin{figure}
	\centerline{\includegraphics[width=10cm]{images/DictionaryLearningCartoon.png}}
	\caption{\label{fig:DLcartoon} A schematic representation of the dictionary learning example. The digit on the right is the non-annotated image $b$, which is approximated in terms of the annotated atoms $w_i$ on the right. The coefficients $\alpha_i$ can then be used to identify the digit.}
\end{figure}

In this example, we run the global hybrid IAS algorithm using the parameters $(r^{(1)},\beta^{(1)},\vartheta^{(1)}) = (1,3/2+10^{-4},10^{-5})$, where all components of the vector $\vartheta^{(1)}$ are assumed equal, as sensitivity is not an issue in this example, and  $(r^{(2)},\beta^{(2)}) = (-1,1)$, with the value $\vartheta^{(2)}$ determined from the compatibility condition (\ref{compatibility}). Furthermore, since the digit images are non-negative, after each update step of the pair $(\alpha,\theta)$, we project the image to the positive cone. A theoretical justification of the projection step was given in \cite{CPrSS}. We switch from the first to the second model in the hybrid IAS scheme when either the relative change in $\theta$ with respect to the $\ell_2$-norm falls below $10^{-3}$ or 80 iterations have been completed.

Figures~\ref{fig:DL1}, \ref{fig:DL2} and \ref{fig:DL3} show the results with different choices of the standard deviation of the likelihood. Observe that here, the noise term $\varepsilon$ represents the discrepancy between the data $b$ and its representation in terms of the dictionary, and can be chosen according to how much fidelity is required. Choosing $\sigma$ large allows a very sparse representation, as the required quality of the approximation is low, however, poor approximation easily leads to a mis-labeling of the digit. On the other hand, decreasing $\sigma$ forces the approximation to be better, and more atoms are required. The labeling can be done using the majority vote principle. In the computed example, Figures~\ref{fig:DL1} and \ref{fig:DL2}, the labeling with majority vote is correct in each case, while in Figure~\ref{fig:DL3} with sparser representation, mislabelings occur.

\begin{figure}
	\centerline{\includegraphics[width=13cm]{images/DictionaryLearningResultHist_1.png}}
	\caption{\label{fig:DL1} Dictionary learning results. The first row shows the test images of the digits to be classified by the dictionary learning algorithm (vector $b$), the true annotation indicated in the figure, the second row the vectors $\theta$ after the IAS iteration with the first hyperprior, and the third row after the iteration with the second hyperprior. The fourth row represents the synthesis $\mW\alpha$ approximating the original digit, and finally, the fifth row gives the histogram of the annotations of the atoms corresponding to coefficients above a threshold $\tau = 0.01$. The annotation is done by majority vote, choosing the largest of the bins. In this example, the standard deviation of the noise representing the mismatch was $\sigma = 0.01$.}
\end{figure}

\begin{figure}
	\centerline{
		\includegraphics[width=13cm]{images/DictionaryLearningResultHist_2.png}
	}
	\caption{\label{fig:DL2} The rows are as in Figure~\ref{fig:DL1}. In this example, the standard deviation of the noise representing the mismatch was $\sigma = 0.05$. Observe that the approximation becomes sparser.}
\end{figure}

\begin{figure}
	\centerline{
		\includegraphics[width=13cm]{images/DictionaryLearningResultHist_3.png}
	}
	\caption{\label{fig:DL3} The rows are as in Figure~\ref{fig:DL1}. In this example, the standard deviation of the noise representing the mismatch was $\sigma = 0.1$. The increased sparsity here is traded with an increased number of misclassifications, such as in the first and the third columns.}
\end{figure}

\section{Conclusions and future work}
\label{sec:concl}
The hierarchical Bayesian framework combined with Krylov subspace iterative solvers for large linear systems is well suited for the design of computationally efficient methods to solve large scale ill posed inverse problems with sparsity constraints. From the point of view of computational efficiency, replacing a whitened Tikhonov-type penalty by a Krylov subspace iteration equipped with early stopping is of crucial importance. The estimate found by this process is not guaranteed to coincide with the MAP estimate, and it was referred to as quasi-MAP estimate in \cite{CPPSV}. The statistical analysis of the early stopping regularization strategy based on Krylov subspace methods is not straightforward, since the estimate depends non-linearly on the data, and it can be seen as an approximate Bayesian computing (ABC) strategy. For further discussion, see, e.g. \cite{SIREV}.
Here we have shown that the framework can be naturally adapted for dealing with overcomplete systems, consisting of,  e.g., combined frames or bases. The approach has significant potential when it may not be known a priori which frame is best suited for representing the unknown, leaving it up to the algorithm to find the most parsimonious representation. In order to avoid that one frame is favored over another, however, it is important that the data are equally sensitive to components in every frame. Fortunately, the sensitivity analysis developed by the authors in \cite{calvetti2019brain,CSS,CPrSS}, provides naturally such scaling. The proposed sensitivity weights are rooted in the very natural Bayesian principle of exchangeability, stating that no set of non-zero components with a given cardinality should be favored over any other. In light of this principle, the scaling guarantees the same explanation power for every sub-frame, so the one leading to most sparse solution is automatically selected. This feature may turn out particularly useful in machine learning, with applications such as MRI fingerprinting (see,e.g., \cite{ma2013magnetic}).
In \cite{CSS}, a connection between the proposed IAS algorithm and the compressed sensing literature  \cite{candes1} was considered, suggesting that when the forward model guarantees perfect sparse recovery, the IAS algorithm effectively finds a good approximation of it. It is reasonable to believe that the results can be extended to overcomplete dictionaries, for which similar recovery results are known \cite{drip}.

The methodology developed in this paper has been tested only with Gaussian noise, leading to a quadratic fidelity term in the optimization problem. The IAS framework has been shown to work well with other noise models, e.g., Poisson distributed noise in connection with low dose X-ray tomography and PET, see \cite{BCS}. The applicability of the approach to non Gaussian noise may be very important for its use for dictionary learning problems where the data consist, e.g., of word counts. The extension of the method to large scale problems, different noise models and nonlinear forward models is the next step and will be addressed in separate future contributions.

config: arxiv


TITLE: example of convexity
QUESTION [3 upvotes]: Let $\mathcal{B}(F)$ the algebra of all bounded linear operators on a complex Hilbert space $F$. Let $(S_1,S_2)\in \mathcal{B}(F)^2$. We define
$$W(S_1,S_2)=\{(\langle S_1 y\; ,\;y\rangle,\langle S_2 y ,\;y\rangle):y \in F,\;\;\|y\|=1\}.$$
Consider the matrices
$$
S_1 = \left[\begin{array}{ccc} 0&1&0 \\ 0& 0&0 \\0&0&0 \end{array}\right] \ \ \textrm{and} \ \ S_2 = \left[\begin{array}{ccc} 0&0&0 \\ 0& 0&0 \\0&1&0 \end{array}\right].
$$
Clearly, $S_1S_2=S_2S_1$. Moreover, we get
$$
W(S_1,S_2)=\{(b\overline{a},b\overline{c});\;(a,b,c) \in \mathbb{C}^3\;\;\hbox{and}\;|a|^2+|b|^2+|c|^2=1\}.
$$

Is $W(S_1, S_2)$ a convex subset of $\mathbb{C}^2$ in this case?

REPLY [3 votes]: In your specific case, $W(S_1,S_2)$ is the polydisk of radius $\frac12$. For
$$|b\bar a|^2+|b\bar c|^2=|b|^2(1-|b|^2)\le\frac14$$
and conversely if $|w|^2+|z|^2\le\frac14$, one takes $b=\sqrt t\in(0,1]$ a root of $t(1-t)=|w|^2+|z|^2$, then $\bar a=w/t$, $\bar c=z/t$.


TITLE: Sum of $\sum\limits_{i=1}^\infty (-1)^i\frac{x^{2i+1}}{2i+1}$
QUESTION [1 upvotes]: Can someone help me with this series? It was on my exam and I don't know how to do it.
For $|x| < 1$ determine the sum of
$$\sum\limits_{i=1}^\infty  (-1)^i\frac{x^{2i+1}}{2i+1}$$

REPLY [5 votes]: The derivative of the given sum is the geometric sum
$$\sum_{n=1}^\infty  (-x^2)^n=-\frac{x^2}{1+x^2}=-1+\frac{1}{1+x^2}$$
so the given sum which vanish at $0$ is
$$\int_0^x\left(-1+\frac{1}{1+t^2}\right)dt=-x+\arctan x$$

REPLY [3 votes]: Hint: Expand $\frac{1}{1+t^2}$ using the familiar $\frac{1}{1-y}=1+y+y^2+y^3+\cdots$, and integrate term by term from $0$ to $x$.  Then adjust for the fact your sum starts at $i=1$.


TITLE: Behaviour of inductor (understanding equation)
QUESTION [0 upvotes]: If an inductor is connected to A.C. source (thinking no resistance is there), and if we travel across circuit there is voltage drop, and it is said
$$V=L\frac{\text dI}{\text dt}$$
How we write this equation? And if the applied voltage equals back emf, how does current flow?

REPLY [0 votes]: The back EMF is created by the current flow, via the associated magnetic field. It rises until it equals the applied voltage. The difference is then zero and the current can no longer increase. But it is still flowing.
The relationship between voltage and current for a given coil is called its impedance $Z$. It is strongly frequency-dependant, but that is another story.


TITLE: Can law of conservation of energy be violated?
QUESTION [1 upvotes]: Consider a positive charge of charge q placed in between plates which are charge with charges opposite in sign and equal in magnitude.A uniform magnetic field is present downwards in the right half of the midway between charged plates then charge slowly accelerates and after reaching the right half,the charge deflects and gets out of the electric field produced due to charged plates but magnetic field is present even outside the plates then charge further deflects and comes in between the plates and it starts accelerating again.this process repeats and the velocity keeps on increasing. Here work done by magnetic field is zero.Charge in the plates remain constant ensuring constant electric field between plates.The ball is gaining kinetic energy continuously.Where does this energy come from?The charged plates are kept between rigid walls to stop movement of plates due to force.The image of this situation is given below.

REPLY [1 votes]: The magnetic field may do no work on the test charge, but the electric field does.
I assume you understand what's happening if you only look in between the plates (because there the electric field is constant).  If the charge gets closer to the negative plate, it gains (kinetic) energy.  If the charge moves further from the negative plate, it loses (kinetic) energy.
Now, when you leave the area between the plates, the specifics of the field are messier, but it's not sharp boundary between a full-strength field inside and zero field outside.  The field becomes weaker, but is still present.  If the particle can gain $E$ energy by moving between two points between the plates, it will require $E$ energy to move backward even if you pick a path that is outside the plates.   The field there might be weaker, but your path will be longer.  Together, it will sum to the same value since the electric field is conservative.

REPLY [0 votes]: You're not correct saying that the charge is constantly gaining kinetic energy. Any motion towards the left reduces its kinetic energy because it gains potential energy. Further, if the charge doesn't come in exactly perpendicular to the left plate, then it will be deflected to the right plate at an angle this time, and eventually the charge will collide with the right plate.


TITLE: A finite alternating sum
QUESTION [18 upvotes]: We have stumbled upon the following finite alternating sum, which we have trouble analyzing. The sum is:
$$
S_n = \sum_{j=0}^n \frac{ (-1)^j e^{-j} }{j!} (n-j)^j
$$
We have observed numerically that $S_n \approx 2 n e^{-n}$. We would like to establish whether this conjecture is true. More precisely, we would like to show that $S_n= \Theta(n e^{-n})$.
The sum is quite hard to evaluate numerically because the summands grow in absolute value initially and then decrease toward zero. The largest term is exponential in $n$ while the sum is conjectured to converge to zero exponentially fast. Any references or ideas are welcome!

REPLY [18 votes]: WhatsUp's generating function allows to get an asymptotics: we have $$\sum S_n t^n=\frac{t}{e^{t/e}-t},\sum S_n e^nt^n=\frac{t}{e^{t-1}-t}.$$
Denote $t-1=x$, then $$\frac{t}{e^{t-1}-t}=\frac{1+x}{e^x-1-x}=\frac{1+x}{x^2/2+x^3/6+\ldots}=\frac{2(1+x)}{x^2(1+x/3+\ldots)}=\frac2{x^2}+\frac{4/3}{x}+O(1)$$
for small $x$. Thus $$\frac{t}{e^{t-1}-t}=\frac2{(1-t)^2}-\frac{4/3}{1-t}+F(t)=\sum \left(2n+\frac23\right)t^n+F(t),$$
where $F$ does not have pole at 1. I claim that there is also no pole with absolute value at most 1, that is, the function $e^{t-1}-t$ does not have zeroes in the closed unit disc other than 1. Indeed, assume that $0\leqslant |r|\leqslant 1$ and $\varphi\in (-\pi,\pi]$ and $e^{re^{i\varphi}-1}=re^{i\varphi}$. Taking the arguments of both sides we get $r\sin \varphi=\varphi$ that yields $\varphi=0$, otherwise $|\varphi|>|\sin \varphi|\geqslant r|\sin \varphi|$. So $r=e^{r-1}\geqslant 1+(r-1)=r$ with equality only for $r=1$. This all implies that $F$ is analytic in a disc with radius greater than 1, therefore the coefficients of $F$ are $O(a^n)$ for certain $a\in (0,1)$ and in your question we get $$S_n=(2n+2/3)e^{-n}(1+O(a^n)).$$


TITLE: Proof that a convergent sequence has a monotones subsequence?
QUESTION [2 upvotes]: I want to prove that every convergent sequence has a monotone subsequence. My idea is to trap some terms of the convergent sequence in successively smaller intervals around the limit.
Let $\lim_{n\to\infty} a_n=A$, and suppose that $(a_n)$ is not eventually constant. There must be infinitely terms of $(a_n)$ in at least one of the intervals $[A-1, A)$, $(A, A+1]$.
Suppose there are infinitely many terms in $(A, A+1]$. Let $I_j$ be the interval $(A+2^{-j}, A+2^{1-j}]$, $j \in \mathbb N$. Then each interval is non-overlapping, and $$\bigcup_{j=1}^\infty I_j = (A, A+1].$$
Let $I_{j_k}$ be the ordered, non-empty intervals (intervals containing terms of $(a_n)$). Let $a_{n_k}$ be any one term contained in $I_{j_k}$. Then $a_{n_k}\gt a_{n_{k+1}}$, so that $(a_{n_k})$ is a monotone decreasing subsequence of $(a_n)$.
I'm sure there's a much simpler, more elegant way to do this. But does my way work? Is there anything I'm overlooking?

REPLY [1 votes]: Suppose there are infinitely many terms in $(A,\infty).$ Choose $n_1$ so that $a_{n_1}\gt A.$ Then choose $n_2\gt n_1$ so that $a_{n_1}\gt a_{n_2}\gt A.$ Then choose $n_3\gt n_2$ so that $a_{n_2}\gt a_{n_3}\gt A.$ And so on.


text:
\begin{document}
\title{Extraordinary numbers and the Riemann hypothesis}
\author{Tatenda Kubalalika}
\email{tatendakubalalika@yahoo.com}
\maketitle
\textbf{ABSTRACT.} In this note, we demonstrate the Riemann Hypothesis via Robin's criterion concerning the magnitude of $\sigma(n):=\sum_{d\mid n} d$. Our approach is particularly concerned with extraordinary $n$.
\\
\\
\textit{2020 MSC Subject classifications: 11M26; 11AXX.}\\
\textit{Keywords and phrases: Extraordinary numbers; Riemann Hypothesis; Robin's criterion.}
\\
\\
\textbf{Introduction.} \textit{Throughout this paper,} we denote by $p_j$ the $j$-th prime. We also: let $N_j =p_{1}p_{2}\cdots p_j$, denote by $\gamma=0.57721\cdots$ the Euler-Mascheroni constant and define, for all $n \in \mathbb{N}_{\geq 2}$, $G(n):=\frac{\sigma(n)}{n\log \log n}$ where $\sigma(n):=\sum_{d\mid n} d$. It is a classical result of Gronwall [7] that $\limsup_{n\rightarrow \infty} G(n)=e^{\gamma}$. Therefore, the choice of the letter $G$ is in honour of Gronwall. Due to Robin [1, Theorem 1], we know that the Riemann Hypothesis (RH) is equivalent to the statement that $G(n) \leq e^{\gamma}$ for every integer $n>5040$. A similar criterion was established by Nicolas [3], who showed that the RH holds if and only if $H_{j}:\frac{N_j}{\varphi(N_j)\log \log N_j} \geq  e^{\gamma}$ for every integer $j\geq 2$, where $\varphi$ denotes the Euler totient function. For a far more extensive discussion of $G, H$ and other related functions, the interested reader is referred to the cited literature.
\begin{center}
\section*{Main results}
\end{center}
\textbf{Definition 1.} An integer $n>1$ is said to be a GA1 number if it is composite and $G(n)\geq G(n/p)$ for every prime $p\mid n$. We say $n$ is a GA2 number if $G(n)\geq G(cn)$ for all $c \in \mathbb{N}$. Finally, $n$ is said to be \textit{extraordinary} if it is both GA1 and GA2.
\\
\\
\textbf{Definition 2.} A \textit{superabundant (SA)} number is a positive integer $s$ for which $\frac{\sigma(s)}{s} \geq \frac{\sigma(t)}{t}$ for every positive integer $t\leq s$, with equality only for $s=t$.
\\
\\
\textbf{Lemma 1 (Theorem 1 of [5]).} \textit{Let $v>36$ be an SA number. Then $v$ must have a prime factorisation of the form $\prod_{j=1}^{k} {p_j}^{a_j}$ where $a_{k}=1$ and $a_{j} \geq a_{j+1}\in \mathbb{N}$ for every positive integer $j \leq k-1$.}
\\
\\
\textbf{Lemma 2 (Theorems 9 and 12 of [2]).} \textit{Suppose $n>500$ is an extraordinary number whose largest prime factor is $P(n)=p$. Then $p<\log n$ and $\gcd(\frac{n}{p}, p)=1$.}
\\
\\
\textbf{Remark on Lemma 2.} Note that Theorems 9 and 12 of [2] were specifically proved for GA1 numbers with at least three prime factors (counting multiplicity). However, we claim that any extraordinary number $n>500$ must satisfy this condition. Indeed, suppose $n>500$ is an extraordinary number with at most two prime factors (counting multiplicity). Since $n$ must be even (see e.g. p.4 of the arxiv version of [2] or p.362 of its journal version), put $n=2w$ where $w>250$ is a prime. Thus $\sigma(n) = 3(w+1)$ hence $G(n)<e^{\gamma}=1.78107\cdots$. But this contradicts Fact 2 of [2] that $G(\alpha)\geq e^{\gamma}$ for every GA2 number $\alpha$, hence our supposition must be false so our claim holds.
\\
\\
\textbf{Lemma 3 (see [1, p.203]).} \textit{Let $v=\prod_{j=1}^{k} {p_j}^{a_j}>6$ be an SA number and $q_v$ be the largest prime for which ${q_v}^2 \mid v$. If $p_{k}>20000$, then $\log N_k < \log v - 0.998q_v$. It is essential to have $v>6$ so as to be sure that the SA number $v$ is divisible by a square of some prime.}
\\
\\
\textbf{Lemma 4 (Lemmas (b) and (c) of [1, p.190]).} \textit{Let $v=\prod_{j=1}^{k} {p_j}^{a_j}>36$ be an SA number and $q_v$ be as defined in Lemma 2. If $p_{k} \geq 1530$, then $q_v \in (\sqrt{p_k}, \sqrt{2p_k})$.}
\\
\\
\textbf{Lemma 5 (see [1, p.206]).} \textit{Let $v=\prod_{j=1}^{k} {p_j}^{a_j}$ be an SA number. If $p_{k} > 20000$, then $v \in [N_k, N_{k+1})$.}
\\
\\
\textbf{Remark on Lemmas 3, 4 and 5.} Note that in [1, p.203], Robin used $N,x$ and $x_2$ in place of $v, p_{k}$ and $q_v$ respectively. He also denoted $\log N_k$ by $\theta(x):=\sum_{p\leq x} \log p$ and he considered \textit{colossally abundant (CA)} numbers \textit{throughout} his paper, instead of SA numbers. Therefore, to be specific, our Lemmas 3, 4 and 5 hold for every CA number $v$ whose largest prime factor is $>20000$. Since every CA number is also an SA number, it follows that our Lemmas 3, 4 and 5 also hold for every SA number $v$ whose largest prime factor is $>20000$. Note that a CA number is a positive integer $U$ for which there exists some fixed $\varepsilon>0$ such that $\frac{\sigma(U)}{U} \geq \frac{\sigma(V)}{V}(U/V)^{\varepsilon}$ for every positive integer $V$.
\\
\\
\textbf{Lemma 6 (Theorem 4 of [2]).} \textit{If the Riemann Hypothesis is false, then there exists some extraordinary number $\nu>5040$ such that $G(\nu)\geq e^{\gamma}$ and $G(\nu) \geq G(\xi)$ for every integer $\xi > 5040$. }
\\
\\
\textbf{Theorem 1.} \textit{The Riemann Hypothesis holds.}
\\
\begin{proof} Suppose that the RH is false. Then by Lemma 6, there should exist some extraordinary number $n>5040$ such that $G(n)\geq e^{\gamma}$ and $G(n)\geq G(y)$ for every integer $y \in (5040, n)$, thus
\begin{equation}
\frac{\sigma(n)}{n} \geq \frac{\sigma(y)\log \log n}{y\log \log y} > \frac{\sigma(y)}{y}
\end{equation}  for every integer $y \in (5040, n)$. It is known [4] that $G(t)<e^{\gamma}$ for all $t \in (5040, 10^{10^{10}}]$ hence we must have $n>10^{10^{10}}$. Thus it follows from $(1)$ that $n$ must be SA since $10080 \in (5040, n)$ is SA. Let $P(n)=p$ be the largest prime factor of $n$. By Lemma 2, we can put $n=mp$ where $\gcd(m, p)=1$ hence $\sigma(n)=\sigma(mp)=\sigma(m)\sigma(p)$. Since $n=mp$ is a GA2 number, note that $G(mp)\geq G(mp^2)$. That is,
$\frac{\sigma(m)\sigma(p)}{mp \log \log mp} \geq \frac{\sigma(m)\sigma(p^2)}{mp^{2} \log \log mp^{2}}$ hence
\begin{equation}
1 - \frac{1}{p^{2} + p+1}=\frac{p\sigma(p)}{\sigma(p^2)} \geq \frac{\log \log mp}{\log \log mp^2}=1+\frac{\log(1-\frac{\log p}{\log (mp^2)})}{\log \log mp^2}.
\end{equation}
By Lemma 2, we have $p<\log n$ thus $\frac{\log p}{\log mp^2}<\frac{\log \log n}{\log n}<0.01$ since $n>10^{10^{10}}$. This entails that $\log \Big(1-\frac{\log p}{\log mp^2}\Big)>-\frac{2\log p}{\log mp^2}$. Plugging this into $(2)$ gives
\begin{equation}
1 - \frac{1}{p^{2} + p+1} > 1 -\frac{2\log p}{(\log mp^2)(\log \log mp^2)}
\end{equation} thus
\begin{equation}
2(p^2 + p + 1)\log p > (\log mp^2)(\log \log mp^2).
\end{equation}
 From $p < \log n = \log mp < \log mp^2$ we obtain $\log p < \log \log mp^2$. Applying this into $(4)$ yields $p^2 + p + 1 > \frac{1}{2}\log mp^2$ hence
\begin{equation}
p > \frac{1}{3}\sqrt{\log mp^2}>\frac{1}{3}\sqrt{\log n}.
\end{equation} Recall that $n>10^{10^{10}}$. Putting this into $(5)$ gives
\begin{equation}
p > 20000.
\end{equation}
 Since $n$ is SA, we know by Lemma 1 that it must have a prime factorisation of the form $2^{a_1}3^{a_2}\cdots p_r$ where $p_r = p$. Combining $(6)$ with Lemma 5 yields $n \in [N_r, N_{r+1})$ thus \begin{equation}
p<\log n < \log N_{r+1} = \log p_{r+1} + \log N_r < \log p_{r+1} + \log n - 0.998q_n
\end{equation} by Lemma 3. Applying Bertrand's Postulate into $(7)$ yields
 \begin{equation}
\log (2p_r)>\log p_{r+1} > 0.998q_n.
\end{equation} But $(8)$ contradicts Lemma 4 by which we know that $q_n > \sqrt{p_r}$ where $p_{r}=p>20000$. Hence we deduce that our supposition must be false, so we are done.
\end{proof}
\textit{Dedicated to my dear parents and all-weather friends. Unfortunately, my father recently departed (25/11/1933-25/01/2021).}


TITLE: Introducing a mathematical structure.
QUESTION [6 upvotes]: Most of the times, introducing a new mathematical structure is done in the following path.

Start with a set/collection, name it as $X$. It is possible that $X$ already have a structure with it, namely the structure of topological  space/manifold/vector space etc.
Define a structure on $X$, denote it by $\mathcal{A}$. So, a structure is a pair $(X,\mathcal{A})$.
Define what does it mean to say a substructure of $(X,\mathcal{A})$. Giving names to well-behaved substructures.
Define what are maps between two structures, say $(X,\mathcal{A})$ and $(Y,\mathcal{B})$. Giving names to well-behaved maps between two structures.
Define what does it mean to say two structures $(X,\mathcal{A})$ and $(Y,\mathcal{B})$ are "equivalent".
Construct new structures from old structures.

This step 6 differs drastically from one structure to another structure. Some ways to produce new structures from old structures are

Quotients.
Pullbacks.
Products (direct).
Sums (direct).
Limits (Injective/Projective).
...
...

For example, defining the notion of structure of a group on a set, we follow the same procedure. Defining a group, defining what does it mean to say a group morphism, what it means to say a subgroup (a normal subgroup), quotients of subgroup, (direct) sum, (direct) product of groups, (Injective/Projective) limit of (a collection) of groups and so on.
In this post, I want to collect this procedure for most of the structures introduced in undergraduate or beginning graduate courses in Mathematics.

REPLY [3 votes]: Definition. Let $G$ be a set and let $(\cdot):G \times G \rightarrow G$ be a binary operation, called multiplication. Then $(G, \cdot)$ is a group if

Asocciativity: $(g_1g_2)g_3 = g_1(g_2g_3)$ for all $g_1, g_2, g_3 \in G$.
Identity: There exists an element $e$ such that $eg = ge = g$ for all $g \in G$.
Inverses: For any $g \in G$ there exists an element $g^{-1} \in G$ such that $gg^{-1} = g^{-1}g = e$.

Substructures. A subgroup of a group $G$ is a subset of $G$ that is itself a group and the inclusion is a group homomorphism (see below).
Maps. A group homomorphism is a map $f: G \rightarrow H$ between groups such that for all $g, h \in G$ we have $f(gh) = f(g)f(h)$. A group isomorphism is a bijective group homomorphism.
Constructions. The product group of a family $\{G_i\}_{i \in I}$ of groups is obtained as follows. The underlying set is the product $\prod_{i \in I} G_i$ and the group operation is coordinatewise multiplication.
The direct sum of a family $\{G_i\}_{i \in I}$ of groups is very similar in construction to the product, but differs in one important detail: The underlying set is the subgroup $\bigoplus_{i \in I}G_i$ of $\prod_{i \in I}G_i$ such that all but finitely many components are equal to the identity.
A normal subgroup of a group is a subset $N$ of $G$ such that $gNg^{-1} = N$ for all $g \in G$. In this case, the quotient group $G/N = \{gN : g \in G\}$ has multiplication $gN \cdot hN = ghN$. The requirement that $N$ is normal is equivalent to the requirement that the multiplication on $G/N$ is well-defined.


TITLE: The relation between quasi-permutation matrix and permutation matrix?
QUESTION [3 upvotes]: We know that  a quasi-permutation matrix is a square matrix over the
complex numbers with non-negative integral trace.
Can anyone tell me why it is called "quasi-permutation matrix"? Is there any relation between quasi-permutation matrix and permutation matrix?

REPLY [2 votes]: Perhaps the best way to appreciate the choice of terminology is by looking back at the 1962 paper Linear groups analogous to permutation groups by W.J. Wong, where the concept is introduced:

If $G$ is a finite linear group of degree $n$, that is, a finite group of automorphisms of an $n$-dimensional complex vector space (or, equivalently,
  a finite group of non-singular matrices of order $n$ with complex coefficients),
  I shall say that $G$ is a quasi-permutation group if the trace of every element of $G$ is a non-negative rational integer. The reason for this terminology is that, if $G$ is a permutation group of degree $n$, its elements, considered as acting on the elements of a basis of an $n$-dimensional complex vector space $V$, induce automorphisms of $V$ forming a group isomorphic to $G$. The trace of the automorphism corresponding to an element $x$ of $G$ is equal to the number of letters left fixed by $x$, and so is a non-negative integer. Thus, a permutation group of degree $n$ has a representation as a quasi-permutation group of degree $n$.

Although stated slightly differently than I did in my Comment, it provides the same essential justification for the terminology, anticipating that this more general framework for finite group representations (linear = using matrices) may share some nice properties with permutation groups.
Wong next states two properties of a permutation group $G$ of degree $n$:

The order $|G|$ is a divisor of $n!$
If $p$ is a prime number exceeding $\sqrt{n}$, then the Sylow $p$-subgroup of $G$ is of elementary Abelian type.

The rest of the paper basically proves analogous statements about quasi-permutation groups:

Theorem $1.\;\;$  The order of a quasi-permutation group $G$ of degree $n$ is a divisor of $n!$
Theorem $2.\;\;$  The Sylow $p$-group of a quasi-permutation group of degree $n$ smaller than $p^2$ is of elementary Abelian type.

The closing paragraph of that paper summarizes the two results as showing that "any $p$-group which can be represented as a quasi-permutation group of degree $n$ can also be represented as a permutation group of degree n, provided that $p^2 \gt n$."  The paragraph then notes that the condition $p^2 \gt n$ cannot be omitted by giving the example of the quaternion group of order $8$, which has a representation as a quasi-permutation group of degree $4$, but none as a permutation group of degree $4$.
One might read this example as suggesting that generalizing to quasi-permutation groups allows for representations of smaller degree than permutation groups do, while still providing concrete opportunities for investigation.  Much of the subsequent research literature has compared the minimal degree $p(G)$ of a (faithful) representation of certain finite groups $G$ by permutation groups with the minimal degree $c(G)$ by quasi-permutation groups.  The introductory fact that permutation groups are quasi-permutation groups implies $c(G) \le p(G)$.
Finally we note that quasi-permutation group has to be understood in a slightly more permissive (artful?) spirit than permutation group is.  The product of two compatible sized permutation matrices is again a permutation matrix, but the product of two quasi-permutation matrices need not be a quasi-permutation matrix (though it will be a posteriori if both are elements of a quasi-permutation group).  For example:
$$ \begin{pmatrix} 0 & i \\ i & 0 \end{pmatrix}
   \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}
 = \begin{pmatrix} i & 0 \\ 0 & i \end{pmatrix} $$
This last has trace $2i$, not an non-negative integer, so it is not a quasi-permutation matrix despite being the product of two quasi-permutation matrices.


TITLE: Diagonalization with orthogonal matrix?
QUESTION [1 upvotes]: Determine an orthogonal amtrix $P$ such that $P^{-1}AP$ is a diagonal matrix, where
$$A= \begin{bmatrix} 7 & 4 & -4 \\ 4 & -8 & -1 \\ -4 & -1 & -8\end{bmatrix}$$
I found the required matrix by using eigen vector approach. But it is not orthogonal.
How to proceed to find orthogonal matrix. Please help

REPLY [1 votes]: Guide:
Suppose you have found two eigenvectors but they are not orthogonal to each other. $\{ u_1, u_2\}$.
$$Au_1 = \lambda u_1$$
$$Au_2 = \lambda u_2$$
Let $v_1 = u_1$ and $v_2 = u_2 - \frac{v_1.u_2}{\|v_1\|^2}v_1$, then $v_2$ and $v_1$ are orthogonal to each other.
To make them orthonormal, divide by their length.
You might want to check out Gram-Schmidt process.


TITLE: Two definitions of the half-derivative of $e^x$ don't match. Who's right?
QUESTION [3 upvotes]: I find two main sources on how to compute the half-derivative of $e^x$. Both make sense to me, but they give different answers.
Firstly, people argue, that
$$\begin{align}
\frac{\mathrm{d}}{\mathrm{d} x} e^{k x} &= k e^{k x} \\[4pt]
\frac{\mathrm{d}^2}{\mathrm{d} x^2} e^{k x} &= k^2 e^{k x} \\[4pt]
\frac{\mathrm{d}^n}{\mathrm{d} x^n} e^{k x} &= k^n e^{k x}
\end{align}$$
Therefore, it seems very reasonable, that
$$\frac{\mathrm{d}^{1/2}}{\mathrm{d} x^{1/2}} e^{k x} = \sqrt{k} e^{k x}$$
But this is not what the usual formula gives
$$
\frac{\mathrm{d}^{1/2}}{\mathrm{d} x^{1/2}} e^{k x} = \frac{1}{\Gamma (1/2)} \frac{\mathrm{d}}{\mathrm{d} x} \int \limits_0^x \mathrm{d} t \frac{e^{k t}}{\sqrt{x-t}} = \frac{1}{\Gamma (1/2)} \frac{\mathrm{d}}{\mathrm{d} x} e^{k x} \int \limits_0^x \mathrm{d} u \frac{e^{- k u}}{\sqrt{u}} = \\
= \frac{1}{\Gamma (1/2)} \frac{\mathrm{d}}{\mathrm{d} x} \frac{e^{k x}}{\sqrt{k}} \int \limits_0^{\sqrt{k x}} \mathrm{d} s \, e^{-s^2}
$$
We can already see that this is not equal to $\sqrt{k} e^{k x}$.
So who is right?
Why do we use the latter formula in almost all cases but somehow we settle for the simpler formula when it comes to the exponential?

Note that both satisfy that if we apply the half-derivative twice, we get the usual first derivative. (In the first case, it's simply because $\sqrt{k} \sqrt{k} = k$; in the second case, there's a proof on Wiki using the properties of beta and gamma functions - plus I verified this numerically even though I couldn't express the integrals in a closed-form.)
I also have hard time accepting this second, complicated, formula, mainly because any integer derivative of the exponential gives exponential, but for the half-derivative we get this weird monstrosity. On the other hand, it should be consistent with the formulas for the half-derivative for all powers of $x$ when it's put together in the infinite series for $e^{k x}$.
Can anyone shed some light on this issue for me, please?

REPLY [3 votes]: Both are right. Think of it like this: what is "the" square root of $i$? Even for $1$ there are two of them, but we forget because of the habit to take the positive root. For $i$ there is no habit available. Things get even more plural when we ask for square roots of $I$, the $n\times n$ identity matrix. There are $2^n$ square roots just among diagonal matrices, and if $A$ is a square root so is $SAS^{-1}$ for any invertible matrix $S$. So we should not expect that square root of $D=\frac{d}{d x}$ is a single something either.
A better analogy here is not derivative but antiderivative, which is, after all $D^{-1}$. There is no "the" antiderivative, there is always a family of them. Again, we tend to overlook that, because they only differ by constants of integration. And we could force uniqueness by declaring the definite integral $\int_0^xf(x)\,dx$ to be "the" antiderivative of $f(x)$. But if we do "the" antiderivative of $e^x$ will be not our favorite $e^x$, but $e^x-1$. To get $e^x$, we should choose $\int_{-\infty}^xf(x)\,dx$ to be "the" antiderivative instead, but then this choice won't work for functions that do not quickly vanish at $-\infty$. There just is no perfect choice. Just as antiderivatives need a reference value to become definitive ($x_0=0$ or $-\infty$ above) so do fractional derivatives more generally. The case of the derivative, and its integer powers, is an exception rather than a rule.
For half-derivatives $D^{\frac12}$ the integral is a bit more involved, in one of the forms
$$D^{\frac12}f(x)=\frac{1}{\Gamma\left(\frac{1}{2}\right)}\frac{d}{dx}\int \frac{f(u)}{(x-u)^\frac12}du,$$
I left out the limits of integration. As shown in Fractional Calculus on Mathpages, if we choose the limits from $0$ to $x$ the answer is $e^x\left(\operatorname{erf}\sqrt{x} + \frac{e^{-x}}{\sqrt{\pi x}}\right)$. On the other hand, with the already familiar choice from $-\infty$ to $x$ the answer is the desired $e^x$. Moreover, then $D^\nu(e^{kx})=k^\nu e^{kx}$ for any $k$ and $\nu$. This is easier to check by using the inverted formula:
$$D^{-\frac12}f(x)=\frac{1}{\Gamma\left(\frac{1}{2}\right)}\int_{-\infty}^x \frac{f(u)}{(x-u)^\frac12}du.$$
The difference in the answers is more visible now, because it is not by an additive constant of integration, but the principle is the same - we get a whole $1$-parameter family of half-derivatives.