text:
\begin{document}
 \begin{abstract} We propose and analyze a mathematical model of Modulated Luminescent Tomography. We show that when single X-rays or  focused X-rays are used as an excitation, the problem is similar  to the inversion of  weighted X-ray transforms. In particular, we give an explicit inversion in the case of Dual Cone X-ray excitation.
 \end{abstract}

 \maketitle
\section{Introduction: A brief description of the method and the model}
\subsection{ A brief description of the modalities} \
The purpose of this work is to build and analyze the mathematical model of several medical multiwave imaging techniques that we call Modulated Luminescent Tomography. Certain phosphorescent contrast agents (for example, nanophosphors) are delivered to the cells and then illuminated with    X-rays. The illuminated particles emit photons which scatter through the tissue and are detected outside the body, see, e.g., \cite{Liu2013}. The goal is to recover the concentration of the contrast agents. The resolution is expected to come from the X-ray excitation. In this respect, those modalities behave differently than Photoacoustic and Thermoacoustic Tomography, where the excitation is highly diffusive but the emitted ultrasound signal allows for high resolution, see, e.g., \cite{Wang2012}.

Two such techniques, X-ray luminescence computed tomography (XLCT),  \cite{Pratx10} and X-ray micro-modulated luminescence tomography (XMLT) \cite{Cong14, St-double_cone_eng}, described below, have been proposed recently.

Microscopy is the principal observational tool that has made fundamental contributions to our understanding of biological systems and engineered tissues. Popular microscopy techniques use visible light and electrons. Sample preparation and imaging with these techniques are relatively simple, being good for in situ or in vivo studies of cultured cell/tissue samples. Inherently, the resolution of optical microscopy is limited by diffraction. With additional sample preparation, stochastic information and innovative interference techniques ~100nm resolution is achievable. Three-dimensional images can be obtained with optical sectioning. Ultimately, multiple scattering prevents these techniques from imaging thick samples. Photoacoustic tomography permits scalable resolution with a depth of up to ~7cm with a depth-to-resolution ratio ~200. Photoacoustic microscopy aims at millimeter depth and micron-scale resolution based on absorption contrast \cite{Wang2013}, which can be used to characterize the structure of a scaffold but does not provide the sensitivity of fluorescence imaging. All these methods are bounded by 1mm imaging depth.

In the in vivo imaging field, fluorescence molecular tomography (FMT) and bioluminescence tomography (BLT) are capable of visualizing biological processes at molecular and cellular levels deep inside living tissue \cite{Ntziachristos2005}. Such optical molecular imaging tools have major advantages in terms of high sensitivity, biological specificity, non-invasiveness, and cost-effectiveness, which are widely used in preclinical investigations and a limited number of clinical applications \cite{Wang2006}. However, resolution and stability of FMT and BLT still remain unsatisfactory. Recently, X-ray luminescence computed tomography (XLCT) was developed \cite{Pratx10} with nanophosphors (NPs) as imaging probes. Excited with a pencil beam of X-rays, NPs simultaneously gives rise to luminescence emission whose photons can be efficiently collected with a CCD camera. This mechanism allows tomographic reconstruction of a NP distribution, similar to the case of X-ray CT. X-ray excitation not only confines luminescence sources within the active beam but also eliminates auto-fluorescence related artifacts in FMT images.  However, due to several technical and physical limitations, the lower bound of XLCT resolution is about one millimeter and seems infeasible to break through.

Recently,  a unique imaging approach called “X-ray micro-modulated luminescence tomography (XMLT)” approach was proposed \cite{Cong14, St-double_cone_eng}, which combines X-ray focusing, nanophosphor excitation, optical sensing, and image reconstruction in a synergistic fashion, and promises significant imaging performance gains in both microscopic and preclinical studies. The XMLT approach is
an exogenous reporter-based imaging system in which spatial resolution is determined by the excitation of the luminescence from nanophosphors with a micro-focused X-ray beam such as a zone-plane or a micro-focus X-ray source coupled with a polycapillary lens. High in vivo spatial resolution can be achieved due to the short wavelength and small spot size of the X-rays.  When X-rays are focused in this way, double cones of X-rays are formed with their shared vertex point inside a sample or a subject.

Our main results state that we can use boundary averaged measurements only and the problem is reduced to an inversion of a weighted linear transform determined by the excitation. For example, if the excitation consist of single X-rays, in Section~\ref{sec5}, we get a weighted X-ray transform with a weight depending on the diffusion and the absorption coefficients $D$ and $\mu_a$ of the medium with respect to the emitted photons, see \r{M1} below.
In the case of double cone excitation (XMLT) with a constant aperture, or in XLCT, we have an explicit and stable inversion under an explicit  if and only if condition \r{stab_cond}, see Theorem~\ref{thm4.1} and Section~\ref{sec4}. We show that the local problem in a region of interest behaves similarly to the corresponding X-ray problems, and we have explicit microlocal inversions of the visible singularities.

\subsection{The model}

Let $\Omega\subset \R^n$ be a smooth domain.
We illuminate the medium in different ways. Each illumination excites the particles at a rate $I_\alpha(x)$, which can be a distribution, where $\alpha$ belongs to some index set $\mathcal{A}$, discrete or not. In the X-ray case, $I_\alpha(x)$ is just a superposition of delta functions along straight lines; which also allows for  attenuated X-rays by introducing an appropriate weight. The structure of $I_\alpha$, XLCT or XMLT or something else, see section~\ref{sec4} and section~\ref{sec5}, is not important at the moment.
The modulated source then is
\be{S}
s_\alpha(x) = I_\alpha(x)f(x),
\ee
where $f$ is the concentration of the active particles, which we want to reconstruct.   We assume that the light from the source propagates according to the diffusion model \cite{Welch-book}
\be{M1}
\left(-\nabla\cdot D\nabla +\mu_a\right) u_\alpha= s_\alpha (x) \quad \text{in $\Omega$}.
\ee
Here, $\mu_a(x)\ge0$ is the absorption, $D(x)>0$ is the diffusion coefficient,
and $u_\alpha$ is the photon density. The coefficients $\mu_a$ and $D$ are given. In the diffusion approximation regime, $D(x)= [3(\mu_a+\mu'_s)]^{-1}>0$, where   $\mu_s'=(1-g)\mu_s$ is the reduced scattering coefficient, $\mu_s$ is the scattering coefficient and $g$ is the scattering anisotropy.

The boundary conditions are of semi-transparent type due to the different index of refraction of the tissue and the air around it. They are of Robin type and have the form
\be{M2}
u_\alpha+2AD\partial_\nu u_\alpha \big|_{\bo}=0,
\ee
where $A>0$ is a given coefficient, well approximated by $A=(1+R)/(1-R)$ with $R$ closely approximated by $R= -1.4399m^{-2} + 0.7099m^{-1}+ 0.6681+0.063m$ \cite{Schweiger95}, where $m$ is the  refractive index   of the tissue.
Let $G$ be the solution operator (the Green's function) of \r{M1}, \r{M2}, i.e., $u_\alpha = Gs_\alpha = GI_\alpha f$.

What we measure is the outgoing photon density
\be{M3}
Q_\alpha f= -D\partial_\nu u_\alpha \big|_{\bo} =\frac1{2A} u_\alpha |_{\bo} =\frac1{2A} (GI_\alpha f)|_{\bo}.
\ee

The data therefore is
\[
\{Q_\alpha f\}_{\alpha\in\mathcal{A}},
\]
and we want to reconstruct $f$.

For each illumination choice, $\alpha$, the data $Q_\alpha$ is a function of $n-1$ variables. When $\alpha$ runs over a continuous space with dimension $m$, we get $m+n-1$ variables. This could make the problem formally overdetermined. On the other hand, because of the diffusion nature of the data for a fixed $\alpha$, we cannot have much resolution hidden in that $Q_\alpha$. The only way to get resolution (i.e., stability) is to have a well chosen set of illuminations $I_\alpha$,  $\alpha\in\mathcal{A}$ which have enough singularities. We show below that we do not actually need to know or measure $Q_\alpha$ pointwise; some average over $\bo$ is enough for a stable recovery. This removes $n-1$ variables from the data.

\section{The first phase: recovery of an averaged intensity.}
The operator $L:= -\nabla\cdot D\nabla +\mu_a$ is symmetric and positive on smooth functions satisfying the the Robin boundary condition \r{M2} because
\be{L}
\begin{split}
\int_\Omega (Lu)\bar u\,\d x &= \int_\Omega\left(  |\nabla u|^2+\mu_a|u|^2 \right)\d x - \int_\bo (\partial_\nu u)\bar u\,\d\sigma\\
 &=  \int_\Omega\left(  |\nabla u|^2+\mu_a|u|^2 \right)\d x +\int_\bo \frac{1}{2AD}  |\bar u|^2\,\d\sigma,
\end{split}
\ee
where $\d \sigma$ is the surface measure.  By
Green's formula, for such $u$ and $v$,
\be{M4}
\begin{split}
\int_\Omega \left( vLu -uLv\right)\,\d x &= -\int_{\bo} \left(  v\partial_\nu u - u\partial_\nu v    \right)D\,\d \sigma\\
 &= -\int_{\bo} \big(  (v +2AD\partial_\nu v) \partial_\nu u - (u  +2AD\partial_\nu u )\partial_\nu v    \big)D\,\d \sigma    =0.
\end{split}
\ee
Therefore, $L$ with the Robin boundary conditions if symmetric and positive by \r{L}.
It is well known that $L$ extends to a self-adjoint operator on $L^2(\Omega)$  with a compact resolvent. In particular, $0$ is in the resolvent set by \r{L}. This shows that the Green's function $G$ is well-defined, and self-adjoint on $L^2(\Omega)$. Moreover, the boundary-value problem \r{M6} below is well posed.

If $Lv_\alpha=0$, and $u_\alpha$ solves  \r{M1}, \r{M2},   we get (we drop the index $\alpha$ in this formula)
\be{M5}
\begin{split}
\int_\Omega   vs\,\d x &= -\int_{\bo} \left(  v\partial_\nu u - u\partial_\nu v     \right)D\,\d \sigma= -\int_\bo \left(   v +2AD \partial_\nu v    \right)D\partial_\nu u \,\d \sigma \\
&= \int_{\bo} \left(   v +2AD \partial_\nu v    \right)Qf \,\d \sigma.
\end{split}
\ee
This suggests the following. Choose $h_\alpha$ and let $v_\alpha:=Vh_\alpha$ solve
\be{M6}
\begin{split}
\left(-\nabla\cdot D\nabla +\mu_a\right) v_\alpha&= 0 \quad \text{in $\Omega$},\\
 v_\alpha+2AD\partial_\nu v_\alpha\big|_{\bo}&=h_\alpha.
\end{split}
\ee
 Then, by \r{M5},
\be{M7}
\int_\Omega v_\alpha  I_\alpha f  \,\d x = \int_{\bo}h_\alpha Q_\alpha f\,\d \sigma.
\ee
In other words, choose $h_\alpha$ somehow; then we can recover $\int_\Omega\! (Vh_\alpha)I_\alpha f\,\d x$ from the data. This requires us to solve \r{M6} first but this can be done numerically in a very efficient way. If $D$ and $\mu_a$ are constants, and $\Omega$ is a circle or a rectangle, it can be done explicitly, as we show below. This gives us a way to compute an averaged value of $I_\alpha f$, with a weight $v_\alpha$ depending on the choice of $h_\alpha $. If   $\mu_a=0$, one can take $h_\alpha=1$ and then $v_\alpha =1$.
 We would want to have positive solutions $v>0$, and we do not want $v\ll1$ because we will divide by it eventually. Positivity can be guaranteed by the maximum principle if $h>0$, since  $0<v<\max_{\bo}h$.
The lower bound of $u$ could be very small if the attenuation $\mu_a$ is large, and that would create stability problems; but this is natural to expect.

Note that we do not need to solve the boundary value problem \r{M6} to find $v$. We could just take any solution $v$ of $(-\nabla\cdot D\nabla +\mu_a)v= 0$ and compute $h=  (v+2AD\partial_\nu v)\big|_{\bo}$ next. If $D$ and $\mu_0$ are constants, we can generate many explicit solutions, regardless of the shape of the domain.

Finally, we want to emphasize that we cannot recover $s_\alpha= I_\alpha f$ from \r{M7} known for all $h_\alpha$; for example,  for any $\phi\in C_0^\infty(\Omega)$, adding $L\phi$ to $I_\alpha f$ would not change the left-hand side of \r{M7}.

We consider the constant coefficient cases below and continue with the general case in the next section.

\subsection{Partial cases: $D$ and $\mu_a$ constants}
Assume now that $D>0$ and $\mu_a\ge0$ are constants.

\subsection{$n=3$, $\Omega$ a ball}
Let $n=3$. The Green's function of $L=-D\Delta+\mu_a$ is then given by
\be{MG}
G_k(x,y) = \frac{e^{-k|x-y|}}{4\pi D|x-y|},
\ee
with $k=\sqrt{\mu_a/D}$, i.e., $(-D\Delta_x+\mu_a) G(x,y)=\delta(x-y)$
Then $G_{-k}$ is a Green's function as well, and the difference is in the kernel of $L$. Therefore,
\[
(-D\Delta_x+\mu_a)\frac{\sinh\left(k|x-y|\right)}{|x-y|}=0.
\]
Therefore, we can chose
\be{Mv}
v(x) = \frac{\sinh\left(k|x-y|\right)}{|x-y|}
\ee
with $y$ arbitrary but fixed. Then compute $h$ from \r{M6}.

If $\Omega$ is the ball $\Omega= \{x|\;|x|\le a\}$, then we can choose $y=0$ to get
\be{M8a}
v(x) = \frac{\sinh{(k|x|)}}{|x|}, \quad
h=  a^{-1}\sinh(ka)+     2AD\left(ka^{-1}\sinh(ka)- a^{-2}\sinh(ka)\right).
\ee
Note that $h$ is just a constant on the boundary $|x|=a$.
Then by \r{M7}, we can recover
\[
\int_{|x|\le a}I_\alpha(x)f(x)\frac{\sinh{(k|x|)}}{|x|}\, \d x.
\]
When $ka\gg1$, the function $v$ is very small near $0$ compared to its boundary values, which means poor recovery of $s_\alpha$ there (we do not actually recover $s_\alpha$ here, just a weighed integral of it). That can be expected though --- signals coming from the center would  be attenuated the most.

\subsection{$n=3$, $\Omega$ arbitrary} One can still use $v$ as in \r{M8a} (or, as in\r{Mv}, see the remark above); then $h$ needs to be calculated  by the second equation in \r{M6}.

\subsection{$n=2$, $\Omega$ a disk} In this case, the fundamental solution involves the Bessel function $K_0$ and is $K_0(k|x-y|)$. Note that there is a weak (logarithmic) singularity at $x=y$. This is equivalent to \r{MG} in the 3D case.  One can then take
\[
v(x) = I_0(k|x|),
\]
where $I_0$ is a Bessel function again, and there is no singularity. This is equivalent to taking $\sinh(k|x|)/|x|$ in 3D. For $h$ on the sphere $|x|=a$ the boundary of the ball), we get
\[
h= I_0(ak)+ 2ADI_1(ak),
\]
which is a constant again.

\section{Reducing the problem to the invertibility of a weighted version of the illumination transform}\label{sec3}
As we showed above, we cannot recover the source $s_\alpha =I_\alpha f$ from \r{M7}. We are not trying to reconstruct $s_\alpha$ however.
 We view \r{M7} as a linear integral transform
\be{M8}

R : f \longmapsto Rf(\alpha) := \int_\Omega v_\alpha I_\alpha f\, \d x, \quad \alpha\in \mathcal{A}.
\ee
We can write equation \r{M7} then as
\be{Q1}
Rf(\alpha) = \int_\bo  h_\alpha Q_\alpha f\,\d\sigma.
\ee
The right-hand side is determined by the data.
The goal then is to invert $R$. Recall that we have some freedom to choose $v_\alpha=Vh_\alpha$, and in particular, we can choose them independently of $\alpha$. Then $v$ just multiples $f$, and we need to invert the transform with kernel $I_\alpha(x)$; and then divide by $v$.

We consider two main examples below: when $R$ is a dual cone transform modeling the XMLT; and when $R$ is a weighted X-ray transform, which models the XLCT.

\section{Double cone excitation (XMLT)}\label{sec4}
\subsection{Formulation}
 We work in $\R^n$ but the interesting case is $\R^3$.
The idea of the XMLT is to focus X-rays, thus forming a double cone,

at each point $x$ in some region of interest ROI at $N$ different directions $\theta_j$, $j=1,\dots, N$. Then
\[
\alpha=(x,j)\in \Omega\times \{1,2,\dots,N\}.
\]
We model the cones with their aperture functions $a_{x,j}(\theta)$, $\theta\in S^{n-1}$, in other words, $I_\alpha(x)$ is a superposition of X-rays with density $a_{x,j}$.
 Then the corresponding  intensity $I_{x,\alpha}(y)$ is given by
\[
I_{x,j}(y) = \frac{  a_{x,j}\big(\frac{x-y}{|x-y|}\big)}{|x-y|^{n-1}}.
\]
For example, if we assume uniform angular density, then $a_{x,j}$ would be constant on the intersection of the interior of the cone with the unit sphere, and zero otherwise. We assume that the cones are double; then $a_{x,j}$ are even functions on the sphere. The simplest case is when $a_{x,j}$ is independent of $x$, and $a_{x_1,j_1}$ is obtained from $a_{x_2,j_2}$ by a translation in the $x$ variable and a rotation in the angular one.

By \r{S}, the transform $R$ becomes
\be{X2a}
Rf(x,j) =\int_\Omega\frac{  a_{x,j}\big(\frac{x-y}{|x-y|}\big)}{|x-y|^{n-1}}v_{j,x}(y)f(y)\,\d y.
\ee
If $a_{x,j}$  are smooth functions, then $R$ is a \PDO\ of order $-1$ with principal symbol
\[
r_j(x,\xi) :=
\pi \int_{S^{n-1}}a_{x,j}(\theta)   v_{j,x}(x)  \delta (\xi\cdot\theta)\,\d\theta,
\]
see, e.g., \cite{FSU}. Here, $\delta$ is the Dirac delta ``function''; therefore, the integral is taken over the sphere of co-dimension two (the grand circle if $n=3$, two points if $n=2$) of the unit sphere intersected with the plane normal to $\xi$.
Note that the symbol above is homogeneous of order $-1$. This shows that the ``visible'' set of singularities $\mathcal{U}\subset T^*\Omega\setminus 0 $ is where the system $\{r_j\}_{j=1}^N$ is elliptic, i.e.,
\be{U}
\mathcal{U} = \{  (x,\xi)\in T^*\Omega\setminus 0;\;  \exists j,   \theta\perp\xi\; \text{so that}\; a_{x,j}(\theta)\not=0   \}.
\ee
\subsection{Recovery of the singularities in a region of interest} If $\Omega_0\subset\Omega$ is an open subset (a region of interest), and if $ a_{x,j}$ are smooth, we can recover the singularities of $f$ there, if the condition for $\mathcal{U}$ in \r{U} is satisfied for any $x\in\Omega_0$, i.e., if
\be{stab_cond}
T^*\Omega_0\setminus 0\subset\mathcal{U}.
\ee
 This can be done constructively as follows.
\begin{itemize}
\item For any $j$ and $x\in \Omega_0$,   choose a smooth function $0<h_{j,x}$ on $\bo$ and compute $v_{j,x}=Vh_{j,x}$ by solving \r{M6} with $h_\alpha =h_{j,x}$.  This determines the transform $R$.
\item Compute the right-hand side of \r{Q1}
 for $x\in\Omega_0$.
\item In \r{Q1}, apply a left parametrix for the matrix operator $f\mapsto Rf(j,\cdot)$ in $\Omega_0$.
\end{itemize}
One way to construct a parametrix of order $1$ is to construct a \PDO\ with the following principal symbol
\[
q_j := \Big (\sum_kr_k^2\Big)^{-1}   r_j,
\]
and then $\sum_j  q_j(x,D)Rf(j,x)=f+Kf$, where $K$ is of order $-1$. If $a_{x,j}\ge0$, one can use instead
\[
 q := \Big (\sum_k r_k\Big)^{-1}
\]
and then $ q(x,D)\sum_j Rf(j,x)=f+Kf$ with $K$ as above.

Note that we can simply choose $h_{j,x}$ independent of ${j,x}$.

The aperture functions $a_{x,j}$ only need to be smooth of finite order to recover the leading singularities, like jumps across surfaces. In practical implementations, this condition can be satisfied approximately by introducing  vignetting near the edges of the cones.  Even without applying the operator $Q$, the $\sum_j Rf(j,x)$ would recover all singularities in the correct places but with varying amplitudes.

When the stability condition \r{stab_cond} does not hold, we can still reconstruct the visible singularities in $\mathcal{U}$ stably  by a microlocal inversion, if $a_{x,j}$ are smooth but  not the ones in the complement of $\bar{\mathcal{U}}$, see, e.g., \cite{SU-JFA09}.

\subsection{Explicit global inversion when the apertures are translation-independent}
Assume now that the aperture functions are non-negative and independent of $x$, i.e., they are given by $a_j(\theta)$.  The stability condition \r{stab_cond} then takes the following form, see \r{U} :
\be{U2}
\text{Any plane through the origin intersects at least one of the cones $\{\theta;\; a_j(\theta)>0\}$.}
\ee
Assume also that $h_{j,x}$ are independent of $j$ and $x$.
The construction above then simplifies as follows. We have
\[
r_j(x,\xi) :=   r^0_j(\xi)v(x), \quad r^0_j(\xi):= \pi
\int_{S^{n-1}}a_j(\theta)    \delta (\xi\cdot\theta)\,\d\theta,
\]
and $r_j(\xi)v(y)$ is actually the full (not just the principal symbol) amplitude of $f\mapsto Rf(j,\cdot)$. Symbols independent of $x$ are  Fourier multipliers, i.e., $r_j(D)= \mathcal{F}^{-1}r_j\mathcal{F}$, where $\mathcal{F}$ is the Fourier transform.
Then $R_j= r^0_j(D)v$ with $v$ regarded as a multiplication, where we changed the notation for $R$ since $\alpha=j$ now. In integral form,
\[
R_jf(x) = \int_\Omega\frac{  a_j\big(\frac{x-y}{|x-y|}\big)}{|x-y|^{n-1}}v(y)f(y)\,\d y.
\]
Note that $r_j^0(D)$ makes sense even if $a_j$ is $L^\infty$ only.

Equation \r{Q1} takes the form
\[
R_jf (x)= \int_\bo h Q_{j,x}\,\d\sigma.
\]
When the stability condition \r{U2} holds,   this equation can be solved explicitly as follows
\be{f}
f = \frac1v \Big(\sum_j r_j(D)  \Big)^{-1}\sum_j \int_\bo h Q_{j,x}\,\d\sigma.
\ee
We summarize this in the following:

\begin{theorem}\label{thm4.1}
Let the aperture functions $0  \le  a_j(\theta)\in L^\infty$ be independent on the point $x$ where we focus.  Choose $0<h\in L^\infty(\bo)$ and let $v=Vh$ be the solution of \r{M6} corresponding to that $h$. If the stability condition \r{U2} holds, then $f$ can be explicitly computed by \r{f}.

Moreover,
\[
\|f\|_{L^2(\Omega)} \le C\frac{\sup_{\bo}h}{\inf_\Omega v}\Big\| \int_\bo \sum_j Q_{j,x}\,\d\sigma  \Big\|_{H^1(\Omega_x)}.
\]
\end{theorem}

The subscript $x$ in $\Omega_x$ indicates that the $H^1$ norm is taken w.r.t.\ to the variable $x$.

\begin{remark}\label{rem1}
One way to satisfy the stability condition \r{U2} is to make sure that union of the open sets $\{a_j>0\}$  where the apertures do not vanish covers the equator of the unit sphere (or any fixed in advance grand circle). This is the situation in the numerical example below. If the apertures $\{a_j>0\}$ are too small, then there will be a large variation between the minimum and the maximum of the symbol $q$ below, which may lead to a mild instability; the most stable singularities would be the vertical ones.
\end{remark}

\section{X-ray excitation (XLCT)} \label{sec5}
In the X-ray excitation case, we send individual  X-rays through the medium. If the direction is $\theta\in S^{n-1}$, and we parameterize lines in this direction by initial points $z$ on $ \theta^\perp =\{x;\; x\cdot\theta=0\}$, we have
\[
\alpha = (\theta,z), \quad \theta\in S^{n-1}, \; z\in \theta^\perp.
\]
In other words, we identify the directed lines in $\R^n$ with such $\alpha$'s. Then $I_\alpha$ is just a delta unction on the line parameterized by $\alpha$. Therefore, $R$ is just the weighted X-ray transform with weight $h_{\theta,z}$. Assume now that $h$ is chosen independently of the line parameterized by $(\theta,z)$. Then
\[
Rf= R_0vf, \quad R_0g(\theta,z) = \int g(z+t\theta)\,\d t.
\]
Equation \r{Q1} takes the form
\[
(R_0vf)(\theta,z)= \int_\bo hQ_{\theta,z}\,\d\sigma.
\]
In other words, each measurement gives us an integral of $vf$ over a single line.
The function $f$ can then be found by inverting the X-ray transform $R_0$ and then dividing by $v$. The problem with a subset of lines and microlocal recovery in a ROI is the same as with the X-ray transform.
Note that in practice, diffraction and other engineering challenges limit the ability to concentrate the radiation too closely along a single line thus limiting the resolution.

\section{Numerical simulations}

A spherical phantom of radius 10 mm was employed for the numerical simulation. The phantom is assigned with optical parameters: absorption coefficient $\mu_a=0.05$ mm$^{-1}$, scattering coefficient $\mu_s=15.$0 mm$^{-1}$, anisotropic coefficient $ g=0.9$, and relative refractive index of $1.37$. The phantom was discretized into $65,775$ tetrahedral elements and $12,044$ nodes. A total of $2,108$ virtual detectors were distributed over the phantom surface to record the photon fluence rates. Then two light sources of radius $1.0$ mm were embedded into the phantom, and filled with nanophosphors of concentrations of $5$ $\mu$g/mL and $10$ $\mu$g/mL, respectively. The centers of the two light sources were at $(2.5, 2.5, 0.0)$ and $(3.5, 0.0, 0.0)$. Two spherical subregions are nearly connected to test the spatial resolution. Using a polycapillary lens, X-ray beams are reshaped to double cone with a cone angle of $19.2$ degrees. The focal point of the X-ray is scanned at grid points in a region of interest  $\text{ROI} = \{ -2<x<5,\; -2<y<5,\; -2<z<2   \} $  along a scanning direction.  Then, the phantom was rotated with $10$ equal angles over a $360$ degrees range to acquire sufficient information for an improved stability of the reconstruction.
The stability condition \r{U2} then holds, see Remark~\ref{rem1}.

At each scanning, the intensity of the photon luminescence on the surface of the phantom was acquired. The intensity describes the optical emission of nanophosphors with double cone excitation.
 Poisson noise was added to the synthetic data for the simulation of measurements. After discretization,  the LSQR method \cite{Paige:1982} was used to solve solve the resulting system, in order to reconstruct the nanophosphor distribution from simulation data. The reconstructed results are in excellent agreement with the true phantom, and the average relative error of the reconstructed nanophosphor concentration was less than $5.37\%$, which was defined as
 \[
 \text{Error} = \frac1{ \#\{i|\; \rho_i^T>\epsilon  \}}   \sum_{\rho_i^T>\epsilon} \frac{\rho_i^r-\rho_i^T }{  \rho_i^T}    ,
  \]
where $ \rho_i^T $ and $\rho_i^r$ are the true and reconstructed nanoparticles concentrations, respectively,  and $\epsilon$ is a background noise level.  Figures~\ref{pic1}--\ref{pic2} presents the comparison between the true and reconstructed nanophosphor distribution, showing the quantification accuracy of the image reconstructions.

\begin{figure}[!ht]
  \centering
  \includegraphics[width = 2in]{Xray_RPI_picA1} \ \ \ \ \ \ \ \ \ \
  \includegraphics[width = 2.0in]{Xray_RPI_picA2}

  \caption{Modulated luminescence tomography simulation. Left: the true nanophosphor distribution in the phantom. Right: the reconstructed nanophosphor distribution.
 The different colors here do not represent different values.  }
  \label{pic1}
\end{figure}

\begin{figure}[!ht]
  \centering
  \includegraphics[width = 2.8in]{Xray_RPI_picB1} \ \
  \includegraphics[width = 2.8in]{Xray_RPI_picB2}

  \caption{Modulated luminescence tomography simulation. Left: The true nanophosphor distribution on the 2D slice  $z=0$. Right:  the reconstructed nanophosphor distributionat on the 2D slice $z=0$.  }
  \label{pic2}
\end{figure}

\bibliographystyle{abbrv}
\bibliography{myreferences,Plamenbiomed}
\end{document}


TITLE: Arbitrary entangled matter-photon emitter
QUESTION [0 upvotes]: I am recently approaching to photonics and integration between matter qubits together with photons. I have an interest in understanding the assumption I can do when abstracting such technologies.
I read in this paper about the generation of an EPR state between a matter qubit (ion) and a photon
$$\frac{1}{\sqrt{2}}(|\uparrow\rangle|H\rangle + |\downarrow\rangle|V\rangle)$$
Where the arrows represent an ion, while H/V labels represent the photon -- see equation (1) from the paper.
This can be therefore be encoded in:
$$\frac{1}{\sqrt{2}}(|00\rangle + |11\rangle)$$
I would like to know whether there is any (theoretical or experimental) reference about the creation of an arbitrary state:
$$\alpha|00\rangle + \beta|11\rangle$$
EDIT 1: I consider a valid answer anyone of the following kinds:

Nobody ever tried to implement this.

Unfortunately, what you ask for is fundamentally impossible to achieve, because...

EDIT 2:
This paper seems like a valid reference to me.

REPLY [3 votes]: The scheme you linked to already allows for significant flexibility by just tweaking the geometry.
As the paper makes clear in its (unnumbered) first displayed equation, the decay scheme they use does not naturally produce a balanced weighting of the two states in question, but rather the state
$$
|\psi⟩ = \sqrt{\frac23}|{↓}⟩|\sigma^+⟩ + \sqrt{\frac13}|{↑}⟩|\pi⟩,
$$
where $|\sigma^+⟩$ and $|\pi⟩$ represent two separate photonic modes with different polarizations and different angular emission profiles with respect to the axis defined by the external magnetic field. Figure 1(a) of the paper shows how their protocol uses the geometry of the problem to even out the superposition, by observing the emission from a direction in which the $|\sigma^+⟩$ mode is weaker than the $|\pi⟩$ mode:

Thus, by simply changing the angle between the observation axis and the applied magnetic field, the scheme should be able to produce mixtures at all sorts of possible relative weights between the even 50:50 superposition in the paper and the fully-$|\sigma^+⟩$ limit.
And, because the two states are eigenstates at different energies, you are able to set their relative phase arbitrarily by just waiting during an appropriate interval.
Of course, this only gets you half of the Bloch sphere, and not a full control, but this is only a minor change to the protocol. It is definitely possible to produce arbitrary entangled states between an atomic emitter and the polarization of the photon it has produced. I can't point you to a reference that exhibits this (and I won't undertake a search, because your question is phrased in a way that indicates that you have multiple requirements about what you're looking for that will lead you to discount alternatives, but you are not stating those requirements upfront), but I can offer you a guarantee that there is no fundamental roadblock to producing this.
What can offer roadblocks is that not everything is equally publishable. So, even if the specific task that you're interested in is possible and has been achieved, it might not have been deemed interesting enough to publish. If you want to look for it in literature, ask yourself what the task would be useful for, and then use those prospective applications as a guideline when doing your literature search.


text:
\appendix

\section{Proofs for Section~\ref{sec:preliminaries}}\label{sec:powers}
\begin{lemma}[Lemma \ref{lem:powers}]
  For any  $\power'\leq \power$,   $\Was{\power}(X,Y)\leq \Was{\power'}(X,Y)$.
\end{lemma}
\begin{proof}
For any matching between diagrams, $\Mch: X\rightarrow Y$, we have

\begin{align*}
&\left( \sum\limits_{\Dpt\in X} \left( \death(\Dpt)- \death(\Mch(\Dpt))\right)^\power + \left(\birth(\Dpt)- \birth(\Mch(\Dpt))\right)^\power \right)^{1/\power}\\
&= \left(
\sum\limits_{\Dpt\in X}
\left( \death(\Dpt)- \death(\Mch(\Dpt))\right)^\power + \left(\birth(\Dpt)- \birth(\Mch(\Dpt))\right)^\power
\right)^{\frac{\power'}{\power'\power}} \\
&\leq \left(
\sum\limits_{\Dpt\in X}
\left( \death(\Dpt)- \death(\Mch(\Dpt))\right)^{\power\power'/\power} + \left(\birth(\Dpt)- \birth(\Mch(\Dpt))\right)^{\power\power'/\power}
\right)^{1/\power'} \\

&= \left( \sum\limits_{\Dpt\in X} \left( \death(\Dpt)- \death(\Mch(\Dpt))\right)^{\power'} + \left(\birth(\Dpt)- \birth(\Mch(\Dpt))\right)^{\power'} \right)^{1/\power'}\\
\end{align*}
To prove the inequality, we show that for $0<a\leq 1$
$$(x+y)^a \leq x^a +y^a$$
This follows from the fact that the above inequality holds if and only if, $(1+t)^a \leq 1 + t^a$ with $0\leq t\leq 1$.  The result follows since the function  $1 + t^a -(1+t)^a$ is positive. Taking the derivative of this function, we obtain $a(t^{a-1} - (1+t))^{a-1}$). Since $0<a<1$, the derivative is positive and at $t=0$, the function is 0, hence the function is positive.
Since the inequality is true for any matching begin with the matching for $\power'$, the resulting matching induces a smaller norm for $\power$, and since the distance is infimum over all matchings the Wasserstein distance is smaller.
\end{proof}

\section{Proofs for Section~\ref{sec:rips}}
 The proofs in this section are straightforward or known. They are included primarily for completeness for the interested reader.

\subsection{Counter-Example}\label{sec:counterexample}
Here we describe the counter example for 1-dimensional homology in $\R^3$. Recall, we refer to the points upon which we build the Rips complex as vertices and classes in the persistence diagram as points.

We construct a configuration of points such that moving one vertex in the configuration moves a linear (in the number of vertex in the configuration) number of points in the corresponding persistence diagram.

For any constant $C$,  set $\alpha = \frac{k\pi}{C}$ for $k\leq\frac{1}{2}$. Note that this restriction is not necessary as we could have set $k<2$, but would introduce some additional corner cases. Consider a vertex $v$ at the center of a sphere of radius $r$. Place two vertices at a distance $a$ on the sphere creates a triangle with edge-lengths $(r,r,a)$ respectively.
Rotating this triangle by $\alpha$ perpendicularly to the chord between the other two vertices yields the construction.

First note that this gives rise to rectangles with points on the sphere where the long edge has a length of   $a$ and the length of the short edge is given by
$$ b = r\left(\sin \frac{\alpha}{2}\right). $$
The rectangle becomes bounded when the diagonal is introduced, which is given by
$\sqrt{a^2 +r^2 \sin^2\frac{\alpha}{2}}$. However, we would like the cone formed by the vertex at the center of the sphere to bound the cycle rather than the diagonal. Therefore,
$$r < a < \sqrt{a^2 +r^2 \sin\frac{\alpha}{2}} $$
or
\begin{align*}
r^2 &<  a^2 +r^2 \sin^2\frac{\alpha}{2} \\
&    < \frac{a^2}{1-\sin^2\frac{\alpha}{2}}
\end{align*}
Hence, if $r\in \left(a, \frac{a}{\sqrt{1-\sin^2\frac{\alpha}{2}}}\right)$, which is always possible for $\alpha>0$, the relation holds and by moving the center vertex, we change the death time of $C$ points in the diagram.

Finally, we note that the birth times for 2-dimensional homology in $\R^3$ can also be achieved by this example, by coning the rectangles off with points outside the sphere, at a distance less than $r$. In this case, the coning points outside the sphere bound the cycles and the center vertex creates a 2-dimensional class for each rectangle.

\subsection{Algebraic Mayer-Vietoris Long Exact Sequence}\label{sec:algebraic_mv}
The following is a well-known construction included here for completeness as we could not find a suitable reference.
 Assume we have the following diagrams
\begin{center}
  \begin{tikzcd}
    \cdots \arrow["\delta_{k+1}"]{r} &   A_k \arrow["\varphi_k"]{r}\arrow["f_k"]{d} &   B_k  \arrow["\psi_k"]{r}\arrow["g_k"]{d} &   C_k  \arrow["\delta_k"]{r}\arrow["h_k"]{d} &  A_{k-1} \arrow["\varphi_{k-1}"]{r} \arrow["f_{k-1}"]{d} & \cdots\\
      \cdots \arrow["\delta'_{k+1}"]{r} &   A'_k \arrow["\varphi'_k"]{r} &   B'_k \arrow["\psi'_k"]{r} &   C'_k \arrow["\delta'_k"]{r}&  A'_{k-1} \arrow["\varphi'_{k-1}"]{r}  &\cdots
  \end{tikzcd}
\end{center}
where the rows are exact and every third vertical morphism is an isomorphism. Without loss of generality, assume $h_k$ is an isomorphism for all $k$. The result says that there exists a long exact sequence
\begin{center}
  \begin{tikzcd}
    \cdots \arrow["d_{k+1}"]{r} &  A_k  \arrow["i_k"]{r}&   B_k\oplus A'_k \arrow["j_k"]{r} &B'_k
\arrow["d_{k}"]{r} & A_{k-1}   \arrow["i_{k-1}"]{r} &\cdots
  \end{tikzcd}
\end{center}
where
\begin{align*}
  i_k &= (\varphi_k,f_k)\\
  j_k &= g_k-\varphi'_k\\
\end{align*}
Before defining the connecting homomorphism $d$, we prove that $\im i_k = \ker j_k$.
First, we show that $\im i_k \subseteq \ker j_k$. Let $i_k(x)\mapsto (a,b)$, then
$$j_k\circ i_k(x) =  j_k(a,b) = g_k(a) - \varphi'_k(b) =  g_k\circ\varphi_k(x) - \varphi'_k\circ f_k(x) = 0 $$
where the last equality follows from commutativity. To show $ \ker j_k \subseteq \im i_k$, take $(a,b)\in \ker j_k$.  We observe that $b\in \ker \psi_k$. This follows from the fact that if $\psi_k(b)\neq 0$, then since $h_k$ is an isomorphism  $h_k\circ \psi_k(b)\neq 0$ and by commutativity $\psi'_k\circ g_k(b)\neq 0$. Since $b\in \ker j_k$, it follows that $g_k(b)=\varphi'_k(a)$. If $g_k(b)=0$, this is a contradiction since $h_k\circ \psi_k(b)\neq 0$. If $g_k(b)\neq 0$, then there exists an element $a$ such that $\psi_k\circ\varphi'_k(a)\neq 0 $, contradicting that the lower row is exact. Hence if $b\neq 0$, then $(a,b)\in \im i_k$. If $b=0$, since $\varphi'_k(a) = g_k(b) = 0$, which implies $a\in \ker f_k$. By exactness, this implies $b\in \im \delta'_{k+1}$. Since $h_{k+1}$ is an isomorphism, there exists an element $c\in C_{k+1}$ such that $\delta'_{k+1}\circ h_{k+1}(c) = a$. By commutativity, this implies $a\in \im f_k\circ\delta_{k+1}$ which implies $a\in \im f_k$, so $(a,0)\in \im i_k$.

Now we are ready to define the connecting homomorphism. Define
$$d = \delta_k \circ h^{-1}_k\circ\psi'_k.$$
Note that $h^{-1}_k$ is well-defined since it is an isomorphism.

Assume $b'\in \im j_k$ and let $j_k(a,b) = b'$. If $a\neq 0$, then $\psi'_k(b') = 0$ by exactness.  If $b\neq 0$, then $h^{-1}_k\circ\psi'_k (b') = \varphi_k(b)$, but $\delta_k\circ \varphi_k(b)=0$ again by exactness. Hence $\im j_k \subset \ker d_k$. The other direction is equivalent, since $b'\in \ker d_k$, either $b'\in
\ker \psi'_k$ in which case it is in $\im \varphi'_k$ and so in $\im j_k$. Alternatively, let $c = h^{-1}_k\circ\psi'_k(b') \neq 0$. If $c\in \ker \delta_k$ then there must exist an element $b$ such that $\psi_k(b)=c$ and by commutativity $g_k(b) = b'$ again implying $b' \im j_k$.

All that remains is to show that $\im d_k = \ker i_{k-1}$. First, we show $\im d_k \subseteq \ker i_{k-1}$. Let $b'\in \im d_k$ and hence $b'\not\in\im j_k$. This implies that $b'\not\in \im \varphi'_k$ and $b'\not\in \im g_k$. Since $b'\not\in \im \varphi'_k$, it follows that $b'\in \im \delta'_k$. Hence, we can choose a lift $c = h_k^{-1}(\delta'_k(b'))$.
 By commutativity, $c\not\in \im \varphi_k$, since $b'\not\in \im g_k$. Hence, $\delta_k(c)\neq 0$.  Now, applying $i_k$ to $\delta_k \circ h^{-1}_k\circ\psi'_k$. We observe that
$f_k\circ d_k = 0$, since $f_k\circ \delta_k \circ h^{-1}_k\circ\psi'_k = \delta'_k\circ\psi'_k =0$  by exactness. Similarly, $\varphi_{k-1}\circ \delta_k \circ h^{-1}_k\circ\psi'_k =0 $ since
 $\varphi_{k-1}\circ \delta_k =0 $.

 Finally, we show that $ \ker i_{k-1}  \subseteq\im d_k$. If $a\in \ker i_{k-1}$, then $a\in \ker \varphi_{k-1}$ and  $a\in \ker f_{k-1}$. Hence, there exists a $c$ such that $\delta_k(c) = a$. Since $h_k$ is an isomorphism, $h_k(c)\neq 0$. By commutativity,  $h_k(c)\in \ker \delta'_k$ so  $h_k(c) \in \im \psi'_k$. Hence there exists a $b'$ such that $\psi'_k(b')=c$ and so $d_k(b')=a$, completing the proof.
 \subsection{Basic Geometric Lemmas}
The following are basic facts of Euclidean planar geometry that we make use of.
 \begin{lemma}\label{lem:cosine}
  Consider a triangle $(xya)$ such that $\angle axy \leq \frac{\pi}{3}$ and $(ax)$ is longer than $(xy)$, then
  $$ ||a-y|| \leq ||x-y||$$
 \end{lemma}
 \begin{proof}
  Applying the cosine rule
  \begin{align*}
   ||a-y||^2 & = ||a-x||^2 + ||x-y||^2 - 2||a-x||\cdot||x-y||\cos(\angle axy)                                                      \\
             & \leq ||a-x||^2\left(1 + \left(\frac{||x-y||}{||a-x||}\right)^2 - 2\frac{||x-y||}{||a-x||}\cos(\angle abc)   \right) \\
             & \leq ||a-x||^2\left(1 + \left(\frac{||x-y||}{||a-x||}\right)^2 - \frac{||x-y||}{||a-x||}  \right)                   \\
             & \leq ||a-x||^2
  \end{align*}
  since  $\cos (\angle axy) \geq \frac{1}{2}$ and $\frac{||x-y||}{||a-x||}\leq 1$.
 \end{proof}
 \begin{figure}\label{fig:cone_picture1}
  \begin{tikzpicture}
   \begin{scope}
    \draw[fill] (0,0) circle (1pt) node[left] {$x$};
    \draw[fill] (4,0) circle (1pt) node[right] {$y$};
    \draw[fill] (5,2) circle (1pt) node[right] {$a$};
    \draw[red,<->] (0:1) arc (0:21:1) node[right,near end] {$\angle axy$};
    \draw[fill=blue,opacity=0.2] (5,2) -- (0,0) -- (4,0);

    \draw (0,0) -- (4,0);
    \draw (0,0) --  (5,2);
    \draw (4,0) --  (5,2);
   \end{scope}
  \end{tikzpicture}
  \caption{The setting for Lemma~\ref{lem:cosine}. }
 \end{figure}

 \begin{figure}\label{fig:cone_picture2}
  \begin{tikzpicture}
   \begin{scope}
    \draw[fill] (0,0) circle (1pt) node[left] {$x$};
    \draw[fill] (4,0) circle (1pt) node[below] {$y$};
    \draw[fill] (1,3) circle (1pt) node[above] {$a$};
    \draw[fill] (5,2) circle (1pt) node[right] {$y'$};

    \draw[fill] (5.39,0) circle (1pt) node[right] {$z$};

    \draw[red,<->] (0:1) arc (0:21:1) node[right,near end] {$\angle yxy'$};
    \draw[red,<->] (22:1) arc (22:70:1) node[right,near end] {$\angle axy'$};
    \draw (0,0) --   (4,0);
    \draw (0,0) --   (5,2);
    \draw[dashed] (0,0) --   (5.39,0);
    \draw[dashed] (1,3) --   (5.39,0);
    \draw[dashed,thick,blue] (0:5.39) arc (0:22:5.39);
    \draw (1,3) --  (4,0);
    \draw[thick,dashed] (5,2) --  (4,0);
    \draw[thick,dashed] (5,2) --  (1,3);

    \draw[fill=blue,opacity=0.2] (0,0) -- (1,3) -- (4,0);
    \draw[fill=blue,opacity=0.2] (5,2) -- (1,3) -- (4,0);

    \draw (0,0) --  (1,3);
   \end{scope}
  \end{tikzpicture}
  \caption{The setting for Lemma~\ref{lem:far_cone}}
 \end{figure}

 \begin{lemma}\label{lem:far_cone}
  Consider a triangle $(axy)$ such that $\angle axy \leq \frac{\pi}{2}$ and $e=(x,y)$ is the longest edge. For any point $y'$ such that $\angle yxy'\leq \frac{\pi}{3}$ and $||x-y'||\geq||x-y||$,
  $$ ||y-y'|| \leq ||y'-x|| \qquad\qquad ||y'-a|| \leq ||y'-x|| $$
 \end{lemma}
 \begin{proof}
  See Figure~\ref{fig:cone_picture2}. $||y-y'|| \leq ||y'-x||$ follows from Lemma~\ref{lem:cosine} applied to the triangle $(yxy')$.

  To show $||y'-a|| \leq ||y'-x||$, we first define the point $z$ which lies on $(xy)$  such that $||z-x|| = ||y'-x||$. It is straightforward to show that  $||z-a||\leq ||y'-a||$, as the distance from $y'$ to $a$ increases as we decrease $\angle yxy'$.
  Hence, if $||z-a||<||z-x||$, it implies the inequality. Consider the above picture
  \begin{center}
   \begin{tikzpicture}
    \begin{scope}
     \draw[fill] (0,0) circle (1pt) node[left] {$x$};
     \draw[fill] (4,0) circle (1pt) node[below] {$y$};
     \draw[fill] (1,3) circle (1pt) node[above] {$a$};
     \draw[fill] (5.39,0) circle (1pt) node[right] {$z$};
     \draw[fill] (1,0) circle (1pt) node[below] {$b$};
     \draw (1,0) --   (1,3);
     \draw (0,0) --   (4,0);
     \draw[dashed] (0,0) --   (5.39,0);
     \draw[dashed] (1,3) --   (5.39,0);
     \draw (1,3) --  (4,0);
     \draw (0,0) --  (1,3);
    \end{scope}
   \end{tikzpicture}
  \end{center}
  By assumption $||a-y||\leq||x-y||$. Using the Pythagorean theorem,
  \begin{align*}
   ||a-z||^2 & = ||a-b||^2 + ||b-z||^2 =  ||a-y||^2- ||b-y||^2 +(||b-y||+||y-z||)^2=  \\
             & = ||a-y||^2 + 2||b-y||\cdot||y-z|| + ||y-z||^2                         \\
             & \leq (||x-b||+||b-y||)^2  + 2||b-y||\cdot||y-z|| + ||y-z||^2           \\
             & =    ||x-b||^2 +||b-y||^2    +2||b-y|| (||x-b|| + ||y-z|| )+ ||y-z||^2 \\
             & \leq (||x-b|| +||b-y|| + ||y-z||)^2 = ||x-z||^2
  \end{align*}
 \end{proof}

 \section{Proof of the Rearrangement Inequality for $p\geq 1$}\label{sec:rearrangement}

\begin{lemma}\label{lem:convex_inequality}
  Given a sequence
  $$a_n \leq a_{n-1} \leq \ldots \leq a_{1} \leq b_{1}\leq b_{2} \leq \ldots \leq b_{n}$$
  and an increasing convex function $f$, the cost of a matching $\Mch: \{b_i\}\mapsto \{a_j\},$
  is
$$f(\Mch) = \sum\limits_{i=1}^n f(b_i-\Mch(b_i)).$$
The identity matching, i.e. $b_i\mapsto a_i$, maximizes this sum.
\end{lemma}
 \begin{proof}
   This follows the proof of the rearrangement inequality in \cite[(p. 78)]{steele2004cauchy}.  Begin with any matching $\Mch=\Mch_0$. Define $\Mch_1$ by performing an inversion:
  \begin{align*}
    (b_i,\Mch_{i-1}(b_i)) &\mapsto (b_i,a_i)\\
    (\Mch^{-1}_{i-1}(a_i),a_i)) &\mapsto (\Mch^{-1}_{i-1}(a_i),Mch_{i-1}(b_i))
  \end{align*}
  That is we work from the middle of the sequence outward.

  Since all pairs with index less than $i$  are paired, it follows that
  $$\Mch(b_i) \leq a_i\leq  b_i \leq \Mch^{-1}(a_i)$$

  Now we show that $f(\Mch_i)\leq f(\Mch_{i+1})$. If the matching is an identity on $i$, the statement is trivial. Let  $j$ denote the index of $\Mch(b_i)$ and $j'$ the index of $\Mch^{-1}(a_i)$. We must show
  $$f(b_{j'} - a_i) + f(b_i - a_j) \leq f(b_{j'} - a_j) + f(b_i - a_i)$$
As we are only concerned with differences, without loss of generality we can set $a_j'=0$ and rearranging terms gives
$$b_i^\power  - (b_i - a_i)^\power\leq b_{j'}^\power - (b_{j'} - a_i)^\power$$
Using the substitutions
$$x_1 = b_i - a_i, \qquad x_2=b_{j'} - a_i \qquad  d = a_i$$
we obtain
$$f(x_1+d)  - f(x_1) \leq f(x_2+d) - x_2$$
For convex functions, the quantity
$$\frac{f(y) - f(z)}{y-z}$$
is non-decreasing in both $y$ and $z$. Since the function is increasing, the quantity is always positive for $y\leq z$. This implies
$$ \frac{f(x_2+d) - f(x_2)}{d} \leq \frac{f(x_1+d) - f(x_2)}{d+x_2-x_1} \leq \frac{f(x_1+d) - f(x_1)}{d} $$
 \end{proof}
 \begin{corollary}\label{cor:general_p}
   For any sequence as above and $\power\geq 1$, the identity matching maximizes
   $$\sum\limits_{i=1}^n (b_i-\Mch(b_i))^\power$$
 \end{corollary}
 \begin{proof}
   As the exponentiation function $x^\power$ is an increasing convex function for non-negative $x$,  Lemma~\ref{lem:convex_inequality} implies the result.
 \end{proof}

 \section{Notation}\label{sec:notation}
 \begin{itemize}
  \item $\SC$ -- a finite CW complex
  \item $\Dgm_k(\SC,f)$ or $\Dgm_k(f)$ -- the $k$-th dimensional persistence diagram for the filtration induced on $\SC$ by a monotone function $f$.
  \item $\Dgm(\SC,f)$ or $\Dgm(f)$ -- $\bigoplus\limits_k \Dgm_k(\SC,f)$
  \item $\Mod{F}$ - a persistence module
  \item $\Dgm(\Mod{F})$ -- the persistence diagram of a persistence module
  \item $\overline{\Dgm}(\Mod{F})$ -- the finite intervals of a persistence diagram
  \item $\Cech(\Xp)$ -- the \v Cech complex
  \item $\VR(\Xp)$ -- the Vietoris-Rips complex, where $\VR_\delta$ refers to the complex at parameter $\delta$
  \item $\Sh$ -- the shadow complex
  \item $\Was{\power}{\cdot}{\cdot}$ - the $\power$-th Wasserstein distance
  \item $\Wasalg{\power}{\cdot}{\cdot}$ - the $\power$-th algebraic Wasserstein distance
 \item $\Cg_k(\cdot)$ -- $k$-dimensional chain complex
  \item $\Hg_k(\cdot)$ -- $k$-dimensional homology
  \item $\Cov$ -- an open cover of a space
  \item $\Nrv$ --  the nerve of an open cover
  \item $\Omega$ -- the inverse map taking points from a persistence diagram to an underlying CW-complex
  \item $\Mch$ -- a bijective matching between sets
  \item $\Csp$ -- a corrspondence
  \item $\St()$ -- the open star of a simplex
  \item $\Cl()$ -- closure of a set
  \item $\Lnk()$ -- the link a simplex, i.e. $\Cl(\St()) - \St()$
  \item $\mono$ -- monomorphism
  \item $\epi$ -- epimorphism
  \item $\rk$ -- the rank
    \item $\Gpt, \Gpta$ -- geometric points
    \item $\Dpt$ -- points in persistence diagrams
    \item $\birth(\Dpt)$ -- birth time of point $\Dpt$
    \item $\death(\Dpt)$ -- death time of point $\Dpt$
    \item $\ell$ -- lifetime of point $\Dpt$, i.e. $\death(\Dpt)-\birth(\Dpt)$
  \end{itemize}

config: arxiv


TITLE: Does this algorithm terminate in all scenarios?
QUESTION [0 upvotes]: Let $x \in \mathbb{R}^p$ denote a $p$-dimensional data point (a vector). I have two sets $A = \{x_1, \dots, x_n\}$ and $B = \{x_{n+1}, \dots, x_{n+m}\}$, so $|A| = n$, and $|B| = m$. Given $k \in \mathbb{N^*}$, let $d_x^{(A, k)}$ denote the mean Euclidean distance from $x$ to its $k$ nearest points in $A$; and $d_x^{(B, k)}$ denote the mean Euclidean distance from $x$ to its $k$ nearest points in $B$.
I have the following algorithm:

$A' = \{ x_i \in A \mid d_{x_i}^{A, k)} > d_{x_i}^{(B, k)} \}$ ... (1)
$A = A \setminus A'$ ... (2)
$B' = \{ x_i \in B \mid d_{x_i}^{A, k)} < d_{x_i}^{(B, k)}$ ... (3)
$B = B \setminus B'$ ... (4)
$A = A \cup B'$ ... (5)
$B = B \cup A'$  ... (6)
Repeat (1), (2), (3), (4), (5) and (6) until: (no element moves from $A$ to $B$ or from $B$ to $A$, that is A' and B' become empty) or (|A| $\leq$ k or |B| $\leq$ k)

Does this algorithm terminate, and if so, is it possible to easily prove it ? Is it also possible to have an upper bound for the number of iterations required to terminate ?
Note: "The $k$ nearest points to $x$ in a set $S$" means: The $k$ points (other than $x$) in $S$, having the smallest Euclidean distance to $x$.

REPLY [2 votes]: Counterexample
This algorithm does not terminate. Here is a counterexample, in dimension $p=2$, with $m=n=2$ points in each set, only looking at the $k=1$ nearest neighbours.
$$
A=\left\{
\begin{pmatrix}-2\\0\end{pmatrix},
\begin{pmatrix}0\\1\end{pmatrix}
\right\}
\qquad
B=\left\{
\begin{pmatrix}3\\0\end{pmatrix},
\begin{pmatrix}0\\-1\end{pmatrix}
\right\}
$$
The two points with $0$ in their first coordinate will change sets in every iteration. After two iterations you are back where you started, so this cannot terminate.
Step by step
Since you doubt my counterexample in a comment below, here is a detailed step-by-step execution of your algo in this case. I'm using the same names for the points as you did:
$$
x_1=\begin{pmatrix}-2\\0\end{pmatrix} \qquad
x_2=\begin{pmatrix}0\\1\end{pmatrix} \qquad
x_3=\begin{pmatrix}3\\0\end{pmatrix} \qquad
x_4=\begin{pmatrix}0\\-1\end{pmatrix} \\
A=\{x_1,x_2\} \qquad B=\{x_3,x_4\}
$$

$A'=\{x_2\}$ since $d_{x_2}^{(A,1)}=\lVert x_1-x_2\rVert=\sqrt5>2=\lVert x_4-x_2\rVert=d_{x_2}^{(B,1)}$
$A=A\setminus A'=\{x_1\}$
$B'=\{x_4\}$ since $d_{x_4}^{(A,1)}=\lVert x_1-x_4\rVert=\sqrt5<\sqrt{10}=\lVert x_3-x_4\rVert=d_{x_4}^{(B,1)}$
$B=B\setminus B'=\{x_3\}$
$A=A\cup B'=\{x_1,x_4\}$
$B=B\cup A'=\{x_2,x_3\}$
$A'=\{x_4\}$ since $d_{x_4}^{(A,1)}=\lVert x_1-x_4\rVert=\sqrt5>2=\lVert x_2-x_4\rVert=d_{x_4}^{(B,1)}$
$A=A\setminus A'=\{x_1\}$
$B'=\{x_2\}$ since $d_{x_2}^{(A,1)}=\lVert x_1-x_2\rVert=\sqrt5<\sqrt{10}=\lVert x_3-x_2\rVert=d_{x_2}^{(B,1)}$
$B=B\setminus B'=\{x_3\}$
$A=A\cup B'=\{x_1,x_2\}$
$B=B\cup A'=\{x_3,x_4\}$


TITLE: Can we get different answers using FTC1 and FTC2?
QUESTION [0 upvotes]: I have just starting studying Calculus II after a very long break in college. In my class we just covered FTC1 and FTC2 in class and I was talked to evaluate two integrals using both of the ways.
I have gotten different answers using FTC1 and FTC2 and I was wondering if it is correct!

REPLY [1 votes]: Your FTCII calculation is correct. You have used linearity of the integral, splitting it into two smaller integrals, which is a perfectly valid approach. I might approach the problem more directly, and simply write
$$
F(x) = \int_1^x t - t^3 \ dt = \left[\frac{1}{2}t^2 - \frac{1}{4}t^4 \right]_1^x = \frac{x^2}{2}-\frac{x^4}{4}-\frac{1}{2} + \frac{1}{4} = \frac{- x^4+2x^2 - 1}{4} = -\frac{(x^2-1)^2}{4},
$$
matching your answer.
But as spaceisdarkgreen pointed out in the comments, if we then wish to compute the derivative, the chain rule gives
$$
F'(x) = -\frac{1}{4}(2x)(2)(x^2-1) = -x(x^2-1) = -x^3 + x,
$$
matching your FTCI answer (before the false manipulation in your last equality, I think you just forgot to include the $-$ sign). All good!


TITLE: How large should n be to guarantee the error in approximating ln(.5) with ln(1+x) centered at a = 0 is less than .0001?
QUESTION [2 upvotes]: I am asked to find the number of terms such that the error of approximating ln(.5) < 0.001 using Taylor's inequality. I can do these using the remainder theorem from the alternating series test, but we are required to use Taylor's inequality.
I am struggling right now to find M such that
$ |f^{n+1}(x)| \le M$.
I have found the derivative
$|f^{n+1}(x)|= \frac{n!}{(1+x)^{n+1}}$. My problem is with finding M. I know that $x \in (-.5, 0)$, and the derivative (absolute value) is decreasing. If I plug -.5 into the derivative to get M, M seems to just blow up and isn't bounded.
Is there any other way to find M in this situation?

REPLY [0 votes]: Making it general, we have (for $|x|<1$)
$$\log(1+x)=\sum_{n=1}^\infty (-1)^{n+1}\,\frac {x^n} n=\sum_{n=1}^p (-1)^{n+1}\,\frac {x^n} n+\sum_{n=p+1}^\infty (-1)^{n+1}\,\frac {x^n} n$$ and you want to know $p$ such that
$$R_p=\frac {x^{p+1}}{p+1} \leq 10^{-k}$$ that is to say
$$(p+1) \geq x^{p+1} 10^k$$
In the real domain, the only analytical solution of the equation is
$$p=-\frac{W\left[-10^k \log (x)\right]}{\log (x)}-1$$ where $W(.)$ is Lambert function whicch is widely used (on the search bar here, just type Lambert and you will find $3392$ entries). For sure, after computing $p$, you will use $\lceil p\rceil$.
For $k=5$ and $x=\frac 12$, this will give $p=11.9183$, then $\lceil p\rceil=12$.
Let us check
$$R_{11}=\frac{1}{49152}\sim 2.03\times 10^{-5} > 10^{-5}$$
$$R_{12}=\frac{1}{106496}\sim 9.39\times 10^{-6} < 10^{-5}$$
If you do not have access to this function, you can use the expansion for large values of $t$ as given in the linked page
$$W(t)= L_1-L_2+\frac{L_2}{L_1}+\frac{L_2(L_2-2)}{2L_1^2}+\frac{L_2(2L_2^2-9L_2+6)}{6L_1^3}+\cdots$$ where $L_1=\log(t)$ and $L_2=\log(L_1)$.
From a coding point of view, there is big advantage to know in advance the number of terms to be added since, in the calculation loop, you will not have any more IF test (this is a very expensive operation).
Suppose that, for the same $x=\frac 12$, we want $k=100$. We should find $p=322.854$  then $\lceil p\rceil=323$. Checking again
$$R_{322}=4.81\times 10^{-100} >10^{-100}$$ $$ R_{323}=9.03\times 10^{-101} <10^{-100} $$
Welcome to the world of the magnificent Lambert function !


TITLE: What is the experimental reason to believe in a Unified Theory?
QUESTION [4 upvotes]: Is there any experimental findings on behalf of the belief that all forces can be unified to a single force? Or the idea has it's inspiration from kind of thought experiment or philosophical belief.

REPLY [4 votes]: The experimental evidence comes from the running coupling constants.

One important clue that they (couplings)  might all be the same comes from the fact that the coupling constants are not constants at all. Instead, they vary with the energy of the phenomenon in which they are measured. The value of as quoted above is only true for phenomena that occur at around 91 GeV of energy, the mass of the Z particle. Similarly, the coupling constant for the electromagnetic force has been measured to be about 1/137 in studies of atomic phenomena. When it is measured at the higher energies of LEP, however, its value has increased to about 1/128! The value of the coupling constant is said to 'run' with energy, and the constants themselves are usually referred to as running coupling constants.

You might read the paragraph comparison to experiment on page 29 in this article.
In this talk the experimentally found running of the strong coupling is seen on a graph on  page 3.
The running of the electromagnetic coupling is studied experimentally here.


TITLE: Is $\{A_{n} \text{ finitely often} \}$ the complement of $\{A_{n} \text{ infinitely often} \}$
QUESTION [1 upvotes]: I tried to think about this using De-Morgan's laws:
$\{A_{n} \text{ i.o} \} = \bigcap_{m\ge 1} \bigcup_{n\ge m} A_n$
so $\{A_{n} \text{ i.o} \}^c = \bigcup_{m\ge 1} \bigcap_{n\ge m} {A_n}^c$
once you take the union of $m=1$ and $m=2$:
$({A_1}^c\cap{A_2}^c\cap{A_3}^c...)\cup({A_2}^c\cap{A_3}^c...)$
this means the element does not have to be in  ${A_1}^c$ so can be in $A_1$. This continues for all $A_n$ so I'm not sure this shows the element has to be in finitely many $A_n$?
Have I thought about this the correct way?

REPLY [1 votes]: Your reasoning is correct. It helps to think of $\cap$ as "and" and $\cup$ as "or". Your last expression is a union, so for an element to be in it, it has to be in one of the sets $(A_m^c\cap A_{m+1}^c\cap\dots)$, meaning that it is not in any of the sets $A_n$ for $n\ge m$. This means it's only in finitely many of the $A_i$.
Another good way to think about it is that you're applying negation rules to logical quantifiers: $\lnot(\forall m\ge1,\exists n\ge m,P_n)\iff \exists m\ge1,\forall n\ge m,\lnot P_n$.


TITLE: "Either A and B is open, then A + B is open" (typo sense-making, Stein Shakarchi Real Analysis)
QUESTION [0 upvotes]: Please advise about the most reasonable way to read this statement. My interpretations are below.
The authors do not define the set operation A + B; I assume A + B = $A \cup B$.
Their statement "Show that if either A and B is open, then A + B is open" is incorrect unless it should read "both A and B are open."
Counterexample when only A is open and B is not: Let $A=B_1(O)\subset \mathbb{R}^2,\ B=\overline{B_1(2)}\subset\mathbb{R}^2$. Then there exist points in B which are not interior points of $A\cup B$ (e.g. (3,0), a boundary point of B).

REPLY [3 votes]: Let $A,B\subseteq G$ be subsets of a topological group (for example $G=\mathbb R^n$).
We have
$$A+B=\bigcup_{b\in B}A+b.$$
If $A$ is open then so is $A+x$ for all $x\in G$. So $A+B$ is a union of open sets, hence open.
In case we know that $A+B$ is open we don't know if one of $A$, $B$ is open. That is already the case in $\mathbb R$. Take for example $A=[0,1]$ and $B=\mathbb Z$. Then $A+B=\mathbb R$ is open but none of $A$ and $B$.


TITLE: Computing order of point at infinity when computing principal divisor
QUESTION [1 upvotes]: Let the elliptic curve $y^2=x^3+1 \pmod{13}$ and the rational function $g = \frac{x^2}{y}$. Compute the principal divisor of $g$ on the above curve.
What I've done. First homogenize both (the curve and the function) and obtain
$\frac{x^2}{yz} = 0$ and $y^2z-x^3-z^3=0$.
Zeros are obtain solving $x^2 = 0$ and $y^2z-x^3-z^3=0$. This leads to the points $\mathcal{O}(0:1:0)$ with order $n_{\mathcal{O}}$ and points $P(0:1:1)$ and $Q(0:-1:1)$ both of order $2$ (i.e note that $x$ is a uniformizer and $g = x^2 \cdot \frac{1}{yz}$ and $\frac{1}{yz}$ neither has a pole nor a zero in points $P$ and $Q$.)
Poles are obtain solving $yz = 0$ and $y^2z-x^3-z^3=0$. This leads to the points points $\mathcal{O}(0:1:0)$ with order $m_{\mathcal{O}}$ and points $R(-1:0:1)$, $S(4:0:1)$ and $T(-3:0:1)$ all three of them with order $-1$ (i.e. note that $y$ is a uniformizer and $g=y^{-1} \cdot \frac{x^2}{z}$ and $\frac{x^2}{z}$ neither has a pole nor a zero in $R,S$ and $T$).
Hence by definition
$$div(g) = n_{\mathcal{O}}(\mathcal{O}) + 2(P) + 2(Q) - m_{\mathcal{O}}(\mathcal{O}) - (R) - (S) - (T)$$
It is well know that $deg(div(g)) = 0$ in this case, so I expect that $n_{\mathcal{O}} = m_{\mathcal{O}} - 1$.
Question: How to compute $n_{\mathcal{O}}$ and $m_{\mathcal{O}}$ using uniformizers ?

REPLY [4 votes]: Let $k=\mathbb F_{13}$. The trick with these kinds of problems is always to work in affine charts. Since we already have the equation in the $xy$-chart, notice that the only way the function $\frac{x^2}{y}$ (or really it's image in the coordinate ring $A=k[x,y]/\langle y^2-x^3-1 \rangle$) vanishes on the $xy$ chart on your curve is if $x=0$ so that $y=\pm 1$ (which are distinct mod $13$.)
Let's start and zoom in at the point $p_1=(0,1)$, i.e. localize the coordinate ring $A$ at the ideal $m_1=(x-0,y-1)=(x,y-1)$. The maximal ideal in the local ring $A_{m_1}$ is $(x,y-1)$ but we twist $y^2-1-x^3=0$ into $y-1=\frac{x^3}{y+1}$ and since $y+1$ does not vanish at our point, we are allowed to invert it in $A_{m_1}$. So the ideal $(x,y-1)$ in $A_{m_1}$ is actually secretly just $(x) $ so $x$ is a uniformizer which is just fancy speak for $\text{ord}_{p_1}(x)=1 $.
So $$\text{ord}_{p_1}(\frac{x^2}{y})=\text{ord}_{p_1}(x^2)-\text{ord}_{p_1}(y)=2\text{ord}_{p_1}(x)-0=2$$
Next up, the same song and dance with $p_2=(0,-1)$ with the main point being $y+1=\frac{x^3}{y+1}$ shows that
$$\text{ord}_{p_2}(\frac{x^2}{y})=2$$
Similarly, for poles we focus on when $y=0$ so $x=-1,-3,4$ so play the same uniformizer game we did above.
Now the only point on the curve missing from our chart must happen when $z=0$ so plug in $z=0$ into your homogeneous equation to get $x^3=0 $ so that $x=0,y=1,z=0$. So lets work in the $x,z$ chart to get $z-x^3-z^3$ and our function $\frac{x^2}{y}$ becomes $\frac{x^2}{z}$ and the only point we haven't dealt with is $(0,0)$ (This is important to keep track of, otherwise you will end up double counting points, which is bad). So localize $k[x,z]/(z-x^3-z^3)$ at $(x,z)$. But since the equation here can be twisted into $$z=\frac{x^3}{(1-z^2)} ,$$
 we have $x$ is the king, I mean uniformizer. So $$\text{ord}_{(0,0)}(\frac{x^2}{z})=2-1=1$$ so zero of order $1$.
This is really wordy but the point is once you learn this, it's pretty routine.
Now I should stop procrastinating and do my own work.


TITLE: Suppose $K$ is compact. What other types of coverings must have finite subcovers?
QUESTION [6 upvotes]: Let $X$ be a topological space. Call $S\subseteq X$ an $\mathcal{O}$-set if there exists an open set $O$ such that $O\subseteq S \subseteq \overline{O}$. Suppose $X$ is compact. Is it true that any covering of $X$ by $\mathcal{O}$-sets has a finite subcovering?
It seems intuitively that it should be true for $X \subseteq \mathbb{R}$, but I'm not sure if this relies on some other topological property of the reals.
Below the fold I have included some background as to why I would ask such a question, if you're interested.

Background: here on Reddit user posted an interpretation of compactness based on Zeno's paradox--the "steps" that the racer takes are considered as open sets, and a compact "racetrack" implies that the the race must be completed in finite time if the velocity function is continuous. I wondered: why are open sets the correct way to view the steps? Then I got to thinking, a set being non-empty and open is one way to guarantee that it really represents a "chunk" of the space...not just a point, or a subspace of lower dimension (if it were a vector space). But still, why open? I'd rather have the freedom to consider the "step" as any set containing an open set, but contained in its closure. If this is possible, then it becomes clear that this holds for any reasonable interpretation of a "step".

REPLY [3 votes]: ADDED: There's actually a counterexample way simpler than the one below: $$[0,1]=\left[\frac12,1\right]\cup\bigcup_{n\in\mathbb N}\left[0,\frac12-\frac{1}{4n}\right).$$

ORIGINAL ANSWER: Let $X=$ the Alexandroff one-point compactification of $\mathbb R$, with the extra point denoted as $\omega$. Let $S_0=\{\omega\}\cup(-\infty,0]\cup[1,\infty)$ and $S_n=(0,(n-1)/n)$ for all $n\in\mathbb N$.
Claim$\phantom{---}$ $S_n$ is an $\mathcal O$-set for all $n\in\mathbb N\cup\{0\}$.
Proof:$\phantom{---}$ First note that $S_0$ is closed, because $S_0^c=(0,1)$ is open. Moreover, $U_0=\{\omega\}\cup(-\infty,0)\cup(1,\infty)$ is open because its complement $[0,1]$ is compact in $\mathbb R$. Additionally, neither $\{\omega\}\cup(-\infty,0]\cup(1,\infty)$ nor $\{\omega\}\cup(-\infty,0)\cup[1,\infty)$ is closed (their complements are $(0,1]$ and $[0,1)$, respectively, which are not open), so $\overline{U_0}$ must be actually equal to $S_0$.  It follows that $S_0$ is an $\mathcal O $-set. If $n\in\mathbb N$, then $S_n$ is open in $\mathbb R$ and hence in $X$, so it is obviously an $\mathcal O$-set. $\blacksquare$
A counterexample is now straightforward. $X$ is compact by construction (and Hausdorff, so it is normal as well—it is even metrizable!), but the conjecture fails: $$X=\bigcup_{n\in\mathbb N\cup\{0\}}S_n,$$
and there is no finite subcover.


TITLE: Overall equal probabilities when replacing choices
QUESTION [3 upvotes]: Explanation:
I want to choose marbles from a bag x times. Once chosen, the marbles are not put back in the bag. After each choice, some number of new marbles are added to the bag. x is less than the total number of marbles, i.e. not all marbles will be chosen. My goal is to make it so every marble has an equal chance of being chosen overall (or as close as possible).
Example 1 (x = 2):
Start with marbles A and B. Choose one at random. Let's say B was chosen. Now marble C is added to the bag and the second marble is chosen. A has already had a 50% chance of being chosen, so giving A and C the same odds would mean A was more likely than C. Is there any weighting possible to make it such that A, B, and C were all equally likely in the beginning? I assume it is not possible because C was not available in the beginning. I believe the next best thing would be to make it so A and C were equally likely from the beginning. My guess is that this can be done by making C twice as likely as A for the second draw (A should have 1/3 chance and C should have 2/3 chance). However, I am not sure how to prove it.
Example 2 (x = 3):
Start with marbles A and B. Again let's say B is chosen. Now Marbles C and D are added and a second is chosen. Again, I want to say C and D should be twice as likely as A. Assuming I was correct before, A should have a 1/5 chance and C and D should both have a 2/5 chance. Let's say C was chosen. Now marbles E and F are added and so the odds of each marble for the third draw should be 0.0909 for A, 0.1818 for C, 0.3636 for E and 0.3636 for F.
Is this correct? How can it be shown mathematically? Does my general approach of marbles added after the last draw being twice as likely as the marbles in the draw before that which are twice as likely as the marbles in the draw before that... work?

REPLY [1 votes]: There are actually three parameters here: the number of starting balls $s$, the number added each time $t$, and the total number of draws to be undertaken $x$
You can say that the total number of balls involved is $xt + s-t$ and if $x$ are drawn then the average overall probability a particular ball is drawn at some stage is $\frac{x}{xt+s-t}$.  For the balls added after the penultimate draw, this must be the probability they are individually drawn in the final round and so the combined probability they are drawn that round is $\frac{xt}{xt+s-t}$, leaving a combined probability of $\frac{s-t}{xt+s-t}$ for the all the other balls in that round.
This puts a constraint on the parameters: you need $s\ge t$ (or $x=1$) to avoid an impossible situation.
I think you can then work backwards, finding weights for the balls introduced after the antepenultimate draw which keep to same probabilities, and so on.  In your examples:

$s=2,t=1,x=2$ imply $\frac{x}{xt+s-t} =\frac{2}{3}$ so that means you want $C$ to have a probability of $\frac{2}{3}$ in the second round and so the older remaining ball to have probability $\frac13$.  This then gives each of the original balls an overall probability of being chosen $\frac12+\frac12 \times \frac13=\frac23$, so this works
$s=2,t=2,x=3$ imply $\frac{x}{xt+s-t} =\frac{1}{2}$ so that means you want $E$ and $F$ to each have a probability of $\frac{1}{2}$ in the third draw and so the older remaining balls to have probability $0$. That shuts out $C$ and $D$ from the third draw so you want them to each have a probability of $\frac{1}{2}$ in the second draw, shutting $A$ and $B$ out from the second and third rounds.  This works, just, since each ball has an overall probability of $\frac12$ of being drawn, but does not allow a finite weighting ratio between new and remaining balls


TITLE: calculation of divergence of a covariant vector
QUESTION [1 upvotes]: $(x,y)$ are rectangular Cartesian coordinates in a plane and
$(\bar{x},\bar{y})$ are coordinates defined by
$$x=\frac{\bar{x}^2+\bar{y}^2}{2} \ \ , \ \ y=\bar{x}\bar{y}$$ If
$A_1,A_2$ are coordinates of a covariant vector in the
$\bar{x},\bar{y}$ coordinates defined by
$A_1=A_2=(\bar{x}^2-\bar{y}^2)^2$ , show that its divergence is
$$\frac{2(\bar{x}-\bar{y})^2}{\bar{x}+\bar{y}}$$

Clearly the metric tensor in the new coordinate is $$g_{11}=g_{22}=\bar{x}^2+\bar{y}^2$$ Hence the conjugate metric tensors are $$g^{11}=g^{22}=\frac{1}{\bar{x}^2+\bar{y}^2}$$ Hence $$\text{div}(A_i)=g^{ij}A_{i \ , \ j}=g^{11}A_{1,1}+g^{22}A_{2,2}$$ "$,$" is the covariant differentiation operator . So $$\text{div}(A_i)=g^{11}\bigg(\frac{\partial A_1}{\partial\bar{x}}-\Gamma^1_{11}A_1-\Gamma^2_{11}A_2\bigg)+g^{22}\bigg(\frac{\partial A_2}{\partial\bar{y}}-\Gamma^1_{22}A_1-\Gamma^2_{22}A_2\bigg)$$ But upon calculation of the Christoffel symbols and substituting in the expression gives $$\text{div}(A_i)=\frac{4(\bar{x}^2-\bar{y}^2)(\bar{x}-\bar{y})}{\bar{x}^2+\bar{y}^2}$$ I don't get what is going possibly wrong there , any help will be appreciated .

REPLY [1 votes]: You have already correctly calculated $g_{11}$ and $g_{22}$ of metric tensor $g_{ik}$. The non-diagonal component $g_{12}=g_{21}=2xy$ (to save time I don't write barr above x and y). $G=det|g_{ik}|=(x^2-y^2)^2$, $A^1=A^2=x^2+y^2-2xy$ $(A^1=g^{11}A_1+g^{12}A_2)$. Then you can calculate Christoffel symbols (from learning point this is a good approach), but let me go directly to the goal: in case of divergence of a vector we may use the formula for $A_{i;j}g^{ij}$ (; means covariant differentiation): $A_{i;j}g^{ij}=\frac{1}{\sqrt{G}}\frac{d}{dx^i}(\sqrt{G}A^i)$(summation over repeated indexes i)$=\frac{1}{\sqrt{G}}\frac{d}{dx}(\sqrt{G}A^1)+\frac{1}{\sqrt{G}}\frac{d}{dy}(\sqrt{G}A^2)=\frac{1}{x^2-y^2}\frac{d}{dx}\left((x^2-y^2)(x-y)^2\right)+\frac{1}{x^2-y^2}\frac{d}{dy}\left((x^2-y^2)(x-y)^2\right)=2\frac{(x-y)^3}{x+y}$
Hopefully this will be helpful.


TITLE: Alternative Proofs: Closed proper subspace of Hilbert space has non-empty orthogonal complement
QUESTION [1 upvotes]: Let $(H,\langle \cdot, \cdot \rangle)$ be a Hilbert space. I am asking for people to present any short/basic/fundamental proofs they know of the following theorem.

Theorem 1: Let $V \subseteq H$ be a closed, proper vector subspace of $H$. Then $V^\perp$ contains a non-zero vector.

Theorem 1 follows pretty easily as a corollary of either of the following theorems.

Theorem 2: A vector subspace $U$ of $H$ is dense if and only if $U^\perp$ is trivial.

Notes on Theorem 2: The proof may be found here, requiring many pages of prior results. These results include $U^\perp = (\overline{U})^\perp$ and $(U^\perp)^\perp = \overline{U}$ for vector subspaces $U$ of $H$, which in turn require the "Hilbert projection theorem" and the continuity of the inner product.

Theorem 3: $H = U \oplus U^\perp$ for any closed vector subspace $U$ of $H$.

Notes on Theorem 3: A proof of this may be found here. I believe this proof requires less machinery than the proof of Theorem 2. In particular, only the construction of the orthogonal projection and some basic properties of the projection are required.

Allow me to explain my thoughts a bit further. Theorem 2 and 3 are definitely not long proofs which make use of important, fundamental facts about orthogonality in Hilbert space. However, Theorem 2 and 3 seem to be much stronger results than Theorem 1. I have a gut feeling that Theorem 1 may be proven without having to resort to the "machinery" required for Theorem 2 or 3 (if you could really call it machinery). I guess another way to put it is that I believe there should be much more basic proofs of Theorem 1.
However, I have not been able to make progress on such a proof myself. If I were to guess, a proof by contradiction making clever use of the continuity of the inner product might work. But, again, I've had no such success. Please feel free to drop your proof approaches below or any reasoning as to why you think Theorem 2 or 3 are the shortest / most fundamental methods.
I hope this post allows for a better understanding for the true difficulty or simplicity of Theorem 1. So, although this is not strictly what I am looking for, also feel free to drop any super sleek proofs of Theorem 1 that use really heavy machinery ahaha Maybe such proofs will enlighten why Theorem 1 cannot be proven in a more fundamental manner. Thanks!

REPLY [1 votes]: You can prove it without heavy machinery, but a bit of computation, using distance-minimising sequences.
Suppose $V \le H$ is closed and proper. Choose any point $x_0 \in H \setminus V$. Let
$$r = \inf_{x \in V} \|x - x_0\| > 0,$$
since $V$ is closed. Suppose $(x_n)_{n \ge 1} \in V$ such that $\lim_{n \to \infty} \|x_n - x_0\| = r$. I claim that $(x_n)$ is Cauchy.
Fix $\varepsilon > 0$. Let
$$\delta = \sqrt{r^2 + \frac{\varepsilon}{4}} - r > 0.$$
Since $\|x_n - x_0\| \to r$, there exists some $N$ such that
$$n \ge N \implies r \le \|x_n - x_0\| < r + \delta = \sqrt{r^2 + \frac{\varepsilon}{4}}.$$
Suppose $m, n \ge N$. Parallelogram law states that:
$$\|(x_m - x_0) - (x_n - x_0)\|^2 + \|(x_m - x_0) + (x_n - x_0)\|^2 = 2\|x_m - x_0\|^2 + 2\|x_n - x_0\|^2,$$
or equivalently,
$$\|x_m - x_n\|^2 = 2\|x_m - x_0\|^2 + 2\|x_n - x_0\|^2 - 4\left\|\frac{x_m + x_n}{2} - x_0\right\|^2.$$
Note that $\frac{x_m + x_n}{2} \in V$, so
$$\left\|\frac{x_m + x_n}{2} - x_0\right\|^2 \ge r^2,$$
by definition of $r$. Therefore,
$$\|x_m - x_n\|^2 < 2\left(r^2 + \frac{\varepsilon}{4}\right) + 2\left(r^2 + \frac{\varepsilon}{4}\right) - 4r^2 = \varepsilon.$$
Thus, $(x_n)$ is Cauchy. It has a limit $x_\infty$, by completeness. I claim that $x_\infty - x_0 \in V^\perp \setminus \{0\}$.
First,
\begin{align*}
\|x_\infty - x_0\| &= \left\|\lim_{n \to \infty} x_n - x_0\right\| \\
&= \lim_{n \to \infty} \|x_n - x_0\| \\
&= r > 0,
\end{align*}
so $x_\infty - x_0 \neq 0$.
To see that $x_\infty - x_0 \in V^\perp$, fix $v \in V \setminus \{0\}$, and suppose $t \in \Bbb{C}$ (or just $\Bbb{R}$, if we're in the reals) to be decided shortly. Then $x_\infty + tv \in V$, so
$$\|x_\infty + tv - x_0\| \ge r = \|x_\infty - x_0\|.$$
Then
\begin{align*}
0 &\le \|x_\infty - x_0 + tv\|^2 - \|x_\infty - x_0\|^2 \\
&= |t|^2\|v\|^2 - 2\operatorname{Re}(t \langle x_\infty - x_0, v \rangle).
\end{align*}
In particular, let
$$t = \frac{\langle v, x_\infty - x_0\rangle}{\|v\|^2},$$
so that the above becomes:
$$0 \le -\frac{|\langle x_\infty - x_0, v\rangle|^2}{\|v\|^2},$$
which implies $\langle v, x_\infty - x_0\rangle = 0$ for all $v \in V \setminus \{0\}$ (and hence all $v \in V$). Thus $x_\infty - x_0 \in V^\perp$.


TITLE: Constrained optimization: equality constraint
QUESTION [3 upvotes]: I have this very general problem (for $n>2$):
$$
\begin{align}
& \max Z = f(x_1,\ldots ,x_n) \\[10pt]
\text{s.t. } & \sum_{i=1}^{n} x_i = B \\[10pt]
& x_i \geq 0
\end{align}
$$
Assume that the function $f(x_1,\ldots,x_n)$ is not a "nice" function, non-linear, non-convex, hardly differentiable. Therefore, in general I cannot use the Lagrangian multipliers method.
For finding the $x_i$'s, I am thinking to numerical optimization approaches.
My questions are:
1) Is there a common, smart and effective way to deal with such an equality constraint in numerical optimization?
2) Am I missing something? Do you have any other direction to suggest?
Thanks

REPLY [2 votes]: if $B > 0$ by dividing on B you can bring all your constraint, both equality and inequalities, to single "unit simplex constraint" There are a lot of papers for optimization on unit simplex. Some of them are using  solving unconstrained problem and projecting it projection on unit simplex each iteration (projected gradient for example, google for it if interested) For projection on simplex algo check wiki:
http://en.wikipedia.org/wiki/Simplex#Projection_onto_the_standard_simplex
If cost function is non-convex you will have to use some form of trust region for those method to work.
Also most(or even all) of the method for optimization on $L_1$ ball are using methods which are also working for unit simplex constraint. You can google on that subject too.
Another (though it could be combined with first) approach is alternating direction method - separate you cost function into smooth and nonsmooth part
$F_0(x) = f(x) + g(x)$
and split them with new variables, after taht descend by x and u intermittengly
$F(x, u) = f(x) + g(u)  + (x-u)^2$
and you can add simplex constraint here as well
$F(x, u, v) = f(x) + g(u) + \delta_\Delta^-(v) + (x-u)^2 + (x-v)^2$
It usually used in the form of augmented Lagrangian or ADMM -Alternating Direction Method of Multipliers
All this assume g is convex and and f can be approximated by convex at least locally


TITLE: m-point-homogeneous, but not (m+1)-point-homogeneous
QUESTION [6 upvotes]: It is straightforward to check that the discrete cube $Q=\{0,1\}^n$ with $\ell^1$-metric is 3-point-homogeneous, but not 4-point-homogeneous (assuming $n$ is large).
In other words, if $A\subset Q$ has at most 3 points, then any distance-preserving map $A\to Q$ can be extended to an isometry $Q\to Q$, but there is a 4-point set for which it does not hold.
I would like to see an example of $m$-point-homogeneous that is not $(m+1)$-point-homogeneous for $m\geqslant 4$.

REPLY [3 votes]: The Schläfli graph, viewed as a metric space, is 4-homogeneous but not 5-homogeneous.


TITLE: Volume vs. Surface Area Integrals
QUESTION [4 upvotes]: In order to find the volume of a sphere radiud $R$, one way is to slice it up into a stack of thin, concentric disks, perpendicular to the $z$-axis. a disk at any point $z$ will have radius $r=\sqrt{R^2-z^2}$ and infinitesimal thickness $dz$, so the volume integral is:
$$ V = \int dV = \int_{-R}^{R} \pi r^2 dz = \int_{-R}^{R} \pi (R^2-z^2) \,dz $$
The above integral of course evaluates to $\frac{4}{3}\pi r^3$, the desired answer. Now, when I tried to do the same thing for the surface area, I treated it as a thin, hollow shell, sliced it up into a stack of concentric rings (or simply the outer surface area of the disks). The area integral in that case is:
$$ A = \int dA = \int_{-R}^{R} 2\pi r \,dz = \int_{-R}^{R} 2\pi \sqrt{R^2-z^2} dz $$
Unfortunately, this integral does not evaluate to $4\pi r^2$. I discovered that instead of taking the thickness of the thin rings as $dz$, I have to take the infinitesimal arclength $dL$ at point $z$ on the cross-section $x = \sqrt{R^2-z^2}$
$$ dL = \sqrt{\left(\frac{z}{\sqrt{R^2-z^2}} \right)^2 +1} \,dz = \frac{R}{\sqrt{R^2-z^2}} dz $$
Putting that in place of $dz$, the integral evaluates perfectly. Conversely, if the thickness in the volume integral was $dL$, in won't work out.
This baffles me. Why do the disks have to have different thickness when you're taking the volume, compared to when you're taking the surface area? Why does it matter? When you take the limit as the number of disk slices goes to infinity and the thickness goes to zero, the sum will approach the sphere, won't it?
Is this a case where the sum of those volumes will approach the sphere, but the sum of their surface areas approaches some different shape, or doesn't approach anything if the thickness is $dz$ and reverse if the thickness is $dL$? If so, why?
The same thing applies for any arbitrary surface/solid of evolution. The surface area integral always involves $dL$, but the volume only has $dx$

REPLY [1 votes]: Try doing what you did to find the circumference of a circle. It's like you approximated the circumference with a bunch of vertical line segments. But no matter how many you make and how small you make them, they still just have a total lenght of $4R$. They're not doing a good job of approximating the curve. To do that, you need to make the ones closer to the top and bottom slanted and therefore longer.

The total length of these segments is just $4R$, not $2\pi R$, and even using more of them at a smaller size would not change this. Try to see the next higher dimension analogy.
Now, if you used the rectangles that these segments define (by adding in horizontal segments) you do get a good approximation of the area, which does get better as the height of these segments gets smaller. The missing area omitted by such rectangles will approach zero as $n\to\infty$. But again, this is not the case with the lengths, which are a constant $4R$ regardless of $n$.


TITLE: if AD ≥ BC, BC ≥ AD then BC = AD
QUESTION [1 upvotes]: I thought it could be solved this way:
(i)$ad-bc \geq 0$ and $bc-ad \geq 0$
Therefore, when adding it, we are left with the inequality
$(ad - bc) + (bc-ad) \geq 0$
when ordering, we are left with $0$, therefore $0=0$
This proves that 0=0, which we already knew, but it does not show that $BC=AD$.
I thought to use antisymmetric but I don't know how it would look

REPLY [0 votes]: In Real Analysis, given $r,s \in \Bbb{R},$ either $r < s, r= s,$ or $r > s.$
Further, $r > s \iff (r - s) > 0$
and $r \geq s \iff (r - s) \geq 0.$
So, you have that $r - s \geq 0.$
Suppose that $s - r > 0.$
Then, you would have that $0 = (r - s) + (s - r) > 0.$  This is a contradiction.
Therefore, it is not the case that $s - r > 0.$
However, you know that $s - r \geq 0.$  Therefore, you can conclude that $s - r = 0 \implies s = r.$

The problem can be attacked a different way.
Since $r \geq s$, it is not the case that $r < s.$
Since $s \geq r,$ and it is not the case that $s > r$, it must be the case that $s = r.$


TITLE: Conjugation of $3$-cycles in $A_5$
QUESTION [6 upvotes]: How do you show that any two 3-cycles are conjugates in $A_5$.

I know we have to take $2$ $3$-cycles say $a$ and $b$ in $A_5$ then we have to show there exists a $c\in A_5$ such that
$a = c b c^{-1}$. But I don't know where to begin. Help me please!

REPLY [7 votes]: Let $(x_1x_2x_3)$ and $(y_1y_2y_3)$ represent the two $3$-cycles. Assume that $x_4,x_5$, respectively $y_4,y_5$ are the remaining two symbols in the five-element set $S_5$ operates on. There is a permutation $\tau$ sending $x_i$ to $\tau x_i=y_i$ for $i=1,...,5$. This will be either a permutation with positive sign or with negative sign. But if it's of negative sign, i.e. an element of $S_5\setminus A_5$, then we can compose $\tau$ with a transposition $(x_4x_5)$ to get an element in $A_5$. Now $\tau(x_1x_2x_3)\tau^{-1}=(τx_1\ τx_2\ τx_3)=(y_1y_2y_3)$

REPLY [0 votes]: You could check it by hand! If you have never done such a thing, it can be instructional to concretely calculate with permutations. Take a particular $3$ cycle, for instance $(123)$, and caculate all $c(123)c^{-1}$ with $c$ running in $A_5$. That would no more than 60 calculations, which shouldn't take more than an hour if you carry them out entirely by hand, maybe even less than half an hour.
You can speed up the process if you know (or prove) that for any cycle $c=(a_1\:a_2\:\cdots a_p)\in S_n$ and $\sigma\in S_n$,
$$\sigma c\sigma^{-1}=\big(\sigma(a_1)\:\sigma(a_2)\:\cdots \sigma(a_p)\big)$$
Actually, once you know this formula, you can create, for every distinct $a,b,c\in\lbrace 1,2,3,4,5\rbrace$, a permutation $\sigma=\sigma_{a,b,c}\in S_5$ with $\sigma(123)\sigma^{-1}=(abc)$, just take
$$\sigma_{a,b,c}=\begin{pmatrix}
1&2&3&4&5\\
a&b&c&d&e
\end{pmatrix}$$
or
$$\sigma_{a,b,c}'=\begin{pmatrix}
1&2&3&4&5\\
a&b&c&e&d
\end{pmatrix}$$
where $d,e$ are the two elements in $\lbrace 1,2,3,4,5\rbrace\smallsetminus\lbrace a,b,c\rbrace$
Now, at first sight, this doesn't solve the question for conjugacy classes in $A_5$. After all, the two choices I gave for $\sigma$ lie in $S_5$, and you want conjugacy in $A_5$. But if you take a closer look at how the two choices of conjugating element are related, you'll find that this actually answers your question. If you want I can add details tomorow, but you should be try to work this out by yourself.

More generally, the conjugacy classes of elements in $A_n$ are the same as those in $S_n$, except for those elements of $A_n$ that are products of even cycles of distinct lengths (fixed points being cycles of length one) : the conjugacy classes of such elements in $S_n$ split into two conjugacy classes in $A_n$.
For instance the $5$-cycles are conjugate in $S_5$ and $S_6$, yet they aren't in $A_5$ and $A_6$, for a $5$-cycle has a cycle decomposition as a $5$-cycle in $A_5$, and a $5$-cycle in $A_6$ has cycle structure a $5$-cycle plus a $1$-cycle. $3$-cycles are conjugate in $A_5$ because their cycle structure is one $3$-cycle plus two $1$-cycles


TITLE: Riemann Integrals
QUESTION [0 upvotes]: If I have a function giving different values when rational or irrational, P is a portition of [0,1]. How do i find the upper and lower Riemann integrals to prove if f is Riemann integrable ?

REPLY [1 votes]: For simplicity, let us assume that $f:[0,1] \to \mathbb{R}$ is such that $f(x)=1$ if $x \in [0,1] \cap \mathbb{Q}$ and $f(x)=0$ if $x \in [0,1] \cap \mathbb{I}$ (here $\mathbb{I}:=\mathbb{R}\backslash \mathbb{Q}$). Since $\mathbb{Q}$ and $\mathbb{I}$ are both dense in $[0,1]$ we have $\displaystyle \sup_{x \in [0,1]}f(x)=1$ and $\displaystyle \inf_{x \in [0,1]}f(x)=0$, so $$\overline{\int_0^1} f(x)dx=1$$ and $$\underline{\int_0^1}f(x)dx=0.$$ In partcular, we see that $f$ is not Riemann-integrable since the value of the upper Darboux-integral does not match the value of the upper Darboux-integral.


TITLE: Trying to find a converse (or counterexample of it) to implicit function theorem
QUESTION [2 upvotes]: The implicit function theorem stated in Spivak's Calculus of Manifolds is as follows:

My question is:

if the above $f$ is continuously differentiable on an open set of
  $(a,b)$ and $f(a,b)=0$
and moreover $\exists$ open set $A$ containing $a$ and open set $B$
  containing $b$ and a continuously differentiable function $g$ such
  that for each $x\in A$ there is a unique point $g(x) \in B$ such that
  $f(x,g(x))=0$
is it necessarily true that
$M$ is invertible?

I ask this because I was talking with my professor about this (and maybe I misheard) but he said that the converse to implicit function theorem is easy to see.
However, I do not know where to begin with this.

REPLY [2 votes]: If $M$ fails to be invertible, then you certainly are in trouble with the theorem. Consider the simplest example: $f(x,y) = x^2+y^2-1$ with $(a,b) = (1,0)$. In no neighborhood of that point is $f=0$ the graph of a function $y=g(x)$.
Now, for your question. If you take $f(x,y) = y^2 = 0$ with $(a,b)=(a,0)$ for any $a$, then we have $g(x)=0$ for all $x$ (which is nicely continuously differentiable) and yet $\partial f/\partial y = 0$ at every point of the subset. So $M$ can fail to be invertible. (This is a standard game to play. Instead of taking the "right" equation $y=0$, you take the "wrong" equation $y^2=0$. This phenomenon can be generalized as much as you want working over the real numbers.)


TITLE: Show that if $X$ is compact, and if $A$ is a closed subset of $X,$ then $A$ is also compact.
QUESTION [0 upvotes]: Show that if $X$ is compact, and if $A$ is a closed subset of $X,$ then $A$(with the subspace topology ) is also compact.
A HINT:
Use the subspace topology to relate open sets in $A$ with open sets in $X;$ somehow you need to find an open cover of $X.$
Still I am unable to solve it using the hints given, could anyone help me in solving it using the hints given?

REPLY [2 votes]: Let $(U_i)_{i\in I}$ be an open cover of $A$, $U_i=V_i\cap A$ where $V_i$ is open in $X$, $(V_i)_{i\in I}\cup (X-A)$ is an open cover of $X$ and has a finite subcover $V_1,...,V_n, X-C$, $A=X\cap A)=(V_1\cup..\cup V_n\cup (X-A))\cap A= U_1\cup..\cup U_n$.