Numerical Methods

This page covers the numerical foundations: how residuals and Jacobians are computed, how the sparse linear algebra works, and what convergence criteria the solver uses.

Analytical Jacobian vs Numerical

The solver computes analytical Jacobians for all 47+ constraint types. Each constraint type has a hand-derived partial derivative implementation in jacobian.rs:

// PointOnCircle: residual = (px-cx)² + (py-cy)² − r²
// ∂r/∂px = 2(px−cx),  ∂r/∂cx = −2(px−cx),  ∂r/∂r = −2r

A numerical fallback using central differences is available for debugging:

∂r_i/∂x_j ≈ (r_i(x_j + h) − r_i(x_j − h)) / (2h),   h = 1e-7

The proptest suite (50 random cases per constraint type) verifies that analytical and numerical Jacobians match within tolerance: absolute < 1e-4, relative < 1e-2.

Residual Formulation

Each constraint contributes one or more scalar residuals to the vector r(x). Examples:

Angular constraints require special treatment due to branch ambiguity (e.g., θ vs θ + π). The solver detects the nearest branch and uses that for the residual.

Normal Equations (LM)

For the Levenberg-Marquardt algorithm, the normal equations are:

(JᵀJ + λ · diag(JᵀJ)) · Δx = −Jᵀr

The damping term λ · diag(JᵀJ) serves two purposes:

1. Positive definiteness — ensures the system is solvable even when JᵀJ is singular
2. Step control — large λ → small steps (gradient descent), small λ → large steps (Gauss-Newton)

The initial λ is set to 1e-3 · max(diag(JᵀJ)) and adapted based on step acceptance.

Sparse Jacobian (CSC Format)

For large problems, the Jacobian is stored in Compressed Sparse Column (CSC) format:

J = [col_ptr, row_indices, values]

Each constraint touches at most 6–8 parameters, giving a very sparse Jacobian. For a 200-parameter problem with 100 constraints, the Jacobian has at most 800 nonzeros out of 20,000 entries (4% density).

SpGEMM (sparse matrix–matrix multiplication) computes A = JᵀJ without forming dense intermediates. The result is another sparse matrix used by the Conjugate Gradient solver.

Convergence Criteria

The solver uses three convergence checks:

1. Residual tolerance: ‖r‖₂ < tol (default: 1e-6)
2. Step tolerance: ‖Δx‖₂ < xtol (default: 1e-10)
3. Gradient tolerance: ‖Jᵀr‖∞ < gtol (default: 1e-8)

Any one of these is sufficient to declare convergence. The solver also enforces a maximum iteration count (default: 100) to prevent infinite loops.

Radius Clamping and Radmin

Circle and arc radii are clamped to radmin = 1e-12 to prevent degenerate geometry:

• During solving, if r < radmin, the radius is clamped and a warning is emitted
• This prevents division by zero in tangent and point-on-circle residuals
• The clamping is reversible: the unclamped value is restored on solve failure

Warm Start and Drag Mode

In interactive CAD applications, the solver is called repeatedly as the user drags geometry. Warm start uses the previous solution as the initial guess:

• Typically converges in 2–5 iterations (vs 20–50 from a cold start)
• The parameter vector is preserved between calls
• Only changed constraints trigger recomputation of the Jacobian structure

Drag mode pins the dragged entity's position as a constraint, ensuring the solver respects the user's input.

Underconstrained Handling

When DOF > 0, the system is underconstrained. Two strategies are available:

1. SoftAnchor — adds weak springs k·(x − x₀) that resist deviation from the initial guess
2. MinimumNorm — uses SVD to find the minimum-norm solution among all solutions satisfying the constraints

MinimumNorm is more expensive (O(n³) for the SVD) but produces smoother motion in interactive scenarios.

Parameter Reduction

The normalize.rs module eliminates constraints that can be expressed as parameter substitutions:

After reduction, the parameter vector is compacted and a mapping is maintained to expand the solution back to the original parameters.

Constraint Propagation

Twenty deduction rules in solver/propagation.rs infer additional constraints from existing ones:

• Horizontal ∧ Vertical → Perpendicular
• Coincident ∧ Horizontal → SameY
• EqualRadius ∧ Concentric → SameCircle

These deductions are applied before solving to tighten the constraint graph and improve convergence.