# Model-Centric Volumetric Point Cloud Attributes

Ricardo L. de Queiroz, Camilo Dorea, Davi R. Freitas, Maja Krivokuća, and
Gustavo P. Sandri  R. L. de Queiroz is with the Computer Science Department at
Universidade de Brasilia, Brasilia, Brazil, e-mail<EMAIL_ADDRESS>Dorea
is with the Computer Science Department at Universidade de Brasilia, Brasilia,
Brazil, e-mail<EMAIL_ADDRESS>Freitas is with the Post-Graduate
Program at the Electrical Engineering Science Department at Universidade de
Brasilia, Brasilia, Brazil, e-mail<EMAIL_ADDRESS>Krivokuća is
with INRIA, Rennes, France, e-mail<EMAIL_ADDRESS>Sandri is with the
Federal Institute of Brasilia, Brasilia, Brazil, e-mail:
<EMAIL_ADDRESS>
###### Abstract

Point clouds have recently gained interest, especially for real-time
applications and for 3D-scanned material, such as is used in autonomous
driving, architecture, and engineering, to model real estate for renovation or
display. Point clouds are associated with geometry information and attributes
such as color. Be the color unique or direction-dependent (in the case of
plenoptic point clouds), it reflects the colors observed by cameras displaced
around the object. Hence, not only are the viewing references assumed, but the
illumination spectrum and illumination geometry is also implicit. We propose a
model-centric description of the 3D object, that is independent of the
illumination and of the position of the cameras. We want to be able to
describe the objects themselves such that, at a later stage, the rendering of
the model may decide where to place illumination, from which it may calculate
the image viewed by a given camera. We want to be able to describe transparent
or translucid objects, mirrors, fishbowls, fog and smoke. Volumetric clouds
may allow us to describe the air, however “empty”, and introduce air
particles, in a manner independent of the viewer position. For that, we rely
on some eletromagnetic properties to arrive at seven attributes per voxel that
would describe the material and its color or transparency. Three attributes
are for the transmissivity of each color, three are for the attenuation of
each color, and another attribute is for diffuseness. These attributes give
information about the object to the renderer, with whom lies the decision on
how to render and depict each object.

##  I Introduction

Point clouds are becoming more popular today for representing 3D objects and
people in scenarios that are captured from the real world, without the
additional complexity of deriving meshes or other connected (surface)
structures [1]. Points in a cloud are by definition unconnected, and are
usually described by a list of points describing the geometry of each point
(i.e. the usual $xyz$ coordinates in 3D space) and any attribute the point
might have. Examples of attributes are colors, in any color space, reflection
intensity, or surface normals. There are voxelized point clouds, wherein the
space is subdivided into a regular grid of “volumetric pixels” or voxels. One
may convert back and forth from voxelized and non-voxelized point clouds, by
snapping points to the voxel grid and by creating a point at the center of a
voxel. Most point clouds, even those that are voxelized, are in fact surface
point clouds, since they do not represent the interior of objects, just the
visible surfaces. Another alternative is a volumetric point cloud, in which
all the voxels inside and on the surface of the 3D object are described, with
information on which voxels are occupied or not, perhaps using functions that
cover the entire space. Fig. 1 illustrates the volumetric point cloud case.

In computer graphics [2], describing the model and rendering it are completely
different processes. The model describes the object, its location and its
properties. The rendering is done via ray-tracing or another method, assuming
a given model, the illumination, and the camera viewpoint. Point clouds, on
the other hand, usually mix both concepts, as the voxel color is typically
captured by cameras in given positions, subject to a given illumination. We
propose a point cloud attribute format that is model-centric, illumination-
independent, and viewpoint-independent, which would describe the objects
themselves independently of the rendering.

The proposed model-centric concept would provide light independence in
describing objects, for example buildings [3]. The same point cloud could be
rendered in daytime, night-time, or even include interior lights. Transparent
materials could be represented, as well as reflective materials (e.g., brass),
or even mirrors. We anticipate that scanning transparent objects may be a
great challenge. Volumetric objects would need to be represented in this case,
instead of just the object surfaces. Interior materials could also be
represented, since visibility is not implied in the proposed point cloud. We
have developed a set of seven attributes that we believe can be used to
reasonably describe the object materials, in terms of their visibility, in a
voxelized volumetric point cloud. Their reasoning is described next.

Figure 1: 3D objects can be represented by a cloud of points or a grid of
voxels, encompassing occupied or empty voxels. Conceptually, all voxels may
receive attributes.

##  II Eletromagnetism-based motivation

Just a quick word on eletromagnetism, pertaining to transversal eletromagnetic
waves in a given medium [4]. Light in most cases can be modeled as such and it
is of interest to understand what happens to the wave as it changes media.

###  II-A Planar waves

The wave propagates due to its varying electric and magnetic fields. From
Maxwell’s equations, variations in the magnetic field induce changes in the
electric field and vice-versa. The result is a wave that propagates in space.
In order to simplify the description, assume a planar wave propagating in one
direction, with fields polarized along the perpendicular directions. Let ${\bf
a}_{x}$, ${\bf a}_{y}$ and ${\bf a}_{z}$ be vectors aligned to the coordinate
system and assume the wave is propagating along ${\bf a}_{z}$ and the electric
and magnetic fields are given by

${\bf E}=E_{x}e^{-\gamma z}{\bf a}_{x},$ (1) ${\bf
H}=\frac{E_{x}}{\eta}e^{-\gamma z}{\bf a}_{y}.$ (2)

where $E_{x}$ is a time-varying component of the electric field over the $x$
axis and $\eta$ is the impedance of the medium. This is a linearly polarized
wave. Both $\eta$ and $\gamma$ are potentially complex numbers. If we express

$\gamma=\alpha+j\beta,$ (3)

so that

${\bf E}=E_{x}e^{-\alpha z}e^{-j\beta z}{\bf a}_{x},$ (4)

which tells us that $\alpha$ is responsible for attenuating the wave while
$\beta$ is its phase propagation term. The same occurs for the magnetic field.

The wave is assumed as having a frequency $\omega=2\pi f=2\pi/\lambda$, where
$\omega$ is given in radians per second, $f$ is the frequency in Hertz and
$\lambda$ is its wavelength.

###  II-B Light in a medium

Before we describe the quantities $\gamma$ and $\eta$ in the previous section,
keep in mind that a medium is associated with three eletromagnetic quantities.
There is the electric permittivity
$\varepsilon=\varepsilon_{r}\varepsilon_{0}$, where $\varepsilon_{r}$ is the
relative permittivity (or dielectric constant) of the medium and
$\varepsilon_{0}$ is the permittivity of the vacuum. There is also the
magnetic permeability $\mu=\mu_{r}\mu_{0}$, where $\mu_{r}$ is the relative
permeability of the medium and $\mu_{0}$ is the permeability of the vacuum.
$\sigma$ is the conductivity of the medium.

The wave propagation term is then

$\gamma=\sqrt{j\omega\mu(\sigma+j\omega\varepsilon)},$ (5)

i.e.

$\alpha=\omega\sqrt{\mu\varepsilon}\sqrt{\frac{1}{2}\left[\sqrt{1+\left(\frac{\sigma}{\omega\varepsilon}\right)^{2}}-1\right]},$
(6)

$\beta=\omega\sqrt{\mu\varepsilon}\sqrt{\frac{1}{2}\left[\sqrt{1+\left(\frac{\sigma}{\omega\varepsilon}\right)^{2}}+1\right]}.$
(7)

The medium impedance is

$\eta=\sqrt{\frac{j\omega\mu}{\sigma+j\omega\varepsilon}}.$ (8)

The wave travels in such a medium with a phase velocity (which, in this case,
is a pretty good approximation to the true velocity) given by

$v=\frac{1}{\sqrt{\frac{\mu\varepsilon}{2}\left[1+\sqrt{1+\left(\frac{\sigma}{\omega\varepsilon}\right)^{2}}\right]}}$
(9)

###  II-C Cases of interest

The above equations impose some algebraic difficulties, but we are most
interested in two cases: the good dielectric and the good conductor. For the
good dielectric, we assume $\sigma$ to be very small, such that $\alpha\approx
0$, $\beta\approx\omega\sqrt{\mu\varepsilon}$, and the impedance is:

$\eta\approx\sqrt{\frac{\mu}{\varepsilon}},$ (10)

with equalities for a perfect dielectric material ($\sigma=0$). For the same
conditions, the speed in the medium becomes

$v\approx\frac{1}{\sqrt{\mu\varepsilon}}.$ (11)

A special case of a good dielectric is the open air, which approximates the
vacuum, wherein $\mu_{r}=\varepsilon_{r}=1$, such that
$\eta=\sqrt{\mu_{0}/\varepsilon_{0}}\approx 377\ \Omega$. The speed is
$c=1/\sqrt{\mu_{0}\varepsilon_{0}}$, which is the causality speed of the
universe or the “speed of light”, and it is the maximum speed any wave can
travel.

In the case of a good conductor and when $\sigma$ is large,

$\eta\approx\sqrt{\frac{\omega\mu}{\sigma}}e^{j\pi/4}\ \ \mbox{and}\ \
v\approx\sqrt{\frac{2\omega}{\mu\sigma}}.$ (12)

One last important note is that when there is conductivity, the propagation is
attenuated. Any imaginary part of the impedance term causes attenuation. A
good conductor, thus, does not allow propagation in its core. In fact the
eletromagnetic wave barely penetrates a conductor, as much as what is named
the skin depth, which is the point where the electric field has dropped by a
factor of $e^{-1}$. The skin depth is $\delta=v/\omega$. For 100 MHz and
copper, for example, $\delta=6.6\mu m$.

###  II-D Light crossing a media interface

Assume, then, light as a planar wave crossing an infinite planar interface
(perpendicular to normal vector $\bf n$) in between materials 1 and 2. Assume
the origin to be where the wave hits the planar surface. Assume that the
incident light comes from direction ${\bf n}_{i}$. It is well known that some
light is reflected at direction ${\bf n}_{r}$ while some is transmitted
(refracted) into medium 2 at direction ${\bf n}_{t}$. All four vectors lie in
the plane defined by vectors ${\bf n}_{i}$ and ${\bf n}$. If all vectors have
unit norm, ${\bf n}_{r}=2({\bf n}_{i}\cdot{\bf n})\ {\bf n}-{\bf n}_{i}$. Note
that $\theta_{1}={\bf n}_{i}\cdot{\bf n}={\bf n}_{r}\cdot{\bf n}$ is the
incidence (and reflectance) angle, while the transmission (refraction) angle
$\theta_{2}$ is given by Snell’s law:

$\frac{\sin\theta_{2}}{\sin\theta_{1}}=\frac{v_{2}}{v_{1}},$ (13)

where the index indicates the medium. If both media are good dielectrics, we
have

$\frac{\sin\theta_{2}}{\sin\theta_{1}}=\frac{\sqrt{\mu_{1}\varepsilon_{1}}}{\sqrt{\mu_{2}\varepsilon_{2}}}=\frac{\beta_{1}}{\beta_{2}}.$
(14)

Snell’s equation is sometimes expressed in terms of the refraction index which
is $c/v$. Such a redundant quantity was left out of our discussion.

The refraction (transmission) vector is:

${\bf n}_{t}=(\sin\theta_{2}\cos\theta_{1}-\cos\theta_{2}){\bf n}-{\bf
s}_{i}.$ (15)

Figure 2: Wave incidence geometry on a planar interface. Top: $E$
perpendicular to the surface; bottom: $E$ parallel to the surface.

There are three waves to consider: the incident (fields ${\bf E}_{i}$ and
${\bf H}_{i}$), the reflected (fields ${\bf E}_{r}$ and ${\bf H}_{r}$), and
the transmitted or refracted (fields ${\bf E}_{t}$ and ${\bf H}_{t}$), which
are depicted in Fig. 2. Since the magnitude of the electric and magnetic
fields are related by the impedance, we are concerned with the propagation
ratio of the electric field into medium 2, i.e., $T=E_{t}/E_{i}$ and the
amount of the field that is reflected $\Gamma=E_{t}/E_{i}$. As also depicted
in Fig. 2, there are two cases to consider, and a more general case which is a
combination of both: when the electric field ${\bf E}_{i}$ is parallel to the
surface (${\bf n}\cdot{\bf E}_{i}=0$)

$\Gamma_{\perp}=\frac{E_{r}}{E_{i}}=\frac{\eta_{2}\cos\theta_{1}-\eta_{1}\cos\theta_{2}}{\eta_{2}\cos\theta_{1}+\eta_{1}\cos\theta_{2}},$
(16)

$T_{\perp}=\frac{E_{t}}{E_{i}}=\frac{2\
\eta_{2}\cos\theta_{1}}{\eta_{2}\cos\theta_{1}+\eta_{1}\cos\theta_{2}},$ (17)

and when the magnetic one is parallel to the surface (${\bf n}\cdot{\bf
H}_{i}=0$)

$\Gamma_{\parallel}=\frac{E_{r}}{E_{i}}=\frac{\eta_{1}\cos\theta_{1}-\eta_{2}\cos\theta_{2}}{\eta_{1}\cos\theta_{1}+\eta_{2}\cos\theta_{2}},$
(18)

$T_{\parallel}=\frac{E_{t}}{E_{i}}=\frac{2\
\eta_{2}\cos\theta_{1}}{\eta_{1}\cos\theta_{1}+\eta_{2}\cos\theta_{2}}.$ (19)

###  II-E Considerations

In an imperfect dielectric, a small but non-zero $\sigma$ may cause a small
but significant attenuation. Attenuation, in general, is caused by
“absorption”, which is dissipation of the energy in the material, being
converted into heat. The reason is that conductivity causes minor currents
which have dissipation losses due to electrical resistance.

As conductivity increases, the refraction angle decreases. For a good
conductor, the field rapidly attenuates but $\theta_{2}\approx 0$, i.e., the
transmitted wave is refracted along the normal no matter where the incident
wave comes from.

For non-ferromagnetic materials, $\mu_{r}\approx 1$ and, for perfect
dielectric materials,

$\eta\approx\sqrt{\frac{\mu_{0}}{\varepsilon_{r}\varepsilon_{0}}}=\frac{1}{\sqrt{\varepsilon_{r}}}\eta_{0}$
(20)

and

$v\approx\sqrt{\frac{1}{\varepsilon_{r}\varepsilon_{0}\mu_{0}}}=\frac{1}{\sqrt{\varepsilon_{r}}}c$
(21)

where $\eta_{0}$ is the impedance of vacuum and $c$, again, is the causality
speed or the light speed in the vacuum, the maximum speed there is. In this
case, $1/\sqrt{\varepsilon_{r}}$ modulates both speed and impedance, which we
use to calculate the refraction angle ($\theta_{2}$) and the transmission
($T$) and reflection ($\Gamma$) coefficients.

In this case, as the relative permittivity $\varepsilon_{r}$ drops to 1, we
obtain the impedance and speed of vacuum. Conversely, for very large
permittivity, both impedance and speed drop. In the limit, if medium 2 has
such large permittivity and medium 1 does not, $\Gamma\rightarrow-1$ and
$T\rightarrow 0$, i.e., there is total reflection.

###  II-F Light Beam Diffusion

The assumption so far was that of a planar wave incident onto a planar surface
separating two uniform media. This is never completely true. The rugosity of a
surface, even at a micron level, may make a light beam appear as having
multiple incident angles since the wavelength of visible light is in the order
of 0.5 micron. This could result in dispersion of both the transmitted and
reflected beams. This is the reason why conductor materials are not
necessarily good mirrors. Commercial mirrors are usually made of a conducting
material, first in liquid form, pressed within smooth plates, one of then
usually made of glass.

Another source of light dispersion is the layering of materials, like coatings
on surfaces, which cause many reflections and refractions from stacked rough
planes to be combined, thus forming diffuse beams.

Another beam deformity can be caused by heterogeneous materials. For example,
a perfect dielectric may be polluted with conductive or high-impedance
materials, in a mix of attenuation, reflection and refraction. One example is
sea water, which can be clean and transparent, or murky and dark, depending on
the sediments in suspension. Similarly, mist or smoke can dampen light
propagation in air.

##  III Attributes

In our use case of a voxelized point cloud, surface refraction and reflection
occurs at voxel boundaries. A voxel represents a uniform medium, but its
boundaries are too large compared to visible light wavelength. If a voxel
represents a cubic centimeter of space, for example, all material
imperfections and rugosity at the millimeter down to micron scale may cause
dispersion and diffusion.

We need measures that are simple, normalized from 0 to 1, for example, without
reverting to difficult properties such as permittivity or conductivity. Yet,
we need to be able to calculate refraction angles, beam dispersal and
transparency.

We want to describe on a per-voxel basis: how much light is reflected and how
it refracts; how much light is attenuated within the medium; how light is
dispersed when crossing rough surfaces at sub-voxel scales. We propose to use
3 parameters, as described below.

Transmissivity - This parameter $P_{t}$ controls how much light is absorbed
into the voxel, and is related to the inverse of the medium eletromagnetic
impedance. It ranges from 0.0 to 1.0. When $P_{t}=0$, we mean that light is
completely propagated into the voxel. Conversely, when $P_{t}=1$, the light is
mostly reflected by the voxel.

One rendering suggestion is to associate $P_{t}$ with $\varepsilon_{r}$, for
example, setting $P_{t}=k_{t}\log(\varepsilon_{r})$ for $k_{t}$ some constant
of comfort. Since $\varepsilon_{r}\geq 1$, the idea is that $P_{t}=0.0$ would
represent the permittivity of the vacuum, while a sufficiently large
permitivitty can be made to correspond to $P_{t}=1.0$. Further biasing of the
correspondence can be made using ”gamma” correction, as depicted in Fig. 3,
where the typical “gamma” parameter is not related to the parameter with the
same symbol from the previous Section. In the junction between two voxels,
from $P_{1}$ one may find $\varepsilon_{r}$, more precisely
$1/\sqrt{\varepsilon_{r}}$ of each medium, from where we can refer to Eq. (21)
in order to estimate velocities, and use Snell’s equation (13) to find
refraction angles. From the $1/\sqrt{\varepsilon_{r}}$ of each medium, one can
also calculate the media impedance using Eq. (20), from which the intensity of
light reflected of refracted can be calculated form Fresnel equations in Eqs.
(16) through (19).

Figure 3: Gamma correction for parameter assignment.

Attenuation - This parameter $P_{a}$ could be useful to express the
attenuation of the medium. It may be related to the conductivity of the medium
$\sigma$, but it may also indicate how much of the incoming (refracted) light
is to be propagated to the next voxel. For example, the quantity $P_{a}=0.0$
may indicate a conductor with (near) total attenuation. Recall that the “skin
depth” of a conductor at visible light frequency (in the order of $10^{15}$
rad/s) is incredibly small compared to a voxel size. $P_{a}=1.0$ may indicate
that all incoming light to a voxel is passed to the next one.

Diffuseness - This parameter $D$ controls the diffusion of the incident light
caused by the non-uniformity of the interface at sub-voxel dimensions. Since
it is caused by the rugosity of the interface, the dispersion affects both the
reflection and the transmission (refraction). In a sense, it is a measure of
specularity of the voxel surface. If $D=0$, we assume there is no dispersion
and light is reflected (refracted) with perfect specularity (transmission). If
$D=1$, the material is assumed perfectly Lambertian and light is reflected
(refracted) equally in all directions. Intermediary values of $D$ assume a
progressive increase in dispersion of the main beam and increasing Lambertian
properties. This is illustrated in Fig. 4. Rendering for those intermediary
values can use, for example, the dispersion of the Phong [5] or Blinn-Phong
[6] models. The Schlick model [7] assumes multiple layers, and that may be
accounted for by imposing diffusion over multiple layers of voxels.

Figure 4: Diffusion parameter relfection examples. When $D=0$, there is
perfect specularity and light is reflected without dispersion. As $D$
increases, reflected light becomes more diffuse with some Lambertian
reflection. For maximum diffusion ($D=1$), reflection is completely Lambertian
with no preferred direction regardless of the incidence geometry. Similar
concepts apply to transmitted (refracted) light.

##  IV Colors and transparency

Both the impedance and the propagation speed are dependent on $\omega$, the
wave frequency. For a very small conductivity, the wave across two dielectric
materials is much more dependent on the permittivity than on the frequency.
The permittivity, however, may also depend on the frequency. Rather than
defining $\varepsilon$ for a spectrum of frequencies, since colors are
commonly described in terms of red, green and blue components, we use these
three color components. They act like filtered averages of the spectrum, or,
in some sense, they are some kind of sampling of the spectrum.

Thus, in order to convey the colors of the materials, as a reflection or
refraction of a white incident color (flat spectrum, or equal red-green-blue
components), we propose to convey the parameters $P_{a}$ and $P_{t}$ for each
component (red, green and blue). It is reasonable to assume that the rugosity
which causes diffusion has dimensions larger than the wavelengths of either
red, green or blue components. Hence, $D$ may not need to be represented for
each color component.

As a result, we propose to use seven parameters as attributes to describe the
materials for each voxel of a volumetric point cloud. The parameters are:
$R_{t}$, $G_{t}$, $B_{t}$, $R_{a}$, $G_{a}$, $B_{a}$, and $D$. The first three
are the $P_{t}$ (transmissivity) for each color component, with similar
notation for the next three in terms of $P_{a}$ (attenuation). In summary, we
propose the attributes in Table I.

TABLE I: Proposed attributes per voxel. | Min. | Meaning | Max. | Meaning
---|---|---|---|---
$R_{t}$ | 0.0 | All red transmitted | 1.0 | All red reflected
$G_{t}$ | 0.0 | All green transmitted | 1.0 | All green reflected
$B_{t}$ | 0.0 | All blue transmitted | 1.0 | All blue reflected
$R_{a}$ | 0.0 | No red attenuation | 1.0 | All red absorbed
$G_{a}$ | 0.0 | No green attenuation | 1.0 | All green absorbed
$B_{a}$ | 0.0 | No blue attenuation | 1.0 | All blue absorbed
$D$ | 0.0 | Perfect specularity | 1.0 | Lambertian diffusion

In Table II we list examples of attributes for a number of materials of
interest. All the parameter numbers are guessed examples just to illustrate
how the 7 parameters may describe the appearance of the diverse materials.
These numbers are not of any quantitative interest, just for qualitative
illustration. We only used simplifying numbers 0.0, 0.2, 0.5, 0.8, and 1.0 to
represent “nothing”, “a little”, “medium”, “a lot”, and “all”. We used shirts
as examples of a Lambertian surface that reflects the light in different
amounts to give the color appearance of the shirt. Note that all shirts have
full attenuation (they are not transparent) and maximum diffusion. The more
they reflect, the lighter they appear. Brass is conductive with large
attenuation. Brass objects are usually shiny objects, i.e., reflective
depending on how polished they are. We assumed brass to be somewhat polished
and with a small but significant diffusion. Glass is a transparent object. We
assumed it not to be perfectly transparent, with $P_{t}=0.2$. We also assumed
a small attenuation and a very flat glass without diffusion. We included a
large diffusion in frosted glass. Similarly for water, where we changed the
values to include more “blue” reflection. Air (clear air) has
$P_{t}=P_{a}=D=0$, i.e., no reflection, no attenuation, no diffusion. Smoke or
mist may be represented by increasing a little some of those parameters. In
any case, these parameters are only illustrative, with no empirical
evaluation.

TABLE II: Example parameters for a few materials. Material | $R_{t}$ | $G_{t}$ | $B_{t}$ | $R_{a}$ | $G_{a}$ | $B_{a}$ | $D$
---|---|---|---|---|---|---|---
White shirt | 0.8 | 0.8 | 0.8 | 1 | 1 | 1 | 1
Dark shirt | 0.2 | 0.2 | 0.2 | 1 | 1 | 1 | 1
Red shirt | 0.8 | 0 | 0 | 1 | 1 | 1 | 1
Green shirt | 0 | 0.8 | 0 | 1 | 1 | 1 | 1
Blue shirt | 0 | 0 | 0.8 | 1 | 1 | 1 | 1
Color shirt | 0.8 | 0.5 | 0.2 | 1 | 1 | 1 | 1
Skin | 0.5 | 0.5 | 0.2 | 0.8 | 1 | 1 | 0.8
Brass | 0.8 | 0.8 | 0.2 | 1 | 1 | 1 | 0.2
Glass | 0.2 | 0.2 | 0.2 | 0.2 | 0.2 | 0.2 | 0
Frosted Glass | 0.2 | 0.2 | 0.2 | 0.2 | 0.2 | 0.2 | 0.8
Water | 0.2 | 0.2 | 0.5 | 0.2 | 0.2 | 0.2 | 0
Mirror | 1 | 1 | 1 | 1 | 1 | 1 | 0
Air | 0 | 0 | 0 | 0 | 0 | 0 | 0
Smoke/mist | 0.5 | 0.5 | 0.5 | 0.2 | 0.2 | 0.2 | 0.5

##  V Conclusions

The purpose of this communication is to start a conversation on a volumetric-
point-cloud model description that is independent of lighting and of camera
viewpoints. Rendering, or capturing/scanning should be independent of the
model description, allowing for example to change lighting. One example is to
be able to represent a building at day or night, even placing illumination
therein. All these options should be rendering choices, not an intrinsic part
of the point cloud.

Our next steps for such point clouds are: (i) synthetic generation; (ii)
rendering; (iii) compression; and (iv) scanning.

This proposal is intended for non-luminous materials only. Internal lights, of
any luminescence, which may include diverse light bulbs, may be represented
under different sets of attributes of the same point cloud.

##  VI Acknowledgment

The authors would like to thank Prof. Leonardo R. A. Xavier, from Universidade
de Brasilia, for the help with eletromagnetism issues.

## References

  * [1] M. Gross and H. Pfister, Point-Based Graphics, Morgan Kaufmann Series on Computer Graphics, Elsevier, 2007.
  * [2] J. D. Foley , A. Van Dam, and S. K. Feiner, Computer Graphics: Principles and Practices, 3rd Edition, Addison-Wesley Professional, 2013.
  * [3] 3D Immersive Virtual Experiences,
https://www.360virtualtour.co/3d-scanning/ .

  * [4] D. T. Paris, and F. K. Hurd, Basic Eletromagnetic Theory, McGraw Hill, New York, NY, US, 1969.
  * [5] B. T. Phong, “Illumination for computer generated pictures,” Coomunications of tha ACM, Graphics and Image Processing, Vol. 18, No. 6, June 1975.
  * [6] J. F. Blinn, “Models of light reflection for computer synthesized pictures,” Proc. 4th Annual Conference on Computer Graphics and Interactive Techniques, pp. 192–-198, 1977.
  * [7] C. Schlick,  “An inexpensive BRDF model for physically-based rendering,”  Proc. Eurographics Computer Graphics Forum, Vol. 13, No. 3, 1994.


# Fisher Matrix Stability

Naren Bhandari<EMAIL_ADDRESS>McWilliams Center for Cosmology,
Department of Physics, Carnegie Mellon University,
Pittsburgh, PA 15213, USA  C. Danielle Leonard<EMAIL_ADDRESS>School of Mathematics, Statistics and Physics, Newcastle University, Newcastle
upon Tyne, NE1 7RU, UK  Markus Michael Rau<EMAIL_ADDRESS>McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon
University,
Pittsburgh, PA 15213, USA  Rachel Mandelbaum<EMAIL_ADDRESS>McWilliams Center for Cosmology, Department of Physics, Carnegie Mellon
University,
Pittsburgh, PA 15213, USA

##  1 Abstract

Fisher forecasts are a common tool in cosmology with applications ranging from
survey planning to the development of new cosmological probes. While
frequently adopted, they are subject to numerical instabilities that need to
be carefully investigated to ensure accurate and reproducible results. This
research note discusses these challenges using the example of a weak lensing
data vector and proposes procedures that can help in their solution.

##  2 Introduction

With the advent of ongoing large photometric surveys like the Dark Energy
Survey (DES; The Dark Energy Survey Collaboration, 2005), the Kilo-Degree
Survey (KiDS; de Jong et al., 2013), the Hyper Suprime-Cam survey (HSC; Aihara
et al., 2018) and upcoming surveys like the Vera C. Rubin Observatory Legacy
Survey of Space and Time (LSST; Ivezić et al., 2019) and Euclid (Laureijs et
al., 2011), the need for computationally efficient forecasting for survey
planning or to investigate the impact of sources of systematic error
increases. Fisher forecasts are a commonly used tool for this purpose, that
has been adopted early in many areas of astrophysics and cosmology (e.g.
Bassett et al., 2011; Cutler & Flanagan, 1994; Jungman et al., 1996a, b;
Tegmark, 1997; Tegmark et al., 1997; Hart & Chluba, 2020) to forecast
parameter constraints without the barrier of the computational expense of e.g.
Markov chain Monte Carlo methods. However, in making use of them, there are
issues of which we must be aware. We highlight the importance of rigorous
investigations of numerical stability in Fisher forecasts while using a weak
lensing convergence power spectrum data vector. We discuss possible pitfalls
in the application of Fisher forecasts and propose methodological tests to
ensure numerical accuracy. Our work builds on that of Yahia-Cherif et al.
(2020), where this question was considered for the complementary case of a
data vector consisting of the galaxy power spectrum measured from a
spectroscopic galaxy survey (with a specific focus on Euclid), and on that of
Blanchard et al. (2020), which performed an analytic estimation of the
relationship between the maximum fractional error on an element of the Fisher
matrix and the error on the elements of the resulting forecast parameter
covariance matrix.

In the course of this investigation we will use two libraries that provide the
capability to model weak lensing convergence power spectra: CCL v2.0.1
(Chisari et al., 2019) and CosmoSIS v1.6 (Zuntz et al., 2015). We note that
the settings for CosmoSIS that we use in this work are optimized for MCMC
sampling and not for high accuracy Fisher matrix calculations. This is not a
limitation of the code itself, since higher accuracy configurations are
supported. The purpose of this work is not to perform a code comparison, but
rather to investigate the importance and impact of numerical derivative
calculations on the reliability of Fisher forecasts, and to provide tests that
can be used to avoid numerical pitfalls. The CosmoSIS code package has been
subject to significant scrutiny in terms of the accuracy of likelihood
evaluations in e.g. Krause et al. (2017).

The structure of this research note is as follows. In the ‘Forecast
Methodology’ section, we outline our fiducial assumptions and how we set up
our Fisher forecast. In the ‘Results’ section, we detail our findings using
both libraries and gauge their numerical stability in the context of Fisher
forecasting. Finally, in the ‘Conclusions’, we summarize and present our
conclusions regarding these forecasts.

##  3 Forecast Methodology

We perform a Fisher forecast for an LSST-like convergence power spectrum data
vector assuming the fiducial cosmology and survey setup summarized in Table 1.
We use $\theta=(\Omega_{c},w_{0},h,10^{9}A_{s},\Omega_{b},n_{s},w_{a})$ as our
cosmological parameter vector and we use a redshift distribution111We use the
zdistri_model_z0=1.100000e-01_beta=6.800000e-01_Y10_source redshift
distribution from the LSST Dark Energy Science Collaboration (DESC) Science
Requirements Document v1 (The LSST Dark Energy Science Collaboration et al.,
2018a, b) which we bin into five tophat tomographic bins that are equally
spaced in redshift. Figure 1 shows the redshift probability density function,
$p(z)$, i.e. the probability that a randomly selected galaxy has redshift $z$,
along with the five tophat selections. Furthermore, we take the fractional sky
coverage, $f_{\rm{sky}}$, from Table 1.1 in LSST Science Collaboration et al.
(2009), and choose the number density, $\overline{n}_{g}$, to be consistent
with Table 1 in Chang et al. (2013).

Table 1: Fiducial Parameter Values and Survey Assumptions Cosmological/Survey Parameter |  | Value
---|---|---
Dark Matter Density | Ω_c | 0.25
Baryon Density | Ω_b | 0.04
Dark Energy Equation of State | w_0 | -0.90
$w(a)=w_{0}+w_{a}(1-a)$ | w_a | 0.00
Reduced Hubble Constant | h | 0.72
Primordial Power Spectrum Amplitude | 10^9 A_s | 2.10
Slope of Primordial Power Spectrum | n_s | 0.96
Shape Noise |  σ_ϵ^2 | 0.23
Fractional Sky Coverage | f_sky | 0.48
Number of Galaxies per Steradian | ¯n_g | 3.1×10^8
Figure 1: Redshift probability density function, $p(z)$, as a function of
redshift, $z$.

### Convergence Power Spectrum Modelling

Using CCL and CosmoSIS, we calculate the shear-shear angular power spectra
which serve as our data vector. The convergence power spectrum
$C^{\kappa}(\ell)$ probes the spatial correlation of the lensing convergence
signal (e.g. Blandford et al., 1991; Kaiser, 1992; Bartelmann & Schneider,
2001)

$C^{\kappa}(\ell)=\frac{9}{4}\Omega_{m}^{2}H_{0}^{4}\int_{0}^{\chi_{H}}\mathrm{d}\,\chi\,\frac{q^{2}(\chi)}{a^{2}(\chi)}P_{\rm\delta}\left(k=\frac{\ell}{f_{K}(\chi)},\chi\right)\,.$
(1)

Here, $\Omega_{m}$, $H_{0}$ denote the matter density and Hubble constant
respectively, $\chi_{H}$ is the comoving distance $\chi$ for which
$\chi(a)\rightarrow\chi_{H}$ as the scale factor $a\rightarrow 0$. Note that
$\Omega_{m}=\Omega_{b}+\Omega_{c}$. From the structure of Eq. 1 we see that
the matter power spectrum
$P_{\rm\delta}\left(k=\frac{\ell}{f_{K}(\chi)},\chi\right)$ is projected along
the line-of-sight using a geometrical factor, the ‘lensing weight’:

$q(\chi)=\int_{\chi}^{\chi_{H}}\mathrm{d}\,\chi^{{}^{\prime}}\,n(\chi^{{}^{\prime}})\left(\frac{f_{K}(\chi^{{}^{\prime}}-\chi)}{f_{K}(\chi^{{}^{\prime}})}\right)\,,$
(2)

where the function $f_{K}(\chi)$ is a curvature-dependent factor in the FLRW
metric. We see that this geometrical factor critically depends on the redshift
distribution $n(\chi^{{}^{\prime}})$ that constitutes an important source of
systematic uncertainty in weak lensing measurements. For a more detailed
treatment of the convergence power spectrum we refer the reader to the
aforementioned references. For the matter power spectrum $P_{\delta}$, we use
the Takahashi et al. (2012) version of the nonlinear power spectrum from
Halofit (Smith et al., 2003).

### Covariance Matrix Modelling

As in Kirk et al. (2015), we use the convergence power spectrum to calculate
the covariance matrix. The Gaussian model for the covariance matrix of the
data vector is given as

$\Sigma^{(k,m)}_{(i,j)}(\ell)=\frac{1}{(2\ell+1)f_{\rm
sky}}\left(\overline{C}^{(i,k)}(\ell)\overline{C}^{(j,m)}(\ell)+\overline{C}^{(i,m)}(\ell)\overline{C}^{(j,k)}(\ell)\right)\,,$
(3)

where

$\overline{C}^{(i,j)}(\ell)=C^{(i,j)}(\ell)+\delta_{i,j}\,\left(\frac{\sigma_{\epsilon}^{2}}{2\,\overline{n}_{g}^{i}}\right)\,.$
(4)

Here, $\overline{n}_{g}^{i}$ denotes the number of galaxies per steradian and
$\sigma_{\epsilon}$ is the shape noise expressed as a standard deviation.

Note that the indices denote samples selected in bins of increasing redshift,
i.e. tomographic bins, whose redshift probability density functions are shown
in Fig. 1. The index combinations $(i,j)$ denote the respective auto- ($i=j$)
and cross-correlations ($i\neq j$) obtained by correlating the samples in the
tomographic selections $i$ and $j$. For example, at fixed $\ell$, the
correlation between the auto-correlation convergence power spectrum of the
first bin with the cross-correlation convergence power spectrum of first and
second tomographic redshift bin, would be $\Sigma^{(k=1,m=1)}_{(i=1,j=2)}$. At
each $\ell$, the data-vector, denoted by $C(\ell)$, is constructed from all
convergence power spectra formed from the available tomographic redshift bins,
i.e., by enumerating the possible combinations we can map the pair of indices
$(i,j)$ to a single index and make a vector. Similarly, at each $\ell$, the
data-covariance matrix, denoted by $\Sigma(\ell)$, is constructed such that
the rows and columns are indexed by the enumeration above.

### Fisher Forecast

Suppose that the negative log-likelihood of our data vector is
$\mathcal{L}(C(\ell);\theta)$. Then the Fisher matrix is defined as the
expectation value of the Hessian matrix of $\mathcal{L}$ with respect to the
parameters $\theta$:

$F_{\alpha\beta}=\left\langle\frac{\partial^{2}\mathcal{L}}{\partial\theta_{\alpha}\partial\theta_{\beta}}\right\rangle$
(5)

Assuming a Gaussian likelihood and that $\Sigma(\ell)$ is independent of
$\theta$, we can calculate the Fisher matrix as

$F_{\alpha\beta}=\sum_{\ell}\frac{\partial
C(\ell)^{T}}{\partial\theta_{\alpha}}\cdot\Sigma(\ell)^{-1}\cdot\frac{\partial
C(\ell)}{\partial\theta_{\beta}}$ (6)

where the derivatives are evaluated at the fiducial parameter values (Kirk et
al., 2015). The sum is taken over $\ell_{\min{}}=76$ to $\ell_{\max{}}=999$ —
this range is motivated by accessible values to the aforementioned surveys and
also to stay in the linear and quasi-linear regime. We refer the reader to
Tegmark (1997); Tegmark et al. (1997); Bassett et al. (2011); Heavens et al.
(2014) for more detailed treatments of the Fisher matrix formalism.

### Numerical Differentiation

We calculate the derivatives above numerically by using centered finite
difference methods outlined in Fornberg (1988):

$\frac{\partial f}{\partial
x}\bigg{|}_{x=x_{0}}=\sum_{j=-n}^{n}\mu_{j}f(x+j\Delta s)$ (7)

where $f$ is some function that depends on the variable $x$, $\mu_{j}$’s are
coefficients given by the aforementioned reference, $\Delta s$ is the step-
size, and $n$ is an integer that determines the number of points surrounding
the reference point where the derivative needs to be approximated. The
geometrical arrangement of these points is referred to as the ‘stencil’. The
number of points in the stencil, namely $2n+1$ for centered finite difference
methods, determines how the numerical error of the stencil method scales with
the step-size. For example, the five-point stencil approximation is given by

$\frac{\partial f}{\partial x}\bigg{|}_{x=x_{0}}=\frac{1}{12\Delta
s}f(x_{0}-2\Delta s)-\frac{8}{12\Delta s}f(x_{0}-\Delta
s)+0f(x_{0})+\frac{8}{12\Delta s}f(x_{0}+\Delta s)-\frac{1}{12\Delta
s}f(x_{0}+2\Delta s).$ (8)

Note that all the centered finite differences we use have a zero coefficient
on the central term. As such, while entering into the theoretical modelling of
the covariance matrix, the central term does not enter the calculation of the
numerical derivative. Thus, the order of accuracy is given by $2n$, i.e.,
evaluating more points in our stencil will lead to more accurate derivatives
at the cost of computational complexity (here, more calls to CCL or CosmoSIS).
We will investigate the numerical performance as a function of the number of
points in the stencil and the step-size. We do so by calculating derivatives
using three-, five-, seven-, and nine-point stencils and seeing how some
figure of merit changes in response to varying the step-size $\Delta s$,
looking for a region where it is stable — henceforth, we will refer to this
process as tuning the derivatives. In this context it is noteworthy that other
derivative tuning methods have been employed by the community. For example
Camera et al. (2017), due the large number of observables and cosmological
parameters in their work, use a linear interpolation constructed on a set of
sampling points around the fiducial value for derivative estimation. They then
decrease the point spacing adaptively until the required numerical accuracy is
met. In this work we will restrict our approach to a visual approach to gain
understanding of the behaviour of different codes as a function of step size.
For the sake of comparisons, it will be convenient to work with a normalized
step-size given by the nominal step-size divided by the absolute value of the
fiducial parameter value instead. For the parameter, $w_{a}$, which has a
fiducial value of 0, we divide the nominal step size by the absolute value of
the fiducial value of $w_{0}$, which is nonzero.

### Performance Evaluation Metrics

In order to gauge how well the measurement we are forecasting constrains dark
energy, we use the Dark Energy Task Force’s Figure of Merit (DETF FoM) given
by

$\text{DETF FoM}=\sqrt{\det{F^{w_{0}w_{a}}}}$ (9)

where $F^{w_{0}w_{a}}$ is the reduced Fisher matrix, obtained by inverting the
Fisher matrix $F$, marginalizing over parameters that are not $w_{0}$ or
$w_{a}$, and inverting back (Albrecht et al., 2006). Note that we can also
define similar quantities for other pairs of parameters, like $\Omega_{c}$ and
$10^{9}A_{s}$ for example. To distinguish which pairs of parameters we are
talking about, we refer explicitly to the $w_{0}-w_{a}$ FoM (as in Eq. 9) or
$\Omega_{c}-10^{9}A_{s}$ FoM where necessary. Since the main purpose of this
note is to study how numerical errors can bias Fisher matrix inference, we are
not marginalizing over intrinsic alignments or other sources of systematic
uncertainty.

##  4 Results

Both CCL and CosmoSIS are well-trusted and mature libraries that we would
expect to produce theoretical predictions and forecasts that are consistent,
so we begin by checking this assumption. As a first test of consistency
between the data vectors produced by CCL and CosmoSIS we investigate the
signal-to-noise ratio (${\rm SNR}$) given by:

${\rm SNR}=\sqrt{\sum_{\ell}C(\ell)\cdot\Sigma(\ell)^{-1}\cdot C(\ell)}$ (10)

These were calculated to be $342.58$ and $342.44$ for CCL and CosmoSIS,
respectively — a relative difference of $-0.04\%$, which we consider to be a
good agreement. This first test is important to ensure that the data vector
and the covariance calculation give reasonable results, before numerical
differentiation introduces numerical sources of error that need to be
isolated.

In order to calculate the Fisher Matrices, we tune the derivatives as
mentioned in the ‘Forecast Methodology’ Section. We select 0.01 as an initial
guess for a starting normalized step-size for all parameters and calculate all
the derivatives using a five-point stencil and arrive at an initial DETF FoM
estimate. To then tune the derivatives, we calculate them using the three-,
five-, seven-, and nine-point stencils, varying the normalized step-size to
look for regions where the DETF FoM does not vary much. In the CCL case, we
were able to find clear regions of stability, wherein the normalized step-
sizes for all parameters was chosen to be 0.01 using a nine-point stencil — a
few examples can be seen in Figure 2. This spacing is of the same order of
magnitude as that determined to be optimal for each cosmological parameter in
Yahia-Cherif et al. (2020) for the spectroscopic galaxy case; this makes sense
as both observables rely on underlying calculations of the power spectrum via
Boltzmann codes. In the CosmoSIS case, we were unable find such a region of
stability for the parameters $\Omega_{c}$ and $w_{0}$, seen in Figure 3.
However, for the sake of further comparison, we can pick values of the
normalized step-sizes and the number of points in the stencil used (like 0.02
for the nine-point stencil for $\Omega_{c}$) that bring the DETF FoM as close
as possible to its CCL counterpart, as seen in the third column of Table 2.
(Despite this, we still see, in the next column, that the analogous FoM for
$\Omega_{c}$ and $10^{9}A_{s}$ are quite different between CCL and CosmoSIS).
While we were unable to successfully tune derivatives for CosmoSIS, there are
other methods of forecasting with the library. We note that our CosmoSIS setup
is mainly optimized for MCMC techniques that may not require a high numerical
precision in the calculation of angular power spectra. Both the internal
Boltzmann solvers and the CosmoSIS modules that perform the integration over
the window functions can be adjusted towards increased numerical accuracy and
are subject to future accuracy improvements. We reiterate that we aren’t
trying to compare the numerical accuracy of the two libraries, but use them to
investigate numerical issues that can arise in Fisher forecasts.

(a)

(b)

Figure 2: Derivative tunings for CCL, focusing on parameters $\Omega_{c}$
(left) and $w_{0}$ (right), using various order stencils. Note that a
normalized step-size of 0.01 appears stable.

(a)

(b)

Figure 3: Derivative tunings for CosmoSIS, focusing on parameters $\Omega_{c}$
(left) and $w_{0}$ (right) using various order stencils. Note that there does
not seem to be a clear region of stability in either case. In the latter case,
this is true even restricting the range of the plot to focus on the larger
normalized step-sizes.

Continuing with the calculation of our Fisher forecast, we get Fisher matrices
from CCL and CosmoSIS that have entry-wise relative differences that are no
larger than $0.6\%$ in magnitude as seen in Table 3. However, these Fisher
matrices have high condition numbers (see Table 2) that lead to large
magnifications in the entry-wise relative differences when we invert the
matrices and look at the corresponding covariance matrices as seen in Table 4.
At worst, the relative difference is almost a factor of 15. This highlights
how small instabilities in the Fisher matrix can translate to much larger
deviations in the parameter covariance matrix. Note that we use the inversion
method, `linalg.inv`, implemented in the NumPy package (Oliphant, 2006; Van
Der Walt et al., 2011) that is based on LU decomposition. In the following
paragraph we will describe the consistency tests we performed to ensure
numerical accuracy in the inversion.

Given the high condition numbers, we gauge how accurate our calculated
inverses are by one of two methods:

  1. 1.

We multiply a Fisher matrix with its corresponding covariance matrix (its
inverse) and check whether we get something that is numerically close to the
identity matrix.

  2. 2.

We follow Köhlinger et al. (2017), which uses a technique from Newman et al.
(1974) to verify matrix inversions in its forecasts. The latter gives upper
and lower bounds on the matrix norm of the difference between the calculated
and actual inverse of a matrix. Let $C$ be our calculated inverse to the
Fisher matrix $F$, let $F^{-1}$ be the actual inverse of $F$, let $I$ be the
Identity matrix, let $R=I-FC$ and let $N$ be a matrix norm (we use the
Frobenius norm). Then

$\frac{N(CR)}{1+N(R)}\leq N(F^{-1}-C)\leq\frac{N(CR)}{1-N(R)}$ (11)

For CCL, the lower bound was calculated to be $1.88\times 10^{-12}$ and the
upper bound differed by $1.70\times 10^{-22}$. For CosmoSIS, the respective
lower bound was $3.09\times 10^{-13}$ and the upper bound differed by
$2.28\times 10^{-23}$. By subtracting the lower bound from the upper bound, we
get the size of the interval that contains the norm in question. This ‘Inverse
Norm Spread’ is shown again in the last column of Table 2. Since these spreads
are very small in all cases, we have another test suggesting that our
calculated inverses are accurate.

We also test the sensitivity of our results with respect to a
reparametrization wherein we use $\Omega_{c}h^{2}$ and $\Omega_{b}h^{2}$
instead of $\Omega_{c}$ and $\Omega_{b}$ respectively, and repeat the above
steps. The last row of Table 2 summarize the condition number, the DETF and
$\Omega_{c}h^{2}-A_{s}$ FoM, and the Inverse Norm Spread for this case using
CCL. While the condition number is still high, the DETF FoM is roughly the
same as earlier and the small Inverse Norm Spread suggests our calculated
covariance matrix is numerically accurate.

In the field of numerical stability, it is well known that when inverting a
matrix, it should be equilibrated, i.e., given a norm, its rows (and columns)
should have similar lengths in that norm. This ensures that small
perturbations in all rows (or columns) are of similar magnitude. For a more
detailed treatment of the subject, we refer the reader to Croz & Higham
(1992); Wilkinson (1961); Bunch (1971). In the context of Fisher forecasting,
if one is using $A_{s}$ as a parameter, equilibration can be used to lower the
condition number. For example, we can scale $A_{s}$ by $10^{9}$ and use
$10^{9}A_{s}$ as a parameter. While this does not affect the relative
differences in the Fisher and covariance matrix, we found that it does make
eigenvalue/eigenvector calculations more accurate. Since doing so comes at no
computational cost, we incorporated this into our fiducial ansatz. Note that
$\log[10^{10}A_{s}]$ has been used as well (Planck Collaboration et al., 2016;
Raveri et al., 2016; Ballardini et al., 2018).

Additionally, we also test how adding priors on $\Omega_{b}$ and $h$ affects
our results. Using Monte Carlo Markov Chains provided in the 2015 release of
data products222We use the chains in the base_plikHM_TT_lowTEB_post_BAO.txt
files from the ‘Planck baseline model: TT low-ell temperature and LFI
polarization (bflike, 2 <= ell <= 29)’ data available at
https://pla.esac.esa.int/pla/aio/product-
action?COSMOLOGY.FILE_ID=COM_CosmoParams_base-plikHM-TT-lowTEB_R2.00.tar.gz by
the Planck Collaboration et al. (2016), we create a covariance matrix
containing the parameters we are interested in. By taking the square root of
the diagonal elements corresponding to $\Omega_{b}$ and $h$, we get the
uncertainties $\sigma_{\Omega_{b}}=0.00066$ and $\sigma_{h}=0.0057$. To add
these priors to our Fisher matrices, we add $\sigma_{h}^{-2}$ and
$\sigma_{\Omega_{b}}^{-2}$ to the corresponding diagonal values in the
original matrices. We therefore neglect, for simplicity, the correlation
between $h$ and $\Omega_{b}$ in the Planck prior. We note that our goal is not
primarily to perform a realistic forecast, but to study the impact of prior
choices on the numerical stability of the Fisher analysis in a simple
framework.

While adding this prior reduces the largest entry-wise relative differences
between the covariance matrices (for the $h-w_{0}$ and $n_{s}-w_{0}$ elements)
from a factor of $\sim 15$ to $\sim 0.5$, it also increases other elements
like the $h-10^{9}A_{s}$ element from a factor of $\sim 0.5$ to $\sim 2$. As
such, we do not find that adding priors solves the problem of high entry-wise
relative differences upon inversion of Fisher matrices.

Lastly, it is worth noting that the high condition numbers appear to stem from
the size of our parameter space. By only varying a subset of our original
parameters, we can reduce the size of the parameter space and generate new
Fisher matrices using CCL and CosmoSIS. After trying various combinations of
parameters, we found that once the parameter space included at least four
parameters, then entry-wise relative differences between the inverted Fisher
matrices started to exceed 10%, despite much smaller entry-wise relative
differences, of similar magnitude to those in Table 4, in the original Fisher
matrices. This agrees with idea that increasing the dimensionality of a matrix
to be inverted results in this inversion being more sensitive to small
deviations in the original matrix, as discussed in a Fisher forecasting in
e.g. Blanchard et al. (2020).

Table 2: Summary table for CCL and CosmoSIS forecasts.  Forecast | Condition Number |  $w_{0}-w_{a}$ FoM |  $\Omega_{c}-10^{9}A_{s}$ FoM |  $\Omega_{c}h^{2}-10^{9}A_{s}$ FoM | Inverse Norm Spread
---|---|---|---|---|---
CCL | $1.07\times 10^{7}$ | $25.24$ | $96.07$ | N/A | $1.70\times 10^{-22}$
CosmoSIS | $8.08\times 10^{6}$ | $29.50$ | $221.5$ | N/A | $2.28\times 10^{-23}$
CCL ($\Omega_{c}h^{2}$ & $\Omega_{b}h^{2}$) | $3.33\times 10^{7}$ | $26.84$ | N/A | $70.89$ | $1.55\times 10^{-21}$
Table 3: Magnitude of relative differences for CCL and CosmoSIS Fisher matrices. The largest differences are highlighted in bold.  | $\Omega_{c}$ | $w_{0}$ | $h$ | $10^{9}A_{s}$ | $\Omega_{b}$ | $n_{s}$ | $w_{a}$
---|---|---|---|---|---|---|---
$\Omega_{c}$ | 4.68E-5 | 1.92E-3 | 3.26E-4 | 5.98E-4 | 1.49E-3 | 1.01E-4 | 1.93E-4
$w_{0}$ | 1.92E-3 | 3.84E-3 | 2.39E-3 | 1.53E-3 | 5.61E-4 | 1.36E-3 | 2.34E-3
$h$ | 3.26E-4 | 2.39E-3 | 7.20E-4 | 1.77E-4 | 1.41E-3 | 4.82E-5 | 5.84E-4
$10^{9}A_{s}$ | 5.98E-4 | 1.53E-3 | 1.77E-4 | 1.05E-3 | 2.47E-3 | 9.57E-4 | 4.11E-4
$\Omega_{b}$ | 1.49E-3 | 5.61E-4 | 1.41E-3 | 2.47E-3 | 5.90E-3 | 1.84E-4 | 2.51E-3
$n_{s}$ | 1.01E-4 | 1.36E-3 | 4.82E-5 | 9.57E-4 | 1.84E-4 | 6.15E-4 | 6.83E-4
$w_{a}$ | 1.93E-4 | 2.34E-3 | 5.84E-4 | 4.11E-4 | 2.51E-3 | 6.83E-4 | 3.32E-5
Table 4: Magnitude of relative differences for CCL and CosmoSIS covariance matrices. The largest differences are highlighted in bold.  | $\Omega_{c}$ | $w_{0}$ | $h$ | $10^{9}A_{s}$ | $\Omega_{b}$ | $n_{s}$ | $w_{a}$
---|---|---|---|---|---|---|---
$\Omega_{c}$ | 7.91E-1 | 1.14E+0 | 7.71E-1 | 7.36E-1 | 7.82E-1 | 6.20E-1 | 1.00E+0
$w_{0}$ | 1.14E+0 | 3.54E-1 | 1.15E+1 | 2.88E-2 | 7.35E-1 | 1.48E+1 | 3.60E-1
$h$ | 7.71E-1 | 1.15E+1 | 6.73E-1 | 4.71E-1 | 6.82E-1 | 4.20E-1 | 1.31E+1
$10^{9}A_{s}$ | 7.36E-1 | 2.88E-2 | 4.71E-1 | 2.00E-1 | 4.66E-1 | 3.17E-1 | 3.34E-2
$\Omega_{b}$ | 7.82E-1 | 7.35E-1 | 6.82E-1 | 4.66E-1 | 6.51E-1 | 5.95E-1 | 5.55E-1
$n_{s}$ | 6.20E-1 | 1.48E+1 | 4.20E-1 | 3.17E-1 | 5.95E-1 | 3.15E-4 | 3.45E-1
$w_{a}$ | 1.00E+0 | 3.60E-1 | 1.31E+1 | 3.34E-2 | 5.55E-1 | 3.45E-1 | 2.99E-1

##  5 Conclusions

In this note, we investigated numerical issues that may arise when Fisher
forecasting LSST-like weak lensing probes to constrain dark energy, focusing
on an angular power spectrum formalism whose performance was gauged by the
DETF FoM. Using two numerical libraries, CCL and CosmoSIS, we generated data-
vectors and data-covariances that were very close numerically, as seen by
their signal-to-noise ratios having a relative difference of $-0.04$%.
However, when trying to calculate the Fisher matrix, we found that numerical
instability can arise when tuning the step-size for the derivatives that go
into Equation 6, which consequently introduces numerical instability in the
DETF FoM as seen in Figure 3. Furthermore, even with two Fisher matrices whose
entry-wise relative differences were small (around 0.6% at most), we found
that these small perturbations can be magnified greatly upon inversion (to
order of magnitude differences) due to these Fisher matrices having high
condition numbers — due, at least in part, to the number of parameters chosen
to investigate. Our work agrees with similar findings in Blanchard et al.
(2020) and Yahia-Cherif et al. (2020), with the former estimating that for the
forecast parameter covariance matrix to be stable to 10%, the maximum error on
an element of the Fisher matrix permitted would be 0.01%. Ensuring all
parameters have a similar order of magnitude may help reduce the condition
number; however, reparameterizing $\Omega_{c}$ by $\Omega_{c}h^{2}$ and
$\Omega_{b}$ and $\Omega_{b}h^{2}$, or even applying priors on certain
parameters does not necessarily yield significant benefits. While the large
condition numbers may not necessarily propagate into large differences for the
DETF FoM, as in our case, it worth carrying out similar numerical stability
tests in future work with Fisher forecasts, as it may impact other quantities
whose derivation involves inverting the Fisher matrix.

We thank the CCL team for making their code available at
https://github.com/LSSTDESC/CCL and the CosmoSIS team for making their code
available at https://bitbucket.org/joezuntz/cosmosis. We also thank Alain
Blanchard, Stefano Camera, Luke Hart, Mark Kamionkowski, Eric Linder, Anz̆e
Slosar and Joe Zuntz, for useful discussions that improved the content of this
paper.

## References

  * Aihara et al. (2018) Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, PASJ, 70, S4, doi: 10.1093/pasj/psx066
  * Albrecht et al. (2006) Albrecht, A., Bernstein, G., Cahn, R., et al. 2006, Report of the Dark Energy Task Force.  https://arxiv.org/abs/astro-ph/0609591
  * Ballardini et al. (2018) Ballardini, M., Finelli, F., Maartens, R., & Moscardini, L. 2018, Journal of Cosmology and Astroparticle Physics, 2018, 044, doi: 10.1088/1475-7516/2018/04/044
  * Bartelmann & Schneider (2001) Bartelmann, M., & Schneider, P. 2001, Phys. Rep., 340, 291, doi: 10.1016/S0370-1573(00)00082-X
  * Bassett et al. (2011) Bassett, B. A., Fantaye, Y., Hlozek, R., & Kotze, J. 2011, International Journal of Modern Physics D, 20, 2559–2598, doi: 10.1142/s0218271811020548
  * Blanchard et al. (2020) Blanchard, A., Camera, S., Carbone, C., et al. 2020, Astronomy & Astrophysics, 642, A191
  * Blandford et al. (1991) Blandford, R. D., Saust, A. B., Brainerd, T. G., & Villumsen, J. V. 1991, MNRAS, 251, 600, doi: 10.1093/mnras/251.4.600
  * Bunch (1971) Bunch, J. R. 1971, J. ACM, 18, 566–572, doi: 10.1145/321662.321670
  * Camera et al. (2017) Camera, S., Harrison, I., Bonaldi, A., & Brown, M. L. 2017, MNRAS, 464, 4747, doi: 10.1093/mnras/stw2688
  * Chang et al. (2013) Chang, C., Jarvis, M., Jain, B., et al. 2013, MNRAS, 434, 2121, doi: 10.1093/mnras/stt1156
  * Chisari et al. (2019) Chisari, N. E., Alonso, D., Krause, E., et al. 2019, The Astrophysical Journal Supplement Series, 242, 2, doi: 10.3847/1538-4365/ab1658
  * Croz & Higham (1992) Croz, J. J. D., & Higham, N. J. 1992, IMA Journal of Numerical Analysis, 12, 1, doi: 10.1093/imanum/12.1.1
  * Cutler & Flanagan (1994) Cutler, C., & Flanagan, É. E. 1994, Phys. Rev. D, 49, 2658, doi: 10.1103/PhysRevD.49.2658
  * de Jong et al. (2013) de Jong, J. T. A., Verdoes Kleijn, G. A., Kuijken, K. H., & Valentijn, E. A. 2013, Experimental Astronomy, 35, 25, doi: 10.1007/s10686-012-9306-1
  * Fornberg (1988) Fornberg, B. 1988, Mathematics of computation, 51, 699
  * Hart & Chluba (2020) Hart, L., & Chluba, J. 2020, MNRAS, 495, 4210, doi: 10.1093/mnras/staa1426
  * Heavens et al. (2014) Heavens, A. F., Seikel, M., Nord, B. D., et al. 2014, Monthly Notices of the Royal Astronomical Society, 445, 1687–1693, doi: 10.1093/mnras/stu1866
  * Ivezić et al. (2019) Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111, doi: 10.3847/1538-4357/ab042c
  * Jungman et al. (1996a) Jungman, G., Kamionkowski, M., Kosowsky, A., & Spergel, D. N. 1996a, Phys. Rev. D, 54, 1332, doi: 10.1103/PhysRevD.54.1332
  * Jungman et al. (1996b) —. 1996b, Phys. Rev. Lett., 76, 1007, doi: 10.1103/PhysRevLett.76.1007
  * Kaiser (1992) Kaiser, N. 1992, ApJ, 388, 272, doi: 10.1086/171151
  * Kirk et al. (2015) Kirk, D., Lahav, O., Bridle, S., et al. 2015, MNRAS, 451, 4424, doi: 10.1093/mnras/stv1268
  * Krause et al. (2017) Krause, E., Eifler, T. F., Zuntz, J., et al. 2017, arXiv e-prints, arXiv:1706.09359.  https://arxiv.org/abs/1706.09359
  * Köhlinger et al. (2017) Köhlinger, F., Viola, M., Joachimi, B., et al. 2017, Monthly Notices of the Royal Astronomical Society, 471, 4412, doi: 10.1093/mnras/stx1820
  * Laureijs et al. (2011) Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv preprint arXiv:1110.3193
  * LSST Science Collaboration et al. (2009) LSST Science Collaboration, Abell, P. A., Allison, J., et al. 2009, arXiv e-prints, arXiv:0912.0201.  https://arxiv.org/abs/0912.0201
  * Newman et al. (1974) Newman, M., et al. 1974, J. Res. Nat. Bur. Stand., 78, 65
  * Oliphant (2006) Oliphant, T. E. 2006, A guide to NumPy, Vol. 1 (Trelgol Publishing USA)
  * Planck Collaboration et al. (2016) Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016, A&A, 594, A1, doi: 10.1051/0004-6361/201527101
  * Raveri et al. (2016) Raveri, M., Martinelli, M., Zhao, G., & Wang, Y. 2016, Information Gain in Cosmology: From the Discovery of Expansion to Future Surveys.  https://arxiv.org/abs/1606.06273
  * Smith et al. (2003) Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, Monthly Notices of the Royal Astronomical Society, 341, 1311–1332, doi: 10.1046/j.1365-8711.2003.06503.x
  * Takahashi et al. (2012) Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, The Astrophysical Journal, 761, 152, doi: 10.1088/0004-637x/761/2/152
  * Tegmark (1997) Tegmark, M. 1997, Physical Review Letters, 79, 3806–3809, doi: 10.1103/physrevlett.79.3806
  * Tegmark et al. (1997) Tegmark, M., Taylor, A. N., & Heavens, A. F. 1997, The Astrophysical Journal, 480, 22–35, doi: 10.1086/303939
  * The Dark Energy Survey Collaboration (2005) The Dark Energy Survey Collaboration. 2005, The Dark Energy Survey.  https://arxiv.org/abs/astro-ph/0510346
  * The LSST Dark Energy Science Collaboration et al. (2018a) The LSST Dark Energy Science Collaboration, Mandelbaum, R., Eifler, T., et al. 2018a, The LSST Dark Energy Science Collaboration (DESC) Science Requirements Document v1 Released Data Products, 1.0.1, Zenodo, doi: 10.5281/zenodo.2662127
  * The LSST Dark Energy Science Collaboration et al. (2018b) —. 2018b, arXiv e-prints, arXiv:1809.01669.  https://arxiv.org/abs/1809.01669
  * Van Der Walt et al. (2011) Van Der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science & Engineering, 13, 22
  * Wilkinson (1961) Wilkinson, J. H. 1961, J. ACM, 8, 281–330, doi: 10.1145/321075.321076
  * Yahia-Cherif et al. (2020) Yahia-Cherif, S., Blanchard, A., Camera, S., et al. 2020, arXiv preprint arXiv:2007.01812
  * Zuntz et al. (2015) Zuntz, J., Paterno, M., Jennings, E., et al. 2015, Astronomy and Computing, 12, 45–59, doi: 10.1016/j.ascom.2015.05.005