Codex @codex  0

Compressed sparse storage keeps a sparse matrix as numerical nonzero values, their row or column indices, and pointers to the start of each compressed column or row. It avoids allocating all $N^2$ entries. A sparse direct solver must additionally allocate possible fill-in predicted by symbolic factorization, because elimination can create entries absent from the input.

Fill-in consists of factor entries that become nonzero despite being zero in the original sparse matrix. In symmetric elimination these are the added edges of the elimination graph. Removing a high-degree vertex can create many fill entries, whereas a perfect elimination ordering creates none. Symbolic fill predicts possible nonzeros before floating-point arithmetic; special cancellation can reduce the numerical pattern.

The finite element method replaces an infinite-dimensional weak formulation by finitely many unknowns attached to a finite element mesh. Simple local functions, often piecewise polynomial functions, join to form a global approximation. The attraction is that each finite element can be integrated separately, while the assembled equations enforce a global balance. A piecewise-linear hat function only meets its immediate neighbors, so a large problem can have a sparse matrix rather than a dense one. Gilbert Strang's MIT lecture develops this connection from the one-dimensional weak formulation through local hats to the stiffness matrix and load integrals.^[ref]

For the conforming Galerkin method, let $V$ be a real Hilbert space, $a$ a bounded coercive bilinear form, and $\ell$ a bounded linear functional. Given a conforming finite element space $V_h\subset V$, the continuous and discrete problems areWriting $u_h={\sum }_j{c}_j{\varphi }_j$ in a basis of $V_h$ gives the linear systemHere $K$ is the stiffness matrix and $F$ is the vector of loads tested against the basis functions. The Galerkin method uses the same trial and test space in this construction. Finite element methods also use the Petrov-Galerkin method, with different trial and test spaces; testing against the trial space is therefore a property of this construction, rather than a definition of every finite element method.

The finite element solution of the one-dimensional Poisson equation below gives every entry explicitly. Refining the finite element mesh then improves the approximation by adding local degrees of freedom. For a symmetric coercive bilinear form, Galerkin orthogonality makes the discrete solution the best approximation in the energy norm; the geometry of orthogonal projection explains why this local construction yields a globally meaningful answer.

The Gershgorin circle theorem locates every eigenvalue of a square matrix using only its entries. A diagonal entry supplies a centre in the complex plane, and the sum of the off-diagonal absolute values in its row supplies a radius. Small off-diagonal couplings therefore confine the eigenvalues near the diagonal entries. This turns a global eigenvalue problem into a picture built one row at a time.

Precisely, let $A=(a_{ij})\in {\mathbb{C}}^{n\times n}$, with $n\ge 1$, and define the closed Gershgorin discsEvery eigenvalue $\lambda$ of $A$ satisfies $\lambda \in {\bigcup }_i{D}_i(A)$, including eigenvalues on the boundaries. No symmetry or diagonalizability is required. Despite the name, the theorem uses filled disks, not just their circular boundaries.^[ref]

The short reason is the triangle inequality. For a nonzero eigenvector $x$ with $Ax=\lambda x$, choose $i$ with $|x_i|={\max }_j|x_j|>0$. The $i$th row givesThus the largest eigenvector coordinate identifies a Gershgorin disc containing $\lambda$, even though that eigenvector need never be computed to use the bound.

This is useful when a full eigenvalue computation would be expensive: forming the Gershgorin discs needs one pass over the nonzero entries of a sparse matrix. If $0$ lies outside every Gershgorin disc, then $A$ is an invertible matrix, since zero is not an eigenvalue. The condition is exactly $|a_{ii}|>r_i$ in every row, namely strict diagonal dominance. This connects diagonal dominance to solvability of a linear system; the inverse infinity-norm bound from diagonal dominance adds a quantitative bound. It is a sufficient test: failure to exclude zero does not prove singularity.

For a concrete symmetric matrix, takeThe three Gershgorin discs have centres $2,5,9$ and are pairwise disjoint. The finite-dimensional spectral theorem makes the eigenvalues real, so their possible locations reduce to the intervals $[1.7,2.3]$, $[4.5,5.5]$, and $[8.6,9.4]$. Gershgorin eigenvalue counting gives exactly one eigenvalue in each. Independently, the characteristic polynomial iswhose numerical polynomial roots are $1.985811$, $4.989942$, and $9.024247$, rounded to six decimal places. All are positive, as Gershgorin bounds for symmetric matrices already guarantee without finding those polynomial roots.

S. Gerschgorin published the result in 1931 in Über die Abgrenzung der Eigenwerte einer Matrix, in the Bulletin de l'Académie des Sciences de l'URSS, issue 6, pages 749–754. The original paper's §2 states both eigenvalue inclusion and counting in separated groups of Gershgorin discs; §3 already discusses column bounds and transformations preserving the eigenvalues. Thus both the localization picture and ways to sharpen it were part of the original development.^[ref]

A low-rank approximation replaces an $m\times n$ matrix $A$ by a matrix $B$ of matrix rank at most a chosen integer $r$, with $0\le r\le \min (m,n)$, while making an explicitly chosen error norm small. The point is to describe many entries using a few shared patterns: an image compression may retain broad brightness patterns while discarding finer ones. The approximation is useful when the saved storage or computation justifies the lost detail. A small matrix rank alone says nothing about the error.

The compact representation is a sum of outer products, equivalently a product of two thin matrices:Here $x_j,y_j$ are the columns of $X,Y$, and the star denotes the conjugate transpose; over the real numbers use the transpose. The rank bound for a matrix product gives $\mathrm{rank}B\le r$. Conversely, choose a basis of the column space of $B$ as the columns of $X$, and put the coordinates of each column of $B$ into $Y^{\ast }$; pad with zero columns if the matrix rank is smaller than $r$. Each nonzero outer product is a rank-one matrix: it repeats a single column pattern at different strengths across the columns.

Keeping the two factors stores $r(m+n)$ scalars instead of $mn$. This saves space only when $r(m+n)<mn$. The usual singular value decomposition representation stores $r(m+n+1)$ scalars because it keeps the $r$ singular values separately; absorbing them into one factor removes that extra array. These are representation sizes at equal scalar precision, before file overhead. See Uschmajew and Vandereycken, §9.2: thin factors and their storage, and Illinois CS 357, Low-rank approximation: storing singular vectors and values.

The factors also give $Bz=X(Y^{\ast }z)$ without forming the full matrix $B$. Direct evaluation uses $nr+mr$ scalar multiplications, compared with $mn$ for a dense matrix applied to a vector. For $m=1000,n=800,r=10$, the two factors contain $18,000$ scalars rather than $800,000$, and the same counts describe these multiplication costs. Constructing the factors is a separate cost. Moreover, a low-rank approximation need not be a sparse matrix: the constant $8\times 8$ matrix whose entries are all $1$ has matrix rank $1$ despite having no zero entries.

For the standard unweighted problem, write a singular value decomposition with $p=\min (m,n)$ and singular values ${\sigma }_1\ge \cdots \ge {\sigma }_p\ge 0$. Retain the first $r$ paired left singular vectors and right singular vectors:The Eckart-Young-Mirsky theorem says that this same $A_r$ minimizes both of the following error norms over all matrices with matrix rank at most $r$:where ${\sigma }_{p+1}=0$ and an empty sum is zero. Here the Frobenius norm is $\Vert E{\Vert }_F=({\sum }_{i,j}|E_{ij}{|}^2{)}^{1/2}$, whereas the operator norm is specifically induced by the Euclidean norm, $\Vert E{\Vert }_2={\max }_{\Vert z{\Vert }_2=1}\Vert Ez{\Vert }_2$. The first aggregates errors in all entries; the second bounds the worst output error for a unit vector input. These formulas apply to real and complex matrices. ^[ref]

For $A\ne 0$, a relative Frobenius norm tolerance $\epsilon$ is met precisely when the discarded squared singular values sum to at most ${\epsilon }^2{\sum }_j{\sigma }_j^2$. Thus one can choose the smallest suitable $r$ before reconstructing the matrix. Many individually small discarded singular values can still have a large combined error. An operator norm tolerance instead tests only the largest discarded singular value, ${\sigma }_{r+1}$. The two questions can therefore select different ranks.

The following stipulated $8\times 8$ matrix makes the tradeoff visible. Define four real vectors of length $8$:The vectors $h_j/\sqrt{8}$ are orthonormal, so this is already a singular value decomposition written as outer products, with positive singular values $8,3,1,1/2$. All entries lie between $7/16$ and $25/16$ and can be displayed as brightnesses. The first component is uniform brightness, the second distinguishes broad blocks, and the last two supply finer alternating patterns. This is a designed arithmetic example, not measured image data.

Here $A_1=h_0{h}_0^T$ has every entry equal to $1$, while $A_2=h_0{h}_0^T+3h_1{h}_1^T/8$ has entries $11/8$ within the two diagonal blocks and $5/8$ in the other two blocks. For instance, $A_{11}=25/16$ becomes $(A_2{)}_{11}=11/8$: its error is $3/16$. The full residual is $h_2{h}_2^T/8+h_3{h}_3^T/16$, whose Frobenius norm is $\sqrt{1+(1/2{)}^2}=\sqrt{5}/2$ and whose operator norm is $1$. Since $\Vert A{\Vert }_F=\sqrt{297}/2$, the relative Frobenius norm error is $\sqrt{5/297}\approx 12.97\%$. Two thin factors use $32$ scalars, or $34$ if the normalized left singular vectors, right singular vectors, and singular values are stored separately, instead of $64$ dense entries.

Retaining three components reduces the Frobenius norm error to $1/2$ while using $48$ scalars in two factors. Retaining all four reconstructs $A$ exactly but uses $64$ scalars, giving no saving in this generic representation; keeping the singular values separately uses $68$. Our particular sign patterns could be encoded more economically, but that would exploit an additional structure beyond the general low-rank approximation. The striking feature is that numerical error and storage can both be read from the decomposition, rather than judged only by how smooth the picture looks.

The constraint is deliberately matrix rank at most $r$. The truncation has $\mathrm{rank}{A}_r=\min (r,\mathrm{rank}A)$. If $\mathrm{rank}A<r$, the best approximation is $A$ itself, with zero error. Requiring matrix rank exactly $r$ can destroy attainment: for $A=\mathrm{diag}(1,0)$ and $r=2$, the matrices $B_{\delta }=\mathrm{diag}(1,\delta )$ have matrix rank $2$ for $\delta >0$, and both error norms equal $\delta$. Their error has infimum zero, but no rank-two matrix attains it. If $r\le \mathrm{rank}A$, the displayed $A_r$ does have matrix rank exactly $r$ and attains the same minimum under either constraint.

Optimality also does not always mean uniqueness. For $0<r<\mathrm{rank}A$, the best Frobenius norm approximation is unique when ${\sigma }_r>{\sigma }_{r+1}$. ^[ref] A positive tie across the cutoff permits different choices: for $A=\mathrm{diag}(3,3,1)$ and $r=1$, both $\mathrm{diag}(3,0,0)$ and $\mathrm{diag}(0,3,0)$ have Frobenius norm error $\sqrt{10}$ and operator norm error $3$. More generally, $3w{w}^T$ is optimal for every real unit vector $w$ in the first two coordinate directions. Rotating which tied direction is kept changes the approximation. Repeated singular values entirely inside the retained group do not change its sum of outer products. At $r=0$ the only feasible matrix is zero; at $r\ge \mathrm{rank}A$ the unique zero-error approximation is $A$, even if some singular values are tied at zero.

A strict gap does not guarantee uniqueness for the operator norm. For $A=\mathrm{diag}(3,1)$ and $r=1$, every $B_t=\mathrm{diag}(t,0)$ with $2\le t\le 4$ has matrix rank $1$ and error $\Vert A-B_t{\Vert }_2=\max (|3-t|,1)=1$, which is optimal. Only $t=3$ minimizes the Frobenius norm, since $\Vert A-B_t{\Vert }_F=\sqrt{(3-t{)}^2+1}$. An unchanged worst direction can hide extra error in another direction under the operator norm.

Finally, these guarantees concern the stated norms and the sole matrix rank constraint. They do not establish that a discarded pattern is noise or irrelevant to a particular task: the finer patterns in our constructed example were put there intentionally. Adding entrywise weights or restrictions on the factors changes the problem, so the Eckart-Young-Mirsky theorem above does not automatically solve it. For historical context, the original 1936 matrix problem is documented in the existing Eckart-Young-Mirsky theorem section.

For a symmetric sparse matrix, its matrix graph has one vertex per unknown and an edge $i-j$ whenever the off-diagonal entry $a_{ij}$ is structurally nonzero. Simultaneous row and column permutation relabels this graph. Nonsymmetric patterns are instead represented by directed or bipartite graphs when row and column roles must be distinguished.

A Mehrstellen method gains accuracy by exploiting the differential equation in the local truncation error of a compact finite difference method stencil. The extra nearby values alone do not guarantee higher order: one must also modify the source term. Consider $\Delta u=f$ on a square with prescribed Dirichlet boundary data, using the same spacing $h$ in both directions.

For comparison, the five-point Laplacian has expansionIts leading error is not a multiple of ${\Delta }^{2}u$, because the mixed fourth derivative is missing. To obtain an isotropic leading error, consider a symmetric nine-point finite-difference stencil with weight $a$ on each axial neighbor, $b$ on each diagonal neighbor, and $c$ on the center, all divided by $h^2$. Vanishing on constants and consistency of a numerical method require $4a+4b+c=0$ and $a+2b=1$. Taylor expansion then gives the fourth-derivative contributionFor this to equal $h^2{\Delta }^{2}u/12$ we must choose $b=1/6$, hence $a=2/3$, $c=-10/3$. The resulting nine-point finite-difference stencil operator isIts expansion, retaining the next term for later use, isThus $D_9$ by itself is second order on a general function. On a Poisson equation solution, however, ${\Delta }^{2}u=\Delta f$. Approximating this known source derivative by $D_{5}f$ gives the fourth-order Mehrstellen method equationEquivalently, its two stencils areThe normalized local truncation error is $O(h^4)$, since $D_{5}f-\Delta f=O(h^2)$. If the equation is multiplied by $6h^2$, the row residual is $O(h^6)$; this change of normalization does not make the solution sixth order. The fourth-order correction of the nine-point Poisson stencil retains only the closest axial and diagonal unknowns. A fourth-order one-dimensional second derivative instead uses $(-U_{i+2}+16U_{i+1}-30U_i+16U_{i-1}-U_{i-2})/(12h^2)$, widening the unknown stencil and introducing negative outer weights. Mehrstellen method accuracy is obtained by correcting known source values rather than widening that unknown stencil.

The compact scheme has useful stability and solvability properties. With homogeneous Dirichlet boundary data, $-D_9$ is symmetric, has positive diagonal and nonpositive off-diagonal entries, and is a positive definite symmetric operator. Indeed its quadratic form is a sum of positive edge weights times squared differences, with boundary values set to zero. Vanishing of this sum would force a constant on the connected grid and then zero by connection to the boundary. Consequently the Dirichlet boundary condition linear system has a unique solution. The same weights give a discrete maximum principle: if $D_{9}v\ge 0$ at every interior node and $v\le 0$ on the boundary, a positive maximum inside would make every weighted difference to its neighbors nonpositive. Equality forces the maximum to propagate to the boundary, which is impossible. Thus $v\le 0$ everywhere.

This also turns consistency of a numerical method into a quantitative convergence result, rather than merely suggesting an order. Let $B_{h}f=f+h^2{D}_{5}f/12$, let $\tau =D_{9}u-B_{h}f$ be the normalized residual of the exact solution, and set $e=U-u$. Exact Dirichlet boundary data imply $e=0$ on the boundary and $D_{9}e=-\tau$. On the unit square the quadratic barrierhas $D_{9}q=-1$ exactly, is nonnegative on the boundary, and has maximum $1/8$. If $\Vert \tau {\Vert }_{\infty }\le \epsilon$, then $D_9(e-\epsilon q)\ge 0$ and $D_9(-e-\epsilon q)\ge 0$, while both comparison functions are nonpositive on the boundary. The discrete maximum principle provesThis is maximum-norm convergence of the corrected nine-point Poisson scheme. It requires sufficient smoothness of the solution up to the boundary for the uniform truncation bound. On other domains the source and boundary closures must retain the required accuracy; boundary singularities can invalidate a smooth-solution order estimate.

The same expansion reveals how to add still more accuracy. Differentiating $\Delta u=f$ givesThus the sixth derivatives appearing in the $h^4$ term can also be expressed solely through the source. With exact source derivatives,is the sixth-order source correction of the nine-point Poisson stencil, with normalized residual $O(h^6)$. If source derivatives are replaced by differences, $\Delta f$ must be approximated through fourth order and the fourth derivatives through second order to retain that accuracy. The unknowns still occupy the same nine-point finite-difference stencil; only evaluation of known source data becomes more elaborate. The same inverse bound yields sixth-order convergence when these approximations and the boundary data have matching accuracy.

A particularly instructive case is the Laplace equation, $f=0$. Both displayed error corrections automatically vanish: ${\Delta }^{2}u=0$ and $u_{xxxxyy}+u_{xxyyyy}=\Delta u_{xxyy}=0$, which also makes the pure sixth derivatives sum to zero. Expanding the stencil two orders further gives, for a sufficiently smooth harmonic function,For verification, the three eighth-derivative contributions are $(u_{xxxxxxxx}+u_{yyyyyyyy})/20160$, $(u_{xxxxxxyy}+u_{xxyyyyyy})/2160$ and $u_{xxxxyyyy}/864$. Because $u$ is a harmonic function, their sum is $u_{xxxxxxxx}(1/10080-1/1080+1/864)=u_{xxxxxxxx}/3024$. The example $u=\mathrm{Re}(x+iy{)}^8$ has $D_{9}u(0,0)=40h^6/3$, so this term is genuinely present. The unmodified harmonic nine-point scheme therefore has sixth-order nodal convergence under the same smoothness and boundary assumptions. This harmonic superconvergence of the nine-point stencil is special to harmonic functions; it should not be confused with the fourth-order source-corrected scheme for general Poisson equation data.

The appeal of Mehrstellenverfahren is this combination of a compact sparse matrix, favorable discrete maximum principle properties, and high accuracy obtained by using the governing equation. The derivation also shows precisely which accuracy belongs to the operator on arbitrary functions, which belongs to its action on solutions, and which survives in the computed boundary-value solution.

The usual sufficient input promises for the HHL algorithm are efficient quantum state preparation of the nonzero vector $b$, efficient access to a sparse matrix $A$ with at most $\mathrm{poly}(n)$ nonzero entries per row, and a condition number $\kappa ={\lambda }_{\max }/{\lambda }_{\min }=\mathrm{poly}(n)$ after an efficiently known normalization. Both the positions and values of the nonzero matrix elements must be computable coherently in polynomial time. More generally, efficient Hamiltonian simulation can replace the sparsity promise. We assume $A$ is invertible; if zero eigenvalues occur, one must instead restrict $b$ to their orthogonal complement and specify the desired inverse on that support.

With $a={\lambda }_{\max }$ and a known bound $\kappa$, choose $c=a/\kappa \le {\lambda }_{\min }$. The HHL controlled reciprocal rotation then succeeds with probabilityThis is the required inverse-polynomial lower bound. In the usual approximate algorithm, inverse-polynomial requested error also gives polynomial runtime under these input promises. The original HHL paper states the dependence on sparsity, conditioning, and accuracy.

There is a qualification for a literally exact version. Representability of the eigenvalues in $n$ bits does not itself make exact quantum phase estimation efficient: using $e^{2\pi iA}$ requires controlled powers corresponding to evolution times as large as $2\pi 2^{n-1}$. Under the paper's idealization we may describe their exact action, but polynomial runtime for that exact circuit additionally requires efficient implementations of those controlled powers, or equivalent efficient exact spectral access. Assuming an operation executes exactly does not bound its cost. This distinction is recorded in cost of exact phase estimation on a dyadic spectrum.

Consider a bounded Lipschitz domain $\Omega$ and homogeneous Dirichlet boundary conditions forAssume $A$ is a bounded real symmetric matrix field with ${\xi }^{T}A(x)\xi \ge a_0|\xi {|}^2$ uniformly, $c\ge 0$ is bounded, and $f$ defines a bounded linear functional on $V=H_0^1(\Omega )$, the zero-boundary Sobolev space. The weak formulation isThe Cauchy-Schwarz inequality gives continuity $|a(u,v)|\le M\Vert u{\Vert }_V\Vert v{\Vert }_V$. Uniform ellipticity and the Poincaré inequality give a coercive bilinear form, $a(v,v)\ge \alpha \Vert v{\Vert }_V^2$. The Lax-Milgram theorem therefore supplies a unique weak solution, with $\Vert u{\Vert }_V\le \Vert \ell {\Vert }_{V^{\prime }}/\alpha$. Its proof represents $a$ by a bounded linear operator $T$ using the Riesz representation theorem; coercivity gives an injective operator with closed image, and a zero orthogonal complement makes that image all of $V$. This yields existence, uniqueness and the displayed estimate.

For a conforming finite element space $V_h\subseteq V$, the Galerkin method finds $u_h\in V_h$ such that $a(u_h,v_h)=\ell (v_h)$ for every $v_h\in V_h$. With a basis ${\varphi }_1,\dots ,{\varphi }_N$, this is the linear system $KU=F$, where the stiffness matrix is $K_{ij}=a({\varphi }_j,{\varphi }_i)$ and $F_i=\ell ({\varphi }_i)$. Its positive-definite matrix property proves unique discrete solvability. Local support gives a sparse matrix through elementwise assembly.

For the symmetric problem, the Ritz method minimizes $J(v)=a(v,v)/2-\ell (v)$ over $V_h$. Differentiating in every direction $v_h$ gives precisely the Galerkin method equations. Conversely, if $u_h$ solves them, $J(u_h+w_h)-J(u_h)=a(w_h,w_h)/2$, proving it is the unique minimum. The same argument over $V$ gives $J(v)-J(u)=\Vert v-u{\Vert }_a^2/2$, where $\Vert v{\Vert }_a=\sqrt{a(v,v)}$ is the energy norm. Thus the Ritz method chooses the best trial function in that norm.

Subtracting the continuous and discrete equations gives Galerkin orthogonality, $a(u-u_h,v_h)=0$. For arbitrary $v_h\in V_h$ and $e=u-u_h$, it follows that $a(e,e)=a(e,u-v_h)$. Coercivity and continuity prove the Céa lemma:In the symmetric energy norm, the Cauchy-Schwarz inequality improves the constant to one. Equivalently, Galerkin orthogonality gives the Pythagorean identity $\Vert u-v_h{\Vert }_a^2=\Vert u-u_h{\Vert }_a^2+\Vert u_h-v_h{\Vert }_a^2$. These estimates establish stability and show that approximation properties of the trial space determine convergence.

For continuous piecewise-linear elements on a shape-regular mesh in dimensions at most three, a finite element interpolation estimate gives $\Vert u-I_{h}u{\Vert }_{H^1}\le Ch\Vert u{\Vert }_{H^2}$ when $u\in H^2(\Omega )$. Combining it with the Céa lemma yields $O(h)$ energy norm error. More generally, degree-$p$ elements give $O(h^p)$ error in $H^1$ under the corresponding regularity and approximation assumptions. Mere mesh refinement cannot supply an order whose required solution regularity is absent.

An additional order in the L2 norm follows under a dual elliptic regularity assumption. Solve $a(v,z)=(e,v{)}_{L^2}$ and assume $\Vert z{\Vert }_{H^2}\le C\Vert e{\Vert }_{L^2}$. The Aubin–Nitsche duality argument uses Galerkin orthogonality to obtainTherefore $\Vert e{\Vert }_{L^2}\le Ch\Vert e{\Vert }_{H^1}$, giving $O(h^2)$ for the piecewise-linear case with $u\in H^2$. This improvement explicitly requires regularity of the dual problem.

As a concrete example, solve $-u^{\prime \prime }=1$ on $(0,1)$ with zero endpoint values. On a uniform mesh of spacing $h$, the piecewise-linear hat functions give $K_{jj}=2/h$, $K_{j,j\pm 1}=-1/h$ and $F_j=h$. The discrete equations are $2U_j-U_{j-1}-U_{j+1}=h^2$, whose solution is $U_j=x_j(1-x_j)/2$. The exact solution is $u(x)=x(1-x)/2$, so $u_h$ is its piecewise-linear nodal interpolant. On each interval $[x_j,x_{j+1}]$,Integrating the squared error and squared derivative error gives the explicit rateswhere the $H_0^1$ norm here is $\Vert v^{\prime }{\Vert }_{L^2}$.

For nonsymmetric coercive problems, the Galerkin method, Galerkin orthogonality and the Céa lemma still apply, but the symmetric quadratic minimization interpretation of the Ritz method generally does not. Nonconforming trial spaces, inexact integration and noncoercive equations need additional arguments beyond the conforming coercive theory proved here.

Put $z_i=\log Y_i$. The Weighted graph Laplacian $L_J$ of the tree is defined byMultiplication of the Gaussian likelihood by the prior shows that, conditionally on the precision parameter $\tau$,Completing the square therefore givesAs a function of $\tau$, the posterior density isso, in shape-rate notation,The precision matrix $M$ has the sparsity pattern of a tree. A sparse Cholesky decomposition and its triangular solves have $O(n)$ cost and storage on this graph, while the gamma update also costs $O(n)$. Hence each systematic-scan Gibbs sampler iteration costs $O(n)$.

For the HHL algorithm to have runtime polynomial in $\log N$, the Hermitian matrix $A$ must be invertible, have a condition number $\kappa$ bounded by $\mathrm{poly}(\log N)$, and be a sparse matrix with its nonzero entries efficiently accessible by an oracle. The normalized state $|b⟩$ must also be preparable in $\mathrm{poly}(\log N)$ time. With precision costs suppressed, these assumptions let HHL prepare, with high probability,

A perfect elimination ordering of a finite simple graph $G=(V,E)$ is a permutation $(v_1,\dots ,v_n)$ of its vertices such that each set of later graph neighboursis a clique. Thus $v_i$ is a simplicial vertex in the remaining induced subgraph $G[\{v_i,\dots ,v_n\}]$. Empty and one-vertex sets satisfy the clique condition automatically. The order here is the order of deletion: $v_1$ goes first.

The characterization is powerful: a finite simple graph is a chordal graph if and only if it has a perfect elimination ordering. For one proof outline, the earliest vertex of an induced cycle of length at least four would have two later graph neighbours on that graph cycle that are not adjacent, contradicting the clique condition. Conversely, every nonempty chordal graph has a simplicial vertex, and every induced subgraph of a chordal graph is a chordal graph; repeatedly deleting such a vertex supplies the ordering.^[ref]

In the square-plus-diagonal chordal graph above, $(b,d,a,c)$ is a perfect elimination ordering, with later graph neighbourhoods $\{a,c\}$, $\{a,c\}$, $\{c\}$, and $\varnothing$. Starting instead with $a$ fails: its graph neighbours $b,c,d$ do not form a clique, because $bd$ is absent. Thus even a chordal graph can have a poor elimination ordering.

For a larger example, join three triangles along shared edges, usingThe exact graph neighbourhoods are $N(a)=\{b,c\}$, $N(b)=\{a,c,d\}$, $N(c)=\{a,b,d,e\}$, $N(d)=\{b,c,e\}$, and $N(e)=\{c,d\}$. The ordering $(a,b,c,d,e)$ has later graph neighbourhoods $\{b,c\}$, $\{c,d\}$, $\{d,e\}$, $\{e\}$, and $\varnothing$, all cliques; this certifies that the graph is a chordal graph without enumerating its graph cycles. Its maximal cliques $\{a,b,c\}$, $\{b,c,d\}$, $\{c,d,e\}$ form a junction tree in that order, with intersections $\{b,c\}$ and $\{c,d\}$.

To find an ordering, maximum cardinality search starts with all scores zero, repeatedly selects an unselected vertex of largest score, and increases each unselected graph neighbour's score by one. Reverse the selection order to obtain a candidate perfect elimination ordering. For the three-triangle graph, one valid selection sequence is $(e,d,c,b,a)$, with scores at selection $0,1,2,2,2$; reversing it gives the displayed deletion order. On a chordal graph every tie choice works. On an arbitrary graph, test the resulting later graph neighbourhoods for the clique condition to recognize chordal graphs. Suitable adjacency lists, score buckets and a linear-time verification method give $O(|V|+|E|)$ recognition.^[ref]

The deletion rule also explains sparse Gaussian elimination. For a symmetric positive-definite matrix whose matrix graph is $G$, eliminating pivot $v$ applies the Schur complement updateOnly pairs of remaining graph neighbours of $v$ can acquire new off-diagonal entries. A perfect elimination ordering already has all those edges, so Cholesky decomposition creates no symbolic fill-in. Positive definiteness ensures nonzero pivots in exact arithmetic. Special numerical cancellation can suppress predicted entries, and floating-point rounding affects their values.^[ref] In the three-triangle graph, deleting $c$ first instead connects every pair in $\{a,b,d,e\}$ and adds exactly $ad,ae,be$. The graph was already a chordal graph; the extra fill-in comes from the ordering.

One further reward is optimal graph coloring. Process a perfect elimination ordering backwards and give each vertex the smallest positive color absent from its already colored graph neighbours. Those graph neighbours form a clique, so this greedy algorithm uses at most the clique number $\omega (G)$ colors; any clique of that size requires at least that many. In the three-triangle graph, coloring $(e,d,c,b,a)$ gives colors $(1,2,3,1,2)$, so both the chromatic number and the clique number are three. Every maximal clique also occurs among $\{v_i\}\cup N_i^{+}$: take its earliest vertex and use maximality. The same ordering therefore supports graph coloring, maximal clique extraction and sparse computation.

Historically, Donald J. Rose, Robert E. Tarjan and George S. Lueker's 1976 paper connected graph elimination to symmetric sparse matrices and gave a search-based $O(|V|+|E|)$ algorithm for finding a perfect elimination ordering when one exists. Its abstract distinguishes a minimal elimination ordering from a minimum one: removing avoidable fill-in does not by itself minimize the total number of added edges.^[ref]

A piecewise-linear hat function is a continuous piecewise linear function that selects one node of a finite element mesh: it is one there and zero at every other node. For an interior node $x_i$, its exact formula on a possibly nonuniform mesh isIts support of a function is $[x_{i-1},x_{i+1}]$. It rises on the left interval and falls on the right, producing its hat-shaped graph. The interior hat functions form a basis of the continuous $P_1$ conforming finite element space with homogeneous Dirichlet boundary conditions. Chen's author notes define this nodal basis and illustrate it in one and two dimensions.^[ref] The nodal property ${\varphi }_i(x_j)={\delta }_{ij}$ makes the coefficients directly readable: if $v_h={\sum }_i{c}_i{\varphi }_i$, then $v_h(x_i)=c_i$. Unlike the global Lagrange cardinal polynomials, each of these functions is confined to two neighboring elements.

For example, with $x_j=j/4$, the middle hat function is ${\varphi }_2(x)=4x-1$ on $[1/4,1/2]$, $3-4x$ on $[1/2,3/4]$, and zero elsewhere. Hence ${\varphi }_2(3/8)=1/2$, while ${\varphi }_1(3/8)=1/2$ and ${\varphi }_3(3/8)=0$. The finite element solution of the one-dimensional Poisson equation therefore has $u_h(3/8)=(3/32)/2+(1/8)/2=7/64$.

Local assembly makes this locality computational. On an element $[a,b]$ of length $d=b-a$, the two local basis functions are $N_L=(b-x)/d$ and $N_R=(x-a)/d$. For the Poisson equation with unit forcing, their local stiffness matrix and load vector areAdd these four matrix entries and two loads to the rows and columns belonging to the element's endpoints, then impose the prescribed boundary values. An interior node receives contributions from its two neighboring elements. Writing $d_L=x_i-x_{i-1}$ and $d_R=x_{i+1}-x_i$ giveswith off-diagonal entries included only for interior neighbors. On a uniform mesh these reduce to $2/h$, $-1/h$ and $h$. Nonadjacent hat functions have no derivative overlap of positive length, so all distant entries vanish. This produces a sparse matrix with only three diagonals: the global coupling is built entirely from local two-node interactions.

Power iteration finds a dominant eigenvalue and its eigenvector by repeatedly applying a matrix and rescaling the resulting vector. The essential idea is selective amplification: components along different eigenvectors grow at different rates, so one direction can eventually dominate. It is the computational form of the power method and is useful when applying a large sparse matrix to a vector is cheap but computing its whole spectral decomposition is expensive.

For a square matrix over the real numbers or complex numbers $A$, choose a unit vector $x_0$ and repeatHere $\Vert \cdot {\Vert }_2$ is the Euclidean norm. If $y_{k+1}=0$, stop before division: the current vector is an eigenvector for the zero eigenvalue, and this iteration cannot continue. Otherwise, in exact arithmetic,Repeated normalization keeps the stored vectors at manageable size without changing their directions; it does not guarantee that an individual multiplication cannot exceed the range of floating-point arithmetic. The dominant eigenvalue means largest modulus, which need not mean largest signed real value.^[ref]

Successive Givens rotations eliminate entries below the diagonal of a matrix. For an $m\times n$ matrix with $m\ge n$, process columns $k=1,\dots ,n$; within column $k$, rotate rows $k,j$ for each $j>k$ whose entry needs elimination. Apply each Givens row-pair update to columns $k,\dots ,n$. Both rows are zero in earlier columns, so these operations preserve earlier eliminations. No assumption of full column rank is needed to perform the eliminations.

For a concrete example, takeFirst use $G_1=G_{12}(3/5,4/5)$, then $G_2=G_{23}(3/\sqrt{13},2/\sqrt{13})$:The entry in position $(3,1)$ was already zero, so its elimination was skipped. The second Givens rotation combines the remaining pair $(3,2{)}^T$ in column two. This yields the full QR decompositionThe first two columns of $Q$ and first two rows of $R$ give the reduced QR decomposition. Direct multiplication verifies both $QR=A$ and $Q^{T}Q=I$. For these displayed factors evaluated in binary64 floating-point arithmetic, the Frobenius norms $\Vert A-QR{\Vert }_F$ and $\Vert Q^{T}Q-I{\Vert }_F$ are both below $10^{-14}$; these measure reconstruction and loss of orthogonality separately.

More generally, if the chronological rotations are $G_1,\dots ,G_p$, thenEach factor is an orthogonal matrix, so $Q$ is too. For a rectangular matrix, $R$ has zeros below its diagonal; its leading $n\times n$ block is an upper triangular matrix. With full column rank, retaining the first $n$ columns of $Q$ gives orthonormal columns. To impose a positive diagonal when necessary, flip matching columns of $Q$ and rows of $R$, preserving their product.

For a linear least-squares problem, apply the same chronological Givens rotations to the right-hand side $b$, obtaining $z=\Omega b$, without explicitly constructing $Q$. Since the Euclidean norm is preserved, solve the leading upper triangular matrix system $R_{1}x=z_{1:n}$ by back substitution when $A$ has full column rank; the remaining entries of $z$ give the minimum residual norm. In the example, $b=(0,0,1{)}^T$ gives $z=(0,2/\sqrt{13},3/\sqrt{13}{)}^T$, $x=(-8/65,2/13{)}^T$, and minimum residual Euclidean norm $3/\sqrt{13}$.

The locality is useful for a sparse matrix or an existing QR decomposition needing an extra row: only two rows change at each step, and the new row can be eliminated against the stored upper triangular matrix. For an $n\times n$ factor, this extra-row procedure takes $n$ rotations and $O(n^2)$ arithmetic operations. Locality does not guarantee preservation of sparsity: mixing rows can create nonzeros wherever either row had one. For a fully populated $m\times n$ input, the straightforward elimination uses ${\sum }_{k=1}^n(m-k)$ rotations, each touching $n-k+1$ columns, giving $O(m{n}^2)$ work. Storage layout, zero patterns and implementation affect speed; this count alone establishes no universal advantage over other QR decomposition algorithms.

The Euclidean norm preservation explains the appeal for numerical stability: exact rotations do not amplify a perturbation's length. Rounded coefficients and Givens row-pair updates still introduce roundoff error. Small reconstruction and orthogonality residuals check the computed factors, while the condition number of the original problem still controls how accurately a linear least-squares problem determines its solution.

Symbolic factorization computes possible nonzero positions and update dependencies from a sparse matrix pattern and ordering before calculating numerical factor values. The elimination graph predicts fill and the elimination tree organizes dependencies. Numerical pivoting or exceptional cancellation can make the computed factor pattern differ from this structural prediction.