https://nhigham.com/2022/10/11/seven-sins-of-numerical-linear-algebra/ Skip to content [logo-blue-on-white_narrower] Nick Higham Applied mathematics, numerical linear algebra and software. Primary Menu * Home * Blog * Papers + Linear Systems and Condition Estimation + Least Squares Problems + Correlation Matrices + Eigenvalue Problems + Matrix Functions and Nonlinear Matrix Equations + Miscellaneous Papers + Articles and Book Reviews * Books + Accuracy and Stability of Numerical Algorithms + Functions of Matrices: Theory and Computation + Handbook of Writing for the Mathematical Sciences + MATLAB Guide + The Princeton Companion to Applied Mathematics + Penguin Dictionary of Mathematics * Conferences * Talks + Videos + Slides * Resources * People * Photos * Miscellanea * Popular Posts * SIAM News + Posts at SIAM News + SIAM News home page * "What Is" Seven Sins of Numerical Linear Algebra In numerical linear algebra we are concerned with solving linear algebra problems accurately and efficiently and understanding the sensitivity of the problems to perturbations. We describe seven sins, whereby accuracy or efficiency is lost or misleading information about sensitivity is obtained. 1. Inverting a Matrix In linear algebra courses we learn that the solution to a linear system Ax = b of n equations in n unknowns can be written x = A^{-1}b , where A^{-1} is the matrix inverse. What is not always emphasized is that there are very few circumstances in which one should compute A^{-1}. Indeed one would not solve the scalar (n=1) system 7x = 21 by computing x = 7^{-1} \times 21, but rather would carry out a division x = 21/7. In the n\times n case, it is faster and more accurate to solve a linear system by LU factorization (Gaussian elimination) with partial pivoting than by inverting A (which has, in any case, to be done by LU factorization). Rare cases where A^{-1} is required are in statistics, where the diagonal elements of the inverse of the covariance matrix are relevant quantities, and in certain algorithms for computing matrix functions. 2. Forming the Cross-Product Matrix A^TA The solution to the linear least squares problem \min_x\| b - Ax \|_2 , where A is a full-rank m\times n matrix with m\ge n, satisfies the normal equations A^T\!A x = A^Tb. It is therefore natural to form the symmetric positive definite matrix A^T\!A and solve the normal equations by Cholesky factorization. While fast, this method is numerically unstable when A is ill conditioned. By contrast, solving the least squares problem via QR factorization is always numerically stable. What is wrong with the cross-product matrix A^T\!A (also known as the Gram matrix)? It squares the data, which can cause a loss of information in floating-point arithmetic. For example, if A = \begin{bmatrix} 1 & 1 \\ \epsilon & 0 \end{bmatrix}, \quad 0 < \ epsilon < \sqrt{u}, where u is the unit roundoff of the floating point arithmetic, then A^T\!A = \begin{bmatrix} 1 + \epsilon^2 & 1 \\ 1 & 1 \end{bmatrix} is positive definite but, since \epsilon^2 0, * |a_{ij}| < \sqrt{a_{ii}a_{jj}} for all i \ne j, but none of these conditions, or even all taken together, guarantees that the matrix has positive eigenvalues. The best way to check definiteness is to compute a Cholesky factorization, which is often needed anyway. The MATLAB function chol returns an error message if the factorization fails, and a second output argument can be requested, which is set to the number of the stage on which the factorization failed, or to zero if the factorization succeeded. In the case of failure, the partially computed R factor is returned in the first argument, and it can be used to compute a direction of negative curvature (as needed in optimization), for example. This sin takes the top spot in Schmelzer and Hauser's Seven Sins in Portfolio Optimization, because in portfolio optimization a negative eigenvalue in the covariance matrix can identify a portfolio with negative variance, promising an arbitrarily large investment with no risk! 5. Not Exploiting Structure in the Matrix One of the fundamental tenets of numerical linear algebra is that one should try to exploit any matrix structure that might be present. Sparsity (a matrix having a large number of zeros) is particularly important to exploit, since algorithms intended for dense matrices may be impractical for sparse matrices because of extensive fill-in (zeros becoming nonzero). Here are two examples of structures that can be exploited. Matrices from saddle point problems are symmetric indefinite and of the form \notag C = \begin{bmatrix} A & B^T \\ B & 0 \end{bmatrix}, with A symmetric positive definite. Much work has been done on developing numerical methods for solving Cx = b that exploit the block structure and possible sparsity in A and B. A second example is a circulant matrix \notag C = \begin{bmatrix} c_1 & c_2 & \dots & c_n \\ c_n & c_1 & \ dots & \vdots \\ \vdots & \ddots & \ddots & c_2 \\ c_2 & \dots & c_n & c_1 \\ \end{bmatrix}. Circulant matrices have the important property that they are diagonalized by a unitary matrix called the discrete Fourier transform matrix. Using this property one can solve Cx = v in O(n \ log_2n) operations, rather than the O(n^3) operations required if the circulant structure is ignored. Ideally, linear algebra software would detect structure in a matrix and call an algorithm that exploits that structure. A notable example of such a meta-algorithm is the MATLAB backslash function x = A\b for solving Ax = b. Backslash checks whether the matrix is triangular (or a permutation of a triangular matrix), upper Hessenberg, symmetric, or symmetric positive definite, and applies an appropriate method. It also allows A to be rectangular and solves the least squares problem if there are more rows than columns and the underdetermined system if there are more columns than rows. 6. Using the Determinant to Detect Near Singularity An n\times n matrix A is nonsingular if and only if its determinant is nonzero. One might therefore expect that a small value for \det(A) indicates a matrix that is nearly singular. However, the size of \det (A) tells us nothing about near singularity. Indeed, since \det(\ alpha A) = \alpha^n \det(A) we can achieve any value for the determinant by multiplying by a scalar \alpha, yet \alpha A is no more or less nearly singular than A for \alpha \ne 0. Another limitation of the determinant is shown by the two matrices \notag T = \begin{bmatrix} 1 & -1 & -1 & \dots & -1\\ & 1 & -1 & \ dots & -1\\ & & 1 & \dots & \vdots\\ & & & \ddots & -1 \\ & & & & 1 \ end{bmatrix}, \quad U = \begin{bmatrix} 1 & 1 & 1 & \dots & 1\\ & 1 & 1 & \dots & 1\\ & & 1 & \dots & \vdots\\ & & & \ddots & 1 \\ & & & & 1 \end{bmatrix} \qquad (1) Both matrices have unit diagonal and off-diagonal elements bounded in modulus by 1. So \det(T) = \det(U) = 1, yet \notag T^{-1} = \begin{bmatrix} 1 & 1 & 2 & \dots & 2^{n-2}\\ & 1 & 1 & \dots & \vdots\\ & & 1 & \ddots & 2\\ & & & \ddots & 1 \\ & & & & 1 \end{bmatrix}, \quad U^{-1} = \begin{bmatrix} 1 & -1 & & & \\ & 1 & -1 & & \\ & & 1 & \ddots & \\ & & & \ddots & -1 \\ & & & & 1 \end {bmatrix}. So T is ill conditioned for large n. In fact, if we change the (n,1) element of T to -2^{n-2} then the matrix becomes singular! By contrast, U is always very well conditioned. The determinant cannot distinguish between the ill-conditioned T and the well-conditioned U. 7. Using Eigenvalues to Estimate Conditioning For any n\times n matrix A and any consistent matrix norm it is true that \|A\| \ge |\lambda_i| for all i, where the \lambda_i are the eigenvalue of A. Since the eigenvalues of A^{-1} are \lambda^{-1}, it follows that the matrix condition number \kappa(A) = \|A\| \, \|A^ {-1}\| is bounded below by the ratio of largest to smallest eigenvalue in absolute value, that is, \notag \kappa(A) \ge \displaystyle\frac{ \max_i |\lambda_i| } { \ min_i |\lambda_i| }. But as the matrix T in (1) shows, this bound can be very weak. It is singular values not eigenvalues that characterize the condition number for the 2-norm. Specifically, \notag \kappa_2(A) = \displaystyle\frac{\sigma_1}{\sigma_n}, where A = U\Sigma V^T is a singular value decomposition (SVD), with U and V orthogonal and \Sigma = \mathrm{diag}(\sigma_i), \sigma_1 \ge \ sigma_2 \ge \cdots \ge \sigma_n \ge 0. If A is symmetric, for example, then the sets \{ |\lambda_i| \} and \{\sigma_i \} are the same, but in general the eigenvalues \lambda_i and singular values \ sigma_i can be very different. Related Blog Posts * What Is a Cholesky Factorization? (2020) * What Is a Circulant Matrix? (2022) * What Is a Condition Number? (2020) * What is a Sparse Matrix? (2020) * What Is a Symmetric Positive Definite Matrix? (2020) * What Is the Singular Value Decomposition? (2020) Share this: * Print * Email * Twitter * More * * Facebook * Reddit * * LinkedIn * Related Posted on October 11, 2022October 11, 2022 by Nick HighamPosted in research Post navigation Previous Previous post: Cleve Moler Wins ICIAM Industry Prize 2023 Leave a Reply Cancel reply Enter your comment here... [ ] Fill in your details below or click an icon to log in: * * * * Gravatar Email (required) (Address never made public) [ ] Name (required) [ ] Website [ ] WordPress.com Logo You are commenting using your WordPress.com account. ( Log Out / Change ) Twitter picture You are commenting using your Twitter account. ( Log Out / Change ) Facebook photo You are commenting using your Facebook account. ( Log Out / Change ) Cancel Connecting to %s [ ] Notify me of new comments via email. [ ] Notify me of new posts via email. [Post Comment] [ ] [ ] [ ] [ ] [ ] [ ] [ ] D[ ] Search for: [ ] [Search] Recent Posts * Seven Sins of Numerical Linear Algebra * Cleve Moler Wins ICIAM Industry Prize 2023 * What Is a Circulant Matrix? * Photos and Videos from NJH60 Conference * What Is Fast Matrix Multiplication? Twitter Feed My Tweets Recent Comments * Cleve Moler on Cleve Moler Wins ICIAM Industry Prize 2023 * Aaron Meurer on Half Precision Arithmetic: fp16 Versus bfloat16 * Jose-Javier Martinez on What Is the Pascal Matrix? * Nick Higham on The Big Six Matrix Factorizations * Pieter Ghysels on The Big Six Matrix Factorizations Categories * books * conferences * Emacs * LaTeX * matrix computations * miscellaneous * people * Princeton Companion * publication peculiarities * publishing * research * software * what-is * writing Tags Adobe_Acrobat Algol ASCII_art Basic Beamer bfloat16 BibTeX blogs bohemian_matrices C CMYK comma Commodore_64 Commodore_Pet correlations creativity dictionary DOI ellipsis error_analysis Forth Fortran fp16 GitHub Helm Householder_symposium IEEE_arithmetic ill_posed_problem IMA Julia Lambert_W_function Lanczos LINPACK Lisp lists logarithm Mac machine_learning Manchester mathematician MATLAB matrix_function Moler NAG Org orthogonal_matrix Pandoc paper PCAM PDF PDF_viewers pencil performance_profile photocopying photography Princeton_University_press punctuation Python rounding scanning SIAM SIAM_News softmax spreadsheets stationery talk Twitter typesetting unwinding_function Vandermonde video wilkinson Windows Subscribe to Blog via Email Enter your email address to subscribe to this blog and receive notifications of new posts by email. Email Address: [ ] Subscribe Join 405 other followers Archives * October 2022 * September 2022 * July 2022 * June 2022 * May 2022 * April 2022 * March 2022 * February 2022 * January 2022 * December 2021 * November 2021 * October 2021 * September 2021 * August 2021 * July 2021 * June 2021 * May 2021 * April 2021 * March 2021 * February 2021 * January 2021 * December 2020 * November 2020 * October 2020 * September 2020 * August 2020 * July 2020 * June 2020 * May 2020 * April 2020 * March 2020 * February 2020 * December 2019 * November 2019 * August 2019 * June 2019 * May 2019 * January 2019 * December 2018 * November 2018 * October 2018 * August 2018 * June 2018 * April 2018 * January 2018 * December 2017 * November 2017 * October 2017 * August 2017 * July 2017 * June 2017 * May 2017 * April 2017 * March 2017 * February 2017 * January 2017 * December 2016 * November 2016 * October 2016 * September 2016 * August 2016 * June 2016 * May 2016 * April 2016 * March 2016 * February 2016 * January 2016 * December 2015 * November 2015 * October 2015 * September 2015 * August 2015 * July 2015 * June 2015 * May 2015 * April 2015 * February 2015 * December 2014 * November 2014 * September 2014 * August 2014 * July 2014 * June 2014 * May 2014 * April 2014 * March 2014 * February 2014 * January 2014 * December 2013 * November 2013 * October 2013 * September 2013 * August 2013 * July 2013 * June 2013 * May 2013 * April 2013 * March 2013 * February 2013 * January 2013 RSS Feed * RSS - Posts * RSS - Comments Search for: [ ] [Search] Recent Posts * Seven Sins of Numerical Linear Algebra * Cleve Moler Wins ICIAM Industry Prize 2023 * What Is a Circulant Matrix? * Photos and Videos from NJH60 Conference * What Is Fast Matrix Multiplication? Recent Comments * Cleve Moler on Cleve Moler Wins ICIAM Industry Prize 2023 * Aaron Meurer on Half Precision Arithmetic: fp16 Versus bfloat16 * Jose-Javier Martinez on What Is the Pascal Matrix? * Nick Higham on The Big Six Matrix Factorizations * Pieter Ghysels on The Big Six Matrix Factorizations Categories * books (18) * conferences (32) * Emacs (9) * LaTeX (16) * matrix computations (7) * miscellaneous (16) * people (17) * Princeton Companion (12) * publication peculiarities (7) * publishing (2) * research (30) * software (30) * what-is (76) * writing (16) Top Posts in Last 2 Days * Seven Sins of Numerical Linear Algebra * What Is a Cholesky Factorization? * The Top 10 Algorithms in Applied Mathematics * What Is Fast Matrix Multiplication? * What Is the Singular Value Decomposition? * What Is a Condition Number? * Better LaTeX Tables with Booktabs * What Is a Symmetric Positive Definite Matrix? * What Is a Circulant Matrix? * What Is an LU Factorization? Recent Posts from NLA Group: Numerical Linear Algebra Group Talks for NJH60 Now Available on YouTube A Logo for the Numerical Linear Algebra Group Our Alumni - Younes Chahlaoui Nick Higham to give International Congress of Mathematicians Lecture 2022 NLA group photo Website Powered by WordPress.com. * Follow Following + [croppe] Nick Higham Join 405 other followers [ ] Sign me up + Already have a WordPress.com account? Log in now. * + [croppe] Nick Higham + Customize + Follow Following + Sign up + Log in + Copy shortlink + Report this content + View post in Reader + Manage subscriptions + Collapse this bar [b]