#  Frozen stars: Black hole mimickers sourced by a string fluid

Ram Brustein(1), A.J.M. Medved(2,3)
${}^{\textrm{\normalsize(1)\ Department of Physics, Ben-Gurion University,
Beer-Sheva 84105, Israel}}$ ${}^{\textrm{\normalsize(2)\ Department of Physics
\& Electronics, Rhodes University, Grahamstown 6140, South Africa}}$
${}^{\textrm{\normalsize(3)\ National Institute for Theoretical Physics
(NITheP), Western Cape 7602, South Africa}}$<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The frozen star is a non-singular, ultracompact object that, to an external
observer, looks exactly like a Schwarzschild black hole, but with a different
interior geometry and matter composition. The frozen star needs to be sourced
by an extremely anisotropic fluid, for which the sum of the radial pressure
and energy density is either vanishing or perturbatively small. Here, we show
that this matter can be identified with the string fluid resulting from the
decay of an unstable $D$-brane or a brane-antibrane system at the end of open-
string tachyon condensation. The string fluid corresponds to flux tubes
emanating from the center and ending at the Schwarzschild radius of the star.
The effective Lagrangian for this fluid can be recast into a Born-Infeld form.
When the fluid Lagrangian is coupled to that of Einstein’s Gravity, the
static, spherically symmetric solutions of the equations of motion are shown
to be the same as those describing the frozen star model. Frozen stars can
therefore be viewed as gravitationally back-reacted BIons. The Born-Infeld
Lagrangian provides a complete set of equations that describe the dynamics of
the frozen star in a generic state, which is not necessarily static nor
spherically symmetric. Additionally, this description provides a new physical
perspective on the structure of the frozen star in terms of the corresponding
electric fields and charges. The electric field is sourced by a point-like
charge at the center of the star, while its outer layer is equal and
oppositely charged. The electric force between the charges is offset because
the mass of the star is fixed.

##  1 Introduction

The frozen star is a type of black hole (BH) mimicker: a static, spherically
symmetric solution of Einstein’s equations that is a regular and horizonless
alternative to the singular Schwarzschild solution [1, 2, 3, 4, 5].
Furthermore, the frozen star is able to reproduce all of the standard
thermodynamic properties of a Schwarzschild BH of the same mass [6]. It is
also possible to incorporate rotation into the frozen star model [7],
resulting in a regular and horizonless mimicker of the Kerr solution.

The motivation for the interior metric (1) follows from the initial impetus
for the model itself. As discussed at length in [1] (also see [2]), this
solution is meant to be an effective classical description of a highly quantum
state for the object’s interior that is known as the polymer model [8, 9]. In
this picture, the BH mimicker consists of an extremely hot collection of long,
closed, interacting, fundamental strings. One can view the polymer model as
the microscopic description of the frozen star. As the story unfolds, this
string picture will come around full circle, as we arrive at yet another
description of the interior that is motivated by string theory but differently
from the polymer model.

We will show that the frozen star is sourced by a fluid of cold strings
resulting from the decay of an unstable $D$-brane or a brane–antibrane system
at the end of open-string tachyon condensation [10, 11, 12, 13, 14]. Gibbons,
Hori and Yi (GHY) [15] (also see [16]) reformulated Sen’s effective action
describing the state of unstable $D$-branes after the process of tachyon
condensation. Many more references to discussions on this process can be found
in [17]. The reformulated Lagrangian is of a specific Born–Infeld form which
describes a fluid of rigid electric-flux tubes. In [16], the possible bending
and stretching of the flux tubes were also considered. As was duly noted in
[15], a similar form of Lagrangian with a two-form field first appeared in
[18]. After that, the same Lagrangian was proposed in the context of the
cloud-of-strings model in [19] and the string-dust model in [20] and was
discussed in a cosmological context as a “hedgehog compactification” in [21,
22].

When gravity is neglected, the Born–Infeld theory of $Dp$-branes has been
shown to give rise to spherically symmetric, static, solitonic solutions of
finite energy that are known as BIons [23]. (See [24] for a later review.)
BIons were shown to have a point-like source in their core with strings, or
flux tubes, emanating from the core and going all the way to infinity. From
this perspective, frozen stars can be viewed as a specific form of
gravitationally back-reacted BIons whose energy and spatial extent are both
finite.

The Born–Infeld Lagrangian completes the Einstein equations and describes the
dynamics of the frozen star in a generic state, which is not necessarily
static nor spherically symmetric. We can thus show that the frozen star
solution is a consistent and complete model which describes an ultracompact
object whose equilibrium state is practically indistinguishable from that of a
Schwarzschild BH.

Although beyond the scope of the current investigation, the framework that is
developed here should also allow us to study extensions of our model to
rotating frozen stars [7] and to the “defrosted star” model, which allows for
deviations away from $\;\rho+p=0\;$ [5]. The latter is a necessary step for
the frozen star interior to support dynamical modes because of the
ultrastability of the undeformed model [1, 3, 5]. The current framework should
eventually enable us to study extensions of our model to cases in which the
departures from equilibrium physics are macroscopic — the importance of which
has been stressed in [25]. Out-of-equilibrium physics could prove to be
important in describing the dynamics of astrophysical BH mergers, which would
help in distinguishing the frozen star from the Schwarzschild solution, as
well as from other BH mimickers.

##  2 The frozen star

The simplest form of frozen star metric, which is the one that will be
considered in this paper, is as follows:

$ds^{2}\;=\;-\varepsilon^{2}dt^{2}+\frac{1}{\varepsilon^{2}}dr^{2}+r^{2}(d\theta^{2}+\sin^{2}\theta
d\phi^{2})\;,$ (1)

where $\varepsilon^{2}$ is a small, constant and dimensionless parameter.
111The same parameter was referred to as $\epsilon$ or $\varepsilon$ in our
earlier papers on this topic. A recent paper that used data from the Event
Horizon Telescope constrained the parameter to be extremely small,
$\;\varepsilon^{2}\lesssim 10^{-22}\;$ [26]. The outermost surface of the star
is pushed out slightly from the location of the would-be horizon, $\;R\sim
2MG(1+\varepsilon^{2})\;$, where $R$ is the star’s radius and $M$ is its mass.
Additionally, just like the horizon-like outer surface, each radial slice of
the interior is a surface of exponentially large but finite redshift and thus
can be viewed, approximately, as a marginally trapped surface.

The stress–energy tensor $T^{a}_{\;\;b}$ that is needed to source this
simplest frozen star geometry is distinguished by having a radial component of
pressure, $\;p\equiv p_{r}\;$, that takes on the most negative value that is
allowed by causality: $\;p=-\rho\;$, where $\rho$ is the energy density.
Meanwhile, the transverse components of pressure $p_{\perp}\;$ are vanishing,

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle-T^{t}_{\;\;t}\;=\;-\frac{1}{8\pi
G}G^{t}_{\;\;t}\;=\;\frac{1}{8\pi G}\frac{1-\varepsilon^{2}}{r^{2}}\;,$ (2)
$\displaystyle p$ $\displaystyle=$ $\displaystyle
T^{r}_{\;\;r}\;=\;\frac{1}{8\pi G}G^{r}_{\;\;r}\;=\;-\frac{1}{8\pi
G}\frac{1-\varepsilon^{2}}{r^{2}}\;,$ (3) $\displaystyle p_{\perp}$
$\displaystyle=$ $\displaystyle
T^{\theta}_{\;\;\theta}\;=\;T^{\phi}_{\;\;\phi}\;=\;\frac{1}{8\pi
G}G^{\theta}_{\;\;\theta}\;=\;0\;.$ (4)

It follows that

$2GM\;=\;2G\int^{R}dr\;4\pi r^{2}\rho\;=\;R(1-\varepsilon^{2})\;.$ (5)

The frozen star geometry can be viewed as a spherically symmetric collection
of straight, radially pointing, rigid strings (i.e., a “hedgehog”), each with
a constant tension of $\frac{1}{8\pi G}$ [4]. The geometry is mildly singular
near the center of the star, as are $\rho$ and $p$, so that a very small
concentric sphere must be regularized to ensure that these densities remain
finite. This process was described in [4]. Also, as detailed in [3], a
matching process is required near the outermost layer of the star so that the
metric in Eq. (1) and its corresponding stress tensor in Eqs. (2-4) match
smoothly to the exterior Schwarzschild geometry. Later in the paper, we will
offer a different perspective on the regularization and smoothing in terms of
a point-like charge at the center of the star and an equal and opposite charge
that is distributed uniformly over the surface of the star.

A perk of the maximally negative radial pressure is that the frozen star model
is able to evade the singularity theorems of Hawking and Penrose [27, 28]. For
a finite value of $\varepsilon^{2}$, a trapped surface is never actually
formed and having $\;p+\rho=0\;$ means that geodesics do not converge. The way
to understand this is to realize that the equation of state $\;p+\rho=0\;$ can
also be viewed as the saturation of the radial component of the null-energy
condition. Thus, the conditions under which the singularity theorems are valid
are not satisfied by the frozen star geometry. Similarly, a frozen star evades
the “Buchdahl-like” bounds which limit the compactness of matter [29, 30, 31,
32, 33] by having a large negative pressure throughout its interior.

An important characteristic of the frozen star geometry is that the deviations
from the Schwarzschild solution are large throughout the interior of the
object; that is, on horizon-length scales. This goes against common lore that
singularity resolution requires some quantum corrections only near the would-
be singularity, but such reasoning has been shown to lead to energy loss by
radiation that far exceeds the original mass of the object [34, 35].

##  3 Born-Infeld effective Lagrangian and BIons

Here, we review for completeness relevant portions of the analysis of Gibbons,
Hori and Yi (GHY) in [15] (also see [16]), where further details and
references can be found. Our conventions are that an index of 0 denotes time,
one of $a$,$b$,$\cdots$ denotes an arbitrary spacetime dimension and one of
$i$,$j$,$\cdots$ denotes an arbitrary spatial dimension. We assume three large
spatial dimensions for concreteness and, for spherical coordinates,
$\;(0,1,2,3)=(t,r,\theta,\phi)\;$.

The starting point for the analysis in [15] is Sen’s effective action for the
decay of unstable $D$-branes; specifically, at the end of tachyon condensation
and near the minimum of the tachyon potential, where it is vanishing [10, 11,
12, 13, 14]. Sen’s effective Lagrangian (density) can be expressed as

${\cal L}\;=\;-V(T)\sqrt{-Det\left(\eta+2\pi\alpha^{\prime}{\cal
F}\right)}\;+\;\sqrt{-\eta}A^{a}J_{a}\;,$ (6)

where $V(T)$ is the tachyon potential or, equivalently, the $D$-brane tension,
$2\pi\alpha^{\prime}$ is the inverse of the fundamental string tension,
$\;\eta^{a}_{\;\;b}=\delta^{a}_{\;\;b}\;$ is the Minkowski background metric,
$\;{\cal F}_{ab}=\partial_{a}A_{b}-\partial_{b}A_{a}\;$ is the field-strength
tensor for the gauge field $A_{a}$ and $J_{a}$ is the source, a 4-vector
current density. We have not included the kinetic terms for the tachyon field
and for the scalars that are associated with the additional transverse
dimensions. These fields are readily restored if need be (see [16] for formal
details) but are not necessary for the current analysis.

The Lagrangian (6) vanishes when $V(T)$ vanishes; however, the Hamiltonian
${\cal H}$ is well defined. In the case of no magnetic sources, this is

$\;{\cal H}\;=\;D^{i}{\cal F}_{0i}-{\cal L}\;=\;E_{i}D^{i}\;,$ (7)

where the last equality has used the fact that $\;{\cal L}=0\;$ at the minimum
of the potential, $\;E_{i}={\cal F}_{0i}\;$ is the electric field and
$\;D^{i}$ is the electric displacement which, in the current case, is the
canonical conjugate of the gauge field, $\;D^{i}=\frac{\delta{\cal
L}}{\delta(\partial_{0}A_{i})}\;$. The displacement $D^{i}$ is naturally
preferred over $E^{i}$ to play the role of the “electric field” because it is
the field in a Born–Infeld theory that always satisfies the Gauss’-law
constraint,

$\partial_{i}D^{i}\;=\;J_{0}\;=\;\rho_{e}\;.$ (8)

More generally, the Hamiltonian can be expressed as

${\cal H}\;=\;\frac{\delta{\cal L}}{\delta(\partial_{0}A_{i})}-{\cal
L}\;=\;\frac{1}{2\pi\alpha^{\prime}}\sqrt{D^{i}D_{i}+{\cal P}^{i}{\cal
P}_{i}}\;,$ (9)

where $\;{\cal P}_{i}=-{\cal
F}_{ij}D^{j}=\left(\vec{D}\times\vec{B}\right)_{i}\;$ is the conserved
momentum associated with spatial translations. The definition of the magnetic
induction is standard, $\;B^{i}=\frac{1}{2}\epsilon^{ijk}{\cal F}_{jk}\;$.

The sources can be included implicitly by imposing the Gauss’-law constraint
(8), Ampere’s law $\;\partial_{i}{\cal
F}^{i}_{\;\;j}-\partial_{0}{E}_{j}=J_{j}\;$, along with the Bianchi
identities, which include Faraday’s law
$\;{\nabla\times}\;\vec{E}+\partial_{0}{\vec{B}}=0\;$ and

$\partial_{i}B^{i}\;=\;0\;.$ (10)

The last equation assumes that there are no sources for magnetic monopoles.
The inline equations above Eq. (10) are of less importance here, as our main
concern is the case for which magnetic sources and time-dependent fields are
absent.

To obtain a useful Lagrangian, GHY follow techniques from [36] and regard the
magnetic fields ${\cal F}_{ij}$ as the conjugates with respect to a new set of
dual variables $K_{ij}$ such that $\;K^{ij}=2\frac{\delta{\cal H}}{\delta{\cal
F}_{ij}}$. The fields ${\cal F}_{ab}$ and ${\cal K}_{ab}$ — the latter being
an extension of $K_{ij}$ as defined in Eq. (14) — should be regarded as
independent variables. This will be relevant when we derive the equations of
motion (EOM) for the tensor field ${\cal K}$.

The resulting Lagrangian is then defined by an appropriate Legendre
transformation,

${\cal L}^{\prime}\;=\;{\cal H}-\frac{1}{2}{\cal
F}_{ij}K^{ij}\;=\;\frac{1}{2\pi\alpha^{\prime}}\sqrt{D^{i}D_{i}-\frac{1}{2}K^{ij}K_{ij}}\;,$
(11)

and is clearly non-vanishing in general. The latter equality can be obtained
using the relation

$\frac{1}{2}K^{ij}K_{ij}\;=\;\frac{1}{{\cal H}^{2}}{\cal P}^{i}{\cal
P}_{i}D^{j}D_{j}\;.$ (12)

GHY then define a two-form field ${\cal
K}_{ab}=\partial_{a}\widetilde{A}_{b}-\partial_{b}\widetilde{A}_{a}\;$ that
acts as an effective field strength,

$\displaystyle{\cal K}_{0i}$ $\displaystyle=$ $\displaystyle D_{i}\;,$ (13)
$\displaystyle{\cal K}_{ij}$ $\displaystyle=$ $\displaystyle K_{ij}\;.$ (14)

The new Lagrangian can now be written in a manifestly Born–Infeld form,

${\cal L}^{\prime}\;=\;\frac{1}{2\pi\alpha^{\prime}}\sqrt{-\frac{1}{2}{\cal
K}^{ab}{\cal K}_{ab}}\;,$ (15)

where the negative sign is a consequence of the time–time component of the
metric appearing in the contraction when $\;a=0\;$ and $\;b=i\;$ (and vice
versa). Here, the (effective) electric-field term is presumed to be the
dominant one. In the case of no magnetic sources, one finds that the new
Lagrangian is the same as the original Hamiltonian, $\;{\cal
L}^{\prime}\;=\;E_{i}D^{i}\;$.

There is a subtlety in this procedure in that Eq. (14) implies that $\;{\cal
K}\wedge{\cal K}=0\;$ but, given that this constraint is in effect, the
canonical analysis of ${\cal L}^{\prime}$ does not lead back to the same
Hamiltonian ${\cal H}$. As explained in [15], this situation can be rectified
by adding a Lagrange-multiplier term to ${\cal L}^{\prime}$ that imposes the
constraint explicitly.

To summarize, the Born–Infeld Lagrangian describing the string fluid is the
following:

${\cal L}^{\prime}\;=\;\frac{1}{2\pi\alpha^{\prime}}\sqrt{-\frac{1}{2}{\cal
K}^{ab}{\cal K}_{ab}}\;+\;\lambda_{1}\epsilon^{abcd}{\cal K}_{ab}{\cal
K}_{cd}\;+\sqrt{-\eta}\;A^{a}J_{a}\;,$ (16)

where $\lambda_{1}$ is the Lagrange multiplier and the constraint is
automatically satisfied (and the constraint term vanishes) by imposing Eq.
(14). The connection with strings follows from the fact that ${\cal K}^{ab}$
can be identified as a surface-forming bivector [15], which can then be
interpreted as a cross-sectional slice of the world sheet of an open string or
a flux tube.

The EOM for this Lagrangian, $\;d{\cal L}^{\prime}=0\;$, are most transparent
when expressed in terms of the original field-strength tensor as these are
equivalent to the original Bianchi identities [15],

$d{\cal F}\;=\;0\;.$ (17)

On the other hand, the Bianchi identities for the new field-strength tensor

$d{\cal K}\;=\;0\;$ (18)

are equivalent to the EOM for the original Lagrangian [15].

There is also a subtlety concerning the source term. The original source term
in Eq. (6) is unaffected by the Legendre transformation; however, $A_{a}$ is
not the gauge field for the field strength ${\cal K}_{ab}$. Let us denote its
gauge field as $\widetilde{A}_{a}$, then
$\;D_{i}=\partial_{i}\widetilde{A}^{0}\;$ whereas
$\;E_{i}=\partial_{i}A^{0}\;$. Fortunately, this tension is resolved because,
as we have noted, the two types of field strengths are independent, as must
also be true of their respective gauge fields. It follows that we can vary the
Lagrangian ${\cal L}^{\prime}$ with respect to either gauge field, and it
happens to be the variation with respect to $A^{a}$ that leads to the expected
Gauss’s-law constraint in Eq. (8). This is most clear in the case of no
magnetic sources, as can be seen from the form of the inline equation for
${\cal L}^{\prime}$ below Eq. (15).

If there are no bulk sources — our case of particular interest — the
stress–energy tensor is given by

$\displaystyle T_{ab}\;=\;2\frac{\delta{\cal
L}^{\prime}}{\delta\eta^{ab}}\;=\;\frac{1}{2\pi\alpha^{\prime}}\frac{{\cal
K}_{a}^{\;\;c}{\cal K}_{bc}}{\sqrt{-\frac{1}{2}{\cal K}^{ab}{\cal
K}_{ab}}}\;.$ (19)

That this is the appropriate definition of $T_{ab}$ is justified in [19]. In
terms of the effective electric fields, magnetic fields and momenta,

$\displaystyle T_{00}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi\alpha^{\prime}}\frac{D^{i}D_{i}}{\sqrt{D^{i}D_{i}-\frac{1}{2}K^{ij}K_{ij}}}\;=\;{\cal
H}\;,$ (20) $\displaystyle T_{0i}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi\alpha^{\prime}}\frac{D^{j}K_{ij}}{\sqrt{D^{i}D_{i}-\frac{1}{2}K^{ij}K_{ij}}}\;=\;-{\cal
P}_{i}\;,$ (21) $\displaystyle T_{ij}$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi\alpha^{\prime}}\frac{-D^{i}D_{j}+K_{i}^{\;\;k}K_{jk}}{\sqrt{\pi^{i}\pi_{i}-\frac{1}{2}K^{ij}K_{ij}}}\;=\;\frac{1}{2\pi\alpha^{\prime}}\frac{-D_{i}D_{j}+{\cal
P}_{i}{\cal P}_{j}}{{\cal H}}\;,$ (22)

where the right-most relations make use of Eq. (12).

The simplest solution of the EOM is a static, spherically symmetric
configuration, which happens when $\;{\cal K}_{01}=-{\cal K}_{10}\;$ and all
other elements of this field strength vanish. In which case, $\;{\cal
K}^{ab}{\cal K}_{ab}=-2D_{1}D^{1}$. Then the only non-vanishing elements of
$T_{ab}$ are

$T_{00}\;=-\;T_{11}\;=\;\frac{1}{2\pi\alpha^{\prime}}\sqrt{D^{1}D_{1}}\;,$
(23)

implying that $\;p=-\rho\;$ and $\;p_{\perp}=0\;$. This is the so-called
string fluid [15], for which

$E_{1}\;=\;\frac{1}{2\pi\alpha^{\prime}}\frac{D_{1}}{\sqrt{(D_{1})^{2}}}\;=\;\frac{1}{2\pi\alpha^{\prime}}\;,$
(24)

and so

$T_{00}\;=\;E_{1}D^{1}\;.$ (25)

Importantly, the bulk portion (or first term) of the Lagrangian (16) is
exactly the same as $T_{00}$,

${\cal
L}^{\prime}_{bulk}\;=\;\frac{1}{2\pi\alpha^{\prime}}\sqrt{D^{1}D_{1}}\;=\;E_{1}D^{1}\;.$
(26)

The static and spherically symmetric case is also related to the BIon solution
[23], for which $\;D_{1}=\frac{q}{4\pi r^{2}}\;$, corresponding to a point-
like charge at the origin $\;\partial_{i}D^{i}=q\leavevmode\nobreak\
\delta^{(3)}(\vec{r})\;$. BIons can therefore be viewed as being sourced by a
fluid of electric flux lines that are emanating radially from a point source
at the center and extending all the way to infinity.

##  4 Frozen stars as gravitationally backreacted BIons

In this section, we couple the BIons to Einstein’s gravity, discuss the
backreaction on the BIons and, finally, show that they correspond to solutions
of the frozen star model.

The complete gravitational and matter action includes, in addition to the
Einstein–Hilbert and Born–Infeld Lagrangians, the constraint and source terms

$S_{GBI}\;=\;\int
d^{4}x\left\\{\sqrt{-g}R+\frac{1}{2\pi\alpha^{\prime}}\sqrt{-\frac{1}{2}{\cal
K}^{ab}{\cal K}_{ab}}+\lambda_{1}\epsilon^{abcd}{\cal K}_{ab}{\cal
K}_{cd}+\sqrt{-g}J_{a}A^{a}\right\\}\;.$ (27)

As discussed at the end of this section, a second constraint which fixes the
mass is still to be included. But, as a surface term, it will not affect the
bulk EOM. We will also be assuming that there are no sources in the bulk.

The gravitational EOM are $\;T^{a}_{\;\;b}=\frac{1}{8\pi G}G^{a}_{\;\;b}\;$,
with $T^{a}_{\;\;b}$ given in Eq. (19), again because there are no bulk
sources. The Born–Infeld EOM are as discussed in the previous section, except
partial derivatives should be replaced by covariant derivatives (although
often inconsequential because of the anti-symmetry properties of the field-
strength tensor). The same applies to the Bianchi identities and source
constraints.

To reproduce the frozen star solution, we will need to choose specific source
terms and impose certain boundary conditions such that the total charge of the
frozen star vanishes and its total mass is fixed (the latter constraint will
be discussed further on). The source current $J_{a}$, for which only $J_{0}$
is non-vanishing, therefore includes a point-like positive electric charge at
the star’s center and a negative charge of equal magnitude which is spread
evenly over the outer surface of the star.

One might be concerned that the attractive electric force between the opposite
charges in the core and the outer layer will endanger its stability. However,
as discussed later, fixing the mass of the star exactly offsets this
attractive force. The reason being that a fixed mass for the solution
translates into a fixed radius, so that the outer layer cannot move inwards
(or outwards) in response to a supplementary force. This was a key ingredient
in the discussion of the stability of the frozen star [1, 3, 5].

The frozen star solution can be related to the BIon solution. The matching can
be made precise for the static, spherically symmetric case by using the
${}^{t}_{\ t}$ component of the Einstein equations along with Eq. (23),

$\rho_{FS}\;=\;\frac{1}{8\pi
G}\frac{1-\varepsilon^{2}}{r^{2}}\;=\;\frac{1}{2\pi\alpha^{\prime}}\sqrt{D^{1}D_{1}}\;,$
(28)

such that $\rho_{FS}$ is the energy density as expressed in Eq. (2). The
result is

$\sqrt{D^{1}D_{1}}\;=\;\frac{\alpha^{\prime}}{4G}\frac{1}{r^{2}}+{\cal
O}[\varepsilon^{2}]\;,$ (29)

corresponding to the existence of an electric point-like charge at the center
of the star, since the solution follows from
$\;\nabla_{i}D^{i}=\nabla_{1}D^{1}=q_{core}\leavevmode\nobreak\
\delta^{(3)}(\vec{r})\;.$ It can then be deduced that

$q_{core}\;=\;\pi\frac{\alpha^{\prime}}{G}\;.$ (30)

An interesting feature of the solution is that the ratio
${\alpha^{\prime}}/{G}$ in weakly coupled string theories scales as $\;\sim
1/g_{s}^{2}\;$, where $g_{s}^{2}$ is the closed string coupling. Moreover, the
description of the solution as radial flux lines fits in well with the
hedgehog picture of the frozen star that was recalled in Section 2.

The total charge of the frozen star needs to vanish, as its energy density,
and thus its electric field, vanishes in the exterior region where the
geometry is described by the vacuum Schwarzschild solution. Correspondingly,
there has to be a charge $q_{out}$ that is spread out uniformly over the
star’s outer surface layer and equal to $-q_{core}$. This leads to a source
term at the outer surface as follows:

$\left(J_{0}\right)_{out}\;=\;\frac{q_{out}}{4\pi
r^{2}}\delta(r-R)\;=\;-\frac{q_{core}}{4\pi r^{2}}\delta(r-R)\;.$ (31)

Figure 1: A frozen star sourced by a fluid of electric flux lines that are
emanating radially from a point-like charge in its core and ending on its
outer layer which is oppositely charged. The electric force between the
charges is offset because the mass of the star is fixed.

Let us now discuss the boundary condition that fixes the mass of the frozen
star. This needs to be formally imposed at infinity, as this is the location
where the mass is defined by way of a surface integral. However, because the
solution outside of the frozen star is strictly a vacuum, the mass can also
evaluated as a surface integral at the outermost radius of the star, or,
alternatively, as a volume integral that enforces the constraint $\;\int
d^{3}x\;\sqrt{g_{22}g_{33}}\;\rho_{FS}=M\;$, as explained in some discussions
about BH mechanics [37, 38]. The latter option entails adding a second
Lagrange-multiplier term of the form (initially expressed as a scalar, not a
density)

$\Delta{L}\;=\;\lambda_{2}\left(\int_{0}^{R}\left[dr\leavevmode\nobreak\ 4\pi
r^{2}\leavevmode\nobreak\ \rho_{FS}\right]\;-\;M\right)\;.$ (32)

Using Eq. (28), we may rewrite the constraint term as

$\Delta{L}\;=\;\lambda_{2}\left(\int_{0}^{R}\left[dr\leavevmode\nobreak\ 4\pi
r^{2}\leavevmode\nobreak\ E_{1}D^{1}\right]\;-\;M\right)\;,$ (33)

which leads to the following Lagrangian density $\Delta{\cal L}$:

$\Delta{\cal L}\;=\;\lambda_{2}\frac{\delta\left(r-R\right)}{4\pi
r^{2}}\left(\int_{0}^{R}\left[dr\leavevmode\nobreak\ 4\pi
r^{2}\leavevmode\nobreak\ E_{1}D^{1}\right]\;-\;M\right)\;.$ (34)

We have, for simplicity, assumed a static and spherically symmetric solution
in expressing the above constraint term, as this will allow us to show
explicitly how the electric force between the charges at the core and the
outer surface is canceled by the mass constraint. However, it is possible to
express the mass constraint covariantly, along the lines of discussions on
isolated and dynamical horizons such as [38]. The outer layer of the frozen
star can be treated as a marginally trapped surface and, therefore, the
formalism for calculating the mass as a covariant surface integral is
applicable. We leave this to a future investigation.

We can solve for $\lambda_{2}$ by varying $\Delta{\cal L}$ with respect to the
original gauge field $A^{0}(R)$, as described in the previous section.
222Note: The integrand factor of $\;4\pi r^{2}=4\pi\sqrt{-g}\;$ is unaffected
by a covariant derivative. The equation of motion that results from varying by
$A^{0}(R)$ at the outer surface, $\;\frac{\delta\Delta{\cal L}}{\delta
A^{0}(R)}+J_{0}(R)=0\;$, along with the Gauss’-law constraint (8) and the
relation between the charges (31), leads to $\;\lambda_{2}=-1\;$, which simply
means that the constraint force exactly cancels the attractive electric force
on the outer charge distribution.

##  5 Summary and outlook

We have shown how our frozen star model for an ultracompact BH mimicker can be
described effectively as the spherically symmetric solution of Einstein’s
gravity coupled to the GHY form of Born–Infeld Lagrangian, which uses Sen’s
effective action for $D$-brane decay as its starting point.

The current framework can be extended in several directions. Including
rotation is straightforward. Let us recall the zero-angular-momentum-
observer’s form of the stress tensor for the rotating star, as presented in
[7]. In this case $\;p=-\rho\;$ and $\;p_{\perp}\sim{\cal
O}[\varepsilon^{2}]\;$, so that the corresponding axially symmetric solution
can be found with only an electric field and vanishing magnetic fields. The
only difference from the static, spherically symmetric case is the dependence
of the electric field on $r$ and $\theta$ rather than on $r$ alone.

To connect this framework with the defrosted star, for which $\;\rho+p\;$ and
$p_{\perp}$ are no longer vanishing but are perturbatively small, will require
more work. In addition to a spherical electric field, a weak spherical
magnetic field will be required. This can be realized by adding a magnetic
monopole source, meaning that the star is sourced by a dyon, as both its
electric and magnetic fluxes are in the radial direction.

It is also of interest that the Born–Infeld model can be thought of as
describing a fluid of open strings, given that the frozen star model has its
conceptual origins as a classical version of our earlier-proposed polymer
model, which contains an extremely hot fluid of closed strings. It is tempting
to conjecture that there is some sort of duality at work here. Nevertheless,
of the pair, it is the Born–Infeld description that has an emphatic advantage;
namely, a well-defined mathematical framework, which we look forward to
exploiting.

## Acknowledgments

We thank Yoav Zigdon for discussions, Suvendu Giri for pointing to us the
similarity of the frozen star model to the clouds of string model and to the
Born-Infeld BIons and Frans Pretorius for insisting that we find the source
matter Lagrangian of the frozen star. We extend our special thanks to Piljin
Yi, for helping us understand his work and its implications to the frozen star
model. The research is supported by the German Research Foundation through a
German-Israeli Project Cooperation (DIP) grant “Holography and the Swampland”
and by VATAT (Israel planning and budgeting committee) grant for supporting
theoretical high energy physics. The research of AJMM received support from an
NRF Evaluation and Rating Grant 119411 and a Rhodes Discretionary Grant
SD07/2022. AJMM thanks Ben Gurion University for their hospitality during past
visits.

## References

  * [1] R. Brustein and A. J. M. Medved, “Resisting collapse: How matter inside a black hole can withstand gravity,” Phys. Rev. D 99, no.6, 064019 (2019) [arXiv:1805.11667 [hep-th]].
  * [2] R. Brustein and A. J. M. Medved, “Non-Singular Black Holes Interiors Need Physics Beyond the Standard Model,” Fortsch. Phys. 67, no.10, 1900058 (2019) [arXiv:1902.07990 [hep-th]].
  * [3] R. Brustein, A. J. M. Medved and T. Simhon, “Black holes as frozen stars,” Phys. Rev. D 105, no.2, 024019 (2022) [arXiv:2109.10017 [gr-qc]].
  * [4] R. Brustein, A. J. M. Medved, T. Shindelman and T. Simhon, “Black Holes as Frozen Stars: Regular Interior Geometry,” Fortsch. Phys. 72, no.1, 2300188 (2024) [arXiv:2301.09712 [gr-qc]].
  * [5] R. Brustein, A. J. M. Medved and T. Shindelman, “Defrosting frozen stars: spectrum of internal fluid modes,” Phys. Rev. D 108, no.4, 044058 (2023) [arXiv:2304.04984 [gr-qc]].
  * [6] R. Brustein, A. J. M. Medved and T. Simhon, “Thermodynamics of frozen stars,” [arXiv:2310.11572 [gr-qc]].
  * [7] R. Brustein and A. J. M. Medved, “Sourcing the Kerr geometry,” [arXiv:2310.16467 [gr-qc]].
  * [8] R. Brustein and A. J. M. Medved, “Black holes as collapsed polymers,” Fortsch. Phys. 65, no. 1, 1600114 (2017) [arXiv:1602.07706 [hep-th]].
  * [9] R. Brustein and A. J. M. Medved,“Emergent horizon, Hawking radia- tion and chaos in the collapsed polymer model of a black hole,” Fortsch. Phys. 65, 0116 (2017) [arXiv:1607.03721 [hep-th]].
  * [10] A. Sen, “Tachyon condensation on the brane anti-brane system,” JHEP 08, 012 (1998) [arXiv:hep-th/9805170 [hep-th]].
  * [11] A. Sen, “BPS D-branes on nonsupersymmetric cycles,” JHEP 12, 021 (1998) [arXiv:hep-th/9812031 [hep-th]].
  * [12] A. Sen, “NonBPS states and Branes in string theory,” [arXiv:hep-th/9904207 [hep-th]].
  * [13] A. Sen, “Supersymmetric world volume action for nonBPS D-branes,” JHEP 10, 008 (1999) [arXiv:hep-th/9909062 [hep-th]].
  * [14] A. Sen, “Universality of the tachyon potential,” JHEP 12, 027 (1999) [arXiv:hep-th/9911116 [hep-th]].
  * [15] G. W. Gibbons, K. Hori and P. Yi, “String fluid from unstable D-branes,” Nucl. Phys. B 596, 136 (2001) [arXiv:hep-th/0009061 [hep-th]].
  * [16] H. U. Yee and P. Yi, “Open / closed duality, unstable D-branes, and coarse grained closed strings,” Nucl. Phys. B 686, 31 (2004) [arXiv:hep-th/0402027 [hep-th]].
  * [17] A. Sen, “Tachyon dynamics in open string theory,” Int. J. Mod. Phys. A 20, 5513 (2005) [arXiv:hep-th/0410103 [hep-th]].
  * [18] H. B. Nielsen and P. Olesen, “Local field theory of the dual string,” Nucl. Phys. B 57, 367 (1973).
  * [19] P. S. Letelier, “Clouds Of Strings In General Relativity,” Phys. Rev. D 20 1294 (1979).
  * [20] J. Stachel, “Thickening The String. I. The String Perfect Dust,” Phys. Rev. D 21, 2171 (1990).
  * [21] E. I. Guendelman and A. Rabinowitz, “The Gravitational field of a hedgehog and the evolution of vacuum bubbles,” Phys. Rev. D 44, 3152 (1991).
  * [22] E. I. Guendelman and A. I. Rabinowitz, “Hedgehog compactification,” Phys. Rev. D 47, 3474 (1993) [erratum: Phys. Rev. D 48, 2961 (1993)].
  * [23] G. W. Gibbons, “Born-Infeld particles and Dirichlet p-branes,” Nucl. Phys. B 514, 603 (1998) [arXiv:hep-th/9709027 [hep-th]].
  * [24] G. W. Gibbons, “Aspects of Born-Infeld theory and string / M theory,” AIP Conf. Proc. 589, no.1, 324 (2001) [arXiv:hep-th/0106059 [hep-th]].
  * [25] R. Brustein and A. J. M. Medved, “Quantum hair of black holes out of equilibrium,” Phys. Rev. D 97, no.4, 044035 (2018) [arXiv:1709.03566 [hep-th]].
  * [26] R. Carballo-Rubio, F. Di Filippo, S. Liberati and M. Visser, “Constraints on thermalizing surfaces from infrared observations of supermassive black holes,” [arXiv:2306.17480 [astro-ph.HE]].
  * [27] R. Penrose, “Gravitational Collapse and Space-Time Singularities,” Phys. Rev. Lett. 14, 57 (1965).
  * [28] S. W. Hawking and R. Penrose, “The singularities of gravitational collapse and cosmology,” Proc. R. Soc. Lond. A 314, 529 (1970).
  * [29] H. Buchdahl, “General Relativistic Fluid Spheres,” Phys. Rev. 116, 1027 (1959).
  * [30] S. Chandrasekhar “Dynamical Instability of Gaseous Masses Approaching the Schwarzschild Limit in General Relativity,” Phys. Rev. Lett. 12, 114 (1964).
  * [31] S. Chandrasekhar, “The Dynamical Instability of Gaseous Masses Approaching the Schwarzschild Limit in General Relativity,” Astrophys. J. 140, 417 (1964).
  * [32] H. Bondi, “Massive spheres in general relativity,” Proc. Roy. Soc. Lond. A 282, 303 (1964).
  * [33] P. O. Mazur and E. Mottola, “Surface tension and negative pressure interior of a non-singular ,” Class. Quant. Grav. 32, no. 21, 215024 (2015) [arXiv:1501.03806 [gr-qc]].
  * [34] V. P. Frolov and A. Zelnikov, “Quantum radiation from an evaporating nonsingular black hole,” Phys. Rev. D 95, no. 12, 124028 (2017) [arXiv:1704.03043 [hep-th]].
  * [35] R. Carballo-Rubio, F. Di Filippo, S. Liberati, C. Pacilio and M. Visser, “On the viability of regular black holes,” JHEP 1807, 023 (2018) [arXiv:1805.02675 [gr-qc]].
  * [36] A. A. Tseytlin, “Selfduality of Born-Infeld action and Dirichlet three-brane of type IIB superstring theory,” Nucl. Phys. B 469, 51 (1996) [arXiv:hep-th/9602064 [hep-th]].
  * [37] J. M. Bardeen, B. Carter and S. W. Hawking, “The Four laws of black hole mechanics,” Commun. Math. Phys. 31, 161 (1973).
  * [38] A. Ashtekar and B. Krishnan, “Isolated and dynamical horizons and their applications,” Living Rev. Rel. 7, 10 (2004) [arXiv:gr-qc/0407042 [gr-qc]].


# A Systematic Bias of Machine Learning Regression Models and Its Correction:
an Application to Imaging-based Brain Age Prediction

Hwiyoung Lee  Maryland Psychiatric Research Center, University of Maryland,
School of Medicine  University of Maryland, Institute for Health Computing
Shuo Chen  Maryland Psychiatric Research Center, University of Maryland,
School of Medicine  University of Maryland, Institute for Health Computing

###### Abstract

Machine learning models for continuous outcomes often yield systematically
biased predictions, particularly for values that largely deviate from the
mean. Specifically, predictions for large-valued outcomes tend to be
negatively biased, while those for small-valued outcomes are positively
biased. We refer to this linear central tendency warped bias as the
“systematic bias of machine learning regression”. In this paper, we first
demonstrate that this issue persists across various machine learning models,
and then delve into its theoretical underpinnings. We propose a general
constrained optimization approach designed to correct this bias and develop a
computationally efficient algorithm to implement our method. Our simulation
results indicate that our correction method effectively eliminates the bias
from the predicted outcomes. We apply the proposed approach to the prediction
of brain age using neuroimaging data. In comparison to competing machine
learning models, our method effectively addresses the longstanding issue of
“systematic bias of machine learning regression” in neuroimaging-based brain
age calculation, yielding unbiased predictions of brain age.

##  1 Introduction

Constructing predictive models with continuous outcomes is a fundamental
aspect of modern data science. Numerous tools have been developed for this
purpose, including statistical methods such as ordinary linear regression,
regression shrinkage, and Generalized Additive Models (GAM), as well as
machine learning methods like random forests, XGBoost, and support vector
regression, among others (Hastie et al.,, 2009). A general objective of these
methods is to minimize the discrepancy between the predicted continuous
outcomes and true values, particularly in independent testing datasets. Using
this heuristic, the predicted outcome is unbiased under the classic linear
regression setting. In contrast, the predicted outcome from a machine learning
model for continuous outcomes is often systematically biased (Zhang and Lu,,
2012; Belitz and Stackelberg,, 2021). This systematic bias is problematic for
the applications of machine learning regression models, leading to inaccurate
conclusions and forecasts.

We first illustrate the machine learning regression introduced bias using a
simple simulation study and multiple machine learning models including Kernel
Ridge Regression (KRR), Lasso, XGBoost, and Random Forest. Specifically, we
first generate synthetic training and testing data, each comprising $n=1,000$
observations. The predictors $\mathbb{X}$ were generated from a multivariate
normal distribution with dimension $p=200$ (i.e., $\mathbb{X}\sim
N(0,\Sigma)$, where $\Sigma=I_{p}$. The response $\mathbb{y}$ is linearly
related to the predictor $\mathbb{X}$, i.e.,
$\mathbb{y}=\mathbb{X}\boldsymbol{\beta}+\epsilon$. Here,
$\boldsymbol{\beta}\in\mathbb{R}^{p}$ is a vector of regression coefficients,
with each coefficient set to 0.1, and $\epsilon$ denotes random noise, which
follows the standard normal distribution.

Figure 1 presents a scatter plot of the predicted values against the observed
responses. A solid line representing the linear regression fit of the
predicted $\widehat{\mathbb{y}}$ to $\mathbb{y}$ is also shown, whereas the
dotted line has a slope of 1. When $\widehat{\mathbb{y}}$ is unbiased, the
slope of the solid line is also 1, otherwise biased. As demonstrated in this
figure, a systematic bias arises, characterized by a tendency to predict
larger values when $y$ is small (overestimation) and to predict lower values
when $y$ is large (underestimation). The sum of these biases for all
observations is approximately zero due to the optimization of the commonly
used objective function. This linear biased tendency generally presents across
all machine learning models and is more apparent in the testing dataset.

Since the systematic bias shows a consistent center-warping and linear
tendency across all machine learning models, we refer to this phenomenon as
the “systematic bias of machine learning regression” (SBMR). We denote the
linear trend of this systematic bias by $c=sin(\theta)$, where $\theta$ is the
angle between the dotted line and the solid line. In the following 1, we show
that the mean squared error of the biased prediction from a machine learning
regression model is smaller than the unbiased prediction. This theoretically
explains why the aforementioned systematic bias is preferred by all machine
learning regression models with the objective function to minimize the mean
squared error.

######  Proposition 1.

Consider the outcome $Y$ with a mean of zero and variance of $\sigma^{2}$,
$\widehat{y}_{i}$ as the unbiased prediction of $y_{i}$ ($i=1,\cdots,n$), and
the systematically biased machine learning regression prediction
$\tilde{y_{i}}=\widehat{y_{i}}-c\widehat{y_{i}}$, where $0<c<1$. Then for $c$,
and the coefficient of determination $R^{2}$ satisfying
$c<\frac{2-2R^{2}}{2-R^{2}}$, we have

$\displaystyle\mathbb{E}\left[\sum_{i=1}^{n}\left(y_{i}-\widehat{y}_{i}\right)^{2}\right]\geq\mathbb{E}\left[\sum_{i=1}^{n}\left(y_{i}-\tilde{y}_{i}\right)^{2}\right].$

This can be straightforwardly checked by using
$\mathbb{E}[\sum_{i=1}^{n}\left(y_{i}-\widehat{y_{i}}\right)^{2}]=(n-1)\sigma^{2}(1-R^{2})$,
and
$\mathbb{E}[\sum_{i=1}^{n}\left(y_{i}-\tilde{y}_{i}\right)^{2}]=(n-1)c^{2}\sigma^{2}+(n-1)(1-c)^{2}\sigma^{2}(1-R^{2})$.

Although this issue frequently arises during the regression of continuous
variables, it has received little attention in general statistics and machine
learning literature. However, in the neuroimaging community, especially in the
context of predicting brain age—calculated by regressing chronological age on
neuroimaging data—this phenomenon is widely acknowledged. Many studies have
endeavored to address this issue; the majority of proposed solutions focus on
the implementation of a post-bias correction step Smith et al., (2019) after
fitting the model, rather than addressing it at the level of the main
objective function. However, post-bias correction is often criticized Butler
et al., (2021) because this additional step relies on certain assumptions
about the distributions or relationships within the data. Thus, if these
assumptions are violated, the correction may not only fail to accurately
address the issue but could potentially worsen it. Nevertheless, to the best
of our knowledge, there are no other methodologies besides the bias correction
approach. To fill this gap and tackle the SBMR issue efficiently, we propose a
novel approach that addresses this problem at the objective function level,
without requiring any subsequent steps. This is achieved by imposing
constraints designed to prevent the fitted values from being forced toward the
mean.

Figure 1: Illustration of regression to the mean: We simulated 1,000
observations for each training and testing set. We compared our proposed
methods (represented in orange panels) with conventional machine learning
methods (shown in grey panels), including Kernel Ridge Regression (KRR), Lasso
Regression, XGBoost, and Random Forest (RF).

##  2 Method

In this section, we propose a solution that addresses the SBMR issue by
imposing a linear constraint on two existing methods: Lasso and KRR. This
additional constraint prevents the fitted values from being forced toward the
overall mean of the data.

Specifically, we partition the data into two groups based on whether the
response values are less than or greater than the overall mean of the training
set. We then adjust the fitted values so that the mean of the fitted values
for observations less than the overall mean equals the mean of these
observations, and similarly for observations greater than the overall mean.
This ensures that the mean of the fitted values corresponds appropriately. We
impose these constraints for each method, which are elaborated in the
following subsections.

###  2.1 Constrained Lasso

We first impose this penalty on the linear model. The regularized linear
regression has the following form:

$\displaystyle\operatornamewithlimits{argmin}_{\beta\in\mathbb{R}^{p}}\sum_{i=1^{n}}\left(\mathbb{y}_{i}-\mathbb{x}_{i}^{\top}\boldsymbol{\beta}\right)^{2}+\lambda\Omega(\boldsymbol{\beta}),$
(1)

where $\Omega$ is a penalty function. One popular choice is $\ell_{1}$
penalty, i.e.,
$\|\boldsymbol{\beta}\|_{1}=\sum_{j=1}^{p}|\boldsymbol{\beta}_{j}|$. We denote
the overall mean of the response as
$\overline{\mathbb{y}}=\frac{1}{n}\sum_{i=1}^{n}\mathbb{y}_{i}$. The response
variable can then be divided into two groups based on the overall training
mean: (i) those less than the overall training mean, i.e.,
$\mathbb{I}_{<}=\\{i\in n|\mathbb{y}_{i}<\overline{\mathbb{y}}\\}$, and (ii)
those greater than the overall training mean, i.e., $\mathbb{I}_{>}=\\{i\in
n|\mathbb{y}_{i}>\overline{\mathbb{y}}\\}$. The corresponding means for each
group are
$\overline{\mathbb{y}}_{<}=\frac{1}{|\mathbb{I}_{<}|}\sum_{i\in\mathbb{I}_{<}}\mathbb{y}_{i}$,
and
$\overline{\mathbb{y}}_{>}=\frac{1}{|\mathbb{I}_{>}|}\sum_{i\in\mathbb{I}_{>}}\mathbb{y}_{i}$.
Accordingly, we can divide the predictors into two sets:
$\mathbb{X}_{<}=[\mathbb{x}_{i}:i\in\mathbb{I}_{<}]\in\mathbb{R}^{|\mathbb{I}_{<}|\times
p}$, and
$\mathbb{X}_{>}=[\mathbb{x}_{i}:i\in\mathbb{I}_{>}]\in\mathbb{R}^{|\mathbb{I}_{>}|\times
p}$.

As illustrated in Figure 1, there is a trend of overestimation to the left of
the overall mean and underestimation to the right. Therefore, we intend to
adjust this systematic bias through the imposition of constraints.
Specifically, predicted values corresponding to the set $\mathbb{I}_{<}$ will
be adjusted to have the same mean as the observed values within their group,
that is $\texttt{Avg}(\widehat{\mathbb{y}}_{<})=\overline{\mathbb{y}}_{<}$,
where Avg($\cdot$) denotes average, i.e.,
$\texttt{Avg}(\widehat{\mathbb{y}}_{<})=\frac{1}{|\mathbb{I}_{<}|}\sum_{i\in\mathbb{I}_{<}}\widehat{\mathbb{y}}_{i}$.
This constraint is applied to correct the overestimation bias associated with
lower response values. Similarly, for $\mathbb{I}_{>}$, we impose the
constraint $\texttt{Avg}(\widehat{\mathbb{y}}_{>})=\overline{\mathbb{y}}_{>}$
to correct the underestimation bias observed to the right of the overall mean.

Under the linear model, where the prediction for the $i$-th subject can be
calculated as $\widehat{y}_{i}=\mathbb{X}_{i}\beta$, the aforementioned
constraints can be expressed as linear equality constraints in matrix form.
Thus, we can formulate the lasso problem with an equality constraint:

$\displaystyle\operatornamewithlimits{argmin}_{\beta\in\mathbb{R}^{p}}\quad$
$\displaystyle\|\mathbb{y}-\mathbb{X}\boldsymbol{\beta}\|^{2}+\lambda\|\boldsymbol{\beta}\|_{1}$
subject to
$\displaystyle\begin{pmatrix}\bf{1}_{|\mathbb{I}_{<}|}^{\top}\mathbb{X}_{<}\\\
\bf{1}_{|\mathbb{I}_{>}|}^{\top}\mathbb{X}_{>}\end{pmatrix}\boldsymbol{\beta}=\begin{pmatrix}|\mathbb{I}_{<}|\cdot\overline{\mathbb{y}}_{<}\\\
|\mathbb{I}_{>}|\cdot\overline{\mathbb{y}}_{>}\end{pmatrix}$

The above form is the constrained lasso problem Gaines et al., (2018) with an
equality-only constraint. It can be solved a quadratic programming.

###  2.2 Constrained KRR

A similar approach can be extended to nonlinear models. A general
nonparametric regression problem can be formulated as
$\mathbb{y}_{i}=f(\mathbb{x}_{i})+\epsilon_{i}$ where the goal is to estimate
a function $\widehat{f}$ from a function space
$\mathcal{F}:\mathcal{X}\rightarrow\mathbb{R}$. One popular choice for the
function space $\mathcal{F}$ is the Reproducing Kernel Hilbert Space (RKHS),
denoted by $\mathcal{H}$. Given that $\mathcal{H}$ is infinite-dimensional,
regularization is essential for estimation. The optimization problem known as
kernel ridge regression (KRR) Schölkopf and Smola, (2002) can be formulated as
follows:

$\displaystyle\widehat{f}=\operatornamewithlimits{argmin}_{f\in\mathcal{H}}\sum_{i=1}^{n}\left(\mathbb{y}_{i}-f(\mathbb{x}_{i})\right)^{2}+\lambda\|f\|_{\mathcal{H}}^{2},$

where the first term is goodness-of-fit to data, and the second term is the
regularization term which controls the trade-off between the fit of the model
to the data and the smoothness of the function $f$. By the representer theorem
Wahba, (1990), the function $f$ can be expressed as a linear combination of
kernel functions, i.e.,
$f(\mathbb{x})=\sum_{i=1}^{n}\alpha_{i}\kappa(\mathbb{x},\mathbb{x}_{i})$. The
function $\kappa:\mathcal{X}\times\mathcal{X}\rightarrow\mathbb{R}$ is called
a kernel function if it is symmetric and positive definite. The most popular
choice of the kernel is Gaussian Radial Basis Function (RBF) Kernel, i.e.,
$\kappa(\mathbb{x}_{i},\mathbb{x}_{j})=\|\mathbb{x}_{i}-\mathbb{x}_{j}\|^{2}/\sigma^{2}$.

By imposing the two previously used equality constraints for adjusting the
mean of the predictive values, the KRR objective function can be expressed as
follows:

$\displaystyle\widehat{f}=$
$\displaystyle\operatornamewithlimits{argmin}_{f\in\mathcal{H}}$
$\displaystyle\sum_{i=1}^{n}\left(\mathbb{y}_{i}-f(\mathbb{x}_{i})\right)^{2}+\lambda\|f\|_{\mathcal{H}}^{2}$
subject to
$\displaystyle\sum_{i\in\mathbb{I}_{<}}\overline{\mathbb{y}}_{<}-f(\mathbb{x}_{i})=0$
(2)
$\displaystyle\sum_{i\in\mathbb{I}_{>}}\overline{\mathbb{y}}_{>}-f(\mathbb{x}_{i})=0$

One of the primary advantages of the kernel method is its ability to apply
nonlinear techniques within a linear algorithmic framework by utilizing the
kernel trick Hofmann et al., (2008). Consequently, unlike other nonlinear
regression methods, KRR benefits from the incorporation of the mean equality
constraints, which can be treated as linear constraints with respect to the
optimization variables. Specifically, by substituting the representation of
$f(\cdot)$ into (2), it can express in a dual form:

$\displaystyle\widehat{\boldsymbol{\alpha}}=$
$\displaystyle\operatornamewithlimits{argmin}_{\alpha\in\mathbb{R}^{n}}$
$\displaystyle\|\mathbb{y}-\mathbb{K}\boldsymbol{\alpha}\|_{2}^{2}+\lambda\boldsymbol{\alpha}^{\top}\mathbb{K}\boldsymbol{\alpha}$
subject to
$\displaystyle\begin{pmatrix}\bf{1}_{|\mathbb{I}_{<}|}^{\top}\mathbb{K}_{<}\\\
\bf{1}_{|\mathbb{I}_{>}|}^{\top}\mathbb{K}_{>}\end{pmatrix}\boldsymbol{\alpha}=\begin{pmatrix}|\mathbb{I}_{<}|\cdot\overline{\mathbb{y}}_{<}\\\
|\mathbb{I}_{>}|\cdot\overline{\mathbb{y}}_{>}\end{pmatrix},$

where $\mathbb{K}\in\mathbb{R}^{n\times n}$ denotes the kernel gram matrix of
the complete data with the $(i,j)$ element is
$\kappa(\mathbb{x}_{i},\mathbb{x}_{j})$, and
$\mathbb{K}_{<}\in\mathbb{R}^{|\mathbb{I}_{<}|\times n}$,
$\mathbb{K}_{>}\in\mathbb{R}^{|\mathbb{I}_{>}|\times n}$ denote the kernel
matrices corresponding to datasets $\mathbb{I}_{<}$, $\mathbb{I}_{>}$,
respectively. Thus, for example the $i$-the row of $\mathbb{K}_{<}$ is
$\mathbb{K}_{{<}(i,\cdot)}=\left(\kappa(\mathbb{x}_{{<}_{i}},\mathbb{x}_{1}),\cdots,\kappa(\mathbb{x}_{{<}_{i}},\mathbb{x}_{n})\right)$,
where ${\mathbb{I}_{<}}_{i}$ denotes the $i$-th subject in the set
$\mathbb{I}_{<}$. Because the main objective function is quadratic and the
constraint is linear in $\boldsymbol{\alpha}$, the above optimization problem
can be readily solved using quadratic programming.

By formulating the above problem in the Lagrangian form:

$\displaystyle\mathcal{L}(\boldsymbol{\alpha},{\boldsymbol{\rho}})=\frac{1}{2}\boldsymbol{\alpha}^{\top}\mathbb{K}\left(\mathbb{K}+\frac{\lambda}{2}{\bf{I}}\right)\boldsymbol{\alpha}-\mathbb{y}^{\top}\mathbb{K}\boldsymbol{\alpha}+{\boldsymbol{\rho}}^{\top}\left(\underbrace{\begin{pmatrix}\bf{1}_{|\mathbb{I}_{<}|}^{\top}\mathbb{K}_{<}\\\
\bf{1}_{|\mathbb{I}_{>}|}^{\top}\mathbb{K}_{>}\end{pmatrix}}_{\mathbf{B}}\boldsymbol{\alpha}-\underbrace{\begin{pmatrix}|\mathbb{I}_{<}|\cdot\overline{\mathbb{y}}_{<}\\\
|\mathbb{I}_{>}|\cdot\overline{\mathbb{y}}_{>}\end{pmatrix}}_{\mathbf{d}}\right)$

we can derive the following Algorithm 1.

Algorithm 1 Algorithm for constrained KRR

1:Update $\boldsymbol{\alpha}^{(t)}$ by finding the stationary point of
$\mathcal{L}(\boldsymbol{\alpha},\boldsymbol{\rho}^{(t-1)}$)
$\displaystyle\nabla_{\boldsymbol{\alpha}}\mathcal{L}(\boldsymbol{\alpha},{\boldsymbol{\rho}^{(t-1)}})=\mathbb{K}\left(\mathbb{K}+\frac{\lambda}{2}{\bf{I}}\right)\boldsymbol{\alpha}-\mathbb{K}\mathbb{y}+\mathbf{B}^{\top}{\boldsymbol{\rho}^{(t-1)}}\equiv
0$

2:Update $\boldsymbol{\rho}^{(t)}$
$\displaystyle{\boldsymbol{\rho}^{(t)}}\leftarrow{\boldsymbol{\rho}^{(t-1)}}+s^{(t)}(\mathbf{B}\boldsymbol{\alpha}^{(t)}-\mathbf{d}),$
where $s^{(t)}$ is a step size, which can be chosen by line search.

3:Iterate updating $\boldsymbol{\alpha}$, and $\boldsymbol{\rho}$ until
converge.

##  3 Simulation

In this section, we benchmark our proposed methods (KRRConst, and LassoConst)
against their unconstrained counterparts and other ML methods including
XGBoost (Chen and Guestrin,, 2016) and Random Forest (RF; Breiman, (2001)). We
consider both linear and nonlinear relationships between predictor
$\mathbb{X}$ and response $\mathbb{y}$. In the linear case, we first generate
$\mathbb{X}$ from the multivariate normal distribution with dimension $p=10$
and then generate the response
$\mathbb{y}=\mathbb{X}\boldsymbol{\beta}+\epsilon$, where $\epsilon\sim
N(0,\sigma_{e}^{2})$. In the nonlinear case, we generate
$\mathbb{y}=f(\mathbb{X})+\epsilon$, and for the form of $f(\mathbb{X})$, we
use two of Friedman’s test functions, which are commonly used for evaluating
machine learning algorithms (Friedman,, 1991):

  1. 1.

$f(\mathbb{X})=0.1\exp(4\mathbb{X}_{1})+4/(1+\exp(-20(\mathbb{X}_{2}-0.5))+3\mathbb{X}_{3}+2\mathbb{X}_{4}+\mathbb{X}_{5}+0\cdot\sum_{p=6}^{10}\mathbb{X}_{p}$

  2. 2.

$f(\mathbb{X})=10\sin(\pi\mathbb{X}_{1}\mathbb{X}_{2})+20(\mathbb{X}_{3}-0.5)^{2}+10\mathbb{X}_{4}+5\mathbb{X}_{5}$.

We fitted the model using only the training set and evaluated its performance
on both the training set and the testing set, which was unseen during the
training stage. For all settings, we use an equal number of samples for the
training and testing sets, with each set containing $n=100$.

To assess whether the methods suffer from or resolve the SBMR in predictions,
we evaluate the bias by comparing the regression slope—obtained by regressing
predicted values on observed values—to the reference line set at 1 (i.e.,
$\widehat{\mathbb{y}}=\mathbb{y}$). In addition, to determine if there is a
linear relationship between the residuals and the predicted values, we
calculated the correlation between them. To measure the overall prediction
accuracy, we also compare the Root Mean Square Error (RMSE). To further
investigate the presence of systematic bias around the tail areas, we
calculated the average bias for observations falling below the lower quartile
(i.e., $<Q_{1}$) and above the upper quartile (i.e., $>Q_{3}$). The results
averaged over 100 replications are provided in Table 1.

Our proposed methods perform well, exhibiting only minimal bias in slope
across all settings and in both the training and testing datasets, with biases
consistently less than 0.1, except for KRRConst in the testing set under the
linear setting.

However, all conventional methods exhibit a biased slope, high linear
correlation between predicted values and residuals, and systematic bias in the
tail areas. Specifically, models tend to produce a positive bias for responses
under the lower quantiles, indicating that the predicted values are generally
larger than the observed values. Conversely, responses above the upper
quantiles typically result in negative biases, suggesting an underestimation
of the observed values. In addition, both XGBoost and random forest exhibit
overfitting, as evidenced by their superior performance on the training set,
particularly in terms of RMSE, compared to their poorer performance on the
testing set. The issue of overfitting with these ML methods is more pronounced
in the linear case, as evidenced by a bias in slope greater than 0.65. This
implies that the predicted values tend to be closer to the overall mean,
indicating that the SBMR occurs.

It should be noted that our proposed method does not achieve optimal results
in terms of RMSE. This is because our objective function does not globally
minimize the squared error loss, as it incorporates constraints designed to
adjust the bias. Consequently, this results in an increase in variance, which
is a trade-off we have to accept to eliminate systematic bias.

|  |  | KRRConst | LassoConst | KRR | Lasso | XGBoost | RF
---|---|---|---|---|---|---|---|---
Linear | Train | Bias (Slope) | 0.0204 | 0.0243 | 0.2962 | 0.5068 | 0.0187 | 0.3227
$\operatornamewithlimits{Cor}(\widehat{\mathbb{y}},\hat{\epsilon})$ | -0.0332 | -0.0357 | -0.5811 | -0.8034 | -0.6516 | -0.9187
RMSE | 0.5765 | 0.6137 | 0.3985 | 0.4945 | 0.0228 | 0.2772
Bias ($<Q_{1}$) | 0.0233 | 0.0260 | 0.2998 | 0.5076 | 0.0201 | 0.3186
Bias ($>Q_{3}$) | -0.0202 | -0.0221 | -0.2959 | -0.5039 | -0.0157 | -0.3206
Test | Bias (Slope) | 0.1178 | 0.0835 | 0.3907 | 0.5599 | 0.6815 | 0.7547
$\operatornamewithlimits{Cor}(\widehat{\mathbb{y}},\hat{\epsilon})$ | -0.1561 | -0.1061 | -0.5825 | -0.8021 | -0.7706 | -0.9281
RMSE | 0.6351 | 0.6820 | 0.5403 | 0.5620 | 0.7142 | 0.6576
Bias ($<Q_{1}$) | 0.1014 | 0.0748 | 0.3852 | 0.5592 | 0.6965 | 0.7584
Bias ($>Q_{3}$) | -0.1349 | -0.0905 | -0.4079 | -0.5774 | -0.6818 | -0.7710
Nonlinear 1 | Train | Bias (Slope) | 0.0292 | 0.0463 | 0.3069 | 0.3854 | 0.0112 | 0.2125
$\operatornamewithlimits{Cor}(\widehat{\mathbb{y}},\hat{\epsilon})$ | -0.0592 | -0.0794 | -0.6380 | -0.7388 | -0.4957 | -0.8560
RMSE | 1.3677 | 1.6200 | 1.2824 | 1.3899 | 0.0635 | 0.6638
Bias ($<Q_{1}$) | 0.0122 | 0.0763 | 0.9644 | 1.2281 | 0.0100 | 0.6812
Bias ($>Q_{3}$) | -0.1195 | -0.1781 | -1.0680 | -1.3374 | -0.0595 | -0.7311
Test | Bias (Slope) | 0.0781 | 0.0906 | 0.3344 | 0.4015 | 0.2943 | 0.4939
$\operatornamewithlimits{Cor}(\widehat{\mathbb{y}},\hat{\epsilon})$ | -0.1376 | -0.1413 | -0.6289 | -0.7275 | -0.5324 | -0.8632
RMSE | 1.5603 | 1.8139 | 1.4328 | 1.4856 | 1.4950 | 1.5449
Bias ($<Q_{1}$) | 0.1544 | 0.1868 | 1.0477 | 1.2675 | 0.7899 | 1.5844
Bias ($>Q_{3}$) | -0.3167 | -0.3578 | -1.1998 | -1.4338 | -1.1879 | -1.7578
Nonlinear 2 | Train | Bias (Slope) | 0.0065 | 0.0123 | 0.3442 | 0.4350 | 0.0124 | 0.2497
$\operatornamewithlimits{Cor}(\widehat{\mathbb{y}},\hat{\epsilon})$ | -0.0149 | -0.0244 | -0.6590 | -0.7790 | -0.4076 | -0.9011
RMSE | 2.9920 | 3.9202 | 2.5692 | 2.7455 | 0.1558 | 1.3704
Bias ($<Q_{1}$) | -0.0622 | -0.0512 | 2.0947 | 2.6742 | 0.0153 | 1.5576
Bias ($>Q_{3}$) | -0.0987 | -0.1143 | -2.2178 | -2.7873 | -0.1281 | -1.5597
Test | Bias (Slope) | 0.0349 | 0.0625 | 0.3590 | 0.4495 | 0.2882 | 0.5220
$\operatornamewithlimits{Cor}(\widehat{\mathbb{y}},\hat{\epsilon})$ | -0.0642 | -0.0932 | -0.6644 | -0.7789 | -0.5669 | -0.9125
RMSE | 3.1257 | 4.1235 | 2.7400 | 2.9182 | 2.5776 | 2.9002
Bias ($<Q_{1}$) | -0.0102 | 0.1719 | 2.1599 | 2.7855 | 1.6632 | 3.3347
Bias ($>Q_{3}$) | -0.3796 | -0.5188 | -2.4279 | -2.9739 | -2.0054 | -3.3655
Table 1: Results of the experiments described in Section 3. The best-ranked
method for each measure is highlighted in bold.

##  4 Real-life example of MTM (Estimating Brain Age)

Regression to the mean has been commonly reported in neuroimaging studies,
particularly in the context of estimating brain age to quantify an
individual’s brain condition by regressing chronological age on neuroimaging
measurements. In such cases, older individuals tend to be estimated with a
younger brain age, while younger individuals tend to be estimated with an
older brain age.

We applied the proposed methods to two neuroimaging studies: The Lifespan
Human Connectome Project Aging Study (HCP-A,
https://www.humanconnectome.org/study/hcp-lifespan-aging), and UK Biobank
(UKBB, https://www.ukbiobank.ac.uk/) data. In both the HCP-A and UKBB
datasets, we utilized fractional anisotropy (FA), a measure of white matter
integrity, measured across various brain regions as the neuroimaging
predictors. The dimensions of the predictors are $p=64$ for HCP-A and $p=39$
for UKBB, respectively. We use chronological age as the response variable. The
age range for the HCP-A dataset is 36-90 years, and for the UKBB dataset, it
is 45-82 years. The number of samples for which neuroimaging measurements are
available is 662 for HCP-A and 36,856 for UKBB.

We trained the model on randomly sampled subsets from the entire dataset—300
samples for HCP-A and 500 for UKBB. We then evaluated model performance using
the unseen testing set, which consists of the remaining data for HCP-A and
another 500 randomly selected samples from the remaining data for UKBB. We
compare the performance of our method with other ML methods. The results from
100 replications are summarized in Table 2.

|  | KRRConst | LassoConst | KRR | Lasso | XGBoost | RF
---|---|---|---|---|---|---|---
HCP-A | Bias (Slope) | 0.0769 | 0.1057 | 0.4995 | 0.4194 | 0.4538 | 0.5046
$\operatornamewithlimits{Cor}(\widehat{\mathbb{y}},\hat{\epsilon})$ | -0.0874 | -0.1525 | -0.7553 | -0.6559 | -0.6086 | -0.7449
RMSE | 13.0182 | 10.0511 | 9.4118 | 9.0879 | 10.6284 | 9.6430
Bias ($<Q_{1}$) | 3.8303 | 3.8248 | 10.4100 | 8.9860 | 9.1444 | 10.3606
| Bias ($>Q_{3}$) | -0.2421 | -0.8128 | -8.5090 | -7.0122 | -8.2300 | -8.6196
UKBB | Bias (Slope) | 0.0659 | 0.1450 | 0.6705 | 0.7098 | 0.6298 | 0.6725
$\operatornamewithlimits{Cor}(\widehat{\mathbb{y}},\hat{\epsilon})$ | -0.0501 | -0.1147 | -0.8354 | -0.8679 | -0.7292 | -0.8390
RMSE | 10.8041 | 9.8240 | 6.1308 | 6.2405 | 6.6023 | 6.1208
Bias ($<Q_{1}$) | 1.5131 | 1.9598 | 7.2653 | 7.6575 | 6.3815 | 7.1399
| Bias ($>Q_{3}$) | 0.1360 | -0.8688 | -6.3362 | -6.7404 | -6.3815 | -6.4899
Table 2: Results of the real data application. The best-ranked method for each
measure is highlighted in bold.

Overall, the constrained KRR outperforms other methods, with the only
exceptions being lower tail bias in the HCP-A data, where the proposed
constrained lasso performs the best with only a small margin, and RMSE for
both datasets, where conventional methods perform well. The small bias in
slope (all less than 0.1) of constrained KRR indicates that the predicted
values are closely aligned with the observed values. Additionally, the
correlation between residuals and response values is significantly weaker for
the proposed methods compared to other competing methods. This is also
illustrated in Figure 2, which displays a scatter plot of one sample result
from 100 replications. This lack of a linear trend indicates that the proposed
method removes systematic biases, leading to more accurate and reliable
predictions across the range of data. This removal of systematic bias is
further evidenced by the small biases observed in the tail areas, below the
first quartile ($<Q_{1}$) and above the third quartile ($>Q_{3}$).

As noted, the proposed method’s reduction in bias may lead to increased
variance in model predictions, potentially raising the RMSE. Consequently, our
method is not optimal in terms of RMSE, illustrating a trade-off between bias
and variance. However, the difference in RMSE between the proposed methods and
other ML methods is not critical, considering the substantial improvement in
tail bias.

Figure 2: This Scatter plot, one result from 100 replications, shows the
observed response $\mathbb{Y}$ ($x$-axis) and the residuals ($y$-axis) from
various methods on the testing set. The reference dashed line represents
$y=0$.

##  5 Conclusion

The bias-variance trade-off is fundamental in machine learning research. In
order to reduce the mean squared error, many machine learning regression
models tolerate biases of predicted outcomes. When the prediction is
systematically biased, as in the scenario “Machine learning towards the mean”,
the foundation of machine learning regression models is undermined. We
proposed a simple remedy to the systematic bias of machine learning regression
models. We developed new algorithms for two regression shrinkage models
handling high-throughput predictors (i.e., $p>n$) by imposing two-sided
constraints to correct the systematic bias. In both simulation analysis and
real-world data applications our methods yield unbiased prediction.

Although our correction methods can avoid systematic bias, they are not exempt
from the bias-variance trade-off. The mean squared error of corrected methods
is generally higher than the classic machine learning regression models. In
future research, we will integrate the proposed corrections into commonly used
machine learning regression models such as random forest, XGBoost, and others
and explore adaptive solutions to address the bias-variance trade-off.

## References

  * Belitz and Stackelberg, (2021) Belitz, K. and Stackelberg, P. (2021).  Evaluation of six methods for correcting bias in estimates from ensemble tree machine learning regression models.  Environmental Modelling & Software, 139:105006.
  * Breiman, (2001) Breiman, L. (2001).  Random forests.  Machine Learning, 45(1):5–32.
  * Butler et al., (2021) Butler, E. R., Chen, A., Ramadan, R., Le, T. T., Ruparel, K., Moore, T. M., Satterthwaite, T. D., Zhang, F., Shou, H., Gur, R. C., Nichols, T. E., and Shinohara, R. T. (2021).  Pitfalls in brain age analyses.  Human Brain Mapping, 42(13):4092–4101.
  * Chen and Guestrin, (2016) Chen, T. and Guestrin, C. (2016).  Xgboost: A scalable tree boosting system.  KDD ’16, page 785–794, New York, NY, USA. Association for Computing Machinery.
  * Friedman, (1991) Friedman, J. H. (1991).  Multivariate Adaptive Regression Splines.  The Annals of Statistics, 19(1):1 – 67.
  * Gaines et al., (2018) Gaines, B. R., Kim, J., and Zhou, H. (2018).  Algorithms for fitting the constrained lasso.  Journal of Computational and Graphical Statistics, 27(4):861–871.
  * Hastie et al., (2009) Hastie, T., Tibshirani, R., and Friedman, J. (2009).  The Elements of Statistical Learning.  Springer New York, NY, 2 edition.
  * Hofmann et al., (2008) Hofmann, T., Schölkopf, B., and Smola, A. J. (2008).  Kernel methods in machine learning.  The Annals of Statistics, 36(3):1171 – 1220.
  * Schölkopf and Smola, (2002) Schölkopf, B. and Smola, A. J. (2002).  Learning with Kernels:Support vector machines, Regularization, Optimization, and Beyond.  MIT Press, Cambridge, MA, USA.
  * Smith et al., (2019) Smith, S. M., Vidaurre, D., Alfaro-Almagro, F., Nichols, T. E., and Miller, K. L. (2019).  Estimation of brain age delta from brain imaging.  NeuroImage, 200:528–539.
  * Wahba, (1990) Wahba, G. (1990).  Spline Models for Observational Data.  Society for Industrial and Applied Mathematics.
  * Zhang and Lu, (2012) Zhang, G. and Lu, Y. (2012).  Bias-corrected random forests in regression.  Journal of Applied Statistics, 39(1):151–160.

##  Appendix A Proof of Proposition 1

###### Proof.

Let $\widehat{y}_{i}$ denote an unbiased prediction $y_{i}$. Recall that
$R^{2}=1-\dfrac{\sum_{i=1}^{n}(y_{i}-\widehat{y}_{i})^{2}}{\sum_{i=1}^{n}(y_{i}-\overline{y})^{2}}$,
and
$\sum_{i=1}^{n}(y_{i}-\widehat{y}_{i})^{2}=\sum_{i=1}^{n}(y_{i}-\overline{y})^{2}\times(1-R^{2})$.
Without loss of generality, we have
$\mathbb{E}\left[\sum_{i=1}^{n}(y_{i}-\widehat{y}_{i})^{2}\right]=(n-1)\sigma^{2}(1-R^{2})$.

Suppose $\widetilde{y}_{i}=\widehat{y}_{i}-c\widehat{y}_{i}$ is the linear
trend of the systematic bias, then following holds.

$\displaystyle\mathbb{E}\left[\sum_{i=1}^{n}(y_{i}-\widetilde{y}_{i})^{2}\right]$
$\displaystyle=\mathbb{E}\left[\sum_{i=1}^{n}\left(y_{i}-\mathbb{E}(\widetilde{y}_{i})+\mathbb{E}(\widetilde{y}_{i})-\widetilde{y}_{i}\right)^{2}\right]$
$\displaystyle=\mathbb{E}\left[\sum_{i=1}^{n}(y_{i}-\mathbb{E}(\widetilde{y}_{i}))^{2}+\sum_{i=1}^{n}(\mathbb{E}(\widetilde{y}_{i})-\widetilde{y}_{i})^{2}+\sum_{i=1}^{n}2(y_{i}-\mathbb{E}(\widetilde{y}_{i}))(\mathbb{E}(\widetilde{y}_{i})-\widetilde{y}_{i})\right]$
$\displaystyle=\mathbb{E}\sum_{i=1}^{n}(cy_{i})^{2}+\mathbb{E}\sum_{i=1}^{n}[(1-c)(y_{i}-\widehat{y}_{i})]^{2}+2\mathbb{E}\sum_{i=1}^{n}(y_{i}-\mathbb{E}(\widetilde{y}_{i}))(\mathbb{E}(\widetilde{y}_{i})-\widetilde{y}_{i})$
$\displaystyle=(n-1)c^{2}\sigma^{2}+(n-1)(1-c)^{2}\sigma^{2}(1-R^{2}).$

It suffices to verify the condition where
$\mathbb{E}\left[\sum_{i=1}^{n}(y_{i}-\widehat{y}_{i})^{2}\right]>\mathbb{E}\left[\sum_{i=1}^{n}(y_{i}-\widetilde{y}_{i})^{2}\right]$:

$\displaystyle(1-R^{2})$ $\displaystyle>c^{2}+(1-c)^{2}(1-R^{2})$
$\displaystyle(1-R^{2})$ $\displaystyle>1-2c+2c^{2}-R^{2}+2cR^{2}-c^{2}R^{2}$
$\displaystyle 2c-2c^{2}-2cR^{2}+c^{2}R^{2}$ $\displaystyle>0$ $\displaystyle
1-c-R^{2}+cR^{2}/2$
$\displaystyle>0\rightarrow\left(1-\frac{R^{2}}{2}\right)c<1-R^{2}$

Thus, we have
$\mathbb{E}\left[\sum_{i=1}^{n}(y_{i}-\widehat{y}_{i})^{2}\right]>\mathbb{E}\left[\sum_{i=1}^{n}(y_{i}-\widetilde{y}_{i})^{2}\right]$,
when $c<\frac{1-R^{2}}{1-R^{2}/2}$. ∎

Therefore, machine learning regression models with the objective function to
minimize mean squared error naturally incorporate the linear bias
$\widetilde{y}_{i}=\widehat{y}_{i}-c\widehat{y}_{i}$ in their prediction.