# GGHead: Fast and Generalizable 3D Gaussian Heads

Tobias Kirschstein  Technical University of MunichGermany ,  Simon Giebenhain
Technical University of MunichGermany ,  Jiapeng Tang  Technical University of
MunichGermany ,  Markos Georgopoulos  Independent ResearcherSwitzerland and
Matthias Nießner  Technical University of MunichGermany

###### Abstract.

Learning 3D head priors from large 2D image collections is an important step
towards high-quality 3D-aware human modeling. A core requirement is an
efficient architecture that scales well to large-scale datasets and large
image resolutions. Unfortunately, existing 3D GANs struggle to scale to
generate samples at high resolutions due to their relatively slow train and
render speeds, and typically have to rely on 2D superresolution networks at
the expense of global 3D consistency. To address these challenges, we propose
Generative Gaussian Heads (GGHead), which adopts the recent 3D Gaussian
Splatting representation within a 3D GAN framework. To generate a 3D
representation, we employ a powerful 2D CNN generator to predict Gaussian
attributes in the UV space of a template head mesh. This way, GGHead exploits
the regularity of the template’s UV layout, substantially facilitating the
challenging task of predicting an unstructured set of 3D Gaussians. We further
improve the geometric fidelity of the generated 3D representations with a
novel total variation loss on rendered UV coordinates. Intuitively, this
regularization encourages that neighboring rendered pixels should stem from
neighboring Gaussians in the template’s UV space. Taken together, our pipeline
can efficiently generate 3D heads trained only from single-view 2D image
observations. Our proposed framework matches the quality of existing 3D head
GANs on FFHQ while being both substantially faster and fully 3D consistent. As
a result, we demonstrate real-time generation and rendering of high-quality
3D-consistent heads at $1024^{2}$ resolution for the first time.

Figure 1. GGHead: Our method can generate diverse 3D head representations
based on 3D Gaussian Splatting. The sampled persons exhibit detailed geometry
and appearance. Both generation and rendering can be conducted in real-time at
full image resolution without the need for 2D super-resolution networks.
††footnotetext: Project Website: https://tobias-kirschstein.github.io/gghead/

##  1\. Introduction

A high-quality 3D generative model that can sample human heads with diverse
geometry and appearance has to satisfy two important constraints: A high
rendering resolution and strict 3D consistency of generated imagery. Once met,
these properties enable exciting use cases such as single-image-to-3D
reconstruction, 3D content creation, or 3D-consistent editing of images.

Unfortunately, obtaining such a high-quality 3D prior over human appearances
is extremely challenging. Not only does such a model have to explain all
diverse human appearances, but there is also a serious lack of large-scale 3D
human datasets due to expensive data capture procedures, rendering any attempt
to learn 3D human priors from 3D datasets alone impractical. Hence, methods
need to be able to learn 3D appearance and geometry from large 2D image
collections, requiring them to render during training. As a result, to meet
the criteria for quality and 3D-consistency, the employed rendering procedure
has to be scalable, both in terms of speed and rendering resolution.

A well-established approach are 3D Generative Adversarial Networks (3D GANs)
that can learn unconditional 3D generation models from 2D images with
corresponding camera poses. These methods have shown remarkable 3D-aware image
synthesis capabilities, rapidly closing the quality gap to methods that
directly generate 2D face images without view control. Under the hood, 3D GANs
combine a differentiable rendering pipeline with a view-conditioned
discriminator architecture to supervise 3D generation solely with 2D images.
To obtain high-quality 3D samples with this approach, a high rendering
resolution is crucial. A commonly chosen and successful differentiable
rendering pipeline is Neural Radiance Fields (NeRFs) (Mildenhall et al.,
2021). However, so far, scaling the rendering resolution of NeRFs has been
bottle-necked by the expensive ray marching procedure necessary to synthesize
the final image. Several works have attempted to alleviate this limitation:
EG3D (Chan et al., 2022) renders images at lower resolutions and then
upsamples them using 2D superresolution, sacrificing 3D consistency in the
process. GRAM (Deng et al., 2022) and GRAM-HD (Xiang et al., 2023) speed up
ray marching by simplifying the 3D representation to a set of 2D manifolds,
allowing cheap ray intersection tests. However, the simplistic 3D
representation worsens rendering quality especially for side-views. Mimic3D
(Chen et al., 2023) generates a pseudo-GT dataset from EG3D to supervise the
3D representation directly with patch-based rendering. Due to the inherent
view-inconsistencies in the generated multi-view images from EG3D, only a
perceptual loss can be used for supervision to gloss over the misalignment.
Most recently, Trevithick et al. (2024) investigate how ray marching can be
sped up with proposal networks, allowing training at $512^{2}$ resolution.
Nevertheless, rendering speeds during both training and inference remain a
significant challenge. For instance, currently, none of the existing methods
are capable of natively generating and rendering images at full $512^{2}$
resolution in real-time. Since the creation of a 3D GAN involves several
million render passes, this imposes a considerable computational burden on
training.

To address these scalability challenges, we propose to replace the commonly
used NeRF 3D representation with the recent 3D Gaussian Splatting (3DGS)
pipeline (Kerbl et al., 2023). Based on rasterization instead of expensive ray
marching, 3DGS has shown remarkable rendering speeds harboring great potential
for more scalable 3D GAN pipelines. To generate a set of 3D Gaussian
primitives, we exploit the UV space of a template mesh, allowing us to
leverage powerful 2D CNN architectures such as StyleGAN2 (Karras et al.,
2020b). From the predicted 2D maps for position, scale, rotation, color, and
opacity, we sample 3D Gaussians and place them in the 3D scene relative to the
template. Finally, the 3D representation is rendered and supervised with a
view-point aware discriminator. The purpose of the template mesh is two-fold:
(i) Simplify the generation process by predicting structured 2D UV maps
instead of unstructured 3D primitives, (ii) Regularize the predicted 3D
representation. The latter is necessary because 3D Gaussians, as opposed to
radiance fields, are highly sensitive to just a few parameters. Therefore, bad
gradients during the early stages of adversarial training can cause blow-up.
To counteract this, we propose a novel UV total variation loss that makes the
3D representation more well-behaved by enforcing continuity of rendered pixels
with respect to their UV coordinates. Our final model can learn 3D head
geometry and appearance priors from 2D data collections alone. Furthermore,
GGHead generates and renders photorealistic 3D heads in real-time thanks to
its efficient 3D representation, allowing it to scale to larger resolutions.

In summary, our contributions are as follows:

  * •

We propose GGHead, a 3D GAN for human heads based on 3D Gaussian Splatting
that is trained only from 2D image observations.

  * •

We parameterize 3D Gaussian Heads as UV maps that can be generated with
efficient 2D CNNs in conjunction with a novel UV total variation
regularization that improves geometric fidelity.

  * •

Our approach is highly scalable, facilitating training at 1k resolutions while
enabling both real-time sample generation and head rendering.

##  2\. Related Work

In this section, we discuss prior related work in the field of generative
modeling of 3D heads. We begin by reviewing existing 3D GANs and their
application to 3D heads. We then discuss 3DGS-based models and their extension
to generalizable representations.

###  2.1. 3D Generative adversarial models

Following the success of their 2D counterpart, 3D GANs employ adversarial
training to obtain a generative model of 3D representations from 2D images.
Earlier works utilise implicit 3D representations and directly render pixels
(Chan et al., 2021; Schwarz et al., 2020) or features that are decoded with a
CNN-based neural renderer (Niemeyer and Geiger, 2021; Xue et al., 2022). The
latter approach has also been adopted by efforts that aim to repurpose
successful 2D GAN architectures (e.g., (Karras et al., 2019a, 2020b, 2020a))
and lift the representation to 3D using implicit neural representations (Chan
et al., 2022; Or-El et al., 2022; Gu et al., 2021). Due to their success,
these approaches have been extended to introduce different modalities of
control, e.g., semantics (Sun et al., 2022).

To obtain higher resolutions, some methods such as EG3D (Chan et al., 2022)
apply a super-resolution network to the rendered features. While efficient,
this approach can introduce unwanted artefacts in the form of 3D
inconsistencies as well as low-resolution geometry. To alleviate this issue,
several works focus on rendering in the pixel-space while maintaining fidelity
and efficiency. Among them, Skorokhodov et al. (2022) propose to render
patches instead of full images to reduce the computational cost. While this
allows supervising a generative 3D representation without super-resolution, it
impacts global consistency due to the discriminator’s limited field of view.
Mimic3D (Chen et al., 2023) proposes an imitation strategy that forces the 3D
representations to mimic the result of the 2D super-resolution branch.
Similarly, GRAM-HD (Xiang et al., 2023) utilizes radiance manifolds and
performs super-resolution on 2D manifolds instead of the rendered images. More
recently, Trevithick et al. (2024) introduce a learning-based sampling
strategy to effectively render every pixel and to model high-resolution
geometry. Orthogonally to these work, our method aims to improve rendering
speed by adopting 3D Gaussian Splatting as a 3D representation instead of
NeRFs.

###  2.2. Generalizable 3D Gaussian Splatting

To mitigate the computational overhead of volumetric rendering imposed by
implicit 3D representations, Kerbl et al. (2023) introduced the seminal work
of 3DGS. In this work, a set of Gaussians is optimized from multi-view images
using volume splatting (Zwicker et al., 2001). The means of the Gaussians are
initialized using points from SfM (Schonberger and Frahm, 2016) while
subsequent efforts try to break this dependency (Fu et al., 2023). Inspired by
its efficient rendering, a plethora of methods have been proposed to extend
3DGS for real-time human avatar modeling utilising 3DMMs (Qian et al., 2023b,
a; Xu et al., 2023; Zielonka et al., 2023; Jiang et al., 2024; Svitov et al.,
2024; Hu and Liu, 2023; Moreau et al., 2023; Saito et al., 2023). Although
these efforts focus on modeling a single scene, 3DGS has been extended to
generalizable methods for scene reconstruction (Xu et al., 2024; Zou et al.,
2023), as well as text-driven synthesis using Score Distillation Sampling (Yi
et al., 2023; Tang et al., 2023; Li et al., 2023; Ren et al., 2023; Ling et
al., 2023). Another line of work investigates the use of 3DGS for
generalizable few-shot 3D reconstruction (Wewer et al., 2024; Charatan et al.,
2023). However, none of these works focus on learning an unconditional 3D
generative model (e.g., a GAN) directly on 3DGS-based representations. Closer
to our work, Gaussian Shell Maps (GSM) (Abdal et al., 2023) introduces a
Gaussian-based 3D GAN for human bodies by relying on shellmaps of fixed
Gaussians. In contrast, we propose a 3D GAN for human heads without the use of
shell maps or an animatable mesh prior, and instead with freely moving
Gaussians and a novel regularization scheme that exploits the UV space of a
template mesh.

##  3\. Method

Figure 2. Method Overview: We adopt 3D Gaussian Splatting inside a 3D GAN
formulation to learn a distribution of 3D heads solely from 2D images. To
build the 3D Gaussian representation, we exploit the UV space of a template
mesh. The CNN generator then takes a normally distributed latent code $z$ and
predicts 2D maps for each Gaussian attribute. To obtain the actual 3D Gaussian
primitives for rasterization, we sample the UV maps uniformly and place the
primitives relative to the template mesh. The generated 3D Gaussian
representation is then rasterized and supervised by the discriminator. To
increase training stability, especially during early stages of adversarial
training, we regularize the predicted position offsets, scales and opacities
via $\mathcal{L}_{reg}^{pos}$, $\mathcal{L}_{reg}^{scale}$ and
$\mathcal{L}_{reg}^{opac}$. Furthermore, we propose a novel UV total variation
loss $\mathcal{L}_{uv}$ to improve the geometric fidelity of generated 3D
heads by enforcing UV renderings to be smooth.

Figure 2 shows an overview of our proposed pipeline. We pursue a 3D GAN
approach, consisting of a generator that predicts 3D representations, a
differentiable rasterizer, and a discriminator that supervises the image
formation process. To make the setup scalable, we employ 3D Gaussian Splatting
as an underlying 3D representation. In the following, we will discuss how we
use a template mesh in conjunction with a 2D CNN to predict 3D Gaussian heads
(section 3.2), and how we regularize the 3D representation to ensure stable
training in an adversarial setting (sections 3.3 and 3.4).

###  3.1. 3D Gaussian Splatting (3DGS)

3D Gaussian Splatting (Kerbl et al., 2023) is a point-based scene
representation that assigns each point five different attributes: The Gaussian
center $\mu\in\mathbb{R}^{3}$, scale $\mathbf{s}\in\mathbb{R}^{3}$, rotation
parameterized as quaternion $\mathbf{q}\in\mathbb{R}^{4}$, color
$\mathbf{c}\in\mathbb{R}^{3}$, and opacity $\sigma\in\mathbb{R}$.

(1)
$\displaystyle\mathcal{G}^{i}=\\{\mathbf{\mu},\mathbf{s},\mathbf{q},\mathbf{c},\sigma\\}$

The resulting annotated point cloud can then be efficiently rendered into an
image $I$ using the differentiable tile-based rasterizer $\mathcal{R}$ and
camera parameters $\pi$:

(2) $\displaystyle I=\mathcal{R}(\mathcal{G},\pi)$

In the following, we will discuss how one can efficiently generate a 3D
Gaussian pointcloud $\mathcal{G}$ using 2D CNN architectures.

###  3.2. Template-based 3D Gaussian Generation

3D Gaussians, being a point-based representation, are inherently unstructured.
This poses a significant challenge for generation tasks due to issues such as
order ambiguity or large regions of empty space. We therefore follow Gaussian
Shell Maps (Abdal et al., 2023) and Relightable Gaussian Codec Avatars (Saito
et al., 2023), and rig the 3D Gaussian representation to a template mesh with
corresponding UV layout. This enables the use of an efficient 2D CNN backbone
$\mathcal{B}$ to predict the Gaussian representation $\mathcal{G}$,
considerably simplifying the generator’s task. Formally, we generate one UV
map $M$ for each Gaussian attribute:

(3) $\displaystyle M_{\star}=\mathcal{B}(z,\pi)$

where
$\star\in\\{\textsc{position},\textsc{scale},\textsc{rotation},\textsc{color},\textsc{opacity}\\}$
and $\mathcal{B}$ is a StyleGAN2 (Karras et al., 2020b) generator, mapping a
normally distributed latent code $z\in\mathbb{R}^{512}$ to UV maps
$M\in\mathbb{R}^{256\times 256\times 14}$.

For predicted positions, we add another learnable CNN layer $\mathcal{Z}$ with
weights initialized as zeros following Zhang et al. (2023):

(4) $\displaystyle M_{position}\leftarrow\mathcal{Z}(M_{position})$

This ensures that predicted position offsets are 0 at the beginning of
training, causing Gaussians to be predicted on the template mesh. This is part
of a set of measures that improve training stability as further elaborated in
sections 3.4 and 3.3.

To obtain the actual 3D Gaussian primitives, we uniformly sample the predicted
2D maps $M$. Each Gaussian $\mathcal{G}^{i}$ has a fixed coordinate
$x_{uv}^{i}\in[0,1]^{2}$ in the template’s UV space and will query its
attributes from $M$:

$\displaystyle\mathcal{G}=\textsc{GridSample}(M,x_{uv})$

where the $\textsc{GridSample}(\cdot)$ operater performs lookup in the
discrete maps $M$ via bilinear interpolation. The sampling scheme governs how
many Gaussians will be created and is shared across all generated persons.
Note that the generator is mostly agnostic to how many Gaussians are sampled
meaning that the sampling scheme may be changed during training, e.g., to
increase the number of Gaussians when rendering at higher resolutions.

###  3.3. Gaussian Attribute Modeling

To ensure that the predicted Gaussian attributes always stay within a
reasonable domain, we employ activation functions. Most notably, for the
predicted positions, we add the Gaussian’s corresponding 3D position $v^{i}$
on the template mesh and limit the predicted offset with a tanh activation
function:

(5) $\displaystyle\mathcal{G}_{position}^{i}$ $\displaystyle\leftarrow
v^{i}+\gamma_{pos}\textsc{tanh}\left(\mathcal{G}_{position}^{i}\right)$

where $\gamma_{pos}$ is the maximum allowed predicted offset from the template
mesh. For our experiments on FFHQ, we use $\gamma_{pos}=0.25$ meaning that
Gaussians can move at most 25cm away from the template. This step is crucial,
as it prevents the Gaussians from leaving the training view frustums during
early stages of adversarial training.

###  3.4. Gaussian Geometry Regularization

Simply predicting 3D Gaussians in an adversarial setting leads to poor
results. The reasons for this are two-fold: (i) The 3D Gaussians react very
sensitively to gradient updates during early training stages where the
discriminator is not calibrated yet. (ii) The 3D Gaussian representation is
quite powerful. Therefore, the generator is capable of predicting scenes that
look plausible when rendered, but the underlying geometry lacks realism.
We find that a set of regularizations that implement simple intuitions
effectively stabilize training as well as improve the underlying geometry. To
tackle training stability, we use a simple L2 regularization on the predicted
position offsets and Gaussian scales. To improve the geometric fidelity of
generated 3D representations, we propose to use a beta regularization on the
predicted offsets as well as a novel total variation loss on UV renderings.

####  3.4.1. Gaussian Position Regularization

A common cause of instabilities during training are Gaussians moving around
too much. We prevent this behavior with a position regularization on predicted
offsets:

(6) $\displaystyle\mathcal{L}_{reg}^{pos}=\|M_{position}\|_{2}$

This term encourages predicted Gaussians to stay close to the template mesh.

####  3.4.2. Gaussian Scale Regularization

As opposed to implicit scene representations, such as NeRFs, 3D Gaussians are
very sensitive to just a few parameters. For example, it is possible for a
single Gaussian to grow tremendously in scale, covering the whole view frustum
with just one primitive. Such cases are detrimental for training stability. We
employ a simple regularization on the predicted Gaussian scales to avoid such
degenerate 3D representations:

(7) $\displaystyle\mathcal{L}_{reg}^{scale}=\|M_{scale}-\gamma_{scale}\|_{2}$

where $\gamma_{scale}=e^{-5}$ is a target scale specified to be just large
enough that Gaussians would cover the template mesh.

####  3.4.3. Gaussian Opacity Regularization

To improve the geometric fidelity of the generated Gaussian representation, we
employ a beta regularization on the opacities, encouraging that Gaussians
should be either fully transparent or fully opaque:

(8) $\displaystyle\mathcal{L}_{reg}^{opac}=\textsc{Beta}(M_{opacity})$

where Beta is the negative log-likelihood of a Beta(0.5, 0.5) distribution
(Lombardi et al., 2019).

####  3.4.4. UV Total Variation Loss

A naive implementation of our method generates Gaussian representations with
many degrees of freedom. As a result, it can happen during training that the
generator “fakes” high-frequency detail by letting colors from further back
shine through. This still generates high-quality images but causes noticeable
artifacts when rendering videos. As the discriminator is not aware of the time
dimension, this behavior is never punished during training. We therefore
design a novel regularization scheme to improve the geometric fidelity of
generated representations. The intuition behind our regularization is that
neighboring pixels in a rendered image should be modelled by Gaussians that
are also close in UV space. This naturally leads to smooth surfaces and
penalizes holes. Formally, we apply a total variation (TV) regularization on
rendered UV images as follows:

(9) $\displaystyle\mathcal{\hat{G}}$ $\displaystyle\leftarrow\mathcal{G}$ (10)
$\displaystyle\mathcal{\hat{G}}_{color}^{i}$
$\displaystyle\leftarrow(u^{i},v^{i},0)$ (11) $\displaystyle R_{uv}$
$\displaystyle=\mathcal{R}(\mathcal{\hat{G}},\pi)$ (12) $\displaystyle
R_{uv}^{\prime}$ $\displaystyle=\frac{R_{uv}-(1-R_{\alpha})}{R_{\alpha}}$ (13)
$\displaystyle\mathcal{L}_{uv}$ $\displaystyle=TV(R_{uv}^{\prime})$

where $u^{i},v^{i}=x^{i}_{uv}$ are a Gaussian’s UV-coordinates, $\pi$ is the
camera, and $R_{\alpha}\in\mathbb{R}^{H\times W}$ is the rendered alpha map
containing the per-pixel accumulations from rasterization. We refer to
$R_{uv}$ as a UV rendering. Equation 12 reverses the alpha compositing that is
performed during rasterization and is necessary in order to disregard
transparency for the regularization. Otherwise, the border region between
foreground and background, which is quite likely to include semi-transparent
pixels due to the fuzzy nature of Gaussians, would contribute consistently to
the regularization term, incentivizing foreground Gaussians to cover the
background. With our novel UV total variation loss in place, holes in the 3D
representation are stitched during training, noticeably improving the
generated geometry. As a positive side effect, the regularization also removes
other types of artifacts such as floating Gaussians that can appear in front
of the face. Figure 8 shows how the generated 3D representation is improved by
our novel regularization term.

###  3.5. Training Pipeline

We generate UV maps $M$ at $256^{2}$ resolution and initially sample one
Gaussian per UV map texel resulting in $65k$ Gaussians. We begin training at
$256\times 256$ rendering resolution and apply progressive growing after
$7000k$ training images. Since the generator is resolution agnostic, we only
need additional layers for the discriminator when increasing the rendering
resolution. For progressive growing, we add a CNN block in front of the
discriminator with a skip connection to ensure a smooth transition and blend
in the contributions from the untrained layers over the course of $1000k$
images. We also increase the number of Gaussians when increasing the rendering
resolution to facilitate modelling more high-frequency details. To achieve
this, we simply increase the sampling density of Gaussians in UV space from
$256^{2}$ to $512^{2}$ for a total of $262k$ Gaussians. This sudden increase
in Gaussians causes the predictions to become temporarily more blurry.
Nevertheless, training stability is not adversely affected by that and the
generator quickly adapts to the presence of more Gaussians. Also note that the
Gaussians are never all active at the same time for any given input latent
$z$. The generator uses the opacity attribute to switch Gaussians on and off
as needed, avoiding topological issues that would be implied by the template
mesh. In practice, around half of the Gaussians are active on average.
The final optimization term for the generator is as follows:

(14) $\displaystyle\mathcal{L}_{G}$
$\displaystyle=\mathcal{L}_{adv}+\lambda_{p}\mathcal{L}_{reg}^{pos}+\lambda_{s}\mathcal{L}_{reg}^{scale}+\lambda_{o}\mathcal{L}_{reg}^{opac}+\lambda_{uv}\mathcal{L}_{uv}$
(15) $\displaystyle\mathcal{L}_{adv}$
$\displaystyle=\textsc{softplus}(-D(\mathcal{R}(\mathcal{G}),\pi))$

where $\mathcal{L}_{adv}$ is the standard non-saturating GAN loss (Goodfellow
et al., 2014) applied on the rendered 3D Gaussian representation. The
discriminator $D$ is trained with R1 gradient regularization (Mescheder et
al., 2018) with weight $1$. For the regularization weights, we choose
$\lambda_{p}=0.1$ and $\lambda_{s}=0.05$. The opacity and UV total variation
loss are initially turned off for training at $256^{2}$ resolution and are set
to $\lambda_{o}=1$ and $\lambda_{uv}=100$ as soon as progressive growing is
applied. We use a batch size of $32$ and learning rates of $0.02$ for both
generator and discriminator using the Adam optimizer (Kingma and Ba, 2014). In
total, we train our model for 25M images following Chan et al. (2022) which
takes roughly 12 days on 2 RTX A6000 GPUs with 48G of VRAM each.

depth |  normals — rendering |  rendering |  |  depth |  normals — rendering |  rendering
---|---|---|---|---|---|---

Figure 3. Generated Geometry: We show samples from our model with
corresponding depth and normal maps. GGHead can generate fine geometric
details such as eyeglasses and hair strands.

188 fps |  42 fps |  25 fps |  3.5 fps |  228 fps
---|---|---|---|---
Gaussian Shell Maps |  EG3D |  GRAM-HD |  Mimic3D |  Ours

Figure 4. Qualitative Comparison: We show uncurated samples (seeds 0-4) of
different 3D GANs trained on FFHQ at $512^{2}$ resolution. Our method matches
the quality of existing approaches while being strictly 3D-consistent and one
order of magnitude faster to render. Times measured on an RTX A6000 GPU with a
batch size of 1.

EG3D (Chan et al., 2022) |  Ours
---|---

Figure 5. Analysis of 3D Consistency: We show spatio-temporal line textures
akin to the Epipolar Line Images (EPI) (Bolles et al., 1987) obtained by
rotating the camera horizontally and stacking the pixels of a fixed horizontal
line segment. EG3D suffers from texture-sticking artifacts due to its use of a
2D super-resolution network, leading to staircases in the epiploar line
images. In contrast, our method provides smooth renderings without any
flickering.

GRAM-HD (Xiang et al., 2023) |  |  Mimic3D (Chen et al., 2023) |  |  Ours
---|---|---|---|---

Figure 6. Analysis of Side-views: We show renderings of GRAM-HD, Mimic3D and
our model from extreme side-views. GRAM-HD exposes the peculiarities of its
underlying 3D representation in this scenario since the predicted 2D manifolds
often run in parallel, causing striped artifacts. Mimic3D suffers from blurred
colors when viewed from the side. In contrast, the 3D Gaussian representation
generated by our method still looks appealing.

##  4\. Experimental Results

###  4.1. Dataset

We train our models on data from the FFHQ dataset (Karras et al., 2019b),
which consists of 70k mostly frontal images of human faces. We use the
preprocessed version of Chan et al. (2022) who aligned the images via facial
landmarks, cropped them into $512^{2}$ resolution, and computed camera poses
by resorting to a 3DMM fitting (Deng et al., 2019). We further follow EG3D’s
processing pipeline to obtain a $1024^{2}$ resolution dataset of cropped and
aligned face images with corresponding camera poses from the original FFHQ in-
the-wild dataset. Finally, we employ the off-the-shelf background matting
network MODNet (Ke et al., 2022) to replace the background with white pixels.
That way, our models can focus their generation power solely on the actual 3D
head. We conduct all our comparisons on this masked 2D face image dataset
derived from FFHQ.

###  4.2. Baselines

We compare our method to several competitive baseline methods from the 3D GAN
literature. The baselines are selected based on three criteria: quality,
speed/resolution, and architectural similarity.

EG3D (Chan et al., 2022) is a state-of-the-art 3D GAN that synthesizes high-
quality face images by relying on a super-resolution network to upsample low-
resolution renderings from generated TriPlanes.

Mimic3D (Chen et al., 2023) is SoTA among fully 3D-consistent GANs on FFHQ. It
distills a TriPlane generator backbone by directly supervising rendered images
with pseudo-GT multi-view images obtained from EG3D.

GRAM-HD (Xiang et al., 2023) proposed an efficient radiance field
representation based on 2D manifolds, enabling faster rendering and training
at $1024^{2}$ resolution.

Gaussian Shell Maps (GSM) (Abdal et al., 2023) shares important architectural
design decisions with our method, most notably that it also predicts 3D
Gaussians in the UV space of a template.

For EG3D and Mimic3D, we use the official implementation and retrain their
methods on our masked FFHQ dataset. Since Gaussian Shell Maps was originally
proposed for animatable full body synthesis, we adapt their implementation and
retrain it on our masked FFHQ dataset. For GRAM-HD, we use official model
checkpoints to conduct qualitative and timing comparisons.

###  4.3. Qualitative Comparison

Figure 3 shows samples from our model alongside extracted depth and normal
maps. The generated 3D representations exhibit a great amount of geometric
detail and diversity. We further compare samples from GGHead to those of EG3D,
Mimic3D, GRAM-HD and GSM in fig. 4. Our method produces renderings of great
quality matching the current state-of-the-art at much higher speed.
Interestingly, GSM synthesizes images with noticeable artifacts, indicating
that generating 3D Gaussians alone is not sufficient for high-quality
renderings. We attribute this to the design of GSM which disallows Gaussians
from moving, severely restricting the expressiveness of the generated 3D
representation.
Next, we analyse the view-consistency of generated images from EG3D and our
method. EG3D generally can synthesize high quality static images, but this
comes at the cost of view-inconsistencies due to the use of a 2D super-
resolution network. In fig. 5, we present the spatio-temporal texture images
similar to Epipolar Line Images (EPI) (Bolles et al., 1987) for EG3D and our
method. Specifically, we smoothly rotate the camera around the generated head
while stacking the pixels from a fixed horizontal line segment for each
timestep. For multi-view consistent images, the resulting spatio-temporal
texture should appear smooth. Conversely, edged artifacts and noise in the
spatio-temporal texture indicate 3D inconsistencies. For EG3D, such artifacts
are prevalent and manifest as flickering in video rendering. In contrast, our
method generates strictly 3D-consistent images.
Furthermore, we investigate the generations of Mimic3D and GRAM-HD in more
detail in fig. 6. We find that both methods generally produce high-quality
renderings but expose idiosyncracies of the underlying 3D representation when
viewed from extreme angles. While GRAM-HD’s use of 2D manifolds becomes
visible in this scenario, Mimic3D produces slightly fuzzy surfaces, which is
typical for NeRFs, causing blurry renderings for side views. Instead, our
method produces more well-behaved renderings for such extreme camera angles.
We attribute this to the fact that generating 3D Gaussian primitives in UV
space forces them to lie on a single 2D manifold. This acts as an inductive
bias towards more surface-like 3D representations which helps GGHead to
generalize to such challenging scenarios.

Table 1. Quantitative comparison: We compare our method to other 3D GANs on FFHQ at $512^{2}$ resolution. Our method fares favorably against other fully 3D-consistent methods and comes close to the performance of EG3D despite not using any 2D super-resolution. GGHead also produces more consistent 3D surfaces while simultaneously being considerably faster. |  | Method | FID $\downarrow$ | $\mathrm{PSNR}_{mv}\uparrow$ | $\mathrm{SSIM}_{mv}\uparrow$
---|---|---|---|---|---
FFHQ@512 | w/ SR | EG3D | 3.28 | 28.79 | 0.886
w/o SR | GSM | 28.19 | 25.53 | 0.895
Mimic3D | 4.27 | 28.07 | 0.921
Ours | 4.06 | 33.55 | 0.959
Table 2. 3D GAN performance benchmark: We measure the time it takes to generate 3D representations (Gen.), render them (Rend.), and perform a full forward-backward pass during training (Forw. + Backw.). ”Res.” denotes the native rendering resolution excluding super-resolution. Note that EG3D and Mimic3D have to lower the native rendering resolution substantially during training. In contrast, our method runs much faster at full resolution without sacrificing rendering quality. Times averaged over 1000 images. | Train | Inference
---|---|---
| Res. | Forw. + Backw.$\downarrow$ | Res. | Gen.$\downarrow$ | Rend.$\downarrow$
EG3D | 128 | 209.61 ms | 128 | 6.1 ms | 23.4 ms
Mimic3D | 64 | 298.73 ms | 512 | 20.3 ms | 286.7 ms
GRAM-HD | - | - | 512 | 157.6 ms | 38.2 ms
GSM | 512 | 120.16 ms | 512 | 8.8 ms | 5.2 ms
Ours | 512 | 127.78 ms | 512 | 10.4 ms | 3.8 ms
Figure 7. Generation and Rendering Resolution Benchmark: We measure the
rendering times (including generation of the 3D representation) on a single
RTX A6000 GPU with a batch size of 1. Our method can effortlessly generate
high-quality 3D representations and render them at 1k native resolution in
real-time while still matching other competitive 3D GANs in terms of FID.

###  4.4. Quantitative Comparison

We conduct quantitative comparisons with EG3D, Mimic3D, and Gaussian Shell
Maps on the masked FFHQ dataset at $512^{2}$ resolution. Since there is no
direct ground truth to evaluate the quality of generated 3D representations,
we resort to measuring the Fréchet Inception Distance (FID) (Heusel et al.,
2017). To do that, we generate 50k samples for each trained model using the
training camera distribution and compare its statistics to the full dataset.
The results are shown in table 1. Our method matches the performance of
current state-of-the-art 3D GANs, obtaining slightly better FID scores than
Mimic3D. It is only surpassed by EG3D, which employs a super-resolution module
to achieve its rendering quality. To quantify the adverse effects of screen-
space super-resolution on 3D consistency, we follow the procedure laid out by
GRAM-HD (Xiang et al., 2023). For each method, we generate 30 images from
different viewpoints of the same person, and then train the surface
reconstruction method NeuS2 (Wang et al., 2023) on the generated images. We
then measure the reconstruction error with PSNR and SSIM scores. The reasoning
behind this is that view-inconsistencies in the images given to NeuS2 will
lead to a higher reconstruction error. In practice, we observe that this
metric also measures “surfaceness” of a 3D representation, since NeuS2 favors
smooth and well-defined surfaces by design. We report the multi-view surface
reconstruction results averaged over 30 randomly generated heads in table 1.
GGHead leads to significantly better surface reconstructions than EG3D,
underscoring its great 3D consistency.

###  4.5. Timing Comparison

A core motivation for our method is to lower the time needed for generating
and rendering 3D heads, ultimately leading to faster training of 3D GANs at
higher resolutions. We benchmark training and inference times of existing 3D
GANs, and compare them to our method in table 2. As expected, we observe a
huge performance gap between the NeRF-based approaches (EG3D, Mimic3D, GRAM-
HD) and methods that employ 3D Gaussian Splatting (GSM, Ours). This is true
not only for the actual rendering of the 3D representation itself, but also
for the required time for a full training iteration where rendering also has
to be performed. Note that due to expensive ray marching, EG3D and Mimic3D
have considerably slower training iterations despite rendering only
$\frac{1}{16}$th of the pixels.
We further study the scalability of existing 3D GAN methods by measuring the
required time for a full forward pass under different rendering resolutions.
The scaling behavior is plotted in fig. 7. Again, our method demonstrates a
considerably better scaling behavior on larger rendering resolutions thanks to
the efficient 3D Gaussian Splatting rasterizer.

###  4.6. Ablations

$R_{uv}$ |  Rendered Image |  |  |  $R_{uv}$ |  Rendered Image |
---|---|---|---|---|---|---
without $\mathcal{L}_{uv}$ |  | Ours

Figure 8. Effect of TV UV regularization: Without $\mathcal{L}_{uv}$, the
predicted Gaussian geometry is flawed. Two common failure cases are the
emergence of floating lines of Gaussians that are especially visible in video
renderings, and improper surfaces where skin texture is created by letting
Gaussians from the back of the head shine through. Both cases are easily
detectable in the UV renderings $R_{uv}$. Our novel UV total variation loss
exploits this fact and effectively addresses both issues by enforcing $R_{uv}$
to be smooth.

In the following, we analyze important design decisions of our method,
including the choice of template mesh, the number of generated Gaussians, and
the effect of our novel UV total variation loss.

####  4.6.1. Effect of Template Mesh

In table 3, we compare multiple versions of GGHead trained with different
template meshes on the masked FFHQ dataset at $256^{2}$ resolution. We
consider different types of template meshes, ranging from simple domain-
agnostic shapes, like planes and spheres, to meshes that are designed for the
face domain. Here, we compare two different head meshes: FLAME (Li et al.,
2017) and an adapted version of the FLAME mesh taken from Abdal et al. (2023)
which has a more efficient UV layout and the back of the head removed. It
turns out that GGHead is fairly robust to the exact choice of template mesh
since all templates produce reasonable results in terms of FID, with the
adjusted FLAME mesh scoring best. We hypothesize that the efficiency of the UV
layout and the density of allocated points in the frontal face region are more
important than the actual shape of the template mesh.

####  4.6.2. Effect of UV Total Variation Loss

Being fundamentally a 3D point representation, the 3D Gaussians lack the
notion of a proper surface. Consequently, it can happen during training that
Gaussians in the face area become too small, allowing other Gaussians from
further back to shine through. This can create high frequency detail that
satisfies the discriminator loss and looks plausible in static renderings.
However, such a 3D representation remains flawed and is easily exposed in
video renderings. Since we cannot supervise the time dimension during
training, we instead employ our novel UV total variation loss. The effect can
be seen in fig. 8. The regularization term effectively reduces the amount of
Gaussians that shine through, considerably improving the surface quality of
the generated 3D representations. As a side effect, the loss also resolves
another failure case where lines of Gaussians would appear that mimic the
appearance of hair strands but occur out of place. In table 3, we also study
the effect of the UV total variation loss quantitatively and find that it
slightly improves FID.

Table 3. Ablation study: We analyze the effect of different template meshes, different numbers of Gaussians, and our novel UV total variation loss. | Method | FID$\downarrow$
---|---|---
256 | Sphere Template | 4.59
Plane Template | 4.27
FLAME Template | 4.74
Adjusted FLAME Template | 4.03
512 | 65k Gaussians | 6.39
262k Gaussians | 4.28
\+ $\mathcal{L}_{uv}$ (Ours) | 4.06

####  4.6.3. Effect of Number of Gaussians

Due to our generator design, we can flexibly sample different numbers of
Gaussians from the same model. We can use this to our advantage when
progressively growing to larger resolutions where more high-frequency detail
needs to be synthesized. As opposed to directly training with all Gaussians
from the beginning, this has two advantages: Training is both faster and more
stable due to fewer degrees of freedom. We analyze the effect of more
Gaussians on the final results in table 3 where increasing the number of
Gaussians (262k) gives a noticeable improvement in quality over keeping the
same number of Gaussians (65k) as in lower resolution training.

##  5\. Limitations and Future Work

We have demonstrated that GGHead can generate and render 3D heads at high
quality. However, the generated 3D representations only provide a user with
viewpoint control. Having additional control over the facial expressions would
enable further use-cases such as 3D-consistent expression editing of single
images. Since our method already employs a coarse template mesh, it stands to
reason to extend it in similar fashion as Qian et al. (2023a); Sun et al.
(2023); Abdal et al. (2023), e.g., by employing FLAME and training on FFHQ
with corresponding expression codes. One step further, the improved
scalability of our method may enable training on large facial video datasets
(Zhu et al., 2022), ultimately paving the way for building a photorealistic
3DMM from 2D data alone.
Another axis of improvement is generality. We have shown how a 3D generative
model can be obtained with our pipeline for the narrow domain of human heads.
However, extending our method to other domains such as ImageNet categories
(Deng et al., 2009) could enable learning more generic high-quality 3D priors.
This could be done by having one template mesh per category and making the
mesh itself learnable. As a result, the model could uncover suitable templates
on its own. Additionally, an approach like that of Skorokhodov et al. (2023)
could be employed to drop the requirement for domain-specific keypoint
detectors and alignment procedure.

##  6\. Conclusion

In this work, we have proposed Generative Gaussian Heads (GGHead), a fast
method for training 3D generative models for human heads on large 2D data
collections. At the core of our approach lies the combination of a 3D GAN
pipeline with efficient rendering from 3D Gaussian Splatting. We achieve this
by parameterizing 3D Gaussian Heads as a set of 2D maps that live in the UV
space of a template mesh and that can be efficiently predicted with powerful
2D CNN architectures. We further propose a novel UV total variation loss that
exploits the design of our generator to improve geometric fidelity of
generated 3D heads. In experiments on FFHQ, we demonstrate that our method
matches the image synthesis quality of the current state of the art while
being both fully 3D-consistent and considerably faster. In timing benchmarks,
we confirm the great scalability of our method giving rise to real-time
generation and rendering of photo-realistic 3D heads at $1024^{2}$ resolution.

### Acknowledgements

This work was supported by the ERC Starting Grant Scan2CAD (804724) and the
German Research Foundation (DFG) Research Unit “Learning and Simulation in
Visual Computing”. We would also like to thank Karla Weighart for proofreading
and Angela Dai for the video voice-over.

## References

  * (1)
  * Abdal et al. (2023) Rameen Abdal, Wang Yifan, Zifan Shi, Yinghao Xu, Ryan Po, Zhengfei Kuang, Qifeng Chen, Dit-Yan Yeung, and Gordon Wetzstein. 2023.  Gaussian shell maps for efficient 3d human generation.  _arXiv preprint arXiv:2311.17857_ (2023).
  * Bolles et al. (1987) Robert C Bolles, H Harlyn Baker, and David H Marimont. 1987.  Epipolar-plane image analysis: An approach to determining structure from motion.  _International journal of computer vision_ 1, 1 (1987), 7–55.
  * Chan et al. (2022) Eric R Chan, Connor Z Lin, Matthew A Chan, Koki Nagano, Boxiao Pan, Shalini De Mello, Orazio Gallo, Leonidas J Guibas, Jonathan Tremblay, Sameh Khamis, et al. 2022\.  Efficient geometry-aware 3d generative adversarial networks. In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_. 16123–16133.
  * Chan et al. (2021) Eric R Chan, Marco Monteiro, Petr Kellnhofer, Jiajun Wu, and Gordon Wetzstein. 2021.  pi-gan: Periodic implicit generative adversarial networks for 3d-aware image synthesis. In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_. 5799–5809.
  * Charatan et al. (2023) David Charatan, Sizhe Li, Andrea Tagliasacchi, and Vincent Sitzmann. 2023.  pixelsplat: 3d gaussian splats from image pairs for scalable generalizable 3d reconstruction.  _arXiv preprint arXiv:2312.12337_ (2023).
  * Chen et al. (2023) Xingyu Chen, Yu Deng, and Baoyuan Wang. 2023.  Mimic3d: Thriving 3d-aware gans via 3d-to-2d imitation. In _2023 IEEE/CVF International Conference on Computer Vision (ICCV)_. IEEE Computer Society, 2338–2348.
  * Deng et al. (2009) Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei. 2009.  Imagenet: A large-scale hierarchical image database. In _2009 IEEE conference on computer vision and pattern recognition_. Ieee, 248–255.
  * Deng et al. (2022) Yu Deng, Jiaolong Yang, Jianfeng Xiang, and Xin Tong. 2022.  Gram: Generative radiance manifolds for 3d-aware image generation. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_. 10673–10683.
  * Deng et al. (2019) Yu Deng, Jiaolong Yang, Sicheng Xu, Dong Chen, Yunde Jia, and Xin Tong. 2019.  Accurate 3d face reconstruction with weakly-supervised learning: From single image to image set. In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition workshops_. 0–0.
  * Fu et al. (2023) Yang Fu, Sifei Liu, Amey Kulkarni, Jan Kautz, Alexei A Efros, and Xiaolong Wang. 2023.  Colmap-free 3d gaussian splatting.  _arXiv preprint arXiv:2312.07504_ (2023).
  * Goodfellow et al. (2014) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. 2014.  Generative adversarial nets.  _Advances in neural information processing systems_ 27 (2014).
  * Gu et al. (2021) Jiatao Gu, Lingjie Liu, Peng Wang, and Christian Theobalt. 2021.  Stylenerf: A style-based 3d-aware generator for high-resolution image synthesis.  _arXiv preprint arXiv:2110.08985_ (2021).
  * Heusel et al. (2017) Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter. 2017.  Gans trained by a two time-scale update rule converge to a local nash equilibrium.  _Advances in neural information processing systems_ 30 (2017).
  * Hu and Liu (2023) Shoukang Hu and Ziwei Liu. 2023.  Gauhuman: Articulated gaussian splatting from monocular human videos.  _arXiv preprint arXiv:2312.02973_ (2023).
  * Jiang et al. (2024) Yujiao Jiang, Qingmin Liao, Xiaoyu Li, Li Ma, Qi Zhang, Chaopeng Zhang, Zongqing Lu, and Ying Shan. 2024.  UV Gaussians: Joint Learning of Mesh Deformation and Gaussian Textures for Human Avatar Modeling.  _arXiv preprint arXiv:2403.11589_ (2024).
  * Karras et al. (2020a) Tero Karras, Miika Aittala, Janne Hellsten, Samuli Laine, Jaakko Lehtinen, and Timo Aila. 2020a.  Training Generative Adversarial Networks with Limited Data. In _Proc. NeurIPS_.
  * Karras et al. (2019a) Tero Karras, Samuli Laine, and Timo Aila. 2019a.  A style-based generator architecture for generative adversarial networks. In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_. 4401–4410.
  * Karras et al. (2019b) Tero Karras, Samuli Laine, and Timo Aila. 2019b.  A style-based generator architecture for generative adversarial networks. In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_. 4401–4410.
  * Karras et al. (2020b) Tero Karras, Samuli Laine, Miika Aittala, Janne Hellsten, Jaakko Lehtinen, and Timo Aila. 2020b.  Analyzing and improving the image quality of stylegan. In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_. 8110–8119.
  * Ke et al. (2022) Zhanghan Ke, Jiayu Sun, Kaican Li, Qiong Yan, and Rynson WH Lau. 2022.  Modnet: Real-time trimap-free portrait matting via objective decomposition. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 36. 1140–1147.
  * Kerbl et al. (2023) Bernhard Kerbl, Georgios Kopanas, Thomas Leimkühler, and George Drettakis. 2023.  3d gaussian splatting for real-time radiance field rendering.  _ACM Transactions on Graphics_ 42, 4 (2023), 1–14.
  * Kingma and Ba (2014) Diederik P Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * Li et al. (2017) Tianye Li, Timo Bolkart, Michael. J. Black, Hao Li, and Javier Romero. 2017.  Learning a model of facial shape and expression from 4D scans.  _ACM Transactions on Graphics, (Proc. SIGGRAPH Asia)_ 36, 6 (2017), 194:1–194:17.  https://doi.org/10.1145/3130800.3130813
  * Li et al. (2023) Xinhai Li, Huaibin Wang, and Kuo-Kun Tseng. 2023.  GaussianDiffusion: 3D Gaussian Splatting for Denoising Diffusion Probabilistic Models with Structured Noise.  _arXiv preprint arXiv:2311.11221_ (2023).
  * Ling et al. (2023) Huan Ling, Seung Wook Kim, Antonio Torralba, Sanja Fidler, and Karsten Kreis. 2023.  Align your gaussians: Text-to-4d with dynamic 3d gaussians and composed diffusion models.  _arXiv preprint arXiv:2312.13763_ (2023).
  * Lombardi et al. (2019) Stephen Lombardi, Tomas Simon, Jason Saragih, Gabriel Schwartz, Andreas Lehrmann, and Yaser Sheikh. 2019.  Neural volumes: Learning dynamic renderable volumes from images.  _arXiv preprint arXiv:1906.07751_ (2019).
  * Mescheder et al. (2018) Lars Mescheder, Andreas Geiger, and Sebastian Nowozin. 2018.  Which training methods for gans do actually converge?. In _International conference on machine learning_. PMLR, 3481–3490.
  * Mildenhall et al. (2021) Ben Mildenhall, Pratul P Srinivasan, Matthew Tancik, Jonathan T Barron, Ravi Ramamoorthi, and Ren Ng. 2021.  Nerf: Representing scenes as neural radiance fields for view synthesis.  _Commun. ACM_ 65, 1 (2021), 99–106.
  * Moreau et al. (2023) Arthur Moreau, Jifei Song, Helisa Dhamo, Richard Shaw, Yiren Zhou, and Eduardo Pérez-Pellitero. 2023.  Human gaussian splatting: Real-time rendering of animatable avatars.  _arXiv preprint arXiv:2311.17113_ (2023).
  * Niemeyer and Geiger (2021) Michael Niemeyer and Andreas Geiger. 2021.  Giraffe: Representing scenes as compositional generative neural feature fields. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_. 11453–11464.
  * Or-El et al. (2022) Roy Or-El, Xuan Luo, Mengyi Shan, Eli Shechtman, Jeong Joon Park, and Ira Kemelmacher-Shlizerman. 2022.  Stylesdf: High-resolution 3d-consistent image and geometry generation. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_. 13503–13513.
  * Qian et al. (2023a) Shenhan Qian, Tobias Kirschstein, Liam Schoneveld, Davide Davoli, Simon Giebenhain, and Matthias Nießner. 2023a.  Gaussianavatars: Photorealistic head avatars with rigged 3d gaussians.  _arXiv preprint arXiv:2312.02069_ (2023).
  * Qian et al. (2023b) Zhiyin Qian, Shaofei Wang, Marko Mihajlovic, Andreas Geiger, and Siyu Tang. 2023b.  3dgs-avatar: Animatable avatars via deformable 3d gaussian splatting.  _arXiv preprint arXiv:2312.09228_ (2023).
  * Ren et al. (2023) Jiawei Ren, Liang Pan, Jiaxiang Tang, Chi Zhang, Ang Cao, Gang Zeng, and Ziwei Liu. 2023.  Dreamgaussian4d: Generative 4d gaussian splatting.  _arXiv preprint arXiv:2312.17142_ (2023).
  * Saito et al. (2023) Shunsuke Saito, Gabriel Schwartz, Tomas Simon, Junxuan Li, and Giljoo Nam. 2023.  Relightable gaussian codec avatars.  _arXiv preprint arXiv:2312.03704_ (2023).
  * Schonberger and Frahm (2016) Johannes L Schonberger and Jan-Michael Frahm. 2016.  Structure-from-motion revisited. In _Proceedings of the IEEE conference on computer vision and pattern recognition_. 4104–4113.
  * Schwarz et al. (2020) Katja Schwarz, Yiyi Liao, Michael Niemeyer, and Andreas Geiger. 2020.  Graf: Generative radiance fields for 3d-aware image synthesis.  _Advances in Neural Information Processing Systems_ 33 (2020), 20154–20166.
  * Skorokhodov et al. (2023) Ivan Skorokhodov, Aliaksandr Siarohin, Yinghao Xu, Jian Ren, Hsin-Ying Lee, Peter Wonka, and Sergey Tulyakov. 2023.  3D generation on ImageNet. In _International Conference on Learning Representations_.  https://openreview.net/forum?id=U2WjB9xxZ9q
  * Skorokhodov et al. (2022) Ivan Skorokhodov, Sergey Tulyakov, Yiqun Wang, and Peter Wonka. 2022.  Epigraf: Rethinking training of 3d gans.  _Advances in Neural Information Processing Systems_ 35 (2022), 24487–24501.
  * Sun et al. (2022) Jingxiang Sun, Xuan Wang, Yichun Shi, Lizhen Wang, Jue Wang, and Yebin Liu. 2022.  Ide-3d: Interactive disentangled editing for high-resolution 3d-aware portrait synthesis.  _ACM Transactions on Graphics (ToG)_ 41, 6 (2022), 1–10.
  * Sun et al. (2023) Jingxiang Sun, Xuan Wang, Lizhen Wang, Xiaoyu Li, Yong Zhang, Hongwen Zhang, and Yebin Liu. 2023.  Next3d: Generative neural texture rasterization for 3d-aware head avatars. In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_. 20991–21002.
  * Svitov et al. (2024) David Svitov, Pietro Morerio, Lourdes Agapito, and Alessio Del Bue. 2024.  HAHA: Highly Articulated Gaussian Human Avatars with Textured Mesh Prior.  _arXiv preprint arXiv:2404.01053_ (2024).
  * Tang et al. (2023) Jiaxiang Tang, Jiawei Ren, Hang Zhou, Ziwei Liu, and Gang Zeng. 2023.  Dreamgaussian: Generative gaussian splatting for efficient 3d content creation.  _arXiv preprint arXiv:2309.16653_ (2023).
  * Trevithick et al. (2024) Alex Trevithick, Matthew Chan, Towaki Takikawa, Umar Iqbal, Shalini De Mello, Manmohan Chandraker, Ravi Ramamoorthi, and Koki Nagano. 2024.  What You See is What You GAN: Rendering Every Pixel for High-Fidelity Geometry in 3D GANs.  _arXiv preprint arXiv:2401.02411_ (2024).
  * Wang et al. (2023) Yiming Wang, Qin Han, Marc Habermann, Kostas Daniilidis, Christian Theobalt, and Lingjie Liu. 2023.  Neus2: Fast learning of neural implicit surfaces for multi-view reconstruction. In _Proceedings of the IEEE/CVF International Conference on Computer Vision_. 3295–3306.
  * Wewer et al. (2024) Christopher Wewer, Kevin Raj, Eddy Ilg, Bernt Schiele, and Jan Eric Lenssen. 2024.  latentSplat: Autoencoding Variational Gaussians for Fast Generalizable 3D Reconstruction.  _arXiv preprint arXiv:2403.16292_ (2024).
  * Xiang et al. (2023) Jianfeng Xiang, Jiaolong Yang, Yu Deng, and Xin Tong. 2023.  Gram-hd: 3d-consistent image generation at high resolution with generative radiance manifolds. In _Proceedings of the IEEE/CVF International Conference on Computer Vision_. 2195–2205.
  * Xu et al. (2023) Yuelang Xu, Benwang Chen, Zhe Li, Hongwen Zhang, Lizhen Wang, Zerong Zheng, and Yebin Liu. 2023.  Gaussian head avatar: Ultra high-fidelity head avatar via dynamic gaussians.  _arXiv preprint arXiv:2312.03029_ (2023).
  * Xu et al. (2024) Yinghao Xu, Zifan Shi, Wang Yifan, Hansheng Chen, Ceyuan Yang, Sida Peng, Yujun Shen, and Gordon Wetzstein. 2024.  Grm: Large gaussian reconstruction model for efficient 3d reconstruction and generation.  _arXiv preprint arXiv:2403.14621_ (2024).
  * Xue et al. (2022) Yang Xue, Yuheng Li, Krishna Kumar Singh, and Yong Jae Lee. 2022.  Giraffe hd: A high-resolution 3d-aware generative model. In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_. 18440–18449.
  * Yi et al. (2023) Taoran Yi, Jiemin Fang, Guanjun Wu, Lingxi Xie, Xiaopeng Zhang, Wenyu Liu, Qi Tian, and Xinggang Wang. 2023.  Gaussiandreamer: Fast generation from text to 3d gaussian splatting with point cloud priors.  _arXiv preprint arXiv:2310.08529_ (2023).
  * Zhang et al. (2023) Lvmin Zhang, Anyi Rao, and Maneesh Agrawala. 2023.  Adding conditional control to text-to-image diffusion models. In _Proceedings of the IEEE/CVF International Conference on Computer Vision_. 3836–3847.
  * Zhu et al. (2022) Hao Zhu, Wayne Wu, Wentao Zhu, Liming Jiang, Siwei Tang, Li Zhang, Ziwei Liu, and Chen Change Loy. 2022.  CelebV-HQ: A large-scale video facial attributes dataset. In _European conference on computer vision_. Springer, 650–667.
  * Zielonka et al. (2023) Wojciech Zielonka, Timur Bagautdinov, Shunsuke Saito, Michael Zollhöfer, Justus Thies, and Javier Romero. 2023.  Drivable 3d gaussian avatars.  _arXiv preprint arXiv:2311.08581_ (2023).
  * Zou et al. (2023) Zi-Xin Zou, Zhipeng Yu, Yuan-Chen Guo, Yangguang Li, Ding Liang, Yan-Pei Cao, and Song-Hai Zhang. 2023.  Triplane meets gaussian splatting: Fast and generalizable single-view 3d reconstruction with transformers.  _arXiv preprint arXiv:2312.09147_ (2023).
  * Zwicker et al. (2001) Matthias Zwicker, Hanspeter Pfister, Jeroen Van Baar, and Markus Gross. 2001.  EWA volume splatting. In _Proceedings Visualization, 2001. VIS’01._ IEEE, 29–538.


# Immunity to Increasing Condition Numbers of Linear Superiorization versus
Linear Programming

Jan Schröder1, Yair Censor2, Philipp Süss1 and Karl-Heinz Küfer1

(1Optimization Department, Fraunhofer - ITWM, Kaiserslautern 67663, Germany
2Department of Mathematics, University of Haifa, Mt. Carmel, Haifa 3498838,
Israel
July 26, 2024)

###### Abstract

Given a family of linear constraints and a linear objective function one can
consider whether to apply a Linear Programming (LP) algorithm or use a Linear
Superiorization (LinSup) algorithm on this data. In the LP methodology one
aims at finding a point that fulfills the constraints and has the minimal
value of the objective function over these constraints. The Linear
Superiorization approach considers the same data as linear programming
problems but instead of attempting to solve those with linear optimization
methods it employs perturbation resilient feasibility-seeking algorithms and
steers them toward feasible points with reduced (not necessarily minimal)
objective function values. Previous studies compared LP and LinSup in terms of
their respective outputs and the resources they use. In this paper we
investigate these two approaches in terms of their sensitivity to condition
numbers of the system of linear constraints. Condition numbers are a measure
for the impact of deviations in the input data on the output of a problem and,
in particular, they describe the factor of error propagation when given wrong
or erroneous data. Therefore, the ability of LP and LinSup to cope with
increased condition numbers, thus with ill-posed problems, is an important
matter to consider which was not studied until now. We investigate
experimentally the advantages and disadvantages of both LP and LinSup on
examplary problems of linear programming with multiple condition numbers and
different problem dimensions.

Keywords: Condition number, ill-posed problem, Linear Programming, Linear
Superiorization, bounded perturbations, feasibility seeking, Scipy, Gurobi.

##  1 Introduction

In this paper we investigate and discuss the behavior of the linear
superiorization (LinSup) methodology for various different condition numbers
and compare it to the behavior of multiple linear optimization solvers. We set
up differently conditioned linear optimization problems in a manner that keeps
them structurally similar in order to keep them comparable to each other. Then
we solve these problems with the introduced superiorization methodology and
with established implementations of the simplex algorithm and an interior
point method. We compare the necessary runtimes in problems with multiple
dimensions and discuss the observed effects. We show experimentally that
superiorization is rather robust in comparison to scipy’s simplex
implementations observing both faster runtimes and better objective function
values.

###  1.1 Context and previous work

The superiorization methodology (SM) is a relatively new method that aims to
improve an existing iterative algorithm by interlacing into it perturbation
steps ([10]). In the context of optimization problems, this concept can be
applied by interlacing into a feasibility-seeking algorithm (for example a
projection method) perturbations of negative gradient steps that reduce the
objective function values. Since its development, the SM has successfully been
applied in various practical applications such as intensity-modulated radio
therapy ([3]), image reconstruction ([11]) and telecommunication networks
([17]). In the particular context of linear optimization problems (LPs), there
is a vast literature for both the feasibility problem (for example the Agmon-
Motzkin-Schönberg algorithm [1, 24],[14, algorithm 5.4.1]) and the
optimization problem (e. g., simplex algorithm, ellipsoid method [25]).
Furthermore, for the LinSup case, the “guarantee problem of the
superiorization method” has been answered positively. This problem is the
question whether superiorization can actually converge to a feasbile point
with a lower or equal objective function value than the point of the
unsuperiorized version of the same feasibility-seeking algorithm. In [13] the
authors employ the principle of concentration of measure to show the result
with high probability. A similar conclusion for the nonlinear case is yet to
be found.
This paper aims to compare the superiorization methodology with optimization
algorithms for LPs in an organized, reproducible and fair manner for problems
of varying difficulty, indicated by the problems’ dimensions and their
condition numbers. Accordingly, our work is an extension of the results in
[9], where a superiorization method for linear problems (called LinSup) was
compared to the simplex method. We augment these results by varying the
condition numbers across multiple problem instances. The condition number has
been known to have a large impact on the performance of certain methods and
often leads to the failure of an algorithm ([28]). As many real-world problems
have an inherently high condition number, immunity against ill-conditioning of
problems is a desirable property for any algorithm. Often it is enough to have
a rough understanding of the order of magnitude of the condition number and
there are efficient algorithms for its estimation ([15],[22]). There is a vast
literature on preconditioning of matrix problems, which has been reviewed in
[5].

###  1.2 Outline

We give a brief overview of the superiorization method and of condition
numbers in Sections 2.1 and 2.2, respectively. We set up multiple linear
optimization problems with varying condition numbers in Section 3.1 and run
the superiorization method and several LP algorithms on them.

For details of our implementation of the superiorization method see Section
3.2. In Section 4 we present the results, followed by a discussion of
remaining challenges and future work in Section 5.

##  2 Preliminaries

###  2.1 The superiorization methodology

The Superiorization Methodology (SM) developed from the investigation of
feasibility-seeking models of some important real-world problems such as image
reconstruction from projections and radiation therapy treatment planning.
Feasibility-seeking algorithms, mainly projection methods, generate iterative
sequences that (under reasonable conditions) converge to a point in the
feasible set. Their main advantage is that they perform projections onto the
individual sets whose intersection is the feasible set and not directly onto
the feasible set and the underlying situation is that such projections onto
the individual sets are more manageable.

When one wishes to find feasible points with a reduced, not necessarily
minimized, value of an imposed objective function then the SM comes into play.
The principle of the SM is to interlace into the iterates of a feasibility-
seeking iterative process perturbations that will steer the iterates toward
superior (meaning smaller or equal) objective function values without losing
the overall convergence of the perturbed iterates to a feasible point. To this
end “bounded perturbations” are used.

How all this is done is rigorously described in earlier papers on the SM,
consult, e.g., [10] for a recent review, read also [20]. A key feature of the
SM is that it does not aim for a constrained optimal function value, but is
content with settling for a feasible point with reduced objective function
value – reduced in comparison to the objective function value of a feasible
point that would be reached by the same feasibility-seeking algorithm without
perturbations. For many applications this is sufficient, in particular,
whenever the introduction of an objective function is only a secondary goal.
Fulfillment of the constraints, in this context, is considered by the modeler
of the real-world problem to be much more important, see, e.g., [10, 11, 3].

Many papers on the SM are cited in [8] which is a Webpage dedicated to
superiorization and perturbation resilience of algorithms that contains a
continuously updated bibliography on the subject. This Webpage is a source for
the wealth of work done in this field to date, including two special issues of
journals [12] and [18] dedicated to research of the SM. Recent work includes
[7, 19, 21, 23, 29].

###  2.2 Condition numbers

The relative condition number is a measure for the impact of deviations in the
input data on the output data of a problem. In particular, it describes the
factor of error propagation when given wrong or erroneous data. Let the
function $f\colon\mathbb{R}^{n}\to\mathbb{R}^{m}$ represent some mathematical
problem and let $x\in\mathbb{R}^{n}$ be its input, where $\mathbb{R}^{n}$
stands for the Euclidean $n$-dimensional space. Denote with
$\tilde{x}\in\mathbb{R}^{n}$ the disturbed input data. Then the relative
condition number of the problem at the point $x$ is defined as (see [31,
equation (12.4)])

$\kappa_{\text{rel}}:=\limsup_{\tilde{x}\rightarrow
x}\frac{\|f(\tilde{x})-f(x)\|}{\|\tilde{x}-x\|}\frac{\|x\|}{\|f(x)\|},$ (1)

as long as $f(x)\neq 0$. Otherwise, it is $\kappa_{\text{rel}}=\infty$. In
particular, the condition number is independent of a chosen numerical
algorithm for solving the problem $f$, but the algorithms convergence speed
may depend on the magnitude of the condition number (see [28, section 4]). In
the following we are interested in the condition number of matrices. When
$f(x)=Ax$ describes the problem of matrix multiplication, where
$A\in\mathbb{R}^{m\times n}$ and with $\|\cdot\|=\|\cdot\|_{2}$, the above
formula becomes

$\kappa_{\text{rel}}=\limsup_{\tilde{x}\rightarrow
x}\frac{\|A(\tilde{x}-x)\|}{\|\tilde{x}-x\|}\frac{\|x\|}{\|Ax\|}.$ (2)

As $f$ is differentiable and writing $\tilde{x}-x=hv$, for some unit vector
$v$ and $h=\|\tilde{x}-x\|$ we get

$\kappa_{\text{rel}}=\lim_{h\rightarrow
0}\frac{\|A(x+hv-x)\|}{h}\frac{\|x\|}{\|Ax\|}=\|A\|\frac{\|x\|}{\|Ax\|}\leq\|A\|\|A^{-1}\|.$
(3)

The term on the right is called the condition number of the matrix $A$ (see
[31, equation (12.15)]), where $A^{-1}$ denotes the inverse or, if $A$ is non-
square, the pseudo-inverse of $A$:

$\kappa(A)=\|A\|\|A^{-1}\|=\frac{\sigma_{\text{max}}}{\sigma_{\text{min}}},$
(4)

where $\sigma_{\text{max}}$ and $\sigma_{\text{min}}$ are the maximal and
minimal singular values of $A$, respectively. The condition number plays a
significant role in the analysis of numerical problems and is subject to
extensive studies in the literature ([16], [31], [15], [28]). Several methods
exist to improve high condition numbers of ill-conditioned problems (these are
the, so-called, pre-conditioning methods, see e. g., [5]) in order to increase
the accuracy of calculated solutions. This is often necessary because many
real-world applications give rise to condition numbers of significant
magnitude. This is the key motivation for the investigation in this paper.

##  3 Problem Setup and Implementation Details

###  3.1 The problem formulation

We consider linear problems of the form

$\displaystyle\\!\min_{x\in\mathbb{R}^{n}}$ $\displaystyle\ \langle
c,x\rangle$ s.t. $\displaystyle Ax\leq b$ (5) $\displaystyle\ell\leq x\leq u,$

where $A\in\mathbb{R}^{m\times
n},b\in\mathbb{R}^{m},c,\ell,u\in\mathbb{R}^{n}$. Write $A=U\Sigma V$ via the
singular value decomposition (cf. [6]) with semi-orthogonal matrices $U$ and
$V$ and diagonal matrix
$\Sigma=\operatorname{diag}(\sigma_{1},\dots,\sigma_{q})$ of singular values.
Without loss of generality, let $\sigma_{1}\geq\ \dots\ \geq\sigma_{q}\geq 0$.
The condition number $\kappa$ of $A$ is (see [16])

$\kappa(A)=\frac{\sigma_{\text{max}}}{\sigma_{\text{min}}}=\frac{\sigma_{1}}{\sigma_{q}}.$
(6)

We want to construct a sequence of matrices $A_{\kappa}$ of a specified
condition number $\kappa$ in such a way that they remain structurally similar
to each other. To this end we reverse the singular value decomposition, that
is, we create exactly one pair of $U$ and $V$ which contains the structure of
the problem and construct, for different values of $\kappa$, diagonal matrices
$\Sigma_{\kappa}$, which impose the condition number of the problem. Then, we
calculate $A_{\kappa}:=U\Sigma_{\kappa}V$. As any matrix has a singular value
decomposition, this makes it possible to define any matrix via this approach
too. In our construction we focus on matrices of full rank, because otherwise
one can remove rows or columns until full rank is achieved. In order to
construct $\Sigma$, setup $q=\min(m,n)=\operatorname{rank}(A)$ singular values
as

$\sigma_{i}=\frac{t}{z_{i}}+\frac{1-t}{s}$ (7)

where $t=\frac{\kappa-1}{q-1}$, $s=10$ and $z_{i}=\frac{si}{q}$. Note that
this choice is somewhat arbitrary. We chose this setup of the $\sigma_{i}$,
because in real-world applications the singular values often behave
approximately proportional to $\frac{1}{i}$ (instead of linear which seems
like an obvious first choice for our problem). The parameter $s$ is used to
control the magnitude of the singular values since in this model we always
have $\sigma_{q}=\frac{1}{s}$. Due to the choice of $t$ it is easy to see that

$\frac{\sigma_{1}}{\sigma_{q}}=\frac{\sigma_{\text{max}}}{\sigma_{\text{min}}}=\kappa$
(8)

as desired. For $U$ and $V$ choose random semi-orthogonal matrices and set
$A=U\cdot\operatorname{diag}(\sigma_{1},\dots,\sigma_{q})\cdot V$.
Furthermore, set $b=A\mathbbm{1}+\mathbbm{1}$, $u=-\ell=100\cdot\mathbbm{1}$
and randomly choose $c$. This choice of parameters implies the feasibility of
$x=\mathbbm{1}$. We set up this problem for multiple dimensions, ranging from
$80\times 100$ to $2000\times 2500$ (cf. [9]).

###  3.2 The superiorization algorithm

We aim to apply separately linear programming and the superiorization
methodology for the data $A=(a^{i})_{i=1}^{m},b,c,\ell,u$ that appears in (5).
In the SM we chose for the feasibility-seeking algorithm the projection method
of Agmon Motzkin and Schönberg (AMS) with overshoot $0<r:=10^{-3}$ as the
“basic algorithm”. This algorithm cyclically projects onto the individual
half-spaces $\langle a^{i},x\rangle\leq b_{i}$ via (as defined in [27, p.
411])

$\mathcal{A}_{i}(x):=\begin{cases}x-\frac{\langle a^{i},x\rangle-
b_{i}+r\|a^{i}\|}{\|a^{i}\|^{2}}a^{i}\text{, if $\langle
a^{i},x\rangle>b_{i}$}\\\ x\text{, otherwise}\end{cases}$ (9)

and a full sweep through all half-spaces is done by the algorithmic operator
$\mathcal{A}$, which is a composition of the individual projections

$\mathcal{A}:=\mathcal{A}_{m}\circ\dots\circ\mathcal{A}_{1}.$ (10)

The parameter $r$ describes, how far into the half-space the projection maps
(i. e., beyond the bounding hyperplane). For the direction vectors in the
perturbations used in the SM we chose the normalized negative gradient of the
objective function in (5), which is constant throughout and equals
$-\frac{c}{\|c\|}$. For the step-sizes $\eta_{k}$ we take an exponentially
decreasing null sequence $(10\cdot 0.99^{n})_{n\in\mathbb{N}}$ with resets
every $\tau_{\text{reset}}=20$ iterations as described in [2, p. 6]. The
starting step-size $\eta_{0}=10$ is decreased by the kernel $\alpha=0.99$ in
each iteration, unless there is a reset. In that case, set
$\eta_{k}=\eta_{0}\alpha^{\rho}$, where $\rho$ is the number of previous
resets during this run.

In other words, the $k$-th iteration consists of a gradient step
$-\eta_{k}\frac{c}{\|c\|}$ followed by a cyclic sweep of projections onto the
half-spaces via $\mathcal{A}$ as given in (9) and (10). This process was
repeated, until the iterate $x^{k}$ was feasible up to a tolerance of
$\epsilon=10^{-8}$ and the relative change
$\frac{\|x^{k}-x^{k-1}\|}{\|x^{k-1}\|}$ from the previous iterate was
negligible, i. e., became smaller than $\epsilon$. Algorithm 1 contains the
pseudocode for the described process. All runs were initialized at the all-
ones vector $x^{0}=\mathbbm{1}$. The parameters $\epsilon$, $\alpha$,
$\eta_{0}$, $r$, $\tau_{\text{reset}}$ can be adjusted for individual
preferences or a particular problem.

Data: $A$, $b$, $c$, $\ell$, $u$, $x^{0}$

Output: Superiorized vector $y$

Parameters : $\epsilon=10^{-8}$, $\alpha=0.99$, $\eta_{0}=10$, $r=10^{-3}$,
$\tau_{\text{reset}}=20$

1 $k\leftarrow 0$;

2 $\tau\leftarrow 0$;

3 $\rho\leftarrow 0$;

4 $x^{-1}\leftarrow x^{0}+\mathbbm{1}$;

5 while  _$\max_{i}(\langle a^{i},x^{k}\rangle-b_{i})\geq\epsilon$ or
$\frac{\|x^{k}-x^{k-1}\|}{\|x^{k-1}\|}\geq\epsilon$_ do

6 $x^{k+1}\leftarrow x^{k}-\eta_{k}\frac{c}{\|c\|}$;

7 $x^{k+1}\leftarrow\mathcal{A}(x^{k+1})$;

8 $x^{k+1}\leftarrow\max(x^{k+1},\ell)$;

9 $x^{k+1}\leftarrow\min(x^{k+1},u)$;

10

11 $\tau\leftarrow\tau+1$;

12 if  _$\tau <\tau_{\text{reset}}$_ then

13 $\eta_{k+1}\leftarrow\eta_{k}\alpha$;

14

15 else

16 $\tau\leftarrow 0$;

17 $\rho\leftarrow\rho+1$;

18 $\eta_{k+1}\leftarrow\eta_{0}\alpha^{\rho}$;

19

20 end if

21

22 $k\leftarrow k+1$;

23

24 end while

25$y\leftarrow x^{k}$;

return $y$

Algorithm 1 The superiorization algorithm

###  3.3 The LP algorithms

We compare the superiorization method with the following LP algorithms:

  1. 1.

scipy.simplex

  2. 2.

scipy.revised simplex

  3. 3.

scipy.interior-point

  4. 4.

gurobi.primal simplex.

Scipy is a library for scientific computing in the programming language
Python. It is freely available at www.scipy.org. Its optimization suite
scipy.optimize contains multiple implementations of common optimization
algorithms like the SQP method, the dogleg method or the conjugate gradient
method. Since we are working with the data of an LP, we employed specialized
LP methods only, in particular the simplex and revised simplex method as well
as the interior-point method of scipy’s linprog function.
Gurobi is a commercial state-of-the-art solver for linear and nonlinear,
continuous and (mixed) integer problems. It is available at www.gurobi.com and
offers a wide range of customization, including multiple algorithms,
globalization strategies, preconditioning, etc. For our experiments we turned
off Gurobi’s automated choice of algorithm and instead forced it to use its
implementation of the primal simplex to ensure consistency of the output data.

##  4 Experimental Results of the Comparative Investigation

Our main finding is, that, based on our computational experiments, the
superiorization methodology is quite robust with respect to increasing
condition numbers. In particular, the runtime remains relatively small when
compared to the scipy implementations of the simplex algorithm. We conjecture
that this is due to the fact, that our superiorization implementation never
considers the full problem at once, but performs individual projections onto
the half-spaces instead. This comes at the cost of sacrificing feasibility
during early iterations. Furthermore, we conjecture, that the bounded
perturbation resilience of the basic algorithm may play a role in absorbing
errors that occur during computation, which would normally be amplified by the
condition number.

Depending on the rate at which the step-sizes of the perturbations converge to
zero, this would mean that superiorization, in general, may be less affected
by high condition numbers, than other current algorithms. A secondary finding
is, that superiorization can produce better results than the scipy simplex
implementations, when both are terminated at a certain time, that is before
the usual termination criterion is met, cf. Figure 1(c). Furthermore, we see
that superiorization reaches its termination criterion considerably faster
than scipy in higher dimensions.

The trend is clear: While in problems with smaller dimensions the simplex
algorithms are considerably faster, with increasing problem dimensions the
scipy implementations take much longer to terminate. We also see a clear
difference in the algorithms styles: Simplex aims to become feasible first and
then starts to improve the objective function. Superiorization, on the other
hand, reduces objective function values in its initial iterations because then
the step-sizes of the perturbations are still large and only as iterations
proceed the effect of feasibility-seeking becomes stronger.

Finally, it should be noted that scipy’s revised simplex performs better than
scipy’s standard simplex method. This can be attributed to the internal $LU$
decomposition of the basic matrices, leading to better numerical stability as
described in [4].

(a) $80\times 100$

(b) $200\times 250$

(c) $400\times 500$. Had both implementations been terminated at an earlier
time (say after 25 seconds), the superiorized solution would be both less
violating the constraints and closer to optimality. As the problem dimensions
increase this trend becomes more and more prominent

(d) $800\times 1000$

Figure 1: Figures 1(a)-1(d) show the behavior of the different algorithms on
problems with increasing dimension for fixed condition number $\kappa=1000$.
On the left we see the maximum violation of the constraints $\max_{i}(\langle
a^{i},x\rangle-b_{i})$ plotted against the runtime. On the right we have the
corresponding objective function values $\langle c,x\rangle$ plotted against
the runtime.

(a) $400\times 500$

(b) $800\times 1000$

(c) $2000\times 2500$

Figure 2: Figures 2(a)-2(c) demonstrate well the robustness with respect to
increasing condition numbers. On the left we see the runtimes of each
algorithm, plotted against the condition numbers, on the right we see the
objective function value at termination, plotted against the condition
numbers. Notice the severely suboptimal objective function value of scipy’s
simplex for high condition numbers. This explains the low runtime as the
algorithm realizes that these problems are hard and quickly “surrenders”. This
trend continues in Figures 2(b) and 2(c) and the revised simplex too starts to
reach its limits. Superiorization on the other hand proves to be quite stable
in terms of its runtime with regards to increasing condition numbers.

We observe that in our experiments the interior point method and the Gurobi
implementation outperformed the other algorithms by a lot. For the interior
point method this matched our expectations and can be explained in the
following way: The condition number of a matrix can be interpreted as a
measure of how linearly dependent its rows or columns are. A well-conditioned
matrix (i. e., $\kappa=1$) only has a single singular value and will be semi-
orthogonal, whereas an ill-conditioned matrix ($\kappa$ “large”) will have
“almost” linearly dependent entries. Consequently, the half-spaces
$\\{x\in\mathbb{R}^{n}\colon\langle a^{i},x\rangle\leq b_{i}\\}$ will be
almost parallel and the resulting polyhedron will consist of many facets and
vertices. A basic simplex implementation, which moves from vertex to vertex,
will consequently face a long runtime. The interior point method, on the other
hand, is not dependent on the vertices. It will take its path through the
interior of the polyhedron regardless of its boundary.
The Gurobi implementation is based on observations of the condition number. In
that sense, our experiments played into one of Gurobis strongest suits. As one
of the biggest commercial state of the art solvers we did not expect to
outperform Gurobi with our implementation, but rather used it as a benchmark
to compare with the other algorithms.

##  5 Conclusions

In this paper we experimentally discussed the superiorization method and
constrained optimization on a set of examplary linear problems with varying
condition numbers and problem sizes with the aim of investigating and
comparing their immunity to increasing condition numbers. Our experimental
results are promising for the observed problem sizes, but more work needs to
be done for larger instances to verify that the trend that we observed
continues.

The superiorization method and constrained optimization use the same input
data which consists of a family of constraints obtained from the modeling
process along with a user-chosen objective function. But the two approaches
aim at different end-points of their iterative processes. The easy
accessibility of the superiorization methodology allows for quick
implementations with the advance knowledge that the aim is not to reach a
constrained optimum. At the same time the SM may compute its solutions at a
lower runtime (in the case of scipy’s simplex and revised simplex).

As superiorization is a relatively new concept, we expect that, with further
tuning of its parameters it will continue to find a place as a computational
model and tool in situations in which users do not wish to invest efforts in
seeking a constrained optimal point but rather wish to find a feasible point
which is “superior” in the sense of having a smaller or equal objective
function value than that of a feasible point reached by the same feasibility-
seeking algorithm.

Another interesting point is that, as is well-known, interior point methods
reach their performance limits for ill-conditioned nonlinear problems (e. g.,
in intensity-modulated radio therapy). It would be interesting to compare the
superiorization methodology in a nonlinear setting, with a different basic
feasibility-seeking algorithm, to the interior point method to see if
superiorization can contribute to solving these problems faster.

Acknowledgements: The work of Y.C. is supported by U.S. National Institutes of
Health Grant Number R01CA266467 and by the Cooperation Program in Cancer
Research of the German Cancer Research Center (DKFZ) and Israel’s Ministry of
Innovation, Science and Technology (MOST).

## References

  * [1] S. Agmon, The relaxation method for linear inequalities, Canadian Journal of Mathematics, 6 (1954), pp. 382–392.
  * [2] F. J. Aragón-Artacho, Y. Censor, A. Gibali, and D. Torregrosa-Belén, The superiorization method with restarted perturbations for split minimization problems with an application to radiotherapy treatment planning, Applied Mathematics and Computation, 440 (2023).  Article 127627.
  * [3] F. Barkmann, Y. Censor, and N. Wahl, Superiorization as a novel strategy for linearly constrained inverse radiotherapy treatment planning, Frontiers in Oncology, 13 (2023), pp. 1–15.
  * [4] R. H. Bartels, A stabilization of the simplex method, Numerische Mathematik, 16 (1971), pp. 414–434.
  * [5] M. Benzi, Preconditioning techniques for large linear systems: A survey, Journal of Computational Physics, 182 (2002), pp. 418–477.
  * [6] J. Bisgard, Analysis and Linear Algebra: The Singular Value Decomposition and Applications, The American Mathematical Society, 2021.
  * [7] E. Bonacker, A. Gibali, K.-H. Küfer, and P. Süss, Speedup of lexicographic optimization by superiorization and its applications to cancer radiotherapy treatment, Inverse Problems, 33 (2017).  Article 044012.
  * [8] Y. Censor, Superiorization and perturbation resilience of algorithms: A bibliography compiled and continuously updated.  Available at: http://math.haifa.ac.il/yair/bib-superiorization-censor.html.
  * [9] Y. Censor, Can linear superiorization be useful for linear optimization problems?, Inverse Problems, 33 (2017).  Article 044006.
  * [10] Y. Censor, Superiorization: The asymmetric roles of feasibility-seeking and objective function reduction, Applied Set-Valued Analysis and Optimization, 5 (2023), pp. 325–346.
  * [11] Y. Censor, R. Davidi, G. T. Herman, R. W. Schulte, and L. Tetruashvili, Projected subgradient minimization versus superiorization, Journal of Optimization Theory and Applications, 160 (2014), pp. 730–747.
  * [12] Y. Censor, G. T. Herman, and M. Jiang, Superiorization: Theory and applications, Inverse Problems, 33 (2017).  Article 040301.
  * [13] Y. Censor and E. Levy, An analysis of the superiorization method via the principle of concentration of measure, Applied Mathematics & Optimization, 83 (2021), pp. 2273–2301.
  * [14] Y. Censor and S. A. Zenios, Parallel Optimization: Theory, Algorithms, and Applications, Oxford University Press, New York, NY, USA, 1997\.
  * [15] A. K. Cline, C. B. Moler, G. W. Stewart, and J. H. Wilkinson, An estimate for the condition number of a matrix, SIAM Journal on Numerical Analysis, 16 (1979), pp. 368–375.
  * [16] P. Deuflhard and A. Hohmann, Numerische Mathematik 1, De Gruyter, Berlin, Germany, 2008.
  * [17] J. Fink, Fixed Point Algorithms and Superiorization in Communication Systems, PhD thesis, Fakultät IV-Elektrotechnik und Informatik, Technische Universität Berlin, Germany, 2022.
  * [18] A. Gibali, G. T. Herman, and C. Schnörr, Superiorization versus constrained optimization: Analysis and applications, Journal of Applied and Numerical Optimization, 2 (2020), pp. 1–2.
  * [19] M. Guenter, S. Collins, A. Ogilvy, W. Hare, and A. Jirasek, Superiorization versus regularization: A comparison of algorithms for solving image reconstruction problems with applications in computed tomography, Medical Physics, 49 (2022), pp. 1065–1082.
  * [20] G. T. Herman, Problem structures in the theory and practice of superiorization, Journal of Applied and Numerical Optimization, 2 (2020), pp. 71–76.
  * [21] G. T. Herman and R. Davidi, Image reconstruction from a small number of projections, Inverse Problems, 24 (2008).  Article 045011.
  * [22] N. J. Higham, A survey of condition number estimation for triangular matrices, SIAM Review, 29 (1987), pp. 575–596.
  * [23] V. Martín-Márquez, S. Reich, and S. Sabach, Iterative methods for approximating fixed points of Bregman nonexpansive operators, Discrete and Continuous Dynamical Systems - Series S, 6 (2012), pp. 1043–1063.
  * [24] T. S. Motzkin and I. J. Schönberg, The relaxation method for linear inequalities, Canadian Journal of MAthematics, 6 (1954), pp. 393–404.
  * [25] K. G. Murty, Linear Programming, Wiley & Sons Ltd., Hoboken, NJ, USA, 1991.
  * [26] J. Nocedal and S. J. Wright, Numerical Optimization, Springer, New York, NY, USA, 1999.
  * [27] B. T. Polyak, Random algorithms for solving convex inequalities, Studies in Computational Mathematics, 8 (2001), pp. 409–422.
  * [28] A. Pyzara, B. Bylina, and J. Bylina, The influence of a matrix condition number on iterative methods’ convergence, 2011 Federated Conference on Computer Science and Information Systems, (2011), pp. 459–464.
  * [29] M. J. Schrapp and G. T. Herman, Data fusion in x-ray computed tomography using a superiorization approach, Review of Scientific Instruments, 85 (2014).  Article 053701.
  * [30] B. Schultze, Y. Censor, P. Karbasi, K. E. Schubert, and R. W. Schulte, An improved method of total variation superiorization applied to reconstruction in proton computed tomography, IEEE Transactions on Medical Imaging, 39 (2018), pp. 294–307.
  * [31] L. N. Trefethen and D. Bau, Numerical Linear Algebra, Twenty-fifth Anniversary Edition, SIAM, Philadelphia, PA, USA, 2022.