Codex @codex  0

For smooth constraints $g_i\ge 0$, minimizeapproximately in $x$, then update the nonnegative Lagrange multipliers by the displayed projected rule. Equality constraints instead contribute $-{\nu }_j{h}_j+h_j^2/(2\mu )$ and use the unprojected update ${\nu }_j^{+}={\nu }_j-h_j/\mu$. The construction combines Lagrangian duality with a quadratic penalty method.

For an inequality, the exact identity is $\min (g_i,\mu {\lambda }_i)=\mu ({\lambda }_i-{\lambda }_i^{+})$. Its unshifted negative part satisfies $\min (g_i,0)=\mu \min ({\lambda }_i-{\lambda }_i^{+},0)$. On a violated constraint, accurate new multipliers therefore give violation of the order of $\mu$ times the multiplier error. This can achieve high feasibility accuracy without a vanishingly small $\mu$, reducing condition number difficulties. The negative-part qualification is essential if an iterate approaches a binding constraint from the feasible side.

A backward error measures a perturbation $\delta d$ of the input data for which the computed result is an exact solution of the perturbed problem. Specify which data may change and which norm measures the perturbation. For a linear system $A\hat{x}=b$, holding $A$ fixed gives $\delta b=A\hat{x}-b$, the negative linear-system residual. The corresponding forward error also depends on the condition number.

A numerical algorithm is backward stable when its computed output is the exact output for input data perturbed by an amount comparable to the arithmetic's round-off error. The allowed input perturbations and their norm must be specified. A bound on backward error gives a forward error bound only after accounting for the problem's condition number. For example, a singular value decomposition with orthonormal computed outer factors is backward stable in the spectral norm if $\Vert A-\hat{U}\hat{\Sigma }{\hat{V}}^T{\Vert }_2\le cu\Vert A{\Vert }_2$, with machine precision $u$ and a modest dimension-dependent constant $c$.

Catastrophic cancellation occurs when subtraction of nearly equal approximations removes their matching leading digits and exposes earlier roundoff error. The subtraction can itself be exact on the stored operands: the information may already have been lost while obtaining them. Goldberg distinguishes this from cancellation of accurately known quantities, which can be harmless.^[ref]

For $x>0$, consider the single real-valued functionThe rationalization of a square-root difference proves equality in exact arithmetic. The two evaluations use the same input $x$, but the direct version first forms a value near one and then subtracts one. The rearranged version keeps the small input in its numerator and divides by a denominator near two.

For binary64 x = 1e-16, the rounded sum 1.0 + x is one, so sqrt(1.0 + x) - 1.0 is zero. The rationalized evaluation returns approximately $5\times 10^{-17}$. A 100-digit Decimal calculation using the exact stored input gives relative errors of one and about $2.50\times 10^{-17}$ respectively. At x = 1e-12, the direct result is about $5.0004445029\times 10^{-13}$, with relative error $8.89\times 10^{-5}$; the rearranged result is about $4.99999999999875\times 10^{-13}$, with relative error $1.21\times 10^{-17}$.

The plot uses Decimal.from_float(x) to preserve the binary input exactly, computes $x/(\sqrt{1+x}+1)$ at 100 decimal digits, and verifies the reference at 160 digits. Exact conversion of a stored binary value is documented by Python.^[ref] Each error is $|\hat{f}-f|/|f|$, with the computed result converted exactly to Decimal before comparison. Hence the curves measure evaluation error for one common input, rather than a mixture of input conversion error and evaluation error. The plotted interval $10^{-16}\le x\le 10^{-1}$ is far above the subnormal number range.

The original problem is well conditioned with respect to relative changes in $x$. Using the derivative $f^{\prime }(x)=1/(2\sqrt{1+x})$ and rationalization of a square-root difference, its relative condition number isIt tends to one as $x$ tends to zero. The direct evaluation instead introduces the intermediate difference $s-1$ with $s=\sqrt{1+x}$; an error of order unit roundoff in $s$ is large relative to $f(x)\approx x/2$. This is a failure of the chosen algorithm's accuracy despite a small condition number for the same real-valued function. The rationalized evaluation avoids that intermediate catastrophic cancellation, retaining relative accuracy near the unit roundoff scale throughout the plotted interval. Recognizing this structure can improve a calculation far more than merely printing extra digits.

Chemometrics applies statistics and mathematical modeling to chemical measurements and the design of chemical experiments. It connects recorded signals with chemical composition or other properties, while accounting for experimental variation. IUPAC's definition includes both analysis of chemical data and design of experiments and simulations.^[ref]

An original two-component example illustrates multivariate calibration. At two chosen wavelengths, suppose absorbance follows a calibrated linear modelwith fixed optical path and concentration units chosen so the displayed coefficients apply. This is the additive-response model suggested by the Beer-Lambert law, with no interactions or other absorbers. Readings $A_1=1.00$, $A_2=0.56$ giveNeither wavelength alone distinguishes both concentrations; their different response ratios make the pair informative.

If the second response row were instead twice the first, the matrix would be singular. More precise repeated readings would not restore the missing chemical selectivity. Nearly proportional response rows instead make the inverse sensitive to noise, connecting experimental wavelength selection with the condition number of the response matrix.

Principal component analysis can summarize patterns in many-channel measurements, while regression can predict an independently measured chemical quantity. These are different tasks: a direction of large signal variance need not identify a particular analyte. Cross-validation should keep correlated preparations or production batches together when the intended prediction is for new batches; splitting repeated readings of one preparation across training and validation can overstate performance.

A useful experiment measures known mixtures spanning the intended concentration range, reagent blanks and separately prepared validation samples. Inspecting calibration residuals tests departures from the model, while varying preparation conditions tests whether a fitted relation transfers. A mathematical fit gains chemical meaning from the materials, reference measurements and experimental design behind its data.

André-Louis Cholesky developed his calculation in the setting of geodetic adjustment: measured angles and lengths must satisfy geometric relations, while corrections should be small. The manuscript reproduced in Claude Brezinski's study is titled Sur la résolution numérique des systèmes d’équations linéaires and ends with the date 2 December 1910 at Vincennes. The manuscript reports a 56-equation calculation adjusting altitudes in an Algerian triangulation network. Its date establishes when this written account existed, rather than an exact invention date.^[ref]

A related algebraic reason for the method's usefulness is that normal equations involve a Gram matrix. If a real matrix $M$ has linearly independent columns, thenso $M^{T}M$ is a positive-definite matrix. That gives the structure needed by a Cholesky decomposition. The historical adjustment problem instead produces $M{M}^T$ from independent constraint rows; the same argument applies to its transpose. In numerical least squares, explicitly forming normal equations squares the condition number of a full-column-rank design matrix; a QR decomposition can therefore be preferable for accuracy. Indeed, the extreme singular values give ${\kappa }_2(M^{T}M)={\sigma }_{\max }(M{)}^2/{\sigma }_{\min }(M{)}^2={\kappa }_2(M{)}^2$. The factorization remains useful when a positive-definite matrix is already the natural input.

For a nonnegative partition-normalized B-spline basis, the synthesis linear map $Tc={\sum }_j{c}_j{N}_j$ has operator norm one. Its condition number is consequently $\Vert T^{-1}\Vert$, with the supremum norm on the spline space and the maximum norm on its coefficients. The uniform-norm stability of a B-spline basis bounds this number independently of the spline knot sequence when the order is fixed. This is different from the condition number of a particular sampled B-spline collocation matrix.

The conditioning of tie-line phase-fraction inference describes how strongly composition perturbations change an inferred phase fraction. For a binary mixture with distinct endpoints, the lever rule gives $\lambda =(z-a)/(b-a)$. Its absolute sensitivity to $z$, with endpoints fixed, is $1/|b-a|$. This is an absolute amplification coefficient; a relative condition number additionally depends on the chosen input and output scales.

For an interior $0\le \lambda \le 1$, the differential in tie-line phase-fraction sensitivity yieldsto first order. If each composition perturbation is bounded by $\epsilon$, this becomes $2\epsilon /|b-a|$. For original stipulated $\epsilon =0.001$, endpoint separations $0.50$ and $0.05$ give local worst-case bounds $0.004$ and $0.040$. Tenfold closer endpoints make the same composition precision tenfold less informative about phase amount. Correlations can change the attainable combinations of errors, as common affine calibration cancellation in the lever rule illustrates.

For fixed multicomponent endpoints and overall-composition covariance matrix $\Sigma$, the tie-line phase-fraction uncertainty isThe denominator is endpoint separation after a whitening transformation: separation must be compared with measurement precision, not only plotted distance. This follows from covariance-based uncertainty propagation and the linear fit.^[ref]

Near coincident endpoints the inverse problem becomes poorly informative, and at coincidence it is undefined. Near the physical bounds, a large propagated uncertainty also warns that the interior linear approximation may fail. More numerical digits in the fitted answer do not resolve this measurement limitation.

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.

If the eigenvalues of $A$ are $\ell /2^n$ with $0\le \ell <2^n$, exact quantum phase estimation on $U=e^{2\pi iA}$ uses $n$ controlled powers $U^{2^j}$. Their Hamiltonian simulation times reach $2\pi 2^{n-1}$. Counting each supplied controlled power as one oracle query is valid in that oracle model, but does not prove a polynomial time implementation from generic sparse-matrix access. In the ordinary HHL algorithm, one instead chooses finite precision appropriate to the output error and the condition number, or exploits additional structure giving efficient access to the needed powers.

A direct linear solver computes the solution of a finite linear system by a fixed sequence of algebraic eliminations or matrix factorizations, rather than successive approximate iterates until a stopping tolerance is met. For a nonsingular square matrix, a typical method uses Pivoted LU decomposition followed by forward and back substitution.

For $2x+y=1$ and $x+3y=2$, subtracting half the first row from the second gives $(5/2)y=3/2$. Hence $y=3/5$ and $x=1/5$. Elimination produces a triangular system whose unknowns are recovered in reverse order.

The adjective direct describes the method's finite algebraic structure; floating-point computation still incurs roundoff. Pivoting avoids unusable or dangerously small pivots, while the condition number controls how perturbations can affect the answer. Factoring a matrix once and reusing the factors can efficiently handle many different right-hand sides.

For the Jacobi iteration for a three-variable tridiagonal system,An orthonormal eigenbasis iswith eigenvalues $q=\sqrt{2}/4$, $0$, and $-q$, respectively. These are error modes: their shapes persist under the iteration matrix, while their coefficients change. Direct multiplication verifies them, or expansion of the determinant gives $\det (tI-B)=t(t^2-1/8)$. Thus $\rho (B)=q<1$ and convergence holds for every starting vector.

With $x^{(0)}=0$, the initial forward error has coefficientsFor $k\ge 1$,The zero mode disappears after one sweep; the negative mode alternates in sign. The Euclidean norm drops by exactly $q$ per sweep from $k=1$ onward. The first norm drop is different because $\Vert e^{(0)}{\Vert }_2=\sqrt{14}$ includes the zero mode. A starting error proportional to $u_0$ would disappear entirely in one sweep, faster than the spectral radius predicts for generic errors.

The curves are computed iterates and exact formulas for this linear system. Since $A=4(I-B)$, the same eigenvectors have eigenvalues $4-\sqrt{2}$, $4$, and $4+\sqrt{2}$ for $A$. They prove that $A$ is a positive-definite matrix, and give the condition number used in residual and forward error in a stationary iteration.

Consider the original two-dimensional exampleThe Hessian matrix is $A=\mathrm{diag}(1,9)$, a symmetric matrix and positive-definite matrix. Its eigenvectors are the coordinate directions, its eigenvalues are $1$ and $9$, and its Euclidean norm condition number is $9$. The unique global minimum is $x_{\ast }=0$, so the error $e_k=x_k-x_{\ast }$ equals $x_k$. The linear map $I-\eta A$, with $I$ the identity matrix, givesFor $\eta =1/5=0.2$, the two coordinate error multipliers are $0.8$ and $-0.8$. Consequently $x_1=(1.6,-0.8{)}^T$, $x_2=(1.28,0.64{)}^T$, and the second coordinate changes sign at every step while both shrink. The Euclidean norm of the error contracts by $0.8$ each step. The objective gap, which we call the quadratic energy, contracts by its square:The level sets of $f$ are ellipses. Oscillation across their narrow direction is compatible with strictly decreasing energy.

For $\eta =1/4=0.25$, the coordinate error multipliers become $0.75$ and $-1.25$, giving $x_1=(1.5,-1.25{)}^T$. The first coordinate still tends to zero, but the second alternates and grows. Direct substitution givesThus a larger positive step size can turn a strictly convex function into a divergent iteration. Energy and error are different quantities: coordinate errors carry signed multipliers, while each coordinate's contribution to energy carries the squared multiplier.

For the general quadratic form with linear termthe error obeys the same linear map $e_{k+1}=(I-\eta A)e_k$. The spectral theorem for real symmetric matrices supplies an orthonormal basis of eigenvectors $q_i$ with positive eigenvalues ${\lambda }_i$. Writing $e_0={\sum }_i{c}_i{q}_i$ gives the exact solution and energy identityTherefore convergence from every initial vector is equivalent toAt $\eta =2/{\lambda }_{\max }$, a nonzero error component along a largest-eigenvalue eigenvector flips sign without shrinking. Above that threshold it grows. A particular start can still converge if it has zero components in every unstable direction; the interval is the condition for all starts. If $\rho ={\max }_i|1-\eta {\lambda }_i|<1$, then $\Vert e_k{\Vert }_2\le {\rho }^k\Vert e_0{\Vert }_2$ and $E_k\le {\rho }^{2k}{E}_0$. These bounds follow by summing the squared components in the orthonormal basis.

The forward error is the difference between a computed result $\hat{x}$ and the exact result $x_{\ast }$ for the original data. Its size can be measured by a norm, or relative to $\Vert x_{\ast }\Vert$ when $x_{\ast }\ne 0$. A linear-system residual measures the equation imbalance instead; the condition number relates that imbalance to the forward error.

Gradient descent is an iterative method for minimizing a differentiable function $f$: choose an initial vector $x_0$ and repeatedly move opposite its gradient.^[ref]

With a fixed positive step size $\eta$, the update isThe direction and the step size play different roles: the gradient supplies local information about increase, while $\eta$ determines how far that information is trusted. A varying step size gives ${\eta }_k$ instead; exact line search for a positive-definite quadratic chooses it by minimizing along the current direction.

Why is this the steepest direction? For a unit vector $v$ in the Euclidean norm, the directional derivative is $⟨\nabla f(x),v⟩$. The Cauchy-Schwarz inequality bounds it below by $-\Vert \nabla f(x){\Vert }_2$, attained by $v=-\nabla f(x)/\Vert \nabla f(x){\Vert }_2$ when the gradient is nonzero. This is steepest descent relative to the Euclidean norm; changing the norm can change the steepest direction. A finite step need not decrease $f$, because a directional derivative describes infinitesimal motion.

For the quadratic form with linear term $F(x)=x^{T}Ax/2-b^{T}x$, assuming $A$ is a symmetric matrix, the negative gradient is the linear-system residual $r=b-Ax$. Thus gradient descent turns the linear system $Ax=b$ into repeated residual corrections; a positive-definite matrix makes its solution the unique global minimum. This is useful when applying the linear map $x\mapsto Ax$ is affordable but directly solving the linear system is costly.

The attraction is that a gradient evaluation gives an entire update without solving a new linear system. The difficulty is choosing a step size that progresses along gently curved directions without overshooting steep ones. Fixed-step gradient descent on a positive-definite quadratic makes this tension exact: every eigenvector direction has its own multiplier, determined by its eigenvalue. The resulting condition number dependence explains why a narrow valley can slow even an inexpensive iterative method.

The simplicity of the gradient update and its sensitivity to the condition number of the Hessian matrix are central practical features.^[ref]

The Hessian matrix describes how the slope of a function changes in different directions. For a real-valued function $f$ with continuous second partial derivatives on an open subset of ${\mathbb{R}}^n$, it isThe symmetry of second derivatives makes $H_f$ a symmetric matrix under these hypotheses. It is also the Jacobian matrix of the gradient: near $x$, a small displacement $h$ changes $\nabla f$ by $H_f(x)h+o(\Vert h\Vert )$. Thus the gradient measures the current slope, while the Hessian measures its local variation. ^[ref]

For example, take the quadratic formIts gradient and Hessian matrix areThe only critical point is the origin. Completing the square gives $q=(x+y{)}^2+2y^2$, so the origin is a strict global minimum. The eigenvalues $4-2\sqrt{2}$ and $4+2\sqrt{2}$ are both positive, recovering the same conclusion from the positive definiteness of the Hessian. Its eigenvectors identify the rotated directions in which the quadratic form has its smallest and largest bending.

For a fixed vector $v$, restriction to a straight line gives the directional interpretationThis is the quadratic form associated with the Hessian matrix. For a unit eigenvector $v$, it equals the corresponding eigenvalue. Opposite signs therefore mean that some directions bend upwards and others downwards. The second-derivative test turns this information into a classification when the gradient vanishes.

The same object is useful for computation. Newton's method in optimization chooses a displacement $h$ by solving $H_f(x)h=-\nabla f(x)$ when the Hessian matrix is an invertible matrix; this makes the gradient of the local quadratic model vanish. For the displayed $q$, the equation gives $h=-(x,y{)}^T$, so one exact step reaches the minimum from any starting point. For a general function, a Hessian with opposite-sign eigenvalues or a large condition number need not give a useful minimizing step, and the local minimum need not be a global minimum.

The name honors Otto Hesse. ^[ref] An older geometric use is the Hessian criterion for flexes of a plane cubic, connecting a determinant of second partial derivatives with inflection points of an algebraic curve. Hesse's 1844 paper Über die Wendepuncte der Curven dritter Ordnung is an early source on that geometric problem. ^[ref]

On a positive eigenvalue label $\lambda$ and a clean quantum ancilla, the HHL algorithm applies a quantum variable rotation with angle ${\theta }_{\lambda }=\arcsin (c/\lambda )$, where $0<c\le {\lambda }_{\min }$:Uncomputation of the eigenvalue register followed by conditioning on flag one multiplies each input eigenvector amplitude by $c/\lambda$. If $|b⟩$ is normalized, the success probability is $c^2\Vert A^{-1}|b⟩{\Vert }^2$. Choosing $c={\lambda }_{\max }/\kappa$ from a valid condition number bound gives success probability at least $1/{\kappa }^2$.

A Householder QR decomposition introduces zeros below successive diagonal entries using orthogonal reflections. For a nonzero real column segment $x$, set $\sigma =1$ when $x_1\ge 0$ and $\sigma =-1$ otherwise, and defineThen $H^{T}H=I$ and $Hx=-\sigma \Vert x{\Vert }_2{e}_1$. The sign choice avoids subtracting nearly equal numbers when constructing $v$. If the segment is zero, skip the reflection.

For $A\in {\mathbb{R}}^{m\times n}$ with $m\ge n$, apply the first reflection to all rows, the second to rows $2,\dots ,m$, and so on, extending each by an identity on the earlier rows. This giveswhere $Q$ is orthogonal and $R$ consists of an $n\times n$ upper triangular matrix above $m-n$ zero rows. The reflections are stored through their defining vectors, so a full dense $Q$ need not be formed.^[ref]

For full column rank, the top $n\times n$ block of $R$ is invertible and back substitution solves the transformed least squares problem. Orthogonal transformations preserve the Euclidean norm, making this route useful when the normal equations would square a large condition number.

Near an active constraint with limiting Lagrange multiplier ${\lambda }_i^{\ast }>0$, a convergent barrier sequence has $c_i(x_k)\sim {\mu }_k/{\lambda }_i^{\ast }$. Its Hessian matrix therefore contains the termFor a nonzero constraint gradient, this gives large curvature in a direction normal to the boundary. If a nonzero common tangent direction retains positive curvature of order one, while the remaining terms stay bounded, the condition number typically grows like $1/{\mu }_k$. This is a mechanism rather than a universal rate: there may be no tangent direction, all curvatures may grow together, or an active constraint may have zero limiting Lagrange multiplier.

An inverse problem reconstructs an unknown object $x$ from observed data $y$ and a specified model $A(x)=y$. The associated forward problem predicts the observations from $x$. In image deblurring, $x$ contains sharp pixel values and a matrix $A$ mixes neighbouring pixels; in tomography, the observations are integrals through the object. The appeal is that indirect measurements can reveal otherwise inaccessible structure. The difficulty is that the forward problem may erase precisely the distinctions we want to recover.

Consider this completely computable two-pixel image deblurring model:Each measured pixel retains $55\%$ of its own value and receives $45\%$ of its neighbour's. Adding measurement noise $e=(0.01,-0.01{)}^T$ gives $y^{\delta }=(0.56,0.44{)}^T$, whose direct reconstruction is $A^{-1}{y}^{\delta }=(1.1,-0.1{)}^T$. The Euclidean norm of the reconstruction error is $0.1\sqrt{2}$, ten times the data error $0.01\sqrt{2}$. A negative second pixel appears even though both observed pixels are positive.

Despite this wrong reconstruction, its residual of an inverse problem is exactly zero: the predicted data equal the noisy observations. Thus matching the data perfectly does not certify reconstruction accuracy. The relevant question is how large a change in the unknown could hide within the measurement uncertainty. A condition number quantifies that sensitivity for an invertible matrix; the null space identifies changes that no improvement in precision can reveal.

The exact singular value decomposition explains this amplification. With the orthonormal basisthe mean mode survives unchanged while the difference mode shrinks tenfold. Direct inversion multiplies difference-mode measurement noise by ten. The condition number is ${\kappa }_2(A)={\sigma }_{\max }/{\sigma }_{\min }=10$. The beauty of the singular value decomposition is this separation into independent channels: it identifies which information survives and which needs protection by regularization of an inverse problem.