The linear least squares problem is where $A$ is an $m$ by $n$ matrix, $b$ is a given $m$ element vector and $x$ is an $n$-element solution vector.

In the most usual case $m\ge n$ and $\mathrm{rank}\left(A\right)=n$, so that $A$ has full rank and in this case the solution to problem (1) is unique; the problem is also referred to as finding a least squares solution to an overdetermined system of linear equations.

When $m<n$ and $\mathrm{rank}\left(A\right)=m$, there are an infinite number of solutions $x$ which exactly satisfy $b-Ax=0$. In this case it is often useful to find the unique solution $x$ which minimizes ${‖x‖}_2$, and the problem is referred to as finding a minimum norm solution to an underdetermined system of linear equations.

In the general case when we may have $\mathrm{rank}\left(A\right)<\min \left(m,n\right)$ – in other words, $A$ may be rank-deficient – we seek the minimum norm least squares solution $x$ which minimizes both ${‖x‖}_2$ and ${‖b-Ax‖}_2$.

This chapter ( F08 class) contains driver methods to solve these problems with a single call, as well as computational methods that can be combined with methods in F07 class to solve these linear least squares problems. (F04 not in this release) also contains Black Box methods to solve these linear least squares problems in standard cases. The next two sections discuss the factorizations that can be used in the solution of linear least squares problems.

A number of methods are provided for factorizing a general rectangular $m$ by $n$ matrix $A$, as the product of an orthogonal matrix (unitary if complex) and a triangular (or possibly trapezoidal) matrix.

A real matrix $Q$ is orthogonal if $Q^{\mathrm{T}}Q=I$; a complex matrix $Q$ is unitary if $Q^{\mathrm{H}}Q=I$. Orthogonal or unitary matrices have the important property that they leave the $2$-norm of a vector invariant, so that if $Q$ is orthogonal or unitary. They usually help to maintain numerical stability because they do not amplify rounding errors.

Orthogonal factorizations are used in the solution of linear least squares problems. They may also be used to perform preliminary steps in the solution of eigenvalue or singular value problems, and are useful tools in the solution of a number of other problems.

The singular value decomposition (SVD) of an $m$ by $n$ matrix $A$ is given by where $U$ and $V$ are orthogonal (unitary) and $\Sigma$ is an $m$ by $n$ diagonal matrix with real diagonal elements, ${\sigma }_i$, such that The ${\sigma }_i$ are the singular values of $A$ and the first $\min \left(m,n\right)$ columns of $U$ and $V$ are the left and right singular vectors of $A$. The singular values and singular vectors satisfy where $u_i$ and $v_i$ are the $i$th columns of $U$ and $V$ respectively.

The computation proceeds in the following stages.

1.

The matrix $A$ is reduced to bidiagonal form $A=U_{1}B{V}_1^{\mathrm{T}}$ if $A$ is real ($A=U_{1}B{V}_1^{\mathrm{H}}$ if $A$ is complex), where $U_1$ and $V_1$ are orthogonal (unitary if $A$ is complex), and $B$ is real and upper bidiagonal when $m\ge n$ and lower bidiagonal when $m<n$, so that $B$ is nonzero only on the main diagonal and either on the first superdiagonal (if $m\ge n$) or the first subdiagonal (if $m<n$).

2.

The SVD of the bidiagonal matrix $B$ is computed as $B=U_2\Sigma V_2^{\mathrm{T}}$, where $U_2$ and $V_2$ are orthogonal and $\Sigma$ is diagonal as described above. The singular vectors of $A$ are then $U=U_1{U}_2$ and $V=V_1{V}_2$.

If $m\gg n$, it may be more efficient to first perform a $QR$ factorization of $A$, and then compute the SVD of the $n$ by $n$ matrix $R$, since if $A=QR$ and $R=U\Sigma V^{\mathrm{T}}$, then the SVD of $A$ is given by $A=\left(QU\right)\Sigma V^{\mathrm{T}}$.

Similarly, if $m\ll n$, it may be more efficient to first perform an $LQ$ factorization of $A$.

This chapter supports two primary algorithms for computing the SVD of a bidiagonal matrix. They are:

(i)	the divide and conquer algorithm;
(ii)	the $QR$ algorithm.

The divide and conquer algorithm is much faster than the $QR$ algorithm if singular vectors of large matrices are required.

The SVD may be used to find a minimum norm solution to a (possibly) rank-deficient linear least squares problem (1). The effective rank, $k$, of $A$ can be determined as the number of singular values which exceed a suitable threshold. Let $\hat{\Sigma }$ be the leading $k$ by $k$ sub-matrix of $\Sigma$, and $\hat{V}$ be the matrix consisting of the first $k$ columns of $V$. Then the solution is given by where ${\hat{c}}_1$ consists of the first $k$ elements of $c=U^{\mathrm{T}}b=U_2^{\mathrm{T}}{U}_1^{\mathrm{T}}b$.

The simple type of linear least squares problem described in [Linear Least Squares Problems] can be generalized in various ways.

1.

Linear least squares problems with equality constraints: $$\text{find }x\text{ to minimize }S={{‖c-Ax‖}_2}^2\text{\hspace{1em} subject to \hspace{1em}}Bx=d\text{,}$$ where $A$ is $m$ by $n$ and $B$ is $p$ by $n$, with $p\le n\le m+p$. The equations $Bx=d$ may be regarded as a set of equality constraints on the problem of minimizing $S$. Alternatively the problem may be regarded as solving an overdetermined system of equations $$\left(\begin{array}{c}A\\ B\end{array}\right)x=\left(\begin{array}{c}c\\ d\end{array}\right)\text{,}$$ where some of the equations (those involving $B$) are to be solved exactly, and the others (those involving $A$) are to be solved in a least squares sense. The problem has a unique solution on the assumptions that $B$ has full row rank $p$ and the matrix $\left(\begin{array}{c}A\\ B\end{array}\right)$ has full column rank $n$. (For linear least squares problems with inequality constraints, refer to E04 class.)

2.

General Gauss–Markov linear model problems: $$\text{minimize }{‖y‖}_2\text{\hspace{1em} subject to \hspace{1em}}d=Ax+By\text{,}$$ where $A$ is $m$ by $n$ and $B$ is $m$ by $p$, with $n\le m\le n+p$. When $B=I$, the problem reduces to an ordinary linear least squares problem. When $B$ is square and nonsingular, it is equivalent to a weighted linear least squares problem: $$\text{find }x\text{ to minimize }{‖B^{-1}\left(d-Ax\right)‖}_2\text{.}$$ The problem has a unique solution on the assumptions that $A$ has full column rank $n$, and the matrix $\left(A,B\right)$ has full row rank $m$. Unless $B$ is diagonal, for numerical stability it is generally preferable to solve a weighted linear least squares problem as a general Gauss–Markov linear model problem.

The symmetric eigenvalue problem is to find the eigenvalues, $\lambda$, and corresponding eigenvectors, $z\ne 0$, such that For the Hermitian eigenvalue problem we have For both problems the eigenvalues $\lambda$ are real.

When all eigenvalues and eigenvectors have been computed, we write where $\Lambda$ is a diagonal matrix whose diagonal elements are the eigenvalues, and $Z$ is an orthogonal (or unitary) matrix whose columns are the eigenvectors. This is the classical spectral factorization of $A$.

The basic task of the symmetric eigenproblem methods is to compute values of $\lambda$ and, optionally, corresponding vectors $z$ for a given matrix $A$. This computation proceeds in the following stages.

1.

The real symmetric or complex Hermitian matrix $A$ is reduced to real tridiagonal form $T$. If $A$ is real symmetric this decomposition is $A=QT{Q}^{\mathrm{T}}$ with $Q$ orthogonal and $T$ symmetric tridiagonal. If $A$ is complex Hermitian, the decomposition is $A=QT{Q}^{\mathrm{H}}$ with $Q$ unitary and $T$, as before, real symmetric tridiagonal.

2.

Eigenvalues and eigenvectors of the real symmetric tridiagonal matrix $T$ are computed. If all eigenvalues and eigenvectors are computed, this is equivalent to factorizing $T$ as $T=S\Lambda S^{\mathrm{T}}$, where $S$ is orthogonal and $\Lambda$ is diagonal. The diagonal entries of $\Lambda$ are the eigenvalues of $T$, which are also the eigenvalues of $A$, and the columns of $S$ are the eigenvectors of $T$; the eigenvectors of $A$ are the columns of $Z=QS$, so that $A=Z\Lambda Z^{\mathrm{T}}$ ($Z\Lambda Z^{\mathrm{H}}$ when $A$ is complex Hermitian).

This chapter supports four primary algorithms for computing eigenvalues and eigenvectors of real symmetric matrices and complex Hermitian matrices. They are:

(i)	the divide-and-conquer algorithm;
(ii)	the $QR$ algorithm;
(iii)	bisection followed by inverse iteration;
(iv)	the Relatively Robust Representation (RRR).

The divide-and-conquer algorithm is generally more efficient than the traditional $QR$ algorithm for computing all eigenvalues and eigenvectors, but the RRR algorithm tends to be fastest of all. For further information and references see Anderson et al. (1999).

This section is concerned with the solution of the generalized eigenvalue problems $Az=\lambda Bz$, $ABz=\lambda z$, and $BAz=\lambda z$, where $A$ and $B$ are real symmetric or complex Hermitian and $B$ is positive definite. Each of these problems can be reduced to a standard symmetric eigenvalue problem, using a Cholesky factorization of $B$ as either $B=L{L}^{\mathrm{T}}$ or $B=U^{\mathrm{T}}U$ ($L{L}^{\mathrm{H}}$ or $U^{\mathrm{H}}U$ in the Hermitian case).

With $B=L{L}^{\mathrm{T}}$, we have Hence the eigenvalues of $Az=\lambda Bz$ are those of $Cy=\lambda y$, where $C$ is the symmetric matrix $C=L^{-1}A{L}^{-\mathrm{T}}$ and $y=L^{\mathrm{T}}z$. In the complex case $C$ is Hermitian with $C=L^{-1}A{L}^{-\mathrm{H}}$ and $y=L^{\mathrm{H}}z$.

Table 1 summarises how each of the three types of problem may be reduced to standard form $Cy=\lambda y$, and how the eigenvectors $z$ of the original problem may be recovered from the eigenvectors $y$ of the reduced problem. The table applies to real problems; for complex problems, transposed matrices must be replaced by conjugate-transposes.

When the generalized symmetric-definite problem has been reduced to the corresponding standard problem $Cy=\lambda y$, this may then be solved using the methods described in the previous section. No special methods are needed to recover the eigenvectors $z$ of the generalized problem from the eigenvectors $y$ of the standard problem, because these computations are simple applications of Level 2 or Level 3 BLAS (see F06 class).

Methods which handle symmetric matrices are usually designed so that they use either the upper or lower triangle of the matrix; it is not necessary to store the whole matrix. If either the upper or lower triangle is stored conventionally in the upper or lower triangle of a two-dimensional array, the remaining elements of the array can be used to store other useful data. However, that is not always convenient, and if it is important to economize on storage, the upper or lower triangle can be stored in a one-dimensional array of length $n\left(n+1\right)/2$; that is, the storage is almost halved.

This storage format is referred to as packed storage; it is described in [Packed storage] in the F07 class Chapter Introduction.

Methods designed for packed storage are usually less efficient, especially on high-performance computers, so there is a trade-off between storage and efficiency.

A band matrix is one whose elements are confined to a relatively small number of subdiagonals or superdiagonals on either side of the main diagonal. Algorithms can take advantage of bandedness to reduce the amount of work and storage required. The storage scheme for band matrices is described in [Band storage] in the F07 class Chapter Introduction.

If the problem is the generalized symmetric definite eigenvalue problem $Az=\lambda Bz$ and the matrices $A$ and $B$ are additionally banded, the matrix $C$ as defined in [Generalized Symmetric-definite Eigenvalue Problems] is, in general, full. We can reduce the problem to a banded standard problem by modifying the definition of $C$ thus: where $Q$ is an orthogonal matrix chosen to ensure that $C$ has bandwidth no greater than that of $A$.

A further refinement is possible when $A$ and $B$ are banded, which halves the amount of work required to form $C$. Instead of the standard Cholesky factorization of $B$ as $U^{\mathrm{T}}U$ or $L{L}^{\mathrm{T}}$, we use a split Cholesky factorization $B=S^{\mathrm{T}}S$, where with $U_{11}$ upper triangular and $L_{22}$ lower triangular of order approximately $n/2$; $S$ has the same bandwidth as $B$.

The nonsymmetric eigenvalue problem is to find the eigenvalues, $\lambda$, and corresponding eigenvectors, $v\ne 0$, such that More precisely, a vector $v$ as just defined is called a right eigenvector of $A$, and a vector $u\ne 0$ satisfying is called a left eigenvector of $A$.

A real matrix $A$ may have complex eigenvalues, occurring as complex conjugate pairs.

This problem can be solved via the Schur factorization of $A$, defined in the real case as where $Z$ is an orthogonal matrix and $T$ is an upper quasi-triangular matrix with $1$ by $1$ and $2$ by $2$ diagonal blocks, the $2$ by $2$ blocks corresponding to complex conjugate pairs of eigenvalues of $A$. In the complex case, the Schur factorization is where $Z$ is unitary and $T$ is a complex upper triangular matrix.

The columns of $Z$ are called the Schur vectors. For each $k$ ($1\le k\le n$), the first $k$ columns of $Z$ form an orthonormal basis for the invariant subspace corresponding to the first $k$ eigenvalues on the diagonal of $T$. Because this basis is orthonormal, it is preferable in many applications to compute Schur vectors rather than eigenvectors. It is possible to order the Schur factorization so that any desired set of $k$ eigenvalues occupy the $k$ leading positions on the diagonal of $T$.

The two basic tasks of the nonsymmetric eigenvalue methods are to compute, for a given matrix $A$, all $n$ values of $\lambda$ and, if desired, their associated right eigenvectors $v$ and/or left eigenvectors $u$, and the Schur factorization.

These two basic tasks can be performed in the following stages.

1.	A general matrix $A$ is reduced to upper Hessenberg form $H$ which is zero below the first subdiagonal. The reduction may be written $A=QH{Q}^{\mathrm{T}}$ with $Q$ orthogonal if $A$ is real, or $A=QH{Q}^{\mathrm{H}}$ with $Q$ unitary if $A$ is complex.
2.	The upper Hessenberg matrix $H$ is reduced to Schur form $T$, giving the Schur factorization $H=ST{S}^{\mathrm{T}}$ (for $H$ real) or $H=ST{S}^{\mathrm{H}}$ (for $H$ complex). The matrix $S$ (the Schur vectors of $H$) may optionally be computed as well. Alternatively $S$ may be postmultiplied into the matrix $Q$ determined in stage $1$, to give the matrix $Z=QS$, the Schur vectors of $A$. The eigenvalues are obtained from the diagonal elements or diagonal blocks of $T$.
3.	Given the eigenvalues, the eigenvectors may be computed in two different ways. Inverse iteration can be performed on $H$ to compute the eigenvectors of $H$, and then the eigenvectors can be multiplied by the matrix $Q$ in order to transform them to eigenvectors of $A$. Alternatively the eigenvectors of $T$ can be computed, and optionally transformed to those of $H$ or $A$ if the matrix $S$ or $Z$ is supplied.

The accuracy with which eigenvalues can be obtained can often be improved by balancing a matrix. This is discussed further in [Balancing and condition for the nonsymmetric eigenproblem] below.

The generalized nonsymmetric eigenvalue problem is to find the eigenvalues, $\lambda$, and corresponding eigenvectors, $v\ne 0$, such that

More precisely, a vector $v$ as just defined is called a right eigenvector of the matrix pair $\left(A,B\right)$, and a vector $u\ne 0$ satisfying is called a left eigenvector of the matrix pair $\left(A,B\right)$.

If $B$ is singular then the problem has one or more infinite eigenvalues $\lambda =\infty$, corresponding to $Bv=0$. Note that if $A$ is nonsingular, then the equivalent problem $\mu Av=Bv$ is perfectly well defined and an infinite eigenvalue corresponds to $\mu =0$. To deal with both finite (including zero) and infinite eigenvalues, the methods in this chapter do not compute $\lambda$ explicitly, but rather return a pair of numbers $\left(\alpha ,\beta \right)$ such that if $\beta \ne 0$  and if $\alpha \ne 0$ and $\beta =0$ then $\lambda =\infty$. $\beta$ is always returned as real and non-negative. Of course, computationally an infinite eigenvalue may correspond to a small $\beta$ rather than an exact zero.

For a given pair $\left(A,B\right)$ the set of all the matrices of the form $\left(A-\lambda B\right)$ is called a matrix pencil and $\lambda$ and $v$ are said to be an eigenvalue and eigenvector of the pencil $\left(A-\lambda B\right)$. If $A$ and $B$ are both singular and share a common null space then so that the pencil $\left(A-\lambda B\right)$ is singular for all $\lambda$. In other words any $\lambda$ can be regarded as an eigenvalue. In exact arithmetic a singular pencil will have $\alpha =\beta =0$ for some $\left(\alpha ,\beta \right)$. Computationally if some pair $\left(\alpha ,\beta \right)$ is small then the pencil is singular, or nearly singular, and no reliance can be placed on any of the computed eigenvalues. Singular pencils can also manifest themselves in other ways; see, in particular, Sections 2.3.5.2 and 4.11.1.4 of Anderson et al. (1999) for further details.

The generalized eigenvalue problem can be solved via the generalized Schur factorization of the pair $\left(A,B\right)$ defined in the real case as where $Q$ and $Z$ are orthogonal, $T$ is upper triangular with non-negative diagonal elements and $S$ is upper quasi-triangular with $1$ by $1$ and $2$ by $2$ diagonal blocks, the $2$ by $2$ blocks corresponding to complex conjugate pairs of eigenvalues. In the complex case, the generalized Schur factorization is where $Q$ and $Z$ are unitary and $S$ and $T$ are upper triangular, with $T$ having real non-negative diagonal elements. The columns of $Q$ and $Z$ are called respectively the left and right generalized Schur vectors and span pairs of deflating subspaces of $A$ and $B$, which are a generalization of invariant subspaces.

It is possible to order the generalized Schur factorization so that any desired set of $k$ eigenvalues correspond to the $k$ leading positions on the diagonals of the pair $\left(S,T\right)$.

The two basic tasks of the generalized nonsymmetric eigenvalue methods are to compute, for a given pair $\left(A,B\right)$, all $n$ values of $\lambda$ and, if desired, their associated right eigenvectors $v$ and/or left eigenvectors $u$, and the generalized Schur factorization.

These two basic tasks can be performed in the following stages.

1.	The matrix pair $\left(A,B\right)$ is reduced to generalized upper Hessenberg form $\left(H,R\right)$, where $H$ is upper Hessenberg (zero below the first subdiagonal) and $R$ is upper triangular. The reduction may be written as $A=Q_{1}H{Z}_1^{\mathrm{T}},B=Q_{1}R{Z}_1^{\mathrm{T}}$ in the real case with $Q_1$ and $Z_1$ orthogonal, and $A=Q_{1}H{Z}_1^{\mathrm{H}},B=Q_{1}R{Z}_1^{\mathrm{H}}$ in the complex case with $Q_1$ and $Z_1$ unitary.
2.	The generalized upper Hessenberg form $\left(H,R\right)$ is reduced to the generalized Schur form $\left(S,T\right)$ using the generalized Schur factorization $H=Q_{2}S{Z}_2^{\mathrm{T}}$, $R=Q_{2}T{Z}_2^{\mathrm{T}}$ in the real case with $Q_2$ and $Z_2$ orthogonal, and $H=Q_{2}S{Z}_2^{\mathrm{H}},R=Q_{2}T{Z}_2^{\mathrm{H}}$ in the complex case. The generalized Schur vectors of $\left(A,B\right)$ are given by $Q=Q_1{Q}_2$, $Z=Z_1{Z}_2$. The eigenvalues are obtained from the diagonal elements (or blocks) of the pair $\left(S,T\right)$.
3.	Given the eigenvalues, the eigenvectors of the pair $\left(S,T\right)$ can be computed, and optionally transformed to those of $\left(H,R\right)$ or $\left(A,B\right)$.

The accuracy with which eigenvalues can be obtained can often be improved by balancing a matrix pair. This is discussed further in [Balancing the generalized eigenvalue problem] below.

The Sylvester equation is a matrix equation of the form where $A$, $B$, and $C$ are given matrices with $A$ being $m$ by $m$, $B$ an $n$ by $n$ matrix and $C$, and the solution matrix $X$, $m$ by $n$ matrices. The solution of a special case of this equation occurs in the computation of the condition number for an invariant subspace, but a combination of methods in this chapter allows the solution of the general Sylvester equation.

Methods are also provided for solving a special case of the generalized Sylvester equations where $\left(A,D\right)$, $\left(B,E\right)$ and $\left(C,F\right)$ are given matrix pairs, and $R$ and $L$ are the solution matrices.

In this section we discuss the effects of rounding errors in the solution process and the effects of uncertainties in the data, on the solution to the problem. A number of the methods in this chapter return information, such as condition numbers, that allow these effects to be assessed. First we discuss some notation used in the error bounds of later sections.

The bounds usually contain the factor $p\left(n\right)$ (or $p\left(m,n\right)$), which grows as a function of the matrix dimension $n$ (or matrix dimensions $m$ and $n$). It measures how errors can grow as a function of the matrix dimension, and represents a potentially different function for each problem. In practice, it usually grows just linearly; $p\left(n\right)\le 10n$ is often true, although generally only much weaker bounds can be actually proved. We normally describe $p\left(n\right)$ as a ‘modestly growing’ function of $n$. For detailed derivations of various $p\left(n\right)$, see Golub and Van Loan (2012) and Wilkinson (1965).

For linear equation (see F07 class) and least squares solvers, we consider bounds on the relative error $‖x-\hat{x}‖/‖x‖$ in the computed solution $\hat{x}$, where $x$ is the true solution. For eigenvalue problems we consider bounds on the error $\left|{\lambda }_i-{\hat{\lambda }}_i\right|$ in the $i$th computed eigenvalue ${\hat{\lambda }}_i$, where ${\lambda }_i$ is the true $i$th eigenvalue. For singular value problems we similarly consider bounds $\left|{\sigma }_i-{\hat{\sigma }}_i\right|$.

Bounding the error in computed eigenvectors and singular vectors ${\hat{v}}_i$ is more subtle because these vectors are not unique: even though we restrict ${‖{\hat{v}}_i‖}_2=1$ and ${‖v_i‖}_2=1$, we may still multiply them by arbitrary constants of absolute value $1$. So to avoid ambiguity we bound the angular difference between ${\hat{v}}_i$ and the true vector $v_i$, so that Here $\arccos \left(\theta \right)$ is in the standard range: $0\le \arccos \left(\theta \right)<\pi$. When $\theta \left(v_i,{\hat{v}}_i\right)$ is small, we can choose a constant $\alpha$ with absolute value $1$ so that ${‖\alpha v_i-{\hat{v}}_i‖}_2\approx \theta \left(v_i,{\hat{v}}_i\right)$.

In addition to bounds for individual eigenvectors, bounds can be obtained for the spaces spanned by collections of eigenvectors. These may be much more accurately determined than the individual eigenvectors which span them. These spaces are called invariant subspaces in the case of eigenvectors, because if $v$ is any vector in the space, $Av$ is also in the space, where $A$ is the matrix. Again, we will use angle to measure the difference between a computed space $\hat{S}$ and the true space $S$:

$\theta \left(S,\hat{S}\right)$ may be computed as follows. Let $S$ be a matrix whose columns are orthonormal and $\mathrm{span} S$. Similarly let $\hat{S}$ be an orthonormal matrix with columns spanning $\hat{S}$. Then Finally, we remark on the accuracy of the bounds when they are large. Relative errors like $‖\hat{x}-x‖/‖x‖$ and angular errors like $\theta \left({\hat{v}}_i,v_i\right)$ are only of interest when they are much less than $1$. Some stated bounds are not strictly true when they are close to $1$, but rigorous bounds are much more complicated and supply little extra information in the interesting case of small errors. These bounds are indicated by using the symbol $\lesssim$, or ‘approximately less than’, instead of the usual $\le$. Thus, when these bounds are close to 1 or greater, they indicate that the computed answer may have no significant digits at all, but do not otherwise bound the error.