A square root algorithm is presented in the following,
where the  covariance update is  computed using matrix
factorization and triangularization \cite{Kailath-encyc}.

Consider the matrix
  \[
	S = \left [ \begin{array}{cc}
			R_{e,t}	&  HP_{t} \Phi ^T \\
			\Phi P_{t}H^{T}	& \Phi P_{t} \Phi ^T+GQG^T
		    \end{array}
	    \right ],
  \]
where some of the subscripts are not explicitly shown.
Define the following factorizations:
  \begin{eqnarray*}
	P_{t}	 &=& L \overline{D} L^{T}			\\
	P_{t+1}	 &=& L_{+}D_{+}L_{+}^{T}	\\
	R_{e,t} &=& L_{e}D_{e}L_{e}^{T}.
  \end{eqnarray*}
 The matrix $S$ can be factored as
  \begin{eqnarray*}
    S &=& ADA^{T} \\
      &=&
	\left [ \begin{array}{ccc}
		I & HL & 0		\\
		0 & \Phi L & G
		\end{array}
	\right ] \! \! \! \!
	\left [ \begin{array}{ccc}
		R & 0 & 0	\\
		0 & \overline{D} & 0	\\
		0 & 0 & Q
		\end{array}
	\right ] \! \! \! \!
	\left [ \begin{array}{cc}
		    I		&       0		\\
		L^{T}H^{T}	& L^{T} \Phi ^{T}	\\
		    0		&	G^{T}		
		\end{array}
	\right ],
  \end{eqnarray*}
where $A$ is an ($n$+$m$) x (2$n$+$m$) matrix
and $D$ is a (2$n$+$m$) x (2$n$+$m$) diagonal matrix.
\footnote{For this paper, matrices $R$ and $Q$ are assumed to be diagonal.
If these matrices are not diagonal, the factorization can be changed
so that $D$ remains diagonal.}
Similarly, it is easily verified that $S$ also can be factored as

  \begin{eqnarray*}
    S &=& A'D'(A')^{T} \\
      &=&
	\left [ \begin{array}{ccc}
		L_{e}		&   0	& 0	\\
		K_tL_{e}		& L_{+} & 0
		\end{array}
	\right ] \! \! \! \!
	\left [ \begin{array}{ccc}
		D_{e}	& 0 	& 0	\\
		0 	& D_{+} & 0	\\
		0 	& 0 	& D_{a}
		\end{array}
	\right ] \! \! \! \!
	\left [ \begin{array}{cc}
		L_{e}^{T}	& L_{e}^{T}K_t^{T}	\\
		     0		& 	L_{+}^{T}	\\
		     0		&	    0	
		\end{array}
	\right ]
  \end{eqnarray*}
where $A'$ is an ($n$+$m$) x (2$n$+$m$) unit lower triangular matrix
and $D'$ is a (2$n$+$m$) x (2$n$+$m$) diagonal matrix.
Because of the nature of this factorization, the matrix $D_{a}$ is
arbitrary.

From these two block factorizations of $S$, it is known that there exists 
an orthogonal transformation matrix $\Theta$ such that
  \begin{eqnarray}
  \lefteqn{
	\left [ \begin{array}{ccc}
		I & HL & 0		\\
		0 & \Phi L & G
		\end{array}
	\right ]
	\left [ \begin{array}{ccc}
		R & 0 & 0	\\
		0 & \overline{D} & 0	\\
		0 & 0 & Q
		\end{array}
	\right ]^{\frac{1}{2}} \Theta =
	  }					\nonumber \\
 & & \hspace*{0.5in}
	\left [ \begin{array}{ccc}
		L_{e}		&   0	& 0	\\
		K_tL_{e}		& L_{+} & 0
		\end{array}
	\right ]
	\left [ \begin{array}{ccc}
		D_{e}	& 0 	& 0	\\
		0 	& D_{+} & 0	\\
		0 	& 0 	& D_{a}
		\end{array}
	\right ]^{\frac{1}{2}}
     \label{pre-ADThetaA'D'}
  \end{eqnarray}
or, for a suitably defined $\Psi$,
  \begin{eqnarray}
	\left [ \begin{array}{ccc}
		I & HL & 0		\\
		0 & \Phi L & G
		\end{array}
	\right ]  \Psi &=& 
	\left [ \begin{array}{ccc}
		L_{e}		&   0	& 0	\\
		K_tL_{e}		& L_{+} & 0
		\end{array}
	\right ].
     \label{pre-APsiA'}
  \end{eqnarray}
It is seen that the right side of Equation~\ref{pre-APsiA'}
has a triangular structure.  
Using fast Givens rotations \cite{Gentleman73}, we can find $\Psi$
as a product of elementary matrices, none of which contain any scalar
square roots.  The fast Givens technique is discussed in the next
section.

Since matrix $A'$ contains $L_{+}$ and matrix $D'$ contains $D_{+}$, the 
matrix $P_{t+1}$ can be explicitly computed as soon as the triangularization
procedure completes.  However, since the recursive algorithm requires the 
factorized version of $P_{t}$, explicit calculation of the
error covariance matrices is unnecessary.

Since matrix $A'$ contains both $L_{e}$ and
$K_t L_{e}$, the gain matrix can be computed by 
  \[ K_t = [K_t L_{e}][L_{e}]^{-1}. \]
Note that this involves computing the inverse of  
unit lower triangular $m$ x $m$ matrix $L_e$.


