Lecture 24: Rank and proof

The last lecture ended with a sketch whose points the solver put somewhere so that every residual is below one part in a billion of the sketch’s size. Why is that not good enough?

Take a rectangle extruded into a box. The solver reports its top line as Horizontal to the solver’s tolerance: the two ends’ heights differ by perhaps $10^{-10}$. Built as though exact, that line tilts a face by a hair where the designer meant square. It is a real feature, and a boolean with another box later finds it: a sliver of tiny faces where two planes that should be one cross at an angle of $10^{-10}$. Enclosures are supposed to save us from this, and an interval around “the height” could. But around which value, and how wide? The answer has to come from the constraints themselves.

So we ask for something stronger than a small residual: a box of numbers, proven to contain a true solution of the constraints. The geometry built from the box then carries exactly the uncertainty the solve left, and the interval tests downstream can see that the two faces may be one. This lecture is about how to obtain that proof.

Newton’s method, in intervals

We want to show that $f(x)=0$ has a solution in a box $X\subset\mathbb{R}^n$, for $f:\mathbb{R}^n\to\mathbb{R}^n$. For the moment let the system be square and regular. The naive idea is to run Newton with interval numbers. It fails, because Newton divides by $J$, and $J$ over a box is a whole family of matrices.

Krawczyk’s idea (1969) is to not invert an interval matrix at all. Pick any matrix $Y$, a free choice that should be good: the inverse of the Jacobian at the middle of the box. Pick a sharp point $\tilde x$ in $X$. Then define

$$g(x)=x-Yf(x).$$

A zero of $f$ is a fixed point of $g$. Now apply the mean value theorem to each component of $f$: for any $x\in X$, there is a matrix $\bar J$, each row the gradient of that row of $f$ at a point of the segment from $\tilde x$ to $x$, with

$$f(x)=f(\tilde x)+\bar J\,(x-\tilde x).$$

Every row of $\bar J$ lies in the row of the interval Jacobian $J(X)$, so $\bar J\in J(X)$. Substituting,

$$g(x)=\tilde x-Yf(\tilde x)+(I-Y\bar J)(x-\tilde x)\ \in\ \underbrace{\tilde x-Yf(\tilde x)+\bigl(I-YJ(X)\bigr)(X-\tilde x)}_{K(X)}.$$

$K(X)$ is the Krawczyk operator. It is computed with plain interval arithmetic: one evaluation of $f$ at a point, one of $J$ over the box, and a matrix product. What we have is a theorem with three parts.

(a) No solution is lost. If $x\in X$ has $f(x)=0$, then $x=g(x)\in K(X)$. So every root in $X$ lies in $K(X)$. In particular, if $K(X)\cap X=\emptyset$, then $X$ holds no root.

(b) Existence. If $K(X)\subseteq X$, then $g$ maps the compact convex box $X$ into itself, and it is continuous. Brouwer’s fixed point theorem gives a fixed point: a zero of $Yf$ in $X$.

(c) Uniqueness and regularity. Suppose moreover that $\lVert I-YJ(X)\rVert_\infty<1$. For every matrix $A$ in $J(X)$, $\lVert I-YA\rVert<1$ makes $YA$ regular, hence $Y$ and $A$ regular. So the zero of $Yf$ from (b) is a zero of $f$, and if $x,y$ were two of them, then $0=f(x)-f(y)=\bar J(x-y)$ with regular $\bar J$, so $x=y$. And “regular over the whole box” says the independent rows stay independent over the whole box.

So $K(X)\subseteq X$ together with the norm bound proves that $X$ contains exactly one solution of $f$. Rounding is no threat: with outward rounding, the computed $K(X)$ contains the true one, so if the computed $K(X)$ fits inside $X$, the true one does too. The proof runs on a computer.

Note what $Y$ is. It is a free choice. It does not enter the conclusions at all; only whether the test passes depends on it. $Y$ need not be accurate, and it should be computed once and carried as an honest enclosure, or as plain numbers treated as exact. And $\tilde x$ is a free choice too.

A worked example

Take the unit circle and the line $y=0.6$:

$$f(x,y)=\begin{pmatrix}x^2+y^2-1\\ y-0.6\end{pmatrix},\qquad J(x,y)=\begin{pmatrix}2x&2y\\0&1\end{pmatrix}.$$

Let $\tilde x=(0.8,0.6)$, a root, and take the box of half width $0.1$. Then

$$J(\tilde x)=\begin{pmatrix}1.6&1.2\\0&1\end{pmatrix},\qquad Y=\begin{pmatrix}0.625&-0.75\\0&1\end{pmatrix},\qquad J(X)=\begin{pmatrix}[1.4,1.8]&[1.0,1.4]\\0&1\end{pmatrix}.$$

Multiplying,

$$I-YJ(X)=\begin{pmatrix}[-0.125,0.125]&[-0.125,0.125]\\0&0\end{pmatrix},$$

with norm $0.25<1$. Since $f(\tilde x)=0$,

$$K(X)=\tilde x+\bigl(I-YJ(X)\bigr)\begin{pmatrix}[-0.1,0.1]\\ [-0.1,0.1]\end{pmatrix}=\begin{pmatrix}0.8+[-0.025,0.025]\\0.6\end{pmatrix}.$$

$K(X)=[0.775,0.825]\times\{0.6\}$, well inside $X$. The test passes. In this box there is one point on both the circle and the line. Notice how much it contracted, by a factor of 4 from the box: the bound $\lVert I-YJ(X)\rVert$ is the contraction. It is proportional to the box’s width, since $J(X)$ deviates from $J(\tilde x)$ by that much.

How wide should the box be?

There is a Goldilocks zone. Too small a box and $K(X)$ falls outside it, because $\tilde x-Yf(\tilde x)$ is a whole Newton step away from $\tilde x$. Too big and the norm exceeds one, and the spread blows $K(X)$ past the box. So: polish first and then widen. Solving::prove in core/geop-core-math/src/solvers/system/enclose.rs does exactly this:

1. Newton first. At most POLISH_STEPS (6) steps, each squaring the error, on the determined variables. Each iterate is sharp, being only a seed for the next. The last pass takes no step: it is the point the box is proven around.
2. Start from what the step still reaches. The first candidate half-width is twice $\lvert Yf(\tilde x)\rvert$, per variable.
3. Grow only as needed. If $K\not\subseteq$ box, the new radius is what $K$ reached. If the radius is no longer finite, the system is singular. At most ENCLOSE_ATTEMPTS (24) tries: how hard it tries, never what a verified box means.
4. Check the norm. If $K\subseteq$ box but the norm bound fails, the rows are not proven independent, and that is an error, not a success.
5. Return $K\cap$ box. Both enclose the solution, so their intersection is honest and tighter. The rule of Lecture 5 again: intersect two enclosures of one value.

In the demo, after “Newton polish” with the line at height 0.6, “Widen until proven” reports a half-width of $10^{-14}$ at the very first attempt. That width says how precisely the constraints pin down the solution at this scale: a few rounding errors. A well-conditioned sketch comes out with vertices of width near the machine precision; a badly conditioned one comes out wider, honestly.

Formally, the kernel writes the operator for a box $\Delta$ of steps around $\tilde x$,

$$K(\Delta)=-Yf(\tilde x)+\bigl(I-YJ(\Delta)\bigr)\Delta,$$

the same thing with the origin moved to $\tilde x$. That move is not cosmetic. Part VIII needs it, since a pose’s step $\Delta$ is not added to the pose but retracted along its own tangent space (Lecture 36). $J(\Delta)$ is the Jacobian at each value the box reaches, times the differential of the retraction there (Variable::differential). For a number the differential is the identity.

Non-square systems: which rows, which variables?

So far $f$ was square. A real sketch has $m$ rows and $n$ unknowns, usually not equal. Two things must be decided before the proof can start: which rows are independent, and which variables they determine. The rest of the variables are free, the designer’s choice, and they stay exactly as they are. The rest of the rows are left out, and have to be checked afterwards.

eliminate in core/geop-core-math/src/solvers/linear.rs decides this by Gauss-Jordan with complete pivoting (Lecture 23). A pivot is the largest entry that is definitely not zero, among the rows and columns not used yet. Elimination stops when every candidate could be zero. Each pivot (column, row) assigns one variable to one row. The result is a square system on the determined variables, with the free ones held.

So everything hangs on a rank decision. And now comes the central difficulty of this lecture.

A redundant row looks independent at a point

Take a rectangle’s fourth right angle, or an arc concentric with a center it already has. These constraints say what the others say, but they say it on the solution only. Where the blocks are now, a solution merely to the solver’s tolerance, the row’s remainder after elimination is not exactly zero. It has a pivot about as small as the tolerance: say $10^{-10}$. As an interval that pivot may very well exclude zero.

Elimination at the point would then count the row as independent. It takes a variable for it. The square system now has too many rows for what the constraints really determine; its solution set is a curve rather than a point, and nothing can prove a curve is a point. The old code rejected such sketches, so the sketch tools had to avoid redundancy by hand. Worse, as AGENTS.md records, it proved an “exact solution” for two lengths $10^{-11}$ apart.

The fix, from Decide a rank over the box you prove, not at a point: the rows are chosen over a box around the point, widened from the point itself through the solver’s tolerance. A pivot that could be zero anywhere in the box is no pivot. Once the box reaches the true solution, a redundant row’s pivot straddles zero, and the row is left out.

Concretely, enclose_group tries in order: the point itself; a box as wide as the tolerance (tolerance, one part in RELATIVE_TOLERANCE of the size); then twice as wide, and so on. Each distinct choice of pivots is put to the Krawczyk test. The first that passes wins.

Why try the point first? Because when no redundancy is present, the point’s pivots are the best ones. And the box must not start too large, or a real pivot of a tiny feature straddles zero too.

The rank must be checked again, over the proven box

Here is the trap in that fix, and it is a very good one. Take a point on the unit circle and on its tangent $y=1$:

$$f(x,y)=\begin{pmatrix}x^2+y^2-1\\ y-1\end{pmatrix}.$$

There is exactly one solution, $(0,1)$, and it is singular: at it, both rows have gradient $(0,\cdot)$, and $J$ has rank 1. Let the solve end at $x=1.9\times10^{-8}$, as AGENTS.md records. Eliminating over a box as wide as the solver’s tolerance: the circle row’s $\partial/\partial y=2$ pivots on $y$, and the tangent row’s remainder is about $-x$ in $x$, a pivot that could be zero over that box. So the row is left out, and $x$ counts as free and is held where the solve left it. The box then proves $y$ by the circle alone: a beautiful proof that $y=\sqrt{1-x^2}$, for the $x$ we are holding.

And the left-out tangent row, checked over the box? Its value is about $y-1\approx-x^2/2$, of size $2\times10^{-16}$, less than the box’s width. It encloses zero. Everything seems consistent, and a “true solution” $(1.9\times10^{-8},1)$ is returned, where there is none.

The catch: the row’s derivative over the proven box is about $-1.9\times10^{-8}$ in $x$, a very small number that is definitely not zero. Over this box, unlike the wider one, the row is independent. A constraint the others truly imply leaves the rank as it is wherever they all hold. This one does not. So:

A left-out row must vanish over the proven box, and it must add nothing to the rank there.

Solving::verify runs this second elimination over the proven box, with the rule that if there are more pivots than picked rows, then “the constraints are singular at their solution”. The test a_tangency_is_not_proven_in_every_scalar in core/geop-core-math/src/solvers/system/tests.rs drives this at $x=0.3,\ 10^{-4},\ 10^{-6},\ 10^{-8}$: no box is claimed, and no conflict either.

The demo’s “Tangent” button shows the deeper reason no box can pass. With $\tilde x=(x,1)$,

$$J(\tilde x)=\begin{pmatrix}2x&2\\0&1\end{pmatrix},\quad Y=\frac{1}{2x}\begin{pmatrix}1&-2\\0&2x\end{pmatrix},\quad I-YJ(X)=\begin{pmatrix}\mp r/x&\mp r/x\\0&0\end{pmatrix}.$$

(using a box of half-width $r$, whose $2(x\pm r)$ and $2(1\pm r)$ are the first row of $J(X)$). So $\lVert I-YJ(X)\rVert=2r/x$, and regularity requires $r<x/2$. But the only root is $x$ away from $\tilde x$, so a box containing it needs $r\ge x$. Those are incompatible: however close the solve ends, no box both contains the root and proves it regular. The demo reproduces this: at $x=0.01$, half-width $0.005$ gives a norm of $1.00$; half-width $0.001$ passes the norm test, and correctly proves that box empty.

Three things that look alike

We can now give precise meanings to the three cases that all present as “an extra row”:

Getting a conflict right takes care of its own. A left-out row that does not vanish over the proven box is not yet a conflict, since the free variables were held: other values of them might satisfy it. So Solving::conflict looks for a real proof. The combination $g$ of the rows that elimination reduced the row to vanishes wherever all the rows do. Over a box $B$ in every variable, by the mean value theorem, $g\in g(x)+\nabla g(B)\,B$. If this is definitely not zero, no point of $B$ is a solution. Only that gets the name of conflict, with the box’s size stated: “cannot hold where the others do within $\rho$ of here”. The test a_line_missing_a_circle_conflicts_in_every_scalar shows the unit circle with $y=1.5$.

And here is the epistemic heart of the matter. Interval arithmetic can prove a conflict, never its absence. A row that encloses zero agrees with the others to every digit the box resolves. It might still be off by $10^{-30}$. So do not say “implied”; say “consistent to what the box resolves”. This is not a weakness to engineer away. In a floating-point world, two lengths equal to the last digit and two lengths that differ in the 40th are the same input.

Enclosing what a value stands for

When the variables are poses, a new trap opens. A pose is a quaternion and a position, seven numbers for six degrees of freedom. A solver’s step sharpens both, so the quaternion comes out a hair off unit. AGENTS.md records what happened when a solved revolute chain was enclosed as $\tilde x\circ\exp(\delta)$ for a box of twists $\delta$: every pose in the box kept that norm. Turning a point by a quaternion off unit scales the point by its norm squared. So no pose in the box was a rigid motion. The true solution was in none of them. A joint’s redundant row came out $10^{-15}$ off zero over the proven box: a “conflict” between constraints that agree.

retract_over in core/geop-core-math/src/primitives/pose.rs now makes the rotation unit honestly, by an interval normalization, before it steps. The general lesson: wherever a representation has more numbers than the thing it stands for (a quaternion for a rotation, homogeneous coordinates for a point), ask whether an enclosure covers the thing, not just the numbers.

Recap

Exercises

1. For the line $y=0.6$ and the unit circle, take $\tilde x=(0.9,0.5)$ and the half-width $r=0.2$. Compute $Y$, $f(\tilde x)$, $K(X)$ and the norm by hand, and confirm that the box is proven. Check it against the demo.

2. Prove part (a) of the theorem without Brouwer: any root in $X$ lies in $K(X)$. Why is this true even when $Y$ is a terrible choice? What happens to the use of this fact when $Y$ is terrible?

3. Two lines are constrained Parallel and then Horizontal each. At the sketch’s solved point the parallel row’s remainder is $3\times10^{-11}$. Describe what elimination at the point does, what elimination over the box of half-width $10^{-9}\ell$ does, and which of them the Krawczyk test can accept.

4. Why is a proof over the box of steps of a pose not simply a proof over a box of seven numbers (quaternion and position)? Give two reasons.

5. Suppose a sketch’s solve ends well, the Krawczyk test passes, and SolveReport::dof says 2. A user adds a dimension that is exactly redundant with those already there, as computed by hand in decimal. What should happen? Why can the solver not distinguish that from a dimension differing by $10^{-20}$?

6. In the demo, set the line to a height of $1.2$. No root exists. Which verdicts can you obtain by moving the centre and the box? Can it ever claim one root? Why not?