text:
\begin{document}

\title{Reducing Basis Mismatch in Harmonic Signal Recovery via Alternating Convex Search}

\author{Jonathan M. Nichols*, Albert K. Oh, Rebecca M. Willett \\

\thanks{J. M. Nichols* (jonathan.nichols@nrl.navy.mil) is with the U. S. Naval Research Laboratory, Washington, D.C. 20375. A. K. Oh (akoh2@wisc.edu) and R. M. Willett (rmwillett@wisc.edu) are with the Department of Electrical and Computer Engineering, University of Wisconsin-Madison, Madison, WI 53706. Code associated with this paper is available at: \url{https://sites.google.com/site/albertktoh/software/acs}.}
}

\maketitle

\begin{abstract}

The theory behind compressive sampling pre-supposes that a given sequence of observations may be exactly
represented by a linear combination of a small number of basis vectors. In practice, however, even small deviations from an exact signal model can result in dramatic increases in estimation error; this is the so-called ``basis mismatch'' problem.  This work provides one possible solution to this
problem in the form of an iterative, biconvex search algorithm.  The approach uses standard $\ell_1$-minimization to find the
signal model coefficients followed by a maximum likelihood estimate of the signal model. The algorithm is illustrated on harmonic signals of varying sparsity and outperforms the current state-of-the-art.

\end{abstract}

\begin{IEEEkeywords} Basis mismatch, compressive sampling, sparsity, biconvex optimization, alternating convex search
\end{IEEEkeywords}

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

One of the most fundamental tasks in signal processing is the estimation of the coefficients associated with a given signal model.  This problem frequently takes the form of a linear inverse problem ${\bf y=\Psi x+\bm{\eta} }$ where ${\bf y}\in\R^N$ are the observed data, ${\bf x}\in\R^N$ are the model coefficients and $\bm{\eta} \in\R^N$ is a ``noise'' vector, chosen from some joint probability distribution, usually Gaussian.  A good model adheres to the principle of parsimony, allowing $S<N$ non-zero coefficients to accurately describe the observations, in which case we say the model is $S$-sparse \cite{Donoho:09}.  For real-valued signals, a common signal model ${\bf \Psi}\in\R^{N\times N}$ consists of sine and cosine vectors (i.e. the Fourier basis), in which case the Discrete Fourier Transform (DFT) of the data yields the estimates ${\bf \hat{x}}$.

The problem becomes more complex, however, if we seek to estimate a high-fidelity model using low-fidelity data.
The compressed sensing (CS) framework introduced by Donoho \cite{Donoho:06}, Candes {\it et al.} \cite{Candes:06a}, and others considers the sampling model
\begin{align}
{\bf y}={\bf A\Psi x}+{\bm\eta}
\label{eqn:1}
\end{align}
where the matrix ${\bf A}\in \R^{M\times N}$ projects the high-dimensional signal ${\bf z=\Psi x}$ onto the (potentially low-dimensional) observations ${\bf y}\in\R^M$, and $\bm{\eta} \in\R^M$.  For the undersampled case, where $M<N$,  conditions for accurate estimation are well-established and mandate that $S<\rho M$ where the fraction
$\rho\in (0,1)$ depends on the degree of undersampling, $\delta=M/N$ \cite{Donoho:10}.  The greater the degree of undersampling, the smaller $\rho$ must be for accurate solutions.  If $\rho$ is small enough, and ${\bf A}$ is chosen appropriately (e.g., entries are random draws from some probability distribution), then (\ref{eqn:1}) can be solved via
\begin{align}
\hat{\bf x}= \underset{\bf x}{\arg \min} \ \|{\bf y-A\Psi x}\|_2^2+\lambda\|{\bf x}\|_1,
\label{eqn:L1}
\end{align}
where $\lambda$ is a positive constant that penalizes non-sparse solutions and where $\|{\bf x}\|_1:=\sum_{k=1}^{N} |x_k|$ is large.   The initial theoretical work on this estimator \cite{Candes:05, Donoho:06} has given rise to new approaches to super-resolution \cite{Candes:12}, analog-to-digital conversion \cite{Tropp:10, Nichols:11, Yenduri:12}, and imaging \cite{Romberg:08,lustig2007sparse,willett2014sparsity,harmany2011spatio}.

The key to accurate estimation is the sparsifying transformation matrix ${\bf \Psi}$, chosen so that only $S<M$ elements of the solution ${\bf x}$ are non-zero.  For example, the sum of $S$ harmonics can be described by exactly $S$ complex Fourier coefficients ($2S$ real coefficients) of the aforementioned Fourier basis.  However, because we observe discrete signals, only certain
frequencies can be represented exactly.  The Fourier basis for an $N$-point signal (sampled at unit intervals) is comprised of sine and cosine vectors at discrete frequencies $f_k=(k-1)/N,~k=1,\ldots,  N/2 + 1$.  So long as the signal frequencies match {\it exactly} those in this set,
the representation (Fourier basis) indeed yields a $S$-sparse representation.  However, if the signal frequencies lie off this frequency grid (e.g., $f=(k+0.5)/N$ for $k$ integer), some number $J>S$ frequencies will be required in the representation. Unfortunately, $J$ may be quite a bit larger than $S$, precluding accurate reconstruction of ${\bf x}$ with only $M = S/\rho$ measurements \cite{Nichols:13}.

This is the crux of the so-called ``basis mismatch'' problem, whereby even a good signal model can yield a poor reconstruction due to seemingly small differences between basis vectors assumed by the reconstruction algorithm and a similar set yielding a far sparser representation of the signal.  This problem has received recent theoretical
\cite{Calderbank:11,Gilbert:12,Tang:12} and experimental \cite{Ekanadham:11,Nichols:13,Rao:12} treatment and solutions are still under development.  One straightforward approach to the basis mismatch problem is to simply oversample the frequency space (e.g., let ${\bf \Psi} \in \R^{N \times QN}$ contain sinusoids with frequencies 1/QN apart instead of 1/N apart for a large integer Q).  Though a higher resolution frequency discretization can shrink errors, it comes at the cost of an
increasingly underdetermined problem with stronger coherence between the vectors in the signal model, both of which work against the conditions needed for successful recovery.

Another approach is to treat both the model coefficients {\it and} the associated frequencies as unknowns to be solved. Boufounos {\it et al.} developed an approach that isolates the unknown frequencies in separate, non-overlapping bins and then solves for their
location and amplitude, though the method is restricted to specific sampling strategies \cite{Gilbert:12}.  For more general sensing matrices, Tang {\it et al.} have shown that both frequency locations and amplitudes can be estimated by
solving the constrained ``atomic norm minimization'' problem via semi-definite programming (SDP) \cite{Tang:12}.  Rao {\it et al.} also formulate the problem as one of atomic norm minimization but propose a greedy ``forward-backward'' (GFB) algorithm as the solution \cite{Rao:12}.  The forward-backward approach lacks the theoretical guarantees of the SDP, however it runs in faster time for larger problems.

In this work, we present an Alternating Convex Search (ACS) algorithm \cite{Gorski:07} that significantly reduces the influence of basis mismatch on the recovery of a linear combination of sinusoids.  The approach iterates between a constrained minimization using the familiar $\ell_1$-norm, and a component-wise minimization
involving the model vectors themselves.  Rather than relying on a random gridding of the frequency space, as was done in \cite{Rao:12}, we estimate the frequency vectors that minimize the mean squared error between the data and model at each step. The method is computationally fast relative to SDP, although slightly less accurate, and appears more inline with the approach described in \cite{Rao:12} as well as the work by Valley and Shaw \cite{Valley:13}.

\section{Algorithm Description \label{sec:algorithm}}

Consider the real harmonic continuous time signal $z(\tau)=\sum_{k=1}^{S} x_k \exp(i2\pi f_k \tau)$ with frequencies $f_k\in [0,1/2]$ and complex coefficients given by $x_k$.  We can accurately model $N$ values of this signal as ${\bf z=\Psi x}$ taking ${\bf \Psi}$ to be the $N\times N$ discrete Fourier basis.  As our sampling model, we consider the aformentioned compressed sensing (CS) framework and write ${\bf y}={\bf Az}+{\bm\eta}={\bf A\Psi x}+{\bm\eta}$ where the entries of the $M\times N$ matrix ${\bf A}$ are chosen in accordance with CS theory to be independent draws from some probability distribution \cite{Donoho:10}.  As we have mentioned, the estimator (\ref{eqn:L1}) will degrade (possibly fail) if the signal frequencies are not chosen from the set $f_k=(k-1)/N,~k=1,\ldots, N/2+1$ (see e.g., \cite{Nichols:13}).

Our solution is to treat the frequencies in the signal model as unknowns to be estimated.  Rather than restrict ourselves to the Fourier basis, we consider the more general overcomplete harmonic ``dictionary'' where each of the columns are referred to as ``atoms'' \cite{Donoho:01}.  We let $Q$ be a factor specifying the degree of overcompleteness, chosen so that $QN/2$ is an integer, and form the $N\times QN$ real-valued matrix ${\bf \Psi}_{\bm\theta}$ where the entries are given by

\small
\begin{align}
&\left({\bf\Psi}_{\bm\theta}\right)_{n,k}=\nonumber \\
&~\left\{\begin{array}{r@{,~}l}
		\sqrt{\frac{2}{N}}\cos\left(2\pi n\left(\frac{k-1}{QN}+\theta_k\right)\right) & k=1,\ldots,  \frac{QN}{2} \\
		-\sqrt{\frac{2}{N}}\sin\left(2\pi n\left(\frac{N-k}{QN}-\theta_{N-k+1}\right) \right) & k=\frac{QN}{2}+1,\ldots,  QN
		\end{array}\right.\nonumber \\
~&n=0,\ldots,  N-1
\label{eqn:basisparam}
\end{align}
\normalsize
where each of the frequency parameters $\theta_k,~k=1,\ldots,  QN/2$ can take on values in the range $[-\frac{1}{2QN},\frac{1}{2QN}]$. Although such a representation provides continuous coverage over the frequency space, for $S\ge 6$ (three or more harmonics) the $f_k$ must be spaced at least $1/(QN)$ apart in order to be resolved.  Note also that though there are $QN$ indices, there are only $QN/2$ elements of $\bm{\theta}$ corresponding to the $QN/2$ unique frequencies in the parameterized basis set. This has implications for the implementation of our approach, as any change to the ``positive'' frequencies must be accompanied by an equal and opposite change to the ``negative'' frequencies (sidestepping complex exponentials this way allows us to leverage off-the-shelf sparse optimization tools for real-valued signals and bases.).
The end result is $QN/2$ unknown frequency perturbations that must be estimated.

Clearly if $Q=1$ and all of the
${\bm\theta}$ are zero, then (\ref{eqn:basisparam}) defines the standard Fourier basis.  However, in the general case we modify (\ref{eqn:L1}) and
propose solving the following minimization problem
\begin{align}
\left(\hat{\bf x},{\hat{\bm\theta}}\right)&=\underset{{\bf x},{\bm\theta}}{\arg \min} \ f({\bf x},{\bm\theta}) \equiv \|{\bf y-A\Psi}_{\bm\theta}{\bf x}\|_2^2+\lambda\|{\bf x}\|_1,
\label{eqn:cost}
\end{align}
which simply states that we would like to find both a sparse signal model ${\bf \Psi}_{\bm\theta}$ by estimating the ${\bm\theta}_k$ and the associated model coefficients ${\bf x}$.

In the first step we solve the familiar $\ell_2 \mbox{-} \ell_1$ problem
\begin{align}
{\bf \hat{x}}^{(t+1)}=\underset{\bf x}{\arg \min} \ \|{\bf y-A\Psi}_{{\hat{\bm\theta}}^{(t)}}{\bf x}\|^2_2+\lambda\|{\bf x}\|_1,
\label{eqn:cost1}
\end{align}
holding perturbations ${\bm \theta}$ fixed. This is a convex optimization problem that is known to yield unique sparse solutions provided that the problem size is chosen in accordance with Figure 5 of  \cite{Donoho:10}; the solution will be sparse with only $S<QN$ non-zero coefficients.

The next step is to hold the model coefficients ${\bf x}$ fixed, and find a better signal model by solving for the frequency perturbations of non-zero coefficients. Note that we declare a coefficient to be ``non-zero'' if its magnitude exceeds a threshold $\kappa$ and we capture the associated indices in the set $\mathcal{S}^{(t+1)} \equiv \{\ell:|x_\ell^{(t+1)}| \ge \kappa\}$. In principle, one could set $\kappa=0$, in which case the algorithm must minimize over all $QN/2-1$ unknown perturbations ${\bm\theta}$; the cost is a significant increase in CPU time.  On the other hand, too large a threshold risks ``missing'' a frequency component that requires adjusting. The updated signal model is found via
\begin{align}
{\hat{\bm\theta}}^{(t+1)}_{\mathcal{S}^{(t)}} &=\underset{{\bm \theta}_{\mathcal{S}^{(t)}} }{\arg \min} \ \|{\bf
y-A\Psi}_{\bm \theta}{\bf \hat{x}}^{(t+1)}\|_2^2.
\label{eqn:cost2}
\end{align}
As implied by the notation, rather than solving for all unknown frequency parameters, we solve only for ${\bm\theta}_{\mathcal{S}^{(t)}}$, which are the indices of ${\bm \theta}$ in the set $\mathcal{S}^{(t)}$. To solve for \eqref{eqn:cost2}, we leverage that each of the atoms in our collection are nearly linearly independent, so rather than searching the $S$-dimensional frequency space, we divide the problem into a series of $S$ one-dimensional searches; this simplification was also noted in \cite{Stoica:89}.
Thus, for each non-zero frequency index $\ell \in \mathcal{S}^{(t)}$ we solve the one-dimensional problem
\begin{align}
\hat{\theta}_\ell^{(t+1)} =  \underset{\theta_\ell}{\arg \min} \ \|{\bf y-A\Psi}_{{\bm \theta}}{\bf \hat{x}}^{(t+1)}\|_2^2
\label{eqn:cost3}
\end{align}
while holding each of the other parameters fixed to their current values.

We repeat these two stages of the algorithm until a convergence criteria is met, as described in Algorithm \ref{alg:ACS_alg}.

\begin{algorithm}
\caption{ACS Algorithm}
\label{alg:ACS_alg}
\begin{algorithmic}

\STATE Input signal $\bf y$ and sensing matrix $\bf A$. Select $Q$.
\STATE Set $\alpha = 0.1$, $\beta = 0.1$, and $TOL = 10^{-5}$
\STATE Set $t \leftarrow 0 $
\STATE Initialize $\hat{{\bm \theta}}^{(t)} = 0$

\REPEAT
\STATE Compute ${\bf \Psi}_{\hat{\bm \theta}^{(t)}}$ via \eqref{eqn:basisparam}
\STATE $\lambda = \alpha \| ({\bf A} {\bf \Psi}_{{\hat{\bm\theta}}^{(t)}})^T {\bf y} \|_{\infty}$
\STATE $\hat{{\bf x}}^{(t+1)} \leftarrow$ solution of \eqref{eqn:cost1}
\STATE $\kappa = \beta \| {\hat{\bf x}}^{(t+1)} \|_2$
\STATE Create $\mathcal{S}^{(t+1)} \equiv \{\ell:|\hat{x}_\ell^{(t+1)}| \ge \kappa\}$

\FOR {$\ell \in \mathcal{S}^{(t+1)}$}
\STATE ${\hat{\theta}_\ell^{(t+1)}} \leftarrow$ solution of \eqref{eqn:cost3}
\ENDFOR

\STATE $\hat{f}^{(t+1)} \triangleq f(\hat{{\bf x}}^{(t+1)}, \hat{{\bm \theta}}^{(t+1)})$
\STATE $t \leftarrow t+1$
\UNTIL  {$| (\hat{f}^{(t+1)} - \hat{f}^{(t)})/\hat{f}^{(t)} | < TOL$}

\STATE Output $\hat{{\bf x}} = \hat{{\bf x}}^{(t)}$, $\hat{{\bm \theta}} = \hat{{\bm \theta}}^{(t)}$

\end{algorithmic}
\end{algorithm}

In the supplementary material, we show that the function (\ref{eqn:cost3}) is convex on the frequency range $\theta_k\in [-1/(2QN),1/(2QN)]$ so long as $Q\ge 3/2$ and the true frequencies are separated by at least $1/QN$.  Hence, our algorithm is appropriately described as an Alternating Convex Search (ACS) \cite{Gorski:07}. A benefit to this strategy is that we may make use of Theorem 4.5 in \cite{Gorski:07}, which states that for a biconvex function $f({\bf x},{\bm\theta})$, the sequence $f(\{ {\bf x}^{(i)},{\bm\theta^{(i)}} \})_{i\in \N}$ generated by ACS converges monotonically.  Although convergence can only be guaranteed for an overcomplete dictionary, we found that setting $Q=1$ is sufficient and even advantageous in practice (when the true tones are sufficiently well-separated) as the resulting problem is not as underdetermined.  Nonetheless, results will be shown for both $Q=1$ and $Q=3/2$.

An interesting feature of this algorithm is that it will sometimes shift two frequency atoms toward the same value.  For example, if the signal frequency lies near $k/2QN$ for some k, the algorithm may set $\theta_{k-1} = 1/2QN$ and $\theta_k = -1/2QN$, in which case the associated coefficients, $\hat{x}_{k-1}$ and $\hat{x}_k$, will sum to the true value (i.e., the signal energy is split between two atoms lying at or near the same frequency line).  From a sparsity standpoint this is somewhat undesirable as a single coefficient is preferred.  However, the resulting model is perfectly valid and produces accurate results at the cost of only two additional real coefficients.  This effect will be clearly seen in the results that follow (see Figure \ref{fig:coeff3}, second frequency line).

\section{Numerical Experiments \label{sec:experiment}}

To study the proposed algorithm, we construct a harmonic signal with $S/2$ non-zero complex
($S$ real) coefficients.  Specifically we set, for $i=1,\ldots, N$,
\begin{align}
z(\tau_i)=\sqrt{2/N}\sum_{s=1}^{S/2} \cos\left(2\pi \left(f_s+\theta_s\right) \tau_i+\varphi_s\right),
\label{eqn:sigmodel}
\end{align}
where each frequency is given a uniformly distributed phase $\varphi_s\in [0,2\pi)$. The frequencies themselves are chosen $f_k=k/QN+\theta_k$ where $\theta_k\sim Uniform(-1/(2QN),1/(2QN))$ and \\ $k\sim Uniform(\{0,1,\ldots,QN/2-1\})$ (without replacement). The sampling interval was assumed to be $\tau_{i+1}-\tau_i=1$.

We consider $M=128$ projections as the measured signal ${\bf y}$.  From these projections, we seek to infer the $N=256$ values of $\bm{x}$ that define the underlying signal ${\bf z}={\bf \Psi x}$. Our method easily scales to much larger values of $N$ and $M$, but we consider smaller problem sizes here to facilitate comparisons with the much slower SDP.  Each entry of the ${\bf A}$ matrix is assumed known and chosen independently from a standard normal distribution.  The sequence of noise values $\bm{\eta}$ are also chosen independently from a normal distribution with standard deviation $\sigma=\left(\frac{1}{N}\sum_{i=1}^{N} z(\tau_i)^2\right)^{1/2}/(10^{(\mathrm{SNR \ dB}/20)})$, where $\mathrm{SNR \ dB}$ is the signal-to-noise ratio (log base 10 scale); we chose $\mathrm{SNR \ dB}=40$.

We note that the Gradient Projections for Sparse Reconstruction (GPSR) algorithm \cite{gpsr} was used to solve the $\ell_2 \mbox{-} \ell_1$ optimization problem in \eqref{eqn:cost1}. It was also used to solve \eqref{eqn:L1} using the standard Fourier basis as a baseline comparison. Note that the GPSR algorithm for both setups used identical termination criteria and we set the recommended $\lambda = \alpha \| ({\bf A} {\bf \Psi})^T {\bf y} \|_{\infty}$ value. As a first example we set $S=6$; Figure \ref{fig:coeff3} shows the coefficients obtained from ACS, the standard GPSR implementation, and true coefficient values. ACS clearly produces a sparser model and correctly identifies the true frequency locations and corresponding magnitudes.

\newlength{\figwidth}

\setlength{\figwidth}{0.4\linewidth}
\begin{figure}[t]
\begin{center}
\includegraphics[width=3in]{CoefficientComparison3-eps-converted-to.pdf}
 \caption{Real and imaginary components of the identified coefficients $x_k$ given in the frequency ranges $[0,0.5]$ and $(0.5,0)$, respectively. Shown are true coefficient values and those identified by standard GPSR and the ACS algorithm.}
   \label{fig:coeff3}
\end{center}
\vspace{-3ex}
\end{figure}

\subsection{Comparison to Overcomplete Fourier Dictionary}

In order to address this basis mismatch problem, perhaps the simplest solution is to oversample the frequencies by creating an overcomplete Fourier dictionary $\bm{\Psi}_{OC} \in \R^{N \times QN}$, which is not parameterized with respect to a frequency perturbation, but rather includes a finer discretization of the frequency components. The atoms in this dictionary are identical to those defined in (\ref{eqn:basisparam}), excluding the frequency perturbations. Our discretized signal ${\bf z}$ is then modeled to be sparse in this dictionary by writing ${\bf z} = {\bf\Psi}_{OC} {\bf w}$, where the coefficients ${\bf w} \in \R^{QN}$ can be solved via (\ref{eqn:L1}).  This approach is also
implemented with GPSR and will be referred to as OC-GPSR.

Figure \ref{fig:oc_comp} compares the reconstruction performance of ACS and OC-GPSR for two measures of model error and algorithm execution time.  Specifically, we analyze the normalized RMSE, defined as $\|{\bf z}-\hat{{\bf z}}\|_2^2/\|{\bf z}\|_2^2$, as well as a measure of the identified model support.  Denote as $\mathcal{J}_i\equiv \{j: ~|\hat{f}_j-f_i|<\epsilon\}$ the set of indices for which the identified tones are within $\epsilon$ Hz of the $i^{th}$ tone and compute
\begin{align}
\mathrm{err}=\sum_{i=1}^S |x_i-\sum_{j\in\mathcal{J}_i} \hat{x}_j|.
\label{eqn:err1}
\end{align}
This measure captures the ability of the approach to correctly identify the support of the true model (see e.g., \cite{Recht:14}).  In this work we always use $\epsilon=1/(5QN)$, i.e., 1/5 of a frequency bin, as the threshold for closeness.

Perhaps the biggest advantage to the ACS approach is correct identification of the frequency support (small $\mathrm{err}$).  The OC-GPSR will typically place many frequencies near the true frequency, resulting in a low RMSE.  However the resulting model is not a parsimonious description of the data and is clearly over-fitting the model.  In terms of execution time the two methods are comparable for $Q>8$, however for mildly overcomplete dictionaries the OC-GPSR approach is a good deal faster which may be essential for some applications.

\setlength{\figwidth}{0.4\linewidth}
\begin{figure}[t]
\begin{center}
\subfloat{\includegraphics[width=\figwidth]{CompareACStoOvercompleteError-eps-converted-to.pdf}} ~
\subfloat{\includegraphics[width=\figwidth]{CompareACStoOvercompleteL1-eps-converted-to.pdf}} \\
\subfloat{\includegraphics[width=\figwidth]{CompareACStoOvercompleteTime-eps-converted-to.pdf}} ~
   \caption{Comparison between the ACS approach and overcomplete dictionaries in their respective ability to recover harmonic signals of varying sparsity.  Shown are the normalized RMSE, the $\mathrm{err}$ metric, and the algorithm runtime. (Averaged over 50 random realizations).}
   \label{fig:oc_comp}
\end{center}
\vspace{-3ex}
\end{figure}

Finally, while the ACS approach is only guaranteed to be biconvex for $Q \ge 3/2$, we see similar RMSE values as for the $Q=1$ case, suggesting that $Q \ge 3/2$ is not necessary in practice and good results can be obtained by perturbing the traditional Fourier basis.   We have found this to be true for a variety of sensing matrices and problem sizes.  In the results of the next section, we therefore use only ACS with $Q=1$.

\subsection{Random Temporal Sampling Strategy}

As a final experiment, we again consider the signal model (\ref{eqn:sigmodel}), however, we use the random temporal sampling strategy used for SDP \cite{Tang:12} and GFB \cite{Rao:12} to form
the compressed samples.  That is to say, the compressed samples are a random subset of all possible $z(\tau_i)$ elements, without replacement. In this case the
sensing matrix ${\bf A}$ is a binary matrix consisting of a single ``$1$'' in each row, corresponding to the selected column $\tau_j$.  This is precisely the problem studied in
\cite{Tang:12} and \cite{Rao:12}, and hence forms a good basis for a comparison of the methods.
We note that it is for this random sampling architecture that a comparison to the SDP and GFB approaches can be easily made, without having to make nontrivial modifications to the provided codes. On the other hand, ACS easily applies to generic ${\bf A}$.

\setlength{\figwidth}{0.4\linewidth}
\begin{figure}[t]
\begin{center}
\subfloat{\includegraphics[width=\figwidth]{CompareToSDPErr-eps-converted-to.pdf}} ~
\subfloat{\includegraphics[width=\figwidth]{CompareToSDPTime-eps-converted-to.pdf}} ~
   \caption{Expected RMSE (left) and runtimes (right) associated with the ACS, SDP, and GFB as a function of signal sparsity. (Averaged over 50 random realizations).}
   \label{fig:avg2}
\end{center}
\vspace{-3ex}
\end{figure}

Comparisons can be seen in Figure \ref{fig:avg2}.
With regard to recovery error, the SDP produces the most accurate recovery with ACS offering comparable performance.  By contrast, the GFB approach tends to produce considerably larger RMSE. However, with regard to computation time, both GFB and ACS are very fast. By contrast, SDP is quite slow, taking roughly 2 minutes for all sparsity levels.

We note that only no-phase signals were considered in this simulation. This was because GFB did not seem robust to phase changes, which was not trivial to fix. Both ACS and SDP, on the other hand, were shown to be robust to phase changes and remain the two main competitors for an arbitrary signal where phases may be unknown.

In short, ACS offers a good balance between computation time and accuracy.  It runs faster than SDP, while nearly matching it in accuracy.  Moreover, the ACS approach does not require known phases, and is flexible in the sense that it is easily implementable with any valid CS matrix, not just one for random temporal sampling.  ACS also produces better overall performance (RMSE, $\mathrm{err}$) than an overcomplete dictionary.

\section{Conclusions \label{sec:conclusions}}

This work describes an Alternating Convex Search (ACS) algorithm for correcting frequency errors in the signal model used in compressive sampling
applications for real harmonic signals.  The algorithm treats both frequencies and model coefficients as unknowns and uses an iterative approach in developing estimates.
Specifically, the approach uses the familiar GPSR algorithm to update the model coefficients, followed by a maximum likelihood estimate of the unknown
frequency locations. The algorithm was demonstrated effective at recovering harmonic signals possessing varying levels of sparsity in competitive computation times.
Finally, we note that the code associated with this paper is available to download at: \url{https://sites.google.com/site/albertktoh/software/acs}.

\bibliographystyle{IEEEtran}

\bibliography{refs}

\clearpage

\begin{appendix} \label{app:A}

\renewcommand{\theequation}{A\arabic{equation}}
\setcounter{equation}{0}
\setcounter{figure}{0}
\section*{Convexity of Surrogate Cost Function}

In this section, we discuss the convexity of the function (\ref{eqn:cost}).  This is the primary requirement for making statements about algorithm convergence.  In particular, we may make use of the following theorem:
\begin{theorem}
(Theorem 4.5 in \cite{Gorski:07}) Let $B\subset\R^{M}\times\R^{N}$, let $f:B\rightarrow\R$ be a biconvex function which is bounded from below, and let the optimization problems (\ref{eqn:cost1},\ref{eqn:cost2}) be
solvable.  Then the sequence $f(\{ {\bf x}^{(i)},{\bm\theta^{(i)}} \})_{i\in \N}$ generated by ACS converges monotonically \cite{Gorski:07}.
\label{theorem:converge}
\end{theorem}
For fixed ${\bm\theta}^*$ the problem (\ref{eqn:cost1}) is known to be convex in ${\bf x}$.

While the function described by (\ref{eqn:cost2}) is not convex, it can be well-approximated by a function that is convex provide each of the frequency parameter is restricted to the range $\theta_k\in[-1/(2QN),1/(2QN)]$, and where $Q$ is determined in what follows.  To begin, denote the true and estimated signals at a given iteration of the algorithm as ${\bf z}$, where each element is $z_n=\sum_k x_k\sin(2\pi f_k n), n=0\cdots N-1,~f_k\in[0,1/2]$ and $\hat{\bf z}({\bm\theta})$, where each element is $\hat{z}({\bm\theta})_n=\sum_k \hat{x}_k\sin(2\pi (f_k+\theta_k)n), n=0\cdots N-1$ respectively.  We are interested in the cost function
\begin{align}
\|{\bf y}-{\bf A}\hat{\bf z}\|_2^2&=\|{\bf Az}-{\bf A}\hat{\bf z}({\bm\theta})+{\bm \eta}\|_2^2\nonumber \\
&=\|{\bf A}({\bf z}-\hat{\bf z}({\bm\theta}))+{\bm \eta}\|_2^2\nonumber \\
&\equiv \|{\bf Av}_{\bm\theta}+{\bm \eta}\|_2^2,
\label{eqn:app1}
\end{align}
where ${\bf A}$ is an $M\times N$ matrix comprised of i.i.d. entries chosen from some probability distribution
function with zero mean and variance $\sigma_A^2$.  The noise vector ${\bm \eta}$ is an $M\times 1$ vector consisting of
zero-mean, i.i.d. Gaussian entries $\eta_m$, each with variance $\sigma^2$. The $N\times 1$ vector ${\bf v}_{\bm\theta}$ is a measure of the discrepancy between signal and estimate where the subscript signifies this vector is a function of our unknown parameters ${\bm\theta}$.  Each member of this vector will be denoted $v_{n\bm\theta},~n=1\cdots N$.  Similarly, each member of ${\bf A}$ will be denoted $a_{mn}$ and each projection of the difference vector will be written ${\bf a}_m{\bf v}_{\bm\theta}$ where we have denoted ${\bf a}_m$ as the $m^{th}$ row of the matrix ${\bf A}$.  Moreover, it should be noted that each $E[a_{mn}v_{n\bm\theta}]=0$, and $E[(a_{mn}v_{n\bm\theta})^2]=\sigma_A^2 v^2_{n\bm\theta}$ where the expectation is taken with respect to the probability distribution governing the $a_{mn}$.

Expanding (\ref{eqn:app1}) and noting that the minimum of the cost function is not affected by re-scaling, we have
\begin{align}
\frac{1}{M}\|{\bf Av}_{\bm\theta}+{\bm \eta}\|_2^2&=\frac{1}{M}\left\{{\bf v}_{\bm\theta}^T{\bf A}^T{\bf Av}_{\bm\theta}+2{\bm \eta}^T{\bf Av}_{\bm\theta}+{\bm \eta}^T{\bm \eta}\right\}\nonumber \\
\label{eqn:expand}
\end{align}
In what follows we will make use of the central limit theorem in approximating this cost function.  Recall that the central limit theorem states that for $M$ independent realizations, $x_m$, of a random variable $X$ with mean $\mu_X$ and variance $\sigma_X^2$, the summation $\frac{1}{M}\sum_{m}x_m\sim \mathcal{N}(\mu_X,\sigma_X^2/M)$ for $M$ large.  In other words, the discrete summation is equal to the true value plus an error term of order $O(\sigma_X/\sqrt{M})$.

With this in mind, we first note that the last term is simply the noise variance
\begin{align}
\frac{1}{M}\sum_{m=0}^{M-1} \eta_m^2=\sigma^2+O(\sqrt{2}\sigma^2/\sqrt{M})
\end{align}
and provides a constant offset to the cost function.  The order of the error is determined by noting that for a normally distributed random variable $X$, with variance $\sigma_X^2$, the variance of $X^2$ is $2\sigma_X^4$.

The cost function is largely governed by ${\bf v}_{\bm\theta}^T{\bf A}^T{\bf Av}_{\bm\theta}$ which is a non-convex function of the unknown parameter vector ${\bm\theta}\in\R^N$. However, over a convex set $\Theta$ which restricts the range of each $\theta_k$, this function can be well approximated by a function that is convex.   To show this, we first re-write this term as
\begin{align}
\frac{1}{M}{\bf v}_{\bm\theta}^T&{\bf A}^T{\bf Av}_{\bm\theta}=\frac{1}{M}\sum_{m=0}^{M-1}({\bf a}_m{\bf v}_{\bm\theta})^2=\frac{1}{M}\sum_{m=0}^{M-1}\left(\sum_{n=0}^{N-1} a_{mn}v_{n\bm\theta}\right)^2\nonumber \\
&=\frac{1}{M}\sum_{m=0}^{M-1} \left[\sum_{n=0}^{N-1} (a_{mn}v_{n\bm\theta})^2+\sum_{n=0}^{N-1} \sum_{q\ne n}^Na_{mn}a_{mq} v_{n\bm\theta}v_{q\bm\theta}\right]\nonumber \\
&=\sum_{n=0}^{N-1} v_{n\bm\theta}^2 \frac{1}{M}\sum_{m=0}^{M-1} a_{mn}^2+\sum_{n=0}^{N-1} \sum_{q\ne n}^Nv_{n\bm\theta}v_{q\bm\theta}\frac{1}{M}\sum_{m=0}^{M-1} a_{mn}a_{mq} \nonumber \\
&=\sum_{n=0}^{N-1} v_{n\bm\theta}^2 \left\{E\left[a_{mn}^2\right]+O\left(\frac{\sqrt{2+\gamma(A)}\sigma_A^2}{\sqrt{M}}\right)\right\}\nonumber \\
&~~+\sum_{n=0}^{N-1} \sum_{q\ne n}^N v_{n\bm\theta}v_{q\bm\theta} \left\{E\left[a_{mn}a_{qn}\right]+O\left(\frac{\sigma_A^2}{\sqrt{M}}\right)\right\}\nonumber \\
&=\sigma_A^2 \sum_{n=0}^{N-1} v_{n\bm\theta}^2\nonumber \\
&~~ + O\left(\frac{\sqrt{2+\gamma(A)}\sigma_A^2}{\sqrt{M}}\sum_{n=0}^{N-1} v_{n\bm\theta}^2+\frac{\sigma_A^2}{\sqrt{M}}\sum_{n=0}^{N-1} \sum_{q\ne n}^N v_{n\bm\theta}v_{q\bm\theta} \right).
\label{eqn:cross}
\end{align}
The constant $\gamma(X)\equiv \frac{E[(X-\mu)^4]}{\sigma_X^4}-3$ is the kurtosis associated with the random variable $X$ and depends on the probability distribution governing $X$.  For example, if the $a_{mn}$ are normally distributed $\gamma(A)=0$ while for entries chosen from a uniform distribution, $\gamma(A)=-6/5$.   Furthermore we have under the assumption that the values $a_{mn}$ are independently chosen, $E[a_{mn}a_{qn}]=0$ while the variance of the second term in (\ref{eqn:cross}) is the product of the variances of the constituent variables.  Thus, by the central limit theorem, approximating this covariance term with an expectation yields an error of order $\sigma_A^4$.

In short, the degree to which the cost function is approximated by the last line in (\ref{eqn:cross}) is entirely predicated on the behavior of the function
\begin{align}
\sum_{n=0}^{N-1}&v_{n\bm\theta}^2\nonumber \\
&=\sum_{n=0}^{N-1}\left[\sum_{k=1}^S x_k\sin(2\pi f_k n)-\sum_{k=1}^S
\hat{x}_k\sin(2\pi (f_k+\theta_k)n)\right]^2\nonumber \\
&=\sum_{n=0}^{N-1} \left[\sum_{k=1}^S \left[x_k\sin(2\pi f_kn)-\hat{x}_k\sin(2\pi(f_k+\theta_k)n)\right]^2\right.\nonumber \\
&\left.+\sum_{k=1}^S\sum_{k\ne k'}x_kx_{k'}
\sin(2\pi f_kn)\sin(2\pi f_{k'}n)\right.\nonumber \\
&\left.-2\sum_{k=1}^S\sum_{k\ne k'} x_k\hat{x}_k\sin(2\pi f_k n)\sin(2\pi f_{k'}n)\right.\nonumber \\
&\left.+\sum_{k=1}^S\sum_{k\ne k'}\hat{x}_k\hat{x}_{k'}\sin(2\pi
f_kn)\sin(2\pi (f_k+\theta_k)n)\right].
\end{align}
The asymptotic (large $N$) properties of the above are easily obtained by noting that
\begin{align}
\lim_{N\rightarrow\infty} \sum_{n=0}^{N-1}\sin(2\pi f_k n)\sin(2\pi f_{k'}n)= 0.
\end{align}
The result is that
\begin{align}
\lim_{N\rightarrow\infty} &\sum_{n=0}^{N-1}v_{n\bm\theta}^2\nonumber \\
&=\lim_{N\rightarrow\infty}\sum_{k=1}^{N}\sum_{n=0}^{N-1} \left[x_k\sin(2\pi f_k n)-\hat{x}_k\sin(2\pi (f_k+\theta_k) n)\right]^2\nonumber \\
&=\lim_{N\rightarrow\infty} \sum_{k=1}^{N} g_k(\theta_k)
\end{align}
where
\begin{align}
g_k(\theta_k)&=\sum_{n=0}^{N-1} \left[x_k\sin(2\pi f_k n)-\hat{x}_k\sin(2\pi (f_k+\theta_k) n)\right]^2
\label{eqn:app3}
\end{align}
Thus, the $N$-dimensional cost can be approximated by a sum of $N$ one-dimensional functions.  We also have that $g_k(\theta_k)\ge 0$, ensuring that the $N$-dimensional minimum of this function can be found by taking each solution component (each $\theta_k$) as the minimizers of (\ref{eqn:app3}).

Returning to Eqn. (\ref{eqn:cross}) it can be seen that for large $N$ the second error term vanishes due to the near orthogonality of the product of sinusoids at different frequencies.  Hence, for even moderately large $N$ we have that
\begin{align}
{\bf v}_{\bm\theta}^T{\bf A}^T{\bf Av}_{\bm\theta}&=\sigma_A^2 \sum_{k=1}^S g_k(\theta_k) + O\left(\frac{\sqrt{2+\gamma(A)}\sigma_A^2}{\sqrt{M}}\sum_{k=1}^S g_k(\theta_k)\right).
\end{align}

Finally, the second term in (\ref{eqn:expand}) can be written
\begin{align}
\frac{2}{M}\sum_{m=0}^{M-1} ({\bf a}_m{\bf v}_{\bm\theta})\eta_m&=\frac{2}{M}\sum_{m=0}^{M-1}\sum_{n=0}^{N-1} a_{mn}\eta_mv_{n{\bm\theta}}\nonumber \\
&=2\sum_{n=0}^{N-1} v_{n{\bm\theta}} \left\{E[a_{mn}]E[\eta_m]+O\left(\frac{\sigma_A\sigma}{\sqrt{M}}\right)\right\}\nonumber \\
&=2\sum_{n=0}^{N-1} v_{n{\bm\theta}} \times O\left(\frac{\sigma_A\sigma}{\sqrt{M}}\right)\nonumber \\
&=2\sum_{k=1}^S h(\theta_k) \times O\left(\frac{\sigma_A\sigma}{\sqrt{M}}\right)
\end{align}
where we have noted that for independent random variables the expectation of the product factors.  We have also made the same orthogonally argument as above in order to write the cost as a function of
\begin{align}
h_k(\theta_k)=\sum_{n=0}^{N-1} \left[x_k\sin(2\pi f_k n)-\hat{x}_k\sin(2\pi (f_k+\theta_k) n)\right].
\end{align}
Thus, for either zero mean noise, or zero mean entries in the matrix ${\bf A}$, this part of the cost function reduces to an error term of order $\frac{2\sigma_A\sigma}{\sqrt{M}}\sum_{k=1}^S h_k(\theta_k)$.

Incorporating all components of the cost function we therefore have
\begin{align}
\frac{1}{M}&\|{y-\bf A}\hat{\bf z}_{\bm\theta}\|_2^2=
\sigma_A^2\sum_{k=1}^{S} g_k(\theta_k)+\sigma^2\nonumber \\
&+O\left(\frac{\sqrt{2+\gamma(A)}\sigma_A^2}{\sqrt{M}}\sum_{k=1}^{S} g_k(\theta_k)+\frac{\sigma_A\sigma}{\sqrt{M}}\sum_{k=1}^S h_k(\theta_k)+\frac{\sigma^2}{\sqrt{M}}\right)
\label{eqn:app2}
\end{align}
so that the quality of the approximation is seen to be dependent on the functions $g_k(\theta),~h_k(\theta)$.    By inspection the function $h_k(\theta)$ is simply a sinusoid with offset governed by $x_k$ and an amplitude governed by $\hat{x}_k$ hence, the order of this error term is at most $O(|x_k+\hat{x}_k|)$ for any value of $\theta_k$ and for any number of data $N$.  As will be shown next, the function $h_k(\theta)$ is $O(1/N)$ relative to the term $g_k(\theta)$ and can safely be neglected in the analysis.

The approximate cost function and the quality of the approximation are therefore determined by the behavior of $g_k(\theta_k)$.  In describing the behavior of this function we also define the difference between the true amplitude and estimated amplitude as $\delta=|x_k-\hat{x}_k|$.  In this case it can be shown that the function possesses a minimum value of $g_k(0)=\frac{N}{2}\delta^2$ and grows to approach a value of $N|\hat{x}_kx_k|+\frac{N}{2}\delta^2$ for $|\theta_k|\ge 1/2/N$.
\begin{figure}[tbc]
  \centerline{
   \begin{tabular}{c}
        \includegraphics[width=3.0in,angle=0]{Hfunc-eps-converted-to.pdf}
   \end{tabular}
   }
   \caption{Behavior of the function $g_k(\theta)$.  Also shown are the limits of a standard frequency bin width of $\pm 1/2/N$.  We are ultimately interested in the behavior of this function on this interval for a given value of $N$.}
   \label{fig:H}
\end{figure}
This function is shown in Figure (\ref{fig:H}) for a signal consisting of a single frequency and for two different sets of signal amplitudes.  In the first example we set the true amplitudes equal to the values $\hat{x}=x=1$, hence $\delta=0$.  In the second we set $\hat{x}_k=4,~x_k=1$ to simulate a very large error in the estimated coefficient value.  The approximation is clearly better when the estimated and true amplitudes are close.  However, even with a large estimated coefficient error the function $g_k(0)$ is still the minimum and the shape of the function is very much the same as the true cost function.  In fact, the approximation remains good for $\theta_k$ extending well beyond a standard Fourier bin width.

Neglecting the higher-order error terms from (\ref{eqn:app2}), we have approximately that
\begin{align}
\frac{1}{M}&\|{y-\bf A}\hat{\bf z}_{\bm\theta}\|_2^2\approx \sigma^2+\sigma_A^2\sum_{k=1}^{S} g_k(\theta_k)+\nonumber \\
&O\left(\frac{\sqrt{2+\gamma(A)}\sigma_A^2SN(|\hat{x}{x}|+\frac{\delta^2}{2})}{\sqrt{M}}\right), \left[-\frac{1}{2N}\le\theta_k\le \frac{1}{2N}\right]\nonumber \\
&\approx \sigma_A^2\sum_{k=1}^{S} g_k(\theta_k)
\label{eqn:appnew}
\end{align}
where $S$ is the number of tones in the signal.  In short, to understand the properties of the stochastic, non-convex function (\ref{eqn:app1}), we may study the properties of the deterministic function (\ref{eqn:appnew}).

As an illustration of the approximate cost (\ref{eqn:appnew}), we choose each $a_{mn}\sim\mathcal{N}(0,1)$, take $M=N=1024$, and plot both (\ref{eqn:app1}) and (\ref{eqn:appnew}) for a single tone ($S=1$) with true frequency $f=0.2$.  We further assumed the case of zero additive noise, i.e., $\sigma=0$.   This comparison is shown in Figure \ref{fig:costcompare}.
\begin{figure}[tbc]
  \centerline{
   \begin{tabular}{ccc}
        \includegraphics[width=1.75in,angle=0]{CostFigZoomCompare-eps-converted-to.pdf} &\hspace*{0.1in} &
         \includegraphics[width=1.75in,angle=0]{CostFigZoomCompare2-eps-converted-to.pdf}
   \end{tabular}
   }
   \caption{Comparison of cost functions (\ref{eqn:app1}) and (\ref{eqn:app2}) for $N=1024,~M=512$ over a standard frequency bin width, $\theta=[-1/2/N,1/2/N]$.  The entries of $A$ were chosen from a standard normal distribution ($\sigma_A^2=1$).  In the first example (left) we set $\hat{x}=x=1$ so that $\delta=0$.  As predicted, the approximation is best in the vicinity of the true frequency ($\theta=0$), rising to be of order $\sqrt{2}N/\sqrt{M}=63$ near the edges of the bin.  For the second example (right) we set $\hat{x}_1=1$, and $x_1=4$, hence the minimum of the cost function attains the value $\frac{N}{2}\delta^2=4608$ while the error is predicted to be of order $\sqrt{2}\times 8.5N/\sqrt{M}=500$.  Although the discrepancy in estimated and true amplitude drives the cost function to larger values, the minimum remains at the true frequency location, i.e., where $\theta=0$.}
   \label{fig:costcompare}
\end{figure}
We have found the approximation to be accurate to the stated order for a variety of probability distributions governing the entries of ${\bf A}$ and for a variety of sampling parameters $N,~M$.

Given this approximation, we are interested in the range over which (\ref{eqn:appnew}) is convex.  In order to show this we need only demonstrate that each of the $g_k(\theta_k)$ is convex.
Differentiating (\ref{eqn:app3}) twice w.r.t. $\theta_k$, we seek the range over which the function $g''(\theta_k)$ is positive.  Under the large sample approximation we see that the zeros of $g''(\theta_k)$ are given by the zeros of
\begin{align}
g_Z''(\theta_k)&=\lim_{N\rightarrow\infty} g''(\theta_k)=2 N^2 \csc (\pi  \theta_k) \sin (\pi  (2 N+1) \theta_k)\nonumber \\
&+2 N \csc ^2(\pi \theta_k) \cos (2 \pi  (N+1) \theta_k)\nonumber \\
&+3 \csc ^2(\pi \theta_k) \cos (2 \pi  (N+1) \theta_k)\nonumber \\
&+\csc ^3(\pi \theta_k) (-\sin (\pi  (2 N+3) \theta_k))\nonumber \\
&-\cot ^2(\pi  \theta_k)+\csc ^2(\pi  \theta_k)
\label{eqn:secondderiv}
\end{align}
which is independent of the true and estimated amplitudes $x,~\hat{x}$.  In the discrete setting we are interested in determining the frequency bin width over which this expression is guaranteed to be positive.  In the worst case, each unknown frequency is precisely at the midpoint of each bin in our discretized frequency space.  We therefore re-write the continuous parameter $\theta_k$ in terms of a discrete
frequency bin size as $\theta_k=\pm \frac{1}{2QN}$ and seek the value of $Q$ for which (\ref{eqn:secondderiv}) is positive.  This substitution yields the transcendental equation
\begin{align}
\tan\left(\frac{\pi}{|Q|}\right)&=\frac{4|Q|\pi}{4|Q|^2-2\pi^2}
\end{align}
for which we require the roots for $|Q|>1$.  It turns out that the only such root is $|Q|=1.509$, so the cost function (\ref{eqn:appnew}) is convex within $\pm \frac{1}{3.018N}$ of the true frequency $f$.  Hence, in order to ensure that the optimization problem being solved is indeed biconvex and that Theorem (\ref{theorem:converge}) applies, we are required to ``tile'' this space with bins of width $\frac{2}{3N}$ which is tantamount to using a Fourier dictionary that is $\frac{3}{2}$ times overcomplete.  In practice, however, we have found that a standard discrete Fourier basis ($Q=1$) is sufficient to produce good convergence.

\end{appendix}

\ifthenelse{\boolean{FIGSINBACK}}{

}

\end{document}


TITLE: Prove that the solution for $y'=y^3(1-\tan^2(\arcsin(y)))$ , $y(0)= {\pi \over 8}$ , is bounded.
QUESTION [0 upvotes]: I got this problem to prove, and I assume I need to use the existence and uniqueness theorem for non-linear ODE's, so I set $y' = f(x,y)$ and differentiating in respect to $y$ gives: $f_y(x,y)$.
And both of them are continuous at $y=\Large  {\pi \over 8}$.
But from this point on I don't really know what I'm doing, how do I get a solution? And how do I check if it is bounded?
I'd be happy for some explanation and a correct solution, I need to know how to prove similar problems for my upcoming test.
Thanks in advance.

REPLY [1 votes]: Let $\sin^{-1}y=\theta\iff \sin \theta=y$
$$\tan\theta=\frac{\sin \theta}{\cos\theta}=\frac{\sin \theta}{\sqrt{1-\sin^2\theta}}=\frac{y}{\sqrt{1-y^2}}$$
$$\tan(\sin^{-1}y)=\frac{y}{\sqrt{1-y^2}}$$
Now, we have $$y'=y^3\left(1-\tan^2(\sin^{-1}y)\right)$$
setting the value of $\tan(\sin^{-1}y)$,
$$\frac{dy}{dx}=y^3\left(1-\left(\frac{y}{\sqrt{1-y^2}}\right)^2\right)$$
$$\frac{dy}{dx}=y^3\left(\frac{1-2y^2}{1-y^2}\right)$$
$$\left(\frac{1-y^2}{y^3(1-2y^2)}\right)\ dy=dx$$
$$\int \left(\frac{1}{y^3}+\frac{1}{y}+\frac{2y}{1-2y^2}\right)\ dy=\int dx$$
$$\int \frac{1}{y^3}\ dy+\int \frac{1}{y}\ dy-\frac{1}{2}\int \frac{d(1-2y^2)}{1-2y^2}\ dy=\int dx$$


text:
\begin{document}
\bibliographystyle{plain-annote}
\title{Cameron-Liebler line classes}

\author{Morgan Rodgers
  \thanks{This research was supported in part by NSF \#0742434 UCD GK12 Transforming Experiences Project.}\\
  University of Colorado Denver\\
  Department of Mathematical and Statistical Sciences\\
  morgan.rodgers@ucdenver.edu
}
\date{July 1, 2013}

\maketitle

\begin{abstract}

New examples of Cameron-Liebler line classes in $\pg(3,q)$ are given with parameter  $\frac{1}{2}(q^2 -1)$.
These examples have been constructed for many odd values of $q$ using a computer search, by forming a union
of line orbits from a cyclic collineation group acting on the space. While there are many equivalent
characterizations of these objects, perhaps the most significant is that a set of lines $\mathcal{L}$ in $\pg(3,q)$
is a Cameron-Liebler line class with parameter $x$ if and only if every spread $\mathcal{S}$ of the space
shares precisely $x$ lines with $\mathcal{L}$. These objects are related
to generalizations of symmetric tactical decompositions of
$\pg(3,q)$, as well as to subgroups of $\pgl(4,q)$ having equally many orbits on points and lines of $\pg(3,q)$.
Furthermore, in some cases the line classes we construct are related to two-intersection sets in $\ag(2,q)$.
Since there are very few known examples of these sets for $q$ odd, any new results in this direction are of
particular interest.

\end{abstract}

\section{Introduction}\label{intro}
Cameron and Liebler~\cite{CL} studied sets of lines in $\pg(3,q)$ having certain nice properties, today known as
\emph{Cameron-Liebler line classes}.  Such a set contains $x(q^2+q+1)$ lines for some integer $x$, and among many
equivalent properties, shares with every spread of $\pg(3,q)$ precisely $x$ lines; $x$ is called the \emph{parameter}
of the set.  There are some trivial examples of Cameron-Liebler line classes; if we take all of the lines through a
common point, or all of the lines in a common plane, we get a line class with parameter 1.  If we take a
non-incident point-plane pair, and take all of the lines through the point along with all of the lines in the plane,
we get a line class having parameter 2.  Also, it is easy to see that the complement of a Cameron-Liebler line class
with parameter $x$ is a line class with parameter $q^2 +1 -x$.

While it was originally conjectured by Cameron and Liebler that only the trivial examples of these line classes
would exist, this was disproved by an example due to Drudge in~\cite{D}, which was shortly thereafter shown to be
part of an infinite family of examples by Bruen and Drudge~\cite{BD}. These examples exist for all odd values of
$q$, and have $x = \frac{1}{2}(q^2+1)$.  The first counterexample to the conjecture for an even value of $q$ was
given by Goevarts and Penttila in~\cite{GP}, having $q = 4$ and $x = 7$.  There has been much work excluding certain
parameters for Cameron-Liebler line classes;
Govaerts and Storme~\cite{GS} first showed that there are no Cameron-Liebler line classes in $\pg(3,q)$ with
parameter $2 < x \leq q$ when $q$ is prime.  Some time later, De Beule, Hallez, Storme~\cite{BHS} excluded
parameters $2 < x \leq q/2$ for all values of $q$.  Most recently, Metsch~\cite{M} eliminated the possibility of
having $2 < x \leq q$ for any prime power $q$.

In this work, we detail new Cameron-Liebler line classes constructed for many odd values of $q$ satisfying
$q \equiv 1 \bmod{4}$ and $q \not\equiv 1 \bmod{3}$, having parameter $x = \frac{1}{2}(q^2 -1)$.  These new examples
are made up of a union of orbits of a cyclic collineation group having order $q^2 +q +1$.  A Cameron-Liebler line
class with parameter $x$ Klein-corresponds to an $x$-tight set of points in $Q^{+}(5,q)$~\cite{DT}; we can use this
model to easily test if a given set of points in $Q^{+}(5,q)$ corresponds to a Cameron-Liebler line class using an
eigenvector condition.  Among these new examples, when $q \equiv 0 \bmod{3}$ we seem to have some especially
interesting behavior; in these cases, the Cameron-Liebler line class constructed gives rise to a symmetric tactical
decomposition of $\pg(3,q)$, which can be used to derive a set of type $(m,n)$ in the affine plane $\ag(2,q)$.  Our
hope is that these new Cameron-Liebler line classes will belong to an infinite family, and that there will also be
an infinite family of related affine sets of type $(m,n)$.  Affine sets of type $(m,n)$ have previously only been
found in planes of order $9$, so the discovery of an infinite family of new examples would hold particular interest.

\section{Cameron-Liebler line classes (tight sets of $Q^{+}(5,q)$)}\label{sec:1}
\begin{definition}
  A \textbf{Cameron-Liebler line class}  $\mathcal{L}$ is a set of lines in $\pg(3,q)$ such that  any line $\ell$ of
  $\pg(3,q)$ is incident with
  \begin{equation*}
    |\{\ m \in \mathcal{L} : \ \mbox{$m$ meets $\ell$, $m \neq \ell$} \}| =
    \left\lbrace
      \begin{array}{ll}
        (q+1)x + (q^2 -1) & \mbox{ if $\ell \in \mathcal{L}$}\\
        (q+1)x & \mbox{ if $\ell \not\in \mathcal{L}$}
      \end{array}
    \right.
  \end{equation*}
  for a fixed integer $x$, called the \emph{parameter} of $\mathcal{L}$.
\end{definition}
There are many other equivalent characterizations of these sets of lines; for an extensive list see~\cite{Pe}.  The
lines of $\pg(3,q)$ Klein-correspond to points of $Q^{+}(5,q)$~\cite{DT}.  The Klein correspondence is a bijection
(which we shall always denote by $\kappa$) between the set of lines of $\pg(3, q)$ and the set points of
$Q^{+}(5, q)$, such that lines $\ell$ and $\ell^{\prime}$ of $\pg(3, q)$ are incident if and only if $\kappa(\ell)$
and $\kappa(\ell^{\prime})$ are collinear.  In this context, we will be interested in \emph{tight sets} of points.
The notion of an $i$-tight set of a finite generalized quadrangle was introduced by Payne in~\cite{P}, and was
extended to polar spaces of higher rank by Drudge in~\cite{DT}.
\begin{definition}
  A set of points $\mathcal{T}$ in $Q^{+}(5,q)$ is said to be \textbf{$i$-tight} if
  \begin{equation*}
    \left| P^{\perp} \cap \mathcal{T} \right| = \left\lbrace
      \begin{array}{ll}
        i (q+1) + q^{2} & \mbox{ if } P \in \mathcal{T} \\
        i (q+1)             & \mbox{ if } P \not\in \mathcal{T} \mbox{.}
      \end{array}
    \right.
  \end{equation*}
\end{definition}
It is easy to see that a Cameron-Liebler line class of $\pg(3,q)$ with parameter $x$ Klein-corresponds to an
$x$-tight set of $Q^{+}(5,q)$.

We now describe the model of the Klein quadric we will be using.
Let $F = GF(q)$, $E = GF(q^3)$, and $T$ be the relative trace function from $E$ to $F$.   We consider the quadric
$Q^{+}(5,q)$ to have the underlying vector space $V = E^2$ considered as a vector space over $F$, and equipped with
the quadratic form $Q : (x,y) \to T(xy)$.  The polar form of $Q$ is given by
$f((u_{1},u_{2}), \ (v_{1}, v_{2})) = T(u_{1}v_{2}) + T(u_{2}v_{1})$.  This form is nondegenerate, and it can be
seen that $\pi_{1} = \{(x,0) : \ x \in E^{*}\}$ and $\pi_{2} = \{(0,y) : \ y \in E^{*}\}$ are totally singular
planes in the quadric.

We will be using the orbits of a cyclic group to construct our Cameron-Liebler line classes.  Take $\mu \in E^{*}$
with $|\mu| = q^2 +q + 1$.  Define the map $g$ on $Q^{+}(5,q)$ by
$g: \langle (x,y) \rangle \to \langle (\mu x, \mu^{-1} y) \rangle$.  Now $g$ is a projective isometry, and the group
$C = \langle g \rangle$  can be seen to have $|C| = q^2 +q +1$. This group also stabilizes the planes $\pi_{1}$ and
$\pi_{2}$.
\begin{theorem}
If $q \not\equiv 1 \bmod{3}$, the group $C$ acts semi-regularly on the points of $Q^{+}(5,q)$.
\end{theorem}
\begin{proof}
Notice that $g^{i}\langle (x,y) \rangle = \langle (x,y) \rangle$ implies that $\mu^{i} \in F$.
But $(q^2 +q +1, q-1) = (q-1, 3) = 1$, since $q \not\equiv 1 \bmod{3}$.  Thus we have that $\mu^{i} = 1$, and so
$g^{i}$ is the identity map.
\end{proof}
This model of the Klein quadric and this cyclic group are very nice to work with algebraically; they were used
in~\cite{PW} to construct cyclic parallelisms.  For us, since each orbit has size $q^2 +q +1$, taking $x$ orbits
gives the right amount of points to have an $x$-tight set, so we will seek to combine these orbits in an appropriate
way so that the proper conditions are satisfied.

\section{New line classes}\label{sec:2}
The main tool we will use to find tight sets of $Q^{+}(5,q)$ is an eigenvector condition given in~\cite{BKLP}.
\begin{theorem}
Let $A$ be the collinearity matrix of $Q^{+}(5,q)$ and let $\mathcal{L}$ be an $x$-tight set with characteristic
vector $\chi$.  Then
\begin{equation}
\left(\chi - \frac{x}{q^2 + 1}\mathbf{j}\right)
\end{equation}
is an eigenvector for $A$ with eigenvalue $q^{2} - 1$, where $\mathbf{j}$ is the vector consisting of all ones.
\end{theorem}
This theorem follows directly from our characterization of $i$-tight sets.

We construct new examples of Cameron-Liebler line classes using the computational software MAGMA~\cite{MAG}.  We
begin by taking a distinguished plane $\pi$ in $\pg(3,q)$ along with a distinguished point $P \not\in \pi$.  Thus
$\mbox{line}(\pi) \cup \mbox{star}(P)$ is a Cameron-Liebler line class of parameter $2$, and the complement of this
set is a Cameron-Liebler line class with parameter $q^2 -1$.  Working in the Klein quadric $Q^{+}(5,q)$ as detailed
in the previous section, we can assume that $\mbox{line}(\pi)$ Klein-corresponds to
$\pi_{1} = \{(x,0) : \ x \in E^{*}  \}$ and
$\mbox{star}(P)$ corresponds to $\pi_{2} = \{(0,y) : \ y \in E^{*}  \}$.  Defining our group $C$ as before, and
requiring that $q \not\equiv 1 \bmod{3}$, we have $\pi_{1}$ and $\pi_{2}$ as orbits of $C$, as well as $q^2 -1$
other orbits each having size $q^2 +q +1$.  Our wish is to split these other orbits in half in such a way that we
obtain two Cameron-Liebler line classes $\mathcal{L}_{1}$ and $\mathcal{L}_{2}$ each having parameter
$\frac{1}{2}(q^2 -1)$.

The difficulty in searching for these line classes is twofold.  First, constructing the collinearity matrix $A$ of
$Q^{+}(5,q)$ is very time consuming.  Second, combining orbits through brute force is not at all computationally
efficient, and searching the eigenspace of $A$ corresponding to $(q^2 -1)$ for eigenvectors of the appropriate form
is not much better.  To circumvent this first problem, we avoid constructing the entire collinearity matrix of
$Q^{+}(5,q)$.  Instead, we use the following result~\cite{God}:
\begin{lemma}
  Suppose $A$ can be partitioned as
  \begin{equation}
    A =
    \left[
      \begin{array}{ccc}
        A_{11}   & \cdots & A_{1k} \\
        \vdots & \ddots & \vdots \\
        A_{k1}   & \cdots & A_{kk}
      \end{array}
    \right]
  \end{equation}
with each $A_{ii}$ square, $1 \leq i \leq k$, and each $A_{ij}$ having constant row sum $b_{ij}$.  Then any
eigenvalue of the matrix $B = (b_{ij})$ is also an eigenvalue of $A$.
\end{lemma}
The nice thing about applying this lemma is that we can easily construct an eigenvector of $A$ from an eigenvector
of $B$.  By using the point orbits of a group which we assume stabilizes the Cameron-Liebler line class we wish to
construct to partition the points of $Q^{+}(5,q)$, we can construct a matrix $B$ more easily than constructing $A$,
and an appropriate eigenvector of $B$ corresponding to the eigenvalue $(q^2-1)$ will give rise to a tight point set
of $Q^{+}(5,q)$.

Let $q \equiv 1 \bmod{4}$ and $q \not\equiv 1 \bmod{3}$.  Assume we have an $\frac{1}{2}(q^2-1)$-tight set
$\mathcal{L}_{1}$ of $Q^{+}(5,q)$, disjoint from $\pi_{1}$ and $\pi_{2}$, which is stabilized by $C$ as well as the
maps $\sigma : (x,y) \mapsto (x^q, y^q)$ and $\theta: (x,y) \mapsto (x, \omega^4 y)$, where
$\langle \omega \rangle = F^{*}$.  We should notice that
$\mathcal{L}_{2} = Q^{+}(5,q) \setminus (\pi_{1} \cup \pi_{2} \cup \mathcal{L}_{1})$ is also a
$\frac{(q^2-1)}{2}$-tight set.  By using the orbits of the group
$G = C\langle \sigma \rangle \langle \theta \rangle$ (which has order $\frac{3}{4}(q-1)(q^2+q+1)$) to partition the
points of $Q^{+}(5,q)$, we obtain a matrix $B$ that is $\frac{(q+1)}{3} \times \frac{(q+1)}{3}$ when
$q \equiv 2 \bmod{3}$, and $(\frac{q}{3}+1)\times(\frac{q}{3}+1)$ when $q \equiv 0 \bmod{3}$.  We have applied this
technique for all prime powers $q < 200$ satisfying $q \equiv 1 \bmod{4}$ and $q \not\equiv 1 \bmod{3}$, and in each
case we are able find eigenvectors corresponding to $\frac{(q^2-1)}{2}$-tight sets (one unique example up to
isomorphism for each $q$). In each case, $G$ ends up being the full stabilizer of the constructed set.\\

\section{Affine two-intersection sets}\label{sec:4}
A set of type $(m,n)$ in a projective or affine plane is a set $\mathcal{K}$ of points such that every line of the
plane contains either $m$ or $n$ points of $\mathcal{K}$;  we require that $m < n$, and we want both values to
occur.  For projective planes, there are many examples of these types of sets with $q$ both even and odd.  However,
the situation is quite different for affine planes.  When $q$ is even, we obtain a set of type $(0,2)$ in $\ag(2,q)$
from a hyperoval of the corresponding projective plane, and similarly a set of type $(0,n)$ from a maximal arc.
Examples of sets of type $(m,n)$ in affine planes of odd order, on the other hand, are extremely scarce.  The only
previously known examples exist in affine planes of order $9$~\cite{PR}, where we have sets of type $(3,6)$
containing either $36$ or $45$ points (these are complementary cases).  These examples were found by a computer
search.

A Cameron-Liebler line class $\mathcal{L}_{1}$ of parameter $40$ in $\pg(3,9)$ constructed in the previous section
induces a symmetric tactical decomposition on the space having four classes of points and lines. The four line
classes are $\mbox{line}(\pi)$, $\mbox{star}(P)$, $\mathcal{L}_{1}$, $\mathcal{L}_{2}$, where
$\mathcal{L}_{2} = Q^{+}(5,q) \setminus (\pi_{1} \cup \pi_{2} \cup \mathcal{L}_{1})$, and the four point classes are
$\pi$, $P$, and two others, $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$.  $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ are
obtained as follows: each point in the space, excluding $P$ and the points of $\pi$, lies on either $30$ or $60$
lines of $\mathcal{L}_{1}$.  We define $\mathcal{P}_{1}$ to be the set of points on $30$ lines of $\mathcal{L}_{1}$,
and $\mathcal{P}_{2}$ to be the set of points on $60$ lines of $\mathcal{L}_{1}$.  Now, if we take a plane
$\pi^{\prime}$ of $\pg(3,9)$ not equal to $\pi$, and not containing $P$, $\pi^{\prime}$ contains precisely one line
of $\mbox{line}(\pi)$, no lines of $\mbox{star}(P)$, and either $30$ or $60$ lines of $\mathcal{L}_{1}$, so $60$ or
$30$ lines of $\mathcal{L}_{2}$.  WLOG we may assume that $\pi^{\prime}$ contains $30$ lines of $\mathcal{L}_{1}$
and $60$ lines of $\mathcal{L}_{2}$.  $\pi^{\prime}$ also contains $(q+1)$ points of $\pi$, and, under our
assumptions, $30$ points of $\mathcal{P}_{1}$ and $60$ points of $\mathcal{P}_{2}$.  In fact, this is a symmetric
tactical decomposition of $\pi^{\prime}$ having $3$ classes on points and lines.  Finally, by taking
$\pi^{\prime} \cap \pi$ to be the line at $\infty$, we derive the affine plane $\ag(2,9)$.  It can be easily verified
that $\pi^{\prime} \cap \mathcal{P}_{1}$ is a set of size $30$ in $\ag(2,9)$ that is of type $(3,6)$.  As the sets of
type $(m,n)$ in $\ag(2,9)$ were completely classified in~\cite{PR}, this set is not new.  However, following this
same procedure with our Cameron-Liebler line class of parameter $3280$ in $\pg(3,81)$ yields an affine set of type
$(36, 45)$ in $\ag(2,81)$, which is new.  Our hope is that, if related Cameron-Liebler line classes exist in
$\pg(3, 3^{2e})$ for all $e$, they will always give rise to sets of type $(m,n)$ in $\ag(2, 3^{2e})$.

\bibliography{CLlineclasses}

\end{document}

config: arxiv


TITLE: How to tackle this squaring of inequality problem
QUESTION [3 upvotes]: If the roots of quadratic equation     $$x^2 − 2ax + a^2 + a – 3 = 0$$
  are real and less than $3$, find the range of $a$.

The roots are   $a \pm \sqrt {3 – a}$
For the roots to be real, we must have  a < 3.
Also, for the roots to be less than 3, we must have $\pm \sqrt {3 – a } \lt 3 – a $
If squaring both sides is allowable, I will get $(a – 2)(a – 3) > 0$. Then the problem is solved.
The question is:- how to convince others that the squaring of both sides of $\pm \sqrt {3 - a } \lt 3 - a $ is allowable?

REPLY [1 votes]: Since $x_1,x_2<3$ we have $3-x_i>0$ so their product is positive:$$ 0<9-3(x_1+x_2)+x_1x_2$$
thus $$0<9-6a+a^2+a-3 = a^2-5a+6=(a-3)(a-2)$$
So $a\in (-\infty ,2)\cup (3,\infty)$. But since the discriminat must be $\geq 0$ we get $a\leq 3$ so we have $a\in (-\infty ,2)$
Now let us prove that all $a<2$ are good.
We have $$x_1= a-\sqrt{3-a}\leq a<3$$ We are left if $x_2= a+\sqrt{3-a}$ is smaller than $3$ if $a<2$ i.e. $$\sqrt{3-a} <3-a$$ which is true since $3-a>1$


TITLE: Are `subclass` and `instance` commands equivalent in defining subclass relation between typeclasses in Isabelle?
QUESTION [8 upvotes]: Given two typeclasses (e.g. group and monoid), I have seen both the subclass and the instance commands used to establish the fact that any group is also a monoid (i.e. the subclass relation), as shown in the example below.
The two commands seem to generate similar proof requirements, but also have some differences so that the proof for one does not work for the other sometimes. For example (Isabelle 2021):
theory test
  imports Main HOL.Set
begin

declare [[show_types]] [[show_sorts]]

instance group ⊆ monoid
proof(intro_classes)
  fix x::"'a::group"
  have "÷ x ⊕ (x ⊕ ) = ÷x ⊕ x" by (simp add: invl neutl assoc_left)
  from this show "x ⊕  = x" by (simp add: left_cancel)
qed

subclass (in group) monoid
proof
  fix x::"'a::group"  (* problem line*)
qed

end

The instance command generates the following proof goal:
proof (state)
goal (1 subgoal):
 1. ⋀x::'a. x ⊕  = x
type variables:
  'a :: group

, while the subclass command generates a similar goal but without the type constraint 'a::group:
proof (state)
goal (1 subgoal):
 1. ⋀x::'a. x ⊕  = x
variables:
  neutral :: 'a
  plus :: 'a ⇒ 'a ⇒ 'a
type variables:
  'a :: type

So the first command fix x::"'a::group"  that worked in the instance case no longer works for the subclass case.
To me, the proof goal of the instance command seems more natural, as one starts from a group and just need to prove that it also respects the laws of monoid.
My questions are:
Are the two usages with instance/subclass supposed to be equivalent or do they express some fundamentally different facts/ideas?
If yes, when is one preferred over the other?
If no, how are they different in terms of intended usage.
(I asked the same question in SO a couple of months ago, but have not received any response. Posting here to hopefully get more attention.)

REPLY [2 votes]: Citing The Isabelle/Isar Reference which is part of the official documentation:

A weakened form of [subclass] is available through a further variant of
instance: instance c1 ⊆ c2 opens a proof that class c2 implies c1
without reference to the underlying locales; this is useful if the properties to prove the logical connection are not sufficient on the locale level
but on the theory level.

We can see this with the same example as yours in Tutorial.Axioms, on line 149 we have a theory level lemma:
theorem left_cancel:
 ((?x::?'a::group) ⊕ (?y::?'a::group) = ?x ⊕ (?z::?'a::group)) = (?y = ?z)

And on line 166 (instance group ⊆ monoid) can fix x :: 'a group as you noticed, but we also get access to left_cancel on line 172.
For subclass (in group) monoid we only get x :: 'a :: type and no access to theory level theorems (eg. left_cancel). It is still provable and will register a sublocale group < monoid which instance won't.

Now, to answer your two questions:

These usages of instance and subclass are not supposed to be equivalent. For example, as you have noticed yourself, the former carries the theory level sort 'a :: group.

If you the locale level is not sufficient to prove the connection, you can go with the former. The latter is useful if you want to access theorems from c2 inside the c1 context.

More references: Haskell-style type classes with Isabelle/Isar and Tutorial to Locales and Locale Interpretation


TITLE: Square root principle value convention
QUESTION [0 upvotes]: Why is the principal square root of a complex number defined as
$\sqrt z = \sqrt r e^{-i \varphi / 2}$
for $\varphi \in (-\pi, \pi]$ ?
Wouldn't it be more natural to let $\varphi \in [0, 2\pi)$ as it is usual for polar coordinates?
If you think that you are free to choose whatever variant you like, then it's not quite the case. Software libraries definitely prefer the first one so if you are to take some assistance from a computer you are forced to use $\varphi \in (-\pi, \pi]$.

REPLY [0 votes]: Not to leave the question unanswered, per @copper.hat 's comment.
The choice is driven by desirability to have continuous square root in a neighborhood of positive real numbers.


TITLE: How do we know quantum entanglement really happens?
QUESTION [1 upvotes]: Quantum entanglement says that we can entangle two particles with an opposite property, separate them, and that they exist in an indeterminate state until we observe one of them.
eg. If we create a left/right pair, separate them in space then observe one of the particles then the other particle of the pair will instantly resolve into the opposite state of the particle we observed.
Question: How do we know this is true? How do we know that something happened to the second particle as a result of observing the first?
A much simpler explanation is that both the particles flipped into fixed left/right states at the moment we started to separate them, we simply hadn't opened the containers yet to see which was which.
I assume that there can be no possible experiment to distinguish between these two explanations since we're not allowed to observe either particle.
What reason(s) do we have to believe the magical 'quantum' explanation when the simple explanation explains it just as well?

REPLY [1 votes]: It's a very good question, one that John Bell asked. In 1964, he proposed an experiment to distinguish between the two cases, one case where a hidden variable is carried along, and another case where the measurement result of one somehow propagates to the other. Bell derived "Bell's inequality" which would be violated if the second case turned out to be correct.
In 1982, Alain Aspect performed the first experiment to test the inequality, and the result confirmed entanglement, not hidden variables. Since then, many other experiments have also confirmed entanglement.
Like you, and most people who first come across this, I used to feel that "a simpler explanation is that both the particles flipped into fixed left/right states at the moment we started to separate them." So I went through the exercise of trying to come up with a function that maps the hidden variables to up/down (or left/right) in a way that violates the inequality like the entanglement experiments do. After enough playing around, I got a feel for how hidden variable cannot do this. I recommend this exercise for those who want to get better feel for how entanglement works.


TITLE: Cardinality of $\mathbb{R}^\mathbb{R}$ by diagonal argument
QUESTION [2 upvotes]: I want to prove that the set of all real functions $\mathbb{R}^\mathbb{R}$ has a higher cardinality than the real numbers $\mathbb R$, by Cantor's diagonal argument.
I'm having difficulties with approaching this problem.
What I'm looking for is a hint in the right direction.
I've seen an example where it is shown that the power set $2^S$ of a countable set $S$ is uncountable, by the diagonal argument.

REPLY [0 votes]: A rough outline of the proof.
Consider a subset $U$ of $\mathbb R ^\mathbb R$. It can be defined by a sequence of characteristic functions $s_i$.
Now, consider a subset $T$ of $\mathbb R$. It can be defined by a single characteristic  function $\chi_j$.
It is clear that a one-to-one correspondence between $U$ and $T$ is not possible; because a map $\chi_j \mapsto s_i$ can never be surjective.
Hence, $\mathbb R ^\mathbb R$ has a higher cardinality than $\mathbb R$.


TITLE: Rice’s theorem explanation
QUESTION [3 upvotes]: I've recently read a article 'Theories of computational complexity'. In the article it states that with Rice's theorem one can prove that, for example, {$i∈[\mathbb N:\phi_i$ is a constant function}, is not recursive, but I can't see how. I already know the theorem, but I don't really get how I can use it. Thank you for your answer!

REPLY [5 votes]: Rice's theorem says that for any subset $F$ of the class $T$ of partial computable functions, the set $\{i\mid \phi_i\in F\}$ is recursive iff $F=\emptyset$ or $F=T$.
Let $F$ be the set of partial computable functions that are constant on their domain. $F$ is certainly non-empty, since the function that always takes the value 0 is certainly computable. $F$ is also not the whole class $T$, since the function $\phi$ with $\phi(0)=0$, $\phi(n)=1$ for $n>0$ is certainly computable and not constant.
It follows that the set you are interested in is not recursive.
If your confusion lies in how to translate the usual statement of Rice's theorem (in terms of decision problems) into the version I wrote, simply remember that formally a decision problem is a set $D\subseteq{\mathbb N}$, and the problem is decidable iff $D$ is recursive. Of course, this is usually stated by saying that there is a Turing machine that, given $n$, decides whether $n\in D$ or not. But this is certainly the same as saying that the machine computes the characteristic function of $D$, which is (the definition of) being recursive.


TITLE: {$f_n(x)$} is defined by $f_{n+1}(x)=\int_{0}^{x}f_n(t)dt$, Prove: $\sum_{1}^{\infty}f_n(x)$ converges and is differentiable in $(0,1)$.
QUESTION [4 upvotes]: Given $f_1(x)$ continuous function in $[0,1]$ and differentiable in $(0,1)$. The series {$f_n(x)$} is defined by $f_{n+1}(x)=\int_{0}^{x}f_n(t)dt$ for $x\in [0,1]$.
I need to prove that $\sum_{1}^{\infty}f_n(x)$ converges and is differentiable in $(0,1)$.
I tried to bound $f_n(x)$ by some other series and to use Weierstrass M-Test to prove uniform convergence which will lead to regular one, but I didn't manage to find a series.
I'd really love your help with this one.
Thanks a lot!

REPLY [8 votes]: Since $f_1$ is continuous on $[0,1]$, there exists $M$ such that $|f_1(x)|\leq M$ for all $x\in [0,1]$. Then
$$|f_{2}(x)|=\left|\int_{0}^{x}f_1(t)dt\right|\leq\int_{0}^{x}|f_1(t)|dt\leq Mx.$$
Similarly,
$$|f_{3}(x)|=\left|\int_{0}^{x}f_2(t)dt\right|\leq\int_{0}^{x}|f_2(t)|dt\leq Mx^2/2,$$
and by induction,
$$|f_{n+1}(x)|=\left|\int_{0}^{x}f_n(t)dt\right|\leq\int_{0}^{x}|f_n(t)|dt\leq Mx^n/n!.$$
Therefore, we have
$$\left|\sum_{n=1}^{\infty}f_n(x)\right|\leq\sum_{n=1}^{\infty}|f_n(x)|\leq\sum_{n=1}^{\infty}Mx^n/n!$$
which converges for all $x\in(0,1)$. By Weierstrass_M-test, $\displaystyle \sum_{n=1}^{\infty}f_n(x)$ converges on $(0,1)$.
On the other hand, Since $f_1$ is differentiable on $(0,1)$ and $f_{n+1}(x)=\int_{0}^{x}f_n(t)dt$ for $n\geq 2$, by induction and the fundamental theorem of calculus, $f_n$ is differentiable on $(0,1)$ for all $n$. This implies that $\displaystyle \sum_{n=1}^{\infty}f_n(x)$ is differentiable on $(0,1)$.


TITLE: Do all singular $n\times n$ matrices form a vector subspace when $n\ge2$?
QUESTION [2 upvotes]: I want to prove or disprove that the set of all $n\times n$ singular matrices form a vector subspace of $M_{nn}$ when $n\geq 2$. So, let:
$$
A_{n,n} =
\begin{pmatrix}
  a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\
  a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\
  \vdots  & \vdots  & \ddots & \vdots  \\
  a_{n,1} & a_{n,2} & \cdots & a_{n,n}
 \end{pmatrix}
 ,\text{ where }a_{i,j}\in \mathbb R$$
I can't just take the case where $A_{n,n} = $ the zero vector and show the space is non-empty, closed under addition and scalar multiplication, can I?
this is only one case of $\det(A) =0$
so how can I prove/disprove this for all cases where $\det(A) =0$
Thanks in advance.

REPLY [6 votes]: Hint.  What do you get when you add
$$\pmatrix{1&0\cr0&0\cr}\quad\hbox{and}\quad \pmatrix{0&0\cr0&1\cr}\ ?$$
Note that

to prove a set is closed under addition you must show that if $M,N$ are in the set then $M+N$ is always in the set;
to show a set is not closed under addition, you must find one example of $M,N$ in the set, such that $M+N$ is not in the set.


TITLE: Simulation of PDF of inverse random varible
QUESTION [0 upvotes]: Let $X$ be exponential random variable with parameter $\beta$
$$
P_X(x)=\frac{1}{\beta}e^{-x/\beta}.
$$
I have found a lemma saying

Let $X$ be a random variable with pdf $p_X (x)$ for all $x\geq 0$ and $p_X (x) = 0$
  for $x < 0$. Then, the pdf of $Y = 1/ X$ is
  $$
P_Y(y)=\frac{1}{y^2}P_X(\frac{1}{y}).
$$

My question is, How we use MATLAB with Monte Carlo simulation to generate or simulate $X$ and compared with theoretical PDF $P_X(x)$ and also simulate $P_Y(y)$ and compared with theoretical expression.

REPLY [1 votes]: The classical way to simulate an $Exp[\beta]$ random variable is by considering random variable
$$X'=-\beta \log(U)$$
where $U$ is uniformly distributed on $[0,1] \ $ ($\log$ is the natural logarithm). It is easy to prove that $X'$ has a  $Exp[\beta]$ distribution.
Thus the simulation program in Matlab is :
 clear all;close all;
 n=100000;
 beta=10;
 E=-beta*log(rand(1,n));
 hist(E,30)

Second part : What can be said about $Z$ defined by
$$\dfrac{1}{Z}=\dfrac{1}{X}+\dfrac{1}{Y}$$
with $X,Y$ following a $Exp[\beta]$ distribution ?
The problem with the inverse of an $Exp[\beta]$ distribution (like $\dfrac{1}{X}$ or $\dfrac{1}{Y}$) is that it has an infinite mean, thus a non-defined variance (https://stats.stackexchange.com/questions/229543/mean-of-inverse-exponential-distribution)
Thus, instead of simulating $1/X$; $1/Y$, adding them and taking the inverse, I have directly done this simulation, for the particular case $\beta=1$ :
 clear all;close all;
 n=1000000;
 U=-log(rand(1,n));V=-log(rand(1,n));
 Z=(U.*V)./(U+V);
 hist(Z,200)
 [mean(Z),var(Z)]

The results I obtain tend to give $E(Z)=\dfrac{1}{3}$ and $var(X)=\dfrac{4}{45}$.
Here is the histogram generated by this program. One notices that after a certain small peak, there is rapid decline.
I have no idea about the way to get these results throretically.


TITLE: Alternative approaches to probability theory
QUESTION [6 upvotes]: I'm undergraduate student in probability theory (and its applications). There are lots of different and definitely good text on standard, functional analysis-based approach, but I'm interested in alternative approaches - maybe, some more algebraic variants. Could you name some ideas/papers/texts about this? I'm especially interested in the ones those can be used in applied problems (such as financial mathematics or something). I've surely saw some approaches in Wikipedia, but I have absolutely no idea of using them in practical problems.

REPLY [1 votes]: Peter Whittle's "Probability via expectation" is a very nice book. In taking expectation (rather than probability measure) as a starting point, the text manages to

Avoid measure theory almost entirely
Get to interesting applied problems very quickly

Though Whittle (probably deliberately) never says this, the approach here is essentially the `Daniell integral' approach to integration theory, cast in a probabilistic context. The text is (in my opinion) spectacularly well written, both in terms of the choice/order of material, and in the quality and warmth of the written English. I cannot think of many books of this type (i.e. one advocating an 'alternative' approach to the development an established theory) which cover so much material quite so masterfully whilst remaining so focused and true to their founding premise.


text:
\section{Open Unit Interval is Proper Subset of Closed Unit Interval}
Tags: Real Intervals

\begin{theorem}
The [[Definition:Open Unit Interval|open unit interval]]:
:$I_o = \left({0 \,.\,.\, 1}\right)$
is a [[Definition:Proper Subset|proper subset]] of the [[Definition:Closed Unit Interval|closed unit interval]]:
:$I_c = \left[{0 \,.\,.\, 1}\right]$
\end{theorem}

\begin{proof}
Let $x \in I_o$.
Then by definition:
:$0 < x < 1$
and so:
:$0 \le x \le 1$
and so:
:$x \in I_c$.
Thus:
:$I_o \subseteq I_c$
Consider:
:$0 \in I_c$
by definition of [[Definition:Closed Real Interval|closed interval]].
But it is not the case that $0 < 0$.
So $0 \notin I_o$ and so $I_c \not \subseteq I_o$.
Hence the result by definition of [[Definition:Proper Subset|proper subset]].
{{qed}}
\end{proof}

config: wiki


TITLE: Not every submanifold is an inverse image of a regular value
QUESTION [5 upvotes]: $\Bbb RP^1 \to \Bbb RP^2$ is an embedding. Show that $\Bbb RP^1$ can't be denoted as an inverse image of a regular value for some smooth map on $\Bbb RP^2$.

REPLY [5 votes]: Suppose $f$ is any smooth function from $\mathbb{R}P^2$ to $\mathbb{R}$ which has $\mathbb{R}P^1$ as a level set.  We will show that every point of $\mathbb{R}P^1$ is singular.
By changing $f$ to $f + c$ for a constant $c$, we may assume wlog that $\mathbb{R}P^1 = f^{-1}(0)$.
Since $\mathbb{R}P^2 - \mathbb{R}P^1$ is connected and since $f$ is never $0$ on it, $f$ must have constant sign on $\mathbb{R}P^2 - \mathbb{R}P^1$.  By replacing $f$ by $-f$ if necessary, we may assume the sign is positive.
But then $0$ is the global minimum of $f$ and so by standard calculus, the derivative of $f$ must vanish at all points in $\mathbb{R}P^1$.  In particular, every point in $\mathbb{R}P^1$ is singular.


TITLE: Poincaré duality and Hodge duality
QUESTION [5 upvotes]: The Poincaré duality is defined in Greub's Multilinear algebra (1967) in Chapter 6, §2 as a isomorphism between $\bigwedge V$ and $\bigwedge V^*$, where $V$ is a finite-dimensional vector space, $V^*$ is its dual, and $\bigwedge V$ is the exterior algebra of $V$. More precisely, the isomorphism maps $\bigwedge^k V$ to $\bigwedge^{n-k} V^*$, $n=\dim V$. Greub's definition has nothing to do with cohomology groups of manifolds and is based entirely on multilinear algebra. I have tried reading Chapter 6 and quickly realised I would need to study the entire book to understand the definition. I am only interested in the cases $V=\mathbb R^n$ where $n\le5$, and I only need to see what the isomorphism looks like on the standard basis of $\bigwedge \mathbb R^n$ and $\bigwedge (\mathbb R^n)^*$.
For example, let $e_1,e_2,e_3$ denote the standard basis of $\mathbb R^3$ and $e_1^*,e_2^*,e_3^*$ the standard basis of $({\mathbb R}^3)^*$. Then the basis of $\bigwedge\mathbb R^3$ consists of 1, $e_1,e_2,e_3$, $e_1\wedge e_2$, $e_2\wedge e_3$, $e_3\wedge e_1$, $e_1\wedge e_2\wedge e_3$ and the basis of $\bigwedge(\mathbb R^3)^*$ consists of 1, $e_1^*,e_2^*,e_3^*$, $e_1^*\wedge e_2^*$, $e_2^*\wedge e_3^*$, $e_3^*\wedge e_1^*$, $e_1^*\wedge e_2^*\wedge e_3^*$. The isomorphism maps, say, $e_1^*$ to $e_2\wedge e_3$ or perhaps $-e_2\wedge e_3$. I am not sure which sign is the right one. That is, I already know what the isomorphism looks like but only up to a sign. So, my question is what the correct signs are in the mapping Greub uses.
My second question is how Poincaré duality as defined by Greub is related to Hodge duality. The latter depends on a bilinar form on $(\mathbb R^n)^*$ but I am only interested in the form defined by $\langle e_i^*,e_i^*\rangle=1$ and $\langle e_i^*,e_j^*\rangle=0$ if $i\ne j$. The Hodge duality maps $\bigwedge^k(\mathbb R^n)^*$ to $\bigwedge^{n-k}(\mathbb R^n)^*$, so it is restricted to the exterior algebra of $(\mathbb R^n)^*$, but its effect is somewhat similar to that of the Poincare duality. For instance, the Hodge duality maps $e_1^*$ to $e_2^*\wedge e_3^*$. If I define a map $\mathcal{R}:\bigwedge (\mathbb R^n)^*\to\bigwedge \mathbb R^n$ in the obvious way by removing $*$ from the relevant basis vectors, e.g. $\mathcal{R}(e_1^*\wedge e_2^*)=e_1\wedge e_2$, would the Poincaré duality $\mathcal{P}:\bigwedge (\mathbb R^n)^*\to\bigwedge \mathbb R^n$ and the Hodge duality $\mathcal{H}$ be related by $\mathcal{P}(a)=\mathcal{R}(\mathcal{H}(a))$ where $a\in \bigwedge (\mathbb R^n)^*$.
EDIT:
Despite what Bruno Joyal says, I think $\mathcal{P}(a)=\mathcal{R}(\mathcal{H}(a))$ is valid without any sign adjustment, provided that $\langle,\rangle$ is defined as above. I define the Hodge dual (or Hodge star *, $\mathcal{H}(a)=*a$) by
$a\wedge(*b)=\langle a,b\rangle \omega$,
where $\omega=e_1^*\wedge\dots\wedge e_n^*$ is the standard element of $\bigwedge^n(\mathbb R^n)^*$, and $\langle a,b\rangle=\det(\langle a_i,b_i\rangle)$ for $a=a_1^*\wedge\dots\wedge a_k^*$ and $b=b_1^*\wedge\dots\wedge b_k^*$, which can be replaced by 1 in calculations with the basis vectors thanks to my choice of the inner product.
I don't understand the definition of Poincaré duality but Greub does give the formula for computing the duals,
$\mathcal{P}(e_{\nu_1}^*\wedge\dots\wedge e_{\nu_p}^*)=(-1)^{\sum_{i=1}^p(\nu_i-i)}e_{\nu_{p+1}}\wedge\dots\wedge e_{\nu_n}$,
where $\nu_1<\nu_2<\dots < \nu_n$, and $(\nu_{p+1},\dotsc,\nu_n)$ is complementary to $(\nu_1,\dots,\nu_p)$. As far as I can see, the factor $(-1)^{\sum_{i=1}^p(\nu_i-i)}$ is the sign of the permutation that takes $(\nu_1,\dots,\nu_n)$ to $(1,\dots,n)$, but $a$ and $*a$ are related in the same way for the basis elements of $\bigwedge(\mathbb R^n)^*$. I just need to use the map $\mathcal{R}$ to go to the primary space. Hence, $\mathcal{P}(a)=\mathcal{R}(\mathcal{H}(a))$. I confirmed that by calculations for $n=2,3,4$. The fact that $*a$ is not an involution doesn't really matter.

REPLY [0 votes]: $
\newcommand\R{\mathbb R}
\newcommand\Ext\bigwedge
\newcommand\form[1]{\langle#1\rangle}
\newcommand\rev\widetilde
\newcommand\hodge\star
$
The punchline is that
$$
  \hodge X = D_\#(Q(X))\qquad \forall X\in\Ext\R^n,
$$
where $\hodge$ is the Hodge star, $D_\# : \Ext(\R^n)^* \to \Ext\R^n$ is one of Greub's Poincaré isomorphisms, and $Q : \Ext\R^n \to \Ext(\R^n)^*$ is the isomorphism induced by the non-degenerate bilinear form used to define $\hodge$.

Choosing a non-degenerate bilinear form on $\R^n$ amounts to a choice of isomorphism $Q : \R^n \to (\R^n)^*$, and upon making this choice $\Ext\R^n$ inherits the Clifford product of the associtated Clifford algebra. By the universal property of Clifford algebras $Q$ extends to a Clifford algebra isomorphism $\Ext\R^n \to \Ext(\R^n)^*$ where we give $\Ext(\R^n)^*$ the dual Clifford product defined by $Q^{-1}$. Writing Clifford multiplication as juxtaposition, for $X, Y \in \Ext\R^n$ we get
$$
  \form{Q(X), Y} = \form{\rev XY}_0,
$$
where $\rev X$ is the reversal (i.e. main anti-automorphism) and $\form{\cdot}_0$ is scalar projection.
One of the Poincaré isomorphisms defined by Greub is
$$
  D_\# : \Ext(\R^n)^* \to \Ext\R^n,\quad D_\#(\xi) = \iota_\xi(I)
$$
where $\iota$ is the interior product and for some choice of $I \in \Ext^n\R^n$, which with our Clifford product we will choose to be a unit $I^2 = \pm1$ to line up with the Hodge star later. For $A \in \Ext^{n-k}\R^n$ and $B \in \Ext^{k}\R^n$, we see that the manifestation of $D_\#$ in the Clifford algebra is
$$\begin{aligned}
  \form{\rev A(D_\#\circ Q(B))}_0
  &= \form{Q(A), D_\#(Q(B))}
\\
  &= \form{Q(A), \iota_{Q(B)}(I)}
\\
  &= \form{Q(B)\wedge Q(A), I}
\\
  &= \form{Q(B\wedge A), I}
\\
  &= \form{\rev{B\wedge A}\:I}_0
\\
  &= \form{\rev A\wedge\rev B\:I}_0
\\
  &= \form{\rev A(\rev BI)}_0
\end{aligned}$$$$
  \implies D_\#\circ Q(B) = \rev BI,
$$
with the last line following since $A$ is arbitrary and $\form{\rev{(\cdot)}\,{\cdot}}_0$ is a non-degenerate bilinear form. The step $Q(B)\wedge Q(A) = Q(B\wedge A)$ is valid since $\wedge$ can be realized with the Clifford product, e.g. $A\wedge B = \form{AB}_n$ with $\form{\cdot}_n$ the projection onto $\Ext^n\R^n$. Since we proved the above for simple $k$-vectors, it holds for all multivectors:
$$
  D_\#\circ Q(X) = \rev XI\qquad\forall X\in\Ext\R^n.
$$
Now consider the Hodge star, which following this is defined by
$$
  A\wedge(\hodge B) = \form{\rev AB}_0I
$$
for $A, B \in \Ext^k\R^n$. Converting the RHS to a $\wedge$ and multiplying by $I$ on the left, we immediately get
$$
  I\:A\wedge(\hodge B) = I\:A\wedge(\rev BI),
$$
where we've also move the reversal onto $B$ since both $A$ and $B$ are $k$-vectors. Then
$$
  \form{(IA)(\hodge B)}_0 = \form{(IA)(\rev BI)}_0,
$$
and it follows by non-degeneracy of the form that
$$
  \hodge B = \rev BI
$$
and this result extends to all multivectors. We have finally proven that
$$
  \hodge X = D_\#(Q(X))\qquad \forall X\in\Ext\R^n.
$$


TITLE: Question about quaternionic conjugation
QUESTION [1 upvotes]: The quaternionic conjugation is defined by
$$\begin{aligned}i &\mapsto -i\\j&\mapsto -j\\k&\mapsto -k\end{aligned}$$
But since $ij=k$, shouldn't we have that $k = ij \mapsto (-i)(-j) = ij = k$?

REPLY [3 votes]: That's because quaternion conjugation, unlike complex conjugation, is not an automorphism. In fact, for quaternions $a$, $b$ we have
$$
(a b)^* = b^* a^*
$$
so in you example
$$
(ij)^* = j^* i^* = (-j)(-i)=ji=-ij=-k=k^*
$$


TITLE: Representation of direct product $G\times H$ on tensor product
QUESTION [2 upvotes]: Let $G$ and $H$ by finite groups.
Let $V$ and $U$ be irreducible representations of $G$ and $H$, respectively.
When the ground field is $\mathbb{C}$, I know how to show that $V\otimes U$ is an irreducible representation of $G\times H$ (using character theory).
What about other ground fields? Is $V\otimes U$ guaranteed to be an irreducible rep of $G\times H$ if the ground field is $F$ and $\operatorname{char}F\nmid|G|$?

REPLY [3 votes]: No. For example, let $G$ be cyclic of order $3$. Then $G$ has an ireducible module $U$ of dimension $2$ over ${\mathbb Q}$. But $U \otimes U$ is not irreducible for $G \times G$; it decomposes into two $2$-dimensional modules.