\documentstyle[11pt]{article}

\begin{document}

\def\bff{{\bf f}}
\def\bfk{{\bf k}}
\def\bfx{{\bf x}}
\def\bfy{{\bf y}}
\def\bfr{{\bf r}}
\def\bfu{{\bf u}}
\def\bfv{{\bf v}}
\def\bfw{{\bf w}}
\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\bfV{{\bf V}}
\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\half{{1\over 2}}

\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}}
\def\eqref#1{\ref{eq:#1}}


\title{Fast Operations Involving the Matrix $|t_i-t_j|$}
\author{George B. Rybicki}
\date{4 March 1998}

\maketitle

\section{Theory}

This note is
concerned with finding the inverse of the symmetric matrix ${\bf
\Psi}$, with components
\bel{1}
        \Psi_{ij}=|t_i-t_j|
\ee
where the values $t_i$, $i=1,2,3,\ldots,n$, are all distinct.
This problem arises in the statistical treatment of time series,
related to the ``random walk'' process. 

We assume that the sequence $t_i$ is given in strict ascending order,
$t_1 < t_2 < \ldots t_n$ (if not true initially, this can be arranged 
by simple relabelling of the indices).  Using the monotonicity
of the $t_i$, the matrix can be written,
\def\matel#1#2{t_{#1}-t_{#2}}
\bel{2}
\Psi =\pmatrix{
    0 & \matel{2}{1} & \matel{3}{1} & \ldots & \matel{n}{1} \cr
    \matel{2}{1} & 0 & \matel{3}{2} & \ldots & \matel{n}{2} \cr
    \matel{3}{1} & \matel{3}{2} & 0 & \ldots & \matel{n}{3} \cr
    \vdots &\vdots &\vdots &\ddots &\vdots \cr
    \matel{n}{1} &\matel{n}{2} & \matel{n}{3} &\ldots & 0\cr}
\ee
It can be shown that the inverse of $\Psi$ may be expressed in the form
\bel{3}
     \Gamma = \Psi\inv = \bfT + \mu \bfL \bfL\tr
\ee
where $\bfT$ is a symmetric tridiagonal matrix, and $\bfL$ is a column
vector consisting of zeroes, except for the first and last elements,
which are unity.  In particular,
\bel{4}
      \bfT =\pmatrix{ 
    d_1 & e_1 & 0 & \ldots &0 & 0 \cr
    e_1 &d_2 & e_2 & \ldots &0 &0 \cr
    0& e_2 &d_3 & \ldots &0 &0 \cr
    \vdots &\vdots &\vdots &\ddots &\vdots &\vdots\cr
    0 &0 &0 &\ldots &d_{n-1} & e_{n-1} \cr
    0 &0 &0 &\ldots & e_{n-1} &d_n\cr}, \qquad
        \bfL =\pmatrix{ 1 \cr
                        0 \cr
                        0 \cr
			\vdots \cr
			0 \cr
			1 \cr}.
\ee
where $e_i$ and $d_i$ are given by 
\bel{5}
e_i = \half(t_{i+1}-t_i)\inv, \qquad 1 \le i \le n-1,
\ee
and
\beal{6}
d_1 &=& -e_1 \\
d_i &=& -e_{i} -e_{i-1}, \qquad 2 \le i \le n-1, \\
d_n &=& -e_n,
\eea
and where $\mu $ is given by
\bel{7}
\mu  = \half(t_n-t_1)\inv.
\ee

The inverse can also be written in the form
\bel{8}
\Gamma=\Psi\inv =\pmatrix{ 
    d_1+\mu  & e_1 & 0 & \ldots &0 & \mu  \cr
    e_1 &d_2 & e_2 & \ldots &0 &0 \cr
    0& e_2 &d_3 & \ldots &0 &0 \cr
    \vdots &\vdots &\vdots &\ddots &\vdots &\vdots\cr
    0 &0 &0 &\ldots &d_{n-1} & e_{n-1} \cr
    \mu  &0 &0 &\ldots & e_{n-1} &d_n+\mu \cr}
\ee
which is seen to be a tridiagonal matrix, modified by the addition of
$\mu $ to each of the four ``corner'' elements.

We remark that these results can be trivially generalized to the case
of the scaled matrix $S_{ij}=\alpha|t_i-t_j|$, where $\alpha$ is a
constant.  All elements of the inverse matrix (i.e., the $d_i$, $e_i$,
and $\mu$) are then simply scaled by the factor $\alpha\inv$.

The above results, especially the representation (\eqref{3}) for the
inverse of $\Psi$, permit a number of ``fast,'' $O(n)$
matrix operations involving $\Psi$.  We present three such operations:

\begin{enumerate}

\item There is a fast solution of the linear system
\bel{9}
     \bfy= \Psi \bfx
\ee
for $\bfx$ given $\bfy$, namely,
\bel{10}
      \bfx= (\bfT +\mu \bfL \bfL\tr)\bfy = \bfT \bfy + \mu (\bfL\tr \bfy) \bfL
\ee
which involves only $O(n)$ operations.

\item There is a fast evaluation of the direct product
$\bfy=\Psi \bfx$ itself, a kind of fast ``convolution.''  In this case we
write
\bel{11}
      \bfy =  (\bfT +\mu \bfL \bfL\tr)\inv \bfx
\ee
Using the Lemma in Appendix A, we have
\bel{12}
      \bfy=\bfu - {\mu \bfL\tr \bfu \over 1+\mu \bfL\tr \bfv}\bfv
\ee
where
\bel{13}
       \bfu=\bfT\inv \bfx, \qquad\qquad \bfv=\bfT\inv \bfL
\ee

\item The following matrix evaluation is often required in practice,
\bel{14}
    \bfw= (\Psi +\bfN)\inv \bfy
\ee
where $\bfN$ is a diagonal matrix and $\bfy$ is a column vector.
This can also be written
\bel{15}
     \bfw= \bfN\inv (\bfN\inv + \Psi\inv )\inv \Psi\inv \bfy
\ee
Noting that $(\bfN\inv +\Psi\inv)\inv (\bfN\inv +\Psi\inv) 
= (\bfN\inv +\Psi\inv)\inv \bfN\inv +(\bfN\inv +\Psi\inv)\inv \Psi\inv = 1$, we
have
\bel{16}
   \bfw= \bfN\inv\left[ 1-(\bfN\inv +\Psi\inv)\inv \bfN\inv \right] \bfy
\ee
or
\bel{17}
   \bfw= \bfN\inv \bfy 
       -\bfN\inv ( \bfN\inv +\bfT +\mu \bfL\bfL\tr)\inv \bfN\inv \bfy
\ee
Since $\bfN\inv +\bfT$ is tridiagonal, this can be evaluated using the
fast method of the Lemma in the Appendix.

\end{enumerate}

\appendix

\section{A Useful Lemma}

\newtheorem{Lemma}{Lemma}    
\begin{Lemma} 
 Let $\bfT$ be a tridiagonal matrix and let $\bfU$, $\bfV$, and $\bfx$
be column vectors.  Then
\bel{A1}
     (\bfT +\bfU\bfV\tr)\inv \bfx 
     = \bfA - {\bfV\tr \bfA \over a+\bfV\tr \bfB} \bfB,
     \quad {\rm where} \quad 
     \bfA = \bfT\inv \bfx, \quad \bfB = \bfT\inv \bfU
\ee
which can be evaluated in $O(n)$ operations.
\end{Lemma}

\noindent
{\bf Proof } The Sherman-Morrison formula is
\bel{A2}
    (\bfT +\bfU\bfV\tr)\inv 
    = \bfT\inv - \bfT\inv \bfU (1+ \bfV\tr \bfT\inv \bfU)\inv \bfV\tr \bfT\inv
\ee
Application of this immediately gives the stated result (\eqref{A1}).
Since this formula involves only two tridiagonal solutions
and two scalar products, this is a ``fast'' process,
involving only $O(n)$ operations. $\bullet$



\end{document}
