# Gradient-based bilevel optimization for multi-penalty Ridge regression
through matrix differential calculus

Gabriele Maroni, Loris Cannelli, Dario Piga

( IDSIA Dalle Molle Institute for Artificial Intelligence USI-SUPSI, Via la
Santa 1, CH-6962 Lugano-Viganello, Switzerland.

November 23 2023)

###### Abstract

Common regularization algorithms for linear regression, such as LASSO and
Ridge regression, rely on a regularization hyperparameter that balances the
tradeoff between minimizing the fitting error and the norm of the learned
model coefficients. As this hyperparameter is scalar, it can be easily
selected via random or grid search optimizing a cross-validation criterion.
However, using a scalar hyperparameter limits the algorithm’s flexibility and
potential for better generalization. In this paper, we address the problem of
linear regression with $\ell_{2}$-regularization, where a different
regularization hyperparameter is associated with each input variable. We
optimize these hyperparameters using a gradient-based approach, wherein the
gradient of a cross-validation criterion with respect to the regularization
hyperparameters is computed analytically through matrix differential calculus.
Additionally, we introduce two strategies tailored for sparse model learning
problems aiming at reducing the risk of overfitting to the validation data.
Numerical examples demonstrate that our multi-hyperparameter regularization
approach outperforms LASSO, Ridge, and Elastic Net regression. Moreover, the
analytical computation of the gradient proves to be more efficient in terms of
computational time compared to automatic differentiation, especially when
handling a large number of input variables. Application to the identification
of over-parameterized Linear Parameter-Varying models is also presented.

##  1 Introduction

###  1.1 Research problem

Linear regression problems often require regularization of the model
parameters for two main reasons: ($i$) to ensure better numerical conditioning
in cases of multicollinearity (i.e., when input variables are highly
correlated); ($ii$) to enhance the generalization capabilities of the
estimated model, as regularization drives the model coefficients towards zero,
thereby reducing model complexity and preventing overfitting.

The most widely used regularization techniques for linear regression are the
LASSO [18] and the Ridge algorithm [9]. These techniques penalize the
$\ell_{1}$ and the $\ell_{2}$ norm of the model coefficients, respectively. In
the simplest case, the importance of the penalty term is governed by a scalar
hyperparameter, which controls the balance between fitting the training data
and enforcing a penalty on the norm of the model coefficients. The use of a
singular scalar hyperparameter simplifies its selection, usually accomplished
through random or grid search within a cross-validation framework. However,
relying solely on a single hyperparameter presents certain limitations.
Specifically, it applies the same regularization strength to all input
variables without considering the varied importance or relevance of different
features. Furthermore, regularization inherently introduces bias by shrinking
model coefficients towards zero, and relying on just one hyperparameter limits
the flexibility in managing the bias-variance tradeoff. In cases involving
inherently sparse models, where only a few input variables are relevant, it
may be more suitable to encourage only the model parameters associated with
non-relevant variables to approach zero, rather than penalizing the other
model coefficients as much or at all.

While addressing these limitations by introducing multiple hyperparameters,
potentially assigning one to each input feature, seems theoretically
appealing, it comes with challenges. The curse of dimensionality becomes
apparent, and the risk of overfitting increases. The search space for
hyperparameters grows significantly, making grid search or random search
strategies practically infeasible. Even more scalable gradient-free
algorithms, such as Bayesian optimization [17], are generally limited to
handling a maximum of roughly 50 hyperparameters. Paradoxically, while
regularization aims to prevent overfitting, introducing numerous
hyperparameters, such as a regularization coefficient for each feature, can
inadvertently lead to overfitting the validation set during hyperparameter
tuning.

###  1.2 Contribution

This paper focuses on linear regression with an $\ell_{2}$ penalty on the norm
of the model coefficients, employing multiple hyperparameters. We consider the
general scenario where each input variable is associated with its unique
hyperparameter. Essentially, our problem can be seen as a generalization of
Ridge regression with multiple hyperparameters. Throughout this paper, we will
refer to this problem as _multi-ridge regression_ , adopting the terminology
in [20]. To select these hyperparameters, we employ a gradient-based approach
that optimizes a cross-validation criterion. The gradient with respect to the
regularization hyperparameters is analytically computed using _matrix
differential calculus_. This computation leverages efficient matrix-vector
multiplication available in modern scientific computing libraries such as
NumPy, PyTorch, and JAX.

As outlined in [1], optimizing a large number of hyperparameters through
cross-validation can lead to overfitting to the validation set. To mitigate
this risk, we introduce two strategies that are particularly suitable for
solving sparse model learning problems.

In the numerical examples, we will demonstrate that using multiple
regularization hyperparameters enhances the performance compared to single-
hyperparameter Ridge regression. Furthermore, it outperforms both LASSO and
Elastic Net, even in the context of sparse model learning problems.

To ensure the reproducibility of our results and to encourage further
contributions to the field, we have made a PyTorch implementation of our
proposed approach. This implementation is tailored with a scikit-learn
compatible interface, ensuring seamless integration into existing workflows.
Both the implementation and the results presented in this paper are freely
available on the GitHub repository at the following link:
https://github.com/gabribg88/Multiridge.

###  1.3 Related works

The problem of Ridge regression with multiple penalties was early proposed in
[9]. A gradient-based approach to optimize the algorithm’s hyperparameter was
later addressed in the seminal work [1]. The author here proposes to solve the
inefficient step of computing the gradient of the problem variables with
respect to the hyperparameters through a combination of Cholesky decomposition
and backpropagation. The approach introduced in [1], later implemented in
[11], has been then extensively explored in subsequent studies. For instance,
the works [6, 8, 13, 15, 3] propose efficient solutions to compute the
computationally expensive gradient step mentioned in [1]. These solutions
employ approximators (e.g., Neumann series as demonstrated in [13, 8]) or
iterative procedures (e.g., conjugate gradient in [15], reverse-hyper gradient
in [6], or proximal non-smooth techniques in [3]). Although these approaches
provide approximate solutions, they have proven to be satisfactory,
particularly for machine learning tasks, as discussed [15, 4].

Another research avenue that has gained momentum from [1] is the exploration
of optimizing regularization hyperparameters by targeting specific machine-
learning problems. By selecting a particular machine learning application,
researchers can frame the task of finding the optimal variable fitting, while
simultaneously determining the best configuration for the hyperparameters, as
a bilevel optimization problem. This allows for tailored approaches to
effectively address the challenge. Interesting tools have been proposed, among
others, in [2] for regularized Support Vector Machines, in [5] for log-linear
models, in [7] for the group-LASSO setting, and in [10] for Ridge and LASSO
regularizations. An alternative approach to nested optimization for selecting
the multi-ridge hyperparameters is based on maximization of the marginal
likelihood and discussed in [20].

In relation to existing literature, we present a custom approach to linear
regression with $\ell_{2}$ regularization. As previously highlighted, this
approach enables us to compute efficiently an analytical expression for the
gradient of a cross-validation-based criterion with respect to the
hyperparameters, using matrix differential calculus.

###  1.4 Paper organization

The paper is organized as follows. Section 2 introduces the mathematical
notation used throughout the paper. Section 3 presents the multi-
hyperparameter regression problem addressed in our contribution, along with an
intuitive motivation of the advantages of using multiple regularization
hyperparameters over a single one, especially in sparse model learning
problems. In Section 4, we present the analytical computation of the gradient
of the cross-validation criterion with respect to the regularization
hyperparameters. Detailed derivations can be found in the Appendix to
facilitate a smooth reading flow. Section 5 discusses the two strategies we
propose to reduce the risk of overfitting to the validation set. Section 6
provides numerical examples showcasing the effectiveness of our approach. One
of these examples is focused on the identification of over-parameterized
Linear Parameter-Varying models, a topic previously tackled in the literature
using LASSO-like and kernel-based algorithms (see, e.g., [19, 16, 12]).
Finally, Section 7 offers conclusions and directions for future research.

##  2 Notation

We denote with $\otimes$ the Kronecker product and with $\odot$ the
elementwise matrix multiplication. Given a column vector $v\in\mathbb{R}^{n}$,
$\mathrm{Diag}(v)$ is an $n\times n$ diagonal matrix having the elements of
$v$ itself on the main diagonal and 0s otherwise; conversely, given a matrix
$V\in\mathbb{R}^{n\times n}$, $\mathrm{diag}(V)$ is the column vector
$v\in\mathbb{R}^{n}$ collecting the entries of the main diagonal of $V$.
$\mathrm{vec}(\cdot)$ denotes the vectorization operator that stacks the
columns of a matrix vertically to form a column vector. Given any matrix
$A\in\mathbb{R}^{m\times n}$, $\left\|A\right\|_{F}$ denotes its Frobenius
norm, and $A^{\prime}$ the transpose of $A$. Finally, we indicate with $I_{N}$
the identity matrix of size $N$.

Given a matrix $U\in\mathbb{R}^{m\times n}$ and a real-valued function
$\phi:\mathbb{R}^{m\times n}\rightarrow\mathbb{R}$, $U\mapsto\phi(U)$, the
gradient of $\phi$ w.r.t. the vectorized version of $U$ ($\mathrm{vec}(U)$) is
denoted as $\nabla_{U}\phi$, with $\nabla_{U}\phi\in\mathbb{R}^{mn}$. If
$U\in\mathbb{R}^{m\times m}$ is a diagonal matrix, with diagonal elements
$u=\mathrm{diag}(U)\in\mathbb{R}^{m}$, $\nabla_{u}\phi\in\mathbb{R}^{m}$ is
the gradient of the function $\phi$ w.r.t. the diagonal elements of the matrix
$U$.

Given a finite-dimensionals set $\mathcal{S}$, $|\mathcal{S}|$ is its
cardinality.

##  3 Problem Formulation

###  3.1 Parametric estimate: setting

In a parametric supervised learning setting, we aim to select a function $f$
that maps input features $x\in\mathcal{X}$ to output targets $y\in\mathcal{Y}$
from a set $\mathcal{F}$ of candidate functions, parameterized by a vector of
parameters $\theta$, such that each value assumed by the vector of parameters
$\theta$ corresponds to a different function in the set $\mathcal{F}$. In this
work we focus on regression problems, and we assume
$\mathcal{X}=\mathbb{R}^{D}$, with $D$ being the number of input features; and
$\mathcal{Y}=\mathbb{R}$ in the case of simple regression with a real-valued
target, or $\mathcal{Y}=\mathbb{R}^{M}$ in the case of multi-task/target
regression, with $M$ denoting the number of real-valued target variables.

In a real-world setting, we are given a dataset of $N$ input-output pairs
$\mathcal{D}=\left\\{\left(x^{(i)},y^{(i)}\right)\right\\}_{i=1}^{N}$ which
are assumed to be independently drawn from the same unknown probability
distribution $\mathbb{P}_{x,y}$. In principle, we would like to select the
values of $\theta$ (and hence the function $f_{\theta}\in\mathcal{F}$) which
minimize the expected risk, defined as the expected value
$\mathbb{E}_{x,y}[\ell(y,f_{\theta})]$ of some scalar loss function
$\ell(y,f_{\theta})$, which is a measure of the discrepancy between the target
and the model’s prediction. In general we do not have access to the data
distribution $\mathbb{P}_{x,y}$, so we approximate the expected risk with an
empirical risk defined as the average loss over the dataset $\mathcal{D}$:

$L(\theta)=\frac{1}{N}\sum_{i=1}^{N}\ell\left(y^{(i)},f_{\theta}\left(x^{(i)}\right)\right).$
(1)

In solving regression problems, a regularization term $R(\theta,\lambda)$ is
typically added to control the complexity of the estimated function
$f_{\theta}$. Such a regularization term is parameterized by a vector of
hyperparameters $\lambda$ that regulate the trade-off between the fit and the
regularization term.

The overall training criterion $C(\theta,\lambda)$ is then defined as:

$C(\theta,\lambda)=L(\theta)+R(\theta,\lambda),$ (2)

and, for a given value of $\lambda$, the problem of finding the optimal
parameterized function is typically formulated as the following optimization
problem:

$\hat{\theta}(\lambda)=\arg\min_{\theta}C(\theta,\lambda).$ (3)

Among different possibilities, the hyperparameters $\lambda$ can be chosen by
optimizing an evaluation criterion $E(\hat{\theta}(\lambda))$, such as
minimizing the average loss over a new dataset that was not used in evaluating
$L(\theta)$. For example, in the celebrated $K$-fold cross-validation
strategy, the original dataset $\mathcal{D}$ is divided into $K$ partitions
$\mathcal{D}_{1},\ldots,\mathcal{D}_{K}$, and, consequently, the evaluation
criterion takes the form:

$\displaystyle
E\left(\hat{\theta}(\lambda)\right)=\frac{1}{K}\sum_{k=1}^{K}L_{k}\left(\hat{\theta}^{(\backslash
k)}(\lambda)\right),$ (4)

where $L_{k}$ is the average loss over the subset $\mathcal{D}_{k}$ (namely,
validation fold). Specifically:

$\displaystyle L_{k}\left(\hat{\theta}^{(\backslash
k)}(\lambda)\right)=\frac{1}{|\mathcal{D}_{k}|}\sum_{(x^{(i)},y^{(i)})\in\mathcal{D}_{k}}\ell\left(y^{(i)},f_{\theta^{(\backslash
k)}(\lambda)}\left(x^{(i)}\right)\right),$ (5)

and $\hat{\theta}^{(\backslash k)}(\lambda)$ represents the estimate of the
model parameters $\theta$ obtained from constructing $C(\theta,\lambda)$ using
all remaining folds $\mathcal{D}_{j}$ (with $j=1,\ldots,K$ and $j\neq k$),
i.e.,

$\displaystyle\hat{\theta}^{(\backslash k)}(\lambda)=$ (6)
$\displaystyle\arg\min_{\theta}\frac{1}{|\bigcup_{j\neq
k}\mathcal{D}_{j}|}\sum\limits_{\begin{subarray}{c}\left(x^{(i)},y^{(i)}\right)\\\
\in\bigcup_{j\neq
k}\mathcal{D}_{j}\end{subarray}}\ell\left(y^{(i)},f_{\theta}\left(x^{(i)}\right)\right)+R(\theta,\lambda).$

Thus, the hyperparameters $\lambda$ are selected by solving the minimization
problem:

$\displaystyle\lambda^{\star}=\arg\min_{\lambda}E\left(\hat{\theta}(\lambda)\right).$
(7)

The hyperparameter optimization problem is thus a bilevel optimization
problem, wherein the objective function of the outer optimization problem (7)
depends on the solution of the inner optimization task (3).

###  3.2 Linear-in-the parameter models and quadratic loss

In this work, we consider:

  * •

as the set of candidate functions, linear-in-the-parameter functions of the
form

$\displaystyle f_{\mathbf{\Theta}}(x)=x^{\prime}\mathbf{\Theta},$ (8)

where $\mathbf{\Theta}\in\mathbb{R}^{D\times M}$ is the coefficient matrix and
$f_{\mathbf{\Theta}}(x)$ is a row vector of size $M$ whose its $i$-th
component represents the model prediction for the $i$-th output;

  * •

a quadratic fitting loss:

$\displaystyle\ell\left(y^{(i)},f_{\mathbf{\Theta}}\left(x^{(i)}\right)\right)=\frac{1}{2}\left\|y^{(i)}-\left(x^{(i)}\right)^{{}^{\prime}}\mathbf{\Theta}\right\|^{2};$
(9)

  * •

a quadratic regularization term:

$\displaystyle
R(\mathbf{\Theta},\lambda)=\frac{1}{2}\left\|\Lambda\mathbf{\Theta}\right\|_{F}^{2},$
(10)

where $\Lambda=Diag(\lambda)$, and $\lambda\in\mathbb{R}^{D}$ is the vector of
hyperparameters.

Thus, the general expressions (6) and (7) of the $K$-fold cross-validation
specialized, respectively, for the above choices (8)-(9)-(10), become:

$\displaystyle\hat{\mathbf{\Theta}}^{(\backslash k)}(\Lambda)=$ (11)
$\displaystyle\arg\min_{\mathbf{\Theta}}\frac{1}{2|\bigcup_{j\neq
k}\mathcal{D}_{j}|}\sum\limits_{\begin{subarray}{c}\left(x^{(i)},y^{(i)}\right)\\\
\in\bigcup_{j\neq
k}\mathcal{D}_{j}\end{subarray}}\left\|y^{(i)}-\left(x^{(i)}\right)^{{}^{\prime}}\mathbf{\Theta}\right\|^{2}+\frac{1}{2}\left\|\Lambda\mathbf{\Theta}\right\|_{F}^{2},$

and

$\displaystyle\Lambda^{\star}=$ (12)
$\displaystyle\arg\min_{\Lambda}\frac{1}{2K}\sum_{k=1}^{K}\frac{1}{|\mathcal{D}_{k}|}\sum\limits_{(x^{(i)},y^{(i)})\in\mathcal{D}_{k}}\left\|y^{(i)}-\left(x^{(i)}\right)^{{}^{\prime}}\hat{\mathbf{\Theta}}^{(\backslash
k)}(\Lambda)\right\|^{2}.$

To ease readability, without any loss of generality, in the rest of the paper
we assume that all cross-validation folds have equal size, namely:
$|\mathcal{D}_{1}|=\ldots=|\mathcal{D}_{K}|=N_{V}=\frac{N}{K}$. This implies
that the amount of samples used in the validation step (see (12)) of each
cycle of the cross-validation is $\frac{N}{K}$, while $N_{T}=(K-1)\frac{N}{K}$
is the number of samples used to estimate the model parameters
$\hat{\mathbf{\Theta}}^{(\backslash k)}(\Lambda)$ in (11).

To leverage the matrix differential calculus approach, presented in Section 4,
and to enhance notation clarity, we introduce the following matrix notation.
Let $X_{\backslash k}\in\mathbb{R}^{N_{T}\times D}$ and $Y_{\backslash
k}\in\mathbb{R}^{N_{T}\times M}$ be matrices stacking in their rows,
respectively, the input samples $x^{(i)}$ and the output samples $y^{(i)}$,
with $\left(x^{(i)},y^{(i)}\right)\in\bigcup_{j\neq k}\mathcal{D}_{j}$.
Similarly, let $X_{k}\in\mathbb{R}^{N_{V}\times D}$ and
$Y_{k}\in\mathbb{R}^{N_{V}\times M}$ be matrices stacking in their rows,
respectively, the input samples $x^{(i)}$ and the output samples $y^{(i)}$,
with $\left(x^{(i)},y^{(i)}\right)\in\mathcal{D}_{k}$.

Based on the introduced notation, the model parameter estimate (11) can be
rewritten in the following compact Tikhonov regularization form:

$\displaystyle\hat{\mathbf{\Theta}}^{(\backslash k)}(\Lambda)$
$\displaystyle=\arg\min_{\mathbf{\Theta}}\,C(\mathbf{\Theta},\Lambda)$
$\displaystyle=\arg\min_{\mathbf{\Theta}}\,L(\mathbf{\Theta})+R(\mathbf{\Theta},\Lambda)$
$\displaystyle=\arg\min_{\mathbf{\Theta}}\,\frac{1}{2N_{T}}\left\|Y_{\backslash
k}-X_{\backslash
k}\mathbf{\Theta}\right\|_{F}^{2}+\frac{1}{2}\left\|\Lambda\mathbf{\Theta}\right\|_{F}^{2},$
(13)

while the optimal hyperparameters $\Lambda$ in (12) can be expressed as:

$\displaystyle\Lambda^{\star}$
$\displaystyle=\arg\min_{\Lambda}\,E\left(\hat{\mathbf{\Theta}}(\Lambda)\right)$
$\displaystyle=\arg\min_{\Lambda}\,\frac{1}{K}\sum_{k=1}^{K}L_{k}\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}(\Lambda)\right)$
$\displaystyle=\arg\min_{\Lambda}\,\frac{1}{K}\sum_{k=1}^{K}\frac{1}{2N_{V}}\left\|Y_{k}-X_{k}\hat{\mathbf{\Theta}}^{(\backslash
k)}(\Lambda)\right\|_{F}^{2}.$ (14)

###  3.3 Motivation

The main motivation behind the optimization of multiple hyperparameters stems
from the well-known fact that conventional regularization techniques, such as
LASSO, Ridge regression, or Elastic Net, reduce the variance of the estimate
while concurrently introducing bias in the estimation of model parameters, by
shrinking them closer to zero. In a intuitive sense, employing multiple
hyperparameters enables variance reduction of the estimate variance, while
upholding an unbiased estimation of the model parameters.

In order to intuitively explain the concept, let us consider an illustrative
example where the output variable $y$ is scalar and generated according to:

$y=x^{\prime}\mathbf{\Theta}_{\textrm{o}}+e,$ (15)

where $e$ is white noise, independent of the input $x$, and
$\mathbf{\Theta}_{\textrm{o}}\in\mathbb{R}^{D}$ denotes the “true” (unknown)
parameter vector. Let us take a scalar $\bar{\lambda}\in\mathbb{R}$,
$\bar{\lambda}\neq 0$, and let $\mathbf{\Theta}_{\textrm{o},j}$ be the $j$-th
element of $\mathbf{\Theta}_{\textrm{o}}$. If the true parameter vector
$\mathbf{\Theta}_{\textrm{o}}$ is assumed to be sparse, then a hyperparameter
matrix $\Lambda\in\mathbb{R}^{D\times D}$ with diagonal elements

$\displaystyle\lambda_{j}=\left\\{\begin{array}[]{ll}\bar{\lambda}&\textrm{ if
}\mathbf{\Theta}_{\textrm{o},j}=0\\\ 0&\textrm{ if
}\mathbf{\Theta}_{\textrm{o},j}\neq 0\end{array}\right.,\quad j=1,\ldots,D,$
(18)

would regularize (thus shrinking towards zero) the model parameters
$\mathbf{\Theta}_{j}$ such that the corresponding true parameters are
$\mathbf{\Theta}_{\textrm{o},j}=0$. On the other hand, it would not shrink
towards zero the model parameters $\mathbf{\Theta}_{j}$ such that
$\mathbf{\Theta}_{\textrm{o},j}\neq 0$. Hence, this type of regularization
does not introduce bias into the model while it reduces the variance of the
estimate compared to non-regularized approaches. Furthermore, for
$\bar{\lambda}\rightarrow\infty$, $\mathbf{\Theta}_{j}\rightarrow 0$, and thus
input variable selection is also performed without introducing model bias.

Nevertheless, it is important to note, as also discussed in [1], that using
multiple hyperparameters $\lambda$ introduces greater adaptability in the
average loss $L_{k}\left(\hat{\theta}^{(\backslash k)}(\lambda)\right)$ in
(5). More flexibility could potentially increase the risk of overfitting to
the validation data employed for hyperparameter selection. In Section 5, we
discuss two strategies aimed at mitigating the risk of overfitting in the case
where the true underlying parameter matrix $\mathbf{\Theta}_{\textrm{o}}$ is
sparse.

##  4 Gradient-based multi-hyperparameter optimization

In this section we show how to compute the gradient $\nabla_{\Lambda}E$ of the
evaluation criterion $E$ (appearing as a cost function of the outer
optimization problem (3.2)) with respect to the hyperparameter matrix
$\Lambda$. This allows us to compute the optimal hyperparameters
$\Lambda^{\star}$ through any gradient-based numerical optimization algorithm,
such as Vanilla (stochastic) gradient descent; ADAM, Adagrad; RMSProp; etc.

The estimation of the model parameters $\hat{\mathbf{\Theta}}^{(\backslash
k)}(\Lambda)$ in the Tikhonov regularized form (3.2) has the following
analytical solution:

$\displaystyle\hat{\mathbf{\Theta}}^{(\backslash
k)}(\Lambda)=\left(X_{\backslash k}^{\prime}X_{\backslash
k}+N_{T}\Lambda\Lambda\right)^{-1}X_{\backslash k}^{\prime}Y_{\backslash k},$
(19)

where the transpose of $\Lambda$ is omitted in the equation above since
$\Lambda$ is diagonal, thus symmetric.

The gradient $\nabla_{\Lambda}E$ can be readily computed via automatic
differentiation (e.g., using the deep learning Python libraries such as
PyTorch or TensorFlow). Nevertheless, employing backpropagation to solve the
Tikhonov regularization problem (3.2) was observed to be inefficient when
handling a growing number of features, as also discussed in the numerical
example in Section 6.2. For this reason, an analytical strategy for gradient
computation is adopted in this paper. The proposed approach leverages the
framework of matrix differential calculus [14], enabling the derivation of
gradients in a matrix form, thereby exploiting the efficient matrix-vector
multiplication features available in scientific computing libraries such as
NumPy, PyTorch, and JAX.

The following expression for the gradient $\nabla_{\Lambda}E$ is derived:
$\displaystyle\nabla_{\Lambda}E$
$\displaystyle=-\frac{N_{T}}{KN_{V}}\sum_{k=1}^{K}\mathrm{vec}\left(\Lambda
B_{k}+B_{k}\Lambda\right),$ (20a) with
$B_{k}=A_{k}^{\prime}X_{k}^{\prime}R_{k}\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}\right)^{\prime}$, $A_{k}=\left(X_{\backslash k}^{\prime}X_{\backslash
k}+N_{T}\Lambda\Lambda\right)^{-1}$, and
$R_{k}=X_{k}\hat{\mathbf{\Theta}}^{(\backslash k)}(\Lambda)-Y_{k}$.

Since $\Lambda$ is constrained to be diagonal, at each update step of
gradient-descent, only the gradient of the evaluation criterion $E$ with
respect to the diagonal elements $\lambda$ of the hyperparameter matrix
$\Lambda$ is of practical interest, and its expression is given by:

$\displaystyle\nabla_{\lambda}E$
$\displaystyle=-\frac{N_{T}}{KN_{V}}\sum_{k=1}^{K}\mathrm{diag}\left(\Lambda
B_{k}+B_{k}\Lambda\right).$ (20b)

All the details behind the derivation of the gradients (20a) and (20b) are
reported in Appendix A.

##  5 Reducing the risk of overfitting

Assuming that the true parameter matrix $\mathbf{\Theta}_{\textrm{o}}$ is
sparse, this section introduces two strategies with the objective of guiding
the hyperparameters $\Lambda$ towards the form outlined in (18). These
strategies aim to mitigate the potential of overfitting the hyperparameters
$\Lambda$ during the cross-validation stage.

###  5.1 Optimal-solution augmentation

The strategy discussed in this subsection is based on the observation that any
non-zero scaling of the hyperparameter matrix $\Lambda$ does not affect its
“optimal” structure in (18). Thus, given: i) true model parameters
$\mathbf{\Theta}_{\textrm{o}}$, and ii) the optimal structure of
$\Lambda^{\star}$ from (18), i) and ii) implies that the term
$\left\|\gamma\Lambda^{\star}\mathbf{\Theta}_{\textrm{o}}\right\|_{F}^{2}$,
equals $0$ for any scaling parameter $\gamma\in\mathbb{R}$ ($\gamma\neq 0$).
Nevertheless, the model parameters $\mathbf{\Theta}_{j}$ associated with zero
elements of the true parameters $\mathbf{\Theta}_{\textrm{o},j}$ (or
equivalently associated with non-zero hyperparameter $\lambda^{\star}_{j}$)
remain subject to regularization towards zero.

To promote the sparse structure of $\Lambda^{\star}$ in (18), we consider a
finite set $\Gamma$ of non-zero scaling parameters $\gamma$. The optimal
hyperparameters $\Lambda^{\star}$ are then computed as follows:

$\displaystyle\Lambda^{\star}=\arg\min_{\Lambda}\frac{1}{2N}\frac{1}{|\Gamma|}\sum_{k=1}^{K}\sum_{\gamma\in\Gamma}\left\|Y_{k}-X_{k}\hat{\mathbf{\Theta}}^{(\backslash
k)}(\gamma\Lambda)\right\|_{F}^{2},$ (21)

with model parameters given by:

$\displaystyle\hat{\mathbf{\Theta}}^{(\backslash k)}(\gamma\Lambda)=$ (22)
$\displaystyle\arg\min_{\mathbf{\Theta}}\frac{1}{2N_{T}}\sum\limits_{\begin{subarray}{c}(x^{(i)},y^{(i)})\\\
\in\bigcup_{j\neq
k}\mathcal{D}_{j}\end{subarray}}\left\|y^{(i)}-\left(x^{(i)}\right)^{\prime}\mathbf{\Theta}\right\|^{2}+\frac{1}{2}\left\|\gamma\Lambda\mathbf{\Theta}\right\|_{F}^{2}.$

This approach bears similarity to data augmentation, a well-established
technique often employed in deep learning tasks such as image recognition.
Data augmentation involves introducing synthetic images (generated through
operations like, for instance, rotation or adjustments in illumination) to the
training dataset. This helps enforce desirable properties such as invariance
to rotations and illumination. Instead of adding artificial data, our approach
ensures that the hyperparameters $\Lambda$ are not just adjusted to match the
data, but rather ensure a small fitting error
$\left\|Y_{k}-X_{k}\hat{\mathbf{\Theta}}^{(\backslash
k)}(\gamma\Lambda)\right\|_{F}^{2}$ in (21), for any $\gamma\in\Gamma$, as
anticipated by the theoretical optimal hyperparameters $\Lambda^{\star}$ in
(18). In other words, the optimal model parameters
$\mathbf{\Theta}^{(\backslash k)}(\gamma\Lambda)$ are expected to be invariant
to any non-zero scaling factor $\gamma$. To further highlight the similarity
with data augmentation, it is interesting to note that newly constructed
artificial model parameters $\hat{\mathbf{\Theta}}^{(\backslash
k)}(\gamma\Lambda)$ (for each $\gamma\in\Gamma$) are used in the computation
of $\Lambda$ in (21).

The practical implementation of the regularization strategy discussed above
requires minor changes to (19), (20a) and (20b), which become, respectively:

$\displaystyle\hat{\mathbf{\Theta}}^{(\backslash
k)}(\gamma\Lambda)=\left(X_{\backslash k}^{\prime}X_{\backslash
k}+N_{T}\gamma^{2}\Lambda\Lambda\right)^{-1}X_{\backslash
k}^{\prime}Y_{\backslash k}\,,$ (23)
$\displaystyle\nabla_{\Lambda}E=-\gamma^{2}\frac{N_{T}}{KN_{V}}\frac{1}{|\Gamma|}\sum_{k=1}^{K}\sum_{\gamma\in\Gamma}\mathrm{vec}\left(\Lambda\tilde{B}_{k}+\tilde{B}_{k}\Lambda\right)\,,$
(24)
$\displaystyle\nabla_{\lambda}E=-\gamma^{2}\frac{N_{T}}{KN_{V}}\frac{1}{|\Gamma|}\sum_{k=1}^{K}\sum_{\gamma\in\Gamma}\mathrm{diag}\left(\Lambda\tilde{B}_{k}+\tilde{B}_{k}\Lambda\right)\,.$
(25)

with
$\tilde{B}_{k}=\tilde{A}_{k}^{\prime}X_{k}^{\prime}R_{k}\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}\right)^{\prime}$, $\tilde{A}_{k}=\left(X_{\backslash
k}^{\prime}X_{\backslash k}+N_{T}\gamma^{2}\Lambda\Lambda\right)^{-1}$, and
$R_{k}=X_{k}\hat{\mathbf{\Theta}}^{(\backslash k)}(\Lambda)-Y_{k}$.

###  5.2 Regularization in validation

The second strategy for reducing the risk of overfitting stems from the same
considerations discussed in the previous paragraph. Indeed, in the case of a
sparse true parameter vector $\mathbf{\Theta}_{\textrm{o}}$, we would aim to
guide the hyperparameter $\Lambda$ to take the structure in (18), rather than
simply minimizing the fitting herror
$\left\|Y_{k}-X_{k}\hat{\mathbf{\Theta}}^{(\backslash
k)}(\Lambda)\right\|_{F}^{2}$ of the $k$-th validation fold. To this aim, the
validation loss function in (3.2) is modified as follow, by introducing a
regularization term:

$\displaystyle\Lambda^{\star}=\arg\min_{\Lambda}$
$\displaystyle\frac{1}{2N}\sum_{k=1}^{K}\left\|Y_{k}-X_{k}\hat{\mathbf{\Theta}}^{(\backslash
k)}(\Lambda)\right\|_{F}^{2}+$
$\displaystyle+\frac{1}{2}\mu\sum_{k=1}^{K}\left\|\Lambda\hat{\mathbf{\Theta}}^{(\backslash
k)}(\Lambda)\right\|_{F}^{2},$ (26)

where $\mu\in\mathbb{R}$, with $\mu>0$, acts as an additional hyperparameter.
Thus, unlike the “optimal-solution augmentation” strategy discussed in Section
5.1, this approach requires the tuning of an additional (albeit scalar)
regularization hyperparameter, hence the requirement for an additional hyper-
validation dataset.

The practical implementation of the regularization strategy discussed in this
subsection requires the derivation of the gradient of the regularization term:

$\displaystyle
Q=\frac{1}{2}\mu\sum_{k=1}^{K}\left\|\Lambda\hat{\mathbf{\Theta}}^{(\backslash
k)}(\Lambda)\right\|_{F}^{2},$ (27)

with respect to $\Lambda$ (or equivalently $\lambda$), whose analytical
expressions are given by:

$\displaystyle\nabla_{\Lambda}Q$
$\displaystyle=\mu\sum_{k=1}^{K}\mathrm{vec}\left(D_{k}\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}\right)^{\prime}-N_{T}\left(\Lambda G_{k}+G_{k}\Lambda\right)\right),$
(28a) $\displaystyle\nabla_{\lambda}Q$
$\displaystyle=\mu\sum_{k=1}^{K}\mathrm{diag}\left(D_{k}\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}\right)^{\prime}-N_{T}\left(\Lambda G_{k}+G_{k}\Lambda\right)\right),$
(28b)

with $D_{k}=\Lambda\hat{\mathbf{\Theta}}^{(\backslash k)}$ and
$G_{k}=A_{k}^{\prime}\Lambda^{\prime}D_{k}\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}\right)^{\prime}$.

Derivations of the gradients $\nabla_{\Lambda}Q$ (resp. $\nabla_{\lambda}Q$)
in (28a) (resp. (28b)) are based on similar differential matrix calculus
arguments used to compute (20a) and (20b), and are reported in Appendix B.

##  6 Examples

In this section, we present numerical examples to demonstrate the
effectiveness of the multi-hyperparameter optimization approach. In this
context, we conduct a comparative analysis between this approach and
established benchmark techniques, such as LASSO, Ridge and Elastic Net
regressions, all of which are readily available through the Python _scikit-
learn_ package. Comparison in terms of final model performance and
computational complexity is discussed in Section 6.1 and 6.2, respectively.
Parametric identification of _Linear Parameter-Varying_ (LPV) dynamical
systems with an unknown model structure is addressed in Section 6.3.

Model performance is measured on a test set (used neither for training nor for
validation) through the R2 index, defined as:

$\displaystyle\textrm{R}^{2}=\max\left\\{0;1-\frac{\sum_{i=1}^{N_{\rm
test}}\left\|y^{(i)}-\hat{y}^{(i)}\right\|^{2}}{\sum_{i=1}^{N_{\rm
test}}\left\|y^{(i)}-\bar{y}\right\|^{2}}\right\\},$ (29)

where $\hat{y}^{(i)}$ is the estimated model output for the $i$-th sample,
$\bar{y}$ is the sample mean of the outputs over the test dataset, and $N_{\rm
test}$ is length of the test dataset. In assessing algorithms’ performance, we
use test datasets larger (up 100x) than the training set. Although this is not
realistic in practice, it allows us to have a fair and statistically
significant evaluation and comparison of different algorithms.

Output data (including test data) is corrupted by an additive noise $e$. The
effect of the noise is measured in terms of the Signal-to-Noise Ratio (SNR),
defined as:

$\mathrm{SNR}=10\log{\left(\frac{\left\|y-e\right\|^{2}}{\left\|e\right\|^{2}}\right)}.$
(30)

For the multi-ridge case discussed in the paper, numerical optimization is
performed through gradient descent, with learning rate chosen manually. In the
first and third examples and in all algorithms discussed in this section,
5-fold cross-validation is implemented to select the regularization
hyperparameters. Once the hyperparameters are chosen, models are re-estimated
on the complete training dataset. In the second example, a simple holdout
validation is used, as the focus is only on the numerical computation of the
gradient $\nabla_{\lambda}E$ and not on assessing model performance.

All computations are performed on a server with two 64-core AMD EPYC 7742
Processors, 256 GB of RAM, and 4 Nvidia RTX 3090 GPUs. The hardware resources
of the server have been limited to 10 CPU threads and 1 GPU.

The software has been developed in PyTorch with a scikit-learn compatible
interface. All the codebase, along with documentation, is available in the
GitHub repository https://github.com/gabribg88/Multiridge.

###  6.1 Example 1: Sensitivity to the number of features

In this numerical example, we investigate the sensitivity of different
regularization strategies with respect to the number of features, while
keeping the number of data points and the information ratio (i.e., the ratio
between the total number of input features available in the dataset and the
number of informative features) constant.

Since the number of features effectively determines the complexity of the
model, we expect that as the ratio of input features to data points grows, the
model’s capacity to capture not only the underlying pattern but also the noise
within the training data, increases, potentially resulting in overfitting. As
the overfitting increases, the impact of the regularization term becomes more
significant. This provides a basis for comparing the effectiveness of various
regularization strategies.

$14$ distinct conditions are simulated, with the number of features $D$
incrementally ranging from $100$ to $1400$ with a step of $100$. To enhance
the statistical significance of the results, we replicated each condition $10$
times, utilizing different random seeds for each replication. In each
individual experiment, a total of $N_{train}=10^{3}$ training data samples and
$N_{test}=10^{5}$ test samples were generated, according to the following
linear data generation procedure:

$\displaystyle y=x^{\prime}\tilde{\mathbf{\Theta}}_{\rm o}+e,$ (31)

where each component of the feature vector $x\in\mathbb{R}^{D}$ is drawn from
a Normal distribution $\mathcal{N}(0,1)$. The noise $e$ is also Normally
distributed, with $e\sim\mathcal{N}(0,\sigma_{e}^{2})$ and variance
$\sigma_{e}^{2}$ chosen to ensure an SNR of approximately $20$ dB. The actual
coefficient vector $\tilde{\mathbf{\Theta}}_{\rm o}$ is a sparse version of a
vector $\mathbf{\Theta}_{\rm o}$, whose components are drawn from a uniform
distribution $\mathcal{U}_{\,[-50,50]}$ between $-50$ and $50$. The relation
between $\tilde{\mathbf{\Theta}}_{\rm o}$ and $\mathbf{\Theta}_{\rm o}$ is
thus expressed as:

$\displaystyle\tilde{\mathbf{\Theta}}_{{\rm{o}},j}=\left\\{\begin{array}[]{ll}\mathbf{\Theta}_{{\rm{o}},j}&\textrm{
if }j\in\mathcal{I}\\\ 0&\textrm{ if }j\notin\mathcal{I}\end{array}\right.,$
(34)

where $\mathcal{I}\subseteq\\{1,2,\dots,D\\}$ represents the index set of
informative features. The cardinality of $\mathcal{I}$ is
$|\mathcal{I}|=\left\lfloor DI_{frac}\right\rfloor$, with
$\left\lfloor\cdot\right\rfloor$ being the floor operator. The fraction of
informative features $I_{frac}$ is set to $0.5$ and kept constant throughout
this example. Each element $j\in\mathcal{I}$ is uniformly drawn at random,
without replacement, in $\\{1,2,\dots,D\\}$.

In each experiment, the benchmark learning algorithms underwent a systematic
training process with the following sequential stages: i) standardization of
both features and target variable; ii) hyperparameter optimization through
5-fold cross-validation;, iii) re-training of the models on the complete
training dataset; iv) evaluation of algorithms’ performances on a test set.
Hyperparameter optimization varies depending on the regularization algorithm,
as detailed below:

  * •

Ridge regression: a grid-search was conducted on a logarithmically spaced grid
of $10^{3}$ values for the $l_{2}$-regularization strength. The grid values
ranged from $10^{-3}$ to $10^{6}$.

  * •

LASSO: a grid-search was performed using a logarithmically spaced grid of
$10^{3}$ values for the $l_{1}$-regularization strength. The grid values
spanned from $10^{-5}$ to $10^{2}$.

  * •

Elastic Net regression: a random-search was executed with $10^{3}$ repetitions
on pairs of: logarithmically spaced values of the regularization strength
ranging from $10^{-5}$ to $10^{3}$, and ratio between $\ell_{1}$ and
$\ell_{2}$ regularization spanning linearly from $0$ to $1$.

For multi-ridge regression, the hyperparameters $\Lambda$ are optimized
through a gradient-based approach, and initialized as an identity matrix
$I_{D}$. The learning rate is set to $350$ across all scenarios, with a decay
rate of $0.999$ applied at each epoch over the course of $300$ training
epochs.

Figure 1: Example 1: model performance _vs_ number of input features for
different regularization algorithms. Median over $10$ repetitions (dots),
along with $5$-th and $95$-th percentile (bars) are depicted.

Results are presented in Fig. 1, illustrating the sensitivity of model
performance with respect to the number of input features for Ridge, LASSO,
Elastic Net, and the proposed multi-ridge regression. For less than $300$
features, all methods deliver performance comparable to the oracle model with
the true parameter vector $\tilde{\mathbf{\Theta}}_{\rm o}$. However, as the
number of input features grows, each method exhibits a gradual decline in
performance due to overfitting.

Multi-ridge regression consistently outperforms the other benchmark
regularization strategies. Notably, simple Ridge regression records the lowest
performance. This outcome is expected for the given problem with a sparse true
parameter vector. Thus, algorithms that promote sparse solutions, such as
LASSO, are more suitable. Nevertheless, even though multi-ridge regression
employs $\ell_{2}$ regularization instead of $\ell_{1}$, it still outperforms
LASSO. This enhanced performance stems from our method’s flexibility in
assigning a distinct regularization strength to each feature, based on its
inherent informativeness.

###  6.2 Example 2: Computational performances

In this numerical example we show the advantages in terms of both
computational efficiency (measured in terms of GPU time) and numerical
accuracy (assessed by quantifying the norm of the difference between the
analytic calculation and the automatic differentiation, as implemented within
the PyTorch framework) during the computation of the gradient in (20b). These
benefits become increasingly pronounced as the number of features in the
simulated dataset grows. More specifically, we constructed $50$ distinct
scenarios, each characterized by a progressively growing number of features,
with logarithmically spaced values. For every scenario, we conducted $10$
experiments, each repeated with varying random seeds to enhance the
statistical significance of the results.

In all the experiments, a fixed number of $N_{train}=10^{3}$ training data
samples were generated. Each data sample was composed of an input-output pair
$(x,y)$, where $x\in\mathbb{R}^{D}$ and $y\in\mathbb{R}^{M}$, with $M=10$.
Generated samples were subsequently partitioned into training fold data and
validation fold data, according to an 80:20 ratio. This is equivalent to
simple holdout validation instead of $K$ fold.

The input features $x$, the targets $y$, as well as the hyperparameters
$\lambda$, are generated according to a Normal distribution, with zero mean
and unitary variance.

(a) Double-precision floating-point format (float64)

(b) Single-precision floating-point format (float32)

Figure 2: Example 2: Results of computational performance. GPU time required
to compute the gradient $\nabla_{\lambda}E$ (upper panels); norm of the
difference between analytic computation of the gradient $\nabla_{\lambda}E$
and its numerical computation through automatic differentiation in PyTorch,
denoted as $\overline{\nabla_{\lambda}E}$ (lower panels).

The experiment is conducted in 2 different scenarios. Firstly, with variables
represented in double-precision floating-point format (float64), and secondly,
with variables represented in single-precision floating-point format
(float32). The number of features ranges from $10$ to $10^{4}$.

Results are illustrated in Fig. 2. The upper panel of the figure displays the
GPU time required for the computation of the gradient
$\overline{\nabla_{\lambda}E}$. The orange line corresponds to the time taken
by the PyTorch autograd, whereas the blue line indicates the time consumed by
the analytic calculation. As the number of features increases, the analytic
calculation showcases a significant enhancement in computational efficiency,
obtaining computational speed more than an order of magnitude faster w.r.t.
automatic differentiation.

The lower panel of the figure displays the norm of the difference between the
analytic calculation and the computation through automatic differentiation. In
the scenario where variables are represented in double-precision floating-
point format, this difference is negligible. Nonetheless, when variables are
represented in single-precision floating-point format, this difference becomes
noteworthy as the number of features increases.

###  6.3 LPV identification

As a last example, we consider the identification of an LPV system described
by the AutoRegressive eXogenous (ARX) form:

$\displaystyle y(k)=$
$\displaystyle\;0.5\cos{(p(k))}y(k-2)-0.1\sin^{2}{(p(k))}y(k-3)$
$\displaystyle+(\cos{(p(k))}-\sin{(p(k))})u(k-2)$
$\displaystyle+3\sin{(p(k))}u(k-3)+e(k),$ (35)

where $u(k)$, $y(k)$, and $p(k)$ represent the measurement of the input,
output, and scheduling variable at time step $k$, respectively, while $e(k)$
is a random Gaussian white noise.

A small-size dataset $\mathcal{D}$ of size $N=50$ is gathered by exciting the
system (6.3) with a white noise input sequence $\\{u(k)\\}_{k=1}^{N}$, with
$u(k)$ drawn from a Gaussian distribution $\mathcal{N}(0,1)$. The scheduling
variable $p(k)$ is also drawn from a Gaussian distribution, with
$p(k)\sim\mathcal{N}(0,\pi)$. The amplitude of the noise $e(k)$ is chosen such
that the SNR is about $14$ dB (namely, the ratio between the power of the
noise and the power of the noise-free output is approximately $4$%).

In estimating the LPV system (6.3), we assume to be in a setting where the
system structure is not known, and we consider an over-parameterized LPV-ARX
model of the form:

$\displaystyle\hat{y}(k)=\sum_{j=1}^{n_{a}}{a}_{j}(p(k))y(k-j)+\sum_{j=1}^{n_{b}}{b}_{j}(p(k))u(k-j),$

with $n_{a}=n_{b}=30$. Furthermore, the dependence of the coefficient
functions ${a}_{j}$ and ${b}_{j}$ ($j=1,\ldots,30$) on the scheduling
parameter $p$ is parameterized in terms of a linear combination of a finite
number of _a-priori_ chosen basis functions, as follows:

$\displaystyle a_{j}(p(k))=$
$\displaystyle\sum_{s=1}^{n_{\alpha}}a_{j,s}\psi_{s}(p(k)),$ $\displaystyle
i=1,\ldots,n_{a},$ (36a) $\displaystyle b_{j}^{\mathrm{}}(p(k))=$
$\displaystyle\sum_{s=1}^{n_{\alpha}}b_{j,s}\psi_{s}(p(k)),$ $\displaystyle
j=1,\ldots,n_{b},$ (36b)

where $n_{\alpha}=8$, $\\{a_{j,s}\in\mathbb{R}\\}_{s=1}^{n_{\alpha}}$, and
$\\{b_{j,s}\in\mathbb{R}\\}_{s=1}^{n_{\alpha}}$ are the unknown constant
parameters to be estimated, and $\psi_{s}(p(k))$ is the $s$-th element of the
vector function:

$\displaystyle\psi_{s}(p)=\left[1\ p\ p^{2}\ p^{3}\ \sin(p)\ \cos(p)\
\sin^{2}(p)\ \cos^{2}(p)\right].$ (37)

Summarizing, we have to estimate $n_{\alpha}(n_{a}+n_{b})=480$ parameters from
a dataset of length $N=50$. In such a setting, regularization is expected to
play a fundamental role for the estimation process.

For multi-ridge regression, the initial conditions of the hyperparameters
$\lambda$ are extracted from the solution of a LASSO problem. More
specifically, if the model parameter $\mathbf{\Theta}_{i}$ estimated by the
LASSO is zero, then the related hyperparameter $\lambda_{j}$ is set to $10$,
otherwise $\lambda_{j}=1$. This means that, in the initialization of
$\lambda$, a 10x higher penalization is given to the model parameters
estimated as zero by the LASSO. The optimal-solution augmentation strategy
discussed in Section 5.1 is implemented to reduce the risk of overfitting to
the validation datasets, by minimizing the evaluation criterion in (21) for
$\gamma\in\Gamma=\\{0.5,\ 1,\ 2\\}$.

The performance of all the tested algorithms is evaluated on a test data of
length $N_{\rm test}=3000$. Figure 3 shows the boxplots of the different
approaches over $200$ Monte Carlo simulations, with: input, noise and
scheduling variable independently regenerated at each run. The figure shows
that LASSO, Elastic Net and the multi-hyperparameter approach proposed in this
paper provide satisfactory performance, achieving a median R2 value of $0.88$,
$0.86$ and $0.91$, respectively. However, this example also highlights how
simple Ridge regression does not have enough flexibility in its regularization
term. In contrast, multi-ridge regression, having a much larger flexibility,
surpasses LASSO’s performances too, even when dealing with a sparse structure
in the “true” parameter vector.

Figure 3: LPV identification: boxplots of the R2 index over 200 Monte Carlo
runs. Medians: Least Squares (0.18); Ridge (0.21); LASSO (0.88); Elastic Net
(0.86); multi-ridge (0.91).

##  7 Conclusions

Multi-penalty Ridge regression increases the degrees of freedom of the
regression algorithm, thus allowing to enhance model’s performance compared to
widely adopted regularization techniques such as LASSO, Ridge, and Elastic Net
regression. The optimization of multiple regularization hyperparameters can be
achieved through gradient descent. Here, the analytical expression of the
gradient for a cross-validation criterion with respect to the regularization
hyperparameters can be derived using matrix differential calculus, avoiding
the need for automatic differentiation and backpropagation. As shown in the
numerical examples, this offers numerical benefits, especially when dealing
with a large number of features.

The approach is specifically tailored to $\ell_{2}$ regularization and takes
advantage of the fact that the optimal model parameters can be determined
analytically. However, ongoing research aims to extend the matrix differential
calculus-based approach presented in this paper to other types of problems
where optimal model parameters cannot be expressed analytically and must be
determined through iterative numerical optimization. These types of problems
include LASSO with multiple hyperparameters or nested optimization problems
arising in meta-learning, where models are trained on the task of learning
themselves, often with the aim of rapidly adapting to new tasks using only a
small-size dataset.

## References

  * [1] Y. Bengio.  Gradient-based optimization of hyperparameters.  Neural computation, 12(8):1889–1900, 2000.
  * [2] K. P. Bennett, J. Hu, X. Ji, G. Kunapuli, and J.-S. Pang.  Model selection via bilevel optimization.  The 2006 IEEE International Joint Conference on Neural Network Proceedings, pages 1922–1929, 2006.
  * [3] Q. Bertrand, Q. Klopfenstein, M. Massias, M. Blondel, S. Vaiter, A. Gramfort, and J. Salmon.  Implicit differentiation for fast hyperparameter selection in non-smooth convex learning.  The Journal of Machine Learning Research, 23(1):6680–6722, 2022.
  * [4] M. Blondel, Q. Berthet, M. Cuturi, R. Frostig, S. Hoyer, F. Llinares-López, F. Pedregosa, and J.-P. Vert.  Efficient and modular implicit differentiation.  Advances in neural information processing systems, 35:5230–5242, 2022.
  * [5] C. B. Do, C.-S. Foo, and A. Ng.  Efficient multiple hyperparameter learning for log-linear models.  Advances in neural information processing systems, 20, 2007.
  * [6] L. Franceschi, M. Donini, P. Frasconi, and M. Pontil.  Forward and reverse gradient-based hyperparameter optimization.  International Conference on Machine Learning, pages 1165–1173, 2017.
  * [7] J. Frecon, S. Salzo, and M. Pontil.  Bilevel learning of the group lasso structure.  Advances in neural information processing systems, 31, 2018.
  * [8] R. Hataya, J. Zdenek, K. Yoshizoe, and H. Nakayama.  Meta approach to data augmentation optimization.  Proceedings of the IEEE/CVF Winter Conference on Applications of Computer Vision, pages 2574–2583, 2022.
  * [9] A. E. Hoerl and R. W. Kennard.  Ridge regression: Biased estimation for nonorthogonal problems.  Technometrics, 12(1):55–67, 1970.
  * [10] K. Kunisch and T. Pock.  A bilevel optimization approach for parameter learning in variational models.  SIAM Journal on Imaging Sciences, 6(2):938–983, 2013.
  * [11] S. Latendresse and Y. Bengio.  Linear regression and the optimization of hyper-parameters.  Web Published Writings, pages 1–7, 2000.
  * [12] V. Laurain, R. Tóth, D. Piga, and M. A. H. Darwish.  Sparse RKHS estimation via globally convex optimization and its application in LPV-IO identification.  Automatica, 115:108914, 2020.
  * [13] J. Lorraine, P. Vicol, and D. Duvenaud.  Optimizing millions of hyperparameters by implicit differentiation.  International conference on artificial intelligence and statistics, pages 1540–1552, 2020.
  * [14] J. R. Magnus and H. Neudecker.  Matrix differential calculus with applications in statistics and econometrics.  John Wiley & Sons, 2019.
  * [15] F. Pedregosa.  Hyperparameter optimization with approximate gradient.  International conference on machine learning, pages 737–746, 2016.
  * [16] D. Piga and R. Tóth.  LPV model order selection in an LS-SVM setting.  In 52nd IEEE Conference on Decision and Control, pages 4128–4133, 2013.
  * [17] J. Snoek, H. Larochelle, and R. P. Adams.  Practical Bayesian optimization of machine learning algorithms.  Advances in neural information processing systems, 25, 2012.
  * [18] R. Tibshirani.  Regression shrinkage and selection via the lasso.  Journal of the Royal Statistical Society: Series B (Methodological), 58(1):267–288, 1996.
  * [19] R. Tóth, H. Hjalmarsson, and C. R. Rojas.  Order and structural dependence selection of LPV-ARX models revisited.  In 51st IEEE Conference on Decision and Control (CDC), pages 6271–6276, 2012.
  * [20] M. A. van de Wiel, M. M. van Nee, and A. Rauschenberger.  Fast cross-validation for multi-penalty high-dimensional ridge regression.  Journal of Computational and Graphical Statistics, 30(4):835–847, 2021.

## Appendix

##  Appendix A Derivation of (20a) and (20b)

By the _first identification theorem_ [14, Chapter 9, Section 5], the
relationship between the differential and the gradient of the evaluation
criterion $E$ is given by:

$\displaystyle dE$
$\displaystyle=\left(\nabla_{\Lambda}E\right)^{\prime}d\mathrm{vec}\Lambda,$
(38a) or equivalently: $\displaystyle dE$
$\displaystyle=\left(\nabla_{\lambda}E\right)^{\prime}d\lambda.$ (38b)

Thus, from (38), both $\nabla_{\Lambda}E$ and $\nabla_{\lambda}E$ can be
readily obtained, as discussed in the following. By differentiating both sides
of (4) and leveraging linearity of the differential operator, we obtain:

$\displaystyle
dE=d\left(\frac{1}{K}\sum_{k=1}^{K}L_{k}\right)=\frac{1}{K}\sum_{k=1}^{K}dL_{k}.$
(39)

To derive $dL_{k}$ (and thus $dE$), $L_{k}$ is rewritten in terms of the trace
operator as follows:

$\displaystyle L_{k}$
$\displaystyle=\frac{1}{2N_{V}}\left\|R_{k}\right\|_{F}^{2}=\frac{1}{2N_{V}}\mathrm{tr}\left(R_{k}^{\prime}R_{k}\right)\,,$
(40)

where $R_{k}=X_{k}\hat{\mathbf{\Theta}}^{(\backslash k)}(\Lambda)-Y_{k}$. By
differentiating (40):

$\displaystyle dL_{k}$
$\displaystyle=\frac{1}{2N_{V}}d\mathrm{tr}\left(R_{k}^{\prime}R_{k}\right)$
$\displaystyle=\frac{1}{2N_{V}}\mathrm{tr}\left[\left(dR_{k}^{\prime}\right)R_{k}+R_{k}^{\prime}dR_{k}\right]$
$\displaystyle=\frac{1}{2N_{V}}\mathrm{tr}\left[\left(dR_{k}\right)^{\prime}R_{k}+R_{k}^{\prime}dR_{k}\right]$
$\displaystyle=\frac{1}{2N_{V}}\mathrm{tr}\left(R_{k}^{\prime}dR_{k}+R_{k}^{\prime}dR_{k}\right)$
$\displaystyle=\frac{1}{N_{V}}\mathrm{tr}\left(R_{k}^{\prime}dR_{k}\right),$

where the following properties of the differential operator are used, for any
matrix $U\in\mathbb{R}^{m\times n}$ and $V\in\mathbb{R}^{m\times p}$:

$\displaystyle d(U^{\prime})=dU^{\prime},$ (42a) $\displaystyle
d\left(UV\right)=\left(dU\right)V+UdV\quad\textit{(Product rule)}.$ (42b)

The computation of $dL_{k}$ in (A) requires to compute the differential
$dR_{k}$. From the definition of $R_{k}$ and the analytical expression of the
model parameters $\hat{\mathbf{\Theta}}$ in (19):

$\displaystyle dR_{k}$
$\displaystyle=d\left(X_{k}\hat{\mathbf{\Theta}}^{(\backslash k)}\right)$ (43)
$\displaystyle=d\left(X_{k}\left(X_{\backslash k}^{\prime}X_{\backslash
k}+N_{T}\Lambda\Lambda\right)^{-1}X_{\backslash k}^{\prime}Y_{\backslash
k}\right)$ $\displaystyle=d\left(X_{k}A_{k}X_{\backslash
k}^{\prime}Y_{\backslash
k}\right)=R_{k}^{\prime}X_{k}\left(dA_{k}\right)X_{\backslash
k}^{\prime}Y_{\backslash k},$

with $A_{k}=\left(X_{\backslash k}^{\prime}X_{\backslash
k}+N_{T}\Lambda\Lambda\right)^{-1}$. By substituting (43) into (A), and from
the cyclic property of the trace operator, we obtain:

$\displaystyle dL_{k}$
$\displaystyle=\frac{1}{N_{V}}\mathrm{tr}\left[R_{k}^{\prime}X_{k}\left(dA_{k}\right)X_{\backslash
k}^{\prime}Y_{\backslash k}\right]$
$\displaystyle=\frac{1}{N_{V}}\mathrm{tr}\left(X_{\backslash
k}^{\prime}Y_{\backslash k}R_{k}^{\prime}X_{k}dA_{k}\right).$ (44)

The differential $dA_{k}$ appearing in (43) can be easily computed starting
from the definition of the definition of $A_{k}$. Specifically, we have:

$\displaystyle dA_{k}$ $\displaystyle=d\left(X_{\backslash
k}^{\prime}X_{\backslash k}+N_{T}\Lambda\Lambda\right)^{-1}$
$\displaystyle=-N_{T}A_{k}d\left(\Lambda\Lambda\right)A_{k},$ (45)

where the following identity is used:

$\displaystyle d\left(W^{-1}\right)=-W^{-1}\left(dW\right)W^{-1}\,,$ (46)

for any invertible square matrix $W\in\mathbb{R}^{m\times m}$. By substituting
(45) into (A), we obtain:

$\displaystyle dL_{k}$
$\displaystyle=-\frac{N_{T}}{N_{V}}\mathrm{tr}\left[X_{\backslash
k}^{\prime}Y_{\backslash
k}R_{k}^{\prime}X_{k}A_{k}d\left(\Lambda\Lambda\right)A_{k}\right]$
$\displaystyle=-\frac{N_{T}}{N_{V}}\mathrm{tr}\left[A_{k}X_{\backslash
k}^{\prime}Y_{\backslash
k}R_{k}^{\prime}X_{k}A_{k}d\left(\Lambda\Lambda\right)\right]$
$\displaystyle=-\frac{N_{T}}{N_{V}}\mathrm{tr}\left[\hat{\mathbf{\Theta}}^{(\backslash
k)}R_{k}^{\prime}X_{k}A_{k}d\left(\Lambda\Lambda\right)\right]\,.$ (47)

Finally, from the product rule and by defining
$B_{k}=A_{k}^{\prime}X_{k}^{\prime}R_{k}\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}\right)^{\prime}$, (A) can be written as:

$\displaystyle dL_{k}$
$\displaystyle=-\frac{N_{T}}{N_{V}}\mathrm{tr}\Big{[}\hat{\mathbf{\Theta}}^{(\backslash
k)}R_{k}^{\prime}X_{k}A_{k}\left(d\Lambda\right)\Lambda$
$\displaystyle\qquad+\hat{\mathbf{\Theta}}^{(\backslash
k)}R_{k}^{\prime}X_{k}A_{k}\Lambda d\Lambda\Big{]}$
$\displaystyle=-\frac{N_{T}}{N_{V}}\mathrm{tr}\left[B_{k}^{\prime}\left(d\Lambda\right)\Lambda+B_{k}^{\prime}\Lambda
d\Lambda\right]$ $\displaystyle=-\frac{N_{T}}{N_{V}}\mathrm{tr}\left[\Lambda
B_{k}^{\prime}\left(d\Lambda\right)+B_{k}^{\prime}\Lambda d\Lambda\right]$
$\displaystyle=-\frac{N_{T}}{N_{V}}\mathrm{tr}\left[\left(\Lambda
B_{k}^{\prime}+B_{k}^{\prime}\Lambda\right)d\Lambda\right]\,.$ (48)

Using the following property: for any matrix $U\in\mathbb{R}^{m\times n}$ and
$V\in\mathbb{R}^{m\times p}$:

$\displaystyle\mathrm{tr}\left(U^{\prime}V\right)=\mathrm{vec}\left(U\right)^{\prime}\mathrm{vec}\left(V\right)\,,$
(49)

(48) becomes:

$\displaystyle dL_{k}$
$\displaystyle=-\frac{N_{T}}{N_{V}}\mathrm{vec}\left(\Lambda
B_{k}+B_{k}\Lambda\right)^{\prime}d\mathrm{vec}\Lambda.$ (50)

By combining (39) and (50), the differential $dE$ can be written as:

$\displaystyle dE=-\frac{N_{T}}{KN_{V}}\sum_{k=1}^{K}\mathrm{vec}\left(\Lambda
B_{k}+B_{k}\Lambda\right)^{\prime}d\mathrm{vec}\Lambda.$ (51)

Thus, from (51) and (38a), we observe:

$\displaystyle\nabla_{\Lambda}E$
$\displaystyle=-\frac{N_{T}}{KN_{V}}\sum_{k=1}^{K}\mathrm{vec}\left(\Lambda
B_{k}+B_{k}\Lambda\right),$ (52)

which proves (20a). In order to prove (20b), we rewrite (48) as:

$\displaystyle dL_{k}$
$\displaystyle=-\frac{N_{T}}{N_{V}}\mathrm{diag}\left(\Lambda
B_{k}+B_{k}\Lambda\right)^{\prime}d\lambda,$ (53)

where we used the following properties:

$\displaystyle\mathrm{tr}\left(W\mathrm{Diag}\left(v\right)\right)=\mathrm{diag}\left(W\right)^{\prime}v\,,$
(54a) $\displaystyle d\mathrm{Diag}(v)=\mathrm{Diag}(dv)\,,$ (54b)

for any square matrix $W\in\mathbb{R}^{m\times m}$ and vector
$v\in\mathbb{R}^{m}$. By combining (39) and (53), the differential $dE$ can be
written as:

$\displaystyle
dE=-\frac{N_{T}}{KN_{V}}\sum_{k=1}^{K}\mathrm{diag}\left(\Lambda
B_{k}+B_{k}\Lambda\right)^{\prime}d\lambda.$ (55)

Thus, from (55) and (38b), we observe:

$\displaystyle\nabla_{\lambda}E$
$\displaystyle=-\frac{N_{T}}{KN_{V}}\sum_{k=1}^{K}\mathrm{diag}\left(\Lambda
B_{k}+B_{k}\Lambda\right),$ (56)

which proves (20b).

##  Appendix B Derivation of (28a) and (28b)

First, let us rewrite the regularization term $Q$ in (27) as:

$\displaystyle Q=\sum_{k=1}^{K}Q_{k},$ (57)

with:

$\displaystyle Q_{k}$
$\displaystyle=\frac{1}{2}\mu\left\|\Lambda\hat{\mathbf{\Theta}}^{(\backslash
k)}\right\|_{F}^{2}=\frac{1}{2}\mu\,\mathrm{tr}\left(\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}\right)^{\prime}\Lambda\Lambda\hat{\mathbf{\Theta}}^{(\backslash
k)}\right)$
$\displaystyle=\frac{1}{2}\mu\,\mathrm{tr}\left(D_{k}^{\prime}D_{k}\right),$
(58)

where $D_{k}=\Lambda\hat{\mathbf{\Theta}}^{(\backslash k)}$ and the dependence
of $\mathbf{\Theta}^{(\backslash k)}$ on $\Lambda$ is omitted to ease reading.
Differentiating both sides of (B) and using the trace and differential
operator properties already adopted in Appendix A we obtain:

$\displaystyle dQ_{k}$
$\displaystyle=\frac{1}{2}\mu\,d\left[\mathrm{tr}\left(D_{k}^{\prime}D_{k}\right)\right]=\frac{1}{2}\mu\,\mathrm{tr}\left[d\left(D_{k}^{\prime}D_{k}\right)\right]$
$\displaystyle=\frac{1}{2}\mu\,\mathrm{tr}\left[\left(dD_{k}^{\prime}\right)D_{k}+D_{k}^{\prime}dD_{k}\right]$
$\displaystyle=\frac{1}{2}\mu\,\mathrm{tr}\left[\left(dD_{k}\right)^{\prime}D_{k}+D_{k}^{\prime}dD_{k}\right]$
$\displaystyle=\frac{1}{2}\mu\,\mathrm{tr}\left(D_{k}^{\prime}dD_{k}+D_{k}^{\prime}dD_{k}\right)$
$\displaystyle=\mu\,\mathrm{tr}\left(D_{k}^{\prime}dD_{k}\right).$ (59)

From the definition of $D_{k}$, (B) can be rewritten as:

$\displaystyle dQ_{k}$
$\displaystyle=\mu\,\mathrm{tr}\left[D_{k}^{\prime}d\left(\Lambda\hat{\mathbf{\Theta}}^{(\backslash
k)}\right)\right]$
$\displaystyle=\mu\,\mathrm{tr}\left[D_{k}^{\prime}\left(d\Lambda\right)\hat{\mathbf{\Theta}}^{(\backslash
k)}+D_{k}^{\prime}\Lambda d\hat{\mathbf{\Theta}}^{(\backslash k)}\right]$
$\displaystyle=\mu\left[\mathrm{tr}\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}D_{k}^{\prime}d\Lambda\right)+\mathrm{tr}\left(D_{k}^{\prime}\Lambda
d\hat{\mathbf{\Theta}}^{(\backslash k)}\right)\right].$ (60)

Let us now focus on the second term in the right side of (B), which is
rewritten by using the same properties adopted in Appendix A, as:

$\displaystyle\mathrm{tr}\left(D_{k}^{\prime}\Lambda
d\hat{\mathbf{\Theta}}^{(\backslash k)}\right)=$
$\displaystyle=\mathrm{tr}\left[D_{k}^{\prime}\Lambda
d\left(\left(X_{\backslash k}^{\prime}X_{\backslash
k}+N_{T}\Lambda\Lambda\right)^{-1}X_{\backslash k}^{\prime}Y_{\backslash
k}\right)\right]$
$\displaystyle=\mathrm{tr}\left[D_{k}^{\prime}\Lambda\left(dA_{k}\right)X_{\backslash
k}^{\prime}Y_{\backslash k}\right]$
$\displaystyle=\mathrm{tr}\left(X_{\backslash k}^{\prime}Y_{\backslash
k}D_{k}^{\prime}\Lambda dA_{k}\right)$
$\displaystyle=-N_{T}\,\mathrm{tr}\left[X_{\backslash k}^{\prime}Y_{\backslash
k}D_{k}^{\prime}\Lambda A_{k}d\left(\Lambda\Lambda\right)A_{k}\right]$
$\displaystyle=-N_{T}\,\mathrm{tr}\left[\hat{\mathbf{\Theta}}^{(\backslash
k)}D_{k}^{\prime}\Lambda A_{k}d\left(\Lambda\Lambda\right)\right]$
$\displaystyle=-N_{T}\,\mathrm{tr}\left[\hat{\mathbf{\Theta}}^{(\backslash
k)}D_{k}^{\prime}\Lambda
A_{k}\left(d\Lambda\right)\Lambda+\hat{\mathbf{\Theta}}^{(\backslash
k)}D_{k}^{\prime}\Lambda A_{k}\Lambda d\Lambda\right].$

By defining
$G_{k}=A_{k}^{\prime}\Lambda^{\prime}D_{k}\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}\right)^{\prime}$, the above equation can be written in the compact form:

$\displaystyle\mathrm{tr}\left(D_{k}^{\prime}\Lambda
d\hat{\mathbf{\Theta}}^{(\backslash k)}\right)$
$\displaystyle=-N_{T}\,\mathrm{tr}\left[G_{k}^{\prime}\left(d\Lambda\right)\Lambda+G_{k}^{\prime}\Lambda
d\Lambda\right]$ $\displaystyle=-N_{T}\,\mathrm{tr}\left[\left(\Lambda
G_{k}^{\prime}+G_{k}^{\prime}\Lambda\right)d\Lambda\right].$

Then, by substituting (B) into (B), we obtain:

$\displaystyle dQ_{k}$
$\displaystyle=\mu\left[\mathrm{tr}\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}D_{k}^{\prime}d\Lambda\right)-N_{T}\,\mathrm{tr}\left(\left(\Lambda
G_{k}^{\prime}+G_{k}^{\prime}\Lambda\right)d\Lambda\right)\right]$
$\displaystyle=\mu\left[\mathrm{tr}\left(\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}D_{k}^{\prime}-N_{T}\left(\Lambda
G_{k}^{\prime}+G_{k}^{\prime}\Lambda\right)\right)d\Lambda\right)\right].$
(63)

Thanks to (49), (B) can be written as:

$\displaystyle dQ_{k}$
$\displaystyle=\mu\,\mathrm{vec}\left(D_{k}\hat{\mathbf{\Theta}}^{(\backslash
k)}-N_{T}\left(\Lambda
G_{k}+G_{k}\Lambda\right)\right)^{\prime}d\mathrm{vec}\Lambda.$ (64)

Based on the same considerations discussed in Appendix A, from (64), and the
definition of $Q$ in (27), we obtain:

$\displaystyle\nabla_{\Lambda}Q$
$\displaystyle=\mu\sum_{k=1}^{K}\mathrm{vec}\left(D_{k}\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}\right)^{\prime}-N_{T}\left(\Lambda G_{k}+G_{k}\Lambda\right)\right),$ (65)

thus proving (28a). Using the matrix properties in (54), (B) can be also
written as:

$\displaystyle dQ_{k}$
$\displaystyle=\mu\,\mathrm{diag}\left(D\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}\right)^{\prime}-N_{T}\left(\Lambda
G_{k}+G_{k}\Lambda\right)\right)^{\prime}d\lambda.$ (66)

Therefore:

$\displaystyle\nabla_{\lambda}Q$
$\displaystyle=\mu\sum_{k=1}^{K}\mathrm{diag}\left(D_{k}\left(\hat{\mathbf{\Theta}}^{(\backslash
k)}\right)^{\prime}-N_{T}\left(\Lambda G_{k}+G_{k}\Lambda\right)\right),$ (67)

thus proving (28b).


# Testing regularity of black holes with X-rays and Gravitational Waves

Swarnim Shashank 111On behalf of: S. Riaz, R. Roy, A. B. Abdikamalov, D.
Ayzenberg, C. Bambi, Z. Zhang, M. Zhou Center for Field Theory and Particle
Physics and Department of Physics, Fudan University,
200438 Shanghai, China

###### Abstract

Physically relevant solutions in general relativity often contain spacetime
singularities, which are typically interpreted as a sign of breakdown of the
theory at high densities/curvatures. Hence, there has been a growing interest
in exploring phenomenological scenarios that describe singularity-free black
holes, gravitational collapses, and cosmological models. We examine the metric
put forth by Mazza, Franzin & Liberati for a rotating regular black hole and
estimate the regularization parameter $l$ based on existing X-ray and
gravitational wave data for black holes. When $l=0$, the solution corresponds
to the singular Kerr solution of general relativity, while a non-zero value of
$l$ yields a regular black hole or a regular wormhole. The analysis reveals
that the available data support a value of $l$ that is close to zero.

##  1 The BH mimicker metric

The metric is stationary, axisymmetric, and asymptotically flat, and its line
element in the Boyer-Lindquist coordinates is written as [2]

${\rm d}s^{2}=-\left(1-\frac{2M\sqrt{r^{2}+l^{2}}}{\Sigma}\right){\rm
d}t^{2}+\frac{\Sigma}{\Delta}{\rm d}r^{2}+\Sigma{\rm d}\theta^{2}\\\
-\frac{4Ma\sin^{2}\theta\sqrt{r^{2}+l^{2}}}{\Sigma}{\rm d}t{\rm
d}\phi+\frac{A\sin^{2}\theta}{\Sigma}{\rm d}\phi^{2}\,$ (1)

where,

$\Sigma=r^{2}+l^{2}+a^{2}\cos^{2}\theta\,,\quad\Delta=r^{2}+l^{2}+a^{2}-2M\sqrt{r^{2}+l^{2}}\,,\quad
A=\left(r^{2}+a^{2}+l^{2}\right)-\Delta a^{2}\sin^{2}\theta\,.$ (2)

The Horizon in the metric is expressed as,

$r_{\rm H\pm}=\sqrt{\left(M\pm\sqrt{M^{2}-a^{2}}\right)^{2}-l^{2}},$ (3)

##  2 Constraints from X-rays

For the case of X-ray data, we analyzed the 2019 NuSTAR observation of the
Galactic black hole EXO 1846–031. The data from NuSTAR observation is quite
good and suitable for such analysis. The constraint on the regularization
parameter we find is,

$l/M<0.49\quad\textnormal{(90\% CL, only statistical error)}.$ (4)

Full details can be found in Ref. [1].

##  3 Gravitational Wave Constraints

Considering the equatorial geodesics of massive particles around a black hole,
one can find the gravitational wave phase in leading order which for the
regularization parameter $l^{2}$ is 2PN.

$\Psi_{\mathrm{GW}}(f)=\Psi_{\mathrm{GW}}^{\mathrm{GR}}(f)-\frac{75}{32}u^{-1/3}\eta^{-4/5}l^{2}+\mathcal{O}[l^{4},u^{0}]\,.$
(5)

From the parameterized post-Einsteinian (ppE) formalism one can find the
relation,

$l^{2}=-\frac{1}{100}\varphi_{4}\delta\varphi_{4}\,,$ (6)

where $\varphi_{4}$ is the PN phase for 2PN and $\delta\varphi_{4}$ is the
correction to the phase used by the LIGO-Virgo Collaboration. Full formalism
is provided in Ref. [1]. Tab. 1 lists the constraints obtained for various GW
events.

Event |  $l/M$ (IMRPhenomPv2) |  $l/M$ (SEOBNRv4P)
---|---|---
GW151226 | $0.50^{+0.26}$ | $0.51^{+0.31}$
GW170608 | $0.49^{+0.26}$ | $0.57^{+0.29}$
GW190408A | $0.58^{+0.39}$ | $0.66^{+0.53}$
GW190412A | $0.57^{+0.27}$ | $0.59^{+0.26}$
GW190630A | $0.48^{+0.38}$ | $0.47^{+0.31}$
GW190707A | $0.44^{+0.24}$ | $0.44^{+0.28}$
GW190708A | $0.47^{+0.28}$ | $0.45^{+0.28}$
GW190720A | $0.45^{+0.29}$ | $0.48^{+0.38}$
GW190728A | $0.41^{+0.36}$ | $--$
GW190828B | $0.46^{+0.33}$ | $0.57^{+0.52}$
Table 1: Median and the $90$th percentile of the parameter $l/M$ from the
significant binary black hole events in GWTC-1 and GWTC-2. We only consider
events with a total redshifted mass below 40 $M_{\odot}$ as we use only the
inspiral phase for our analysis. Note the lower limit for the data is set to
zero.

## Acknowledgments

I would like to thank C. Bambi and S. Riaz for their comments. I acknowledge
the support from the China Scholarship Council (CSC), Grant No. 2020GXZ016646.
The trip for Moriond was partially supported by the China Postdoctoral Science
Foundation, Grant No. 2022M720035.

## References

## References

  * [1] Riaz, S., Shashank, S., Roy, R., Abdikamalov, A., Ayzenberg, D., Bambi, C., Zhang, Z. & Zhou, M. Testing regular black holes with X-ray and GW data. JCAP. 10 pp. 040 (2022)
  * [2] Mazza, J., Franzin, E. & Liberati, S. A novel family of rotating black hole mimickers. JCAP. 4 pp. 082 (2021)