\documentstyle[12pt]{article} \begin{document} \def\bff{{\bf f}} \def\bfk{{\bf k}} \def\bfx{{\bf x}} \def\bfy{{\bf y}} \def\bfr{{\bf r}} \def\bfA{{\bf A}} \def\bfB{{\bf B}} \def\bfC{{\bf C}} \def\bfD{{\bf D}} \def\bfE{{\bf E}} \def\bfI{{\bf I}} \def\bfGamma{{\bf \Gamma}} \def\bfK{{\bf K}} \def\bfL{{\bf L}} \def\bfN{{\bf N}} \def\bfPhi{{\bf \Phi}} \def\bfT{{\bf T}} \def\bfU{{\bf U}} \def\bfone{{\bf 1}} \def\bfxi{{\bf \xi}} \def\tb#1{{\tilde {\bf #1}}} \def\tbf{\tb{f}} \def\tbk{\tb{k}} \def\tbA{\tb{A}} \def\tbB{\tb{B}} \def\tbC{\tb{C}} \def\tbE{\tb{E}} \def\tbI{\tb{I}} \def\xbar{{\bar x}} \def\ybar{{\bar y}} \def\Nbar{{\bar N}} \def\nbar{{\bar n}} \def\chibar{{\bar \chi}} \def\bfxbar{{\bar \bfx}} \def\bfybar{{\bar \bfy}} \def\bfNbar{{\bar \bfN}} \def\bfPhibar{{\bar \bfPhi}} \def\tr{^T} \def\inv{^{-1}} \def\tchi{{\tilde \chi}} \def\chiprh{\chi_{\rm PRH}} \def\tQ{{\tilde Q}} \def\Ahat{{\hat {\bf A}}} \def\be{\begin{equation}} \def\bel#1{\begin{equation}\label{eq:#1}} \def\ee{\end{equation}} \def\bea{\begin{eqnarray}} \def\beal#1{\begin{eqnarray}\label{eq:#1}} \def\eea{\end{eqnarray}} \title{Reduction Method for Duplicated Times} \author{George B. Rybicki} \date{3 May 1995} \maketitle Consider the linear equation \bel{1} (\bfPhi + \bfN)\bfx = \bfy, \ee which plays an important role in the fast methods proposed by Rybicki and Press (1995; hereafter RH). Here $\bfx$ and $\bfy$ are $n$-vectors,and $\bfN$ is the $n\times n$ diagonal ``noise'' matrix, with components \bel{2} N_{ij} = N_i \delta_{ij}, \ee and $\bfPhi$ is the $n\times n$ correlation matrix of the ``signal'', with components \bel{3} \Phi_{ij}=\Phi(t_i,t_j)=\langle s(t_i)s(t_j) \rangle. \ee For convenience, the times $t_i$, $i=1,2,\ldots,n$ are assumed to be given in (or sorted into) non-decreasing order, $t_1 \le t_2 \le \ldots \le t_n$. The object is to solve for the unknown $\bfx$ for a given $\bfy$ A problem arises in solving Eq.~(\ref{eq:1}) when some of the times $t_i$ are duplicated. Although the inverse $(\bfPhi+\bfN)\inv$ can still exist, providing the elements of $\bfN$ are nonzero, the inverse $\bfT=\bfPhi^{-1}$ does not exist, and the formal solution $\bfx=\bfT(\bfone+\bfN\bfT)^{-1}\bfy$, suggested in RH, cannot be used to effect a fast solution. The purpose of this note is to show how to solve Eq.~(\ref{eq:1}) by reducing it to one of the same form, but one of lower dimensionality in which all duplicate times are eliminated. This reduction works for matrices $\bfPhi$ of the general form (\ref{eq:3}), not just for the exponential matrices treated by RH. \section{The Reduction} There can be a number of sets of duplicated times, each with its own multiplicity. Let us concentrate on one such set consisting of $m$ duplicated times, say, $t_{k}=t_{k+1}=\ldots=t_{k+m-1}$. Then the $m$ rows $k,k+1,\ldots,k+m-1$ of $\bfPhi$ are identical, and so are the corresponding $m$ columns (this shows that $\bfPhi$ is singular). Eq.~(\ref{eq:1}) can now be substantially simplified by use of simple row and column operations. Subtracting the $k$th row from rows $k+1$ through $k+m-1$, all elements of $(\bfPhi+\bfN)$ in these rows are eliminated, except for the $k$th through $k+m-1$st columns. Other rows of the matrix are not affected. Rows $k+1$ through $k+m-1$ of the right hand side are changed from $y_i$ to $y_i-y_k$. A similar set of column operations can be effected by replacing the variable $x_k$ by the new variable \bel{5} \xbar_k= x_k + x_{k+1} + \ldots + x_{k+m-1}, \ee keeping all other variables the same. This eliminates all elements of the matrix in columns $k+1$ through $k+m-1$, except for rows $k$ through $k+m-1$. Other columns are not affected, nor is the right hand side. The effect of these operations is to partially decouple the system, so that the set of variables outside the range $k$ to $k+m-1$ couple only with themselves and with $\xbar_k$. The main coupling to be resolved is between the new variable $\xbar_k$ and the set of variables $x_{k+1}$ through $x_{k+m-1}$. The relevant set of equations consists of the system \bel{6} \matrix{ F_{k+1}x_{k+1} & + N_k x_{k+2} &+\ldots & + N_k x_{k+m-1} &= & N_k \xbar_k+y_{k+1}-y_k \cr N_k x_{k+1} & +F_{k+2} x_{k+2} &+\ldots & + N_k x_{k+m-1} &= & N_k \xbar_k+y_{k+2}-y_k \cr \vdots &\vdots &\ddots &\vdots & = & \vdots \cr N_k x_{k+1} & +N_k x_{k+2} &+\ldots & +F_{k+m-1}x_{k+m-1} &= & N_k \xbar_k+y_{k+m-1}-y_k, \cr} \ee where $F_l=N_k+N_l$, plus the single equation \bel{7} \ldots + (\Phi_{kk} +N_k)\xbar_k -N_k x_{k+1} \ldots -N_k x_{k+m-1} +\ldots = y_k. \ee The dots here refer to terms containing unknowns other than $\xbar_k$ and $x_{k+1}$ through $x_{k+m-1}$. The system (\ref{eq:5}) and Eq.~(\ref{eq:6}) can be solved for the variables $x_l$ in terms of $\xbar_k$ and the quantities $y_l$ and $N_l$ for $k\le l \le k+m-1$. Defining the quantities $\Nbar_k$ and $\ybar_k$ by \beal{8} \Nbar_k\inv &=& \sum_{l=k}^{k+m-1} N_l\inv \nonumber\\ \Nbar_k\inv\ybar_k &=& \sum_{l=k}^{k+m-1} N_l\inv y_l, \eea the result can be expressed in the simple form \bel{9} x_l=(\Nbar_k \xbar_k + y_l-\ybar_k)/N_l, \qquad k \le l \le k+m-1. \ee Using this equation, we can eliminate variables $x_{k+1}$ through $x_{k+m-1}$ from Eq.~(\ref{eq:7}), and we obtain \bel{10} \ldots + (\Phi_{kk} +\Nbar_k)\xbar_k +\ldots = \ybar_k. \ee We note that Eq.~(\ref{eq:10}) plus all the other equations in rows outside the range $k$ to $k+m-1$ form a system of $n-m+1$ equations of precisely the original form without the duplicated times, except for the replacements $x_k \rightarrow \xbar_k$, $y_k \rightarrow \ybar_k$, and $N_k \rightarrow \Nbar_k$. Once this system is solved, the original variables $x_k$ through $x_{k+m-1}$ can be determined by use of Eq.~(\ref{eq:9}). Since the altered system is of the same general form as the original, the same reduction can be applied to it to eliminate a second set of duplicate times, and a third, etc., until all duplicate times are eliminated. By relabeling the variables, it is possible to use the same method to solve the reduced system as is used to solve a system without duplicate times. Since the additional manipulations implied in Eqs.~(\ref{eq:8}) and (\ref{eq:9}) add at most only $O(n)$ operations, the reduction does not spoil any fast methods. \section{Further Developments} \subsection{Determinants} For some statistical applications it is necessary to compute the determinant of the matrix $\bfPhi+\bfN$. We show here how this can be done in terms of the reduced problem. After a little consideration of the reduced problem, it can be seen that this determinant differs from that of the reduced problem by the determinant of the $(m-1)\times (m-1)$ matrix $\bfK_k$ on the left hand side of Eq.~(\ref{eq:6}), that is, \bel{11} \det (\bfPhi+\bfN) = \det\bfK_k \det(\bfPhi+\bfN)_{\rm reduced}. \ee The matrix $\bfK_k$ can be written \bel{12} \bfK_k = \bfD_k + N_k \bfE \bfE\tr = (\bfone + N_k \bfE\bfE\tr\bfD_k\inv)\bfD_k, \ee where $\bfD_k$ is the $(m-1)\times (m-1)$ diagonal matrix $(D_k)_{ij}=N_{k+j-1}\delta_{ij}$ and $\bfE$ is the column vector with unity elements $E_i=1$. Using the general result $\det(\bfone+\bfA\bfB)=\det(\bfone+\bfB\bfA)$, valid for recangular matrices $\bfA$ and $\bfB$ for which $\bfA\bfB$ and $\bfB\bfA$ are square, we obtain the final result \beal{14} \det\bfK_k &=&\det(1+N_k \bfE\tr\bfD_k\inv\bfE) \det\bfD_k \nonumber\\ &=& \left(1+N_k\sum_{l=k+1}^{k+m-1}N_l\inv \right) N_{k+1}N_{k+2} \cdots N_{k+m-1} \nonumber\\ &=& \Nbar_k\inv \prod_{l=k}^{k+m-1} N_l. \eea When there are several sets of duplicate times, the correction to the determinant will be a product of determinants of the form (\ref{eq:14}), one for each set. In practice, it is more convenient to work with the logarithms of the determinant, since these are less likely to cause overflow. Thus we write Eq.~(\ref{eq:14}) in the equivalent form \bel{14.1} \log\det\bfK_k = -\log\Nbar_k +\sum_{l=k}^{k+m-1} \log N_l \ee \subsection{Evaluation of $\chi^2$} One of the important statistical uses of the fast methods is to calculate the quantity $\chi^2$, defined by \bel{15} \chi^2 = \bfy\tr (\bfPhi+\bfN)\inv \bfy. \ee By Eq.~(\ref{eq:1}) $\chi^2$ can be written simply \bel{16} \chi^2=\bfy\tr\bfx=\bfx\tr\bfy=\sum_{i=1}^{n} x_iy_i. \ee The evaluation can be done straightforwardly by solving the system (\ref{eq:1}) by the reduction method to get $\bfx$. However, it is of some interest to see how this $\chi^2$ relates to the reduced problem. Let us now adopt the notation that the reduced problem is \bel{17} (\bfPhibar + \bfNbar)\bfxbar = \bfybar, \ee which is of dimension $\nbar$. The reduced quantity analogous to $\chi^2$ is \bel{18} \chibar^2 = \sum_{i=1}^{\nbar} \xbar_i \ybar_i. \ee Let us assume that there is one set of duplicated times. Then $\chi^2$ and $\chibar^2$ are related by \bel{19} \chi^2=\chibar^2 + \sum_{l=k}^{k+m-1} x_l y_l. \ee Using Eq.~(\ref{eq:9}) to substitute for $x_l$ and with some manipulation, we obtain \bel{20} \chi^2=\chibar^2 + \chi^2_k \ee where \bel{21} \chi^2_k=\sum_{l=k}^{k+m-1} N_l\inv (y_l-\ybar_k)^2. \ee That is, the original $\chi^2$ is equal to the reduced $\chibar^2$ plus a ``local'' contribution $\chi^2_k$, due to the dispersion of the $y_l$ values about their local mean $\ybar_k$, inversely weighted by the appropriate mean square noise values $N_l$. In general, there will be a sum of local correction terms of this sort from each set of duplicated times. It is interesting to note that for Gaussian processes the expectation value of $\chibar^2$ is equal to the reduced number of variables, and it is the local contributions that bring the expectation value of $\chi^2$ up to the full number of variables $n$. \subsection{Reduced Data Sets} The analysis of large data sets (i.e., with large dimension $n$) can be computationally expensive. It is therefore of considerable importance to find ways to reduce the dimensionality of a data set, while retaining as much of its essential information as possible. In this subsection we consider how the preceding formulas can be used to construct reduced data sets. We restrict the discussion to cases where the data is to be used only for calculating $\chi^2$ and $\log\det(\bfPhi + \bfN)$. For example, these are the essential quantities in the method of Press, Rybicki and Hewitt (1992a,b) for determining the time delay of gravitational lenses. The reduction proposed here is based on the assumption that some of the times $t_i$ associated with the data, while not strictly at duplicate times, can be grouped into ``clusters'' having almost the same times. This might occur, for example, if astronomical observations of an object are taken over many years, but a number of the observations may have been taken within relatively short intervals, short in comparison to the correlation time of the signal. A reduced data set might be constructed which contains only one representative observation from each night (or maybe each hour; or maybe each week). The question then arises, how best to characterize this reduced, representative data. The formulas of the preceding subsections give the answer, namely, for each group of times the new data set consists of the reduced values $\ybar_k$ and $\Nbar_k$, defined by Eqs.~(\ref{eq:8}), as well as the associated local contribution $\chi^2_k$, defined by Eq.~(\ref{eq:21}) and the local contribution to the logarithm of the determinant, given by Eq.~(\ref{eq:14.1}). This constitutes the desired reduced data set. This reduced data set is very much like the old one in structure, except it now has the new quantities of local $\chi^2$ and log-determinant $\log\det \bfK_k$ attached to each point. We adopt the convention that the original set also had such local quantities associated with each point, call them $\chi^2_l$ and $\log\det\bfK_l$, but they were all zero. Then we can rewrite Eq.~(\ref{eq:21}) as \bel{22} \chi^2_k=\sum_{l=k}^{k+m-1} N_l\inv (y_l-\ybar_k)^2 +\sum_{l=k}^{k+m-1}\chi^2_l \ee and Eq.~(\ref{eq:14.1}) as \bel{23} \log\det\bfK_k = -\log\Nbar_k +\sum_{l=k}^{k+m-1} \log N_l +\sum_{l=k}^{k+m-1}\log\det\bfK_l \ee With this slight change, it can be shown that the reduction proceedure is ``associative,'' that is, if one goes through a series of reductions of the data, the result is independent of the intermediate stages and depends only on the final time groupings. It is easy to see that this is true for the quantities $\ybar_k$ and $\Nbar_k$, owing to the form of Eqs.~(\ref{eq:8}). It is also true for the formulas (\ref{eq:22}) and (\ref{eq:23}), but we shall not give the proof here. Note that the reduction achieved here is not necessarily in the storage requirement, since two extra quantities are now included per time. Of much greater importance is the reduction in the number of independent times $n$, since the computational times for the data analysis depend strongly on this. \end{document}