TITLE: Algorithm design for enumerating pairs of noncommuting elements up to conjugacy
QUESTION [5 upvotes]: I am trying to write some Magma code that, given a group $G$, returns a list of pairs $(x,y)$ with $x,y\in G$ such that $[x,y]\neq 1$ and such that every pair $(z,w)$ in the group with $[z,w]\neq 1$ is conjugate (disregarding order) to some $(x,y)$ output by the algorithm. I realized when I began that I have no idea of the computational cost of the various steps in the process. I would love to learn a lot about this, but to stay focused, I will ask a restricted question:

Which of the following three algorithm designs will more efficiently solve the above problem? And how will the answer to this be affected by the group order? By number of conjugacy classes?

Algorithm 1a

Compute the action by conjugation of $G$ on the set $Sym^2 G$ of unordered pairs of elements of $G$.
Restrict attention to the pairs that don't commute.
Partition $Sym^2 G$ into orbits and select a representative from each orbit. Return these selections.

Algorithm 1b
Same as algorithm 1 except allow redundancy by working with the ordered pairs $G\times G$ instead of $Sym^2 G$.
Algorithm 2

Decompose $G$ into conjugacy classes.
Ignoring the singleton classes, select a representative $x$ from each class.
For each representative $x$, compute the centralizer $H=Z_G(x)$.
Compute the action of $H$ by conjugation on the complement $G\setminus H$.
From each orbit of this action, select a $y$, and output $(x,y)$.

Thoughts
From a mathematical aesthetics point of view, for example if I were writing a proof, algorithm 1a would be the way to go, clearly. Both other algorithms will produce redundant lists containing $(x,y)$ and $(y',x')$ as separate outputs, where $x,x'$ and $y,y'$ are simultaneously conjugate (unless the group has a very special property). More fundamentally, the set of conjugacy orbits of unordered pairs is actually the thing I'm trying to get a hold of, so that algorithm 1a aims directly for it while the other two come at it from the side.
But OTOH, if I were to try to generate this list by hand, I would absolutely choose algorithm 2. Both other algorithms force me to take hold of the set of pairs (ordered or unordered) all at once, which is quadratic in the group order, whereas algorithm 2 never requires me to contemplate any set bigger than the group itself at any one time.
How would a computer fare with these questions? Is computing orbits of an action on a big set much harder than on a small set? How about computing all those centralizers? Does the need in algorithm 2 to partition $G$ into orbits of numerous different groups cause extra difficulties? What are the relevant issues?
Thanks in advance.

REPLY [4 votes]: There are a number of things you could optimize on, runtime, your own work time implementing code, or memory use. Most likely memory use will become the first obstacle when looking at larger examples (as you can't wait to hope for more memory). Thus in the end an algorithm like your #2 will get you further.
(If you want a quick hack for small examples use #1.)
What I would do is however the following, which uses group conjugation to reduce the problem:

Determine Conjugacy classes
For each class representative $x$, let $C=C_G(x)$. Run (second loop) through all conjugacy classes, let $y$ be the class representative, $D=C_G(y)$.
Compute double cosets $C\setminus G/D$.
For each double coset representative $r$ let $z=y^r$. (This is the $C_i$ orbits on the class of $X_j$ Derek Holt's answer is already mentioning.)
Test the pair $x,z$ for your property.

These pairs $x,z$ are now enumerated only up to $G$-conjugacy. In most cases this means that you never need to store, or run through, something as large as a whole class, but the number of double cosets is significantly smaller.
If you need help with implementing this with GAP -- I cannot comment on Magma -- let me know.


TITLE: Are these two versions of identity theorem of power series equivalent?
QUESTION [2 upvotes]: I have two versions of identity theorem of power series:

From PMA - Rudin

From Analysis I - Amann

It seems to me that the one from PMA - Rudin is more general since it does not require the limit point to be $0$.
My question: How to generalize the theorem from Analysis I - Amann to not require that $(y_j)$ converges to $0$?
Thank you for your help!

REPLY [3 votes]: Lemma. Let the power series $h(x) = \sum_kc_kx^k$ converge for $x\in G := (-R,R)$ and set $A := \{x\in G\,|\,\forall n\in\mathbb N_0 : h^{(n)}(x) = 0\}$. Then either $A=\emptyset$ or $A = G$.

Proof. We only have to prove that $A$ is both open and closed in $G$. The closedness of $A$ follows immediately from $A = \bigcap_n (h^{(n)})^{-1}(\{0\})$. Now, let $x_0\in A$. Then $h(x) = \sum_k\frac{h^{(k)}(x_0)}{k!}(x-x_0)^k$ vanishes completely for $x$ in a neighborhood $U$ of $x_0$. In particular $h^{(n)}(x) = 0$ for $x\in U$ and thus $U\subset A$. $\quad\square$
Now, let $h(x) = \sum_kc_kx^k$ be as above and assume the existence of a sequence $(z_j)\subset G$ such that $z_j\to x_0$, where $x_0\in G$, and $h(z_j) = 0$ for all $j$. Write $h(x) = \sum_kd_k(x-x_0)^k$ for $x$ in a neighborhood $U$ of $x_0$ and $z_j = x_0+e_j$. Then $e_j\to 0$ and $\sum_kd_ke_j^k = 0$ for all $j$. By Amann's result, $d_k = 0$ for all $k$. In particular $h(x) = 0$ for $x\in U$ and thus the set $A$ from the lemma is not empty, which implies $A = G$. Once more in particular, $0\in A$, which implies that $c_k = \tfrac{h^{(k)}(0)}{k!} = 0$ for all $k$.
Let me remark that Amann's result does not help that much here because one can easily show by hand that the existence of $(z_j)$ leads to $x_0\in A$.


text:
\chapter{Conclusion} \label{conclusion}

\section{Closing remarks}

This area of mathematics has proven for us to be a very intriguing
one, uniting elements of algebra and combinatorics from ``pure''
mathematics with ``applied'' ideas from spectral analysis, algorithms,
and computational linear algebra.  We hope that our exposition and
results can spur, in some small way, further interest in the field.

\section{Future work}

Our research in this area has raised many more questions than it has
answered.  We will list several problems which we feel are worthy of
future investigation.

\subsection{Conjecture \ref{the-big-conjecture}}

A proof of Conjecture \ref{the-big-conjecture} would be very nice to
have, especially if it is constructive.  Armed with separating sets
for all $W_n$, we would immediately have an isotypic projection
algorithm.

It would also be useful to find a bound for the sizes of minimal
separating sets for the permutation representation of $W_n$, as
described in Section \ref{perm-sepset-sec}.

\subsection{Other separating sets}

We considered only separating sets consisting of class sums.  Although
these have many nice properties, they are not necessarily optimal.
Other possibilities should be considered.  In particular, thinking of
the Jucys-Murphy elements for $S_n$ (see Section \ref{jucys-murphy}),
one could consider conjugacy classes intersected with subgroups,
or some similar construction.  It is especially suggestive that $S_1
\subgroup S_2 \subgroup \dots \subgroup S_n$ and $W_1 \subgroup W_2
\subgroup \dots \subgroup W_n$ both have a strongly recursive
structure.  Also, the Jucys-Murphy elements separate representations
into finer pieces than isotypic submodules, and can even be used to
compute a genuine discrete Fourier transform; it would be very helpful
to be able to duplicate these properties for $W_n$.

\subsection{Computational bounds}

We concentrated on finding minimal-size separating sets.  However, as
we saw in Section \ref{cyclic-example}, minimal size is not always
best when we actually want to compute projections.  In fact, in order
to say anything about the computational properties of the separating
sets we found, we would have to look at how the eigenspaces of our
elements interact, considering the dimensions of their intersections
as they decompose the space.  Once computational bounds are
established for isotypic projections using our separating sets, we
could evaluate them with respect to other possible separating sets to
find one with the best computational properties.

\subsection{Greedy algorithm}

In Section \ref{separating-sets}, we described a greedy algorithm for
quickly finding separating sets from a character table.  It would be
useful to know how optimal its results are.  As mentioned in Appendix
\ref{reduction}, a greedy algorithm for a related problem (MINIMUM
TEST COLLECTION) has been well studied, and it seems likely that these
results could be brought to bear on the separating set problem.  Since
this method bounds separating set sizes, it is possible that we could
thus obtain a (nonconstructive) proof for Conjecture \ref{the-big-conjecture}.

\subsection{Extensions to iterated wreath products of cyclic groups}

We have only examined the group $W_n = Z_2 \wr \dots \wr Z_2$.  More
generally, the groups $W_{n,r} = Z_r \wr \dots \wr Z_r$ are also of
great interest.  Many of our results about the groups extend to this
case (especially when $r$ is prime); see also \cite{orrison:root} and
\cite{foote:image1}.  In particular, the case $r=4$ gives rise to a
so-called ``wreath product transform'' with useful applications in
image processing (see \cite{foote:image1}).

config: arxiv


text:
\begin{document}

\maketitle

\begin{abstract}
In many inverse problems it is essential to use regularization methods that
preserve edges in the reconstructions, and many reconstruction models have been
developed for this task, such as the Total Variation (TV) approach.
The associated algorithms are complex and require a good knowledge of large-scale
optimization algorithms, and they involve certain tolerances that the user must choose.
We present a simpler approach that relies only on standard computational building
blocks in matrix computations, such as orthogonal transformations, preconditioned
iterative solvers,
Kronecker products, and the discrete cosine transform\,---\,hence the term
``plug-and-play.''
We do not attempt to improve on TV reconstructions, but rather provide an easy-to-use
approach to computing reconstructions with similar properties.
\end{abstract}

\begin{keywords}
Image deblurring, inverse problems, $p$-norm regularization, projection algorithm
\end{keywords}

\begin{AMS}
65F22, 65F30
\end{AMS}

\pagestyle{myheadings}
\thispagestyle{plain}
\markboth{D. CHEN, M. E. KILMER AND P. C. HANSEN}{``PLUG-AND-PLAY'' EDGE-PRESERVING REGULARIZATION}

\section{Introduction}\label{eppsec:intr}

This paper is concerned with  discretizations of linear ill-posed problems,
which arise in many technical and scientific applications such as astronomical
and medical imaging, geoscience, and non-destructive testing
\cite{hansenbook,MuellerSiltanen}.
The underlying model is $b = A\, \bar{x} + \eta$,
where $b$ is the noisy data, the matrix $A$ (which is often structured or sparse)
represents the forward operator, $\bar{x}$ is the exact solution, and $\eta$
denotes the unknown noise.
We present a new large-scale regularization algorithm which is able
to reproduce sharp gradients and edges in the solution.
Our algorithm uses only standard linear-algebra building blocks and is
therefore easy to implement and to tune to specific applications.

For ease of exposition, we focus on image deblurring problems involving
$m\times m$ images $B$ (the blurred and noisy image) and $X$ (the reconstruction).
With $b = \text{vec}(B)$ and $x = \text{vec}(X)$, both of length $n=m^2$,
the $n\times n$ matrix $A$ is determined by the point-spread function (PSF)
and corresponding boundary conditions~\cite{hanaol06}.
This matrix is very ill-conditioned (or rank deficient), and computing the ``naive solution''
$A^{-1}b = \bar{x}+A^{-1}\eta$ (or, in the rank-deficient case, the minimum norm solution)
results in a reconstruction that is completely dominated by the inverted noise~$A^{-1}\eta$.

Classical regularization methods, such as Tikhonov regularization or truncated SVD,
damp the noise component in the solution by suppressing high-frequency components
at the expense of smoothing the edges in the reconstruction.
The same is true for regularizing iterations (such as CGLS or GMRES) based on computing
solutions in a low-dimensional Krylov subspace.
The underlying characteristic in these methods is that regularization is achieved
by projecting the solution onto a low-dimensional \textit{signal subspace} $\mathcal{S}_k$
spanned by $k$, low-frequency basis vectors, with the result that the high-frequency components
are missing, hindering the reconstruction of sharp edges.

The projection approach is a powerful paradigm that can often be tailored to
the particular problem.
While these projected solutions may not always have satisfactory accuracy or details,
they still contain a large component of the desired solution, namely, the
low-frequency component which can be reliably determined from the noisy data.
What is missing is the high-frequency component, spanned by high-frequency basis vectors,
and this component must be determined via our prior information about the desired solution.

This work describes an easy-to-use large-scale method for computing the needed
high-frequency component,
related to the prior information that image must have smooth regions while the gradient
of the reconstructed image is allowed to have some (but not too many) large values.
This idea is similar in spirit to Total Variation regularization, where the gradient
is required to be sparse;
but by relaxing this constraint we arrive at problems that are simpler to solve.
The work can be considered as a continuation of earlier work \cite{PPTSVD2,pptsvd96,mtsvd92}
by one of us; it is also related to the decomposition approach in~\cite{BaglamaReichel}.

The remainder of this paper is organized as follows. In Section~\ref{eppsec:epp} we present
the new edge-preserving algorithm and the convergence analysis. Section~\ref{eppsec:alg} discusses the efficient numerical implementation issues. Section~\ref{eppsec:exp} presents numerical experiments of the new deblurring algorithm and comparisons with other state-of-art deblurring algorithms. The conclusions are presented in Section~\ref{eppsec:con}.

\section{The Projection-Based Edge-Preserving Algorithm}
\label{eppsec:epp}

This section pre\-sents the main ideas of the algorithm, while the implementation details
for large-scale problems are discussed in the next section.

\subsection{Mathematical Model}

Throughout, the matrix $L$ defines a discrete derivative or gradient of the solution
(to be precisely defined later), and $\| \cdot \|_p$ denotes the vector $p$-norm.
The underlying prior information is then that the solution's seminorm $\| L\, x \|_p$,
with $1< p<2$,
is not large (which allows some amount of large gradients or edges in the reconstruction).
The choice of the combination of $L$ and $p$ is important and, of course, somewhat problem dependent;
the matrix $L$ used here is the $2m(m-1)\times n$ matrix given by
  \[
    L = \begin{pmatrix} \!\!\!\! L_1 \otimes I \\ \ \ I \otimes L_1 \end{pmatrix}, \quad \mbox{where} \quad
    L_1 = \begin{pmatrix} -1 & 1 & 0& \cdots & 0 \\ 0 & -1 & 1& \cdots & 0 \\
          \vdots  & \vdots  & \ddots  & \ddots & \vdots \\ 0 & 0 & \cdots & -1& 1 \end{pmatrix}
    \in \mathbb{R}^{m-1\times m},
  \]
where $\otimes$ is the Kronecker product \cite{PPTSVD2}.
The one-dimensional null space $\mathcal{N}(L)$ of this matrix is spanned by
the $n$-vector $e$ of all ones.

Assume $W_k\in \mathbb{R}^{n\times k}$ is a matrix with orthonormal columns that span the
signal subspace $\mathcal{S}_k$, and let $W_0$ be the matrix containing the orthonormal
basis vectors for the complementary space $\mathcal{S}_k^\perp$.
The fundamental assumption is that the columns of $W_k$ represent ``smooth" modes
in which it is possible to distinguish a substantial component of the signal from the noise.
In other words, with the model from Section~1, we assume that
  \begin{equation}
  \label{eq:separate}
    \| W_k^T \bar{x} \|_2 \gg \| W_k^T (A^{-1} \eta) \|_2 \ .
  \end{equation}
This ensures that we can compute a good, but smooth, approximation to $\bar{x}$
as
  \[
    x_k = W_k\, y_k \ , \qquad y_k = \mathrm{argmin}_y \| (A\, W_k)\, y - b \|_2 \ ,
  \]
and we refer to the minimization problem for $y_k$ as the \textit{projected problem},
which we assume is easy to solve.
To obtain a reconstruction with the desired features,
our strategy is then to compute the solution of the following modified projection problem
  \be
  \label{eppeq:epp}
    \min_{x\in \mathcal{B}} \| L \, x \|_p \quad \text{with} \quad
    \mathcal{B} = \{ x: x = \argminz \| (A\,W_k W_k^T)z - b \|_2 \} \ ,
  \ee
with  $L$ defined above.
As we shall see, we can express the solution to (\ref{eppeq:epp}) as the low-frequency
solution $x_k$ plus a high-frequency correction.

\subsection{Uniqueness Analysis}

Eld\'en \cite{elden82} provides an explicit solution of~\eqref{eppeq:epp} for the case
$p=2$ and proves the uniqueness condition for the minimizer.
The MTSVD algorithm \cite{mtsvd92} corresponds the case where $p=2$ and $W_k$ consists of
the first $k$ right singular vectors, while the PP-TSVD algorithm \cite{pptsvd96}
and its 2D extension \cite{PPTSVD2} correspond to the same $W_k$ and $p=1$.
In this work, we extend these results by solving~\eqref{eppeq:epp} for $1<p<2$ and for
different choices of $W_k$.
Below we present results that give conditions for the existence and uniqueness of the
solution to~\eqref{eppeq:epp}.

\begin{lemma} \label{lem:epp}
The linear $p$-norm problem
  \be
      \argminx \| A\, x-b \|_p^p, \qquad p>1, \label{eppeq:pmin}
  \ee
has a unique minimizer if and only if $A$ has full column rank.
\end{lemma}

\begin{proof}

The function $\| x \|_p^p$ is strictly convex for $1 < p$, or equivalently,
the Hessian $H(x)$ of $\|x\|_p^p$ is positive definite for all $x$.
This implies that $\|Ax-b\|_p^p$ is strictly convex (or equivalently, its Hessian $A^* H(x)\, A$
is positive definite for all $x$) if and only if $A$ has full column rank.
It follows from strict convexity that the minimizer is unique.\footnote{We thank Martin S.
Andersen for help with this proof.}
$\qquad$
\end{proof}

\begin{theorem} \label{thm:epp}
The modified projection problem \eqref{eppeq:epp} has a unique minimizer if and only if
$\mN(AW_kW_k^T) \cap \mN(L) = \{ 0 \}$.
\end{theorem}

\begin{proof}
From \cite{elden82}, the constraint set $\mathcal{B}$ in~\eqref{eppeq:epp} can be written as
  \[
    \mathcal{B} = \{ x: x = (AW_kW_k^T)^\dagger b + Px', \ x' \mbox{ arbitrary} \},
  \]
where $\dagger$ denotes the Moore-Penrose pseudoinverse~\cite{bjorck96} and
  \[
    P = I-(AW_kW_k^T)^\dagger (AW_kW_k^T)
  \]
is the orthogonal projector onto $\mN(AW_kW_k^T)$.
Let $\tilde{b} = (AW_kW_k^T)^\dagger b$. Solving the constrained minimization~\eqref{eppeq:epp}
is equivalent to solving the following unconstrained problem
  \[
    \min_{\tilde{x}} \| LP\tilde{x} - (-L\tilde{b}) \|_p.
  \]
By Lemma~\ref{lem:epp}, the above minimization problem has a unique solution if and only if
$\mN(LP) = \{ 0 \}$.
This is true for $P = I-(AW_kW_k^T)^\dagger (AW_kW_k^T)$, the projection onto
$\mN(AW_kW_k^T)$, if and only if $\mN(AW_kW_k^T) \cap \mN(L) = \{ 0 \}$. $\qquad$
\end{proof}

\subsection{Algorithm}

It follows from the proof of Theorem~\ref{thm:epp} that we can solve the modified projection
problem \eqref{eppeq:epp} in two steps.
We first compute an approximate solution $x_k \in \mathcal{S}_k$ that contains the smooth components,
and then we compute the edge-correction component $x_0$ in the orthogonal complement $\mS_k^\perp$.
As a result,
  \[
    x = x_k + x_0 = W_k y_k + W_0 y_0 \ ,
  \]
where $y_k$ is the solution to the projected problem, and $y_0$ is the solution to an
associated $p$-norm problem.
These two solutions are computed sequentially, as shown in the EPP Algorithm~\ref{eppalg:eppalg}.

\begin{algorithm}[h]
\caption{Edge-Preserving Projection (EPP) Algorithm}
\begin{algorithmic}[1]
\State Compute the smooth component $x_k = W_k y_k$ using the projected problem
  \be
  \label{eppeq:eppy}
    y_k = \argminy \| (AW_k)y-b \|_2.
  \ee
\vspace*{-4mm}
\State Compute the correction component $x_0 = W_0 y_0$ using the $p$-norm problem
  \be
    y_0 = \argminy \| (LW_0)y - (-LW_k y_k)  \|_p. \label{eppeq:eppz}
  \ee
\vspace*{-4mm}
\State The regularized solution is then
\vspace*{-2mm}
  \be
    x = W_k y_k + W_0 y_0 .
  \ee
\vspace*{-6mm}
\end{algorithmic}\label{eppalg:eppalg}
\end{algorithm}

\subsection{Choosing Projection Spaces}

From Lemma~\ref{lem:epp}, a sufficient condition for the uniqueness of $x$ is that
both $AW_k$ and $LW_0$ have full column rank, for then \eqref{eppeq:eppy} and \eqref{eppeq:eppz}
in the EPP Algorithm have unique solutions $y_k$ and $y_0$, correspondingly.
In principle, we can choose any subspace $\mS_k$ and its orthogonal complement $\mS_k^\perp$
with corresponding $W_k$ and $W_0$.
But in practice, however, in order to have a useful and efficient numerical implementation,
we must choose suitable basis vectors for $\mS_k$ with the following requirements:
\begin{itemize}
\item The matrix $W_k$ must separate signal from noise according to \eqref{eq:separate}.
\item The matrices $AW_k$ and $LW_0$ must have full column rank.
\item There are efficient algorithms to compute multiplications with $W_k$ and $W_0$
      and their transpose.
\end{itemize}

\subsubsection{Singular Vectors}

The MTSVD and PP-TSVD algorithms proposed in \cite{PPTSVD2,pptsvd96,mtsvd92} use the first
$k$ singular vectors as the basis vectors for $\mathcal{S}_k$.
In this case we have the following result.

\begin{theorem}\label{thm:svd}
Assume that $W_k = [v_1, v_2, \cdots,  v_k]$, where $ v_i$ are right singular vectors of $A$
corresponding to nonzero singular values.
Then the modified projection problem \eqref{eppeq:epp} has a unique solution if and only if
$\id\notin \mathcal{N}(A) = \mathrm{range}(W_0)$.
\end{theorem}

\begin{proof}
Since $W_k W_k^T$ is the orthogonal projector onto $\mathrm{range}(W_k)$ it
follows that
 \[
   \null( A W_k W_k^T) = \mathrm{range}(W_0) = \mathrm{span}\{ v_{k+1},\ldots,v_n \},
 \]
and the requirement from Theorem \ref{thm:epp} becomes
$\mathrm{range}(W_0) \cap \{ e \} \notin \{ 0 \}$, which is clearly satisfied
if $\id\notin \mathrm{range}(W_0)$. $\qquad$
\end{proof}

For blurring operators, the SVD-based subspace $\mS_k$ contains low-frequency components,
while $\mS_k^\perp$ contains relatively high-frequency components.
It is therefore very likely that the projection of $\id$ onto $\mS_k$ is not zero, and in fact
this is easy to check.

\subsubsection{Discrete Cosine Vectors}

Another suitable set of basis vectors for $\mathcal{S}_k$ are those associated with
spectral transforms such as the discrete sine or cosine transforms (DST or DCT)
and their multidimensional extensions~\cite{haje08,hanaol06}.
Recall that for 1D signals of length $m$, the orthogonal DCT matrix $C$ has elements
\[
  c_{ij} = \left\{
  \begin{array}{l l}
    \sqrt{\nicefrac{1}{m}} & \quad \text{if $i=0$}\\[1mm]
    \sqrt{\nicefrac{2}{m}}\,\cos\!\left(\frac{(2j+1)i\pi}{2m}\right) & \quad \text{if $i>0$}\\
  \end{array} \right.
  \quad \mbox{for} \quad i,j = 0,1,2,\cdots, m-1 \ .
\]
The 2-dimensional DCT matrix is the Kronecker product $C \otimes C$ of the
above matrix~\cite{vapi93}.
The DCT basis vectors, which are the \textit{rows} of the DCT matrix, have the
desired spectral properties.
Multiplications with $W_k$ and $W_0$ and their transposes are equivalent to computing
either a DCT transform or its inverse, which is done by fast algorithms similar to the FFT\@.
For this basis we have the following result.

\begin{theorem}\label{thm:dct}
Let the columns of $W_k$ be the first $k$ 2D DCT basis vectors.
Then the modified projection problem \eqref{eppeq:epp} has a unique solution if and only if
$\id\notin \mN(A)$.
\end{theorem}

\begin{proof}
From the definition of the DCT it follows that $w_1 = e / \| e \|_2$ and hence
$A W_k W_k^T e = A e$, and therefore
$\mathcal{N}(A W_k W_k^T) \cap \{ e \} = \{ 0 \} \Leftrightarrow
 Ae \neq 0 \Leftrightarrow e \notin \mathcal{N}(A)$.  $\qquad$
\end{proof}

\input{OneDimEx.tex}

\section{Computational Issues and Numerical Implementations}
\label{eppsec:alg}

While the above analysis guarantees the existence and uniqueness of the solution to
\eqref{eppeq:epp}, it is critical to develop an efficient numerical implementation
for large-scale problems, which must take the following three issues into account:
\begin{itemize}
\item efficiently construct, or perform operations with, the basis vectors
      for the subspace~$\mS_k$,
\item robustly choose the optimal dimension $k$ of $\mS_k$, and
\item efficiently solve the $p$-norm minimization problem \eqref{eppeq:eppz}.
\end{itemize}
The optimal subspace dimension can be computed standard parameter-choice algorithms
\cite{hansenbook}; here we use the GCV method as explained in Section~\ref{eppsec:exp}.

\subsection{Working with the Projection Spaces}

As discussed above, the singular vectors and the 2D DCT matrix can be used
as the basis vectors for $\mS_k$ and $\mS_k^\perp$.
Here we will address numerical implementation issues with these choices.

For large-scale deblurring problems it is impossible to obtain $W_k = [v_1, \cdots, v_k]$
by computing the SVD of the blurring matrix $A$ without utilizing its structure.
Fortunately, in many problems the underlying point-spread function is separable or can be
approximated by a separable one \cite{hanaol06,kana99,nangpe03,vapi93}.
Hence, the blurring matrix $A$ can be represented as a Kronecker product $A_1 \otimes A_2$.
Given the SVDs of the two matrices $A_1 = U_1 \Sigma_1 V_1^T$ and $A_2 = U_2 \Sigma_2 V_2^T$,
the right singular matrix of $A$ is (or can be approximated by)
$V = \Pi (V_1 \otimes V_2)$, where the permutation matrix $\Pi$ ensures that the
ordering of the singular vectors is in accordance with decreasing singular values
(i.e., the diagonal elements of $\Pi(\Sigma_1 \otimes \Sigma_2)$).

For the DCT basis, multiplication with the $m\times m$ DCT matrix $C$ is implemented
in an very efficient way using the FFT algorithm, requiring only $O(m\log m)$ operations,
and a similar fast algorithm is available for the 2D DCT\@.
The multiplications with $W_k$ and its transpose are equivalent to applying
either the DCT or its inverse.
Therefore, it is unnecessary to form the matrix $W_k$ explicitly.

\subsection{Iteratively Reweighted Least Squares and AMG Preconditioner}

The key to the success of the EPP Algorithm is an efficient solver for the
$p$-norm minimization problem~\eqref{eppeq:eppz}.
A standard and robust approach is to use the iteratively reweighted least squares (IRLS)
method \cite{bjorck96,osborne85,irls88}, which is identical to Newton's method
with line search.
This approach reduces the $p$-norm problem to the solution of a sequence of
weighted least squares problems, which can be solved using iterative solvers.
Osborne \cite{osborne85} shows that the IRLS method is convergent for $1<p<3$.

\begin{algorithm}[t]
\caption{Iterative Reweighted Least Squares (IRLS) Algorithm
for minimizing $f(\hx)=\|\hA\,\hx -\hb\|_p$}
\label{alg:irls}
\begin{algorithmic}[1]
\State $\hx_0 = 0$ \quad (starting vector)
\For{$j = 0, 1, 2, \ldots$ }
\State $\widehat{r}^j = \hb - \hA\,\hx^j $
\State $D_j = \diag(|r^j|^{(p-2)/2})$
\State $z^j = \mathrm{argmin}_z \| D_j ( \hA\, z - \widehat{r}^j )  \|_2$ \quad (solved iteratively)

\State $\alpha_j = \mathrm{argmin}_{\alpha} f(\hx^j+\alpha z^j)$ \quad (line search)
\State $\hx^{j+1} = \hx^j + \alpha^j z^j$
\EndFor
\end{algorithmic}
\end{algorithm}

For convenience, we briefly summarize the IRLS algorithm for solving the linear
$p$-norm problem $\min_{\hx} \| \hA \, \hx  - \hb  \|_p$.
We denote the $j$th iteration vector by $\hx^j$, and we introduce
the diagonal matrix $D_j$ determined by $j$th residual vector $\widehat{r}^j=\hb - r\hA \, \hx^j $ as
  \[
    D_j = \diag\Bigl( \bigl| \hA \, \hx^j - \hb \bigr|^{\frac{p-2}{2}} \Bigr).
  \]
The Newton search direction $z^j$ is identical to the solution of the weighted least squares problem
  \be
  \label{eppeq:irlsls}
    \min_{z} \| D_j \bigl( \hA \, z - \widehat{r}^j) \bigr \|_2.
  \ee
For $1<p<2$, as the iteration vector $\hx ^j$ gets close to the solution, the diagonal elements
in $D^j$ increase to infinity, and this tendency increases as $p$ approach~1.
Hence, the matrix $D_j\hA $ in \eqref{eppeq:irlsls} becomes increasingly ill-conditioned
as the iterations converge.
It is therefore difficult to find a suitable preconditioner for the least squares
problem~\eqref{eppeq:irlsls}.

Consider the corresponding normal equations
  \[
    \hA ^TD_j^2\hA \, z^j = \hA ^TD_j^2 \, \widehat{r}^j = \hA ^TD_j^2(\hb - \hA \, \hx^j ) ,
  \]
and define the new variable $q^{j} = z^j+\hx^j$. The normal equations can then be rewritten
  \be
  \label{eppeq:irls}
    \hA ^TD_j^2\hA \, q^{j}= \hA ^TD_j^2 \, \hb .
  \ee
The benefit of the above transformation is that the right-hand side in the new
system \eqref{eppeq:irls} depends on iteration $j$ only through $D^j$,
which is known in the $j$th iteration.\footnote{We thank Eric de Sturler for pointing this out.}

For our algorithm, it follows from \eqref{eppeq:eppz} that $\hA = LW_0$ and
$\hb = -LW_k y_k$, so \eqref{eppeq:irls} can be rewritten as
  \be
  \label{eppeq:irls2}
    W_0^T (L^TD_j^2L) W_0 \, q^j = - W_0^T (L^T D_j^2L) W_k y_k.
  \ee
Since the condition number increases as the IRLS algorithm converges to the solution,
preconditioning is necessary in solving \eqref{eppeq:irls2}.
Recall that $L$ is a gradient operator, and hence $L^TD_j^2L$ represents a diffusion
operator with large discontinuities in the diffusion coefficients.
Algebraic multi-grid (AMG) methods are robust when the diffusion coefficients are
discontinuous and vary widely~\cite{rust87, troo01}.
Therefore, we employ an AMG method to develop a right preconditioner $M$ for~\eqref{eppeq:irls2}.
The right-preconditioned problem is
  \be
  \label{eppeq:amgirls}
    [W_0^T (L^TD_j^2L) W_0 M] \tilde{q}^j = - W_0^T (L^T D_j^2L) W_k y_k,
  \ee
where $q^j = M\tilde{q}^j$.
In our implementation, given a vector $v$, the matrix-vector multiplication $w=Mv$
is implemented in three steps:
\begin{enumerate}
\item Compute $\tilde{v} = W_0 v$.
\item Use the AMG method to solve $(L^TD_j^2L)u = \tilde{v}$ for $u$.
\item Compute the result $w=W_0^Tu$.
\end{enumerate}

The matrix $W_0^T (L^TD_j^2L) W_0$ is symmetric positive definite if $D_j^2$ is positive definite.
If not, positive definiteness of $D_j^2$ can be guaranteed by adding a small positive number
to the diagonal elements.
A first thought may be to solve \eqref{eppeq:amgirls} with the conjugate gradient (CG)
method; but this requires that the preconditioner $M$ is also symmetric and positive definite.
In our implementation we use the Gauss-Seidel method in the pre- and post-relaxations of
the AMG method, and hence the AMG residual reduction operator is not symmetric~\cite{rust87},
and consequently the preconditioner is not symmetric.
Instead we solve \eqref{eppeq:irlsls} with the GMRES algorithm with right AMG
preconditioning~\cite{gmres86}.

\section{Numerical Results}
\label{eppsec:exp}

We present numerical experiments using the EPP algorithm, and we
perform a brief comparison with Total Variation deblurring.
To better visualize the impact of the high-frequency correction we use Matlab's
colormap \texttt{Hot} for the first example, which varies smoothly from black
through shades of red, orange, and yellow, to white, as the intensity increases.
Throughout we use the $256 \times 256$ ``cameraman'' test image. All the
numerical simulations are performed using Matlab R2009b on Windows 7 x86
32-bit system. The C compiler used to build AMG preconditioner MEX-files is
Microsoft Visual Studio 2008.

\subsection{Image Quality, PSFs, and Algorithm Parameters}

 \begin{figure}[t]
\centering
\includegraphics[width= 0.25\textwidth]{psfGauss.png}
\hspace{.25in}
\includegraphics[width= 0.25\textwidth]{psfDisk.png}
\caption{Gaussian PSF with $\sigma=5$ (left) and out-of-focus PSF with $r=5$ (right).}
\label{fig:psfs}
\end{figure}

The ``noise level" of a test image is defined as $\| \eta \|_2/\| b \|_2$.
The quality of the restored images is measured by the relative error
$\| x_\text{restored} - \bar{x} \|_2 / \| \bar{x} \|_2$ and by the MSSIM \cite{ssim04}
(for which a larger value is better).
In our experiments the test images are generated with two common types of PSFs,
Gaussian blur and out-of-focus blur, and we use reflexive boundary conditions
in the restorations.
The elements of the \textbf{Gaussian PSF} are
  \[
    p_{ij} = \exp\!\left(-\frac{\sigma^2}{2} \Bigl( (i-k)^2 + (j-\ell)^2 \Bigr) \right),
  \]
and the elements of the \textbf{out-of-focus PSF} are
  \[
    p_{ij} = \left\{ \begin{array}{ll}
      \nicefrac{1}{\pi r^2} & \quad \text{if} \quad (i-k)^2+(j-l)^2 \leq r^2 \\[1mm]
    0 & \quad \text{otherwise} ,
  \end{array} \right.
\]
where $(k,\ell)$ is the center of the PSF, and $\sigma$ and $r$ are parameters
that determine the amount of blurring; see Fig.~\ref{fig:psfs}.
Both are doubly symmetric, but the latter is not separable, and therefore it is not
possible to efficiently compute the exact SVD of the corresponding matrix~$A$.

To compute the subspace dimension $k$ we use the GCV method \cite{hansenbook}, which can
be implemented very efficiently when the singular vectors or the DCT basis are used.
The GCV function can be expressed as
  \[
    G(k) = \frac{\sum_{i=k+1}^n \beta_i^2}{(n-k)^2}, \quad \mbox{for} \qquad k = 1, 2, \cdots, n-1,
  \]
where $\beta_i = w_i^T b$ ($w_i$ being either the left singular vectors $u_i$ or
the DCT basis vectors).
As noted in \cite{chnaol08}, the GCV method very often provides a parameter
that is too large.
Also, in some of our experiments we assume that the singular vectors are approximated by a
Kronecker product, which might be not accurate.
Hence, we choose $k$ to be equal to $2/3$ of the output from GCV algorithm,
where the factor of $2/3$ was chosen on the basis of numerous experiments.

The stopping criteria used in the iterative methods were chosen based on exhaustive experiments
(see \cite{DonghuiThesis} for details) to balance computational time against the quality of the reconstruction.
Results computed with smaller tolerances than those used here are qualitatively similar to those computed
with the chosen tolerances, but the computational time is much longer.

\subsection{Performance of the EPP Algorithm}

In the EPP algorithm the norm parameter $p$ can be any number between $1$ and $2$.
For smaller $p$, the solution tends to have sharper edges, but as $p$ gets closer to $1$
the $p$-norm minimization in \eqref{eppeq:eppz} becomes more ill-conditioned requiring
much more computational work, while there are very little improvement of the restored images.
Hence, we show computed results with $p=1.01, 1.05, 1.1$, and $1.2$.

 \begin{figure}[t]
\centering
\includegraphics[width= 0.24\textwidth]{dctDcameraman_5_5.png}

\includegraphics[width= 0.24\textwidth]{dctDcameraman_5_5_Xk.png}

\includegraphics[width= 0.24\textwidth]{dctDcameraman_5_5_X0.png}

\includegraphics[width= 0.24\textwidth]{dctDcameraman_5_5_X.png}
\vspace*{-2mm}
\caption{\label{fig:dctdisk}
DCT-EPP algorithm, out-of-focus blur with $r=5$, noise level $5\%$.
Left to right:\ blurred noisy image, low-frequency component $x_k$, high-frequency
component $x_0$, and reconstruction.}
\vspace{5mm}

\includegraphics[width= 0.24\textwidth]{dctcameraman_5_1.png}

\includegraphics[width= 0.24\textwidth]{dctcameraman_5_1_Xk.png}

\includegraphics[width= 0.24\textwidth]{dctcameraman_5_1_X0.png}

\includegraphics[width= 0.24\textwidth]{dctcameraman_5_1_X.png}
\vspace*{-2mm}
\caption{\label{fig:dctgauss}
DCT-EPP algorithm, Gaussian blur with $\sigma=5$, noise level $1\%$.}
\vspace{5mm}

\centering
\includegraphics[width=0.23\textwidth]{svdcameraman_5_1.png}
\hspace{.01in}
\includegraphics[width=0.23\textwidth]{svdcameraman_5_1_Xk.png}
\hspace{.01in}
\includegraphics[width=0.23\textwidth]{svdcameraman_5_1_X0.png}
\hspace{.01in}
\includegraphics[width=0.23\textwidth]{svdcameraman_5_1_X.png}
\vspace*{-2mm}
\caption{\label{fig:svdgauss}
SVD-EPP algorithm, Gaussian blur with $\sigma=5$, noise level is $1\%$.}
\end{figure}

Table~\ref{tab:dct} shows the results of the restored out-of-focus blurred images
using the DCT-EPP algorithm.
The blur radius $r$ varies from $5$ to $15$ pixels, and the noise level
varies from 1\% to 10\%.
The table reports the computed truncation parameter $k$, the relative errors, and the MSSIM
for both $x_k$ and the final restored image.
Compared to the restored quality of $x_k$, the latter image has larger MSSIM
and smaller relative error, demonstrating that the correction step~\eqref{eppeq:eppz}
improves the image quality.
This is illustrated by the example in Figure~\ref{fig:dctdisk}.
The restored images computed using smaller $p$ are generally better than the results using larger $p$.
The corresponding results for Gaussian blur with $\sigma=5,10,15$, still using the DCT-EPP algorithm,
are also shown in Table~\ref{tab:dct}; see Fig.~\ref{fig:dctgauss} for an example.

Table \ref{tab:svd} summarizes the results for the SVD-EPP algorithm, again for out-of-focus
and Gaussian blur; see also Fig.~\ref{fig:svdgauss}.
For the Gaussian blur, the performance is similar to the DCT case.
The out-of-focus blur, however, is not separable. Therefore, we feed the SVD-EPP algorithm the

approximate singular vectors obtained from a Kronecker-product approximation of~$A$
(with Toeplitz blocks). As this particular Kronecker approximation to the true SVD is not sufficiently good,
the algorithm performs poorly, whereas use of the DCT basis gives better quality reconstructions.

\begin{landscape}
\renewcommand{\arraystretch}{1.1}
\begin{table}
\caption{Comparison of the quality of images restored by the DCT-EPP Algorithm
for $p=1.01,1.05,1.1,1.2$.
As expected, $p=1.01$ gives the best results.}~\label{tab:dct}
\centering
\begin{tabular}{|c|c|c|cc|cccc|cccc|} \hline
$r$ & noise & $k$ & \multicolumn{2}{c|}{$x_k$} & \multicolumn{8}{c|}{$x_{\mathrm{restored}}$} \\
\cline{4-13}
or & level &  & relative  & MSSIM & \multicolumn{4}{c|}{relative error} & \multicolumn{4}{c|}{MSSIM} \\
\cline{6-13}
$\sigma$ & \% &  & error &  & $p=1.01$ & $p=1.05$ & $p=1.1$ & $p=1.2$ &  $p=1.01$& $p=1.05$ &  $p=1.1$& $p=1.2$\\
\hline\hline
\multicolumn{13}{|c|}{\textbf{Out-of-focus PSF}} \\
\hline
5 & 1 & 2519 & 0.148 & 0.617& 0.135& 0.135& 0.135& 0.136& 0.707& 0.706& 0.704& 0.704  \\
5 & 5 & 1555 & 0.161 & 0.600& 0.151& 0.151& 0.151& 0.151& 0.668& 0.667& 0.666& 0.665  \\
5 & 10 & 1254 & 0.168 & 0.568& 0.159& 0.158& 0.158& 0.159& 0.637& 0.644& 0.644& 0.642  \\
\hline
10 & 1 & 1310 & 0.177 & 0.523& 0.158& 0.160& 0.159& 0.159& 0.640& 0.629& 0.634& 0.633  \\
10 & 5 & 427 & 0.204 & 0.493& 0.192& 0.192& 0.192& 0.193& 0.579& 0.578& 0.577& 0.572  \\
10 & 10 & 418 & 0.205 & 0.487& 0.194& 0.193& 0.194& 0.195& 0.570& 0.573& 0.565& 0.561  \\
\hline
15 & 1 & 772 & 0.195 & 0.488& 0.174& 0.174& 0.173& 0.175& 0.609& 0.611& 0.614& 0.604  \\
15 & 5 & 204 & 0.233 & 0.467& 0.223& 0.223& 0.223& 0.224& 0.531& 0.529& 0.527& 0.522  \\
15 & 10 & 203 & 0.234 & 0.463& 0.224& 0.224& 0.224& 0.225& 0.524& 0.523& 0.525& 0.519  \\
\hline\hline
\multicolumn{13}{|c|}{\textbf{Gaussian PSF}} \\
\hline
5 & 1 & 1188 & 0.168 & 0.577& 0.158& 0.158& 0.158& 0.158& 0.657& 0.656& 0.656& 0.654  \\
5 & 5 & 678 & 0.188 & 0.535& 0.176& 0.177& 0.177& 0.177& 0.613& 0.606& 0.610& 0.608  \\
5 & 10 & 560 & 0.196 & 0.508& 0.185& 0.186& 0.185& 0.186& 0.587& 0.580& 0.583& 0.579  \\
\hline
10 & 1 & 337 & 0.214 & 0.490& 0.202& 0.202& 0.203& 0.203& 0.563& 0.565& 0.560& 0.556  \\
10 & 5 & 242 & 0.228 & 0.463& 0.219& 0.219& 0.219& 0.220& 0.524& 0.524& 0.522& 0.516  \\
10 & 10 & 174 & 0.240 & 0.476& 0.231& 0.231& 0.231& 0.232& 0.531& 0.531& 0.530& 0.523  \\
\hline
15 & 1 & 167 & 0.241 & 0.471& 0.232& 0.232& 0.232& 0.233& 0.529& 0.525& 0.526& 0.522  \\
15 & 5 & 117 & 0.252 & 0.461& 0.240& 0.241& 0.240& 0.241& 0.525& 0.522& 0.524& 0.518  \\
15 & 10 & 100 & 0.257 & 0.460& 0.246& 0.246& 0.246& 0.247& 0.518& 0.517& 0.514& 0.510  \\
\hline
\end{tabular}
\end{table}
\end{landscape}

\begin{landscape}
\renewcommand{\arraystretch}{1.1}
\begin{table}
\caption{Comparison of the quality of images restored by the SVD-EPP Algorithm; similar to
Table~\ref{tab:dct}.}~\label{tab:svd}
\centering
\begin{tabular}{|c|c|c|cc|cccc|cccc|} \hline
$r$ & noise & $k$ & \multicolumn{2}{c|}{$x_k$} & \multicolumn{8}{c|}{$x_{\mathrm{restored}}$} \\
\cline{4-13}
or & level &  & relative  & MSSIM & \multicolumn{4}{c|}{relative error} & \multicolumn{4}{c|}{MSSIM} \\
\cline{6-13}
$\sigma$  & \% &  & error &  & $p=1.01$ & $p=1.05$ & $p=1.1$ & $p=1.2$ &  $p=1.01$& $p=1.05$ &  $p=1.1$& $p=1.2$\\
\hline\hline
\multicolumn{13}{|c|}{\textbf{Out-of-focus PSF}} \\
\hline
5 & 1 & 5648 &  0.239 &  0.448 & 0.246 & 0.246& 0.245& 0.245& 0.530 & 0.531& 0.530& 0.530  \\
5 & 5 & 2102 &  0.218 &  0.499 & 0.223 & 0.224& 0.223& 0.222& 0.565 & 0.564& 0.565& 0.565  \\
5 & 10 & 1311 &  0.219 &  0.492 & 0.219 & 0.218& 0.217& 0.216& 0.578 & 0.579& 0.578& 0.577  \\
\hline
10 & 1 & 2483 &  0.261 &  0.380 & 0.267 & 0.267& 0.267& 0.266& 0.443 & 0.444& 0.445& 0.444  \\
10 & 5 & 835 &  0.234 &  0.454 & 0.237 & 0.236& 0.236& 0.235& 0.525 & 0.525& 0.526& 0.523  \\
10 & 10 & 468 &  0.234 &  0.450 & 0.234 & 0.233& 0.233& 0.232& 0.535 & 0.533& 0.532& 0.527  \\
\hline
15 & 1 & 2173 &  0.323 &  0.202 & 0.334 & 0.335& 0.331& 0.331& 0.235 & 0.237& 0.234& 0.234  \\
15 & 5 & 486 &  0.254 &  0.433 & 0.253 & 0.254& 0.253& 0.252& 0.509 & 0.512& 0.509& 0.502  \\
15 & 10 & 336 &  0.250 &  0.421 & 0.249 & 0.250& 0.249& 0.248& 0.493 & 0.495& 0.494& 0.489  \\
\hline\hline
\multicolumn{13}{|c|}{\textbf{Gaussian PSF}} \\
\hline
5 & 1 & 1188 &  0.168 &  0.577 & 0.157 & 0.157& 0.157& 0.158& 0.664 & 0.662& 0.662& 0.657  \\
5 & 5 & 679 &  0.188 &  0.536 & 0.175 & 0.176& 0.175& 0.177& 0.621 & 0.616& 0.618& 0.611  \\
5 & 10 & 561 &  0.196 &  0.508 & 0.184 & 0.184& 0.185& 0.186& 0.594 & 0.589& 0.588& 0.580  \\
\hline
10 & 1 & 338 &  0.213 &  0.489 & 0.202 & 0.201& 0.202& 0.203& 0.564 & 0.570& 0.568& 0.559  \\
10 & 5 & 242 &  0.228 &  0.463 & 0.219 & 0.219& 0.219& 0.219& 0.527 & 0.523& 0.520& 0.512  \\
10 & 10 & 174 &  0.240 &  0.476 & 0.231 & 0.230& 0.231& 0.232& 0.535 & 0.536& 0.531& 0.524  \\
\hline
15 & 1 & 168 &  0.241 &  0.471 & 0.232 & 0.232& 0.232& 0.233& 0.535 & 0.533& 0.530& 0.522  \\
15 & 5 & 124 &  0.250 &  0.470 & 0.238 & 0.239& 0.239& 0.241& 0.534 & 0.532& 0.529& 0.519  \\
15 & 10 & 100 &  0.257 &  0.460 & 0.245 & 0.245& 0.245& 0.247& 0.524 & 0.520& 0.518& 0.508  \\
\hline
\end{tabular}
\end{table}
\end{landscape}

\subsection{Comparison with Total Variation Deblurring\label{eppsec:tv}}

We conclude by briefly comparing the performance of the EPP algorithm with the TV deblurring algorithm,
using the algorithm proposed in~\cite{hungwe08}.
In order to avoid giving our algorithm an advantage, the parameters of the TV algorithm
were chosen to optimize the MSSIM (which obviously requires the true image).
As shown in Table \ref{tab:TV}, the images restored by the TV method qualitatively
have better quality as those computed by EPP algorithm as measured by both the
relative error and the MMSIM\@.
This demonstrates that the EPP algorithm can be a computationally attractive
alternative to TV.

\begin{table}[t]
\renewcommand{\arraystretch}{1.1}
\caption{Comparison of the restored images by the TV and DCT-EPP algorithms with $p=1.01$.
There is no dramatic difference between the performance of the two algorithms.}~\label{tab:TV}
\centering
\begin{tabular}{|c|c|cc|cc|cc|cc|} \hline
$r$ & noise & \multicolumn{2}{c|}{relative error}  & \multicolumn{2}{c|}{MSSIM}
   & \multicolumn{2}{c|}{relative error}  & \multicolumn{2}{c|}{MSSIM} \\
\cline{3-10}
$\sigma$ & level \% & EPP & TV & EPP & TV & EPP & TV & EPP & TV \\
\hline\hline
 & & \multicolumn{4}{|c|}{\textbf{Out-of-focus PSF}} & \multicolumn{4}{|c|}{\textbf{Gaussian PSF}} \\
\hline
 5 & 1 & 0.135 & 0.150 & 0.707 & 0.677 & 0.158 & 0.167& 0.657 & 0.624  \\
 5 & 5 & 0.151 & 0.153 & 0.668 & 0.656 & 0.176 & 0.176& 0.613 & 0.587  \\
 5 & 10 &0.159 & 0.159 & 0.637 & 0.626 & 0.185 & 0.182& 0.587 & 0.559  \\
\hline
10 & 1 & 0.158 & 0.172 & 0.640 & 0.606 & 0.202 & 0.205& 0.563 & 0.528  \\
10 & 5 & 0.192 & 0.180 & 0.579 & 0.571 & 0.219 & 0.213& 0.524 & 0.507  \\
10 & 10 &0.194 & 0.187 & 0.570 & 0.495 & 0.231 & 0.219& 0.531 & 0.496  \\
\hline
15 & 1 & 0.174 & 0.184 & 0.609 & 0.562 & 0.232 & 0.232& 0.529 & 0.489  \\
15 & 5 & 0.223 & 0.195 & 0.531 & 0.520 & 0.240 & 0.236& 0.525 & 0.479  \\
15 & 10 &0.224 & 0.205 & 0.524 & 0.462 & 0.246 & 0.249& 0.518 & 0.455  \\
\hline
\end{tabular}
\end{table}

 \begin{figure}[h]
\centering
\includegraphics[width= 0.24\textwidth]{rice.png}

\includegraphics[width= 0.24\textwidth]{svd_rice_3_3_n.png}

\includegraphics[width= 0.24\textwidth]{dct_rice_3_3_x.png}

\includegraphics[width= 0.24\textwidth]{TV_rice_3_3.png}
\vspace*{-3mm}
\caption{Left to right:\ original rice grain image, blurred noisy image (Gaussian blur
with $\sigma = 3$ and noise level 0.03), DCT-EPP reconstruction
(PSNR = 24.3, MSSIM = 0.72), and TV reconstruction (PSNR = 24.9, MSSIM = 0.74).}
\label{fig:rice}
\end{figure}

To illustrate that the EPP and TV reconstructions have different features
(due to the different reconstruction models)
we consider MATLAB's $256 \times 256$
``rice grain'' image shown in Fig.~\ref{fig:rice} together with a Gaussian-blurred version
and the DCT-EPP and TV reconstructions.
The TV reconstruction has sharper edges, which comes at the expense of a more
complicated algorithm with much larger computing time.

\section{Conclusions}
\label{eppsec:con}

We developed a new computational framework for projection-based edge-preserving regularization,
and proved the existence and uniqueness of the solution.
Our algorithm uses standard computational building blocks and is therefore easy to implement
and tune to specific applications.
Our experimental results for image deblurring show that the reconstructions are better than
those from standard projection algorithms, and they are competitive with those from other
edge preserving restoration techniques.

\bibliographystyle{siam}
\bibliography{epp_ref}

\end{document}


TITLE: Intuition behind the implicit function theorem.
QUESTION [0 upvotes]: Is there any intuition behind the implicit function theorem? Take $F(x,y)=0$ where $y=f(x)$. Then,
$$
\frac{dy}{dx}=-\frac{F_x}{F_y}
$$
I see that the derivative is equal to the reciprocal of the partial derivatives, what is the intuition behind this?

REPLY [0 votes]: The intuition for the implicit function theorem as I understand it is: the solutions to $F(x,y)=0$ are a union of graphs of functions.
In your case, it looks like you are talking about one function you've already picked, $y=f(x)$. Sometimes, you might need to have $x=g(y)$ to cover the whole solution set.
One of the typical examples given is $F(x,y)=x^2+y^2 -1$.
Since $dF/dy=2y$ is nonzero whenever you are not on the $x$-axis, you're guaranteed that the points outside of the $x$-axis can be represented as a union of graphs.
In this case, the graph of $f(x)=\sqrt{1-x^2}$ and $g(x)=-\sqrt{1-x^2}$ for $x\in [-1, 1]$ might be one way to express it.
However, some authors require that the domains of these functions be open sets. In that case, it would be impossible to pick up the points of the graph above using functions of $x$. But, thankfully, we could also use the equivalent functions of $y$ to cover those points:
$f'(y)=\sqrt{1-y^2}$ and $g'(y)=-\sqrt{1-y^2}$
They are able to cover everything where $dF/dx$ is nonzero (which is everything off the $y$-axis.)

From your example it looks like you are rather looking for an explanation of why this applies to implicit differentiation, to which I would say: since you can "zoom in" to one of these functions no matter where you are in the solution set, you'll always be able to learn about the tangent (locally) using regular differentiation (since they are functions.)


text:
\newcommand{\spl}{\mathrm{spl}}
\newcommand{\frT}{\mathfrak{T}}
\newcommand{\Tens}{\mathbf{T}}
\renewcommand{\theta}{\vartheta}

\section{Algebraic dynamics arising from automorphism groups of geometric structures}
\label{sec:invariance}

In this section, we explore the dynamics arising from automorphism groups of multilinear
forms on manifolds of finite volume.

\subsection{Preliminaries:  isotropic and split subgroups of algebraic groups}

\subsubsection{Solvable groups}
For a solvable linear algebraic group $\ag{A}$ defined over $\RR$,
let $\ag{A}_\spl$ denote the maximal $\RR$-split connected subgroup of
$\ag{A}$. This means that $\ag{A}_\spl$ is the maximal connected subgroup
trigonalizable over the reals.
For the real algebraic group $A = \ag{A}_{\RR}$, its maximal trigonalizable
subgroup  is $A_{\spl} =  A \cap \ag{A}_\spl$.
Let $\aT$ be a torus  defined over $\RR$. Then $\aT$ is called
\emph{anisotropic}  if $\ag{T}_\spl=\one$.  Otherwise, $\aT$ is called \emph{isotropic}.
Equivalently,  $\aT$ is anisotropic if its group of real points
$T=\aT_{\RR}$ is compact.  Every torus defined over $\RR$ has a
decomposition into subgroups

$ \aT=\aT_\spl\cdot \aT_{\rm c}$,

where $ \aT_{\rm c}$ is a maximal anisotropic torus defined over
$\RR$ and $\aT_\spl \cap \aT_{\rm c}$ is finite.
Moreover, if $\aT \subseteq \ag{A}$ is a maximal torus defined over $\RR$
and $\aU$ is the unipotent radical of $\ag{A}$,
then there is a semidirect product decomposition
\[
\ag{A}_{\spl} = \aU \cdot \aT_{\spl}.
\]
Note also that the split part $\ag{A}_{\spl}$ is preserved under
$\RR$-defined morphisms of algebraic groups.
See Borel \cite[\textsection 15]{borel_book} for more background.

\subsubsection{Semisimple groups}
A simple linear algebraic group $\aL$ defined over $\RR$ is called
\emph{isotropic} if any maximal $\RR$-defined torus is isotropic.
A semisimple linear algebraic group $\aL$ defined over $\RR$ will be called
\emph{isotropic} if all of its simple factors are isotropic.

\subsubsection{The isotropic kernel}
If $\aG$ is a linear algebraic group defined over $\RR$
with solvable radical $\ag{R}$ and a maximal semisimple subgroup $\ag{L}$,
then we define the \emph{isotropic part}  of $\aG$ to be the normal subgroup
\begin{equation}
\aG_{\rm is} = \ag{S}\cdot\ag{R}_{\spl} \, ,
\label{eq:splitpart}
\end{equation}
where  $\ag{S}$ is the maximal normal isotropic subgroup
of $\aL$.  Note that $\aG_{\rm is}$ is a normal subgroup and that there is an almost
semidirect product decomposition
\[
\aG  = \ag{B} \cdot   \aG_{\rm is} \, ,
\]
where  $\ag{B}$ is a reductive anisotropic subgroup.
The normal subgroup $ \aG_{\rm is}$ will be called the \emph{isotropic kernel}
of $\aG$. Any $\RR$-defined homomorphism of $\aG$ to an anisotropic
group contains $\aG_{\rm is}$ in its kernel, that is, it factors over
$\aG/\aG_{\rm is}$.

\begin{lem}\label{lem:split0}
Let $\aG$ and $\aH$ be $\RR$-defined linear algebraic groups and
$\psi:\aG\to\aH$ an $\RR$-defined algebraic homomorphism.
Then $\psi(\aG_{\rm is})\subseteq\aH_{\rm is}$.
\end{lem}
\begin{proof}
The split part of the solvable radical of $\aG$ is mapped to the split part
of the solvable radical of $\aH$ (see Borel \cite[15.4 Theorem]{borel_book}).
Moreover, $\psi$ maps the isotropic simple subgroups of $\aG$
to isotropic simple subgroups of $\aH$.
\end{proof}

\subsubsection{Real algebraic groups}

Let $G = \aG_{\RR}$ be the group of real points of $\aG$.
Then $G_{\rm is} =  G \cap  \aG_{\rm is}$
is called the isotropic kernel  of $G$ and $B = \ag{B} \cap G$ is a compact group. We note
that  $$    G =  B \cdot G_{\rm is}  \;. $$
We conclude that  $G_{\rm is}$ is a normal uniform subgroup of $G$, and, in fact, every morphism of
real algebraic groups from  $G$ to a compact group factors over
$G/G_{\rm is}$.
Recall that a closed subgroup $H$ of $G$ is called \emph{uniform}  if $G/H$ is compact.

\subsection{Dynamics of linear maps, invariant probability measures}

Let  $V$ be a real vector space and $G \subseteq  \GL(V)$.
Let  $A= \ac{G}_{\RR}$ denote the real Zariski closure of $G$.

\subsubsection{Quasi-invariant vectors} \label{sec:quasi-invariant}
For $v \in V$, put
$G_{v} = \{ g \in G \mid g \cdot v = v \}$. Recall that  $v$ is
called \emph{invariant} (by $G$) if $G = G_{v}$.

\begin{definition} \label{def:quasi-invariant}
We call $v \in V$ \emph{quasi-invariant by $G$} if the closure of the orbit $G \cdot v$
is compact. Moreover, $v$ is called  \emph{strongly quasi-invariant} if $\ac{G}_{\RR} \cdot v$ has compact closure in $V$.
\end{definition}

Clearly, for any group $G$, strong quasi-invariance implies quasi-invariance.

If $G \subseteq \GL(V)$ is a real linear algebraic group, then quasi-invariance and
strong quasi-invariance of $v \in V$ are equivalent.

\begin{example} \label{ex:trig}
Let $T  \subseteq \GL(V)$ be a trigonalizable subgroup.
Then $v \in V$ is quasi-invariant  by $T$ if and only if $v$ is
invariant by $T$.
\end{example}

\begin{lem} \label{lem:Gis}
Let $G \subseteq \GL(V)$ be a real linear algebraic group and $v \in V$.
Put $H  = G_{v}$. Suppose that $v$ is quasi-invariant by $G$. Then $H$ contains
the isotropic kernel  $G_{\rm is}$ of $G$. In particular, $H$ is a uniform subgroup of $G$.
\end{lem}
\begin{proof}  By Example \ref{ex:trig},  $H = G_{v}$ must contain every trigonalizable subgroup of $G$.
Therefore,  $H$ contains $G_{\rm is}$.
 \end{proof}

 \subsubsection{Invariant probability measures}

We now consider  $G$-invariant probability measures on $V$.

\begin{example}[Invariant measures are supported on fixed points] \label{ex:trio}
Let $A \subseteq \GL(V)$ be a one-dimensional trigonalizable group and $\mu$ an
$A$-invariant probability measure.   If $\mu$ is ergodic then, by Zimmer \cite[Proposition 2.1.10]{zimmer},
$\mu$ is supported on some orbit $A \cdot v$. Since the Haar-measure on $A$ is
infinite, we must have $A \cdot v = v$. By the ergodic decomposition theorem
\cite[Chapter 5]{viana-oliveira}, any invariant
probability measure $\mu$ may be approximated by finite convex combinations of ergodic invariant
measures, which (as we just remarked) are all supported on the fixed point set $V^{A}$. Therefore,  any $A$-invariant probability measure $\mu$ is supported on  $V^A$.
\end{example}

 The following is a key observation in the context of invariant probability measures
 on representation spaces:

 \begin{prop} \label{prop:measure}
 Let $G \subseteq  \GL(V)$ and $\mu$ a $G$-invariant probability measure on $V$.
 Put  $A = \ac{G}_{\! \RR}$. Then the following hold:
\begin{enumerate}
\item $\mu$ is invariant by $A $.
\item
The action of $A$ on the support of $\mu$ factors over a compact group.

\item
There exists a $G$-invariant subspace
$W \subseteq V$ containing the support of $\mu$,
such that the action of $G$ on $W$ factors over a compact group.
\item
If $A =  A_{\rm is}$ then $\mu$ is supported on $V^A$.
\end{enumerate}
\end{prop}
\begin{proof}
Embed  $V$ as an affine subspace into some real
projective space $\mathrm{P}^n\RR$.
For any probability measure $\nu$ on $\mathrm{P}^n\RR$,
let  $\PGL(n)_{\nu}$ denote the stabilizer of $\nu$  in the
projective linear group. By \cite[Theorem 3.2.4]{zimmer}, $\PGL(n)_{\nu}$
is a real algebraic subgroup of $\PGL(n)$. Since $\mu$ pushes forward to a $G$-invariant measure
on $\mathrm{P}^n\RR$, this implies (1).

To prove (4), observe that  $ A_{\rm is}$ is generated by its one-dimensional split subgroups.
In light of Example \ref{ex:trio}, we deduce that $\supp \mu$ is contained
in the  fixed point set of $A_{\rm is}$.
In particular, $A_{\rm is}$ acts trivially on the support of $\mu$, so (2) follows.
Let $W$ be the smallest $G$-invariant subspace of $V$
that contains $\supp \mu$.
Then $A_{\rm is}$ acts trivially on $W$. Hence (3) holds.
\end{proof}

Note that the action of $G$ on $W$ factors over a compact group and the elements of $W$ are
therefore strongly quasi-invariant by $G$.

\subsection{The characteristic map of a geometric $G$-manifold}

Let $M$ be a differentiable manifold, and $\theta$ an $s$-multilinear
form on  $M$.

The group of diffeo\-morphisms of $M$ that
preserve $\theta$ will be  denoted by $\Aut(M,\theta)$. Let
$$  G\subseteq\Aut(M,\theta)$$ be a Lie group.

We may identify the Lie algebra $\frg$ of $G$ with the subalgebra of
vector fields on $M$ whose flows are one-parameter subgroups of
diffeomorphisms contained in $G$.
We let $L^s(\frg)$ denote the space of $s$-linear forms on $\frg$,
and consider the associated  map
\[
\Phi:M \to L^s(\frg),\quad
p \mapsto \Phi_{p}
\]
where, for $X_1,\ldots,X_s\in\frg$,
the multilinear form  $\Phi_{p}$ on  $\frg$ is declared by
\[
\Phi_{p}(X_1,\ldots,X_s)
= \theta_{p}(X_{1,p},\ldots, X_{s,p}) \; .
\]
The adjoint representation $\Ad$ of $G$ on $\frg$ induces a representation
$$ \rho:G\to\GL(L^s(\frg))   \; .  $$

The map $\Phi$ is equivariant with respect to the
action of $G$ on $M$ and the representation $\rho$.

\begin{thm}\label{thm:isom_is_quasi-invariant}
Let $M$ be a differentiable manifold  and  $G\subseteq\Aut(M,\theta)$.
Suppose that the Lie group $G$ preserves a finite measure on $M$. Then the action
of $G$ on $$ \Phi(M) \, \subseteq  \, L^s(\frg)$$  factors over a compact group.
Moreover, the forms $\Phi_{p}$, $p \in M$, are strongly quasi-invariant by $G$.
\end{thm}
\begin{proof}
Put $V=L^s(\frg)$ and
$A={\ac{\rho(G)}}_{\! \RR}$.
Since $\Phi$ is $G$-equivariant, the finite $G$-invariant measure $\mu$ on $M$ pushes forward
to a finite $G$-invariant measure $\mu'$ on $V$, which is supported on $\Phi(M)$. By Proposition
\ref{prop:measure} (3), there exists an $A$-invariant linear subspace $W$ containing
$\Phi(M)$ on which the action of $A$ factors over a compact group. In particular,
$\Phi_{p}$ is strongly quasi-invariant for all $p \in M$.
\end{proof}

\subsection{Nil-invariant and quasi-invariant bilinear forms on Lie algebras}
\label{sec:nilinvariant}

Let $\frg$ be a Lie algebra and consider its  inner auto\-morphism group  $$G = \Inn(\frg) \, . $$
By definition,
$G$ is the connected Lie subgroup of $\GL(\frg)$ that is tangent to the Lie algebra  of linear
transformations
$$ \ad(\frg) = \{ \ad(x)  \mid x \in \frg \} \; \subseteq \; \gl(\frg) \, . $$
We also consider the real Zariski closure $$\ac{\Inn(\frg)}  \, \subseteq \, \Aut(\frg) \, \subseteq \,  \GL(\frg)$$ of $\Inn(\frg)$, which is an algebraic group of automorphisms of $\frg$ with  Lie algebra $$\ac{\ad(\frg)} \,  \subseteq \der(\frg) \, \subseteq \,  \gl(\frg) \, .$$

We are interested in the  action of $G$ on the space of bilinear forms on $\frg$. The following notions are fundamental (\cf Section \ref{sec:quasi-invariant}):

\begin{definition} \label{def:quasi-invariant_beta}
Let $\beta$ be a  bilinear form on $\frg$. Then
$\beta$ is called \emph{invariant} if $G$ is contained in the orthogonal group  $\OO(\beta)$ of
$\beta$. Likewise,  $\beta$ will be called \emph{quasi-invariant}  if
$$ \ac{\Inn(\frg)} \cap \OO(\beta)   \; \,  \subseteq \; \,  \ac{\Inn(\frg)}$$
is the inclusion of a uniform subgroup.
\end{definition}

According to Lemma \ref{lem:Gis}, $\beta$ is  quasi-invariant  if and only if
\begin{equation} \label{eq:quasi}
 \ac{\Inn(\frg)} \cap \OO(\beta) \, \supseteq  \; \,   \ac{\Inn(\frg)}_{\rm is} \, .
\end{equation}
That is,
the stabilizer of $\beta$ contains the maximal normal subgroup
$$       \ac{\Inn(\frg)}_{\rm is} $$ that is  generated by linear maps in  $\ac{\Inn(\frg)}$
that are trigonalizable over $\RR$.

\subsection{The quasi-invariance condition} \label{sec:quasi_invariance}
The condition \eqref{eq:quasi}  may be formulated in Lie algebra terms  only:

\begin{lem} \label{lem:Gisbeta}
The bilinear form $\beta$ is quasi-invariant by $G$ if and only if the Lie algebra $$ \ac{\ad(\frg)} \cap  \fro(\beta)$$  contains all elements of\/  $\ac{\ad(\frg)}$ that are trigonalizable over $\RR$.
\end{lem}

The Lie algebra  $\ac{\ad(\frg)}_{\rm is}$ generated by all elements trigonalizable over $\RR$ is tangent to
$  \ac{\Inn(\frg)}_{\rm is}$. In particular,
$\ac{\ad(\frg)}_{\rm is}$ is an ideal in  $\ac{\ad(\frg)}$, since
$\ac{\Inn(\frg)}_{\rm is}$ is a normal subgroup of  $\Inn(\frg)_{\rm is}$.
The following  is a noteworthy  consequence:

\begin{lem} \label{lem:acis}
Suppose that $\beta$ is a quasi-invariant bilinear form. Then $$  \ac{\Inn(\frg)}_{\rm is} (\frg)  \subseteq \frg_{\beta}. $$
 That is, for any $\varphi \in  \ac{\ad(\frg)}_{\rm is}$ and $a \in \frg$, we have $\varphi(a) \in
\frg_{\beta}$.
 \end{lem}
 \begin{proof}
 Since  $\ac{\ad(\frg)}_{\rm is}$ is an ideal of $\ac{\ad(\frg)}$, we have
 $$ [\ac{\ad(\frg)}_{\rm is}, \ad(\frg)] \,  \subseteq \,   \ac{\ad(\frg)}_{\rm is} . $$
 But $\ad(\frg)$ is also an ideal
 in the Lie algebra $\der(\frg)$ of all derivations of $\frg$. This implies that, for any
 $\varphi \in \ac{\ad(\frg)}_{\rm is}$,
$$ [\varphi, \ad(a)] = \ad (\varphi(a)) \,  \in  \, \ad(\frg) \cap \ac{\ad(\frg)}_{\rm is} \; . $$
If $\beta$ is  quasi-invariant,
we have  $\ac{\ad(\frg)}_{\rm is} \subseteq \fro(\beta)$. Therefore, for any $\varphi \in \ac{\ad(\frg)}_{\rm is}$,
$\ad(\varphi(a)) \in \ac{\ad(\frg)}_{\rm is} \subseteq  \fro(\beta)$. In particular, $\varphi(a) \in \frg_{\beta}$.
\end{proof}

\subsection{The nil-invariance condition} \label{sec:nil_invariance}

The following condition is somewhat weaker than quasi-invariance since it takes into account only the nilpotent elements in the endomorphism Lie algebra $\ac{\ad(\frg)}$:

\begin{definition} \label{def:nilinvariant_beta}
A bilinear form $\beta$ on $\frg$ is called \emph{nil-invariant} if
$$     \varphi \in   \fro(\beta)   $$
for all nilpotent elements $\varphi$ of the Lie algebra  $\ac{\ad(\frg)}$.
\end{definition}

We fix  a Levi decomposition
\[
\frg = (\frk\times\frs)\ltimes\frr ,
\]
where $\frk$ is semisimple of compact type, $\frs$ is semisimple without
factors of compact type, and $\frr$ is the solvable radical of $\frg$.

\begin{lem} If $\beta$ is nil-invariant, then
\begin{equation} \label{eq:nilinvariant}
	\frs + \frn \, \subseteq \, \frg_\beta \, .
\end{equation}
\end{lem}
To see that $\frs$ is contained in $\frg_\beta$ recall that the subalgebra $\frs$ is generated by the set of elements $x \in \frs$ such that $\ad_\frg(x)$ is nilpotent (for example \cite[Lemma 5.1]{BGZ}).

\subsection{Quasi-invariant and nil-invariant symmetric bilinear forms}

In the following, we recall important structural properties of
Lie algebras with nil-invariant symmetric bilinear form.

\begin{thm}[\mbox{\cite[Theorem A, Corollary C]{BGZ}}]\label{thm:BGZ_invariance}
Let $\frg$ be a finite-dimensional real Lie algebra with
nil-invariant symmetric bilinear form $\met$.
Let $\met_{\gs}$ denote the restriction of $\met$ to $\gs = \frs \ltimes \frr$.
Then:
\begin{enumerate}
\item
$\met_{\gs}$ is invariant by the adjoint action of $\frg$ on $\gs$.
\item
$\met$ is invariant by the adjoint action of $\gs$.
\end{enumerate}
Furthermore, if $(\frg, \met)$ is effective,

then
\begin{enumerate}
\item[(3)] $ \frg^{\perp} \, \subseteq \, \frk \ltimes \frz(\gs)$  and $\; [\frg^{\perp}, \gs] \, \subseteq \, \frz(\gs) \cap \frg^{\perp}$. If  $\frg^\perp \subseteq \frg_{\met}$ then $[\frg^\perp , \gs] = \zsp$.
\end{enumerate}
\end{thm}

By definition, quasi-invariance always implies nil-invariance. It is somewhat surprising and non-trivial that for a \emph{symmetric} bilinear form also the converse holds and an even stronger implication is satisfied in the solvable case:

\begin{cor}  \label{cor:nil_and_quasi}
Let $\beta$ be a symmetric bilinear form on $\frg$. Then $\beta$ is quasi-invariant if and only if $\beta$ is nil-invariant.  Moreover, if $\frg$ is solvable then any nil-invariant symmetric bilinear form $\beta$ is
invariant.
\end{cor}

Corollary  \ref{cor:nil_and_quasi} is directly implied by Theorem \ref{thm:BGZ_invariance}, which streng\-thens
\eqref{eq:nilinvariant} for symmetric bilinear forms considerably.

\subsection{Nil-invariant and quasi-invariant skew-symmetric bilinear forms}

As noted above for symmetric bilinear forms on Lie algebras, nil-invariance already implies quasi-invariance.
The following simple example illustrates
that, in general,  this is not the case for skew-symmetric bilinear forms:

\begin{example} \label{ex:nil_noquasi}
Let $\frg$ be the two-dimensional solvable Lie algebra with generators
$x,y$, satisfying $[x,y]=y$.
Consider the skew-symmetric form $$\omega = \d x \wedge \d y$$  on $\frg$
(which satisfies

$\omega(x,y)=1$).

It is easily verified that $\omega$ is invariant by $\ad(y)$, but not by
$\ad(x)$.
This means $\omega$ is nil-invariant without being quasi-invariant.
\end{example}

In the following Section \ref{sec:skew_Lie}, we will develop the structure theory of quasi-invariant and nil-invariant skew-symmetric forms in more detail.


TITLE: Evaluating integral involving Bessel function.
QUESTION [1 upvotes]: Evaluate
$$\int_0^{\infty } \frac{2^{\frac{r}{\delta }} \left(2^{\frac{r}{\delta }}-1\right) r\ e^{-\frac{\alpha ^2+\left(2^{\frac{r}{\delta }}-1\right)^2}{2 \beta ^2}}\log 2 }{\beta ^2 \delta } I_0\left(\frac{\alpha  \left(2^{\frac{r}{\delta }}-1\right)}{\beta ^2}\right) \, dr$$
I arrived at the above equation while trying to obtain the Expectation of $r$.
I'm at my wits end with this problem. I don't know how to continue so either some hints, pushes in the right direction or the full evaluation would be highly appreciated. Thanks

REPLY [0 votes]: What is $I_0$? A constant? Without this, make a change of variables
$$r\to r\ln 2/\delta$$to obtain up to a multiplicative constant
$$\int_0^{\infty }   e^{r} \left(e^{ {r} }-1\right) r\ e^{-\frac{\alpha ^2+\left(e^{ {r} }-1\right)^2}{2 \beta ^2}}   I_0\left(\frac{\alpha  \left(e^{ {r} }-1\right)}{\beta ^2}\right)  dr$$
Now make a change of variables $u=e^r-1$, $du = e^rdr$ to obtain (again, up to a multiplicative constant)
$$\int_0^{\infty }    u \ln(u+1)\ e^{-\frac{\alpha ^2+u^2}{2 \beta ^2}}   I_0\left(\frac{\alpha  u}{\beta ^2}\right)   du,$$
which seems to be more or less doable depending on the form of $I_0$.


TITLE: Help understanding centre of mass
QUESTION [2 upvotes]: So basically the question I posed myself a few days ago was for a torus shaped planet, what would the gravitational field be like inside the torus? After viewing some answers on the internet, (Gravity on a doughnut-shaped/Möbius planet, and also having my friend code http://ideone.com/s9tlyM for me to check), it seems that anywhere inside the torus, you would experience a force towards the .centre (ie the centre of mass). TORUS RING misinterpreted results oops
So what does the centre of mass really mean for an object, if when I'm in the torus planet, and I'll fly to the edges, instead of being attracted to the centre of mass (which seems logical to me)
Or does centre of mass have no relation to gravity at all?
http://en.wikipedia.org/wiki/Barycentric_coordinates_(astronomy)
As best as I understand this article, it says that a system of stars will orbit around their centre of mass, so if you treat me and the torus planet as two stars shouldn't I be attracted to the centre of mass?
Thanks for any help clearing this up.

REPLY [2 votes]: When considering the gravitational field at a point, you need to sum the effects of all mass based on the quantity, direction, and distance of those masses.   In the general case, there may be no simplification for this summation.
But there are several cases where simplifications are possible.  In particular, spheres and shells with radial symmetry can have their contribution replaced by a point mass of equal mass at their center of gravity for regions outside their surface.  Because many real bodies such as stars and planets are a very close fit for this distribution, center of mass is used quite often.  The barycenter article that you linked to is using the assumption that the masses involved are radially symmetric.
But as you've seen in your torus, it cannot be used as a simplistic replacement in all situations.  Even for a spherical shell of mass, objects inside the shell are not attracted to the center.