{'fraction_non_alphanumeric': 0.054336835626676044, 'fraction_numerical': 0.01754077460764057, 'mean_word_length': 4.264522101365425, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Learning Linear Complementarity Systems

This paper investigates the learning, or system identification, of a class of piecewise-affine dynamical systems known as linear complementarity systems (LCSs). We propose a violation-based loss which enables efficient learning of the LCS parameterization, without prior knowledge of the hybrid mode boundaries, using gradient-based methods. The proposed violation-based loss incorporates both dynamics prediction loss and a novel complementarity - violation loss. We show several properties attained by this loss formulation, including its differentiability, the efficient computation of first- and second-order derivatives, and its relationship to the traditional prediction loss, which strictly enforces complementarity. We apply this violation-based loss formulation to learn LCSs with tens of thousands of (potentially stiff) hybrid modes. The results demonstrate a state-of-the-art ability to identify piecewise-affine dynamics, outperforming methods which must differentiate through non-smooth linear complementarity problems.


Introduction
Many physical systems of interest are well captured by multi-modal or hybrid representations. For example, robotics problems which treat contact with the environment (Stewart and Trinkle, 2000;Brogliato, 1999), optimal control problems (Bemporad et al., 2000), and control of networks (Heemels et al., 2002), all exhibit switching or hybrid properties.
In this work, we are interested in system identification/model learning of multi-modal systems. We focus on piecewise-affine (PWA) models as they can sufficiently describe the multi-modal nature of dynamics due to the approximation properties of affine functions (Breiman, 1993;Lin and Unbehauen, 1992) but are tractable enough for control tasks due to their simple (affine) structure over polyhedral regions (Bemporad and Morari, 1999). Even though PWA models are widely used, it is well-known that PWA regression is NPhard in general (see (Lauer, 2015) for a detailed analysis), because it requires simultaneous classification of the data points into modes and the regression of a submodel for each mode.
In this paper, we consider PWA models in the context of linear complementarity systems (LCSs) (Heemels et al., 2000). We focus on a subclass of LCS models (with P-matrix assumption) that are equivalent to continuous piecewise affine models (Heemels et al., 2001;Camlibel et al., 2007). LCS models are efficient representations of PWA models and an LCS has the ability to represent/approximate a system with large number of hybrid modes compactly, with only few complementarity variables. In some cases, an LCS with n λ complementarity variables is equivalent to a PWA model with 2 n λ modes. Many robotics problems that involve contact can be efficiently locally approximated via LCS models, e.g., we have exploited the LCS representation to enable contact-aware  and realtime control of robotic tasks . In this work, we propose an approach that learns an LCS from state-input data of a hybrid system, which does not contain any prespecified number of modes. The approach is able to identify LCS models by proposing an implicit loss function.

Related Work
Many successful approaches in identifying piecewise models have been proposed over the years. See the survey paper (Paoletti et al., 2007) for a detailed overview. Mixed integer formulations that mainly focus on hinging hyperplanes and piecewise affine Wiener models have been proposed (Roll et al., 2004) but as the number of integer variables scale with number of data points such approaches are only applicable in small data regime. On the contrary, researchers have also focused on convex formulations where first they estimate a set of submodels and then select few of them that explains the data (Elhamifar et al., 2014), but the approach relies on restricting the parameter space and can be overly conservative.
Many alternate approaches that enable PWA system identification from data exist such as (Ferrari-Trecate et al., 2003;Nakada et al., 2005;Bemporad et al., 2005;Hartmann et al., 2015;Du et al., 2020). Researchers also suggested recursive PWA identification algorithms (Bako et al., 2011), (Breschi et al., 2016). Most of the above methods are clustering-based where a predetermined number of models are identified and each training data point is associated with one of the models. Then, linear separation techniques are used to compute the polyhedral partitions. This iterative nature of the methods can lead to overly conservative, suboptimal solutions. Approaches that simultaneously cluster, PWL-separate and fit (Bemporad, 2021) rely heavily on initial assignment of data points to clusters. Unlike our approach, none of the methods above have been tested on identifying PWA functions with thousands of partitions, and most of them have been only tested on functions with less than 30 pieces.
For more expressive models such as deep neural networks, researchers have explored the positive effect of imposing structured knowledge to capture the multi-modality (de Avila Belbute-Peres et al., 2018;Li et al., 2018;Battaglia et al., 2016). Particularly in robotics, special emphasis has been on multi-body systems with frictional contact (Geilinger et al., 2020) and it has been shown that imposing structure leads to accurate, sample efficient strategies (Pfrommer et al., 2020). Similar work has demonstrated the difficulty inherent in learning non-smooth dynamical systems without exploiting particular structures (Parmar et al., 2021). Researchers have also explored learning models as functionals of signeddistance fields (Driess et al., 2021). These methods lead to rich, accurate but complex models that are not amenable to techniques of model-based control. On the contrary, here our focus is on simple models such as PWA models that enable model-based control while sufficiently capturing the hybrid dynamics.
Notation: In this paper, regular and bold lowercase letters represent scalar and vectors, respectively. Uppercase letters represent matrices. For vector v ∈ R n , v[i] is the i-th entry, i = 1, 2, ..., n. diag(v) is to diagonalize a vector v into a matrix. vec (A) denotes the vectorlization of a matrix A into a column vector; ⊗ is the Kronecker product. I n denotes the identity matrix with size of n × n. A 0 means symmetric A is positive definite.

Problem Statement
Consider the following discrete-time linear complementarity system (LCS), where the state evolution is governed by a linear dynamics in (1a) and a linear complementarity problem (LCP) in (1b): Here, x t ∈ R nx and u t ∈ R nu are the system state and input at time step t, respectively; and λ t ∈ R r is the complementarity variable at time step t. A ∈ R nx×nx , B ∈ R nx×nu , C ∈ R nx×n λ , d ∈ R nx , D ∈ R n λ ×nx , E ∈ R n λ ×nu , F ∈ R n λ ×n λ , and c ∈ R n λ are system matrix/vector parameters. At (x t , u t ), λ t is solved from the LCP in (1b), written as It is well-known that LCP(F, q t ) has a unique solution λ t for every q t if and only if F is P -matrix (Cottle et al., 2009). We will discuss this in the next section.
In this paper, we consider to learn a LCS from a dataset . Specifically, we aim to find the system parameter by minimizing a loss L(θ, D). Thus, the problem of interest in this paper is to solve min θ∈Θ L(θ, D) + R(θ).
Here, R(θ) can be any regularization term imposed on θ which will be discussed later.

LCP and Prediction-based Formulation
This section will discuss the solution to LCP in (2), and then describe a prediction-based loss formulation L(θ, D) for (4). To start, we make the following assumption on F in (2).
The set of F satisfying Assumption 1 contains all positive matrices of feasible dimension and any asymmetric matrices with definite-positive symmetric part. In robotics applications, F with Assumption 1 has been widely used in soft contact dynamics problems such as [TBD]. Any F satisfying Assumption 1 can be shown to be a P -matrix (Tsatsomeros, 2002), thus leading to the existence and uniqueness of λ t . In fact, under Assumption 1, λ t = LCP(F, q t ) can be solved by the following convex optimization due to the fact λ T F λ = 1 2 λ T (F + F T )λ, With the above assumption, one natural loss in (4) can be Here, x θ t+1 is the predicted next state, implicitly depending on θ. We call (6) predictionbased loss, as it evaluates the difference between the predicted x θ t+1 and observed x * t+1 . One can minimize (6) via any gradient-based method by differentiating through LCP (de Avila Belbute-Peres et al., 2018). This requires differentiblity of a LCP, given below.
A sketch of a proof of Lemma 1 is given in Appendix 6. The above prediction-based loss in (6) will serve as a benchmark in the following method development.

Proposed Method for Learning LCS
This section will develop a new method for learning LCS. As we will show this section and the experiments in next section, the proposed method attains several advantages over the prediction-based loss (6) both in theoretical property and implementation.

Violation-based Loss
To start, we give the following lemma stating an equivalence of a LCP.
Lemma 2 Given any q t ∈ R n λ and F satisfying Assumption 1, solving λ t = LCP(F, q t ) is the equivalent to solving the following strongly-convex quadratic program: with any constant 0 < γ < σ min (F T +F ) ( σ min (·) denotes the smallest singular value) . (7). To do so, we compute the Hessian of f (λ, φ), Due to Schur complement, . This is true because γ<σ min (F T +F ). This completes the proof.
In (7), we have introduced a proxy variable φ ≥ 0 to represent LCP constraint F λ+q t ≥ 0. Compared to other equivalences of LCP, such as (5), we emphasize the following benefits of (7) with the introduced proxy variable φ. First, (7) now only has box constraints which are independent from θ; this will facilitate the learning process because one does not need to explicitly track the active and inactive constraints and differentiate through constraints (which usually leads to numerical difficulty as shown in (Jin et al., 2021)). Second, compared to (5), (7) turns hard constraint F λ + q t ≥ 0 into a soft penalty; this may smooth the landscape of the proposed loss, facilitating the optimization process over θ. With Lemma 2, we are now in a position to propose the following loss function for learning LCS, includes two parts: the violation of dynamics (1a) and the violation of the LCP, as stated in Lemma 2. We have introduced parameter > 0 to control the weight of penalties on the two violations. (9) is to minimize data's violation to both dynamics (1a) and complementarity constraints (1b), thus we name it violation-based loss. In what follows, we will show that the violation-based loss attains some good properties both for theoretical analysis and algorithmic implementation, in comparison with the prediction-based loss (6).

Properties of Violation-based Loss
The first lemma shows that the violation-based loss (9) is a strongly-convex quadratic program w.r.t. (λ, φ) and allows much easier computation of the gradient w.r.t. θ.
(b) Let (λ ,θ t , φ ,θ t ) be the solution to (9b). L (θ, D) is differentiable with respect to θ if the strict complementarity holds for both λ t ≥0 and φ t ≥0 at (λ ,θ t , φ ,θ t ), t = 1, 2, ..., N . The gradient is given by with e dyn the strict complementarity condition holds for constrained optimization in (9b). This is a direct result from the well-known sensitivity analysis theory (see Theorem 2.1 in (Fiacco, 1976)). The gradient of L (θ, D) can be obtained directly applying the envelope theorem (Afriat, 1971). For example, the gradient of l (θ, Writing the above into the matrix form leads to ∇ A L = N t=1 e dyn t (x * t ) T . Similar derivations also apply to ∇ B L , ∇ B L , ∇ C L , ∇ D L , ∇ E L , and ∇ F L . This completes the proof.
In addition to the strongly-convex quadratic problem with bound constraints in (9b), Lemma 3 state that L (θ, D) allows for much simpler differentiation, as stated in claim (b). Note that differentiation of L (θ, D) in (10) does not involve any matrix inverse. This is in stark contrast with the prediction-based loss (6), whose differentiation (de Avila Belbute-Peres et al., 2018) is based on the implicit function theorem (Rudin et al., 1976) and requires the inverse of Jacobian matrix of KKT equations (which is computationally expensive).
Another implication of Lemma 3 is that the Lipschitz constant of L (θ, D) with respect to the LCP matrices (D, E, F, c) can be controlled by the choice of . Specifically, the second line of (10) shows that one can always choose a large to produce a small Lipschitz constant of the loss landscape with respect to (D, E, F, c). This property can facilitate the learning process by controlling the smoothness of the loss landscape, and also helpful in the generalization of learned results as analyzed in a concurrent work (Bianchini et al., 2021). However, we also need to note that the large choice of could lead to the bias learning results, as shown and analyzed in the later simulation examples.
We further have the following result, which states the second-order derivative of L (θ, D).
Proof To prove Lemma 4, we need first to show diag ∂L is invertible. We only provide the sketch of the proof due to page limits. First, the KKT conditions at z t = (λ ,θ t , φ ,θ t ) can be written as the following LCP The strict complementarity (differentiblity of L ) stated in claim (b) of Lemma 3 is equivalent to say the above LCP in (12) is strictly complementarity. By claim (a) in Lemma 3, we have known ∂ 2 L ∂z t ∂z t 0, which is a P-matrix. Following the same proof in Appendix 6, one can show that diag ∂L Now we prove Lemma 4. By applying envelop theorem (Afriat, 1971) to l (θ, x * t , u * t , x * t+1 ) in (9b), one can write ∇ θ L = ∂L ∂θ T . When taking the second-order derivative, one has Here, by differentiating through the LCP in (12), one can obtain Plugging the above to (13) leads to (11).
Lemma 4 states that first, Hessian of the violation-based loss with respect to the system parameter θ can also be analytically obtained. Such Hessian is important both for algorithmic implementation and theoretical analysis. Arithmetically, the above Hessian can be used to develop second-order methods for optimizing (4). Analytically, the Hessian can be used to analyze the convexity of the problem. Specifically, if (9b) is convex jointly with respect to (λ t , φ t ) and θ, one can show that L (θ, D) will be convex (also see Section 3.2.5 in (Boyd et al., 2004)). This holds for all other system matrices/vectors except matrices C and F , which imposes the challenges for learning process.
Finally, we give the following result showing the violation-based loss L (θ, D) in (9) can be controlled to approximate the prediction-based loss L pred (θ, D) (6) in terms of both the loss itself and its differentiability.
Lemma 5 Given F satisfying Assumption 1 and any constant 0 < γ < σ min (F T + F ), there exists ∆ > 0 such that for any ∈ (0, ∆], Proof We here only provide the sketch for the proof of the above lemma due to page limits. In the proof of claim (a), first, we can show that the strict complementarity in Lemma (1) is equivalently to say the strict complementarity for (7) in Lemma 2, i.e., the following LCP is strict complementarity. Further, one can show that (12) will converge to (15) as → 0.
Since the strict complementarity preserves as falls in a small neighborhood around 0 (this is similar to the proof of Theorem 2.1 in (Fiacco, 1976)), one can say L (θ, D) is differentiable with any small > 0 in the neighborhood around 0. The proof of claim (b) can directly follow the standard proof in penalty-based optimization (Fiacco and McCormick, 1990).
The above lemma has shown that the proposed violation-based L (θ, D) and predictionbased loss L pred (θ, D) are essentially related to each other in term of function itself and its differentiability with respect to θ. Specifically, the differentiability of prediction-based loss L pred (θ, D) in Lemma 1, i.e., strict complementarity for the LCP, always implies the differentiability of the violation-based loss L (θ, D) for any choice of small > 0. Second, by controlling → 0, the violation-based formulation approximates to prediction-based one.
In light of all properties stated above, we now summarize the advantage of the proposed violation-based loss L (θ, D) over the prediction-based loss L pred (θ, D). First, L (θ, D) allows for more efficient computation of gradient, while differentiation of L (θ, D) always requires the matrix inverse. Second, L (θ, D) permits analytical Hessian information, which is important both for the algorithmic implementation and theoretical analysis, such Hessian matrix is difficult to obtain for L pred (θ, D). Finally, the violation-based loss L (θ, D) is flexible to approximate L pred (θ, D) in terms of both the function itself and its differentiability with respect to θ, by controlling the weight parameter .

Examples
In implementation, one way to enforce Assumption 1 is using re-parameterization tricks: by re-parameterizing F = GG T + H − H T with any matrices G ∈ R n λ ×n λ and H ∈ R n λ ×n λ , one can easily see that F + F T 0. Also note that F and C in a LCS (1) are permutationand scaling-invariant with respect to D it also satisfies the following LCS withλ * t = P Sλ t .
Here P ∈ R n λ ×n λ , is any permutation matrix and S is any diagonal matrix with positive diagonal entries. To mitigate this ambiguity, we add a regularizing cost R(θ) = ω C 2 F , where · F is the matrix Frobenius norm and ω is the weighting parameter. We set ω = 10 −5 in our following experiments for both methods.
We randomly define a ground-truth LCS with θ * , where all parameters are selected from a uniform distribution in range [−1, 1]. To generate training data D = {(x * t , u * t , x * t+1 )} N train t=1 , we sample x * t and u * t from uniform distributions over [−10, 10] and [−5, 5], respectively, and then solve x * t+1 based on θ * . We also add zero-mean Gaussian noise with standard deviation σ=10 −2 to D. We generate similar, but noiseless, testing data T = {(x t ,ū t ,x t+1 )} Ntest t=1 . To evaluate the learned LCS on T , we define the following mean relative prediction error, with x θ t the state predicted by the learned LCS at (x t ,ū t ). The size of T is N test =1000. The following simulations will evaluate different aspects of the proposed violation-based learning formulation (9), in comparison with the prediction-based learning formulation (6). For any learning formulation, each evaluation includes a total of 30 training rounds (otherwise stated), and each training round uses a random ground-truth LCS of θ * to generate D and T and uses a random initialization θ to initialize the training. Each training round uses the Adam algorithm (Kingma and Ba, 2014) with the mini-batch size 200 and the learning rate 10 −3 (other Adam parameters: β 1 = 0.9, β 2 = 0.9, = 10 −6 ). Because we randomize the LCS systems, and some systems may be easier to identify than others, we expect fairly high variance of the results.

Results and Analysis
The evaluations of different aspects of the proposed violation-based formulation (9) versus the prediction-based formulation (6) are shown in Fig. 1. In summary of all evaluations, one can conclude that the proposed violation-based learning outperforms the prediction-based learning, specifically when (a) handling high numbers of system modes, e.g., 16k modes at n λ = 20 as shown in Fig. 1(a); (b) dealing with high system dimensions, e.g., n x = 128 as shown in Fig. 1(b); and (c) learning high-stiff LCS system as shown in Fig. 1(c). Further, Fig. 1(d ) and Fig. 1(d ) show that parameters γ and in the violation-based loss (9) are not sensitive to the performance, and thus finding proper γ and is not difficult in practice.   (a) In Fig. 1(a), we vary the number of complementarity constraints, i.e., n λ , with fixed n x = 10 and n u = 4. D has the size of N train = 50k for all evaluations, and σ min (F +F T ) = 1. Note that the maximum achievable number of modes depends on both n λ and n x , specifically is 2 n λ if n x = n λ . In our above evaluations, at n λ = 20, D contains around 16k modes. In the violation-based loss, γ = 10 −2 and = 10 −4 . Fig. 1(a) shows an obvious advantage of the violation-based method, especially for learning LCS with high number of modes, e.g., 16k modes. This ability of the violation-based method is also in contrast to the prior PWA work which only handles tens of modes, as discussed in introduction.
(b) In Fig. 1(b), we vary system state dimension n x with n λ = 10 and n u = 4. Here, each evaluation includes 15 training rounds, and other training settings follow (a). The results show the advantage of the proposed violation-based method over the predictionbased method in handling high dimension of system states, such as n x ≥ 64.
(c) Fig. 1(c) varies the system stiffness indicated by σ min (F +F T ), i.e., smaller σ min (F +F T ) means a stiffer system. Here, n x =8, n u =2, n λ =10, and others follow (a). The results show a signifcant advantage of the violation-based formulation (9) over prediction-based learning (6). By weighting the prediction error against the LCP violation, loss (9) generalizes better in the presence of stiffness and can be connected to graph distance (Bianchini et al., 2021).
(d-e) Fig. 1(d ) and 1(e) plot the performance of the violation-based method for different choices of γ and , respectively. Here n x =4, n λ =4, n u =2, N train =5k, and others follow (a). Fig. 1(b) confirms Lemma 2 that any 0<γ<σ min (F +F T ) does not significantly influence the results. Fig. 1(e), investigating the choice of , shows, across many orders of magnitude, the results are largely invariant to this choice. For small , following Lemma 5, we see the violation-based formulation perform similarly to the prediction-based form. Large can lead to results which grossly violate complemenatrity, and thus perform poorly. Empirically, it is not difficult to find a in the middle range which does not much influence the performance.

Conclusion
In this paper, we have proposed a violation-based loss formulation which enables to learn a LCS using gradient-based methods. The violation-based loss is a sum of dynamics prediction loss and a novel complementarity violation loss. We have shown several some properties attained by this loss formulation. The numerical results demonstrate a state-of-the-art ability to identify piecewise-affine dynamics, outperforming methods which must differentiate through non-smooth linear complementarity problems.