Numerical linear algebra

Numerical linear algebra,sometimes calledapplied linear algebra,is the study of howmatrix operationscan be used to createcomputer algorithmswhichefficientlyand accurately provide approximate answers to questions incontinuousmathematics. It is a subfield ofnumerical analysis,and a type oflinear algebra.Computersusefloating-point arithmeticand cannot exactly representirrationaldata, so when a computer algorithm is applied to a matrix of data, it can sometimesincrease the differencebetween a number stored in the computer and the true number that it is an approximation of. Numerical linear algebra uses properties of vectors and matrices to develop computer algorithms that minimize the error introduced by the computer, and is also concerned with ensuring that the algorithm is as efficient as possible.

Numerical linear algebra aims to solve problems of continuousmathematicsusing finite precision computers, so its applications to thenaturalandsocial sciencesare as vast as the applications of continuous mathematics. It is often a fundamental part ofengineeringandcomputational scienceproblems, such asimageandsignal processing,telecommunication,computational finance,materials sciencesimulations,structural biology,data mining,bioinformatics,andfluid dynamics.Matrix methods are particularly used infinite difference methods,finite element methods,and the modeling ofdifferential equations.Noting the broad applications of numerical linear algebra,Lloyd N. Trefethenand David Bau, III argue that it is "as fundamental to the mathematical sciences as calculus and differential equations",even though it is a comparatively small field.Because many properties of matrices and vectors also apply to functions and operators, numerical linear algebra can also be viewed as a type offunctional analysiswhich has a particular emphasis on practical algorithms.

Common problems in numerical linear algebra include obtaining matrix decompositions like thesingular value decomposition,theQR factorization,theLU factorization,or theeigendecomposition,which can then be used to answer common linear algebraic problems like solving linear systems of equations, locating eigenvalues, or least squares optimisation. Numerical linear algebra's central concern with developing algorithms that do not introduce errors when applied to real data on a finite precision computer is often achieved byiterativemethods rather than direct ones.

History

Numerical linear algebra was developed by computer pioneers likeJohn von Neumann,Alan Turing,James H. Wilkinson,Alston Scott Householder,George Forsythe,andHeinz Rutishauser,in order to apply the earliest computers to problems in continuous mathematics, such as ballistics problems and the solutions to systems ofpartial differential equations.The first serious attempt to minimize computer error in the application of algorithms to real data is John von Neumann andHerman Goldstine's work in 1947.The field has grown as technology has increasingly enabled researchers to solve complex problems on extremely large high-precision matrices, and some numerical algorithms have grown in prominence as technologies likeparallel computinghave made them practical approaches to scientific problems.

Matrix decompositions

Further information:Matrix decomposition

Partitioned matrices

Main article:Block matrix

For many problems in applied linear algebra, it is useful to adopt the perspective of a matrix as being a concatenation of column vectors. For example, when solving the linear system$x=A^{-1}b$,rather than understandingxas the product of$A^{-1}$withb,it is helpful to think ofxas the vector ofcoefficientsin the linear expansion ofbin thebasisformed by the columns ofA.Thinking of matrices as a concatenation of columns is also a practical approach for the purposes of matrix algorithms. This is because matrix algorithms frequently contain two nested loops: one over the columns of a matrixA,and another over the rows ofA.For example, for matrices$A^{m\times n}$and vectors$x^{n\times 1}$and$y^{m\times 1}$,we could use the column partitioning perspective to computey:=Ax+yas

forq=1:n
forp=1:m
y(p)=A(p,q)*x(q)+y(p)
end
end

Singular value decomposition

Main article:singular value decomposition

The singular value decomposition of a matrix$A^{m\times n}$is$A=U\Sigma V^{\ast }$whereUandVareunitary,and$\Sigma$isdiagonal.The diagonal entries of$\Sigma$are called thesingular valuesofA.Because singular values are the square roots of theeigenvaluesof$AA^{\ast }$,there is a tight connection between the singular value decomposition and eigenvalue decompositions. This means that most methods for computing the singular value decomposition are similar to eigenvalue methods;perhaps the most common method involvesHouseholder procedures.

QR factorization

Main article:QR decomposition

The QR factorization of a matrix$A^{m\times n}$is a matrix$Q^{m\times m}$and a matrix$R^{m\times n}$so thatA = QR,whereQisorthogonalandRisupper triangular.The two main algorithms for computing QR factorizations are theGram–Schmidt processand theHouseholder transformation.The QR factorization is often used to solvelinear least-squaresproblems, and eigenvalue problems (by way of the iterativeQR algorithm).

LU factorization

Main article:LU decomposition

An LU factorization of a matrixAconsists of a lower triangular matrixLand an upper triangular matrixUso thatA = LU.The matrixUis found by an upper triangularization procedure which involves left-multiplyingAby a series of matrices$M_{1},\ldots,M_{n-1}$to form the product$M_{n-1}\cdots M_{1}A=U$,so that equivalently$L=M_{1}^{-1}\cdots M_{n-1}^{-1}$.

Eigenvalue decomposition

Main article:Eigendecomposition of a matrix

The eigenvalue decomposition of a matrix$A^{m\times m}$is$A=X\Lambda X^{-1}$,where the columns ofXare the eigenvectors ofA,and$\Lambda$is a diagonal matrix the diagonal entries of which are the corresponding eigenvalues ofA.There is no direct method for finding the eigenvalue decomposition of an arbitrary matrix. Because it is not possible to write a program that finds the exact roots of an arbitrary polynomial in finite time, any general eigenvalue solver must necessarily be iterative.

Algorithms

Gaussian elimination

Main article:Gaussian elimination

From the numerical linear algebra perspective, Gaussian elimination is a procedure for factoring a matrixAinto itsLUfactorization, which Gaussian elimination accomplishes by left-multiplyingAby a succession of matrices$L_{m-1}\cdots L_{2}L_{1}A=U$untilUis upper triangular andLis lower triangular, where$L\equiv L_{1}^{-1}L_{2}^{-1}\cdots L_{m-1}^{-1}$.Naive programs for Gaussian elimination are notoriously highly unstable, and produce huge errors when applied to matrices with many significant digits.The simplest solution is to introducepivoting,which produces a modified Gaussian elimination algorithm that is stable.

Solutions of linear systems

Further information:System of linear equations § Other methods

Numerical linear algebra characteristically approaches matrices as a concatenation of columns vectors. In order to solve the linear system$x=A^{-1}b$,the traditional algebraic approach is to understandxas the product of$A^{-1}$withb.Numerical linear algebra instead interpretsxas the vector of coefficients of the linear expansion ofbin the basis formed by the columns ofA.

Many different decompositions can be used to solve the linear problem, depending on the characteristics of the matrixAand the vectorsxandb,which may make one factorization much easier to obtain than others. IfA=QRis a QR factorization ofA,then equivalently$Rx=Q^{\ast }b$.This is as easy to compute as a matrix factorization.If$A=X\Lambda X^{-1}$is an eigendecompositionA,and we seek to findbso thatb=Ax,with$b'=X^{-1}b$and$x'=X^{-1}x$,then we have$b'=\Lambda x'$.This is closely related to the solution to the linear system using the singular value decomposition, because singular values of a matrix are the absolute values of its eigenvalues, which are also equivalent to the square roots of the absolute values of the eigenvalues of the Gram matrix$X^{*}X$.And ifA=LUis anLUfactorization ofA,thenAx=bcan be solved using the triangular matricesLy=bandUx=y.

Least squares optimisation

Main article:Numerical methods for linear least squares

Matrix decompositions suggest a number of ways to solve the linear systemr=b−Axwhere we seek to minimizer,as in theregression problem.The QR algorithm solves this problem by computing the reduced QR factorization ofAand rearranging to obtain${\widehat {R}}x={\widehat {Q}}^{\ast }b$.This upper triangular system can then be solved forx.The SVD also suggests an algorithm for obtaining linear least squares. By computing the reduced SVD decomposition$A={\widehat {U}}{\widehat {\Sigma }}V^{\ast }$and then computing the vector${\widehat {U}}^{\ast }b$,we reduce the least squares problem to a simple diagonal system.The fact that least squares solutions can be produced by the QR and SVD factorizations means that, in addition to the classicalnormal equationsmethod for solving least squares problems, these problems can also be solved by methods that include the Gram-Schmidt algorithm and Householder methods.

Conditioning and stability

Main article:Numerical analysis § Numerical stability and well-posed problems

Allow that a problem is a function$f:X\to Y$,whereXis a normed vector space of data andYis a normed vector space of solutions. For some data point$x\in X$,the problem is said to be ill-conditioned if a small perturbation inxproduces a large change in the value off(x). We can quantify this by defining acondition numberwhich represents how well-conditioned a problem is, defined as $${\widehat {\kappa }}=\lim _{\delta \to 0}\sup _{\|\delta x\|\leq \delta }{\frac {\|\delta f\|}{\|\delta x\|}}.$$

Instability is the tendency of computer algorithms, which depend onfloating-point arithmetic,to produce results that differ dramatically from the exact mathematical solution to a problem. When a matrix contains real data with manysignificant digits,many algorithms for solving problems like linear systems of equation or least squares optimisation may produce highly inaccurate results. Creating stable algorithms for ill-conditioned problems is a central concern in numerical linear algebra. One example is that the stability of householder triangularization makes it a particularly robust solution method for linear systems, whereas the instability of the normal equations method for solving least squares problems is a reason to favour matrix decomposition methods like using the singular value decomposition. Some matrix decomposition methods may be unstable, but have straightforward modifications that make them stable; one example is the unstable Gram–Schmidt, which can easily be changed to produce the stablemodified Gram–Schmidt.Another classical problem in numerical linear algebra is the finding that Gaussian elimination is unstable, but becomes stable with the introduction of pivoting.

Iterative methods

Main article:Iterative methods

There are two reasons that iterative algorithms are an important part of numerical linear algebra. First, many important numerical problems have no direct solution; in order to find the eigenvalues and eigenvectors of an arbitrary matrix, we can only adopt an iterative approach. Second, noniterative algorithms for an arbitrary$m\times m$matrix require$O(m^{3})$time, which is a surprisingly high floor given that matrices contain only$m^{2}$numbers. Iterative approaches can take advantage of several features of some matrices to reduce this time. For example, when a matrix issparse,an iterative algorithm can skip many of the steps that a direct approach would necessarily follow, even if they are redundant steps given a highly structured matrix.

The core of many iterative methods in numerical linear algebra is the projection of a matrix onto a lower dimensionalKrylov subspace,which allows features of a high-dimensional matrix to be approximated by iteratively computing the equivalent features of similar matrices starting in a low dimension space and moving to successively higher dimensions. WhenAis symmetric and we wish to solve the linear problemAx=b,the classical iterative approach is theconjugate gradient method.IfAis not symmetric, then examples of iterative solutions to the linear problem are thegeneralized minimal residual methodandCGN.IfAis symmetric, then to solve the eigenvalue and eigenvector problem we can use theLanczos algorithm,and ifAis non-symmetric, then we can useArnoldi iteration.

Software

Main article:List of numerical analysis software

Several programming languages use numerical linear algebra optimisation techniques and are designed to implement numerical linear algebra algorithms. These languages includeMATLAB,Analytica,Maple,andMathematica.Other programming languages which are not explicitly designed for numerical linear algebra have libraries that provide numerical linear algebra routines and optimisation;CandFortranhave packages likeBasic Linear Algebra SubprogramsandLAPACK,Pythonhas the libraryNumPy,andPerlhas thePerl Data Language.Many numerical linear algebra commands inRrely on these more fundamental libraries likeLAPACK.More libraries can be found on theList of numerical libraries.

Further reading

• Dongarra, Jack;Hammarling, Sven (1990). "Evolution of Numerical Software for Dense Linear Algebra". In Cox, M. G.; Hammarling, S. (eds.).Reliable Numerical Computation.Oxford: Clarendon Press. pp. 297–327.ISBN0-19-853564-3.

External links

• Freely available software for numerical algebra on the web,composed by Jack Dongarra and Hatem Ltaief, University of Tennessee
• NAG Library of numerical linear algebra routines
• Numerical Linear Algebra GrouponX(Research group in the United Kingdom)
• siaglaonX(Activity group about numerical linear algebra in theSociety for Industrial and Applied Mathematics)
• The GAMM Activity Group on Applied and Numerical Linear Algebra