# On the limitations of low-rank approximations in contact mechanics problems

Kiran Sagar Kollepara, José M Navarro-jiménez, Yves Le Guennec, Luisa Silva, José V Aguado

► To cite this version:

Kiran Sagar Kollepara, José M Navarro-jiménez, Yves Le Guennec, Luisa Silva, José V Aguado. On the limitations of low-rank approximations in contact mechanics problems. International Journal for Numerical Methods in Engineering, 2023, 124 (1), pp.217-234. 10.1002/nme.7118. hal-03918970

HAL Id: hal-03918970

https://hal.science/hal-03918970

Submitted on 25 Oct 2023

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>

Distributed under a Creative Commons Attribution 4.0 International License
# On the limitations of low rank approximations in contact mechanics problems

Kiran Sagar Kollepara
Yves Le Guennec

José V. Aguado
José M. Navarro

Luisa Silva

## Abstract

Typical strategies for reducing the computational cost of contact mechanics models use low-rank approximations. The underlying hypothesis is the existence of a low-dimensional subspace for the displacement field and a non-negative low-dimensional subcone for the contact pressure. However, given the local nature of contact, it seems natural to wonder whether low-rank approximations are a good fit for contact mechanics or not. In this paper, we investigate some of their limitations and provide numerical evidence showing that contact pressure is linearly inseparable in many practical cases. To this end, we consider various mechanical problems involving non-adhesive frictionless contacts and analyse the performance of the low-rank models in terms of three different criteria, namely compactness, generalization and specificity.

**Keywords:** contact mechanics, variational inequalities, reduced order models, low-rank approximation

## 1 Introduction

Reduced order models (ROM) have been the key to the application of numerical models in near real-time simulations. ROMs reduce the computational complexity of discretized partial differential equations, which permits their usage in a multi-query context. ROMs have found successful applications in a wide variety of problems such as design optimization [1], uncertainty quantification [2], inverse problems [3], optimal control [4], real-time monitoring [5], and others. At the core of these applications, is the idea of solving a parametrized problem, where one or more input parameters of the numerical model influence the state of the system. ROMs are used to generate parametric models that are able to deal with multi-query situations efficiently.

The underlying principle behind many ROMs is the low-rank hypothesis that states that many real-life mathematical/numerical models defined in high-dimensional spaces can be fitted using low-dimensional spaces. Consequently, models can be expressed using only a few degrees of freedom that have been intelligently selected. Typically, this is done using dimensionality reduction
techniques that involve computing reduced bases (RB) for the unknowns of the problem. The reduced bases contain the underlying physical structure, thus spanning the low-dimensional space where the model lies.

ROMs have been applied extensively to unconstrained problems such as thermal and mechanics problems [6, 7] and also to mixed problems with equality constraints such as the incompressibility condition in Stokes and Navier-Stokes flow problems [8]. Application of ROMs to variational problems with inequality constraints (also referred to as variational inequalities) has been more recent [9–15]. Although the scope of this article is limited to contact mechanics, several other applications of variational inequalities are found in porous media flow problems [16], cavitation problems in lubrication systems [17], anti-plane frictional problems [18] and even in financial trading problems [10]. Inequality constraints appear in mechanical problems with obstacles or multi-body mechanical problems where there is a possibility of contact between bodies and obstacles, or with other bodies [19]. Moreover, the region of contact is unknown a-priori.

Inequality constraint problems are often posed in mixed form with the so-called Karush–Kuhn–Tucker (KKT) conditions that force the Lagrange Multipliers to be non-negative. Enforcing this condition in ROMs is non-trivial as traditional methods of computing RBs do not preserve the non-negative nature of the input information.

One of the first works on reducing the contact mechanics problem [9], proposed the idea of using the contact pressure snapshots directly to define a non-negative subspace. Both displacement and contact pressure snapshots were generated in a greedy fashion using error estimators. Fundamental aspects of the reduced problem, such as the existence and uniqueness of the solution and inf-sup stability were also explored.

Compression of snapshots to create a reduced basis was studied in [11], where Non-Negative Matrix Factorization (NNMF) was used to compute a basis with user-specified cardinality, but does not provide any means to specify truncation tolerance for the reduced basis. Also, error estimators for greedy sampling of parametric space are developed.

Algorithms to sort a precomputed set of snapshots, in order of importance, to create a compact basis are explored using projection methods in [10] and [14], using an Angle-Greedy and Cone-Projected Greedy (CPG) procedures, respectively. The former does not take into account non-negative restrictions, whereas the latter uses a cone projection involving non-negative coefficients. This makes CPG more efficient in capturing the contact pressure subcone.

The hyper-reduction approach in [13] defines a subdomain of the contact problem containing the most important points. The reduction is achieved with the usage of a POD basis for displacement and resolution of the weak form on the reduced integration domain. Reconstruction of displacement solution on the full domain is relatively straightforward using the POD basis. On the other hand, the reconstruction of contact pressure involves solving a non-negative least square problem using the snapshots.

A Craig-Bampton based resolution of the contact problem was discussed
in [15] where the reduction of the displacement field was achieved using the Krylov subspace method. The Lagrange Multiplier method was not reduced under the assumption that the number of contact dofs remains small.

All of these works are based on a Lagrange Multiplier approach to solve the inequality constrained reduced problem. An exception to this trend is [12], where the penalty approach is used. The Empirical Interpolation Method (EIM) is employed for efficient computation of barrier functions that restrict the solution to the feasible domain.

Most contributions build reduced models on problems with static contact pairs, meaning contact pairs do not change with the state of the system. However, [14] considers dynamic contact pairs, with node-to-segment formulation. To efficiently evaluate non-linearities due to dynamic pairing, they use EIM to define an affine decomposition of distance functions.

Most ROMs discussed above involve an implicit assumption that contact pressure space can be approximated by a low-dimensional subspace, i.e., it satisfies the low-rank hypothesis. However, given the local nature of contact, the applicability of the low-rank hypothesis may be questionable. As already pointed out in Reference [13], singular values for contact pressure decay slowly compared to those of displacement field, the first hint on lack of low-dimensional subspace. This article focuses on studying the validity of the low-rank hypothesis to the variational inequality problems, specifically contact mechanics problems and its contact pressure field. We show that in many practical cases, it seems difficult to justify the existence of low-rank representation. We provide a quantitative assessment of low-rank approximations' robustness through validation metrics [20] like compactness, generalization ability and specificity.

The paper is organized as follows: The variational inequalities for contact mechanics are first introduced in Sec. 2.1 and 2.2. Then a generic methodology of reduced modelling in variational inequalities is discussed in Sec. 2.3. In Sec. 3, a reduced model of a simple 1D obstacle problem is demonstrated, followed by the main argument of this article on the linearly inseparable nature of Lagrange multipliers as a limitation of the low-rank approach. To endorse this argument quantitatively, validation metrics are defined in Sec. 4. The arguments about linear inseparability are supported using validation metrics for specific numerical examples in Sec. 5. Finally, concluding arguments about linear inseparability and a few perspectives on circumventing the limitations due to lack of low-rank structure are given.

## 2 Model problem and basic low-rank formulation of contact problems

In this section, mechanical problems with constraints on the displacement field are discussed. Linear elastic and small deformation problems are considered, keeping the evaluation of internal energy simple. The inequality constraints may be a result of the presence of either an obstacle or a second body in the
<figure><img src="image_2.png" /><figcaption>Figure 1: Kinematic description of mechanical problems with possibility of contact. Left: Type 1: Body-obstacle problem and Right: Type 2: Two body problem</figcaption></figure>

domain, with the possibility of contact between different surfaces in the deformed configuration. For convenience, these problems will be referred to as Type-1 problems, involving a deformable body and a fixed obstacle, and Type-2 problems, involving two deformable bodies. A schematic of the two types of contact problems is shown in Fig. 1. Contact phenomena like friction and adhesion are neglected. The resolution of contact mechanics problems depends on computing the distance between the surfaces of bodies involved. As the evaluation of distance is not a trivial task, different strategies exist based on underlying simplifications. For the generic case of non-conforming meshes, this paper uses node-to-segment formulation along with the closest point projection. The details to these approaches can be found in [19, 21].

## 2.1 Model problem

The generic weak form of a parametrized mechanical problem involving contacts is described by the following inequality constrained minimization problem:

$$ 
\begin{array}{ll}
\mathbf{u} = \arg \min_{\mathbf{v} \in \mathcal{V}} \frac{1}{2} a(\mu; \mathbf{v}, \mathbf{v}) - f(\mu; \mathbf{v}) & (1) \\
\text{s.t.} \quad k(\mu, \mathbf{v}; \mathbf{v})(\mathbf{x}) \le g(\mu, \mathbf{v})(\mathbf{x}) & \text{on } \Gamma^2
\end{array}
$$

where

* $a(\mu; \mathbf{v}, \mathbf{v})$ is strain energy function and $f(\mu; \mathbf{v})$ is the work contribution by the external forces. The operators are defined $\forall \mathbf{v} \in \mathcal{V}$ and $\forall \mu \in \mathcal{P}$, where $\mathcal{V}$ be an appropriate function space for displacement and $\mathcal{P}$ is the
parametric space. Note that $\mu$ may not be necessarily a scalar and the parametric space $\mathcal{P}$ can be possibly multidimensional.

* $k(\mu, v; v)(x)$ and $g(\mu, v)(x)$ are the distance functions that indicate separation between bodies, defined conveniently on surface $\Gamma^2$. The first term indicates the contribution of displacement field $v$ to the distance and the second term indicates the distance in the undeformed configuration.

## NOTATION

Strong and weak forms are distinguished using notations adopted from [14]. Linear and bilinear operators that appear in weak formulations are expressed in the format $o(\cdot; \cdot, \cdot)$, where the non-linear dependencies and linear dependencies are separated by the semi-colon ';'. On the other hand, strong form functions are expressed in the format $o(\cdot; \cdot, \cdot)(\cdot)$, distinguished from weak forms using an additional argument which is the geometric position. The same applies to linear forms expressed as $o(\cdot; \cdot)$ and $o(\cdot; \cdot)(\cdot)$. Note that the notation $\mathbf{u}$ and $\mathbf{v}$ are used for both continuous and discrete versions of displacement fields.

Both Type-1 and Type-2 problems can be expressed in the above form. As a result of the node-to-segment formulation and closest point projection, the distance function for Type-2 can be expressed as:

$$ 
\begin{aligned}
k(\mu, v; v)(\mathbf{x}) &= [v^2(\mathbf{x}^2) - v^1(\bar{\mathbf{x}}^1)] \cdot \mathbf{n}^1 \\
g(\mu, v)(\mathbf{x}) &= [\mathbf{x}^2 - \bar{\mathbf{x}}^1] \cdot \mathbf{n}^1
\end{aligned} 
\quad (2) 
$$

where $\bar{\mathbf{x}}^1$ is the nearest to point $\mathbf{x}^2$, on the surface $\Gamma^2$, as shown in Fig. 1. The pairing $(\bar{\mathbf{x}}^1, \mathbf{x}^2)$ is evaluated in the deformed configuration and hence is dependent on the displacement field. This forces the distance function $k$ to have an implicit non-linear dependence on displacement, apart from the explicit linear dependence that is evident in the expression. For Type-1 problem, $v^1(\bar{\mathbf{x}}^1)$ in the above expressions vanishes as the obstacle is fixed.

## 2.2 Mixed form

The constrained minimization problem can be expressed in a mixed form by introducing Lagrange multipliers and KKT conditions. The resulting optimization problem can be expressed as a saddle point problem of the Lagrangian functional:

$$ (\mathbf{u}, \lambda) = \arg\min_{\mathbf{v} \in \mathcal{V}, \eta \in \mathcal{W}_+} \mathcal{L}(\mathbf{v}, \eta) \quad (3) $$

where

$$ \mathcal{L}(\mathbf{v}, \eta) = \frac{1}{2}a(\mu; \mathbf{v}, \mathbf{v}) - f(\mu; \mathbf{v}) + b(\mu, \mathbf{u}; \eta, \mathbf{u}) - d(\mu, \mathbf{u}; \eta) \quad (4) $$
and the distance function is expressed by following weak forms:

$$b(\mu, \mathbf{u}; \eta, \mathbf{u}) = \int_{\Gamma^2} \eta k(\mu, \mathbf{u}; \mathbf{u})(\mathbf{x}) \partial \Gamma$$

$$d(\mu, \mathbf{u}; \eta) = \int_{\Gamma^2} \eta g(\mu, \mathbf{u})(\mathbf{x}) \partial \Gamma$$

$W_+$ indicates a non-negative function cone for contact pressure defined on the surface $\Gamma^2$. The condition $\lambda \ge 0$ is implied naturally as $W_+$ is a function cone and admits a physical meaning because $\lambda$ is equivalent to contact pressure. A negative contact pressure implies traction, in other words, adhesion between contact surfaces which contradicts the simplifying assumptions.

The KKT conditions associated to the optimizer $(\mathbf{u}, \lambda) \in \mathcal{V} \times W_+$ of (4) can be expressed as follows:

$$ \begin{aligned} a(\mu; \mathbf{v}, \mathbf{u}) - f(\mu; \mathbf{v}) + b(\mu, \mathbf{u}; \lambda, \mathbf{v}) &= 0, & \mathbf{v} \in \mathcal{V} \\ b(\mu, \mathbf{u}; \eta, \mathbf{u}) - d(\mu, \mathbf{u}; \eta) &\le 0, & \eta \in W_+ \\ b(\mu, \mathbf{u}; \lambda, \mathbf{u}) - d(\mu, \mathbf{u}; \lambda) &= 0 \end{aligned} \quad (5) $$

The discrete equivalent of the Lagrangian functional in (4) can be derived by introducing finite element spaces $\mathcal{V}^h \subset \mathcal{V}$ and $W_+^h \subset W_+$:

$$ \mathcal{L}^h(\mathbf{v}^h, \eta^h) = \frac{1}{2} \mathbf{v}^{hT} \mathbf{K}(\mu) \mathbf{v}^h - \mathbf{v}^{hT} \mathbf{f}(\mu) + \eta^h T [\mathbf{C}(\mu, \mathbf{v}^h) \mathbf{v}^h - \mathbf{g}(\mu, \mathbf{v}^h)] \quad (6) $$

where $\mathbf{v}^h \in \mathbb{R}^{N_u^h}$ and $\eta^h \in \mathbb{R}^{N_\lambda^h}$. Here, $N_u^h$ and $N_\lambda^h$ indicate the number of respective finite element degrees of freedoms (dofs). Given the finite element basis $\{N_i\}$ for displacement, a vector field, and $\{M_j\}$ for contact pressure, a scalar field defined on potential contact surface, the discrete operators can be defined as:

$$ \begin{aligned} \mathbf{K}(\mu)_{ij} &= a(\mu; N_i, N_j) \\ \mathbf{C}(\mu, \mathbf{v}^h)_{ij} &= b(\mu, \mathbf{v}^h; M_i, N_j) \end{aligned} $$

and the vectors $\mathbf{f}$ and $\mathbf{g}$ are defined as:

$$ \begin{aligned} \mathbf{f}(\mu)_i &= f(\mu; N_i) \\ \mathbf{g}(\mu, \mathbf{v}^h)_i &= d(\mu, \mathbf{v}^h; M_i) \end{aligned} $$

The bilinear and linear terms $b$ and $d$ are non-linear w.r.t. the displacement field $\mathbf{v}$ and the parameter $\mu$, which is reflected in discrete form as $\mathbf{C}(\mu, \mathbf{v}^h)$ and $\mathbf{g}(\mu, \mathbf{v}^h)$ respectively. Using the above operators, the discrete form of the KKT conditions in (5) for the optimal solution $(\mathbf{u}^h, \mathbf{\lambda}^h) \in \mathbb{R}^{N_u^h} \times \mathbb{R}^{N_\lambda^h}$ can be written
as:

$$
\begin{align}
& \mathbf{K}(\mu)\mathbf{u}^h - \mathbf{f}(\mu) + \mathbf{C}^T(\mu, \mathbf{u}^h)\boldsymbol{\lambda}^h = \mathbf{0} \nonumber \\
& \qquad \boldsymbol{\lambda}^h \ge \mathbf{0} \nonumber \\
& \mathbf{C}(\mu, \mathbf{u}^h)\mathbf{u}^h - \mathbf{g}(\mu, \mathbf{u}^h) \le \mathbf{0} \tag{7} \\
& \boldsymbol{\lambda}^h \mathbf{T} [\mathbf{C}(\mu, \mathbf{u}^h)\mathbf{u}^h - \mathbf{g}(\mu, \mathbf{u}^h)] = \mathbf{0} \nonumber
\end{align}
$$

Various strategies to linearize the above problems exist, including fixed point and Newton-Raphson methods. The details of above formulations and lineariza- tion strategies are not detailed here, but can be found in references [19, 22, 23]. Finite Element computations for the contact problems in this paper were per- formed using a Newton-Raphson approach.

2.3 Low-rank approach to parametrized contact mechan- ics

Resolution of parametric contact mechanics problems is not only hindered by the curse of dimensionality, but also by the non-linear inequality constraints. Low-rank approaches can potentially address this issue by reducing the number of constraints [11, 14]. The idea behind ROMs is to split the cost of computation into two stages. The offline stage, where most of the computational complexity is resolved, consists of extracting the underlying structure of the system. This is done by generating a RB from a set of solutions from the high-fidelity solver. These solutions are usually referred to as snapshots in ROM literature. This step is typically performed only once to create the reduced model. The online step is performed every time a new query is submitted, in which the on-demand solution is computed using the reduced model, which is cheaper to evaluate than the high-fidelity model.

Using notation from the previous section, the high-fidelity (finite element)
space is denoted by $\mathbf{V}^h$. Let the reduced subspace be denoted by $\mathbf{V}^r = \operatorname{span}(\mathbf{\Phi})$,
where $\mathbf{\Phi}$ denotes the corresponding RB. In the online stage, solutions in the
subspace $\mathbf{V}^r$ are sought to approximate the high-fidelity solution. In discrete
sense, $\mathbf{\Phi}$ is a matrix with each column corresponding to a basis vector. Then,
the approximation $\mathbf{u}^r$ in the subspace can be expressed as:

$$
\mathbf{u}^r \approx \mathbf{\Phi} \hat{\mathbf{u}} \qquad \hat{\mathbf{u}} \in \mathbb{R}^{N_u^r} \tag{8}
$$

All reduced degrees of freedom and reduced operators are indicated by a hat $\hat$.

For constructing the RB Φ in the offline stage, a finite subset of the paramet-
ric space P<sub>tr</sub> ∈ P is explored and an approximate subspace V<sup>r</sup> is extracted from
the set of high-fidelity snapshots from the training set. The number of degrees
of freedom associated to the reconstruction problem in the online stage is same
as the cardinality of basis Φ, and it influences the efficiency of the online phase.
A low-rank basis Φ will, therefore, lead to a lower number of unknowns resulting
in an efficient RB. As this assumption holds true for many physical models, RB
methods are widely applicable. Otherwise, it may be difficult to generate a RB that approximates well the system behaviour in the entire parametric space.

Reduction of variational inequality problems poses additional complications that other mixed problems do not face. In mixed problems like incompressible Stokes flow or incompressible elasticity, the conservation of mass equation is essentially seen as an equality constraint and the pressure is seen as a Lagrange multiplier or penalization. This constraint is always active and, therefore, must always be included in the linearized discrete system. Also, negative pressures are not explicitly prohibited, in neither the full nor the reduced solvers [8], as the momentum equation only includes the gradient of pressure and is unaffected even if the whole pressure field is globally lifted by a constant value. On the other hand, for inequality constrained problems, only a subset of inequality constraints satisfy the equality and the inactive ones cannot be included in the linearized system matrices. Resolution of active and inactive constraints becomes part of the linearization scheme, along with the elimination of negative Lagrange multipliers.

As seen in 2.2, variational inequalities involve primal and dual unknowns, $u^h$ and $\lambda^h$. Therefore, in the reduced sense, one might introduce the terms primal and dual bases corresponding to the two unknown variables. Reduction of variational inequalities typically involves computing the reduced primal and dual bases, and then finding a solution in the spaces spanned by these bases. Such approaches are discussed in [9, 11, 14].

### 2.3.1 Offline stage

As two finite element subspaces were defined for the mixed problem in Section 2.2, similarly two reduced spaces must be defined for the reduced problem. For the primal variable, an orthogonal basis is generated using Proper Orthogonal Decomposition (POD). In POD setting, high-fidelity snapshots are generated over the training set $\mathcal{P}_{tr}$ are collected. Let $\mathbf{S}_{tr} = \{\mathbf{u}_s^h\}_{s=1}^{NS}$ be a matrix containing the training set snapshots. Left singular vectors of the snapshot matrix are taken as the RB, as shown in (9). This can be done by arranging snapshots as columns of a matrix and using the Singular Value Decomposition (SVD) or the Principal Component Analysis (PCA) methods. The most important singular vectors are selected using the singular values and the truncation tolerance $\delta$.

$$ \Phi \leftarrow svd(\mathbf{S}_{tr}, \delta) \qquad (9) $$

**Need for a non-negative basis:** Computing the dual basis is more complicated than the primal basis. The dual field i.e. the Lagrange multiplier field must satisfy the non-negativity constraints. In case of the contact problems without adhesive and cohesive surfaces, this constraint admits a physical meaning as negative contact pressure cannot be admitted, as discussed in Section 2.2. An orthogonal basis generated using projection based methods cannot satisfy such constraints, as it is surely bound to contain negative entries. To ensure such constraints, one way is to define a subcone for the dual field instead of
a subspace. A function cone, unlike a function space, is spanned by a set of non-negative basis functions and non-negative coefficients. In fact, $W_+$ and $W_+^h$ in Section 2.1 are also cones. By extension, a subcone, unlike a subspace, must be equipped by a non-negative RB and must be spanned by non-negative coefficients.

$$W_+^r = \operatorname{span}^+(\Theta) = \sum_{i=1}^{N} \hat{\lambda}_i \theta_i, \quad \hat{\lambda}_i \ge 0 \qquad (10)$$

where $\Theta$ and $\{\theta_i\}_{i=1}^N$ are the dual RB and dual RB functions. $\hat{\lambda}_i$ is the dual reduced dofs.

## Methods without a Lagrange Multiplier

The non-negativity condition appears because of the Lagrange Multiplier formulation, thereby forcing the necessity of a non-negative dual basis. In a penalty formulation, such as the barrier method in [12], the contact pressure variable does not appear explicitly and therefore, there is no need for computing a dual basis. However, this comes at the cost of additional approximation, as the penalty approach permits a small but non-zero penetration. Any attempt at minimizing the penetrations by increasing the penalty parameter leads to ill conditioning of the numerical problem.

The Augmented Lagrange Multiplier Method (ALMM) circumvents these limitations. In this approach, the Lagrange Multiplier is initialized at zero value, and the contact pressure contributed by penalty is added to the multiplier progressively. Despite the intervention of the penalty term, the non-negativity condition appears indirectly when the Lagrange Multipliers are updated. In a low-rank setting, the update step will potentially need a non-negative basis to avoid reconstruction of the full Lagrange Multiplier. Thus, the necessity of a non-negative dual basis prevails.

In the Nitsche method[24], the contact pressure is derived as surface traction from the constitutive model. Thus, the Lagrange Multiplier does not appear explicitly. However, the operators that enforce the non-penetration constraint might show dependence on the constitutive model, as mentioned in [19, Chapter 6], which may lead to additional difficulties in creating a reduced model in case of complex constitutive models.

Once the reduced subspace $V^r$ and subcone $W_+^r$ are available, the reduced problem can be generated from (5), replacing the continuous spaces $V$ and $W_+$ with the reduced spaces $V^r$ and $W_+^r$. The reduced KKT conditions on the solution $(\hat{u}, \hat{\lambda}) \in \mathbb{R}^{N_u^r} \times \mathbb{R}^{N_u^r}$ can be expressed as:
$$
\begin{align}
\hat{\mathbf{K}}(\mu)\hat{\mathbf{v}} - \hat{f}(\mu) + \hat{\mathbf{C}}^T(\mu, \hat{\mathbf{v}})\hat{\lambda} &= \mathbf{0} \nonumber \\
\hat{\lambda} &\ge \mathbf{0} \nonumber \\
\hat{\mathbf{C}}(\mu, \hat{\mathbf{u}})\hat{\mathbf{u}} - \hat{g}(\mu, \hat{\mathbf{u}}) &\le \mathbf{0} \nonumber \\
\hat{\lambda}^T(\hat{\mathbf{C}}(\mu, \hat{\mathbf{u}})\hat{\mathbf{u}} - \hat{g}(\mu, \hat{\mathbf{u}})) &= \mathbf{0} \tag{11}
\end{align}
$$

where the discrete reduced operators are built by introducing the reduced basis functions $\phi$ and $\theta$

$$
\begin{align*}
\hat{\mathbf{K}}(\mu)_{ij} &= a(\mu; \phi_i, \phi_j) \\
\hat{\mathbf{C}}(\mu, \hat{\mathbf{u}})_{ij} &= b(\mu, \mathbf{u}^r; \theta_i, \phi_j) \\
\hat{f}(\mu)_i &= f(\mu; \phi_i) \\
\hat{g}(\mu, \hat{\mathbf{u}})_i &= d(\mu, \mathbf{u}^r; \theta_i)
\end{align*}
$$

with $\mathbf{u}^r = \mathbf{\Phi}\hat{\mathbf{u}}$

**Generation of dual basis:** Non-negative nature of the dual RB $\Theta$ prevents us from using orthogonal decompositions of Lagrange Multiplier snapshot matrix $\Lambda$. To this end, non-negativity preserving decompositions have been explored by various authors. Non-negative matrix factorization (NNMF) method [25] can be used to decompose a non-negative matrix, such as $\Lambda$ into two low-rank non-negative matrices $\mathbf{W}$, $\mathbf{H}$, such that $\Lambda \approx \mathbf{W}\mathbf{H}$. The left-hand matrix $\mathbf{W}$ is used as the dual RB in [11].

On the other hand, one could directly use the snapshot vectors as the basis vectors. This was the approach of [9], where the possibility of non-unique dual solution is also discussed, as snapshot vectors are not guaranteed to be linearly independent. Also, the number of dual dofs increase with the number of snapshots, preventing optimal reduction of the system. A greedy snapshot selection method based on the criteria of maximizing the volume of the reduced cone, namely the Cone-Projected Greedy (CPG) algorithm, was proposed in [14]. This algorithm creates a dual RB by greedily selecting snapshots from the snapshot matrix, attempting to create a more compact basis in comparison to using the full snapshot matrix. The algorithm is based on projection $\Pi_{\Theta}$ of a vector $\lambda$ on a vector cone $\Theta$,

$$ \Pi_{\Theta}(\lambda) \approx \Theta\alpha, \text{ where } \alpha \text{ is } \arg\min_{\gamma} \|\lambda - \Theta\gamma\| \quad \forall \gamma \ge \mathbf{0} \quad (12) $$

The greedy algorithm evaluates the error between each snapshot and its projection on the cone spanned by previously selected snapshots, and then adds that snapshot with maximum cone projection error. The process continues until the cone projection error is within a set tolerance. [10] had also proposed a similar algorithm called Angle-Greedy algorithm, but the projection error is calculated based on computation of the angle between the candidate snapshot vector and the space, and not the cone, spanned by the previously selected snapshots i.e.
it does not place non-negativity constrain on the coefficients $\gamma$ in (12). This is not the best way of selecting snapshots for variational inequality problems since solutions are sought in the reduced cone, and not the entire reduced space.

### 2.3.2 Online Stage

The non-linear reduced problem in (11) is solved using the fixed-point iterations. The linearized problem, for a fixed point method, at an iteration level $p$ can be expressed as follows:

$$ \begin{bmatrix} \hat{\mathbf{K}} & \hat{\mathbf{C}}_A^T(\mu, \hat{\mathbf{u}}^p) \\ \hat{\mathbf{C}}_A(\mu, \hat{\mathbf{u}}^p) & \mathbf{0} \end{bmatrix} \begin{bmatrix} \hat{\mathbf{u}}^{p+1} \\ \hat{\mathbf{\lambda}}_A^{p+1} \end{bmatrix} = \begin{bmatrix} \hat{\mathbf{f}} \\ \hat{\mathbf{g}}_A(\mu, \mathbf{u}^p) \end{bmatrix} \quad (13) $$

The subscript $A$ indicates the active set of reduced contact constraints, i.e. the constraints that satisfy the equality and force the solution to lie on the boundary of the feasible region. With all reduced operators defined, the fixed point algorithm is stated in Alg. 1.

**Algorithm 1 Online phase**

1. Input: Queryed value of parameter $\mu$

2. Given: Primal basis $\mathbf{\Phi}$ and dual basis $\mathbf{\Theta}$
Reduced operators $\hat{\mathbf{K}}, \hat{\mathbf{f}}$ (possible to build offline)

3. Initiate boolean array $\mathbf{mask_active}$ with one random element set to True.

4. **while** $\hat{\mathbf{u}}, \hat{\mathbf{\lambda}}$ not converge **do**

5. Build constraint operators $\mathbf{C}(\mu, \hat{\mathbf{u}}^p), \mathbf{g}(\mu, \hat{\mathbf{u}}^p)$ using FEM

6. Project constraint operators on RBs

$$ \hat{\mathbf{C}}(\mu, \hat{\mathbf{u}}^p) = \mathbf{\Theta}^T \mathbf{C}(\mu, \hat{\mathbf{u}}^p) \mathbf{\Phi} $$

$$ \hat{\mathbf{g}}(\mu, \hat{\mathbf{u}}^p) = \mathbf{\Theta}^T \mathbf{g}(\mu, \hat{\mathbf{u}}^p) $$

7. Filter-out rows that are not in active set:

$$ \hat{\mathbf{C}}(\mu, \hat{\mathbf{u}}^p)_A = \hat{\mathbf{C}}(\mu, \hat{\mathbf{u}}^p)[\mathbf{mask_active},:] $$

$$ \hat{\mathbf{g}}(\mu, \hat{\mathbf{u}}^p)_A = \hat{\mathbf{g}}(\mu, \hat{\mathbf{u}}^p)[\mathbf{mask_active}] $$

8. Solve system (13)

9. Set $\hat{\mathbf{\lambda}}^{p+1}[\mathbf{mask_active}] \leftarrow \hat{\mathbf{\lambda}}_A^{p+1}$ and $\hat{\mathbf{\lambda}}^{p+1}[\mathbf{NOT mask_active}] \leftarrow \mathbf{0}$

10. Update active constraints set

$$ \mathbf{mask_active}[i] = \begin{cases} \text{False} & \text{if } \hat{\mathbf{\lambda}}_i^{p+1} < 0 \\ \text{True} & \text{if } (\hat{\mathbf{C}}\hat{\mathbf{u}}^{p+1} - \hat{\mathbf{g}} \geq 0)_i \end{cases} $$

11. **end while**

12. Reconstruct $\mathbf{u} = \mathbf{\Phi}\hat{\mathbf{u}}$ and $\mathbf{\lambda} = \mathbf{\Theta}\hat{\mathbf{\lambda}}$

13. Output: $\mathbf{u}, \mathbf{\lambda}$

**Construction of operators:** For an efficient ROM, the construction of non-linear operators must also be inexpensive. This is not the case in Alg. 1, where
the full order operators $C(\mu, \hat{u})$ and $g(\mu, \hat{u})$ are constructed in each iteration (Step 5). Affine decompositions of the nonlinear constraint operators using the EIM are discussed in [14]. EIM decomposition of the constraint operators further splits computational complexity into offline and online stages, by selecting the so-called “magic points” that are a small subset of all potential contact nodes where the nonlinear terms are evaluated in the online phase. The hyper-reduction method also permits efficient computation of operators by selecting fewer quadrature points [13]. The distance functions are then computed only on the selected quadrature points. This is crucial to solve problems where the contact pairs are strongly dependent on the displacement field. However, the efficient construction of constraint operators is not covered in this paper, and focus is maintained on the non-negative dual space and its low-rank approximation.

## 3 Illustrative case

In this section, the Type-1 problem of a 1D elastic rope-obstacle problem from [9, 11, 12] is considered. The geometry of the problem enjoys further simplification of static contact pairs. In other words, each point on the elastic rope can come in contact with a unique point on the obstacle, thereby stripping the distance functions $k$ and $g$ of their non-linear dependence on the displacement field $u$.

The parametrized model of the elastic rope with obstacle function can be expressed as:

$$ 
\begin{array}{ll}
\nu \nabla^2 u(x) = f & \text{on } x \in [0, 1] \\
u(0) = u(1) = 0 & (14) \\
u(x) \ge g(x, \gamma) & \text{on } \gamma \in [-0.5, 0.5]
\end{array}
$$

where the parameterized obstacle function is defined as:

$$ g(x, \gamma) = -0.2(\sin(\pi x) - \sin(3\pi x)) - 0.5 + \gamma x \quad (15) $$

and the quantities $\nu = 30$ and $f = 250$ are independent of the parameter $\gamma$.

Snapshots are generated in the training set, $P_{tr} \in P$ consisting of 10 equidistant points in the parametric space $P$ defined in (15). The deformation and contact pressure snapshots are shown in Fig. 2. The primal RB is built using the POD approach and dual RB is created using the cone-projected greedy algorithm of [14].

The reconstruction errors are shown in Fig. 3 for the two cases of full and truncated dual RB. As expected, the points in the training set show very low reconstruction error, and the points outside the training set show a moderate error. The only exception at $\gamma = -0.35$ where reconstruction with full RB is highly accurate because the contact area does not change in the regime $\gamma \in [-0.4, -0.3]$ and the training set has two snapshots in this region. When dual RB is truncated, some points in the training set achieve the same error level as that of points outside the training set. These points correspond to the same
<figure><img src="image_3.png" /><figcaption>Figure 2: Sample snapshots for $\gamma \in \mathcal{P}_{\text{tr}} \subset [-0.5, 0.5]$. Left: Deformation $u(x, \gamma)$ and obstacle $g(x, \gamma)$ snapshots. Right: Contact pressure snapshots</figcaption></figure>

<table><caption>Table 1: Truncation tolerances and ranks for each basis for the rope-ostacle problem. $\epsilon$ indicates numerical precision</caption><thead><tr><th></th><th></th><th>Primal basis<br>(POD)</th><th>Dual basis<br>(CP-greedy)</th></tr></thead><tbody><tr><td rowspan="2">Full</td><td>Tolerance</td><td>&epsilon;</td><td>&epsilon;</td></tr><tr><td>Rank</td><td>11</td><td>11</td></tr><tr><td rowspan="2">Truncated</td><td>Tolerance</td><td>&epsilon;</td><td>2 &times; 10<sup>-1</sup></td></tr><tr><td>Rank</td><td>11</td><td>7</td></tr></tbody></table>

dual snapshots that were eliminated by the truncation procedure of the CPG algorithm. Another relevant observation is that the number of active dofs of the reduced dual solution $\hat{\lambda}$ is small for all reconstruction cases. For points in the training set, $\hat{\lambda}$ has exactly 1 active dof, whereas for the points outside the training set, it has a maximum of 2 active dofs. This is due to the fact that the dual basis is composed primarily of snapshots and does not undergo any compression like in the case of the primal basis generated by POD.

In general, any efficiently generated reduced basis has its vectors arranged in order of importance. In case the basis consists of left singular vectors, generated by SVD, the decreasing importance of basis vectors is indicated by the singular values. This is a very useful indication of the rank of the subspace. Even though the dual basis is not computed using SVD, it can still be useful to compute singular values and assess their decay. This evolution of singular values, seen in Fig. 4, shows that the decay is slow compared to the primal variable.

Also, the discussions in previous sections highlight the fact that dual basis $\Theta$ cannot be truncated without a significant loss of accuracy. This is also evident from Table 1, even a truncation tolerance as high as 0.2 leads to truncation of just four vectors in dual RB.

The above observations are explainable by the hypothesis that the dual snap-
<figure><img src="image_4.png" /><figcaption>Figure 3: Primal reconstruction ($H_1$-norm) errors for points in and outside the training set for the rope-obstacle problem. Left: Full Rank. Right: Full rank primal and truncated dual RB</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Decay of singular values (cumulative and normalized) for the rope-obstacle problem. A similar figure is also given in Reference [13, Figure 4]</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>Figure 5: Sparsity pattern of dual snapshots for the rope-obstacle problem. Each row corresponds to a snapshot</figcaption></figure>

shots are highly “inseparable” and do not lie in a low-dimensional subcone. Each contact pressure snapshot in Fig. 2 is linearly independent of the other due to difference in the contact area. Hence, if any snapshot is picked randomly from this figure, it cannot be reasonably approximated by any linear combination of the rest of the snapshots. In other words, each snapshot vanishes outside the contact zone, and this zone is almost “unique” for each snapshot. For a particular snapshot to be reasonably approximated, other snapshots with similar contact zones are needed. Therefore, the elimination of snapshots from the dual basis leads to a drastic increase in the reconstruction error for the corresponding parametric values in the training set. The dual variable, therefore, not only needs special treatment due to its positivity constraints but is also highly inseparable because of sensitivity to contact position and area. The inseparability is also evident in the non-zero pattern (“spy” plot) of dual snapshots visualized in Fig. 5, where most rows have an almost unique sparsity pattern. It is unlikely that a subspace whose members show varying sparsity patterns will admit a low-dimensional behaviour.

More numerical examples are explored to study the lack of low-rank behaviour in the Sec. 5.

# 4 Metrics for low-rank methods

In this section, we introduce the validation metrics [20], that provide a quantitative evaluation of the robustness of low-rank models. Validation metrics defined in [20] identify three different measures for this purpose:

*   *Compactness:* As the name suggests, this metric is a measure of dimensionality of the reduced space. Reference [20] defines compactness as the squared sum of the first *m* singular values of the snapshot matrix. This metric is similar to the decay of singular values studied in the previous section in Fig. 4. However, an alternative definition given by (16) is used here. The new definition allows extending the concept of compactness to CPG reduced bases, where equivalents of singular values are undefined.
Also, in case of an orthogonal basis, this definition is equivalent to the normalized singular values used in Fig. 4.

$$C(m) := \frac{\|\mathbf{S}_{\text{tr}} - \Pi_{\Psi_{[m]}}(\mathbf{S}_{\text{tr}})\|_{F}}{\|\mathbf{S}_{\text{tr}}\|_{F}} \qquad (16)$$

where $\mathbf{S}_{\text{tr}}$ is a matrix containing (either primal or dual) snapshots in $\mathcal{P}_{\text{tr}}$ arranged column-wise. $\Pi_{\Psi_{[m]}}(\mathbf{S}_{\text{tr}})$ is a projection operators that projects each column in $\mathbf{S}_{\text{tr}}$ on the first $m$ vectors of the basis $\Psi$ (indicated in python notation by $\Psi[:m]$). The snapshot matrix $\mathbf{S}_{\text{tr}}$, projection operator $\Pi$ and the basis $\Psi$ are generic symbols and are used in the following combinations, indicated by the following labels:

<table><thead><tr><th>Label in figures</th><th>Snapshots in <b>S</b><sub>tr</sub></th><th>Projection operator <i>&Pi;</i></th><th>Basis <i>&Psi;</i></th></tr></thead><tbody><tr><td>primal orth</td><td>Primal</td><td>Orthogonal</td><td>POD primal basis <b>&Phi;</b></td></tr><tr><td>dual orth</td><td>Dual</td><td>Orthogonal</td><td>POD dual basis <b>&Upsilon;</b></td></tr><tr><td>dual cone</td><td>Dual</td><td>Cone</td><td>CPG dual basis <b>&Theta;</b></td></tr></tbody></table>

The Frobenius norm is chosen in (16), because it is a natural choice for a POD/SVD basis that contains the best rank-wise approximations in the sense of Frobenius norm, but other norms can be also be used. This is actually equivalent to the original definition of compactness in [20], because in the case $\Psi$ consists of the left singular vectors, compactness corresponds to cumulative energy of singular values.

For the dual field, the metric is given using two kinds of projections: orthogonal and CPG. Though the dual CPG metric is more relevant to the inequality constrained problems, the dual orthogonal metric is also included as it is more intuitive due to correspondence with the singular values.

* **Generalization Ability:** This metric measures the ability of the reduced basis to approximate parametric instances that are not in the training set. It is computed by using a leave-one-out approach on the training set to generate a $m$-rank basis. Then, the reconstruction error of the eliminated snapshot is used to compute generalization ability as given by (17),

$$G(m) = \frac{1}{N} \sum_{i=1}^{N} \frac{\| \mathbf{s}_i - \Pi_{\Psi(\mathbf{S}_{\text{tr}} \setminus \mathbf{s}_i)}(\mathbf{s}_i) \|}{\| \mathbf{s}_i \|} \qquad (17)$$

where $\Psi(\mathbf{S}_{\text{tr}} \setminus \mathbf{s}_i)$ is the basis created using the snapshot subset $\mathbf{S}_{\text{tr}}$ with $\mathbf{s}_i$ removed (i.e. $\mathbf{S}_{\text{tr}} \setminus \mathbf{s}_i$). The symbols $\mathbf{S}_{\text{tr}}$, $\Pi$ and $\Psi$ follow the definition in (16).

* **Specificity:** The subspace spanned by the reduced basis is expected to contain elements similar to that of the training set. This metric measures
the extent of dissimilarity between randomly picked elements from the subspace and the snapshots in the training set. To compute specificity, a set of random shapes in the associated subspace are generated using random coefficients $R = \{\Psi\alpha_i\}_{i=1}^N$, where each vector $\alpha_i$ is drawn randomly within a predefined range, and compared to the closest snapshot in the training set. The mean value of this error is defined as Specificity.

$$S(m) = \frac{1}{N_s} \sum_{i=1}^{N_s} \min_{\mathbf{s}_j \in \mathbf{S}_{\text{tr}}} \frac{\|\mathbf{s}_j - \Psi[:m]\alpha_i\|}{\|\mathbf{s}_j\|} \quad (18)$$

The vectors $\{\alpha_i | \alpha_i \in \mathbb{R}^m\}_{i=1}^{N_s}$ are drawn randomly in the range defined by $\mu_\alpha \pm \sigma_\alpha \subset \mathbb{R}^m$. The vectors $\mu_\alpha$ and $\sigma_\alpha$ contain the element-wise mean and standard deviation of vectors $\{\gamma_k\}_{k=1}^m$, where $\gamma_k$ is the reduced coordinate of the training set snapshot $k$. (e.g. if $\Psi$ is sourced from SVD of a centered $\mathbf{S}_{\text{tr}}$, the vectors $\gamma_k$ are the columns of matrix $\Gamma := \Psi^T\mathbf{S}_{\text{tr}}$, and also the standard deviations $\sigma_\alpha$ can be computed as $\sigma_s/\sqrt{N}$ where $\sigma_s$ are the singular values). In the specific case where the basis $\Psi$ is the $\Theta$, the CPG basis, only non-negative entries are allowed in $\{\alpha_i\}$. The symbols $\mathbf{S}_{\text{tr}}$, $\Psi$ follows the definition in (16).

# 5 Numerical Examples

## 5.1 Hertz problem

This section deals with the reduced model of the Hertz contact problem of two half-cylinders loaded against each other, as shown in Fig. 6a. A loading parametrization problem is solved, where the imposed displacement $d \in (0, 0.3)$ is the parametric space, with $R_1 = R_2 = 1.0$. The contact pressure snapshots for various imposed displacement values are shown in Fig. 6b.

**High-fidelity model:** A finite element model of the geometry in Fig. 6a is created using a quad mesh. Each half-cylinder has 513 nodes and 466 elements of which 78 elements lie of the potential contact surface i.e. the semi-circular edges. Linear shape functions are used to discretize the displacement field, while contact pressure is discretized using piece-wise constant shape functions centred at the surface nodes of the quad mesh (collocation method, as described in [14]). For surface integrals, a single-point gauss quadrature centred on the node is used. This formulation gives mostly smooth contact pressure profiles, except near the peak pressure.

**Reduced Model:** Snapshots are generated in the training set, $\mathcal{P}_{\text{tr}} \in \mathcal{P}$ consisting of 12 equidistant points in the parametric space. The full rank and truncated dual RB (using CPG) is considered, as given in Table 2, whereas the primal RB is not truncated. The full RB has 12 basis vectors for both primal

<footnote n="1">¹The mean of snapshot matrix is removed from each column</footnote>
<figure><img src="image_7.png" /><figcaption>Figure 6: Left: Hertz problem: Two half cylinders loaded against each other. Displacement *d*, imposed on the top cylinder, is treated as parameter in the reduced model. Right: Contact pressure snapshots</figcaption></figure>

<table><caption>Table 2: Truncation tolerances and ranks for each basis for Hertz problem. $\epsilon$ indicates numerical precision</caption><thead><tr><th></th><th></th><th>Primal basis<br>(POD)</th><th>Dual basis<br>(CP-greedy)</th></tr></thead><tbody><tr><td rowspan="2">Full</td><td>Tolerance</td><td>&epsilon;</td><td>&epsilon;</td></tr><tr><td>Rank</td><td>12</td><td>12</td></tr><tr><td rowspan="2">Truncated</td><td>Tolerance</td><td>&epsilon;</td><td>5 &times; 10<sup>-2</sup></td></tr><tr><td>Rank</td><td>12</td><td>8</td></tr></tbody></table>

and dual fields, whereas the truncated dual RB has 8 vectors. Note that dual RB truncates only 4 out of 12 basis vectors for a high truncation tolerance of 0.05.

The reconstruction errors are shown in Fig. 7 for full and truncated dual RBs. As expected, the points in the training set show very low reconstruction error, and the points outside the training set show a relatively high error. When the dual RB is truncated, points corresponding to truncated dual snapshots have a high error level, in the same order of error for points outside the training set. Like the illustrative case in Sec. 3, the non-zero pattern in snapshots in Fig. 8 also show a unique sparsity pattern for each snapshot. Also, the comptactness shown in Fig. 9 displays a slow decay of the dual orth and dual cone, indicating a high-rank behaviour of the dual variable. The slightly slower dual cone curve than the dual orth was expected because of the non-negativity constraint that appears in approximation using the cone of the dual basis. These observations again reinforce the proposition that the dual basis does not admit a low-rank behaviour.

Generalization ability and Specificity metrics introduced in Sec. 4 are plotted in Fig. 10. The generalization ability for primal field is much better than the
<figure><img src="image_8.png" /><figcaption>Figure 7: Primal reconstruction ($H_1$-norm) errors for the Hertz problem. Left: Full Rank and Right: Full rank primal and truncated dual RB</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: Sparsity pattern of dual snapshots for Hertz problem</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Compactness of reduced bases for Hertz problem</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 10: Metrics for Hertz problem. Left: Generalization ability and Right: Specificity. The bars indicate ±1σ interval of corresponding metrics</figcaption></figure>

dual field, another reflection of separability issues with the dual field. On the other hand, primal and dual specificity are about the same order, still the primal specificity is consistently lower.

## 5.2 Ironing problem

A Type-2 problem where the contact area changes quite significantly, namely the ironing problem [21], is considered. The ironing problem consists of an iron block pressed against a flat slab and moved along the length of the slab (Fig. 11). The problem is simplified with two more assumptions: the first is that iron moves slowly enough that the problem can be considered to be *quasi-static* and the second is that surfaces are frictionless. The horizontal position $d_x$ of the iron is taken as the parameter for the reduced model. In this problem, as the potential contact surface on the slab is quite larger than the actual position of contact, the contact pressure snapshots display large changes in contact position, shown in Fig. 12. The same is reflected in *spy* pattern of the contact pressure snapshot matrix, shown in Fig. 13

**Details of FE model:** A finite element model of the ironing problem is created using a structured quad mesh for both the iron and the slab. Like the Hertz problem, the displacement field is approximated using linear shape functions and contact pressure using piecewise-constant shape functions centred at the surface nodes. To better demonstrate the inseparability issues, two meshes: coarse and fine, are considered for iron and slab.

<table><thead><tr><th></th><th>Iron</th><th>Slab</th></tr></thead><tbody><tr><td>Coarse Mesh</td><td>30 × 30</td><td>20 × 100</td></tr><tr><td>Fine Mesh</td><td>60 × 60</td><td>40 × 200</td></tr></tbody></table>

The compactness metric, which is equivalent to decay of truncation error as discussed before, is shown in Fig. 14 for a coarse and fine mesh, created using 128 snapshots. The contrasting part of the two graphs is the green curve representing the dual orth compactness, as this curve vanishes at rank 100 for
<figure><img src="image_12.png" /><figcaption>Figure 11: Ironing problem: 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>

<figure><img src="image_13.png" /><figcaption>Figure 12: Sample snapshots of contact pressure for ironing problem</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 13: A typical sparsity pattern of dual snapshots for ironing problem</figcaption></figure>
the coarse mesh, but not for the fine mesh. This is natural, since coarse mesh has only 100 nodes on the surface of the slab and thus allowing a maximum rank of 100, whereas the fine mesh has around 200 nodes. Since the orthogonally computed dual (dual orth) compactness is also a measure for determining the true rank of the snapshot matrix, the green curve naturally vanishes at a rank of 100. However, the most interesting part of this figure is the blue curve representing the truncation computed using cone-projection. Even for the coarse mesh, whose true rank is bounded by the number of surface nodes, i.e. 100, the dual cone compactness doesn't decay well beyond this rank. This is another interesting behaviour of the low-rank approach since the dual solution is sought in a subcone that involves positivity constraints, further degrading compactness.

Apart from compactness, another way to demonstrate the dimensionality of the dual subcone is by using the projection error of nested training sets, which will be referred to as nested error subsequently. To compute the nested error, snapshots are computed in a nested set of points in the parametric space, so that n-th level training set is a subset of the (n + 1)-th level training set. In each nested level, 2^n + 1 points are uniformly distributed in the parametric space [0, L]. This means, every alternate point in (n + 1)-th level is the mid-point of two consecutive points in n-th level. The nested error H(n) can then be defined as:

$$H(n) = \frac{\|\mathbf{S}_{\mathrm{tr}}^{n+1} - \Pi_{\mathbf{S}_{\mathrm{tr}}^n}(\mathbf{S}_{\mathrm{tr}}^{n+1})\|_F}{\|\mathbf{S}_{\mathrm{tr}}^{n+1}\|_F} \quad (19)$$

where the symbols Π and S_{tr} carry same meaning defined in Sec. 4 and ||·||_F indicates the Frobenius norm. The snapshot matrix for n-th level is indicated by S_{tr}^n. The nested error is computed by projecting the (n + 1)-th level snapshots on the subspace/subcone spanned by n-th level snapshots. The slope of nested errors computed using orthogonal and cone projections indicate that the offline stage will require a very large number of snapshots to explore the dual subcone, which is another indication of its high-dimensionality, as shown in Fig. 15. The effects of mesh size are similar to those seen in compactness, as nested error using dual cone projection does not decay even after training set is full rank in case of coarse mesh.

Generalization ability and Specificity metrics are plotted in Fig. 16. The generalization ability shows a similar trend as that of nested error, as they are similar quantities measuring reconstruction errors outside the training set. The primal basis outperforms its dual counterparts in this metric, as expected. The specificity metric is quite higher for dual quantities than the primal ones, which is because random shapes generated using the RB are quite different.

# Conclusions and Perspectives

**Conclusions:** In this article, the limitations of low-rank ROMs to contact mechanics problems are discussed. Demonstrations of this approach show moderate reconstruction errors for parametric points outside the training set. To
<figure><img src="image_15.png" /><figcaption>Figure 14: Compactness for reduced bases of ironing problem. Left: Coarse Mesh and Right: Fine Mesh</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 15: Projection error of nested level $n + 1$ snapshots on dual RB of level $n$. Left: Coarse Mesh and Right: Fine Mesh</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 16: Metrics for ironing problem (for fine mesh). Left: Generalization ability and Right: Specificity. The fill regions indicate $\pm 1\sigma$ interval of corresponding metrics.</figcaption></figure>
explore the possible sources of this error, we focus on the key assumption of the low-rank hypothesis i.e. the solution subspace is low-dimensional. It turns out that the Lagrange Multiplier associated with the inequality constraints, i.e. the contact pressure, does not admit a low-rank subspace. Qualitative and quantitative assessments are provided to support this argument. The high-dimensional nature of dual subspace/subcone is demonstrated using various metrics that show high-dimensionality of contact pressure compared to the relatively low-dimensionality of the displacement field.

**Perspectives:** To circumvent the limitations due to the lack of low-rank structure, several potential approaches need to be explored. One possibility is the usage of dictionary-based approximation with a relatively large number of snapshots spanning the dual subspace efficiently. Another possibility is resorting to either a non-linear transformation, like the change of variable approach; potentially transforming the dual snapshots in a way that they lie in a low-dimensional subspace.

## 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 [26] and SciPy [27] for scientific computations and Matplotlib [28] for creating plots. Meshes were created using Gmsh [29].

## References

[1] S. McBane and Y. Choi. Component-wise reduced order model lattice-type structure design. *Computer Methods in Applied Mechanics and Engineering*, 381:113813, 2021. doi:10.1016/j.cma.2021.113813.

[2] F. D. Freitas, R. Pulch, and J. Rommes. Fast and accurate model reduction for spectral methods in uncertainty quantification. *International Journal for Uncertainty Quantification*, 6(3):271–286, 2016.

[3] D. González, F. Masson, F. Poulhaon, A. Leygue, E. Cueto, and F. Chinesta. Proper generalized decomposition based dynamic data driven inverse identification. *Mathematics and Computers in Simulation*, 82(9): 1677–1695, 2012. doi:10.1016/j.matcom.2012.04.001.

[4] M. Kärcher, Z. Tokoutsi, M. A. Grepl, and K. Veroy. Certified reduced basis methods for parametrized elliptic optimal control problems with distributed controls. *Journal of Scientific Computing*, 75(1):276–307, Aug 2017. doi:10.1007/s10915-017-0539-z.
[5] D. Modesto, S. Zlotnik, and A. Huerta. Proper generalized decomposition for parameterized helmholtz problems in heterogeneous and unbounded domains: Application to harbor agitation. *Computer Methods in Applied Mechanics and Engineering*, 295:127 – 149, 2015. doi:10.1016/j.cma.2015.03.026.

[6] A. Quarteroni, A. Manzoni, and F. Negri. *Reduced Basis Methods for Partial Differential Equations*, volume 92 of UNITEXT. Springer International Publishing, Cham, 2016. ISBN 978-3-319-15430-5. doi:10.1007/978-3-319-15431-2.

[7] J. S. Hesthaven, G. Rozza, and B. Stamm. *Certified Reduced Basis Methods for Parametrized Partial Differential Equations*. SpringerBriefs in Mathematics. Springer International Publishing, Cham, 2016. ISBN 978-3-319-22469-5. doi:10.1007/978-3-319-22470-1.

[8] A. Quarteroni and G. Rozza. Numerical solution of parametrized navier–stokes equations by reduced basis methods. *Numerical Methods for Partial Differential Equations*, 23(4):923–948, 2007. doi:10.1002/num.20249.

[9] 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.

[10] O. Burkovska, B. Haasdonk, J. Salomon, and B. Wohlmuth. Reduced basis methods for pricing options with the black–scholes and heston models. *SIAM Journal on Financial Mathematics*, 6(1):685–712, 2015. doi:10.1137/140981216.

[11] 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.

[12] 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.

[13] 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.

[14] 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.

[15] 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.
[16] J. Oden and N. Kikuchi. Theory of variational inequalities with applications to problems of flow through porous media. *International Journal of Engineering Science*, 18(10):1173-1284, 1980. doi:10.1016/0020-7225(80)90111-1.

[17] D. Kinderlehrer and G. Stampacchia. *An Introduction to Variational Inequalities and Their Applications*. Society for Industrial and Applied Mathematics, 2000. doi:10.1137/1.9780898719451.

[18] M. Sofonea and A. Matei. *Variational Inequalities with Applications*, volume 18 of *Advances in Mechanics and Mathematics*. Springer New York, New York, NY, 2009. doi:10.1007/978-0-387-87460-9.

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

[20] 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.

[21] 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.

[22] 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.

[23] 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.

[24] 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.

[25] D. D. Lee and H. S. Seung. Learning the parts of objects by non-negative matrix factorization. *Nature*, 401:788-791, oct 1999. doi:10.1038/44565.

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

[27] 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.
[28] J. D. Hunter. Matplotlib: A 2d graphics environment. *Computing in Science & Engineering*, 9(3):90–95, 2007. doi:10.1109/MCSE.2007.55.

[29] 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.

## About the authors

Kiran Sagar Kollepara, José V. Aguado, Luisa Silva, José M. Navarro

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

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

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

Kiran Sagar Kollepara, Yves Le Guennec

IRT Jules Verne

Chemin du Chaffault, 44340 Bouguenais, France

e-mail: {, yves.le-guennec}@irt-jules-verne.fr

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