13.5. History of the Eigenvalue Problem

John G.F. Francis developed two related algorithms in London, England, between 1959 and 1961 that finally solved the eigenvalue problem. Francis found success that eluded several mathematicians since the eighteenth century. Francis’s insights were certainly brilliant. However, Francis can not claim full credit for his fantastic achievement. Issai Schur established that cleverly applying unitary transformation matrices could solve the eigenvalue problem (Schur Triangularization Theorem) fifty years before Francis’s efforts. Then, in the years just before Francis’s efforts, Wallace Givens and Alston Householder provided the methods for creating the needed unitary transformation matrices (Unitary Transformation Matrices).

13.5.1. Early Explorations

Eigenvalues and eigenvectors were first discovered in the eighteenth century by Leonhard Euler and Joseph-Louis Lagrange while studying problems related to the rotational motion of rigid bodies. Several mathematicians discovered properties of eigensystems in the eighteenth and nineteenth centuries. In particular, two nineteenth-century mathematicians established important facts that provided clues relative to the eigenvalue problem. French mathematician Augustin-Louis Cauchy determined in 1827 that symmetric matrices would be the simplest type of matrix to diagonalize [STEWART93]. He also discovered that similarity transformations allow changes to matrices while preserving the eigenvalues. German mathematician Carl Gustav Jacob Jacobi proposed an iterative method in 1846 for calculating the eigenvalues and eigenvectors of real symmetric matrices [STRANG16]. His method used similarity transformations and rotation matrices. Jacobi was also the first to note that the diagonal values of a triangular matrix are the matrix’s eigenvalues.

13.5.2. Schur’s Insight

Despite the early discoveries, the eigenvalue problem was an unsolved conundrum at the beginning of the twentieth century. It was unclear then that the solution would not come from factoring the characteristic polynomial, rather than operations on the matrix. Issai Schur hit on the strategy that led to the ultimate solution in 1909. He published a paper that established the existence of a matrix decomposition satisfying the similarity requirement with factors of a unitary matrix, its transpose, and a triangular matrix [SCHUR09]. He developed the triangularization theorem (Schur Triangularization Theorem in Schur Decomposition) before suitable unitary transformation matrices or algorithms were discovered that could bring about the proposed factorization, the Schur decomposition. The advantages of unitary matrices were recently discovered in 1902 [AUTONNE02]. A review of linear algebra literature from the beginning of the twentieth century until the 1970s shows that systems of linear equations and the eigenvalue problem were paramount to the researchers. Many researchers proposed algorithms to solve the eigenvalue problem in the fifty years following Schur’s seminal paper, but none were successful.

Note

Issai Schur was the last Jewish professor employed at the University of Berlin during the late 1930s. Although not directly killed by the events of the Holocaust, he was greatly persecuted. He died in 1941 on his sixty-sixth birthday [OCONNOR98].

13.5.3. Anticipation of Unitary Matrices

Chapter VIII of Herbert Turnbull and Alexander Aitken’s 1932 text [TURNBULL32] discusses the utility of unitary transformation matrices. Section 2 covers the topic of vector rotations from multiplication by a unitary matrix. They prove the existence of a unitary matrix that matches what we now call a Householder reflection matrix. They begin with a lemma and its proof for using a unitary transformation to rotate a vector to another vector of the same length.

Lemma

Any two nonzero vectors of order \(n\) having equal norms may be transformed into each other by unitary transformations. For vectors \(\bm{p}\) and \(\bm{q}\) if \(\norm{\bm{p}} = \norm{\bm{q}}\), then there exists a unitary matrix \(\bf{R}\) such that \(\bm{p} = \mathbf{R}\bm{q}\).

They extend the lemma to establish the existence of the Householder reflection transform.

Corollary

A unitary matrix \(\bf{R}\) exists, which transforms a given vector of norm \(\alpha^H \alpha\) into {\(\alpha, 0, 0, \ldots, 0\)}.

We recognize \(\bf{R}\) from the corollary as an early form of Householder reflection because the resulting vector contains only zero elements in positions other than the first. Turnbull and Aitken project an application of the transform to the Schur decomposition. The form of the Schur decomposition [SCHUR09] was already established when Turnbull and Aitken wrote, although it was not yet so named. They give a theorem that restates the Schur triangularization theorem and proves that a transform such as what Householder later developed meets the unitary requirements of Schur Triangularization Theorem. Although Turnbull and Aitken establish the existence of the modern Householder reflection and its utility for finding the Schur decomposition, they do not show algorithms for finding either. That must wait for future generations of mathematicians and for computers to run the algorithms.

13.5.4. Oak Ridge and Unitary Matrices

Wallace Givens and Alston Householder developed the Givens rotation transform and the Householder reflection transform during the 1950s at Oak Ridge National Laboratory [1] (ORNL) in Oak Ridge, Tennessee. The set of technologies that researchers at ORNL study includes mathematics and computing. Alston Householder became the head of a new department called the Mathematics Panel in 1948. Householder attended the Mark I unveiling at Harvard and became convinced that Oak Ridge needed a large computing facility. After receiving bids from several vendors, ORNL purchased a copy of the Institute for Advanced Study (IAS) machine that John von Neumann designed. They called their computer ORACLE, an acronym for Oak Ridge Automatic Computer and Logical Engine. It was built at Argonne National Laboratory and installed at ORNL in February 1954 [MERTZ70].

What ORNL had in ORACLE, and the IBM machines that later replaced it, was a machine that required a staff of engineers to maintain and operate it. The facilities for input/output and programming were limited. The Semiannual Progress Reports of the Mathematics Panel [2] for the 1950s and 1960s document the progression of developing a library of essential programming functions, including mathematical functions. They wrote ORACLE’s early programs in assembly language. The programming capabilities improved considerably after they developed a compiler for the ALGOL programming language.

The Mathematics Panel worked on and wrote computer programs to solve problems in several branches of mathematics. Householder and a few other researchers, including Wallace Givens, had a strong interest in linear algebra and numerical analysis. Their highest priorities included unitary matrices, algorithms for orthogonal matrix factoring, and the eigenvalue problem.

13.5.4.1. Origins of Givens Rotation Transformations

Since geometric rotation matrices were known to be orthogonal, the strategy of using multiplication with a unitary matrix to change a matrix element to zero using rotation about another matrix element, as described in Givens Rotation Matrices, appears to have been a fairly obvious starting point for Alston Householder and Wallace Givens. Jacobi successfully used similar rotations more than one hundred years earlier to find the eigenvalues of the simplest matrix type (real and symmetric). Householder reports in his 1953 text [HOUSEHOLDER53] that Givens wrote about the construction of the rotation transformations in an ORNL internal memorandum in 1951. It is unknown to what extent Householder was involved in Givens’ work. He would have certainly been an advisor since Givens then reported to Householder.

Beyond just developing the unitary transformation, Givens developed programs to apply it to numerical analysis problems. He authored an ORNL Technical Report in March of 1954 on using unitary rotations to calculate the eigenvalues of symmetric matrices [GIVENS54]. Then, he reported in the December 1954 Semiannual Progress Report for the Mathematics Panel that he could also find the eigenvectors of symmetric matrices [ORNL54]. He reported using unitary rotations to find the inverse of a matrix in the Semiannual Progress Report of June 1956. However, he also reported only partial success in finding the eigenvalues of non-symmetric matrices [ORNL56].

After leaving ORNL for a faculty position at Wayne State University, Givens published a paper describing difficulties associated with the eigenvalue problem for non-symmetric matrices [GIVENS57]. This report does not go into the same level of detail about the algorithms he used as the technical report he wrote at ORNL. He presented a paper at a conference in late 1957 that was later published as a full JACM journal article in 1958 giving more details about unitary rotations and their application to the eigenvalue problem [GIVENS58]. The JACM article is the most frequently cited reference for Givens rotation transformations.

13.5.4.2. Origins of Householder Reflection Transformations

Owing likely to the humility of Dr. Householder, reflection transformations (Householder Reflection Matrices) were introduced with minimal publicity. The paper in which they were introduced [HOUSEHOLDER58] is only slightly over three pages long, and the only suggested application is for QR factoring to solve systems of equations or to calculate the inverse of a matrix. His stated reason for developing the unnamed transform is to provide an alternative to Givens rotators for some applications. He points out that \(n(n-1)/2\) rotators, each with a computationally expensive square root, are required for a \(n{\times}n\) QR factorization. The Householder transform only needs \(n - 1\) transforms to accomplish the same task. Householder transforms are also faster on integer processors because they do not require square root calculations. He also briefly described the transformation in his 1964 book on numerical analysis [HOUSEHOLDER64].

It is now widely accepted that the Householder reflector is perhaps the most useful tool available in numerical analysis [DUBRULLE00]. The Givens rotator certainly has applications, but the Householder transform provides the foundation for orthogonal matrix factoring. Householder mentioned applying it to the QR decomposition, but it now has a role in every orthogonal matrix factorization.

Householder does not refer to the transform as a reflection matrix. Others gave it that name. The general equation for a reflection was known, but Householder adjusted the application. A reflection problem is often viewed from the perspective of finding the reflection of a vector relative to a fixed hyperplane. Imagine the hyperplane as a reflective surface like a mirror [HIGHAM02]. In the Householder transform equation, we know what the vector is and the target where it should reflect. The location of the reflective hyperplane is found in the form of its normal vector (\(\bm{u}\) in equation (12.2)). Then the reflector matrix is found from \(\bm{u}\).

13.5.5. The LR Algorithm

During the 1950s, while Alston Householder and his colleagues in the Mathematics Panel 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].

Rutishauser developed what was called the qd algorithm [3] 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 \(\m{L}\) and \(\m{R}\), and proposed the LR algorithm as a method to find the eigenvalues of a matrix or the roots of a polynomial. Rutishauser’s LR algorithm was the inspirational strategy that led John Francis to the QR algorithm. He proposed the LR algorithm based on the following similarity transform.

\[\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}\]

Today, we would call his LR factoring the LU decomposition, which is solved with an elimination procedure (LU Decomposition). LU decomposition requires pivoting for numerical stability, which would complicate a general eigenvalue finding algorithm. Still, it provides a way to find the eigenvalues of some matrices. Consider the following example: a close approximation to the eigenvalues of a \(3{\times}3\) symmetric matrix is found in MATLAB’s Command Window with ten iterations of the similarity transformation using factors from LU decomposition.

>> 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

13.5.6. John Francis

As mentioned, John G.F. Francis solved the eigenvalue problem between 1959 and 1961. Working for the National Research Development Corporation (NRDC), Francis was assisting the British Ministry of Aircraft Production with stability analysis of airplane designs. He was given matrices of size \(27{\times}27\) representing models of aircraft flutter for which he needed to find the eigenvalues and eigenvectors [UHLIG09]. He met James H. Wilkinson, a linear algebra expert from the National Physical Laboratory in London [UHLIG09], [IEEE13]. He decided to try a strategy that combined the latest developments from three sources—Wilkinson, Heinz Rutishauser from the Institute of Applied Mathematics at ETH Zurich, and colleagues Wallace Givens and Alston Householder from Oak Ridge National Laboratory in Tennessee. He wrote a program in assembly language that implemented an algorithm called the QR algorithm (Iterative QR Algorithm). Later, he improved his first algorithm, creating what he called the QR algorithm with implicit shifts [4] (Iterative Francis Algorithm). He published two papers [FRANCIS61], [FRANCIS62] about his algorithms.

Certainly, Francis’s success can be partially attributed to working on the problem at the right time [UHLIG09]. Heinz Rutishauser’s LR algorithm (The LR Algorithm) was close to solving it [GUTKNECHT10]. The primary change that Francis made to the LR algorithm was to follow the advice of Wilkinson and use the unitary transformation matrices of Givens and Householder [GIVENS58], [HOUSEHOLDER58], [WILKINSON59], [FRANCIS61] instead of an algorithm using elimination procedures. Francis credits Rutishauser in his first publication, stating that his QR algorithm is similar to the LR algorithm except that the transformations are all unitary, which should improve the numerical stability [FRANCIS61]. The switch to unitary transformations also allowed shifts and deflation, as suggested by Wilkinson [WILKINSON59], [WILKINSON68shift] to be more effective [WILKINSON65]. Indeed, the accumulated knowledge from all of the past efforts made the time right for Francis.

13.5.7. Vera Kublanovskaya

Vera Nikolaevna Kublanovskaya independently developed the QR algorithm to solve the eigenvalue problem. She is usually not given the same level of recognition for doing so as John Francis. There are a few reasons for this.

1. She published her algorithm after Francis did. There is no doubt that her research was independent of Francis’s. Like Francis, she cited Rutishauser’s LR algorithm (The LR Algorithm) as similar to her algorithm, except that she uses unitary transformations. Nevertheless, more credit is usually given to the first one to publish.

2. She had no computer access to write a program to run her algorithm. She verified her algorithm with hand-written manual calculations.

3. Her algorithm was the same as Francis’s first QR algorithm. Francis’s second algorithm, QR with implicit shifts (now called the Francis algorithm), is the better algorithm.

4. Being in the Soviet Union, her publications were in Russian and may have been more obscure to researchers in other countries.

Footnotes:

[1]

After the conclusion of the Manhattan Project and World War II, the laboratory at Oak Ridge took on the task of developing nuclear reactor technology for the U.S. Department of Energy in 1943. The laboratory still operates as a research institution. As their website (https://www.ornl.gov/) states, “Oak Ridge National Laboratory delivers scientific discoveries and technical breakthroughs needed to realize solutions in energy and national security and provide economic benefit to the nation.”

[2]

Scanned copies of these reports are available on the website of the Department of Energy’s Office of Scientific and Technical Information (https://www.osti.gov/).

[3]

Although not capitalized, qd is an abbreviation for quotient-difference.

[4]

The implicit algorithm is now also called the Francis algorithm.