\section{Motivation and introduction}
\label{sec:Mot}
The accuracy and precision of a measurement are reduced by drifts in the  offset and sensitivity of the instrument used to collect the data. The offset or zero-point  is additive  and  independent of the magnitude of the value being measured. If uncorrected,  it directly biases the result of the measurement. Sensitivity on the other hand is multiplicative and  affects the output of the instrument more for larger quantities than for smaller ones. Sensitivity can be measured and any variation corrected by repeated or even continuous calibration of the instrument. The offset can be eliminated  using modulation, i.e. switching the measurement setup such that the quantity being measured alternates in sign with  measurement states $A$ and $B$.  As long as  the zero-point does not change  over the time it takes to make the measurement,  the result determined from the difference of these states   is independent of the offset. As an example, a  technique  employed by many laboratories  uses a single pan balance  for making mass measurements. The sample $A$ and a known calibrated mass  $B$ are alternately placed on the balance's pan and their values read out. The mass difference, obtained by subtracting two subsequent readings, and the known mass $B$ give the mass of $A$. An offset that is constant in time does not bias the mass difference. The result can be made less dependent on the sensitivity of the balance by choosing a calibration weight whose mass is nearly identical to that of sample. Alternately, this can be viewed as modulation where the measurement result is half the difference of the measured masses and the zero-point is their average.

If the time required to make an $AB$ comparison is sufficiently long compared to the drift rate, the assumption of an unchanging zero-point  is no longer  valid. In this case the result includes a  systematic bias from the zero-point drift. This problem is more acute for measurements taken over very long times due to the $1/f$ behavior of noise. The term $1/f$ noise reflects the fact that the power spectral density of the noise  is inversely proportional to the frequency~\cite{Milotti02}. Measurements made over long time intervals would experience larger fluctuations in the zero-point than those over short time intervals. As a consequence, an ensemble of mean values taken at different times has more scatter than the errors of the means themselves. The presence of $1/f$ noise has been shown to be fully consistent with a stationary distribution~\cite{stoisiek76} and this extra noise can be added to the standard deviation of the distribution.

The technique described in the above example is not limited to mass metrology. Many precision experiments change one critical control parameter or a part of the geometry to achieve modulation in the experiment's output. Examples include measurements of the gravitational constant~\cite{Sch06}, the search for parity violation in nucleon-nucleon interactions~\cite{Swa09,Bass09}, the determination of Planck's constant using a Watt balance~\cite{Ste05} and searches for equivalence-principle violations~\cite{Sch08}. As the examples indicate square wave modulation is a powerful tool used for Null and non-Null experiments. The lock-in amplifier is able to recover noisy signals using similar techniques.

Consider a series of measurements,${u_i}$, spaced equally in time  where an observable $\eta$ changes sign with each measurement, $z_0(i)$ denotes the zero-point that varies slowly with $i$, and $\varepsilon_i$ is a random variable that adds noise to the measurement. The expectation value of $\varepsilon$ is $E(\varepsilon)=0$ and the noise is assumed to be uncorrelated and stationary, $\Cov{(\varepsilon_i,\varepsilon_j)}=\delta_{ij} \sigma^2$.  The time series can be written as
\begin{equation}
u_i = z_0(i) + (-1)^i \eta + \varepsilon_i,
\label{eqn:def1}
\end{equation}
where $i$ runs from 1 to $N$, which is assumed to be even.
An estimator $S$ of the signal can be obtained by multiplying the time series by $(-1)^i$ and subsequent averaging,
\begin{equation}
S = \frac{1}{N} \sum_i ^N \left( (-1)^i z_0(i)  + \eta + (-1)^i \varepsilon_i \right).
\end{equation}
The expectation value and the variance of the estimator are
\begin{eqnarray}
E(S) &=& \eta + \frac{1}{2} E (z_0\prime ) \\
\Var(S) &=& E (z_0 \prime^2) -  (E {z_0}\prime)^2 +\sigma^2
\end{eqnarray}
The signal $\eta$ is recovered, however it is biased by the mean value of the slope of the varying background, $z_0\prime\equiv z_0(i+1)-z_0(i)$. The variance of the estimator contains the variance of the measurement noise $\sigma^2$ plus variations about the mean slope of $z_0(i)$. This later part limits the ultimate statistical precision of the measurement.

Various patterns for taking data are known to completely cancel linear and even higher order drifts in calculations of the mean~\cite{Sut93,Gla00}.  For example, the sequences $ABA$ and $ABBA$ both remove linear drift terms.
$ABBA$  with 4 measurements ${u_1,u_2,u_3,u_4}$, recovers the signal by computing the mean according to the expression  $S_4 = \frac{1}{4} (u_1-u_2 - (u_3-u_4))$. Assuming uncorrelated data ($\Cov{(u_i,u_j)}=\delta_{ij} \sigma^2$), the variance of this estimator is given by $\Var{(S_4)}=\sigma^2/4$.  Averaging $N/4$ of these $S_4$ results gives an overall variance of  $\sigma^2/N$.
In the $ABA$ sequence, the estimator is given by $S_3 = \frac{1}{4} (u_1 -2u_2 + u_3)$ with $\Var{(S_3)}=3\sigma^2/8$. Averaging $N/3$ of these $S_3$ results leads to a variance of  $(9/8)\sigma^2/N$.
In general, the smallest possible variance is achieved when all measurements are weighted with equal absolute values. When comparing different drift cancellation methods a useful figure of merit is the relative increase in the uncertainty of the result $S_k$ over the minimum variance $\sigma/\sqrt{N}$. Thus  we define the statistical cost of $S_k$ as
\begin{equation*}
c(N) \equiv  \sqrt{ \frac{N} {\sigma^2 }\Var(S_k) } -1.
\nonumber
\end{equation*}
Comparing the statistical cost values, 0.06 for $ABA$ and 0 for $ABBA$, shows the latter to be better. Strictly speaking, the data set size $N$ should be divisible by both $3$ and $4$ when directly comparing these two methods.  Although by this metric $S_4$ is the best unbiased estimator in the presence of linear drift, the $ABBA$ pattern has  several disadvantages: (1) The series is asymmetric. The first $B$ occurs after a configuration change, whereas no such change occurs for the second $B$ measurement. It is reasonable to assume the two $B$ measurements have different uncertainties and drift properties.   (2) The measurement pattern is fixed in advance. If the data exhibits additional non-linear drift, removing it would require a new measurement with a more appropriate pattern. (3) Drift cancellation only occurs if all measurements of the series are used. If one measurement is missing or affected by an out-lier, the remaining measurements in that sequence would have to be discarded.

In this article we propose an alternative approach  that avoids the above shortcomings. Data is taken in an alternating (square wave) sequence $ABAB...$ without regard for any inherent drift. An estimator $F_N^p$  is then constructed that removes drift up to any desired polynomial order $p$ where $p<N$. When applied to the data it yields an expectation value for the mean independent of the  drift orders selected. This plays the same roll as the data taking sequences discussed above where linear drift is removed.   The statistical cost associated with $F_N^p$ increases with increasing $p$ and decreases with increasing $N$.   A cost -- benefit analysis can help choose  appropriate values for $p$ and $N$.

The next two sections contain the mathematical justification for the methods used. The construction of the estimator $F_N^p$ is described in two steps. First, we introduce a filter that removes  polynomial drift from the input data. The output of this filter is a shorter time series without drift to the order  selected. Next, the filtered data are averaged to a single number which is  the best estimate of the signal. In practice, the estimator $F_N^p$ is used to generate a set of weights that are  applied to the entire data set. These products are then summed together for the result.

\section{Design of a filter that removes polynomial drift from data}
\label{sec:filt_design}

The filter is a linear combination of consecutive terms in the original data series ${u_1,u_2,..,u_N}$. A new series ${y_1,y_2,...,y_{N-p-1}}$ is calculated as
\begin{equation}
y_i = \sum_{k=0}^{p} C_k u_{i+k}
\label{eq:lin-comb}
\end{equation}
The $C_k$s are chosen to cancel variations across the included points that are less than a given polynomial order $p$ and  preserve a signal that changes sign with each successive term.
 These design requirements for the filter can be met by solving the matrix equation
\begin{eqnarray}
 \left(
\begin{array}{ccccccc}
1  	&1 	&1     	& 1 	&\cdots	 &1  \\
0 	&1  	&2   	& 3 	&\cdots	 &p  \\
0    	&1 	&4   	& 9 	&\cdots	 &p^2 \\
 	&   	& 	&\vdots	&	 &   \\
0   & 1    	& 2^{p-1} & 3^{p-1} &\cdots	 &p^{p-1} \\
1 	&-1   	&1 	&-1	&\cdots	 &(-1)^p \\
\end{array}
\right)
 \left(
\begin{array}{c}
C_0 \\ C_1 \\ C_2 \\ \vdots \\ C_{p-1} \\ C_p \\
\end{array}
\right)=
\left(
\begin{array}{c}
0 \\ 0 \\ 0 \\ \vdots \\ 0 \\1 \\
\end{array}
\right),
\label{eq:matrix}
\end{eqnarray}
for $C_k$. The first $p$ rows cancel variations with polynomial orders up to $p-1$ including the zeroth order, i.e. offset. The last row normalizes the part of $u_i$ that changes sign with each successive $i$. The $p+1$ equations uniquely determine $p+1$ coefficients of the sum in equation~(\ref{eq:lin-comb}).

An alternate analytical expression for these coefficients can be obtained from the generating function $F(z) = (1-z^{-1})^{p}$ where $z$ is complex. In the limit $|z|\rightarrow 1$ it can be seen that the derivatives of $F(z)$ less than order $p$ vanish. This has a correspondence with equation~(\ref{eq:matrix}), where in the first $p$ rows the sum of each successive polynomial order is required to vanish. The binomial expansion of $F(z)$ gives

\begin{equation}
F(z) = (1-z^{-1})^p = \sum_{k=0}^p (-1)^k \left( \begin{array}{c}
p \\ k
\end{array} \right) z^{-k} \label{eqn:zt}
\end{equation}
$F(z)$  is therefore the Z-transform of a discrete time series $f[k]$ with the property that the sum of terms $k^jf[k]$ equal zero for $j < p$. By taking the Z-transform of both sides of equation~(\ref{eq:lin-comb}) its transfer function $H(z)$ can be computed,

\begin{equation}
H(z) = \frac{Y(z)}{U(z)} = \sum_{k=0}^p C_k z^{-k}.
\end{equation}

Equating $H(z)$ to $F(z)$ and normalizing with $1/2^p$ to preserve the signal amplitude, the coefficients
\begin{equation}
C_k = \frac{(-1)^k}{2^p}
\left( \begin{array}{c}
p \\ k
\end{array} \right)
\end{equation}
are obtained.

Substituting these coefficients into equation~(\ref{eq:lin-comb}) one obtains
\begin{equation}
y_i(p) = \frac{1}{2^p}\sum_{k=0}^p  (-1)^{i+k}  \left( {p \atop k} \right) u_{i+k}.
\label{eq:lin}
\end{equation}
An additional factor of $(-1)^{i}$ is included here to demodulate the alternating signal $\eta$ found in equation~(\ref{eqn:def1}).

Equation~(\ref{eq:lin}) is a discrete convolution of the input data $u_i$ with a finite input response (FIR) filter kernel. The filter described here has the property of removing polynomial drifts. In principle, any (FIR) filter could be used, and another  may be more efficient in environments where the functional form of the drift is known.

To properly account for uncertainties it is advantageous to express the convolution as a linear transformation of the input data set $\bi{U}={u_1,..,u_N}$ to a filtered data set $\bi{Y}={y_1,..,y_{N-p}}$:
\begin{equation}
\bi{Y}=\bi{A}\bi{U}
\label{eq:trans}
\end{equation}
The matrix $\bi{A}$ has dimensions $(N-p) \times N$ as incompletely defined terms in the convolution are discarded. The elements of each row are taken from the coefficients in equation~(\ref{eq:lin}):
\begin{equation}
a_{ij} = \frac{1}{2^p} \left( {p \atop j-i} \right)(-1)^{j}  \mbox{~~~for $0\le j-i \le  p$, zero otherwise}.
\label{eq:aij}
\end{equation}
Adjacent values in the filtered data set are not statistically independent, since filtering introduces correlations between the elements of $\bi{Y}$. The covariance matrix of the filtered data can be calculated from the covariance matrix of the input data $ \Cov{(U,U)} $ using,
\begin{equation}
\bi{C} \equiv \Cov{(Y,Y)}  =\bi{A} \Cov{(U,U)} \bi{A}^T.
\end{equation}
We make the assumption that the $u_i$ values stem from an uncorrelated and stationary distribution with sample variance $\sigma^2$ (See the $1/f$ noise discussion in section~\ref{sec:Mot}).  The covariance matrix $\Cov{(U,U)}$ can then be written as the $N\times N$ identity matrix scaled by $\sigma^2$, however, the calculations presented here do not require a priori knowledge of $\sigma$. The elements of $\bi{C}$ are given by,
\begin{eqnarray}
c_{ij} &=& \frac{\sigma^2}{4^p}\sum_{r=0}^{N-p} a_{ir} a_{jr} \\
 &=& \frac{\sigma^2}{4^p}\sum_{r=0}^{N-p}(-1)^{2r}\left( {p \atop r-j} \right)  \left( {p \atop r-i} \right),
\label{eq:AAT1}
\end{eqnarray}
which simplifies using Vandermonde's identity and a change of variables to
\begin{equation}
c_{ij} = \frac{\sigma^2}{4^p}\left( {2p \atop p+i-j} \right) \mbox{~~~for $-p\le i-j \le  p$, zero otherwise}.
\label{eq:AAT2}
\end{equation}
Every covariance matrix is positive semi definite and symmetric, i.e. $\bi{C} = \bi{C}^T$~\cite{Strang}. Here, the latter is apparent from the inherent symmetry of the binomial coefficients. The covariance matrix $\bi{C}$ is no longer diagonal and the degree of correlation between the filtered data points is given by the off-diagonal elements.

\section{Extracting the signal from drift-free data}
\label{sec:drift_removal}

The filtered data set, $\bi{Y}$, contains the demodulated difference signal $\eta$, which can be retrieved using linear regression to estimate its mean value $\mu$. Given the form of the covariance matrix, the best linear unbiased estimator (BLUE) is obtained using Aitken's generalized least squares procedure~\cite{Aitken35, Rao}.

\begin{equation}
\mu = \frac{\bi{X}^T {\bi C}^{-1} } {\bi{X}^T  {\bi C}^{-1} \bi{X}} \bi{Y},
\end{equation}
where the design matrix $\bi{X}$, is a column vector of length $N-p$ with all elements equal to one. It should be noted that the matrix expressions in the denominators evaluate to scalar values. For notational purposes, we define the estimator of our signal $\mu$ to be $F_N^p$ in order to reflect its two parameter nature ($N$ data points and drift reduction up to order $p$). Since the result is independent of the data uncertainty $\sigma^2$, it can be expressed in terms of known quantities,
\begin{equation}
F_N^p= \frac{\bi{X}^T (\bi{A}\bi{A}^T )^{-1} } {\bi{X}^T  (\bi{A}\bi{A}^T )^{-1} \bi{X}} \bi{Y}.
\label{eqn:estm}
\end{equation}
Replacing $\bi{Y}$ by $\bi{AU}$, a weight vector $\bi{V}$ can be defined which allows calculation of the result directly from the input data set.
The estimate for the signal can be obtained in one simple step, the scalar product,
\begin{equation}
F_N^p = \bi{V} \bi{U} \;\;\;\mbox{with}\;\;\;\bi{V} =  \frac{\bi{X}^T (\bi{A}\bi{A}^T )^{-1} } {\bi{X}^T  (\bi{A}\bi{A}^T )^{-1} \bi{X}} \bi{A}.
\end{equation}

The weight vector $\bi{V}$ is unique to the values of $N$ and $p$ but can be reused as long as the length of the data set and the desired order of drift suppression remain unchanged.

\section{The variance of the estimator}
\label{sec:variance}

The variance of $F_N^p$ is calculated using error propagation. Let the factor multiplying $\bi{Y}$ in equation~(\ref{eqn:estm}) be
\begin{equation}
\bi{W} = \frac{\bi{X}^T (\bi{A}\bi{A}^T )^{-1} } {\bi{X}^T  (\bi{A}\bi{A}^T )^{-1} \bi{X}}.
\label{eqn:w}
\end{equation}
It then follows that the variance of $F_N^p$ is given by
\begin{equation*}
\Var(F_N^p) = \bi{W} \Cov{(Y,Y)} \bi{W}^T,
\end{equation*}
which simplifies to
\begin{equation}
\Var(F_N^p) = \frac{\sigma^2}{\bi{X}^T  (\bi{A}\bi{A}^T )^{-1} \bi{X}}  \label{eq:variance_wm}.
\end{equation}
The variance of the estimator is the inverse of the sum of all elements of the inverse correlation matrix scaled by $\sigma^2$.
If the value of $\sigma^2$ is not known, one estimate for $\Var(F_N^p)$ can be obtained by calculating the scatter in the estimator for a number of different data segments each of length $N$. If only a single set of $N$ measurements is available, $\sigma^2$ can be estimated from the sum of squares of the residuals about the regression mean $\mu$ as shown.
The quantity
\begin{equation*}
(\bi{Y}-\mu \bi{X})^T \bi{C}^{-1}(\bi{Y}-\mu \bi{X})
\nonumber
\end{equation*}
follows a $\chi^2$-distribution which has a mean value equal to the number of degrees of freedom $\nu = N - p -1$. Hence,
\begin{equation}
s^2 = (\bi{Y}-\mu \bi{X})^T (\bi{A}\bi{A}^T )^{-1}(\bi{Y}-\mu \bi{X})/(N-p-1) \label{eq:s}
\end{equation}
is an estimator of $\sigma^2$.
Using equation~(\ref{eq:variance_wm})  the variance in the mean of the filtered data can be compared to that of the un-filtered input data.  The corresponding variance in the mean of $N$ un-correlated data points is $\sigma^2 N^{-1}$  where the data variability $\sigma^2$ is common to both. Using the ratio of their variances, the statistical cost function becomes
\begin{equation}
c(N) =  \sqrt{\frac{N}{\bi{X}^T  (\bi{A}\bi{A}^T )^{-1} \bi{X}}}-1 .
\label{eq:cost}
\end{equation}
This gives the relative increase in statistical error resulting from the smaller correlated data set $(N-p)$ obtained in equation~(\ref{eq:trans}).

Figure \ref{fig:weightedmean} shows the statistical cost for three different filters labeled by their order $p$ versus the number of data points $N$. Cost values for a fixed value of $p$ oscillate between the upper and lower edges of their respective shaded bands with increasing $N$. The upper cost value for a filter of order $p$ is the same as the lower value for a filter with $p+1$.
\begin{figure}[ht!]
\begin{center}

\includegraphics[height=8cm, clip]{fig1}
\end{center}
\caption{The statistical cost is shown as a percentage versus the number of data points. The shaded bands are included to guide the eye and correspond to three different order filters. The cost values for a given filter order $p$ oscillate between the margins of their respective band with increasing $N$. This trajectory is shown using arrows for the first few points of the $p = 3$ filter.}
\label{fig:weightedmean}
\end{figure}

This oscillation can be understood as follows:
Each $y_i$ is a particular linear combination of the raw data as given in equation~(\ref{eq:lin}). If they have the same weighting factor, we can sum two at a time where their $i$ values differ by $m$.
\begin{equation}
y_i(p)  + y_{i+m}(p) = \frac{1}{2^p}\sum_{k=0}^p  (-1)^{i+k}  \left( {p \atop k} \right)( u_{i+k} + (-1)^m u_{i+k+m}).
\end{equation}
If $m$ is odd this is equivalent to applying the filter to the difference of the two points. The filter removes variations across the points that go as $i^q$ where $q$ is some power less than $p$. The difference $(i+m)^q - i^q$ is always one order less than $q$. In the process of taking the mean, for those points whose index differs by an odd number, the data is effectively filtered to the next higher order than the one chosen. In computing the mean, the  $y_i$ values are weighted by the vector $\bi{W}$ given in equation~(\ref{eqn:w}). Elements of this weight vector are each a sum over the corresponding column of the inverse correlation matrix $(\bi{A}\bi{A}^T )^{-1}$. Since the correlation matrix is symmetric about its diagonal it follows that the inverse correlation matrix is symmetric as well.\footnote{This follows from a property of invertible matrices $(\bi{C^{-1}})^T = (\bi{C^T})^{-1}$}  The weight vector therefore has this same symmetry about its center. If the length of the weight vector  is even, then the positions of its elements that are equal  in value differ by an odd number and the polynomial order filtered increases by 1. The oscillating behavior in the variance reflects this additional cancellation. When choosing either the filter order or the length of the data set, the statistical error is smallest when $N-p$ is odd.

The statistical precision of the mean of a measurement series scales as $\sigma/\sqrt{N}$ where $N$ is the number of measurements. If all measurements come from a population with variance  known a priori to be small enough, the desired statistical precision can be achieved with just a few measurements. In this case drift suppression is best achieved with sequence patterns as discussed in section~\ref{sec:Mot}.  In all other cases  precision measurements  require $N$  sufficiently large that the increase in  error from the filter would be minimal.
In  measurement scenarios with two different outcomes, equal amounts of data are taken in each state to insure  an unbiased result. The lowest order filter consistent with even values of $N$ has $p = 3$. It can be seen from figure~\ref{fig:weightedmean} that after  8  measurements the increase in statistical error is below a few percent. If only  linear drift  were removed, this same  statistical cost  could be achieved  with even fewer measurement. If 20 or more measurements are available the increase in error is minimal for all filter orders shown.

This analysis method can also be used to characterize the zero-point drift of the instrument. We start with the null hypothesis, i.e. that no significant drift exists so the results would be the same whether or not the filter is used. Without loss in generality we assume an even $N$. Hence, we choose to compare this with a filter of order $p = 3$ to achieve the smallest error.  An estimator for un-filtered data can be obtained by letting $p \rightarrow 0$. In this case elements of the $\bi{A}$ matrix ($N \times N $ ) are given by $a_{i,j} = (-1)^i\delta _{i,j}$ and the design matrix $\bi{X}$ has length $N$. Since both the filtered and un-filtered estimators are calculated from the same raw data, the null hypothesis can be checked. Using equation~(\ref{eq:s}) an estimator of the standard deviation for both cases can be obtained, $s^2_{p=0}$ and $s^2_{p=3}$. Multiplying these estimators  by their degrees of freedom $\nu = N-p-1$, yields variables that follow  $\chi^2$-distributions with $N-1$ and $N-4$ degrees of freedom respectively. In each case the $\chi^2$ probability is calculated and the null hypothesis is accepted if the $\chi^2$ probability of the filtered data is not significantly larger than that of the un-filtered data.

\section{Conclusion}

We have shown how zero-point drift can introduce a bias into the outcome of a measurement. Difference measurements in particular are  sensitive to drifts, especially when the time required for a measurement is long. If there are two states $A$ and $B$, choosing specific measurement sequences like $ABBA$ can completely eliminate this  bias from the difference.

The alternative method proposed here acquires data with a square wave modulation sequence $ABAB$... and recovers the signal with a post-hoc analysis. It is equally effective in removing linear drift and, in addition, can be configured to suppress zero-point drift of any desired polynomial order. Besides the simplicity of the measurement sequence, it offers the following advantages: (1) The measurement series is symmetric as every measurement is performed after the state has changed, (2) the highest order drift  suppressed can be selected after the data has been successfully taken, and (3) a missing data point has only a minor effect on  drift suppression. Its main disadvantage  is that it results in a slightly higher uncertainty in the end result when compared with schemes like $ABBA$ where data points are weighted with equal absolute values. The relative amount by which the uncertainty is increased depends on the number of data points $N$ and the order of the polynomial drift suppressed ($<p$) but is typically less than a few percent as shown in figure~\ref{fig:weightedmean}.

In the analysis process the set of weights (weight vector) used to average the $ABAB$... data series is calculated according to the prescription given for the estimator $F_N^p$.
Applying the method is straight forward due to the use of matrix equations. Typically the filter would be chosen to have an order $p$ of 2 or 3 depending on whether the data set size $N$ is odd or even, respectively. With $p$ and $N$ known one need only construct the appropriate $\bi{A}$ matrix and  design matrix $\bi{X}$.
There are software packages available such as MatLab, Mathematica, Wavemetric's IGOR, and Python with SciPi and NumPi to perform the matrix operations. The estimator  of the result is then the scalar product of the data vector with the weight vector.

We have chosen a polynomial expansion for the drift because the lowest order term contributes directly as a bias offset of the measurement. There may be specific environments where a different expansion would lead to smaller errors with fewer terms.
The matrix equations have been kept sufficiently general as to allow incorporating other functional forms of the drift. We note that the formalism presented here can be adopted to accommodate irregularly-spaced data. In this case each row of matrix, $\bf{A}$, has to be changed to meet the irregularly-spaced analog of equation~(\ref{eq:matrix}).

In conclusion, we believe that the filter presented here in combination with an $ABAB$... data sequence is an elegant and powerful tool for obtaining accurate measurements in the presence of drift.

\section*{References}
\bibliographystyle{unsrt}


\section{Introduction}

Quantization is born with quantum mechanics, as the fundamental attempt to establish a correspondence between the mathematical structures of classical and quantum mechanics, presented in the following table.

\begin{center}
  \begin{tabular}{c|c|c}
  & classical & quantum\\
  \hline
Phase space & symplectic manifold  $(\cM,\om)$ & Hilbert space $\calH$\\
\hline
Observables & Poisson algebra  $A\subset\Cinfty(\cM)$ & associative algebra $\cA\subset\ccL(\calH)$\\
\hline
Symmetries & Lie subalgebra $\fkg\subset\mathrm{ham}(\cM,\om)$ & Lie subalgebra $\fkg\subset\mathrm{u}(\calH)$.
  \end{tabular}\\
  \end{center}
One of the most celebrated quantization procedure is the geometric quantization \cite{Kos70,Sou70}, whose main draw-back is its too small set of quantizable observables. Equivariant quantization \cite{LOv99,DLO99} aims to overcome this issue for systems admitting a configuration space $M$ with a large enough group $G$ of (local) symmetries, as the projective or conformal group. More precisely, the inverse of the obtained quantization map is a $G$-equivariant symbol map on $M$, from differential operators to symmetric tensors.

We present here the equivariant quantization of spin systems whose configuration space is a spin manifold $M$ endowed with a metric $g$ of signature $(p,q)$.   This suppose to introduce a framework for classical mechanics of spin systems, namely the supercotangent bundle of $(M,g)$ \cite{BMa77,Get83} endowed with its canonical symplectic structure \cite{Rot90}. We recover some of the main objects of spin geometry via its geometric quantization, and the privileged status of conformal transformations of $(M,g)$ is highlighted. Restricting us to a conformally flat manifold, we describe explicitly the action of $\ro(p+1,q+1)$ on the space of classical and quantum observables, plus on a related space of tensors. We state then our main results on the existence and uniqueness of conformally equvariant quantization and superization, which are isomorphisms between these three $\ro(p+1,q+1)$-modules. Some applications are given. We refer to \cite{Mic09,Mic10a} for more details and proofs.

\section{From classical to quantum spin systems}
The quantum framework for spin systems with configuration space $(M,g)$ is well-known: the state space $\calH$ is obtained by completion of the space of sections of the spinor bundle $\sS\rightarrow M$, and we choose the algebra $\D(M,\sS)$ of spinor differential operators as space of quantum observables. Its usual algebra of symbols is $\Pol(T^*M)\otimes\Ga(\bbCl(M,g))$, the tensor product of the space of functions on $T^*M$, which are polynomial in the fiber variables, with the space of sections of the complex Clifford bundle of $(M,g)$. Replacing $\Ga(\bbCl(M,g))$ by its graded counterpart, namely the algebra of complex differential forms $\Om_\bbC(M)$, we end up with a superalgebra of functions on the supercotangent bundle of $M$. That provides us with the algebra of symbols for $\D(M,\sS)$ w.r.t. its bifiltration \cite{Get83}, as well as with the classical setting for a spin system on $M$.

\subsection{Supercotangent bundle and pseudomechanics}
The supercotangent bundle of the manifold $M$ is $\cM=T^*M\oplus\Pi TM$, i.e. the direct sum of the cotangent bundle and the tangent bundle with reverse parity. Thus, its superalgebra of functions is $\Cinfty(T^*M)\otimes\Om(M)$, generated locally by coordinates $(x^i,p_i,\xi^i)$, where $\xi^i$ identifies with $dx^i$. The general study of symplectic supermanifolds by Rothstein \cite{Rot90} proves that a symplectic structure on $\cM$ is equivalent to the data of a metric and a compatible connexion on $M$. As a consequence, to any pseudo-Riemannian manifold $(M,g)$ corresponds a canonical symplectic form $\om$ on $\cM$, given by
\begin{equation}
\om=d\a \quad \text{and} \quad \a=p_idx^i+\frac{\hbar}{2\bi}g_{ij}\xi^id^\nabla\xi^j,
\end{equation}
where $d^\nabla$ is the covariant differential w.r.t. the Levi-Civita connexion of $g$. The spin components are defined as $S^{ij}=\frac{\hbar}{\bi}\xi^i\xi^j$ so that, together with the Poisson bracket associated to $\om$, they generate a Lie algebra isomorphic to $\ro(p,q)$. The equations of motion of a spin particle in an exterior electromagnetic field can be easily recovered in that framework \cite{BMa77,Mic09}, and also the coupling of the spin with the gravitation \cite{Rav80,Mic09}. Thus, the Hamiltonian flows of the kinetic energy $g^{ij}p_ip_j$ leads to the Papapetrou's equations \cite{Pap51},
\begin{eqnarray*}
&&\dot{x}^j\nabla_j\dot{x}^i=-\half g^{ik}R(S)_{jk}\dot{x}^j,\\
&&\dot{x}^k\nabla_k S^{ij}=0,
\end{eqnarray*}
where the spin happens to be coupled with the curvature via $R(S)_{jk}=g_{im}R^i_{ljk}S^{jk}$, with $(R^i_{ljk})$ the components of the Riemann tensor.

\subsection{Spin geometry via geometric quantization}

Thanks to geometric quantization, we can built the main objects of spin geometry from $(\cM,\om)$ endowed with a polarization, i.e. a Lagrangian distribution. Upon topological restrictions on $M$, the vertical polarization of the cotangent bundle of $M$ can be completed by a maximal isotropic complex distribution for $g$ on $\Pi TM$, to give a polarization on $(\cM,\om)$. In the simplest case of a Riemannian metric, the geometric quantization leads then to the construction of Hitchin \cite{Hit74} for the spinor bundle $\sS$ of a pseudo-Hermitian manifold $(M,g)$, where spinor fields identify with antiholomorphic differential forms tensorized with square root of the volume form of $(M,g)$.

Besides, for Darboux coordinates $(x^i,\tp_i,\txi^i)$ of $(\cM,\om)$, the quantum map\footnote{we consider here that $\QG$ takes its values in usual spinor fields rather than spinor half-densities.} satisfies
\begin{equation}
\QG(x^i)=x^i,\qquad\QG(\tp_i)=\frac{\hbar}{\bi}\partial_i \quad \text{and}\quad \QG(\txi^i)=\frac{\tg^i}{\sqrt{2}},
\end{equation}
where $\tg^i$ is a Clifford matrix for the flat metric given by $(\eta_{ij})=\mathbb{I}_p\oplus-\mathbb{I}_q$. Let us remind that the vector fields on $M$ can be lifted to Hamiltonian vector fields on $T^*M$, giving rise to a tautological momentum map $J$, such that $J_X=p_iX^i$ for $X=X^i\partial_i\in\Vect(M)$. This map can be lifted to $\cM$ and quantized,
\begin{equation}
\QG(J_X)=\frac{\hbar}{\bi}\nabla_X,
\end{equation}
giving rise to an essential object of the spin geometry: the covariant derivative of spinors.

\subsection{Conformal geometry of the supercotangent and spinor bundles}

The symplectic structure of the supercotangent bundle $\cM$ depends on a metric $g$ on~$M$, as a consequence the natural lift to $\cM$ of $X\in\Vect(M)$ does not preserve the potential $1$-form $\a$. Moreover, the condition $L_{\tilde{X}}\a=0$ does not uniquely determine the lift $\tilde{X}$ of~$X$. One way out is to impose the further condition $L_{\tilde{X}}\beta\propto\beta$, where $d\b=g_{ij}d^\nabla\xi^i\wedge dx^j$ is the odd symplectic form on $\Pi TM$. The lift is then unique but exists only if $X$ is a conformal Killing vector field. We refer to it later as $\tilde{X}$.
Let us notice that, as expected, the momentum of an infinitesimal rotation $X_{ij}$ is,
\begin{equation}
\cJ_{X_{ij}}=p_ix_j-x_jp_i+S_{ij}.
\end{equation}
Remarkably, the quantization of the momentum map $\cJ$ leads to
\begin{equation}
\QG(\cJ_X)=\frac{\hbar}{\bi}\sL_X,
\end{equation}
the Lie derivative of spinors introduced by Kosmann \cite{Kos72}, well-defined precisely for conformal Killing vector fields. Via geometric quantization, we get thus the conformal geometry of the spinor bundle of $(M,g)$ out of the one of its supercotangent bundle.

\section{Conformally equivariant quantization}

We suppose from now on that $(M,g)$ is a conformally flat spin manifold, i.e. $g_{ij}=F\eta_{ij}$ for some positive function $F$. The conformal Killing vector fields on $(M,g)$ generate then a Lie algebra isomorphic to $\ro(p+1,q+1)$, which is the one of infinitesimal conformal transformations of $(\bbR^n,\eta)$. The aim is to compare the action of those vector fields on $\sD(M,\sS)$, the algebra of spinor differential operators, with those on its algebras of symbols, i.e. to compare their $\ro(p+1,q+1)$-module structures.

\subsection{Conformal geometry of spinor differential operators and of their symbols}
Let us define the space of $\l$-densities by $\cF^\l=\Ga(|\Lambda^n T^*M|^{\otimes\l})$, with $\l\in\bbR$ and $n=\dim M$.
Instead of $\sD(M,\sS)$, we will rather study the two parameters family of $\ro(p+1,q+1)$-modules $(\Dslm)$. Each of those modules is defined as the space of differential operators $D:\Ga(\sS)\otimes\cF^\l\rightarrow\Ga(\sS)\otimes\cF^\m$, endowed with the adjoint action
\begin{equation}
\LD D=\sLl D-\sLm D,
\end{equation}
where $\sLl=\sL_X+\l\Div(X)$ is the action of $X$ on $\Ga(\sS)\otimes\cF^\l$. The corresponding ${\ro(p+1,q+1)}$-module of classical observables is naturally the space $\cS^\d[\xi]=\Pol(T^*M)\otimes\Om_\bbC(M)\otimes\cF^\d$, with $\d=\mu-\l$, endowed with the Hamiltonian action
\begin{equation}
\Ld =\tilde{X}+\d\Div(X).
\end{equation}
The explicit expressions of these both classical and quantum actions have been computed in \cite{Mic10a}, showing that $\Dslm$ and $\cS^\d[\xi]$ are filtered modules, by the order of differential operators and the degree in the $p$ variables respectively. Their common associated graded module is the module of tensorial symbols $\cT^\d[\xi]=\bigoplus_{\k=0}^n\Pol(T^*M)\otimes\Omega^{\k}_\bbC\otimes\cF^{\d-\frac{\k}{n}}$, endowed with the natural action on weighted tensors. It identifies to $\Pol(T^*M)\otimes\Om_\bbC(M)$ as an algebra and to the usual space of symbols $\Pol(T^*M)\otimes\Ga(\bbCl(M,g))$ as a module.

\subsection{Main results}
Let us begin with few definitions. A map is called conformally equivariant if it is an isomorphism of $\ro(p+1,q+1)$-modules. Besides, as $\Pol(T^*M)$ is a submodule of $\cT^0[\xi]$,  a linear isomorphism $\cT^\d[\xi]\rightarrow\cS^\d[\xi]$ which preserves the principal symbol\footnote{that is the higher order term of an element w.r.t. a filtration, here it identifies to a tensor over $M$.} is named a superization, whereas a linear isomorphism $\cS^\d[\xi]\rightarrow \Dslm$ preserving the principal symbol is a quantization, since it relates classical and quantum observables.

\begin{thm}
There exists $ I^e_{\rS}\subset I_{\rS}\subset\bbQ_+^*$ such that the conformally equivariant superization
$
\rST^\d:\cT^\delta[\xi]\rightarrow \cS^\delta[\xi],
$
exists if $\d\not\in I^e_{\rS}$ and is unique if $\d\not\in I_{\rS}$.
\end{thm}
\begin{thm}
Let $\d=\mu-\l\in\bbR$. There exists $I^e_{\cQ}\subset I_{\cQ}\subset\bbQ_+^*$ such that the conformally equivariant quantization
$
\cQlm:\cS^\delta[\xi]\rightarrow \Dslm,
$
exists and is unique if $\d\not\in I_{\cQ}$, and exists for at least one value of $\l\in\bbR$ if $\d\not\in I^e_{\cQ}$.
\end{thm}
The values of $\d$ for which existence or uniqueness of $\cQlm$ is lost are called resonances. The fact that $\d=0$ is not a resonance is crucial, thus the conformally equivariant quantization $\cQ^{\half,\half}$  extends uniquely the quantization map provided by geometric quantization.
Let us remark that $\cQlm\circ\rST^\d$ is a particular case of AHS-equivariant quantization \cite{CSi09}, but each single map deserves interest, at least from a physical point of view.

The idea of the proofs is the same than in the spinless case \cite{DLO99}, and relies on the use of the Casimir operators of each module. To be concrete, we name $\tC$ and $\Cd$ those of $\cT^\d[\xi]$ and $\cS^\d[\xi]$. If the conformally equivariant superization exists, then $\Cd\rST^\d=\rST^\d\tC$, and in particular every eigenvectors of $\tC$ is send to an eigenvector of $\Cd$ with the same eigenvalue and the same principal symbol. The main point is to prove that the eigenvectors of $\tC$ are uniquely determined by their eigenvalue and their principal symbol. Then, if this is the case for those of $\Cd$ too, we get the uniqueness of the superization. As $\Cd$ is equal to $\tC$ plus an operator lowering the degree in the $p$ variable, this is simply checked by the resolution of triangular systems, and the resonant values of $\d$ are precisely those leading to degenerated systems.
If unique, the map constructed in this way is easily proved to be conformally equivariant.

\subsection{Some applications}

We give now two applications that we hope to investigate further in forthcoming papers. The first one relies on the explicit formulas for the conformally equivariant superization that we determine in \cite{Mic09} for symbols of degree $1$ in $p$. Let us recall that a (conformal) Killing-Yano tensor on $M$ is a skew-symmetric tensor describing higher symmetries of $(M,g)$. As Killing tensors, it generates constant of motion but for spin particles  \cite{GRV93}. In fact, a correspondence has been obtained in \cite{Tan95} between Killing-Yano tensors and classical supercharges, which happens to be generalized by the conformally equivariant superization.
 \begin{thm}
Let $f$ be a skew-symmetric tensor and $P_f=f^i_{j_1\ldots j_{\k-1}}\xi^{j_1}\ldots\xi^{j_{\k-1}}p_i$ the associated tensor symbol. Denoting by $\Delta=p_i\xi^i$, we have
\begin{equation}
\{\Delta,\rST^0(P_f)\}=0\;(\propto\Delta ) \Longleftrightarrow f \text{ is a (conformal) Killing-Yano tensor}.
\end{equation}
\end{thm}

The second application deals with the conformal invariants of the $\ro(p+1,q+1)$-modules that we have introduced. Let $R=g^{ij}p_ip_j$. The Weyl theory of invariants together with the explicit actions of $\ro(p+1,q+1)$ on each module leads to the following Theorem.
\begin{thm}
The conformal invariants of each family of modules are
\begin{enumerate}
\item $\Delta^a\, R^s \in\cT^{\frac{2s+a}{n}}[\xi],\quad$ where $s\in\bbN$ and $a=0,1$.
\item $\Delta\, R^s \in\cS^{\frac{2s+a}{n}}[\xi],\quad$ where $s\in\bbN$.
\item $\cQlm(\Delta\, R^s) \in\sD^{\frac{n-2s-1}{2n},\frac{n+2s+1}{2n}},\quad$ where $s\in\bbN$.
\end{enumerate}
\end{thm}
In particular, we get the Dirac operator as the lower order conformal invariant of $(\Dslm)$. Since the maps $\rST^\d$ and $\cQlm$ preserve conformal invariance,  we deduce from the last Theorem that $\rST^\frac{2s}{n}$ does not exist and that the module $\sD^{\frac{n-2s-1}{2n},\frac{n+2s+1}{2n}}$ is exceptional: this the only one in the family $(\sD^{\l,\l+\frac{2s+1}{n}})$ to be isomorphic to its space of symbols $\cS^\frac{2s+1}{n}[\xi]$. As a consequence, $\d=\frac{2s+1}{n}$ are resonances. In fact, every resonances correspond to some conformal invariants, that will be the matter of our next paper.

 \subsection*{Acknowledgements}
 It is a pleasure to acknowledge Christian Duval for his essential guidance, and invaluable discussions. Special thanks are due to Valentin Ovsienko for his constant interest in this work.

\bibliographystyle{plain}


\section{Introduction}\label{sec1}
To introduce the main notions of the
triangular factorization (see  \cite{3, 5, 8, 14, 15, 20})
consider a Hilbert space $L^{2}(a,\, b)\, $
$(-\infty{\leq}a<b{\leq}\infty)$. The orthogonal projectors $P_{\xi}$ in $L^{2}(a,\, b)$
are defined by the relations
\[
\big(P_{\xi}f\big)(x)=f(x)  \, \, \mathrm{for} \, \,  a<x<\xi,  \,\,
\big(P_{\xi}f\big)(x)=0 \, \, \mathrm{for} \, \,  \xi <x<b
\, \, \big( f \in L^{2}(a,b)\big).
\]
Denote the identity operator by $I$.
\begin{defn}\label{Definition 1.1.}
A bounded operator $S_{-}$ on
 $L^{2}(a,b)$  is called lower triangular if for
every
 $\xi$ the relations \begin{equation}\label{1.1}
  S_{-}Q_{\xi}=Q_{\xi}S_{-}Q_{\xi} , \end{equation} where
  $Q_{\xi}=I-P_{\xi}$, are true. The operator  $S_{-}^{\ast}$ is called
  upper triangular.
\end{defn}

\begin{defn} \label{Definition 1.2.}
 A bounded, positive definite and invertible
operator $S$ on $L^{2}(a,b)$ is said to admit a left (right)
triangular factorization if it can be represented in the form
\begin{equation}\label{1.2}
S=S_{-}S_{-}^* \quad (S=S_{-}^*S_{-}),\end{equation} where $S_{-}$ and $S_{-}^{-1}$ are
bounded and lower triangular operators.
\end{defn}
Further, we often write factorization meaning a left triangular factorization.

In paper \cite{20} (see p. 291)  we formulated necessary and sufficient conditions
under which the positive definite operator $S$ admits  a triangular
factorization. The factorizing operator $S_{-}^{-1}$ was
constructed in the explicit form.
We proved that a wide class of  operators admits a triangular factorization \cite{20}.

 D. Larson proved  \cite{8} the \emph{existence}  of positive definite and invertible
 but non-factorable operators. In the present article we construct  \emph{concrete examples} of
 such operators. In particular, the following operator
 \begin{equation}\label{1.3}   Sf=f(x)-{\mu}\int_{0}^{\infty}\frac{\mathrm{sin}{\pi}(x-t)}{{\pi}(x-t)}f(t)dt,\quad
 f(x){\in}L^{2}(0,\infty),\quad 0<\mu<1
 \end{equation}
 is positive definite and invertible but non-factorable.
 Using  positive definite and invertible but non-factorable operators  we have managed to substitute pure  existence  theorems \cite{8} by concrete examples in the
 well-known problems posed by J.R. Ringrose \cite{13}, R.V. Kadison and I.M. Singer \cite{5}.
 We note that Kadison-Singer problem was posed independently by
 I. Gohberg and M.G. Krein \cite{4}.

 The non-factorable operator $S$, which is defined by formula \eqref{1.3}, is used in a number of theoretical and applied problems (in optics \cite{7},
random matrices \cite{23}, generalized stationary processes \cite{11,12}, and Bose gas theory \cite{10}). The results obtained in  the paper are interesting
 from this point of view too.
 \section{A special class of operators \\ and corresponding differential systems}\label{sec2}

In this section we consider operators $S$ of the form
 \begin{equation}\label{2.1}
 Sf=f(x)-{\mu}\int_{0}^{\infty}h(x-t)f(t)dt,\quad
 f(x){\in}L^{2}(0,\infty),
 \end{equation}
 where $\mu=\overline{\mu}$ and  $h(x)$ admits representation
 \begin{equation}\label{2.2}
   h(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{ix\lambda}\rho(\lambda)d\lambda.
 \end{equation}
 We suppose that the function $\rho(\lambda)$  satisfies the following conditions

1. The function $\rho(\lambda)$ is real and bounded
\begin{align}\label{d0}&
|\rho(\lambda)|{\leq}U^{2},\quad U>0 \quad  (-\infty<\lambda<\infty).
\end{align}

2.
  $\rho(\lambda)=\rho(-\lambda){\in}L(-\infty,\infty)$.\\
  Hence,  the function $h(x)\quad
  (-\infty<x<\infty)$ is continuous and real. The corresponding operator
  \begin{equation}\label{2.3}
   Hf=\int_{0}^{\infty}h(x-t)f(t)dt
  \end{equation}
  is self-adjoint and bounded, where $||H||{\leq}U.$
  We introduce the operators
   \begin{equation}\label{2.4}   S_{\xi}f=f(x)-{\mu}\int_{0}^{\xi}h(x-t)f(t)dt,\quad
 f(x){\in}L^{2}(0,\xi),\quad 0<\xi<\infty.\end{equation}
 The following statement is true.

 \begin{prop}\label{Proposition 2.1.}
 {If $-1/U<\mu<1/U$, then the operator}
 $S_{\xi}$, which is defined by formula \eqref{2.4}, is positive definite, bounded
 and invertible.
 \end{prop}

 Hence, we have
 \begin{equation}\label{2.5}
   S_{\xi}^{-1}f=f(x)+\int_{0}^{\xi}R_{\xi}(x,t,\mu)f(t)dt.
   \end{equation}
 The function $R_{\xi}(x,t,\mu)$ is jointly continuous in $x,t,\xi , \mu$.
 M.G. Krein (see \cite{4}, Ch. IV, Section 7) proved that
 \begin{equation}\label{2.6}
  S_{b}^{-1}=(I+V_{+})(I+V_{-}),\quad
 0<b<\infty,
 \end{equation}
 where the operators $V_{+}$ and $V_{-}$ are defined in $L^{2}(0,b)$ by the
 relations
\begin{equation}\label{2.7}
 \big(V_{+}^{\ast}f\big)(x)=\big(V_{-}f\big)(x)=\int_{0}^{x}R_{x}(x,t,\mu)f(t)dt.
\end{equation}
 The Krein's formula \eqref{2.6} is true for the Fredholm class of operators. The operator $S_{b}$ belongs to this class. The kernel of the operator
$V_{-}$ does not depend of $b$. Hence, if the operator $S$ admits the factorization,
then formula \eqref{2.6} holds for the case $b=\infty$ too, i.e.
\begin{equation}\label{2.8}   S^{-1}=(I+V_{+})(I+V_{-}).
\end{equation}

\begin{rem}\label{Remark 2.1.}
Relation \eqref{2.8} also follows  from Theorem 2.1 in the  paper \cite{20}.
\end{rem}

Let us introduce the function
\begin{equation}\label{2.9}
  q_{1}(x)=1+\int_{0}^{x}R_{x}(x,t,\mu)dt.
\end{equation}
Using the relation  $R_{x}(x,t,\mu)=R_{x}(x-t,0,\mu)$ (see \cite{4}, formula (8.12)), we obtain
\begin{equation}\label{2.10}
  q_{1}(x)=1+\int_{0}^{x}R_{x}(u,0,\mu)du.
  \end{equation}
According to the
well-known Krein's formula (\cite{4}, Ch. IV, formulas (8.3) and (8.14)) we have
\begin{equation}\label{2.11}
 q_{1}(x)=\exp\left\{\int_{0}^{x}R_{t}(t,0,\mu)dt\right\} .
 \end{equation}
Together with $q_{1}(x)$ we shall consider the function
\begin{equation}\label{2.12}
  q_{2}(x)=M(x)+\int_{0}^{x}M(t)R_{x}(x,t,\mu)dt,
\end{equation}
where
\begin{equation}\label{2.13}
  M(x)=\frac{1}{2}-{\mu}\int_{0}^{x}h(s)ds.
  \end{equation}
The functions $q_{1}(x)$ and $q_{2}(x)$ generate the $2{\times}2$ differential system
\begin{equation}\label{2.14}
  \frac{dW}{dx}=izJH(x)W,\quad W(0,z)=I_{2}.
  \end{equation}
  Here $W(x,z)$ and $H(x)$
are   $2{\times}2$ matrix functions and $J$ is a $2{\times}2$ matrix:
\begin{equation}\label{2.15}
  H(x)=\left[\begin{array}{cc}
                       q_{2}^{2}(x) & 1/2 \\
                       1/2 & q_{1}^{2}(x)
                     \end{array}\right], \quad
                     J=\left[\begin{array}{cc}
                        0 & 1 \\
                        1 & 0
                      \end{array}\right].
\end{equation}
Note that according to  \cite{19} (see formulas (53) and (56) therein) we have:
\begin{equation}\label{2.16}
  q_{1}(x)q_{2}(x)=1/2.
  \end{equation}
It is easy to see that
\begin{equation}\label{2.17}
 JH(x)=T(x)PT^{-1}(x),
 \end{equation}
where
\begin{equation}\label{2.18}
  T(x)=\left[\begin{array}{cc}
                             q_{1}(x) & -q_{1}(x) \\
                             q_{2}(x) & q_{2}(x)
                           \end{array}\right], \quad
P=\left[\begin{array}{cc}
    1 & 0 \\
    0 & 0
  \end{array}\right].
\end{equation}
Consider the matrix function
\begin{equation}\label{2.19}   V(x,z)=e^{-ixz/2}T^{-1}(x)W(x,z)T(0).\end{equation}
Due to \eqref{2.14}-\eqref{2.19} we get
\begin{equation}\label{2.20}   \frac{dV}{dx}=(iz/2)jV+\Gamma(x)V, \quad V(0)=I_{2},\end{equation}
where
\begin{equation}\label{2.21}   \Gamma(x)=\left[\begin{array}{cc}
                             0 & B(x) \\
                             B(x)& 0
                           \end{array}\right], \quad
j=\left[\begin{array}{cc}
          1 & 0 \\
          0 & -1
        \end{array}\right],
\end{equation}
\begin{equation}\label{2.22}
  B(x)=\frac{q_{1}^{\prime}(x)}{q_{1}(x)}=R_{x}(x,0,\mu).
  \end{equation}
Let us introduce the functions
\begin{equation}\label{2.23}   \Phi_{n}(x,z)=v_{1n}(x,z)+v_{2n}(x,z)\quad (n=1,2),\end{equation}
\begin{equation}\label{2.24}   \Psi_{n}(x,z)=i[v_{1n}(x,z)-v_{2n}(x,z)] \quad (n=1,2),\end{equation}
where $v_{in}(x,z)$ are elements of the matrix function $V(x,z)$.
It follows from \eqref{2.20} that
\begin{align}\label{2.25} & \frac{d\Phi_{n}}{dx}=(z/2)\Psi_{n}-B(x)\Phi_{n},\quad \Phi_{1}(0,z)=\Phi_{2},(0,z)=1, \\
& \label{2.26}   \frac{d\Psi_{n}}{dx}=-(z/2)\Phi_{n}+B(x)\Psi_{n},
\quad \Psi_{1}(0,z)=-\Psi_{2}(0,z)=i .\end{align}
Consider again the differential system \eqref{2.14} and the  solution $W(x,z)$ of this system.
The element $w_{1,2}(\xi,z)$ of the matrix function $W(\xi,z)$
can be represented in the form  (see \cite{17}, p. 54, formula \eqref{2.5})
\begin{equation}\label{2.27}   w_{1,2}(\xi,z)=iz\left(\left( I-zA\right)^{-1}1,S_{\xi}^{-1}1\right)_{\xi},\end{equation}
where the operator $A$ has the form
\begin{equation}\label{2.28}   Af=i\int_{0}^{x}f(t)dt.\end{equation}
It is well-known that
\begin{equation}\label{2.29}   (I-zA)^{-1}1=e^{izx}.\end{equation}
We can obtain a representation of $W(\xi,z)$ without using the operator $S_{\xi}^{-1}$.
Indeed, it follows from
\eqref{2.19},  \eqref{2.23}, and \eqref{2.24} that
\begin{equation}\label{2.30}   W(x,z)=(1/2)e^{ixz/2}T(x)\left[\begin{array}{cc}
                                            \Phi_{1}-i\Psi_{1}& \Phi_{2}-i\Psi_{2} \\
                                         \Phi_{1}+i\Psi_{1} & \Phi_{2}+i\Psi_{2}
                                          \end{array}\right]T^{-1}(0).\end{equation}
According to equality \eqref{2.10} we have $q_{1}(0)=1$. Due to \eqref{2.18} we infer
\begin{equation}\label{2.31}   T(0)=\left[\begin{array}{cc}
                             1& -1 \\
                             1/2 & 1/2
                           \end{array}\right],\quad T^{-1}(0)=\left[\begin{array}{cc}
                             1/2& 1 \\
                             -1/2 & 1
                           \end{array}\right].\end{equation}
Further we plan to use a Krein's result from \cite{6}. For that purpose we introduce the functions
\begin{equation}\label{2.32}
P(x,z)=e^{ixz/2}[\Phi(x,z)-i\Psi(x,z)]/2,
\end{equation}
\begin{equation}\label{2.33}
P_{\ast}(x,z)=e^{ixz/2}[\Phi(x,z)+i\Psi(x,z)]/2,
\end{equation}
where
\begin{equation}\label{2.34}
\Phi(x,z)=\Phi_{1}(x,z)+\Phi_{2}(x,z),\quad \Psi(x,z)=\Psi_{1}(x,z)+\Psi_{2}(x,z).
\end{equation}
Using \eqref{2.25}, \eqref{2.26} and  \eqref{2.32},   \eqref{2.33} we see that the pair $P(x,z)$ and $P_{\ast}(x,z)$ is a solution
of the following Krein system
\begin{equation}\label{2.35}   \frac{dP}{dx}=izP-B(x)P_{\ast},\quad  \frac{dP_{\ast}}{dx}=-B(x)P,
\end{equation} where
\begin{equation}\label{2.36}   P(0,z)=P_{\ast}(0,z)=1.\end{equation}
It follows from \eqref{2.32} and \eqref{2.33} that
\begin{equation}\label{2.37}   P(x,z)-P_{\ast}(x,z)=-ie^{ixz/2}\Psi(x,z).\end{equation}

\section{Non-factorable positive definite operators, a sufficient condition}\label{sec3}

We assume that the following relation  is true:
\begin{equation}\label{3.1}   M(x)=(1-\mu)/2+q(x),\quad q(x){\in}L^{2}(0,\infty),\end{equation}
where the function $M(x)$ is defined by \eqref{2.13}.
Condition \eqref{3.1} can be rewritten in an equivalent form:
\begin{equation}\label{3.2}   \int_{0}^{\infty}h(x)dx=1/2,\quad \int_{x}^{\infty}h(x)dx{\in}L^{2}(0,\infty).
\end{equation}
Now, we need the relations (see \cite{16}, Ch. 1, formulas (1.37) and (1.44)):
\begin{equation}\label{3.3}   S_{\xi}1=M(x)+M(\xi-x),\quad S_{\xi}=U_{\xi}S_{\xi}U_{\xi}, \end{equation}
where $U_{\xi}f(x)=\overline{f(\xi-x)},\quad 0{\leq}x{\leq}\xi.$
  It follows from \eqref{3.1} and \eqref{3.3} that
\begin{equation}\label{3.4}   S_{\xi}1=1-\mu +q(x)+U_{\xi}q(x).\end{equation}Hence the relation
\begin{equation}\label{3.5}   S_{\xi}^{-1}1=\frac{1}{(1-\mu)}[1-r_{\xi}(x)-U_{\xi}r_{\xi}(x)]\end{equation}
is true. Here  $r_{\xi}(x)=S_{\xi}^{-1}q(x).$ Using formulas \eqref{2.27},  \eqref{3.1}, and \eqref{3.5},
we obtain the following representation of  $w_{1,2}(\xi,z)$.
\begin{lem}\label{Lemma 3.1.}
{The function $w_{1,2}(\xi,z)$ has the form}
\begin{equation}\label{3.6}   w_{1,2}(\xi,z)=e^{iz\xi}G(\xi,z)-\overline{G(\xi,\overline{z})},\end{equation}
{where}
\begin{equation}\label{3.7}    G(\xi,z)=\frac{1}{1-\mu}\left[ 1-iz\int_{0}^{\xi}e^{-izx}r_{\xi}(x)dx\right].\end{equation}
\end{lem}
Note that the operator $S$ is positive definite, bounded and invertible.
 According to \eqref{2.6} we have
\begin{equation}\label{3.8}   Q(x)=(I+V_{-})q(x){\in}L^{2}(0,\infty).\end{equation}
Hence, there exists a sequence $x_{n}$ such that
\begin{equation}\label{3.9}   Q(x_{n}){\to}0,\quad x_{n}{\to}\infty.\end{equation}
Now, we prove the following statement.
\begin{lem}\label{Lemma 3.2.}
{Let relation \eqref{3.9} be true. Then we have}
\begin{equation}\label{3.10}
  \lim_{x_{n}{\to}\infty}{q_{1}(x_{n})}=\frac{1}{\sqrt{1-\mu}} .
\end{equation}
\end{lem}
\begin{proof}
In view of \eqref{2.9}, \eqref{2.12}, and \eqref{3.1} we get
\begin{equation}\label{3.11}
  q_{2}(x)=q_{1}(x)(1-\mu)/2+Q(x).
  \end{equation}
Taking into account the relation
 $q_{1}(x)q_{2}(x)=1/2$ (see \cite{19}, formulas (53) and  (56)),  we obtain   the equality
\begin{equation}\label{3.12}   1/2=q_{1}^{2}(x)(1-\mu)/2+q_{1}(x)Q(x).\end{equation}
 Formula  \eqref{3.10} follows directly from \eqref{3.9}, \eqref{3.12}, and inequality $q_{1}(x)>0$.
 \end{proof}

It follows from \eqref{2.18} and \eqref{3.10} that
 \begin{equation}\label{3.15}   T(x_{n}){\to}\left[\begin{array}{cc}
                            C &-C \\
                            1/2C & 1/2C
                          \end{array}\right],\quad x_{n}{\to}\infty,\quad C=1/\sqrt{(1-\mu)}.
                          \end{equation}
Hence, in view of \eqref{2.31}, \eqref{2.32}, \eqref{2.34}, and \eqref{3.15} the following assertion is true.
\begin{lem}\label{Lemma 3.3.} Let $x_n$ tend to $\infty$. Then,
$w_{1,2}$ has the following asymptotics
\begin{equation}\label{3.16}   w_{1,2}(x_{n},z)=-iCe^{ix_{n}z/2}\Psi(x_{n},z)\big(1+o(1)\big).
\end{equation}
\end{lem}

\begin{lem} \label{Lemma 3.4.}
Suppose that the operator $S$ admits a factorization.
Then we have
\begin{equation}\label{3.17}    \lim_{\xi{\to}\infty}e^{-iz\xi}w_{1,2}(\xi,z)=G(z),\quad {\Im}z<0,\end{equation}
\begin{equation}\label{3.18}    \lim_{\xi{\to}\infty}w_{1,2}(\xi,z)=-\overline{G(\overline{z})},\quad {\Im}z>0.\end{equation}
{where}
\begin{equation}\label{3.19}    G(z)=\frac{1}{1-\mu}[1-iz\int_{0}^{\infty}e^{-izx}r(x)dx], \quad
r(x)=S^{-1}q(x).\end{equation}
\end{lem}
\begin{proof}
According to \eqref{2.8} we have
$S_{-}^{-1}=I+V_{-}$, where $V_{-}$ is defined by \eqref{2.7}.  Hence,  the operator function
 $S_{\xi}^{-1}$ strongly converges to the operator  $S^{-1}$ when $\xi{\to}\infty$. Then
 the function $r_{\xi}(x)=S_{\xi}^{-1}q(x)$ strongly converges to $r(x)=S^{-1}q(x)$, when
 $\xi{\to}\infty$, and
 $r(x){\in}L^{2}(0,\infty)$. Using \eqref{3.6} and \eqref{3.7} we obtain relations \eqref{3.17} and \eqref{3.18}. The lemma is proved.
 \end{proof}

 From Lemma \ref{Lemma 3.4.} we derive the following important assertion.

 \begin{prop}\label{Proposition 3.1.}
 {If at least one of the equalities \eqref{3.17} and  \eqref{3.18} is not true, then
the corresponding operator $S$ does not admit factorization.}
\end{prop}

Note that a new approach to the notion of the limit of a function was used in Lemma \ref{Lemma 3.2.}. Namely,
we introduce a continuous function $F(x)$, which belongs to $L(0,\infty)$, and consider
sequences $x_{n}{\to}\infty$, such that
\begin{equation}\label{d1} F(x_{n}){\to}0.\end{equation}
\begin{defn} \label{Definition 3.1.} We  say that the function $f(x)$ tends to $A$
almost sure $(a.s.)$ if relation \eqref{d1} implies
\begin{equation}\label{d2} f(x_{n}){\to}A,\qquad x_{n}{\to}\infty.\end{equation}
\end{defn}
Equality (3.10) can be written in the form
\begin{equation}\label{d3} \lim_{x{\to}\infty}{q_{1}(x)}=\frac{1}{\sqrt{1-\mu}},\quad  a.s. \end{equation}
\begin{rem}\label{Remark 3.1.}
From heuristic point of view "almost all" sequences $x_{n}{\to}\infty$
satisfy relation \eqref{d1}. This is the reason of using the probabilistic term
"almost sure".
\end{rem}

\section{A class of non-factorable positive definite operators}\label{sec4}

Introduce a partition
\begin{equation}\label{4.1}
0=a_{0}<a_{1}<...<a_{n}=a,
\end{equation}
and consider the function
$\rho(\lambda)=\rho(-\lambda)$ such that
\begin{equation}\label{4.2}
\rho(\lambda)= \left\{
\begin{array}{rl}
0, & a \leq \lambda, \\
b_{k-1}, & a_{k-1}{\leq}\lambda<a_{k},
\end{array}
\right.
\end{equation}
where
\begin{equation}\label{4.3}
b_{0}=1;\quad -1{\leq}b_{k}{\leq}1\quad (0<k{\leq}n-1).\end{equation}
In the case of $\rho$ given by \eqref{4.2} and  \eqref{4.3} we can put  $U=1$ in \eqref{d0}.
 Further we investigate the operators $S$, which are defined by formulas \eqref{2.1}, \eqref{2.2}, and
\eqref{4.2}. The spectral function $\sigma(\lambda)$ of the corresponding system
\eqref{2.35} is absolutely continuous and such that (see \cite{6}):
\begin{equation}\label{4.4}   \sigma^{\prime}(\lambda)=[1-\mu\rho(\lambda)]/(2\pi).\end{equation}
\begin{rem}\label{Remark 4.1.}
The operators $S$, which are defined by formulas \eqref{2.1}, \eqref{2.2}, and
\eqref{4.2}, appear in the theory of generalized stationary processes of
white noise type (see \cite{11,12}). If $n=1$ and $a_{1}=\pi$,  then the corresponding operator $S$ has  the form \eqref{1.3}.
\end{rem}

It follows from \eqref{2.2} and \eqref{4.2} that
\begin{equation}\label{4.5}   h(x)=\frac{1}{\pi}\sum_{k=1}^{n}b_{k}\frac{{\sin}a_{k}x-{\sin}a_{k-1}x}{x}.
\end{equation}
According to \eqref{4.4} we have
\begin{equation}\label{4.6}   \int_{-\infty}^{\infty}\frac{{\log}\sigma^{\prime}(u)}{1+u^{2}}du<\infty.
\end{equation}
It follows from \eqref{4.7}  (see \cite{6}) that
\begin{equation}\label{4.7}   \int_{0}^{\infty}|P(x,z_{0})|^{2}dx<\infty,\quad {\Im}z_{0}>0.\end{equation}
Hence, there exists a sequence $x_{n}$ such that
\begin{equation}\label{4.8}   |P(x_{n},z_{0})|^{2}{\to}0,\quad x_{n}{\to}\infty.\end{equation}
Now, we use the  corrected form of Krein's theorem (see \cite{6,21}):

\begin{prop}\label{Proposition 4.1.}
$1)$ { There exists the limit}
\begin{equation}\label{4.9}   \Pi(z)=\lim_{x_{n}{\to}\infty} P_{\ast}(x_{n},z),\end{equation}
{where the convergence is uniform at any bounded closed set  of the upper half-plane
$\Im z>0$. }\\
$2)$ {The function $\Pi(z)$ can be represented in the form}
\begin{equation}\label{4.10}    \Pi(z)=\frac{1}{\sqrt{2\pi}}\exp\left\{\frac{1}{2i\pi}\int_{-\infty}^{\infty}
\frac{1+tz}{(z-t)(1+t^{2})}(\log{\sigma}^{\prime}(t))dt+i\alpha\right\} ,\end{equation}
{where} $\alpha=\overline{\alpha}$. Here  $\sigma$ is the spectral function of system \eqref{2.35},
which corresponds to $\rho$ given by \eqref{4.2} and \eqref{4.3}, that is, this $\sigma$
is defined by \eqref{4.4}.
\end{prop}
\begin{rem}\label{Remark 4.2.}
The function $|Q(x)|^{2}+|P(x,z_{0})|^{2}$ belongs to
the space $L(0,\infty)$. Hence, there exists a sequence $x_{n}$ such that
relations \eqref{3.9} and  \eqref{4.8} are true simultaneously.
\end{rem}

If \eqref{4.5} holds, then the following conditions are fulfilled:
\begin{equation}\label{4.11}    0<\delta{\leq}||S||{\leq}\Delta<\infty,\qquad \int_{0}^{\infty}|h(x)|^{2}dx<\infty.\end{equation}
Therefore,
 in formula \eqref{4.10} we have (see \cite{19}, Proposition 1):
\begin{equation}\label{4.12}    \alpha =0.\end{equation}
One can easily  see that
\begin{equation}\label{4.13}   -\frac{1}{2i\pi}\int_{-\infty}^{\infty}
\frac{1+tz}{(z-t)(1+t^{2})}\log(2\pi)dt=\frac{1}{2}\log(2\pi). \end{equation}
It follows from \eqref{4.10}, \eqref{4.12}, and \eqref{4.13} that $\Pi(z)$ has the form
\begin{equation}\label{4.14}   \Pi(z)=
\prod_{k=0}^{n-1}\left[\left(\frac{a_{k+1}+z}{a_{k+1}-z}\right)\left(\frac{a_{k}-z}{a_{k}+z}\right)\right]^{\log(1-b_{k}\mu)/2i\pi},\quad {\Im}z>0.
\end{equation}
Next, we prove  the main result of this paper.
 \begin{thm}\label{Theorem 4.1.}
 { The bounded  positive definite and invertible operator $S$, which is defined by formulas \eqref{2.1} and \eqref{4.5},
does not admit a left triangular factorization.}
\end{thm}
\begin{proof}
Taking into account  Lemma 3.3 and relations  \eqref{2.37},  \eqref{4.8},  and \eqref{4.9} we have
 \begin{equation}\label{4.24}    \lim_{x_{n}{\to}\infty}w_{1,2}(x_{n},z)=-C\Pi(z),\quad {\Im}z>0,\quad C=1/\sqrt{(1-\mu)}.\end{equation}
Now, we use the following relations
\begin{equation}\label{4.25}   \lim_{y{\to}+0}\left(\frac{a_{k+1}-iy}{a_{k+1}+iy}\right)\left(\frac{a_{k}+iy}{a_{k}-iy}\right)=1,\quad k>0,
 \end{equation}
 \begin{equation}\label{4.26}   \lim_{y{\to}+0}\left(\frac{a_{k+1}-iy}{a_{k+1}+iy}\right)\left(\frac{a_{k}+iy}{a_{k}-iy}\right)=-1,\quad k=0.
 \end{equation}
 Formulas \eqref{4.14},  \eqref{4.25}, and \eqref{4.26} imply that
  \begin{equation}\label{4.27}   \lim_{y{\to}+0}\Pi(iy)=\sqrt{(1-\mu)}.\end{equation}

 Suppose that the operator $S$ admits a factorization.
It follows from  the asymptotics of sinus integral (see \cite{2}, Ch. 9, formulas (2) and (10)),
that the kernel $h(x)$, defined by formula \eqref{4.5}, satisfies conditions \eqref{3.2}.
Hence, the conditions of Lemma 3.4 are fulfilled.
 Comparing formulas \eqref{3.18} and \eqref{4.24},  we see that
\begin{equation}\label{4.29}   -\lim_{y{\to}+0}\overline{G(-iy)}=-1/(1-\mu){\ne}-C\lim_{y{\to}+0}\Pi(iy)=-1.\end{equation}
Hence, the relation
 \eqref{3.18} is not true. According to Proposition 3.1 the operator $S$ does not admit a factorization.
 The theorem is proved.
 \end{proof}

\section{Examples instead of existence theorems} \label{sec5}

Let the nest $N$ be the family of subspaces $Q_{\xi}L^{2}(0,\infty).$
The corresponding \emph{nest algebra} $Alg(N)$ is the algebra of all linear bounded operators
 in the space $L^{2}(0,\infty)$ for which every subspace from $N$ is an invariant subspace. Put $D_{N}=Alg(N){\bigcap}Alg(N)^{\ast}$.
 The set $N$ has \emph{ multiplicity one} if the diagonal $D_{N}$ is abelian, that is,
 $D_{N}$ is a commutative algebra.
 We can see that the lower triangular operators $S_{-}$ form the algebra
 $Alg(N)$, the corresponding diagonal $D_{N}$ is abelian, and it consists
 of the commutative  operators
 \begin{equation}\label{5.1}   T_{\varphi}f=\varphi(x)f,\quad f{\in}L^{2}(0,\infty),\end{equation}
 where $\varphi(x)$ is bounded. Hence, the introduced nest $N$ has the multiplicity $1$.\\

\noindent {\textbf{Ringrose Problem.}}
 \textit{Let} $N$ \textit{be a multiplicity one nest and} $T$ \textit{be a bounded invertible operator. Is} $TN$ \textit{necessarily multiplicity one nest?}\\

We obtain a concrete counterexample to Ringrose's hypothesis.

  \begin{prop}\label{Proposition 5.1.}
{Let the positive definite,
invertible operator $S$ be defined by the relations \eqref{2.1} and \eqref{4.5}.
 The  set $S^{1/2}N$ fails to have  multiplicity 1.}
 \end{prop}

 \begin{proof}
We use the well-known result (see \cite{3}, p. 169):\\
\emph{The following assertions are equivalent:}\\
\emph{1. The  positive definite, invertible operator T admits factorization.\\
2. $T^{1/2}$ preserves  the multiplicity of $N$.}\\
We stress that in our case the set $N=Q_{\xi}L^{2}(0,\infty)$ is fixed.)
The operator $S$  does not admit the factorization. Therefore,
the  set $S^{1/2}N$ fails to have  multiplicity $1$. The proposition is proved.
\end{proof}
Next, consider the operator
\begin{equation}\label{5.2}    Vf=\int_{0}^{x}e^{-(x+y)}f(y)dy, \quad f(x){\in}L^{2}(0,\infty).\end{equation}
An operator is said to be \emph{hyperintransitive} if its lattice of invariant
subspaces contains a multiplicity one nest. Note that the lattice of invariant
subspaces of the operator  $V$ coincides with $N$, see \cite{9} and \cite{22} (Ch. 11, Theorem 150). Hence we  deduce the answer
to Kadison-Singer \cite{5} and to Gohberg-Krein \cite{4} question.

\begin{cor} \label{Corollary 5.1.}
{The operator $W=S^{1/2}VS^{-1/2}$  is a non-hyperintransitive
compact operator.}
\end{cor}

Indeed, the lattice of the invariant subspaces of
 the operator $W$ coincides with $S^{1/2}N.$

\begin{rem}\label{Remark 5.1.}
The existence parts of Theorem \ref{Theorem 4.1.}, Proposition \ref{Proposition 5.1.}, and
Corollary \ref{Corollary 5.1.}
 are proved by
D.R. Larson \cite{8}.
\end{rem}