<jupyter_start><jupyter_text>Lecture 10:  Sparse direct solvers  Recap of the previous part- Distributed memory for huge dense matrices- Sparse matrix formats (COO, CSR, CSC)- Matrix-by-vector product- Inefficiency of sparse matrix processing- Approaches to reduce cache misses  Plan for todaySparse direct solvers: - LU decomposition of sparse matrices- fill-in of $L$ and $U$ factors- nested dissection- spectral clustering in details  LU decomposition of sparse matrix- Why sparse linear systems can be solved faster, what is the technique? - In the LU factorization of the matrix $A$ the factors $L$ and $U$ can be also sparse:$$A = L U$$- And solving linear systems with **sparse** triangular matrices is very easy.   Note that the inverse matrix is not sparse!<jupyter_code>import numpy as np
import scipy.sparse as spsp
n = 7
ex = np.ones(n);
a = spsp.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); 
b = np.array(np.linalg.inv(a.toarray()))
print(a.toarray())
print(b)
np.linalg.svd(b[:3, 4:])[1]<jupyter_output>[[-2.  1.  0.  0.  0.  0.  0.]
 [ 1. -2.  1.  0.  0.  0.  0.]
 [ 0.  1. -2.  1.  0.  0.  0.]
 [ 0.  0.  1. -2.  1.  0.  0.]
 [ 0.  0.  0.  1. -2.  1.  0.]
 [ 0.  0.  0.  0.  1. -2.  1.]
 [ 0.  0.  0.  0.  0.  1. -2.]]
[[-0.875 -0.75  -0.625 -0.5   -0.375 -0.25  -0.125]
 [-0.75  -1.5   -1.25  -1.    -0.75  -0.5   -0.25 ]
 [-0.625 -1.25  -1.875 -1.5   -1.125 -0.75  -0.375]
 [-0.5   -1.    -1.5   -2.    -1.5   -1.    -0.5  ]
 [-0.375 -0.75  -1.125 -1.5   -1.875 -1.25  -0.625]
 [-0.25  -0.5   -0.75  -1.    -1.25  -1.5   -0.75 ]
 [-0.125 -0.25  -0.375 -0.5   -0.625 -0.75  -0.875]]<jupyter_text>And the factors...$L$ and $U$ are typically sparse. In the tridiagonal case they are even bidiagonal!<jupyter_code>from scipy.sparse.linalg import splu
T = splu(a.tocsc(), permc_spec="NATURAL")
print(T.L.toarray())<jupyter_output>[[ 1.          0.          0.          0.          0.          0.
   0.        ]
 [-0.5         1.          0.          0.          0.          0.
   0.        ]
 [ 0.         -0.66666667  1.          0.          0.          0.
   0.        ]
 [ 0.          0.         -0.75        1.          0.          0.
   0.        ]
 [ 0.          0.          0.         -0.8         1.          0.
   0.        ]
 [ 0.          0.          0.          0.         -0.83333333  1.
   0.        ]
 [ 0.          0.          0.          0.          0.         -0.85714286
   1.        ]]<jupyter_text>Interesting to note that ```splu``` without ```permc_spec``` will produce permutations which will not yield the bidiagonal factor:<jupyter_code>from scipy.sparse.linalg import splu
T = splu(a.tocsc())
print(T.L.todense())
print(T.perm_c)<jupyter_output>[[ 1.          0.          0.          0.          0.          0.
   0.        ]
 [-0.5         1.          0.          0.          0.          0.
   0.        ]
 [ 0.         -0.66666667  1.          0.          0.          0.
   0.        ]
 [ 0.          0.         -0.75        1.          0.          0.
   0.        ]
 [ 0.          0.          0.          0.          1.          0.
   0.        ]
 [ 0.          0.          0.         -0.8        -0.5         1.
   0.        ]
 [ 0.          0.          0.          0.         -0.5        -0.71428571
   1.        ]]
[0 1 2 3 5 4 6]<jupyter_text>2D-caseFrom a matrix that comes as a discretization of  a two-dimensional problem everything is much worse:<jupyter_code>import scipy as sp
import scipy.sparse as spsp
import scipy.sparse.linalg
import matplotlib.pyplot as plt
%matplotlib inline
n = 20
ex = np.ones(n);
lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); 
e = sp.sparse.eye(n)
A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)
A = spsp.csc_matrix(A)
T = spsp.linalg.splu(A)
plt.spy(T.L, marker='.', color='k', markersize=8)
#plt.spy(A, marker='.', color='k', markersize=8)<jupyter_output><empty_output><jupyter_text>For correct permutation in 2D case the number of nonzeros in  $L$ factor grows as $\mathcal{O}(N \log N)$. But  complexity is $\mathcal{O}(N^{3/2})$.  Sparse matrices and graph ordering- The number of nonzeros in LU decomposition has a deep connection to the graph theory.- ``networkx package`` can be used to visualize graphs, given only the adjacency matrix. - It may even recover to some extend the graph structure.<jupyter_code>import networkx as nx
n = 10
ex = np.ones(n);
lp1 = sp.sparse.spdiags(np.vstack((ex,  -2*ex, ex)), [-1, 0, 1], n, n, 'csr'); 
e = sp.sparse.eye(n)
A = sp.sparse.kron(lp1, e) + sp.sparse.kron(e, lp1)
A = spsp.csc_matrix(A)
G = nx.Graph(A)
nx.draw(G, pos=nx.spectral_layout(G), node_size=50)<jupyter_output><empty_output><jupyter_text>Fill-inThe fill-in of a matrix are those entries which change from an initial zero to a nonzero value during the execution of an algorithm.The fill-in is different for different permutations. So, before factorization we need to find reordering which produces the smallest fill-in.**Example**$$A = \begin{bmatrix} * & * & * & * & *\\ * & * & 0 & 0 & 0 \\ * & 0  & * & 0 & 0 \\ * & 0 & 0& * & 0 \\ * & 0 & 0& 0 & * \end{bmatrix} $$If we eliminate elements from the top to the bottom, then we will obtain dense matrix.However, we could maintain sparsity if elimination was done from the bottom to the top.  Fill-in minimizationReordering the rows and the columns of the sparse matrix in order to reduce the number of nonzeros in $L$ and $U$ factors is called **fill-in** minimization.  Examples of algorithms for fill-in minimization:- **Minimum degree ordering** - order by the degree of the vertex- **Cuthill–McKee algorithm** (and reverse Cuthill-McKee) - reorder to minimize the bandwidth (does not exploit graph representation).- **Nested dissection**: split the graph into two with minimal number of vertices on the separator (set of vertices removed after we separate the graph into two distinct connected graphs).  Complexity of the algorithm depends on the size of the graph separator. For 1D Laplacian separator contains only 1 vertex, in 2D - $\sqrt{N}$ vertices.  Gaussian elimination for sparse matricesGiven matrix $A=A^*>0$ we calculate its Cholesky decomposition $A = LL^*$.Factor $L$ can be dense even if $A$ is sparse:$$\begin{bmatrix} * & * & * & * \\ * & * &  &  \\ * &  & * &  \\ * &  &  & * \end{bmatrix} = \begin{bmatrix} * &  &  &  \\ * & * &  &  \\ * & * & * &  \\ * & * & * & * \end{bmatrix}\begin{bmatrix} * & * & * & * \\  & * & * & * \\  &  & * & * \\  &  &  & * \end{bmatrix}$$How to make factors sparse, i.e. to minimize the **fill-in**?  Gaussian elimination and permutationWe need to find a permutation of indices so that factors are sparse, i.e. we build Cholesky factorisation of $PAP^\top$, where $P$ is a permutation matrix.For the example from the previous slide$$P \begin{bmatrix} * & * & * & * \\ * & * &  &  \\ * &  & * &  \\ * &  &  & * \end{bmatrix} P^\top = \begin{bmatrix} * &  &  & *  \\  & * &  & * \\  &  & * & * \\ * & * & * & * \end{bmatrix} = \begin{bmatrix} * &  &  &  \\  & * &  &  \\  &  & * &  \\ * & * & * & * \end{bmatrix}\begin{bmatrix} * &  &  & *  \\  & * &  & * \\  &  & * & * \\  & &  & * \end{bmatrix}$$where$$P = \begin{bmatrix}  &  &  & 1 \\  &  & 1 &  \\  & 1 &  &  \\ 1 &  &  &  \end{bmatrix}$$- Arrowhead form of the matrix gives sparse factors in LU decomposition<jupyter_code>import numpy as np
import scipy.sparse as spsp
import scipy.sparse.linalg as spsplin
import scipy.linalg as splin
import matplotlib.pyplot as plt
%matplotlib inline

A = spsp.coo_matrix((np.random.randn(10), ([0, 0, 0, 0, 1, 1, 2, 2, 3, 3], 
                                           [0, 1, 2, 3, 0, 1, 0, 2, 0, 3])))
print("Original matrix")
plt.spy(A)
plt.show()
lu = spsplin.splu(A.tocsc(), permc_spec="NATURAL")
print("L factor")
plt.spy(lu.L)
plt.show()
print("U factor")
plt.spy(lu.U)
plt.show()
print("Column permutation:", lu.perm_c)
print("Row permutation:", lu.perm_r)<jupyter_output>Original matrix<jupyter_text>Block arrowhead structure$$PAP^\top = \begin{bmatrix} A_{11} &  & A_{13} \\  & A_{22} & A_{23} \\ A_{31} & A_{32} & A_{33}\end{bmatrix}$$then$$PAP^\top = \begin{bmatrix} A_{11} & 0 & 0 \\ 0 & A_{22} & 0 \\ A_{31} & A_{32} & A_{33} - A_{31}A_{11}^{-1} A_{13} - A_{32}A_{22}^{-1}A_{23} \end{bmatrix} \begin{bmatrix}  I & 0 & A_{11}^{-1}A_{13} \\ 0 & I & A_{22}^{-1}A_{23} \\ 0 & 0 & I\end{bmatrix}$$- Block $ A_{33} - A_{31}A_{11}^{-1} A_{13} - A_{32}A_{22}^{-1}A_{23}$ is Schur complement for block diagonal matrix $\begin{bmatrix} A_{11} & 0 \\ 0 & A_{22} \end{bmatrix}$- We reduce problem to solving smaller linear systems with $A_{11}$ and $A_{22}$   How can we find permutation?- Key idea comes from graph theory- Sparse matrix can be treated as an **adjacency matrix** of a certain graph:the vertices $(i, j)$ are connected, if the corresponding matrix element is non-zero.  ExampleGraphs of $\begin{bmatrix} * & * & * & * \\ * & * &  &  \\ * &  & * &  \\ * &  &  & * \end{bmatrix}$ and $\begin{bmatrix} * &  &  & *  \\  & * &  & * \\  &  & * & * \\ * & * & * & * \end{bmatrix}$ have the following form: and * Why the second ordering is better than the first one?  Graph separator**Definition.** A **separator** in a graph $G$ is a set $S$ of vertices whose removal leaves atleast two connected components.Separator $S$ gives the following ordering for an $N$-vertex graph $G$:- Find a separator $S$, whose removal leaves connected components$T_1$, $T_2$, $\ldots$, $T_k$- Number the vertices of $S$ from $N − |S| + 1$ to $N$- Recursively, number the vertices of each component: $T_1$ from $1$ to$|T_1|$, $T_2$ from $|T_1| + 1$ to $|T_1| + |T_2|$, etc- If a component is small enough, enumeration in this component is arbitrarily  Separator and block arrowhead structure: exampleSeparator for the 2D Laplacian matrix $$ A_{2D} = I \otimes A_{1D} + A_{1D} \otimes I, \quad A_{1D} = \mathrm{tridiag}(-1, 2, -1),$$is as follows  Once we have enumerated first indices in $\alpha$, then in $\beta$ and separators indices in $\sigma$ we get the following matrix$$PAP^\top = \begin{bmatrix} A_{\alpha\alpha} &  & A_{\alpha\sigma} \\  & A_{\beta\beta} & A_{\beta\sigma} \\ A_{\sigma\alpha} & A_{\sigma\beta} & A_{\sigma\sigma}\end{bmatrix}$$which has arrowhrad structure.- Thus, the problem of finding **permutation** was reduced to the problem of finding **graph separator**!  Nested dissection- For blocks $A_{\alpha\alpha}$, $A_{\beta\beta}$ we continue splitting recursively.- When the recursion is done, we need to eliminate blocks $A_{\sigma\alpha}$ and $A_{\sigma\beta}$. - This makes block in the position of $A_{\sigma\sigma}\in\mathbb{R}^{n\times n}$ dense.Calculation of Cholesky of this block costs $\mathcal{O}(n^3) = \mathcal{O}(N^{3/2})$, where $N = n^2$ is the total number of nodes.So, the complexity is $\mathcal{O}(N^{3/2})$  Packages for nested dissection- MUltifrontal Massively Parallel sparse direct Solver ([MUMPS](http://mumps.enseeiht.fr/))- [Pardiso](https://www.pardiso-project.org/)- [Umfpack as part of SuiteSparse](http://faculty.cse.tamu.edu/davis/suitesparse.html)All of them have interfaces for C/C++, Fortran and Matlab   Nested dissection summary- Enumeration: find a separator.- Divide-and-conquer paradigm- Recursively process two subsets of vertices after separation- In theory, nested dissection gives optimal complexity. - In practice, it beats others only for very large problems.  Separators in practice- Computing separators is not a **trivial task**.- Graph partitioning heuristics has been an active research area for many years, often motivated by partitioning for parallel computation.Existing approaches:- Spectral partitioning (uses eigenvectors of Laplacian matrix of graph) - more details below- Geometric partitioning (for meshes with specified vertex coordinates) [review and analysis](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.31.4886&rep=rep1&type=pdf)- Iterative-swapping ([(Kernighan-Lin, 1970)](http://xilinx.asia/_hdl/4/eda.ee.ucla.edu/EE201A-04Spring/kl.pdf), [(Fiduccia-Matheysses, 19820](https://dl.acm.org/citation.cfm?id=809204))- Breadth-first search [(Lipton, Tarjan 1979)](http://www.cs.princeton.edu/courses/archive/fall06/cos528/handouts/sepplanar.pdf)- Multilevel recursive bisection (heuristic, currently most practical) ([review](https://people.csail.mit.edu/jshun/6886-s18/lectures/lecture13-1.pdf) and [paper](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.499.4130&rep=rep1&type=pdf)). Package for such kind of partitioning is called METIS, written in C, and available [here](http://glaros.dtc.umn.edu/gkhome/views/metis)  Spectral graph partitioningThe idea of spectral partitioning goes back to <NAME>, who studied connectivity of graphs ([paper](https://dml.cz/bitstream/handle/10338.dmlcz/101168/CzechMathJ_23-1973-2_11.pdf)).We need to split the vertices into two sets.Consider +1/-1 labeling of vertices and **the cost**$$E_c(x) = \sum_{j} \sum_{i \in N(j)} (x_i - x_j)^2, \quad N(j) \text{ denotes set of neighbours of a node } j. $$We need a balanced partition, thus $$\sum_i x_i =  0 \quad \Longleftrightarrow \quad x^\top e = 0, \quad e = \begin{bmatrix}1 & \dots & 1\end{bmatrix}^\top,$$and since we have +1/-1 labels, we have$$\sum_i x^2_i = n \quad \Longleftrightarrow \quad \|x\|_2^2 = n.$$  Graph LaplacianCost $E_c$ can be written as (check why)$$E_c = (Lx, x)$$where $L$ is the **graph Laplacian**, which is defined as a symmetric matrix with$$L_{ii} = \mbox{degree of node $i$},$$$$L_{ij} = -1, \quad \mbox{if $i \ne j$  and there is an edge},$$and $0$ otherwise.- Rows of $L$ sum to zero, thus there is an eigenvalue $0$ and gives trivial eigenvector of all ones.- Eigenvalues are non-negative (why?).  Partitioning as an optimization problemMinimization of $E_c$ with the mentioned constraints leads to a partitioning that tries to minimize number of edges in a separator, while keeping the partition balanced. We now relax the integer quadratic programming to the continuous quadratic programming$$E_c(x) = (Lx, x)\to \min_{\substack{x^\top e =0, \\ \|x\|_2^2 = n}}$$  Fiedler vector- The solution to the minimization problem is given by the eigenvector (called Fiedler vector) corresponding  to the second smallest eigenvalue of the graph Laplacian. Indeed,$$    \min_{\substack{x^\top e =0, \\ \|x\|_2^2 = n}} (Lx, x) = n \cdot \min_{{x^\top e =0}} \frac{(Lx, x)}{(x, x)} = n \cdot \min_{{x^\top e =0}} R(x), \quad R(x) \text{ is the Rayleigh quotient}$$- Since $e$ is the eigenvector, corresponding to the smallest eigenvalue, on the space $x^\top e =0$ we get the second minimal eigevalue.- The sign $x_i$ indicates the partitioning.- In computations, we need to find out, how to find this second minimal eigenvalue –– we at least know about power method, but it finds the largest. We will discuss iterative methods for eigenvalue problems later in our course.- This is the main goal of the iterative methods for large-scale linear problems, and can be achieved via few matrix-by-vector products.<jupyter_code>import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
import networkx as nx
kn = nx.read_gml('karate.gml')
print("Number of vertices = {}".format(kn.number_of_nodes()))
print("Number of edges = {}".format(kn.number_of_edges()))
nx.draw_networkx(kn, node_color="red") #Draw the graph
import scipy.sparse.linalg as spsplin
Laplacian = nx.laplacian_matrix(kn).asfptype()
plt.spy(Laplacian, markersize=5)
plt.title("Graph laplacian")
plt.axis("off")
plt.show()
eigval, eigvec = spsplin.eigsh(Laplacian, k=2, which="SM")
print("The 2 smallest eigenvalues =", eigval)
plt.scatter(np.arange(len(eigvec[:, 1])), np.sign(eigvec[:, 1]))
plt.show()
print("Sum of elements in Fiedler vector = {}".format(np.sum(eigvec[:, 1].real)))
nx.draw_networkx(kn, node_color=np.sign(eigvec[:, 1]))<jupyter_output><empty_output><jupyter_text>Summary on demo- Here we call SciPy sparse function to find fixed number of eigenvalues (and eigenvectors) that are smallest (other options are possible)- Details of the underlying method we will discuss soon- Fiedler vector gives simple separation of the graph- To separate graph on more than two parts you should use eigenvectors of laplacian as feature vectors and run some clustering algorithm, e.g. $k$-means  Fiedler vector and algebraic connectivity of a graph**Definition.** The algebraic connectivity of a graph is the second-smallest eigenvalue of the Laplacian matrix.**Claim.** The algebraic connectivity of a graph is greater than 0 if and only if a graph is a connected graph.  Practical problemsComputing bisection recursively is expensive. As an alternative, one typically computes **multilevel bisection** that consists of 3 phases.- Graph coarsening: From a given graph, we join vertices into larger nodes, and get sequences of graphs $G_1, \ldots, G_m$. - At the coarse level, we do high-quality bisection- Then, we do **uncoarsening**: we propagate the splitting from $G_k$ to $G_{k-1}$ and improve the quality of the split by local optimization algorithms (refinement).   Practical problems (2)Once the permutation has been computed, we need to implement the elimination, making use of efficient computational kernels. If in the elemination we will be able to get the elements into **blocks**, we will be able to use BLAS-3 computations. It is done by **supernodal** data structures: If adjacent rows have the same sparsity structure, they can be stored in **blocks**:Also, we can use such structure in efficient computations![Details in this survey](https://arxiv.org/pdf/1907.05309.pdf)  Minimal degree orderings- The idea is to eliminate rows and/or columns with fewer non-zeros, update fill-in and then repeat- Efficient implementation is an issue (adding/removing elements).- Current champion is "approximate minimal degree" by Amestoy, Davis, Duff.- It is **suboptimal** even for 2D problems- In practice, it often wins for medium-sized problems- SciPy sparse package [uses](https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.splu.html) minimal ordering approach for different matrices ($A^{\top}A$, $A + A^{\top}$)   Take home message- Sparse matrices & graphs ordering- Ordering is important for LU fill-in: more details in the next lecture- Separators and how do they help in fill-in minimization- Nested dissection idea- Fiedler vector and spectral bipartitioning- Other orderings from SciPy sparse package  Plan for the next part- Basic iterative methods for solving large linear systems- Convergence - Acceleration  Questions?<jupyter_code>from IPython.core.display import HTML
def css_styling():
    styles = open("./styles/custom.css", "r").read()
    return HTML(styles)
css_styling()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14307108092207305, 'fraction_numerical': 0.035460234262181096, 'mean_word_length': 3.105841018884497, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 31, 'https://': 6, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 2, 'substrings_counts': 3, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '322628', 'license': 'MIT', 'n_tokens_mistral': 6713, 'n_tokens_neox': 6141, 'n_words': 2174, 'path': 'lectures/lecture10/lecture-10.ipynb', 'repo_name': 'mikhailpautov/nla2020', 'size': '18707'}