# Sparse approximations for contact mechanics

Kiran Sagar Kollepara, José V Aguado, Yves Le Guennec, Luisa Silva, Domenico Borzacchiello

► To cite this version:

Kiran Sagar Kollepara, José V Aguado, Yves Le Guennec, Luisa Silva, Domenico Borzacchiello. Sparse approximations for contact mechanics. 2024. hal-04615501

HAL Id: hal-04615501

https://hal.science/hal-04615501

Preprint submitted on 18 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Sparse approximations for contact mechanics

Kiran Sagar Kollepara
Luisa Silva

José V. Aguado
Domenico Borzacchiello

Yves Le Guennec

**Abstract**

Low-rank model order reduction strategies for contact mechanics show limited dimension-
ality reduction due to linear inseparability of contact pressure field. Therefore, a dictionary
based strategy is explored for creating efficient models for frictionless non-adhesive contact.
A large dictionary of contact pressure trajectories is generated using a high-fidelity finite
element model, while approximating the online query with a small number of dictionary en-
tries. This is achieved by inducing sparsity in the approximation. Accuracy, computational
effort and limitations of such methods are demonstrated on few numerical examples.

**Keywords:** contact mechanics, variational inequalities, reduced order models, sparse approximations, dictionary based approximation

# 1 Introduction

Contact mechanics deals with multi-body mechanical problems with the possibility of the outer surfaces coming to contact. The occurrence and position of contact between bodies is generally unknown a priori, which adds non-linearity to the problem. Such problems have a wide range of applications; including tribology [1], crash simulations [2], manufacturing [3] and biomechanics [4] to name a few. Important contributions to numerical treatment of contact problems have been made in the past few decades with advances in algorithms as well as the interface physics [5–9]. The complexity of the contact model can have varying degrees, with as simple as hard contact with only compressive surface traction, to as complex as frictional and adhesive contact. Also, the search for contact pairs can be a computationally intensive task.

In the context of finite element methods (FEM), major bottlenecks include solving contact detection and associated non-linearities. Contact mechanics problems are usually variational inequalities, in other words, they are inequality constrained minimization problems. In the simplest case of non-adhesive frictionless contact, the inequality constraint is essentially a non-penetration condition. Such problems involve the computation of mainly two physical unknowns: the displacement field (primal) and a non-negative contact pressure (dual). Various methods exist to resolve the mixed problem, such as the Lagrange Multiplier method (LMM), Augmented LMM, and the mortar method, among others. In the LMM, the so-called Karush-Kuhn-Tucker (KKT) conditions impose the non-negativity of contact pressure, non-penetration of contact surfaces and finally a complimentary slackness condition which forces either the contact pressure or the gap between surfaces to be null at every position on the potential contact surface [8, 10].

As the importance of real time mechanical simulations across industries grows rapidly [11, 12], it is important to address reduced order modelling (ROM) methods applied to contact mechanics. Usually, ROM methods typically consist of an *offline* and *online* stage. In the offline stage, a high-fidelity model is evaluated for different parametric instances to obtain training data, referred to as *snapshots*. The reduced model is then built using the offline data and the governing equations. In
Sparse approximations for contact mechanics

Kiran Sagar Kollepara
Luisa Silva

José V. Aguado
Domenico Borzacchiello

Yves Le Guennec

Abstract

Low-rank model order reduction strategies for contact mechanics show limited dimension-
ality reduction due to linear inseparability of contact pressure field. Therefore, a dictionary
based strategy is explored for creating efficient models for frictionless non-adhesive contact.
A large dictionary of contact pressure trajectories is generated using a high-fidelity finite
element model, while approximating the online query with a small number of dictionary en-
tries. This is achieved by inducing sparsity in the approximation. Accuracy, computational
effort and limitations of such methods are demonstrated on few numerical examples.

**Keywords:** contact mechanics, variational inequalities, reduced order models, sparse approximations, dictionary based approximation

# 1 Introduction

Contact mechanics deals with multi-body mechanical problems with the possibility of the outer surfaces coming to contact. The occurrence and position of contact between bodies is generally unknown a priori, which adds non-linearity to the problem. Such problems have a wide range of applications; including tribology [1], crash simulations [2], manufacturing [3] and biomechanics [4] to name a few. Important contributions to numerical treatment of contact problems have been made in the past few decades with advances in algorithms as well as the interface physics [5–9]. The complexity of the contact model can have varying degrees, with as simple as hard contact with only compressive surface traction, to as complex as frictional and adhesive contact. Also, the search for contact pairs can be a computationally intensive task.

In the context of finite element methods (FEM), major bottlenecks include solving contact
detection and associated non-linearities. Contact mechanics problems are usually variational
inequalities, in other words, they are inequality constrained minimization problems. In the
simplest case of non-adhesive frictionless contact, the inequality constraint is essentially a non-
penetration condition. Such problems involve the computation of mainly two physical unknowns:
the displacement field (primal) and a non-negative contact pressure (dual). Various methods
exist to resolve the mixed problem, such as the Lagrange Multiplier method (LMM), Augmented
LMM, and the mortar method, among others. In the LMM, the so-called Karush-Kuhn-Tucker
(KKT) conditions impose the non-negativity of contact pressure, non-penetration of contact
surfaces and finally a complimentary slackness condition which forces either the contact pressure
or the gap between surfaces to be null at every position on the potential contact surface [8, 10].

As the importance of real time mechanical simulations across industries grows rapidly [11, 12],
it is important to address reduced order modelling (ROM) methods applied to contact mechanics.
Usually, ROM methods typically consist of an *offline* and *online* stage. In the offline stage, a high-
fidelity model is evaluated for different parametric instances to obtain training data, referred to as
*snapshots*. The reduced model is then built using the offline data and the governing equations. In
the online phase, the user can evaluate the reduced model very efficiently, permitting exploration of the parametric space.

In an LMM-based approach to reduced modelling of contacts, two reduced basis are sought. This can be seen in one of the earliest contributions from Haasdonk et al. [13], where the reduced contact pressure space is defined as the span of contact pressure snapshots. Compression of contact pressure data was explored by Balajewicz et al. [14] and [15, 16] using Non-Negative Matrix Factorization and greedy algorithms with a non-negativity preserving projection, respectively. A reduced integration domain (RID) based approach is presented in Fauque et al. [17], Le Berre et al. [18], where the surface integrals and contact pressure degrees of freedom (dofs) are evaluated only on the surface of the RID in a reduced setup, followed by recovery of contact pressure through a non-negative least-squares interpolation from a training set.

Some works have bypassed the creation of dual reduced bases. For example, Bader et al. [19], Niakh et al. [20] use the penalty approach instead of LMM which eliminates the dual variable at the cost of allowing small penetration values. On the other hand, [21] puts forward the argument that reducing the dual dofs is unnecessary as the number of surface dofs is much less than the number of internal dofs. However, such assumptions may not always hold true in examples such as membrane-like or shell-like geometries that experience contact [13, 22].

Most works mentioned above assume the presence of a low-dimensional contact pressure space for building the reduced model. However, as shown recently in our recent work Kollepara et al. [23], contact problems have a peculiar feature that renders this assumption untrue in many cases. This feature is the contact zone, that exhibits a local nature. The contact zone is a function of the state of the system, and hence, may strongly depend on the loading, geometry, and other physical parameters.

In general, for variational inequalities, only some of the constraints are usually active. An inequality constraint is referred to as active only if the current state of the problem satisfies the equality. The so-called complementarity slack KKT condition implies that the Lagrange multipliers assume a non-zero value only where constraints are active. In other words, contact pressure is non-zero only at parts of the surface where contact is established. The strong dependence of the active zone on the state of the system leads to the so-called *linear inseparability* of contact pressure. In such a situation, contact pressure field cannot be approximated well using linear combinations of other snapshots. Thus, the snapshots can neither be compressed nor do they *generalise*<fnref n="1" /> well [23]. This has two implications from an ROM point of view:

* It becomes necessary to generate a large training set of contact pressure snapshots.

* The low-rank contact pressure basis will have approximately the same dimension as the high-fidelity solution space; thus, offering negligible reduction.

The first problem cannot be avoided, as a large database of snapshots is required for a satisfactory interpolation of the non-linear manifold. However, it is possible to mitigate the second one, the ideas for which are discussed using the visual aid of Fig. 1. The figure illustrates the high-dimensional solution space, inside which snapshots lie in a low-dimensional solution manifold. If the solution manifold is significantly non-linear, a low-rank approximation will generally not be satisfactory (Fig. 1a). To capture the manifold, a better strategy would be to generate an approximation based on snapshots nearby the queried point. This can be done in two ways:

1. Multiple localized solution subspaces will better approximate the curvature of the manifold [25], as depicted in Fig. 1b. However, to construct localized subspaces, a clustering technique is needed to group snapshots into sets based on distance.

<footnote n="1">¹Generalisation ability is a measure of model error outside the training set [24]</footnote>
the online phase, the user can evaluate the reduced model very efficiently, permitting exploration
of the parametric space.

In an LMM-based approach to reduced modelling of contacts, two reduced basis are sought. This can be seen in one of the earliest contributions from Haasdonk et al. [13], where the reduced contact pressure space is defined as the span of contact pressure snapshots. Compression of contact pressure data was explored by Balajewicz et al. [14] and [15, 16] using Non-Negative Matrix Factorization and greedy algorithms with a non-negativity preserving projection, respectively. A reduced integration domain (RID) based approach is presented in Fauque et al. [17], Le Berre et al. [18], where the surface integrals and contact pressure degrees of freedom (dofs) are evaluated only on the surface of the RID in a reduced setup, followed by recovery of contact pressure through a non-negative least-squares interpolation from a training set.

Some works have bypassed the creation of dual reduced bases. For example, Bader et al. [19], Niakh et al. [20] use the penalty approach instead of LMM which eliminates the dual variable at the cost of allowing small penetration values. On the other hand, [21] puts forward the argument that reducing the dual dofs is unnecessary as the number of surface dofs is much less than the number of internal dofs. However, such assumptions may not always hold true in examples such as membrane-like or shell-like geometries that experience contact [13, 22].

Most works mentioned above assume the presence of a low-dimensional contact pressure space for building the reduced model. However, as shown recently in our recent work Kollepara et al. [23], contact problems have a peculiar feature that renders this assumption untrue in many cases. This feature is the contact zone, that exhibits a local nature. The contact zone is a function of the state of the system, and hence, may strongly depend on the loading, geometry, and other physical parameters.

In general, for variational inequalities, only some of the constraints are usually active. An inequality constraint is referred to as active only if the current state of the problem satisfies the equality. The so-called complementarity slack KKT condition implies that the Lagrange multipliers assume a non-zero value only where constraints are active. In other words, contact pressure is non-zero only at parts of the surface where contact is established. The strong dependence of the active zone on the state of the system leads to the so-called *linear inseparability* of contact pressure. In such a situation, contact pressure field cannot be approximated well using linear combinations of other snapshots. Thus, the snapshots can neither be compressed nor do they *generalise*<fnref n="1" /> well [23]. This has two implications from an ROM point of view:

* It becomes necessary to generate a large training set of contact pressure snapshots.

* The low-rank contact pressure basis will have approximately the same dimension as the high-fidelity solution space; thus, offering negligible reduction.

The first problem cannot be avoided, as a large database of snapshots is required for a
satisfactory interpolation of the non-linear manifold. However, it is possible to mitigate the
second one, the ideas for which are discussed using the visual aid of Fig. 1. The figure illustrates
the high-dimensional solution space, inside which snapshots lie in a low-dimensional solution
manifold. If the solution manifold is significantly non-linear, a low-rank approximation will
generally not be satisfactory (Fig. 1a). To capture the manifold, a better strategy would be to
generate an approximation based on snapshots nearby the queried point. This can be done in
two ways:

1. Multiple localized solution subspaces will better approximate the curvature of the manifold [25], as depicted in Fig. 1b. However, to construct localized subspaces, a clustering technique is needed to group snapshots into sets based on distance.

<footnote n="1">¹Generalisation ability is a measure of model error outside the training set [24]</footnote>
2. Sparse regression methods can be used to select the best snapshots from the training set to render a close approximation (Fig. 1c). In other words, the sparse algorithm is allowed to select the locally significant snapshots, without the need for clustering. Indeed, this approach is equivalent to finding the best “cluster” for the queried point on-the-go. Thus, there is no need of “discretizing” the manifold into local subspaces. Also, the reconstruction quality in regions near intersection of two local subspaces need not be compromised.

3. Another possibility is to use non-linear approximation techniques [26], which is outside the scope of this article (and will be explored in a future article).

In this article, we propose to use sparse approximations from the domain of dictionary learning [27, 28]. Dictionary-based sparse approximations have already been applied to various numerical models including transport, acoustics [29] and biomedical [30] applications. In contrast to a reduced basis, an over-complete dictionary contains a large number of snapshots, some of which would essentially contain redundant information. For an approximation using linear combinations of the dictionary elements, multiple solutions would likely exist. The idea is to induce sparsity in the coefficients of the online problem, encouraging selection of a few snapshots that are the best candidates for the problem. The reader should beware that the underlying linear inseparability is not resolved by this approach. While the use of large dictionaries provide a better resolution of the non-linear manifold, the primary focus of this work is overcoming the added complexity resulting from use of large dictionaries.

Our work has some similarities as well as differences with earlier approaches. While Haasdonk et al. [13] uses snapshots to define solution subspaces, sparse regularization is not imposed for the dual reduced variable. Moreover, our approach does not attempt to compute a dual reduced basis like the works of Balajewicz et al. [14], Niakh et al. [15], Benaceur et al. [16]. Both Haasdonk et al. [13], Fauque et al. [17] tackle problems with node-to-node interaction, whereas our work includes problems with node-to-segment interactions and large “slip” between contact surfaces.

The article is structured as follows: first a very brief introduction to the discretized contact problem is given, followed by a brief introduction to sparse regression methods. Then, two proposed formulations of dictionary-based sparse methods to contact mechanics are introduced with numerical examples. Finally, some observations on the two methods are discussed and analyzed, prior to concluding the article.

## 2 Model Problem

Typical contact problems can be framed as inequality constrained minimization problems. Consider two elastic bodies with frictionless and non-adhesive surfaces, that might come into contact in a deformed state (Fig. 2). The contact problem can be framed as a minimization of the internal elastic energy while respecting the non-penetration constraints. To remain concise, the discretized version of this problem is given here. For a more detailed formulation, the reader may refer to the preceding publication [23] as well as the literature [5, 6, 8].

$$ \begin{gathered} \arg\min_{\mathbf{v}^h} \frac{1}{2} \mathbf{v}^h T \mathbf{K}(\mu) \mathbf{v}^h - \mathbf{v}^h T \mathbf{f}(\mu) \\ \mathbf{C}(\mu, \mathbf{v}^h) \mathbf{v}^h - \mathbf{g}(\mu, \mathbf{v}^h) \le \mathbf{0} \end{gathered} \tag{1} $$

where

$$ \bullet \mathbf{v}^h = \begin{bmatrix} \mathbf{v}_1^h \\ \mathbf{v}_2^h \end{bmatrix} \text{ is the stacked nodal displacement vector} $$
2. Sparse regression methods can be used to select the best snapshots from the training set to render a close approximation (Fig. 1c). In other words, the sparse algorithm is allowed to select the locally significant snapshots, without the need for clustering. Indeed, this approach is equivalent to finding the best “cluster” for the queried point on-the-go. Thus, there is no need of “discretizing” the manifold into local subspaces. Also, the reconstruction quality in regions near intersection of two local subspaces need not be compromised.

3. Another possibility is to use non-linear approximation techniques [26], which is outside the scope of this article (and will be explored in a future article).

In this article, we propose to use sparse approximations from the domain of dictionary learning [27, 28]. Dictionary-based sparse approximations have already been applied to various numerical models including transport, acoustics [29] and biomedical [30] applications. In contrast to a reduced basis, an over-complete dictionary contains a large number of snapshots, some of which would essentially contain redundant information. For an approximation using linear combinations of the dictionary elements, multiple solutions would likely exist. The idea is to induce sparsity in the coefficients of the online problem, encouraging selection of a few snapshots that are the best candidates for the problem. The reader should beware that the underlying linear inseparability is not resolved by this approach. While the use of large dictionaries provide a better resolution of the non-linear manifold, the primary focus of this work is overcoming the added complexity resulting from use of large dictionaries.

Our work has some similarities as well as differences with earlier approaches. While Haasdonk et al. [13] uses snapshots to define solution subspaces, sparse regularization is not imposed for the dual reduced variable. Moreover, our approach does not attempt to compute a dual reduced basis like the works of Balajewicz et al. [14], Niakh et al. [15], Benaceur et al. [16]. Both Haasdonk et al. [13], Fauque et al. [17] tackle problems with node-to-node interaction, whereas our work includes problems with node-to-segment interactions and large “slip” between contact surfaces.

The article is structured as follows: first a very brief introduction to the discretized contact problem is given, followed by a brief introduction to sparse regression methods. Then, two proposed formulations of dictionary-based sparse methods to contact mechanics are introduced with numerical examples. Finally, some observations on the two methods are discussed and analyzed, prior to concluding the article.

## 2 Model Problem

Typical contact problems can be framed as inequality constrained minimization problems. Consider two elastic bodies with frictionless and non-adhesive surfaces, that might come into contact in a deformed state (Fig. 2). The contact problem can be framed as a minimization of the internal elastic energy while respecting the non-penetration constraints. To remain concise, the discretized version of this problem is given here. For a more detailed formulation, the reader may refer to the preceding publication [23] as well as the literature [5, 6, 8].

$$ \begin{gathered} \arg\min_{\mathbf{v}^h} \frac{1}{2} \mathbf{v}^h T \mathbf{K}(\mu) \mathbf{v}^h - \mathbf{v}^h T \mathbf{f}(\mu) \tag{1} \\ \mathbf{C}(\mu, \mathbf{v}^h) \mathbf{v}^h - \mathbf{g}(\mu, \mathbf{v}^h) \le \mathbf{0} \end{gathered} $$

where

$$ \bullet \mathbf{v}^h = \begin{bmatrix} \mathbf{v}_1^h \\ \mathbf{v}_2^h \end{bmatrix} \text{ is the stacked nodal displacement vector} $$
<figure><img src="image_1.png" /><figcaption>Figure 1: Illustration of low-rank, localized low-rank and dictionary-based approximation. Snapshots (green circles) in an *n*-dimensional system lie on a manifold (green-curve) of dimension lower than *n*. Low-rank subspaces have a dimension ≪ *n*. The approximation to the queried solution (⊕) and its closeness to the manifold can be visualized for the three methods.</figcaption></figure>

• $\mathbf{K}(\mu) = \begin{bmatrix} \mathbf{K}_1(\mu) & 0 \\ 0 & \mathbf{K}_2(\mu) \end{bmatrix}$ is the block elastic stiffness matrix

• $\mathbf{f}(\mu) = \begin{bmatrix} \mathbf{f}_1(\mu) \\ \mathbf{f}_2(\mu) \end{bmatrix}$ is the stacked nodal loading vector

• $\mathbf{g}(\mu, \mathbf{v}^h)$ is the node-to-segment distance between the surfaces of two bodies in reference configuration 2

• $\mathbf{C}(\mu, \mathbf{v}^h)\mathbf{v}^h$ gives the change in distance between the surfaces in the deformed configuration. with the subscripts 1 and 2 referring to the two bodies.

<footnote n="2">²This quantity depends on $\mathbf{v}^h$ because contact pairs are established in deformed configuration. See [23] for details.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Illustration of low-rank, localized low-rank and dictionary-based approximation. Snapshots (green circles) in an *n*-dimensional system lie on a manifold (green-curve) of dimension lower than *n*. Low-rank subspaces have a dimension ≪ *n*. The approximation to the queried solution (⊕) and its closeness to the manifold can be visualized for the three methods.</figcaption></figure>

• $\mathbf{K}(\mu) = \begin{bmatrix} \mathbf{K}_1(\mu) & 0 \\ 0 & \mathbf{K}_2(\mu) \end{bmatrix}$ is the block elastic stiffness matrix

• $\mathbf{f}(\mu) = \begin{bmatrix} \mathbf{f}_1(\mu) \\ \mathbf{f}_2(\mu) \end{bmatrix}$ is the stacked nodal loading vector

• $\mathbf{g}(\mu, \mathbf{v}^h)$ is the node-to-segment distance between the surfaces of two bodies in reference configuration 2

• $\mathbf{C}(\mu, \mathbf{v}^h)\mathbf{v}^h$ gives the change in distance between the surfaces in the deformed configuration. with the subscripts 1 and 2 referring to the two bodies.

<footnote n="2">²This quantity depends on $\mathbf{v}^h$ because contact pairs are established in deformed configuration. See [23] for details.</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 2: Kinematic description of two body contact problems with possibility of contact</figcaption></figure>

Solution to such an inequality constrained problem is given by the following Karush-Kuhn-Tucker (KKT) conditions, where an additional variable i.e. the Lagrange Multiplier $\lambda^h$ is introduced:

$$ K(\mu)u^h - f(\mu) + C^T(\mu, u^h)\lambda^h = 0 \quad (2a) $$

$$ \lambda^h \ge 0 \quad (2b) $$

$$ C(\mu, u^h)u^h - g(\mu, u^h) \le 0 \quad (2c) $$

$$ \lambda^h [C(\mu, u^h)u^h - g(\mu, u^h)] = 0 \quad (2d) $$

In this case, the Lagrange Multiplier $\lambda^h$ corresponds to the contact pressure between the two surfaces.

# 3 Dictionary learning with sparse methods: a brief overview

Before getting started with the contact problem, sparse regression techniques and corresponding literature are briefly introduced. This is followed by a brief discussion on their application to partial differential equations (PDEs).

## 3.1 Sparse regression methods

Sparse regression techniques have been under development in recent decades across various fields, especially in signal processing and statistics[31, 32]. The main idea behind these techniques is to approximate a signal with the few most suitable elements of a precomputed dictionary of signals [33], or to develop the simplest predictive model [34]. A typical sparse regression can be written as the following optimization statement:

$$ \begin{array}{ll} \min & ||\alpha||_p \\ \text{s.t.} & ||x - D\alpha|| < \varepsilon \end{array} \quad (3) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Kinematic description of two body contact problems with possibility of contact</figcaption></figure>

Solution to such an inequality constrained problem is given by the following Karush-Kuhn-Tucker (KKT) conditions, where an additional variable i.e. the Lagrange Multiplier $\lambda^h$ is introduced:

$$ K(\mu)u^h - f(\mu) + C^T(\mu, u^h)\lambda^h = 0 \qquad (2a) $$

$$ \lambda^h \ge 0 \qquad (2b) $$

$$ C(\mu, u^h)u^h - g(\mu, u^h) \le 0 \qquad (2c) $$

$$ \lambda^h [C(\mu, u^h)u^h - g(\mu, u^h)] = 0 \qquad (2d) $$

In this case, the Lagrange Multiplier $\lambda^h$ corresponds to the contact pressure between the two surfaces.

# 3 Dictionary learning with sparse methods: a brief overview

Before getting started with the contact problem, sparse regression techniques and corresponding literature are briefly introduced. This is followed by a brief discussion on their application to partial differential equations (PDEs).

## 3.1 Sparse regression methods

Sparse regression techniques have been under development in recent decades across various fields, especially in signal processing and statistics[31, 32]. The main idea behind these techniques is to approximate a signal with the few most suitable elements of a precomputed dictionary of signals [33], or to develop the simplest predictive model [34]. A typical sparse regression can be written as the following optimization statement:

$$ \begin{array}{ll} \min & ||\alpha||_p \\ \text{s.t.} & ||x - D\alpha|| < \varepsilon \end{array} \qquad (3) $$
where $x$ is the signal being approximated, $D$ could either be a dictionary of signals or a basis that possibly admits sparse representation, $\alpha$ is a vector of coefficients, $||\cdot||_p$ is a sparsity inducing norm, and $\varepsilon$ is the tolerance on the approximation error. It is also possible to write other forms of sparse regression problems where the approximation error is minimized and $p$-norm term is constrained to stay under a specific tolerance of sparsity or the minimization of a penalty form with a weighted sum of the two terms. A primer on sparse regression methods can be found in Mairal et al. [35].

It is common knowledge in the communities using sparse regression methods that the $l_0$-norm is a measure of cardinality, hence (3) with $p = 0$ is the natural choice in a theoretical sense. However, this problem is NP-hard in general, and therefore, it is more common to use $l_p$-norm with $p \in ]0, 1]$. Examples of such methods include the LASSO [36] and Dantzig Selector [37], which use $p = 1$. Other possibilities include the use of greedy algorithms, such as Orthogonal Matching Pursuit (OMP) [38] and weight $l_2$-norm minimization such as FOCUSS [39]. OMP and FOCUSS methods, which will be referred to in upcoming sections, are briefly detailed in Appendix A.

## 3.2 Numerical modelling using sparse methods

Sparse regression methods, though originally developed for application on pure data approximation or data compression problems, have seen recent interest in applications towards physics problems [34, 40, 41]. The methods that have been discussed previously are also meant for data approximation. At this point, the application of such methods to unconstrained PDE is discussed, which can be formulated as a residual minimization problem involving dictionary-based approximation and sparse regression. Consider a discretized problem with the state vector $v$ and residual $r(v)$. Given a dictionary $D$, the basic idea is to find an approximation $v \approx D\hat{v}$ that minimizes a sparsity inducing $l_p$-norm of $\hat{v}$ but also keeping the residual $r(D\hat{v})$ below a certain tolerance, as given in (4).

$$ \begin{array}{ll} \min & ||\hat{v}||_p \\ \text{s.t.} & ||r(D\hat{v})|| < \varepsilon \end{array} \qquad (4) $$

Dictionary-based sparse approximation for PDEs has been discussed in [29], along with the application of random sketching to residual evaluation which provides an efficient way of solving the minimization problem. This is similar to the approach of randomized SVD proposed in [42], where random combinations of rows/columns are used. The idea is to project the residual on a low-dimensional subspace in $colsp(D)$, where $colsp$ denotes the column-space. The matrix $B$, containing a basis corresponding to a random linear combination of dictionary columns, is used for this projection in (5). This projection essentially reduces the complexity of enforcing the residual constraint. Given a dictionary $D$ of size $m \times n$, a smaller matrix that has information from all columns of $D$ can be generated by multiplying a matrix $R$ of size $n \times L$ with $L \ll \min(m, n)$. After orthogonalization, the output matrix $B$ can be used for projecting the residual $r$.

$$ \min ||\hat{v}||_p \qquad (5a) $$

$$ \text{s.t.} \quad ||\mathbf{B}^T \mathbf{r}(\mathbf{D}\hat{\mathbf{v}})|| < \varepsilon \qquad (5b) $$

$$ \text{where } \mathbf{B} = \text{orth}(\mathbf{DR}), \quad \mathbf{R} \in \mathcal{U}_{[0,1]}^{n \times L} \qquad (5c) $$

where $\text{orth}(\cdot)$ denotes an orthogonalizing operation and $\mathcal{U}_{[0,1]}^{n \times L}$ denotes a random matrix of size $n \times L$ whose entries are drawn on the interval $[0, 1]$ using a uniform distribution. The number $L$
where $\boldsymbol{x}$ is the signal being approximated, $\mathbf{D}$ could either be a dictionary of signals or a basis that possibly admits sparse representation, $\boldsymbol{\alpha}$ is a vector of coefficients, $||\cdot||_p$ is a sparsity inducing norm, and $\varepsilon$ is the tolerance on the approximation error. It is also possible to write other forms of sparse regression problems where the approximation error is minimized and $p$-norm term is constrained to stay under a specific tolerance of sparsity or the minimization of a penalty form with a weighted sum of the two terms. A primer on sparse regression methods can be found in Mairal et al. [35].

It is common knowledge in the communities using sparse regression methods that the $\ell_0$-norm is a measure of cardinality, hence (3) with $p = 0$ is the natural choice in a theoretical sense. However, this problem is NP-hard in general, and therefore, it is more common to use $\ell_p$-norm with $p \in ]0, 1]$. Examples of such methods include the LASSO [36] and Dantzig Selector [37], which use $p = 1$. Other possibilities include the use of greedy algorithms, such as Orthogonal Matching Pursuit (OMP) [38] and weight $\ell_2$-norm minimization such as FOCUSS [39]. OMP and FOCUSS methods, which will be referred to in upcoming sections, are briefly detailed in Appendix A.

## 3.2 Numerical modelling using sparse methods

Sparse regression methods, though originally developed for application on pure data approximation or data compression problems, have seen recent interest in applications towards physics problems [34, 40, 41]. The methods that have been discussed previously are also meant for data approximation. At this point, the application of such methods to unconstrained PDE is discussed, which can be formulated as a residual minimization problem involving dictionary-based approximation and sparse regression. Consider a discretized problem with the state vector $\boldsymbol{v}$ and residual $\boldsymbol{r}(\boldsymbol{v})$. Given a dictionary $\mathbf{D}$, the basic idea is to find an approximation $\boldsymbol{v} \approx \hat{\boldsymbol{v}}$ that minimizes a sparsity inducing $\ell_p$-norm of $\hat{\boldsymbol{v}}$ but also keeping the residual $\boldsymbol{r}(\mathbf{D}\hat{\boldsymbol{v}})$ is below a certain tolerance, as given in (4).

$$ \begin{array}{ll} \min & ||\hat{\boldsymbol{v}}||_p \\ \text{s.t.} & ||\boldsymbol{r}(\mathbf{D}\hat{\boldsymbol{v}})|| < \varepsilon \end{array} \qquad (4) $$

Dictionary-based sparse approximation for PDEs has been discussed in [29], along with the application of random sketching to residual evaluation which provides an efficient way of solving the minimization problem. This is similar to the approach of randomized SVD proposed in [42], where random combinations of rows/columns are used. The idea is to project the residual on a low-dimensional subspace in $colsp(\mathbf{D})$, where $colsp$ denotes the column-space. The matrix $\mathbf{B}$, containing a basis corresponding to a random linear combination of dictionary columns, is used for this projection in (5). This projection essentially reduces the complexity of enforcing the residual constraint. Given a dictionary $\mathbf{D}$ of size $m \times n$, a smaller matrix that has information from all columns of $\mathbf{D}$ can be generated by multiplying a matrix $\mathbf{R}$ of size $n \times L$ with $L \ll \min(m, n)$. After orthogonalization, the output matrix $\mathbf{B}$ can be used for projecting the residual $\boldsymbol{r}$.

$$ \min ||\hat{\boldsymbol{v}}||_p \qquad (5a) $$

$$ \text{s.t.} \quad ||\mathbf{B}^T \boldsymbol{r}(\mathbf{D}\hat{\boldsymbol{v}})|| < \varepsilon \qquad (5b) $$

$$ \text{where } \mathbf{B} = \text{orth}(\mathbf{DR}), \quad \mathbf{R} \in \mathcal{U}_{[0,1]}^{n \times L} \qquad (5c) $$

where $\text{orth}(\cdot)$ denotes an orthogonalizing operation and $\mathcal{U}_{[0,1]}^{n \times L}$ denotes a random matrix of size $n \times L$ whose entries are drawn on the interval $[0, 1]$ using a uniform distribution. The number $L$
can be tuned to not only control the complexity of the minimization problem but also to crudely constrain the sparsity of the solution. For $L \ll n$, $colsp(B) \subset colsp(D)$ is low-dimensional and hence, the projection $B^T r$ is a “weaker” evaluation of $r$, which allows for a sparser $\hat{v}$ at the price of permitting residuals orthogonal to $B$.

Another related approach in which snapshots are directly used for reducing the model without any explicit compression is the so-called CUR decomposition [43]. In this method, not only columns but also rows from the snapshot matrix are used for developing a parametric regression. It works by algorithmically selecting a column and row subset from the snapshot matrix, and using them for developing a regression model. A contrasting feature of the CUR approach is the use of regression without invoking the partial differential equation or its weak form, resulting in a non-intrusive method.

In the proposed dictionary-based approach for contact mechanics, unlike the low-rank approach [14, 16], no compression of the contact pressure snapshots is applied. The dual solution can simply be restricted to the column-cone<fnref n="3" /> of the dictionary $D$, which might be reasonably rich to approximate the true contact pressure subcone. In the online phase, only a few snapshots will be selected to estimate the contact pressure field. This is justified because not all information in the dictionary is necessary for the reconstruction of a particular instance and hence, the dual reduced dofs will admit a sparse solution. Therefore, enough motivation exists to use sparsity-enforcing methods to choose a few columns of a large dictionary to approximate the dual solution.

# 4 Dictionary-based sparse methods for contact problems

This section introduces a dictionary-based sparse approximation algorithm for inequality-constrained problems. For the contact pressure dictionaries, the term “over-complete” is used loosely. It is simply used to imply that the dictionary has numerous snapshots, possibly larger than the number of contact dofs in the high-fidelity problem.

**Remark 1.** The authors investigated another idea to solve the dictionary-based sparse problem. It involved strong assumptions about the convexity of the feasible region defined by the non-penetration constraint. Therefore, the method is not included in the main article because of its limited scope. However, it is briefly explained in the Appendix B with an illustration.

A greedy method inspired by the Orthogonal Matching Pursuit (OMP) (see Appendix A) and the active set method is devised, referred to as greedy active-set method. The primal and dual reduced spaces are defined using a reduced basis and a dictionary, as detailed below:

* Primal reduced basis $\Phi$ is defined using the left singular vectors of primal snapshots, i.e. $\Phi \leftarrow svd(D_p, \delta)$ where $\delta$ is the truncation tolerance. The first singular vectors that contribute to $1 - \delta$ fraction of total energy (computed using squared singular values) are kept.

* The dual dictionary $D_d$ contains all the snapshots from the training phase.

The above choice is justified since the displacement field is known to be relatively separable compared to the contact pressure[23]. The unknowns of the online problem are the reduced displacement dofs $\hat{u}$ and the sparse contact pressure dofs $\hat{\lambda}_I$. Using the defined notations, a projection-based, sparse problem can be stated as:

<footnote n="3">³column-cone can be thought of as the subset of the column-space restricted to non-negative coefficients.</footnote>
can be tuned to not only control the complexity of the minimization problem but also to crudely constrain the sparsity of the solution. For $L \ll n$, $colsp(B) \subset colsp(D)$ is low-dimensional and hence, the projection $B^T r$ is a “weaker” evaluation of $r$, which allows for a sparser $\hat{v}$ at the price of permitting residuals orthogonal to $B$.

Another related approach in which snapshots are directly used for reducing the model without any explicit compression is the so-called CUR decomposition [43]. In this method, not only columns but also rows from the snapshot matrix are used for developing a parametric regression. It works by algorithmically selecting a column and row subset from the snapshot matrix, and using them for developing a regression model. A contrasting feature of the CUR approach is the use of regression without invoking the partial differential equation or its weak form, resulting in a non-intrusive method.

In the proposed dictionary-based approach for contact mechanics, unlike the low-rank approach [14, 16], no compression of the contact pressure snapshots is applied. The dual solution can simply be restricted to the column-cone<fnref n="3" /> of the dictionary $D$, which might be reasonably rich to approximate the true contact pressure subcone. In the online phase, only a few snapshots will be selected to estimate the contact pressure field. This is justified because not all information in the dictionary is necessary for the reconstruction of a particular instance and hence, the dual reduced dofs will admit a sparse solution. Therefore, enough motivation exists to use sparsity-enforcing methods to choose a few columns of a large dictionary to approximate the dual solution.

# 4 Dictionary-based sparse methods for contact problems

This section introduces a dictionary-based sparse approximation algorithm for inequality-constrained problems. For the contact pressure dictionaries, the term “over-complete” is used loosely. It is simply used to imply that the dictionary has numerous snapshots, possibly larger than the number of contact dofs in the high-fidelity problem.

**Remark 1.** The authors investigated another idea to solve the dictionary-based sparse problem. It involved strong assumptions about the convexity of the feasible region defined by the non-penetration constraint. Therefore, the method is not included in the main article because of its limited scope. However, it is briefly explained in the Appendix B with an illustration.

A greedy method inspired by the Orthogonal Matching Pursuit (OMP) (see Appendix A) and the active set method is devised, referred to as greedy active-set method. The primal and dual reduced spaces are defined using a reduced basis and a dictionary, as detailed below:

* Primal reduced basis $\Phi$ is defined using the left singular vectors of primal snapshots, i.e. $\Phi \leftarrow svd(D_p, \delta)$ where $\delta$ is the truncation tolerance. The first singular vectors that contribute to $1 - \delta$ fraction of total energy (computed using squared singular values) are kept.

* The dual dictionary $D_d$ contains all the snapshots from the training phase.

The above choice is justified since the displacement field is known to be relatively separable compared to the contact pressure[23]. The unknowns of the online problem are the reduced displacement dofs $\hat{u}$ and the sparse contact pressure dofs $\hat{\lambda}_I$. Using the defined notations, a projection-based, sparse problem can be stated as:

<footnote n="3">³column-cone can be thought of as the subset of the column-space restricted to non-negative coefficients.</footnote>
$$ \min \left\| \hat{\lambda} \right\|_{0} \qquad (6a) $$

$$ \text{s.t.} \quad \left[ \begin{array}{cc} \mathbf{\Phi}^T \mathbf{K} \mathbf{\Phi} & \mathbf{\Phi}^T \mathbf{C}^T \mathbf{D}_{d,I} \\ \mathbf{D}_{d,I}^T \mathbf{C} \mathbf{\Phi} & \mathbf{0} \end{array} \right] \left[ \hat{\mathbf{u}} \right] = \left[ \begin{array}{c} \mathbf{\Phi}^T \mathbf{f} \\ \mathbf{D}_{d,I}^T \mathbf{g} \end{array} \right] \qquad (6b) $$

$$ \text{and} \quad \mathbf{D}_{d}^{T} (\mathbf{C} \mathbf{\Phi} \hat{\mathbf{u}} - \mathbf{g}) \leq \tau \qquad (6c) $$

$$ \hat{\lambda} \geq 0 \qquad (6d) $$

where $\mathcal{I}$ indicates the active set and $\mathbf{D}_{d,I} = \mathbf{D}_{d}[\cdot, \mathcal{I}]$

(6b), (6c) and (6d) are based on the same principle as the KKT conditions (2), but applied over solution spaces defined by the reduced basis $\mathbf{\Phi}$ and the dictionary $\mathbf{D}_d$ (rather than a low-rank dual basis). Comparing (6) with the generic form of sparse regression statement (4), the balance equation is imposed as an equality constraint, instead of an inequality constraint where the residual norm is limited to a certain tolerance. Imposing limits on the residual norm is usually more relevant when using random sketching, like in (5) when the system is under-determined. Also, it is not straightforward to use standard sparse regression algorithms to solve inequality constraints (6c) and (6d).

The truncated primal basis $\mathbf{\Phi}$ cannot contain the entire information learnt from the training set. The solutions generated using the resulting subspace tend to have some false penetrations, which tend to drive the selection of some additional but irrelevant contact pressure snapshots from the dictionary $\mathbf{D}_d$ (see Appendix C). To avoid this issue, a small penetration $\tau$ is permitted in (6c) to avoid this spurious selection of snapshots. Since the unwanted spurious effects were observed to increase with truncation, the permitted penetration $\tau$ must be linked to the tolerance $\delta$. Unless specified, the $\tau := \delta$ is used in the rest of the article.

**Remark 2.** As a consequence of (6c), the non-penetration condition is not satisfied exactly by the sparse model. Thus, one might question the interest of using an LMM-based model over a penalty-based unconstrained model that is simpler to reduce. While this is a valid argument, the community remains interested in reduction of LMM-based models [13, 14, 16–18], possibly for reasons such as avoiding ill-conditioning problems of penalty methods.

The primary notion in greedy algorithms is to enrich the approximation by a rank-1 term in each step. For standard sparse regression problems like (4), the OMP algorithm starts with an empty coefficient vector and adds a single non-zero entry in one step. The position of the new non-zero entry corresponds to the vector in dictionary with the highest correlation to the residual. Using the same idea, the greedy active-set algorithm devised here performs enrichment by correlating the contact pressure snapshots with the violation of inequality constraint (6c). More precisely, the algorithm starts with a zero $\hat{\lambda}$ and in each enrichment step, the dictionary column with the highest correlation to the violations of the non-penetration condition is added in the greedy enrichment. Thus, enrichment index of the dictionary $p_{\text{add}}$ can be computed as:

$$ p_{\text{add}} := \arg \max_{p} \mathbf{D}_{d}[\cdot, p]^{T} (\mathbf{C} \mathbf{\Phi} \hat{\mathbf{u}}^{k} - \mathbf{g})^{+}, \tau \qquad (7) $$

where $(\mathbf{z})_i^{+}, \tau := \begin{cases} z_i & \text{if } z_i > \tau \\ 0 & \text{otherwise} \end{cases}$, is a hard-thresholding operator giving the part greater than $\tau$.
$$ \min \left\| \hat{\lambda} \right\|_0 \qquad (6a) $$

$$ \text{s.t.} \quad \left[ \begin{array}{cc} \mathbf{\Phi}^T \mathbf{K} \mathbf{\Phi} & \mathbf{\Phi}^T \mathbf{C}^T \mathbf{D}_{d,I} \\ \mathbf{D}_{d,I}^T \mathbf{C} \mathbf{\Phi} & \mathbf{0} \end{array} \right] \left[ \hat{\mathbf{u}} \right] = \left[ \begin{array}{c} \mathbf{\Phi}^T \mathbf{f} \\ \mathbf{D}_{d,I}^T \mathbf{g} \end{array} \right] \qquad (6b) $$

$$ \text{and} \quad \mathbf{D}_d^T (\mathbf{C}\mathbf{\Phi}\hat{\mathbf{u}} - \mathbf{g}) \leq \tau \qquad (6c) $$

$$ \hat{\lambda} \geq 0 \qquad (6d) $$

where $\mathcal{I}$ indicates the active set and $\mathbf{D}_{d,I} = \mathbf{D}_d[\cdot, \mathcal{I}]$

(6b), (6c) and (6d) are based on the same principle as the KKT conditions (2), but applied over solution spaces defined by the reduced basis $\mathbf{\Phi}$ and the dictionary $\mathbf{D}_d$ (rather than a low-rank dual basis). Comparing (6) with the generic form of sparse regression statement (4), the balance equation is imposed as an equality constraint, instead of an inequality constraint where the residual norm is limited to a certain tolerance. Imposing limits on the residual norm is usually more relevant when using random sketching, like in (5) when the system is under-determined. Also, it is not straightforward to use standard sparse regression algorithms to solve inequality constraints (6c) and (6d).

The truncated primal basis $\mathbf{\Phi}$ cannot contain the entire information learnt from the training set. The solutions generated using the resulting subspace tend to have some false penetrations, which tend to drive the selection of some additional but irrelevant contact pressure snapshots from the dictionary $\mathbf{D}_d$ (see Appendix C). To avoid this issue, a small penetration $\tau$ is permitted in (6c) to avoid this spurious selection of snapshots. Since the unwanted spurious effects were observed to increase with truncation, the permitted penetration $\tau$ must be linked to the tolerance $\delta$. Unless specified, the $\tau := \delta$ is used in the rest of the article.

**Remark 2.** As a consequence of (6c), the non-penetration condition is not satisfied exactly by the sparse model. Thus, one might question the interest of using an LMM-based model over a penalty-based unconstrained model that is simpler to reduce. While this is a valid argument, the community remains interested in reduction of LMM-based models [13, 14, 16–18], possibly for reasons such as avoiding ill-conditioning problems of penalty methods.

The primary notion in greedy algorithms is to enrich the approximation by a rank-1 term in each step. For standard sparse regression problems like (4), the OMP algorithm starts with an empty coefficient vector and adds a single non-zero entry in one step. The position of the new non-zero entry corresponds to the vector in dictionary with the highest correlation to the residual. Using the same idea, the greedy active-set algorithm devised here performs enrichment by correlating the contact pressure snapshots with the violation of inequality constraint (6c). More precisely, the algorithm starts with a zero $\hat{\lambda}$ and in each enrichment step, the dictionary column with the highest correlation to the violations of the non-penetration condition is added in the greedy enrichment. Thus, enrichment index of the dictionary $p_{\text{add}}$ can be computed as:

$$ p_{\text{add}} := \arg \max_{p} \mathbf{D}_{d}[\cdot, p]^{T} (\mathbf{C}\mathbf{\Phi}\hat{\mathbf{u}}^{k} - \mathbf{g})^{+}, \tau \qquad (7) $$

where $(\mathbf{z})_i^{+}, \tau := \begin{cases} z_i & \text{if } z_i > \tau \\ 0 & \text{otherwise} \end{cases}$, is a hard-thresholding operator giving the part greater than $\tau$.
(7) can be computed more efficiently using the reduced operators $\hat{C}(\mu, \hat{u}^{k-1})$ and $\hat{g}(\mu, \hat{u}^{k-1})$:

$$p_{\text{add}} := \arg \max_{p} (\hat{C}\hat{u}^k - \hat{g})^{+,\tau} [p] \quad (8)$$

where

$$\hat{C}(\mu, \hat{u}^{k-1}) = D_d^T C(\mu, \hat{u}^{k-1}) \Phi$$

$$\hat{g}(\mu, \hat{u}^{k-1}) = D_d^T g(\mu, \hat{u}^{k-1})$$

However, to satisfy the inequality constraint (6d) that prohibits non-negative pressures, it is necessary to not only enrich but also to eliminate terms that violate the constraints in each greedy step. This is done by eliminating the largest negative $\hat{\lambda}$ i.e. setting the largest negative entry to zero. In the absence of a negative entry, elimination is not performed. The elimination index $p_{\text{rem}}$ is computed as:

$$p_{\text{rem}} = \arg \min_{p} (\hat{\lambda}^k)^- [p] \quad (9)$$

where $(z)_i^- := \begin{cases} z_i & \text{if } z_i < 0 \\ 0 & \text{otherwise} \end{cases}$

Although the motivation for this algorithm was based on the OMP, the greedy active-set algorithm ultimately resembles a fixed-point active-set method. However, a maximum of one dual dictionary column is activated in a given iteration, thereby maintaining a small active set in intermediate iterations. The algorithm is given in Alg. 1, where Steps 8 to 10 contain the enrichment and the elimination process. In each iteration, only one of the two operations is performed. Due to the architecture of the algorithm, the complexity of operations does not scale with the dictionary size, except enrichment/elimination Steps 8 to 10. The small sizes of active set $\mathcal{I}$, which is directly related to the sparsity of $\hat{\lambda}$, prevent the size of mixed system in Step 6 from growing significantly during the iterative process.

Note that at any given greedy step, $p_{\text{rem}}$ will never be same as $p_{\text{add}}$. For any active dual dof indexed $\tilde{p} \in \mathcal{I}$, (6b) implies $\hat{C}[\tilde{p}]\hat{u}^k - \hat{g}[\tilde{p}] = 0 \implies p_{\text{add}} \neq \tilde{p} \implies p_{\text{add}} \neq p_{\text{rem}}$, since $p_{\text{rem}} \in \mathcal{I}$. Moreover, to avoid any possibility of an endless loop between same $p_{\text{add}}$ and $p_{\text{rem}}$, an upper limit to number of iterations $k_{\text{max}}$ is added.

## 4.1 Application to the Hertz problem

The Hertz contact problem between two half cylinders (Fig. 3) is considered, with $R_1 = R_2 = 1$. A finite element model with 466 quad elements in each half-cylinder is used as a high-fidelity model. Bilinear shape functions describe the displacement field, while node-centered piece-wise constant shape functions at the boundary describe the contact pressure field. Surface nodes also serve as Gauss quadrature points for boundary integration. The imposed displacement on the top cylinder $d \in (0, 0.3)$ is considered as the parameter. Each half-cylinder has 513 nodes and 466 elements of which 78 elements lie of the potential contact surface, that is, the semi-circular edges. Computational time for each snapshot is around 0.4-1.0s.

Details of the offline phase:

* Training set: Primal and dual snapshots were computed for four training sets of sizes 12, 30, 60, and 120 in the offline stage. The training points $P_{\text{tr}}$ are distributed uniformly in parametric space $\mathcal{P} = (0, 0.3]$.
(7) can be computed more efficiently using the reduced operators $\hat{C}(\mu, \hat{u}^{k-1})$ and $\hat{g}(\mu, \hat{u}^{k-1})$:

$$p_{\text{add}} := \arg \max_{p} (\hat{C}\hat{u}^k - \hat{g})^{+,\tau} [p] \quad (8)$$

where

$$\begin{aligned}
\hat{C}(\mu, \hat{u}^{k-1}) &= D_d^T C(\mu, \hat{u}^{k-1}) \Phi \\
\hat{g}(\mu, \hat{u}^{k-1}) &= D_d^T g(\mu, \hat{u}^{k-1})
\end{aligned}$$

However, to satisfy the inequality constraint (6d) that prohibits non-negative pressures, it is necessary to not only enrich but also to eliminate terms that violate the constraints in each greedy step. This is done by eliminating the largest negative $\hat{\lambda}$ i.e. setting the largest negative entry to zero. In the absence of a negative entry, elimination is not performed. The elimination index $p_{\text{rem}}$ is computed as:

$$p_{\text{rem}} = \arg \min_{p} (\hat{\lambda}^k)^{-} [p] \quad (9)$$

where $(z)_i^{-} := \begin{cases} z_i & \text{if } z_i < 0 \\ 0 & \text{otherwise} \end{cases}$

Although the motivation for this algorithm was based on the OMP, the greedy active-set algorithm ultimately resembles a fixed-point active-set method. However, a maximum of one dual dictionary column is activated in a given iteration, thereby maintaining a small active set in intermediate iterations. The algorithm is given in Alg. 1, where Steps 8 to 10 contain the enrichment and the elimination process. In each iteration, only one of the two operations is performed. Due to the architecture of the algorithm, the complexity of operations does not scale with the dictionary size, except enrichment/elimination Steps 8 to 10. The small sizes of active set $\mathcal{I}$, which is directly related to the sparsity of $\hat{\lambda}$, prevent the size of mixed system in Step 6 from growing significantly during the iterative process.

Note that at any given greedy step, $p_{\text{rem}}$ will never be same as $p_{\text{add}}$. For any active dual dof indexed $\tilde{p} \in \mathcal{I}$, (6b) implies $\hat{C}[\tilde{p}]\hat{u}^k - \hat{g}[\tilde{p}] = 0 \implies p_{\text{add}} \neq \tilde{p} \implies p_{\text{add}} \neq p_{\text{rem}}$, since $p_{\text{rem}} \in \mathcal{I}$. Moreover, to avoid any possibility of an endless loop between same $p_{\text{add}}$ and $p_{\text{rem}}$, an upper limit to number of iterations $k_{\text{max}}$ is added.

## 4.1 Application to the Hertz problem

The Hertz contact problem between two half cylinders (Fig. 3) is considered, with $R_1 = R_2 = 1$. A finite element model with 466 quad elements in each half-cylinder is used as a high-fidelity model. Bilinear shape functions describe the displacement field, while node-centered piece-wise constant shape functions at the boundary describe the contact pressure field. Surface nodes also serve as Gauss quadrature points for boundary integration. The imposed displacement on the top cylinder $d \in (0, 0.3)$ is considered as the parameter. Each half-cylinder has 513 nodes and 466 elements of which 78 elements lie of the potential contact surface, that is, the semi-circular edges. Computational time for each snapshot is around 0.4-1.0s.

Details of the offline phase:

* Training set: Primal and dual snapshots were computed for four training sets of sizes 12, 30, 60, and 120 in the offline stage. The training points $P_{\text{tr}}$ are distributed uniformly in parametric space $\mathcal{P} = (0, 0.3]$.
<caption>**Algorithm 1 Greedy active-set algorithm**</caption>

1: Input: Queried value of parameter $\mu$
2: Given: Primal basis $\Phi$ and dual Dictionary $D_d$
Reduced operators $\Phi^T K \Phi$ and $\Phi^T f$
$\triangleright$ can be built offline
3: Initialize: $k = 0, \mathcal{I} = \emptyset$
4: **while** $\hat{u}$ and $\mathcal{I}$ not converge AND $k < k_{\max}$ **do**
5: Build reduced constraint operators $\hat{C}(\mu, \hat{u}^{k-1})$ and $\hat{g}(\mu, \hat{u}^{k-1})$
6: Solve the linear system (6b)
7: Set $\hat{\lambda}^k[\mathcal{I}] \leftarrow \hat{\lambda}_{\mathcal{I}}$ and $\hat{\lambda}^k[\mathcal{I}^c] \leftarrow 0$
$\triangleright$ $\mathcal{I}^c$ is the complementary set of $\mathcal{I}$
Enrich/Eliminate decision:
8: **if** $\hat{\lambda}^k \ge 0$ **then**
$\triangleright$ All multipliers are non-negative, so enrich
$$p_{\text{add}} = \arg \max_p (\hat{C}\Phi\hat{u}^k - \hat{g})^{+,\tau}$$
$$\mathcal{I} \leftarrow \mathcal{I} \cup \{p_{\text{add}}\}$$
9: **else:**
$\triangleright$ Eliminate negative multiplier
$$p_{\text{rem}} = \arg \min_p (\hat{\lambda}^k)^- [p]$$
$$\mathcal{I} \leftarrow \mathcal{I} \setminus \{p_{\text{rem}}\}$$
10: **end if**
11: $k \leftarrow k + 1$
12: **end while**
13: Reconstruct $\mathbf{u} = \Phi\hat{u}$ and $\lambda = D_d\hat{\lambda}$
14: Output: $\mathbf{u}, \lambda$

<figure><img src="image_3.png" /><figcaption>Figure 3: Hertz problem: two half cylinders loaded against each other and the resulting contact pressure snapshots. Displacement $d$, imposed on the top cylinder, is treated as a parameter in the reduced model.</figcaption></figure>
<caption>**Algorithm 1 Greedy active-set algorithm**</caption>

1: Input: Queried value of parameter $\mu$
2: Given: Primal basis $\Phi$ and dual Dictionary $D_d$
Reduced operators $\Phi^T K \Phi$ and $\Phi^T f$
$\triangleright$ can be built offline
3: Initialize: $k = 0, \mathcal{I} = \emptyset$
4: **while** $\hat{u}$ and $\mathcal{I}$ not converge AND $k < k_{\max}$ **do**
5: Build reduced constraint operators $\hat{C}(\mu, \hat{u}^{k-1})$ and $\hat{g}(\mu, \hat{u}^{k-1})$
6: Solve the linear system (6b)
7: Set $\hat{\lambda}^k[\mathcal{I}] \leftarrow \hat{\lambda}_{\mathcal{I}}$ and $\hat{\lambda}^k[\mathcal{I}^c] \leftarrow 0$
$\triangleright$ $\mathcal{I}^c$ is the complementary set of $\mathcal{I}$
Enrich/Eliminate decision:
8: **if** $\hat{\lambda}^k \ge 0$ **then**
$\triangleright$ All multipliers are non-negative, so enrich
$$p_{\text{add}} = \arg \max_{\mathcal{I}} (\hat{C}\Phi\hat{u}^k - \hat{g})^{+,\tau}$$
$$\mathcal{I} \leftarrow \mathcal{I} \cup \{p_{\text{add}}\}$$
9: **else:**
$\triangleright$ Eliminate negative multiplier
$$p_{\text{rem}} = \arg \min_{\mathcal{I}} (\hat{\lambda}^k)^{-} [p]$$
$$\mathcal{I} \leftarrow \mathcal{I} \setminus \{p_{\text{rem}}\}$$
10: **end if**
11: $k \leftarrow k + 1$
12: **end while**
13: Reconstruct $u = \Phi\hat{u}$ and $\lambda = D_d\hat{\lambda}$
14: Output: $u, \lambda$

<figure><img src="image_4.png" /><figcaption>Figure 3: Hertz problem: two half cylinders loaded against each other and the resulting contact pressure snapshots. Displacement $d$, imposed on the top cylinder, is treated as a parameter in the reduced model.</figcaption></figure>
* Validation Set: $P_{\text{val}} \subset P$ consisting 119 points, which are also the mid-points of the fourth (and largest) training set points, is used to study reconstruction errors.

The primal basis is computed using SVD with truncation parameter $\delta$, while the dual dictionary is composed of contact pressure snapshots. A spy plot of the dual dictionary is shown in Fig. 4.

<figure><img src="image_4.png" /><figcaption>Figure 4: A spy pattern of dual dictionary with 30 elements for the Hertz problem. The 79 points on the y-axis are the dofs on the potential contact surface. The dictionary columns are arranged in the increasing order of loading parameter $d$.</figcaption></figure>

Using the greedy active-set method, two instances of contact pressure reconstructions are shown in Fig. 5. The dictionary snapshots selected by the algorithm for reconstruction are also displayed. As expected, only a few dictionary elements, two in this case, are chosen. These two selected snapshots match the contact position required for the specific instance, indicating that the enrichment and elimination process is effectively selecting the correct snapshots for reconstruction.

Mean reconstruction errors over the validation set improve by more than an order of magnitude as dictionary size is increased from 12 to 120 for the case $\delta = 10^{-10}$ (Fig. 6a). However, such gains with over-complete dictionaries diminish for the case $\delta = 10^{-6}$, probably because the primal truncation starts to become the dominant source of error, rather than the lack of contact pressure separability.

Mean computational time on the validation set (Fig. 6c) initially increases with the dictionary size, but later stabilises. However, it is evident that the computational time is primarily influenced primarily by the number of iterations (Fig. 6b). This observation is further supported by Fig. 6d, where the computational time per iteration remains relatively constant w.r.t. the dictionary size. The small increases in computational time per iteration with dictionary size can be explained due to the increase in the dimension of primal basis $\Phi$, especially in the case of $\delta = 10^{-10}$ where the primal rank is strongly dependent on the training set size.

Detailed plots of reconstruction errors over the validation set $P_{\text{val}}$ are shown in Fig. 7 for primal truncation tolerance $\delta \in \{10^{-10}, 10^{-8}, 10^{-6}\}$; and each curve corresponds to one of the training sets. An unexpected observation in the case of $\delta = 10^{-6}$ is that certain regions have higher dual errors for larger training sets. This happens due to the spurious penetrations discussed previously (Section 4) and in Appendix C. Even after relaxation of the non-penetration constraint by the value of $\tau$ (see (6c)), spurious selections are still possible for higher values of $\delta$. Moreover, the errors tend to be larger in general near $d = 0.0$, due to higher sensitivity w.r.t.
* Validation Set: $P_{\text{val}} \subset P$ consisting 119 points, which are also the mid-points of the fourth (and largest) training set points, is used to study reconstruction errors.

The primal basis is computed using SVD with truncation parameter $\delta$, while the dual dictionary is composed of contact pressure snapshots. A spy plot of the dual dictionary is shown in Fig. 4.

<figure><img src="image_5.png" /><figcaption>Figure 4: A spy pattern of dual dictionary with 30 elements for the Hertz problem. The 79 points on the y-axis are the dofs on the potential contact surface. The dictionary columns are arranged in the increasing order of loading parameter $d$.</figcaption></figure>

Using the greedy active-set method, two instances of contact pressure reconstructions are shown in Fig. 5. The dictionary snapshots selected by the algorithm for reconstruction are also displayed. As expected, only a few dictionary elements, two in this case, are chosen. These two selected snapshots match the contact position required for the specific instance, indicating that the enrichment and elimination process is effectively selecting the correct snapshots for reconstruction.

Mean reconstruction errors over the validation set improve by more than an order of magnitude as dictionary size is increased from 12 to 120 for the case $\delta = 10^{-10}$ (Fig. 6a). However, such gains with over-complete dictionaries diminish for the case $\delta = 10^{-6}$, probably because the primal truncation starts to become the dominant source of error, rather than the lack of contact pressure separability.

Mean computational time on the validation set (Fig. 6c) initially increases with the dictionary size, but later stabilises. However, it is evident that the computational time is primarily influenced primarily by the number of iterations (Fig. 6b). This observation is further supported by Fig. 6d, where the computational time per iteration remains relatively constant w.r.t. the dictionary size. The small increases in computational time per iteration with dictionary size can be explained due to the increase in the dimension of primal basis $\Phi$, especially in the case of $\delta = 10^{-10}$ where the primal rank is strongly dependent on the training set size.

Detailed plots of reconstruction errors over the validation set $P_{\text{val}}$ are shown in Fig. 7 for primal truncation tolerance $\delta \in \{10^{-10}, 10^{-8}, 10^{-6}\}$; and each curve corresponds to one of the training sets. An unexpected observation in the case of $\delta = 10^{-6}$ is that certain regions have higher dual errors for larger training sets. This happens due to the spurious penetrations discussed previously (Section 4) and in Appendix C. Even after relaxation of the non-penetration constraint by the value of $\tau$ (see (6c)), spurious selections are still possible for higher values of $\delta$. Moreover, the errors tend to be larger in general near $d = 0.0$, due to higher sensitivity w.r.t.
<figure><img src="image_5.png" /><figcaption>Figure 5: Greedy active-set reconstruction of contact pressure for certain parametric instances in the validation set of the Hertz problem. Elements of the dictionary chosen by the algorithm are shown and the reconstruction errors are given. Dictionary of size 12 is used in these examples.</figcaption></figure>

$d$ in this region. The cross (×) markers, indicating cases where algorithm does not converge due to oscillations between $p_{\text{add}}$ and $p_{\text{rem}}$.

The sparsity pattern of the dual reduced dofs $\hat{\lambda}$ for the training set of size 30 is shown in Fig. 8. As the snapshots in $D_d$ are arranged in increasing order of $d$, most of the columns of the sparsity pattern have two close non-zeros, indicating that the nearest two snapshots to the targeted reconstruction were chosen by the algorithm. Spurious selections are seen in some rows, but the values of corresponding coefficients are quite small and do not have any significant impact on the solution.

## 4.2 Application to the Ironing problem

Unlike the Hertz problem, other contact problems might exhibit larger changes in the contact zone, like the ironing problem, where an iron block is pressed against a flat slab and moved along its length (Fig. 9). The contact zone varies strongly with the horizontal displacement parameter $d_x$ (Fig. 10), which worsens the separability [23]. Therefore, over-complete dictionaries might be more relevant in generating an effective reduced model.

The ironing problem is considered with two simplifications: the surfaces are assumed to be frictionless and the problem is *quasi-static*, i.e. the iron is moved very slowly on the slab surface. The two bodies and surfaces are discretized using same type of elements as the Hertz problem. The iron is discretized using a structured mesh of 15 × 15 nodes, while the slab is discretized using a Cartesian mesh of 10 × 50 nodes.

The parametric space is defined as $\mathcal{P} := \{d_x \mid 0 \le d_x \le L\}$, given length of the slab $L$. The vertical displacement is fixed at $d_y = 0.3$. The training sets $\mathcal{P}_{\text{tr}}^n$ are defined in a nested manner, where $n$ is the nested level index. Each training set contains uniformly distributed points in $\mathcal{P}$, such that $\mathcal{P}_{\text{tr}}^n \subset \mathcal{P}_{\text{tr}}^{n+1}$ and $\#\mathcal{P}_{\text{tr}}^n = 2^n + 1$. Training sets from nested levels $n = 3$ to $n = 7$ are used, whereas the validation set is defined as $\mathcal{P}_{\text{val}} = \mathcal{P}_{\text{tr}}^8 \setminus \mathcal{P}_{\text{tr}}^7$. Computational time for each snapshot is around 0.85 seconds to 1.1 seconds.

The dual dictionary $D_d$ is composed of contact pressure snapshots at the $2^n + 1$ training points, while the primal basis consists of the left singular vectors truncated by filtering the
<figure><img src="image_6.png" /><figcaption>Figure 5: Greedy active-set reconstruction of contact pressure for certain parametric instances in the validation set of the Hertz problem. Elements of the dictionary chosen by the algorithm are shown and the reconstruction errors are given. Dictionary of size 12 is used in these examples.</figcaption></figure>

$d$ in this region. The cross (×) markers, indicating cases where algorithm does not converge due to oscillations between $p_{\text{add}}$ and $p_{\text{rem}}$.

The sparsity pattern of the dual reduced dofs $\hat{\lambda}$ for the training set of size 30 is shown in Fig. 8. As the snapshots in $D_d$ are arranged in increasing order of $d$, most of the columns of the sparsity pattern have two close non-zeros, indicating that the nearest two snapshots to the targeted reconstruction were chosen by the algorithm. Spurious selections are seen in some rows, but the values of corresponding coefficients are quite small and do not have any significant impact on the solution.

## 4.2 Application to the Ironing problem

Unlike the Hertz problem, other contact problems might exhibit larger changes in the contact zone, like the ironing problem, where an iron block is pressed against a flat slab and moved along its length (Fig. 9). The contact zone varies strongly with the horizontal displacement parameter $d_x$ (Fig. 10), which worsens the separability [23]. Therefore, over-complete dictionaries might be more relevant in generating an effective reduced model.

The ironing problem is considered with two simplifications: the surfaces are assumed to be frictionless and the problem is *quasi-static*, i.e. the iron is moved very slowly on the slab surface. The two bodies and surfaces are discretized using same type of elements as the Hertz problem. The iron is discretized using a structured mesh of 15 × 15 nodes, while the slab is discretized using a Cartesian mesh of 10 × 50 nodes.

The parametric space is defined as $\mathcal{P} := \{d_x | 0 \le d_x \le L\}$, given length of the slab $L$. The vertical displacement is fixed at $d_y = 0.3$. The training sets $\mathcal{P}_{\text{tr}}^n$ are defined in a nested manner, where $n$ is the nested level index. Each training set contains uniformly distributed points in $\mathcal{P}$, such that $\mathcal{P}_{\text{tr}}^n \subset \mathcal{P}_{\text{tr}}^{n+1}$ and $\#\mathcal{P}_{\text{tr}}^n = 2^n + 1$. Training sets from nested levels $n = 3$ to $n = 7$ are used, whereas the validation set is defined as $\mathcal{P}_{\text{val}} = \mathcal{P}_{\text{tr}}^8 \setminus \mathcal{P}_{\text{tr}}^7$. Computational time for each snapshot is around 0.85 seconds to 1.1 seconds.

The dual dictionary $D_d$ is composed of contact pressure snapshots at the $2^n + 1$ training points, while the primal basis consists of the left singular vectors truncated by filtering the
<figure><img src="image_6.png" /><figcaption>(a) Mean relative errors (bar plot, left y-axis) and primal rank (line plot, right y-axis). + labels indicate cases where $\#P_{tr} > \#FOM$ contact dofs.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b) Mean number of iterations</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) Mean computation time per iteration (s)</figcaption></figure>

Figure 6: Evolution of (a) relative reconstruction errors and primal rank, (b) number of iterations, (c) total and (d) per iteration computation time with training set size for the Hertz problem. Mean values of these quantities over the validation set $P_{val}$ are shown. Primal and dual errors are computed using $\mathcal{H}^1$ and $\mathcal{L}^2$ norm, respectively. Quantities are shown for different primal truncation tolerance $\delta \in \{10^{-6}, 10^{-8}, 10^{-10}\}$. Note that computation times exclude the construction time of non-linear operators.
<figure><img src="image_7.png" /><figcaption>(a) Mean relative errors (bar plot, left y-axis) and primal rank (line plot, right y-axis). + labels indicate cases where $\#P_{tr} > \#FOM$ contact dofs.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(b) Mean number of iterations</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d) Mean computation time per iteration (s)</figcaption></figure>

Figure 6: Evolution of (a) relative reconstruction errors and primal rank, (b) number of iterations, (c) total and (d) per iteration computation time with training set size for the Hertz problem. Mean values of these quantities over the validation set $P_{val}$ are shown. Primal and dual errors are computed using $\mathcal{H}^1$ and $\mathcal{L}^2$ norm, respectively. Quantities are shown for different primal truncation tolerance $\delta \in \{10^{-6}, 10^{-8}, 10^{-10}\}$. Note that computation times exclude the construction time of non-linear operators.
<figure><img src="image_9.png" /><figcaption>(a) Primal and dual errors for $\delta = 10^{-10}$</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Primal and dual errors for $\delta = 10^{-8}$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(c) Primal and dual errors for $\delta = 10^{-6}$</figcaption></figure>

Figure 7: Detailed primal and Dual reconstruction errors for the Hertz problem using the training sets 1,2 and 3. The validation set is common to all curves, consisting of 60 equispaced points in $\mathcal{P}$. Each subfigure corresponds to a different truncation tolerance $\delta$ for the primal basis. Primal and dual errors are computed using $\mathcal{H}^1$ and $\mathcal{L}^2$ norm, respectively. Cross-marks $\times$ indicate points where the greedy active-set algorithm reached max the number of iterations (50) and did not converge to the defined tolerance $10^{-5}$
<figure><img src="image_10.png" /><figcaption>(a) Primal and dual errors for $\delta = 10^{-10}$</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Primal and dual errors for $\delta = 10^{-8}$</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) Primal and dual errors for $\delta = 10^{-6}$</figcaption></figure>

Figure 7: Detailed primal and Dual reconstruction errors for the Hertz problem using the training sets 1,2 and 3. The validation set is common to all curves, consisting of 60 equispaced points in $\mathcal{P}$. Each subfigure corresponds to a different truncation tolerance $\delta$ for the primal basis. Primal and dual errors are computed using $\mathcal{H}^1$ and $\mathcal{L}^2$ norm, respectively. Cross-marks $\times$ indicate points where the greedy active-set algorithm reached max the number of iterations (50) and did not converge to the defined tolerance $10^{-5}$
<figure><img src="image_12.png" /><figcaption>Figure 8: Sparsity of dual reduced dofs $\hat{\lambda}$ selected by greedy active-set method for reconstructions in the validation set of the Hertz problem using a dual dictionary of size 30. Dots indicate the non-zero positions.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 9: Ironing problem: The iron block is pressed against the flat slab by a displacement $d_y$ and moved horizontally. The horizontal displacement of the iron block $d_x \in [0, L]$ is treated as the parameter in the reduced model</figcaption></figure>

singular values w.r.t. a fixed tolerance $\delta$.

Figs. 11 to 13 contain results of the greedy active-set method applied to the ironing problem. The following are the observations that contrast with the results of Hertz problem:

1. **Inseparability:** Fig. 11 shows the contrast between approximations by dictionaries of different sizes. Its easy to see that the effects of inseparability are more pronounced in the case of ironing problem.

2. **Spurious selections:** Comparing Fig. 8 and Fig. 12, the sparsity pattern for ironing problem does not exhibit any spurious snapshot selections. This difference can be attributed to the stronger variations in the contact position; unlike the Hertz problem where contact
<figure><img src="image_13.png" /><figcaption>Figure 8: Sparsity of dual reduced dofs $\hat{\lambda}$ selected by greedy active-set method for reconstructions in the validation set of the Hertz problem using a dual dictionary of size 30. Dots indicate the non-zero positions.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 9: Ironing problem: The iron block is pressed against the flat slab by a displacement $d_y$ and moved horizontally. The horizontal displacement of the iron block $d_x \in [0, L]$ is treated as the parameter in the reduced model</figcaption></figure>

singular values w.r.t. a fixed tolerance $\delta$.

Figs. 11 to 13 contain results of the greedy active-set method applied to the ironing problem. The following are the observations that contrast with the results of Hertz problem:

1. **Inseparability:** Fig. 11 shows the contrast between approximations by dictionaries of different sizes. Its easy to see that the effects of inseparability are more pronounced in the case of ironing problem.

2. **Spurious selections:** Comparing Fig. 8 and Fig. 12, the sparsity pattern for ironing problem does not exhibit any spurious snapshot selections. This difference can be attributed to the stronger variations in the contact position; unlike the Hertz problem where contact
<figure><img src="image_14.png" /><figcaption>Figure 10: Sample snapshots of contact pressure for ironing problem</figcaption></figure>

pressure curves were centred at the same position. The activation-deactivation process is,
thus, less ambiguous and less sensitive to primal truncations in the case of the ironing
problem.

3. **Reconstruction errors** in Fig. 13a are quite similar for $\delta = 10^{-10}$ and $\delta = 10^{-8}$ possi-
bly because the bottlenecks are the inseparability and finite, although large, size of dual
dictionary.

4. **Number of iterations and computational time** seem to increase significantly with dictionary size (see Fig. 13b and Fig. 13c). A possible explanation for this behaviour could be that the larger dictionaries of the ironing problem, coupled with stronger parametric dependence of contact zone, force the greedy active-set to perform more iterations to find the appropriate dual dictionary candidates and contact pairs.

<figure><img src="image_15.png" /><figcaption>Figure 11: Contact pressure reconstructions for certain parametric instances of the ironing problem using greedy active-set method. The dictionary nested level, dictionary elements chosen by the algorithm and the reconstruction errors are also shown.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure 10: Sample snapshots of contact pressure for ironing problem</figcaption></figure>

pressure curves were centred at the same position. The activation-deactivation process is,
thus, less ambiguous and less sensitive to primal truncations in the case of the ironing
problem.

3. **Reconstruction errors** in Fig. 13a are quite similar for $\delta = 10^{-10}$ and $\delta = 10^{-8}$ possi-
bly because the bottlenecks are the inseparability and finite, although large, size of dual
dictionary.

4. **Number of iterations and computational time** seem to increase significantly with dictionary size (see Fig. 13b and Fig. 13c). A possible explanation for this behaviour could be that the larger dictionaries of the ironing problem, coupled with stronger parametric dependence of contact zone, force the greedy active-set to perform more iterations to find the appropriate dual dictionary candidates and contact pairs.

<figure><img src="image_16.png" /><figcaption>Figure 11: Contact pressure reconstructions for certain parametric instances of the ironing problem using greedy active-set method. The dictionary nested level, dictionary elements chosen by the algorithm and the reconstruction errors are also shown.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 12: Sparsity of dual reduced dofs $\hat{\lambda}$ selected by greedy active-set method for various reconstructions of ironing problem using the dual dictionary of nested level 4</figcaption></figure>

## 4.3 Application to a two-parameter ironing problem

To validate the same approach on a multidimensional parametric domain, the ironing problem with an additional parameter, i.e. the vertical displacement is considered. Thus, the new parametric space is defined as $\mathcal{P} := \{(d_x, d_y) | 0 \le d_x \le L, 0.1 \le d_y \le 0.3\}$. In the 2D parametric space, the nested training sets are defined such that the number of points is doubled in each parametric direction with increasing levels. The resulting error profiles and computational times are given in Fig. 14, with the same trends as the single parameter case.

# 5 Discussions

The application of dictionary methods was proposed to mitigate the challenges posed by linear inseparability of the contact pressure. Unlike a classical reduced basis, an over-complete dictionary has the potential to preserve local information computed in the offline stage, which would otherwise be lost by truncation of a reduced basis. The same has been demonstrated by building reasonably accurate reconstructions using dictionary methods.

The reconstruction errors of the Hertz problem decay consistently when richer training sets are used, and the computational time per iteration seems to be fairly constant, but change slightly with the primal rank. For the ironing problem with large contact variations, the reconstruction errors decays comparatively slower than Hertz problem. Though the dictionary size didn't have any significant influence on the computational time per iteration, the number of iterations (and total time) increased significantly for larger dictionaries in case of ironing problem. This can be explained by the inherently high non-linear behaviour that comes with the large contact surface and could possibly be reduced by using smarter initialization of the contact pairs or by increasing the hierarchical level of the training set gradually.

As mentioned earlier, the resolution of the inseparabilities and reducing the subsequent errors is not within the scope of this article. The primary objective is to focus on efficient algorithms to utilise over-complete dictionaries for contact problems. A comparison of the greedy active set method with its non-sparse version is shown in Table. 1. The sparse method provides considerable gains in computational time but not in terms of reconstruction accuracy.
<figure><img src="image_17.png" /><figcaption>Figure 12: Sparsity of dual reduced dofs $\hat{\lambda}$ selected by greedy active-set method for various reconstructions of ironing problem using the dual dictionary of nested level 4</figcaption></figure>

## 4.3 Application to a two-parameter ironing problem

To validate the same approach on a multidimensional parametric domain, the ironing problem with an additional parameter, i.e. the vertical displacement is considered. Thus, the new parametric space is defined as $\mathcal{P} := \{(d_x, d_y) | 0 \le d_x \le L, 0.1 \le d_y \le 0.3\}$. In the 2D parametric space, the nested training sets are defined such that the number of points is doubled in each parametric direction with increasing levels. The resulting error profiles and computational times are given in Fig. 14, with the same trends as the single parameter case.

# 5 Discussions

The application of dictionary methods was proposed to mitigate the challenges posed by linear inseparability of the contact pressure. Unlike a classical reduced basis, an over-complete dictionary has the potential to preserve local information computed in the offline stage, which would otherwise be lost by truncation of a reduced basis. The same has been demonstrated by building reasonably accurate reconstructions using dictionary methods.

The reconstruction errors of the Hertz problem decay consistently when richer training sets are used, and the computational time per iteration seems to be fairly constant, but change slightly with the primal rank. For the ironing problem with large contact variations, the reconstruction errors decays comparatively slower than Hertz problem. Though the dictionary size didn't have any significant influence on the computational time per iteration, the number of iterations (and total time) increased significantly for larger dictionaries in case of ironing problem. This can be explained by the inherently high non-linear behaviour that comes with the large contact surface and could possibly be reduced by using smarter initialization of the contact pairs or by increasing the hierarchical level of the training set gradually.

As mentioned earlier, the resolution of the inseparabilities and reducing the subsequent errors is not within the scope of this article. The primary objective is to focus on efficient algorithms to utilise over-complete dictionaries for contact problems. A comparison of the greedy active set method with its non-sparse version is shown in Table. 1. The sparse method provides considerable gains in computational time but not in terms of reconstruction accuracy.
<figure><img src="image_17.png" /><figcaption>(a) Mean reconstruction errors (bar plot, left y-axis) and primal rank (line plot, right y-axis). + labels indicate cases where $\#\mathcal{P}_{\text{tr}} > \#\text{FOM}$ contact dofs.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 13: Evolution of (a) mean reconstruction (relative) errors and primal rank, (b) number of iterations, (c) total and (d) per iteration computation time with training set size for the Ironing problem. Mean values of these quantities over validation set $\mathcal{P}_{\text{val}}$ are shown. Primal and dual errors are computed using $\mathcal{H}^1$ and $\mathcal{L}^2$ norm, respectively. Quantities are shown for different primal truncation tolerance $\delta \in \{10^6, 10^{-8}, 10^{-10}\}$. Note that computation time excludes the construction time of non-linear operators.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>(a) Mean reconstruction errors (bar plot, left y-axis) and primal rank (line plot, right y-axis). + labels indicate cases where $\#\mathcal{P}_{\text{tr}} > \#\text{FOM}$ contact dofs.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Figure 13: Evolution of (a) mean reconstruction (relative) errors and primal rank, (b) number of iterations, (c) total and (d) per iteration computation time with training set size for the Ironing problem. Mean values of these quantities over validation set $\mathcal{P}_{\text{val}}$ are shown. Primal and dual errors are computed using $\mathcal{H}^1$ and $\mathcal{L}^2$ norm, respectively. Quantities are shown for different primal truncation tolerance $\delta \in \{10^6, 10^{-8}, 10^{-10}\}$. Note that computation time excludes the construction time of non-linear operators.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 14: Evolution of (a) mean reconstruction (relative) errors and rank of primal basis, (b) mean computation time and (c) mean computation time per iteration with training set size for the two-parameter ironing problem.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 14: Evolution of (a) mean reconstruction (relative) errors and rank of primal basis, (b) mean computation time and (c) mean computation time per iteration with training set size for the two-parameter ironing problem.</figcaption></figure>
<table><caption>(a) Hertz problem with 30 snapshots in training set.</caption><thead><tr><th></th><th><em>Non-sparse</em><br>Active set with CPG<br>tolerance = 1 × 10<sup>-2</sup></th><th><em>Sparse</em><br>Greedy active set</th></tr></thead><tbody><tr><td>Primal RB size</td><td>22</td><td>22</td></tr><tr><td>Dual RB/Dictionary size</td><td>21</td><td>30</td></tr><tr><td>Time (seconds)</td><td>1.7 × 10<sup>-2</sup>s</td><td>4 × 10<sup>-3</sup> s</td></tr><tr><td>Primal reconstruction errors</td><td>1.4 × 10<sup>-3</sup></td><td>2 × 10<sup>-3</sup></td></tr><tr><td>Dual reconstruction errors</td><td>5 × 10<sup>-2</sup></td><td>5 × 10<sup>-2</sup></td></tr></tbody></table>

<table><caption>(b) Ironing problem with 129 snapshots in training set.</caption><thead><tr><th></th><th><em>Non-sparse</em><br>Active set with CPG<br>tolerance = 5 × 10<sup>-2</sup></th><th><em>Sparse</em><br>Greedy active set</th></tr></thead><tbody><tr><td>Primal RB size</td><td>99</td><td>99</td></tr><tr><td>Dual RB/Dictionary size</td><td>52</td><td>129</td></tr><tr><td>Time (seconds)</td><td>0.08s</td><td>0.02 s</td></tr><tr><td>Primal reconstruction errors</td><td>3.7 × 10<sup>-2</sup></td><td>3 × 10<sup>-2</sup></td></tr><tr><td>Dual reconstruction errors</td><td>2.3 × 10<sup>-2</sup></td><td>1.5 × 10<sup>-2</sup></td></tr></tbody></table>

<caption>Table 1: Performance comparison of non-sparse and sparse methods. Mean of the quantities (time and reconstruction errors) observed on the validation set are given.</caption>

# 6 Conclusions and Perspectives

*Conclusions*: The key idea of using over-complete dictionaries for contact mechanics problems is introduced, whose necessity is justified using the known limitations of current reduced methods. Based on the concepts of sparse regressions, an algorithm is proposed to solve dictionary-based approximations of contact problems. Approximations based on large dictionaries were shown to be limit the inseparability issues to a certain extent. Performance of the proposed method seems to be correlated to the contact position variability, with satisfactory performance for small variabilities.

*Perspectives*: The primary bottleneck in the proposed method is the evaluation of non-linear operators that describe the inequality constraint i.e. the non-penetration condition. Current methods also run into issues caused by local properties of contact, as contact pairs change significantly for problems with large contact variations. Moreover, non-linear methods to interpolate the solution manifold should be explored.

# Acknowledgements

The research has been carried out as part of a doctoral thesis funded by IRT Jules Verne under the PERFORM programme. All computations were done on **python** using open source packages, namely NumPy [44] and SciPy [45] for scientific computations and Matplotlib [46] for creating
<table><caption>(a) Hertz problem with 30 snapshots in training set.</caption><thead><tr><th></th><th><em>Non-sparse</em><br>Active set with CPG<br>tolerance = 1 &times; 10<sup>-2</sup></th><th><em>Sparse</em><br>Greedy active set</th></tr></thead><tbody><tr><td>Primal RB size</td><td>22</td><td>22</td></tr><tr><td>Dual RB/Dictionary size</td><td>21</td><td>30</td></tr><tr><td>Time (seconds)</td><td>1.7 &times; 10<sup>-2</sup>s</td><td>4 &times; 10<sup>-3</sup> s</td></tr><tr><td>Primal reconstruction errors</td><td>1.4 &times; 10<sup>-3</sup></td><td>2 &times; 10<sup>-3</sup></td></tr><tr><td>Dual reconstruction errors</td><td>5 &times; 10<sup>-2</sup></td><td>5 &times; 10<sup>-2</sup></td></tr></tbody></table>

<table><caption>(b) Ironing problem with 129 snapshots in training set.</caption><thead><tr><th></th><th><em>Non-sparse</em><br>Active set with CPG<br>tolerance = 5 &times; 10<sup>-2</sup></th><th><em>Sparse</em><br>Greedy active set</th></tr></thead><tbody><tr><td>Primal RB size</td><td>99</td><td>99</td></tr><tr><td>Dual RB/Dictionary size</td><td>52</td><td>129</td></tr><tr><td>Time (seconds)</td><td>0.08s</td><td>0.02 s</td></tr><tr><td>Primal reconstruction errors</td><td>3.7 &times; 10<sup>-2</sup></td><td>3 &times; 10<sup>-2</sup></td></tr><tr><td>Dual reconstruction errors</td><td>2.3 &times; 10<sup>-2</sup></td><td>1.5 &times; 10<sup>-2</sup></td></tr></tbody></table>

<caption>Table 1: Performance comparison of non-sparse and sparse methods. Mean of the quantities (time and reconstruction errors) observed on the validation set are given.</caption>

# 6 Conclusions and Perspectives

*Conclusions*: The key idea of using over-complete dictionaries for contact mechanics problems is introduced, whose necessity is justified using the known limitations of current reduced methods. Based on the concepts of sparse regressions, an algorithm is proposed to solve dictionary-based approximations of contact problems. Approximations based on large dictionaries were shown to be limit the inseparability issues to a certain extent. Performance of the proposed method seems to be correlated to the contact position variability, with satisfactory performance for small variabilities.

*Perspectives*: The primary bottleneck in the proposed method is the evaluation of non-linear operators that describe the inequality constraint i.e. the non-penetration condition. Current methods also run into issues caused by local properties of contact, as contact pairs change significantly for problems with large contact variations. Moreover, non-linear methods to interpolate the solution manifold should be explored.

# Acknowledgements

The research has been carried out as part of a doctoral thesis funded by IRT Jules Verne under the PERFORM programme. All computations were done on **python** using open source packages, namely NumPy [44] and SciPy [45] for scientific computations and Matplotlib [46] for creating
References

[1] B. Li, P. Li, R. Zhou, X.-Q. Feng, and K. Zhou. Contact mechanics in tribological and contact damage-related problems: A review. *Tribology International*, 171:107534, 2022. ISSN 0301-679X. doi:10.1016/j.triboint.2022.107534.

[2] K. Brown, S. Attaway, S. Plimpton, and B. Hendrickson. Parallel strategies for crash and impact simulations. *Computer Methods in Applied Mechanics and Engineering*, 184(2):375-390, 2000. ISSN 0045-7825. doi:10.1016/S0045-7825(99)00235-2.

[3] S. T. Button. *Tribology in Metal Forming Processes*, pages 103–120. Springer Berlin Heidelberg, Berlin, Heidelberg, 2013. ISBN 978-3-642-31683-8. doi:10.1007/978-3-642-31683-8_3.

[4] G. A. Ateshian, C. R. Henak, and J. A. Weiss. Toward patient-specific articular contact mechanics. *Journal of Biomechanics*, 48(5):779–786, 2015. ISSN 0021-9290. doi:10.1016/j.jbiomech.2014.12.020.

[5] P. Wriggers and J. C. Simo. A note on tangent stiffness for fully nonlinear contact problems. *Communications in Applied Numerical Methods*, 1(5):199–203, 1985. doi:10.1002/cnm.1630010503.

[6] J. Hallquist, G. Goudreau, and D. Benson. Sliding interfaces with contact-impact in large-scale lagrangian computations. *Computer Methods in Applied Mechanics and Engineering*, 51(1):107–137, 1985. doi:10.1016/0045-7825(85)90030-1.

[7] K. A. Fischer and P. Wriggers. Mortar based frictional contact formulation for higher order interpolations using the moving friction cone. *Computer Methods in Applied Mechanics and Engineering*, 195(37-40):5020–5036, 2006. doi:10.1016/j.cma.2005.09.025.

[8] P. Wriggers. *Computational Contact Mechanics*. Springer Berlin Heidelberg, 2006. ISBN 978-3-540-32608-3. doi:10.1007/978-3-540-32609-0.

[9] F. Chouly, M. Fabre, P. Hild, R. Mlika, J. Pousin, and Y. Renard. An Overview of Recent Results on Nitsche's Method for Contact Problems. pages 93–141. 2017. doi:10.1007/978-3-319-71431-8_4.

[10] V. A. Yastrebov. *Numerical Methods in Contact Mechanics*. John Wiley & Sons, Inc., Hoboken, NJ USA, feb 2013. ISBN 9781118647974. doi:10.1002/9781118647974.

[11] U. Meier, O. López, C. Monserrat, M. Juan, and M. Alcañiz. Real-time deformable models for surgery simulation: a survey. *Computer Methods and Programs in Biomedicine*, 77 (3):183–197, 2005. ISSN 0169-2607. doi:https://doi.org/10.1016/j.cmpb.2004.11.002. URL https://www.sciencedirect.com/science/article/pii/S0169260704002093.

[12] X. Tan, W. Chen, T. Zou, J. Yang, and B. Du. Real-time prediction of mechanical behaviors of underwater shield tunnel structure using machine learning method based on structural health monitoring data. *Journal of Rock Mechanics and Geotechnical Engineering*, 15(4): 886–895, 2023. ISSN 1674-7755. doi:https://doi.org/10.1016/j.jrmge.2022.06.015.

[13] B. Haasdonk, J. Salomon, and B. Wohlmuth. A Reduced Basis Method for Parametrized Variational Inequalities. *SIAM J. Numer. Anal.*, 50(5):2656–2676, jan 2012. doi:10.1137/110835372.
plots. Meshes were created using Gmsh [47].

## References

[1] B. Li, P. Li, R. Zhou, X.-Q. Feng, and K. Zhou. Contact mechanics in tribological and contact damage-related problems: A review. *Tribology International*, 171:107534, 2022. ISSN 0301-679X. doi:10.1016/j.triboint.2022.107534.

[2] K. Brown, S. Attaway, S. Plimpton, and B. Hendrickson. Parallel strategies for crash and impact simulations. *Computer Methods in Applied Mechanics and Engineering*, 184(2):375-390, 2000. ISSN 0045-7825. doi:10.1016/S0045-7825(99)00235-2.

[3] S. T. Button. *Tribology in Metal Forming Processes*, pages 103-120. Springer Berlin Heidelberg, Berlin, Heidelberg, 2013. ISBN 978-3-642-31683-8. doi:10.1007/978-3-642-31683-8_3.

[4] G. A. Ateshian, C. R. Henak, and J. A. Weiss. Toward patient-specific articular contact mechanics. *Journal of Biomechanics*, 48(5):779-786, 2015. ISSN 0021-9290. doi:10.1016/j.jbiomech.2014.12.020.

[5] P. Wriggers and J. C. Simo. A note on tangent stiffness for fully nonlinear contact problems. *Communications in Applied Numerical Methods*, 1(5):199-203, 1985. doi:10.1002/cnm.1630010503.

[6] J. Hallquist, G. Goudreau, and D. Benson. Sliding interfaces with contact-impact in large-scale lagrangian computations. *Computer Methods in Applied Mechanics and Engineering*, 51(1):107-137, 1985. doi:10.1016/0045-7825(85)90030-1.

[7] K. A. Fischer and P. Wriggers. Mortar based frictional contact formulation for higher order interpolations using the moving friction cone. *Computer Methods in Applied Mechanics and Engineering*, 195(37-40):5020-5036, 2006. doi:10.1016/j.cma.2005.09.025.

[8] P. Wriggers. *Computational Contact Mechanics*. Springer Berlin Heidelberg, 2006. ISBN 978-3-540-32608-3. doi:10.1007/978-3-540-32609-0.

[9] F. Chouly, M. Fabre, P. Hild, R. Mlika, J. Pousin, and Y. Renard. An Overview of Recent Results on Nitsche's Method for Contact Problems. pages 93-141. 2017. doi:10.1007/978-3-319-71431-8_4.

[10] V. A. Yastrebov. *Numerical Methods in Contact Mechanics*. John Wiley & Sons, Inc., Hoboken, NJ USA, feb 2013. ISBN 9781118647974. doi:10.1002/9781118647974.

[11] U. Meier, O. López, C. Monserrat, M. Juan, and M. Alcañiz. Real-time deformable models for surgery simulation: a survey. *Computer Methods and Programs in Biomedicine*, 77 (3):183-197, 2005. ISSN 0169-2607. doi:https://doi.org/10.1016/j.cmpb.2004.11.002. URL https://www.sciencedirect.com/science/article/pii/S0169260704002093.

[12] X. Tan, W. Chen, T. Zou, J. Yang, and B. Du. Real-time prediction of mechanical behaviors of underwater shield tunnel structure using machine learning method based on structural health monitoring data. *Journal of Rock Mechanics and Geotechnical Engineering*, 15(4): 886-895, 2023. ISSN 1674-7755. doi:https://doi.org/10.1016/j.jrmge.2022.06.015.

[13] B. Haasdonk, J. Salomon, and B. Wohlmuth. A Reduced Basis Method for Parametrized Variational Inequalities. *SIAM J. Numer. Anal.*, 50(5):2656-2676, jan 2012. doi:10.1137/110835372.
[14] M. Balajewicz, D. Amsallem, and C. Farhat. Projection-based model reduction for contact problems. *Int. J. Numer. Methods Eng.*, 106(8):644–663, may 2016. doi:10.1002/nme.5135.

[15] I. Niakh, G. Drouet, V. Ehrlacher, and A. Ern. Stable model reduction for linear variational inequalities with parameter-dependent constraints. *ESAIM: Mathematical Modelling and Numerical Analysis*, 57(1):167–189, Jan. 2023. ISSN 2822-7840, 2804-7214. doi:10.1051/m2an/2022077.

[16] A. Benaceur, A. Ern, and V. Ehrlacher. A reduced basis method for parametrized variational inequalities applied to contact mechanics. *Int. J. Numer. Methods Eng.*, 121(6):1170–1197, mar 2020. doi:10.1002/nme.6261.

[17] J. Fauque, I. Ramière, and D. Ryckelynck. Hybrid hyper-reduced modeling for contact mechanics problems. *Int. J. Numer. Methods Eng.*, 115(1):117–139, Jul 2018. doi:10.1002/nme.5798.

[18] S. Le Berre, I. Ramière, J. Fauque, and D. Ryckelynck. Condition Number and Clustering-Based Efficiency Improvement of Reduced-Order Solvers for Contact Problems Using Lagrange Multipliers. *Mathematics*, 10(9):1495, Apr. 2022. ISSN 2227-7390. doi:10.3390/math10091495.

[19] E. Bader, Z. Zhang, and K. Veroy. An empirical interpolation approach to reduced basis approximations for variational inequalities. *Math. Comput. Model. Dyn. Syst.*, 22(4):345–361, 2016. doi:10.1080/13873954.2016.1198388.

[20] I. Niakh, G. Drouet, V. Ehrlacher, and A. Ern. A reduced basis method for frictional contact problems formulated with Nitsche's method. 2023. doi:10.48550/ARXIV.2307.11541. Publisher: arXiv Version Number: 1.

[21] D. Manvelyan, B. Simeon, and U. Wever. An efficient model order reduction scheme for dynamic contact in linear elasticity. *Computational Mechanics*, 68(6):1283–1295, dec 2021. doi:10.1007/s00466-021-02068-4.

[22] P. Mulye, J. Hemmer, L. Morançay, C. Binetruy, A. Leygue, S. Comas-Cardona, P. Pichon, and D. Guillon. Numerical modeling of interply adhesion in composite forming of viscous discontinuous thermoplastic prepregs. *Composites Part B: Engineering*, 191:107953, 2020. ISSN 1359-8368. doi:10.1016/j.compositesb.2020.107953.

[23] K. S. Kollepara, J. M. Navarro-Jiménez, Y. Le Guennec, L. Silva, and J. V. Aguado. On the limitations of low-rank approximations in contact mechanics problems. *International Journal for Numerical Methods in Engineering*, page nme.7118, Sept. 2022. ISSN 0029-5981, 1097-0207. doi:10.1002/nme.7118.

[24] R. H. Davies. *Learning shape: optimal models for analysing natural variability*. PhD thesis, University of Manchester, 2002. URL http://www.bmva.org/theses/2002/2002-davies.pdf.

[25] D. Amsallem, M. J. Zahr, and C. Farhat. Nonlinear model order reduction based on local reduced-order bases. *International Journal for Numerical Methods in Engineering*, 92(10): 891–916, 2012. doi:10.1002/nme.4371.

[26] D. Amsallem, J. Cortial, K. Carlberg, and C. Farhat. A method for interpolating on manifolds structural dynamics reduced-order models. *International Journal for Numerical Methods in Engineering*, 80(9):1241–1258, 2009. doi:10.1002/nme.2681.
[14] M. Balajewicz, D. Amsallem, and C. Farhat. Projection-based model reduction for contact problems. *Int. J. Numer. Methods Eng.*, 106(8):644–663, may 2016. doi:10.1002/nme.5135.

[15] I. Niakh, G. Drouet, V. Ehrlacher, and A. Ern. Stable model reduction for linear variational inequalities with parameter-dependent constraints. *ESAIM: Mathematical Modelling and Numerical Analysis*, 57(1):167–189, Jan. 2023. ISSN 2822-7840, 2804-7214. doi:10.1051/m2an/2022077.

[16] A. Benaceur, A. Ern, and V. Ehrlacher. A reduced basis method for parametrized variational inequalities applied to contact mechanics. *Int. J. Numer. Methods Eng.*, 121(6):1170–1197, mar 2020. doi:10.1002/nme.6261.

[17] J. Fauque, I. Ramière, and D. Ryckelynck. Hybrid hyper-reduced modeling for contact mechanics problems. *Int. J. Numer. Methods Eng.*, 115(1):117–139, Jul 2018. doi:10.1002/nme.5798.

[18] S. Le Berre, I. Ramière, J. Fauque, and D. Ryckelynck. Condition Number and Clustering-Based Efficiency Improvement of Reduced-Order Solvers for Contact Problems Using Lagrange Multipliers. *Mathematics*, 10(9):1495, Apr. 2022. ISSN 2227-7390. doi:10.3390/math10091495.

[19] E. Bader, Z. Zhang, and K. Veroy. An empirical interpolation approach to reduced basis approximations for variational inequalities. *Math. Comput. Model. Dyn. Syst.*, 22(4):345–361, 2016. doi:10.1080/13873954.2016.1198388.

[20] I. Niakh, G. Drouet, V. Ehrlacher, and A. Ern. A reduced basis method for frictional contact problems formulated with Nitsche's method. 2023. doi:10.48550/ARXIV.2307.11541. Publisher: arXiv Version Number: 1.

[21] D. Manvelyan, B. Simeon, and U. Wever. An efficient model order reduction scheme for dynamic contact in linear elasticity. *Computational Mechanics*, 68(6):1283–1295, dec 2021. doi:10.1007/s00466-021-02068-4.

[22] P. Mulye, J. Hemmer, L. Morançay, C. Binetruy, A. Leygue, S. Comas-Cardona, P. Pichon, and D. Guillon. Numerical modeling of interply adhesion in composite forming of viscous discontinuous thermoplastic prepregs. *Composites Part B: Engineering*, 191:107953, 2020. ISSN 1359-8368. doi:10.1016/j.compositesb.2020.107953.

[23] K. S. Kollepara, J. M. Navarro-Jiménez, Y. Le Guennec, L. Silva, and J. V. Aguado. On the limitations of low-rank approximations in contact mechanics problems. *International Journal for Numerical Methods in Engineering*, page nme.7118, Sept. 2022. ISSN 0029-5981, 1097-0207. doi:10.1002/nme.7118.

[24] R. H. Davies. *Learning shape: optimal models for analysing natural variability*. PhD thesis, University of Manchester, 2002. URL http://www.bmva.org/theses/2002/2002-davies.pdf.

[25] D. Amsallem, M. J. Zahr, and C. Farhat. Nonlinear model order reduction based on local reduced-order bases. *International Journal for Numerical Methods in Engineering*, 92(10): 891–916, 2012. doi:10.1002/nme.4371.

[26] D. Amsallem, J. Cortial, K. Carlberg, and C. Farhat. A method for interpolating on manifolds structural dynamics reduced-order models. *International Journal for Numerical Methods in Engineering*, 80(9):1241–1258, 2009. doi:10.1002/nme.2681.
[27] M. Aharon, M. Elad, and A. Bruckstein. K-svd: An algorithm for designing overcomplete dictionaries for sparse representation. *IEEE Transactions on Signal Processing*, 54(11): 4311–4322, 2006. doi:10.1109/tsp.2006.881199.

[28] R. Rubinstein, M. Zibulevsky, and M. Elad. Efficient Implementation of the K-SVD Algorithm Using Batch Orthogonal Matching Pursuit. *CS Techion*, (January 2008), 2014.

[29] A. Balabanov, O.and Nouy. Randomized linear algebra for model reduction—part II: minimal residual methods and dictionary-based approximation. *Advances in Computational Mathematics*, 47(2):26, apr 2021. doi:10.1007/s10444-020-09836-5.

[30] N. Lauzeral. *Reduced order and sparse representations for patient-specific modeling in computational surgery*. PhD thesis, École centrale de Nantes, dec 2019. URL http://www.theses.fr/2019ECDN0062.

[31] E. Simoncelli, W. Freeman, E. Adelson, and D. Heeger. Shiftable multiscale transforms. *IEEE Transactions on Information Theory*, 38(2):587–607, 1992. doi:10.1109/18.119725.

[32] E. J. Candès and D. L. Donoho. Recovering edges in ill-posed inverse problems: optimality of curvelet frames. *The Annals of Statistics*, 30(3):784–842, 2002. doi:10.1214/aos/1028674842. URL https://doi.org/10.1214/aos/1028674842.

[33] M. Elad. *Sparse and Redundant Representations*. Springer New York, New York, NY, 2010. ISBN 978-1-4419-7010-7. doi:10.1007/978-1-4419-7011-4.

[34] S. L. Brunton, J. L. Proctor, and J. N. Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. *Proceedings of the National Academy of Sciences*, 113(15):3932–3937, apr 2016. doi:10.1073/pnas.1517384113.

[35] J. Mairal, F. Bach, and J. Ponce. Sparse modeling for image and vision processing. *Foundations and Trends in Computer Graphics and Vision*, 8(2-3):85–283, 2014. doi:10.1561/0600000058.

[36] R. Tibshirani. Regression Shrinkage and Selection Via the Lasso. *Journal of the Royal Statistical Society: Series B (Methodological)*, 58(1):267–288, jan 1996. doi:10.1111/j.2517-6161.1996.tb02080.x.

[37] E. Candes and T. Tao. The Dantzig selector: Statistical estimation when p is much larger than n. *The Annals of Statistics*, 35(6):2313–2351, dec 2007. doi:10.1214/009053606000001523.

[38] Y. Pati, R. Rezaiifar, and P. Krishnaprasad. Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition. In *Proceedings of 27th Asilomar Conference on Signals, Systems and Computers*, pages 40–44 vol.1, 1993. doi:10.1109/ACSSC.1993.342465.

[39] I. Gorodnitsky and B. Rao. Sparse signal reconstruction from limited data using FOCUSS: a re-weighted minimum norm algorithm. *IEEE Transactions on Signal Processing*, 45(3): 600–616, mar 1997. doi:10.1109/78.558475.

[40] W.-X. Wang, R. Yang, Y.-C. Lai, V. Kovanis, and C. Grebogi. Predicting catastrophes in nonlinear dynamical systems by compressive sensing. *Phys. Rev. Lett.*, 106:154101, Apr 2011. doi:10.1103/PhysRevLett.106.154101.
[27] M. Aharon, M. Elad, and A. Bruckstein. K-svd: An algorithm for designing overcomplete dictionaries for sparse representation. *IEEE Transactions on Signal Processing*, 54(11): 4311–4322, 2006. doi:10.1109/tsp.2006.881199.

[28] R. Rubinstein, M. Zibulevsky, and M. Elad. Efficient Implementation of the K-SVD Algorithm Using Batch Orthogonal Matching Pursuit. *CS Techion*, (January 2008), 2014.

[29] A. Balabanov, O.and Nouy. Randomized linear algebra for model reduction—part II: minimal residual methods and dictionary-based approximation. *Advances in Computational Mathematics*, 47(2):26, apr 2021. doi:10.1007/s10444-020-09836-5.

[30] N. Lauzeral. *Reduced order and sparse representations for patient-specific modeling in computational surgery*. PhD thesis, École centrale de Nantes, dec 2019. URL http://www.theses.fr/2019ECDN0062.

[31] E. Simoncelli, W. Freeman, E. Adelson, and D. Heeger. Shiftable multiscale transforms. *IEEE Transactions on Information Theory*, 38(2):587–607, 1992. doi:10.1109/18.119725.

[32] E. J. Candès and D. L. Donoho. Recovering edges in ill-posed inverse problems: optimality of curvelet frames. *The Annals of Statistics*, 30(3):784–842, 2002. doi:10.1214/aos/1028674842. URL https://doi.org/10.1214/aos/1028674842.

[33] M. Elad. *Sparse and Redundant Representations*. Springer New York, New York, NY, 2010. ISBN 978-1-4419-7010-7. doi:10.1007/978-1-4419-7011-4.

[34] S. L. Brunton, J. L. Proctor, and J. N. Kutz. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. *Proceedings of the National Academy of Sciences*, 113(15):3932–3937, apr 2016. doi:10.1073/pnas.1517384113.

[35] J. Mairal, F. Bach, and J. Ponce. Sparse modeling for image and vision processing. *Foundations and Trends in Computer Graphics and Vision*, 8(2-3):85–283, 2014. doi:10.1561/0600000058.

[36] R. Tibshirani. Regression Shrinkage and Selection Via the Lasso. *Journal of the Royal Statistical Society: Series B (Methodological)*, 58(1):267–288, jan 1996. doi:10.1111/j.2517-6161.1996.tb02080.x.

[37] E. Candes and T. Tao. The Dantzig selector: Statistical estimation when p is much larger than n. *The Annals of Statistics*, 35(6):2313–2351, dec 2007. doi:10.1214/009053606000001523.

[38] Y. Pati, R. Rezaiifar, and P. Krishnaprasad. Orthogonal matching pursuit: recursive function approximation with applications to wavelet decomposition. In *Proceedings of 27th Asilomar Conference on Signals, Systems and Computers*, pages 40–44 vol.1, 1993. doi:10.1109/ACSSC.1993.342465.

[39] I. Gorodnitsky and B. Rao. Sparse signal reconstruction from limited data using FOCUSS: a re-weighted minimum norm algorithm. *IEEE Transactions on Signal Processing*, 45(3): 600–616, mar 1997. doi:10.1109/78.558475.

[40] W.-X. Wang, R. Yang, Y.-C. Lai, V. Kovanis, and C. Grebogi. Predicting catastrophes in nonlinear dynamical systems by compressive sensing. *Phys. Rev. Lett.*, 106:154101, Apr 2011. doi:10.1103/PhysRevLett.106.154101.
[41] H. Schaeffer, R. Caflisch, C. D. Hauck, and S. Osher. Sparse dynamics for partial differential equations. *Proceedings of the National Academy of Sciences*, 110(17):6634-6639, 2013. doi:10.1073/pnas.1302752110.

[42] N. Halko, P. G. Martinsson, and J. A. Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. *SIAM Review*, 53 (2):217-288, 2011. doi:10.1137/090771806.

[43] Y. Le Guennec, J. Brunet, F. Daim, M. Chau, and Y. Tourbier. A parametric and non-intrusive reduced order model of car crash simulation. *Computer Methods in Applied Mechanics and Engineering*, 338:186-207, 2018. doi:10.1016/j.cma.2018.03.005.

[44] C. R. Harris, K. J. Millman, S. J. van der Walt, and et. al. Array programming with NumPy. *Nature*, 585(7825):357-362, Sept. 2020. doi:10.1038/s41586-020-2649-2.

[45] P. Virtanen, R. Gommers, T. Oliphant, and et. al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nature Methods*, 17:261-272, 2020. doi:10.1038/s41592-019-0686-2.

[46] J. D. Hunter. Matplotlib: A 2d graphics environment. *Computing in Science & Engineering*, 9(3):90-95, 2007. doi:10.1109/MCSE.2007.55.

[47] C. Geuzaine and J.-F. Remacle. Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities. *International Journal for Numerical Methods in Engineering*, 79:1309-1331, 2008.

[48] S. Mallat and Z. Zhang. Matching pursuits with time-frequency dictionaries. *IEEE Transactions on Signal Processing*, 41(12):3397-3415, 1993. doi:10.1109/78.258082.

[49] M. Manucci, J. V. Aguado, and D. Borzacchiello. Sparse data-driven quadrature rules via $l^p$-quasi-norm minimization. *Computational Methods in Applied Mathematics*, 2022. doi:10.1515/cmam-2021-0131.

[50] C. L. Lawson and R. J. Hanson. *Solving Least Squares Problems*. Society for Industrial and Applied Mathematics, jan 1995. ISBN 978-0-89871-356-5. doi:10.1137/1.9781611971217.

[51] K. Veroy and A. T. Patera. Certified real-time solution of the parametrized steady incompressible navier-stokes equations: rigorous reduced-basis a posteriori error bounds. *International Journal for Numerical Methods in Fluids*, 47(8-9):773-788, 2005. doi:10.1002/fld.867.

[52] E. Liberge and A. Hamdouni. Reduced order modelling method via proper orthogonal decomposition (pod) for flow around an oscillating cylinder. *Journal of Fluids and Structures*, 26(2):292-311, 2010. ISSN 0889-9746. doi:10.1016/j.jfluidstructs.2009.10.006.
[41] H. Schaeffer, R. Caflisch, C. D. Hauck, and S. Osher. Sparse dynamics for partial differential equations. *Proceedings of the National Academy of Sciences*, 110(17):6634-6639, 2013. doi:10.1073/pnas.1302752110.

[42] N. Halko, P. G. Martinsson, and J. A. Tropp. Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions. *SIAM Review*, 53 (2):217-288, 2011. doi:10.1137/090771806.

[43] Y. Le Guennec, J. Brunet, F. Daim, M. Chau, and Y. Tourbier. A parametric and non-intrusive reduced order model of car crash simulation. *Computer Methods in Applied Mechanics and Engineering*, 338:186-207, 2018. doi:10.1016/j.cma.2018.03.005.

[44] C. R. Harris, K. J. Millman, S. J. van der Walt, and et. al. Array programming with NumPy. *Nature*, 585(7825):357-362, Sept. 2020. doi:10.1038/s41586-020-2649-2.

[45] P. Virtanen, R. Gommers, T. Oliphant, and et. al. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. *Nature Methods*, 17:261-272, 2020. doi:10.1038/s41592-019-0686-2.

[46] J. D. Hunter. Matplotlib: A 2d graphics environment. *Computing in Science & Engineering*, 9(3):90-95, 2007. doi:10.1109/MCSE.2007.55.

[47] C. Geuzaine and J.-F. Remacle. Gmsh: A 3-d finite element mesh generator with built-in pre-and post-processing facilities. *International Journal for Numerical Methods in Engineering*, 79:1309-1331, 2008.

[48] S. Mallat and Z. Zhang. Matching pursuits with time-frequency dictionaries. *IEEE Transactions on Signal Processing*, 41(12):3397-3415, 1993. doi:10.1109/78.258082.

[49] M. Manucci, J. V. Aguado, and D. Borzacchiello. Sparse data-driven quadrature rules via $\ell^p$-quasi-norm minimization. *Computational Methods in Applied Mathematics*, 2022. doi:10.1515/cmam-2021-0131.

[50] C. L. Lawson and R. J. Hanson. *Solving Least Squares Problems*. Society for Industrial and Applied Mathematics, jan 1995. ISBN 978-0-89871-356-5. doi:10.1137/1.9781611971217.

[51] K. Veroy and A. T. Patera. Certified real-time solution of the parametrized steady incompressible navier-stokes equations: rigorous reduced-basis a posteriori error bounds. *International Journal for Numerical Methods in Fluids*, 47(8-9):773-788, 2005. doi:10.1002/fld.867.

[52] E. Liberge and A. Hamdouni. Reduced order modelling method via proper orthogonal decomposition (pod) for flow around an oscillating cylinder. *Journal of Fluids and Structures*, 26(2):292-311, 2010. ISSN 0889-9746. doi:10.1016/j.jfluidstructs.2009.10.006.
# A Sparse methods

As discussed in Section 3, solving the problem (3) with $l_0$-norm is NP-hard. Thus, practically usable sparse regression methods have been devised by the community, such as, sparsity inducing surrogate norms in place of the $l_0$-norm, or greedy methods. Some examples include LASSO [36] and Dantzig Selector [37] that uses $l_0$-norm. The methods used in this work, OMP and FOCUSS are detailed below:

## A.1 Orthogonal Matching Pursuit

The OMP method [38] aims to approximate the solution of (3) with $p = 0$ using a greedy technique. At each greedy step $k$, the OMP algorithm searches for the column of dictionary $D$ that is the least orthogonal to the current residual $(x - D\alpha^{k-1})$ and adds to the set of indices previously selected. This process is carried out until the residual is below the tolerance $\epsilon$. This is done by projecting the residual vector on each column and selecting the largest projection.

$$ \mathcal{I} \leftarrow \mathcal{I} \cup \arg \max_j |\mathbf{D}[\cdot, j]^T (x - D\alpha^{k-1})| $$

where $\mathcal{I}$ is the current set of selected indices and columns of dictionary $D$ are normalized. Then the new set of coefficients $\alpha^k$ are computed using the least square solution:

$$ \alpha_{\mathcal{I}}^k = (\mathbf{D}_{\mathcal{I}}^T \mathbf{D}_{\mathcal{I}})^{-1} \mathbf{D}_{\mathcal{I}}^T x $$

with $\alpha_{\mathcal{I}} = \alpha[\mathcal{I}]$ and $\mathbf{D}_{\mathcal{I}} = \mathbf{D}[\cdot, \mathcal{I}]$.

Update of $\mathcal{I}$ is locally optimal, but the above least-square update of $\alpha^k$, the OMP solution is optimal w.r.t to the currently selected subset $\mathbf{D}_{\mathcal{I}}$ of the dictionary. This is in contrast to the predecessor of OMP, the Matching Pursuit [48], which used locally optimal updates for $\alpha^k$ at the current iteration.

## A.2 FOCUSS

The FOCUSS method [39] finds a sparse approximation using an iterative process that starts from a fully dense solution and progresses towards “localized energy solutions”. The FOCUSS method is broadly based on the following ideas:

* For a full rank dictionary of size, say $m \times n$ ($m < n$), the system $D\alpha = x$ is under-determined. For such a system, the closed form solution for (3) with $p = 2$ is given by $\alpha = D^\dagger x$, where $D^\dagger$ is the Moore-Penrose pseudo-inverse. Note that this solution is not sparse in general, as $l_2$-norm does not have sparsity-inducing properties.

* A modified minimization problem can be defined by replacing the $||\alpha||_2$ with a weighted norm $||W^{-1}\alpha||_2$, or more generically $||W^\dagger\alpha||_2$ if $W$ is singular. The authors of [39] argue that by changing $W$, every possible solution of the under-determined system $D\alpha = x$ can be obtained.

* The basis of the algorithm lies in iterating towards a sparse solution using a weight $W^k$ which induces sparsity. For the iteration $k$, the weight is chosen as $W_k = \text{diag}(\alpha^{k-1})$.
# A Sparse methods

As discussed in Section 3, solving the problem (3) with $l_0$-norm is NP-hard. Thus, practically usable sparse regression methods have been devised by the community, such as, sparsity inducing surrogate norms in place of the $l_0$-norm, or greedy methods. Some examples include LASSO [36] and Dantzig Selector [37] that uses $l_0$-norm. The methods used in this work, OMP and FOCUSS are detailed below:

## A.1 Orthogonal Matching Pursuit

The OMP method [38] aims to approximate the solution of (3) with $p = 0$ using a greedy technique. At each greedy step $k$, the OMP algorithm searches for the column of dictionary $\mathbf{D}$ that is the least orthogonal to the current residual $(\mathbf{x} - \mathbf{D}\alpha^{k-1})$ and adds to the set of indices previously selected. This process is carried out until the residual is below the tolerance $\epsilon$. This is done by projecting the residual vector on each column and selecting the largest projection.

$$ \mathcal{I} \leftarrow \mathcal{I} \cup \arg \max_{j} \left| \mathbf{D}[\cdot, j]^T (\mathbf{x} - \mathbf{D}\alpha^{k-1}) \right| $$

where $\mathcal{I}$ is the current set of selected indices and columns of dictionary $\mathbf{D}$ are normalized. Then the new set of coefficients $\alpha^k$ are computed using the least square solution:

$$ \alpha_{\mathcal{I}}^{k} = (\mathbf{D}_{\mathcal{I}}^{T} \mathbf{D}_{\mathcal{I}})^{-1} \mathbf{D}_{\mathcal{I}}^{T} \mathbf{x} $$

with $\alpha_{\mathcal{I}} = \alpha[\mathcal{I}]$ and $\mathbf{D}_{\mathcal{I}} = \mathbf{D}[\cdot, \mathcal{I}]$.

Update of $\mathcal{I}$ is locally optimal, but the above least-square update of $\alpha^k$, the OMP solution is optimal w.r.t to the currently selected subset $\mathbf{D}_{\mathcal{I}}$ of the dictionary. This is in contrast to the predecessor of OMP, the Matching Pursuit [48], which used locally optimal updates for $\alpha^k$ at the current iteration.

## A.2 FOCUSS

The FOCUSS method [39] finds a sparse approximation using an iterative process that starts from a fully dense solution and progresses towards “localized energy solutions”. The FOCUSS method is broadly based on the following ideas:

* For a full rank dictionary of size, say $m \times n$ ($m < n$), the system $\mathbf{D}\alpha = \mathbf{x}$ is under-determined. For such a system, the closed form solution for (3) with $p = 2$ is given by $\alpha = \mathbf{D}^{\dagger}\mathbf{x}$, where $\mathbf{D}^{\dagger}$ is the Moore-Penrose pseudo-inverse. Note that this solution is not sparse in general, as $l_2$-norm does not have sparsity-inducing properties.

* A modified minimization problem can be defined by replacing the $||\alpha||_2$ with a weighted norm $||\mathbf{W}^{-1}\alpha||_2$, or more generically $||\mathbf{W}^{\dagger}\alpha||_2$ if $\mathbf{W}$ is singular. The authors of [39] argue that by changing $\mathbf{W}$, every possible solution of the under-determined system $\mathbf{D}\alpha = \mathbf{x}$ can be obtained.

* The basis of the algorithm lies in iterating towards a sparse solution using a weight $\mathbf{W}^k$ which induces sparsity. For the iteration $k$, the weight is chosen as $\mathbf{W}_k = \text{diag}(\alpha^{k-1})$.
The trick lies in the fact that the algorithm, effectively, minimizes the following weighted norm:

$$ \|\mathbf{W}_k^\dagger \alpha\|_2^2 = \sum_i \left( \frac{\alpha_i^k}{\alpha_i^{k-1}} \right)^2 $$

due to which smaller entries of $\alpha^{k-1}$ tend to diminish further.

On combining these ideas, the FOCUSS algorithm ends up with quite a simple implementation, given in Alg. 2. Usually, the algorithm is initialized with the $\alpha^0$ containing all non-zeros. In [39], the $||\alpha||_2$ minimizing solution $\alpha^0 = \mathbf{D}^\dagger \mathbf{x}$ was used for initialization.

**Algorithm 2 FOCUSS**

1: Inputs: $\mathbf{D}, \mathbf{x}, \alpha^0$

2: Initialize $k = 1$

3: **while** $\alpha$ not converge **do**

4: $\mathbf{W}^k = \text{diag}(\alpha^{k-1})$

5: $\alpha^k = \mathbf{W}^k (\mathbf{D}\mathbf{W}^k)^\dagger \mathbf{x}$

6: $k \leftarrow k + 1$

7: **end while**

The nnFOCUSS Algorithm: The FOCUSS algorithm can be modified to generate a non-negative solution [30, 49]. The nnFOCUSS algorithm works by computing an appropriate relaxation parameter for each iterative update that maintains non-negativity. After computation of new coefficients at Step 5 of the Alg. 2, the relaxation step of nnFOCUSS algorithm is performed as follows:

if $\min(\alpha^k) < 0$ then

$$ \Delta\alpha = -(\alpha^k - \alpha^{k-1})^{-} $$

$$ \alpha^k \leftarrow \alpha^{k-1} + \min \left( \frac{\alpha^{k-1}}{\Delta\alpha} \right) \Delta\alpha $$

▷ Element wise division

end if

The algorithm must be initialized using non-negative coefficients, $\alpha^0$, which is computed using the non-negative least squares [50] solution (which is not sparse, in general).

## B A non-penetrating convex hull approach with mono-lithic dictionaries

Another novel sparse approximation method was formulated, namely the convex-hull approximation approach. This section briefly addresses the formulation and preliminary results. The convex hull method differs from the greedy active-set method in several aspects. Before this method is formulated, a broad overview of these differences is provided in Table 2.

The inspiration for this approach comes from the reduced modelling of incompressible flow problems, where the velocity field is inherently divergence-free. Consequently, the reduced incompressible flow problem is unconstrained, as all candidates in velocity subspace naturally satisfy the constraints [51, 52]. It is not straightforward to extend the same idea to contact mechanics problems, because linear combinations of the penetration-free primal snapshots do not satisfy the non-penetration condition.
The trick lies in the fact that the algorithm, effectively, minimizes the following weighted norm:

$$ \|\mathbf{W}_k^\dagger \alpha\|_2^2 = \sum_i \left( \frac{\alpha_i^k}{\alpha_i^{k-1}} \right)^2 $$

due to which smaller entries of $\alpha^{k-1}$ tend to diminish further.

On combining these ideas, the FOCUSS algorithm ends up with quite a simple implementation, given in Alg. 2. Usually, the algorithm is initialized with the $\alpha^0$ containing all non-zeros. In [39], the $||\alpha||_2$ minimizing solution $\alpha^0 = \mathbf{D}^\dagger \mathbf{x}$ was used for initialization.

<caption>**Algorithm 2 FOCUSS**</caption>

1: Inputs: $\mathbf{D}, \mathbf{x}, \alpha^0$
2: Initialize $k = 1$
3: **while** $\alpha$ not converge **do**
4: &nbsp;&nbsp;&nbsp;&nbsp; $\mathbf{W}^k = \text{diag}(\alpha^{k-1})$
5: &nbsp;&nbsp;&nbsp;&nbsp; $\alpha^k = \mathbf{W}^k (\mathbf{D}\mathbf{W}^k)^\dagger \mathbf{x}$
6: &nbsp;&nbsp;&nbsp;&nbsp; $k \leftarrow k + 1$
7: **end while**

The nnFOCUSS Algorithm: The FOCUSS algorithm can be modified to generate a non-negative solution [30, 49]. The nnFOCUSS algorithm works by computing an appropriate relaxation parameter for each iterative update that maintains non-negativity. After computation of new coefficients at Step 5 of the Alg. 2, the relaxation step of nnFOCUSS algorithm is performed as follows:

if $\min(\alpha^k) < 0$ then
$$ \Delta\alpha = -(\alpha^k - \alpha^{k-1})^\top $$
$$ \alpha^k \leftarrow \alpha^{k-1} + \min\left(\frac{\alpha^{k-1}}{\Delta\alpha}\right) \Delta\alpha $$
▷ Element wise division
end if

The algorithm must be initialized using non-negative coefficients, $\alpha^0$, which is computed using the non-negative least squares [50] solution (which is not sparse, in general).

## B A non-penetrating convex hull approach with mono-lithic dictionaries

Another novel sparse approximation method was formulated, namely the convex-hull approximation approach. This section briefly addresses the formulation and preliminary results. The convex hull method differs from the greedy active-set method in several aspects. Before this method is formulated, a broad overview of these differences is provided in Table 2.

The inspiration for this approach comes from the reduced modelling of incompressible flow problems, where the velocity field is inherently divergence-free. Consequently, the reduced incompressible flow problem is unconstrained, as all candidates in velocity subspace naturally satisfy the constraints [51, 52]. It is not straightforward to extend the same idea to contact mechanics problems, because linear combinations of the penetration-free primal snapshots do not satisfy the non-penetration condition.
<table><caption>Table 2: Characteristics of the two dictionary-based approximation methods for contact mechanics problems</caption><thead><tr><th></th><th>Convex hull approximation</th><th>Greedy active-set</th></tr></thead><tbody><tr><td>Primal Basis type</td><td>Dictionary</td><td>Truncated low-rank</td></tr><tr><td>Dual Basis type</td><td>Dictionary</td><td>Dictionary</td></tr><tr><td>Monolithic dictionary</td><td>Yes</td><td>No</td></tr><tr><td>Iteration method</td><td>Fixed point</td><td>Fixed point</td></tr><tr><td>Constraint enforcement</td><td>Convex combinations of non-penetrating snapshots</td><td>Active-set</td></tr><tr><td>Sparsity induction</td><td>Non-negative FOCUSS</td><td>Greedy activation and deactivation of dual dofs</td></tr><tr><td>Additional assumptions</td><td>Feasible region is convex and parameter-independent</td><td>None</td></tr></tbody></table>

For simplicity, let us first consider a contact problem where the contact pairs are constant. The discrete inequality constraint can be expressed using matrices that do not change with the displacement field. The feasible region $\mathcal{K}$ in the displacement solution space can be defined as:

$$ \mathcal{K} := \{u \mid Cu - g \le 0\} \qquad (10) $$

<figure><img src="image_20.png" /><figcaption>Figure 15: Illustration of non-penetrating property of convex hull. The convex feasible region $\mathcal{K}$ is separated from non-feasible zone by the hyperplanes defined by $Cu = g$, illustrated schematically by segments. The convex hull of the displacement snapshots is a convex subset of low-dimensional subspace defined by the linear combination of snapshots. It is also a subset of the convex feasible zone $\mathcal{K}$.</figcaption></figure>

As the operators $\mathbf{C}$ and $\mathbf{g}$ are independent of displacement field $u$, it is straightforward to demonstrate the convexity of $\mathcal{K}$. Thus, a convex combination<fnref n="4" /> of non-penetrating snapshots would naturally satisfy the inequality constraint. The implication here is that all the candidates in convex hull<fnref n="5" /> of the training set snapshots satisfy the inequality constraint. This notion is visually shown in Fig. 15, where the convex hull of training set snapshots does not violate the

<footnote n="4">⁴A convex combination is a linear combination with non-negative coefficients that sum up to unity</footnote>

<footnote n="5">⁵Convex hull of a set of points (or snapshots) is the set containing all convex combinations of the points. In other words, it is the smallest convex set that envelopes the set of points.</footnote>
<table><caption>Table 2: Characteristics of the two dictionary-based approximation methods for contact mechanics problems</caption><thead><tr><th></th><th>Convex hull approximation</th><th>Greedy active-set</th></tr></thead><tbody><tr><td>Primal Basis type</td><td>Dictionary</td><td>Truncated low-rank</td></tr><tr><td>Dual Basis type</td><td>Dictionary</td><td>Dictionary</td></tr><tr><td>Monolithic dictionary</td><td>Yes</td><td>No</td></tr><tr><td>Iteration method</td><td>Fixed point</td><td>Fixed point</td></tr><tr><td>Constraint enforcement</td><td>Convex combinations of non-penetrating snapshots</td><td>Active-set</td></tr><tr><td>Sparsity induction</td><td>Non-negative FOCUSS</td><td>Greedy activation and deactivation of dual dofs</td></tr><tr><td>Additional assumptions</td><td>Feasible region is convex and parameter-independent</td><td>None</td></tr></tbody></table>

For simplicity, let us first consider a contact problem where the contact pairs are constant. The discrete inequality constraint can be expressed using matrices that do not change with the displacement field. The feasible region $\mathcal{K}$ in the displacement solution space can be defined as:

$$ \mathcal{K} := \{u \mid Cu - g \le 0\} \qquad (10) $$

<figure><img src="image_21.png" /><figcaption>Figure 15: Illustration of non-penetrating property of convex hull. The convex feasible region $\mathcal{K}$ is separated from non-feasible zone by the hyperplanes defined by $Cu = g$, illustrated schematically by segments. The convex hull of the displacement snapshots is a convex subset of low-dimensional subspace defined by the linear combination of snapshots. It is also a subset of the convex feasible zone $\mathcal{K}$.</figcaption></figure>

As the operators $\mathbf{C}$ and $\mathbf{g}$ are independent of displacement field $u$, it is straightforward to demonstrate the convexity of $\mathcal{K}$. Thus, a convex combination$^4$ of non-penetrating snapshots would naturally satisfy the inequality constraint. The implication here is that all the candidates in convex hull$^5$ of the training set snapshots satisfy the inequality constraint. This notion is visually shown in Fig. 15, where the convex hull of training set snapshots does not violate the

<footnote>$^4$A convex combination is a linear combination with non-negative coefficients that sum up to unity</footnote>

<footnote>$^5$Convex hull of a set of points (or snapshots) is the set containing all convex combinations of the points. In other words, it is the smallest convex set that envelopes the set of points.</footnote>
convex non-penetration constraint. This property of convexity can be exploited to build efficient reduced models, if the following hypothesis holds:

**Hypothesis 1** (Convex Subset hypothesis). Given a parametrized high-dimensional inequality-constrained problem with a convex feasible set $\mathcal{K}$, its solutions lie in a low-dimensional convex subset inside $\mathcal{K}$.

To explore the convex hull of the dictionary, the displacement field will be approximated using convex combinations of the primal dictionary $\mathbf{D}_p$:

$$ \begin{array}{l} u \approx \mathbf{D}_p \hat{u} \\ \text{s.t.} \quad 1^T \hat{u} = 1 \\ \phantom{\text{s.t.} \quad} \hat{u} \ge 0 \end{array} $$

As the non-penetration constraint will be satisfied naturally by all candidates in the convex hull of dictionary snapshots, a dictionary-based approximation can be computed using coefficients $\hat{u}$ that satisfy equilibrium equations. Therefore, the objective is to minimize the following residual:

$$ r(\alpha) = \mathbf{K}(\mu)\mathbf{D}_p\alpha + \mathbf{C}^T(\mu, u)\mathbf{D}_d\alpha - f(\mu) $$

Note that the residual has been defined using $\hat{u} = \hat{\lambda} = \alpha$, which is equivalent to using a monolithic dictionary, in which each column consists of a displacement snapshot stacked above the corresponding contact pressure snapshot. Usage of a monolithic dictionary not only reduces the number of unknowns but also makes the problem of minimizing the residual $||r(\alpha)||$ well-defined. Moreover, the residual can be projected on a low-dimensional subspace of the colsp($\mathbf{D}_p$), described by the matrix $\mathbf{B}$ containing its basis, while also adding sparsity constraints to the unknown $\alpha$, like in (5). Thus, the problem can be stated as:

$$ \min ||\alpha||_p \qquad (11a) $$

$$ \text{s.t.} \quad ||\mathbf{B}^T r(\alpha)|| < \varepsilon \qquad (11b) $$

$$ 1^T \alpha = 1 \qquad (11c) $$

$$ \alpha \ge 0 \qquad (11d) $$

The nnFOCUSS algorithm [30, 49], based on the FOCUSS algorithm [39], is a weighted-norm minimization technique that progressively induces sparsity to the solution(see Appendix A for details). The algorithm takes three inputs: the vector $\mathbf{x}$ that will be approximated, the dictionary $\mathbf{D}$ that serves as the approximation space, and a non-sparse non-negative initial guess $\alpha^0$ of approximation coefficients. It can be used to solve the convex hull approximation problem (11), by plugging the inputs defined as:

$$ \mathbf{D} = \begin{bmatrix} \mathbf{B}^T \mathbf{K}_{\text{mono}} \\ 1^T \end{bmatrix}, \quad \mathbf{x} = \begin{bmatrix} \mathbf{B}^T f(\mu) \\ 1 \end{bmatrix}, \quad \alpha^0 = \text{nnls}(\mathbf{D}, \mathbf{x}) $$

where $\mathbf{K}_{\text{mono}}$ is the operator for the monolithic residual $r(\alpha)$:

$$ \mathbf{K}_{\text{mono}} = \mathbf{K}(\mu)\mathbf{D}_p + \mathbf{C}^T(\mu, u)\mathbf{D}_d $$

nnls indicates the non-negative least squares. An open-source implementation of nnls provided in the SciPy library [45] is used.
convex non-penetration constraint. This property of convexity can be exploited to build efficient reduced models, if the following hypothesis holds:

**Hypothesis 1** (Convex Subset hypothesis). Given a parametrized high-dimensional inequality-constrained problem with a convex feasible set $\mathcal{K}$, its solutions lie in a low-dimensional convex subset inside $\mathcal{K}$.

To explore the convex hull of the dictionary, the displacement field will be approximated using convex combinations of the primal dictionary $\mathbf{D}_p$:

$$ 
\begin{array}{l}
u \approx \mathbf{D}_p \hat{u} \\
\text{s.t.} \quad 1^T \hat{u} = 1 \\
\phantom{\text{s.t.} \quad} \hat{u} \ge 0
\end{array}
$$

As the non-penetration constraint will be satisfied naturally by all candidates in the convex hull of dictionary snapshots, a dictionary-based approximation can be computed using coefficients $\hat{u}$ that satisfy equilibrium equations. Therefore, the objective is to minimize the following residual:

$$ r(\alpha) = \mathbf{K}(\mu)\mathbf{D}_p\alpha + \mathbf{C}^T(\mu, u)\mathbf{D}_d\alpha - f(\mu) $$

Note that the residual has been defined using $\hat{u} = \hat{\lambda} = \alpha$, which is equivalent to using a monolithic dictionary, in which each column consists of a displacement snapshot stacked above the corresponding contact pressure snapshot. Usage of a monolithic dictionary not only reduces the number of unknowns but also makes the problem of minimizing the residual $||r(\alpha)||$ well-defined. Moreover, the residual can be projected on a low-dimensional subspace of the colsp($\mathbf{D}_p$), described by the matrix $\mathbf{B}$ containing its basis, while also adding sparsity constraints to the unknown $\alpha$, like in (5). Thus, the problem can be stated as:

$$ \min ||\alpha||_p \qquad (11a) $$

$$ \text{s.t.} \quad ||\mathbf{B}^T r(\alpha)|| < \varepsilon \qquad (11b) $$

$$ 1^T \alpha = 1 \qquad (11c) $$

$$ \alpha \ge 0 \qquad (11d) $$

The nnFOCUSS algorithm [30, 49], based on the FOCUSS algorithm [39], is a weighted-norm minimization technique that progressively induces sparsity to the solution(see Appendix A for details). The algorithm takes three inputs: the vector $\mathbf{x}$ that will be approximated, the dictionary $\mathbf{D}$ that serves as the approximation space, and a non-sparse non-negative initial guess $\alpha^0$ of approximation coefficients. It can be used to solve the convex hull approximation problem (11), by plugging the inputs defined as:

$$ \mathbf{D} = \begin{bmatrix} \mathbf{B}^T \mathbf{K}_{\text{mono}} \\ 1^T \end{bmatrix}, \quad \mathbf{x} = \begin{bmatrix} \mathbf{B}^T f(\mu) \\ 1 \end{bmatrix}, \quad \alpha^0 = \text{nnls}(\mathbf{D}, \mathbf{x}) $$

where $\mathbf{K}_{\text{mono}}$ is the operator for the monolithic residual $r(\alpha)$:

$$ \mathbf{K}_{\text{mono}} = \mathbf{K}(\mu)\mathbf{D}_p + \mathbf{C}^T(\mu, u)\mathbf{D}_d $$

nnls indicates the non-negative least squares. An open-source implementation of nnls provided in the SciPy library [45] is used.
A test is proposed which checks if Hypothesis 1 hold within a reasonable error using the so-called leave-one-out approach on the training set snapshots. The Hypothesis 1 can be verified by approximating the left-out snapshot with the convex hull of the rest of the snapshots using least-square criteria. The approximation error, referred to here as the convex hull least square (CHLS) error, can be used to assess the existence of a hypothetical low-dimensional convex set. Each dictionary element $d$ is approximated using the convex hull of $\tilde{D}_p$ and the error is recorded:

$$ \epsilon_{CHLS}[k] = \min_{\alpha} \frac{\|\tilde{D}_p\alpha - d\|_2}{\|d\|_2} \qquad (12) \\ \text{s.t. } \mathbf{1}^T \alpha = 1 $$

A high value of $\epsilon_{CHLS}$ might occur due to a lack of low-dimensional convexity or due to scattered training set data.

## B.1 Illustrative example: elastic rope-obstacle problem

The method is demonstrated on a problem where the convexity of feasible space is assured. The elastic rope-obstacle problem from references [13, 14, 19] is considered, where a 1D elastic rope, fixed at the two ends comes into contact with a fixed obstacle (Fig. 16). The elasticity $\nu(x, \gamma)$ of the rope is a function of parameter $\gamma$.

$$ \begin{aligned} \nu(x, \gamma)\nabla^2 u(x) &= f \quad \text{on } x \in [0, 1] \\ u(0) &= u(1) = 0 \\ u(x) &\geq -0.2(\sin(\pi x) - \sin(3\pi x)) - 0.5 \end{aligned} \qquad (13) $$

where

$$ \nu(x, \gamma) = \begin{cases} \gamma & \text{if } x < 0.5 \\ 30 & \text{otherwise} \end{cases} \quad \text{on } \gamma \in [10, 50] $$

Training set: $\gamma \in \mathcal{P}_{tr} = \{10, 15, 20 \dots 45, 50\}$

Validation set: $\gamma \in \mathcal{P}_{val} = \{12.5, 17.5, 22.5 \dots 42.5, 47.5\}$

The monolithic dictionary with 9 snapshots corresponding to the training set defined above is computed in the offline stage. The test results for Hypothesis 1 are shown in Fig. 17, indicating that the hypothesis holds to a reasonable level.

The low-rank matrix $\mathbf{B}$ on which the residual will be projected is computed using left singular vectors of the truncated SVD $\mathbf{B} \leftarrow \text{svd}(\mathbf{D}_p, \delta)$, with truncation tolerance $\delta = 10^{-7}$. The reconstruction errors and resultant sparsity of the solution in the training set and a validation set are shown in Fig. 18. The training set is reconstructed within near numerical precision. In the validation set, the algorithm chooses the training set points that are nearest to $\gamma$ to create the best possible reconstruction with the given dictionary.

As discussed earlier, the non-penetration condition and complementary slack are not explicitly imposed. However, as the feasible region is convex, one would expect that at least the non-penetration condition would be satisfied exactly. Nonetheless, as seen in the Table 3, reconstructions in validation set show a finite but small penetration and complementary slack. These violations appear because the condition (11c) is satisfied only to a certain precision. The violation to convex condition, evaluated as $|\mathbf{1}^T\alpha - 1|$ in Table 3, is of the similar order as penetration and complementary slack.
A test is proposed which checks if Hypothesis 1 hold within a reasonable error using the so-called leave-one-out approach on the training set snapshots. The Hypothesis 1 can be verified by approximating the left-out snapshot with the convex hull of the rest of the snapshots using least-square criteria. The approximation error, referred to here as the convex hull least square (CHLS) error, can be used to assess the existence of a hypothetical low-dimensional convex set. Each dictionary element $d$ is approximated using the convex hull of $\tilde{D}_p$ and the error is recorded:

$$ \epsilon_{CHLS}[k] = \min_{\alpha} \frac{\|\tilde{D}_p\alpha - d\|_2}{\|d\|_2} \qquad (12) \\ \text{s.t. } \mathbf{1}^T\alpha = 1 $$

A high value of $\epsilon_{CHLS}$ might occur due to a lack of low-dimensional convexity or due to scattered training set data.

## B.1 Illustrative example: elastic rope-obstacle problem

The method is demonstrated on a problem where the convexity of feasible space is assured. The elastic rope-obstacle problem from references [13, 14, 19] is considered, where a 1D elastic rope, fixed at the two ends comes into contact with a fixed obstacle (Fig. 16). The elasticity $\nu(x, \gamma)$ of the rope is a function of parameter $\gamma$.

$$ \begin{aligned} \nu(x, \gamma)\nabla^2 u(x) &= f \quad \text{on } x \in [0, 1] \\ u(0) &= u(1) = 0 \\ u(x) &\geq -0.2(\sin(\pi x) - \sin(3\pi x)) - 0.5 \end{aligned} \qquad (13) $$

where

$$ \nu(x, \gamma) = \begin{cases} \gamma & \text{if } x < 0.5 \\ 30 & \text{otherwise} \end{cases} \quad \text{on } \gamma \in [10, 50] $$

Training set: $\gamma \in \mathcal{P}_{tr} = \{10, 15, 20, \dots, 45, 50\}$

Validation set: $\gamma \in \mathcal{P}_{val} = \{12.5, 17.5, 22.5, \dots, 42.5, 47.5\}$

The monolithic dictionary with 9 snapshots corresponding to the training set defined above is computed in the offline stage. The test results for Hypothesis 1 are shown in Fig. 17, indicating that the hypothesis holds to a reasonable level.

The low-rank matrix $\mathbf{B}$ on which the residual will be projected is computed using left singular vectors of the truncated SVD $\mathbf{B} \leftarrow \text{svd}(\mathbf{D}_p, \delta)$, with truncation tolerance $\delta = 10^{-7}$. The reconstruction errors and resultant sparsity of the solution in the training set and a validation set are shown in Fig. 18. The training set is reconstructed within near numerical precision. In the validation set, the algorithm chooses the training set points that are nearest to $\gamma$ to create the best possible reconstruction with the given dictionary.

As discussed earlier, the non-penetration condition and complementary slack are not explicitly imposed. However, as the feasible region is convex, one would expect that at least the non-penetration condition would be satisfied exactly. Nonetheless, as seen in the Table 3, reconstructions in validation set show a finite but small penetration and complementary slack. These violations appear because the condition (11c) is satisfied only to a certain precision. The violation to convex condition, evaluated as $|\mathbf{1}^T\alpha - 1|$ in Table 3, is of the similar order as penetration and complementary slack.
<figure><img src="image_21.png" /><figcaption>Figure 16: Elastic rope-obstacle problem: Snapshots of the elastic rope fixed at the two ends, while it deforms under a uniform load and establishes contact with the obstacle. Snapshots for various values of $\gamma$ are shown.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 17: Testing the Hypothesis 1 for the given training set of elastic rope-obstacle problem (13)</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Figure 16: Elastic rope-obstacle problem: Snapshots of the elastic rope fixed at the two ends, while it deforms under a uniform load and establishes contact with the obstacle. Snapshots for various values of $\gamma$ are shown.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Figure 17: Testing the Hypothesis 1 for the given training set of elastic rope-obstacle problem (13)</figcaption></figure>
<table><caption>Table 3: Observed values of various compliances for the elastic rope-obstacle problem with parameter $\gamma$.</caption><thead><tr><th></th><th>γ ∈ P<sub>tr</sub></th><th>γ ∈ P<sub>val</sub></th></tr></thead><tbody><tr><td>|1<sup>T</sup>α − 1|</td><td>O(10<sup>−13</sup>)</td><td>O(10<sup>−3</sup>)</td></tr><tr><td>Penetration</td><td>O(10<sup>−16</sup>)</td><td>O(10<sup>−4</sup>)</td></tr><tr><td>Complementary Slack</td><td>O(10<sup>−15</sup>)</td><td>O(10<sup>−3</sup>)</td></tr></tbody></table>

<figure><img src="image_23.png" /><figcaption>Figure 18: Convex hull reconstruction errors for elastic rope-obstacle problem with parameter $\gamma$ using $\mathbf{B} \leftarrow svd(\mathbf{D}_p, \delta = 10^{-7})$. Crosses × indicate achieved sparsity (w.r.t. right y-axis). Primal and dual sparsity are naturally equal as the dictionary is monolithic.</figcaption></figure>

## C Spurious effects of POD truncation

In the greedy active-set method, the snapshots of the dual dictionary are selected based on the violations of the non-penetration constraints. The dictionary vector with the highest correlation with the current state of penetration is chosen. However, instead of zero penetration, a small value of penetration $\tau$ is allowed (see (6c) and (7)). This is done to avoid spurious selection of dictionary elements, which happens due to a truncated primal basis. Such a reconstruction example of Hertz problem for parameter value $d = 0.25$ is shown in Fig. 19, where a spurious peak in the contact pressure is evident. In this case, a dual dictionary of size 30 and primal basis truncated at $10^{-6}$ are used. The algorithm selected the appropriate snapshots from the dual dictionary, but also selected some unnecessary snapshots: the first ($\mathbf{D}_d[1]$) and twelfth ($\mathbf{D}_d[12]$) dictionary vectors. Importantly, the $\mathbf{D}_d[1]$ vector contributes to large error in reconstruction as evident in the figure. The false selection of this snapshot is linked to the truncation of primal basis, and therefore, can be avoided by setting the value of $\tau$ same as $\delta$.

## About the authors

Kiran Sagar Kollepara, José V. Aguado, Luisa Silva, Domenico Borzacciello

Institut de Recherche en Génie Civil et Mécanique (GeM) at Ecole Centrale de Nantes
1 rue de la Noë, BP 92101, 44321 Nantes cedex 3, France
<table><caption>Table 3: Observed values of various compliances for the elastic rope-obstacle problem with parameter $\gamma$.</caption><thead><tr><th></th><th>γ ∈ P<sub>tr</sub></th><th>γ ∈ P<sub>val</sub></th></tr></thead><tbody><tr><td>|1<sup>T</sup>α − 1|</td><td>O(10<sup>−13</sup>)</td><td>O(10<sup>−3</sup>)</td></tr><tr><td>Penetration</td><td>O(10<sup>−16</sup>)</td><td>O(10<sup>−4</sup>)</td></tr><tr><td>Complementary Slack</td><td>O(10<sup>−15</sup>)</td><td>O(10<sup>−3</sup>)</td></tr></tbody></table>

<figure><img src="image_24.png" /><figcaption>Figure 18: Convex hull reconstruction errors for elastic rope-obstacle problem with parameter $\gamma$ using $\mathbf{B} \leftarrow svd(\mathbf{D}_p, \delta = 10^{-7})$. Crosses $\times$ indicate achieved sparsity (w.r.t. right y-axis). Primal and dual sparsity are naturally equal as the dictionary is monolithic.</figcaption></figure>

## C Spurious effects of POD truncation

In the greedy active-set method, the snapshots of the dual dictionary are selected based on the violations of the non-penetration constraints. The dictionary vector with the highest correlation with the current state of penetration is chosen. However, instead of zero penetration, a small value of penetration $\tau$ is allowed (see (6c) and (7)). This is done to avoid spurious selection of dictionary elements, which happens due to a truncated primal basis. Such a reconstruction example of Hertz problem for parameter value $d = 0.25$ is shown in Fig. 19, where a spurious peak in the contact pressure is evident. In this case, a dual dictionary of size 30 and primal basis truncated at $10^{-6}$ are used. The algorithm selected the appropriate snapshots from the dual dictionary, but also selected some unnecessary snapshots: the first ($\mathbf{D}_d[1]$) and twelfth ($\mathbf{D}_d[12]$) dictionary vectors. Importantly, the $\mathbf{D}_d[1]$ vector contributes to large error in reconstruction as evident in the figure. The false selection of this snapshot is linked to the truncation of primal basis, and therefore, can be avoided by setting the value of $\tau$ same as $\delta$.

## About the authors

Kiran Sagar Kollepara, José V. Aguado, Luisa Silva, Domenico Borzacciello

Institut de Recherche en Génie Civil et Mécanique (GeM) at Ecole Centrale de Nantes
1 rue de la Noë, BP 92101, 44321 Nantes cedex 3, France
<figure><img src="image_24.png" /><figcaption>Figure 19: Reconstruction instance of Hertz problem with $d = 0.25$ using $\#\mathcal{P}_{tr} = 30$ and primal basis truncated at $\delta = 10^{-6}$, where few dictionary snapshots were spuriously selected by the greedy active-set algorithm if non-penetration condition is applied “hardly”. The falsely selected snapshots are indicated by their index in the dual dictionary. The snapshot which contributes to most of the error is $\mathbf{D}_d[1]$, corresponds to loading parameter $d = 0.01$</figcaption></figure>

e-mail: {kiran-sagar.kollepara, jose.aguado-lopez, luisa.rocha-da-silva, domenico.borzacchiello}@ec-nantes.fr

web: http://gem.ec-nantes.fr

Yves Le Guennec

IRT Jules Verne
1 Mail des 20 000 Lieues, 44340 Bouguenais, France

web: http://irt-jules-verne.fr
<figure><img src="image_25.png" /><figcaption>Figure 19: Reconstruction instance of Hertz problem with $d = 0.25$ using $\#\mathcal{P}_{tr} = 30$ and primal basis truncated at $\delta = 10^{-6}$, where few dictionary snapshots were spuriously selected by the greedy active-set algorithm if non-penetration condition is applied “hardly”. The falsely selected snapshots are indicated by their index in the dual dictionary. The snapshot which contributes to most of the error is $\mathbf{D}_d[1]$, corresponds to loading parameter $d = 0.01$</figcaption></figure>

e-mail: {kiran-sagar.kollepara, jose.aguado-lopez, luisa.rocha-da-silva, domenico.borzacchiello}@ec-nantes.fr

web: http://gem.ec-nantes.fr

Yves Le Guennec

IRT Jules Verne

1 Mail des 20 000 Lieues, 44340 Bouguenais, France

web: http://irt-jules-verne.fr