Paper deep dive
Directional Influence Function: Estimating Training Data Influence in Constrained Learning
Xin Wang, R. Tyrrell Rockafellar, Xuegang, Ban
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:As constrained learning becomes increasingly common, models are trained under explicit feasibility requirements to enforce fairness, safety, robustness, regulariza- tion, and physics or logic constraints. Understanding how training samples in- fluence the model solution (e.g., learned parameters) is crucial for interpretability and robustness. The classical influence function (IF) estimates sample contribu- tions via local sensitivity analysis, measuring how the solution changes when a specific training sample is perturbed or removed. However, IF becomes unreli- able in constrained settings: data perturbations can reshape both the objective and the feasible region, leading to estimates that violate feasibility. In response, we propose the Directional Influence Function (DIF), a novel estimator that explicitly incorporates these constraints into influence estimation. DIF formulates the opti- mality conditions of constrained learning as a variational inequality (VI) and ana- lyzes how perturbing training data affects this VI. We validate DIF on constrained linear regression and demonstrate that it recovers leave-one-out retraining results, whereas IF and penalty-based IF exhibit significant bias. We further apply DIF to fairness-constrained CNNs, where DIF accurately predicts test loss changes under data removal and aligns closely with actual retraining. Our results establish DIF as an efficient and reliable tool for data attribution in constrained learning.
Tags
Links
- Source: https://arxiv.org/abs/2607.23388v1
- Canonical: https://arxiv.org/abs/2607.23388v1
Trouble viewing inline? Open PDF directly →
Full Text
95,769 characters extracted from source content.
Expand or collapse full text
Directional Influence Function: Estimating Training Data Influence in Constrained Learning Xin Wang Department of Civil and Environmental Engineering University of Washington Seattle, WA, USA &R. Tyrrell Rockafellar Department of Mathematics University of Washington Seattle, WA, USA &Xuegang (Jeff) Ban Department of Civil and Environmental Engineering University of Washington Seattle, WA, USA banx@uw.edu Corresponding author. Abstract As constrained learning becomes increasingly common, models are trained under explicit feasibility requirements to enforce fairness, safety, robustness, regularization, and physics or logic constraints. Understanding how training samples influence the model solution (e.g., learned parameters) is crucial for interpretability and robustness. The classical influence function (IF) estimates sample contributions via local sensitivity analysis, measuring how the solution changes when a specific training sample is perturbed or removed. However, IF becomes unreliable in constrained settings: data perturbations can reshape both the objective and the feasible region, leading to estimates that violate feasibility. In response, we propose the Directional Influence Function (DIF), a novel estimator that explicitly incorporates these constraints into influence estimation. DIF formulates the optimality conditions of constrained learning as a variational inequality (VI) and analyzes how perturbing training data affects this VI. We validate DIF on constrained linear regression and demonstrate that it recovers leave-one-out retraining results, whereas IF and penalty-based IF exhibit significant bias. We further apply DIF to fairness-constrained CNNs, where DIF accurately predicts test loss changes under data removal and aligns closely with actual retraining. Our results establish DIF as an efficient and reliable tool for data attribution in constrained learning. 1 Introduction Understanding how individual training samples influence model predictions, also known as data attribution (Pruthi et al., 2020; Feldman and Zhang, 2020; Lin et al., 2024), is fundamental for interpreting model behavior, model correction, and data error debugging. Although retraining the model after removing a sample and comparing the resulting change in model solution (learned parameters) provides ground-truth data influence, this approach is computationally expensive for modern deep learning models. A widely adopted alternative is to approximate sample influence by studying how the solution responds to small data perturbations. As a representative method, the influence function (IF), originally developed in robust statistics (Hampel, 1974) and later adapted to machine learning (Koh and Liang, 2017; Feldman and Zhang, 2020; Zhang et al., 2024), instantiates this idea by differentiating the solution with respect to (w.r.t.) small data perturbations. While these influence estimators work for unconstrained learning problems, they often fail when applied to constrained learning. The latter refers to learning tasks constrained by domain knowledge, such as in physics-informed neural networks (PINNs), or those constrained by embedding requirements, including fairness, safety, and robustness (Hounie et al., 2023; Li et al., 2023; Garcıa and Fernández, 2015; Xia et al., 2025). The empirical formulation of constrained learning is (Chamon and Ribeiro, 2020): θ¯:= θ= argminθ∈ℝd1N0∑i=1N0ℓ0(zi(0),θ) _θ ^d 1N_0 _i=1^N_0 _0 (z_i^(0),θ ) (1) s.t. 1Nj∑i=1Njℓj(zi(j),θ)≤τj,j=1,…,m, 1N_j _i=1^N_j _j (z_i^(j),θ )≤ _j, j=1,…,m, where θ∈ℝdθ ^d denotes the model parameters to be optimized, and θ¯ θ is the optimal solution. The function ℓ0(zi(0),θ) _0 (z_i^(0),θ ) represents the primary loss evaluated on the training data point zi(0)z_i^(0). Each constraint j∈1,…,mj∈\1,…,m\ is represented by an auxiliary loss function ℓj(zi(j),θ) _j (z_i^(j),θ ) evaluated over a separate dataset zi(j)i=1Nj \z_i^(j) \_i=1^N_j. The constant τj _j specifies the upper bound allowed for the jjth constraint. The influence estimators aim to quantify the change in solution, denoted by Δθ θ, when specific training data points participating in the objective function or constraints (e.g., z1(1)z^(1)_1) are removed. However, current IF-based approaches fail in constrained learning for two primary reasons. First, constrained problems (1) require that the Δθ θ remain within the feasible region, while IF methods ignore the constraints, not to mention that the feasible region itself may also be altered by data perturbation. Applying existing IF cannot guarantee that the estimated solution change is feasible. Second, IF primarily relies on the gradient of the optimal solution θ¯ θ w.r.t. data perturbations. In constrained learning, the solution can be only directionally differentiable. The full derivatives do not necessarily exist, rendering IF estimators invalid. This issue arises because the solution updates must adhere to the feasible region when constraints are active, often requiring projections to maintain feasibility. Additionally, data removal can change the active constraint set, leading to sudden shifts in the optimal solution. See the problem (3) for a concrete example. To overcome the above limitations, this paper introduces the Directional Influence Function (DIF), an influence function designed for constrained learning, to estimate the impact of training points participating in either the loss function or the constraints on the model solution. The DIF estimates the change in model solution through a directional derivative approach, aiming to address the following question: How does the solution of constrained learning change when data points are removed from either the objective loss function or the constraints? Our main contributions are as follows: (1) We formalize data attribution for constrained learning by casting the optimality conditions as a variational inequality (VI) and performing local sensitivity analysis of this VI. (2) DIF quantifies the effect of data perturbations on model solutions via directional derivatives, thereby addressing the non-smoothness of solution changes induced by constraints. (3) We propose an approach that computes DIF by solving a quadratic program (QP) and prove that, when all constraints are inactive (i.e., their KKT multipliers are zero), DIF reduces to the classical IF. The remainder of this paper is organized as follows. Section 2 formally defines our problem and illustrates the failure of IF in a constrained learning setting through a linear regression example. Section 3 introduces the proposed DIF. Section 4 analyzes the impact of data perturbations on the VI-formulated optimality conditions, introduces a QP approach to compute DIF, and presents our main theoretical results. Section 5 evaluates DIF on a constrained linear regression task and a constrained CNN model. Throughout this paper, we assume that all functions involved are twice continuously differentiable (C2C^2). 2 Problem Formulation Let Z(0)=zi(0)i=1N0Z^(0)= \z_i^(0) \_i=1^N_0 be the dataset for the objective function and Z(j)=zi(j)i=1NjZ^(j)= \z_i^(j) \_i=1^N_j the dataset for the j-th constraint. We begin by examining how the solution θ¯ θ changes when a selected subset Zr⊂⋃j=0mZ(j)Z^r⊂ _j=0^mZ^(j) is removed. This removal is modeled via a perturbation formulation, where each point in ZrZ^r is assigned a small weight εk:ε0 _k: _0 for points contributing to the objective, and εj,j∈1,…,m _j,j∈\1,…,m\ for the points contributing to the j-th constraint. This formulation, referred to as Perturbed Constrained Learning (PCL), is defined as follows111To simplify our discussion, this constrained learning formulation does not explicitly represent equality constraints, as equalities can be transformed into inequalities.: minθ _θ 1N0∑i=1N0ℓ0(zi(0),θ)+ε0∑zi(0)∈Zrℓ0(zi(0),θ) 1N_0 _i=1^N_0 _0 (z_i^(0),θ )+ _0 _z_i^(0)∈ Z^r _0 (z_i^(0),θ ) (2) s.t. 1Nj∑i=1Njℓj(zi(j),θ)+εj∑zi(j)∈Zrℓj(zi(j),θ)≤τj,j=1,…,m. 1N_j _i=1^N_j _j (z_i^(j),θ )+ _j _z_i^(j)∈ Z^r _j (z^(j)_i,θ )≤ _j,j=1,…,m. Let ε=[ε0,…,εm] = [ _0,…, _m ] denote the perturbation vector. When ε= =0, problem (2) reduces to the problem (1). Shifting ε from ε¯=[0,0,…,0] =[0,0,…,0] to ε^= = [−1N0,−1N1,…,−1Nm] [- 1N_0,- 1N_1,…,- 1N_m ] corresponds to removing ZrZ^r from the constrained learning problem (1). We treat θ as a function of ε , denoted by θ(ε)θ( ), and define the solution change as Δθ=θ(ε^)−θ(ε¯) θ=θ( )-θ( ). By definition, θ(ε¯)=θ¯θ( )= θ is the solution to problem (1). 2.1 Failure of IF in Constrained Learning Toy example. We utilize a ℓ1 _1-constrained least-squares example to demonstrate the failure of the IF in constrained learning. Consider the following toy example, where θ=[θ1,θ2]θ= [ _1, _2 ] represents the model parameters. minθ _θ 13(θ⊤x1−y1)2+13(θ⊤x2−y2)2+13(θ⊤x3−y3)2 13(θ x_1-y_1)^2+ 13(θ x_2-y_2)^2+ 13(θ x_3-y_3)^2 (3) +ε⋅(θ⊤x2−y2)2 + ·(θ x_2-y_2)^2 s.t. ‖θ‖1≤1. \|θ\|_1≤ 1. The data is provided as follows: [x1x2x3]=[101001],[y1y2y3]=[101/2] [ array[]lx_1\\ x_2\\ x_3 array ]= [ array[]l1&0\\ 1&0\\ 0&1 array ], [ array[]ly_1\\ y_2\\ y_3 array ]= [ array[]l1\\ 0\\ 1/2 array ] The explicit solution of this problem is: θ(ε)=Proj∥⋅∥1≤1[(12+3ε,0.5)]θ( )=Proj_\|·\|_1≤ 1 [ ( 12+3 ,0.5 ) ] (4) Removing (x2,y2)(x_2,y_2) is equivalent to shift ε from 0 to −13- 13. The IF defines the impact caused by data removal using the derivative of θ(ε)θ( ): IF(z2)=dθ(ε)dε|ε=0=−Hθ¯−1∇θℓ0(x2,y2,θ¯),IF(z_2)= . dθ( )d |_ =0=-H_ θ^-1 _θ _0(x_2,y_2, θ), where z2=(x2,y2)z_2=(x_2,y_2), Hθ¯=13∑i=13∇θ2ℓ0(xi,yi,θ¯)H_ θ= 13 _i=1^3 _θ^2 _0(x_i,y_i, θ) is the Hessian, and ℓ0(x2,y2,θ¯)=(θ¯⊤x2−y2)2 _0(x_2,y_2, θ)=( θ x_2-y_2)^2. IF then estimates the change Δθ θ ignoring the constraint by: θ(−13)−θ(0)≈−13dθ(ε)dε|ε=0.θ(- 13)-θ(0)≈- 13 . dθ( )d |_ =0. Figure 1: Loss landscape of the toy example when ε=0 =0 and ε=−13 =- 13. The left plot illustrates the non-existence of the derivative of θ(ε)θ( ) at ε=0 =0 as dθ(ε)dε|ε=0+≠dθ(ε)dε|ε=0− . dθ( )d |_ =0^+≠ . dθ( )d |_ =0^-. The right plot shows that the IF estimation of Δθ θ fails to account for the feasible region, whereas our DIF accurately predicts the solution change direction and closely matches the ground truth. However, as shown in the left panel of figure 1, the dθ(ε)dε|ε=0 . dθ( )d |_ =0 is not well defined because the left derivative and right derivative are inconsistent. When ε>0 >0, the point (12+3ε,0.5) ( 12+3 ,0.5 ) satisfies ∥⋅∥1<1\|·\|_1<1 and thus lies in the interior of the feasible region. In contrast, when ε<0 <0, its ℓ1 _1-norm exceeds 11, so the solution is projected onto the ℓ1 _1-boundary, which results in two different one-sided derivatives at ε=0 =0. Moreover, the right panel of figure 1 compares the ground truth Δθ θ with the change estimated by the IF. When ε shifts from 0 to −13- 13, θ moves along the boundary of the feasible region, whereas the estimated change by the IF significantly deviates from the ground truth due to its failure to consider the constraint. In contrast, our estimator DIF (see Definition 2), which accounts for the geometry of the feasible region, closely matches the ground truth Δθ θ along the boundary. The derivation of DIF for this toy example is provided in Appendix C. 3 Directional Influence Function (DIF) To address the non-smoothness of the solution map θ(ε)θ( ) w.r.t. perturbation ε and the inability of IF to handle constrained scenarios, this section proposes an alternative approach called the Directional Influence Function (DIF). We start with the definition of the classical IF. Definition 1 (Influence Function). The classical IF applies only to unconstrained problems. The IF is defined as the derivative of θ w.r.t. the ε0 _0, i.e., IF(θ¯;Zr)=dθ(ε0)dε0|ε0=0.IF( θ;Z^r)= . dθ( _0)d _0 |_ _0=0. (5) Definition 2 (Directional Influence Function). The DIF of the data subset ZrZ^r is defined as the limiting directional derivative of the solution map θ(ε):ℝm+1→ℝdθ( ):R^m+1 ^d at ε¯= =0, if the limit exists222This generalized definition accounts for potential non-smoothness of θ(ε)θ( ) by allowing the direction to vary in a neighborhood of Δε¯ , as is common in nonsmooth analysis(Clarke, 1990; Rockafellar and Wets, 1998). DIF(θ¯;Zr) ( θ;Z^r ) :=Dθ(ε¯;Δε¯)=limt↘0+Δε→Δε¯θ(ε¯+tΔε¯)−θ(ε¯)t, :=Dθ( ; )= _ subarrayct 0^+\\ → subarray θ( +t )-θ( )t, (6) where θ¯=θ(ε¯) θ=θ( ) denotes the model solution before perturbation. Δε¯=ε^−ε¯=[−1N0,−1N1,…,−1Nm] = - = [- 1N_0,- 1N_1,…,- 1N_m ] is the perturbation direction, where Δε¯ corresponds to removing the data subset ZrZ^r from the constrained learning (1). Recall that a function f:X→Yf:X→ Y is called positively homogeneous if f(αx)=αf(x)∀α≥0,∀x∈X.f(α x)=α f(x) ∀α≥ 0,∀ x∈ X. Corollary 3. DIF is positively homogeneous,i.e., αDθ(ε¯;Δε¯)=Dθ(ε¯;αΔε¯),α∈ℝ+α Dθ( ; )=Dθ( ;α ), α ^+ Proof. See Appendix B.4 for the full proof. ∎ When ε is perturbed from ε¯ by a finite Δε¯ , IF estimates Δθ θ via a first-order Taylor expansion. IF pertains to unconstrained learning and ignores constraints; thus Δθ≈dθ(ε0)dε0|ε=ε¯⋅Δε¯0. θ≈ . dθ( _0)d _0 |_ = · _0. (7) This approximation, as demonstrated in the previous section, fails in the presence of constraints. Instead, DIF utilizes a directional derivative approach: Δθ θ ≈Dθ(ε¯;Δε¯‖Δε¯‖)⋅‖Δε¯‖=Dθ(ε¯;Δε¯). ≈ Dθ ( ; \| \| )·\| \|=Dθ ( ; ). (8) where Dθ(ε¯;Δε¯‖Δε¯‖)Dθ( ; \| \|) is the directional derivative, and the second equation follows from Corollary 3. In what follows, we will demonstrate how to compute Dθ(ε¯;Δε¯)Dθ ( ; ) by solving a linearized variational inequality system. We use Δθ θ to denote the DIF-based estimate of Δθ θ. 4 Deriving DIF via Sensitivity Analysis of Optimality Conditions In this section, we show how DIF naturally arises from the sensitivity analysis of the optimality conditions of the constrained learning problem. Specifically, we first formulate the optimality system as a VI, then linearize this VI to obtain an auxiliary VI, and finally convert it into a QP whose solution yields the desired directional change in the model solution. Recall that the constrained learning is typically solved using Lagrangian dual approaches, such as the augmented Lagrangian method, primal-dual methods (Chamon and Ribeiro, 2020; Ahmed et al., 2022; Nandwani et al., 2019). To quantify the impact of data downweighting on the solution of ( 2), we introduce the perturbation ε into the Lagrangian function, enabling us to analyze how the perturbation affects the optimality system. The Lagrangian function of problem (2) is: L(ε,θ,λ)=1N0∑i=1N0ℓ0(zi(0),θ)+ε0∑zi(0)∈Zrℓ0(zi(0),θ) L( ,θ,λ)= 1N_0 _i=1^N_0 _0 (z_i^(0),θ )+ _0 _z_i^(0)∈ Z^r _0 (z_i^(0),θ ) (9) +∑j=1mλj[1Nj∑i=1Njℓj(zi(j),θ)+εj∑zi(j)∈Zrℓj(zi(j),θ)−τj], + _j=1^m _j [ 1N_j _i=1^N_j _j (z_i^(j),θ )+ _j _z_i^(j)∈ Z^r _j (z_i^(j),θ )- _j ], where λ=[λ1,λ2,…,λm],λj≥0λ=[ _1, _2,…, _m], _j≥ 0 denotes the dual variables. Let (θ(ε),λ(ε))(θ( ),λ( )) be the primal–dual solution for a given ε . From Section 2, we use (θ¯,λ¯)=(θ(ε¯),λ(ε¯))( θ, λ)=(θ( ),λ( )) to denote the optimal solution of the constrained learning problem (1). At θ=θ¯θ= θ, we partition the constraints into the active set IActive I_Active and the inactive set IInactiveI_Inactive, where IActive =IBinding∪INon-bindingI_Active =I_Binding∪ I_Non-binding. IActive I_Active :=j|1Nj∑i=1Njℓj(zi(j),θ¯)=τj, := \j\, |\, 1N_j _i=1^N_j _j(z_i^(j), θ)= _j \, IInactive I_Inactive :=j|1Nj∑i=1Njℓj(zi(j),θ¯)<τj, := \j\, |\, 1N_j _i=1^N_j _j(z_i^(j), θ)< _j \, IBinding I_Binding :=j∈IActive∣λ¯j>0, :=\j∈ I_Active λ_j>0\, INon-binding I_Non-binding :=j∈IActive∣λ¯j=0. :=\j∈ I_Active λ_j=0\. As discussed earlier, varying ε from ε¯ to ε corresponds to removing ZrZ^r. In what follows, we study the sensitivity of the optimality conditions with respect to ε . We begin by examining the optimality condition of problem (2). 4.1 Optimality Condition of Problem (2) Rather than relying on the Karush-Kuhn-Tucker (KKT) conditions, we adopt a VI formulation to characterize the optimality conditions. This VI form allows us to analyze solution perturbations through the lens of generalized differential calculus (Rockafellar and Wets, 1998). The optimality condition of problem (2) is: −∇θL(ε,θ,λ)∈Nℝd(θ),∇λL(ε,θ,λ)∈Nℝ+m(λ),- _θL( ,θ,λ)∈ N_R^d(θ), _λL( ,θ,λ)∈ N_R_+^m(λ), (10) where Nℝd(θ),Nℝ+m(λ)N_R^d(θ),N_R_+^m(λ) denote the Normal Cone. Definition 4 (Normal Cone;Dontchev and Rockafellar (2009)). Let C⊆ℝdC ^d be a closed convex set. The normal cone to C at a point x∈Cx∈ C is defined as: NC(x)=v∈ℝd:⟨v,y−x⟩≤0,∀y∈C.N_C(x)=\v ^d: v,y-x ≤ 0,\ ∀ y∈ C\. 4.2 Linear Approximation of Variational Inequality The classical derivation of the IF relies on the Taylor expansion of the first-order optimality condition, which does not naturally handle constraints. In constrained learning, a natural choice for the optimality condition is the KKT system. Yet, estimating how ε affects the KKT system is challenging due to its mixed system of equations and inequalities. DIF adopts an elegant alternative: representing the optimality condition as a VI (see formulation (10)) and applying a first-order approximation directly to the VI. The VI (10) can be compactly written as: f(ε,θ,λ)+NE(θ,λ)∋, f( ,θ,λ)+N_E(θ,λ) 0, (11) wheref(ε,θ,λ)=(∇θL(ε,θ,λ)−∇λL(ε,θ,λ)),E=ℝd×ℝ+m f( ,θ,λ)= pmatrix _θL( ,θ,λ)\\ - _λL( ,θ,λ) pmatrix,E=R^d×R_+^m By definition, the solution (ε¯,θ¯,λ¯)( , θ, λ) satisfies the equation (11). When ε is perturbed by Δε¯ , estimating Δθ θ and Δλ λ amounts to finding Δθ θ and Δλ λ such that f(ε¯+Δε¯,θ¯+Δθ,λ¯+Δλ)+NE(θ¯+Δθ,λ¯+Δλ)∋f( + , θ+ θ, λ+ λ)+N_E( θ+ θ, λ+ λ) 0 (12) The linear approximation of (12) is f(ε¯,θ¯,λ¯)+∇εf(ε¯,θ¯,λ¯)Δε¯+∇(θ,λ)f(ε¯,θ¯,λ¯)[Δθ^Δλ^] f( , θ, λ)+ _ f( , θ, λ)\, + _(θ,λ)f( , θ, λ)\, [ array[]l θ\\ λ array ] (15) +NE((θ¯,λ¯)+(Δθ^,Δλ^))∋. +N_E(( θ, λ)+( θ, λ)) 0. (16) Here, (Δθ^,Δλ^)( θ, λ), which satisfies (16), serves as an approximation of (Δθ,Δλ)( θ, λ). Denote A:=∇(θ,λ)f(ε¯,θ¯,λ¯)=[∇θ2L(ε¯,θ¯,λ¯),∇θλL(ε¯,θ¯,λ¯)−∇θλL(ε¯,θ¯,λ¯),−∇λ2L(ε¯,θ¯,λ¯)], A= _(θ,λ)f( , θ, λ)= [ array[]l _θ^2L( , θ, λ), _θλL( , θ, λ)\\ - _θλL( , θ, λ),-∇^2_λL( , θ, λ) array ], μ:=∇εf(ε¯,θ¯,λ¯)=[∇θεL(ε¯,θ¯,λ¯)−∇λεL(ε¯,θ¯,λ¯)] μ= _ f( , θ, λ)= [ array[]l _θ L( , θ, λ)\\ - _λ L( , θ, λ) array ] Δη^=[Δθ^Δλ^] η= [ array[]l θ\\ λ array ] The linearized VI (16) can be expressed in the form: f(ε¯,θ¯,λ¯)+μΔε¯+AΔη^+NE((θ¯,λ¯)+Δη^)∋0.f( , θ, λ)+μ +A η+N_E(( θ, λ)+ η) 0. (17) Note that all terms in the linearized VI (17), except for Δη η—which is the variable we aim to estimate—can be precomputed from the known solution θ¯ θ and λ¯ λ. The linearized VI (17) can be simplified into the following Auxiliary VI (18). Proposition 5 (Auxiliary VI). The solution Δη η to the linearized VI (17) satisfies μΔε¯+AΔη^+NK(Δη^)∋0,μ +A η+N_K( η) 0, (18) where K=ℝd×DK=R^d× D. D is a space defined as D:=Δλ^∈ℝm|Δλ^j∈ℝfor j∈IBinding,Δλ^j≥0for j∈INon-binding,Δλ^j=0for j∈IInactive.D:= \ λ ^m\; |\; aligned & λ_j j∈ I_Binding,\\ & λ_j≥ 0 j∈ I_Non-binding,\\ & λ_j=0 j∈ I_Inactive aligned \. Proof. See Appendix B.3 ∎ We refer to (18) as the Auxiliary VI, as it admits the same solution as the linearized VI (17), but has a simpler form. Proposition 6. Let Δη^=[Δθ^;Δλ^] η=[ θ; λ] denote the solution to Auxiliary VI (18). Then, Δθ θ exactly recovers the directional derivative of the solution mapping θ(ε)θ( ) at ε¯ along Δε¯ , i.e., Δθ^=Dθ(ε¯;Δε¯). θ=Dθ( ; ). (19) Proof. See Appendix B.5. ∎ Remark 7. By Proposition 6, computing the DIF reduces to solving the Auxiliary VI (18). The following theorem establishes the existence of DIF. Theorem 8. Assume the following regularity conditions hold: Assumption 1. (LICQ: Linear Independence Constraint Qualification) The gradients 1Nj∑i=1Nj∇θℓj(zi(j),θ¯) 1N_j _i=1^N_j _θ _j(z_i^(j), θ) associated with the active constraints (j∈IActivatej∈ I_Activate) are linearly independent. Assumption 2. (SOSC: Second-Order Sufficient Condition) For any Δθ≠0 θ≠ 0 such that Δθ⟂1Nj∑i=1Nj∇θℓj(zi(j),θ¯) θ 1N_j _i=1^N_j _θ _j(z_i^(j), θ) for all j∈IBindingj∈ I_Binding, it holds that ⟨Δθ,∇θθ2L(ε¯,θ¯,λ¯)Δθ⟩>0 θ,∇^2_θL( , θ, λ) θ >0. Then the directional derivative Dθ(ε¯;Δε¯)Dθ( ; ) exists. Proof. See Appendix B.2. ∎ Proposition 9. Δθ=θ(ε¯+Δε¯)−θ(ε¯) θ=θ( + )-θ( ) is the ground truth and Δθ θ is the estimation. For any perturbation Δε¯ sufficiently small, there exists an M s.t. ‖Δθ−Δθ^‖≤M‖Δε¯‖\| θ- θ\|≤ M\| \| (20) 4.3 Computing DIF via Quadratic Programming Proposition 6 implies that computing the DIF reduces to solving the Auxiliary VI (18). In this section, we construct a QP 21 whose optimality is exactly the Auxiliary VI (18). By solving this QP, we can compute DIF and estimate the directional change in solution, i.e., Δθ θ. Figure 2 summarizes the derivation in this section, illustrating the connections between the problem (2), VI (11), Auxiliary VI (18), and QP (21) formulations. PCL (2)VI (11)Auxiliary VI (18) QP (21) Optimality Condition Equivalent Theorem 10 Equivalent LinearizationProposition 5Δθ θΔλ λ Figure 2: The relationship between (2) and (21). We define the following QP: min ω L(ε¯,θ¯,λ¯)+⟨∇θεL(ε¯,θ¯,λ¯)Δε¯,ω⟩+12⟨ω,∇θθ2L(ε¯,θ¯,λ¯)ω⟩ L( , θ, λ)+ _θ L( , θ, λ) ,\;ω + 12 ω,\; _θ^2L( , θ, λ)\,ω (21) s.t. 1Nj∑i=1Njℓj(zi(j),θ¯)+1Nj∑i=1Nj∇θℓj(zi(j),θ¯)⋅ω+∑zi(j)∈ZrΔε¯jℓj(zi(j),θ¯)−τj=0,j∈Ibinding,≤0,j∈Inon-binding,free,j∈IInactive. 1N_j _i=1^N_j _j(z_i^(j), θ)+ 1N_j _i=1^N_j _θ _j(z_i^(j), θ)\!·\!ω+ _z_i^(j)∈ Z^r _j _j(z_i^(j), θ)- _j Theorem 10 (Auxiliary problem). Let (ω⋆,ζ⋆) (ω ,ζ ) denote the optimal primal-dual solution to the QP (21). Then (w⋆,ζ⋆) (w ,ζ ) also satisfies the VI (18). In particular, the QP (21) and the auxiliary VI (18) admit the same solution pair under substitution (w∗,ζ∗)=(Δθ^,Δλ^)(w^*,ζ^*)=( θ, λ). Proof. See Appendix B.6. ∎ Proposition 11. If no constraints are active at the solution θ¯ θ, i.e., IActive=∅I_Active= , then the DIF coincides with the classical IF. 5 Validation 5.1 DIF Validating via Constrained Linear Regression Constrained linear regression is a widely applicable instance of learning with hard constraints. It arises in several domains: in portfolio optimization, asset weights must be non-negative and sum to one, i.e., θi≥0,∑iθi=1 _i≥ 0, _i _i=1; in traffic flow modeling, flow variables must satisfy conservation laws and remain below capacity, e.g., Aθ=0,θ≤cAθ=0,θ≤ c; and in fair machine learning, additional linear conditions are imposed to enforce equity across groups, such as Afair θ≤bfair A_fair θ≤ b_fair . This section leverages a constrained linear regression to evaluate the DIF estimator. We repeatedly remove one data point (100 trials) and compare the solution change Δθ θ predicted by DIF with the ground-truth retraining solution. We include the classical IF estimator and the penalty-based IF approximation as baselines. Data generation. We generate a synthetic regression dataset. Specifically, we draw n=1000n=1000 samples with d=5d=5 features: X∈ℝn×d,Xij∼(0,1),θ∗∼(0,Id),y=Xθ∗+ε,ε∼(0,0.12In). X ^n× d,X_ij (0,1), θ^* (0,I_d ), y=Xθ^*+ , (0,0.1^2I_n ). Constrained Linear Regression. The regression parameters are obtained by solving the constrained least-squares problem: θ^=argminθ∈ℝ512n‖Xθ−y‖22 s.t.Aeq θ=beq ,Aineq θ≤bineq θ= _θ ^5 12n\|Xθ-y\|_2^2 s.t. A_eq θ=b_eq , A_ineq θ≤ b_ineq where Aeq=[1,1,1,1,1],beq=−4.0,Aineq=[1100000−101],bineq=[1.51.5]. gatheredA_eq=[1,1,1,1,1], b_eq=-4.0, A_ineq= [ array[]c1&1&0&0&0\\ 0&0&-1&0&1 array ], b_ineq= [ array[]c1.5\\ 1.5 array ]. gathered Experimental results. We perform 100 single-point removals. For each trial, we estimate the resulting change Δθ θ and compare it with the ground-truth leave-one-out (LOO) retraining. We leverage three estimators: i) The IF estimates the solution change by ignoring the constraints, yielding ΔθIF _IF as in equation (7). i) The penalty-based IF adds soft penalties for the constraints to the objective and applies the IF estimator on the penalized surrogate. i) The DIF enforces feasibility by solving the QP (21), yielding ΔθDIF _DIF as in equation (8). Appendix D provides the full derivations for all three methods in the constrained linear regression setting. All experiments are solved with CVXPY. As shown in Fig. 3, DIF nearly coincides with LOO (points align with y=xy=x), whereas IF and penalty IF exhibit noticeable bias. (a) (b) Figure 3: Comparison of solution changes on a constrained linear regression task. (a) IF and the penalty-based IF estimators versus the ground truth leave-one-out retraining results. (b) Estimation of the proposed DIF versus the same ground-truth values. DIF aligns almost perfectly with the y=xy=x line, demonstrating its ability to accurately capture parameter changes under constraints. 5.2 DIF Validation via Constrained CNN Model We adopt the following constrained learning formulation (Shen et al., 2022): minθ∈ΘR¯(fθ)s.t.Ri(fθ)−R¯(fθ)−τ≤0,i=1,…,m, _θ∈ \; R(f_θ) .t. R_i(f_θ)- R(f_θ)-τ≤ 0,\;\;i=1,…,m, (22) where Ri(fθ)=(x,y)∼i[ℓ(fθ(x),y)]R_i(f_θ)=E_(x,y) _i[ (f_θ(x),y)] is the risk of the client (or group) i, and R¯(fθ)=1C∑i=1CRi(fθ) R(f_θ)= 1C _i=1^CR_i(f_θ) is the average risk. The constraints ensure that no group’s risk exceeds the global average risk by more than τ, thereby limiting the performance disparity across groups. We use a network with seven convolutional layers and tanh(⋅) (·) non-linearities, modeled after the all convolutional network of Springenberg et al. (2014). We train it on the MNIST training set (LeCun et al., 1998). We solve this problem using a primal–dual method, jointly updating the model parameters θ and the Lagrange multipliers. Heterogeneous Data Partitioning. We create non-IID group-wise data partitions following Shen et al. (2022) . We split the dataset into C groups by allocating an α-fraction of the data uniformly at random and distributing the remaining (1−α)(1-α)-fraction in a label-skewed way, where samples are sorted by class labels and assigned consecutively to groups. Unless otherwise specified, we set the number of groups to m=3m=3 and α=0.5α=0.5. This setup creates label-imbalanced distributions across different groups. Influence Estimation. We employ DIF to estimate the effect of removing training samples under the fairness constraint. We select the 100 most influential training samples and remove one sample at a time (100 trials in total). Since CNNs are typically non-convex, the solution may shift to another valley in the loss landscape during retraining. Therefore, instead of directly comparing the predicted solution change ΔθDIF _DIF with the ground-truth change ΔθLOO _LOO, we evaluate DIF indirectly by comparing their resulting loss changes on a misclassified test point. Specifically, we approximate the updated model as θ′≈θ+ΔθDIFθ ≈θ+ _DIF and compute the predicted loss difference Δℓ=R(fθ′(xtest),ytest)−R(fθ(xtest),ytest). =R (f_θ (x_test),y_test )-R (f_θ(x_test),y_test ). Figure 4 compares these DIF-predicted loss differences with the actual loss differences obtained by retraining the model. The points align closely with the y=xy=x line (Pearson r=0.90r=0.90), indicating that DIF accurately predicts the influence of individual training samples on this test loss. Figure 4: Actual vs. DIF predicted loss differences on a misclassified test sample. 6 Conclusion References K. Ahmed, T. Li, T. Ton, Q. Guo, K. Chang, P. Kordjamshidi, V. Srikumar, G. Van den Broeck, and S. Singh (2022) Pylon: a pytorch framework for learning with constraints. In NeurIPS 2021 Competitions and Demonstrations Track, p. 319–324. Cited by: §4. L. Chamon and A. Ribeiro (2020) Probably approximately correct constrained learning. Advances in Neural Information Processing Systems 33, p. 16722–16735. Cited by: §1, §4. F. H. Clarke (1990) Optimization and nonsmooth analysis. SIAM. Cited by: footnote 2. A. L. Dontchev and R. T. Rockafellar (2009) Implicit functions and solution mappings. Vol. 543, Springer. Cited by: Definition 12, Lemma 16, Remark 19, Definition 4. V. Feldman and C. Zhang (2020) What neural networks memorize and why: discovering the long tail via influence estimation. Advances in Neural Information Processing Systems 33, p. 2881–2891. Cited by: §1. J. Garcıa and F. Fernández (2015) A comprehensive survey on safe reinforcement learning. Journal of Machine Learning Research 16 (1), p. 1437–1480. Cited by: §1. F. R. Hampel (1974) The influence curve and its role in robust estimation. Journal of the american statistical association 69 (346), p. 383–393. Cited by: §1. I. Hounie, A. Ribeiro, and L. F. Chamon (2023) Resilient constrained learning. Advances in Neural Information Processing Systems 36, p. 71767–71798. Cited by: §1. P. W. Koh and P. Liang (2017) Understanding black-box predictions via influence functions. In International conference on machine learning, p. 1885–1894. Cited by: §1. Y. LeCun, L. Bottou, Y. Bengio, and P. Haffner (1998) Gradient-based learning applied to document recognition. Proceedings of the IEEE 86 (11), p. 2278–2324. Cited by: §5.2. Z. Li, B. Liu, Z. Yang, Z. Wang, and M. Wang (2023) Double duality: variational primal-dual policy optimization for constrained reinforcement learning. Journal of Machine Learning Research 24 (385), p. 1–43. Cited by: §1. J. Lin, L. Tao, M. Dong, and C. Xu (2024) Diffusion attribution score: evaluating training data influence in diffusion models. arXiv preprint arXiv:2410.18639. Cited by: §1. Y. Nandwani, A. Pathak, and P. Singla (2019) A primal dual formulation for deep learning with constraints. Advances in neural information processing systems 32. Cited by: §4. G. Pruthi, F. Liu, S. Kale, and M. Sundararajan (2020) Estimating training data influence by tracing gradient descent. Advances in Neural Information Processing Systems 33, p. 19920–19930. Cited by: §1. R. T. Rockafellar and R. J. Wets (1998) Variational analysis. Springer. Cited by: §4.1, Definition 13, Definition 14, footnote 2. Z. Shen, J. Cervino, H. Hassani, and A. Ribeiro (2022) An agnostic approach to federated learning with class imbalance.. In ICLR, Cited by: §5.2, §5.2. J. T. Springenberg, A. Dosovitskiy, T. Brox, and M. Riedmiller (2014) Striving for simplicity: the all convolutional net. arXiv preprint arXiv:1412.6806. Cited by: §5.2. J. Xia, Y. Yang, J. Shen, S. Wang, and J. Cao (2025) FairTP: a prolonged fairness framework for traffic prediction. In Proceedings of the AAAI Conference on Artificial Intelligence, Vol. 39, p. 26391–26399. Cited by: §1. Y. Zhang, J. Shen, X. Xiong, and Y. Kwon (2024) Timeinf: time series data contribution via influence functions. arXiv preprint arXiv:2407.15247. Cited by: §1. Appendix A Preliminary Definition 12 (Polar Cone; cf. Dontchev and Rockafellar (2009)). Let K⊆ℝnK ^n be a closed convex cone. The polar cone of K is defined as K∗=y∈ℝn∣⟨x,y⟩≤0∀x∈K.K^*=\\,y ^n\; \; x,y ≤ 0\;\;∀ x∈ K\,\. (23) Moreover, the normal vectors to K and K∗K^* satisfy y∈NK(x)⟺x∈NK∗(y)⟺x∈K,y∈K∗,⟨x,y⟩=0.y∈ N_K(x)\;\; \;\;x∈ N_K^*(y)\;\; \;\;x∈ K,\;y∈ K^*,\; x,y =0. (24) Definition 13 (Tangent Cone, cf. Rockafellar and Wets (1998)). For a set C⊂ℝnC ^n (not necessarily convex) and a point x∈Cx∈ C, a vector v is said to be tangent to C at x if 1τk(xk−x)→vfor some xk→x,xk∈C,τk↓0. 1τ^k(x^k-x)\;→\;v some x^k→ x,\;x^k∈ C,\;τ^k 0. (25) The set of all such vectors v is called the tangent cone to C at x and is denoted by TC(x)T_C(x). For x∉Cx∉ C, we take TC(x)=∅T_C(x)= . Definition 14 (Critical Cone, cf. Rockafellar and Wets (1998)). For a convex set C, any x∈Cx∈ C and any v∈NC(x)v∈ N_C(x), the critical cone to C at x for v is KC(x,v)=w∈TC(x)∣w⟂v.K_C(x,v)=\\,w∈ T_C(x) w v\,\. (26) Definition 15 (Critical Subspaces). For a convex set C⊂ℝnC ^n and (x,v)(x,v) with v∈NC(x)v∈ N_C(x), let KC(x,v)K_C(x,v) denote the critical cone at (x,v)(x,v). Then the associated critical subspaces are defined as KC+(x,v)=KC(x,v)−KC(x,v)=w−w′∣w,w′∈KC(x,v),K_C^+(x,v)\;=\;K_C(x,v)-K_C(x,v)\;=\;\\,w-w w,w ∈ K_C(x,v)\,\, (27) KC−(x,v)=KC(x,v)∩[−KC(x,v)]=w∈KC(x,v)∣−w∈KC(x,v).K_C^-(x,v)\;=\;K_C(x,v)∩[-K_C(x,v)]\;=\;\\,w∈ K_C(x,v) -w∈ K_C(x,v)\,\. (28) Here KC+(x,v)K_C^+(x,v) is the smallest linear subspace containing KC(x,v)K_C(x,v), and KC−(x,v)K_C^-(x,v) is the largest linear subspace contained in KC(x,v)K_C(x,v). Lemma 16 (Reduction Lemma;cf. 2E.4 Dontchev and Rockafellar (2009)). Let C⊂ℝnC ^n be a convex set and x¯∈C x∈ C with v¯∈NC(x¯) v∈ N_C( x). Define the critical cone KC:=KC(x¯,v¯)K_C:=K_C( x, v). Then, for all sufficiently small w,u∈ℝnw,u ^n, we have the equivalence v¯+Δv∈NC(x¯+Δx)⟺Δv∈NKC(Δx). v+ v∈ N_C( x+ x)\ \ v∈ N_K_C( x). (29) Appendix B Proofs B.1 Auxiliary Lemmas For clarity, we first introduce some notation and rewrite the generalized equation in a standard form. Notation. Let η=(θ,λ)η=(θ,λ) denote the parameter vector, η¯=(θ¯,λ¯) η=( θ, λ) be its reference value, and Δη^=(Δθ^,Δλ^) η=( θ, λ) be the estimation to Δη=(Δθ,Δλ) η=( θ, λ). Solution mapping. We define the solution mapping as S(ε):=(θ,λ)∣f(ε,θ,λ)+NE(θ,λ)∋0=η∣f(ε,η)+NE(η)∋0,S( ):=\\,(θ,λ) f( ,θ,λ)+N_E(θ,λ) 0\,\=\\,η f( ,η)+N_E(η) 0\,\, (30) where NE(η)N_E(η) denotes the normal cone to E at η. Generalized equation form. To study the local behavior of S, we rewrite the system as the generalized equation G(η):=f(ε¯,η¯)+∇ηf(ε¯,η¯)(η−η¯)+NE(η),G(η):=f( , η)+ _ηf( , η)(η- η)+N_E(η), (31) with ∇ηf(ε¯,η¯)=∇(θ,λ)f(θ¯,λ¯,x¯)=:A. _ηf( , η)= _(θ,λ)f( θ, λ, x)=:A. (32) At the reference point, the optimality condition (11) ensures that G(η¯)∋0.G( η) 0. (33) Linearization. We then consider the linearized generalized equation (the first-order approximation of G): G0(Δη):=∇ηf(ε¯,η¯)Δη+NK(Δη)=AΔη+NK(Δη),G0(0)∋0,G_0( η):= _ηf( , η) η+N_K( η)=A η+N_K( η), G_0(0) 0, (34) where K:=KE(η¯,−f(ε¯,η¯))K:=K_E( η,-f( , η)) is the critical cone. This coincides with K=ℝd×DK=R^d× D as introduced in Proposition 5. It is convenient to define the inverse-type mapping s¯:=G0−1=(A+NK)−1, s:=G_0^-1=(A+N_K)^-1, (35) which solves the linearized inclusion (34). Given perturbation uΔε¯u , s¯(−uΔε¯)=(A+NK)−1(−uΔε¯)=Δη^∣AΔη^+NK(Δη^)+uΔε¯∋0. s(-u )=(A+N_K)^-1(-u )=\ η A η+N_K( η)+u 0\. (36) In other words, s¯(−uΔε¯) s(-u ) yields the solution to the auxiliary VI (18). Critical subspace of E. Following the definition 15, the critical subspaces of E at (θ¯,λ¯,−f(ε¯,θ¯,λ¯))( θ, λ,-f( , θ, λ)) are defined as KE+(θ¯,λ¯,−f(ε¯,θ¯,λ¯)) K_E^+( θ, λ,-f( , θ, λ)) =ℝd×Kℝm+(λ¯,∇λL(ε¯,θ¯,λ¯)), =R^d× K_R^m^+( λ, _λL( , θ, λ)), (37) KE−(θ¯,λ¯,−f(ε¯,θ¯,λ¯)) K_E^-( θ, λ,-f( , θ, λ)) =ℝd×Kℝm−(λ¯,∇λL(ε¯,θ¯,λ¯)). =R^d× K_R^m^-( λ, _λL( , θ, λ)). (38) Δη^=(Δθ^,Δλ^)satisfiesΔη^∈KE+⟺Δλ^j=0,∀j∈IInactive,and Δθ^∈ℝd. η=( θ, λ)\ satisfies\ η∈ K_E^+\ \ λ_j=0,\ ∀ j∈ I_Inactive,\ and θ ^d. (39) Lemma 17. Under the assumptions of Theorem 8, for any Δη η satisfying the linearized VI (18), if Δη^∈KE+,Δη^≠0,AΔη^⟂KE−, η∈ K_E^+, η≠ 0, A η K_E^-, (40) then ⟨Δη^,AΔη^⟩>0. η,\;A η >0. Proof. For any Δη η satisfying the linearized VI (18), it follows from the definition of the space D that Δλ^j=0 λ_j=0 for all j in the inactive set. Therefore, any Δη η satisfying VI (15) also satisfies Δη^∈KE+ η∈ K_E^+. We next analyze the orthogonality condition AΔη^⟂KE−A η K_E^-. By expanding A(Δθ^,Δλ^)A( θ, λ) and splitting the θ- and λ-parts, we obtain: AΔη^⟂KE− A η K_E^- ⟺A(Δθ^,Δλ^)⟂KE− \ \ A( θ, λ) K_E^- ⟺(HΔθ^+∇θλL(ε¯,θ¯,λ¯)Δλ^,−∇λθL(ε¯,θ¯,λ¯)Δθ^)⟂KE− \ \ (H θ+ _θλL( , θ, λ)\, λ,\;- _λθL( , θ, λ)\, θ ) K_E^- ⟺HΔθ^+∇θλL(ε¯,θ¯,λ¯)Δλ^⟂ℝd,−∇λθL(ε¯,θ¯,λ¯)Δθ^⟂Kℝm−(λ¯,∇λL(ε¯,θ¯,λ¯)) \ \ H θ+ _θλL( , θ, λ)\, λ ^d, - _λθL( , θ, λ)\, θ K_R^m^-\! ( λ, _λL( , θ, λ) ) ⟺HΔθ^+∇θλL(ε¯,θ¯,λ¯)Δλ^=0,−∇λθL(ε¯,θ¯,λ¯)Δθ^⟂Kℝm−(λ¯,∇λL(ε¯,θ¯,λ¯)). \ \ H θ+ _θλL( , θ, λ)\, λ=0, - _λθL( , θ, λ)\, θ K_R^m^-\! ( λ, _λL( , θ, λ) ). Moreover, by the definition of the critical subspace, we have Δλ^∈Kℝm−(λ¯,∇λL(ε¯,θ¯,λ¯))⟺Δλ^j=0,∀j∈Inon-binding∪IInactive, λ∈ K_R^m^-\! ( λ, _λL( , θ, λ) )\ \ λ_j=0,\ ∀ j∈ I_non-binding∪ I_Inactive, (41) where IInactiveI_Inactive and Inon-bindingI_non-binding denote the sets of inactive and non-binding constraints, respectively. Consequently, the condition −∇λθL(ε¯,θ¯,λ¯)Δθ^⟂Kℝm−(λ¯,∇λL)- _λθL( , θ, λ)\, θ\ \ K_R^m^-\! ( λ, _λL ) only requires that ∇λjθLΔθ^=0,∀j∈Ibinding. _ _jθL\, θ=0, ∀ j∈ I_binding. That is, AΔη^⟂KE− A η K_E^- ⟺HΔθ^+∇θλL(ε¯,θ¯,λ¯)Δλ^=0,∇λjθLΔθ^=0,∀j∈Ibinding. \ \ H θ+ _θλL( , θ, λ)\, λ=0, _ _jθL\, θ=0,\ ∀ j∈ I_binding. (42) Notice that ∇λjθLΔθ^=0⟺Δθ^⟂1Nj∑i=1Nj∇θℓj(zi(j),θ¯),∀j∈Ibinding. _ _jθL\, θ=0\ \ θ 1N_j _i=1^N_j _θ _j (z_i^(j), θ ), ∀ j∈ I_binding. (43) Assumption 2 in Theorem 8 implies that for all Δθ θ satisfying (42), we have ⟨Δθ^,∇θθ2L(ε¯,θ¯,λ¯)Δθ^⟩>0. θ, _θ^2L( , θ, λ) θ >0. (44) Moreover, substituting A:=∇(θ,λ)f(ε¯,θ¯,λ¯)=[∇θ2L(ε¯,θ¯,λ¯)∇θλL(ε¯,θ¯,λ¯)−∇θλL(ε¯,θ¯,λ¯)−∇λ2L(ε¯,θ¯,λ¯)],A:= _(θ,λ)f( , θ, λ)= bmatrix _θ^2L( , θ, λ)& _θλL( , θ, λ)\\ - _θλL( , θ, λ)&- _λ^2L( , θ, λ) bmatrix, to (44), we obtain ⟨Δη^,AΔη^⟩=⟨(Δθ^,Δλ^),A(Δθ^,Δλ^)⟩=⟨Δθ^,∇θθ2L(ε¯,θ¯,λ¯)Δθ^⟩. η,\;A η = ( θ, λ),\,A( θ, λ) = θ,\,∇^2_θL( , θ, λ)\, θ . (45) Combining the two results yields ⟨Δη^,AΔη^⟩>0 η,A η >0. ∎ Lemma 18. Under the assumptions of Theorem 8, s¯:=(A+NK)−1 s:=(A+N_K)^-1 is everywhere single-valued. Equivalently, there exists a unique solution Δη η to the auxiliary VI (18) given (ε¯,θ¯,λ¯)( , θ, λ) and Δε . Proof. We prove that G0−1:=(A+NK)−1G_0^-1:=(A+N_K)^-1 is single-valued everywhere. Assume, for contradiction, that there exist two distinct solutions Δη^1 η_1 and Δη^2 η_2 such that r=G0−1(Δη^1)=G0−1(Δη^2)r=G_0^-1( η_1)=G_0^-1( η_2). Then, Δη^1,Δη^2∈K,r−AΔη^1∈NK(Δη^1),r−AΔη^2∈NK(Δη^2). η_1, η_2∈ K, r-A η_1∈ N_K( η_1), r-A η_2∈ N_K( η_2). This implies, by the definition of the normal cone, that ⟨Δη^1,r−AΔη^2⟩≤0,⟨Δη^2,r−AΔη^1⟩≤0. η_1,\,r-A η_2 ≤ 0, η_2,\,r-A η_1 ≤ 0. (46) Note that K is the critical cone. By condition (24), we have Δη^1 η_1 ∈K,r−AΔη^1∈K∗,⟨Δη^1,r−AΔη^1⟩=0, ∈ K, r-A η_1∈ K^*, η_1,\;r-A η_1 =0, (47) Δη^2 η_2 ∈K,r−AΔη^2∈K∗,⟨Δη^2,r−AΔη^2⟩=0. ∈ K, r-A η_2∈ K^*, η_2,\;r-A η_2 =0. By the definition of the polar cone K∗=(K−)⟂K^*=(K^-) , we also have −A(Δη^1−Δη^2)∈K∗−K∗=(K−)⟂,Δη^1−Δη^2∈K−K=K+.-A( η_1- η_2)∈ K^*-K^*=(K^-) , η_1- η_2∈ K-K=K^+. Consider the following inner product: ⟨Δη^1−Δη^2,A(Δη^1−Δη^2)⟩ η_1- η_2,\;A( η_1- η_2) =⟨Δη^1−Δη^2,[r−AΔη^2]−[r−AΔη^1]⟩ = η_1- η_2,\;[r-A η_2]-[r-A η_1] =⟨Δη^1,r−AΔη^2⟩−⟨Δη^1,r−AΔη^1⟩ = η_1,\;r-A η_2 - η_1,\;r-A η_1 −⟨Δη^2,r−AΔη^2⟩+⟨Δη^2,r−AΔη^1⟩ - η_2,\;r-A η_2 + η_2,\;r-A η_1 ≤0. ≤ 0. Following Lemma 17, Assumption 2 of Theorem 8 ensures that for any nonzero Δη^∈K+ η∈ K^+ with AΔη^⟂K−A η K^-, we must have ⟨Δη^,AΔη^⟩>0 η,\;A η >0. This contradicts the above inequality unless Δη^1=Δη^2 η_1= η_2. Therefore, the solution must be unique, which proves that G0−1=(A+NK)−1G_0^-1=(A+N_K)^-1 is everywhere single-valued. ∎ Remark 19. According to Dontchev and Rockafellar (Dontchev and Rockafellar, 2009, Chapter 2E), consider the generalized equation G0=A+NKG_0=A+N_K, where A is a linear operator and NKN_K is the normal cone mapping of a closed convex cone K. If the inverse mapping G0−1G_0^-1 is everywhere single-valued, then G0−1G_0^-1 is (globally) Lipschitz continuous. Lemma 20 (Relation between G−1G^-1 and G0−1G_0^-1). Then for all Δη η sufficiently close to 0, the inverse mappings G−1G^-1 and G0−1G_0^-1 satisfy the following relationship: G−1(γ)=G0−1(γ)+η¯,γnear 0.G^-1(γ)=G_0^-1(γ)+ η, γ\ near 0. (48) Proof. Consider any Δη η close to 0. We first expand the generalized equation G(η¯+Δη)G( η+ η) as G(η¯+Δη) G( η+ η) =f(ε¯,η¯)+∇ηf(ε¯,η¯)Δη+NE(η¯+Δη) =f( , η)+ _ηf( , η)\, η+N_E( η+ η) =f(ε¯,η¯)−∇ηf(ε¯,η¯)η¯+∇ηf(ε¯,η¯)(η¯+Δη)+NE(η¯+Δη). =f( , η)- _ηf( , η)\, η+ _ηf( , η)( η+ η)+N_E( η+ η). Let A:=∇ηf(ε¯,η¯)A:= _ηf( , η). Then γ∈G(η¯+Δη)γ∈ G( η+ η) if and only if γ∈f(ε¯,η¯)−Aη¯+A(η¯+Δη)+NE(η¯+Δη).γ∈ f( , η)-A η+A( η+ η)+N_E( η+ η). (49) Rearranging the terms, this is equivalent to γ−f(ε¯,η¯)+Aη¯∈A(η¯+Δη)+NE(η¯+Δη),γ-f( , η)+A η∈ A( η+ η)+N_E( η+ η), (50) which can be further rewritten as γ−f(ε¯,η¯)∈AΔη+NE(η¯+Δη).γ-f( , η)∈ A η+N_E( η+ η). (51) Note that −f(ε¯,η¯)∈NE(η¯)-f( , η)∈ N_E( η). By the Reduction Lemma 16, this holds if and only if γ∈AΔη+NK(Δη),γ∈ A η+N_K( η), (52) which is exactly γ∈G0(Δη).γ∈ G_0( η). (53) ∎ Lemma 21. Define σ(ε):=G−1(−∇εf(ε¯,η¯)(ε−ε¯)).σ( ):=G^-1\! (- _ f( , η)\,( - ) ). If G−1G^-1 admits a single-valued Lipschitz localization around 0 for η¯ η, then σ(ε)σ( ) is a first-order approximation of S(ε)S( ) at ε¯ . Proof. Since f is strictly differentiable at (ε¯,η¯)( , η), for any small perturbations (Δε,Δη)( , η) we have f(ε¯+Δε,η¯+Δη)=f(ε¯,η¯)+∇εf(ε¯,η¯)Δε+∇ηf(ε¯,η¯)Δη+o(Δε,Δη),f( + , η+ η)=f( , η)+ _ f( , η)\, + _ηf( , η)\, η+o( , η), (54) where ‖o(Δε,Δη)‖\|o( , η)\| denotes the higher-order remainder term satisfying ‖o(Δε,Δη)‖=o(‖Δε‖+‖Δη‖)\|o( , η)\|=o(\| \|+\| η\|). By (30), S(ε¯+Δε) S( + ) =η∣f(ε¯+Δε,η)+NE(η)∋0 =\η f( + ,η)+N_E(η) 0\ (55) =η∣f(ε¯,η¯)+∇εf(ε¯,η¯)Δε+∇ηf(ε¯,η¯)(η−η¯)+o(Δε,η−η¯)+NE(η)∋0. =\η f( , η)+ _ f( , η) + _ηf( , η)(η- η)+o( ,η- η)+N_E(η) 0\. Then, η∈S(ε¯+Δε)η∈ S( + ) satisfieså that f(ε¯,η¯)+∇εf(ε¯,η¯)Δε+∇ηf(ε¯,η¯)(η−η¯)+o(Δε,η−η¯)+NE(η)∋ 0.f( , η)+ _ f( , η)\, + _ηf( , η)(η- η)+o( ,η- η)+N_E(η)\; \;0. (56) Rearranging the terms, we obtain the equivalent inclusion −∇εf(ε¯,η¯)Δε−o(Δε,η−η¯)∈f(ε¯,η¯)+∇ηf(ε¯,η¯)(η−η¯)+NE(η)=G(η),- _ f( , η)\, -o( ,η- η)\;∈\;f( , η)+ _ηf( , η)(η- η)+N_E(η)\;=\;G(η), (57) where G(η):=∇ηf(ε¯,η¯)(η−η¯)+f(ε¯,η¯)+NE(η)G(η):= _ηf( , η)\,(η- η)+f( , η)+N_E(η) denotes the linearized generalized equation in η. Therefore, any solution η∈S(ε¯+Δε)η∈ S( + ) can be written as: η=G−1(−∇εf(ε¯,η¯)Δε−o(Δε,η−η¯)),η\;=\;G^-1\! (- _ f( , η)\, -o( ,η- η) ), (58) which expresses S(ε¯+Δε)S( + ) implicitly via the inverse mapping G−1G^-1. We use κ to denote the Lipschitz constant of G−1G^-1. Since η¯=G(0) η=G(0), we have: we have ‖η−η¯‖ \|η- η\| =‖G−1(−∇εf(ε¯,η¯)Δε−o(Δε,η−η¯))−G−1(0)‖ =\|G^-1\! (- _ f( , η)\, -o( ,η- η) )-G^-1(0)\| (59) ≤k‖−∇εf(ε¯,η¯)Δε−o(Δε,η−η¯)‖. ≤ k\, \|- _ f( , η)\, -o( ,η- η) \|. Applying the triangle inequality gives ‖η−η¯‖≤k‖∇εf(ε¯,η¯)‖‖Δε‖+k‖o(Δε,η−η¯)‖.\|η- η\|≤ k\,\| _ f( , η)\|\,\| \|+k\,\|o( ,η- η)\|. (60) Since o(Δε,η−η¯)=o(‖Δε‖+‖η−η¯‖)o( ,η- η)=o(\| \|+\|η- η\|) as (Δε,η−η¯)→(0,0)( ,η- η)→(0,0), for any δ>0δ>0 there exists a neighborhood of (0,0)(0,0) such that ‖o(Δε,η−η¯)‖≤δ(‖Δε‖+‖η−η¯‖).\|o( ,η- η)\|≤δ (\| \|+\|η- η\| ). (61) Combining (59) and (61) yields ‖η−η¯‖≤k‖∇εf(ε¯,η¯)‖‖Δε‖+kδ‖Δε‖+kδ‖η−η¯‖.\|η- η\|≤ k\,\| _ f( , η)\|\,\| \|+kδ\,\| \|+kδ\,\|η- η\|. (62) Rearranging terms gives (1−kδ)‖η−η¯‖≤k(‖∇εf(ε¯,η¯)‖+δ)‖Δε‖.(1-kδ)\,\|η- η\|≤ k (\| _ f( , η)\|+δ )\,\| \|. (63) Finally, choosing δ>0δ>0 small enough such that kδ<12kδ< 12, we obtain ‖η−η¯‖≤k 1−kδ(‖∇εf(ε¯,η¯)‖+δ)‖Δε‖=O(‖Δε‖),\|η- η\|≤ k\,1-kδ\, (\| _ f( , η)\|+δ )\,\| \|=O(\| \|), (64) Since G−1G^-1 is Lipschitz continuous, S(ε¯+Δε)−σ(ε¯+Δε) S( + )-σ( + ) =‖G−1(−∇εf(ε¯,η¯)Δε−o(Δε,η−η¯))−G−1(−∇εf(ε¯,η¯)Δε)‖ = \|G^-1(- _ f( , η) -o( ,η- η))-G^-1(- _ f( , η) ) \| (65) ≤k‖o(Δε,η−η¯)‖. ≤ k\|o( ,η- η)\|. Substituting (64) to (66) yields: S(ε¯+Δε)−σ(ε¯+Δε) S( + )-σ( + ) ≤k‖o(Δε,η−η¯)‖=o(‖Δε‖). ≤ k\|o( ,η- η)\|=o(\| \|). (66) This proves that σ(ε)σ( ) is a first-order approximation of S(ε)S( ) at ε¯ . ∎ B.2 Proof of Theorem 8 Following Lemma 18 and Remark 19, we establish that the mapping G0−1G_0^-1 is locally Lipschitz continuous. Moreover, by Lemma 20, it holds that G−1(v)=G0−1(v)+η¯G^-1(v)=G_0^-1(v)+ η which implies that G−1G^-1 inherits the Lipschitz continuity and single-valuedness of G0−1G_0^-1. In this case, the Lipschitz continuity of G−1G^-1 ensures the applicability of Lemma 21, which shows that σ(ε)σ( ) serves as a first-order local approximation of the solution mapping S(ε)S( ).. Indeed, for small Δε , σ(ε¯+Δε¯) σ( + ) =G−1(−∇εf(ε¯,η¯)Δε¯) =G^-1 (- _ f( , η)\, ) (67) =G0−1(−∇εf(ε¯,η¯)Δε¯)+η¯. =G_0^-1 (- _ f( , η)\, )+ η. (68) Therefore, S(ε¯+Δε¯)−S(ε¯) S( + )-S( ) =σ(ε¯+Δε¯)−η¯+o(‖Δε‖) =σ( + )- η+o(\| \|) (69) =G0−1(−∇εf(ε¯,η¯)Δε¯)+o(‖Δε¯‖) =G_0^-1\! (- _ f( , η)\, )+o(\| \|) =s¯(−∇εf(ε¯,η¯)Δε¯)+o(‖Δε¯‖). = s (- _ f( , η)\, )+o(\| \|). (70) We obtain limt↓0S(ε¯+tΔε¯)−S(ε¯)t _t 0 S( +t )-S( )t =limt↓0s¯(−∇εf(ε¯,η¯)tΔε¯)+o(‖Δε¯‖)t = _t 0 s (- _ f( , η)\,t )+o(\| \|)t (71) Since NKN_K is the normal cone mapping of the convex cone K, it satisfies NK(αw)=αNK(w)N_K(α w)=α N_K(w) for all α>0α>0. Together with the linearity of A, this gives (A+NK)(αw)=α(A+NK)(w) (A+N_K )(α w)=α (A+N_K )(w). Hence s¯=(A+NK)−1 s= (A+N_K )^-1 is positively homogeneous, limt↓0s¯(−∇εf(ε¯,η¯)tΔε¯)+o(‖Δε¯‖)t=s¯(−∇εf(ε¯,η¯)Δε¯). _t 0 s (- _ f( , η)\,t )+o(\| \|)t= s(- _ f( , η)\, ). (73) Hence, S is directionally differentiable at ε¯ with DS(ε¯)(Δε¯)=s¯(−∇εf(ε¯,η¯)Δε¯).DS( )( )= s (- _ f( , η)\, ). Equation (36) ensures that s¯(−∇εf(ε¯,η¯)Δε¯) s (- _ f( , η) ) is the solution of the auxiliary VI (18). This implies that, under the assumptions of Theorem 8, the solution mapping S is directionally differentiable at ε¯ , and its directional derivative is given by the solution of the auxiliary VI (18). Since S(ε)=(θ,λ)S( )=(θ,λ), the DIF Dθ(ε¯;Δε¯)Dθ( ; ) defined in Definition 2 corresponds to the θ-component ofDS(ε¯)(Δε¯).DS( )( ). This guarantees the existence of the DIF. B.3 Proof of Proposition 5 Proof. Given (ε¯,θ¯,λ¯)( , θ, λ) satisfy the optimality condition (11), we have ∇θL(ε¯,θ¯,λ¯)+Nℝd(θ¯)∋0,−∇λL(ε¯,θ¯,λ¯)+Nℝ+m(λ¯)∋0. _θL( , θ, λ)+N_R^d( θ) 0, - _λL( , θ, λ)+N_R^m_+( λ) 0. (74) Following the definition of the normal cone (Def. 4), condition (74) implies that: θ∈ℝd,∇θL(ε¯,θ¯,λ¯)=0,θ ^d, _θL( , θ, λ)=0, (75) λ∈ℝ+m,−∇λL(ε¯,θ¯,λ¯)≤0.λ ^m_+,- _λL( , θ, λ)≤ 0. (76) Note that R+mR^m_+ is a closed, convex cone. Following condition (24), we have ∇λL(ε¯,θ¯,λ¯)⋅λ¯=0, _λL( , θ, λ)· λ=0, (77) which is consistent with the complementary slackness condition in the KKT system. Now we derive the equivalent of VI (16) and VI (18). We assume the Δε¯,Δθ^,Δλ , θ, λ are sufficiently close to 0. VI (16) requires that ∇θL(ε¯,θ¯,λ¯)+∇θεL(ε¯,θ¯,λ¯)Δε¯+∇θ2L(ε¯,θ¯,λ¯)Δθ^+∇θλL(ε¯,θ¯,λ¯)Δλ^+Nℝd(θ¯+Δθ^)∋0 _θL( , θ, λ)+ _θ L( , θ, λ) + _θ^2L( , θ, λ) θ+ _θλL( , θ, λ) λ+N_R^d( θ+ θ) 0 (78) ⇕ θ¯+Δθ^∈ℝd,∇θL(ε¯,θ¯,λ¯)+∇θεL(ε¯,θ¯,λ¯)Δε¯+∇θ2L(ε¯,θ¯,λ¯)Δθ^+∇θλL(ε¯,θ¯,λ¯)Δλ^=0 θ+ θ ^d, _θL( , θ, λ)+ _θ L( , θ, λ) + _θ^2L( , θ, λ) θ+ _θλL( , θ, λ) λ=0 and −∇λL(ε¯,θ¯,λ¯)−∇λεL(ε¯,θ¯,λ¯)Δε¯−∇λθL(ε¯,θ¯,λ¯)Δθ^−∇λ2L(ε¯,θ¯,λ¯)Δλ^+Nℝ+m(λ¯+Δλ^)∋0 - _λL( , θ, λ)- _λ L( , θ, λ) - _λθL( , θ, λ) θ- _λ^2L( , θ, λ) λ+N_R^m_+( λ+ λ) 0 (79) ⇕ λ¯+Δλ^∈ℝ+m,∇λL(ε¯,θ¯,λ¯)+∇λεL(ε¯,θ¯,λ¯)Δε¯+∇λθL(ε¯,θ¯,λ¯)Δθ^+∇λ2L(ε¯,θ¯,λ¯)Δλ^≤0 λ+ λ ^m_+, _λL( , θ, λ)+ _λ L( , θ, λ) + _λθL( , θ, λ) θ+ _λ^2L( , θ, λ) λ≤ 0 It is straightforward to verify that ∇λ2L(ε¯,θ¯,λ¯)=0 _λ^2L( , θ, λ)=0. By equation (75), ∇θL(ε¯,θ¯,λ¯)=0 _θL( , θ, λ)=0. Substituting ∇θL(ε¯,θ¯,λ¯)=0 _θL( , θ, λ)=0 to (78) and ∇λ2L(ε¯,θ¯,λ¯)=0 _λ^2L( , θ, λ)=0 to (79) yields θ¯+Δθ^∈ℝd,∇θεL(ε¯,θ¯,λ¯)Δε¯+∇θ2L(ε¯,θ¯,λ¯)Δθ^+∇θλL(ε¯,θ¯,λ¯)Δλ^=0, θ+ θ ^d, _θ L( , θ, λ) + _θ^2L( , θ, λ) θ+ _θλL( , θ, λ) λ=0, (80) and λ¯+Δλ^∈ℝ+m,∇λL(ε¯,θ¯,λ¯)+∇λεL(ε¯,θ¯,λ¯)Δε¯+∇λθL(ε¯,θ¯,λ¯)Δθ^≤0 λ+ λ ^m_+, _λL( , θ, λ)+ _λ L( , θ, λ) + _λθL( , θ, λ) θ≤ 0 (81) Since ℝ+mR_+^m is a closed convex cone, by (24), we have (∇λL(ε¯,θ¯,λ¯)+∇λεL(ε¯,θ¯,λ¯)Δε¯+∇λθL(ε¯,θ¯,λ¯)Δθ^)⋅(λ¯+Δλ^)=0 ( _λL( , θ, λ)+ _λ L( , θ, λ)\, + _λθL( , θ, λ)\, θ )·( λ+ λ)=0 (82) In summary, VI (16) is equivalent to the system consisting of (80), (81), and (82). Note that λ=[λ1,…,λj]λ=[ _1,…, _j], we now further discuss the formulation (81) and (82) for j∈IInactive,j∈IBindingj∈ I_Inactive,j∈ I_Binding, and j∈INon-bindingj∈ I_Non-binding. Case 1. If j∈IInactivej∈ I_Inactive, then ∇λjL(ε¯,θ¯,λ¯)<0 _ _jL( , θ, λ)<0 and λ¯j=0 λ_j=0. Since Δε¯ , Δθ θ, and Δλ λ are all sufficiently close to 0 and ∇λjL(ε¯,θ¯,λ¯)<0 _ _jL( , θ, λ)<0, condition (81) ∇λjL(ε¯,θ¯,λ¯)+∇λjεL(ε¯,θ¯,λ¯)Δε¯+∇λjθL(ε¯,θ¯,λ¯)Δθ^≤0 _ _jL( , θ, λ)+ _ _j L( , θ, λ) + _ _jθL( , θ, λ) θ≤ 0 is automatically satisfied. Therefore, the term ∇λjεL(ε¯,θ¯,λ¯)Δε¯+∇λjθL(ε¯,θ¯,λ¯)Δθ^is free. _ _j L( , θ, λ)\, + _ _jθL( , θ, λ)\, θ\;\;is free. (83) Moreover, condition (82) implies that Δλ^j=0for j∈IInactive. λ_j=0 j∈ I_Inactive. (84) Case 2. If j∈INon-bindingj∈ I_Non-binding, then ∇λjL(ε¯,θ¯,λ¯)=0 _ _jL( , θ, λ)=0 and λ¯j=0 λ_j=0. Substituting ∇λjL(ε¯,θ¯,λ¯)=0 _ _jL( , θ, λ)=0 and λ¯j=0 λ_j=0 into condition (81) yields Δλ^j λ_j ≥0 ≥ 0 for j∈INon-binding, j∈ I_Non-binding, (85) ∇λjεL(ε¯,θ¯,λ¯)Δε¯+∇λjθL(ε¯,θ¯,λ¯)Δθ _ _j L( , θ, λ)\, + _ _jθL( , θ, λ)\, θ ≤0 ≤ 0 for j∈INon-binding. j∈ I_Non-binding. Case 3. If j∈IBindingj∈ I_Binding, then ∇λjL(ε¯,θ¯,λ¯)=0 _ _jL( , θ, λ)=0 and λ¯j≥0 λ_j≥ 0. Since Δλ λ is close to 0, λ¯+Δλ^∈ℝ+m λ+ λ _+^m is automatically satisfied. To satisfy conditions (81) and (82), we must have ∇λjεL(ε¯,θ¯,λ¯)Δε¯+∇λjθL(ε¯,θ¯,λ¯)Δθ^=0. _ _j L( , θ, λ)\, + _ _jθL( , θ, λ)\, θ=0. In summary, we obtain Δλ^jis free. λ_j\;\;is free. for j∈IBinding, j∈ I_Binding, (86) ∇λjεL(ε¯,θ¯,λ¯)Δε¯+∇λjθL(ε¯,θ¯,λ¯)Δθ^=0 _ _j L( , θ, λ)\, + _ _jθL( , θ, λ)\, θ=0 for j∈IBinding. j∈ I_Binding. Taken together, Cases 1–3 show that VI (16) is equivalent to the system consisting of (80) and (83–86).On the other hand, the following argument shows that VI (18) is also equivalent to this system. The VI (18) implies that: ∇θεL(ε¯,θ¯,λ¯)Δε¯+∇θ2L(ε¯,θ¯,λ¯)Δθ^+∇θλL(ε¯,θ¯,λ¯)Δλ^+Nℝd(Δθ^)∋0 _θ L( , θ, λ) + _θ^2L( , θ, λ) θ+ _θλL( , θ, λ) λ+N_R^d( θ) 0 (87) ⇕ Δθ^∈ℝd,∇θεL(ε¯,θ¯,λ¯)Δε¯+∇θ2L(ε¯,θ¯,λ¯)Δθ^+∇θλL(ε¯,θ¯,λ¯)Δλ^=0 θ ^d, _θ L( , θ, λ) + _θ^2L( , θ, λ) θ+ _θλL( , θ, λ) λ=0 and −∇λεL(ε¯,θ¯,λ¯)Δε¯−∇λθL(ε¯,θ¯,λ¯)Δθ^−∇λ2L(ε¯,θ¯,λ¯)Δλ^+ND(Δλ^)∋0 - _λ L( , θ, λ) - _λθL( , θ, λ) θ- _λ^2L( , θ, λ) λ+N_D( λ) 0 (88) ⇕ Δλ^∈D,(−∇λεL(ε¯,θ¯,λ¯)Δε¯−∇λθL(ε¯,θ¯,λ¯)Δθ^−∇λ2L(ε¯,θ¯,λ¯)Δλ^)(Δλ^′−λ^)≤0,∀Δλ^′∈D λ∈ D, (- _λ L( , θ, λ) - _λθL( , θ, λ) θ- _λ^2L( , θ, λ) λ)( λ - λ)≤ 0,∀ λ ∈ D Since D is a closed convex cone, by (24), we have (∇λεL(ε¯,θ¯,λ¯)Δε¯+∇λθL(ε¯,θ¯,λ¯)Δθ^)⋅Δλ^=0 ( _λ L( , θ, λ)\, + _λθL( , θ, λ)\, θ )· λ=0 (89) Case 1. If j∈IInactivej∈ I_Inactive, then by the definition of the space D (see Proposition 5), we have Δλ^j=0,∀Δλ^∈D λ_j=0\;\;,∀\, λ∈ D (90) Thus, condition (88) is automatically satisfied. Therefore, the term ∇λjεL(ε¯,θ¯,λ¯)Δε¯+∇λjθL(ε¯,θ¯,λ¯)Δθ^is free. _ _j L( , θ, λ)\, + _ _jθL( , θ, λ)\, θ\;\;is free. (91) Case 2. If j∈INon-bindingj∈ I_Non-binding, then by the definition of the space D (see Proposition 5), we have Δλ^j≥0,∀Δλ^∈D λ_j≥ 0,\;\;∀\, λ∈ D (92) By the condition (89), ∇λjεL(ε¯,θ¯,λ¯)Δε¯+∇λjθL(ε¯,θ¯,λ¯)Δθ^≤0 for j∈INon-binding . _ _j L( , θ, λ) + _ _jθL( , θ, λ) θ≤ 0 for j∈ I_Non-binding . (93) Case 3. If j∈Ibindingj∈ I_binding, then by the definition of the space D (see Proposition 5), we have Δλ^jis free λ_j free (94) By the condition (89), ∇λjεL(ε¯,θ¯,λ¯)Δε¯+∇λjθL(ε¯,θ¯,λ¯)Δθ^=0 for j∈IBinding _ _j L( , θ, λ) + _ _jθL( , θ, λ) θ=0 for j∈ I_Binding (95) Since (87,90–95) and (80,83–86) are equivalent formulations, VI (18) is equivalent to VI (16). ∎ B.4 Proof of corollary 3 Proof. Recall the directional derivative Dθ(ε¯;v):=limt↓0,v′→vθ(ε¯+tv′)−θ(ε¯)t,Dθ( ;v):= _t 0,v → v θ ( +tv )-θ( )t, whenever the limit exists and is finite. Let α≥0α≥ 0. - If α=0α=0 : by the definition with v≡0v≡ 0, Dθ(ε¯;0)=limt↓0,v′→0θ(ε¯+tv′)−θ(ε¯)t=0,Dθ( ;0)= _t 0,v → 0 θ ( +tv )-θ( )t=0, hence 0⋅Dθ(ε¯;v)=Dθ(ε¯;0)0· Dθ( ;v)=Dθ( ;0). - If α>0α>0 : using the definition with the direction αvα v, Dθ(ε¯;αv)=limt↓0,u′→αvθ(ε¯+tu′)−θ(ε¯)t.Dθ( ;α v)= _t 0,u →α v θ ( +tu )-θ( )t. Choose u′=αv′u =α v with v′→v → v. Then Dθ(ε¯;αv) Dθ( ;α v) =limt↓0,v′→vθ(ε¯+t(αv′))−θ(ε¯)t=limt↓0,v′→vθ(ε¯+(αt)v′)−θ(ε¯)t = _t 0,v → v θ ( +t (α v ) )-θ( )t= _t 0,v → v θ ( +(α t)v )-θ( )t =lims↓0,v′→vθ(ε¯+sv′)−θ(ε¯)s/α( set s=αt) = _s 0,v → v θ ( +sv )-θ( )s/α ( set s=α t) =αlims↓0,v′→vθ(ε¯+sv′)−θ(ε¯)s=αDθ(ε¯;v). =α _s 0,v → v θ ( +sv )-θ( )s=α Dθ( ;v). Combining the two cases yields αDθ(ε¯;Δε¯)=Dθ(ε¯;αΔε¯),∀α∈ℝ+α Dθ( ; )=Dθ( ;α ), ∀α ^+ ∎ B.5 Proof of proposition 6 Section B.2 has proved that the directional derivative of the solution mapping S is characterized by the auxiliary VI (18). Specifically, DS(ε¯)(Δε¯)=Δη^whereμΔε¯+AΔη^+NK(Δη^)∋0.DS( )( )= η μ\, +A\, η+N_K( η) 0. (96) Since the DIF Dθ(ε¯;Δε¯)Dθ( ; ) defined in Definition 2 is the θ-component of DS(ε¯)(Δε¯)DS( )( ), we have Δθ^=Dθ(ε¯;Δε¯). θ\;=\;Dθ( ; ). (97) B.6 Proof of Theorem 10 Define the Lagrangian of QP (21) as ℒQP(ω,ζ)= _QP(ω,ζ)= L(ε¯,θ¯,λ¯)+⟨∇θεL(ε¯,θ¯,λ¯)⋅Δε¯,ω⟩+12⟨ω,∇θθ2L(ε¯,θ¯,λ¯)ω⟩ \;L( , θ, λ)+ _θ L( , θ, λ)\!·\! ,\;ω + 12 ω,\; _θ^2L( , θ, λ)\,ω (98) +∑j∈Ibindingζj[1Nj∑i=1Njℓj(zi(j),θ¯)+1Nj∑i=1Nj∇θℓj(zi(j),θ¯)⋅ω+∑zi(j)∈ZrΔε¯jℓj(zi(j),θ¯)−τj] \;+ _j∈ I_binding _j\! [ 1N_j _i=1^N_j _j(z_i^(j), θ)+ 1N_j _i=1^N_j _θ _j(z_i^(j), θ)\!·\!ω+ _z_i^(j)∈ Z^r _j\, _j(z_i^(j), θ)- _j ] +∑j∈Inon-bindingζj[1Nj∑i=1Njℓj(zi(j),θ¯)+1Nj∑i=1Nj∇θℓj(zi(j),θ¯)⋅ω+∑zi(j)∈ZrΔε¯jℓj(zi(j),θ¯)−τj] \;+ _j∈ I_non-binding _j\! [ 1N_j _i=1^N_j _j(z_i^(j), θ)+ 1N_j _i=1^N_j _θ _j(z_i^(j), θ)\!·\!ω+ _z_i^(j)∈ Z^r _j\, _j(z_i^(j), θ)- _j ] +∑j∈IInactiveζj[1Nj∑i=1Njℓj(zi(j),θ¯)+1Nj∑i=1Nj∇θℓj(zi(j),θ¯)⋅ω+∑zi(j)∈ZrΔε¯jℓj(zi(j),θ¯)−τj]. \;+ _j∈ I_Inactive _j\! [ 1N_j _i=1^N_j _j(z_i^(j), θ)+ 1N_j _i=1^N_j _θ _j(z_i^(j), θ)\!·\!ω+ _z_i^(j)∈ Z^r _j\, _j(z_i^(j), θ)- _j ]. A pair (ω⋆,ζ⋆)(ω ,ζ ) satisfies the Karush–Kuhn–Tucker (KKT) conditions if: (i) Primal feasibility. 1Nj∑i=1Njℓj(zi(j),θ¯)+1Nj∑i=1Nj∇θℓj(zi(j),θ¯)⋅ω⋆+∑zi(j)∈ZrΔε¯jℓj(zi(j),θ¯)−τj=0, 1N_j _i=1^N_j _j(z_i^(j), θ)+ 1N_j _i=1^N_j _θ _j(z_i^(j), θ)\!·\!ω + _z_i^(j)∈ Z^r _j\, _j(z_i^(j), θ)- _j=0, ∀j∈Ibinding, ∀ j∈ I_binding, (99) 1Nj∑i=1Njℓj(zi(j),θ¯)+1Nj∑i=1Nj∇θℓj(zi(j),θ¯)⋅ω⋆+∑zi(j)∈ZrΔε¯jℓj(zi(j),θ¯)−τj≤0, 1N_j _i=1^N_j _j(z_i^(j), θ)+ 1N_j _i=1^N_j _θ _j(z_i^(j), θ)\!·\!ω + _z_i^(j)∈ Z^r _j\, _j(z_i^(j), θ)- _j≤ 0, ∀j∈Inon-binding. ∀ j∈ I_non-binding. (100) For j∈IInactivej∈ I_Inactive, there is no constraint. (i) Stationarity. ∇θεL(ε¯,θ¯,λ¯)⋅Δε¯+∇θθ2L(ε¯,θ¯,λ¯)ω⋆+∑j∈Ibinding∪Inon-binding∪IInactiveζj⋆(1Nj∑i=1Nj∇θℓj(zi(j),θ¯))=0. _θ L( , θ, λ)\!·\! + _θ^2L( , θ, λ)\,ω + _j∈ I_binding∪ I_non-binding∪ I_Inactive _j \! ( 1N_j _i=1^N_j _θ _j(z_i^(j), θ) )=0. (101) (i) Dual feasibility. ζj⋆≥0,∀j∈Inon-binding. _j ≥ 0, ∀ j∈ I_non-binding. (102) (iv) Complementary slackness. ζj⋆[1Nj∑i=1Njℓj(zi(j),θ¯)+1Nj∑i=1Nj∇θℓj(zi(j),θ¯)⋅ω⋆+∑zi(j)∈ZrΔε¯jℓj(zi(j),θ¯)−τj]=0,∀j∈Inon-binding. _j \! [ 1N_j _i=1^N_j _j(z_i^(j), θ)+ 1N_j _i=1^N_j _θ _j(z_i^(j), θ)\!·\!ω + _z_i^(j)∈ Z^r _j\, _j(z_i^(j), θ)- _j ]=0, ∀ j∈ I_non-binding. (103) We now show that any (ω⋆,ζ⋆)(ω ,ζ ) satisfying the above KKT conditions also satisfies the auxiliary VI (18). Recall that in Section B.3, we have proved that the auxiliary VI (18) is equivalent to the system (87,90−95)( 6d, 6e-- 6f). Note that L(ε,θ,λ) L( ,θ,λ) =1N0∑i=1N0ℓ0(zi(0),θ)+ε0∑zi(0)∈Zrℓ0(zi(0),θ) = 1N_0 _i=1^N_0 _0(z_i^(0),θ)+ _0\!\! _z_i^(0)∈ Z^r\! _0(z_i^(0),θ) +∑j=1mλj[1Nj∑i=1Njℓj(zi(j),θ)+εj∑zi(j)∈Zrℓj(zi(j),θ)−τj], + _j=1^m _j\! [ 1N_j _i=1^N_j _j(z_i^(j),θ)+ _j\!\! _z_i^(j)∈ Z^r\! _j(z_i^(j),θ)- _j ], we have 1Nj∑i=1Nj∇θℓj(zi(j),θ¯)=∇θλjL(ε¯,θ¯,λ¯). 1N_j _i=1^N_j _θ _j(z_i^(j), θ)= _θ _jL( , θ, λ). (104) The stationarity condition can be rewritten as ∇θεL(ε¯,θ¯,λ¯)Δε¯+∇θ2L(ε¯,θ¯,λ¯)ω⋆+∇θλL(ε¯,θ¯,λ¯)ζ⋆=0. _θ L( , θ, λ)\, + _θ^2L( , θ, λ)\,ω + _θλL( , θ, λ)\,ζ =0. (105) By the definition of the normal cone (Definition 4), this is equivalent to ∇θεL(ε¯,θ¯,λ¯)Δε¯+∇θ2L(ε¯,θ¯,λ¯)ω⋆+∇θλL(ε¯,θ¯,λ¯)ζ⋆∈Nℝd(Δθ^). _θ L( , θ, λ)\, + _θ^2L( , θ, λ)\,ω + _θλL( , θ, λ)\,ζ ∈ N_R^d( θ). (106) We distinguish the three index sets. Case 1 (j∈IInactivej∈ I_Inactive). By complementary slackness, ζj⋆=0 _j =0. Case 2 (j∈Inon-bindingj∈ I_non-binding). By complementary slackness, ζj⋆≥0 _j ≥ 0. Moreover, since 1Nj∑i=1Njℓj(zi(j),θ¯)=τj 1N_j _i=1^N_j _j(z_i^(j), θ)= _j, substituting ∇λjεL(ε¯,θ¯,λ¯)Δε¯=∑zi(j)∈ZrΔε¯jℓj(zi(j),θ¯) _ _j L( , θ, λ)\, = _z_i^(j)∈ Z^r _j _j(z_i^(j), θ) and ∇λjθL(ε¯,θ¯,λ¯)ω⋆=1Nj∑i=1Nj∇θℓj(zi(j),θ¯)ω⋆ _ _jθL( , θ, λ)\,ω = 1N_j _i=1^N_j _θ _j(z_i^(j), θ)ω into the primal feasibility condition gives ∇λjεL(ε¯,θ¯,λ¯)Δε¯+∇λjθL(ε¯,θ¯,λ¯)ω∗≤0. _ _j L( , θ, λ)\, + _ _jθL( , θ, λ)\,ω^*≤ 0. (107) Case 3 (j∈Ibindingj∈ I_binding). Here ζj⋆ _j is free. Since 1Nj∑i=1Njℓj(zi(j),θ¯)=τj 1N_j _i=1^N_j _j(z_i^(j), θ)= _j, using the same substitutions into of case 3 into primal feasibility as above yields ∇λjεL(ε¯,θ¯,λ¯)Δε¯+∇λjθL(ε¯,θ¯,λ¯)ω∗=0. _ _j L( , θ, λ)\, + _ _jθL( , θ, λ)\,ω^*=0. (108) Combining all cases, ζ⋆ζ satisfies ∇λεL(ε¯,θ¯,λ¯)Δε¯−∇λθL(ε¯,θ¯,λ¯)ω⋆−∇λ2L(ε¯,θ¯,λ¯)ζ⋆+ND′(ζ⋆)∋0, _λ L( , θ, λ)\, - _λθL( , θ, λ)\,ω - _λ^2L( , θ, λ)\,ζ +N_D (ζ ) 0, (109) where D′:=ζ∗∈ℝm|ζj∗≥0(j∈Inon-binding),ζj∗=0(j∈IInactive).D := \ζ^* ^m\; |\;ζ^*_j≥ 0\ (j∈ I_non-binding),\ ζ^*_j=0\ (j∈ I_Inactive) \. (109) together with (105) shows that (ω⋆,ζ⋆)(ω ,ζ ) satisfies 0∈A(ω⋆,ζ⋆)+μΔε¯+NK(ω⋆,ζ⋆),0∈ A(ω ,ζ )+μ +N_K(ω ,ζ ), which is exactly the auxiliary VI (18). ∎ Appendix C Details in Toy example We consider the toy example (3) minθ∈ℝ213(θ⊤x1−y1)2+13(θ⊤x2−y2)2+13(θ⊤x3−y3)2+ε(θ⊤x2−y2)2s.t.‖θ‖1≤1, _θ ^2\; 13(θ x_1-y_1)^2+ 13(θ x_2-y_2)^2+ 13(θ x_3-y_3)^2+ \,(θ x_2-y_2)^2 .t. \|θ\|_1≤ 1, with data x1=(1,0),x2=(1,0),x3=(0,1)x_1=(1,0),\,x_2=(1,0),\,x_3=(0,1) and y1=1,y2=0,y3=12y_1=1,\,y_2=0,\,y_3= 12. At ε=0 =0, the constrained solution is θ¯=(0.5, 0.5) θ=(0.5,\,0.5). The corresponding dual variable is λ¯=0 λ=0. Gradients and Hessian. For a single squared loss ℓi(θ)=(θ⊤xi−yi)2 _i(θ)=(θ x_i-y_i)^2, ∇θℓi(θ)=2(θ⊤xi−yi)xi,∇θθ2ℓi(θ)=2xixi⊤. _θ _i(θ)=2(θ x_i-y_i)x_i, ∇^2_θ _i(θ)=2x_ix_i . Hence, at ε=0 =0 the (unconstrained) Hessian of 13∑i=13ℓi 13 _i=1^3 _i is Hθ¯=13∑i=132xixi⊤=23(x1x1⊤+x2x2⊤+x3x3⊤)=(430023).H_ θ= 13\! _i=1^32x_ix_i = 23 (x_1x_1 +x_2x_2 +x_3x_3 )= pmatrix 43&0\\[2.0pt] 0& 23 pmatrix. (110) The mixed derivative of the perturbed part is ∇θεL(ε¯,θ¯)=∇θℓ2(θ¯)=2(θ¯⊤x2−y2)x2=2⋅0.5⋅(1,0)=(1,0). _θ L( , θ)= _θ _2( θ)=2( θ x_2-y_2)x_2=2· 0.5·(1,0)=(1,0). (111) C.1 Classical IF (ignoring constraints) Differentiating the stationarity condition ∇θL(ε,θ)=0 _θL( ,θ)=0 w.r.t. ε at (ε¯,θ¯)( , θ) gives Hdθdε+∇θεL(ε¯,θ¯)=0⇒dθ(ε)dε|ε=0=−Hθ¯−1∇θεL(ε¯,θ¯).H\, dθd + _θ L( , θ)=0 . dθ( )d |_ =0=-H^-1_ θ _θ L( , θ). Using (110–111), dθ(ε)dε|ε=0=−(340032)(1,0)=(−34, 0). . dθ( )d |_ =0=- pmatrix 34&0\\[2.0pt] 0& 32 pmatrix(1,0)= (- 34,\,0 ). (112) Removing sample z2z_2 corresponds to Δε=−13 =- 13, so the IF estimate is ΔθIF≈dθdεΔε=(−34, 0)⋅(−13)=(14, 0). _IF≈ dθd \, = (- 34,\,0 )\!·\! (- 13 )= ( 14,\,0 ). (113) Then θ¯+ΔθIF=(0.75, 0.5) θ+ _IF=(0.75,\,0.5) whose ℓ1 _1-norm equals 1.251.25, i.e., the IF step is infeasible. C.2 DIF (feasible, sensitivity on the linearized VI) DIF linearizes the KKT/VI system at (ε¯,θ¯)( , θ) and searches Δθ θ in the tangent subspace. The QP (21) corresponding to the toy problem (3) is minΔθ∈ℝ2⟨∇θεL(ε¯,θ¯)Δε,Δθ⟩+12Δθ⊤HΔθs.t.Δθ1+Δθ2=0. _ θ ^2\; _θ L( , θ)\, ,\; θ + 12\, θ H\, θ .t. _1+ _2=0. (114) With Δε=−13 =- 13 and b:=∇θεLΔε=(−13,0)b:= _θ L\, =(- 13,0), (114) is minΔθb⊤Δθ+12Δθ⊤HΔθs.t.Δθ1+Δθ2=0. _ θ\;b θ+ 12\, θ H\, θ .t. _1+ _2=0. The KKT conditions (with multiplier μ) are HΔθ+b+μ(1,1)⊤=0,Δθ1+Δθ2=0.H θ+b+μ(1,1) =0, _1+ _2=0. Expanding coordinates yields 43Δθ1−13+μ=0,23Δθ2+μ=0,Δθ1+Δθ2=0. 43 _1- 13+μ=0, 23 _2+μ=0, _1+ _2=0. From the second and third equations μ=−23Δθ2μ=- 23 _2 and Δθ2=−Δθ1 _2=- _1. Substituting into the first gives 43Δθ1−13−23Δθ2=0⇒6Δθ1=1 43 _1- 13- 23 _2=0 6 _1=1, hence ΔθDIF=(16,−16)\; _DIF=( 16,\,- 16)\; (115) which lies on the tangent subspace and thus preserves feasibility to first order. Appendix D Details in Constrained Linear Regression We consider per-sample ℓi(θ)=12(xi⊤θ−yi)2, _i(θ)= 12\,(x_i θ-y_i)^2, The constrained optimum θ¯ θ solves minθ _θ 1n∑i=1nℓi(θ) 1n _i=1^n _i(θ) (116) s.t. Aeqθ=beq,Aineqθ≤bineq. A_eqθ=b_eq, A_ineqθ≤ b_ineq. Assume we remove the data sample (xi,yi)(x_i,y_i). We use H=∇2L(θ¯)=1nX⊤X,gi=∇θℓi(θ¯)=(xi⊤θ¯−yi)xi.H=∇^2L( θ)= 1nX X, g_i= _θ _i( θ)=(x_i θ-y_i)\,x_i. Classical Influence Function (IF). Ignoring feasibility constraints, the first-order effect of removing sample (xi,yi)(x_i,y_i) is ΔθIF=−H−11ngi. _IF\;=\;-\,H^-1\, 1n\,g_i. (117) Penalty-based IF. We approximate feasibility via a penalized surrogate L~(θ)=12n‖Xθ−y‖22+ρ2‖Aeqθ−beq‖22+∑j=1mk((aj⊤θ−bj)+)3, L(θ)= 12n\|Xθ-y\|_2^2+ ρ2\|A_eqθ-b_eq\|_2^2+ _j=1^mk\, ((a_j θ-b_j)_+ )^3, (118) where (t)+=maxt,0(t)_+= \t,0\, ρ,k>0ρ,k>0, and aj⊤a_j is the j-th row of AineqA_ineq. The Hessian at θ¯ θ is Hpen=1nX⊤X+ρAeq⊤Aeq+∑j∈6ktjajaj⊤,tj=aj⊤θ¯−bj,=j:tj>0.H_pen= 1nX X+ρ\,A_eq A_eq+ _j 6k\,t_j\,a_ja_j , t_j=a_j θ-b_j, =\j:\ t_j>0\. (119) The corresponding update is ΔθpenIF=−Hpen−11ngi. _penIF\;=\;-\,H_pen^-1\, 1n\,g_i. (120) Directional Influence Function (DIF). We enforce feasibility with a one-step constrained QP (linearized KKT): ΔθDIF=argminΔθ∈ℝd _DIF= _ θ ^d 12Δθ⊤HΔθ−1ngi⊤Δθ 12\, θ H\, θ\;-\; 1n\,g_i θ (121) s.t. AeqΔθ=0, A_eq θ=0, AineqΔθ=0(active constraints only). A_ineq θ=0\;\;(active constraints only).