Direct methods for sparse matrix solution
Curator: Iain Duff
Evgeni Sergeev
Jack Dongarra
Barbara Zubik-Kowal
Prof. Iain Duff , Rutherford Appleton Laboratory, Chilton, Oxfordshire, UK
Bora Uçar , CNRS and ENS Lyon, France
Direct methods for sparse matrix solutions are characterized by using a matrix factorization to solve a set of equations of the form
where \(b\) is a given vector, \(x\) is the vector of unknowns and \(A\) is a given sparse matrix representing the coefficients of unknowns in each equation. In contrast to iterative methods , direct methods obtain the solution to the above system in a finite and fixed number of steps.
In direct methods, the matrix \(A\) is factorized into a product of other sparse matrices. Systems of equations with the factor matrices as a coefficient matrix are very easy to solve.
Factorization or decomposition
Depending on the properties of the matrix \(A\), different factorizations are used:
- For an \(n\times n\) symmetric positive definite matrix, the Cholesky factorization \(A=LL^T\) is usually computed, where \(L\) is a lower triangular (sparse) matrix. Often an \(LDL^T\) decomposition is used to avoid taking square roots and to avoid problems with matrices that are nearly singular.
- For a \(n\times n\) unsymmetric matrix \(A\), its LU decomposition \(A = LU\) is computed where \(L\) is a unit lower triangular matrix, and \(U\) is an upper triangular matrix.
- For an \(n\times n\) symmetric indefinite matrix, its \(LDL^T\) decomposition \(A = LDL^T\) is computed where \(D\) is a block diagonal matrix (with blocks of order 1 or 2), and \(L\) is a unit lower triangular matrix.
- For an \(m\times n\) rectangular matrix, with \(m\ge n\), its QR decomposition \(A=QR\) is computed. Here, \(Q\) is an \(m\times m\) orthogonal matrix, and \(R\) is an \(m\times n\) upper trapezoidal matrix.
The linear systems with \(L\), \(U\), \(D\), and \(Q\), as defined above, are easy to solve. Those with \(L\) are solved by forward substitution, those with \(U\) are solved by back substitution, those with \(D\) are solved by multiplying by the explicit inverse of the \(1\times 1\) and \(2\times 2\) blocks, and those with \(Q\) are solved by multiplying the right-hand side by \(Q^T\).
The factorizations discussed above and their use in solving the systems are mathematically equivalent to their dense counterparts. The added difficulty in the sparse case is that the factorization should try to avoid operating on the zeros of the matrix and should keep the factors as sparse as possible. Meeting these requirements usually entails first performing a symbolic analysis on the matrix in order to predict and reduce the memory and run time requirements for the subsequent numerical factorization.
The symbolic analysis can usually be viewed using the graph models for sparse matrices. Below some tasks to be performed during the symbolic analysis are discussed.
Cholesky factorization
Consider the sparse symmetric positive definite matrix
\(A = \left [ \begin{array}{rrrr} 4 & -1 & -1 & -1 \\ -1 & 2 & & \\ -1 & & 2 & \\ -1 & & & 2\end{array}\right] \)
and its Cholesky factor \(L_A\) is given by (up to four units of accuracy)
\(L_A = \left [ \begin{array}{rrrr} 2.0000 & & &\\ -0.5000 & 1.3229 & & \\ -0.5000 & -0.1890 & 1.3093 & \\ -0.5000 & -0.1890 & -0.2182 & 1.2910\end{array}\right ]. \)
As is seen, the lower triangular matrix \(L_A\) is full. The entries of \(L_A\) that were zero in \(A\) are called fill-in. In the matrix \(L_A\), the \((3,2)\), \((4,2)\), and \((4,3)\) entries are fill-ins.
Consider the permutation matrix
\(P = \left [ \begin{array}{rrrr} 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1\\ 1 & 0 & 0 & 0\end{array}\right ] . \)
If \(A\) is multiplied by \(P\) from the left and by \(P^T\) from the right, the symmetric positive definite matrix \(B = PAP^T\) is obtained, viz.
\(B = \left [\begin{array}{rrrr} 2 & & & -1 \\ & 2 & & -1 \\ & & 2 & -1\\ -1 & -1 & -1 & 4\end{array}\right ] \)
with its Cholesky factor \(L_B\)
\(L_B = \left [ \begin{array}{rrrr} 1.4142 & & &\\ & 1.4142 & & \\ & & 1.4142 &\\ -0.7071 & -0.7071 & -0.7071 & 1.5811\end{array}\right ] \)
having no fill-ins.
The factorization of \(B = L_BL_B^T\) is clearly preferable to that of \(A = L_AL_A^T\). The solution \(x\) of \(Ax=b\) can be computed by first finding the solution \(y\) of the linear system \(By=c\) with \(c = Pb\), and then by setting \(x=P^Ty\). Since there are fewer entries in \(L_B\) than in \(L_A\), the time and memory requirements will be less for the second alternative. Matrices of the form of \(A\) are called arrowhead matrices and clearly if the dimension is very large there will be a very significant difference in costs between working with the form of \(A\) (first row and column dense) and working with the form of \(B\) (last row and column dense).
Methods used to find a permutation matrix \(P\) to reduce fill-in so that a factorization of the matrix \(PAP^T\) has much less fill-in and requires fewer operations during the Cholesky factorization are called ordering methods. Ordering methods normally employ graph theoretical tools. There is a correspondence between Cholesky factorization and an elimination process in the undirected graph associated with \(A\). In this model, there is a vertex for each row/column of \(A\), and there is an edge for each off-diagonal nonzero entry \(a_{ij}=a_{ji}\neq 0\). One step of the factorization, where a pivot \(a'_{ii}\), for \(i=1, ..., n-1\) is used to reduce all entries below the \(i\)th diagonal entry to zero, corresponds to making the neighbours of vertex \(i\) a clique, and then removing the vertex \(i\) as described by Parter (1961) and Rose (1970) . Using this equivalence, ordering methods can be viewed as ordering the vertices of the graph \(G\) so that the same ordering determines the pivotal sequence for the matrix factorization. Minimizing the fill-in in Cholesky factorization has been shown to be NP-complete by Yannakakis (1981) . Fortunately, there are many heuristics to get an ordering for reducing the fill-in significantly, including the minimum degree (Tinney and Walker, 1967) and the nested dissection (George, 1973) heuristics and their variants.
Other symbolic analysis operations include obtaining main characteristics of the pattern of the Cholesky factor (such as the number of nonzeros for each row and column), and determining the pattern of the Cholesky factor using the elimination tree (Schreiber, 1982) .
LU decomposition
In the case of Cholesky factorization, the symbolic and numerical concerns are nicely separated. When performing an LU decomposition on the other hand, one may need to perform pivoting to maintain numerical stability. Although this means that the symbolic analysis cannot give as accurate a prediction as in the Cholesky case, there are effective symbolic analysis methods for the LU decomposition.
Consider the LU decomposition of \(A\) with partial pivoting, i.e., \(PA=LU\), where \(P\) is a permutation matrix describing the numerical pivoting. Consider a symbolic analysis for the Cholesky factorization on the matrix \(B=A^TA\), where \(B\) is ordered as \(Q^TBQ\) for fill reduction. The pattern of the Cholesky factor \(L_B\) of \(Q^TBQ\) will be a superset of the \(L\) factor of \(AQ\) for any permutation matrix \(P\). Therefore, the symbolic analysis can be performed on \(A^TA\) (either with or without forming the product). In general, the symbolic analysis on \(B\) will greatly overestimate the requirement for \(A\).
An alternative to the above approach is to perform a symbolic analysis on the symmetric matrix \(B= |A| + |A^T|\). This time, row pivoting operations during numerical factorization may greatly perturb the ordering found and hence the symbolic analysis will not be a good forecast of the work and storage involved in the subsequent factorization. If it is possible to choose the pivots from the diagonal during the numerical factorization then it should be possible to keep closer to the forecast of the symbolic phase. In order to encourage such an outcome, the matrix \(A\) can be permuted at the beginning of the symbolic analysis phase to have large diagonal entries, i.e., a matrix \(C =|AQ|+|AQ^T|\) is found whose diagonal entries are large. A permutation \(Q\) that puts large entries into the diagonal of a matrix can be found by using a variant of the maximum weighted bipartite matching algorithm on the bipartite graph model of \(A\) as shown by Duff and Koster (2001) .
\(LDL^T\) decomposition
In order to avoid taking square roots in Cholesky factorization, a symmetric positive definite matrix can be factored as \(A=LDL^T\) with \(D\) being diagonal and positive, and \(L\) with unit diagonal entries.
The same principles can be carried over to symmetric indefinite matrices by permitting \(D\) be a block diagonal matrix; in general with blocks of size \(1\times 1\) and \(2\times 2\). The blocks of \(D\) corresponds to the pivots. A useful approach is to determine those \(1\times 1\) and \(2\times 2\) pivots during the analysis phase with the hope that the predetermined pivots will be numerically favorable during the actual numerical factorization. In current methods, weighted bipartite matching based algorithms are used to determine those potential pivots. Once such pivots are determined, the fill reducing ordering methods are constrained to keep the \(2\times 2\) pivots together.
QR decomposition
Consider an \(m\times n\) sparse rectangular matrix \(A\) with a full column rank, having a QR factorization \(A=QR\). In a thin QR factorization \(R\) is written down as an upper triangular matrix, instead of an upper trapezoidal one. If \(R\) comes from such a factorization, then \(R\) is also the Cholesky factor of \(A^TA=R^TR\). Many of the sparsity oriented issues in the QR decomposition case uses this correspondence and harness the methods developed for sparse Cholesky factorization.
Some other sparsity issues
A common sparsity oriented technique is to permute a sparse matrix into block triangular (BTF) form using a matching of maximum cardinality in the bipartite graph of \(A\) as described by Pothen and Fan (1990) . A matrix can be put into BTF form, if it is reducible. Given a matrix in this form viz.
\( A_{BTF} = \left [\begin{array}{rrrr} A_{11}& A_{12} & \cdots& A_{1K}\\ & A_{22} & \cdots& A_{2K}\\ & &\ddots & \vdots \\ & & &A_{KK}\end{array}\right], \)
it is possible to exploit this to solve the linear system \(Ax=b\) without factorizing the whole matrix. One can factor the last block \(A_{KK}\) to find the corresponding entries in the unknown vector \(x\), and then substitute those in the equations involving \(A_{K-1,K-1}\) and \(A_{K-1,K}\) and so on to solve the linear system \(Ax=b\) using a block back substitution.
Computation schemes
The basic algorithms to compute \(LU\), \(LDL^T\) and Cholesky factorizations are called up-looking, left-looking, and right-looking methods. If \(A\) is very sparse and would have only little fill-in, the up-looking algorithm is preferable. Otherwise, supernodal (an efficient formulation of the left-looking method) and multifrontal (an efficient formulation of the right-looking method) methods are preferred.
The two computational schemes to compute QR factorization are based on Householder reflections and Givens rotations. There are left-looking supernodal, and right-looking, multifrontal methods to compute the QR factorization using Householder reflections. The Givens rotations follow the up-looking method and are preferable when there is only little fill-in.
- I. S. Duff and J. Koster, The design and use of algorithms for permuting large entries to the diagonal of sparse matrices, SIAM Journal on Matrix Analysis and Applications, 20(4), 889-901, 2001 ( doi ).
- A. George, Nested dissection of a regular finite element mesh, SIAM Journal on Numerical Analysis 10 (2): 345–363, 1973 ( doi ).
- S. Parter, The use of linear graphs in Gauss elimination, SIAM Review, 3(2), pp. 119-130, 1961 ( doi ).
- A. Pothen and C.-J. Fan, Computing the block triangular form of a sparse matrix, ACM Transactions on Mathematical Software, 16, pp. 303-324, 1990 ( doi ).
- D. J. Rose, Triangulated graphs and the elimination process, Journal of Mathematical Analysis and Applications, 32 (3), pp. 597-609, 1970 ( doi ).
- R. Schreiber, A New implementation of sparse Gaussian elimination, ACM Transactions on Mathematical Software, 8 (3), pp. 256--276, 1982 ( doi ).
- W. F. Tinney and J. W. Walker, Direct solutions of sparse network equations by optimally ordered triangular factorization, Proceedings of the IEEE, 55, pp. 1801–1809, 1967 ( doi ).
- M. Yannakakis, Computing the minimum fill-in is NP-complete, SIAM Journal on Algebraic and Discrete Methods, 2 (1), pp. 77-79, 1981 ( doi ).
Recommended reading
- T. A. Davis, Direct Methods for Sparse Linear Systems , SIAM, Philadelphia, PA, 2006.
- I. S. Duff, A. M. Erisman, and J. K. Reid, Direct Methods for Sparse Matrices, Oxford University Press, London, 1986.
- I. S. Duff and B. Uçar, Combinatorial problems in solving linear systems, in Combinatorial Scientific Computing, U. Naumann and O. Schenk, eds., CRC Press, Boca Raton, FL, 2012, pp. 21–68.
- A. George and J. W. H. Liu, Computer Solution of Large Sparse Positive Definite Systems, Prentice-Hall, Englewood Cliffs, N. J., 1981.
- J. W. H. Liu, The role of elimination trees in sparse factorization, SIAM Journal on Matrix Analysis and Applications, 11 (1), pp. 134--172, 1990 ( doi )
- Numerical Analysis
- Applied Mathematics
Personal tools
- View source
- View history
- Propose a new article
- Instructions for Authors
- Random article
Focal areas
- Astrophysics
- Celestial mechanics
- Computational neuroscience
- Computational intelligence
- Dynamical systems
- More topics
- Recently published articles
- Recently sponsored articles
- Recent changes
- All articles
- List all Curators
- List all users
- Scholarpedia Journal
- What links here
- Related changes
- Special pages
- Printable version
- Permanent link
- This page was last modified on 31 March 2016, at 07:39.
- This page has been accessed 43,696 times.
- "Direct methods for sparse matrix solution" by Iain Duff and Bora Uçar is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 3.0 Unported License . Permissions beyond the scope of this license are described in the Terms of Use
- Privacy policy
- About Scholarpedia
- Disclaimers
- MATLAB Answers
- Community Home
- File Exchange
- Communities
- Treasure Hunt
- Virtual Badges
- MATLAB FAQs
- Contributors
- Recent Activity
- Flagged Content
- Manage Spam
- Trial software
You are now following this question
- You will see updates in your followed content feed .
- You may receive emails, depending on your communication preferences .
How to solve a sparse matrix efficiently?
Direct link to this question
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently
5 Comments Show Hide 4 older comments
Direct link to this comment
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently#comment_200967
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently#comment_200968
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently#comment_201025
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently#comment_201034
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently#comment_201232
Sign in to comment.
Sign in to answer this question.
Answers (2)
Direct link to this answer
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently#answer_127682
0 Comments Show Hide -1 older comments
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently#answer_127683
6 Comments Show Hide 5 older comments
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently#comment_200966
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently#comment_200978
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently#comment_200987
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently#comment_201024
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently#comment_201068
https://www.mathworks.com/matlabcentral/answers/120850-how-to-solve-a-sparse-matrix-efficiently#comment_201153
- Thermal conduction (central differences)
- Enthalpy flow by mass flux (upwind scheme)
- Boundary conditions
- 3D geometry
- convergence
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
- Switzerland (English)
- Switzerland (Deutsch)
- Switzerland (Français)
- 中国 (English)
You can also select a web site from the following list:
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
- América Latina (Español)
- Canada (English)
- United States (English)
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 简体中文 Chinese
- 日本 Japanese (日本語)
- 한국 Korean (한국어)
Contact your local office
Sparse linear algebra ( scipy.sparse.linalg ) #
Abstract linear operators #, matrix operations #, matrix norms #, solving linear problems #.
Direct methods for linear equation systems:
Iterative methods for linear equation systems:
Iterative methods for least-squares problems:
Matrix factorizations #
Eigenvalue problems:
Singular values problems:
The svds function supports the following solvers:
- svds(solver=’arpack’)
- svds(solver=’lobpcg’)
- svds(solver=’propack’)
Complete or incomplete LU factorizations
Exceptions #
IMAGES
VIDEO
COMMENTS
The six steps of problem solving involve problem definition, problem analysis, developing possible solutions, selecting a solution, implementing the solution and evaluating the outcome. Problem solving models are used to address issues that...
When multiplying or dividing different bases with the same exponent, combine the bases, and keep the exponent the same. For example, X raised to the third power times Y raised to the third power becomes the product of X times Y raised to th...
The four steps for solving an equation include the combination of like terms, the isolation of terms containing variables, the isolation of the variable and the substitution of the answer into the original equation to check the answer.
Data Structure for Sparse Matrix Reordering (cont.) ... (4) Solve (D-1 - H)z = -y.
One can factor the last block A_{KK} to find the corresponding entries in the unknown vector x, and then substitute those in the equations
Solving Linear Systems: Sparse Matrices, Iterative Methods ... Fill-in is a major problem for certain sparse matrices and leads to.
Problem definition. • Input. Problem definition. – Matrix A ... Takes advantage of sparse structure.
where A is an m x n matrix, b is given m-vector, and x is unknown solution n-vector to be determined. Such a system of equations asks question “Can vector b be.
Richard Peng (Georgia Tech)https://kyng.inf.ethz.ch/acseminar/2020-10-22_peng.htmlOctober 22, 2020.
The first of a series of 42 lectures on direct methods for sparse ... The Numerics of Solving Sparse Linear Systems Faster than Matrix
The solution to representing and working with sparse matrices is to use an alternate data structure to represent the sparse data. The zero
schemes for sparse matrices and some simple linear algebra operations using them.
As the upwind scheme is somewhat asymmetric, the sparse Matrix that looks like the one below. So far I use x = A\b, but when the problem starts to get
Compute a lower bound of the 1-norm of a sparse matrix. Solving linear problems#. Direct methods for linear equation systems: spsolve (A