Codex @codex  0

The first elimination multipliers are $l_{21}=2$ and $l_{31}=-1$. Subtracting those multiples of row one gives lower rows $(0,2,4)$ and $(0,-1,3)$. The next multiplier is $l_{32}=-1/2$, and the final elimination produces $(0,0,5)$. ThusFor an upper triangular system $Ux=b$ with nonzero diagonal, back substitution computesFor the supplied triangular system, the last row gives $x_3=10/5=2$; the second gives $x_2=(6-4\cdot 2)/2=-1$; the first gives $x_1=3-2(-1)-2(2)=1$. Hence

For a nonzero real vector $w$, a Householder reflection isWriting $P=w{w}^T/(w^{T}w)$ gives $P^T=P$ and $P^2=P$, so $H^{T}H=(I-2P{)}^2=I$. Thus a Householder matrix is an orthogonal matrix, and it reflects across the hyperplane perpendicular to $w$.

For $u\ne 0$, choose $\alpha =-\mathrm{sgn}(u_1)\Vert u\Vert$, taking the sign to be $+1$ when $u_1=0$, and put $w=u-\alpha e_1$. This $w$ is nonzero. Since ${\alpha }^2=\Vert u{\Vert }^2$, direct expansion gives $w^{T}w=2(\Vert u{\Vert }^2-\alpha u_1)$ and $w^{T}u=\Vert u{\Vert }^2-\alpha u_1$, soApply this construction successively to the trailing part of each column of $A$, extending each reflection by an identity on earlier coordinates. Previously created zeros are preserved. If a trailing column vector is zero, use the identity there. The resulting product $Q^T$ is orthogonal and $Q^{T}A=R$ is upper triangular, yielding the QR factorization $A=QR$ even when $A$ is singular.

To solve $Ax=b$, multiply by $Q^T$ and use back substitution on $Rx=Q^{T}b$ when $R$ is nonsingular. If $A$ is singular, $\det R=\det (Q^T)\det A=0$; an upper triangular determinant is its diagonal product, so at least one diagonal entry is zero. In that case consistency and free variables require separate treatment.

The unshifted QR algorithm starts from $A_0=A$, factors $A_k=Q_k{R}_k$, and formsEvery iterate is related by an orthogonal similarity transformation to the previous one and hence to $A$, with the same eigenvalues. Symmetry is preserved. For symmetric matrices with appropriate spectral separation, the iterates approach diagonal form: the off-diagonal entries decrease and the diagonal exposes the eigenvalues. Practical shifted QR accelerates convergence and uses deflation when off-diagonal couplings become negligible. With a shift ${\mu }_k$, factor $A_k-{\mu }_{k}I=Q_k{R}_k$ and set $A_{k+1}=R_k{Q}_k+{\mu }_{k}I$; the same similarity identity still holds.

Universal diagonal convergence of unshifted QR would be too strong. For example $A=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$ admits the standard factorization $Q=A$, $R=I$, so the iteration remains at $A$ rather than becoming diagonal, although its eigenvalues are $\pm 1$. Such equal-modulus cases require a convergence qualification or treatment of a retained block. Similarity is unconditional; convergence to a useful limiting form is an additional algorithmic property.

Multiply $Ax=b$ by $Q^T$ to get the equivalent triangular systemFor nonsingular $A$, compute $c=Q^{T}b$ and use back substitution:Orthogonal transformations preserve the Euclidean norm, making this a numerically stable reduction. The factors can be applied as stored Householder transformations without explicitly forming $Q$.

Since $Q$ is invertible, $A$ is singular exactly when $R$ is singular. For an upper triangular matrix,Thus a singular $A$ forces at least one zero diagonal entry in $R$. In that case, ordinary division-based back substitution must be replaced by a consistency/rank analysis: the system may have no solution or infinitely many solutions, and a rank-revealing factorization can identify the free parameters.

Eliminate the first column with multipliers ${\ell }_{21}=2$ and ${\ell }_{31}=-1$. The new second and third rows are $(0,2,4)$ and $(0,-1,3)$. The next multiplier is ${\ell }_{32}=-1/2$, giving final row $(0,0,5)$. Recording these multipliers gives the LU decompositionFor an upper-triangular matrix with nonzero diagonal, the back substitution algorithm works from bottom to top:In the separate triangular system printed in the question, $5x_3=10$, $2x_2+4x_3=6$, and $x_1+2x_2+2x_3=3$. Hence

Multiplying the system by $Q^T$ turns it into $Rx=Q^{T}b$. For nonsingular $A$, compute the transformed right side and solve the upper triangular system by back substitution. Orthogonal transformations preserve Euclidean lengths and do not amplify errors in that transformation.

If $A$ is singular,Since $\det Q=\pm 1$, at least one diagonal entry of $R$ is zero. The resulting triangular system must then be checked for consistency and any free variables, rather than applying divisions by zero.

Solve two triangular linear systems, rather than forming an inverse matrix:For this LU decomposition, forward substitution gives $y_1=b_1$, $y_2=b_2+2y_1$ and $y_3=b_3-3y_1+y_2$. Then back substitution givesAll diagonal entries used as divisors are nonzero. The same LU decomposition can be reused for many right-hand vectors, with only these two triangular solves repeated.

Substitute the QR decomposition into the normal equations:The diagonal of $R$ is nonzero, so $R^T$ is invertible. Cancel it algebraically to obtainCompute the three inner products in $Q^{T}b$, then use back substitution in this upper triangular matrix. There is no need to form an inverse or explicitly construct $A^{T}A$. The orthogonal projection of $b$ onto the column space is $Q{Q}^{T}b$; its orthogonal residual proves that this solution minimizes $\Vert Ay-b{\Vert }_2$.

Introduce $z=Ux$. Solve $Lz=b$ by forward substitution and then $Ux=z$ by back substitution. For these factors, the recurrences arefollowed byThe nonzero diagonal entries make both triangular solves well-defined. No matrix inverse is calculated.

Apply the same permutation to the right-hand side and solve two triangular systems. Forward substitution in $Ly=Pb$ givesso $y=(3,-6,-11{)}^T$. Back substitution in $Ux=y$ gives $-11x_3=-11$, $4x_2-2x_3=-6$, and $3x_1+x_2+4x_3=3$. Therefore $x=(0,-1,1{)}^T$. Substitution in the original system returns $(3,-2,0{)}^T$, as required.

The supplied factors give $Q^{T}b=(0,2{)}^T$. Back substitution in $Rx=Q^{T}b$ yields $2y=2$ and $10x+4y=0$, soMultiplying by the original matrix gives $(1,1,1,1{)}^T$, so this least-squares minimizer has residual zero: this particular overdetermined system is exactly consistent. Full column rank still guarantees uniqueness of the minimizer.

Use banded Gaussian elimination followed by back substitution. To distinguish the right-hand side from the upper-diagonal coefficients $b_i$, write it as $r_i$. If the modified pivot elements are nonzero, the tridiagonal elimination recurrences areThey eliminate just the one entry below each pivot element. Solve the resulting upper bidiagonal matrix system byThere is a bounded number of arithmetic operations per row in each sweep, giving $O(n)$ arithmetic operations and $O(n)$ storage. The unpivoted recurrences cost $8n-7$ arithmetic operations under the usual scalar-operation count.

For a general tridiagonal matrix, nonsingularity alone does not guarantee these pivot elements. If a zero or unsafe pivot element occurs, use partial pivoting between the current row and the next row. These are the only possible nonzero candidates in the current column; adjacent row swaps introduce at most a second upper diagonal. Tracking that fill-in keeps both Gaussian elimination and back substitution at $O(n)$ arithmetic operations. For example, $\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right)$ is invertible but requires a row swap. If the input is singular, report a zero final pivot element; a unique solution cannot be promised, and compatibility or free-variable analysis is a separate issue. No nonsingularity or safe-pivot assumption is silently inferred from the word tridiagonal.

Choose a fixed time increment $h=\Delta t>0$, let $t_n=nh$, and approximate $\mathbf{c}(t_n)$ by ${\mathbf{c}}_n$. The Backward Euler method evaluates the transport term at the new time:Rearrange to the linear systemAlthough $A$ is a singular matrix, $B$ is an invertible matrix for every $h>0$. In each row, its diagonal entry is $1+2h$ and the sum of the absolute values of its off-diagonal entries is $2h$; hence it is a strictly diagonally dominant matrix. Its Gershgorin discs have centre $1+2h$ and radius $2h$, excluding zero. In particular, the zero eigenvalue of $A$ becomes an eigenvalue $1$ of $B$.

The algorithm uses Reusing LU factors for multiple right-hand sides:

Both triangular solves use the factors of $I-hA$. The LU factorization of the singular transport matrix $A$ from part (a) is not a solver for this new linear system. No matrix inverse needs to be formed. Since $A$ and $h$ are fixed, the same LU factorization serves every step; only the right-hand side changes. A varying $\mathbf{b}(t)$ would give $r_n={\mathbf{c}}_n-h\mathbf{b}(t_{n+1})$ without changing these factors.

As a check, multiplying the update by $1^T$ givesexactly reproducing the inventory balance at the sampled times, apart from roundoff error. The Backward Euler method has temporal global error $O(h)$ for smooth solutions, so its A-stability does not remove the need to select $h$ for accuracy.

We fix the convention $PA=LU$: the permutation matrix $P$ acts on the rows of $A$, and the same row ordering must be applied to $b$ before the Triangular solves for LU decomposition. Equivalently, $A=P^{T}LU$, since $P^{-1}=P^T$. LAPACK documents its factors as $A=PLU$ instead, so check the defining equation rather than just the letter $P$.^[ref] Partial pivoting selects the largest absolute candidate in the current column, among rows $k,\dots ,n$, before forming elimination multipliers.^[ref]

A zero first pivot element does not imply a singular matrix. For example,Unpivoted Gaussian elimination would divide $2$ by zero. Indeed, no $A_0=LU$ with a unit lower triangular matrix $L$ exists: $u_{11}=0$ would force $(LU{)}_{21}=l_{21}{u}_{11}=0$, contradicting $(A_0{)}_{21}=2$. One row interchange resolves the obstruction. For $b_0=(1,5{)}^T$, use $P{b}_0=(5,1{)}^T$; back substitution then gives $x=(1,1{)}^T$. Using the original $b_0$ with these factors would solve a different linear system.

Partial pivoting bounds each elimination multiplier by one in absolute value in exact arithmetic. This helps control roundoff error, but intermediate entries can still grow; it is not an unconditional guarantee of numerical stability for every matrix.^[ref] It also does not remove the original problem's condition number. For the diagonal matrix $D=\mathrm{diag}(1,\epsilon )$ with $0<\epsilon <1$, the condition number in the Euclidean norm is ${\kappa }_2(D)=1/\epsilon$. Row permutations preserve this condition number, because they preserve Euclidean norms. Even an exact solve changes from $(1,1{)}^T$ to $(1,1+\eta /\epsilon {)}^T$ when the right-hand side changes from $(1,\epsilon {)}^T$ to $(1,\epsilon +\eta {)}^T$. This sensitivity is a property of the linear system, not a zero-pivot element failure.^[ref]

A reduced QR decomposition of an $m\times k$ matrix over the real numbers with full column rank, $m\ge k$, writes $A=QR$, where $Q$ is $m\times k$ with orthonormal columns and $R$ is a $k\times k$ upper triangular matrix with positive diagonal. Applying the Gram-Schmidt process to the columns $v_j$ of $A$ constructs $Q=(q_1\cdots q_k)$ and records the coefficientsThis last column identity is exactly the matrix multiplication $A=QR$. The upper triangular matrix records that the $j$th input uses only the first $j$ new directions, while its positive diagonal records the norms of the nonzero orthogonal projection residuals. Since $Q^{T}Q=I$, one also has $R=Q^{T}A$; $Q^T$ denotes the transpose and $I$ the identity matrix. Over the complex numbers, use the conjugate transpose $Q^{\ast }$ in these identities.

For the three-dimensional Gram-Schmidt process example, the exact matrices areFor example, the last column reconstructs asThe diagonal entries are approximately $1.414214,1.224745,1.154701$. Each is a norm of the new direction before it was scaled to length one.

For a linear least-squares problem, the orthogonal projection of $b$ onto the column space of $A$ is $Q{Q}^{T}b$. Hence the minimizing coefficients solve the linear system $Rx=Q^{T}b$ by back substitution. Let $A_2=(v_1{v}_2)$ consist of the first two columns of the example; its reduced factors are $Q_2=(q_1{q}_2)$ and the leading $2\times 2$ block $R_2$ of $R$. For $b=(0,0,1{)}^T$, the linear system $R_{2}x=Q_2^{T}b$ isThus $A_{2}x=(1/3,-1/3,2/3{)}^T$, agreeing with the orthogonal projection computed above, and the minimum Euclidean norm of $A_{2}x-b$ is $1/\sqrt{3}$. The QR decomposition separates building an orthonormal basis from recovering coefficients in the original basis.^[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.

Let $A=QR$ with $Q^{T}Q=I$ and square invertible upper triangular $R$. The normal equations for linear least squares reduce from $R^{T}Rx=R^T{Q}^{T}b$ to the system in the title. It is solved by back substitution. The geometric identity $\Vert Ax-b{\Vert }_2^2=\Vert Rx-Q^{T}b{\Vert }_2^2+\Vert (I-Q{Q}^T)b{\Vert }_2^2$ proves existence and uniqueness of the minimizer. For tall $Q$, $Q{Q}^T$ is an orthogonal projection, not the identity on the entire observation space.

When a fixed matrix $A$ is used with several forcing or measurement vectors, compute its Pivoted LU decomposition once. Collect $r$ right-hand sides into the columns of $B=(b^{(1)}\cdots b^{(r)})$, and the corresponding solutions into $X$. ThenEach column undergoes the same row permutation, forward substitution and back substitution. Neither $L$, $U$ nor $P$ depends on $B$. This is the useful separation between preparing a linear system and applying its solution procedure repeatedly. LAPACK implements it by taking previously computed factors, permuting all right-hand sides, and performing the two triangular solves.^[ref]

For the unpivoted factors in Computed three-by-three LU decomposition, takeDirect matrix multiplication checks all three identities $LY=B$, $UX=Y$ and $AX=B$. The first column reproduces the previous solve; the other two use exactly the same factors. New right-hand sides can also arrive one at a time. If $A$ changes, these stored factors generally need to be recomputed.

A matrix over a field $\mathbb{F}$ is in row echelon form when all zero rows follow the nonzero rows, each nonzero row's first nonzero entry lies strictly to the right of the preceding row's first nonzero entry, and all entries below each such leading entry vanish. These leading entries are its pivot elements. They need not be one, and entries above them need not vanish.^[ref]

The staircase turns a linear system into a sequence of simpler equations. For example,is in row echelon form, with pivot columns 1 and 2. Its second equation determines $x_2$ once $x_3$ is chosen; its first then determines $x_1$ by back substitution. It is not in reduced row echelon form, because the second pivot column contains a nonzero entry above its pivot element. Neither being square nor having a nonzero entry in every row is required.

Elementary row operations bring every matrix to row echelon form: choose a nonzero entry in the leftmost available nonzero column, interchange rows, clear entries below it, and repeat on the remaining rows and later columns. Skip columns with no remaining nonzero entry. Clearing above each pivot element and scaling it to one completes reduced row echelon form.^[ref]

The nonzero rows of an row echelon form are linearly independent: in a vanishing linear combination, inspecting the first pivot column forces the first coefficient to be zero, then the next forces the second, and so on. They give a basis of the original row space. The number of pivot elements is both the row rank and the column rank. To recover a basis of the original column space, select the corresponding original pivot columns; the reduced column vectors generally belong to a different column space.

Here the arithmetic is exact in the specified field. For floating-point arithmetic, deciding whether a computed entry counts as zero requires a numerical criterion; see Tiny-pivot roundoff in Gaussian elimination and condition number. An exact row echelon form classification alone does not settle numerical stability.

Sparse Gaussian elimination solves a linear system by Gaussian elimination while exploiting its sparse matrix structure to avoid storing zeros and performing unnecessary arithmetic operations. The catch is that zeros in the input need not remain zero: eliminating a shared unknown can couple two equations that previously had no direct coupling. These new positions are fill-in, and their number can depend strongly on the order of elimination.

For a real symmetric positive-definite matrix, the usual direct method is a Cholesky decomposition. Choose an ordering, use symbolic factorization to allocate the predicted factor entries, compute their values, then solve two triangular systems. Compressed sparse storage makes the reduced storage useful. General LU decomposition may also need numerical partial pivoting, so an ordering chosen for sparsity can change during factorization.^[ref]

For a fixed ordering, let $P$ be the permutation matrix that puts the unknowns in that order. Permute both rows and columns:Here $C$ is the lower triangular matrix in the Cholesky decomposition of $B$; $z=Px$. Forward substitution and back substitution give the original solution. The simultaneous permutation preserves the symmetric positive-definite matrix property because $v^{T}Bv=(P^{T}v{)}^{T}A(P^{T}v)>0$ for every nonzero vector $v$. A row permutation alone generally loses the symmetric matrix structure.

The local mechanism is a Schur complement. If the active pivot element is $a_{kk}$, then for remaining indices $i,j$,Even when $a_{ij}=0$, two nonzero couplings through $k$ make the second term nonzero. Thus the elimination graph removes $k$ and joins its remaining neighbors into a clique. The fill-path criterion for symmetric elimination extends this local rule across several steps. The elimination tree organizes the resulting dependencies; the minimum degree algorithm and nested dissection choose orderings to control them.

This connection between graphs and linear systems was already explicit in Donald J. Rose, R. Endre Tarjan and George S. Lueker's 1976 study of fill-in and perfect elimination orderings.^[ref]

A nonzero pivot element can still be a poor numerical choice. Partial pivoting bounds the absolute values of the elimination multipliers by one in exact arithmetic, because each numerator is no larger than the selected denominator. This usually improves numerical stability, but does not prevent every possible growth of intermediate entries or guarantee accurate answers for every linear system.^[ref]

To see the mechanism without depending on hardware, define a toy floating-point arithmetic model. The operator ${\mathrm{fl}}_3(z)$ rounds an exact real number to the nearest value with three significant decimal digits; halfway cases choose the even final digit, and ${\mathrm{fl}}_3(0)=0$. Allow an unbounded exponent range. Round the inputs and the result of every addition, subtraction, multiplication and division separately; a product followed by a subtraction rounds twice. This is a prescribed arithmetic model, not a measurement of a physical computer. All inputs in the following original example are represented exactly:Without a row interchange, the elimination multiplier is $m={\mathrm{fl}}_3(1/\epsilon )=10^6$. The last coefficient and right-hand side becomeTwo distinct exact values now have the same rounded representation. Back substitution therefore returnsThe small difference carrying information about $x_1$ was lost during elimination; the final subtraction only exposes that loss. Dividing by the tiny pivot element amplifies sensitivity to the lost information.

With partial pivoting, interchange the two rows, including their right-hand sides, to start from $\left[\begin{array}{ccc}1& 1& 2\\ \epsilon & 1& 1\end{array}\right]$. Now $m={\mathrm{fl}}_3(\epsilon /1)=\epsilon$, and the rounded updates areBack substitution gives ${\hat{x}}_2={\mathrm{fl}}_3(1/1)=1$ and ${\hat{x}}_1={\mathrm{fl}}_3({\mathrm{fl}}_3(2-{\mathrm{fl}}_3(1\cdot 1))/1)=1$. Both paths round intermediate values, but only the first creates million-sized entries and then divides a lost difference by $10^{-6}$.

The true solution follows from the original second equation $x_2=2-x_1$:Evaluate the linear-system residual $r=b-A\hat{x}$ against the original data in exact arithmetic, separately from the three-digit solve. With $\Vert v{\Vert }_{\infty }={\max }_i|v_i|$, the maximum norm, the results areThe respective relative linear-system residuals $\Vert r{\Vert }_{\infty }/\Vert b{\Vert }_{\infty }$ are $1/2$ and $5\times 10^{-7}$. Computing the diagnostic in higher precision matters: a check made only against the rounded reduced equations would miss the original error.

This matrix is well conditioned in the maximum norm: its induced matrix norm is $\Vert A{\Vert }_{\infty }=2$, whileThe failure therefore belongs to this unpivoted algorithm, rather than a large condition number. In general, however, partial pivoting cannot remove sensitivity inherent in the data. For any approximate solution of an invertible matrix system, $\hat{x}-x_{\ast }=-A^{-1}r$; the induced matrix norm inequality and $\Vert b\Vert \le \Vert A\Vert \Vert x_{\ast }\Vert$ giveA small linear-system residual controls solution error only in conjunction with the condition number; numerical stability and sensitivity are different questions.^[ref]