By 苏剑林 | October 03, 2026
For unconstrained smooth functions, the universal method for finding extrema is naturally to take the derivative (gradient) and set it to zero. What about extrema problems with equality constraints? Many readers are likely no stranger to this; the standard method is the Lagrange multiplier method. However, Lagrange multipliers often feel like a pure algebraic trick, lacking the intuitive feel of perturbation analysis found in unconstrained optimization.
In this article, let's try this: if we apply the perturbation thinking of unconstrained extrema to constrained optimization, what results will we get?
Stationary Point Conditions
First, let's revisit why a gradient equal to zero is a necessary condition (stationary point condition) for unconstrained extrema. Consider the objective function $f(\boldsymbol{x})$, where $\boldsymbol{x}\in\mathbb{R}^d$. We want to find the minimum point of $f(\boldsymbol{x})$. We make a first-order approximation near $\boldsymbol{x}$:
\begin{equation}f(\boldsymbol{x} + \delta\boldsymbol{x}) = f(\boldsymbol{x}) + \delta\boldsymbol{x}^{\top} \nabla_{\boldsymbol{x}} f(\boldsymbol{x})\end{equation}
If $f(\boldsymbol{x})$ is already a minimum, then for any perturbation $\delta\boldsymbol{x}\neq \boldsymbol{0}$, it must hold that $f(\boldsymbol{x} + \delta\boldsymbol{x})\geq f(\boldsymbol{x})$, which implies $\delta\boldsymbol{x}^{\top} \nabla_{\boldsymbol{x}} f(\boldsymbol{x})\geq 0$. However, $\delta\boldsymbol{x}^{\top} \nabla_{\boldsymbol{x}} f(\boldsymbol{x})$ is a linear function with respect to $\delta\boldsymbol{x}$. For it to be non-negative for all $\delta\boldsymbol{x}$, the only choice is:
\begin{equation}\nabla_{\boldsymbol{x}} f(\boldsymbol{x})=\boldsymbol{0}\end{equation}
The same logic applies to maxima, so this is the common stationary point condition for both maxima and minima. When the stationary point condition holds, the next leading term in the expansion is:
\begin{equation}f(\boldsymbol{x} + \delta\boldsymbol{x}) = f(\boldsymbol{x}) + \frac{1}{2}\delta\boldsymbol{x}^{\top}\nabla_{\boldsymbol{x}}^2 f(\boldsymbol{x}) \,\delta\boldsymbol{x}\end{equation}
According to the same principle, for $f(\boldsymbol{x} + \delta\boldsymbol{x})\geq f(\boldsymbol{x})$ to hold consistently, then $\delta\boldsymbol{x}^{\top}\nabla_{\boldsymbol{x}}^2 f(\boldsymbol{x}) \,\delta\boldsymbol{x}\geq 0$ must always hold, meaning the Hessian matrix $\nabla_{\boldsymbol{x}}^2 f(\boldsymbol{x})$ is positive semi-definite. This is a necessary condition for a minimum point. Conversely, a necessary condition for a maximum point is that $\nabla_{\boldsymbol{x}}^2 f(\boldsymbol{x})$ is negative semi-definite.
Of course, this discussion is not a rigorous mathematical proof; it is more of a heuristic derivation.
Equality Constraints
Now let's see what happens to the perturbation process when equality constraints are introduced. Consider:
\begin{equation}\min_{\boldsymbol{x}} f(\boldsymbol{x})\qquad\text{s.t.}\qquad \boldsymbol{g}(\boldsymbol{x})=\boldsymbol{0}\end{equation}
where $\boldsymbol{g}(\boldsymbol{x})\in\mathbb{R}^k$, representing $k$ equality constraints. Pick any $\boldsymbol{x}$ that satisfies $\boldsymbol{g}(\boldsymbol{x})=\boldsymbol{0}$. Similarly, considering the first-order approximation of the perturbation gives:
\begin{equation}f(\boldsymbol{x} + \delta\boldsymbol{x}) = f(\boldsymbol{x}) + \delta\boldsymbol{x}^{\top} \nabla_{\boldsymbol{x}} f(\boldsymbol{x})\end{equation}
Similarly, the necessary condition for $f(\boldsymbol{x})$ to be a minimum is still that $\delta\boldsymbol{x}^{\top} \nabla_{\boldsymbol{x}} f(\boldsymbol{x})\geq 0$ holds consistently. But the difference is that now it doesn't need to hold for any perturbation $\delta\boldsymbol{x}$, but only for all $\delta\boldsymbol{x}$ that satisfy $\boldsymbol{g}(\boldsymbol{x} + \delta\boldsymbol{x})=\boldsymbol{0}$. Applying a first-order approximation to $\boldsymbol{g}(\boldsymbol{x} + \delta\boldsymbol{x})=\boldsymbol{0}$, we get:
\begin{equation}\boldsymbol{J}\delta\boldsymbol{x} = \boldsymbol{0},\qquad \boldsymbol{J} = \nabla_{\boldsymbol{x}} \boldsymbol{g}(\boldsymbol{x})\in\mathbb{R}^{k\times d}\end{equation}
Note that $\boldsymbol{J}\delta\boldsymbol{x} = \boldsymbol{0}$ is a system of linear equations with respect to $\delta\boldsymbol{x}$. We can write its general solution using the "pseudoinverse":
\begin{equation}\delta\boldsymbol{x} = (\boldsymbol{I} - \boldsymbol{J}^{\dagger}\boldsymbol{J})\delta\boldsymbol{z}\end{equation}
where $\boldsymbol{J}^{\dagger}$ is the pseudoinverse of $\boldsymbol{J}$, and $\delta\boldsymbol{z}$ represents any unconstrained perturbation vector. When $\boldsymbol{J}\boldsymbol{J}^{\top}$ is invertible, $\boldsymbol{J}^{\dagger}=\boldsymbol{J}^{\top}(\boldsymbol{J}\boldsymbol{J}^{\top})^{-1}$, and from this it is clear that $\boldsymbol{J}^{\dagger}\boldsymbol{J}$ is a symmetric matrix. Substituting this into $\delta\boldsymbol{x}^{\top} \nabla_{\boldsymbol{x}} f(\boldsymbol{x})\geq 0$, and by the arbitrariness of $\delta\boldsymbol{z}$, we obtain the stationary point condition:
\begin{equation}(\boldsymbol{I} - \boldsymbol{J}^{\dagger} \boldsymbol{J}) \nabla_{\boldsymbol{x}} f(\boldsymbol{x}) = \boldsymbol{0} \label{eq:stable-g}\end{equation}
where $(\boldsymbol{I} - \boldsymbol{J}^{\dagger} \boldsymbol{J}) \nabla_{\boldsymbol{x}} f(\boldsymbol{x})$ is called the projected gradient of $f(\boldsymbol{x})$ on the manifold $\boldsymbol{g}(\boldsymbol{x})=\boldsymbol{0}$. It is the orthogonal projection of the Euclidean gradient $\nabla_{\boldsymbol{x}} f(\boldsymbol{x})$ onto the tangent space. In the language of Riemannian geometry, it is also called the Riemannian gradient.
Lagrange
Now let's compare this with the Lagrange multiplier method. The Lagrange multiplier method introduces undetermined coefficients $\boldsymbol{c}\in\mathbb{R}^k$ and modifies the objective function to:
\begin{equation}F(\boldsymbol{x}, \boldsymbol{c}) = f(\boldsymbol{x}) - \boldsymbol{c}^{\top}\boldsymbol{g}(\boldsymbol{x})\end{equation}
Then we seek the stationary point of the unconstrained objective $F(\boldsymbol{x}, \boldsymbol{c})$. Setting the gradients on both sides to zero gives:
\begin{equation}\nabla_{\boldsymbol{x}} F(\boldsymbol{x}, \boldsymbol{c}) = \nabla_{\boldsymbol{x}} f(\boldsymbol{x}) - \boldsymbol{J}^{\top}\boldsymbol{c} = \boldsymbol{0},\qquad \nabla_{\boldsymbol{c}} F(\boldsymbol{x}, \boldsymbol{c}) = -\boldsymbol{g}(\boldsymbol{x}) = \boldsymbol{0}\end{equation}
Note that for the pseudoinverse, $\boldsymbol{J}^{\dagger} \boldsymbol{J}=\boldsymbol{J}^{\top}(\boldsymbol{J}^{\top})^{\dagger}$ holds. Therefore, equation $\eqref{eq:stable-g}$ can also be written as $\nabla_{\boldsymbol{x}} f(\boldsymbol{x}) - \boldsymbol{J}^{\top}(\boldsymbol{J}^{\top})^{\dagger}\nabla_{\boldsymbol{x}} f(\boldsymbol{x}) = \boldsymbol{0}$. If we denote $(\boldsymbol{J}^{\top})^{\dagger}\nabla_{\boldsymbol{x}} f(\boldsymbol{x})$ as $\boldsymbol{c}$, it becomes exactly the same as the above equation. That is, equation $\eqref{eq:stable-g}$ can derive the same conditions as the Lagrange multiplier method.
Conversely, can equation $\eqref{eq:stable-g}$ be derived from the Lagrange multiplier method? Naturally, yes. Assuming $k < d$ and that the rank of $\boldsymbol{J}$ is $k$, then viewing $\nabla_{\boldsymbol{x}} f(\boldsymbol{x}) - \boldsymbol{J}^{\top}\boldsymbol{c} = \boldsymbol{0}$ as an equation for $\boldsymbol{c}$, it is an overdetermined system. If it has a solution, it must be $\boldsymbol{c}=(\boldsymbol{J}^{\top})^{\dagger}\nabla_{\boldsymbol{x}} f(\boldsymbol{x})$. Substituting this back into the equation yields equation $\eqref{eq:stable-g}$.
This shows that the conditions provided by both are mathematically equivalent. Of course, the Lagrange multiplier method looks cleaner in its writing, but the author believes the perturbation idea is more consistent with the unconstrained case and thus prefers this understanding.
Furthermore, the simplicity of the Lagrange multiplier method is only formal; understanding it truly is no easier than the perturbation perspective. Readers can view the process from the perturbation idea to the projected gradient as a derivation of the Lagrange multiplier method.
Two Examples
We deepen our understanding through the following two examples.
Sphere Constraint
For the first example, consider a simple extrema on a sphere:
\begin{equation}\min_{\boldsymbol{x}} f(\boldsymbol{x})\qquad\text{s.t.}\qquad \Vert\boldsymbol{x}\Vert_2=1\end{equation}
It is easy to find that $\nabla_{\boldsymbol{x}}\Vert\boldsymbol{x}\Vert_2 = \boldsymbol{x}/\Vert\boldsymbol{x}\Vert_2 = \boldsymbol{x}$, so $\boldsymbol{J} = \boldsymbol{x}^{\top}$. Since $\boldsymbol{x}$ is a unit vector, $\boldsymbol{J}^{\dagger} = \boldsymbol{x}$, therefore the stationary point condition is:
\begin{equation}\boldsymbol{0} = (\boldsymbol{I} - \boldsymbol{x}\boldsymbol{x}^{\top})\nabla_{\boldsymbol{x}} f(\boldsymbol{x}) = \nabla_{\boldsymbol{x}} f(\boldsymbol{x}) - \boldsymbol{x}(\boldsymbol{x}^{\top}\nabla_{\boldsymbol{x}} f(\boldsymbol{x}))\end{equation}
The final equation is precisely the projection of the gradient onto the tangent space of the sphere.
Orthogonal Constraint
Another example is the optimization problem on the Stiefel manifold:
\begin{equation}\min_{\boldsymbol{X}} f(\boldsymbol{X})\qquad\text{s.t.}\qquad \boldsymbol{X}^{\top}\boldsymbol{X} = \boldsymbol{I}_m,\quad\boldsymbol{X}\in\mathbb{R}^{n\times m}\,(n\geq m)\end{equation}
If we just want to find the stationary point condition, using the Lagrange multiplier method would be much faster, but here we will derive the projected gradient in full. The difficulty here lies in the fact that $\boldsymbol{X}$ itself is a matrix, so $\boldsymbol{J}$ would be a high-order matrix, which increases the mental burden.
To address this, we use the linear operator perspective. Taking the first-order approximation of $\boldsymbol{X}^{\top}\boldsymbol{X} - \boldsymbol{I}_m$ gives $\delta\boldsymbol{X}^{\top}\boldsymbol{X} + \boldsymbol{X}^{\top}\delta\boldsymbol{X}$. According to the previous derivation, this should equal $\boldsymbol{J}\delta\boldsymbol{x}$. We directly define the linear operator:
\begin{equation}\boldsymbol{J}[\boldsymbol{V}] = \boldsymbol{V}^{\top}\boldsymbol{X} + \boldsymbol{X}^{\top}\boldsymbol{V}\end{equation}
Next, we need to derive $\boldsymbol{J}^{\top}[\boldsymbol{V}]$. This is understood according to $\langle\boldsymbol{U}, \boldsymbol{J}^{\top}[\boldsymbol{V}]\rangle_F=\langle\boldsymbol{V}, \boldsymbol{J}[\boldsymbol{U}]\rangle_F$ (analogous to $\boldsymbol{u}^{\top}\boldsymbol{J}^{\top}\boldsymbol{v}=\boldsymbol{v}^{\top}\boldsymbol{J}\boldsymbol{u}$). We already know $\boldsymbol{J}[\boldsymbol{U}]$ and $\langle\cdot,\cdot\rangle_F$, so we have:
\begin{equation}\newcommand{tr}{\mathop{\text{tr}}}\langle\boldsymbol{U}, \boldsymbol{J}^{\top}[\boldsymbol{V}]\rangle_F = \tr(\boldsymbol{V}^{\top}(\boldsymbol{U}^{\top}\boldsymbol{X} + \boldsymbol{X}^{\top}\boldsymbol{U})) = \tr(\boldsymbol{U}^{\top}\boldsymbol{X}(\boldsymbol{V} + \boldsymbol{V}^{\top})) = \langle\boldsymbol{U}, \boldsymbol{X}(\boldsymbol{V} + \boldsymbol{V}^{\top})\rangle_F \end{equation}
From this, we read $\boldsymbol{J}^{\top}[\boldsymbol{V}] = \boldsymbol{X}(\boldsymbol{V} + \boldsymbol{V}^{\top})$. According to $\boldsymbol{J}^{\dagger}\boldsymbol{J}=\boldsymbol{J}^{\top}(\boldsymbol{J}\boldsymbol{J}^{\top})^{-1}\boldsymbol{J}$, we also need $\boldsymbol{J}\boldsymbol{J}^{\top}$ and $(\boldsymbol{J}\boldsymbol{J}^{\top})^{-1}$:
\begin{equation}\boldsymbol{J}\boldsymbol{J}^{\top}[\boldsymbol{V}] = \boldsymbol{J}[\boldsymbol{X}(\boldsymbol{V} + \boldsymbol{V}^{\top})] = (\boldsymbol{V} + \boldsymbol{V}^{\top})\boldsymbol{X}^{\top}\boldsymbol{X} + \boldsymbol{X}^{\top}\boldsymbol{X}(\boldsymbol{V} + \boldsymbol{V}^{\top}) = 2(\boldsymbol{V} + \boldsymbol{V}^{\top})\end{equation}
If $\boldsymbol{V}$ is a symmetric matrix, then $\boldsymbol{J}\boldsymbol{J}^{\top}[\boldsymbol{V}]=4\boldsymbol{V}$, so $(\boldsymbol{J}\boldsymbol{J}^{\top})^{-1}[\boldsymbol{V}]=\boldsymbol{V}/4$. Since $\boldsymbol{J}[\boldsymbol{V}]$ is coincidentally symmetric, we only need this specific case. Summing it up, we have:
\begin{equation}\boldsymbol{J}^{\dagger}\boldsymbol{J}[\boldsymbol{V}] = \boldsymbol{J}^{\top}(\boldsymbol{J}\boldsymbol{J}^{\top})^{-1}\boldsymbol{J}[\boldsymbol{V}] = \boldsymbol{J}^{\top}[(\boldsymbol{V}^{\top}\boldsymbol{X} + \boldsymbol{X}^{\top}\boldsymbol{V})/4] = \boldsymbol{X}(\boldsymbol{V}^{\top}\boldsymbol{X} + \boldsymbol{X}^{\top}\boldsymbol{V})/2\end{equation}
Let $\boldsymbol{G} = \nabla_{\boldsymbol{X}}f(\boldsymbol{X})$. Then its projected gradient on the Stiefel manifold is:
\begin{equation}\boldsymbol{G} - \boldsymbol{J}^{\dagger}\boldsymbol{J}[\boldsymbol{G}] = \boldsymbol{G} - \boldsymbol{X}[\boldsymbol{X}^{\top}\boldsymbol{G}]_{\text{sym}}\end{equation}
where $[\boldsymbol{Z}]_{\text{sym}} = (\boldsymbol{Z} + \boldsymbol{Z}^{\top})/2$ is the symmetrization operator. Setting it to zero and multiplying both sides by $\boldsymbol{X}^{\top}$ from the left gives $\boldsymbol{X}^{\top}\boldsymbol{G} = [\boldsymbol{X}^{\top}\boldsymbol{G}]_{\text{sym}}$. This means the stationary condition is that $\boldsymbol{X}^{\top}\boldsymbol{G}$ itself must be a symmetric matrix!
Second-Order Conditions
Is there a way to judge the minimum/maximum under equality constraints similar to checking if $\nabla_{\boldsymbol{x}}^2 f(\boldsymbol{x})$ is positive/negative definite? Naturally, there is, and there is a very clever derivation. Similar to the Lagrange multiplier method, we define:
\begin{equation}F(\boldsymbol{x}, \boldsymbol{c}) = f(\boldsymbol{x}) - \boldsymbol{c}^{\top}\boldsymbol{g}(\boldsymbol{x})\end{equation}
For any stationary point $\boldsymbol{x}$, we automatically have $\boldsymbol{g}(\boldsymbol{x})=\boldsymbol{0}$, so for any $\boldsymbol{c}\in\mathbb{R}^k$, it holds that $F(\boldsymbol{x}, \boldsymbol{c})\equiv f(\boldsymbol{x})$. Thus, expanding to the second order gives:
\begin{equation}f(\boldsymbol{x} + \delta \boldsymbol{x}) = F(\boldsymbol{x} + \delta\boldsymbol{x}, \boldsymbol{c}) = \underbrace{F(\boldsymbol{x}, \boldsymbol{c})}_{\equiv f(\boldsymbol{x})} + \delta \boldsymbol{x}^{\top}(\nabla_{\boldsymbol{x}}f(\boldsymbol{x}) - \boldsymbol{J}^{\top} \boldsymbol{c}) + \frac{1}{2}\delta \boldsymbol{x}^{\top}\boldsymbol{H}(\boldsymbol{c})\delta \boldsymbol{x}\end{equation}
Here $\boldsymbol{H}[\boldsymbol{c}] \triangleq \nabla_{\boldsymbol{x}}^2 (f(\boldsymbol{x}) - \boldsymbol{c}^{\top}\boldsymbol{g}(\boldsymbol{x}))$. Due to the arbitrariness of $\boldsymbol{c}$, we can directly take $\boldsymbol{c}=(\boldsymbol{J}^{\top})^{\dagger}\nabla_{\boldsymbol{x}} f(\boldsymbol{x})$, making the first-order term zero, hence:
\begin{equation}f(\boldsymbol{x} + \delta \boldsymbol{x}) = f(\boldsymbol{x}) + \frac{1}{2}\delta \boldsymbol{x}^{\top}\boldsymbol{H}[(\boldsymbol{J}^{\top})^{\dagger}\nabla_{\boldsymbol{x}} f(\boldsymbol{x})]\delta \boldsymbol{x}\end{equation}
Now only the second-order term remains, so the approximation accuracy of $\delta \boldsymbol{x}$ only needs to be first-order, i.e., $\delta\boldsymbol{x} = (\boldsymbol{I} - \boldsymbol{J}^{\dagger}\boldsymbol{J})\delta\boldsymbol{z}$ is sufficient. Substituting this in gives:
\begin{equation}f(\boldsymbol{x} + \delta \boldsymbol{x}) = f(\boldsymbol{x}) + \frac{1}{2}\delta \boldsymbol{z}^{\top}(\boldsymbol{I} - \boldsymbol{J}^{\dagger}\boldsymbol{J})\boldsymbol{H}[(\boldsymbol{J}^{\top})^{\dagger}\nabla_{\boldsymbol{x}} f(\boldsymbol{x})](\boldsymbol{I} - \boldsymbol{J}^{\dagger}\boldsymbol{J})\delta\boldsymbol{z}\end{equation}
From the arbitrariness of $\delta \boldsymbol{z}$, the necessary condition for $f(\boldsymbol{x})$ to be a minimum is that $(\boldsymbol{I} - \boldsymbol{J}^{\dagger}\boldsymbol{J})\boldsymbol{H}[(\boldsymbol{J}^{\top})^{\dagger}\nabla_{\boldsymbol{x}} f(\boldsymbol{x})](\boldsymbol{I} - \boldsymbol{J}^{\dagger}\boldsymbol{J})$ is positive semi-definite. Similarly, for a maximum, the matrix must be negative semi-definite.
Of course, whether in unconstrained or constrained optimization, judging the maximum/minimum through the positive/negative definiteness of the Hessian matrix usually only has theoretical value. In practice, one typically filters for all positions that meet the stationary point condition and then analyzes each situation individually.
Summary
This article applied the perturbation analysis idea of unconstrained extrema to equality constrained optimization, derived the projected gradient in the tangent space, and established that the stationary point condition is setting the corresponding projected gradient to zero. Finally, we proved its equivalence to the conclusions of the Lagrange multiplier method.