# RANSAC Algorithms for
Subspace Recovery and Subspace Clustering

Ery Arias-Castro  Jue Wang

(University of California, San Diego)

###### Abstract

We consider the RANSAC algorithm in the context of subspace recovery and
subspace clustering. We derive some theory and perform some numerical
experiments. We also draw some correspondences with the methods of Hardt and
Moitra (2013) and Chen and Lerman (2009b).

††The present project was initiated in the context of an Independent Study for
Undergraduates (Math 199). We acknowledge support from the US National Science
Foundation (DMS 1513465).

##  1 Introduction

The Random Sample Consensus (with acronym RANSAC) algorithm of Fischler and
Bolles (1981), and its many variants and adaptations, are well-known in
computer vision for their robustness in the presence of gross errors
(outliers). In this paper we focus on the closely related problems of subspace
recovery and subspace clustering in the presence of outliers, where RANSAC-
type methods are believed to be optimal, yet too costly in terms of
computations when the fraction of inliers is small. Although this is a well-
understood limitation of the RANSAC, we nevertheless establish this rigorously
in the present context. In particular, we derive the performance and
computational complexity of RANSAC for these two problems, and perform some
numerical experiments corroborating our theory and comparing the RANSAC with
other methods proposed in the literature.

###  1.1 The problem of subspace recovery

Consider a setting where the data consist of $n$ points in dimension $p$,
denoted $x_{1},\dots,x_{n}\in\mathbb{R}^{p}$. It is assumed that $m$ of these
points lie on a $d$-dimensional linear subspace $L$ and that the points are
otherwise in general position, which means that the following assumption is in
place:

######  Assumption 1.

A $q$-tuple of data points (with $q\leq p$) is linearly independent unless it
includes at least $d+1$ points from $L$.

(We say that points are linearly dependent if they are so when seen as
vectors.)

The points on $L$ are called inliers and all the other points are called
outliers. This is the setting of subspace recovery without noise. When there
is noise, the points are not exactly on the underlying subspace but rather in
its vicinity. In any case, the goal is to recover $L$, or said differently,
distinguish the inliers from the outliers. See Figure 1(a) for an illustration
in a setting where the subspace is of dimension $d=2$ in ambient dimension
$p=3$. The goal is to recover the subspace $L$ and/or identify the inlier
points.

This problem is intimately related to the problem of robust covariance
estimation, which dates back decades (Maronna, 1976; Huber and Ronchetti,
2009; Tyler, 1987), but has attracted some recent attention. We refer the
reader to the introduction of (Zhang and Lerman, 2014) for a comprehensive
review of the literature, old and new. Subspace recovery in the presence of
outliers, as we consider the problem here, is sometimes referred to a robust
principal components analysis, although there are other meanings in the
literature more closely related to matrix factorization with a low-rank
component (Wright et al., 2009; Candès et al., 2011).

(a) Subspace recovery problem.

(b) Subspace clustering problem.

Figure 1: An illustration of the two settings considered in the paper.

###  1.2 The problem of subspace clustering

Consider a setting where the data consist of $n$ points in dimension $p$,
denoted $x_{1},\dots,x_{n}\in\mathbb{R}^{p}$. It is assumed that $m_{k}$ of
these points lie on a $d_{k}$-dimensional linear subspace $L_{k}$, where
$k=1,\dots,K$ (so that there are $K$ subspaces in total). The remaining points
are in general position:

######  Assumption 2.

A $q$-tuple of data points is linearly independent unless it includes at least
$d_{k}+1$ points from $L_{k}$ for some $k\in\\{1,\dots,K\\}$.

In this setting all the points on one of the subspaces are inliers, and all
the other points are outliers. This is the setting of subspace clustering
without noise. When there is noise, the inliers are not exactly on the
subspaces but in their vicinity. See Figure 1(b) for an illustration in a
setting where there is one subspace of dimension $d_{1}=1$ and two subspaces
of dimension $d_{2}=d_{3}=2$, in ambient dimension $p=3$. The goal is to
cluster $m_{k}$ points to their corresponding $L_{k}$ for all
$k\in\\{1,\dots,K\\}$.

The problem of subspace clustering has applications in computer vision, in
particular, movement segmentation (Vidal, 2011; Vidal et al., 2005).

###  1.3 Contents

In Section 2, we consider the problem of subspace recovery. In Section 3, we
consider the problem of subspace clustering. In both cases, we study a
‘canonical’ RANSAC algorithm, deriving some theory and comparing it with other
methods in numerical experiments. We briefly discuss our results in Section 4.

######  Remark 1 (linear vs affine).

Throughout, we consider the case where the subspaces are linear, although some
applications may call for affine subspaces. (This is for convenience.) Because
of this, we are able to identify a point $x\in\mathbb{R}^{p}$ with the
corresponding vector (sometimes written $x-0$).

##  2 Subspace recovery

We consider the setting of Section 1.1 and use the notation defined there. In
particular, we work under Assumption 1. We consider the noiseless setting for
simplicity.

###  2.1 RANSAC for subspace recovery

We propose a simple RANSAC algorithm for robust subspace recovery. In the
present setting, in particular under Assumption 1, the underlying linear
subspace $L$ (which we assumed is of dimension $d$) is determined by any
$(d+1)$-tuple that comes from that subspace. The algorithm starts by randomly
selecting a $(d+1)$-tuple and checking if this tuple forms a linear subspace
of dimension $d$. If so, the subspace is recovered and the algorithm stops.
Otherwise, the algorithm continues, repeatedly sampling a $(d+1)$-tuple at
random until the subspace is discovered. (Optionally, the algorithm can be
made to stop when a maximum number of tuples has been sampled.) In this
formulation, detailed in Algorithm 1, $d$ is known.

Input : data points $x_{1},\dots,x_{n}\in\mathbb{R}^{p}$; dimension $d$

Output : a linear subspace of dimension $d$ containing at least $d+1$ points

1 repeat

2 randomly select a $(d+1)$-tuple of data points

3until  _the tuple is linearly dependent_ ;

return _the subspace spanned by the tuple_

Algorithm 1 RANSAC (Subspace Recovery)

By design, the procedure is exact. (Again, we are in the noiseless setting. In
a noisy setting, the method can be shown to be essentially optimal.) However,
researchers have shied away from a RANSAC approach because of its time
complexity. We formalize what is in the folklore in the following.

######  Proposition 1.

Algorithm 1 is exact and the number of iterations has the geometric
distribution111 Here we consider the variant of the geometric distribution
that is supported on the positive integers. with success probability
$\theta_{1}:=\binom{m}{d+1}/\binom{n}{d+1}$. Thus the expected number of
iterations is $1/\theta_{1}=\binom{n}{d+1}/\binom{m}{d+1}$, which is of order
$O(n/m)^{d+1}$ when $d$ is held fixed.

Note that each iteration requires on the order of $O(pd^{2})$ operations as it
requires computing the rank of a $p$-by-$(d+1)$ matrix.

###### Proof.

The algorithm sample a $(d+1)$-tuple independently and uniformly at random
until the tuple is linearly dependent. Because of Assumption 1, a
$(d+1)$-tuple is linearly dependent if and only if all the points in the tuple
are from $L$. While there are $\binom{n}{d+1}$ $(d+1)$-tuples in total, only
$\binom{m}{d+1}$ fit the bill, so that the probability of drawing a suitable
tuple is $\theta_{1}=\binom{m}{d+1}/\binom{n}{d+1}$. Because the draws are
independent, the total number of draws until the algorithm stops has the
geometric distribution with success probability $\theta_{1}$.

We know that the mean of this distribution is
$1/\theta_{1}=\binom{n}{d+1}/\binom{m}{d+1}$, and when $d$ is assumed fixed,
while $n$ and $m$ are large, we have

$1/\theta_{1}\sim\frac{n^{d+1}/(d+1)!}{m^{d+1}/(d+1)!}=(n/m)^{d+1}.$

(The reader is invited to verify that this still holds true as long as
$d=o(m^{2})$.) ∎

In applications where the number of outliers is a non-negligible fraction of
the sample (meaning that $n/m$ is not close to 1), the RANSAC’s number of
iterations depends exponentially on the dimension of subspace. This confirms
the folklore, at least in such a setting.

######  Remark 2.

For simplicity, we analyzed the variant of the algorithm where the tuples are
drawn with replacement, so that the worst-case number of iterations is
infinite. However, in practice one should draw the tuples without replacement
(which is equally easy to do in the present setting), as recommended in
(Schattschneider and Green, 2012). For this variant, the worst-case time
complexity is $\binom{n}{d+1}-\binom{n}{d+1}+1$. Moreover, Proposition 1 still
applies if understood as an upper bound. (The number of iterations has a so-
called negative hypergeometric distribution in this case.)

######  Remark 3.

If the dimension $d$ is unknown, a possible strategy is to start with $d=1$,
run the algorithm for a maximum number of iterations, and if no pair of points
is found to be aligned with the origin, move to $d=2$, and continue in that
fashion, increasing the dimension. If no satisfactory tuple is found, the
algorithm would start again at $d=1$. The algorithm will succeed eventually.

###  2.2 The algorithm of Hardt and Moitra for subspace recovery

As we said above, researchers have avoid RANSAC procedures because of the
running time, which as we saw can be prohibitive. Recently, however, Hardt and
Moitra (2013) have proposed a RANSAC-type algorithm that strikes an
interesting compromise between running time and precision.

Their algorithm is designed for the case where the sample size is larger than
the ambient dimension, namely $n>p$. It can be described as follows. It
repeatedly draws a $p$-tuple at random until the tuple is found to be linearly
dependent. When such a tuple is found, the algorithm returns a set of linearly
dependent points in the tuple. See the description in Algorithm 2. A virtue of
this procedure is that it does not require knowledge of the dimension $d$ of
the underlying subspace.

Input : data points $x_{1},\dots,x_{n}\in\mathbb{R}^{p}$

Output : a linear subspace

1 repeat

2 randomly select a $p$-tuple of data points

3until  _the tuple is linearly dependent_ ;

return _the subspace spanned by any subset of linearly dependent points in the
tuple_

Algorithm 2 Hardt-Moitra (Subspace Recovery)

######  Proposition 2.

When $n>p$, Algorithm 2 is exact and its number of iterations has the
geometric distribution with success probability $\theta_{2}:=\sum_{k\geq
d+1}\binom{m}{k}\binom{n-m}{p-k}/\binom{n}{p}$. Thus the expected number of
iterations is $1/\theta_{2}$.

Note that each iteration requires on the order of $O(p^{3})$ operations as it
requires computing the rank of a $p$-by-$p$ matrix.

###### Proof.

With Assumption 1 in place, a $p$-tuple is linearly dependent if and only if
it contains at least $d+1$ points from the subspace $L$. Thus the Repeat
statement above stops exactly when it found a $p$-tuple that contains at least
$d+1$ points. Moreover, also because of Assumption 1, the points within that
tuple that are linear dependent must belong to $L$. Therefore, the algorithm
returns $L$, and is therefore exact.

We now turn to the number of iterations. The number of iterations is obviously
geometric and the success probability is the probability that a $p$-tuple
drawn uniformly at random contains at least $d+1$ points from $L$.
$\theta_{2}$ is that probability. Indeed, it is the probability that, when
drawing $p$ balls without replacement from an urn with $m$ red balls out of
$n$ total, the sample contains at least $d+1$ red balls. In the present
context, the balls are of course the points and the red balls are the points
on the linear subspace. ∎

Hardt and Moitra (2013) analyze their algorithm in a slightly different
setting and with the goal of finding the maximum fraction of outliers that can
be tolerated before the algorithm breaks down in the sense that it does not
run in polynomial time. In particular, they show that, if $m/n\geq d/p$, then
their algorithm has a number of iterations with the geometric distribution
with success probability at least $1/(2p^{2}n)$, so that the expected number
of iterations is bounded by $2p^{2}n$, which is obviously polynomial in
$(p,n)$. In fact, it can be better than that. The following is a consequence
of Proposition 2.

######  Corollary 1.

If, in addition to $n>p$, it holds that $m/n\geq d/p$, with $d/p\leq\tau$ and
$p/n\leq\tau$, for some fixed $\tau<1$, then $\theta_{2}$ is bounded from
below by a positive quantity that depends only on $\tau$. Consequently,
Algorithm 2 has expected number of iterations of order $O(1)$.

###### Proof.

Let $U$ denote a random variable with the hypergeometric distribution with
parameters $(p,m,n)$ described above. Then $\theta_{2}=\mathbb{P}(U\geq d+1)$,
and it depends on $(d,p,m,n)$. We show that $\theta_{2}$ is bounded from below
irrespective of these parameters as long as the conditions are met. Noting
that $\theta_{2}$ is increasing in $d$ and $n$, and decreasing in $p$ and $m$,
it suffices to consider how $\theta_{2}$ varies along a sequence where
$n\to\infty$ and $(d,p,m)$ all varying with $n$ in such a way that
$m/n\to\tau$ and $d/p\to\tau$, as this makes the expected number of iterations
largest. Define

$\displaystyle\mu$ $\displaystyle:=\operatorname{\mathbb{E}}(U)=p(m/n);$
$\displaystyle\sigma^{2}$
$\displaystyle=\operatorname{Var}(U)=p(m/n)(1-m/n)(n-p)/(n-1).$

The condition $m/n\geq d/p$ implies that $\mu\geq d$, and along the sequence
of parameters under consideration, $\sigma\to\infty$. Moreover, along such a
sequence, $Z:=(U-\mu)/\sigma$ is standard normal in the limit, so that

$\displaystyle\theta_{2}=\mathbb{P}(U\geq d+1)$
$\displaystyle=\mathbb{P}\big{(}Z\geq(d+1-\mu)/\sigma\big{)}$
$\displaystyle\geq\mathbb{P}\big{(}Z\geq 1/\sigma\big{)}$
$\displaystyle\to\mathbb{P}\big{(}\mathcal{N}(0,1)\geq 0)=1/2,$

using Slutsky’s theorem in the last line. ∎

###  2.3 Numerical experiments

We performed some small-scale numerical experiments comparing RANSAC (in the
form of Algorithm 1), the Hardt-Moitra (HM) procedure (Algorithm 2), and the
Geometric Median Subspace (GMS) of (Zhang and Lerman, 2014), which appears to
be one of the best methods on the market. (We used the code available on Teng
Zhang’s website.)

Each inlier is uniformly distributed on the intersection of the unit sphere
with the underlying subspace. Each outlier is simply uniformly distributed on
the unit sphere. The result of each algorithm is averaged over 1000 repeats.
Performance is measured by the (first principal) angle between the returned
subspace and the true subspace. (This is to be fair to GMS, as the other two
algorithms are exact.) The results are reported in Table 1.

We performed another set of experiments to corroborate the theory established
in Proposition 1 for the complexity of RANSAC. The results are shown in Figure
2, where each setting has been repeated 1000 times. As expected, as the
dimensionality of the problem increases, RANSAC’s complexity becomes quickly
impractical.

parameters |  | average system time |  | difference in angle
---|---|---|---|---
$(d,p,m,m_{0})$ |  | RANSAC | HM | GMS |  | RANSAC | HM | GMS
$(8,10,100,50)$ |  | .0051 | .0009 | .0947 |  | 0 | 0 | .0341
$(4,10,100,50)$ |  | .0011 | .0006 | .1912 |  | 0 | 0 | 0
$(8,20,100,50)$ |  | .0064 | .0002 | .5008 |  | 0 | 0 | .0184
$(6,10,100,20)$ |  | .0001 | .0006 | .1918 |  | 0 | 0 | 0
$(9,10,100,50)$ |  | .0076 | .0093 | .0117 |  | 0 | 0 | .0411
$(18,20,100,50)$ |  | .6123 | .0013 | .0160 |  | 0 | 0 | .2285
Table 1: Numerical experiments comparing RANSAC, HM, and GMS for the problem
of subspace recovery. As in the text, $d$ is the dimension of the subspace,
$p$ is the ambient dimension, $m$ is the number of inliers, $m_{0}$ is the
number of outliers (so that $n=m+m_{0}$ is the sample size). Figure 2: Average
number of iterations for RANSAC (in the form of Algorithm 1) as a function of
the subspace dimension $d$ and the ratio of sample size $n$ to number of
inliers $m$. The dashed lines are the averages from our simulation while the
lines are derived from theory (Proposition 1).

##  3 Subspace clustering

We consider the setting of Section 1.2 and use the notation defined there. In
particular, we work under Assumption 2. We consider the noiseless setting for
simplicity. We also assume that all subspaces are of same dimension, denoted
$d$ (so that $d_{k}=d$ for all $k$).

###  3.1 RANSAC for subspace clustering

We propose a simple RANSAC algorithm for subspace clustering. As before, any
of the linear subspaces is determined by any $(d+1)$-tuple that comes from
that subspace. The algorithm starts by randomly selecting a $(d+1)$-tuple and
checking if this tuple forms a linear subspace of dimension $d$. If so, one of
the subspaces is recovered and all the points on the subspace are extracted
from the data. Otherwise, the algorithm continues, repeatedly sampling a
$(d+1)$-tuple at random until that condition is met. The algorithm continues
in this fashion until all the $K$ subspaces have been recovered. In this
formulation, detailed in Algorithm 3, both $d$ and $K$ are assumed known.

Input : data points $x_{1},\dots,x_{n}\in\mathbb{R}^{p}$; dimension $d$;
number of subspace $K$

Output : $K$ linear subspaces of dimension $d$, each containing at least $d+1$
points

1 for  _$k=1,\dots,K$_ do

2 repeat

3 randomly select a $(d+1)$-tuple of data points

4 until  _the tuple is linearly dependent_ ;

5 return _the subspace spanned by the tuple_

6 remove the points on that subspace from the data

7 end for

Algorithm 3 RANSAC (Subspace Clustering)

Again, the procedure is exact by design, since we are in the noiseless
setting. Here too, researchers have not embraced RANSAC approaches because of
their running time. We confirm this folklore in the following, where we assume
for simplicity that all subspaces have the same number of points $m$ (so that
$m_{k}=m$ for all $k$).

######  Proposition 3.

Algorithm 3 is exact and the number of iterations is has the distribution of
$I_{1}+\dots+I_{K}$, where the $I$’s are independent and $I_{j}$ has the
geometric distribution with success probability
$(K-j+1)\binom{m}{d+1}/\binom{n-(j-1)m}{d+1}$. This is stochastically bounded
by the negative binomial with parameters $(K,\theta_{1})$. Thus the expected
number of iterations is bounded by $K/\theta_{1}$, which is of order
$O(n/m)^{d+1}$ when $d$ and $K$ are held fixed.

The proof is very similar to that of Proposition 1 and is omitted.

######  Remark 4.

When the dimensions of the subspaces are unknown, a strategy analogous to that
described in Remark 3 is of course possible. When the number of subspaces is
unknown, a stopping rule can help decide whether there remains a subspace to
be discovered. Details are omitted as such an approach, although natural,
could prove complicated.

###  3.2 Adapting the algorithm of Hardt and Moitra for subspace clustering

Algorithm 3 consists in applying Algorithm 1 until a subspace is recovered,
removing the points on that subspace, and then continuing, until all $K$
subspaces are recovered. An algorithm for subspace clustering can be based on
the algorithm of Hardt and Moitra (2013) (Algorithm 2) instead. The resulting
algorithm is suited for the case where $n-(m_{1}+\dots+m_{K})>p$. Based on the
fact that Algorithm 2 has expected number of iterations bounded by $2p^{2}n$,
the resulting algorithm for subspace clustering has expected number of
iterations bounded by $2Kp^{2}n$. See Algorithm 4, where we assume that the
number of subspaces is known, but do not assume that the dimensions of the
subspaces are known (and they do not need to be the same).

Input : data points $x_{1},\dots,x_{n}\in\mathbb{R}^{p}$; number of subspace
$K$

Output : $K$ linear subspaces each with a number of points exceeding its
dimension

1 for  _$k=1,\dots,K$_ do

2 repeat

3 randomly select a $p$-tuple of data points

4 until  _the tuple is linearly dependent_ ;

5 repeat

6 find the smallest number of linearly dependent points in the tuple

7 return _the subspace spanned by these points_

8 remove the points on that subspace from the data

9 until  _there are no more linearly dependent points in the tuple_ ;

10 end for

Algorithm 4 Subspace Clustering based on the Hardt-Moitra Algorithm

The reason why we extract the smallest number of linearly dependent points at
each step is to avoid a situation where a $p$-tuple contains $d_{j}+1$ points
from $L_{j}$ and $d_{k}+1$ points from $L_{k}$ (with $j\neq k$), in which
case, assuming $d_{j}+d_{k}<p$, these points are linearly dependent but do not
span one of the subspaces. This particular step is, however, computationally
challenging as it amounts to finding the sparsest solution to a $p$-by-$p$
linear system, a problem known to be challenging (Tropp and Wright, 2010, Eq
1). One possibility is to replace this will finding the solution with minimum
$\ell_{1}$ norm (Tropp and Wright, 2010, Eq 8). The use of the $\ell_{1}$
constraint is central to the method proposed by Elhamifar and Vidal (2009).

###  3.3 The algorithm of Chen and Lerman

The Spectral Curvature Clustering (SCC) algorithm of Chen and Lerman (2009b)
is in fact of RANSAC type. The method was designed for the noisy setting and
is therefore more sophisticated.222 Chen and Lerman (2009b) consider the case
where the subspaces are affine, but we adapt their method to the case where
they area linear. It is based on a function
$\mathcal{A}:(\mathbb{R}^{p})^{d+1}\to[0,1]$ that quantifies how close a
$(d+1)$-tuple is from spanning a subspace of dimension $d$ or less. It is
equal to 1 when this is the case and is strictly less than 1 when this is not
the case. The algorithm draws a number, $c$, of $d$-tuples at random, where
the $s$-th tuple is denoted $(x_{1,s},\dots,x_{d,s})$, and computes the matrix
$\mathbf{W}=(W_{ij})$, where

$W_{ij}=\sum_{s=1}^{c}\mathcal{A}(x_{i},x_{1,s},\dots,x_{d,s})\mathcal{A}(x_{j},x_{1,s},\dots,x_{d,s}).$

It then applies a form of spectral graph partitioning algorithm to
$\mathbf{W}$ closely related to method of Ng et al. (2002). (The method
assumes all subspaces are of same dimension $d$, and both $d$ and $K$ are
assumed known.)

In the noiseless setting, one could take $\mathcal{A}$ to return 1 if the
tuple is linearly dependent and 0 otherwise. In that case, $W_{ij}$ is simply
the number of $d$-tuples among the $c$ that were drawn with whom both $x_{i}$
and $x_{j}$ are linearly dependent. Chen and Lerman (2009a) analyzes their
method in a setting that reduces to this situation and show that the method is
exact in this case.

###  3.4 Numerical experiments

We performed some numerical experiments to compare various methods for
subspace clustering, specifically, RANSAC, Sparse Subspace Clustering (SSC)
(Elhamifar and Vidal, 2009), Spectral Curvature Clustering (SCC) (Chen and
Lerman, 2009b), and Thresholding-based Subspace Clustering (TSC) (Reinhard
Heckel, 2015).

Each inlier is uniformly distributed on the intersection of the unit sphere
with its corresponding subspace. Each outlier is simply uniformly distributed
on the unit sphere. The result of each algorithm is averaged over 500 repeats.
Performance is measured by the Rand index. The results are reported in Table
2.

parameters | average system time |  | rand index
---|---|---|---
$(d,p,K,m,m_{0})$ | RANSAC | SSC | SCC | TSC |  | RANSAC | SSC | SCC | TSC
$(4,8,3,50,50)$ | .0622 | .3188 | .3731 | .1767 |  | 1 | .8407 | .9997 | .7237
$(6,8,3,50,50)$ | .8810 | .2813 | .8257 | .1815 |  | 1 | .6490 | .8322 | .5829
$(4,8,3,50,100)$ | .2108 | .3974 | .4563 | $.1757$ |  | 1 | .7619 | .9689 | $.2849$
$(4,8,5,50,50)$ | .3425 | .6163 | .9395 | .3282 |  | 1 | .9206 | .9548 | .7578
$(8,10,3,50,50)$ | 17.1827 | .3185 | 1.4141 | .2019 |  | 1 | .6148 | .6904 | .5688
Table 2: Numerical experiments comparing RANSAC, SSC, SCC, and TSC for the
problem of subspace clustering. As in the text, $d$ is the dimension of the
subspaces (assumed to be the same), $p$ is the ambient dimension, $K$ is the
number of subspaces, $m$ is the number of inliers per subspace (assumed to be
the same), $m_{0}$ is the number of outliers (so that $n=Km+m_{0}$ is the
sample size).

##  4 Discussion and conclusion

In our small scale experiments, RANSAC is seen to be competitive with other
methods, at least when the intrinsic dimensionality is not too large and when
there are not too many outliers (or too many underlying subspaces) present in
the data. This was observed both in the context of subspace recovery and in
the context of subspace clustering.

## References

  * Candès et al. (2011) Candès, E. J., X. Li, Y. Ma, and J. Wright (2011).  Robust principal component analysis?  Journal of the ACM (JACM) 58(3), 11.
  * Chen and Lerman (2009a) Chen, G. and G. Lerman (2009a).  Foundations of a multi-way spectral clustering framework for hybrid linear modeling.  Foundations of Computational Mathematics 9(5), 517–558.
  * Chen and Lerman (2009b) Chen, G. and G. Lerman (2009b).  Spectral curvature clustering (scc).  International Journal of Computer Vision 81(3), 317–330.
  * Elhamifar and Vidal (2009) Elhamifar, E. and R. Vidal (2009).  Sparse subspace clustering.  In Computer Vision and Pattern Recognition, 2009. CVPR 2009. IEEE Conference on, pp. 2790–2797. IEEE.
  * Fischler and Bolles (1981) Fischler, M. A. and R. C. Bolles (1981).  Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography.  Communications of the ACM 24(6), 381–395.
  * Hardt and Moitra (2013) Hardt, M. and A. Moitra (2013).  Algorithms and hardness for robust subspace recovery.  In Conference on Learning Theory (COLT), Volume 30, pp. 354–375.
  * Huber and Ronchetti (2009) Huber, P. J. and E. M. Ronchetti (2009).  Robust Statistics (2nd edition).  Wiley.
  * Maronna (1976) Maronna, R. A. (1976).  Robust M-estimators of multivariate location and scatter.  The Annals of Statistics 4(1), 51–67.
  * Ng et al. (2002) Ng, A. Y., M. I. Jordan, and Y. Weiss (2002).  On spectral clustering: Analysis and an algorithm.  In Advances in neural information processing systems, pp. 849–856.
  * Reinhard Heckel (2015) Reinhard Heckel, Helmut Bolcskei, D. o. I. . E. E. Z. S. (2015).  Robust subspace clustering via thresholding.  IEEE Transactions on Information Theory 61.
  * Schattschneider and Green (2012) Schattschneider, R. and R. Green (2012).  Enhanced ransac sampling based on non-repeating combinations.  In Proceedings of the 27th Conference on Image and Vision Computing New Zealand, pp. 198–203. ACM.
  * Tropp and Wright (2010) Tropp, J. A. and S. J. Wright (2010).  Computational methods for sparse solution of linear inverse problems.  Proceedings of the IEEE 98(6), 948–958.
  * Tyler (1987) Tyler, D. E. (1987).  A distribution-free m-estimator of multivariate scatter.  The Annals of Statistics, 234–251.
  * Vidal (2011) Vidal, R. (2011).  Subspace clustering.  IEEE Signal Processing Magazine 28(2), 52–68.
  * Vidal et al. (2005) Vidal, R., Y. Ma, and S. Sastry (2005).  Generalized principal component analysis (gpca).  IEEE Transactions on Pattern Analysis and Machine Intelligence 27(12), 1945–1959.
  * Wright et al. (2009) Wright, J., A. Ganesh, S. Rao, Y. Peng, and Y. Ma (2009).  Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization.  In Advances in neural information processing systems, pp. 2080–2088.
  * Zhang and Lerman (2014) Zhang, T. and G. Lerman (2014).  A novel M-estimator for robust PCA.  The Journal of Machine Learning Research 15(1), 749–808.


# Analysis of Black Hole Solutions in Parabolic Class Using Neural Networks

Ehsan<EMAIL_ADDRESS><EMAIL_ADDRESS>University of
Alcalá, Department of Signal Theory and Communications, Research group GRAM,
Alcalá de Henares, Spain.  Armin<EMAIL_ADDRESS>Department of
Mathematics and Statistics, Memorial University of Newfoundland, St John’s,
NL, Canada.  Roberto J<EMAIL_ADDRESS>University of
Alcalá, Department of Signal Theory and Communications, Research group GRAM,
Alcalá de Henares, Spain.

###### Abstract

In this paper, we introduce a numerical method based on Artificial Neural
Networks (ANNs) for the analysis of black hole solutions to the Einstein-
axion-dilaton system in a high dimensional parabolic class. Leveraging a
profile root-finding technique based on General Relativity we describe an ANN
solver to directly tackle the system of ordinary differential equations.
Through our extensive numerical analysis, we demonstrate, for the first time,
that there is no self-similar critical solution for the parabolic class in the
high dimensions of space-time. Specifically, we develop 95% ANN-based
confidence intervals for all the solutions in their domains. At the 95%
confidence level, our ANN estimators confirm that there is no black hole
solution in higher dimensions, hence the gravitational collapse does not
occur. Results provide some doubts about the universality of the Choptuik
phenomena. Therefore, we conclude that the fastest-growing mode of the
perturbations that determine the critical exponent does not exist for the
parabolic class in the high dimensions.

##  1 Introduction

It is well known that three parameters play key roles in describing black
holes. These parameters include the mass, the angular momentum and the charge
of the black holes. According to [1], there still may be another parameter
that may be relevant to black holes and their solutions. Christodoulou [2]
studied the existence of the spherically symmetric collapse solutions to the
Einstein scalar field. Choptuik [1] then proposed a distinguished property of
the collapse solutions illustrating the discrete self-similarity of the real
scalar field. Indeed, one can argue that various gravitational collapse
solutions represent space-time self-similarity; hence dilations happen.

Note that all the critical solutions illustrate the scaling law. The initial
condition of the scalar field can be characterized by the parameter $p$
relating to the field amplitude. Let $p=p_{\text{crit}}$ represent the
critical solution. The black hole is then formed when $p$ takes on a value
larger than $p_{\text{crit}}$. When $p$ becomes bigger than the critical
value, the following scaling law determines the mass of the black hole or the
Schwarzschild radius in $d\geq 4$ by

$r_{S}(p)\propto(p-p_{\text{crit}})^{\gamma}\,,\quad
M_{\text{bh}}(p)\sim(p-p_{\text{crit}})^{(D-3)\gamma}\,,$ (1)

where $\gamma\simeq 0.37$ for four dimensions and for a real scalar field [1,
3, 4, 5].

Several approaches investigated the numerical simulations for studying the
spherically symmetry of the scalar field for some matter content [6, 7, 8, 9,
10]. Moreover, in [11, 12, 4, 13] we find studies that evaluated the
gravitational collapse solutions of the perfect fluid. The critical exponent
$\gamma\simeq 0.36$ was first found in [12]. Strominger and Thorlacius [14]
discussed that $\gamma$ may have a universal value for all fields that can be
coupled to gravity in four dimensions.

Maison [13] has already shown the matter-dependence of the mass-scaling
exponent. Hatefi and Kuntz [15] recently showed that the critical exponent
depends on the space-time dimension and also they confirmed that the critical
exponent depends an the matter content and the ansätze of the various
solutions of self-similar collapse. These results thus provide doubts
concerning the universality of the Choptuik exponent. It should be mentioned
that one can derive the exponent based on various perturbations of self-
similar solutions [4, 13, 16]. In fact, one can start perturbing field $h$ as

$h=h_{0}+\varepsilon\,h_{-\kappa}\;,$ (2)

where $h_{-\kappa}$ includes the scaling $-\kappa\in\mathbb{C}$ related to
different modes. Note that the minus sign indicates a growing mode near the
black-hole formation time as $t\rightarrow 0$. When one requires to find the
most relevant mode $\kappa^{*}$ , corresponding to the highest value of
$\real(\kappa)$, $\gamma$ can be obtained from [4, 13, 16] by

$\gamma=\frac{1}{\real\kappa^{*}}\,.$ (3)

Abrahams and Evans [17] explored numerically the axisymmetric gravitational
collapse. Alvarez-Gaume et al. [18] found the other critical solutions as well
as the shock waves. Hirschmann and Eardley [19] derived the original results
for the axion-dilaton in four dimensions and concluded that $\gamma\simeq
0.2641$.

There is a relationship between the axion-dilaton system and the so-called
formation of black holes, particularly in analyzing their holographic
descriptions in diverse dimensions, see [5]. The axion-dilaton system is
relevant to the gauge/gravity correspondence [20], relating the critical
exponent and Choptuik exponent, the imaginary part of quasinormal modes and
the dual conformal field theory [21]. This axion-dilaton system has found
applications in various areas, including the physics of black holes [22, 23]
as well as the S-duality transformation in this particular class [24]. Another
application of the axion-dilaton system appears in cases where one is
interested in the collapse of spaces that can asymptotically tend to
$AdS_{5}\times S^{5}$. In these cases, the matter content is given by the
axion-dilaton system and the self-dual 5-form field.

Alvarez-Gaume and Hatefi [25] originally investigated the gravitational
collapse of the axion-dilaton system in the elliptic case. They computed the
critical collapse of the system by analyzing the hyperbolic ansatz [26].
Subsequently, the continuous self-similar solutions were calculated in [27] in
various dimensions. In a recent paper, Antonelli and Hatefi [28] derived the
perturbations and reproduced the correct value of $\gamma\sim 0.2641$ in four
dimensions [19]. Later on, in [15] a computation of the perturbation equations
and the extraction of other critical exponents was studied.

Recently, in [29, 30] Fourier-based regression models and nonlinear
statistical spline smothers are explored for estimating the gravitational
collapse functions of the axion-dilaton system under 4-dimensional elliptic
and hyperbolic spaces. While these statistical methods perform well enough in
estimating the nonlinear pattern of the critical collapse functions in four
dimensions case, the Fourier-based estimators may change significantly in the
case of slight changes in the estimates of the regression parameters. Another
challenge associated with the spline regression methods is that the methods
require realizations of the underlying population. In other words, one has to
solve the equations of motion iteratively to observe the response from the
critical collapse functions on their entire domain. This may be challenging
when finding a solution to the equations of motion results in numerical
instability, or there may not be a clear solution to distinguish the domains
of the critical collapse functions.

Previous works focus the analysis of the solutions for the axion-dilaton
system in low dimensional space-time. Our first contribution is that in this
work we extend the study of these solutions, under the parabolic class, to
higher dimensional space-time, ranging from four to nine. Given the recent
advances in Artificial Neural Networks (ANNs) for solving differential
equations [31, 32, 33], we propose, for the first time, an ANN-based solver to
estimate the self-similar gravitational collapse solutions of Einstein axion-
dilaton system under parabolic class. Leveraging a profile root-finding
technique based on General Relativity, our approach is able to directly tackle
the system of ordinary differential equations associated to the axion-dilaton
system. Unlike previous works in the literature on critical solutions in
parabolic class, we develop a stochastic approach to present the universal
estimator, modeled by the ANN, of the critical collapse functions to deal with
the sampling variability involved in neural networks for black hole solutions
in parabolic class. In other words, we construct 95% empirical intervals for
the estimators of all the critical collapse functions in parabolic case.
Through extensive numerical studies, at a 95% confidence level, we show that
there is no black hole solution in higher dimensions for the parabolic class
using our ANN-based solver.

This paper is organized as follows. Section 2 describes the self-similar
axion-dilaton system in $d$ dimensions and the equations of motion for the
third conjugacy class of SL(2,R) transformation. Section 3 explains our
approach by using ANNs to analyze the solutions to the the self-similar axion-
dilaton system. In Section 4, we investigate the performance of the ANN-based
method in estimating the critical collapse functions under parabolic class for
various space-time high dimensions. Section 5, finally, presents a summary and
concluding remarks.

##  2 Equations of Motion of Black Holes

In this section, we discuss the self-similar axion-dilaton system and then
develop the equations of motion for the third conjugacy class of SL(2,R)
transformation.

###  2.1 The Relevant System in Parabolic Class

The effective action of the Einstein axion-dilaton system that can be coupled
to $d$ dimensional gravity [34, 35] is given by

$S=\int
d^{d}x\sqrt{-g}\left(R-\frac{1}{2}\frac{\partial_{a}\tau\partial^{a}\bar{\tau}}{(\imaginary\tau)^{2}}\right)\;.$
(4)

The axion-dilaton can be combined by defining $\tau\equiv a+ie^{-\phi}$. The
SL(2,R) symmetry of action in Eq. (4) is shown by

$\tau\rightarrow M\tau\equiv\frac{a\tau+b}{c\tau+d},\quad\quad ad-bc=1\;,$ (5)

where $a,b,c$ and $d$ are real parameters.

When the quantum effects are taken into account, the SL(2,R) symmetry will be
related to $\mathrm{SL}(2,\mathbb{Z})$ and this duality also holds for the
case of the non-perturbative symmetry [36, 37, 38]. One can obtain the
equations of motion by taking the variations from the metric and the $\tau$
field, respectively. Hence the equations of motion are given by

$R_{ab}=\frac{1}{4(\imaginary\tau)^{2}}(\partial_{a}\tau\partial_{b}\bar{\tau}+\partial_{a}\bar{\tau}\partial_{b}\tau)\;,$
(6)
$\nabla^{a}\nabla_{a}\tau+\frac{i\nabla^{a}\tau\nabla_{a}\tau}{\imaginary\tau}=0\,.$
(7)

The metric in $d$ dimensions for the spherical symmetry can be written as

$ds^{2}=(1+u(t,r))(-b(t,r)^{2}dt^{2}+dr^{2})+r^{2}d\Omega_{d-2}^{2}\,,$ (8)

where $\tau$ field represents a function of both coordinates of $t$ and $r$,
that is $\tau=\tau(t,r)$, and $d\Omega_{d-2}^{2}$ denotes the angular part for
the metric.

If we use the scale-invariant $(t,r)\rightarrow(\Lambda t,\Lambda r)$, then
$ds^{2}$ must scale as $ds^{2}\rightarrow\Lambda^{2}ds^{2}$. Thus all the
functions in the metric should be scale-invariant, namely $u(t,r)=u(z)$,
$b(t,r)=b(z)$, $z\equiv-r/t$. The effective action in Eq. (4) is invariant
under the SL(2,R) transformation (5). This implies that $\tau$ should also be
invariant up to an SL(2,R) transformation as

$\tau(\Lambda t,\Lambda r)=M(\Lambda)\tau(t,r)\;.$ (9)

When the system of $(g,\tau)$ enjoys the above invariant properties, then the
system is called a continuous self-similar (CSS) solution. Note that the
different assumptions can be related to different classes of
$\evaluated{\derivative{M}{\Lambda}}_{\Lambda=1}$ [27], where $\tau(t,r)$ can
take the elliptic, hyperbolic and parabolic assumptions. From (9), the general
form of the ansatz for the elliptic, hyperbolic and parabolic class is given
by

$\tau(t,r)=\begin{dcases}i\frac{1-(-t)^{i\omega}f(z)}{1+(-t)^{i\omega}f(z)}\,,&\quad\text{elliptic}\\\\[5.0pt]
\frac{1-(-t)^{\omega}f(z)}{1+(-t)^{\omega}f(z)}\,,&\quad\text{hyperbolic}\\\\[5.0pt]
f(z)+\omega\log(-t)\,,&\quad\text{parabolic}\end{dcases}$ (10)

where $\omega$ is an unknown real parameter for all the classes. Also function
$f(z)$ satisfies the condition $\absolutevalue{f(z)}<1$ for the elliptic case;
however, function $f(z)$ satisfies the condition $\imaginary f(z)>0$ for
hyperbolic and parabolic classes.

It is worth mentioning that in Appendix A we show that the proposed ANNs-based
approach can confirm that there exists a solution in the four-dimensional
elliptic case which is compatible with findings of [27]. Here, we just study
the parabolic self-similar solutions in detail. From the equations of motion
(6) and (7), it can be shown that the following symmetry is present. In other
words, if

$f(z)\rightarrow f(z)+a\,,$ (11)

then all the equations of motion remain invariant.

Applying CSS ansätze (10) to equations of motion (6) and (7), one can develop
the ordinary differential equations for $u(z)$, $b(z)$, $f(z)$. If we take the
spherical symmetry, $u(z)$ can be eliminated in terms of $b(z)$ and $f(z)$,
employing

$u(z)=\frac{zb^{\prime}(z)}{(q-1)b(z)}\,.$ (12)

Note that the first derivative of $u(z)$ can also be removed from all
equations of motion as

$\frac{qu^{\prime}(z)}{2(1+u(z))}=\frac{w\bar{f}(z)f^{\prime}(z)+wf(z)\bar{f}^{\prime}(z)-2z\bar{f}^{\prime}(z)f^{\prime}(z)}{(f(z)-\bar{f}(z))^{2}}\,.$
(13)

Finally, all the ordinary differential equations (ODEs) are given by

$\displaystyle b^{\prime}(z)$ $\displaystyle=B(b(z),f(z),f^{\prime}(z))\,,$
(14) $\displaystyle f^{\prime\prime}(z)$
$\displaystyle=F(b(z),f(z),f^{\prime}(z)).$ (15)

From (6) and (7), the equations of motion for self-similar solutions for the
parabolic case in $d$ dimensions $(d=4,..,10)$ are given by the following
system

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
b^{\prime}-{2z(b^{2}-z^{2})\over(d-2)b(f-\bar{f})^{2}}f^{\prime}\bar{f}^{\prime}+{2\omega(b^{2}-z^{2})\over(d-2)b(f-\bar{f})^{2}}(\bar{f}^{\prime}+f^{\prime})+{2\omega^{2}z\over(d-2)b(f-\bar{f})^{2}},$
$\displaystyle 0$ $\displaystyle=$
$\displaystyle-f^{\prime\prime}-{2z(b^{2}+z^{2})\over(d-2)b^{2}(f-\bar{f})^{2}}f^{\prime
2}\bar{f}^{\prime}+{2\over(f-\bar{f})}\left(1+{\omega(b^{2}+z^{2})\over(d-2)b^{2}(f-\bar{f})}\right)f^{\prime
2}$ (16)
$\displaystyle+{2\omega(b^{2}+2z^{2})\over(d-2)b^{2}(f-\bar{f})^{2}}f^{\prime}\bar{f}^{\prime}+{2\over
z}\left(-\frac{(z^{2}-\frac{(d-2)b^{2}}{2})}{(z^{2}-b^{2})}+{2\omega
z^{2}\over(b^{2}-z^{2})(f-\bar{f})}\right.$
$\displaystyle+\left.{2\omega^{2}z^{4}\over(d-2)b^{2}(b^{2}-z^{2})(f-\bar{f})^{2}}\right)f^{\prime}-{2\omega^{2}z\over(d-2)b^{2}(f-\bar{f})^{2}}\bar{f}^{\prime}+$
$\displaystyle{2\omega\over(b^{2}-z^{2})}\left(-\frac{1}{2}-{\omega\over(f-\bar{f})}\right.-\left.{\omega^{2}z^{2}\over(d-2)b^{2}(f-\bar{f})^{2}}\right).$

These equations are invariant under some shifts of $f(z)$ by a real number.
Let $f(z)=u(z)+iv(z),\bar{f}(z)=u(z)-iv(z)$ where $u(z),v(z)$ are real
functions. Solving the system of equations (34) in Cartesian coordinates, the
equations of motion, for $d$ dimensions $(d=4,..,10)$ can be written as
follows:

$-2(d-2)b^{\prime}bv^{2}+w^{2}z+(z^{2}-b^{2})(-2wu^{\prime}+zu^{\prime
2}+zv^{\prime 2})=0\,.$ (17)

The other equation for $f^{\prime\prime}$ is given by

$\displaystyle
2(d-2)b^{2}v^{2}(u^{\prime\prime}+iv^{\prime\prime})z(z^{2}-b^{2})+w^{3}z^{3}-2(d-2)v^{2}wzb^{2}-3w^{2}z^{4}u^{\prime}+4(d-2)v^{2}z^{2}b^{2}u^{\prime}$
(18)
$\displaystyle+w^{2}z^{2}b^{2}u^{\prime}-2(d-2)^{2}v^{2}b^{4}u^{\prime}+3wz^{5}u^{\prime
2}-wz^{3}b^{2}u^{\prime 2}-2wzb^{4}u^{\prime 2}$ $\displaystyle-z^{6}u^{\prime
3}+z^{2}b^{4}u^{\prime
3}+4(d-2)vwz^{2}b^{2}v^{\prime}-4(d-2)vz^{3}b^{2}u^{\prime}v^{\prime}$
$\displaystyle+4(d-2)vzb^{4}u^{\prime}v^{\prime}+wz^{5}v^{\prime
2}-wz^{3}b^{2}v^{\prime 2}-z^{6}u^{\prime}v^{\prime
2}+z^{2}b^{4}u^{\prime}v^{\prime 2}$
$\displaystyle+i\bigg{(}2(d-2)vw^{2}zb^{2}-4(d-2)vwz^{2}b^{2}u^{\prime}+2(d-2)vz^{3}b^{2}u^{\prime
2}-2(d-2)vzb^{4}u^{\prime 2}-w^{2}z^{4}v^{\prime}$
$\displaystyle+4(d-2)v^{2}z^{2}b^{2}v^{\prime}-w^{2}z^{2}b^{2}v^{\prime}-2(d-2)^{2}v^{2}b^{4}v^{\prime}+2wz^{5}u^{\prime}v^{\prime}-2wzb^{4}u^{\prime}v^{\prime}$
$\displaystyle-z^{6}u^{\prime 2}v^{\prime}+z^{2}b^{4}u^{\prime
2}v^{\prime}-2(d-2)vz^{3}b^{2}v^{\prime 2}+2(d-2)vzb^{4}v^{\prime
2}-z^{6}v^{\prime 3}+z^{2}b^{4}v^{\prime 3}\bigg{)}=0\,.$

The regularity of $\tau$ and residual symmetries in the equations of motions
(34) entails the initial boundary conditions as

$\displaystyle\left\\{\begin{array}[]{lc}b(0)=1,&f^{\prime}(0)=0\,,\\\
f(0)=i\imaginary f(0),&\imaginary f(0)=x_{0}\,,\end{array}\right.$ (21)

where $x_{0}$ is a real parameter and $\imaginary f(0)>0$.

Hence, the problem is reduced to determine two real parameters, including
$\omega$ and $f(0)$. The equations of motion in (34) have five singularities
(see [25]) located at $z=\pm 0$, $z=\infty$ and $z=z_{\pm}$. The last two
singularities are expressed by $b(z_{\pm})=\pm z_{\pm}$. They are related to
the homothetic horizon where $z=z_{+}$ is just a mere coordinate singularity
[19, 25]. Thus, $\tau$ must be regular across it, representing that
$f^{\prime\prime}(z)$ must be finite as $z\rightarrow z_{+}$. Therefore, the
vanishing of the divergent part of $f^{\prime\prime}(z)$ generates a complex-
valued constraint at $z_{+}$, denoted by
$G(b(z_{+}),f(z_{+}),f^{\prime}(z_{+}))=0$. The explicit form of the $G$
function for the parabolic case is then given by

$\displaystyle G(f(z_{+}),f^{\prime}(z_{+}))=$
$\displaystyle\,-\frac{1}{f(z_{+})}\Bigg{(}-2(d-2)\omega\bar{f}(z_{+})\left(\omega-2z_{+}f^{\prime}(z_{+})\right)$
(22)
$\displaystyle+(d-2)\bar{f}(z_{+})^{2}\left((d-4)z_{+}f^{\prime}(z_{+})+\omega\right)$
$\displaystyle-2(d-2)f(z_{+})\left(\bar{f}(z_{+})\left((d-4)z_{+}f^{\prime}(z_{+})+\omega\right)-\omega\left(\omega-2z_{+}f^{\prime}(z_{+})\right)\right)$
$\displaystyle+(d-2)f(z_{+})^{2}\left((d-4)z_{+}f^{\prime}(z_{+})+\omega\right)+2\omega^{2}\left(\omega-2z_{+}f^{\prime}(z_{+})\right)\Bigg{)}\,.$

Owing to the vanishing of the real and imaginary parts of $G$, one can rewrite
the equations in (34) as a system of equations based on two unknown parameters
$(\omega,\imaginary f(0))$. Once the system is obtained, using $f(0)=ix_{0}$
and the boundary condition at $z=0$, one can find the discrete solutions by
integrating numerically the equations of motion (34). For instance, the
discrete solution in four dimensions for the elliptic case was first found in
[39] and then was confirmed in [25, 27].

###  2.2 Scaling Symmetry

One can show that there is an extra symmetry for the parabolic case from the
ansatz as

$\tau(t,r)=f(z)+\omega\log(-t)\,.$ (23)

If we scale $\omega$ and $f(z)$ by a constant $L$, then $\tau$ also scales by
$L$. In other words, if we have

$\omega\rightarrow L\omega\,,\quad f(z)\rightarrow Lf(z)\,,\quad
L\in\mathbb{R}_{+}\,,$ (24)

then $\tau$ also transforms as $\tau\rightarrow L\tau$. Hence, if
$(\omega,\imaginary f(0))$ is a solution, then $(L\omega,L\imaginary f(0))$
leads to the same solution.

Note that $G$ equation (22) scales homogeneously as $L^{2}$ under
$(\omega,f)\rightarrow(L\omega,Lf)$, because all equations of motion and the
constraint $G(\omega,\imaginary f(0))$ are invariant under this new scaling.
According to this degeneracy, the only real unknown parameter becomes the
ratio $\omega/\imaginary f(0)$, and the two-dimensional constraint should be
solved for it.

Due to the fact that $G(\omega,Imf(0))$ only depends on $\omega/\imaginary
f(0)$, one can apply a profile root-finding method to determine the solutions
to the equation of $G$. To do that, we first set $\imaginary f(0)=1$ and
treated equation (22) as a function of $\omega$, and next found the root of
$G$ on $\omega$ coordinate. If there exists a root, then it generates a
continuous ray of solutions $(\omega,\imaginary f(0))=(K\omega^{*},K)$.
Different solutions in the same ray must be related to an SL(2,R)
transformation where they are not physically distinguishable. When
$\tau(t,r)=\tau(z)$, we can show that this trivial solution can be obtained
from any of the previous solutions by $\omega\rightarrow 0$, using the initial
boundary conditions $b(0)=1$, $\tau^{\prime}(0)=0$. Therefore,
$\tau(0)=\tau_{0}$ should be an arbitrary constant. These are, in fact, the
initial conditions that reproduce flat space-time with a constant axion-
dilaton of $\tau_{0}$.

Figure 1: The contours of the real and imaginary parts of $G(\omega,\imaginary
f(0))$ in four-dimensional (top panels) and five-dimensional (bottom panels)
parabolic class. The blue and orange curves represent the real and imaginary
functions, respectively.

The top panel of Figure 1 shows two-dimensional contours of the zeroes of the
real part (with blue colour) and imaginary part (with orange colour) of
$G(\omega,\imaginary f(0))$, as functions of $(\omega,\log\imaginary f(0))$ in
four-dimensional parabolic space. We observe that the curves corresponding to
the real and imaginary parts of the $G$ function do not intersect in the
entire region. The top panel of Figure 1 demonstrates the degeneracy
associated with the extra scaling invariance (24) in four dimensions, as
described earlier. According to the wide region of $(\omega,\log\imaginary
f(0))$ in the top-left panel of Figure 1, one may argue that there may exist a
solution for the small values of $\omega$ and $\log\imaginary f(0)$. To
elaborate on this argument, we show the contours of the real and imaginary
parts of $G(\omega,\imaginary f(0))$ but now on smaller values of $\omega$ and
$\log\imaginary f(0))$ in the top-right panel of Figure 1. From the later
plot, there is clearly no intersection between the real and imaginary parts of
the $G$ function. This verifies that no solution exists to equations in (34)
in the four-dimensional parabolic space for $\omega>0$.

###  2.3 Solutions in Higher Dimensions

In this subsection, we investigate the solutions to equations of motion in
(34) for parabolic space in various space-time dimensions ranging from five to
nine. Similarly to the four dimensions, one can use the constraints on
vanishing the real and imaginary parts of $G(\omega,\imaginary f(0))$ and re-
write equation (34) as a function of two unknown real parameters $\omega$ and
$\imaginary f(0)$. Using an appropriate starting point for $\imaginary f(0)$,
one can apply the profiling method and find the root of $G$ function (22).

Figure 1 (bottom panel) and Figure 2 demonstrate the contours using the zeroes
of the real and imaginary parts of $G(\omega,\imaginary f(0))$ function (22)
for different dimensions ranging from five to nine. In each plot, blue and
orange curves indicate the real and imaginary parts of the
$G(\omega,\imaginary f(0))$ function, respectively. From Figures 1 and 2
demonstrate that the two parts do not intersect throughout the entire domain
of the parameters $(\omega,\log\imaginary f(0))$. This concludes that no
parabolic black hole solution exists even in higher space-time dimensions
$d\in\\{5,\ldots,9\\}$. This leads to the fact that gravitational collapse
will not take place in high-dimensional parabolic cases and it leads to the
Minkowski flat space-time.

Figure 2: The contours of the real and imaginary parts of $G(\omega,\imaginary
f(0))$ in parabolic case for various dimensions including $d=6$ (top left),
$d=7$ (top right), $d=8$ (bottom left) and $d=9$ (bottom right).

Note that, for the five-dimensional parabolic case, one can search for the
zeroes of the real and imaginary parts of $G(\omega,1)$ and may find a
spurious solution ray. This solution ray may stem from the numerical
instability of the equations of motion for a very small value of $\omega$. For
further remarks in this case, readers are referred to [40]. While the
numerical accuracy is insufficient to assess the solution with certainty in
five dimensions, using the root-finding method, one may identify a very small
value for $|G|$ and a value of $\omega$ as

$\displaystyle\absolutevalue{G}\sim 0.006,\quad\omega\sim 1.65\;.$ (25)

It is evident from Figure 2 that the higher the space-time grows, the more
divergence is observed between the real and imaginary parts of the
$G(\omega,\imaginary f(0))$ function. In other words, in the higher
dimensions, the real part of the $G$ function diverges so quickly and
disappears in the entire domain, that it is more likely that there is no
interaction between the real and imaginary curves. Therefore, it confirms that
no parabolic black hole solution is expected in higher dimensions.

##  3 Search for Solutions Using Neural Networks

In this section, we describe an ANN-based solver for the equations of motion
in (17) and (18). We intend to provide an alternative approach to the method
described, in order to shed more light on the problem, confirming that no
black hole solution exists in higher dimensional parabolic cases.

In the simplest form, an artificial neural network consists in a multi-layer
perceptron [41]. The neural network utilizes linear and nonlinear
transformations to dispatch data from one layer to another, ultimately from
the input layer to the output layer. Many variations of neural networks have
been developed, such as fully connected neural networks, convolutional neural
networks, and recursive neural networks, to name a few [41, 42]. Artificial
neural networks have found applications in almost all scientific fields
ranging from bioinformatics [43] and material sciences [44], to natural
language processing [45]. In this work, we focus on the fully connected neural
networks which are sufficient to deal with the non-linearity of the Einstein-
axion-dilaton system in the high-dimensional parabolic classes.

Let $\mathcal{N}^{L}({\bf x},t,{\bf\phi})$ represent a deep neural network
with $L$ layers mapping
$\mathbb{R}^{d_{\text{in}}}\rightarrow\mathbb{R}^{d_{\text{out}}}$.
$d_{\text{in}}$ and $d_{\text{out}}$ encode the dimensionality of the input
and output of the network. Let $n_{l}$ denote the number of neurons for the
$l$-th layer of the network, for $l=1,\ldots,L-1$. $n_{0}=d_{\text{in}}$ and
$n_{L}=d_{\text{out}}$ are the number of neurons in the input and output
layers, respectively. Suppose ${\bf W}^{l}$ for $l=1,\ldots,L-1$ denotes the
weight matrix of size $(n_{l}\times n_{l-1})$ associated with the $l$-th layer
of the neural network, whose $(j,k)$-th entry, namely ${W}^{l}_{j,k}$,
represents the weight connecting the $j$-th neuron in layer $l$ with the
$k$-th neuron in layer $l-1$. Let ${\bf b}^{l}$ denote the bias vector in the
$l$-th layer where ${b}_{j}^{l}$ gives the bias associated with the $j$-th
neuron in the $l$-th layer. Let ${\bf Z}^{l}$ be the response observed from
the $l$-th layer of the network before the activation. The response variable
${\bf Z}^{l}$ is then given by

$\displaystyle{\bf Z}^{l}={\bf W}^{l}~{}\text{NN}^{l-1}({\bf
x},t,{\bf\phi})+{\bf b}^{l},~{}~{}~{}l=1,\ldots,L-1.$ (26)

Suppose $\sigma$ is a nonlinear activation function applied element-wise to
the response variable in the layers. The input layer is $\text{NN}^{0}({\bf
x},t,{\bf\phi})={\bf x}$. From (26), the output from the $l$-th layer of the
neural network is given by

$\displaystyle\text{NN}^{l}({\bf x},t,{\bf\phi})=\sigma({\bf
Z}^{l}),~{}~{}~{}l=1,\ldots,L-1$ (27)

and the output in the last layer is as follows

$\displaystyle\text{NN}^{L}({\bf x},t,{\bf\phi})={\bf Z}^{L}\;.$ (28)

Overall, the neural network architecture is characterized by its parameters
(or weights) which can be encoded as follows ${\bf\phi}=({\bf
W}^{1},\ldots,{\bf W}^{L},{\bf b}^{1},\ldots,{\bf b}^{L})$.

To train a neural network, one requires the definition of a loss function,
which measures the performance of the trained machine. Let
$\mathcal{L}\left(\mathcal{N}^{L}({\bf x},t,{\bf\phi})\right)$ denote the loss
function of the neural network $\mathcal{N}^{L}({\bf x},t,{\bf\phi})$. In
addition to the above forward propagation steps (see equations (27) and (28)),
one needs to consider the back-propagation to make sure that the responses in
the network make sense. The back-propagation algorithm evaluates the gradients
of the loss function with respect to the response variables in the previous
layers. It is easy to show that the back-propagation boils down to

$\displaystyle\frac{\partial\mathcal{L}\left(\mathcal{N}^{L}({\bf
x},t,{\bf\phi})\right)}{\partial{\bf Z}^{l}}=\left[\left({\bf
W}^{l+1}\right)^{\top}\left(\frac{\partial\mathcal{L}\left(\mathcal{N}^{L}({\bf
x},t,{\bf\phi})\right)}{\partial{\bf
Z}^{l+1}}\right)\right]\odot\frac{\partial\sigma({\bf Z}^{l})}{\partial{\bf
Z}^{l}}\;,$ (29)

where $\odot$ stands for the Hadamard product calculating the element-wise
product between two matrices. Finally, one can compute the back-propagation of
the last layer by

$\displaystyle\frac{\partial\mathcal{L}\left(\mathcal{N}^{L}({\bf
x},t,{\bf\phi})\right)}{\partial{\bf
Z}^{L}}=\frac{\partial}{\partial\text{NN}^{L}({\bf
x},t,{\bf\phi})}\mathcal{L}\left(\mathcal{N}^{L}({\bf
x},t,{\bf\phi})\right)\odot\frac{\partial\sigma({\bf Z}^{L})}{\partial{\bf
Z}^{L}}.$ (30)

For more details about the deep neural network algorithms, readers are
referred to [41, 42].

The Universal Approximation Theorem declares that neural networks are
universal approximators: no matter the nature of a function, there is a neural
network that can approximately approach its result. Because of this property
of neural networks, various researchers have recently investigated ways to
employ neural networks to solve ordinary and partial differential equations,
e.g. [31].

The basic concept in solving differential equations using an ANN is to re-
frame the problem as an optimization problem in which the objective consists
of minimizing the squared residual of the differential equations. To do so,
one can generalize any differential equation as

$\mathcal{D}u({\bf x},t)-f=0\;,$

where $\mathcal{D}$ is the differential operator, $u({\bf x},t)$ is the
solution of interest, and $f$ is a known forcing function.

As described earlier, ANNs can be considered universal approximates for any
function. We consider the output of a neural network as $\mathcal{N}^{L}({\bf
x},t,{\bf\phi})$. $\mathcal{N}^{L}({\bf x},t,{\bf\phi})$ provides us with the
approximation for our solution for the differential equation. Given
$\mathcal{N}^{L}({\bf x},t,{\bf\phi})$ as a solution to the differential
equation, we can construct the loss function by the squared residuals in the
form of

$\mathcal{L}(\mathcal{N}^{L}({\bf x},t,{\bf\phi}))=\left(\mathcal{N}^{L}({\bf
x},t,{\bf\phi})-f\right)^{2}\;.$ (31)

The loss function (31) now can be incorporated into the forward and backward
propagation steps (26)-(30) to measure the loss corresponding to the set of
learned parameters.

It is easily seen that equation (31) translates the problem of finding
solutions to a system of differential equations, to an optimization process of
the neural network loss function by

$\arg\min_{{\bf\phi}}\mathcal{L}(\mathcal{N}^{L}({\bf x},t,{\bf\phi})).$ (32)

Equation (32) indicates that the neural network is designed to minimize the
squared residuals of the trained $\mathcal{N}^{L}({\bf x},t,{\bf\phi})$.
Consequently, this guarantees that $\mathcal{N}^{L}({\bf x},t,{\bf\phi})$
represents the functional form of the solution of the differential equations
when the neural network learning process converges.

We follow this approach in this work to tackle the problem of solving the
system of ordinary differential equations defined by the equations of motion
of black holes in the parabolic class, as we described in Section 2 (Equations
(17) and (18)). This system of equations can be approached with a neural
network-based solver like the one described. Our goal here is to provide an
alternative approach to General Relativity-based profile root-fining,
described in Section 2, so that we can explore whether there are self-similar
solutions for high dimensions in the parabolic class.

##  4 Numerical Studies Based on ANNs

In this section, we explore the properties of the fully connected neural
networks, described in Section 3, in estimating the critical collapse
functions corresponding to equations of motion in higher dimensional parabolic
class. Using the neural network estimates for the critical collapse functions,
we show that no self-similar critical solution exists for the parabolic class
in higher dimensions of space-time.

Due to the invariance of $f(z)-{\bar{f}}(z)$, it is clear from the equations
of motion that all the equations are invariant under a shift of $f(z)$. In
other words, the entire equations do not depend on $u(z)$, implying that
$u(0)=0$. On the other hand, the regularity of $\tau$ (10) at $z=0$ implies
$f^{\prime}(0)=u^{\prime}(0)+iv^{\prime}(0)=0$. Hence, all boundary conditions
are given by

$b(0)=1,~{}u^{\prime}(0)=v^{\prime}(0)=u(0)=0.$ (33)

Similar to Section 2, from the regularity of $\tau$ and residual symmetries,
the boundary conditions of the required critical functions are derived by
(21). In a similar vein, as the behaviours of the real and imaginary parts of
$f$ do not depend on the value of $\imaginary f(0)$, we can set the boundary
conditions of the two real parameters as $\imaginary f(0)=1$ and $\omega=1.65$
in the ANNs estimating method for five dimensional space-time. From (25) we
also set $\omega=1.65$ and $\imaginary f(0)=0.321$ for all the other
dimensions, see [40, 26].

Figure 3: The ANN architecture implemented to optimize the parabolic equations
of motion. It is a multi-layer perceptron of four hidden layers with 16 units
each.

To implement the ANN-based method, described in Section 3, we used the python
package NeuroDiffEq [31]. We have implemented a specific solver to obtain the
estimates of the critical collapse functions $b_{0}(z)$, $Re(f(z))$ and
$Im(f(z))$ corresponding to the parabolic class for $d=4,\ldots,9$ space-time
dimensions. We integrate in this solver the NN network architecture shown in
Figure 3. It consists of a multilayer perceptron for each function, with an
input layer, a set of four hidden fully connected layers with 16 units each,
and the output layer. Every hidden unit is followed by a rectified linear unit
(ReLU) activation function [46, 41]. This architecture provides the necessary
stability to find the solutions of the system of differential equations we are
dealing with.

According to the profile root-finding studies discussed in Section 2, if there
was a parabolic solution to the equations of motion, the solution would appear
in a neighbourhood of small values of space-time, namely $z<3$. For this
reason, using the above ANN configuration, we validated the estimates of the
critical collapse functions at 1000 equally spaced points in the interval
$z_{j}\in(0,3],~{}j=1,\ldots,1000$. We then computed the ANN estimates using
3000 epochs to reach a good level of convergence of the loss functions in
learning all the critical functions for $d=4,\ldots,9$ dimensions.

Due to the non-linear patterns of the parabolic equations of motion, and the
trained ANN solvers’ sampling variability, the critical collapse function
estimates slightly change from one trial to another. To cope with this
stochastic property of the proposed ANN estimates, we independently replicated
20 times the above ANN estimation procedure. That is, the ANN solver was
independently applied to the differential equations for 20 times. We then
validated the performance of the trained ANN models at points
$z_{j},~{}j=1,\ldots,1000$. Skewness in the performance of the ANN estimations
is unavoidable because of the non-linear pattern of the critical functions in
their domain and the stochastic property of the learned parameters. Therefore,
the critical collapse functions in the domain are considered our statistical
parameters of interest. We then constructed the 95% empirical confidence
intervals for all the critical collapse functions. To construct the confidence
intervals, we obtained the 2.5, 50 and 97.5 percentiles of the ANN-based
estimates at each space-time point $z_{j},j=1,\ldots,1000$. We finally
replicated the above ANN-based method to estimate the parabolic critical
collapse functions in all dimensions $d=4,\ldots,9$.

Figure 4: The 95% ANN-based confidence intervals for the critical collapse
functions in the four-dimensional parabolic class. The ANN lower (dotted
lines), upper (dashed lines) and median (solid lines) bands were computed from
20 replicates with 3000 epochs each.  Figure 5: The 95% ANN-based confidence
intervals for the critical collapse functions in the five-dimensional
parabolic class. The ANN lower (dotted lines), upper (dashed lines) and median
(solid lines) bands were computed from 20 replicates with 3000 epochs each.
Figure 6: The differences between train and validation loss functions of the
developed ANNs in estimating the critical collapse functions. The top, middle
and bottom panels show the loss differences of the estimates in 4D, 5D and 6D
parabolic spaces, respectively. The panels show the line plot and box plot of
loss differences of the ANN estimates for 20 replicates with 3000 epochs each.

Figures 4 \- 10 show the 95% ANN-based confidence intervals of the critical
collapse functions $b_{0}(z)$ (green colour), $Re(f(z))$ (red colour) and
$Im(f(z))$ (blue colour) corresponding to parabolic class for $d=4,\ldots,9$
space-time dimensions. The lower, upper and middle bands of the confidence
intervals are represented by 2.5 (dotted line), 97.5 (dashed line) percentiles
and median (solid line), respectively. We observe from Figures 4 \- 10 that no
intersection occurs between 95% ANN-based confidence regions of functions
$Re(f(z))$ and $Im(f(z))$ in parabolic class of $d=4,\ldots,9$ dimensions. As
the dimension grows, in addition, more divergence is observed between the ANN-
based confidence regions of the real and imaginary parts of
$G(\omega,\imaginary f(0))$. Therefore, Figures 4 \- 10 confirm that, at 95%
level of confidence, there exists no self-similar black hole solution for
parabolic class in $d=4,\ldots,9$ dimensions. This ANN-based estimation method
confirms the numerical results of the profile root-finding of Section 2, but
from a different perspective, directly tackling the system of differential
equations defined by the equations of motion of black holes.

We assess the convergence of the ANN estimators for the critical collapse
functions in Figures 6 and 11. With respect to the effect of the learned
parameters (for the ANN) on the behaviour of the loss functions, we analyzed
the loss differences observed in the 20 replicates of the ANN estimation
method. To do so, we computed the difference between training and test loss
functions in the last epoch of the ANN learning process. Figures 6 and 11 show
the line and box plots of this loss difference of the proposed ANN method in
the estimation of the critical collapse functions. We clearly observe that the
median difference between training and test losses is close to zero. Despite
the sampling variability, this confirms, on average, the convergence of the
ANN-based method in solving the equations of motion in $d=4,\ldots,9$
dimensions for parabolic class.

Figure 7: The 95% ANN-based confidence intervals for the critical collapse
functions in the six-dimensional parabolic class. The ANN lower (dotted
lines), upper (dashed lines) and median (solid lines) bands were computed from
20 replicates with 3000 epochs each.

In Figure 5 there is a short intersection between the 95% ANN-based confidence
intervals of the real and imaginary parts of $G(\omega,\imaginary f(0))$ in
five dimensions. This spurious black hole solution is compatible with the
finding of Hatefi and Vanzan [40] where equations of motion in the five-
dimensional parabolic class may sometimes lead to a solution due to the
numerical instability of the solvers.

To elaborate on this argument, we investigated another numerical study where
more estimates are considered. Here, we followed our ANN method, as described
earlier, in estimating the critical functions in five dimensions with 3000
epochs. We replicated the ANN estimates 100 times instead of just 20 times.
Figures 12 and 13 show the results of this extensive numerical study. From
Figure 12, we observe that there is no intersection between 95% ANN-confidence
intervals of the real and imaginary parts of $G(\omega,\imaginary f(0))$ in 5
dimensions. Therefore, at a 95% confidence level, no parabolic black hole
solution exists even in 5 dimensions. We believe that the spurious solution is
due to the complexity of the equations of motion and the numerical instability
of the neural network solvers. The ANN solver sometimes converges to the local
optimum rather than the global optimum. Therefore the variance of the ANN
statistic increases. This numerical instability consequently increases the
length of the ANN-confidence intervals of the critical collapse functions.
Figure 13 shows clearly this numerical instability such that the ANN solver
converges to a local optimum of the equations of motion for 3-4 replicates
(out of 100 replicates) and results in a significant loss difference.

##  5 Summary and Concluding Remarks

In this paper, we proposed a numerical approach to find the self-similar
solutions to the axion-dilaton system in higher dimensions of the parabolic
class. Specifically, we have described an ANN-based approach that allows us to
directly solve the differential equations proposed in the axion-dilaton model.
Unlike the statistical models in [29, 30], the ANN-based estimators take
directly the equations of motion as input and translate the estimation problem
into an optimization process. The proposed numerical method is quite generic
and can be applied to all SL(2, R) transformation classes in any dimension and
any matter content. In fact, we have shown in Appendix A that the proposed
ANN-based approach could confirm that there exists a solution which is
compatible with the findings in the literature for the four-dimensional
elliptic case. We found that the proposed ANNs estimation approach detected a
solution in the equation of motions with probability 1 in the numerical
experiment. This empirically confirms a solution in the domain of the critical
collapse functions in the elliptic class of four dimensions.

Unlike elliptic and hyperbolic classes, for the parabolic class in higher
dimensional space-time we have shown that there exists no self-similar
critical solution due to an extra scaling symmetry. According to the sampling
variability of neural networks, we constructed 95% empirical confidence
intervals based on our fully connected ANNs for all the critical collapse
functions in parabolic class for various dimensions ranging from four to nine.
Through extensive numerical studies, at a 95% confidence level, the ANN-based
solvers confirm that there is no black hole solution in higher dimensions for
the parabolic class.

The gravitational collapse does not thus take place for parabolic ansatz. Our
results provide some doubts about the universality of the Choptuik phenomena.
Therefore, the fastest-growing mode of the perturbations determining the
Choptuik exponent does not exist in these cases, and no gravitational collapse
occurs. It is also worth mentioning that our analysis provides reasons not to
expect to transfer the standard derivations of statistical mechanics to the
critical gravitational collapse.

## Acknowledgments

The authors would like to thank the editor and the anonymous referee for their
constructive comments which improved the quality and presentation of the
manuscript. Ehsan Hatefi would like to thank A. Kuntz, E. Hirschmann, C.
Gutiérrez-Álvarez and L. Álvarez-Gaumé for fruitful conversations. The
research of Ehsan Hatefi is supported by the María Zambrano Grant of the
Ministry of Universities of Spain. Armin Hatefi acknowledges the support from
the Natural Sciences and Engineering Research Council of Canada (NSERC).

## References

  * [1] M.W. Choptuik, Universality and Scaling in Gravitational Collapse of a Massless Scalar Field, Phys. Rev. Lett. 70, 9 (1993).
  * [2] D. Christodoulou, “The Problem of a Self-gravitating Scalar Field,” Commun. Math. Phys. 105 (1986) 337; “Global Existence of Generalized Solutions of the Spherically Symmetric Einstein Scalar Equations in the Large,” Commun. Math. Phys. 106 (1986) 587;“The Structure and Uniqueness of Generalized Solutions of the Spherically Symmetric Einstein Scalar Equations,” Commun. Math. Phys. 109 (1987) 591.
  * [3] R. S. Hamade and J. M. Stewart, “The Spherically symmetric collapse of a massless scalar field,” Class. Quant. Grav. 13 (1996) 497 [arXiv:gr-qc/9506044].
  * [4] T. Koike, T. Hara, and S. Adachi, “Critical Behavior in Gravitational Collapse of Radiation Fluid: A Renormalization Group (Linear Perturbation) Analysis,” Phys. Rev. Lett. 74 (1995) 5170 [gr-qc/9503007].
  * [5] L. Alvarez-Gaume, C. Gomez and M. A. Vazquez-Mozo, “Scaling Phenomena in Gravity from QCD,” Phys. Lett. B 649 (2007) 478 [hep-th/0611312].
  * [6] M. Birukou, V. Husain, G. Kunstatter, E. Vaz and M. Olivier, “Scalar field collapse in any dimension,” Phys. Rev. D 65 (2002) 104036 [gr-qc/0201026].
  * [7] E. Sorkin and Y. Oren, “On Choptuik’s scaling in higher dimensions,” Phys. Rev. D 71, 124005 (2005) [arXiv:hep-th/0502034].
  * [8] J. Bland, B. Preston, M. Becker, G. Kunstatter and V. Husain, “Dimension-dependence of the critical exponent in spherically symmetric gravitational collapse,” Class. Quant. Grav. 22 (2005) 5355 [gr-qc/0507088].
  * [9] E. W. Hirschmann and D. M. Eardley, “Universal scaling and echoing in gravitational collapse of a complex scalar field,” Phys. Rev. D 51 (1995) 4198 [gr-qc/9412066].
  * [10] J. V. Rocha and M. Tomašević, “Self-similarity in Einstein-Maxwell-dilaton theories and critical collapse,” Phys. Rev. D 98 (2018) no.10, 104063 [arXiv:1810.04907 [gr-qc]].
  * [11] L. Alvarez-Gaume, C. Gomez, A. Sabio Vera, A. Tavanfar and M. A. Vazquez-Mozo,“Critical gravitational collapse: towards a holographic understanding of the Regge region,” Nucl. Phys. B 806 (2009) 327 [arXiv:0804.1464 [hep-th]].
  * [12] C. R. Evans and J. S. Coleman, “Observation of critical phenomena and selfsimilarity in the gravitational collapse of radiation fluid,” Phys. Rev. Lett. 72 (1994) 1782 [gr-qc/9402041].
  * [13] D. Maison, “Non-Universality of Critical Behaviour in Spherically Symmetric Gravitational Collapse,” Phys. Lett. B 366 (1996) 82 [gr-qc/9504008].
  * [14] A. Strominger and L. Thorlacius, “Universality and scaling at the onset of quantum black hole formation,” Phys. Rev. Lett. 72 (1994) 1584 [hep-th/9312017].
  * [15] E. Hatefi and A. Kuntz, “On Perturbation Theory and Critical Exponents for Self-Similar Systems,” Eur. Phys. J. C 81 (2021) no.1, 15 [arXiv:2010.11603 [hep-th]].
  * [16] E. W. Hirschmann and D. M. Eardley, ‘Critical exponents and stability at the black hole threshold for a complex scalar field,” Phys. Rev. D 52 (1995) 5850 [gr-qc/9506078].
  * [17] A. M. Abrahams and C. R. Evans, “Critical behavior and scaling in vacuum axisymmetric gravitational collapse,” Phys. Rev. Lett. 70 (1993) 2980.
  * [18] L. Alvarez-Gaume, C. Gomez, A. Sabio Vera, A. Tavanfar and M. A. Vazquez-Mozo, “Critical formation of trapped surfaces in the collision of gravitational shock waves,” JHEP 0902, 009 (2009) [arXiv:0811.3969 [hep-th]].
  * [19] E. W. Hirschmann and D. M. Eardley, “Criticality and bifurcation in the gravitational collapse of a selfcoupled scalar field,” Phys. Rev. D 56 (1997) 4696 [gr-qc/9511052].
  * [20] J. M. Maldacena, “The Large N limit of superconformal field theories and supergravity,”Int. J. Theor. Phys. 38 (1999), 1113-1133, Adv. Theor. Math. Phys. 2, arXiv:hep-th/9711200, E. Witten, “Anti-de Sitter space and holography,” Adv. Theor. Math. Phys. 2 (1998), 253-291,hep-th/9802150, S. Gubser, I. R. Klebanov and A. M. Polyakov, “Gauge theory correlators from noncritical string theory,” Phys. Lett. B 428 (1998), 105-114, hep-th/9802109.
  * [21] D. Birmingham, “Choptuik scaling and quasinormal modes in the AdS / CFT correspondence,” Phys. Rev. D 64 (2001), 064024 [arXiv:hep-th/0101194 [hep-th]].
  * [22] E. Hatefi, A. Nurmagambetov and I. Park, “ADM reduction of IIB on $\mathcal{H}^{p,q}$ to dS braneworld,” JHEP 04 (2013), 170, arXiv:1210.3825 , “$N^{3}$ entropy of $M5$ branes from dielectric effect,” Nucl. Phys. B 866 (2013), 58-71, arXiv:1204.2711, S. de Alwis, R. Gupta, E. Hatefi and F. Quevedo, “Stability, Tunneling and Flux Changing de Sitter Transitions in the Large Volume String Scenario,” JHEP 11 (2013), 179, arXiv:1308.1222,
  * [23] A. Ghodsi and E. Hatefi, “Extremal rotating solutions in Horava Gravity,” Phys. Rev. D 81 (2010) 044016 [arXiv:0906.1237 [hep-th]].
  * [24] R. S. Hamade, J. H. Horne and J. M. Stewart, “Continuous Self-Similarity and $S$-Duality,” Class. Quant. Grav. 13 (1996) 2241 [arXiv:gr-qc/9511024].
  * [25] L. Álvarez-Gaumé and E. Hatefi, “Critical Collapse in the Axion-Dilaton System in Diverse Dimensions,” Class. Quant. Grav. 29 (2012) 025006 [arXiv:1108.0078 [gr-qc]].
  * [26] L. Álvarez-Gaumé and E. Hatefi, “More On Critical Collapse of Axion-Dilaton System in Dimension Four,” JCAP 1310 (2013) 037 [arXiv:1307.1378 [gr-qc]].
  * [27] R. Antonelli and E. Hatefi, “On self-similar axion-dilaton configurations,” , JHEP 03 (2020), 074 [arXiv:1912.00078 [hep-th]].
  * [28] R. Antonelli and E. Hatefi, “On Critical Exponents for Self-Similar Collapse,” JHEP 03 (2020), 180 [arXiv:1912.06103 [hep-th]].
  * [29] E. Hatefi and A. Hatefi, “Nonlinear statistical spline smoothers for critical spherical black hole solutions in 4-dimension,” Annals Phys. 446 (2022), 169112 [arXiv:2201.00949 [gr-qc]].
  * [30] E. Hatefi and A. Hatefi, “Estimation of Critical Collapse Solutions to Black Holes with Nonlinear Statistical Models,” Mathematics 10 (2022) no.23, 4537 [arXiv:2110.07153 [gr-qc]].
  * [31] F. Chen, D. Sondak, P. Protopapas, P.Mattheakis, M. Liu, S.Agarwal, D. & Di Giovanni, M., “NeuroDiffEq: A Python package for solving differential equations with neural networks,” Journal Of Open Source Software. 5, 1931 (2020)
  * [32] Lu. Lu, Xuhui. Meng, Zhiping. Mao, and George Em. Karniadakis, “DeepXDE: A Deep Learning Library for Solving Differential Equations,” SIAM Review, 63, 208-228 (2021)
  * [33] Alexander Koryagin, Roman Khudorozkov, and Sergey Tsimfer “PyDEns framework for solving differential equations with deep learning,” 2019
  * [34] A. Sen, “Strong - weak coupling duality in four-dimensional string theory,” Int. J. Mod. Phys. A 9 (1994) 3707 [hep-th/9402002].
  * [35] J. H. Schwarz, “Evidence for nonperturbative string symmetries,” Lett. Math. Phys. 34 (1995) 309 [hep-th/9411178].
  * [36] M.B. Green, J.H. Schwarz and E. Witten, 1987 Superstring Theory Vols I,II, Cambridge University Press,
  * [37] J. Polchinski, 1998 String Theory, Vols I,II, Cambridge University Press
  * [38] A. Font, L. E. Ibanez, D. Lust and F. Quevedo, “Strong - weak coupling duality and nonperturbative effects in string theory,” Phys. Lett. B 249 (1990) 35.
  * [39] D. M. Eardley, E. W. Hirschmann and J. H. Horne, “S duality at the black hole threshold in gravitational collapse,” Phys. Rev. D 52 (1995) 5397 [arXiv:gr-qc/9505041].
  * [40] E. Hatefi and E. Vanzan, “On Higher Dimensional Self-Similar Axion-Dilaton Solutions,” Eur. Phys. J. C 80 (2020), 10 [arXiv:2005.11646 [hep-th]].
  * [41] Goodfellow, I., Bengio, Y., & Courville, A. “Deep learning”, MIT press. (2016)
  * [42] Ramsundar, B., and Zadeh, R. B. “TensorFlow for deep learning: from linear regression to reinforcement learning,” O Reilly Media, Inc., (2018).
  * [43] Min, S., Lee, B., & Yoon, S., “Deep learning in bioinformatics.” Briefings in bioinformatics. (2017), 851–869.
  * [44] Choudhary, K., et al. “Recent advances and applications of deep learning methods in materials science,” NPJ Computational Materials. (2022)
  * [45] Collobert, R., & Weston, J., “A unified architecture for natural language processing: Deep neural networks with multitask learning,” In Proceedings of the 25-th international conference on Machine learning, (2008), 160–167.
  * [46] Hahnloser, et al. “Digital selection and analogue amplification coexist in a cortex-inspired silicon circuit.” nature (2000): 947-951.

Figure 8: The 95% ANN-based confidence intervals for the critical collapse
functions in the seven-dimensional parabolic class. The ANN lower (dotted
lines), upper (dashed lines) and median (solid lines) bands were computed from
20 replicates with 3000 epochs each. Figure 9: The 95% ANN-based confidence
intervals for the critical collapse functions in the eight-dimensional
parabolic class. The ANN lower (dotted lines), upper (dashed lines) and median
(solid lines) bands were computed from 20 replicates with 3000 epochs each.
Figure 10: The 95% ANN-based confidence intervals for the critical collapse
functions in the nine-dimensional parabolic class. The ANN lower (dotted
lines), upper (dashed lines) and median (solid lines) bands were computed from
20 replicates with 3000 epochs each. Figure 11: The differences between train
and validation loss functions of the developed ANNs in estimating the critical
collapse functions. The top, middle and bottom panels show the loss
differences of the estimates in 7d, 8d and 9d parabolic spaces, respectively.
The panels show the line plot and box plot of loss differences of the ANN
estimates for 20 replicates with 3000 epochs each. Figure 12: The 95% ANN-
based confidence intervals for the critical collapse functions in the five-
dimensional parabolic class. The ANN lower (dotted lines), upper (dashed
lines) and median (solid lines) bands were computed from 100 replicates with
3000 epochs each. Figure 13: The differences between train and validation loss
functions of the developed ANNs in estimating the critical collapse functions
in 5d parabolic class. The panels show the line plot and box plot of loss
differences of the ANN estimates for 100 replicates with 3000 epochs each.

##  Appendix A Solution for the Elliptic Class in 4 dimension

Although we planned to develop the ANNs throughout this manuscript to show no
black hole solution exists in higher dimensional space-time of the parabolic
class, it is instructive to assess if the proposed ANN-based method can detect
if there exists a solution to equations of motion in the case where the
solution is already known in the literature. To do so, we focus on the
elliptic space-time in 4 dimension. In the following, we briefly describe the
equations of motion corresponding to the elliptic class of 4 dimension.

In a similar vein to Section 2, applying CSS ansätze (10) to equations of
motion (6) and (7), one can derive the ordinary differential equations for
$u(z)$, $b(z)$, $f(z)$. Due to the spherical symmetry, $u(z)$ and
$u^{\prime}(z)$ can be eliminated. From (6) and (7), the equations of motion
for self-similar solutions can be re-written as:

$\displaystyle 0$ $\displaystyle=$ $\displaystyle
b^{\prime}+{z(b^{2}-z^{2})\over
b(-1+|f|^{2})^{2}}f^{\prime}\bar{f}^{\prime}-{i\omega(b^{2}-z^{2})\over
b(-1+|f|^{2})^{2}}(f\bar{f}^{\prime}-\bar{f}f^{\prime})-{\omega^{2}z|f|^{2}\over
b(-1+|f|^{2})^{2}},$ $\displaystyle 0$ $\displaystyle=$ $\displaystyle
f^{\prime\prime}-{z(b^{2}+z^{2})\over b^{2}(-1+|f|^{2})^{2}}f^{\prime
2}\bar{f}^{\prime}+{2\over(1-|f|^{2})}\left(1-{i\omega(b^{2}+z^{2})\over
2b^{2}(1-|f|^{2})}\right)\bar{f}f^{\prime 2}$ (34)
$\displaystyle+{i\omega(b^{2}+2z^{2})\over
b^{2}(-1+|f|^{2})^{2}}ff^{\prime}\bar{f}^{\prime}+{2\over z}\left(1+{i\omega
z^{2}(1+|f|^{2})\over(b^{2}-z^{2})(1-|f|^{2})}\right.$
$\displaystyle+\left.{\omega^{2}z^{4}|f|^{2}\over
b^{2}(b^{2}-z^{2})(1-|f|^{2})^{2}}\right)f^{\prime}+{\omega^{2}z\over
b^{2}(-1+|f|^{2})^{2}}f^{2}\bar{f}^{\prime}+$
$\displaystyle{2i\omega\over(b^{2}-z^{2})}\left(\frac{1}{2}-{i\omega(1+|f|^{2})\over
2(1-|f|^{2})}\right.-\left.{\omega^{2}z^{2}|f|^{2}\over
2b^{2}(-1+|f|^{2})^{2}}\right)f.$

The above equations are invariant under a global redefinition of the phase of
$f(z)$. One can solve them in polar coordinates, that is
$f(z)=f_{m}(z)e^{if_{a}(z)}$. The vanishing of the divergent part of
$f^{\prime\prime}(z)$ generates a complex-valued constraint at $z_{+}$,
denoted by $G(b(z_{+}),f(z_{+}),f^{\prime}(z_{+}))=0$. The explicit form of
the $G$ function for the elliptic case in 4d is then given by:

$\displaystyle G(f(z_{+}),f^{\prime}(z_{+}))=$
$\displaystyle\,2z\bar{f}(z_{+})\left(-2\omega^{2}\right)f^{\prime}(z_{+})$
(35)
$\displaystyle+f(z_{+})\bar{f}(z_{+})\left(2z_{+}\bar{f}(z_{+})(-2+2i\omega+2)f^{\prime}(z_{+})+2i\omega\left(2+\omega^{2}\right)\right)$
$\displaystyle-\frac{2z_{+}(2+2i\omega-2)f^{\prime}(z_{+})}{f(z_{+})}+2\omega(\omega-i)f(z_{+})^{2}\bar{f}(z_{+})^{2}-2\omega(\omega+i)\,.$

The regularity of $\tau$ and residual symmetries in the equations of motions
(34) result in the initial boundary conditions as
$b(0)=1,f_{m}^{\prime}(0)=f_{a}^{\prime}(0)=f_{a}(0)=0$.

The behaviours of equations in (34) do not depend on the value of phase of
$f(z)$; hence, one can easily set the boundary conditions as $f_{a}(0)=0$.
Similar to Section 2, applying the profile root-finding, one can easily obtain
$\omega=1.176$ and $f_{m}(0)=0.892$ for four-dimensional space-time. Like in
Section 2, if there is an elliptic solution to the equations of motion, the
solution should appear in a neighbourhood of small space-time values, namely
$z<3$. In this Appendix, we present a numerical study to show if the proposed
ANNs approach can detect the solution in the equations of motion for the
elliptic case in 4d.

Similar to numerical studies of Section 4, we validated the estimates of the
critical collapse functions at 1000 equally spaced points in the interval
$z_{j}=(0,4]$ , $j=1,\ldots,1000$, using the same ANN architecture. We then
computed the ANN estimates using 1000 epochs to reach a good level of
convergence for the loss functions in learning all the critical functions for
$d=4$ dimension.

As discussed in Section 4, the ANNs require a training step to propagate and
back-propagate the gradients for training the parameters; hence the ANN-based
critical collapse functions change slightly from one trial to another,
depending on the trained parameters. To account for this stochastic behaviour,
we ran the ANN method with 1000 epochs in estimating the critical collapse
functions, and then we replicated ANN estimation for $N=20$ times (as in our
previous study for the parabollic class).

We used three statistical measures to evaluate the performance of the ANNs in
estimating the critical collapse functions. Like in Section 4, we first
developed the 95% ANN-based confidence interval for the critical collapse
functions in elliptic space. Figure 14 represents the ANN-based confidence
intervals based on 20 replicates. We also computed the difference between
training and validation losses in the last epoch of the ANNs. Figure 15
represents the line and box plots of the loss differences over 20 replicates.
Finally, we computed the coverage probability to empirically measure the
likelihood that the ANNs approach captures a solution in the domain of the
critical collapse functions. To compute the coverage probability, when ANNs
detect a solution in replicate $i$, we set the coverage index $I_{i}=1$;
otherwise, $I_{i}=0$ for $i=1,\ldots,N$. Finally, the coverage probability is
calculated by $\frac{1}{N}\sum_{i=1}^{N}I_{i}$.

Figure 15 shows that the difference between training and validation loss
functions is almost $10^{-3}$. This indicates that with 1000 epochs the ANNs
have converged in estimating the critical collapse functions in the elliptic
space-time of 4d. From Figure 14, it is seen that the ANNs confirm, at the 95%
confidence level, that there is a solution in the elliptic class of 4d roughly
between $z_{+}\in[0.5,3]$, which is compatible with the analytical solution
$z_{+}=2.605$ from [28]. Finally, we observed that the proposed ANNs
estimation approach detected a solution in the equation of motions with
coverage probability 1 in the numerical experiment. This empirically confirms
the solution in the domain of the critical collapse functions in the elliptic
class of 4d.

Figure 14: The 95% ANN-based confidence intervals for the critical collapse
functions in the four-dimensional elliptic class. The ANN lower (dotted
lines), upper (dashed lines) and median (solid lines) bands were computed from
20 replicates with 1000 epochs each. Figure 15: The differences between train
and validation loss functions of the developed ANNs in estimating the critical
collapse functions in 4d elliptic class. The panels show the line plot and box
plot of loss differences of the ANN estimates for 20 replicates.