.. _LR:

The LR Algorithm
----------------

During the 1950s, when Alston Householder and his colleagues in the
Mathematics Division at Oak Ridge National Laboratory were first
learning how to operate, maintain, and program a computer, Eduard
Stiefel and his staff were going through the same learning experience at
the Institute of Applied Mathematics at ETH Zurich, Switzerland. Stiefel
hired Heinz Rutishauser to research numerical methods that could be
programmed to run on a computer [GUTKNECHT10]_.

.. index:: qd algorithm

Rutishauser developed what was called the
*qd algorithm* [1]_ to find the roots of a particular type of
polynomial. He observed that the work of the qd algorithm was the same
as factoring a tridiagonal matrix into a product of lower and upper
triangular matrices and then forming the product in reverse order. He
called the two triangular matrices :math:`\m{L}` and :math:`\m{R}`.
Thus, he proposed the LR algorithm based on the following similarity
transform.

.. math::

   \m{A} = \m{L\,R}, \quad \m{R} = \m{L}^{-1}\m{A}, 
       \quad \m{A}_1 = \m{R\,L} = \m{L}^{-1}\,\m{A\,L}

.. index:: LR algorithm

The LR algorithm was proposed as a method to find the eigenvalues of a matrix
or the roots of a polynomial. Today, we would call his LR factoring the LU
decomposition, which is solved with an elimination procedure (:ref:`LUdecomp`).
For numerical stability, LU decomposition requires pivoting, which complicates
the algorithm.

Rutishauser’s LR algorithm was the inspirational strategy that led John
Francis to the QR algorithm. Although LU decomposition uses elimination
instead of orthogonal factoring methods, it still provides a way to find
the eigenvalues of some matrices. Consider the following example where a
close approximation to the eigenvalues of a :math:`3{\times}3` positive
definite matrix is found in MATLAB’s Command Window with ten iterations
of the similarity transformation.

::

   >> A = 5*randn(3);
   >> A = A*A';
   >> eig(A)
   ans =
       0.1294
       9.9568
      44.5828 
   >> for i=1:10
   [L,U] = lu(A);
   A = U*L;
   end
   >> sort(diag(A))
   ans =
       0.1294
       9.9567
      44.5829

.. [1]
   Abbreviation qd stands for quotient-difference.