# Reconstruction of Refractive Indices from Spectral Measurements of
Monodisperse Aerosols

Tobias Kyrion  RWTH Aachen University, Aachen, Germany

###### Abstract

For the investigation of two-component aerosols one needs to know the
refractive indices of the two aerosol components. One problem is that they
depend on temparature and pressure, so one needs for their determination a
robust measurement instrument such as the FASP device, which can cope with
rigid environmental conditions. In this article we show that the FASP device
is capable of measuring the needed refractive indices, if monodisperse
aerosols of the pure components are provided. We determine the particle radii
of the monodisperse aerosols needed for this task and investigate how accurate
the measurements have to be in order to retrieve refractive indices in a
sufficient quality, such that they are suitable for investigations of two-
component aerosols.

key words:  nonlinear inverse problem; refractive index reconstruction; Mie
theory; nonlinear Tikhonov regularization; nonlinear regression

AMS Subject Classifications:  15A29; 34K29; 34M50; 45Q05; 70F17

##  1 Setup of the Experiment

This paper provides an algorithm for the reconstruction of refractive indices
from spectral measurements of monodisperse aerosols. The experiments we are
conducting are similar to the experiments presented in [1] with the difference
that we are using air as surrounding medium for the aerosol particles and that
temperature and pressure may approach $200^{\circ}\text{C}$ and $8$ bar
respectively. For these rigid conditions reliable databases for refractive
indices do not exist up to now. These refractive index databases are needed
for the measurement of particle size distributions of polydisperse aerosols
using the FASP.

As outlined in [2] the FASP measures light intensities $I_{long}(l)$ and
$I_{short}(l)$ having passed a long and a short measurement path length
$L_{long}$ and $L_{short}$ respectively. The evaluations of the FASP
measurements are based on the relation

$\int_{0}^{\infty}k(r,l)n(r)dr=e(l)\quad\text{with}\quad
e(l)=-\frac{\log(I_{long}(l))-\log(I_{short}(l))}{L_{long}-L_{short}},$ (1.1)

where $k(r,l):=\pi r^{2}Q_{ext}(m_{med}(l),m_{part}(l),r,l)$ is the so-called
kernel function, $l$ is the wavelength of the incident light, $r$ is the
radius of the spherical scattering particle and $m_{med}(l)$ and $m_{part}(l)$
are the refractive indices of the surrounding medium and the particle material
depending on the wavelength $l$. The function
$Q_{ext}(m_{med}(l),m_{part}(l),r,l)$ is the Mie extinction efficiency from
[3]. The function $n(r)$ is the size distribution of the scattering particels.
The right-hand side $e(l)$ in (1.1) is denoted as the spectral extinction.

Now if $n(r)$ is the size distribution of a monodisperse aerosol, where all
particles possess the same radius $r_{m}$, it is given by
$n(r)=n\delta(r-r_{m})$, where $n$ is the total number of particles and
$\delta(r-r_{m})$ is a Dirac delta distribution truncated on the positive
half-axis. Inserting this into (1.1) gives

$n\pi r_{m}^{2}Q_{ext}(m_{med}(l),m_{part}(l),r_{m},l)=e(l),$ (1.2)

hence the Mie extinction efficiency is measured directly at the radius
$r_{m}$.

The Mie extinction efficiency is given as an infinite series, i.e.

$Q_{ext}(m_{med}(l),m_{part}(l),r,l)=\sum_{n=1}^{\infty}q_{n}(m_{med}(l),m_{part}(l),r,l).$

The computation of the coefficient functions
$q_{n}(m_{med}(l),m_{part}(l),r,l)$ will be discussed in Section 2.

It is clear that in practical computations
$Q_{ext}(m_{med}(l),m_{part}(l),r,l)$ can only be approximated by a truncated
series, because only the computation of a finite number of the
$q_{n}(m_{med}(l),m_{part}(l),r,l)$’s is practically feasible.

We now fix a wavelength $l$. The complex refractive index $m_{part}(l)$ for
the wavelength $l$ is reconstructed from FASP measurements of several
monodisperse aerosols with particle radii $r_{1}$, …, $r_{N}$. Let
$q(r_{1},l)$, …, $q(r_{N},l)$ denote the measured spectral extinctions $e(l)$
corresponding to the particle radii $r_{1}$, …, $r_{N}$. We assume that they
are contaminated by additive Gaussian noise, i.e.
$q(r_{i},l)=q_{true}(r_{i},l)+\delta_{i}$ with
$\delta_{i}\sim\mathcal{N}(0,s_{i}^{2})$ for $i=1,...,N$. Furthermore we
assume that the standard deviations $s_{i}$ can be estimated from measurements
sufficiently accurately, such that we can regard them as known. We have

$q_{true}(r_{i},l)=n_{i}\pi
r_{i}^{2}\sum_{n=1}^{\infty}q_{n}(m_{med}(l),m_{part}(l),r_{i},l),$

where $n_{i}$ is the number of particles having the same radius $r_{i}$. Then
a reconstrution of $m_{part}(l)$ is obtained from the set of solutions $M(l)$
of the nonlinear regression problem

$M(l):=\underset{m\in\mathbb{C}}{\mathrm{argmin}}\;\sum_{i=1}^{N}\frac{1}{2\left(\frac{s_{i}}{n_{i}}\right)^{2}}\left(\pi
r_{i}^{2}\sum_{n=1}^{N_{tr}}q_{n}(m_{med}(l),m,r_{i},l)-\frac{q(r_{i},l)}{n_{i}}\right)^{2}.$
(1.3)

Note that $M(l)$ contains in general more than one solution, especially when
$q(r_{i},l)$ is perturbed by measurement noise. We discuss nonlinear
regression problems with truncated series expansions such as (1.3) in Section
4.

For solving (1.3) we use a global optimization strategy presented in Section 5
to generate reasonable candidates for start values for a local solver for a
regularized version of (1.3). Section 8 provides a selection method to find a
unique start value out of the candidates. In order to apply a gradient-based
local solver we must know the derivatives of the Mie extinction efficiency
series, which are discussed in Section 3.

##  2 Mie Theory

We recapitulate Mie theory in an absorbing medium as presented in [3]. Our
first step is to introduce the complex-valued Riccati-Bessel-functions
$\xi_{n}:\mathbb{C}\rightarrow\mathbb{C}$ and
$\psi_{n}:\mathbb{C}\rightarrow\mathbb{C}$ given by

$\xi_{n}(z)=\sqrt{\textstyle{\frac{\pi}{2}}}\sqrt{z}J_{n+\frac{1}{2}}(z)$
(2.1)

and

$\psi_{n}(z)=\sqrt{\textstyle{\frac{\pi}{2}}}\sqrt{z}J_{n+\frac{1}{2}}(z)+i\sqrt{\textstyle{\frac{\pi}{2}}}\sqrt{z}Y_{n+\frac{1}{2}}(z),$
(2.2)

with the Bessel functions $J_{n+\frac{1}{2}}:\mathbb{C}\rightarrow\mathbb{C}$
and $Y_{n+\frac{1}{2}}:\mathbb{C}\rightarrow\mathbb{C}$ of order
$n+\frac{1}{2}$ of first and second kind. We define the size parameter
$\rho=2\pi\frac{r}{l}$. Then we set $z_{med}:=\rho\cdot m_{med}$ and
$z_{part}:=\rho\cdot m_{part}$. Here and in the following we omit the
wavelength dependence of $m_{med}$ and $m_{part}$ for better readability. We
introduce the notation $m_{med}=n_{med}+ik_{med}$ and
$m_{part}=n_{part}+ik_{part}$.

We introduce the so-called Mie coefficients:

$\begin{split}a_{n}&:=\frac{m_{part}\dot{\xi_{n}}(z_{med})\xi_{n}(z_{part})-m_{med}\xi_{n}(z_{med})\dot{\xi_{n}}(z_{part})}{m_{part}\dot{\psi_{n}}(z_{med})\xi_{n}(z_{part})-m_{med}\psi_{n}(z_{med})\dot{\xi_{n}}(z_{part})}\\\
&\\\
b_{n}&:=\frac{m_{part}\xi_{n}(z_{med})\dot{\xi_{n}}(z_{part})-m_{med}\dot{\xi_{n}}(z_{med})\xi_{n}(z_{part})}{m_{part}\psi_{n}(z_{med})\dot{\xi_{n}}(z_{part})-m_{med}\dot{\psi_{n}}(z_{med})\xi_{n}(z_{part})}\\\
&\\\
c_{n}&:=\frac{m_{part}\psi_{n}(z_{med})\dot{\xi_{n}}(z_{med})-m_{part}\dot{\psi_{n}}(z_{med})\xi_{n}(z_{med})}{m_{part}\psi_{n}(z_{med})\dot{\xi_{n}}(z_{part})-m_{med}\dot{\psi_{n}}(z_{med})\xi_{n}(z_{part})}\\\
&\\\
d_{n}&:=\frac{m_{part}\dot{\psi_{n}}(z_{med})\xi_{n}(z_{med})-m_{part}\psi_{n}(z_{med})\dot{\xi_{n}}(z_{med})}{m_{part}\dot{\psi_{n}}(z_{med})\xi_{n}(z_{part})-m_{med}\psi_{n}(z_{med})\dot{\xi_{n}}(z_{part})}\end{split}$
(2.3)

With the Mie coefficients we can express the coefficient functions

$\displaystyle A_{n}(\rho,m_{med},m_{part})$ $\displaystyle:=\frac{l}{2\pi
m_{part}}\left(\left|c_{n}\right|^{2}\xi_{n}(z_{part})\overline{\dot{\xi_{n}}(z_{part})}-\left|d_{n}\right|^{2}\dot{\xi_{n}}(z_{part})\overline{\xi_{n}(z_{part})}\right)$
(2.4) and $\displaystyle B_{n}(\rho,m_{med},m_{part})$
$\displaystyle:=\frac{l}{2\pi
m_{med}}\left(\left|a_{n}\right|^{2}\dot{\psi_{n}}(z_{med})\overline{\psi_{n}(z_{med})}-\left|b_{n}\right|^{2}\psi_{n}(z_{med})\overline{\dot{\psi_{n}}(z_{med})}\right),$
(2.5)

which finally occur in the series expansion

$Q_{ext}(r,l,m_{med},m_{part})=\frac{l}{2cI(r,l)}\sum_{n=1}^{\infty}(2n+1)\mathrm{Im}\big{(}A_{n}(\rho,m_{med},m_{part})+B_{n}(\rho,m_{med},m_{part})\big{)}.$
(2.6)

Here the quantity $I(r,l)$ is the average incident intensity of light with
wavelength $l$ for a spherical particle with radius $r$ and $c$ denotes the
speed of light in vacuum. The function $I(r,l)$ is given by

$\begin{matrix}I(r,l)=\displaystyle{\frac{l^{2}}{8\pi(k_{med})^{2}}\frac{n_{med}}{2c}}\bigg{(}1+\left(4\pi
k_{med}\frac{r}{l}-1\right)e^{\displaystyle{4\pi
k_{med}\frac{r}{l}}}\bigg{)},&\;\mathrm{for}\;\;k_{med}\neq 0\\\ &\\\
I(r,l)=\pi
r^{2}\displaystyle{\frac{n_{med}}{2c}},&\;\;\mathrm{for}\;\;k_{med}=0.\end{matrix}$
(2.7)

Obviously we cannot evaluate (2.6) exactly, because we cannot compute an
infinite sum due to limited processing resources. Therefore we have to
truncate this series expansion. In [4] a commonly used truncation index
$N_{trunc}$ is presented, which is given by

$\begin{split}N_{trunc}&=\big{\lceil}\left|M\right.+4.05\cdot
M^{\frac{1}{3}}+\left.2\right|\big{\rceil},\\\
\mathrm{with}\;\;M&=\max\lceil\left|\rho\right|,|\rho\cdot m_{med}|,|\rho\cdot
m_{part}|\rceil.\end{split}$ (2.8)

##  3 Derivatives of the Truncated Mie Efficiency Series

The Bessel functions $J_{\alpha}(z)$ and $Y_{\alpha}(z)$ for an arbitrary
weight $\alpha$ fulfill the recurrence relations

$\frac{d}{dz}\big{(}z^{\alpha}J_{\alpha}(z)\big{)}=z^{\alpha}J_{\alpha-1}(z)\quad\text{and}\quad\frac{d}{dz}\big{(}z^{\alpha}Y_{\alpha}(z)\big{)}=z^{\alpha}Y_{\alpha-1}(z),$
(3.1)

cf. [5]. For the Bessel functions occurring in the Riccati-Bessel-functions
$\xi_{n}(z)$ and $\psi_{n}(z)$ follows from this with the weight
$\alpha=n+\frac{1}{2}$ that

$\displaystyle\dot{\xi_{n}}(z)$
$\displaystyle=\sqrt{\textstyle{\frac{\pi}{2}}}\sqrt{z}\left(J_{n-\frac{1}{2}}(z)-\frac{n}{z}J_{n+\frac{1}{2}}(z)\right)$
(3.2) $\displaystyle\text{and}\quad\dot{\psi_{n}}(z)$
$\displaystyle=\sqrt{\textstyle{\frac{\pi}{2}}}\sqrt{z}\left(J_{n-\frac{1}{2}}(z)-\frac{n}{z}J_{n+\frac{1}{2}}(z)\right)+\sqrt{\textstyle{\frac{\pi}{2}}}\sqrt{z}\left(Y_{n-\frac{1}{2}}(z)-\frac{n}{z}Y_{n+\frac{1}{2}}(z)\right)i$
(3.3)

for $z\neq 0$.

We apply (3.1) a second time to get $\ddot{\xi_{n}}(z)$, which yields

$\ddot{\xi_{n}}(z)=\sqrt{\textstyle{\frac{\pi}{2}}}\frac{\sqrt{z}}{z^{2}}\left(n(n+1)J_{n+\frac{1}{2}}(z)+(1-2n)J_{n-\frac{1}{2}}(z)+z^{2}J_{n-\frac{3}{2}}(z)\right).$

For an arbitrary weight $\alpha$ we have the recurrence relation

$J_{\alpha-1}(z)=\frac{2\alpha}{z}J_{\alpha}(z)-J_{\alpha+1}(z),$

see [5], and we use it to eliminate the term $J_{n-\frac{3}{2}}(z)$ in the
expression for $\ddot{\xi_{n}}(z)$. Then also $J_{n-\frac{1}{2}}(z)$ cancels
out, such that we obtain the representation

$\ddot{\xi_{n}}(z)=\sqrt{\textstyle{\frac{\pi}{2}}}\frac{\sqrt{z}}{z^{2}}\left(n(n+1)-z^{2}\right)J_{n+\frac{1}{2}}(z)$
(3.4)

only involving $J_{n+\frac{1}{2}}(z)$.

We recapitulate the Cauchy-Riemann equations in its complex form. For a
holomorphic function $f:\mathbb{C}\rightarrow\mathbb{C}$ with
$f(z)=f(x+iy)=u(x,y)+iv(x,y)$ holds

$\dot{f}(z)=\frac{d}{dz}f(z)=\frac{\partial}{\partial
x}f(x+iy)=-i\frac{\partial}{\partial y}f(x+iy).$ (3.5)

From this follows

$u_{x}=\mathrm{Re}\big{(}\dot{f}(z)\big{)},\quad
u_{y}=-\mathrm{Im}\big{(}\dot{f}(z)\big{)},\quad
v_{x}=\mathrm{Im}\big{(}\dot{f}(z)\big{)}\quad\text{and}\quad
v_{y}=\mathrm{Re}\big{(}\dot{f}(z)\big{)}.$ (3.6)

Using the latter relations we can deduce using
$z_{part}=\rho\left(n_{part}+ik_{part}\right)$

$\displaystyle\frac{\partial}{\partial
n_{part}}\mathrm{Re}\big{(}\xi_{n}(z_{part})\big{)}$
$\displaystyle=\rho\mathrm{Re}\big{(}\dot{\xi_{n}}(z_{part})\big{)},$
$\displaystyle\frac{\partial}{\partial
k_{part}}\mathrm{Re}\big{(}\xi_{n}(z_{part})\big{)}$
$\displaystyle=-\rho\mathrm{Im}\big{(}\dot{\xi_{n}}(z_{part})\big{)},$
$\displaystyle\frac{\partial}{\partial
n_{part}}\mathrm{Im}\big{(}\xi_{n}(z_{part})\big{)}$
$\displaystyle=\rho\mathrm{Im}\big{(}\dot{\xi_{n}}(z_{part})\big{)},$
$\displaystyle\frac{\partial}{\partial
k_{part}}\mathrm{Im}\big{(}\xi_{n}(z_{part})\big{)}$
$\displaystyle=\rho\mathrm{Re}\big{(}\dot{\xi_{n}}(z_{part})\big{)}$

and analogously

$\displaystyle\frac{\partial}{\partial
n_{part}}\mathrm{Re}\big{(}\dot{\xi_{n}}(z_{part})\big{)}$
$\displaystyle=\rho\mathrm{Re}\big{(}\ddot{\xi_{n}}(z_{part})\big{)},$
$\displaystyle\frac{\partial}{\partial
k_{part}}\mathrm{Re}\big{(}\dot{\xi_{n}}(z_{part})\big{)}$
$\displaystyle=-\rho\mathrm{Im}\big{(}\ddot{\xi_{n}}(z_{part})\big{)},$
$\displaystyle\frac{\partial}{\partial
n_{part}}\mathrm{Im}\big{(}\dot{\xi}_{n}(z_{part})\big{)}$
$\displaystyle=\rho\mathrm{Im}\big{(}\ddot{\xi_{n}}(z_{part})\big{)},$
$\displaystyle\frac{\partial}{\partial
k_{part}}\mathrm{Im}\big{(}\dot{\xi}_{n}(z_{part})\big{)}$
$\displaystyle=\rho\mathrm{Re}\big{(}\ddot{\xi_{n}}(z_{part})\big{)}.$

The Bessel function values

$\displaystyle
J_{0+\frac{1}{2}}(z_{med}),...,J_{N_{trunc}+\frac{1}{2}}(z_{med}),\quad
Y_{0+\frac{1}{2}}(z_{med}),...,Y_{N_{trunc}+\frac{1}{2}}(z_{med})$ and
$\displaystyle
J_{0+\frac{1}{2}}(z_{part}),...,J_{N_{trunc}+\frac{1}{2}}(z_{part}).$

already computed for a function evaluation of the truncated Mie extinction
efficiency can be reused for their derivatives.

Now everything is prepared to differentiate the squared magnitude of the Mie
coefficient $a_{n}$ with respect to $n_{part}$ and $k_{part}$. It is
sufficient only to discuss $\left|a_{n}\right|^{2}$ in the following, because
the differentiation of $\left|b_{n}\right|^{2}$, $\left|c_{n}\right|^{2}$ and
$\left|d_{n}\right|^{2}$ works analogously.

First we write the squared norm of the Mie coefficient $a_{n}$ as
$\left|a_{n}\right|^{2}=a_{n}\overline{a_{n}}$, which gives

$\displaystyle\frac{\partial}{\partial n_{part}}\left|a_{n}\right|^{2}$
$\displaystyle=\left(\frac{\partial}{\partial
n_{part}}a_{n}\right)\overline{a_{n}}+a_{n}\left(\frac{\partial}{\partial
n_{part}}\overline{a_{n}}\right)$
$\displaystyle\text{and}\quad\frac{\partial}{\partial
k_{part}}\left|a_{n}\right|^{2}$ $\displaystyle=\left(\frac{\partial}{\partial
k_{part}}a_{n}\right)\overline{a_{n}}+a_{n}\left(\frac{\partial}{\partial
k_{part}}\overline{a_{n}}\right).$

We write

$\displaystyle a_{n}=\frac{E}{D}\quad\text{with}\quad E$
$\displaystyle:=m_{part}\dot{\xi_{n}}(z_{med})\xi_{n}(z_{part})-m_{med}\xi_{n}(z_{med})\dot{\xi_{n}}(z_{part})$
$\displaystyle\text{and}\quad D$
$\displaystyle:=m_{part}\dot{\psi_{n}}(z_{med})\xi_{n}(z_{part})-m_{med}\psi_{n}(z_{med})\dot{\xi_{n}}(z_{part}),$

which yields

$\displaystyle\frac{d}{dm_{part}}a_{n}$
$\displaystyle=\frac{1}{D^{2}}\left(\left(\frac{d}{dm_{part}}E\right)D-E\left(\frac{d}{dm_{part}}D\right)\right)$
$\displaystyle\text{with}\quad\frac{d}{dm_{part}}E$
$\displaystyle=\dot{\xi_{n}}(z_{med})\xi_{n}(z_{part})+m_{part}\dot{\xi_{n}}(z_{med})\rho\dot{\xi_{n}}(z_{part})-m_{med}\xi_{n}(z_{med})\rho\ddot{\xi_{n}}(z_{part})$
$\displaystyle\text{and}\quad\frac{d}{dm_{part}}D$
$\displaystyle=\dot{\psi_{n}}(z_{med})\xi_{n}(z_{part})+m_{part}\dot{\psi_{n}}(z_{med})\rho\dot{\xi_{n}}(z_{part})-m_{med}\psi_{n}(z_{med})\rho\ddot{\xi_{n}}(z_{part}).$

Furthermore follow from (3.5) the relations

$\displaystyle\frac{\partial}{\partial n_{part}}a_{n}$
$\displaystyle=\frac{d}{dm_{part}}a_{n}$
$\displaystyle\text{and}\quad\frac{\partial}{\partial k_{part}}a_{n}$
$\displaystyle=\left(\frac{d}{dm_{part}}a_{n}\right)i.$

Although $\overline{a_{n}}$ is not holomorphic with respect to $m_{part}$, we
can still compute the partial derivatives
$\displaystyle{\frac{\partial}{\partial n_{part}}\overline{a_{n}}}$ and
$\displaystyle{\frac{\partial}{\partial k_{part}}\overline{a_{n}}}$. We obtain
using (3.6) the relations

$\displaystyle\frac{\partial}{\partial n_{part}}\overline{a_{n}}$
$\displaystyle=\overline{\frac{\partial}{\partial n_{part}}a_{n}}$
$\displaystyle\text{and}\quad\frac{\partial}{\partial
k_{part}}\overline{a_{n}}$ $\displaystyle=\overline{\frac{\partial}{\partial
k_{part}}a_{n}}.$

This completes the computations of $\displaystyle{\frac{\partial}{\partial
n_{part}}\left|a_{n}\right|^{2}}$ and $\displaystyle{\frac{\partial}{\partial
k_{part}}\left|a_{n}\right|^{2}}$.

For a holomorphic function $f(x+iy)$ we can easily deduce from (3.6) that

$\displaystyle\frac{\partial}{\partial
x}\mathrm{Im}\left(f(x+iy)\right)=\mathrm{Im}\left(\frac{\partial}{\partial
x}f(x+iy)\right)\quad\text{and}\quad\frac{\partial}{\partial
y}\mathrm{Im}\left(f(x+iy)\right)=\mathrm{Im}\left(\frac{\partial}{\partial
y}f(x+iy)\right).$

This gives with respect to (3.5)

$\displaystyle\frac{\partial}{\partial n_{part}}\mathrm{Im}\left(A_{n}\right)$
$\displaystyle=\mathrm{Im}\left(\frac{\partial}{\partial
n_{part}}A_{n}\right)$ $\displaystyle\text{with}\quad\frac{\partial}{\partial
n_{part}}A_{n}$ $\displaystyle=\frac{l}{2\pi}\left(\frac{\partial}{\partial
n_{part}}\left(\left|c_{n}\right|^{2}\right)U_{1}+\left|c_{n}\right|^{2}\frac{\partial}{\partial
n_{part}}U_{1}\right.$ $\displaystyle\hskip
142.26378pt\left.-\frac{\partial}{\partial
n_{part}}\left(\left|d_{n}\right|^{2}\right)U_{2}-\left|d_{n}\right|^{2}\frac{\partial}{\partial
n_{part}}U_{2}\right),$ $\displaystyle\text{where}\quad U_{1}$
$\displaystyle=\frac{\xi_{n}(z_{part})\overline{\dot{\xi_{n}}(z_{part})}}{m_{part}},$
$\displaystyle\frac{\partial}{\partial n_{part}}U_{1}$
$\displaystyle=\frac{1}{m_{part}^{2}}\left(\left(\rho\dot{\xi_{n}}(z_{part})\overline{\dot{\xi_{n}}(z_{part})}+\xi_{n}(z_{part})\rho\overline{\ddot{\xi_{n}}(z_{part})}\right)m_{part}\right.$
$\displaystyle\hskip
227.62204pt\left.-\xi_{n}(z_{part})\overline{\dot{\xi_{n}}(z_{part})}\right),$
$\displaystyle\text{and}\quad U_{2}$
$\displaystyle=\frac{\dot{\xi}_{n}(z_{part})\overline{\xi_{n}(z_{part})}}{m_{part}},$
$\displaystyle\frac{\partial}{\partial n_{part}}U_{2}$
$\displaystyle=\frac{1}{m_{part}^{2}}\left(\left(\rho\ddot{\xi_{n}}(z_{part})\overline{\xi_{n}(z_{part})}+\dot{\xi}_{n}(z_{part})\rho\overline{\dot{\xi_{n}}(z_{part})}\right)m_{part}\right.$
$\displaystyle\hskip
227.62204pt\left.-\dot{\xi}_{n}(z_{part})\overline{\xi_{n}(z_{part})}\right).$

Analogously we obtain

$\displaystyle\frac{\partial}{\partial k_{part}}\mathrm{Im}\left(A_{n}\right)$
$\displaystyle=\mathrm{Im}\left(\frac{\partial}{\partial
k_{part}}A_{n}\right)$ $\displaystyle\text{with}\quad\frac{\partial}{\partial
k_{part}}A_{n}$ $\displaystyle=\frac{l}{2\pi}\left(\frac{\partial}{\partial
k_{part}}\left(\left|c_{n}\right|^{2}\right)U_{1}+\left|c_{n}\right|^{2}\frac{\partial}{\partial
k_{part}}U_{1}\right.$ $\displaystyle\hskip
142.26378pt\left.-\frac{\partial}{\partial
n_{part}}\left(\left|d_{n}\right|^{2}\right)U_{2}-\left|d_{n}\right|^{2}\frac{\partial}{\partial
k_{part}}U_{2}\right),$
$\displaystyle\text{where}\quad\frac{\partial}{\partial k_{part}}U_{1}$
$\displaystyle=\frac{1}{m_{part}^{2}}\bigg{(}\left(\rho\left(\dot{\xi_{n}}(z_{part})\right)i\>\>\overline{\dot{\xi_{n}}(z_{part})}+\xi_{n}(z_{part})\rho\overline{\left(\ddot{\xi_{n}}(z_{part})\right)i}\right)m_{part}$
$\displaystyle\hskip
213.39566pt-\left(\xi_{n}(z_{part})\overline{\dot{\xi_{n}}(z_{part})}\right)i\bigg{)},$
$\displaystyle\text{and}\quad\frac{\partial}{\partial k_{part}}U_{2}$
$\displaystyle=\frac{1}{m_{part}^{2}}\bigg{(}\left(\rho\left(\ddot{\xi_{n}}(z_{part})\right)i\>\>\overline{\xi_{n}(z_{part})}+\dot{\xi}_{n}(z_{part})\rho\overline{\left(\dot{\xi_{n}}(z_{part})\right)i}\right)m_{part}$
$\displaystyle\hskip
213.39566pt-\left(\dot{\xi}_{n}(z_{part})\overline{\xi_{n}(z_{part})}\right)i\bigg{)}.$

The derivatives of $\mathrm{Im}\left(B_{n}\right)$ are much easier to compute,
since the dependence on $n_{part}$ and $k_{part}$ lies only in
$\left|a_{n}\right|^{2}$ and $\left|b_{n}\right|^{2}$ here. So we can deduce

$\displaystyle\frac{\partial}{\partial n_{part}}\mathrm{Im}\left(B_{n}\right)$
$\displaystyle=\mathrm{Im}\left(\frac{\partial}{\partial
n_{part}}B_{n}\right)$ $\displaystyle\text{with}\quad\frac{\partial}{\partial
n_{part}}B_{n}$ $\displaystyle=\frac{l}{2\pi}\left(\frac{\partial}{\partial
n_{part}}\left(\left|a_{n}\right|^{2}\right)\frac{\dot{\psi_{n}}(z_{med})\overline{\psi_{n}(z_{med})}}{m_{med}}\right.$
$\displaystyle\hskip 142.26378pt\left.-\frac{\partial}{\partial
n_{part}}\left(\left|b_{n}\right|^{2}\right)\frac{\psi_{n}(z_{med})\overline{\dot{\psi_{n}}(z_{med})}}{m_{med}}\right).$

Analogously we get

$\displaystyle\frac{\partial}{\partial k_{part}}\mathrm{Im}\left(B_{n}\right)$
$\displaystyle=\mathrm{Im}\left(\frac{\partial}{\partial
k_{part}}B_{n}\right)$ $\displaystyle\text{with}\quad\frac{\partial}{\partial
k_{part}}B_{n}$ $\displaystyle=\frac{l}{2\pi}\left(\frac{\partial}{\partial
k_{part}}\left(\left|a_{n}\right|^{2}\right)\frac{\dot{\psi_{n}}(z_{med})\overline{\psi_{n}(z_{med})}}{m_{med}}\right.$
$\displaystyle\hskip 142.26378pt\left.-\frac{\partial}{\partial
k_{part}}\left(\left|b_{n}\right|^{2}\right)\frac{\psi_{n}(z_{med})\overline{\dot{\psi_{n}}(z_{med})}}{m_{med}}\right).$

##  4 Nonlinear Regression using Truncated Series Expansions

We wish to reconstruct the refractive indices of a particle material from
spectral measurements by solving a nonlinear regression problem of the form

$X_{t,\delta}:=\underset{\boldsymbol{x}\in\mathbb{R}^{D}}{\mathrm{argmin}}\sum_{i=1}^{N}\frac{1}{2\sigma_{i}^{2}}\left(\sum_{n=1}^{t}a_{n}^{i}(\boldsymbol{x})-\sum_{n=1}^{\infty}a_{n}^{i}(\boldsymbol{x}_{true})-\delta_{i}\right)^{2},$
(4.1)

where $t\in\mathbb{N}$ is a finite truncation index and
$\delta_{i}\sim\mathcal{N}(0,s_{i}^{2})$. Remember that $N$ represents the
number of particle radii $r_{i}$ of the different monodisperse aerosols we are
investigating. We still assume that for each radius $r_{i}$ the standard
deviations $s_{i}$ are determined well enough from a set of experiments, such
that they can be regarded as known. We have to confine ourselves to a finite
truncation index $t$, because it is practically not feasible to compute all
coefficient functions $a_{n}^{i}(\boldsymbol{x})$ for $i=1,...,N$. Throughout
this paper we assume that the feasible set $\Omega$ is compact.

We define the functions
$\boldsymbol{f}_{t}:\mathbb{R}^{D}\rightarrow\mathbb{R}^{N}$ and
$\boldsymbol{f}:\mathbb{R}^{D}\rightarrow\mathbb{R}^{N}$ by

$\displaystyle\boldsymbol{f}_{t}(\boldsymbol{x})$
$\displaystyle:=\left(\sum_{n=1}^{t}a^{1}_{n}(\boldsymbol{x}),\;...,\;\sum_{n=1}^{t}a^{N}_{n}(\boldsymbol{x})\right)^{T}$
$\displaystyle\text{and}\quad\boldsymbol{f}(\boldsymbol{x})$
$\displaystyle:=\left(\sum_{n=1}^{\infty}a^{1}_{n}(\boldsymbol{x}),\;...,\;\sum_{n=1}^{\infty}a^{N}_{n}(\boldsymbol{x})\right)^{T}.$

We set
$\boldsymbol{e}:=\boldsymbol{f}(\boldsymbol{x}_{true})+\boldsymbol{\delta}$
with $\boldsymbol{\delta}:=\left(\delta_{1},...,\delta_{N}\right)^{T}$. Then
the observed probability density is given by

$p_{observed}(\boldsymbol{e}|\boldsymbol{x}):=(2\pi)^{-\frac{N_{l}}{2}}\big{|}\mathrm{det}(\boldsymbol{\Sigma_{\sigma}})\big{|}^{-\frac{1}{2}}\exp(-\textstyle{\frac{1}{2}}\|\boldsymbol{\Sigma_{\sigma}}^{-\frac{1}{2}}(\boldsymbol{f}_{t}(\boldsymbol{x})-\boldsymbol{e})\|_{2}^{2})$

with the covariance matrix
$\boldsymbol{\Sigma_{\sigma}}:=\mathrm{diag}\left(\sigma_{1}^{2},...,\sigma_{N}^{2}\right)$.
We know a priori that the vector $\boldsymbol{x}$ specifying our model
$\boldsymbol{f}_{t}(\boldsymbol{x})$ lies within the set $\Omega$. This
knowledge can be expressed with the prior probability density

$p_{prior}(\boldsymbol{x}):=\left(\mathrm{vol}(\Omega)\right)^{-1}I_{\Omega}(\boldsymbol{x}),$

where $I_{\Omega}$ is the indicator function of $\Omega$. Now $X_{t,\delta}$
is the set of MAP-estimators of the posterior probability density, i.e.

$\begin{split}X_{t,\delta}:=\underset{\boldsymbol{x}}{\mathrm{argmax}}&\;p_{posterior}(\boldsymbol{x}|\boldsymbol{e})\\\
\text{with}\quad p_{posterior}(\boldsymbol{x}|\boldsymbol{e})\propto
p_{observed}(\boldsymbol{e}|\boldsymbol{x})p_{prior}(\boldsymbol{x})&\propto\exp(-\textstyle{\frac{1}{2}}\|\boldsymbol{\Sigma_{\sigma}}^{-\frac{1}{2}}(\boldsymbol{f}_{t}(\boldsymbol{x})-\boldsymbol{e})\|_{2}^{2})I_{\Omega}(\boldsymbol{x}).\end{split}$
(4.2)

We carry out all the following investigations under the next assumption on the
covariance matrix:

######  Assumption 4.1.

The covariance matrix $\boldsymbol{\Sigma_{\sigma}}$ has the simple form

$\boldsymbol{\Sigma_{\sigma}}=\delta^{2}\cdot\mathrm{diag}(\sigma_{1}^{2},...,\sigma_{N}^{2})=:\delta^{2}\cdot\boldsymbol{\Sigma},$

where $\delta\geq 0$ is an arbitrary but fixed noise level and $\sigma_{1}$,
…, $\sigma_{N}$ are fixed.

To simplify notations we introduce the two functions
$\boldsymbol{f}_{t}:\mathbb{R}^{D}\rightarrow\mathbb{R}^{N}$ and
$\boldsymbol{g}_{t}:\mathbb{R}^{D}\rightarrow\mathbb{R}^{N}$ depending on the
truncation index $t$ and defined by

$\displaystyle\left(\boldsymbol{f}_{t}(\boldsymbol{x})\right)_{i}$
$\displaystyle:=\sum_{n=1}^{\lfloor
t\rfloor}a_{n}^{i}(\boldsymbol{x})+\big{(}t-\lfloor t\rfloor\big{)}a_{\lfloor
t\rfloor+1}^{i}(\boldsymbol{x})$
$\displaystyle\text{and}\quad\left(\boldsymbol{g}_{t}(\boldsymbol{x})\right)_{i}$
$\displaystyle:=\left(\boldsymbol{f}(\boldsymbol{x})\right)_{i}-\left(\boldsymbol{f}_{t}(\boldsymbol{x})\right)_{i},\quad\text{for}\quad
i=1,...,N.$

In the following we will investigate how an element
$\boldsymbol{x}_{t,\delta}$ of the set $X_{t,\delta}$ depends on the
truncation index $t$. We change to a continuous truncation index here, i.e. we
change from now on from (4.1) to the new regression problem

$\begin{split}X_{t,\delta}&:=\underset{\boldsymbol{x}\in\mathbb{R}^{D}}{\mathrm{argmin}}\;F_{t,\delta}(\boldsymbol{x})\quad\text{s.t.}\quad\boldsymbol{x}\in\Omega,\\\
\text{with}\quad
F_{t,\delta}(\boldsymbol{x})&:=\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t}(\boldsymbol{x})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}^{2}\end{split}$
(4.3)

where the truncation index $t\geq 0$ is allowed to be non-integer.

As a preparation we prove the following technical lemma, which will form the
basis of our continuity and convergence results.

######  Lemma 4.2.

Let the twice continuously differentiable function
$F:\mathbb{R}^{N}\rightarrow\mathbb{R}$ have a strict local minimum
$\boldsymbol{x}_{0}$ inside a compact set $S\subset\mathbb{R}^{N}$. Let the
function $h:\mathbb{R}^{N}\times\mathbb{R}\rightarrow\mathbb{R}$ have the
property $\lim_{\varepsilon\to 0}h(\boldsymbol{x},\varepsilon)=0$ for all
$\boldsymbol{x}\in S$ and let $h(\boldsymbol{x},\varepsilon)$ be twice
continuously differentiable with respect to $\boldsymbol{x}$ and continuous in
$\varepsilon$. Furthermore we assume that the local minima
$\boldsymbol{x}_{\varepsilon}$ of
$F_{\varepsilon}(\boldsymbol{x}):=F(\boldsymbol{x})+h(\boldsymbol{x},\varepsilon)$
are strict for any $\varepsilon>0$. Then there exists a sequence of local
minima $\boldsymbol{x}_{\varepsilon}$ of $F_{\varepsilon}(\boldsymbol{x})$
with $\lim_{\varepsilon\to 0}\boldsymbol{x}_{\varepsilon}=\boldsymbol{x}_{0}$.

######  Proof.

The strategy of the proof is to construct for given $\varepsilon$ a
neighborhood of $\boldsymbol{x}_{0}$ which must contain a local minimizer
$\boldsymbol{x}_{\varepsilon}$ of the perturbed function
$F_{\varepsilon}(\boldsymbol{x})$. By sending $\varepsilon$ to $0$, this
neighborhood shrinks down to the local minimum $\boldsymbol{x}_{0}$ itself,
thus yielding the convergence of $\boldsymbol{x}_{\varepsilon}$ to
$\boldsymbol{x}_{0}$. To have this neighborhood shrink down to
$\boldsymbol{x}_{0}$, it is crucially important that $\boldsymbol{x}_{0}$ must
be a strict local minimum.

We define $d(\varepsilon):=\sup_{\boldsymbol{x}\in
S}|h(\boldsymbol{x},\varepsilon)|$. From $\lim_{\varepsilon\to
0}h(\boldsymbol{x},\varepsilon)=0$ for all $\boldsymbol{x}\in S$ follows
$\lim_{\varepsilon\to 0}d(\varepsilon)=0$. Let us now introduce the function
$F^{-}(\boldsymbol{x}):=F(\boldsymbol{x})-d(\varepsilon)$. Obviously
$\boldsymbol{x}_{0}$ is also a local minimum of $F^{-}(\boldsymbol{x})$, so
for $\varepsilon$ sufficiently small there exists a neighborhood
$U_{2d(\varepsilon)}(\boldsymbol{x}_{0})\subset S$ of $\boldsymbol{x}_{0}$
with

$F^{-}(\boldsymbol{x})\geq F^{-}(\boldsymbol{x}_{0})\quad\text{and}\quad
F^{-}(\boldsymbol{x})-F^{-}(\boldsymbol{x}_{0})\leq
2d(\varepsilon)\quad\text{for all}\quad\boldsymbol{x}\in
U_{2d(\varepsilon)}(\boldsymbol{x}_{0}).$

In particular we have

$\forall\boldsymbol{x}\in\partial
U_{2d(\varepsilon)}(\boldsymbol{x}_{0}):\;F^{-}(\boldsymbol{x})=F^{-}(\boldsymbol{x}_{0})+2d(\varepsilon)=F(\boldsymbol{x}_{0})+d(\varepsilon).$

Let us assume that there exists an $\boldsymbol{x}\in\partial
U_{2d(\varepsilon)}(\boldsymbol{x}_{0})$ with

$F_{\varepsilon}(\boldsymbol{x})<F(\boldsymbol{x}_{0})+d(\varepsilon)=F^{-}(\boldsymbol{x})=F(\boldsymbol{x})-d(\varepsilon).$

Then
$F_{\varepsilon}(\boldsymbol{x})=F(\boldsymbol{x})+h(\boldsymbol{x},\varepsilon)$
implies $-d(\varepsilon)>h(\boldsymbol{x},\varepsilon)$, hence
$-h(\boldsymbol{x},\varepsilon)>d(\varepsilon)\geq-h(\boldsymbol{x},\varepsilon)$
by definition of $d(\varepsilon)$, contradiction. Therefore we conclude

$\forall\boldsymbol{x}\in\partial
U_{2d(\varepsilon)}(\boldsymbol{x}_{0}):\;F_{\varepsilon}(\boldsymbol{x})\geq
F(\boldsymbol{x}_{0})+d(\varepsilon).$ (4.4)

Since $F_{\varepsilon}(\boldsymbol{x})$ is continuous and
$\overline{U}_{2d(\varepsilon)}(\boldsymbol{x}_{0})$ is compact for
$\varepsilon$ small enough, there exists an
$\boldsymbol{x}_{\varepsilon}\in\overline{U}_{2d(\varepsilon)}(\boldsymbol{x}_{0})$
with

$F_{\varepsilon}(\boldsymbol{x}_{\varepsilon})=\min_{\boldsymbol{x}\in\overline{U}_{2d(\varepsilon)}(\boldsymbol{x}_{0})}F_{\varepsilon}(\boldsymbol{x}).$

Let us assume
$F_{\varepsilon}(\boldsymbol{x}_{\varepsilon})>F(\boldsymbol{x}_{0})+d(\varepsilon)$.
Then by definition of $\boldsymbol{x}_{\varepsilon}$ we get in particular

$\displaystyle
F(\boldsymbol{x}_{0})+h(\boldsymbol{x}_{0},\varepsilon)=F_{\varepsilon}(\boldsymbol{x}_{0})\geq
F_{\varepsilon}(\boldsymbol{x}_{\varepsilon})>F(\boldsymbol{x}_{0})+d(\varepsilon),$

i.e. $h(\boldsymbol{x}_{0},\varepsilon)>d(\varepsilon)\geq
h(\boldsymbol{x}_{0},\varepsilon)$, contradiction. It follows

$F_{\varepsilon}(\boldsymbol{x}_{\varepsilon})\leq
F(\boldsymbol{x}_{0})+d(\varepsilon)\quad\text{and}\quad
F_{\varepsilon}(\boldsymbol{x}_{0})\leq F(\boldsymbol{x}_{0})+d(\varepsilon),$
(4.5)

where the latter follows with a proof by contradiction as well.

If it happens to hold that
$F_{\varepsilon}(\boldsymbol{x}_{\varepsilon})=F(\boldsymbol{x}_{0})+d(\varepsilon)$,
then we also have
$F_{\varepsilon}(\boldsymbol{x}_{0})=F(\boldsymbol{x}_{0})+d(\varepsilon)$.
Otherwise we have
$F_{\varepsilon}(\boldsymbol{x}_{\varepsilon})<F(\boldsymbol{x}_{0})+d(\varepsilon)$
and then (4.4) implies that $\boldsymbol{x}_{\varepsilon}$ cannot lie on
$\partial U_{2d(\varepsilon)}(\boldsymbol{x}_{0})$, thus it must lie within
the interior of $U_{2d(\varepsilon)}(\boldsymbol{x}_{0})$. So in any case
(4.5) gives that $U_{2d(\varepsilon)}(\boldsymbol{x}_{0})$ must contain a
local minimizer $\boldsymbol{x}_{\varepsilon}$ of
$F_{\varepsilon}(\boldsymbol{x})$.

Now $\lim_{\varepsilon\to 0}d(\varepsilon)=0$ gives $\lim_{\varepsilon\to
0}\boldsymbol{x}_{\varepsilon}=\boldsymbol{x}_{0}$. The existence of the last
limit is guaranteed by the fact that $\boldsymbol{x}_{0}$ is strict and the
claim is proved.

$\square$

######  Proposition 4.3.

Let all coefficient functions $a_{n}^{i}(\boldsymbol{x})$ be twice
continuously differentiable and bounded on $\Omega$. We assume that each local
minimum $\boldsymbol{x}_{t,\delta}$ of the right hand side function
$F_{t,\delta}(\boldsymbol{x})$ in (4.3) is strict and lies in the interior of
$\Omega$. Then each local minimum depends continuously on the truncation index
$t$.

######  Proof.

To prove the claim, one could be tempted to apply the implicit function
theorem on the equation $\boldsymbol{d}(t,\boldsymbol{x}_{t,\delta})=0$ with
$\boldsymbol{d}(s,\boldsymbol{x}):=\nabla\boldsymbol{f}_{s}(\boldsymbol{x})$.
This would give that the local minima are parameterized by a function
$\boldsymbol{m}(s)$ with the property
$\boldsymbol{m}(t)=\boldsymbol{x}_{t,\delta}$, where $s$ is from an
environment $U(t)$ of $t$. The problem with this approach is that it requires
continuous differentiability of $\boldsymbol{d}(s,\boldsymbol{x})$ in the
truncation parameter $s$. Thus the continuous truncation we are using would
need more complicated methods such as spline interpolation of the partial
sums, which would increase the overall computational effort.

Therefore we use in the following a more direct approach to prove the claim.
Let $\varepsilon>0$ be arbitrary. First we consider an integer truncation
index $t\in\mathbb{N}$, i.e. we have $t=\lfloor t\rfloor$. Now for
$\varepsilon$ small enough, we get $\lfloor t+\varepsilon\rfloor=t$ and
$\lfloor t-\varepsilon\rfloor=t-1$. This gives

$\displaystyle\left(\boldsymbol{f}_{t+\varepsilon}(\boldsymbol{x})\right)_{i}$
$\displaystyle=\left(\boldsymbol{f}_{t}(\boldsymbol{x})\right)_{i}+\varepsilon
a_{t+1}^{i}(\boldsymbol{x})$
$\displaystyle\text{and}\quad\left(\boldsymbol{f}_{t-\varepsilon}(\boldsymbol{x})\right)_{i}$
$\displaystyle=\left(\boldsymbol{f}_{t-1}(\boldsymbol{x})\right)_{i}+(1-\varepsilon)a_{t}^{i}(\boldsymbol{x})$
$\displaystyle=\left(\boldsymbol{f}_{t}(\boldsymbol{x})\right)_{i}-\varepsilon
a_{t}^{i}(\boldsymbol{x}).$

As next step we turn to an noninteger truncation index $t$. In this case, we
can always select $\varepsilon$ small enough such that $\lfloor
t+\varepsilon\rfloor=\lfloor t\rfloor$ and $\lfloor
t-\varepsilon\rfloor=\lfloor t\rfloor$ respectively hold. This yields

$\displaystyle\left(\boldsymbol{f}_{t+\varepsilon}(\boldsymbol{x})\right)_{i}$
$\displaystyle=\left(\boldsymbol{f}_{t}(\boldsymbol{x})\right)_{i}+\varepsilon
a_{\lfloor t\rfloor+1}^{i}(\boldsymbol{x})$
$\displaystyle\text{and}\quad\left(\boldsymbol{f}_{t-\varepsilon}(\boldsymbol{x})\right)_{i}$
$\displaystyle=\left(\boldsymbol{f}_{t}(\boldsymbol{x})\right)_{i}-\varepsilon
a_{\lfloor t\rfloor+1}^{i}(\boldsymbol{x}).$

Now we introduce the function

$\displaystyle\boldsymbol{a}(\boldsymbol{x}):=\begin{cases}\big{(}a_{t}^{1}(\boldsymbol{x}),...,a_{t}^{N}(\boldsymbol{x})\big{)}^{T},&\text{for}\quad
t-\varepsilon,\;t\in\mathbb{N}\\\ \big{(}a_{\lfloor
t\rfloor+1}^{1}(\boldsymbol{x}),...,a_{\lfloor
t\rfloor+1}^{N}(\boldsymbol{x})\big{)}^{T},&\text{else}.\end{cases}$

For
$F_{t,\delta}(\boldsymbol{x})=\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t}(\boldsymbol{x})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}^{2}$
this yields

$\displaystyle F_{t+\varepsilon,\delta}(\boldsymbol{x})$
$\displaystyle=F_{t,\delta}(\boldsymbol{x})+2\varepsilon\big{\langle}\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{a}(\boldsymbol{x}),\;\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t}(\boldsymbol{x})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\big{\rangle}+\varepsilon^{2}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{a}(\boldsymbol{x})\|_{2}^{2}$
$\displaystyle\text{and}\quad F_{t-\varepsilon,\delta}(\boldsymbol{x})$
$\displaystyle=F_{t,\delta}(\boldsymbol{x})-2\varepsilon\big{\langle}\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{a}(\boldsymbol{x}),\;\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t}(\boldsymbol{x})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\big{\rangle}+\varepsilon^{2}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{a}(\boldsymbol{x})\|_{2}^{2}.$

Therefore we obtain both for $F_{t+\varepsilon,\delta}(\boldsymbol{x})$ and
$F_{t-\varepsilon,\delta}(\boldsymbol{x})$ a decomposition of the form
$F_{t+\varepsilon,\delta}(\boldsymbol{x})=F_{t,\delta}(\boldsymbol{x})+h^{\varepsilon}_{t,\delta}(\boldsymbol{x})$
and
$F_{t-\varepsilon,\delta}(\boldsymbol{x})=F_{t,\delta}(\boldsymbol{x})+h^{\varepsilon}_{t,\delta}(\boldsymbol{x})$
respectively, where the function $h^{\varepsilon}_{t,\delta}(\boldsymbol{x})$
is appropriately selected according to above findings. We can readily check
$\lim_{\varepsilon\to 0}|h^{\varepsilon}_{t,\delta}(\boldsymbol{x})|=0$ for
all $\boldsymbol{x}\in\Omega$ from the boundedness of the
$a_{n}^{i}(\boldsymbol{x})$’s. Then the result follows from Lemma 4.2.

$\square$

######  Corollary 4.4.

Let $t_{1}$ and $t_{2}$ be truncation indices with $t_{1}<t_{2}$. Let
$\boldsymbol{x}_{t_{1},\delta}$ be a local minimizer of (4.1). Let
$\gamma\in[0,1]$ and define
$t_{\gamma}:=t_{1}+\gamma\left(t_{2}-t_{1}\right)$. Then beginning at
$\gamma=0$ one can successively find local minimizers
$\boldsymbol{x}_{t_{\gamma},\delta}$ for the truncation index $t_{\gamma}$
using numerical continuation, see [6]. Here for $\gamma_{1}<\gamma_{2}$ the
minimizer $\boldsymbol{x}_{t_{\gamma_{1}},\delta}$ is used as a start vector
to compute the next minimizer $\boldsymbol{x}_{t_{\gamma_{2}},\delta}$. The
next parameter $\gamma_{2}$ has to be sufficiently close to $\gamma_{1}$, such
that the start vector $\boldsymbol{x}_{t_{\gamma_{1}},\delta}$ still lies
within the domain of convergence for Newton’s method. $\square$

We use Corollary 4.4 to compute $\boldsymbol{x}_{t,\delta}$ for increasing
truncation index $t$ in a stable way. If we would keep $t$ as integer and
increase it in integer steps, we might leave the domain of convergence in the
continuation method. Therefore we increase them using a smaller step width.

In the following we investigate how well the minimizers
$\boldsymbol{x}_{t,\delta}$ of the noise-contaminated regression problem (4.3)
with truncated series expansions approximate the minimizers
$\boldsymbol{x}_{\infty,0}$ of the noise-free and untruncated problem

$\boldsymbol{x}_{\infty,0}:=\underset{\boldsymbol{x}\in\mathbb{R}^{D}}{\mathrm{argmin}}\sum_{i=1}^{N}\frac{1}{2\sigma_{i}^{2}}\left(\sum_{n=1}^{\infty}a_{n}^{i}(\boldsymbol{x})-\sum_{n=1}^{\infty}a_{n}^{i}(\boldsymbol{x}_{true})\right)^{2}\quad\text{s.t.}\quad\boldsymbol{x}\in\Omega.$
(4.6)

######  Proposition 4.5.

Let the noise vector $\boldsymbol{\delta}$ fulfill $\lim_{\delta\to
0}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{\delta}\|_{2}^{2}=0$ and let
the functions $\boldsymbol{f}(\boldsymbol{x})$ and
$\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x})$ be bounded on $\Omega$. Assume
$\lim_{\delta\to
0}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x})\|_{2}^{2}=0$
for all $\boldsymbol{x}\in\Omega$. Then for any strict minimizer
$\boldsymbol{x}_{\infty,0}$ of the right hand side function of (4.6) in the
interior of $\Omega$ exist minimizers $\boldsymbol{x}_{t_{\delta},\delta}$ of
(4.3) with $\lim_{\delta\to
0}\boldsymbol{x}_{t_{\delta},\delta}=\boldsymbol{x}_{\infty,0}$. Here we also
assume the $\boldsymbol{x}_{t_{\delta},\delta}$’s to be strict for all
$\delta>0$.

######  Proof.

With the notation introduced before we can write

$\displaystyle\boldsymbol{x}_{\infty,0}\in X_{\infty,0}$
$\displaystyle:=\underset{\boldsymbol{x}\in\mathbb{R}^{D}}{\mathrm{argmin}}\;F_{\infty,0}(\boldsymbol{x})\quad\text{s.t.}\quad\boldsymbol{x}\in\Omega$
$\displaystyle\text{with}\quad F_{\infty,0}(\boldsymbol{x})$
$\displaystyle:=\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}(\boldsymbol{x})-\boldsymbol{f}(\boldsymbol{x}_{true}))\|_{2}^{2}.$

From the decomposition
$\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x})=\boldsymbol{f}(\boldsymbol{x})\>-\>\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x})$
we obtain

$F_{t_{\delta},\delta}(\boldsymbol{x})=\;\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\left(\boldsymbol{f}(\boldsymbol{x})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta}\right)\|_{2}^{2}-2\big{\langle}\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x}),\;\boldsymbol{\Sigma}^{-\frac{1}{2}}\left(\boldsymbol{f}(\boldsymbol{x})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta}\right)\big{\rangle}\\\
+\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x})\|_{2}^{2}.$

Then a further decomposition of the first term on the right hand side yields

$\displaystyle F_{t_{\delta},\delta}(\boldsymbol{x})=\;$ $\displaystyle
F_{\infty,0}(\boldsymbol{x})-2\big{\langle}\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{\delta},\;\boldsymbol{\Sigma}^{-\frac{1}{2}}\left(\boldsymbol{f}(\boldsymbol{x})-\boldsymbol{f}(\boldsymbol{x}_{true})\right)\big{\rangle}+\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{\delta}\|_{2}^{2}$
$\displaystyle-$ $\displaystyle
2\big{\langle}\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x}),\;\boldsymbol{\Sigma}^{-\frac{1}{2}}\left(\boldsymbol{f}(\boldsymbol{x})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta}\right)\big{\rangle}+\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x})\|_{2}^{2}$
$\displaystyle=:$
$\displaystyle\;F_{\infty,0}(\boldsymbol{x})+H_{t_{\delta},\delta}(\boldsymbol{x}).$

From the limit $\lim_{\delta\to
0}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x})\|_{2}^{2}=0$,
the limit $\lim_{\delta\to
0}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{\delta}\|_{2}^{2}=0$ and the
boundedness of $\boldsymbol{f}(\boldsymbol{x})$ follows $\lim_{\delta\to
0}|H_{t_{\delta},\delta}(\boldsymbol{x})|=0$ for arbitrary but fixed
$\boldsymbol{x}\in\Omega$.

Then the existence of the $\boldsymbol{x}_{t_{\delta},\delta}$’s follows from
Lemma 4.2.

$\square$

At last we study how the minimizers $\boldsymbol{x}_{t_{\delta},\delta}$ of
(4.3) behave for $\delta\rightarrow 0$. We begin with a preparing corollary.

######  Corollary 4.6.

Let the assumptions of Proposition 4.5 hold. Then we have for any local
minimizer $\boldsymbol{x}_{t_{\delta},\delta}$ of (4.3) approximating a local
minimizer $\boldsymbol{x}_{\infty,0}$ of (4.6) for $\delta\rightarrow 0$ with
$\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}(\boldsymbol{x}_{\infty,0})-\boldsymbol{f}(\boldsymbol{x}_{true}))\|_{2}=0$
that

$\lim_{\delta\to
0}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}=0.$

######  Proof.

The assumptions of Proposition 4.5 give

$\lim_{\delta\to
0}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{\delta}\|_{2}^{2}=0\quad\text{and}\quad\lim_{\delta\to
0}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta})\|_{2}^{2}=0.$

We have by continuity of $\boldsymbol{f}(\boldsymbol{x})$ that
$\lim_{\delta\to
0}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}(\boldsymbol{x}_{t_{\delta},\delta})-\boldsymbol{f}(\boldsymbol{x}_{true}))\|_{2}=0$.
Then

$\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}(\boldsymbol{x}_{t_{\delta},\delta})-\boldsymbol{f}(\boldsymbol{x}_{true}))\|_{2}+\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta})\|_{2}+\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{\delta}\|_{2}\\\
\geq\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}$

gives the desired result. $\square$

######  Proposition 4.7.

Let the assumptions of Proposition 4.5 hold. Assume that the local minimizers
$\boldsymbol{x}_{\infty,0}$ of (4.6) with
$\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}(\boldsymbol{x}_{\infty,0})-\boldsymbol{f}(\boldsymbol{x}_{true}))\|_{2}=0$
form a discrete set $S_{\infty,0}$. Then the set $L_{\infty,0}$ consisting of
the limits $\lim_{\delta\to 0}\boldsymbol{x}_{t_{\delta},\delta}$ of local
minimizers $\boldsymbol{x}_{t_{\delta},\delta}$ of (4.3) with $\lim_{\delta\to
0}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}=0$
coincides with $S_{\infty,0}$ and there exists a noise level $\delta_{max}$
such that all minimizers $\boldsymbol{x}_{t_{\delta},\delta}$ approximating
$S_{\infty,0}$ are isolated for all $\delta\leq\delta_{max}$.

######  Proof.

On the one hand from Proposition 4.5 we know that there exists a sequence
$\boldsymbol{x}_{t_{\delta},\delta}$ of minimizers of (4.3) with
$\lim_{\delta\to
0}\boldsymbol{x}_{t_{\delta},\delta}=\boldsymbol{x}_{\infty,0}$. Then
$\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}(\boldsymbol{x}_{\infty,0})-\boldsymbol{f}(\boldsymbol{x}_{true}))\|_{2}=0$
and Corollary 4.6 give $\lim_{\delta\to
0}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}=0$,
which implies $S_{\infty,0}\subseteq L_{\infty,0}$.

On the other hand holds for $\boldsymbol{x}_{t_{\delta},\delta}$ with
$\lim_{\delta\to
0}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}=0$
that

$\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}+\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta})\|_{2}+\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{\delta}\|_{2}\\\
\geq\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}(\boldsymbol{x}_{t_{\delta},\delta})-\boldsymbol{f}(\boldsymbol{x}_{true}))\|_{2}$

which implies $\lim_{\delta\to
0}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}(\boldsymbol{x}_{t_{\delta},\delta})-\boldsymbol{f}(\boldsymbol{x}_{true}))\|_{2}=0$.
In particular this means by continuity of $\boldsymbol{f}(\boldsymbol{x})$
that the vector $\lim_{\delta\to 0}\boldsymbol{x}_{t_{\delta},\delta}$ must be
a local minimizer of
$\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}(\boldsymbol{x})-\boldsymbol{f}(\boldsymbol{x}_{true}))\|_{2}$.
Thus we have also shown $L_{\infty,0}\subseteq S_{\infty,0}$.

In the following we number all elements of $S_{\infty,0}$ with the index $k$,
i.e. we write $\boldsymbol{x}_{\infty,0}^{k}$ for $k=1,...,|S_{\infty,0}|$.
Similarly we number all minimizers $\boldsymbol{x}_{t_{\delta},\delta}$ with
$\lim_{\delta\to
0}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}=0$
approximating the $\boldsymbol{x}_{\infty,0}^{k}$’s with
$\boldsymbol{x}_{t_{\delta},\delta}^{k}$, i.e. $\lim_{\delta\to
0}\boldsymbol{x}_{t_{\delta},\delta}^{k}=\boldsymbol{x}_{\infty,0}^{k}$ for
$k=1,...,|S_{\infty,0}|$. Define

$D_{min}:=\min_{i\neq
j}\|\boldsymbol{x}_{\infty,0}^{i}-\boldsymbol{x}_{\infty,0}^{j}\|_{2}.$

Since $\lim_{\delta\to
0}\boldsymbol{x}_{t_{\delta},\delta}^{k}=\boldsymbol{x}_{\infty,0}^{k}$, we
can find an error levels $\delta_{max}^{k}$ such that

$\|\boldsymbol{x}_{t_{\delta},\delta}^{k}-\boldsymbol{x}_{\infty,0}^{k}\|_{2}<\textstyle{\frac{1}{2}}D_{min}\quad\text{for}\quad
k=1,...,|S_{\infty,0}|,$

which holds for all $0\leq\delta\leq\delta_{max}^{k}$ for each $k$. Then for
all $0\leq\delta\leq\delta_{max}:=\min_{k}\\{\delta_{max}^{k}\\}$ the
$\boldsymbol{x}_{t_{\delta},\delta}^{k}$’s must have pairwise mutual distances
greater than zero.

$\square$

Now Proposition 4.7 gives that the number of local minima
$\boldsymbol{x}_{t_{\delta},\delta}^{k}$ remains constant if the noise level
$\delta$ is small enough. It also yields that these local minima then form a
set of separated continuous curves parametrized in $\delta$.

At last we wish to have an estimate of the convergence of the local minima
$\boldsymbol{x}_{t_{\delta},\delta}^{k}$ of the truncated and noise
contaminated problem to the local minima $\boldsymbol{x}_{\infty,0}^{k}$ of
the noise-free and untruncated problem, which is useful for practical
computations.

######  Proposition 4.8.

For the noise level $\delta$ small enough, we can bound for any local minimum
$\boldsymbol{x}_{t_{\delta},\delta}^{k}$ the approximation error
$\|\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\|_{2}$
with a positively weighted linear combination of the residual
$\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta}^{k})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}$,
the truncation error
$\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}$
and the noise estimate
$\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{\delta}\|_{2}$.

######  Proof.

The first order necessary conditions for a local minimum of
$F_{t_{\delta},\delta}(\boldsymbol{x})=F_{\infty,0}(\boldsymbol{x})+H_{t_{\delta},\delta}(\boldsymbol{x})$
at $\boldsymbol{x}_{t_{\delta},\delta}^{k}$ and a local minimum of
$F_{\infty,0}(\boldsymbol{x})$ at $\boldsymbol{x}_{\infty,0}^{k}$ yield in
particular

$\displaystyle\big{\langle}\nabla
F_{\infty,0}(\boldsymbol{x}_{t_{\delta},\delta}^{k})+\nabla
H_{t_{\delta},\delta}(\boldsymbol{x}_{t_{\delta},\delta}^{k}),\;\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\big{\rangle}$
$\displaystyle\geq 0$ $\displaystyle\text{and}\quad\big{\langle}\nabla
F_{\infty,0}(\boldsymbol{x}_{\infty,0}^{k}),\;\boldsymbol{x}_{t_{\delta},\delta}^{k}-\boldsymbol{x}_{\infty,0}^{k}\big{\rangle}$
$\displaystyle\geq 0,$

Adding the last two inequalities yields

$\big{\langle}\nabla
H_{t_{\delta},\delta}(\boldsymbol{x}_{t_{\delta},\delta}^{k}),\;\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\big{\rangle}\geq\big{\langle}\nabla
F_{\infty,0}(\boldsymbol{x}_{\infty,0}^{k})-\nabla
F_{\infty,0}(\boldsymbol{x}_{t_{\delta},\delta}^{k}),\;\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\big{\rangle}.$
(4.7)

Since $\nabla F_{\infty,0}(\boldsymbol{x})$ is totally differentiable at
$\boldsymbol{x}_{\infty,0}^{k}$ we obtain

$\nabla F_{\infty,0}(\boldsymbol{x}_{t_{\delta},\delta}^{k})=\nabla
F_{\infty,0}(\boldsymbol{x}_{\infty,0}^{k})+\mathrm{Hess}_{F_{\infty,0}}(\boldsymbol{x}_{\infty,0}^{k})\left(\boldsymbol{x}_{t_{\delta},\delta}^{k}-\boldsymbol{x}_{\infty,0}^{k}\right)+\boldsymbol{w}_{\infty,0}(\boldsymbol{x}_{t_{\delta},\delta}^{k},\boldsymbol{x}_{\infty,0}^{k}),$

where
$\boldsymbol{w}_{\infty,0}(\boldsymbol{x},\boldsymbol{x}_{\infty,0}^{k})$
fulfills

$\displaystyle\|\boldsymbol{w}_{\infty,0}(\boldsymbol{x},\boldsymbol{x}_{\infty,0}^{k})\|_{2}\leq\|\boldsymbol{x}-\boldsymbol{x}_{\infty,0}^{k}\|_{2}\epsilon_{\infty,0}(\boldsymbol{x},\boldsymbol{x}_{\infty,0}^{k})\quad\text{with}\quad\lim_{\boldsymbol{x}\to\boldsymbol{x}_{\infty,0}^{k}}\epsilon_{\infty,0}(\boldsymbol{x},\boldsymbol{x}_{\infty,0}^{k})=0.$

Since $\mathrm{Hess}_{F_{\infty,0}}(\boldsymbol{x}_{\infty,0}^{k})$ is
positive definite, the expression
$\left(\big{\langle}\boldsymbol{x},\mathrm{Hess}_{F_{\infty,0}}(\boldsymbol{x}_{\infty,0}^{k})\>\boldsymbol{x}\big{\rangle}\right)^{\frac{1}{2}}$
gives a norm on $\mathbb{R}^{D}$. Because of the equivalence of all norms in
$\mathbb{R}^{D}$, there exists a constant $C_{\infty,0}^{k}>0$ with

$\displaystyle\left(\big{\langle}\boldsymbol{x},\mathrm{Hess}_{F_{\infty,0}}(\boldsymbol{x}_{\infty,0}^{k})\>\boldsymbol{x}\big{\rangle}\right)^{\frac{1}{2}}\geq
C_{\infty,0}^{k}\|\boldsymbol{x}\|_{2}\quad\text{for
all}\quad\boldsymbol{x}\in\mathbb{R}^{D}.$

Since $\lim_{\delta\to
0}\boldsymbol{x}_{t_{\delta},\delta}^{k}=\boldsymbol{x}_{\infty,0}^{k}$ we can
find a noise level $\rho_{max}^{k}$ such that
$|\epsilon_{\infty,0}(\boldsymbol{x}_{t_{\delta},\delta}^{k},\boldsymbol{x}_{\infty,0}^{k})|\leq
d_{\infty,0}^{k}$ for all $\delta\leq\rho_{max}^{k}$, where $d_{\infty,0}^{k}$
is a constant with $0\leq d_{\infty,0}^{k}<\left(C_{\infty,0}^{k}\right)^{2}$.
Then using

$\big{\langle}\boldsymbol{w}_{\infty,0}(\boldsymbol{x}_{t_{\delta},\delta}^{k},\boldsymbol{x}_{\infty,0}^{k}),\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\big{\rangle}\leq\epsilon_{\infty,0}(\boldsymbol{x}_{t_{\delta},\delta}^{k},\boldsymbol{x}_{\infty,0}^{k})\|\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\|_{2}^{2}$

and (4.7) we can estimate

$\displaystyle\|\nabla
H_{t_{\delta},\delta}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}\|\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\|_{2}$
$\displaystyle\geq\;$ $\displaystyle\big{\langle}\nabla
F_{\infty,0}(\boldsymbol{x}_{\infty,0}^{k})-\nabla
F_{\infty,0}(\boldsymbol{x}_{t_{\delta},\delta}^{k}),\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\big{\rangle}$
$\displaystyle=\;$
$\displaystyle\big{\langle}\mathrm{Hess}_{F_{\infty,0}}(\boldsymbol{x}_{\infty,0}^{k})\left(\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\right)-\boldsymbol{w}_{\infty,0}(\boldsymbol{x}_{t_{\delta},\delta}^{k},\boldsymbol{x}_{\infty,0}^{k}),\;\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\big{\rangle}$
$\displaystyle\geq\;$
$\displaystyle\left(\big{(}C_{\infty,0}^{k}\big{)}^{2}-\epsilon_{\infty,0}(\boldsymbol{x}_{t_{\delta},\delta}^{k},\boldsymbol{x}_{\infty,0}^{k})\right)\|\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\|_{2}^{2},$

i.e. this gives

$\|\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\|_{2}\leq\left(\big{(}C_{\infty,0}^{k}\big{)}^{2}-d_{\infty,0}^{k}\right)^{-1}\|\nabla
H_{t_{\delta},\delta}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}$ (4.8)

for all $\delta\leq\rho_{max}^{k}$.

We have

$\displaystyle\nabla H_{t_{\delta},\delta}(\boldsymbol{x})=\;2\Big{(}$
$\displaystyle\mathrm{Jac}_{\boldsymbol{g}_{t_{\delta}}}^{T}(\boldsymbol{x})\boldsymbol{\Sigma}^{-1}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x})-\mathrm{Jac}_{\boldsymbol{g}_{t_{\delta}}}^{T}(\boldsymbol{x})\boldsymbol{\Sigma}^{-1}(\boldsymbol{f}(\boldsymbol{x})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})$
$\displaystyle-$
$\displaystyle\mathrm{Jac}_{\boldsymbol{f}}^{T}(\boldsymbol{x})\boldsymbol{\Sigma}^{-1}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x})-\mathrm{Jac}_{\boldsymbol{f}}^{T}(\boldsymbol{x})\boldsymbol{\Sigma}^{-1}\boldsymbol{\delta}\Big{)},$

i.e. we can find constants $K_{1}^{k}$ and $K_{2}^{k}$ such that

$\displaystyle\|\nabla
H_{t_{\delta},\delta}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}\leq\;$
$\displaystyle
2\Big{(}K_{1}^{k}\|\mathrm{Jac}_{\boldsymbol{g}_{t_{\delta}}}^{T}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\boldsymbol{\Sigma}^{-\frac{1}{2}}\|_{\infty}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}$
$\displaystyle+K_{1}^{k}\|\mathrm{Jac}_{\boldsymbol{g}_{t_{\delta}}}^{T}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\boldsymbol{\Sigma}^{-\frac{1}{2}}\|_{\infty}\big{(}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta}^{k})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}$
$\displaystyle+\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}\big{)}$
$\displaystyle+K_{2}^{k}\|\mathrm{Jac}_{\boldsymbol{f}}^{T}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\boldsymbol{\Sigma}^{-\frac{1}{2}}\|_{\infty}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}$
$\displaystyle+K_{2}^{k}\|\mathrm{Jac}_{\boldsymbol{f}}^{T}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\boldsymbol{\Sigma}^{-\frac{1}{2}}\|_{\infty}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{\delta}\|_{2}\Big{)},$

which gives the result. $\square$

######  Corollary 4.9.

Let the truncation indices $t_{\delta}$ depend on the vector of independent
Gaussian random variables $\boldsymbol{\delta}$ with $\lim_{\delta\to
0}\mathbb{E}\big{(}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{\delta}\|_{2}^{2}\big{)}=0$
such that $\lim_{\delta\to
0}\mathbb{E}\big{(}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x})\|_{2}^{2}\big{)}=0$
holds for all arbitrary but fixed $\boldsymbol{x}\in\Omega$. Then we have for
all minimizers $\boldsymbol{x}_{\infty,0}^{k}$ of (4.6) with
$\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}(\boldsymbol{x}_{\infty,0}^{k})-\boldsymbol{f}(\boldsymbol{x}_{true}))\|_{2}=0$
that for $\delta$ sufficiently small there exist minimizers
$\boldsymbol{x}_{t_{\delta},\delta}^{k}$ of (4.3) with

$\lim_{\delta\to
0}\mathbb{E}\big{(}\|\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\|_{2}\big{)}=0.$

######  Proof.

Proposition 4.5 establishes the existence of the
$\boldsymbol{x}_{t_{\delta},\delta}^{k}$’s. Corollary 4.6 gives

$\displaystyle\lim_{\delta\to
0}\mathbb{E}\big{(}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta}^{k})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}\big{)}=0.$

Proposition 4.7 gives that for $\delta$ sufficiently small there exists a
constant $K_{\infty,0}^{k}$ with

$\mathbb{E}\big{(}\|\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\|_{2}\big{)}\leq
K_{\infty,0}^{k}\mathbb{E}\big{(}\|\nabla
H_{t_{\delta},\delta}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}\big{)}.$

Set
$S_{1}^{k}:=\sup_{\boldsymbol{x}\in\Omega}\|\mathrm{Jac}_{\boldsymbol{g}_{t_{\delta}}}^{T}(\boldsymbol{x})\boldsymbol{\Sigma}^{-\frac{1}{2}}\|_{\infty}<\infty$
and
$S_{2}^{k}:=\sup_{\boldsymbol{x}\in\Omega}\|\mathrm{Jac}_{\boldsymbol{f}}^{T}(\boldsymbol{x})\boldsymbol{\Sigma}^{-\frac{1}{2}}\|_{\infty}<\infty$.
Then the estimate for $\|\nabla
H_{t_{\delta},\delta}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}$ in the
proof of Proposition 4.7 gives

$\displaystyle\mathbb{E}\big{(}\|\nabla
H_{t_{\delta},\delta}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}\big{)}\leq\;$
$\displaystyle
2\Big{(}K_{1}^{k}S_{1}^{k}\mathbb{E}\big{(}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}\big{)}$
$\displaystyle+K_{1}^{k}S_{1}^{k}\Big{(}\mathbb{E}\big{(}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta}^{k})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}\big{)}$
$\displaystyle+\mathbb{E}\big{(}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}\big{)}\Big{)}$
$\displaystyle+K_{2}^{k}S_{2}^{k}\mathbb{E}\big{(}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}\big{)}$
$\displaystyle+K_{2}^{k}S_{2}^{k}\mathbb{E}\big{(}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{\delta}\|_{2}\big{)}\Big{)},$

which proves the claim since Assumption 4.1 gives
$\mathbb{E}\big{(}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{\delta}\|_{2}\big{)}\leq\sqrt{N}\delta$.

$\square$

The strategy for our retrieval algorithm is to start with an initial guess for
the truncation index $t_{start}$ and try to find all local minima
$\boldsymbol{x}_{t_{start},\delta}^{k}$. Then the truncation index is
gradually increased and starting from $\boldsymbol{x}_{t_{start},\delta}^{k}$
the continuation method is applied to find finally the local minima
$\boldsymbol{x}_{t_{\delta},\delta}^{k}$. Motivated by Propositions 4.8 and
4.9 only those local minima are considered to be possible approximations to
our sought-after refractive index, where the residual
$\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta}^{k})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}$
and an estimate of the truncation error
$\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\boldsymbol{g}_{t_{\delta}}(\boldsymbol{x}_{t_{\delta},\delta}^{k})\|_{2}$
are both reasonably small. The latter serves also as a stopping criterion for
the continuation method

The initial guess $t_{start}$ has to be selected with care. On the one hand if
it is to small, the model is to inaccurate and the retrieval of the sought-
after local minima can not be guaranteed. On the other hand if it is too big,
computational effort is wasted, since too many Mie coefficient functions with
almost vanishing magnitudes and thus essentially not changing the local minima
are computed.

##  5 The Reconstruction Algorithm

We now return to our regression problem (1.3). For $i=1,...,N$ we see that the
measured extinctions normalized by the number of particles $n_{i}$ with radius
$r_{i}$, i.e. the quantity $\frac{e_{i}}{n_{i}}$, is Gaussian-distributed with
mean $\frac{1}{n_{i}}q_{true}(r_{i},l)$ and standard deviation
$\sigma_{i}:=\frac{s_{i}}{n_{i}}$. In the following we fix a wavelength $l$,
i.e. we reconstruct the sought-after particle refractive index $m_{part}(l)$
wavelength by wavelength. In the following the unit both for particle radii
and light wavelengths is $\mu$m.

We make use of the function
$\boldsymbol{q}_{N_{tr}}:\mathbb{R}^{2}\rightarrow\mathbb{R}^{N}$ defined by

$\boldsymbol{q}_{N_{tr}}(\boldsymbol{x}):=\left(\pi
r_{1}^{2}\sum_{n=1}^{N_{tr}}q_{n}(\boldsymbol{x},r_{1},l),\;...,\;\pi
r_{N}^{2}\sum_{n=1}^{N_{tr}}q_{n}(\boldsymbol{x},r_{N},l)\right)^{T},$ (5.1)

where $R=\\{r_{1},...,r_{N}\\}$ is the particle radius grid and $N_{tr}$ the
truncation index to be used. We allow non-integer truncation indices $N_{tr}$
as well, where the non-integer truncation is done like in Proposition 4.3.
Here the expression $q_{n}(\boldsymbol{x},r_{k},l)$ is a short notation of
$q_{n}(m_{med}(l),(\boldsymbol{x})_{1}+(\boldsymbol{x})_{2}i,r_{k},l)$ from
Section 1, where the sought-after refrative index $m_{part}(l)$ is identified
with the vector $\boldsymbol{x}$ here, i.e.
$m_{part}(l)=(\boldsymbol{x})_{1}+(\boldsymbol{x})_{2}i$. So its computation
follows Section 2.

In the following the refractive index search area is given by the rectangle
$\Omega:=[0,20]\times[0,40]$, which means that we only consider refractive
indices of particle materials whose real parts lie in the interval $[0,20]$
and its imaginary parts in the interval $[0,40]$. This rather large search
area makes the algorithm suitable for a wide range of aerosol materials.

Algorithm 1 Reconstruction of Refractive Indices

1:$b_{real}=20$

2:$b_{imag}=40$

3:$N_{real}=81$

4:$N_{imag}=161$

5:$c_{i}=(i-1)\frac{b_{real}}{N_{real}-1}$ for $i=1,...,N_{real}$

6:$d_{i}=(i-1)\frac{b_{imag}}{N_{imag}-1}$ for $i=1,...,N_{imag}$

7:$R=\\{0.1,0.2,0.3\\}$

8:$N=3$

9:$N_{tr}=3$

10:$S_{start}=\\{\\}$

11:estimate $\sigma_{1}^{2}$, …, $\sigma_{N}^{2}$ from the sample means
approximating the standard deviations of $\frac{e_{1}}{n_{1}}$, …,
$\frac{e_{N}}{n_{N}}$.

12:$\delta^{2}:=\mathrm{max}\big{\\{}\sigma_{1}^{2},...,\sigma_{N}^{2}\big{\\}}$

13:$\boldsymbol{\Sigma}:=\delta^{-2}\cdot\mathrm{diag}\big{(}\sigma_{1}^{2},...,\sigma_{N}^{2}\big{)}$

14:for $i=1\;\textbf{to}\;N_{real}$ do

15: for $j=1\;\textbf{to}\;N_{imag}$ do

16: compute the Hessian $\boldsymbol{H}(c_{i},d_{j})$ of
$F(\boldsymbol{x}):=\frac{1}{2}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\left(\boldsymbol{q}_{N_{tr}}(\boldsymbol{x})-\boldsymbol{e}_{real}\right)\|_{2}^{2}$
at $\boldsymbol{x}=(c_{i},d_{j})^{T}$

17: if $\boldsymbol{H}(c_{i},d_{j})$ is positive definite then

18: use $(c_{i},d_{j})^{T}$ as start vector to compute

19:
$\boldsymbol{x}_{new}=\underset{\boldsymbol{x}\in\mathbb{R}^{2}}{\mathrm{argmin}}\frac{1}{2}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\left(\boldsymbol{q}_{N_{tr}}(\boldsymbol{x})-\boldsymbol{e}\right)\|_{2}^{2}\quad\text{s.t.}\quad\boldsymbol{x}\in[0,b_{real}]\times[0,b_{imag}]$

20: if $S_{start}$ is empty $\lor$
$\frac{\|\boldsymbol{x}-\boldsymbol{x}_{new}\|_{2}}{\|\boldsymbol{x}\|_{2}}\geq
10^{-2}\quad\forall\boldsymbol{x}\in S_{start}$ then

21: $S_{start}=S_{start}\cup\\{\boldsymbol{x}_{new}\\}$

22: end if

23: end if

24: end for

25:end for

26:$N_{start}=\left|S_{start}\right|$

27:$S_{out}=\\{\\}$

28:$\tau=3$

29:for $i=1\;\textbf{to}\;N_{start}$ do

30: $c=N_{tr}$

31: $Tol_{rel}=10^{-3}$

32: $D_{rel}=\infty$

33: while $D_{rel}>Tol_{rel}$ do

34: for $p=1\;\textbf{to}\;10$ do

35: use the vector $S_{start}(i)$ as start vector to compute

36:
$\boldsymbol{x}_{new}=\underset{\boldsymbol{x}\in\mathbb{R}^{2}}{\mathrm{argmin}}\frac{1}{2}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\left(\boldsymbol{q}_{c+\frac{p}{10}}(\boldsymbol{x})-\boldsymbol{e}\right)\|_{2}^{2}\quad\text{s.t.}\quad\boldsymbol{x}\in[0,b_{real}]\times[0,b_{imag}]$

37:
$Res_{cur}=\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\left(\boldsymbol{q}_{c+\frac{p-1}{10}}(S_{start}(i))-\boldsymbol{e}\right)\|_{2}^{2}$

38:
$Res_{new}=\|\boldsymbol{\Sigma}^{-\frac{1}{2}}\left(\boldsymbol{q}_{c+\frac{p}{10}}(\boldsymbol{x}_{new})-\boldsymbol{e}\right)\|_{2}^{2}$

39: $D_{rel}=\frac{\left|Res_{cur}-Res_{new}\right|}{Res_{cur}}$

40: $S_{start}(i)=\boldsymbol{x}_{new}$

41: end for

42: $c=c+1$

43: end while

44: if $Res_{new}<\tau N\delta^{2}$ then

45: if $S_{out}$ is empty $\lor$
$\frac{\|\boldsymbol{x}-\boldsymbol{x}_{new}\|_{2}}{\|\boldsymbol{x}\|_{2}}\geq
10^{-2}\quad\forall\boldsymbol{x}\in S_{out}$ then

46: $S_{out}=S_{out}\cup\\{\boldsymbol{x}_{new}\\}$

47: end if

48: end if

49:end for

In the first loop from lines 14 \- 25 a search for local minima of the fit
function $F(\boldsymbol{x})$ defined in line 16 for the truncation index
$N_{tr}=3$ is performed. The loop runs through all grid points
$(c_{i},d_{j})^{T}$ of the search grid defined in lines 5 \- 6. If the Hessian
of $F(\boldsymbol{x})$ at some grid point $(c_{i},d_{j})^{T}$ is positive
definite, this point might lie in the vicinity of a local minimum. The Hessian
is computed exactly, where the second partial derivatives of the Mie
extinction efficiency with respect to the real and imaginary part of the
scattering material needed here are computed using the product rule approach
from Section 3. So we use $(c_{i},d_{j})^{T}$ as start point for a local
solver in this case. In line 20 we only accept a new local minimum if it is
sufficiently different from the local minima already found. Then it is stored
in the container $S_{start}$. This simple global search strategy can find all
local minima if the search grid is fine enough.

######  Remark 5.1.

Of course other well established global search heuristics can be applied here
as well. In test runs we compared genetic algorithms with our sequential
search strategy on the two-dimensional refractive index search area, but their
computational effort and reliability remained the same or were even worse. In
[7] the technique of simulated annealing was used to retrieve aerosol
refractive indices, which could be a promising alternative here. In our study
the measurement noise was so high, that a unique global minimum of our fit
function could not be determined. Instead our focus lied on effectively
finding all local minima with small values of the fit function and we regarded
them all as possible approximations to a thought-after refractive index.

The second loop from lines 29 \- 49 uses the local minima found in the first
loop as start points for the continuation method following Proposition 4.3 and
Corollary 4.4. We found that a step width of $0.1$ is for our problem a well-
balanced choice between too big step widths rendering the continuation method
unstable and too small step widths making it computationally inefiicient. With
the stopping criterion $D_{rel}\leq Tol_{rel}$ of the while-loop it is
approximately checked if the magnitude of the remainder term is small enough.
Finally in line 44 it is checked if the residual is small enough. In our
implementation we did another run of lines 44 \- 48 with $\tau=5$ and $\tau=7$
respectively, if none of the reconstructions had a squared residual smaller
than $\tau N_{r}\delta^{2}$ for the previous $\tau$. This had to be done,
because the parameter $\tau$ has to be selected carefully in order to estimate
the bound on
$\mathbb{E}\big{(}\|\boldsymbol{x}_{\infty,0}^{k}-\boldsymbol{x}_{t_{\delta},\delta}^{k}\|_{2}\big{)}$
derived in the proof of Corollary 4.9 correctly.

##  6 Comparison of the Numerical Continuation Approach with Established
Truncation Index Heuristics

As solution of the forward problem we generated for a discrete set of
wavelengths $l_{1}$, …, $l_{N_{l}}$ unperturbed spectral extinctions
normalized with the number of particles of the monodisperse aerosol by
computing

$\left(\boldsymbol{e}_{true}\right)_{i,j}:=\pi
r_{i}^{2}\sum_{n=1}^{N_{tr}}q_{n}(m_{med}(l_{j}),m_{part}(l_{j}),r_{i},l_{j}),\quad\text{for}\quad
i=1,...,N,\quad j=1,...,N_{l}$

with $m_{part}(l_{i})$ taken as the refractive indices of Ag,
$\mathrm{H}_{2}$O and CsI. Here we used the truncation index

$\begin{split}&\rho=2\pi\frac{r}{l},\\\ &M=\max\\{|\rho|,|\rho\cdot
m_{part}(l)|,|\rho\cdot m_{med}(l)|\\},\\\ &N_{tr}:=\lceil|M+4.05\cdot
M^{\frac{1}{3}}+2|\rceil\end{split}$ (6.1)

introduced in [4].

For particle size distribution reconstructions as outlined in [2] we need
particle refractive indices for five optical windows, see [8], so the
wavelength grid of interest consists of five ranges. These ranges are given by
$8$ linearly spaced wavelengths from $0.6-0.8\;\mu$m, $8$ from
$1.1-1.3\;\mu$m, $8$ from $1.6-1.8\;\mu$m, $16$ from $2.1-2.5\;\mu$m and $8$
from $3.1-3.3\;\mu$m, so we have in total $N_{l}=48$ wavelengths.

For each of the $48$ wavelengths we generated noisy spectral extinctions
$\boldsymbol{e}$ by adding zero-mean Gaussian noise to
$\boldsymbol{e}_{true}$, i.e.

$\left(\boldsymbol{e}\right)_{i,j}=\left(\boldsymbol{e}_{true}\right)_{i,j}+\delta_{i,j}\quad\text{
with
}\;\delta_{i,j}\sim\mathcal{N}(0,(0.05\cdot\left(\boldsymbol{e}_{true}\right)_{i,j})^{2}),\quad
i=1,...,N\quad j=1,...,N_{l}.$

Here the standard deviations were taken to be $5\%$ of the original extinction
values. We computed each mean $\left(\boldsymbol{e}_{real}\right)_{i,j}$ of
the noisy spectral extinctions with a sample size of $N_{s}=300$ .

In the following Algortihm 1 is referred to as method 1. On the same simulated
spectral extinctions we let Algorithm 1 run up to line 25, but with the
difference that at each evaluation of
$\boldsymbol{q}_{N_{tr}}(\boldsymbol{x})$ we directly took the trunction index
from (6.1). We denote this approach with method 2. We now display the average
run times of method 1 and method 2 for $10$ sweeps through all $48$
wavelenghs.

####  6.0.1 Results for Ag

####  6.0.2 results for CsI

####  6.0.3 Results for $\mathrm{H}_{2}$O

###  6.1 Maximal Relative Deviations

For the $10$ simulation runs we list the maximal relative deviations

$100\cdot\frac{\big{\|}\left(n_{part}^{1}(l_{j}),k_{part}^{1}(l_{j})\right)^{T}-\left(n_{part}^{2}(l_{j}),k_{part}^{2}(l_{j})\right)^{T}\big{\|}_{2}}{\big{\|}\left(n_{part}^{1}(l_{i}),k_{part}^{1}(l_{j})\right)^{T}\big{\|}_{2}}$

of the refractive index reconstructions
$\left(n_{part}^{2}(l_{j}),k_{part}^{2}(l_{j})\right)^{T}$ from method 2 from
$\left(n_{part}^{1}(l_{j}),k_{part}^{1}(l_{j})\right)^{T}$ of method 1 for
$j=1,...,48$. At each wavelength, multiple local minima can be detected by
both methods. For the relative deviations we always selected the local minima
forming the smoothest reconstructions on each optical window in the sense of
Section 8.

####  6.1.1 Results for Ag

####  6.1.2 results for CsI

####  6.1.3 Results for $\mathrm{H}_{2}$O

###  6.2 Conclusion

For Ag, the average total run time over all $48$ wavelengths for method 1 was
$44.0167\%$ less than for method 2, for $\mathrm{H}_{2}$O $43.9808\%$ and for
CsI $44.7322\%$, i.e. method 1 is almost two times faster than method 2. The
results are of the same quality, since their relative deviations are just
small fractions of percentages.

The continuation method approach saves run time significantly with the same
quality of the results compared to using the truncation index (6.1) all the
time.

##  7 Nonlinear Tikhonov Regularization

So far we have solved the regression problem (4.3) without any regularization,
thus the obtained refractive index reconstructions might still be too error-
contaminated to be of practical use. A widely used regularization strategy for
nonlinear regression problems is Tikhonov regularization, which yields the
regularized regression problem

$\boldsymbol{x}_{\gamma}:=\underset{\boldsymbol{x}\in\mathbb{R}^{D}}{\mathrm{argmin}}\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t}(\boldsymbol{x})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}^{2}+\gamma\|\boldsymbol{x}-\boldsymbol{x}^{*}\|_{2}^{2}\quad\text{s.t.}\quad\boldsymbol{x}\in\Omega$
(7.1)

when we apply it on (4.3), cf. [9]. Here $\gamma$ is a regularization
parameter and $\boldsymbol{x}^{*}$ is an estimate of the sought-after true
solution. In many cases the unregularized problem has a whole set of
minimizers, thus the vector $\boldsymbol{x}^{*}$ works also as a selection
criterion. Now if a reasonable $\boldsymbol{x}^{*}$ is found, the
regularization parameter $\gamma$ can be determined with the discrepancy
principle, i.e. $\gamma$ is computed such that

$\displaystyle\|\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t}(\boldsymbol{x}_{\gamma})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\|_{2}=R(\delta),$

where $R(\delta)$ is an estimate of the residual of the “true” solution which
depends on the noise level $\delta$. For this task monotonicity in the
residual of $\boldsymbol{x}_{\gamma}$ is established in [9].

The problem of finding a good estimate $\boldsymbol{x}^{*}$ still remains. In
[10] an alternative implementable parameter choice strategy without the need
of an $\boldsymbol{x}^{*}$ is derived. Applied on our problem it gives

$\displaystyle\gamma\big{\langle}\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t}(\boldsymbol{x}_{\gamma})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta}),\;\mathrm{J}_{\gamma}^{-1}(\boldsymbol{\Sigma}^{-\frac{1}{2}}(\boldsymbol{f}_{t}(\boldsymbol{x}_{\gamma})-\boldsymbol{f}(\boldsymbol{x}_{true})-\boldsymbol{\delta})\big{\rangle}=R(\delta),$
with
$\displaystyle\mathrm{J}_{\gamma}:=\gamma\mathrm{I}+\boldsymbol{\Sigma}^{-\frac{1}{2}}\mathrm{Jac}_{\boldsymbol{f}_{t}}(\boldsymbol{x}_{\gamma})\mathrm{Jac}_{\boldsymbol{f}_{t}}(\boldsymbol{x}_{\gamma})^{T}\boldsymbol{\Sigma}^{-\frac{1}{2}}.$

This method has the drawback that the matrix $\mathrm{J}_{\gamma}$ needs to be
inverted, which may lead to instabilities.

Nevertheless the quality of the regularized solutions still depends strongly
on the start values for solving (7.1). We know about our sought-after
refractive indices that they form smooth curves on each of the five optical
windows. The complex refractive index curves of most materials can be
described using the so-called Lorentz-oscillator-model, cf. [11]. Here points
with bigger curvatures only occur at so-called resonance frequencies
corresponding to some isolated resonance wavelengths. Motivated by these facts
we derive in the following a method to find reasonable start values for
Phillips-Twomey-regularization out of the results of Algorithm 1, which will
be outlined in Section 9.

##  8 Finding the Smoothest Coupled Solutions

We have the problem of identifying the best approximation to the sought-after
true particle material refractive index $\boldsymbol{x}_{true}$ out of a set
of multiple solutions obtained with Algorithm 1. This problem can be solved by
coupling the solutions, which means that we combine solutions from neighboring
wavelengths $l$ in each of the five optical windows in order to obtain a
unique solution for every optical window. We know about the complex refractive
index curves to be retrieved that they are smooth, hence we expect their sum
of the squared second finite differences both in the real and imaginary parts
to be small.

Let $l_{1}$, …, $l_{s}$ denote the wavelengths of any of our five wavelength
ranges. Let $N_{1}$, …, $N_{s}$ be the number of solutions found for all the
wavelengths. We denote with $\boldsymbol{x}^{i}_{j}$ the $j$-th solution found
for wavelength $l_{i}$ for $i=1,...,s$ and $j=1,...,N_{i}$. Now we wish to
find the smoothest combined solution from all possible combinations
$\boldsymbol{x}^{1}_{j_{1}}$, …, $\boldsymbol{x}^{s}_{j_{s}}$ for
$j_{i}=1,...,N_{i}$, hence we have a total number of $\prod_{i=1}^{s}N_{i}$
combinations. Here we measure smoothness of a combination
$\boldsymbol{x}^{1}_{j_{1}}$, …, $\boldsymbol{x}^{s}_{j_{s}}$ with the sum

$S:=\sum_{i=2}^{s-1}\left(\left(\big{(}\boldsymbol{x}^{i-1}_{j_{i-1}}\big{)}_{1}-2\big{(}\boldsymbol{x}^{i}_{j_{i}}\big{)}_{1}+\big{(}\boldsymbol{x}^{i+1}_{j_{i+1}}\big{)}_{1}\right)^{2}+\left(\big{(}\boldsymbol{x}^{i-1}_{j_{i-1}}\big{)}_{2}-2\big{(}\boldsymbol{x}^{i}_{j_{i}}\big{)}_{2}+\big{(}\boldsymbol{x}^{i+1}_{j_{i+1}}\big{)}_{2}\right)^{2}\right)$

of its second finite differences both in the real parts
$(\boldsymbol{x}^{i}_{j_{i}}\big{)}_{1}$ and its imaginary parts
$(\boldsymbol{x}^{i}_{j_{i}}\big{)}_{2}$, which means that we regard a
combination the smoother the smaller its sum $S$ is.

We encounter the problem that the total number of possible combinations
$\prod_{i=1}^{s}N_{i}$ might get too big to iterate through all combinations
in the search for the smoothest one in acceptable time. Therefore we propose a
greedy algorithm, which uses each second finite difference as start point to
find a smooth combination.

Algorithm 2 Detection of the Smoothest Combination

1:$S_{min}=\infty$

2:$S_{cur}=0$

3:$Comb=\\{\\}$

4:$SmoothestCombination=\\{\\}$

5:for $z=2\;\textbf{to}\;s-1$ do

6: for $c1=1\;\textbf{to}\;N_{z-1}$ do

7: for $c2=1\;\textbf{to}\;N_{z}$ do

8: for $c3=1\;\textbf{to}\;N_{z+1}$ do

9:
$S_{cur}=\Big{(}\big{(}\boldsymbol{x}^{z-1}_{c1}\big{)}_{1}-2\big{(}\boldsymbol{x}^{z}_{c2}\big{)}_{1}+\big{(}\boldsymbol{x}^{z+1}_{c3}\big{)}_{1}\Big{)}^{2}+\Big{(}\big{(}\boldsymbol{x}^{z-1}_{c1}\big{)}_{2}-2\big{(}\boldsymbol{x}^{z}_{c2}\big{)}_{2}+\big{(}\boldsymbol{x}^{z+1}_{c3}\big{)}_{2}\Big{)}^{2}$

10: $Comb(z-1)=\boldsymbol{x}^{z-1}_{c1}$

11: $Comb(z)=\boldsymbol{x}^{z}_{c2}$

12: $Comb(z+1)=\boldsymbol{x}^{z+1}_{c3}$

13: for $k=z-2\;\textbf{to}\;1$ do

14: $D_{min}=\infty$

15: $D_{cur}=\infty$

16: $\boldsymbol{x}_{min}=(0,0)^{T}$

17: $\boldsymbol{x}_{mid}=Comb(k+1)$

18: $\boldsymbol{x}_{right}=Comb(k+2)$

19: for $j=1\;\textbf{to}\;N_{k}$ do

20:
$D_{cur}=\Big{(}\big{(}\boldsymbol{x}^{k}_{j}\big{)}_{1}-2\big{(}\boldsymbol{x}_{mid})_{1}+\big{(}\boldsymbol{x}_{right}\big{)}_{1}\Big{)}^{2}\newline
\hskip
120.92421pt+\Big{(}\big{(}\boldsymbol{x}^{k}_{j}\big{)}_{2}-2\big{(}\boldsymbol{x}_{mid}\big{)}_{2}+\big{(}\boldsymbol{x}_{right}\big{)}_{2}\Big{)}^{2}$

21: if $D_{cur}<D_{min}$ then

22: $D_{min}=D_{cur}$

23: $\boldsymbol{x}_{min}=\boldsymbol{x}^{k}_{j}$

24: end if

25: end for

26: $S_{cur}=S_{cur}+D_{min}$

27: $Comb(k)=\boldsymbol{x}_{min}$

28: end for

29: for $k=z+2\;\textbf{to}\;s$ do

30: $D_{min}=\infty$

31: $D_{cur}=\infty$

32: $\boldsymbol{x}_{min}=(0,0)^{T}$

33: $\boldsymbol{x}_{mid}=Comb(k-1)$

34: $\boldsymbol{x}_{left}=Comb(k-2)$

35: for $j=1\;\textbf{to}\;N_{k}$ do

36:
$D_{cur}=\Big{(}\big{(}\boldsymbol{x}_{left}\big{)}_{1}-2\big{(}\boldsymbol{x}_{mid})_{1}+\big{(}\boldsymbol{x}^{k}_{j}\big{)}_{1}\Big{)}^{2}\newline
\hskip
120.92421pt+\Big{(}\big{(}\boldsymbol{x}_{left}\big{)}_{2}-2\big{(}\boldsymbol{x}_{mid}\big{)}_{2}+\big{(}\boldsymbol{x}^{k}_{j}\big{)}_{2}\Big{)}^{2}$

37: if $D_{cur}<D_{min}$ then

38: $D_{min}=D_{cur}$

39: $\boldsymbol{x}_{min}=\boldsymbol{x}^{k}_{j}$

40: end if

41: end for

42: $S_{cur}=S_{cur}+D_{min}$

43: $Comb(k)=\boldsymbol{x}_{min}$

44: end for

45: if $S_{cur}<S_{min}$ then

46: $S_{min}=S_{cur}$

47: $SmoothestCombination=Comb$

48: end if

49: end for

50: end for

51: end for

52:end for

The main loop spanning over the lines 5 \- 52 iterates through all positions
$z=2,...,s-1$ of a middle point for a second finite difference both in the
real and imaginary part. At each position $z$ the inner loops beginning in
lines 6 \- 8 iterate through all possible second finite differences which can
be formed out of the vectors $\boldsymbol{x}^{z-1}_{c1}$,
$\boldsymbol{x}^{z}_{c2}$ and $\boldsymbol{x}^{z+1}_{c3}$ for
$c1=1,...,N_{z-1}$, $c2=1,...,N_{z}$ and $c3=1,...,N_{z+1}$, i.e. they loop
through all of their middle points and left and right neighbors at position
$z$. In lines 9 \- 12 the variable $S_{cur}$ is initialized with the sum of
the squared second finite differences in the real and imaginary parts of the
current vectors $\boldsymbol{x}^{z-1}_{c1}$, $\boldsymbol{x}^{z}_{c2}$ and
$\boldsymbol{x}^{z+1}_{c3}$ and the positions $z-1$, $z$ and $z+1$ of the
array $Comb$ are filled with the current vectors. For $z\geq 3$ the loop in
lines 13 \- 28 successively fills the positions $k=z-2,z-3,...,1$ of the array
$Comb$. At each new position $k$ the minimal sum of the two squared second
finite differences in the real and imaginary part $D_{min}$ is determined in
lines 19 \- 25, where the middle and right point are fixed and taken as the
leftmost two vectors from the array $Comb$ and the right point runs through
all $\boldsymbol{x}^{k}_{j}$ for $j=1,...,N_{k}$. After the vector
$\boldsymbol{x}_{min}$ giving out of all of the $\boldsymbol{x}^{k}_{j}$’s
giving the minimal sum $D_{min}$ is found, the $D_{min}$ is added to $S_{cur}$
and $\boldsymbol{x}_{min}$ is stored in $k$-th entry $Comb(k)$. In a similar
way the loop in lines 29 \- 44 succesively fills the positions
$k=z+2,z+3,...,s$ for $z\leq s-2$. This time the left and middle point are
fixed and taken as the rightmost two points of the array $Comb$, whereas the
left point iterates through all $\boldsymbol{x}^{k}_{j}$ for $j=1,...,N_{k}$.
Again the vector $\boldsymbol{x}_{min}$ is that one of the
$\boldsymbol{x}^{k}_{j}$’s giving the minimal sum $D_{min}$ and it is stored
in $Comb(k)$. As well the sum $D_{min}$ is added to $S_{cur}$.

In the above procedure every triple of neighboring vectors from the results of
Algorithm 1 is considered to possibly lie on the sought-after smoothest
combination with the smallest sum of all squared second differences. The three
vectors are used as start points to find a smooth combination with a greedy
strategy, where only a vector is added to the current combination, if it gives
the smallest sum $D_{min}$ at the left or right end of the growing set with
vectors already added, until the first and last position are reached.

Finally from all of the combinations constructed this way the smoothest one
with the smallest sum $S_{min}$ out of all $S_{cur}$’s is selected in lines 45
\- 48 to be the final output $SmoothestCombination$.

Define $N_{total}:=\sum_{j=1}^{s}N_{j}$. Then the total number of operations
needed for Algorithm 2 can be estimated by
$\mathcal{O}\big{(}N_{total}\sum_{j=2}^{s-1}N_{j-1}N_{j}N_{j+1}\big{)}$ which
is considerably less than the $\mathcal{O}\big{(}\prod_{j=1}^{s}N_{j}\big{)}$
operations needed by the naive method of iterating through all possible
combinations.

##  9 Further Regularization of Coupled Solutions

Not only for determining the smoothest refractive index curve reconstructions
formed from the results of Algorithm 1 the coupled view on the solutions is
beneficial - it also leads to further improvement of the results by Twomey-
regularization. Let us investigate the coupled approach in a probability
theoretical setting. Here we reuse the notations introduced in Section 4, i.e.
we let $\boldsymbol{x}^{1}$, …, $\boldsymbol{x}^{s}$ denote a set of solution
for any of the five optical windows. Then the joint posterior probability
density of $\boldsymbol{x}^{1}$, …, $\boldsymbol{x}^{s}$ is given by

$p(\boldsymbol{x}^{1},...,\boldsymbol{x}^{s}|\boldsymbol{e}^{1},...,\boldsymbol{e}^{s})=\prod_{j=1}^{s}p(\boldsymbol{x}^{j}|\boldsymbol{e}^{j})\propto\exp\left(-\textstyle{\frac{1}{2}}\sum_{j=1}^{s}\|\boldsymbol{\Sigma}_{j}^{-\frac{1}{2}}(\boldsymbol{f}_{t}^{j}(\boldsymbol{x}^{j})-\boldsymbol{e}^{j})\|_{2}^{2}\right)\prod_{j=1}^{s}I_{\Omega}(\boldsymbol{x}^{j}),$
(9.1)

where $\boldsymbol{e}^{j}$ is the data vector for the $j$-th wavelength
$l_{j}$ having $N$ entries with $N$ being the size of the radius grid.
Moreover $\boldsymbol{\Sigma}_{j}$ is the scaled covariance matrix for $l_{j}$
and $\boldsymbol{f}_{t}^{j}(\boldsymbol{x})$ is the applied model depending on
$l_{j}$. Note that we initially have differing truncation indices $t_{1}$, …,
$t_{s}$. Since the coefficient functions of the truncated model function
$\boldsymbol{f}_{t_{j}}(\boldsymbol{x})$ are decaying fast for each $t_{j}$,
it is convenient to change to the same truncation index
$t:=\max\\{t_{1},...,t_{s}\\}$ for all wavelengths $l_{1}$, …, $l_{s}$. The
errors introduced by doing so are negliglible. It is easy to show that
maximizing the joint density (9.1) is equivalent to maximize all single
densities $p(\boldsymbol{x}^{j}|\boldsymbol{e}^{j})$ independently, i.e. a
joint MAP-estimator

$\boldsymbol{x}^{1}_{opt},...,\boldsymbol{x}^{s}_{opt}=\underset{\boldsymbol{x}^{1},...,\boldsymbol{x}^{s}}{\mathrm{argmax}}\;p(\boldsymbol{x}^{1},...,\boldsymbol{x}^{s}|\boldsymbol{e}^{1},...,\boldsymbol{e}^{s})$

consists of the single MAP-estimators

$\boldsymbol{x}^{j}_{opt}=\underset{\boldsymbol{x}}{\mathrm{argmax}}\;p(\boldsymbol{x}|\boldsymbol{e}^{j})$

for $j=1,...,s$. This means the the results of Algorithm 1 can be used to
construct MAP-estimators for the joint posterior probabilty density.

This behavior changes when we replace the joint prior probabilty density

$p_{prior}(\boldsymbol{x}^{1},...,\boldsymbol{x}^{s})=\left(\mathrm{vol}(\Omega)\right)^{-s}\prod_{j=1}^{s}I_{\Omega}(\boldsymbol{x}^{j})$

with

$\displaystyle
p_{prior}(\boldsymbol{x}^{1},...,\boldsymbol{x}^{s})\propto\exp\left(-\textstyle{\frac{1}{2}}\gamma
S(\boldsymbol{x}^{1},...,\boldsymbol{x}^{s})\right)\prod_{j=1}^{s}I_{\Omega}(\boldsymbol{x}^{j}),$

where

$\displaystyle S(\boldsymbol{x}^{1},...,\boldsymbol{x}^{s}):=$
$\displaystyle\sum_{i=2}^{s-1}\left(\Big{(}\big{(}\boldsymbol{x}^{i-1}\big{)}_{1}-2\big{(}\boldsymbol{x}^{i}\big{)}_{1}+\big{(}\boldsymbol{x}^{i+1}\big{)}_{1}\Big{)}^{2}+\Big{(}\big{(}\boldsymbol{x}^{i-1}\big{)}_{2}-2\big{(}\boldsymbol{x}^{i}\big{)}_{2}+\big{(}\boldsymbol{x}^{i+1}\big{)}_{2}\Big{)}^{2}\right)$
$\displaystyle+\rho$
$\displaystyle\sum_{i=1}^{s}\Big{(}\big{(}\boldsymbol{x}^{i}\big{)}_{1}^{2}+\big{(}\boldsymbol{x}^{i}\big{)}_{2}^{2}\Big{)},$

where $\gamma$ is a regularization parameter and $\rho$ is a parameter
specifying the amount Tikhonov regularization.

In the new prior distribution we use a combination of Tikhonov and Phillips-
Twomey-regularization both in the real and imaginary parts. Here we apply a
small amount of Tikhonov-regularization by setting $\rho=10^{-8}$, such that
the resulting regularization operator gets regular. This means that the
regularized regression problem (9.2) can be transformed into standard Tikhonov
form and that the monotonicity results from [9] are still valid. Each second
finite difference is clearly a function of three neighboring points, therefore
a decoupled computation of the joint MAP-estimator

$\boldsymbol{x}^{1}_{opt},...,\boldsymbol{x}^{s}_{opt}:=\underset{\boldsymbol{x}^{1},...,\boldsymbol{x}^{s}}{\mathrm{argmax}}\;p(\boldsymbol{x}^{1},...,\boldsymbol{x}^{s}|\boldsymbol{e}^{1},...,\boldsymbol{e}^{s})$
(9.2)

with

$p(\boldsymbol{x}^{1},...,\boldsymbol{x}^{s}|\boldsymbol{e}^{1},...,\boldsymbol{e}^{s})\propto\exp\left(-\textstyle{\frac{1}{2}}\sum_{j=1}^{s}\|\boldsymbol{\Sigma}_{j}^{-\frac{1}{2}}(\boldsymbol{f}_{t}^{j}(\boldsymbol{x}^{j})-\boldsymbol{e}^{j})\|_{2}^{2}-\textstyle{\frac{1}{2}}\gamma
S(\boldsymbol{x}^{1},...,\boldsymbol{x}^{s})\right)\prod_{j=1}^{s}I_{\Omega}(\boldsymbol{x}^{j})$

for each wavelength seperately is not possible anymore after changing to the
new prior density. However the result vectors $\boldsymbol{x}^{1}$, …,
$\boldsymbol{x}^{s}$ from Algorithm 2 form a good start vector to solve the
nonlinear regression problem (9.2).

We selected the regularization parameter $\gamma$ using the discrepancy
principle, i.e. we compute $\gamma$ such that the regularized solution

$\displaystyle\boldsymbol{x}^{1}_{\gamma},...,\boldsymbol{x}^{s}_{\gamma}:=\underset{\boldsymbol{x}^{1},...,\boldsymbol{x}^{s}}{\mathrm{argmin}}\;\sum_{j=1}^{s}\|\boldsymbol{\Sigma}_{j}^{-\frac{1}{2}}(\boldsymbol{f}_{t}^{j}(\boldsymbol{x}^{j})-\boldsymbol{e}^{j})\|_{2}^{2}+\gamma
S(\boldsymbol{x}^{1},...,\boldsymbol{x}^{s})\quad\text{s.t.}\quad\boldsymbol{x}^{j}\in\Omega,\quad
j=1,...,s$

fulfills a relation of the form

$\displaystyle\sum_{j=1}^{s}\|\boldsymbol{\Sigma}_{j}^{-\frac{1}{2}}(\boldsymbol{f}_{t}^{j}(\boldsymbol{x}^{j}_{\gamma})-\boldsymbol{e}^{j})\|_{2}^{2}=R(\delta),$

where $R(\delta)$ is a proposed residual value depending on the noise level
$\delta$. In [2] several different residual values are proposed for a fixed
model discretization and a set of regularization parameters is obtained from
those using the discrepancy principle. The pairings of model discretizations
and regularization parameters obtained this way are compared by their Bayesian
posterior probabilities. For these probabilities a set of integrals of the
different model posterior densities is needed to be computed, which can be
done with Monte Carlo integration methods. Due to the highly nonlinear
behavior of our model $\boldsymbol{f}_{t}(\boldsymbol{x})$ such integration
methods are not available here. Therefore we simplified the posterior
exploration in such a way that only one residual value is proposed.

Since each observed probability density
$p(\boldsymbol{e}^{j}|\boldsymbol{x}^{j})$ for $j=1,...,s$ is Gaussian, the
joint observed density
$p(\boldsymbol{e}^{1},...,\boldsymbol{e}^{1}|\boldsymbol{x}^{1},...,\boldsymbol{x}^{s})=\prod_{j=1}^{s}p(\boldsymbol{e}^{j}|\boldsymbol{x}^{j})$
is Gaussian as well. We have $\boldsymbol{x}^{j}\in\mathbb{R}^{2}$, thus the
sum of residuals
$\sum_{j=1}^{s}\|\boldsymbol{\Sigma}_{j}^{-\frac{1}{2}}(\boldsymbol{f}_{t}^{j}(\boldsymbol{x}^{j})-\boldsymbol{e}^{j})\|_{2}^{2}$
running through all wavelengths in the optical window is
$\chi^{2}(2s)$-distributed. This yields

$\displaystyle\mathbb{E}\Big{(}\sum_{j=1}^{s}\|\boldsymbol{\Sigma}_{j}^{-\frac{1}{2}}(\boldsymbol{f}_{t}^{j}(\boldsymbol{x}^{j})-\boldsymbol{e}^{j})\|_{2}^{2}\Big{)}=2s.$

Now a widely proposed residual value for the discrepancy principle is
$\tau\cdot 2s$, where $\tau=1.1$ is the so-called Morozov safety factor. This
choice is prone to under- or overregularization since the residual value
corresponding to the “true” solution might be much smaller or bigger than
$2\tau s$. Therefore we proposed a residual value which depends more
dynamically on the observed behovior of the residual.

Let $\boldsymbol{x}^{1}_{0}$, …, $\boldsymbol{x}^{s}_{0}$ denote the
unregularized solutions, i.e. the results of Algorithm 2. Then their squared
residual is given by
$R_{0}:=\sum_{j=1}^{s}\|\boldsymbol{\Sigma}_{j}^{-\frac{1}{2}}(\boldsymbol{f}_{t}^{j}(\boldsymbol{x}^{j}_{0})-\boldsymbol{e}^{j})\|_{2}^{2}$.
We first proposed

$R(\delta)=\max\left\\{2\tau_{1}s,\;\tau R_{0}\right\\},$

where we selected $\tau_{1}=1.1$. This means that the residual of the
regularized solution is beginning at $R_{0}$ at least increased by the factor
$\tau_{1}$, which avoids underregularization. If it then happens that
$\frac{R(\delta)}{R_{0}}>\theta$ with $\theta=1.5$ the proposed residual is
most likely too big and overregularization occurs. In this case we corrected
$R(\delta)$ by setting

$R(\delta)=\max\left\\{2\tau_{2}s,\;\theta R_{0}\right\\}$

with $\tau_{2}=0.9$.

##  10 Numerical Results

To see how reliable our proposed reconstruction algortihm is, we performed for
each of the scatterer materials Ag, $\mathrm{H}_{2}$O a numerical study with
$100$ sweeps through all $48$ wavelenghts of the five optical windows with the
same settings as in Section 6. We found out that the radii $r_{1}:=0.1\;\mu$m,
$r_{2}:=0.2\;\mu$m and $r_{3}:=0.3\;\mu$m contain the most information about
the refractive indices. This was found by keeping our $48$ wavelengths fixed
and comparing the quality of inversion results under varying aerosol particle
radii. Bigger radii did not improve the results in our simulations and
refractive index reconstructions only using bigger radii even turned out to be
too unstable. A more thorough treatment of this problem can be found in [7],
where a covariance eigenvalue analysis is used. Although not directly
comparable with our study of uncoated particles, the coated radii
$0.0975\;\mu$m, $0.2305\;\mu$m and $0.11\;\mu$m carrying the most information
content found in this study are roughly comparable to our radii.

We computed original spectral extinctions

$\left(\boldsymbol{e}_{true}\right)_{i,j}:=\pi
r_{i}^{2}\sum_{n=1}^{N_{tr}}q_{n}(m_{med}(l_{i}),m_{part}(l_{i}),r_{j},l_{i}),\quad
i=1,...,48,\quad j=1,...,3$

for Ag, $\mathrm{H}_{2}$O and CsI and added zero-mean Gaussian noise to it in
order to obtain the simulated noisy spectral extinctions

$\left(\boldsymbol{e}\right)_{i,j}=\left(\boldsymbol{e}_{true}\right)_{i,j}+\delta_{i,j}\quad\text{
with
}\;\delta_{i,j}\sim\mathcal{N}(0,(0.05\cdot\left(\boldsymbol{e}_{true}\right)_{i,j})^{2}),\quad
i=1,...,48,\quad j=1,...,3.$

The standard deviations were taken to be $5\%$ of the true spectral
extinctions. Real experiments using $500$ wavelengths were contaminated by
Gaussian noise with $30\%$ of the true spectral extinctions as standard
deviations. We expect that switching to $48$ wavelengths and thus increasing
the time resolution of the measurements will lower the standard deviations to
a small percentage. We used a sample size of $N_{s}=300$ to compute each mean
$\left(\boldsymbol{e}_{real}\right)_{i,j}$ of noisy spectral extinctions.

In the following the results are presented separately for each of the three
materials. The uppermost plot presents the relative errors of the
unregularized solutions obtained with Algorithm 2 from the original scatterer
refractive indices. The next plot displays the run times of Algorithm 1, which
returned all local minima of (4.3). These candidate solutions served as input
for Algorithm 2. Then the relative errors of the regularized solutions are
presented. Finally the relative errors of the average of the regularized
solutions are shown.

###  10.1 Results for Ag

####  10.1.1 Results of Algorithms 1 and 2

####  10.1.2 Relative Errors of the Regularized Solutions

####  10.1.3 Relative Errors of the Average of the Regularized Solutions

###  10.2 Results for CsI

####  10.2.1 Results of Algorithms 1 and 2

####  10.2.2 Relative Errors of the Regularized Solutions

####  10.2.3 Relative Errors of the Average of the Regularized Solutions

###  10.3 Results for $\mathrm{H}_{2}$O

####  10.3.1 Results of Algorithms 1 and 2

####  10.3.2 Relative Errors of the Regularized Solutions

####  10.3.3 Relative Errors of the Average of the Regularized Solutions

###  10.4 Conclusion

The severest relative errors can be observed for Ag. For the initial
unregularized solutions they lie between $1$ and $5\%$ on average and can go
up to ca. $53\%$ in the extreme cases as one can see in the leftmost subplot
for the first optical window. The run times of Algorithm 1 lie between $30$
and $50$ seconds in the average case and can rise up to $200$ seconds in the
extreme cases. A typical sweep through all $48$ wavelengths needed ca. $30$
minutes in total and this value was very much the same for all three
materials. For Ag the regularization procedure effectively reduced the
relative errors such that they are in the range between $0.5$ and $2.2\%$ on
average and are below $10\%$ in the extreme cases. Finally one can see in the
last plot that the relative errors of the average of all $100$ regularized
solutions are all below $0.4\%$.

For CsI the relative errors of the unregularized solutions are already quite
small and lie between $0.03$ and $0.065\%$ on average and rise only up to
$0.3\%$ in the extreme cases. The run times of Algorithm 1 are typically in
the range from $25$ to $55$ seconds and are always below $95$ seconds. The
regularization of the solutions brought only a small improvement of the
results here such that the relative errors did not change much. They are still
in the same range from $0.03$ and $0.065\%$ on average but only reach up to
ca. $0.2\%$ now. The relative errors of the average of the $100$ regularized
solutions are between $0.01$ and $0.055\%$.

Also for $\mathrm{H}_{2}$O the relative errors of the unregularized solutions
are comparably small and are below $0.35\%$ on average and still below $1.3\%$
in the extreme cases. Especially the rightmost subplot for the last optical
window shows the biggest relative errors, whereas for all the other optical
windows the relative errors are below $0.03\%$ on average and below $0.15\%$
in the extreme cases. A similar behavior can be observed for the run times of
Algorithm 1. For the first four optical windows they are between $20$ and $45$
seconds on average and below $100$ seconds in the extreme cases, whereas for
the last optical window they are between $30$ and $170$ seconds on average and
can even rise up to $350$ seconds. For $\mathrm{H}_{2}$O the regularization
procedure improves the relative errors only slightly for the first four
optical windows and even increases them for the last optical window such that
they can rise up to ca. $0.06\%$ on average and $1.4\%$ in the extreme cases.
The relative errors of the average of the $100$ regularized solutions are
virtually zero for the first four optical windows and below $0.55\%$ for the
last optical window.

##  11 Higher Noise Levels

To see how our proposed reconstruction algortihm behaves for higher noise
levels, we performed for each of the scatterer materials Ag, $\mathrm{H}_{2}$O
and CsI two numerical studies with $10$ sweeps through all $48$ wavelenghts of
the five optical windows with the same settings as in Section 6. We computed
original spectral extinctions

$\left(\boldsymbol{e}_{true}\right)_{i,j}:=\pi
r_{i}^{2}\sum_{n=1}^{N_{tr}}q_{n}(m_{med}(l_{i}),m_{part}(l_{i}),r_{j},l_{i}),\quad
i=1,...,48,\quad j=1,...,3$

for all wavelengths and added zero-mean Gaussian noise to it in order to
obtain the simulated noisy spectral extinctions

$\left(\boldsymbol{e}\right)_{i,j}=\left(\boldsymbol{e}_{true}\right)_{i,j}+\delta_{i,j}\quad\text{
with
}\;\delta_{i,j}\sim\mathcal{N}(0,(0.15\cdot\left(\boldsymbol{e}_{true}\right)_{i,j})^{2}),\quad
i=1,...,48,\quad j=1,...,3.$

for the first study and

$\left(\boldsymbol{e}\right)_{i,j}=\left(\boldsymbol{e}_{true}\right)_{i,j}+\delta_{i,j}\quad\text{
with
}\;\delta_{i,j}\sim\mathcal{N}(0,(0.3\cdot\left(\boldsymbol{e}_{true}\right)_{i,j})^{2}),\quad
i=1,...,48,\quad j=1,...,3.$

for the second. The standard deviations were taken to be $15\%$ and $30\%$
respectively of the true spectral extinctions. We used a sample size of
$N_{s}=300$ to compute each means $\left(\boldsymbol{e}_{real}\right)_{i,j}$
of noisy spectral extinctions.

For brevity we only present the relative errors of the average of the $10$
regularized solutions.

###  11.1 Results for Ag

####  11.1.1 Standard Deviation of $15\%$

####  11.1.2 Standard Deviation of $30\%$

###  11.2 Results for CsI

####  11.2.1 Standard Deviation of $15\%$

####  11.2.2 Standard Deviation of $30\%$

###  11.3 Results for $\mathrm{H}_{2}$O

####  11.3.1 Standard Deviation of $15\%$

####  11.3.2 Standard Deviation of $30\%$

###  11.4 Conclusion

Whereas the relative errors for CsI and $\mathrm{H}_{2}$O are still below
$1\%$, they can rise up to ca. $53\%$ for Ag. Therefore the reconstructed
refractive indices for Ag under this noise level are most likely not of
practical use. This shows that the FASP measurements of monodisperse aerosols
must be sufficiently accurate in order to retrieve the scatterer refractive
indices from them.

##  12 Numerical Study

We performed four numerical studies for two-component aerosols with log-
normal, RRSB and Hedrih model size distributions as outlined in [2]. The
aerosol particles were assumed to be homogeneously internally mixed, such that
only one effective refractive index was retrieved. One component of the
simulated aerosols was $\mathrm{H}_{2}$O with volume fractions of $0$, $11$,
$22$, $33$, $44$, $56$, $67$, $78$, $89$ and $100\%$. In the first two studies
we simulated mixtures of $\mathrm{H}_{2}$O and CsI, where we used the original
aerosol component refractive indices for the first study. For the second study
we used the average of the $100$ regularized solutions from Section 10. We did
the same for the third and fourth study, but here we simulated mixtures of
$\mathrm{H}_{2}$O and Ag. In the third study we utilized the original aerosol
component refractive indices and for the fourth the average of the $100$
regularized solutions from Section 10.

We applied the same reconstruction methods described in [2] under the same
settings, i.e. for each reconstruction we generated $300$ artificial noisy
measurements for all $48$ wavelengths, where the measurement error was
simulated as additive zero-mean Gaussian noise. For each wavelength, the
standard deviations were taken as $5\%$ of the solutions of the forward
problem. In [2] three different regularization methods, namely Tikhonov,
minimal first differences and Twomey regularization, were compared and their
results turned out to be very similar. Therefore we only used Tikhonov
regularization in the following. The results for the first study were directly
adopted from [2].

For every inversion we computed the $L^{2}$-error of the obtained
reconstruction relative to the original size distribution and measured the
total run time needed for the inversion. The computations were performed on a
notebook with a $2.27$ GHz CPU and $3.87$ GB accessible primary memory.

##  13 Results for Mixtures of $\mathrm{H}_{2}$O and CsI

####  13.0.1 Noise-free Refractive Indices

original $\mathrm{H}_{2}$O |
---|---
average $L^{2}$-errors (%)
Log-Normal |  RRSB |  Hedrih
volume percent |  Distribution |  Distribution |  Distribution
0 % |  33.5018 |  38.4939 |  17.2492
11 % |  30.0691 |  31.7809 |  17.5805
22 % |  28.8329 |  30.4269 |  17.3665
33 % |  24.9249 |  27.7925 |  15.0919
44 % |  23.4809 |  24.2137 |  16.9073
56 % |  22.6134 |  21.5396 |  15.5153
67 % |  20.4546 |  19.5181 |  14.5155
78 % |  18.9590 |  16.9927 |  16.7670
89 % |  18.5494 |  14.4005 |  13.3015
100 % |  17.9452 |  12.0441 |  11.3083
original $\mathrm{H}_{2}$O |
---|---
average fraction deviation (%)
Log-Normal |  RRSB |  Hedrih
volume percent |  Distribution |  Distribution |  Distribution
0 % |  10.7150 |  6.5800 |  5.0800
11 % |  7.3700 |  5.3000 |  6.3050
22 % |  6.1750 |  4.7100 |  4.2600
33 % |  4.3700 |  3.7300 |  4.4350
44 % |  3.9550 |  3.7300 |  3.2100
56 % |  3.2000 |  3.1650 |  3.3750
67 % |  2.6050 |  2.3650 |  1.8700
78 % |  2.2350 |  1.8050 |  2.7750
89 % |  2.2000 |  1.3750 |  1.9150
100 % |  1.4150 |  0.4650 |  1.0050

####  13.0.2 Noisy Refractive Indices

original $\mathrm{H}_{2}$O |
---|---
average $L^{2}$-errors (%)
Log-Normal |  RRSB |  Hedrih
volume percent |  Distribution |  Distribution |  Distribution
0 % |  32.8763 |  42.3379 |  18.1280
11 % |  31.3243 |  31.8214 |  17.1997
22 % |  27.6456 |  32.6273 |  17.8791
33 % |  24.5193 |  27.1816 |  15.4878
44 % |  23.4746 |  24.2689 |  16.0616
56 % |  22.8027 |  21.3374 |  15.1618
67 % |  19.6712 |  17.9423 |  14.4443
78 % |  18.4366 |  16.5406 |  17.5635
89 % |  18.4948 |  14.1970 |  13.8089
100 % |  18.4893 |  11.2708 |  12.4506

original $\mathrm{H}_{2}$O |
---|---
average fraction deviation (%)
Log-Normal |  RRSB |  Hedrih
volume percent |  Distribution |  Distribution |  Distribution
0 % |  10.9350 |  7.6900 |  4.8150
11 % |  6.6400 |  4.8600 |  6.5900
22 % |  5.6600 |  4.8650 |  5.0150
33 % |  4.1400 |  3.9450 |  4.1950
44 % |  3.8600 |  3.4400 |  3.4650
56 % |  3.7550 |  2.8200 |  3.3550
67 % |  2.3500 |  2.3900 |  1.9550
78 % |  2.0250 |  1.9600 |  2.8300
89 % |  2.1350 |  1.2600 |  1.9050
100 % |  1.6100 |  0.3350 |  1.0200

###  13.1 Results for Mixtures of $\mathrm{H}_{2}$O and Ag

####  13.1.1 Noise-free Refractive Indices

original $\mathrm{H}_{2}$O |
---|---
average $L^{2}$-errors (%)
Log-Normal |  RRSB |  Hedrih
volume percent |  Distribution |  Distribution |  Distribution
0 % |  61.2258 |  67.4898 |  54.9872
11 % |  47.6091 |  59.5032 |  49.4157
22 % |  37.1819 |  67.8529 |  36.8106
33 % |  57.4058 |  69.0461 |  46.9500
44 % |  54.5266 |  71.4585 |  51.6323
56 % |  45.6020 |  67.5288 |  36.6557
67 % |  36.0602 |  53.0986 |  22.8959
78 % |  33.9722 |  43.2794 |  16.1026
89 % |  24.1114 |  27.5773 |  16.0464
100 % |  21.0693 |  11.1298 |  11.5123
original $\mathrm{H}_{2}$O |
---|---
average fraction deviation (%)
Log-Normal |  RRSB |  Hedrih
volume percent |  Distribution |  Distribution |  Distribution
0 % |  0 |  0 |  0
11 % |  0.2650 |  0.3500 |  0.2950
22 % |  0.5350 |  3.6650 |  1.8050
33 % |  18.5000 |  16.4750 |  10.4350
44 % |  13.1200 |  13.5400 |  14.5050
56 % |  8.6300 |  9.6450 |  8.0700
67 % |  6.5950 |  6.8300 |  5.7750
78 % |  3.4650 |  3.2900 |  1.9750
89 % |  1.3150 |  1.5450 |  1.0600
100 % |  0.4750 |  0.0600 |  0.1600

####  13.1.2 Noisy Refractive Indices

original $\mathrm{H}_{2}$O |
---|---
average $L^{2}$-errors (%)
Log-Normal |  RRSB |  Hedrih
volume percent |  Distribution |  Distribution |  Distribution
0 % |  62.0576 |  67.5179 |  57.9628
11 % |  46.9321 |  59.9884 |  51.1218
22 % |  38.6799 |  65.8397 |  36.0770
33 % |  56.2737 |  69.8686 |  47.1554
44 % |  55.2457 |  70.4010 |  53.9243
56 % |  44.6657 |  67.1041 |  37.0255
67 % |  37.8969 |  55.6161 |  23.9009
78 % |  39.5606 |  45.8995 |  16.6339
89 % |  24.1072 |  28.0702 |  16.2226
100 % |  18.3892 |  9.8479 |  11.4723
original $\mathrm{H}_{2}$O |
---|---
average fraction deviation (%)
Log-Normal |  RRSB |  Hedrih
volume percent |  Distribution |  Distribution |  Distribution
0 % |  0 |  0 |  0
11 % |  0.2700 |  0.3950 |  0.3200
22 % |  0.4450 |  2.9900 |  1.0750
33 % |  16.0650 |  16.4650 |  10.6400
44 % |  13.1950 |  14.2150 |  15.4550
56 % |  8.2350 |  9.0350 |  8.8350
67 % |  6.4650 |  7.1100 |  5.5450
78 % |  3.8000 |  3.1700 |  1.8200
89 % |  1.3250 |  1.5600 |  1.0900
100 % |  0.3000 |  0.0200 |  0.0800

###  13.2 Conclusion

The resuts of the first and second study only differ by ca. $4\%$ at most and
behave very similarly. The same is for the third and fourth study. These
numerical results indicate that $100$ FASP measurement sweeps consisting of
$300$ single measurements with an accuracy as in Section 10 are sufficient to
determine aerosol refractive indices in such a quality, that they are suitable
for particle size distribution reconstructions for two-component homogeneously
internally mixed aerosols using the FASP. The particle radii of the three
monodisperse aerosols generated for the refractive indices retrieval need to
be $0.1\;\mu$m, $0.2\;\mu$m and $0.3\;\mu$m respectively.

##  14 Outlook

It is of interest to investigate if the methods derived in this study can be
extended to the case of core-plus-shell aerosols.

##  15 Acknowledgement

I thank Graham Alldredge, Ph.D. for proofreading the manuscript, useful
recommendations and fruitful discussions on this topic.

This work is sponsored by the German Federal Ministry of Education and
Research (BMBF) under the contract number 02NUK022A.

Responsibility for the content of this report lies with the authors.

## References

  * [1] A. A. Riziq, C. Erlick, E. Dinar, and Y. Rudich, “Optical properties of absorbing and non-absorbing aerosols retrieved by cavity ring down (crd) spectroscopy,” Atmos. Chem. Phys., Vol. 7, 2007.
  * [2] G. Alldredge and T. Kyrion, “Robust inversion methods for aerosol spectroscopy,” Inverse Problems in Science and Engineering, 2016.
  * [3] Q. Fu and W. Sun, “Mie theory for light scattering by a spherical particle in an absorbing medium,” Applied Optics, Vol. 40, Issue 9, 2001.
  * [4] W. Wiscombe, “Improved mie scattering algorithms,” Applied Optics, Vol. 19, Issue 9, 1980.
  * [5] M. Abromovitz and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.  Dover Publications, 1972.
  * [6] P. Deuflhard and A. Hohmann, Numerical Mathematics I: an algorithmically oriented introduction.  de Gruyter, 2002.
  * [7] C. Erlick, M. Haspel, and Y. Rudich, “Simultaneous retrieval of the complex refractive indices of the core and shell of coated aerosol particles from extinction measurements using simulated annealing,” Appl. Opt., 2011.
  * [8] R. M. Goody and Y. L. Yung, Atmospheric Radiation. Theoretical Basis, second edition.  Oxford University Press, New York, 1989.
  * [9] H. W. Engl, K. Kunisch, and A. Neubauer, “Convergence rates for Tikhonov regularisation of non-linear ill-posed problems,” Inverse Problems, vol. 5, pp. 523–540, 1989.
  * [10] O. Scherzer, H. W. Engl, and K.Kunisch, “Optimal a Posteriori Parameter Choice for Tikhonov Regularization for Solving Nonlinear Ill-Posed Problems,” SIAM J. Numer. Anal., vol. 6, pp. 1796–1838, 1993.
  * [11] M. Quinten, Optical Properties of Nanoparticle Systems.  WILEY-VCH Verlag, 2010.