Codex @codex  0

The starting point of a conforming finite element method is a weak formulation in a Hilbert space $X$. For example, a homogeneous Dirichlet boundary condition on a bounded Lipschitz domain gives $X=H_0^1(\Omega )$. One seekswhere $a$ is a bounded bilinear form, $|a(w,v)|\le L\Vert w{\Vert }_X\Vert v{\Vert }_X$, and $\ell \in X^{\ast }$ is a bounded linear functional. A conforming finite element space $X_h\subset X$ is finite dimensional and is typically built from functions that are polynomial on each element of a finite element mesh. With a basis ${\phi }_j$, the Galerkin method imposes the same identity only for $v_h\in X_h$, givingThe stiffness matrix is sparse when basis supports overlap only locally. Local element integrals assemble its entries; the mass matrix represents an $L^2$ term. This separates the choice of trial functions from the variational principle determining their coefficients.

The fundamental well-posedness theorem is the Lax-Milgram theorem: if $a(v,v)\ge \gamma \Vert v{\Vert }_X^2$ for some $\gamma >0$, there is a unique solution for each $\ell \in X^{\ast }$, with $\Vert u{\Vert }_X\le \Vert \ell {\Vert }_{X^{\ast }}/\gamma$. The same theorem on $X_h$ gives a unique discrete solution. The constants must be independent of the mesh for uniform stability of a numerical method. Boundary conditions are built into the trial space when they are essential constraints; flux conditions arise through boundary terms in integration by parts and are natural constraints in the weak formulation.

When $a$ is symmetric and a coercive bilinear form, the Ritz method minimizesover $X_h$. Its first variation is exactly the Galerkin method, so Ritz minimization and Galerkin testing coincide for symmetric coercive problems. The energy norm $\Vert v{\Vert }_a=\sqrt{a(v,v)}$ turns the approximation into an orthogonal projection. Indeed, subtracting the exact and discrete weak equations gives Galerkin orthogonality, $a(u-u_h,v_h)=0$. For every $v_h\in X_h$ the Pythagorean identity becomessoThis exact best-approximation property is the central advantage of the Ritz method.

Symmetry is unnecessary for the Galerkin method. For any bounded, coercive form, the Céa lemma givesTo prove it, put $e=u-u_h$. Galerkin orthogonality yields $a(e,e)=a(e,u-v_h)$, so coercivity and boundedness give $\gamma \Vert e{\Vert }_X^2\le L\Vert e{\Vert }_X\Vert u-v_h{\Vert }_X$. Divide when $e\ne 0$ and take the infimum. Thus approximation capability plus uniform coercivity gives convergence. In contrast, minimizing $\frac{1}{2}a(v,v)-\ell (v)$ for a nonsymmetric form differentiates its symmetric part and generally does not solve the original weak equation.

A concrete symmetric example is the one-dimensional reaction-diffusion problem from Question 5. Continuous piecewise-linear hat functions give the diffusion stiffness matrix $S$ with diagonal $2/h$ and neighbouring entries $-1/h$, and the mass matrix $M$ with diagonal $2h/3$ and neighbouring entries $h/6$. Therefore $K=S+M$ and $F_i=i{h}^2$. Its positive-definite matrix is the coordinate form of the continuous energy. The exact smooth solution $u(x)=x-\sinh (x)/\sinh 1$ makes this an explicit test of the method, not just an abstract existence result.

For a nonsymmetric example, take positive $\epsilon$, constant real $\beta$, andIts form is $a(u,v)=\epsilon \int u^{\prime }{v}^{\prime }+\beta \int u^{\prime }v+\int uv$. Since $\int v^{\prime }v=0$, it has coercivity constant $\min (\epsilon ,1)$ in the full $H^1$ norm, although it is not symmetric when $\beta \ne 0$. The Galerkin method remains well posed. On the uniform hat basis, $K=\epsilon S+\beta C+M$, where $C_{i,i+1}=1/2$, $C_{i,i-1}=-1/2$ and $C_{ii}=0$. The skew-symmetric matrix $C$ cancels from $c^{T}Kc$, proving nonsingularity. A Ritz energy would omit precisely this advection contribution. Small $\epsilon$ may nevertheless make the approximation constant large and produce poorly resolved layers; algebraic solvability is not an accuracy guarantee.

Quantitative convergence uses a finite element interpolation estimate. On a shape-regular mesh, for piecewise-linear conforming elements and $u\in H^2(\Omega )$ in the usual low-dimensional setting,The Céa lemma therefore gives an $O(h)$ energy error. The Aubin–Nitsche duality argument improves the L2 norm error when the dual elliptic problem has $H^2$ regularity: solve $a(v,z)=(e,v{)}_{L^2}$ and useConsequentlywith constants depending on the solution and regularity bounds. The second estimate is conditional on dual regularity; corners or insufficient data regularity can reduce the rate.

The same Galerkin method handles evolution equations. For the homogeneous heat equation, the semidiscrete finite element heat equation is $M\dot{c}+Sc=0$. Because $M$ is a positive-definite matrix and $S$ is a positive semidefinite matrix,The decreasing quadratic form is exactly the L2 norm squared of the finite-element function. This separates spatial stability of a numerical method from the subsequent time integrator, just as in the Schrödinger equation example. Conformity, coercivity and approximation estimates together explain convergence.

An $s$-stage Runge-Kutta method for $y^{\prime }=f(t,y)$ is defined byThe Butcher tableau records $A=(a_{ij})$, $b$ and $c$; internal consistency usually sets $c=Ae$, and first-order consistency requires $b^{T}e=1$. A strictly lower triangular $A$ gives an explicit method, while other stage dependencies generally require an implicit Runge-Kutta method. The stage solvability of an implicit Runge-Kutta method must be checked separately: if $f$ is globally Lipschitz continuous in $y$ with constant $L$, the stage fixed-point map is a contraction in the maximum stage norm whenever $hL{\max }_i{\sum }_j|a_{ij}|<1$. This is a sufficient small-step condition, not a restriction intrinsic to the definitions of A-stability or B-stability.

Linear stability begins with the Dahlquist test equation $y^{\prime }=\lambda y$. Solving the stages givesThe linear stability domain consists of test values for which the stage system is well defined and $|R(z)|\le 1$. A-stability means this includes the closed left half-plane. For $\mathrm{Re}\lambda <0$, this gives stability with no scalar decay-mode step restriction. For a normal matrix $L$ in $y^{\prime }=Ly$, unitary diagonalization of a normal matrix reduces the norm estimate to these scalar factors. For a non-normal matrix, eigenvalues alone do not establish a uniform norm bound: eigenvector conditioning and transient amplification matter. Mesh-uniform estimates for discretized PDEs must control operator powers, not only their spectra.

Every consistent explicit Runge-Kutta method has a nonconstant polynomial stability function and therefore cannot be A-stable, since that polynomial is unbounded on the negative real axis. For example, forward Euler method has $R(z)=1+z$, stability disk $|1+z|\le 1$, and restriction $0\le hq\le 2$ for $\lambda =-q<0$. For a pure imaginary mode $\lambda =i\omega \ne 0$, $|1+ih\omega |>1$, so it is unstable for every positive step. The classical fourth-order Runge-Kutta method instead hasIt is stable on the imaginary axis exactly for $|q|\le 2\sqrt{2}$, an illustrative conditional stability range for oscillatory evolution.

Backward Euler method has $R(z)=(1-z{)}^{-1}$ and is A-stable. L-stability adds $R(z)\to 0$ as $|z|\to \infty$ in the left half-plane; backward Euler satisfies this and strongly damps unresolved rapidly decaying modes in a stiff differential equation. The implicit midpoint rule and the trapezoidal rule both have $R(z)=(1+z/2)/(1-z/2)$, so they are A-stable but not L-stable. They preserve the modulus of a pure imaginary test mode, but their stiff-decay limit is $-1$, leaving oscillatory numerical remnants. The three-stage Lobatto IIIA method has the fourth-order rational function found in Question 2 and likewise lacks L-stability. Thus high order and A-stability do not themselves imply efficient stiff damping.

Nonlinear stability measures differences between solutions of the same equation. A one-sided Lipschitz condition isDifferentiating the squared difference and applying the Gronwall inequality gives $\Vert y(t)-\tilde{y}(t)\Vert \le e^{\nu (t-t_0)}\Vert y(t_0)-\tilde{y}(t_0)\Vert$. A dissipative vector field has $\nu \le 0$ and is contractive. A B-stable method preserves this property for every positive step size for which the stages are well defined:This includes time-dependent vector fields when dissipativity holds at each common time argument. Applying it to scalar linear dissipative fields shows that B-stability implies A-stability when the linear stages are well defined; the converse fails.

The standard sufficient theorem is algebraic stability implies B-stability, subject to stage solvability. Define $B=\mathrm{diag}(b)$ and $M=BA+A^{T}B-b{b}^T$. Algebraic stability of a Runge-Kutta method means $b_i\ge 0$ and $M$ is a positive semidefinite matrix. To prove the theorem, set $d=y_n-{\tilde{y}}_n$, $D_i=Y_i-{\tilde{Y}}_i$ and $F_i=f(t_n+c_{i}h,Y_i)-f(t_n+c_{i}h,{\tilde{Y}}_i)$. Expanding the output difference and using $D_i=d+h{\sum }_j{a}_{ij}{F}_j$ yields the Runge-Kutta contractivity identityDissipativity makes the first sum nonpositive, while positive semidefiniteness makes the last quadratic form nonnegative. The latter follows by factoring $M=Q^{T}Q$ and writing it as a sum of squared norms of linear combinations of the $F_i$. This proves the claimed contraction. It is a sufficient theorem; absence of algebraic stability is not, by itself, a proof of failure of B-stability for every representation.

Several examples make the distinctions concrete. For backward Euler, $b=1$, $A=(1)$, so $M=(1)$: it is both L-stable and B-stable. For implicit midpoint, $b=1$, $A=(1/2)$, so $M=(0)$: it is B-stable and preserves squared norms for a skew-Hermitian linear equation, but lacks stiff damping. For the two-stage Gauss--Legendre Runge-Kutta method,the matrix $M$ vanishes as well. The collocation theorem gives order four, and its stability function is exactly the same as the three-stage Lobatto IIIA function. Yet the Gauss method is algebraically stable, whereas that Lobatto method is not. Identical scalar linear stability functions need not imply identical nonlinear stability properties.

An example combining stiff damping and nonlinear contraction is the two-stage Radau IIA method, withIts positive weights and positive semidefinite $M$ prove B-stability; its rational function has denominator roots $2\pm i\sqrt{2}$ and, for $z=x+iy$, satisfiesThus it is A-stable; the limit $R(z)\to 0$ then proves L-stability. The standard Radau IIA collocation theorem gives order $2s-1$, here three. These properties explain its usefulness for stiff nonlinear equations.

Finally, A-stability alone does not imply B-stability. The trapezoidal rule fails B-stability even for the scalar dissipative equation $y^{\prime }=-y^3$. With $h=1$, its step map $T$ is uniquely defined byAt $y=\sqrt{2}$, $T(y)=0$. Implicit differentiation givesThus nearby starting values are separated by approximately twice their original distance after one step, although the exact flow is contractive. This example and the contractivity theorem show why scalar linear analysis, stiff damping and genuinely nonlinear norm estimates answer different stability questions.

Direct multiplication and the identities ${\cos }^2\theta +{\sin }^2\theta =1$ and ${\cosh }^2\theta -{\sinh }^2\theta =1$ giveThus $R$ is an orthogonal matrix, while $S$ preserves the indefinite quadratic form $x_1^2-x_2^2$.

The Hasse theorem for elliptic curves states that every elliptic curve $E/{\mathbb{F}}_q$ satisfiesLet $\pi$ be the Frobenius isogeny and put $a=\mathrm{tr}(\pi )$. Its fixed points are exactly $E({\mathbb{F}}_q)$, so separability of $1-\pi$ givesThe degree of an isogeny is a nonnegative quadratic form on the endomorphism ring, and for all integers $m,n$,If $a^2>4q$, this quadratic polynomial has two real roots and takes a negative value at some rational $m/n$ between them, hence after clearing denominators at some integer pair $(m,n)$. Therefore $a^2\le 4q$, and substituting $a=q+1-\#E({\mathbb{F}}_q)$ proves the bound. This is the degree-form proof of the Hasse bound.

Since $\hat{h}-h$ is bounded, a bound on $\hat{h}(P)$ bounds the projective height of $x(P)$. The Northcott theorem gives only finitely many possible rational $x$-coordinates, and each has at most two points above it. Hence

By the Mordell-Weil theorem, $E(\mathbb{Q})\cong {\mathbb{Z}}^r\oplus E(\mathbb{Q}{)}_{\mathrm{tors}}$. The canonical height is a positive-definite quadratic form on the lattice ${\mathbb{Z}}^r$, so canonical-height lattice-point growth givesIf $r=0$, this count is bounded; if $r=1$, it is $O(\sqrt{B})$. Consequently

The Horn copositive matrix is invariant under cyclic permutation of its five coordinates: its negative entries correspond precisely to neighboring indices on the five-cycle. For any $x\ge 0$, cyclically relabel the coordinates so that $x_5$ is a smallest coordinate. In particular, $x_4-x_5\ge 0$.

Using the stated identity in this coordinate order,The square is nonnegative, and both remaining terms are nonnegative since $x_i\ge 0$ and $x_4\ge x_5$. Cyclic invariance means the relabeling has not changed the quadratic form. Thus $H$ is a copositive matrix for every original ordering of $x$.

Let $w=(1,2,1,0,0{)}^T$. Direct multiplication gives $Hw=(0,0,0,2,2{)}^T$, so $w^{T}Hw=0$.

Suppose $H=P+N$, where $P$ is a positive semidefinite matrix and $N$ is a symmetric nonnegative matrix. Since $w\ge 0$, both $w^{T}Pw$ and $w^{T}Nw$ are nonnegative. Their sum is zero, so both vanish. In particular,Every coefficient $w_i{w}_j$ in this sum is strictly positive and every $N_{ij}$ is nonnegative, forcing $N_{ij}=0$ for $1\le i,j\le 3$.

Now apply the same argument to all five cyclic shifts of $w$. The Horn copositive matrix is cyclically invariant, so every shifted vector also has zero quadratic form. It follows that the entries of $N$ vanish on every cyclic block of three consecutive indices. Every pair of indices on a five-cycle lies in such a block, hence $N=0$.

This would imply $P=H$. But the zero quadratic form of a positive semidefinite matrix would then force $Hw=0$, contradicting $Hw=(0,0,0,2,2{)}^T$. Therefore $H$ lies outside the positive-semidefinite-plus-nonnegative cone. By the previous equivalence, its quartic form is a nonnegative polynomial that is not a sum of squares polynomial.

For algebraic stability, the weights must be nonnegative and the real symmetric matrix $M=BA+A^{T}B-b{b}^T$, where $B=\mathrm{diag}(b)$, must be a positive semidefinite matrix.

The weights are positive, but $a_{11}=0$ and $b_1=1/6$, soA positive semidefinite matrix cannot have a negative diagonal entry, as its quadratic form at the first standard basis vector would be negative. Thus this Lobatto IIIA method is not algebraically stable, despite its A-stability.

A Runge-Kutta method with stage matrix $A$, weights $b$ and step size $h$ isIts order measures accuracy as $h\to 0$, while stability concerns how errors or stiff components propagate. Linear A-stability, strong stiff damping through L-stability, and nonlinear B-stability address different questions.

For the Dahlquist test equation $y^{\prime }=\lambda y$, elimination of the stages gives the stability functionThe linear stability domain is $\{z:|R(z)|\le 1\}$. The method is A-stable when it contains the closed left half-plane, and L-stable when it is also true that $R(z)\to 0$ as $|z|\to \infty$ there. These conditions describe, respectively, stability for every non-growing scalar linear mode and damping of very rapidly decaying modes.

The Forward Euler method has $R(z)=1+z$, so it fails A-stability, for example at $z=-3$. In fact no nontrivial explicit Runge-Kutta method is A-stable: its stability function is a nonconstant polynomial, which is unbounded along the negative real axis. The Backward Euler method has $R(z)=1/(1-z)$. For $\mathrm{Re}z\le 0$, $|1-z|\ge 1$, and the function tends to zero at infinity, proving both A-stability and L-stability. The implicit midpoint rule has $R(z)=(1+z/2)/(1-z/2)$. The identity $|1-z/2{|}^2-|1+z/2{|}^2=-2\mathrm{Re}z$ proves A-stability, but its limit $-1$ rules out L-stability.

For nonlinear equations, a dissipative vector field satisfies $⟨f(x)-f(y),x-y⟩\le 0$. The corresponding exact solutions contract because the time derivative of their squared distance is twice that inner product. A B-stable method preserves this contraction for every $h>0$, whenever its stages are well defined.

Let $d=y_n-{\tilde{y}}_n$, $D_i=Y_i-{\tilde{Y}}_i$, and $F_i=f(Y_i)-f({\tilde{Y}}_i)$. Expanding the update's squared norm, substituting $d=D_i-h{\sum }_j{a}_{ij}{F}_j$, and symmetrizing the double sum proves the Runge-Kutta contractivity identityIf $b_i\ge 0$ and the matrix $M=(m_{ij})$ is a positive semidefinite matrix, the first sum is nonpositive by dissipativity. The second double sum is nonnegative: sum the quadratic form of $M$ over each coordinate of the $F_i$. Thus algebraic stability implies B-stability. Applied to complex scalar linear equations with the real part of the Hermitian inner product, the same argument also implies A-stability, where the stages are defined.

The Backward Euler method has $b_1=a_{11}=1$, hence $M=(1)$, and is algebraically stable. One can prove its nonlinear contraction directly: $d_{n+1}=d+h(f(y_{n+1})-f({\tilde{y}}_{n+1}))$ gives $\Vert d_{n+1}{\Vert }^2\le ⟨d,d_{n+1}⟩\le \Vert d\Vert \Vert d_{n+1}\Vert$. The implicit midpoint rule has $a_{11}=1/2$, $b_1=1$ and $M=0$, so it is B-stable even though it is not L-stable.

Linear A-stability does not by itself imply nonlinear B-stability. The ODE trapezoidal rule has the same stability function as the implicit midpoint rule, but take the dissipative scalar field $f(y)=-\max (y,0)$ and $h=6$. For any positive starting value, its trapezoidal update is $y_{n+1}=y_n+3(f(y_n)+f(y_{n+1}))=-2y_n<0$. Thus two positive starting values have their distance doubled, violating B-stability. Its algebraic stability matrix also has a negative first diagonal entry. The field is globally Lipschitz continuous, so this is a well-defined nonlinear counterexample.

These examples distinguish accuracy, scalar stiff-mode stability, nonlinear contraction, and rapid stiff-mode damping. For implicit methods, solvability of the stage equations remains an additional requirement; a formal stability function alone does not establish it for arbitrary nonlinear problems.

The symmetric matrix of the quadratic form in the standard basis isbecause the off-diagonal entries contribute twice to $x^{T}Ax$.

Apply Gram-Schmidt orthogonalization for a symmetric bilinear form toThese vectors are pairwise orthogonal for $B(u,v)=u^{T}Av$, andThus the form is diagonal in the basis $(u_1,u_2,u_3)$, with matrixIt follows that its rank is $3$ and its signature is

By the real spectral theorem, the real symmetric matrix $A$ has an orthonormal basis of eigenvectors. In that basis the associated quadratic form isRescaling the coordinates belonging to nonzero eigenvalues changes every positive coefficient to $+1$ and every negative coefficient to $-1$. Hence the diagonal normal form has one positive square for each positive eigenvalue and one negative square for each negative eigenvalue. By Sylvester's law of inertia, these counts do not depend on the diagonalizing basis, so

The numerical eigenvalues are not invariant under a general change of basis, because the matrix changes by congruence rather than similarity. For example, on a one-dimensional space let $q(e)=1$. Its matrix in the basis $(e)$ is $(1)$, whereas in the basis $(2e)$ it is $(4)$. The eigenvalue changes from $1$ to $4$, although its sign, and therefore the signature, is unchanged.

Forthe chord-and-tangent group law gives, for distinct nonopposite points $P=(x_1,y_1)$ and $Q=(x_2,y_2)$,For doubling, replace the slope byThe point at infinity is the identity and $-(x,y)=(x,-y)$.

Let the $x$-coordinates of $P,Q,P+Q,P-Q$ be $x_1,x_2,x_3,x_4$. Applying the addition formula with $Q$ and $-Q$ givesSubstituting $y_i^2=x_i^3+a{x}_i+b$ and simplifying yieldsMultiplying the two addition formulas and eliminating $y_1{y}_2$ in the same way gives

These identities are the algebraic source of two parallelogram laws. Applied to pullbacks of the pole divisor of $x$, they implyfor isogenies of elliptic curves. Together with $\mathrm{deg}(n\phi )=n^2\mathrm{deg}\phi$, this makes the degree a quadratic form. Applied to the Absolute logarithmic Weil height of the four $x$-coordinates, with bounded terms removed by passage to the limit defining the canonical height of an elliptic curve, they similarly give

Let $f_1,f_2\in V^{\ast }$ satisfy $f_1(v)f_2(v)=0$ for every vector $v$. If $f_1=0$ there is nothing to prove. Otherwise choose $u$ with $f_1(u)\ne 0$. For every $v\in \ker f_1$ and every real number $t$,so the hypothesis gives $f_2(u+tv)=0$. Taking two values of $t$ shows that $f_2(u)=f_2(v)=0$. Sincewe obtain $f_2=0$. Thus one of the two linear functionals is zero.

Sylvester's law of inertia says that, in a suitable basis, every real quadratic form isThe rank is $r=p+m$, and the signature is $\sigma =p-m$; both are independent of the chosen basis.

If $q=f_1{f}_2$, its polar form has image contained in $\mathrm{span}\{f_1,f_2\}$, so $r\le 2$. If $r=2$, the two functionals must be independent and, after taking their sum and difference, $q$ is a difference of two squares. Hence $\sigma =0$. If $r=1$, then $|\sigma |=1$. In every case,Conversely, this inequality leaves onlyup to a linear change of coordinates. These factor respectively as $0$, $(\pm x_1)x_1$, and $(x_1-x_2)(x_1+x_2)$. Therefore

Finally suppose that $q$ takes both positive and negative values. Sylvester's law supplies normalized basis vectorswith $q(e_i)=1$, $q(f_j)=-1$, and the $z_l$ spanning the radical. The $p+m+k$ vectorsare linearly independent and all satisfy $q(v)=0$. They form the required basis of isotropic vectors.

The likelihood function and prior distribution give the posterior densityCompleting the square in the quadratic form givesThus Gaussian conjugacy for a normal linear model yields

This symmetric finite-rank integral operator is not positive. For $f(t)=1-\frac{3}{2}t$, set $A_j={\int }_0^1{t}^{j}f(t)dt$. Directly,The associated quadratic form isEvery covariance operator is a positive operator, so this operator is not a covariance operator.

The symmetric bilinear form associated with a real quadratic form $Q$ is obtained by polarization identity:In coordinates $Q(x)=x^{T}Ax$, replacing $A$ by its symmetric part does not change $Q$, and the formula gives $phi(u,v)=u^{T}Av$; this proves existence. A symmetric form is positive semidefinite when $phi(v,v)\ge 0$ for every $v$, and positive definite when the inequality is strict for every $v\ne 0$.

The diagonalization theorem for real quadratic forms says that some basis puts any symmetric form intowith $r$ further zero coordinates. Sylvester's law of inertia says that $(p,q,r)$ is independent of the diagonalizing basis. To prove this, let $P$ be the span of the positive coordinate vectors and let $N\oplus Z$ be the span of the negative and zero vectors. If $P^{\prime }$ is positive definite for another diagonalization and $\dim P^{\prime }>p$, thenso the two spaces intersect nontrivially. A vector in the intersection would have both positive and nonpositive square, a contradiction. Thus $p^{\prime }\le p$; symmetry gives $p=p^{\prime }$. Applying the same argument to $-\varphi$ gives $q=q^{\prime }$, and then $r=r^{\prime }$.

For the nondegenerate form on $V$, write its inertia as $(p,q,0)$, so $p+q=2n$. The restriction vanishes identically on $E$: polarization gives $phi(u,v)=0$ for $u,v\in E$. Projection of $E$ to the positive coordinate space is injective, since a vector with zero positive projection cannot be isotropic unless it is zero. Hence $k\le p$; projection to the negative space likewise gives $k\le q$. Therefore $k\le \min (p,q)\le n$.

The matrix of $l^2$ is $l{l}^T$. If $l\ne 0$, it has rank one and inertia $(1,0,n-1)$, hence signature one; if $l=0$, its rank and signature are zero. The coefficientwise product $(l^2,s^2)$ has matrix $t{t}^T$, where $t_i=l_i{s}_i$, so it has the same rank-one conclusion when $t\ne 0$ and is zero otherwise.

Finally, diagonalization writes every positive semidefinite form as a sum of squares, say $f={\sum }_a{l}_a^2$ and $g={\sum }_b{s}_b^2$. Bilinearity of coefficientwise multiplication givesa sum of positive semidefinite rank-at-most-one forms. Thus $(f,g)$ is positive semidefinite.

For $h(x)=e_p(q(x))$, the product in the $U^3$ cube average isThe expression in parentheses is the third additive derivative of the quadratic form $q$, so it vanishes. Every cube contributes one and therefore the quadratic phase satisfies

Apply the coordinate-replacement proof from part (ii) directly to the quadratic formThe zero diagonal makes $f$ multilinear. For coordinate $i$,and the row and column assumptions implyEvery hybrid vector appearing during replacement has independent centered variance-one coordinates with fourth moments at most $9$. The degree-one case of part (i) therefore givesThe fourth-order Taylor remainder from replacing coordinate $i$ is at most $(M/12)\mathbb{E}(D_{i}f{)}^4\le 12M$. Summing the $n$ replacement errors yields the stronger estimate

Write $h=x-x_0$, $k=y-y_0$, and let $a=f_{xx}$, $b=f_{xy}$ and $c=f_{yy}$ be evaluated at the critical point. Assuming $f$ is twice continuously differentiable nearby, the Taylor expansion isThe linear terms vanish because the gradient is zero. For $a>0$, completing the square givesConsequently the requested constant is the Hessian determinantThe displayed identity in the question is an identity for the quadratic part; for a general $f$, its right-hand side needs the remainder $o(h^2+k^2)$.