text:
\begin{document}

\title{Point Spread Function Engineering for 3D Imaging of Space Debris using a Continuous Exact $\ell_0$ Penalty (CEL0) Based Algorithm\thanks{The research of the first, third and fourth authors was supported by the US Air Force Office of Scientific Research under grant FA9550-15-1-0286. The work of the second author was supported by HKRGC Grant CUHK14306316, HKRGC Grant CUHK14301718, HKRGC CRF Grant C1007-15G, HKRGC AoE Grant AoE/M-05/12, CUHK DAG No. 4053211, and CUHK FIS Grant No. 1907303.
The work of the third author was also supported by HKRGC Grant CUHK14306316.  }}

\author[1]{Chao Wang}
\author[1]{Raymond H. Chan}
\author[2]{Robert J. Plemmons}
\author[3,4]{Sudhakar Prasad}
\affil[1]{Department of Mathematics, The Chinese University of Hong Kong, Shatin, Hong Kong}
\affil[2]{Departments of Mathematics and Computer Science, Wake Forest University, NC 27109, USA}
\affil[3]{Department of Physics and Astronomy, The University of New Mexico, NM 87131, USA}
\affil[4]{School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA}
\date{ }

\maketitle

\begin{abstract}
 We consider three-dimensional (3D) localization and imaging  of space debris from only one two-dimensional (2D) snapshot image. The technique involves an optical imager that exploits off-center image rotation to encode both the lateral and depth coordinates of point sources, with the latter being encoded in the angle of rotation of the PSF. We formulate 3D localization into a large-scale sparse 3D inverse problem in discretized form.  A recently developed penalty called  continuous exact $\ell_0$ (CEL0)  is applied in this problem for the Gaussian noise model. Numerical experiments and comparisons illustrate the efficiency of the algorithm.

  \end{abstract}

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

The area of 3D imaging and localization

has been getting increasing attention in recent years. The use of 3D localization in single-molecule super-resolution microscopy can obtain a complete picture of subcellubar structures   \cite{Book_SR_micro2017,tilted20183d, 3dsml2017review}. The molecules are labeled by some specific fluorescent proteins or oligonucleotides, which can be regarded as a collection of point sources. Another application of 3D imaging is for space situational awareness (SSA). Currently, there are more than 20,000 objects in orbit around earth \cite{nature_news2018quest}, including operational satellites, dead ones and other human-made debris.  3D localization of micro-scale space debris that become increasingly abundant with decreasing size can be vital for SSA systems responsible for the overall protection of space assets.  Radar systems can sometimes detect such space debris objects, but can at best localize them with lower precision than short-wavelength optical systems.
A stand-alone optical system based on the use of a light-sheet illumination and scattering concept \cite{englert2014optical}
for spotting debris within meters of a spacecraft has also been proposed.
A second system can localize all three coordinates of an unresolved, scattering debris \cite{hampf2015optical,wagner2016detection} by
utilizing either the parallex between two observations, or a pulsed laser ranging system, or a hybrid system.
 However, to the best of our knowledge there is no other proposal of either an optical or an integrated optical-radar system to perform full 3D debris localization and tracking in the range of tens to hundreds of meters.
Prasad \cite{Internal_report2016prasad} has proposed engineering point spread functions for 3D localization by
the use of an optical imager that exploits off-center image rotation.  This system encodes in a single image snapshot both the range $z$ and transverse ($x,y$) coordinates of a swarm of unresolved sources such as
small, sub-centimeter class space debris, which when actively illuminated can scatter a fraction of laser irradiance back into the imaging sensor.  2D image data taken with a specially designed point spread function (PSF) that encodes, via a simple rotation, changing source distance can be employed to acquire a three dimensional (3D) field of unresolved sources like space debris. Here, we propose the 3D localization and tracking of space debris at optical wavelengths by PSF engineering and employing a space-based telescope.

 PSF engineering is widely used in single-molecule super-resolution \cite{lew2011corkscrew,DH2008pavani,DH2009pavani,tetrapods2014shechtman,tetrapods2015shechtman,Rice2016generalized}. It is based on choosing a phase pattern that makes the defocused image of a point source depth-dependent without blurring it excessively. For space-surveillance and SSA, PSF engineering is just beginning to be considered, with \cite{kumar2013psf,prasad2013rotating} proposing the optical theory and simulation of a single lobe rotating PSF. Such PSFs have obvious advantages over multi-lobe PSFs when dealing with high source densities at low light levels. In \cite{kl_nc}, we proposed a mathematical formulation of the 3D localization problem employing such a PSF in the Poisson-noise case. This noise model characterizes an EMCCD sensor operated in the photon-counting (PC) regime. However, when conventional CCD sensors operate at low per-pixel photon fluxes and large read-out noise, a Gaussian noise model describes more accurately this kind of data noise.  In this paper, we consider the latter case.

The following forward model based on the rotating PSF image describes the spatial distribution of image brightness for $M$ point sources in the observed 2D image:
\begin{equation}
	G (x,y) = \sum_{i=1}^M \mathcal{H}_{{z_i}}(x-x_i,y-y_i)f_i + b + \mathcal{N}(x,y),
	\label{equ:forward_model}
\end{equation}
where  $\mathcal{N}$ is the Gaussian noise operator which is data-independent and $b$ is the uniform background value. Here
$\mathcal{H}_{z_i}(x-x_i,y-y_i)$ is the rotating PSF for the $i$-th point source of flux $f_i$ and  3D position coordinates $(x_i, y_i, z_i)$ with the depth information $z_i$ encoded in  $\mathcal{H}_{z_i}$, and $(x, y)$ is  the position in the image plane.
In the Fourier optics model \cite{Goodman17} of image formation, the  incoherent PSF for a clear aperture containing a phase mask that imposes an optical phase retardation, $\psi(\mathbf{s})$, on the imaging wavefront is given by
\begin{equation}
	\label{equ:A}
{\mathcal{H}_{z}}(\mathbf{s} ) =
{1\over \pi}\left|\int P(\mathbf{u} )\mathrm{exp} \left[ \iota( 2\pi \mathbf{u}\cdot\mathbf{s} +  \zeta u^2 - \psi(\mathbf{u}))  \right] d \mathbf{u} \right|^2,
\end{equation}
where
$\zeta = \frac{\pi  (l_0-z) R^2}{\lambda l_0 z}$  is defocus parameter and the imaging wavelength is denoted by $\lambda$.  Here $\iota = \sqrt{-1} $, $\l_0$ is  the distance between the lens and the best focus point, and $P(\mathbf{u} )$ is  the indicator function for the pupil of radius $R$. We use $\mathbf{s}$ with polar coordinates $(s, \phi_{\mathbf{s}}) $ to denote a scaled version of the image-plane position vector, $\mathbf{r}$, namely $\mathbf{s} =  \frac{\mathbf{r}}{\lambda z_I/R} $.
Here $\mathbf{r}$ is measured from the center of the geometric {(Gaussian)} image point located at $\mathbf{r}_I=(x,y)$, and $z_I$ is the distance between the image plane and the lens. The pupil-plane position vector $\boldsymbol{\rho}$ is normalized by the pupil radius, $\mathbf{u} = \frac{\boldsymbol{\rho}}{R}$.
For the single-lobe rotating PSF, $\psi(\mathbf{u}) $ is chosen to be the spiral phase
 profile defined as $$\psi(\mathbf{u}) = l\phi_{\mathbf{u}}, \ \ \text{for } \sqrt{\frac{l-1}{L}}\leq u\leq \sqrt{\frac{l}{L}}, \ l = 1,\cdot\cdot\cdot, L, $$  in which $L$ is the number of concentric annular zones in the phase mask. We evaluate \eqref{equ:A} by using the fast Fourier transform. With such spiral phase retardation, PSF (\ref{equ:A}) performs a complete rotation about the geometrical image center, as $\zeta$ changes between $-L\pi$ and $L\pi$, before it begins to degrade significantly for values of $\zeta$ outside this range.

		Next, we discuss  the problem of 3D localization of closely spaced point sources from simulated noisy image data obtained by using such a rotating-PSF imager. The localization problem is discretized on a cubical lattice where the coordinates and values of its nonzero entries represent the 3D locations and fluxes of the sources, respectively. Finding the locations and fluxes of a few point sources on a large lattice is evidently a large-scale sparse 3D inverse problem. Based on the Gaussian statistical noise model, we describe the results of simulation using a recently developed regularization tool called the continuous exact $\ell_0$ (CEL0) penalty term \cite{CEL02015soubies}, which when added to a least-squares data fitting term constitutes an $\ell_0$-sparsity non-convex optimization protocol with promising results. We use an iteratively reweighted $\ell_1$ (IRL1) algorithm to solve this optimization problem.

The rest of the paper is organized as follows. In Section \ref{sec:non-convex}, we  describe the CEL0-based non-convex optimization model for solving the point source localization problem. In Section \ref{sec:alg}, our non-convex optimization algorithm is developed. Numerical experiments, including comparisons with other optimization methods, are discussed in Section \ref{sec:numerical}.  Some concluding remarks are made in Section \ref{sec:Conclusions}.

\section{CEL0-based Optimization Model}\label{sec:non-convex}

Here, we  build  a forward model for the problem based  on the approach developed in  \cite{Rice2016generalized}.
In order to estimate the 3D locations of the point sources, we assume {that they are distributed on }

a discrete lattice  $\mathcal{X}\in \mathds{R}^{m \times n \times d}$.
The indices of the nonzero entries of $\mathcal{X}$ are the 3-dimensional locations of the point sources and the values at these entries   correspond to the fluxes, {\it i.e.}, the energy emitted by the illuminated point source. The observed 2D image $G \in \mathds{R}^{m \times n}$ can be approximated  as
\begin{equation*}
	G \approx \mathcal{T}(\mathcal{A} \ast \mathcal{X} )+b 1 + \mathcal{N},
\end{equation*}
where
 $b$ is background signal, $1$ is a matrix of $1$s of size the same as the size of $G$ and
$\mathcal{N}$ is the Gaussian noise. Here,   $\mathcal{A} \ast \mathcal{X}$ is the convolution of $\mathcal{X}$ with the 3D PSF $\mathcal{A}$. This 3D PSF ${\mathcal{A}}$ is a cube which is constructed by a sequence of images with respect to different depths of the points. Each horizontal slice is the image corresponding to a point source at the origin in the $(x,y)$ plane and at depth $z$.  This tensor $\mathcal{A}$ {is constructed} by sampling depths at regular intervals in the range, $\zeta_i\in [-\pi L, \ \pi L]$, over which the PSF performs one complete rotation about the geometric image center before it begins to break apart. The $i$-th slice of the dictionary is {$\mathcal{H}_{z_i}$ with certain depth $z_i$}.  Here $\mathcal{T}$ is an operator for extracting the last slice of the cube $\mathcal{A} \ast \mathcal{X}$ since the observed information is a snapshot.

 In order to recover $\mathcal{X}$,  we need to solve a large-scale sparse 3D inverse problem with data-fitting term and regularization term.
Since the Gaussion noise is data-independent, it leads to the use of  least squares for the data-fitting term, {\it i.e.,} $ \frac{1}{2} \left\|\mathcal{T}(\mathcal{A}\ast \mathcal{X})+b 1 -G\right\|_F^2, $
where  $\|Y\|_F$ is the Frobenius norm of $Y$, which is equal to the $\ell_2$ norm of the vectorized $Y$.
For the regularization term, we choose the continuous exact $\ell_0$ (CEL0) penalty, as described in  \cite{cel0_cell2017high,CEL02015soubies,cel0_theory2017unified}.  It is a non-convex term approaching the $\ell_0$
 norm for linear least squares data fitting problems and is constructed as

 $$\mathcal{R}(\mathcal{X}):=\Phi_{\mathrm{CEL0}}(\mathcal{X})=\sum_{u,v,w = 1}^{m,n,d}\phi(\|\mathcal{T}(\mathcal{A} \ast \delta_{uvw} )\|_F, \mu; \mathcal{X}_{uvw}), $$
  where $\phi(a, \mu; u) = \mu  - \frac{a^2}{2}\left(|u|-\frac{\sqrt{2\mu}}{a}\right)^2 \mathds{1}_{\left\{|u| \leq \frac{\sqrt{2 \mu  }}{a} \right\}}$; see Figure~\ref{fig:curve_regul}, and
$\mathds{1}_{\{ u \in E\}} := \begin{cases}
	1 & \text{if} \ u\in E;\\
	0 & \text{others}.
\end{cases}   $
Here $\delta_{uvw}$ is a 3D tensor whose only nonzero entry is at $(u, v, w)$ with value 1 and $\mu $ is the regularization parameter.

\begin{figure}[htbp]
\centering
\includegraphics[width=0.46\textwidth, height=0.27\textheight]{images//cel0_curve.pdf}
\caption{The function $\phi(a, \mu; u)$ for $\ell_2$-CEL0  with $\mu=1$.  }
\label{fig:curve_regul}
\end{figure}

The minimization problem may be stated as
	 \begin{equation}
	 	\min\limits_{\mathcal{X}\geq 0 }\left\{  \frac{1}{2} \left\|\mathcal{T}(\mathcal{A}\ast \mathcal{X})+b-G\right\|_F^2 +\sum_{u,v,w = 1}^{m,n,d}\phi(\|\mathcal{T}(\mathcal{A} \ast \delta_{uvw} )\|_F,\mu; \mathcal{X}_{uvw}) \right\}. \label{equ:min_fun_gaussian}
	 \end{equation}
	 To emphasize that our non-convex optimization model is based on the CEL0 regularization term, we simply designate our optimization model \eqref{equ:min_fun_gaussian} as CEL0.

When combined with a least-squares data fitting term, CEL0 has many good properties and it does not place any strict requirements on the former.

{The global minimizers of the $\ell_0$ penalty model with a least squares data-fitting term ($\ell_2$-$\ell_0$)  are, in fact, contained in the set of global minimizers of CEL0 \eqref{equ:min_fun_gaussian}. A minimizer of \eqref{equ:min_fun_gaussian} can be transformed into a minimizer of  $\ell_2$-$\ell_0$. Moreover, some local minimizers of $\ell_2$-$\ell_0$ are not critical points of CEL0, which means CEL0 can avoid some local minimizers of $\ell_2$-$\ell_0$.
}

\section{ Development of the Algorithm  }\label{sec:alg}

 Note that our optimization model for the Gaussian noise case is  non-convex, due to the regularization term.
We first consider an iterative reweighted $\ell_1$ algorithm (IRL1) \cite{IRL1_2015} to solve the optimization problem. This is a majorization-minimization method which solves a series of convex optimization problems with a weighted-$\ell_1$ regularization term.
 It considers the problem (see Algorithm 3, in \cite{IRL1_2015})
\begin{equation*}
	\min_{x\in X} F(x):=F_1(x)+ F_2(G(x)),
\end{equation*}
 where $X$ is the constraint set, $F$ is a lower semicontinuous (lsc) function, extended, real-valued, proper, while $F_1$ is proper, lower-semicontinous, and convex and $F_2$ is coordinatewise nondecreasing, {\it i.e.} $F_2(x)\leq F_2(x+t e_i)$ with $x, x+t e_i\in G(X)$ and $t >0,$ where  $e_i$ is the $i$-th canonical basis unit vector. The function $F_2$ is concave on $G(X)$.   The IRL1 iterative scheme \cite[Algorithm 3]{IRL1_2015} is
 \begin{equation*}
 	\begin{cases}
 		W^{(k)} = \partial F_2(y), \ y = G(x^{(k)}), \\
 		x^{(k+1)} = \argmin\limits_{x\in X} \left\{F_1(x) + \langle W^{(k)}, G(x)\rangle \right\},
 	\end{cases}
 \end{equation*}
where $\partial$ stands for subdifferential.

For the Gaussian noise problem \eqref{equ:min_fun_gaussian}, we choose
\begin{equation*}
	\begin{split}
		F_1(\mathcal{X}) = & \frac{1}{2} \left\|\mathcal{T}(\mathcal{A} \ast \mathcal{X} )+b 1 -G\right\|_F^2;   \\
		F_2(\mathcal{X}) = & \  \mu  - \frac{a_i^2}{2}\left(\mathcal{X}_{uvw}-\frac{\sqrt{2\mu}}{a_i}\right)^2 \mathds{1}_{\left\{\mathcal{X}_{uvw} \leq \frac{\sqrt{2 \mu  }}{a_i} \right\}}; \\
		G(\mathcal{X}) = & \  |\mathcal{X}|; \\
		X = & \ \{\mathcal{X} \ | \ \mathcal{X}_{uvw} \geq 0 \text{ for all } u,v,w  \}.
	\end{split}
\end{equation*}
Here $a_i = \|\mathcal{T}(\mathcal{A} \ast \delta_{uvw} )\|_F$ and $i = (w-1)mn+(v-1)m+u. $
The algorithm is summarized as Algorithm~\ref{alg:outer}.
\begin{algorithm}[h]
\begin{algorithmic}[1]
\Require $\mathcal{X}^{(0)}\in \mathds{R}^{m\times n \times d}$ and $G\in \mathds{R}^{m \times n}$. Set $ \mu$.
\Ensure The solution $\mathcal{X}^\ast$ which is the minimizer in the last outer iteration.
\Repeat
    \State Compute $W^{(k)}_{uvw} = \left( a_i \sqrt{2\mu} - a_i^2\mathcal{X}^{(k)}_{uvw} \right) \mathds{1}_{\left\{\mathcal{X}^{(k)}_{uvw} \leq \frac{\sqrt{2 \mu  }}{a_i} \right\}};$
    \State Given $G$,  $W^{(k)}_{uvw}$,  obtain $\mathcal{X}^{(k)}$ by solving
    \begin{equation}
        \mathcal{X}^{(k+1)} = \argmin_{\mathcal{X}\geq 0  } \left\{ \frac{1}{2} \left\|\mathcal{T}(\mathcal{A} \ast \mathcal{X} )+b 1 -G\right\|_F^2 + \sum_{u,v,w=1}^{m,n,d} W^{(k)}_{uvw}  |\mathcal{X}_{uvw} | \right\}.\label{equ:subproblem}
    \end{equation}
   \Until{convergence}
\end{algorithmic}
\caption{Iterative reweighted $\ell_1$ algorithm (IRL1) for the rotating PSF problem}
\label{alg:outer}
\end{algorithm}

{\bf Remark: }The minimization problem in \eqref{equ:subproblem} of IRL1 is a weighted  $\ell_1$ model with nonnegative constraints. In  \cite{Rice2016generalized}, the $\ell_1$ model ``without'' nonnegative constraints is solved by the alternating direction method of multipliers (ADMM).

\section{Numerical Results}\label{sec:numerical}

In this section, we apply our optimization  approach to solving  simulated rotating PSF problems for point source localization and compare it to some other competing optimization methods. The codes of our algorithm and the others with which we compared our method were written in $\mathrm{MATLAB \ 9.0 \ (R2016a)}, $ and all our numerical experiments were conducted on a typical personal computer with a standard CPU (Intel i7-6700, 3.4GHz).

The fidelity of localization is assessed in terms of the {\bf recall rate}, defined as the {\it ratio of the number of identified true positive point sources over the number of true positive point sources}, and the {\bf precision rate}, defined as the {\it ratio of the number of identified true positive point sources over the number of all point sources} obtained by the algorithm; see \cite{Book_SR_micro2017}.

To distinguish  true positives from false positives for the estimated point sources, we need to determine the minimum total distance between them and the true point sources. Here all 2D simulated observed images are described by 96-by-96 matrices. We set the number of zones of the spiral phase mask responsible for the rotating PSF at $L=7$  and the aperture-plane side length as 4 which sets the pixel resolution in the 2D image (FFT) plane as 1/4 in units of $\lambda z_I/R$.

The dictionary corresponding to our discretized 3D space contains 21 slices in the axial direction, with the corresponding values of the defocus parameter, $\zeta$, distributed uniformly over the range, $[-21, \  21]$.

According to the Abbe-Rayleigh resolution criterion,  two point sources that are within $(1/2)\lambda z_I/R$ of each other and lying in the same transverse plane cannot be separated in the limit of low intensities.
In view of this criterion and our choice of the aperture-plane side length and if we assume conservatively that our algorithm does not yield any significant super-resolution, we must regard two point sources that  are within 2 image pixel units of each other as a single point source.
Analogously, two point sources along the same line of sight ({\it i.e.,} with the same $x,y$ coordinates) that are axially separated from each other within a single unit of $\zeta$  must also be regarded as a single point source.

As in real problems, our simulation does not assume that the point sources are on the grid points. Rather, a number of point sources are randomly generated in a 3D continuous image space with certain fluxes. We consider a variety of source densities, from 5 point sources to 40 point sources in the same size space.  For each case, we randomly generate 20 observed images and use them for training the parameters in our algorithm, and then test 50 simulated images with the well-selected parameters.
 The number of photons emitted by each point source follows a Poisson distribution with a mean of 2000 photons.

 For adding the Gaussian noise, we use the MATLAB command
  \begin{equation*}
  	\verb|G = I0  + b + sigma*randn(Np)|,
  \end{equation*}
where $\verb|b|$ is the uniform background noise which we set to a typical value 5. Here, $\verb|I0|$ is the 2D original image formed by  adding all the images of the point sources without noise, and  $\verb|Np| = 96$ is the size of the images. The noise level is denoted as $\verb|sigma|$ and we choose  it  to be   10\% of the highest pixel value in original image $\verb|I0|$. Here, \verb|randn| is the $\mathrm{MATLAB}$ command for  the Gaussian distribution with the mean as 0 and standard deviation as 1.

We test our CEL0 based algorithm for several point-source densities. Figure~\ref{fig:15_cel0}-\ref{fig:30_cel0} consider examples of 15 point sources and 30 point sources, respectively.

\begin{figure}[htbp]
\centering
\subfloat[Observed image]{\includegraphics[width=0.3\textwidth]{images//15obs_cel0_1}}
\hspace{0.01mm}
\subfloat[Estimated locations in 2D]{\includegraphics[width=0.3\textwidth]{images//15effect_cel0}}
\hspace{0.01mm}
\subfloat[Estimated locations in 3D]{\includegraphics[width=0.3\textwidth]{images//15x_cel0_1}}
\caption{Localizations for the 15 point sources case. ``$\circ$'' denotes the location of the ground truth point source  and ``+''  the  location of the estimated point source. }\label{fig:15_cel0}
\end{figure}

\begin{figure}[htbp]
\centering
\subfloat[Observed image]{\includegraphics[width=0.3\textwidth]{images//30obs_cel0}}
\hspace{0.01mm}
\subfloat[Estimated locations in 2D]{\includegraphics[width=0.3\textwidth]{images//30effect_cel0}}
\hspace{0.01mm}
\subfloat[Estimated locations in 3D]{\includegraphics[width=0.3\textwidth]{images//30x_cel0}}
\caption{Localizations for the 30 point sources case. }\label{fig:30_cel0}
\end{figure}

From both observed images (see Figure~\ref{fig:15_cel0}(a) and Figure~\ref{fig:30_cel0}(a)), {neither the secondary rings nor the angle of rotation of the PSF are easily identifed}, which means we cannot use a calibration method \cite{huang2008cylendricalstorm,DH2009pavani}. In Figure~\ref{fig:15_cel0}(a), there are two overlapping PSF images corresponding to two different point sources, and our optimization approach is still able to distinguish and estimate them; see the arrows in Figure~\ref{fig:15_cel0}(a),(c).

In Figure~\ref{fig:30_cel0},  many PSF  images are overlapping corresponding to point sources that are very close. Our algorithm estimates the clusters of these point sources but {gives} more point sources than their ground-truth number.

Next, we compare our algorithm with  two other regularization methods: an $\ell_1$ regularization model \cite{Rice2016generalized} as well as a new non-convex model \cite{mila2008nonconvex,mila2010nonconvex,tranformed2014minimization} called  transformed $\ell_1$ (TL1). We use $\ell_1$ to denote the $\ell_1$ regularization model whose regularization term is $$\mathcal{R}(\mathcal{X}):= \mu \left\| \mathcal{X} \right\|_1 = \mu \sum_{u,v,w = 1}^{m,n,d}|\mathcal{X}_{uvw}|. $$
Following \cite{Rice2016generalized}, we solve the optimization problem by ADMM.
 For TL1, the regularization term is $$\mathcal{R}(\mathcal{X}):= \mu \sum_{u,v,w = 1}^{m,n,d} \theta(a;\mathcal{X}) = \mu \sum_{u,v,w = 1}^{m,n,d} \frac{|\mathcal{X}_{uvw}|}{a+|\mathcal{X}_{uvw}|}, $$
where $a$ is fixed parameter and determines the degree of non-convexity. We use IRL1 to solve this model with a similar scheme as Figure~\ref{alg:outer}. The only difference is $$W_{uvw}^{(k)} = \partial F_2(\mathcal{X}^{(k)}_{uvw}) = \frac{a\mu}{\left(a+\mathcal{X}^{(k)}_{uvw}\right)^2}, $$
where $F_2(\mathcal{Y}) = \mu \sum\limits_{u,v,w=1}^{m,n,d} \frac{\mathcal{Y}_{uvw}}{a+\mathcal{Y}_{uvw}}$.

In Figure~\ref{fig:compare_distr_30 cel0},  we again consider the 30 point sources case.  We see that $\ell_2$-$\ell_1$ has more false positives than other algorithms although it detects all the ground truth point sources. TL1 and our algorithm have different but fewer false positives.

\begin{figure}[htbp]
\centering
\subfloat[$\ell_1$]{\includegraphics[width=0.32\textwidth, height=0.22\textheight]{images//30x_l1_1}}
\subfloat[TL1]{\includegraphics[width=0.32\textwidth, height=0.22\textheight]{images//30x_tl1_1}}
\subfloat[CEL0]{\includegraphics[width=0.32\textwidth, height=0.22\textheight]{images//30x_cel0_1}}
\caption{Localizations from three regularization models (30 point sources).  }
\label{fig:compare_distr_30 cel0}
\end{figure}

For further comparisons, we tested 50 different random images and computed the average of recall and precision rates in each density case for both algorithms; see
Figure~\ref{Tab: cel0}.

\begin{table}[htbp]
\begin{center}
\caption{Comparisons of  $\ell_2$-$\ell_1$ with our $\ell_2$-CEL0.  All the results are with post-processing.}
\begin{tabular}{c|cc|cc|cc}
\hline  & \multicolumn{2}{c}{$\ell_1$} & \multicolumn{2}{c}{TL1} & \multicolumn{2}{c}{CEL0}\\
\hline No.  Sources & Recall &  Prec. & Recall &  Prec.   & Recall & Prec.   \\
\hline 10  & 94.80\% & 64.04\% & 89.60\% & 68.79\% & {\bf95.80}\% & {\bf79.72}\% \\
\hline 15  & 90.80\%& 61.68\%  & 87.07\% & 64.67\% &{\bf 93.20}\% &{\bf77.68}\% \\
\hline 20  & 86.60\%& 57.72\% & 83.30\% & 60.78\% & {\bf 89.30}\% &{\bf72.12}\% \\
\hline 30  & {\bf88.80}\% & 47.51\% & 79.80\% & 56.06\% &87.20\% &{\bf58.77}\%  \\
\hline 40  &  {\bf 81.50}\% & 42.03\% & 71.15\% & 48.81\% & 77.40\% &{\bf52.87}\%  \\
\hline
\end{tabular}\label{Tab: cel0}
\medskip
\end{center}
\end{table}

In Figure~\ref{Tab: cel0}, we see that our algorithm is better than $\ell_1$ and TL1  for almost all cases especially in precision rates. For example, in the cases of 10 and 15 point sources, the precision rate in our algorithm is over 10\% higher than the one in $\ell_1$.  In  the higher-density cases, like those with 30 and 40 sources, all methods have more than 5 false positives. We were able to mitigate the latter by further post-processing based on machine learning techniques, as in \cite{Rice2016generalized}. We set the maximum number of iterations for $\ell_1$ at 800, which guaranteed its convergence, and for CEL0 regularization and TL1 we set the maximum number of inner and outer iterations at 400 and 2, respectively.
 Here we  emphasize the advantage of our algorithm in providing a better initial guess than $\ell_1$ and TL1 with a similar cost time.

\section{Conclusions and Future Work}\label{sec:Conclusions}
We have proposed an optimization algorithm, based on a CEL0 penalty term, for the 3D localization of a swarm of randomly spaced point sources using a rotating PSF which has a single lobe in the image of each point source. This has distinct advantages over a double-lobe rotating PSF, e.g.  \cite{moerner2015single,DH2008pavani,DH2009pavani,Rice2016generalized}, especially in cases where the point source density is high and the photon number per source is small. This research focuses on the Gaussian noise case which describes  conventional CCD sensors in low per-pixel photon fluxes and large read-out noise. We note that at high source densities, the optimization can lead to false positives.

We employed a post-processing step based both on centroiding the locations of recovered sources that are tightly clustered  and thresholding the recovered flux values to eliminate obvious false positives from our recovery sources.  These techniques can be applied to other rotating PSFs as well as other depth-encoding PSFs for accurate 3D localization and flux recovery of point sources in a scene from its image data under the Poisson noise model.  Applications include not only 3D localization of  space debris, but also super-resolution 3D single-molecule localization microscopy, e.g.  \cite{Book_SR_micro2017,3dsml2017review}.

Applying recently developed machine learning techniques for removing false positives instead of logistic regression models \cite{Rice2016generalized} will be considered. Tests of this algorithm based on real data collected using phase masks fabricated for both applications are currently being planned.  In addition,  work involving
snapshot multi-spectral/hyperspectral \cite{KK2020hypersp} imaging, which will permit accurate  material characterization, as well as higher 3D
resolution and localization of space microdebris via a sequence of snapshots is under way.

\bibliographystyle{plain}
\bibliography{ref_rPSF}

\end{document}

config: arxiv


TITLE: Spin operators in QM
QUESTION [8 upvotes]: In a text (Introduction to Quantum Mechanics by Griffiths) I am using it states without motivation that spin angular momentum has the same commutations relations as orbital angular momentum (these relations with the ladder operators were used to find the eigenvalue equations of orbital angular momentum) These are the spin angular momentum commutation relations: $$[\hat{S}_x, \hat{S}_y] = i \hbar \hat{S}_z,~~~[\hat{S}_y, \hat{S}_z] = i \hbar \hat{S}_x, ~~~ [\hat{S}_z, \hat{S}_x] = i \hbar \hat{S}_y$$ it follows then that spin angular momentum has the same eigenvalue equations as orbital angular momentum:
$$\hat{S}^2 |s m \rangle = \hbar^2 s (s+1)|s m \rangle;~~~\hat{S}_z |s m \rangle = \hbar m |s m \rangle;$$
Further in the text we consider a two particle system of two spin-$\frac{1}{2}$ particles-for example the electron and proton of a hydrogen atom in ground state where we define the spin operator as $$\hat{S} := \hat{S}^{(1)} + \hat{S}^{(2)}.$$ It then states that in order to confirm eigenvectors belonging to this operator, we have to ensure that the eigenvalue equations above are satisfied. I just want to know if the commutation relations holds for any spin operator, even for multi-particle systems and is that then why the same eigenvalue equations arise? If this is true then we should also be able to apply the ladder operator in each case in about the same way as orbital angular momentum to derive the eigenvectors of a certain spin operator?

REPLY [5 votes]: Coupling two non-interacting quantum systems $\:\alpha,\beta\:$ with angular momenta $\:j_{\alpha},j_{\beta}\:$ (no matter if orbital or spin) we reach to the following equation for the angular momentum $\:j\:$ of the composite system $\:f\:$
\begin{equation}
J_{\boldsymbol{n}}=\bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{n}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\bigr)+ \left(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{n}}\right)
\tag{A-01}
\end{equation}
which, for the $\:\mathbf{n}\:$-components to be more clear, can be expressed as
\begin{equation}
\mathbf{n}\boldsymbol{\cdot}\mathbf{J}=\Bigl[\bigl(\mathbf{n}\boldsymbol{\cdot}\mathbf{J}^{\boldsymbol{\alpha}}  \bigr) \boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}+ \mathbb{I}_{\boldsymbol{\alpha}} \boldsymbol{\otimes} \bigl(\mathbf{n}\boldsymbol{\cdot}\mathbf{J}^{\boldsymbol{\beta}} \bigr)\Bigr]
\tag{A-02}
\end{equation}
On above equations the symbol $''\boldsymbol{\otimes}''$ is used for the product of state vectors, spaces or operators. The vector $\:\mathbf{n}=\left(n_{1},n_{2},n_{3}\right) \:$ is of unit norm. The operators $\:\mathbf{J}^{\boldsymbol{\alpha}},\, J^{\boldsymbol{\alpha}}_{\boldsymbol{n}},\,\mathbb{I}_{\boldsymbol
{\alpha}}\:$ act on the $(2j_{\alpha}+1)$-dimensional Hilbert space $\: \mathsf{H}_{\alpha}\:$ of system $\:\alpha\:$ and on the same footing the operators $\:J^{\boldsymbol{\beta}}_{\boldsymbol{n}},\,\mathbb{I}_{\boldsymbol
{\beta}}\:$ act on the $(2j_{\beta}+1)$-dimensional Hilbert space $\: \mathsf{H}_{\beta}\:$ of system $\:\beta\:$, the symbol $\:\mathbb{I}\:$ being used for the identity. Finally the operators $\:\mathbf{J},\, J_{\boldsymbol{n}}\:$ act on the $(2j_{\alpha}+1)\cdot (2j_{\beta}+1)$-dimensional Hilbert space $\: \mathsf{H}_{f}=\mathsf{H}_{\alpha}\boldsymbol{\otimes}\mathsf{H}_{\beta}\:$ of the composite system $\:f\:$.
We write equation (A-02) for the three axes of a coordinate system
\begin{align}
J_{\boldsymbol{1}}=\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{1}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr)+ \left(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{1}}\right)
\tag{A-03a}\\
J_{\boldsymbol{2}}=\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{2}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr)+ \left(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{2}}\right)
\tag{A-03b}\\
J_{\boldsymbol{3}}=\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{3}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr)+ \left(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{3}}\right)
\tag{A-03c}
\end{align}
These three component equations can be expressed symbolically in one vector equation
\begin{equation}
\mathbf{J}=\bigl(\mathbf{J}^{\boldsymbol{\alpha}} \boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\bigr)+\left(\mathbb{I}_{\boldsymbol{\alpha}} \boldsymbol{\otimes} \mathbf{J}^{\boldsymbol{\beta}} \right)
\tag{A-04}
\end{equation}
Now we must check if this so constructed quantity $\:\mathbf{J}=\left({J}_{1},{J}_{2},{J}_{3}\right) \:$ of the composite system is a consistent angular momentum and the criterion for this is the validation of the equation
\begin{equation}
  \mathbf{J}\boldsymbol{\times}\mathbf{J}= i \, \mathbf{J}
\tag{A-05}
\end{equation}
or by components
\begin{align}
J_{\boldsymbol{1}}J_{\boldsymbol{2}}-J_{\boldsymbol{2}}J_{\boldsymbol{1}}= i \, J_{\boldsymbol{3}}
\tag{A-06a}\\
J_{\boldsymbol{2}}J_{\boldsymbol{3}}-J_{\boldsymbol{3}}J_{\boldsymbol{2}}= i \, J_{\boldsymbol{1}}
\tag{A-06b}\\
J_{\boldsymbol{3}}J_{\boldsymbol{1}}-J_{\boldsymbol{1}}J_{\boldsymbol{3}}= i \, J_{\boldsymbol{2}}
\tag{A-06c}
\end{align}
To prove equations (A-06), let find a general expression for $\:J_{\boldsymbol{n}}J_{\boldsymbol{k}}\:$, where $\:J_{\boldsymbol{n}},\,J_{\boldsymbol{k}}\:$ the components of $\:\mathbf{J}\:$ parallel to the unit vectors $\:\mathbf{n}\:$ and $\:\mathbf{k}\:$ respectively. From equation (A-01) and the following multiplication rule
\begin{equation}
\left(\mathrm{A}_{2} \boldsymbol{\otimes} \mathrm{B}_{2}\right)\left(\mathrm{A}_{1} \boldsymbol{\otimes} \mathrm{B}_{1}\right)= \left( \mathrm{A}_{2}\mathrm{A}_{1}\right)  \boldsymbol{\otimes} \left( \mathrm{B}_{2}\mathrm{B}_{1}\right)
\tag{A-07}
\end{equation}
we have
\begin{align}
  J_{\boldsymbol{n}}J_{\boldsymbol{k}} & = \Bigl[\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{n}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr)+ \Bigl(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{n}}\Bigr)\Bigr] \left[\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{k}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr)+ \left(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{k}}\right)\right]
  \nonumber\\
& =\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{n}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr)\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{k}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr)+\Bigl(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{n}}\Bigr)\left(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{k}}\right)
  \nonumber\\
& + \Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{n}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr)\left(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{k}}\right)
+\Bigl(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{n}}\Bigr)\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{k}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr)
\tag{A-08}
\end{align}
so
\begin{equation}
  J_{\boldsymbol{n}}J_{\boldsymbol{k}} =
   \Bigl[\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{n}}J^{\boldsymbol{\alpha}}_{\boldsymbol{k}}\Bigr)\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr]+\Bigl[\mathbb{I}_{\boldsymbol
{\alpha}}\boldsymbol{\otimes}\Bigl( J^{\boldsymbol{\beta}}_{\boldsymbol{n}}J^{\boldsymbol{\beta}}_{\boldsymbol{k}}\Bigr)\Bigr]
+\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{n}}\boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{k}}\Bigr)
+\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{k}}\boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{n}}\Bigr)
\tag{A-09}
\end{equation}
Permutation of $\:n\:$ and $\:k\:$ yields
\begin{equation}
  J_{\boldsymbol{k}}J_{\boldsymbol{n}} =
   \Bigl[\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{k}}J^{\boldsymbol{\alpha}}_{\boldsymbol{n}}\Bigr)\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr]+\Bigl[\mathbb{I}_{\boldsymbol
{\alpha}}\boldsymbol{\otimes}\Bigl( J^{\boldsymbol{\beta}}_{\boldsymbol{k}}J^{\boldsymbol{\beta}}_{\boldsymbol{n}}\Bigr)\Bigr]
+\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{k}}\boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{n}}\Bigr)
+\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{n}}\boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{k}}\Bigr)
\tag{A-10}
\end{equation}
Subtracting (A-10) from (A-09)
\begin{equation}
  J_{\boldsymbol{n}}J_{\boldsymbol{k}}-J_{\boldsymbol{k}}J_{\boldsymbol{n}}=
  \Bigl[\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{n}}J^{\boldsymbol{\alpha}}_{\boldsymbol{k}}-J^{\boldsymbol{\alpha}}_{\boldsymbol{k}}J^{\boldsymbol{\alpha}}_{\boldsymbol{n}}\Bigr)\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr]+\Bigl[\mathbb{I}_{\boldsymbol
{\alpha}}\boldsymbol{\otimes}\Bigl(J^{\boldsymbol{\beta}}_{\boldsymbol{n}}J^{\boldsymbol{\beta}}_{\boldsymbol{k}}- J^{\boldsymbol{\beta}}_{\boldsymbol{k}}J^{\boldsymbol{\beta}}_{\boldsymbol{n}}\Bigr)\Bigr]
\tag{A-11}
\end{equation}
For $\:n=1\:$ and $\:k=2\:$ above equation (A-11) gives
\begin{align}
  J_{\boldsymbol{1}}J_{\boldsymbol{2}}-J_{\boldsymbol{2}}J_{\boldsymbol{1}} & =
  \Bigl[\overbrace{\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{1}}J^{\boldsymbol{\alpha}}_{\boldsymbol{2}}-J^{\boldsymbol{\alpha}}_{\boldsymbol{2}}J^{\boldsymbol{\alpha}}_{\boldsymbol{1}}\Bigr)}^{i \,J^{\boldsymbol{\alpha}}_{\boldsymbol{3}} }\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr]+\Bigl[\mathbb{I}_{\boldsymbol
{\alpha}}\boldsymbol{\otimes} \overbrace{\Bigl(J^{\boldsymbol{\beta}}_{\boldsymbol{1}}J^{\boldsymbol{\beta}}_{\boldsymbol{2}}- J^{\boldsymbol{\beta}}_{\boldsymbol{2}}J^{\boldsymbol{\beta}}_{\boldsymbol{1}}\Bigr)}^{i \,J^{\boldsymbol{\beta}}_{\boldsymbol{3}}}\Bigr]
\nonumber\\
& = \Bigl[\Bigl(i \,J^{\boldsymbol{\alpha}}_{\boldsymbol{3}}\Bigr)\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol{\beta}}\Bigr]
+\Bigl[\mathbb{I}_{\boldsymbol{\alpha}}\boldsymbol{\otimes}\Bigl(i\,J^{\boldsymbol{\beta}}_{\boldsymbol{3}}\Bigr)\Bigr]
\nonumber\\
& = i\,\Bigl[\Bigl(J^{\boldsymbol{\alpha}}_{\boldsymbol{3}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol{\beta}}\Bigr)
+\Bigl(\mathbb{I}_{\boldsymbol{\alpha}}\boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{3}}\Bigr)\Bigr]
\nonumber\\
& = i \, J_{\boldsymbol{3}}
\tag{A-12}
\end{align}
so proving (A-06a). By cyclic permutation (A-06b) and (A-06c) are proved also.
For the treatment of the angular momentum we make use of equation (A-03c), repeated here for convenience:
\begin{equation}
  J_{\boldsymbol{3}}=\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{3}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr)+ \left(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{3}}\right)
\tag{A-03c}
\end{equation}
This relation has the advantage that if the matrices representing the components $\:J^{\boldsymbol{\alpha}}_{\boldsymbol{3}}\:$ and  $\:J^{\boldsymbol{\beta}}_{\boldsymbol{3}}\:$ of the component systems are diagonal, then  the matrix representing the component $\:J_{\boldsymbol{3}}\:$ of the composite system is diagonal too(1). But for the full treatment of the angular momentum we need the matrix representing the quantity $\:\mathbf{J}^{\boldsymbol{2}}=J^{\boldsymbol{2}}_{\boldsymbol{1}}+J^{\boldsymbol{2}}_{\boldsymbol{2}}+J^{\boldsymbol{2}}_{\boldsymbol{3}}\:$ also. We'll find an expression of $\:\mathbf{J}^{\boldsymbol{2}}\:$ convenient for the determination of its matrix, which isn't from the beginning diagonal as $\:J_{\boldsymbol{3}}\:$ does.
So, inserting in equation  (A-09)  the pair of values $\:(n,k)=(1,1)\:$, $\:(n,k)=(2,2)\:$ and  $\:(n,k)=(3,3)\:$ we have respectively
\begin{align}
 J_{\boldsymbol{1}}^{\boldsymbol{2}} =
   \Bigl[\bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{1}}\bigr)^{\boldsymbol{2}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr]+\Bigl[\mathbb{I}_{\boldsymbol {\alpha}}\boldsymbol{\otimes}\bigl( J^{\boldsymbol{\beta}}_{\boldsymbol{1}}\bigr)^{\boldsymbol{2}}\Bigr]
+2\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{1}}\boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{1}}\Bigr)
\tag{A-13a}\\
J_{\boldsymbol{2}}^{\boldsymbol{2}} =
   \Bigl[\bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{2}}\bigr)^{\boldsymbol{2}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr]+\Bigl[\mathbb{I}_{\boldsymbol {\alpha}}\boldsymbol{\otimes}\bigl( J^{\boldsymbol{\beta}}_{\boldsymbol{2}}\bigr)^{\boldsymbol{2}}\Bigr]
+2\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{2}}\boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{2}}\Bigr)
\tag{A-13b}\\
J_{\boldsymbol{3}}^{\boldsymbol{2}} =
   \Bigl[\bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{3}}\bigr)^{\boldsymbol{2}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr]+\Bigl[\mathbb{I}_{\boldsymbol {\alpha}}\boldsymbol{\otimes}\bigl( J^{\boldsymbol{\beta}}_{\boldsymbol{3}}\bigr)^{\boldsymbol{2}}\Bigr]
+2\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{3}}\boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{3}}\Bigr)
\tag{A-13c}
\end{align}
Having in mind that
\begin{align}
 \bigl(\mathbf{J}^{\boldsymbol{\alpha}}\bigr)^{\boldsymbol{2}}  & =\bigl(      J^{\boldsymbol{\alpha}}_{\boldsymbol{1}}\bigr)^{\boldsymbol{2}} +\bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{2}}\bigr)^{\boldsymbol{2}}+\bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{3}}\bigr)^{\boldsymbol{2}} = j_{\alpha}(j_{\alpha}+1)\mathbb{I}_{\alpha}
\tag{A-14}\\
 \bigl( \mathbf{J}^{\boldsymbol{\beta}}\bigr)^{\boldsymbol{2}}   &=\bigl(      J^{\boldsymbol{\beta}}_{\boldsymbol{1}}\bigr)^{\boldsymbol{2}} +\bigl( J^{\boldsymbol{\beta}}_{\boldsymbol{2}}\bigr)^{\boldsymbol{2}}+\bigl( J^{\boldsymbol{\beta}}_{\boldsymbol{3}}\bigr)^{\boldsymbol{2}} = j_{\beta}(j_{\beta}+1) \mathbb{I}_{\beta}
\tag{A-15}\\
  \mathbb{I}_{\alpha} \boldsymbol{\otimes}\mathbb{I}_{\beta} & \equiv
\mathbb{I}_{f}=\text{identity in  }  \mathsf{H}_{f}=\mathsf{H}_{\alpha}\boldsymbol{\otimes}\mathsf{H}_{\beta}
\tag{A-16}
\end{align}
addition of equations (A-13) yields
\begin{equation}
  \mathbf{J}^{\boldsymbol{2}} =\bigl[ j_{\alpha}(j_{\alpha}+1)+ j_{\beta}(j_{\beta}+1) \bigr] \mathbb{I}_{f} +2\sum_{q=1}^{q=3}\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{q}}\boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{q}}\Bigr)
\tag{A-17}
\end{equation}

(1) More precisely : from the definition of the product of operators and given that  $\:J^{\boldsymbol{\alpha}}_{\boldsymbol{3}}\:$ is represented by the $(2j_{\alpha}+1)$-square matrix
\begin{equation}
   J^{\alpha}_{3} =
   \begin{bmatrix}
     j_{\alpha} & 0 & \cdots & 0 \\
     0 & j_{\alpha}-1 & \cdots & 0 \\
     \vdots & \vdots & m_{\alpha} &  \vdots \\
     0 & 0 & \cdots & -j_{\alpha}
   \end{bmatrix}
\tag{foot-01}
\end{equation}
and $\:J^{\boldsymbol{\beta}}_{\boldsymbol{3}}\:$ is represented by the $(2j_{\beta}+1)$-square matrix
\begin{equation}
   J^{\beta}_{3} =
   \begin{bmatrix}
     j_{\beta} & 0 & \cdots & 0 \\
     0 & j_{\beta}-1 & \cdots & 0 \\
     \vdots & \vdots & m_{\beta} &  \vdots \\
     0 & 0 & \cdots & -j_{\beta}
   \end{bmatrix}
\tag{foot-02}
\end{equation}
equation (A-03c) gives that $\:J_{\boldsymbol{3}}\:$ is represented by the following  $(2j_{\alpha}+1)\cdot (2j_{\beta}+1)$-square diagonal matrix
\begin{equation}
  J_{\boldsymbol{3}}=\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{3}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr)+ \left(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{3}}\right)=
\nonumber\\
\end{equation}
\begin{equation}
 \begin{bmatrix}
     \begin{matrix}
   j_{\alpha}+ j_{\beta} & 0 & \cdots & 0 \\
     0 &   j_{\alpha}+ j_{\beta}-1 & \cdots & 0 \\
     \vdots & \vdots & \ddots &  \vdots \\
     0 & 0 & \cdots &   j_{\alpha} -j_{\beta}
\end{matrix} &  &  & \\
      & \begin{matrix}
       j_{\alpha}-1+ j_{\beta} & 0 & \cdots & 0 \\
     0 & j_{\alpha}-1+ j_{\beta}-1 & \cdots & 0 \\
     \vdots & \vdots & \ddots &  \vdots \\
     0 & 0 & \cdots & j_{\alpha}-1-j_{\beta}
\end{matrix} &  &  \\
     & &\ddots &   \\
  &   &  & -j_{\alpha}-j_{\beta}
\end{bmatrix}
\end{equation}
\begin{equation}
\tag{foot-03}
\end{equation}
Example : for $\:j_{\alpha}=\tfrac{1}{2}\:$ and $\:j_{\beta}=1\:$
\begin{equation}
   J^{\alpha}_{3} =
    \begin{bmatrix}
      \begin{array}{cc}
      +\frac{1}{2}&0\\
      &\\
      0&-\frac{1}{2}
      \end{array}
    \end{bmatrix}
\:,\qquad
   J^{\beta}_{3} =
    \begin{bmatrix}
      \begin{array}{ccc}
      +1&0&0\\
      0&0&0\\
      0&0&-1
      \end{array}
    \end{bmatrix}
\tag{foot-04}
\end{equation}
so
\begin{equation}
   \Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{3}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr)=
    \begin{bmatrix}
      \begin{array}{cc}
      +\frac{1}{2}\cdot\mathbb{I}_{\boldsymbol
{\beta}}&0\cdot\mathbb{I}_{\boldsymbol
{\beta}}\\
&\\
      0\cdot\mathbb{I}_{\boldsymbol
{\beta}}&-\frac{1}{2}\cdot\mathbb{I}_{\boldsymbol
{\beta}}
      \end{array}
    \end{bmatrix}
=
\begin{bmatrix}
      \begin{array}{cccccc}
      +\frac{1}{2}&0&0&0&0&0\\
      0&+\frac{1}{2}&0&0&0&0\\
      0&0&+\frac{1}{2}&0&0&0\\
      0&0&0&-\frac{1}{2}&0&0\\
      0&0&0&0&-\frac{1}{2}&0\\
      0&0&0&0&0&-\frac{1}{2}
      \end{array}
    \end{bmatrix}
\tag{foot-05}
\end{equation}
\begin{equation}
   \left(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{3}}\right)=
    \begin{bmatrix}
      \begin{array}{cc}
      1\cdot J^{\boldsymbol{\beta}}_{\boldsymbol{3}}&0\cdot J^{\boldsymbol{\beta}}_{\boldsymbol{3}}\\
      &\\
      0\cdot  J^{\boldsymbol{\beta}}_{\boldsymbol{3}}&1\cdot J^{\boldsymbol{\beta}}_{\boldsymbol{3}}
      \end{array}
    \end{bmatrix}
=
\begin{bmatrix}
      \begin{array}{cccccc}
      +1&0&0&0&0&0\\
      0&0&0&0&0&0\\
      0&0&-1&0&0&0\\
      0&0&0&+1&0&0\\
      0&0&0&0&0&0\\
      0&0&0&0&0&-1
      \end{array}
    \end{bmatrix}
\tag{foot-06}
\end{equation}
Adding (foot-05), (foot-06) we have
\begin{equation}
  J_{\boldsymbol{3}} =\Bigl( J^{\boldsymbol{\alpha}}_{\boldsymbol{3}}\boldsymbol{\otimes}\mathbb{I}_{\boldsymbol
{\beta}}\Bigr)+ \left(\mathbb{I}_{\boldsymbol
{\alpha}} \boldsymbol{\otimes}J^{\boldsymbol{\beta}}_{\boldsymbol{3}}\right)=
    \begin{bmatrix}
      \begin{array}{cccccc}
      +\frac{3}{2}&0&0&0&0&0\\
      0&+\frac{1}{2}&0&0&0&0\\
      0&0&-\frac{1}{2}&0&0&0\\
      0&0&0&+\frac{1}{2}&0&0\\
      0&0&0&0&-\frac{1}{2}&0\\
      0&0&0&0&0&-\frac{3}{2}
      \end{array}
    \end{bmatrix}
\tag{foot-07}
\end{equation}
which after rearrangement of rows and columns becomes
\begin{equation}
 \widehat{J}_{\boldsymbol{3}}    =
  \begin{bmatrix}
      \begin{array}{cccccc}
      +\frac{3}{2}&0&0&0&0&0\\
      0&+\frac{1}{2}&0&0&0&0\\
      0&0&-\frac{1}{2}&0&0&0\\
      0&0&0&-\frac{3}{2}&0&0\\
      0&0&0&0&+\frac{1}{2}&0\\
      0&0&0&0&0&-\frac{1}{2}
      \end{array}
    \end{bmatrix}
\tag{foot-08}
\end{equation}
recognized later on as the direct sum of $\:j_{1}=\tfrac{1}{2}\:$ and $\:j_{2}=\tfrac{3}{2}\:$
\begin{equation}
\boldsymbol{2}\boldsymbol{\otimes}\boldsymbol{3}=\boldsymbol{2}\boldsymbol{\oplus}\boldsymbol{4}
\tag{foot-09}
\end{equation}
a special case of the more general expression of the product space as the direct sum of mutually orthogonal and invariant under SU(2) subspaces
\begin{equation}
(2j_{\alpha}+1)\boldsymbol{\otimes}(2j_{\beta}+1)=\bigoplus_{j=\vert j_{\beta}-j_{\alpha} \vert }^{j=\left(j_{\alpha}+j_{\beta}\right)}(2j+1)
\tag{foot-10}
\end{equation}

For a more detailed treatment see my answers here : Total spin of two spin-1/2 particles.


text:
In many biological systems of practical interest the applicability of the
fluid approximations we have introduced in Section~\ref{Sec:mat} is limited
by several factors: (a) the presence of boundaries in the state space that
are visited with non-negligible probability, (b) the explicit violation of
the density dependent property, (c) the low population levels that make the
fluid approximation inaccurate, and (d) the presence of ``slow components''
that renders the fluid approximation ineffective. The aim of this section
is to introduce a new model that relies on hybrid switching jump diffusion
(HSJD) processes and that is able to give a reliable approximation of
Markov models for reaction networks even in the presence of these factors.

\subsection{From CTMCs to jump diffusion processes\label{ctmc2jd}}

At the boundary of the state space, indeed, some of the state changes that
are possible in the interior are not enabled any longer (e.g. an enzymatic reaction cannot occur if there are no enzyme molecules available). If the Markov chain spends some time at the boundary, it may happen that special behaviors, which do not occur or occur
rarely in the interior, become significant and cause the appearance of
different modes in the joint probability distribution.  Such multi-modal
behavior cannot be captured by a deterministic approximation which at most
can keep track of the mean value of the marginal distributions of each
component (cf. \cite{srivastava}).  Not even the diffusion
approximation given in eq.~\eqref{eq:sde_ctmc} solves the problem since it is defined
only in the interior of the state space and when the boundary is reached
its validity ceases unless suitable boundary conditions that are not
specified in the original literature are imposed.  In \cite{ATPN14} the
problem was solved in the context of density dependent Markov models that
are described by SPNs whose transitions fire with rates proportional to their enabling degrees (i.e., depend on the number of tokens in their input places)
and whose places are all covered by P-invariants (i.e., the state space is
bounded; for more details on the SPN formalism see \cite{BOABCDF95}).  In
particular, in \cite{ATPN14} a jump diffusion approximation has been
introduced that we now recast in the more general setting of
nearly density dependent models.

Let us consider a nearly density dependent Markov chain $X \BN(u)$ and its
normalized version $Z \BN(u)$. $Z \BN$ has state space
$S\BN \subseteq \mathbb{Q}^{r}$ and it is such that the possible values of
some of its components can be bounded by a minimal and/or by a maximal
value.  Such minimal and maximal values are present, for example, in a
chemical reaction model where the number of molecules of a given chemical
compound cannot become negative, or in a population model, where any given
sub--population can neither become negative nor can exceed the total
population size.

The instantaneous transition rates of the normalized
Markov chain $Z \BN$ will be denoted by $p\BN_{x,x+l/N},
x \in \mathbb{Q}^{r}, l \in \mathbb{Z}^{r}$.  We will use the same notation
for the transition rates extended in the natural way to the state space of
the diffusion approximation where $x \in \mathbb{R}^{r}$.  Let us recall
that the natural state space of the diffusion approximation $Y \BN(u)$ is
the convex hull of the state space of the Markov chain,
$H(S\BN) \subset \mathbb{R}^{r}$.  Notice however that $Y \BN$ is defined
only in the interior of such space and up to the first time the process
visit its boundary.  In order to extend the diffusion approximation to include
the boundary of the state space, we introduce a new approximating process
$\tilde{Y}\BN(u)$ which behaves like the diffusion approximation ${Y}\BN$
in the interior of the state space and that, when the boundary is reached,
mimics the jump behavior of the original Markov chain.  We describe the
process $\tilde{Y}\BN(u)$ more formally in the following paragraphs.

In order to identify those components that are actually at their minimal or
maximal values, we define a map $B:
H(S\BN) \rightarrow \mathcal{P}(\{1,\cdots, r\})$, from the state space of
the diffusion to the power set of the set of indexes, which, given a state
$x$, provides the set of indexes of the extremal components.  Hence,
$B(\tilde{Y}^{\BN}(u))$ is the set of indexes of the components that are on
the boundary at time $u$ \footnote{Let us suppose to have a process with two components, $X_t$ and $Y_t$, both bounded between $0$ and $1$. The map $B$ applied to the vector $(X_t,Y_t)$ returns the set of the indexes of the components at the boundary. If at time $t=7$ the process is in $(X_7=0.34,Y_7=1)$, the map $B$ returns $\{2\}$, meaning that the second component, namely $Y_t$, lies on one of the boundaries.}.  Let us notice that if no components of the
actual state $\tilde{Y}^{\BN}(u)$ attain the boundaries,
then $B(\tilde{Y}^{\BN}(u))=\varnothing$ and $\tilde{Y}^{\BN}(u)$ has the
same behavior as the diffusion $Y^{\BN}(u)$.  As soon as
$B(\tilde{Y}^{\BN}(u))$ becomes non-empty, the set of possible state changes
$C_{x}=\{l \in \mathbb{Z}^{r} : l \neq 0 ,\; p_{x,x+l}\BN \neq 0\}$
is split dynamically (depending on the current state $x=\tilde{Y} ^{\BN}(u)$)
into two subsets, $\stella{C}_{x}$ and $\anello{C}_{x}$. The former
contains those state changes that either modify a component which is
extremal or which are such that the corresponding rate depends on an extremal
component, i.e.
\[ \stella{C}_{x}= \left\{ l : \exists \; i \in B(x) \text{ such that either } l_{i} \neq 0 \text{ or } \frac{\partial p \BN_{x,x+l/N}}{\partial x_{i}} \neq 0 \right\},\]
the latter is the complement set, i.e. $\anello{C}_{x}=C_{x}-\stella{C}_{x}$.

As long as the state changes included in $\stella{C}_{\tilde{Y}^{\BN}(u)}$
do not occur, the subsystem made of the components with indexes not
included in $B(\tilde{Y}^{\BN}(u))$ can still be approximated by
diffusion which are analogous to eq.~\eqref{eq:sde_ctmc} except
that the sums are restricted to the state changes in
$\anello{C}_{\tilde{Y}^{\BN}(u)}$.

The events included in $\stella{C}_{\tilde{Y}^{\BN}(u)}$ encode the behavior at the boundary. We keep them discrete and we treat them as a jump process
which is responsible for all the events of the type ``the $i$-th component
leaves the boundary''. The amplitudes and the intensities of the jumps are
formally taken from the original CTMC and depend on the entire state of the
process, i.e., on all its components, no matter whether they are at the
boundary or not.  The approximating jump diffusion $\tilde{Y} ^{\BN}(u)$
which embodies both the fluid evolution and the discrete events solves the
following system of SDEs
\begin{align}\label{eq:Jsde}
d \tilde{Y} ^{\BN}(u) = \hspace{-2mm}
\sum_{l \in \anello{C}_{\tilde{Y} ^{\BN}(u)}} \hspace{-1mm} l f \left( \tilde{Y} ^{\BN}(u),l \right) du
+ \frac{1}{\sqrt{N}} \sum_{l \in \anello{C}_{\tilde{Y} ^{\BN}(u)}} \hspace{-1mm} l \sqrt{f \left( \tilde{Y} ^{\BN}(u),l \right) }dW_{l}(u)
+ \hspace{-2mm}\sum_{l \in \stella{C}_{\tilde{Y} ^{\BN}(u)}} \hspace{-1mm}\frac{l}{N}\;  dM^{\BN}_{l}(u)
\end{align}
where  $M^{\BN}_{l}(u)$ is the counting process that describes how many events with state change $l$ occurred in the time interval $(0,u]$ and whose intensity is given by
\[
\mu_{l}(\tilde{Y} ^{\BN}(u^-))= p\BN_{\tilde{Y} ^{\BN}(u^-),\tilde{Y} ^{\BN}(u^-)+\frac{l}{N}}
\]
which depends on the actual state of the process $\tilde{Y} ^{\BN}$ right
before the jump\footnote{A function $f$ evaluated in $u^-$ is defined as
the left-sided limit $\displaystyle {f(u^-)= \lim_{x\rightarrow u^-}
f(x)}$.}. Let us recall that a counting process is a stochastic process
with positive, integer and increasing values. It is the natural model for
the number of outcomes in a system over time.  As the diffusion given
in \eqref{eq:sde_ctmc}, also \eqref{eq:Jsde} has its proper mathematical
definition in its integral form, where the integrals with respect to the
Brownian motions $W_l$ and with respect to the counting processes
$M^{\BN}_l$ have to be defined in the general It\^o theory of integration
with respect to semimartingales \cite{klebaner,rogers}.

Equation \eqref{eq:Jsde} is a system of equations, one for every component
of the process.  This might seem contradictory with our previous
description according to which only the components not included in
$B(\tilde{Y} ^{\BN}(u))$ are fluidized.  Let us however remark that the
fluid increments in the first sum of \eqref{eq:Jsde} do not affect the
component at the boundary since if $l \in \anello{C}_{x}$ then the $i$th
entry of $l$ is $0$ for any $i \in B(x)$.  On the other hand, a component
that is not at the boundary at a given time is moved by the continuous
compounding of the fluid increments that sums up with the effect of the
jumps.

\subsection{From CTMCs to hybrid switching diffusions with jumps} \label{subsec:hyb}

Many real systems violate the basic assumption that justifies the fluid
approximation, i.e., that events are very frequent and cause very small
state changes.  An example is a Markov chain in which some of the
transition rates have the nearly density dependent form of
\eqref{eq:qnearly} and others do not. It is reasonable to
``fluidize'' the system partially, making the effect of the density
dependent state changes continuous, but keeping the other state changes as
discrete events. Another example is a Markov model whose components have
values ranging in different scales.  For example, when some components of
the system model resources with limited availability.  Even if the rates
are all nearly density dependent, it is natural to consider a family of
initial conditions where the components not related to the resources grow
large with the indexing parameter of the family, while those modeling
resources are kept fixed.  Accordingly, we propose to fluidize only those
components whose increase is not in contrast with the modeling assumptions.
Moreover, if the transition rates follow the law of mass action given
in eq.~\eqref{eq:lme}, the reactions that involve reactants of which there are
only a few are typically slower than those that involve the abundant
species, introducing different time scales in the reaction rates.  A better
approximation is achieved if the components with small values and the
reactions with law rates are kept discrete, while the rest is fluidized.

Different hybrid approaches of this kind have been adopted in the
literature and we list a few here for comparison. A special mention has to
be done to Piecewise Deterministic Markov Processes (PDMPs, cf. \cite{pdmp}
for the mathematical theory and \cite{hybridLimit, pola2003stochastic,
caravagna,menz, verena,luca2010} for applications in computer science and biology)
that, as we shall see, are the ODE based counterpart to our HSJD approach.
Even hybrid approaches based on diffusion have already appeared in the
literature, mainly using switching diffusions (see \cite{yin} for the
mathematical theory and \cite{intep, haseltine, salis, pola2003stochastic}
for applications) but, as far as we know, no treatment of the behavior at
the boundaries has been included into such models. Other general and
popular approaches to hybrid systems where the fluid components are
evolving according to a diffusion have been studied
in \cite{gshs,lyg1,pola2003stochastic} under the name of (Generalized)
Stochastic Hybrid Models. Such models have a different origin (aircraft
dynamics) and have not been introduced as approximations of Markov
chains. In their dynamics, reaching of the boundaries of the state place
has a role, but it is a different one. Moreover, the partitioning of the
state space between fluid and discrete components is static, independent of
the system evolution.  The specific contribution of this paper is to adopt
the jump diffusion paradigm for the partially fluidized system that allows
the fluid components to visit the boundaries and to jump back in the
interior of the state space in almost the same way as the original chain
does.

For the sake of clarity, we will first introduce the process disregarding
the behavior at the boundary. The machinery that allows the visit of the
boundaries will be added later on.

Let us hence introduce a partial fluidization of a family of
$r$-dimentional Markov chains $X \BN$.  The formal requirements on $X \BN$
for the approximation to hold will be made clear after having introduced
the necessary notation.  Intuitively, we need that a subsystem of the
original model containing a subset of its components and a subset of its
possible state changes is nearly density dependent.  This subsystem is
fluidized and all the rest of the system (the ``problematic'' part) is
left discrete.

We assume that we are given the partition of the set of indexes $\{
1,\dots, r \}$ defined by two disjoint subsets: $F$, with cardinality $n$,
and $D$, with cardinality $m$ such that $n+m=r$, and $F \cup D = \{
1,\dots, r \}$.  The components $X \BN_{i}, i \in F,$ will be fluidized and
their initial condition will be proportional to the indexing parameter $N$.
The initial condition for the components $X \BN_{i}, i \in D,$ will remain
instead constant as $N$ increases.  The partition of the components is
supposed to be explicitly given as part of the model itself.  The process
we aim to approximate is the partially scaled Markov chain $Z \BN$ in which
only the components in $F$ are divided by $N$ (and will be approximated by
a diffusion) while the components in $D$ are kept as they are.
Consequently, $Z_{i} \BN = X_{i} \BN$ for all $i \in D$ and
$Z_{j}\BN=\frac{X_{j}\BN}{N}$ for all $j \in F$. Let us denote the state
space of $Z_{j}\BN$ by ${S'}\BN $.  A transition of the chain $X \BN$
causes a state change from state $k$ to state $k+l$ with a rate
$q_{k,k+l}\BN$. The corresponding transition of the partially scaled chain
$Z\BN$ moves the chain from state $x$ to state $x+l'$ with rate
$p\BN_{x,x+l'}=q_{k,k+l}\BN$ where the state change $l'$ is such that
$l'_{i}=l_{i}$ for all $i \in D$ and $l'_{j}=\frac{l_{j}}{N}$ for all
$j \in F$.

The partition of the components requires to partition the possible state
changes $C=\{ l \in \mathbb{Z}^{r} : \; l \neq 0 \text{ and } \exists\; x :
p\BN_{x,x+l'} \neq 0 \}$ into two disjoint sets $C^{F}$ and $C^{D}$, which
distinguish the fluid events that happen continuously by diffusion and the
discrete ones that happen by jumps.  This partitioning is made on the basis
of the following criteria.  Those transitions that modify the discrete
components necessarily cause discrete increments and hence belong to
$C^{D}$.  The events whose rates depend on a discrete component are also
treated as discrete because their rates could be too slow and the number of
involved object could be too few to justify a continuous approximation.
The remaining events, which are with rates in the nearly density dependent
form, are in $C^{F}$.  Formally we have
\[
C^{F}= \left\{ l \in C : \;   l_{i}=0\ \text{ and } \frac{\partial p_{x,x+l'}\BN}{\partial x_{i}} =0, \; \forall i \in D, x \in {S'} \BN\right\}.
\]
Let us remark that the condition $l_{i}=0, \forall i \in D,$ is strictly
necessary for preserving the discreteness of the components in $D$ while
the condition regarding the transition rates is introduced mainly to
increase the quality of the approximation and might be dropped in the case
one needs to simulate the model faster (but obviously with less
accuracy).  The remaining state changes are included in $C^{D}=C-C^{F}$. Let
us stress that the above partition of the state changes is static,
i.e., it does not depend on the current state of the process.

We approximate the partially scaled process $Z\BN$ with the process
$\Upsilon \BN$ whose components evolve with hybrid switching jump
diffusion dynamics described by the following equations:
\begin{align}
d\Upsilon \BN_{j}(u) &= \sum_{l \in C^{F}} l_{j} f\left( \Upsilon\BN (u),l \right) du
+  \sum_{l \in C^{F}} \hspace{-1mm}  \frac{l_{j}}{\sqrt{N}} \sqrt{f \left( \Upsilon \BN(u),l \right) }dW_{l}(u)
+ \hspace{-2mm}\sum_{l \in C^{D}} \hspace{-1mm}\frac{l_{j}}{N}\;  dJ^{\BN}_{l}(u), \quad \quad j \in F, \nonumber\\
d\Upsilon \BN_{i}(u) &= \sum_{l \in C^{D}} \hspace{-1mm} l_{i}\;  dJ^{\BN}_{l}(u), \quad \quad i \in D
\end{align}
where $J\BN_{l}(u)$ is the process counting the events that causes a
partially scaled state change $l'$ (corresponding to a state change $l$ in
the original chain) which occurred in the time interval $(0,u]$ and whose
intensity is given by
\begin{equation}
\nu_{l}(\Upsilon \BN(u^-))= p\BN_{\Upsilon \BN(u^-),\Upsilon \BN(u^-)+l'} \label{intensity}
\end{equation}

In order to include the behavior on the boundary of the state space for the
components in $F$, we proceed as explained in the previous section.  The
set $C^{F}$ is split dynamically (depending on the current state of the
process $x$) into two subsets $\stella{C}^{F}_{x}$ and $\anello{C}^{F}_{x}$
defined as
\[
\stella{C}^{F}_{x}= \left\{ l \in C^{F} : \exists \; i \in B(x) \cap F \text{ such that either } l_{i} \neq 0 \text{ or } \frac{\partial p \BN_{x,x+l'}}{\partial x_{i}} \neq 0 \right\},
\]
and $\anello{C}^{F}_{x}=C^{F}-\stella{C}^{F}_{x}$, respectively.  Finally,
the hybrid switching jump diffusion process $\tilde{\Upsilon}\BN(u)$
accounting both for the hybrid definition of the model and for the boundary
conditions is given as
\begin{align}
d\tilde{\Upsilon} \BN_{j}(u) &= \sum_{l \in \anello{C}^{F}_{\tilde{\Upsilon}\BN(u)}} l_{j} f\left( \tilde{\Upsilon}\BN (u),l \right) du
+  \sum_{l \in \anello{C}^{F}_{\tilde{\Upsilon}\BN(u)}} \hspace{-1mm}  \frac{l_{j}}{\sqrt{N}} \sqrt{f \left( \tilde \Upsilon \BN(u),l \right) }dW_{l}(u)
+ \hspace{-2mm}\sum_{l \in C^{D} \cup \stella{C}^{F}_{\tilde{\Upsilon}\BN(u)}} \hspace{-1mm}\frac{l_{j}}{N}\;  dJ^{\BN}_{l}(u), \quad \quad j \in F, \nonumber\\
d\tilde \Upsilon \BN_{i}(u) &= \sum_{l \in C^{D}} \hspace{-1mm} l_{i}\;  dJ^{\BN}_{l}(u), \quad \quad i \in D,\label{hsjd}
\end{align}
where, again, $J\BN_{l}(u)$ is the process counting the events that cause
a partially scaled state change $l'$ (corresponding to a state change $l$
in the original chain) which occurred in the time interval $(0,u]$ and
whose intensity is given by \begin{equation}
\nu_{l}(\Upsilon \BN(u^-))= p\BN_{\Upsilon \BN(u^-),\Upsilon \BN(u^-)+l'}. \label{intensity2}
\end{equation}

The resulting process is a \emph{hybrid switching jump diffusion} in which
the jumps account both for the jumps of intrinsically discrete
components and for the jumps out of the boundaries of the fluid
components. Let us remark that \emph{hybrid} jump diffusions are not a
generalization of jump diffusion processes. On the contrary they can be
seen as a special case, since the discrete component can simply be
described as a jump-diffusion process that does only have jumps and no
diffusion.

Let us note that in the limit when $N\rightarrow \infty$ the noise in the
first equation of \eqref{hsjd} vanishes and the model becomes a PDMP in
which the discrete components modulates the ODE and the fluid levels have
an influence on the intensity \eqref{intensity2} of the jumps.  The
mathematical theory of hybrid switching diffusions together with the
special case of PDMP can be found in the recent monograph \cite{yin} where
also some hybrid switching jump-diffusions are briefly sketched in the
appendix.

\subsection{Simulation of the proposed hybrid switching jump diffusion} \label{algoritmoAparole}

From (\ref{hsjd}) an Euler scheme based, approximate simulation algorithm
can be derived directly.  Assume that the state of the process is known at
time $u$ and that we aim to simulate the process in a short time interval.
The length of this interval is $\Delta t$ if no discrete state change
occurs in $[u,u+\Delta t]$ and it is shorter otherwise.  First, we check if
a discrete state change occurs in the time interval $[u,u+\Delta t]$.  For
this purpose, we generate an exponential random variable whose rate,
$\nu$, is equal to the sum of the intensities of the involved counting
processes, i.e., $\nu=\sum_{l \in
C^{D} \cup \stella{C}^{F}_{\tilde{\Upsilon}\BN(u)}} \nu_{l}(\Upsilon \BN(u))$.
Let us denote this random variable by $r$.  If $r>\Delta t$ then no
discrete state change occurs in the considered interval.  If $r\leq
\Delta t$ then we have to simulate which counting process generates the
event.  This is done by generating a discrete random variable according to
a discrete probability distribution with mass function
$m_l=\nu_{l}(\Upsilon \BN(u))/\nu$ with $l \in
C^{D} \cup \stella{C}^{F}_{\tilde{\Upsilon}\BN(u)}$. We denote the
resulting random variable by $k$.  In order to simulate the diffusion we
need to generate $|\anello{C}^{F}_{\tilde{\Upsilon}\BN(u)}|$ random
variables according to the standard normal distribution.  Let us denote
these random variables by $W'_l$ with
$l \in \anello{C}^{F}_{\tilde{\Upsilon}\BN(u)}$.  Having generated the
above random numbers the state of the process is updated first according to
\begin{align}
\nonumber
\tilde{\Upsilon} \BN_{j}(u+\min(r,\Delta t)) = & \tilde{\Upsilon}  \BN_{j}(u)+
\sum_{l \in \anello{C}^{F}_{\tilde{\Upsilon}\BN(u)}} l_{j} f\left( \tilde{\Upsilon}\BN (u),l \right) \min(r,\Delta t)
+  \sum_{l \in \anello{C}^{F}_{\tilde{\Upsilon}\BN(u)}}   \frac{l_{j}}{\sqrt{N}} \sqrt{f \left( \tilde \Upsilon \BN(u),l \right)}
W'_{l} \sqrt{\min(r,\Delta t)}
+ \\ & I(r\leq \Delta t) \frac{l_{k}}{N} \quad \quad ~~~~~~~~~~~~~~ j \in F \nonumber\\
\tilde \Upsilon \BN_{i}(u+ \min(r,\Delta t))
= &  \Upsilon \BN_{i}(u) + I(r\leq \Delta t) l_{k} \quad \quad i \in D \nonumber
\end{align}
where $I$ is the indicator function.  Due to the presence of the diffusion,
it is possible that the previous update results in a state in which some
components are smaller than their possible minimal or larger than their
possible maximal values.  If a component is smaller (greater) than its
minimal (maximal) value then we set it to its minimal (maximal) value.  In
other words, we simulate the process in $[u,u+\Delta t]$ as if it was
unbounded and subsequently we test whether the process left its valid state
space.  If $\Delta t$ is sufficiently small then the above simulation
procedure results in traces that reflect properly the behavior of the
original process.  The choice of $\Delta t$ can be made on the basis of the
actual state of the process by taking into account the maximum among the
drifts.  We found that for our examples $\Delta t=\max\left\{ l_{j}
f\left( \tilde{\Upsilon}\BN (u),l \right)~\mbox{with}~
l \in \anello{C}^{F}_{\tilde{\Upsilon}\BN(u)} \right\}^{-1}/100$ is a
reasonable compromise between execution time and precision.  A more
implementation oriented description of the simulation algorithm is provided
in \ref{app:algo}.

Note that with this simulation approach the borders are reached with delay.
For instance, if a component is positive at time $u$ and it is updated to a
negative value at time $u+\Delta t$, then we assume that the process
reached level 0 at time $u+\Delta t$.  It can also happen that we miss the
fact that the process reached a border.  Indeed there is a positive
probability that an unbounded diffusion process is positive at time $u$ and
at time $u+\Delta t$ and it is negative somewhere inside $[u,u+\Delta]$.
With sufficiently small $\Delta t$ this causes insignificant imprecision.
A more precise treatment of the borders, with which the time when the
process reaches the border is estimated more precisely, is possible and we
plan to experiment with it in the future.

Let us remark that the simulation algorithm described above is similar to
the hybrid stochastic simulation methods based on dynamic partitioning
(e.g., \cite{griffith, salis2}).  The main difference is that, while these
hybrid simulation methods were derived heuristically through clever
simulation speed ups, our simulation algorithm arises by applying a simple
Euler discretization scheme to the SDEs \eqref{hsjd} that describe the
approximating HSJD model.


TITLE: Isn't there a contradiction?
QUESTION [0 upvotes]: In special relativity Einstein used Pythagorean theorem for proving Lorenz transformations. But in general relativity we discovered that space-time has curvature near massive objects, so the geometry near them isn't euclidean. The question is, how can we use Pythagorean theorem in special relativity, if general relativity proves that geometry in our world isn't fully Euclidean ? Or am I missing something?

REPLY [0 votes]: The reason that Special Relativity is called "special" is that it only deals with the special case where there are no large masses nearby--and therefore no gravity. With no gravity, spacetime is flat and the Pythagorean theorem holds.


TITLE: How to get from $\frac{n(n+1)}{2}$ to $n + \binom{n}{2}$?
QUESTION [0 upvotes]: How to get from $\frac{n(n+1)}{2}$ to $n + \binom{n}{2}$?
I tried some transforming, but my main problem is that I don't know how to get to $n!$ and $(n-2)!$ from the original term.

REPLY [3 votes]: Here's a combinatorial proof:
Suppose we have n boys and 1 girl in a room and we have to select 2 people out of them.
Then we can do that in $\binom{n+1}{2}=\frac{n(n+1)}{2}$ ways since we have n+1 people to choose from.
Counting in another way. Case 1: There is a girl in our selection. Then we can do that n ways. Case 2: There is no girl in our selection. We can do that in $\binom{n}{2}$ ways since we have n boys to choose from.
Since these two must be equal. We have:
$$\frac{n(n+1)}{2}=n+\binom{n}{2}$$


TITLE: How to figure out this definite integral
QUESTION [1 upvotes]: How to calculate this definite integral?
$$\int_0^1\frac{x\cos x}{\sqrt{1+x^2}}\mathrm{d}x$$
I've been thinking about this for a few days but still have no ideas.

REPLY [0 votes]: I got an answer with Kampé de Fériet Function.
\begin{aligned}
&\int_0^1{\frac{x\cos x}{\sqrt{1+x^2}}}\text{d}x\\
&=\int_0^{\ln \left( 1+\sqrt{2} \right)}{\text{sinh}x\cos \left( \text{sinh}x \right)}\text{d}x  \ \ \ \ \ (x\rightarrow \mathrm{sinh}x)
\\
&=\sum_{n=0}^{\infty}{\frac{\left( -1 \right) ^n}{\left( 2n \right) !}}\int_0^{\ln \left( 1+\sqrt{2} \right)}{\text{sinh}^{2n+1}x\text{d}x}
\\
&=-\frac{\sqrt{\pi}}{2}\sum_{n=0}^{\infty}{\frac{\Gamma \left( n+1 \right)}{\Gamma \left( 2n+1 \right) \Gamma \left( n+\frac{3}{2} \right)}}\\
&\ \ \ -\sqrt{\frac{\pi}{2}}\sum_{n=0}^{\infty}{\sum_{k=0}^{\infty}{\frac{\left( -1 \right) ^n2^k}{\Gamma \left( 2n+1 \right) \Gamma \left( \frac{3}{2}+k \right) \Gamma \left( \frac{3}{2}+n \right) \Gamma \left( \frac{3}{2}+k+n \right)}}}\ \ \ \ \ (\text{with the help of MMA})
\\
&=-\frac{\pi}{2}\mathbf{L}_{-1}\left( 1 \right) -\sqrt{\frac{\pi}{2}}\sum_{n=0}^{\infty}{\sum_{k=0}^{\infty}{\frac{\left( -1 \right) ^n2^k}{\Gamma \left( 2n+1 \right) \Gamma \left( \frac{3}{2}+k \right) \Gamma \left( \frac{3}{2}+n \right) \Gamma \left( \frac{3}{2}+k+n \right)}}}
\\
&=-\frac{\pi}{2}\mathbf{L}_{-1}\left( 1 \right) -\frac{4\sqrt{2}}{\pi}\sum_{n=0}^{\infty}{\sum_{k=0}^{\infty}{\frac{\left( 1 \right) \!_n\left( 1 \right) \!_k}{\left( 1 \right) \!_n\left( \frac{1}{2} \right) \!_n\left( \frac{3}{2} \right) \!_n\left( \frac{3}{2} \right) \!_k\left( \frac{3}{2} \right) \!_{n+k}}}}\frac{\left( -\frac{1}{4} \right) \!^n}{n!}\frac{2^k}{k!}
\end{aligned}
where $\mathbf{L}_{\nu}\left( z \right)$ is modified Struve function and by the definition of Kampé de Fériet Function
$$\displaystyle \mathbf{F}_{o,m,n}^{r,p,q}\left( \left. \begin{array}{c}
 \mathbf{A}_r;\mathbf{a }_p;\mathbf{a }_{q}^{'}\\[5pt]
 \mathbf{B}_o;\mathbf{b }_m;\mathbf{b }_{n}^{'}\\
\end{array} \right|x,y \right) =\sum_{i,j=0}^{\infty}{\frac{\prod_{s=1}^r{\left( A_s \right) \!_{i+j}}\prod_{k=1}^p{\left( a _k \right) \!_i}\prod_{\ell =1}^q{\left( a _{\ell}^{'} \right) \!_j}}{\prod_{v=1}^o{\left( B_v \right) \!_{i+j}}\prod_{t=1}^m{\left( b _t \right) \!_i}\prod_{u=1}^n{\left( b _{u}^{'} \right) \!_j}}}\frac{x^i}{i!}\frac{y^j}{j!}$$
where $\mathbf{a }_p=\alpha _1,\alpha _2,...,\alpha _p.$ Hence we have
$$\int_0^1{\frac{x\cos x}{\sqrt{1+x^2}}}\text{d}x=-\frac{\pi}{2}\mathbf{L}_{-1}\left( 1 \right) -\frac{4\sqrt{2}}{\pi}\mathbf{F}_{1,3,1}^{0,1,1}\left( \left. \begin{array}{c}
 1;1\\
 \dfrac{3}{2};1,\dfrac{1}{2},\dfrac{3}{2};\dfrac{3}{2}\\
\end{array} \right|-\frac{1}{4},2 \right) $$


TITLE: why cannot we say squareroot of $-1$ as $-i$
QUESTION [3 upvotes]: When $i^2=-1$ why can't we write $i=\pm \sqrt{-1}$ that is why can't we say $$\sqrt{-1}=-i$$

REPLY [0 votes]: For the purposes of traditional analysis, I don't know if it actually even matters which root you define $i$ to be.  The properties of that quantity are all the same regardless of whether $i=-\sqrt{-1}$ or $i=\sqrt{-1}$.  For example
$i^{0+4k}=1 \quad k=0,1,2,...$
$i^{1+4k}=i \quad k=0,1,2,...$
$i^{2+4k}=-1 \quad k=0,1,2,...$
$i^{3+4k}=-i \quad k=0,1,2,...$
$\exp(i\theta) = \cos(\theta) + i\sin(\theta)$
The quaternions give an example of THREE quantities $i^2=j^2=k^2=-1$.  In Clifford algebras there are actually many such general quantities when squared give -1.


TITLE: Finding the Big-O notation with logarithms involved.
QUESTION [0 upvotes]: I am trying to understand algorithms and especially the Big-O notation and I came across this question that included logs:
The question wants me to prove that the Big-O notations for $\log_3 x$ is $\log_2 x$.
I know how to solve the ones with no logs, but I am very confused on how to approach the ones with logs. Thank you so much in advance!
Do I start by breaking it down?
$\log_3 x$ is $O(\log_2 x)$, where $3$ and $2$ are the log bases.

REPLY [3 votes]: Hint: Use the change of base formula for logs
\begin{align}
\frac{\log_2 y}{\log_2 3} = \log_3 y.
\end{align}


TITLE: What physical properties of materials trigger the capacitive touch screen?
QUESTION [5 upvotes]: I have watched some youtube videos about capacitive touch screen
http://www.youtube.com/watch?v=ZHzaVzYEZbw
http://www.youtube.com/watch?v=GmCE18RMEeQ
The videos show that a conductor(finger), in which charges can be induced, can be used to trigger the capacitive touch screen. However, I tried to wear a plastic glove to hold a metal rod to touch the screen and no detection is triggered.
I really want to know the physics theory about the capacitive sensing.
Can someone tell me what physical properties or shape trigger the capacitive touch screen?
Is there any examples that something WITHOUT conducting to human body can trigger the screen?
(P.S. I have read other posts about capacitive touch screen here and do not find a good answer so this is not a duplicate)

REPLY [3 votes]: Capacitive touch screens use a grid of transparent conductors on the inside surface of the glass to sense position.  Each conductor has some series resistance connected so that the capacitance from each conductor to ground can be measured by applying a current pulse to that conductor and measuring the rise and fall times of the voltages at these conductors.  By scanning the x conductors one at a time, the point of contact can be interpolated between the conductors with the slowest rise and fall times to get a precise y position, and the y conductors are similarly scanned to get the x position.  Because capacitance is proportional to the area of intersection between two plates (one being the conductor on the back of the glass, the other being the user's finger), pointed objects, even if good conductors, do not make good surrogates for fingers.  Try the same experiment with a rod of larger diameter (best if it's close to the size of a fingertip, but some touchscreens respond to smaller contact area than others), or some conductive material that deforms as it is pressed against the glass to cover more area, and it should work.  The conductivity of the material is much less important than the area it covers.


TITLE: Does the line passing through $(3,4,-1)$ which is normal to $x+4y-z = -2$, intersect any of the coordinate axes?
QUESTION [0 upvotes]: Does the line passing through $(3,4,-1)$ which is normal to $x+4y-z = -2$, intersect any of the coordinate axes?

I'm not sure how to go about this question. Any help would be greatly appreciated. Thank you in advance!

REPLY [3 votes]: The line passing through $(3,4,-1)$ and normal to the plane $x+4y-z = -2$ has the following parametric equation:
$$\begin{pmatrix}x\\y\\z\end{pmatrix}=\begin{pmatrix}3\\\color{blue}{4}\\\color{blue}{-1}\end{pmatrix}+\lambda\begin{pmatrix}1\\\color{blue}{4}\\\color{blue}{-1}\end{pmatrix}$$
It intersects a coordinate axis if two coordinates are (simultaneously) zero. Can you find a value of $\lambda$ for which this happens? Hint: look at the blue coordinates.


TITLE: Inequality of arithmetic means of $a_1,a_2,\dots,a_n$ and $a_1,2a_2,\dots,na_n$
QUESTION [0 upvotes]: Today while reading some text I stumbled upon an inequality which basically boiled down to the validity of the following:

If $0\le a_1\le a_2 \le \dots \le a_n$ then
  $$\frac{a_1+a_2+\dots+a_n}{a_1+2a_2+\dots+na_n}\le\frac{n}{\frac{n(n+1)}2}.\tag{1}$$
  The above can be also rewritten as
  $$\frac{a_1+a_2+\dots+a_n}{a_1+2a_2+\dots+na_n} \le \frac{n}{1+2+\dots+n},\tag{1'}$$
  maybe this looks a bit nicer.

I managed to prove it in some way. (See below - I hope I have not made a mistake there.) But I would still be interested to see whether there are some other approaches.

My attempt. We can see that the above inequality is equivalent to
$$(n+1)(a_1+a_2+\dots+a_n) \le 2(a_1+2a_2+\dots+na_n). \tag{2}$$
This inequality already seems to be intuitively clear. (The same number of summands on both sides, but on the right hand side we have larger summands more often.)
But it can be also obtained simply by adding the following two inequalities:
\begin{align*}
a_1+2a_2+\dots+na_n &\le a_1+2a_2+\dots+na_n\\
na_1+(n-1)a_2+\dots+a_n &\le a_1+2a_2+\dots+na_n
\end{align*}
where the second inequality is a special case of rearrangement inequality.

REPLY [2 votes]: From what we are given, since $a_1,a_2,\ldots,a_{k-1} \leq a_k$ we have
$$a_1 + a_2 + \cdots +a_k \leq k a_k$$
Hence,
$$\sum_{k=1}^n \left(\sum_{l=1}^k a_l\right) \leq \sum_{k=1}^n k a_k$$
Interchanging the sum on the left hand side gives us
$$\sum_{k=1}^n (n+1-k)a_k \leq \sum_{k=1}^n ka_k$$
Hence,
$$(n+1) \sum_{k=1}^n a_k \leq 2 \sum_{k=1}^n k a_k$$