Least Squares and Quadratic Programs

Fitting a model to data and computing contact forces in a physics engine are the same problem wearing different clothes. Both minimize a convex quadratic, and the machinery that solves one solves the other. Once you see the shape, a least squares fit, a control step, and a stack of boxes resting on the floor all read as instances of a single object.

Least squares is the smallest quadratic program

Overdetermined systems have no exact solution, so we settle for the $x$ that minimizes the residual

$$\min_{x}\ \lVert Ax - b\rVert^{2} = \min_{x}\ x^{\top}\!A^{\top}\!A\,x - 2\,b^{\top}\!A\,x + b^{\top}b .$$

Setting the gradient to zero gives the normal equations

$$A^{\top}\!A\,x = A^{\top}b .$$

This is already a quadratic in $x$ with Hessian $A^{\top}A$, which is symmetric positive semidefinite for any $A$. Least squares is an unconstrained quadratic program, solvable in closed form. Add inequalities and the closed form is gone, but the object stays convex.¹

The general quadratic program

A quadratic program (QP) minimizes a quadratic objective subject to linear constraints:

$$\min_{x}\ \tfrac{1}{2}\,x^{\top}\!Q\,x + c^{\top}x \quad\text{s.t.}\quad G x \le h .$$

The single fact that decides everything is the definiteness of $Q$. If $Q$ is symmetric positive semidefinite, the objective is convex, the feasible set is a convex polyhedron, and any local minimum is global. If $Q$ is indefinite the problem is NP-hard. Physics and control live entirely in the convex case, where $Q = Q^{\top} \succeq 0$ by construction.

A contact solve is a QP

Stack the active contact constraints into a Jacobian $J$ mapping generalized velocity to relative velocity at the contacts. The impulses $\lambda$ that resolve the contacts minimize a quadratic in exactly the QP form

$$\min_{\lambda \in \mathcal{K}}\ \tfrac{1}{2}\,\lambda^{\top}\!A\,\lambda + b^{\top}\lambda, \qquad A = J\,M^{-1}J^{\top}.$$

Here $A$ is the Delassus operator, the inverse inertia seen in contact space. It is symmetric positive semidefinite, so this is a convex QP. The role of $Q$ is played by $A$, and the cone $\mathcal{K}$ carries the physics: a frictionless contact can push but not pull, so $\lambda \ge 0$, a box (orthant) constraint, and this reproduces the linear complementarity condition $0 \le \lambda \perp A\lambda - b \ge 0$. With friction $\mathcal{K}$ becomes a second-order (Coulomb) cone. Modern engines lean into this: MuJoCo solves a convex program each step,² and the SAP solver recasts compliant contact as an unconstrained convex minimization so a Newton method applies directly.³ The same $\tfrac{1}{2}\lambda^{\top}A\lambda + b^{\top}\lambda$ is also the inner problem of a model-predictive control step, which is why one QP solver serves both the simulator and the controller.

A box-constrained QP by hand

Take the separable case first. With

$$Q = \begin{bmatrix} 2 & 0 \\ 0 & 1 \end{bmatrix}, \qquad c = \begin{bmatrix} -2 \\ -1 \end{bmatrix},$$

the unconstrained minimum is $x^{*} = -Q^{-1}c = [\,1,\ 1\,]$, with objective $f(x^{*}) = -1.5$. Now impose the box $0 \le x \le 0.5$. Because $Q$ is diagonal, the objective splits into independent single-variable parabolas, so the constrained minimum is just the per-coordinate clamp:

$$\hat{x} = \operatorname{clamp}\!\big([1,1],\,0,\,0.5\big) = [\,0.5,\ 0.5\,], \qquad f(\hat{x}) = -1.125 .$$

The catch: clamping is exact only because $Q$ is diagonal. Give $Q$ an off-diagonal term and the coordinates couple, so pushing $x_0$ to its bound shifts the optimal $x_1$. The clamp of $x^{*}$ is then feasible but suboptimal, and there is no closed form. You must iterate.

KKT and the active set

At the optimum the Karush-Kuhn-Tucker conditions hold: stationarity $Q\hat{x} + c + G^{\top}\mu = 0$, primal feasibility $G\hat{x} \le h$, dual feasibility $\mu \ge 0$, and complementary slackness $\mu_i (G\hat{x}-h)_i = 0$. The last line says each constraint is either tight with $\mu_i \ge 0$ or slack with $\mu_i = 0$. The set of tight constraints is the active set. If you knew it, the QP would collapse to an equality-constrained problem with a closed-form solution; active-set methods guess it and repair the guess. Projected gradient descent is the first-order cousin: step downhill, then project back onto the feasible box,

$$x \leftarrow \operatorname{clamp}\!\big(x - \alpha\,(Q x + c),\ l,\ u\big),$$

and the clamp automatically discovers which bounds are active as the iterate settles onto them.

Projected gradient descent on a box QPinteractive

c₀ = -3.00c₁ = -1.00

PauseReset

Q

[2 0.6 ; 0.6 1]

iterate k

0 / 40

x* (unconstrained)

[1.46, 0.12]

x̂ (constrained)

[1.00, 0.40]

f(x̂)

-2.08

clamp(x*)

[1.00, 0.12]

active box constraints: x₀ = 1 (upper). When any are active, clamp(x*) differs from x̂: the coupling Q₀₁ makes clamping wrong.

Footnotes

1. Boyd & Vandenberghe, Convex Optimization (2004). ↩

2. MuJoCo docs, Computation. ↩

3. Castro, Permenter, Han, An Unconstrained Convex Formulation of Compliant Contact (2021). ↩