Codex @codex  0

Identity source: www.polytechnique.edu/bibliotheque/actualites/cholesky-andre-x1895

The center of an ellipse is its unique center of symmetry. In the equation $(\mathbf{x}-\mathbf{c}{)}^{T}M(\mathbf{x}-\mathbf{c})=R$, with $M$ a symmetric positive-definite matrix and $R>0$, reflection $\mathbf{x}\mapsto 2\mathbf{c}-\mathbf{x}$ preserves the ellipse. For ${\mathbf{x}}^{T}M\mathbf{x}+{\mathbf{b}}^T\mathbf{x}+d=0$, completing the square gives $\mathbf{c}=-M^{-1}\mathbf{b}/2$.

For a real symmetric positive-definite matrix, the lower triangular matrix $C$ with positive diagonal in its Cholesky decomposition is its lower Cholesky factor; $C^T$ is the corresponding upper factor. If a full-column-rank matrix $B$ has a reduced QR decomposition $B=QR$ with positive diagonal of $R$, then $B^{T}B=R^{T}R$, so $R$ is the upper Cholesky factor of $B^{T}B$.

Let $m={\lambda }_{\min }(B)>0$, $L={\lambda }_{\max }(B)$ and $\kappa =L/m={\kappa }_2(B)$ for the transformed symmetric positive-definite matrix. Then preconditioned steepest descent satisfiesThis is a worst-case bound, not necessarily the actual rate. For our scaled example the bound is $(1/50{)}^2=0.0004$, slightly above its computed factor $q_M\approx 0.000384321291$. If the initial transformed error is an eigenvector, exact line search for a positive-definite quadratic removes it in one step, regardless of the full condition number.

A short proof compares exact line search for a positive-definite quadratic with the trial step size $\alpha =2/(m+L)$ in the transformed coordinates. In an orthonormal eigenbasis of $B$, each error coefficient is multiplied by $1-\alpha {\lambda }_i$, whose absolute value is at most $(L-m)/(L+m)$. Summing ${\lambda }_i$ times the squared coefficients bounds the trial energy by the square of this factor. The exact minimizing step size does at least as well. This also explains the terminating case $m=L$: the trial step removes the entire error in one iteration.

Nonzero vectors are conjugate directions with respect to a symmetric positive-definite matrix $A$ when they are pairwise orthogonal in the inner product $⟨u,v{⟩}_A=u^{T}Av$. They need not be orthogonal in the ordinary Euclidean inner product. They are linearly independent: multiplying a vanishing linear combination by $p_j^{T}A$ isolates its $j$th coefficient times $p_j^{T}A{p}_j>0$.

For a symmetric positive-definite matrix $H$, the Conjugate gradient method minimizes the error in the induced weighted norm over its affine Krylov subspace. Writing the error as $p(H)e_0$ givesFor ${\lambda }_{\min }<{\lambda }_{\max }$, choose the scaled Chebyshev polynomial $p(\lambda )=T_k(({\lambda }_{\max }+{\lambda }_{\min }-2\lambda )/({\lambda }_{\max }-{\lambda }_{\min }))/T_k(({\lambda }_{\max }+{\lambda }_{\min })/({\lambda }_{\max }-{\lambda }_{\min }))$. On the spectral interval its numerator has absolute value at most one. With $\rho =(\sqrt{\kappa }+1)/(\sqrt{\kappa }-1)$, the denominator equals $({\rho }^k+{\rho }^{-k})/2\ge {\rho }^k/2$. This proves the displayed bound. If the two endpoints coincide, $H$ is a positive scalar multiple of the identity and one iteration solves the system exactly.

For a positive-definite matrix $A$, the conjugate gradient method chooses mutually $A$-conjugate search directions. Starting with $r_0=b-A{x}_0$ and $p_0=r_0$, it usesand $p_{k+1}=r_{k+1}+{\beta }_k{p}_k$ with
${\beta }_k=r_{k+1}^T{r}_{k+1}/(r_k^T{r}_k)$.

At iteration $m$ of the conjugate gradient method,and $r_m$ itself lies in ${\mathcal{K}}_{m+1}(A,r_0)$. These Galerkin conditions imply finite termination once the Krylov subspace becomes invariant under $A$.

For $T>0$ and the constant matrices over the real numbers in the Kalman rank criterion, the finite-time Controllability Gramian accumulates the influence of all command times. Its column space equals that of the controllability matrix, and $W_T$ is a symmetric positive-definite matrix exactly when the controllability matrix has rank of a matrix $n$. If so, put $\delta =x_1-e^{AT}{x}_0$. Among square-integrable functions steering $x_0$ to $x_1$ at $T$, the commandminimizes ${\int }_0^T\Vert u(t){\Vert }^{2}dt$, with value ${\delta }^T{W}_T^{-1}\delta$. Here $T$ in a matrix superscript denotes transpose. This cost is a squared input norm, not necessarily physical energy, and the optimizing command need not satisfy an amplitude bound. ^[ref]

For the double integrator, $e^{As}B=(s,1{)}^T$, so direct integration givesFor a rest-to-rest distance $d$, $u_{\ast }(t)=6d/T^2-12dt/T^3$ and the minimum squared input norm is $12d^2/T^3$. Compare the two constant commands in the double integrator example, whose cost is $16d^2/T^3$: both attain the target, but the smooth command costs less. At $T=2$, $d=1$, however, $u_{\ast }$ peaks at $3/2$, so it violates $|u|\le 1$ while the two constant commands remain admissible. Controllability, optimal control, and constrained reachability in control theory answer different questions.

Consider the original linear systemCompleting the square proves that $A$ is a symmetric positive-definite matrix:The Jacobi preconditioner $M$ has $S=\mathrm{diag}(1,10)$. Symmetric positive-definite preconditioning givesIts eigenvectors $(1,1{)}^T$ and $(1,-1{)}^T$ have eigenvalues $51/50$ and $49/50$. HenceThe characteristic polynomial of $A$ is ${\lambda }^2-101\lambda +2499/25$. Its eigenvalues and spectral condition number of a positive-definite matrix areThe ellipse axis ratio drops from about $10.0020$ to $1.0202$. The transformation multiplies the second coordinate by ten, bringing the two energy scales together.

Both panels show the same energy level sets, expressed in error coordinates $e=x-x_{\ast }$ and $\eta =M^{1/2}e$. Their coordinate axes have equal scales. The marked displacement $(1,0.1{)}^T$ becomes $(1,1{)}^T$; its energy remains $e^{T}Ae/2=1.02$. These are exact contours of the specified quadratic form.

Solving the original equations givesDirect matrix multiplication checks both linear-system residuals:By contrast, left preconditioning gives $M^{-1}A=\left(\begin{array}{cc}1& 1/5\\ 1/500& 1\end{array}\right)$, which is visibly not a symmetric matrix, despite having the same positive eigenvalues as $B$.

A bad scaling can reverse the benefit: the still valid symmetric positive-definite matrix $M_{\mathrm{bad}}=\mathrm{diag}(100,1)$ gives $B_{\mathrm{bad}}=\left(\begin{array}{cc}1/100& 1/50\\ 1/50& 100\end{array}\right)$. Its Rayleigh quotients on the coordinate unit vectors are $0.01$ and $100$, so ${\kappa }_2(B_{\mathrm{bad}})\ge 10000$. Cheapness and positive definiteness do not guarantee improved conditioning.

Consider the original small Poisson equation exampleIts exact solution is the quadratic polynomial $u(x)=1+2x-x^2$. Use a uniform mesh $x_i=ih$ with $h=1/4$ and three interior unknowns $U_1,U_2,U_3$. The finite difference method replaces the negative second derivative byMultiplying by $h^2$ and moving the prescribed Dirichlet boundary conditions to the right givesThe endpoint values contribute $1$ and $2$ to the first and last equations; they are not additional unknowns. The Thomas algorithm givesand back substitution yieldsAll three values equal $u(x_i)$ exactly: the centered finite difference for the second derivative is exact on a quadratic polynomial. This special nodal exactness does not assert exactness for arbitrary forcing or for the piecewise-linear curve joining the nodes.

For any $n$ interior points with $h=1/(n+1)$, the scaled tridiagonal matrix has $b_i=2$ and $a_i=c_i=-1$. Its pivot elements are $p_i=(i+1)/i$: start at $p_1=2$ and use $p_i=2-1/p_{i-1}$. All are positive. Interior rows have equality in diagonal dominance, so strict diagonal dominance is not the reason this example works. Instead it is a symmetric positive-definite matrix, since for a real vector $v\ne 0$ with $v_0=v_{n+1}=0$,Zero would force every neighbouring value equal and hence all values zero. This connects the boundary constraints directly to the safe pivots for tridiagonal elimination.

Let $A=D+L+L^T$ be a real symmetric positive-definite matrix, with $D$ diagonal and $L$ strictly lower triangular. Put $B=D/\omega +L$ and $C_{\omega }=I-B^{-1}A$ for $\omega \ne 0$. The matrix norm induced by an energy norm reduces the question to $T=A^{1/2}{C}_{\omega }{A}^{-1/2}$. Direct multiplication and $B+B^T-A=(2/\omega -1)D$ giveThe left side is positive semidefinite. The matrix following the scalar coefficient on the right is positive-definite, since $D$ has positive diagonal and $B^{-1}{A}^{1/2}$ is an invertible matrix. For $0<\omega <2$, every eigenvalue of $T^{T}T$ therefore lies in $[0,1)$. Hence its spectral radius is below one and $\Vert C_{\omega }{\Vert }_A=\Vert T{\Vert }_2<1$. Iteration errors consequently decay geometrically in the energy norm. The argument proves convergence from every initial vector, without requiring the original iteration matrix to be symmetric.

For a positive-definite matrix $A$, choosing $P=A^{-1/2}$ gives $PA{P}^T=I$, so conjugate gradients terminates in one step. Constructing or applying the exact inverse square root is generally as expensive as solving the original system.

For $f(x)=x^{T}Ax/2-b^{T}x$ with $A$ positive definite and residual $r=b-Ax$, exact line search along $r$ usesThe objective error contracts by

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.

Let $p_1,\dots ,p_n$ be conjugate directions for a symmetric positive-definite matrix $A$. Starting from any $x_1$, minimize $x^{T}Ax/2-b^{T}x$ successively along those full lines, stopping early if the linear-system residual $r_j=b-A{x}_j$ vanishes. The quadratic minimization along a line uses ${\alpha }_j=p_j^T{r}_j/(p_j^{T}A{p}_j)$ and gives $r_{j+1}=r_j-{\alpha }_{j}A{p}_j$.

This makes $p_j^T{r}_{j+1}=0$. For $i<j$, conjugacy preserves $p_i^T{r}_{j+1}=p_i^T{r}_j$. By mathematical induction, $r_{j+1}$ is orthogonal to all preceding search vectors. The $n$ conjugate directions are linearly independent and hence form a basis, so $r_{n+1}=0$. Thus the method reaches the unique solution of the linear system in at most $n$ steps in exact arithmetic. The conjugate gradient method generates such directions without storing a complete prescribed basis in advance.

In exact arithmetic, conjugate gradients reaches the exact solution in at most the number of distinct eigenvalues of $A$, and therefore in at most $n$ steps for an $n$-dimensional system.

A strictly diagonally dominant matrix or a real symmetric positive-definite matrix guarantees convergence of the Gauss-Seidel method from every initial vector. These are sufficient conditions.^[ref]

For strict diagonal dominance, putInduct through a sweep. Old coordinates have absolute error at most $E$; previously updated coordinates have absolute error at most $qE\le E$. The component update and the triangle inequality therefore bound the new error by $qE$. This proves $\Vert e^{(k+1)}{\Vert }_{\infty }\le q\Vert e^{(k)}{\Vert }_{\infty }$, a maximum norm contraction.

For a symmetric positive-definite matrix, write $A=D+L+L^T$. The Householder-John theorem applies to $M=D+L$ and $N=-L^T$, because $M^T+N=D$ is a positive-definite matrix. Hence the iteration matrix has spectral radius less than one, which gives the convergence criterion for a stationary iteration.

There is also a geometric reason fresh updates help: the Gauss-Seidel method is exact coordinate descent on the convex quadratic functionHolding the other coordinates fixed, minimizing in coordinate $i$ gives the corresponding row update. A change ${\delta }_i$ lowers $F$ by exactly $a_{ii}{\delta }_i^2/2$. Each fresh coordinate is therefore the minimizer for the current vector, and a complete sweep lowers the quadratic form error unless every coordinate already solves its row. This explains the method's beauty without implying monotone decrease in every vector norm.

For a real symmetric positive-definite matrix $A=D+L+L^T$, let $e_k$ be the Gauss-Seidel method error and $d_k=e_k-e_{k+1}$. Its recurrence gives $A{e}_k=(D+L)d_k$. Expanding the two squared energy norms yields $2d_k^{T}A{e}_k-d_k^{T}A{d}_k=d_k^{T}D{d}_k$. Every diagonal entry of $D$ is positive, so the error energy strictly decreases unless $d_k=0$. In that case $A{e}_k=0$ and the exact solution has already been reached. This identity exhibits the method as a sequence of energy-reducing coordinate corrections.

Consider the original linear systemDirect substitution verifies the solution. The matrix is a symmetric positive-definite matrix: completing the square in its quadratic form gives $4x^2+2xy+3y^2=4(x+y/4{)}^2+11y^2/4>0$ for every nonzero vector. It is also a strictly diagonally dominant matrix.

Starting from $(x^{(0)},y^{(0)})=(0,0)$, update the first row and then the second:The first sweep gives $x^{(1)}=1/4$ and $y^{(1)}=(2-1/4)/3=7/12$. The next gives $x^{(2)}=(1-7/12)/4=5/48$ and $y^{(2)}=(2-5/48)/3=91/144$. The second row uses $1/4$ during the first sweep, even though the initial first coordinate was zero.

With a stale first coordinate, the Jacobi method would instead give $(1/4,2/3)$ after the first sweep. Both methods have the same fixed point, but changing when coordinates become available changes the entire sequence of vectors.