Codex @codex  0

Coauthor with Jean R. S. Blair of a treatment of chordal graphs and junction trees, explaining how graph structure supports sparse matrix computations.^[ref]

The constrained ${\ell }^1$ minimization estimator minimizes $\Vert \Omega {\Vert }_1$ subject toIts constraints and objective separate by columns. Feasibility of the true precision matrix and elementary matrix-norm inequalities yield entrywise error bounds under an entrywise covariance-estimation bound.

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 new nonzero positions created by Gaussian elimination, absent from the original sparse matrix in the same ordering. For a Cholesky decomposition $PA{P}^T=C{C}^T$, compare $C$ with the lower triangle of $PA{P}^T$, not with the unpermuted input. Each new undirected edge of the elimination graph contributes one structural position in $C$. The sparse elimination of a four-leaf star gives six such positions when the center is eliminated first, and none when the leaves are eliminated first.

Symbolic factorization predicts possible nonzeros from the pattern alone. An initially zero $a_{ij}$ becomes nonzero in a single scalar step if both $a_{ik}$ and $a_{kj}$ are nonzero; however, at later steps several contributions can cancel. For an exact example, letEliminate indices $1,2$ in that order. The first Schur complement, in order $(2,3,4)$, isThe second pivot element updates the $(3,4)$ entry to $-1/3-(-1/2)(1)/(3/2)=0$. The remaining block is $\mathrm{diag}(7/3,7/3)$, so all four pivot elements are positive and $H$ is a symmetric positive-definite matrix. The elimination graph retains the structural edge $3-4$, but this exact choice of coefficients makes the numerical entry vanish. A small change to $H_{34}=H_{43}$ generally removes the cancellation.

Such an exact zero is different from discarding a small value computed with floating-point arithmetic. The symbolic factorization deliberately reserves space for possible entries without depending on special numerical identities. A perfect elimination ordering creates no structural fill-in at all; an ordering with cancellations need not have that property.

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.

With $A=D+L+U$, the Gauss-Seidel method solves $(D+L)x^{(k+1)}=b-U{x}^{(k)}$. It uses each newly computed coordinate immediately in the remaining triangular solve. It is successive over-relaxation with parameter one, and converges from every start when $A$ is a symmetric positive-definite matrix.

To solve a linear system $Ax=b$, split its matrix into the diagonal matrix $D$ and the strictly lower and upper parts $L,U$, retaining their actual signs. Every diagonal entry must be nonzero. A sweep updates coordinates in order:The Jacobi method uses the old vector throughout; the Gauss-Seidel method reuses each fresh coordinate. Thus the update order is part of the method.^[ref]

The matrix splitting has $M=D+L$, $N=-U$, so its iteration matrix isCompute the update by forward substitution with the lower triangular matrix $D+L$, rather than forming a matrix inverse. For a sparse matrix, a sweep only visits stored entries and can overwrite the working vector. Dependencies on fresh coordinates constrain parallel execution. This makes the Gauss-Seidel method an inexpensive solver component, including in a multigrid method.^[ref]

George S. Lueker coauthored research on graph representations of vertex elimination. The work connects a combinatorial ordering problem with the computational fill produced by sparse matrix factorization.^[ref]

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]

The graphical Lasso estimates a sparse precision matrix by minimizing a Gaussian negative log-likelihood plus an entrywise ${\ell }^1$ penalty.

Writing $A=D+L+U$, where $D$ is diagonal and $L,U$ are strictly triangular, the Jacobi iteration is

The Jacobi method solves a linear system by updating every coordinate from the same old vector. It requires every diagonal entry $a_{ii}$ to be nonzero. Isolating coordinate $i$ in row $i$ gives the simultaneous ruleHere $B$ is the iteration matrix. A sweep computes a complete new vector before replacing the old one, so its coordinate updates can run independently.^[ref]

For a sparse matrix, a sweep needs work proportional to the stored nonzero entries, while Gaussian elimination may create fill-in. The Jacobi method is therefore a simple starting point for large linear systems, and the weighted Jacobi method supplies a useful building block of a multigrid method. Its beauty is that a local row correction becomes a global eigenvalue question: the convergence criterion for a stationary iteration is exactly $\rho (B)<1$.

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.

Numerical linear algebra develops stable finite algorithms for matrix computations.

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,

Construct a fresh sparse matrix, transform each stored value, and omit products that become zero.Filtering preserves column order. It handles s == 0.0 and any floating point underflow producing zero. Sharing immutable list structure elsewhere would be safe; this code constructs new transformed rows. Within the class, Scala's ordinary private permits access to out.rows. The result isThe time complexity is $O(N+M)$ for $M$ stored entries. As usual for the supplied Double model, arithmetic means machine arithmetic; for a finite-real interpretation assume finite operands and products, excluding infinities and NaNs.

Entry $(j,x)$ in source row i becomes entry $(i,x)$ in destination row j. Visit source rows in descending order and prepend, producing increasing destination column order without a sorting pass.Before processing source row i, each destination row contains exactly the contributions from source rows with index greater than i, in increasing index order. Prepending index i preserves that order, and each source position is unique. Hence the representation invariant holds andThe source sparse matrix is unchanged.