diff --git a/desc/objectives/utils.py b/desc/objectives/utils.py
index 4609c86a43..27fe6fd15d 100644
--- a/desc/objectives/utils.py
+++ b/desc/objectives/utils.py
@@ -80,13 +80,14 @@ def factorize_linear_constraints(objective, constraint, x_scale="auto"):  # noqa
 
     if isinstance(objective, ProximalProjection):
         # remove cols of A corresponding to ["R_lmn", "Z_lmn", "L_lmn", "Ra_n", "Za_n"]
-        # see desc.optimize._constraint_wrappers.ProximalProjection._set_eq_state_vector
+        # see desc.optimize._constraint_wrappers.ProximalState._set_eq_state_vector
         c = 0
         cols = np.array([], dtype=int)
+        eq = objective._state.eq
         for t in objective.things:
-            if t is objective._eq:
-                for arg, dim in objective._eq.dimensions.items():
-                    if arg in objective._args:  # these Equilibrium args are kept
+            if t is eq:
+                for arg, dim in eq.dimensions.items():
+                    if arg in objective._state.args:  # these Equilibrium args are kept
                         cols = np.append(cols, np.arange(c, c + dim))
                     c += dim  # other Equilibrium args are removed
             else:  # non-Equilibrium args are always included
diff --git a/desc/optimize/__init__.py b/desc/optimize/__init__.py
index 7093f03c4a..b7013593cb 100644
--- a/desc/optimize/__init__.py
+++ b/desc/optimize/__init__.py
@@ -1,7 +1,11 @@
 """Functions for minimization and wrappers for scipy methods."""
 
 from . import _desc_wrappers, _scipy_wrappers
-from ._constraint_wrappers import LinearConstraintProjection, ProximalProjection
+from ._constraint_wrappers import (
+    LinearConstraintProjection,
+    ProximalProjection,
+    ProximalState,
+)
 from .aug_lagrangian import fmin_auglag
 from .aug_lagrangian_ls import lsq_auglag
 from .fmin_scalar import fmintr
diff --git a/desc/optimize/_constraint_wrappers.py b/desc/optimize/_constraint_wrappers.py
index 1041b2a1fb..233d36335e 100644
--- a/desc/optimize/_constraint_wrappers.py
+++ b/desc/optimize/_constraint_wrappers.py
@@ -5,7 +5,6 @@
 import numpy as np
 
 from desc.backend import jit, jnp, put
-from desc.batching import batched_vectorize
 from desc.objectives import (
     BoundaryRSelfConsistency,
     BoundaryZSelfConsistency,
@@ -572,6 +571,9 @@ class ProximalProjection(ObjectiveFunction):
     perturb_options, solve_options : dict
         dictionary of arguments passed to Equilibrium.perturb and Equilibrium.solve
         during the projection step.
+    state: ProximalState, optional
+        State manager for the equilibrium constraints. Default is None, in which
+        case the state is created.
     name : str
         Name of the objective function.
     """
@@ -579,109 +581,27 @@ class ProximalProjection(ObjectiveFunction):
     def __init__(
         self,
         objective,
-        constraint,
-        eq,
+        constraint=None,
+        eq=None,
         perturb_options=None,
         solve_options=None,
+        state=None,
         name="ProximalProjection",
     ):
         assert isinstance(objective, ObjectiveFunction), (
             "objective should be instance of ObjectiveFunction." ""
         )
-        assert isinstance(constraint, ObjectiveFunction), (
-            "constraint should be instance of ObjectiveFunction." ""
-        )
-        for con in constraint.objectives:
-            errorif(
-                not con._equilibrium,
-                ValueError,
-                "ProximalProjection method cannot handle general "
-                + f"nonlinear constraint {con}.",
-            )
-            # can't have bounds on constraint bc if constraint is satisfied then
-            # Fx == 0, and that messes with Gx @ Fx^-1 Fc etc.
-            errorif(
-                con.bounds is not None,
-                ValueError,
-                "ProximalProjection can only handle equality constraints, "
-                + f"got bounds for constraint {con}",
-            )
         self._objective = objective
-        self._constraint = constraint
-        solve_options = {} if solve_options is None else solve_options
-        self._solve_during_proximal_build = solve_options.pop(
-            "solve_during_proximal_build", True
-        )  # If user does not want the solve during build, mainly for debug purposes
-        perturb_options = {} if perturb_options is None else perturb_options
-        perturb_options.setdefault("verbose", 0)
-        perturb_options.setdefault("include_f", False)
-        solve_options.setdefault("verbose", 0)
-        self._perturb_options = perturb_options
-        self._solve_options = solve_options
+        if state is None:
+            self._state = ProximalState(eq, constraint, perturb_options, solve_options)
+        else:
+            self._state = state
         self._built = False
         # don't want to compile this, just use the compiled objective and constraint
         self._use_jit = False
         self._compiled = False
-        self._eq = eq
         self._name = name
 
-    def _set_eq_state_vector(self):
-        full_args = self._eq.optimizable_params.copy()
-        self._args = self._eq.optimizable_params.copy()
-        # the eq optimizable variables for proximal are the Rb, Zb and profile
-        # coefficients. Once these are chosen, we will solve the equilibrium to
-        # find the R_lmn, Z_lmn, L_lmn, Ra_n, Za_n. That is why we remove them
-        # from the list of optimizable variables. This is accompanied by not including
-        # self-consistency constraints (see get_combined_constraint_objectives in
-        # desc.optimize.optimizer) and also removing columns corresponding to these
-        # variables from the constraint matrix A in
-        # desc.objectives.utils.factorize_linear_constraints.
-        for arg in ["R_lmn", "Z_lmn", "L_lmn", "Ra_n", "Za_n"]:
-            self._args.remove(arg)
-
-        self._eq_Z, self._eq_D, self._eq_unfixed_idx = (
-            self._eq_solve_objective._Z,
-            self._eq_solve_objective._D,
-            self._eq_solve_objective._unfixed_idx,
-        )
-
-        dxdc = []
-        xz = {arg: np.zeros(self._eq.dimensions[arg]) for arg in full_args}
-
-        for arg in self._args:
-            if arg not in ["Rb_lmn", "Zb_lmn"]:
-                x_idx = self._eq.x_idx[arg]
-                dxdc.append(np.eye(self._eq.dim_x)[:, x_idx])
-            if arg == "Rb_lmn":
-                c = get_instance(self._eq_linear_constraints, BoundaryRSelfConsistency)
-                # We have A @ R_lmn = Rb_lmn
-                A = c.jac_unscaled(xz)[0]["R_lmn"]
-                Ainv = np.linalg.pinv(A)
-                # Once this is multipled by Rb_lmn, we get the full eq state vector
-                # with the R_lmn but rest is 0
-                dxdRb = np.eye(self._eq.dim_x)[:, self._eq.x_idx["R_lmn"]] @ Ainv
-                dxdc.append(dxdRb)
-            if arg == "Zb_lmn":
-                c = get_instance(self._eq_linear_constraints, BoundaryZSelfConsistency)
-                A = c.jac_unscaled(xz)[0]["Z_lmn"]
-                Ainv = np.linalg.pinv(A)
-                dxdZb = np.eye(self._eq.dim_x)[:, self._eq.x_idx["Z_lmn"]] @ Ainv
-                dxdc.append(dxdZb)
-        # dxdc is a matrix that when multiplied by the optimization variables (only
-        # Rb_lmn, Zb_lmn) gives the full state vector of the equilibrium (Rb_lmn and
-        # Zb_lmn part will be 0, but they will be represented by the equivalent
-        # R_lmn and Z_lmn). For example, let's say the eq optimization variables are
-        # ceq = [Rb_lmn, Zb_lmn, p_l, i_l].T                      # noqa : E800
-        # Then, we will use dxdc for the following:
-        # xeq = dxdc @ ceq                                        # noqa : E800
-        # And xeq will be,
-        # xeq = [                                                 # noqa : E800
-        #     R_lmn, Z_lmn, jnp.zeros_like(L_lmn)                 # noqa : E800
-        #     jnp.zeros_like(Rb_lmn), jnp.zeros_like(Zb_lmn),     # noqa : E800
-        #     p_l, i_l,                                           # noqa : E800
-        # ]                                                       # noqa : E800
-        self._dxdc = jnp.hstack(dxdc)
-
     def build(self, use_jit=None, verbose=1):  # noqa: C901
         """Build the objective.
 
@@ -697,105 +617,48 @@ def build(self, use_jit=None, verbose=1):  # noqa: C901
         timer = Timer()
         timer.start("Proximal projection build")
 
-        self._eq_linear_constraints = get_fixed_boundary_constraints(eq=self._eq)
-        self._eq_linear_constraints = maybe_add_self_consistency(
-            self._eq, self._eq_linear_constraints
-        )
-
         # we don't always build here because in ~all cases the user doesn't interact
         # with this directly, so if the user wants to manually rebuild they should
         # do it before this wrapper is created for them.
         if not self._objective.built:
             self._objective.build(use_jit=use_jit, verbose=verbose)
-        if not self._constraint.built:
-            self._constraint.build(use_jit=use_jit, verbose=verbose)
 
-        for constraint in self._eq_linear_constraints:
-            constraint.build(use_jit=use_jit, verbose=verbose)
-
-        # Here we create and build the LinearConstraintProjection
-        # for the equilibrium subproblem using the self._constraint as objective
-        # and our fixed-bdry constraints we just made. This will
-        # be passed as the objective for the eq subproblem, which saves
-        # some time as by building it here we can avoid re-computing the
-        # constraint matrix A and its SVD for the feasible direction method
-        self._eq_solve_objective = LinearConstraintProjection(
-            self._constraint,
-            ObjectiveFunction(self._eq_linear_constraints),
-            name="Eq Update LinearConstraintProjection",
-        )
-        self._eq_solve_objective.build(use_jit=use_jit, verbose=verbose)
-
-        errorif(
-            self._constraint.things != [self._eq],
-            ValueError,
-            "ProximalProjection can only handle constraints on the equilibrium.",
-        )
-
-        self._objectives = [self._objective, self._constraint]
+        self._state.build(use_jit=use_jit, verbose=verbose)
+        self._objectives = [self._objective, self._state.constraint]
         self._set_things()
 
-        self._eq_idx = self.things.index(self._eq)
-
         self._dim_f = self._objective.dim_f
         if self._dim_f == 1:
             self._scalar = True
         else:
             self._scalar = False
 
-        self._set_eq_state_vector()
-
-        # the full state vector includes all the parameters from all the things
-        # however, sub-objectives only need the part for their thing. We will
-        # use this to split the state vector into its components
-        self._dimx_per_thing = [t.dim_x for t in self.things]
-        # we remove the R_lmn, Z_lmn, L_lmn, Ra_n, Za_n from the equilibrium params
-        # dimc_per_thing accounts for that, don't confuse it with reduced state vector
-        self._dimc_per_thing = [t.dim_x for t in self.things]
-        self._dimc_per_thing[self._eq_idx] = np.sum(
-            [self._eq.dimensions[arg] for arg in self._args]
-        )
-
-        # equivalent matrix for A[unfixed_idx] @ D @ Z == A @ feasible_tangents
-        self._feasible_tangents = jnp.eye(self._objective.dim_x)
-        self._feasible_tangents = jnp.split(
-            self._feasible_tangents, np.cumsum(self._dimx_per_thing), axis=-1
-        )
-        # dg/dxeq_reduced = dg/dx_eq_unscaled @ dx_eq_unscaled/dxeq_reduced # noqa: E800
-        # x_eq_unscaled = Deq(xp_eq + Zeq @ xeq_reduced)                    # noqa: E800
-        # So, the feasible tangents (aka. dx_eq_unscaled/dx_reduced) is Deq@Zeq
-        # Since here we are setting the feasible direction for eq parameters only,
-        # we need to add 0 rows for eq fixed parameters and non-eq parameters which we
-        # handle by below operation
-        self._feasible_tangents[self._eq_idx] = self._feasible_tangents[self._eq_idx][
-            :, self._eq_unfixed_idx
-        ] @ (self._eq_Z * self._eq_D[self._eq_unfixed_idx, None])
-        self._feasible_tangents = jnp.concatenate(
-            [np.atleast_2d(foo) for foo in self._feasible_tangents], axis=-1
-        )
-
-        ## history and caching
-        # first, ensure equilibrium is solved to the
-        # specified tolerances, necessary as we assume
-        # eq is solved when taking the derivatives later
-        if self._solve_during_proximal_build:
-            self._eq.solve(
-                objective=self._eq_solve_objective,
-                constraints=None,
-                **self._solve_options,
-            )
-        # then store the now-solved eq state as the initial state
-        self._x_old = self.x(self.things)
-        self._allx = [self._x_old]
-        self._allxopt = [self._objective.x(*self.things)]
-        self._allxeq = [self._eq.pack_params(self._eq.params_dict)]
-        self.history = [[t.params_dict.copy() for t in self.things]]
-
         self._built = True
         timer.stop("Proximal projection build")
         if verbose > 1:
             timer.disp("Proximal projection build")
 
+    def _set_things(self, things=None):
+        """Assign "things" to the wrapper and underlying objectives.
+
+        Parameters
+        ----------
+        things: list of optimizable objects, optional
+            If None, uses "things" of self._objectives.
+
+        """
+        super()._set_things(things)
+
+        # Sync "things" between the wrapper and objective.
+        # Does not include self._constraint
+        self._objective._set_things(self.things)
+
+        self._eq_idx = self.things.index(self._state.eq)
+        self._dimx_per_thing = [t.dim_x for t in self.things]
+        dimc_per_thing = [t.dim_x for t in self.things]
+        dimc_per_thing[self._eq_idx] = self._state.dim_ceq
+        self._dimc_per_thing = dimc_per_thing
+
     def unpack_state(self, x, per_objective=True):
         """Unpack the state vector into its components.
 
@@ -824,17 +687,23 @@ def unpack_state(self, x, per_objective=True):
                 + f"{self.dim_x} got {x.size}."
             )
 
-        xs = jnp.split(x, np.cumsum(self._dimc_per_thing))
+        xs = jnp.split(x, np.cumsum(self._dimc_per_thing)[:-1])
         params = []
         for t, xi in zip(self.things, xs):
-            if t is self._eq:
-                xi_splits = np.cumsum([self._eq.dimensions[arg] for arg in self._args])
-                p = {arg: xis for arg, xis in zip(self._args, jnp.split(xi, xi_splits))}
+            if t is self._state.eq:
+                xi_splits = np.cumsum(
+                    [self._state.eq.dimensions[arg] for arg in self._state.args]
+                )
+                p = {
+                    arg: xis
+                    for arg, xis in zip(self._state.args, jnp.split(xi, xi_splits))
+                }
                 p.update(  # add in dummy values for missing parameters
                     {
                         arg: jnp.zeros_like(xis)
                         for arg, xis in t.params_dict.items()
-                        if arg not in self._args  # R_lmn, Z_lmn, L_lmn, Ra_n, Za_n
+                        if arg
+                        not in self._state.args  # R_lmn, Z_lmn, L_lmn, Ra_n, Za_n
                     }
                 )
                 params += [p]
@@ -862,10 +731,10 @@ def x(self, *things):
         assert [type(t1) is type(t2) for t1, t2 in zip(things, self.things)]
         xs = []
         for t in self.things:
-            if t is self._eq:
+            if t is self._state.eq:
                 xs += [
                     jnp.concatenate(
-                        [jnp.atleast_1d(t.params_dict[arg]) for arg in self._args]
+                        [jnp.atleast_1d(t.params_dict[arg]) for arg in self._state.args]
                     )
                 ]
             else:
@@ -880,13 +749,7 @@ def dim_x(self):
         Note that we remove the R_lmn, Z_lmn, L_lmn, Ra_n, Za_n from the equilibrium
         params.
         """
-        s = 0
-        for t in self.things:
-            if t is self._eq:
-                s += sum(self._eq.dimensions[arg] for arg in self._args)
-            else:
-                s += t.dim_x
-        return s
+        return np.sum(self._dimc_per_thing)
 
     def _update_equilibrium(self, x, store=False):
         """Update the internal equilibrium with new boundary, profile etc.
@@ -897,67 +760,65 @@ def _update_equilibrium(self, x, store=False):
             New values of the state vector of equilibrium (except R_lmn, Z_lmn,
             L_lmn, Ra_n, Za_n) and all the parameters of the other things.
         store : bool
-            Whether the new x should be stored in self.history
+            Whether the new x is stored as the next accepted iterate.
 
         Notes
         -----
-        After updating, if store=False, self._eq will revert back to the previous
-        solution when store was True
+        After updating, if store=False, self._state.eq will revert back to the previous
+        solution when store was True.
 
         """
         # xopt is the full state vector of all the things
         # xeq is the full state vector of the equilibrium only
 
-        # TODO (#1720): We don't need to check the whole state vector, just the
-        # equilibrium parameters should be enough.
-        # first check if its something we've seen before, if it is just return
+        # first check if it's something we've seen before, if it is just return
         # cached value, no need to perturb + resolve
-        xopt = f_where_x(x, self._allx, self._allxopt)
-        xeq = f_where_x(x, self._allx, self._allxeq)
-        if xopt.size > 0 and xeq.size > 0:
+        xs = np.split(x, np.cumsum(self._dimc_per_thing)[:-1])
+        ceq = xs[self._eq_idx]
+        xeq = f_where_x(ceq, self._state.allceq, self._state.allxeq)
+        if xeq.size > 0:
             pass
         else:
-            # After unpack_state, R_lmn, Z_lmn, L_lmn, Ra_n and Za_n in below lists
-            # will be 0s
-            x_list = self.unpack_state(x, False)
-            x_list_old = self.unpack_state(self._x_old, False)
-            xeq_dict = x_list[self._eq_idx]
-            xeq_dict_old = x_list_old[self._eq_idx]
-            deltas = {str(key): xeq_dict[key] - xeq_dict_old[key] for key in xeq_dict}
+            # build a dictionary of the deltas between xeq and xeq_old,
+            # restricted to state.args
+            ceq_split = jnp.split(ceq, self._state.idx_ceq)
+            ceq_dict = dict(zip(self._state.args, ceq_split))
+            deltas = {
+                arg: ceq_dict[arg] - self._state.xeq_old[arg]
+                for arg in self._state.args
+            }
+            # clear cache to reduce memory
+            self._state._tangents = {}
+            self._state._tangent_xf = None
             # We pass in the LinearConstraintProjection object to skip some redundant
             # computations in the perturb and solve methods
-            self._eq = self._eq.perturb(
-                objective=self._eq_solve_objective,
+            self._state.eq = self._state.eq.perturb(
+                objective=self._state.eq_solve_objective,
                 constraints=None,
                 deltas=deltas,
-                **self._perturb_options,
+                **self._state.perturb_options,
             )
-            self._eq.solve(
-                objective=self._eq_solve_objective,
+            self._state.eq.solve(
+                objective=self._state.eq_solve_objective,
                 constraints=None,
-                **self._solve_options,
+                **self._state.solve_options,
             )
-            xeq = self._eq.pack_params(self._eq.params_dict)
-            x_list[self._eq_idx] = self._eq.params_dict.copy()
-            xopt = jnp.concatenate(
-                [t.pack_params(xi) for t, xi in zip(self.things, x_list)]
-            )
-            self._allx.append(x)
-            self._allxopt.append(xopt)
-            self._allxeq.append(xeq)
+            xeq = self._state.eq.pack_params(self._state.eq.params_dict)
+            self._state.allceq.append(ceq)
+            self._state.allxeq.append(xeq)
+            self._state.eq_is_current = False
 
         if store:
-            self._x_old = x
-            x_list = self.unpack_state(x, False)
-            xeq_dict = self._eq.unpack_params(xeq)
-            self._eq.params_dict = xeq_dict
-            x_list[self._eq_idx] = xeq_dict
-            self.history.append(x_list)
-        else:
+            eq_params = self._state.eq.unpack_params(xeq)
+            self._state.eq.params_dict = eq_params
+            self._state.xeq_old = eq_params
+        elif not self._state.eq_is_current:
             # reset to last good params
-            self._eq.params_dict = self.history[-1][self._eq_idx]
-            self._eq_solve_objective.update_constraint_target(self._eq)
+            self._state.eq.params_dict = self._state.xeq_old
+            self._state.eq_solve_objective.update_constraint_target(self._state.eq)
+        self._state.eq_is_current = True
 
+        xopt = jnp.concatenate([*xs[: self._eq_idx], xeq, *xs[self._eq_idx + 1 :]])
         return xopt, xeq
 
     def compute_scaled(self, x, constants=None):
@@ -1058,24 +919,11 @@ def grad(self, x, constants=None):
         # We are looking for the gradient of L = 0.5 * G.T @ G
         # Then, the gradient is ∇L = G.T @ J_of_G
         # where J_of_G is the Jacobian of G with respect to the optimization variables
-        # We explained getting J_of_G in the _jvp method. It is basically,
-        # J_of_G = ∇G @ [dc_tangents - (∇F @ dx_tangents) ^ -1 @ (∇F @ dc_tangents)]
-        # where ∇G is the Jacobian of G with respect to full state vector
-        # and ∇F is the Jacobian of F with respect to full state vector. Then,
-        # ∇L = G.T @ ∇G @ [dc_tangents - (∇F @ dx_tangents) ^ -1 @ (∇F @ dc_tangents)]
-        # We get the part in [] using the _get_tangent method.
-        v = jnp.eye(x.shape[0])
+        # This is a vjp with G serving as the cotangents.
         constants = setdefault(constants, [None, None])
-        xg, xf = self._update_equilibrium(x, store=True)
-        jvpfun = lambda u: self._get_tangent(u, xf, constants, op="scaled_error")
-        tangents = batched_vectorize(
-            jvpfun,
-            signature="(n)->(k)",
-            chunk_size=self._constraint._jac_chunk_size,
-        )(v)
+        xg, _ = self._update_equilibrium(x, store=True)
         g = self._objective.compute_scaled_error(xg, constants[0])
-        g_vjp = self._objective.vjp_scaled_error(g, xg, constants[0])
-        return tangents @ g_vjp
+        return self._vjp(g, x, constants, "scaled_error")
 
     def hess(self, x, constants=None):
         """Compute Hessian of self.compute_scalar.
@@ -1099,6 +947,10 @@ def hess(self, x, constants=None):
         J = self.jac_scaled_error(x, constants)
         return J.T @ J
 
+    def _jac(self, x, constants=None, op="scaled"):
+        # passing v=None corresponds to jvp in all directions
+        return self._jvp(None, x, constants, op).T
+
     def jac_scaled(self, x, constants=None):
         """Compute Jacobian of self.compute_scaled.
 
@@ -1115,8 +967,7 @@ def jac_scaled(self, x, constants=None):
             Jacobian matrix.
 
         """
-        v = jnp.eye(x.shape[0])
-        return self.jvp_scaled(v, x, constants).T
+        return self._jac(x, constants, "scaled")
 
     def jac_scaled_error(self, x, constants=None):
         """Compute Jacobian of self.compute_scaled_error.
@@ -1134,8 +985,7 @@ def jac_scaled_error(self, x, constants=None):
             Jacobian matrix.
 
         """
-        v = jnp.eye(x.shape[0])
-        return self.jvp_scaled_error(v, x, constants).T
+        return self._jac(x, constants, "scaled_error")
 
     def jac_unscaled(self, x, constants=None):
         """Compute Jacobian of self.compute_unscaled.
@@ -1152,8 +1002,7 @@ def jac_unscaled(self, x, constants=None):
         J : ndarray
             Jacobian matrix.
         """
-        v = jnp.eye(x.shape[0])
-        return self.jvp_unscaled(v, x, constants).T
+        return self._jac(x, constants, "unscaled")
 
     def jvp_scaled(self, v, x, constants=None):
         """Compute Jacobian-vector product of self.compute_scaled.
@@ -1220,23 +1069,16 @@ def _jvp(self, v, x, constants=None, op="scaled_error"):
         # and the Jacobian we want is dG/dc - dG/dx * (dF/dx)^-1 * dF/dc
 
         # Note: This Jacobian can be obtained using JVPs in proper tangent directions.
-        # First we will compute the tangent direction (see _get_tangent for details),
-        # then we will compute the Jacobian.
+        # First we will compute the tangent direction (see _proximal_get_tangents
+        # for details), then we will compute the Jacobian.
         v = v[0] if isinstance(v, (tuple, list)) else v
         constants = setdefault(constants, [None, None])
         xg, xf = self._update_equilibrium(x, store=True)
-
-        # we don't need to divide this part into blocked and batched because
-        # self._constraint._deriv_mode will handle it
-        jvpfun = lambda u: self._get_tangent(u, xf, constants, op=op)
-        tangents = batched_vectorize(
-            jvpfun,
-            signature="(n)->(k)",
-            chunk_size=self._constraint._jac_chunk_size,
-        )(v)
-
+        tangents = self._state.get_tangents(
+            xf, self._eq_idx, self._dimc_per_thing, op, v, constants[1]
+        )
         if self._objective._deriv_mode == "batched":
-            # objective's method already know about its jac_chunk_size
+            # objective's method already knows about its jac_chunk_size
             return getattr(self._objective, "jvp_" + op)(tangents, xg, constants[0])
         else:
             return _proximal_jvp_blocked_pure(
@@ -1246,52 +1088,75 @@ def _jvp(self, v, x, constants=None, op="scaled_error"):
                 op,
             )
 
-    def _get_tangent(self, v, xf, constants, op):
-        # Note: This function is vectorized over v. So, v is expected to be 1D array
-        # of size self.dim_x.
-
-        # v contains self._args DoFs from eq and other objects (like coils, surfaces
-        # etc), we want jvp_f to only get parts from equilibrium, not other things
-        vs = jnp.split(v, np.cumsum(self._dimc_per_thing))
-        # This is (dF/dx)^-1 * dF/dc  # noqa : E800
-        dfdc = _proximal_jvp_f_pure(
-            self._constraint,
-            xf,
-            constants[1],
-            vs[self._eq_idx],
-            self._eq_solve_objective._feasible_tangents,
-            self._dxdc,
-            op,
-        )
-        # broadcasting against multiple things
-        dfdcs = [jnp.zeros(dim) for dim in self._dimc_per_thing]
-        dfdcs[self._eq_idx] = dfdc
-        # note that dfdc.size != vs[self._eq_idx].size
-        # dfdc has the size of reduced state vector of the equilibrium
-        # but vs[self._eq_idx] has the size of self._args DoFs
-        dfdc = jnp.concatenate(dfdcs)
-
-        # We try to find dG/dc - dG/dx * (dF/dx)^-1 * dF/dc
-        # where G is the objective function. Since DESC stores x and c in the same
-        # vector, instead of multiple JVP calls, we will just find a tangent direction
-        # that will give us the same result.
-        # For making the explanation clear, assume J is the Jacobian of the objective
-        # function with respect to the full state vector (both x and c). Then,
-        # dG/dc = J @ (tangent vectors in c direction)
-        # dG/dx = J @ (tangent vectors in x direction)
-        # So, dG/dc - dG/dx * (dF/dx)^-1 * dF/dc can be written as
-        # J @ [(tangent vectors in c direction) - (tangent vectors in x direction)@dfdc]
-        # Note: We will never form full Jacobian J, we will just compute the above
-        # expression by JVPs.
-        dxdcv = jnp.concatenate(
-            [
-                *vs[: self._eq_idx],
-                self._dxdc @ vs[self._eq_idx],  # Rb_lmn, Zb_lmn to full eq state vector
-                *vs[self._eq_idx + 1 :],
-            ]
+    def _vjp(self, v, x, constants=None, op="scaled"):
+        constants = setdefault(constants, [None, None])
+        xg, xf = self._update_equilibrium(x, store=True)
+        tangents = self._state.get_tangents(
+            xf, self._eq_idx, self._dimc_per_thing, op, constants=constants[1]
         )
-        tangent = dxdcv - self._feasible_tangents @ dfdc
-        return tangent
+        v_vjp = getattr(self._objective, "vjp_" + op)(v, xg, constants[0])
+        return tangents @ v_vjp
+
+    def vjp_scaled(self, v, x, constants=None):
+        """Compute vector-Jacobian product of self.compute_scaled.
+
+        Parameters
+        ----------
+        v : ndarray or tuple of ndarray
+            Vectors to left-multiply the Jacobian by.
+        x : ndarray
+            Optimization variables.
+        constants : list
+            Constant parameters passed to sub-objectives. (Deprecated)
+
+        """
+        return self._vjp(v, x, constants, "scaled")
+
+    def vjp_scaled_error(self, v, x, constants=None):
+        """Compute vector-Jacobian product of self.compute_scaled_error.
+
+        Parameters
+        ----------
+        v : ndarray or tuple of ndarray
+            Vectors to left-multiply the Jacobian by.
+        x : ndarray
+            Optimization variables.
+        constants : list
+            Constant parameters passed to sub-objectives. (Deprecated)
+        """
+        return self._vjp(v, x, constants, "scaled_error")
+
+    def vjp_unscaled(self, v, x, constants=None):
+        """Compute vector-Jacobian product of self.compute_unscaled.
+
+        Parameters
+        ----------
+        v : ndarray or tuple of ndarray
+            Vectors to left-multiply the Jacobian by.
+        x : ndarray
+            Optimization variables.
+        constants : list
+            Constant parameters passed to sub-objectives. (Deprecated)
+        """
+        return self._vjp(v, x, constants, "unscaled")
+
+    def history(self, allx):
+        """Builds list of params for each proximal iterate."""
+        out = []
+        for x in allx:
+            xs = np.split(x, np.cumsum(self._dimc_per_thing)[:-1])
+            ceq = xs[self._eq_idx]
+
+            # read the full equilibrium parameters corresponding to a given ceq
+            xeq = f_where_x(ceq, self._state.allceq, self._state.allxeq)
+
+            params = self.unpack_state(x, False)
+
+            # unpack_state leaves the solve-output equilibrium params as 0s
+            eq_params = self._state.eq.unpack_params(xeq)
+            params[self._eq_idx] = eq_params
+            out.append(params)
+        return out
 
     @property
     def constants(self):
@@ -1304,23 +1169,289 @@ def constants(self):
             "of their objective compute methods. Instead declare all the "
             "constants in the build method and use as obj._constants.",
         )
-        return [self._objective.constants, self._constraint.constants]
+        return [self._objective.constants, self._state.constraint.constants]
 
     def __getattr__(self, name):
         """For other attributes we defer to the base objective."""
         return getattr(self._objective, name)
 
 
-# in ProximalProjection we have an explicit state that we keep track of (and add
-# to as we go) meaning if we jit anything with self static it doesn't update
-# correctly, while if we leave self unstatic then it recompiles every time because
-# the pytree structure of ProximalProjection is changing. To get around that we
-# define these helper functions that are stateless so we can safely jit them
+class ProximalState:
+    """State manager for objectives and constraints wrapped by ProximalProjection.
+
+    Provides a single source of equilibrium information, which can be shared
+    between different ProximalProjection instances. Stores the equilibrium
+    parameters, history, and caches tangents.
+
+    Parameters
+    ----------
+    eq: Equilibrium:
+        Equilibrium that is subject to the given constraint at each Proximal step.
+    constraint: ObjectiveFunction
+        Equilibrium constraint to enforce. Should be an ObjectiveFunction with one or
+        more of the following objectives: {ForceBalance, CurrentDensity,
+        RadialForceBalance, HelicalForceBalance}
+    perturb_options, solve_options : dict
+        Dictionary of arguments passed to Equilibrium.perturb and Equilibrium.solve
+        during the projection step.
+    cache_tangents : bool
+        Whether to compute and store the full Equilibrium tangents by default.
+        If True, this applies even to callers asking for fewer than dim(xeq)
+        directions. This is useful when the state is shared by multiple
+        ProximalProjection wrappers, which happens in augmented Lagrangian solvers.
+        Default is False.
+    """
+
+    def __init__(
+        self,
+        eq,
+        constraint,
+        perturb_options=None,
+        solve_options=None,
+        cache_tangents=False,
+    ):
+
+        assert isinstance(constraint, ObjectiveFunction), (
+            "constraint should be instance of ObjectiveFunction." ""
+        )
+        for con in constraint.objectives:
+            errorif(
+                not con._equilibrium,
+                ValueError,
+                "ProximalState cannot handle general " + f"nonlinear constraint {con}.",
+            )
+            # can't have bounds on constraint bc if constraint is satisfied then
+            # Fx == 0, and that messes with Gx @ Fx^-1 Fc etc.
+            errorif(
+                con.bounds is not None,
+                ValueError,
+                "ProximalState can only handle equality constraints, "
+                + f"got bounds for constraint {con}",
+            )
+
+        self.eq = eq
+        self.constraint = constraint
+
+        perturb_options = dict(setdefault(perturb_options, {}))
+        solve_options = dict(setdefault(solve_options, {}))
+        self._solve_during_proximal_build = solve_options.pop(
+            "solve_during_proximal_build", True
+        )  # If user does not want the solve during build, mainly for debug purposes
+        perturb_options.setdefault("verbose", 0)
+        perturb_options.setdefault("include_f", False)
+        solve_options.setdefault("verbose", 0)
+
+        self.perturb_options = perturb_options
+        self.solve_options = solve_options
+        self.allxeq = []
+        self.allceq = []
+        self.xeq_old = None
+        self.eq_is_current = True
+
+        # full equilibrium parameters at which tangents are computed
+        self._tangent_xf = None
+        # tangent cache
+        self._cache_tangents = cache_tangents
+        self._tangents = {}
+
+        self._built = False
+
+    def build(self, use_jit=None, verbose=1):
+        """Build the object.
+
+        Parameters
+        ----------
+        use_jit : bool, optional
+            Whether to just-in-time compile the objective and derivatives.
+        verbose : int, optional
+            Level of output.
+        """
+        if self._built:
+            return
+
+        self.eq_linear_constraints = get_fixed_boundary_constraints(eq=self.eq)
+        self.eq_linear_constraints = maybe_add_self_consistency(
+            self.eq, self.eq_linear_constraints
+        )
+
+        # we don't always build here because in ~all cases the user doesn't interact
+        # with this directly, so if the user wants to manually rebuild they should
+        # do it before this wrapper is created for them.
+        if not self.constraint.built:
+            self.constraint.build(use_jit=use_jit, verbose=verbose)
+
+        for constraint in self.eq_linear_constraints:
+            constraint.build(use_jit=use_jit, verbose=verbose)
+
+        # Here we create and build the LinearConstraintProjection
+        # for the equilibrium subproblem using the self._constraint as objective
+        # and our fixed-bdry constraints we just made. This will
+        # be passed as the objective for the eq subproblem, which saves
+        # some time as by building it here we can avoid re-computing the
+        # constraint matrix A and its SVD for the feasible direction method
+        self.eq_solve_objective = LinearConstraintProjection(
+            self.constraint,
+            ObjectiveFunction(self.eq_linear_constraints),
+            name="Eq Update LinearConstraintProjection",
+        )
+        self.eq_solve_objective.build(use_jit=use_jit, verbose=verbose)
+
+        errorif(
+            self.constraint.things != [self.eq],
+            ValueError,
+            "ProximalState can only handle constraints on the equilibrium.",
+        )
+
+        self._set_eq_state_vector()
+
+        if self._solve_during_proximal_build:
+            self.eq.solve(
+                objective=self.eq_solve_objective,
+                constraints=None,
+                **self.solve_options,
+            )
+
+        dims = [self.eq.dimensions[arg] for arg in self.args]
+        self.dim_ceq = int(np.sum(dims))
+        self.idx_ceq = np.cumsum(dims)[:-1]
+        self.allceq = [
+            jnp.concatenate(
+                [jnp.atleast_1d(self.eq.params_dict[arg]) for arg in self.args]
+            )
+        ]
+        self.allxeq = [self.eq.pack_params(self.eq.params_dict)]
+        self.xeq_old = self.eq.params_dict.copy()
+        self.eq_is_current = True
+        self._built = True
+
+    def _set_eq_state_vector(self):
+        """Removes equilibrium DOF which become dependent under Proximal."""
+        full_args = self.eq.optimizable_params.copy()
+        self.args = self.eq.optimizable_params.copy()
+
+        # the eq optimizable variables for proximal are the Rb, Zb and profile
+        # coefficients. Once these are chosen, we will solve the equilibrium to
+        # find the R_lmn, Z_lmn, L_lmn, Ra_n, Za_n. That is why we remove them
+        # from the list of optimizable variables. This is accompanied by not including
+        # self-consistency constraints (see get_combined_constraint_objectives in
+        # desc.optimize.optimizer) and also removing columns corresponding to these
+        # variables from the constraint matrix A in
+        # desc.objectives.utils.factorize_linear_constraints.
+        for arg in ["R_lmn", "Z_lmn", "L_lmn", "Ra_n", "Za_n"]:
+            self.args.remove(arg)
+
+        dxdc = []
+        xz = {arg: np.zeros(self.eq.dimensions[arg]) for arg in full_args}
+
+        for arg in self.args:
+            if arg not in ["Rb_lmn", "Zb_lmn"]:
+                x_idx = self.eq.x_idx[arg]
+                dxdc.append(np.eye(self.eq.dim_x)[:, x_idx])
+            if arg == "Rb_lmn":
+                c = get_instance(self.eq_linear_constraints, BoundaryRSelfConsistency)
+                A = c.jac_unscaled(xz)[0]["R_lmn"]
+                Ainv = np.linalg.pinv(A)
+                dxdRb = np.eye(self.eq.dim_x)[:, self.eq.x_idx["R_lmn"]] @ Ainv
+                dxdc.append(dxdRb)
+            if arg == "Zb_lmn":
+                c = get_instance(self.eq_linear_constraints, BoundaryZSelfConsistency)
+                A = c.jac_unscaled(xz)[0]["Z_lmn"]
+                Ainv = np.linalg.pinv(A)
+                dxdZb = np.eye(self.eq.dim_x)[:, self.eq.x_idx["Z_lmn"]] @ Ainv
+                dxdc.append(dxdZb)
+        # dxdc is a matrix that when multiplied by the optimization variables (only
+        # Rb_lmn, Zb_lmn) gives the full state vector of the equilibrium (Rb_lmn and
+        # Zb_lmn part will be 0, but they will be represented by the equivalent
+        # R_lmn and Z_lmn). For example, let's say the eq optimization variables are
+        # ceq = [Rb_lmn, Zb_lmn, p_l, i_l].T                      # noqa : E800
+        # Then, we will use dxdc for the following:
+        # xeq = dxdc @ ceq                                        # noqa : E800
+        # And xeq will be,
+        # xeq = [                                                 # noqa : E800
+        #     R_lmn, Z_lmn, jnp.zeros_like(L_lmn)                 # noqa : E800
+        #     jnp.zeros_like(Rb_lmn), jnp.zeros_like(Zb_lmn),     # noqa : E800
+        #     p_l, i_l,                                           # noqa : E800
+        # ]                                                       # noqa : E800
+        self.dxdc = jnp.hstack(dxdc)
+
+    def get_tangents(self, xf, eq_idx, dimc_per_thing, op, v=None, constants=None):
+        """Computes tangent directions for the ProximalProjection wrapper.
+
+        Checks if (xf, op) has been seen in the current iteration; if
+        so, returns the tangent vector/matrix. Otherwise, calls a given
+        function to compute tangents.
+
+        Parameters
+        ----------
+        xf : ndarray
+            Equilibrium state vector to compute the tangents at.
+        eq_idx: int
+            index of the equilibrium in the full set of things. Comes
+            from the ProximalProjection wrapper.
+        dimc_per_thing: list[int]
+            Number of optimizable params per thing. Comes from the
+            ProximalProjection wrapper.
+        op : str
+            One of ``scaled``, ``scaled_error``, or ``unscaled``.
+        v : ndarray, optional
+            Directions in the optimization variables. If None, the
+            identity directions are used and the result is cached.
+        constants : list
+            Constant parameters passed to the constraint.
 
+        Returns
+        -------
+        tangents : ndarray
+            Tangent directions in the full state vector of all the things.
 
-def jit_if_possible(func):
+        """
+        key = "scaled" if op in ["scaled", "scaled_error"] else "unscaled"
+        xf = jnp.asarray(xf)
+        if (self._tangent_xf is None) or (not np.array_equal(self._tangent_xf, xf)):
+            self._tangents = {}
+            self._tangent_xf = xf
+
+        v = jnp.eye(sum(dimc_per_thing)) if v is None else jnp.asarray(v)
+        vs = jnp.split(v, np.cumsum(dimc_per_thing)[:-1], axis=-1)
+
+        # If caller is already asking for at least as many directions
+        # as dimc of the equilibrium, then might as well compute and
+        # store tangents.
+        full_tangents = (
+            self._cache_tangents or vs[eq_idx].shape[0] >= dimc_per_thing[eq_idx]
+        )
+        v_eq = jnp.eye(dimc_per_thing[eq_idx]) if full_tangents else vs[eq_idx]
+        if key in self._tangents:
+            eq_tangents = vs[eq_idx] @ self._tangents[key]
+        else:
+            eq_tangents = _proximal_get_tangents(
+                self.constraint,
+                xf,
+                v_eq,
+                constants,
+                self.eq_solve_objective._feasible_tangents,
+                self.dxdc,
+                op,
+            )
+            if full_tangents:
+                self._tangents[key] = eq_tangents
+                eq_tangents = vs[eq_idx] @ eq_tangents
+
+        return jnp.concatenate([*vs[:eq_idx], eq_tangents, *vs[eq_idx + 1 :]], axis=-1)
+
+
+# ProximalState holds explicit state that we keep track of (and add to as we go),
+# meaning if we jit anything with it static it doesn't update correctly, while if we
+# leave it unstatic then it recompiles every time because the pytree structure is
+# changing. To get around that we define these helper functions that are stateless
+# so we can safely jit them.
+
+
+def jit_if_possible(func=None, *, static_argnames=("op",)):
     """Jit a function if use_jit."""
-    jitted_func = functools.partial(jit, static_argnames=["op"])(func)
+    if func is None:
+        return functools.partial(jit_if_possible, static_argnames=static_argnames)
+    jitted_func = functools.partial(jit, static_argnames=list(static_argnames))(func)
 
     @functools.wraps(func)
     def wrapper(*args, **kwargs):
@@ -1335,26 +1466,33 @@ def wrapper(*args, **kwargs):
 
 
 @jit_if_possible
-def _proximal_jvp_f_pure(constraint, xf, constants, dc, eq_feasible_tangents, dxdc, op):
-    # Note: This function is called by _get_tangent which is vectorized over v
-    # (v is called dc in this function). So, dc is expected to be 1D array
-    # of same size as full equilibrium state vector. This function returns a 1D array.
-
-    # here we are forming (dF/dx)^-1 @ dF/dc
-    # where Fxh is dF/dx and Fc is dF/dc
-    Fxh = getattr(constraint, "jvp_" + op)(eq_feasible_tangents.T, xf, constants).T
+def _proximal_eq_tangents(
+    constraint, xf, constants, eq_feasible_tangents, dxdcv, op="scaled_error"
+):
+    # Note: dxdcv holds the directions in c, mapped to the full eq state vector, as
+    # rows. It is either dxdc.T or v @ dxdc.T, the return has the same shape.
+
+    # here Fxh is dF/dx in the reduced (feasible) eq coordinates and Fc is dF/dc. A
+    # single batched JVP gives both, so the SVD below is computed once by
+    # construction, instead of relying on the compiler to hoist it out of a loop.
     # Our compute functions never include variables like Rb_lmn, Zb_lmn etc. So,
     # taking the JVP in just dc direction will give 0. To prevent this, we use dxdc
     # which is the dx/dc matrix and convert the Rb_lmn to R_lmn entries etc.
     # For example, if we want the derivative wrt Rb_023, we should take the derivative
     # wrt all R_lmn coefficients that contribute to Rb_023. See BoundaryRSelfConsistency
     # for the relation between Rb_lmn and R_lmn.
-    Fc = getattr(constraint, "jvp_" + op)(dxdc @ dc, xf, constants)
+    dim_x_reduced = eq_feasible_tangents.shape[-1]
+    tangents = jnp.concatenate([eq_feasible_tangents.T, dxdcv], axis=0)
+    J = getattr(constraint, "jvp_" + op)(tangents, xf, constants)
+    Fxh, Fc = J[:dim_x_reduced].T, J[dim_x_reduced:].T
     cutoff = jnp.finfo(Fxh.dtype).eps * max(Fxh.shape)
     uf, sf, vtf = jnp.linalg.svd(Fxh, full_matrices=False)
     sf += sf[-1]  # add a tiny bit of regularization
     sfi = jnp.where(sf < cutoff * sf[0], 0, 1 / sf)
-    return vtf.T @ (sfi * (uf.T @ Fc))
+    # this is (dF/dx)⁻¹ @ dF/dc for all the directions at once  # noqa : E800
+    dfdc = vtf.T @ (sfi[:, None] * (uf.T @ Fc))
+    # feasible_tangents maps the reduced eq state vector back to the full one
+    return dxdcv - (eq_feasible_tangents @ dfdc).T
 
 
 @jit_if_possible
@@ -1389,3 +1527,42 @@ def _proximal_jvp_blocked_pure(objective, vgs, xgs, op):
             outi = getattr(obj, "jvp_" + op)([_vi for _vi in vi], xi).T
             out.append(outi)
     return jnp.concatenate(out).T
+
+
+@jit_if_possible(static_argnames=("op",))
+def _proximal_get_tangents(
+    constraint,
+    xf,
+    veq,
+    constants,
+    eq_feasible_tangents,
+    dxdc,
+    op="scaled_error",
+):
+    # We try to find dG/dc - dG/dx * (dF/dx)⁻¹ * dF/dc
+    # where G is the objective function. Since DESC stores x and c in the same
+    # vector, instead of multiple JVP calls, we will just find a tangent direction
+    # that will give us the same result.
+    # For making the explanation clear, assume J is the Jacobian of the objective
+    # function with respect to the full state vector (both x and c). Then,
+    # dG/dc = J @ (tangent vectors in c direction)
+    # dG/dx = J @ (tangent vectors in x direction)
+    # So, dG/dc - dG/dx * (dF/dx)⁻¹ * dF/dc can be written as
+    # J @ [(tangent vectors in c direction) - (tangent vectors in x direction)@dfdc]
+    # Note: We will never form full Jacobian J, we will just compute the above
+    # expression by JVPs.
+
+    # veq contains prox._args DoFs from eq. This is the only block which changes
+    # when the equilibrium is re-solved.
+    if veq.ndim == 2 and veq.shape[0] > dxdc.shape[1]:
+        eq_tangents = veq @ _proximal_eq_tangents(
+            constraint, xf, constants, eq_feasible_tangents, dxdc.T, op
+        )
+    else:
+        dxdcv = veq @ dxdc.T
+        # atleast_2d and reshape are to also handle a single (1D) direction
+        eq_tangents = _proximal_eq_tangents(
+            constraint, xf, constants, eq_feasible_tangents, jnp.atleast_2d(dxdcv), op
+        )
+        eq_tangents = eq_tangents.reshape(dxdcv.shape)
+    return eq_tangents
diff --git a/desc/optimize/aug_lagrangian_ls.py b/desc/optimize/aug_lagrangian_ls.py
index 740b8584cd..8fe678f3f2 100644
--- a/desc/optimize/aug_lagrangian_ls.py
+++ b/desc/optimize/aug_lagrangian_ls.py
@@ -25,7 +25,7 @@
     inequality_to_bounds,
     print_header_nonlinear,
     print_iteration_nonlinear,
-    scale_columns,
+    scale_matrix,
     solve_triangular_regularized,
 )
 
@@ -288,7 +288,7 @@ def lagjac(z, y, mu, *args):
     # we don't need unscaled J anymore, so we overwrite it with J_h = J * d to avoid
     # carrying so many J-sized matrices in memory, which can be large. The buffer is
     # donated so the scaling doesn't allocate a second copy of J.
-    J_h = scale_columns(J, d)
+    J_h = scale_matrix(J, d)
     del J
     g_norm = jnp.linalg.norm(
         (g * v * scale if scaled_termination else g * v), ord=jnp.inf
@@ -546,6 +546,7 @@ def lagjac(z, y, mu, *args):
 
             # updating augmented lagrangian params
             if g_norm < gtolk:
+                mu_old = mu
                 y = jnp.where(jnp.abs(c) < ctolk, y - mu * c, y)
                 mu = jnp.where(jnp.abs(c) >= ctolk, tau * mu, mu)
                 if constr_violation < ctolk:
@@ -554,12 +555,14 @@ def lagjac(z, y, mu, *args):
                 else:
                     ctolk = max(eta / (jnp.mean(mu) ** alpha_eta), ctol)
                     gtolk = max(omega / (jnp.mean(mu) ** alpha_omega), gtol)
-                # if we update lagrangian params, need to recompute L and J
+                # if we update lagrangian params, need to recompute L and J.
                 L = lagfun(f, c, y, mu)
                 Lcost = 0.5 * jnp.dot(L, L)
-                del J
-                J = lagjac(z, y, mu, *args)
-                njev += 1
+                row_scale = jnp.concatenate(
+                    [jnp.ones_like(f), jnp.sqrt(mu) / jnp.sqrt(mu_old)]
+                )
+                J = scale_matrix(J, row_scale[:, None])
+
                 g = jnp.dot(L, J)
 
                 if jac_scale:
@@ -577,7 +580,7 @@ def lagjac(z, y, mu, *args):
             d = v**0.5 * scale
             diag_h = g * dv * scale
             g_h = g * d
-            J_h = scale_columns(J, d)
+            J_h = scale_matrix(J, d)
             del J
 
             if g_norm < gtol and constr_violation < ctol:
@@ -611,7 +614,7 @@ def lagjac(z, y, mu, *args):
     active_mask = find_active_constraints(z, zbounds[0], zbounds[1], rtol=xtol)
     # after overwriting J_h with J*d, we have to revert back and store the
     # unscaled version
-    J_h = scale_columns(J_h, 1 / d)
+    J_h = scale_matrix(J_h, 1 / d)
     result = OptimizeResult(
         x=x,
         s=s,
diff --git a/desc/optimize/fmin_scalar.py b/desc/optimize/fmin_scalar.py
index 4c5a85f2b6..1e6004d18d 100644
--- a/desc/optimize/fmin_scalar.py
+++ b/desc/optimize/fmin_scalar.py
@@ -24,7 +24,7 @@
     compute_hess_scale,
     print_header_nonlinear,
     print_iteration_nonlinear,
-    scale_columns,
+    scale_matrix,
 )
 
 
@@ -243,9 +243,9 @@ def fmintr(  # noqa: C901
     # we don't need unscaled H anymore, so we overwrite it with H_h = d[:, None] * H * d
     # to avoid carrying so many H-sized matrices in memory, which can be large. The
     # buffer is donated so the scaling doesn't allocate a second copy of H.
-    H_h = scale_columns(H, d)
+    H_h = scale_matrix(H, d)
     del H
-    H_h = scale_columns(H_h, d[:, None])
+    H_h = scale_matrix(H_h, d[:, None])
 
     g_norm = jnp.linalg.norm(
         (g * v * scale if scaled_termination else g * v), ord=jnp.inf
@@ -440,9 +440,9 @@ def fmintr(  # noqa: C901
 
             g_h = g * d
 
-            H_h = scale_columns(H, d)
+            H_h = scale_matrix(H, d)
             del H
-            H_h = scale_columns(H_h, d[:, None])
+            H_h = scale_matrix(H_h, d[:, None])
 
             x_norm = jnp.linalg.norm(
                 ((x * scale_inv) if scaled_termination else x), ord=2
@@ -474,8 +474,8 @@ def fmintr(  # noqa: C901
     active_mask = find_active_constraints(x, lb, ub, rtol=xtol)
     # after overwriting H_h with the scaled version, we have to revert back and
     # store the unscaled one
-    H_h = scale_columns(H_h, 1 / d)
-    H_h = scale_columns(H_h, 1 / d[:, None])
+    H_h = scale_matrix(H_h, 1 / d)
+    H_h = scale_matrix(H_h, 1 / d[:, None])
     result = OptimizeResult(
         x=x,
         success=success,
diff --git a/desc/optimize/least_squares.py b/desc/optimize/least_squares.py
index a73d14079c..9043f49e84 100644
--- a/desc/optimize/least_squares.py
+++ b/desc/optimize/least_squares.py
@@ -24,7 +24,7 @@
     compute_jac_scale,
     print_header_nonlinear,
     print_iteration_nonlinear,
-    scale_columns,
+    scale_matrix,
     solve_triangular_regularized,
 )
 
@@ -209,7 +209,7 @@ def lsqtr(  # noqa: C901
     # we don't need unscaled J anymore, so we overwrite it with J_h = J * d to avoid
     # carrying so many J-sized matrices in memory, which can be large. The buffer is
     # donated so the scaling doesn't allocate a second copy of J.
-    J_h = scale_columns(J, d)
+    J_h = scale_matrix(J, d)
     del J
     g_norm = jnp.linalg.norm(
         (g * v * scale if scaled_termination else g * v), ord=jnp.inf
@@ -427,7 +427,7 @@ def lsqtr(  # noqa: C901
             diag_h = g * dv * scale
 
             g_h = g * d
-            J_h = scale_columns(J, d)
+            J_h = scale_matrix(J, d)
             del J
             x_norm = jnp.linalg.norm(
                 ((x * scale_inv) if scaled_termination else x), ord=2
@@ -458,7 +458,7 @@ def lsqtr(  # noqa: C901
     active_mask = find_active_constraints(x, lb, ub, rtol=xtol)
     # after overwriting J_h with J*d, we have to revert back and store the
     # unscaled version
-    J_h = scale_columns(J_h, 1 / d)
+    J_h = scale_matrix(J_h, 1 / d)
     result = OptimizeResult(
         x=x,
         success=success,
diff --git a/desc/optimize/optimizer.py b/desc/optimize/optimizer.py
index 2d4b114a55..d7e7caa4e8 100644
--- a/desc/optimize/optimizer.py
+++ b/desc/optimize/optimizer.py
@@ -28,7 +28,11 @@
     warnif,
 )
 
-from ._constraint_wrappers import LinearConstraintProjection, ProximalProjection
+from ._constraint_wrappers import (
+    LinearConstraintProjection,
+    ProximalProjection,
+    ProximalState,
+)
 
 
 class Optimizer(IOAble):
@@ -335,7 +339,7 @@ def optimize(  # noqa: C901
             # reset eq params to initial
             if eq is not None:
                 eq.params_dict = eq_params_init
-            result["history"] = objective.history
+            result["history"] = objective.history(result["allx"])
             objective = objective._objective
         else:
             result["history"] = [
@@ -460,10 +464,11 @@ def _project_x_scale(x_scale, objective):
         # Project equilibrium part: remove excluded parameters
         excluded_params = ["R_lmn", "Z_lmn", "L_lmn", "Ra_n", "Za_n"]
         included_idx = []
-        for arg in prox_obj._eq.optimizable_params:
+        for arg in prox_obj._state.eq.optimizable_params:
             if arg not in excluded_params:
-                included_idx.extend(prox_obj._eq.x_idx[arg])
-        x_scale[prox_obj._eq_idx] = x_scale[prox_obj._eq_idx][jnp.array(included_idx)]
+                included_idx.extend(prox_obj._state.eq.x_idx[arg])
+        eq_idx = prox_obj._eq_idx
+        x_scale[eq_idx] = x_scale[eq_idx][jnp.array(included_idx)]
         x_scale = jnp.concatenate(x_scale)
 
     if isinstance(objective, LinearConstraintProjection):
@@ -530,6 +535,8 @@ def _combine_constraints(constraints):
         Otherwise returns None.
 
     """
+    if len(constraints) == 1 and isinstance(constraints[0], ObjectiveFunction):
+        return constraints[0]
     if len(constraints):
         objective = ObjectiveFunction(constraints)
     else:
@@ -576,6 +583,32 @@ def _parse_constraints(constraints):
     return linear_constraints, nonlinear_constraints
 
 
+def _parse_nonlinear_constraints(eq, nonlinear_constraints):
+    """Split nonlinear constraints for Proximal.
+
+    Parameters
+    ----------
+    eq: Equilibrium
+    nonlinear constraints : tuple of Objective
+        Nonlinear constraints to parse.
+
+    Returns
+    -------
+    eq_constraints : tuple of Objective
+        Constraints which can be handled directly by Proximal,
+        i.e. ones involving the equilibrium and without bounds.
+    other_constraints : tuple of Objective
+        Remaining nonlinear constraints. Can be empty.
+    """
+    eq_constraints, other_constraints = [], []
+    for con in nonlinear_constraints:
+        if con._equilibrium and con.bounds is None and con.things == [eq]:
+            eq_constraints.append(con)
+        else:
+            other_constraints.append(con)
+    return eq_constraints, other_constraints
+
+
 def _maybe_wrap_nonlinear_constraints(
     eq, objective, nonlinear_constraints, method, options
 ):
@@ -599,16 +632,38 @@ def _maybe_wrap_nonlinear_constraints(
                 """))
         wrapper = "proximal"
     if wrapper is not None and wrapper.lower() in ["prox", "proximal"]:
+        eq_constraints, other_constraints = _parse_nonlinear_constraints(
+            eq, nonlinear_constraints
+        )
+        errorif(
+            not len(eq_constraints),
+            ValueError,
+            f"Method {wrapper}-{method} requires at least one equilibrium constraint "
+            + "(e.g. ForceBalance), but none were given.",
+        )
+        if len(other_constraints) and not optimizers[method]["equality_constraints"]:
+            raise ValueError(
+                f"Method {wrapper}-{method} cannot accept nonlinear constraints. "
+                + "Consider an Augmented Lagrangian solver such as proximal-lsq-auglag."
+            )
+
         perturb_options = options.pop("perturb_options", {})
         solve_options = options.pop("solve_options", {})
-        objective = ProximalProjection(
-            objective,
-            constraint=_combine_constraints(nonlinear_constraints),
+        state = ProximalState(
+            eq,
+            _combine_constraints(eq_constraints),
             perturb_options=perturb_options,
             solve_options=solve_options,
-            eq=eq,
+            cache_tangents=bool(other_constraints),
         )
-        nonlinear_constraints = ()
+
+        objective = ProximalProjection(objective, state=state)
+        nonlinear_constraints = (
+            (ProximalProjection(_combine_constraints(other_constraints), state=state),)
+            if len(other_constraints)
+            else ()
+        )
+
     return objective, nonlinear_constraints
 
 
@@ -718,7 +773,7 @@ def get_combined_constraint_objectives(  # noqa: C901
         objective.build(verbose=verbose)
         if nonlinear_constraint is not None:
             nonlinear_constraint = LinearConstraintProjection(
-                nonlinear_constraint, linear_constraint
+                nonlinear_constraint, linear_constraint, **linear_constraint_options
             )
             nonlinear_constraint.build(verbose=verbose)
 
diff --git a/desc/optimize/utils.py b/desc/optimize/utils.py
index d3b938eae0..aa05c47ce3 100644
--- a/desc/optimize/utils.py
+++ b/desc/optimize/utils.py
@@ -504,11 +504,11 @@ def compute_jac_scale(A, prev_scale_inv=None):
 
 
 @functools.partial(jit, donate_argnums=0)
-def scale_columns(A, d):
+def scale_matrix(A, d):
     """Compute `A * d` reusing `A`'s buffer instead of allocating a second copy.
 
     `A` is invalid after this call, so callers must rebind, ie
-    `A = scale_columns(A, d)`.
+    `A = scale_matrix(A, d)`.
     """
     return A * d
 
diff --git a/tests/test_optimizer.py b/tests/test_optimizer.py
index 1a933f6485..9086adbef1 100644
--- a/tests/test_optimizer.py
+++ b/tests/test_optimizer.py
@@ -35,6 +35,7 @@
     AspectRatio,
     BoundaryRSelfConsistency,
     BoundaryZSelfConsistency,
+    CoilCurvature,
     CoilLength,
     Energy,
     FixBoundaryR,
@@ -1188,6 +1189,158 @@ def test_constrained_AL_scalar():
     np.testing.assert_array_less(-Dwell, ctol)
 
 
+@pytest.mark.slow
+@pytest.mark.regression
+@pytest.mark.optimize
+def test_proximal_constrained_AL_lsq():
+    """Test proximal-lsq-auglag with a non-equilibrium nonlinear constraint."""
+    eq = desc.examples.get("SOLOVEV")
+    solve_options = {"ftol": 1e-8, "xtol": 1e-8, "gtol": 1e-8, "maxiter": 50}
+    volume_target = 0.95 * float(eq.compute("V")["V"])
+
+    coil = FourierPlanarCoil(
+        r_n=[1.0], center=[4.0, 0, 0], normal=[0, 1, 0], current=1e6
+    )
+    length_target = 1.1 * float(coil.compute("length")["length"])
+
+    R_modes = np.vstack(
+        (
+            [0, 0, 0],
+            eq.surface.R_basis.modes[
+                np.max(np.abs(eq.surface.R_basis.modes), 1) > 1, :
+            ],
+        )
+    )
+    Z_modes = eq.surface.Z_basis.modes[
+        np.max(np.abs(eq.surface.Z_basis.modes), 1) > 1, :
+    ]
+
+    objective = ObjectiveFunction(
+        (
+            CoilCurvature(coil, target=0.5, weight=1e-2),
+            Volume(eq=eq, target=volume_target),
+        )
+    )
+    constraints = (
+        ForceBalance(eq=eq),  # absorbed by proximal
+        CoilLength(coil, target=length_target),  # auglag
+        FixBoundaryR(eq=eq, modes=R_modes),
+        FixBoundaryZ(eq=eq, modes=Z_modes),
+        FixPressure(eq=eq),
+        FixIota(eq=eq),
+        FixPsi(eq=eq),
+        FixCoilCurrent(coil),
+    )
+
+    objective.build(verbose=0)
+    prox = ProximalProjection(objective, ObjectiveFunction(ForceBalance(eq=eq)), eq)
+    prox.build(verbose=0)
+    state = prox._state
+
+    for arg in ["R_lmn", "Z_lmn", "L_lmn", "Ra_n", "Za_n"]:
+        assert arg not in state.args
+    dim_eq = sum(eq.dimensions[arg] for arg in state.args)
+    assert dim_eq < eq.dim_x
+    assert prox.dim_x == dim_eq + coil.dim_x
+    assert prox._dimc_per_thing[prox._eq_idx] == dim_eq
+    assert prox._dimx_per_thing[prox._eq_idx] == eq.dim_x
+    assert prox._dimc_per_thing[0] == prox._dimx_per_thing[0] == coil.dim_x
+
+    (eq_opt, coil_opt), _ = Optimizer("proximal-lsq-auglag").optimize(
+        (eq, coil),
+        objective=objective,
+        constraints=constraints,
+        maxiter=30,
+        verbose=0,
+        copy=True,
+        options={"solve_options": solve_options},
+    )
+
+    np.testing.assert_allclose(
+        float(coil_opt.compute("length")["length"]), length_target, rtol=1e-8
+    )
+    np.testing.assert_allclose(
+        float(eq_opt.compute("V")["V"]), volume_target, rtol=1e-6
+    )
+
+    force = ObjectiveFunction(ForceBalance(eq=eq_opt))
+    force.build(verbose=0)
+    force_final = np.linalg.norm(force.compute_scaled_error(force.x(eq_opt)))
+    assert force_final < 1e-7
+
+
+@pytest.mark.slow
+@pytest.mark.regression
+@pytest.mark.optimize
+def test_proximal_constrained_AL_scalar():
+    """Test proximal-fmin-auglag with a non-equilibrium nonlinear constraint."""
+    eq = desc.examples.get("SOLOVEV")
+    solve_options = {"ftol": 1e-8, "xtol": 1e-8, "gtol": 1e-8, "maxiter": 50}
+    volume_target = 0.95 * float(eq.compute("V")["V"])
+
+    coil = FourierPlanarCoil(
+        r_n=[1.0], center=[4.0, 0, 0], normal=[0, 1, 0], current=1e6
+    )
+    length_target = 1.1 * float(coil.compute("length")["length"])
+
+    R_modes = np.vstack(
+        (
+            [0, 0, 0],
+            eq.surface.R_basis.modes[
+                np.max(np.abs(eq.surface.R_basis.modes), 1) > 1, :
+            ],
+        )
+    )
+    Z_modes = eq.surface.Z_basis.modes[
+        np.max(np.abs(eq.surface.Z_basis.modes), 1) > 1, :
+    ]
+
+    # dummy objective, targets are expressed as constraints for the aug lagrangian
+    objective = ObjectiveFunction(GenericObjective("0", thing=eq))
+    constraints = (
+        ForceBalance(eq=eq),  # absorbed by proximal
+        CoilLength(coil, target=length_target),  # auglag
+        Volume(eq=eq, target=volume_target),  # auglag
+        FixBoundaryR(eq=eq, modes=R_modes),
+        FixBoundaryZ(eq=eq, modes=Z_modes),
+        FixPressure(eq=eq),
+        FixIota(eq=eq),
+        FixPsi(eq=eq),
+        FixCoilCurrent(coil),
+    )
+
+    ctol = 1e-4
+    (eq_opt, coil_opt), _ = Optimizer("proximal-fmin-auglag").optimize(
+        (eq, coil),
+        objective=objective,
+        constraints=constraints,
+        maxiter=30,
+        ctol=ctol,
+        verbose=0,
+        copy=True,
+        options={
+            "solve_options": solve_options,
+            "initial_penalty_parameter": 1e3,
+            # high penalty parameter helps coil length converge quickly
+        },
+    )
+
+    np.testing.assert_allclose(
+        float(coil_opt.compute("length")["length"]),
+        length_target,
+        rtol=ctol,
+        atol=ctol,
+    )
+    np.testing.assert_allclose(
+        float(eq_opt.compute("V")["V"]), volume_target, rtol=ctol, atol=ctol
+    )
+
+    force = ObjectiveFunction(ForceBalance(eq=eq_opt))
+    force.build(verbose=0)
+    force_final = np.linalg.norm(force.compute_scaled_error(force.x(eq_opt)))
+    assert force_final < 1e-7
+
+
 @pytest.mark.unit
 @pytest.mark.optimize
 def test_optimize_multiple_things_different_order():
@@ -1360,6 +1513,10 @@ def test_proximal_jacobian():
     prox2.build()
     prox3.build()
 
+    unfixed_idx = prox1._state.eq_solve_objective._unfixed_idx
+    Z = prox1._state.eq_solve_objective._Z
+    dxdc = prox1._state.dxdc
+
     x = prox1.x(eq)
     v = np.random.default_rng(1138).random(x.shape)
 
@@ -1369,10 +1526,10 @@ def test_proximal_jacobian():
     # for scaled jacobian
     Fx = con1.jac_scaled(xf)
     Gx = obj1.jac_scaled(xg)
-    Fxh = Fx[:, prox1._eq_unfixed_idx] @ prox1._eq_Z
-    Gxh = Gx[:, prox1._eq_unfixed_idx] @ prox1._eq_Z
-    Fc = Fx @ prox1._dxdc
-    Gc = Gx @ prox1._dxdc
+    Fxh = Fx[:, unfixed_idx] @ Z
+    Gxh = Gx[:, unfixed_idx] @ Z
+    Fc = Fx @ dxdc
+    Gc = Gx @ dxdc
     cutoff = np.finfo(Fxh.dtype).eps * np.max(Fxh.shape)
     uf, sf, vtf = jnp.linalg.svd(Fxh, full_matrices=False)
     sf += sf[-1]  # add a tiny bit of regularization
@@ -1382,10 +1539,10 @@ def test_proximal_jacobian():
     # for unscaled jacobian
     Fx = con1.jac_unscaled(xf)
     Gx = obj1.jac_unscaled(xg)
-    Fxh = Fx[:, prox1._eq_unfixed_idx] @ prox1._eq_Z
-    Gxh = Gx[:, prox1._eq_unfixed_idx] @ prox1._eq_Z
-    Fc = Fx @ prox1._dxdc
-    Gc = Gx @ prox1._dxdc
+    Fxh = Fx[:, unfixed_idx] @ Z
+    Gxh = Gx[:, unfixed_idx] @ Z
+    Fc = Fx @ dxdc
+    Gc = Gx @ dxdc
     cutoff = np.finfo(Fxh.dtype).eps * np.max(Fxh.shape)
     uf, sf, vtf = jnp.linalg.svd(Fxh, full_matrices=False)
     sf += sf[-1]  # add a tiny bit of regularization