Reconstruction of refractive indices from spectral measurements of monodisperse aerosols

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 temperature 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.

Keywords:

• Nonlinear inverse problem
• refractive index reconstruction
• Mie theory
• nonlinear Tikhonov regularization
• nonlinear regression

AMS Subject Classifications:

• 15A29
• 34K29
• 34M50
• 45Q05
• 70F17

1. Set-up 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_{\mathrm{long}}(l)$ and $I_{\mathrm{short}}(l)$ having passed a long and a short measurement path length $L_{\mathrm{long}}$ and $L_{\mathrm{short}}$, respectively. The evaluations of the FASP measurements are based on the relation(1.1) $$\begin{array}{c}\hfill {\int }_0^{\infty }k(r,l)n(r)\mathrm{d}r=e(l)\text{with}e(l)=-\frac{log(I_{\mathrm{long}}(l))-log(I_{\mathrm{short}}(l))}{L_{\mathrm{long}}-L_{\mathrm{short}}},\end{array}$$(1.1)

where $k(r,l):=\mathit{\pi }{r}^2{Q}_{ext}(m_{\mathrm{med}}(l),m_{\mathrm{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_{\mathrm{med}}(l)$ and $m_{\mathrm{part}}(l)$ are the refractive indices of the surrounding medium and the particle material depending on the wavelength l. The function $Q_{\mathrm{ext}}(m_{\mathrm{med}}(l),m_{\mathrm{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(1.1) $$\begin{array}{c}\hfill {\int }_0^{\infty }k(r,l)n(r)\mathrm{d}r=e(l)\text{with}e(l)=-\frac{log(I_{\mathrm{long}}(l))-log(I_{\mathrm{short}}(l))}{L_{\mathrm{long}}-L_{\mathrm{short}}},\end{array}$$(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\mathit{\delta }(r-r_m)$, where n is the total number of particles and $\mathit{\delta }(r-r_m)$ is a Dirac delta distribution truncated on the positive half-axis. Inserting this into (1.1(1.1) $$\begin{array}{c}\hfill {\int }_0^{\infty }k(r,l)n(r)\mathrm{d}r=e(l)\text{with}e(l)=-\frac{log(I_{\mathrm{long}}(l))-log(I_{\mathrm{short}}(l))}{L_{\mathrm{long}}-L_{\mathrm{short}}},\end{array}$$(1.1) ) gives(1.2) $$\begin{array}{c}\hfill n\mathit{\pi }{r}_m^2{Q}_{\mathrm{ext}}(m_{\mathrm{med}}(l),m_{\mathrm{part}}(l),r_m,l)=e(l),\end{array}$$(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.$$\begin{array}{c}\hfill Q_{\mathrm{ext}}(m_{\mathrm{med}}(l),m_{\mathrm{part}}(l),r,l)=\sum _{n=1}^{\infty }{q}_n(m_{\mathrm{med}}(l),m_{\mathrm{part}}(l),r,l).\end{array}$$

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

It is clear that in practical computations $Q_{\mathrm{ext}}(m_{\mathrm{med}}(l),m_{\mathrm{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_{\mathrm{med}}(l),m_{\mathrm{part}}(l),r,l)$’s is practically feasible.

We now fix a wavelength l. The complex refractive index $m_{\mathrm{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_{\mathrm{true}}(r_i,l)+{\mathit{\delta }}_i$ with ${\mathit{\delta }}_i\sim \mathcal{N}(0,s_i^2)$ for $i=1,\dots ,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$$\begin{array}{c}\hfill q_{\mathrm{true}}(r_i,l)=n_i\mathit{\pi }{r}_i^2\sum _{n=1}^{\infty }{q}_n(m_{\mathrm{med}}(l),m_{\mathrm{part}}(l),r_i,l),\end{array}$$

where $n_i$ is the number of particles having the same radius $r_i$. Then, a reconstrution of $m_{\mathrm{part}}(l)$ is obtained from the set of solutions M(l) of the nonlinear regression problem(1.3) $$\begin{array}{c}\hfill M(l):=\underset{m\in \mathbb{C}}{\mathrm{a}rgmin}\sum _{i=1}^N\frac{1}{2{\left(\frac{s_i}{n_i}\right)}^2}{\left(\mathit{\pi }{r}_i^2\sum _{n=1}^{N_{tr}}{q}_n(m_{\mathrm{med}}(l),m,r_i,l)-\frac{q(r_i,l)}{n_i}\right)}^2.\end{array}$$(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(1.3) $$\begin{array}{c}\hfill M(l):=\underset{m\in \mathbb{C}}{\mathrm{a}rgmin}\sum _{i=1}^N\frac{1}{2{\left(\frac{s_i}{n_i}\right)}^2}{\left(\mathit{\pi }{r}_i^2\sum _{n=1}^{N_{tr}}{q}_n(m_{\mathrm{med}}(l),m,r_i,l)-\frac{q(r_i,l)}{n_i}\right)}^2.\end{array}$$(1.3) ) in Section 4.

For solving (1.3(1.3) $$\begin{array}{c}\hfill M(l):=\underset{m\in \mathbb{C}}{\mathrm{a}rgmin}\sum _{i=1}^N\frac{1}{2{\left(\frac{s_i}{n_i}\right)}^2}{\left(\mathit{\pi }{r}_i^2\sum _{n=1}^{N_{tr}}{q}_n(m_{\mathrm{med}}(l),m,r_i,l)-\frac{q(r_i,l)}{n_i}\right)}^2.\end{array}$$(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(1.3) $$\begin{array}{c}\hfill M(l):=\underset{m\in \mathbb{C}}{\mathrm{a}rgmin}\sum _{i=1}^N\frac{1}{2{\left(\frac{s_i}{n_i}\right)}^2}{\left(\mathit{\pi }{r}_i^2\sum _{n=1}^{N_{tr}}{q}_n(m_{\mathrm{med}}(l),m,r_i,l)-\frac{q(r_i,l)}{n_i}\right)}^2.\end{array}$$(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 ${\mathit{\xi }}_n:\mathbb{C}\to \mathbb{C}$ and ${\mathit{\psi }}_n:\mathbb{C}\to \mathbb{C}$ given by(2.1) $$\begin{array}{c}\hfill {\mathit{\xi }}_n(z)=\sqrt{\frac{\mathit{\pi }}{2}}\sqrt{z}{J}_{n+\frac{1}{2}}(z)\end{array}$$(2.1)

and(2.2) $$\begin{array}{c}\hfill {\mathit{\psi }}_n(z)=\sqrt{\frac{\mathit{\pi }}{2}}\sqrt{z}{J}_{n+\frac{1}{2}}(z)+i\sqrt{\frac{\mathit{\pi }}{2}}\sqrt{z}{Y}_{n+\frac{1}{2}}(z),\end{array}$$(2.2)

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

We introduce the so-called Mie coefficients:(2.3) $$\begin{array}{cc}\hfill a_n& :=\frac{m_{\mathrm{part}}\dot{{\mathit{\xi }}_n}(z_{\mathrm{med}}){\mathit{\xi }}_n(z_{\mathrm{part}})-m_{\mathrm{med}}{\mathit{\xi }}_n(z_{\mathrm{med}})\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})}{m_{\mathrm{part}}\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}}){\mathit{\xi }}_n(z_{\mathrm{part}})-m_{\mathrm{med}}{\mathit{\psi }}_n(z_{\mathrm{med}})\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})}\hfill \\ \hfill b_n& :=\frac{m_{\mathrm{part}}{\mathit{\xi }}_n(z_{\mathrm{med}})\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})-m_{\mathrm{med}}\dot{{\mathit{\xi }}_n}(z_{\mathrm{med}}){\mathit{\xi }}_n(z_{\mathrm{part}})}{m_{\mathrm{part}}{\mathit{\psi }}_n(z_{\mathrm{med}})\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})-m_{\mathrm{med}}\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}}){\mathit{\xi }}_n(z_{\mathrm{part}})}\hfill \\ \hfill c_n& :=\frac{m_{\mathrm{part}}{\mathit{\psi }}_n(z_{\mathrm{med}})\dot{{\mathit{\xi }}_n}(z_{\mathrm{med}})-m_{\mathrm{part}}\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}}){\mathit{\xi }}_n(z_{\mathrm{med}})}{m_{\mathrm{part}}{\mathit{\psi }}_n(z_{\mathrm{med}})\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})-m_{\mathrm{med}}\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}}){\mathit{\xi }}_n(z_{\mathrm{part}})}\hfill \\ \hfill d_n& :=\frac{m_{\mathrm{part}}\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}}){\mathit{\xi }}_n(z_{\mathrm{med}})-m_{\mathrm{part}}{\mathit{\psi }}_n(z_{\mathrm{med}})\dot{{\mathit{\xi }}_n}(z_{\mathrm{med}})}{m_{\mathrm{part}}\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}}){\mathit{\xi }}_n(z_{\mathrm{part}})-m_{\mathrm{med}}{\mathit{\psi }}_n(z_{\mathrm{med}})\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})}\hfill \end{array}$$(2.3)

With the Mie coefficients, we can express the coefficient functions(2.4) $$\begin{array}{cc}\hfill A_n(\mathit{\rho },m_{\mathrm{med}},m_{\mathrm{part}})& :=\frac{l}{2\mathit{\pi }{m}_{\mathrm{part}}}\left({\left|c_n\right|}^2{\mathit{\xi }}_n(z_{\mathrm{part}})\overline{\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})}-{\left|d_n\right|}^2\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})\overline{{\mathit{\xi }}_n(z_{\mathrm{part}})}\right)\hfill \end{array}$$(2.4)

and(2.5) $$\begin{array}{cc}\hfill B_n(\mathit{\rho },m_{\mathrm{med}},m_{\mathrm{part}})& :=\frac{l}{2\mathit{\pi }{m}_{\mathrm{med}}}\left({\left|a_n\right|}^2\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}})\overline{{\mathit{\psi }}_n(z_{\mathrm{med}})}-{\left|b_n\right|}^2{\mathit{\psi }}_n(z_{\mathrm{med}})\overline{\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}})}\right),\hfill \end{array}$$(2.5)

which finally occur in the series expansion(2.6) $$\begin{array}{cc}\hfill Q_{\mathrm{ext}}(r,l,m_{\mathrm{med}},m_{\mathrm{part}})& =\frac{l}{2cI(r,l)}\sum _{n=1}^{\infty }(2n+1)\mathrm{I}m(A_n(\mathit{\rho },m_{\mathrm{med}},m_{\mathrm{part}})\hfill \\ \hfill & +B_n(\mathit{\rho },m_{\mathrm{med}},m_{\mathrm{part}})).\hfill \end{array}$$(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(2.7) $$\begin{array}{ccc}\hfill I(r,l)& {\displaystyle =\frac{l^2}{8\mathit{\pi }{(k_{\mathrm{med}})}^2}\frac{n_{\mathrm{med}}}{2c}(1+\left(4\mathit{\pi }{k}_{\mathrm{med}}\frac{r}{l}-1\right)e^{{\displaystyle 4\mathit{\pi }{k}_{\mathrm{med}}\frac{r}{l}}}),}\hfill & \hfill \mathrm{f}or{k}_{\mathrm{med}}\ne 0\\ \hfill I(r,l)& {\displaystyle =\mathit{\pi }{r}^2\frac{n_{\mathrm{med}}}{2c},}\hfill & \hfill \mathrm{f}or{k}_{\mathrm{med}}=0.\end{array}$$(2.7)

Obviously we cannot evaluate (2.6(2.6) $$\begin{array}{cc}\hfill Q_{\mathrm{ext}}(r,l,m_{\mathrm{med}},m_{\mathrm{part}})& =\frac{l}{2cI(r,l)}\sum _{n=1}^{\infty }(2n+1)\mathrm{I}m(A_n(\mathit{\rho },m_{\mathrm{med}},m_{\mathrm{part}})\hfill \\ \hfill & +B_n(\mathit{\rho },m_{\mathrm{med}},m_{\mathrm{part}})).\hfill \end{array}$$(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_{\mathrm{trunc}}$ is presented, which is given by(2.8) $$\begin{array}{cc}\hfill N_{\mathrm{trunc}}& =\lceil \left|M\right.+4.05·M^{\frac{1}{3}}+\left.2\right|\rceil ,\hfill \\ \hfill \mathrm{w}ithM& =max\lceil \left|\mathit{\rho }\right|,|\mathit{\rho }·m_{\mathrm{med}}|,|\mathit{\rho }·m_{\mathrm{part}}|\rceil .\hfill \end{array}$$(2.8)

3. Derivatives of the truncated Mie efficiency series

The Bessel functions $J_{\mathit{\alpha }}(z)$ and $Y_{\mathit{\alpha }}(z)$ for an arbitrary weight $\mathit{\alpha }$ fulfil the recurrence relations(3.1) $$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}z}(z^{\mathit{\alpha }}{J}_{\mathit{\alpha }}(z))=z^{\mathit{\alpha }}{J}_{\mathit{\alpha }-1}(z)\text{and}\frac{\mathrm{d}}{\mathrm{d}z}(z^{\mathit{\alpha }}{Y}_{\mathit{\alpha }}(z))=z^{\mathit{\alpha }}{Y}_{\mathit{\alpha }-1}(z),\end{array}$$(3.1)

cf. [5]. For the Bessel functions occurring in the Riccati–Bessel functions, ${\mathit{\xi }}_n(z)$ and ${\mathit{\psi }}_n(z)$ follow from this with the weight $\mathit{\alpha }=n+\frac{1}{2}$ that(3.2) $$\begin{array}{cc}\hfill \dot{{\mathit{\xi }}_n}(z)& =\sqrt{\frac{\mathit{\pi }}{2}}\sqrt{z}\left(J_{n-\frac{1}{2}}(z)-\frac{n}{z}{J}_{n+\frac{1}{2}}(z)\right)\hfill \end{array}$$(3.2) (3.3) $$\begin{array}{cc}\hfill \text{and}\dot{{\mathit{\psi }}_n}(z)& =\sqrt{\frac{\mathit{\pi }}{2}}\sqrt{z}\left(J_{n-\frac{1}{2}}(z)-\frac{n}{z}{J}_{n+\frac{1}{2}}(z)\right)+\sqrt{\frac{\mathit{\pi }}{2}}\sqrt{z}\left(Y_{n-\frac{1}{2}}(z)-\frac{n}{z}{Y}_{n+\frac{1}{2}}(z)\right)i\hfill \end{array}$$(3.3)

for $z\ne 0$.

We apply (3.1(3.1) $$\begin{array}{c}\hfill \frac{\mathrm{d}}{\mathrm{d}z}(z^{\mathit{\alpha }}{J}_{\mathit{\alpha }}(z))=z^{\mathit{\alpha }}{J}_{\mathit{\alpha }-1}(z)\text{and}\frac{\mathrm{d}}{\mathrm{d}z}(z^{\mathit{\alpha }}{Y}_{\mathit{\alpha }}(z))=z^{\mathit{\alpha }}{Y}_{\mathit{\alpha }-1}(z),\end{array}$$(3.1) ) a second time to get $\ddot{{\mathit{\xi }}_n}(z)$, which yields$$\begin{array}{c}\hfill \ddot{{\mathit{\xi }}_n}(z)=\sqrt{\frac{\mathit{\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).\end{array}$$

For an arbitrary weight $\mathit{\alpha }$, we have the recurrence relation$$\begin{array}{c}\hfill J_{\mathit{\alpha }-1}(z)=\frac{2\mathit{\alpha }}{z}{J}_{\mathit{\alpha }}(z)-J_{\mathit{\alpha }+1}(z),\end{array}$$

see [5], and we use it to eliminate the term $J_{n-\frac{3}{2}}(z)$ in the expression for $\ddot{{\mathit{\xi }}_n}(z)$. Then also $J_{n-\frac{1}{2}}(z)$ cancels out, such that we obtain the representation(3.4) $$\begin{array}{c}\hfill \ddot{{\mathit{\xi }}_n}(z)=\sqrt{\frac{\mathit{\pi }}{2}}\frac{\sqrt{z}}{z^2}\left(n(n+1)-z^2\right)J_{n+\frac{1}{2}}(z)\end{array}$$(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}\to \mathbb{C}$ with $f(z)=f(x+iy)=u(x,y)+iv(x,y)$ holds(3.5) $$\begin{array}{c}\hfill \dot{f}(z)=\frac{\mathrm{d}}{\mathrm{d}z}f(z)=\frac{\mathit{\partial }}{\mathit{\partial }x}f(x+iy)=-i\frac{\mathit{\partial }}{\mathit{\partial }y}f(x+iy).\end{array}$$(3.5)

From this follows(3.6) $$\begin{array}{c}\hfill u_x=\mathrm{R}e(\dot{f}(z)),u_y=-\mathrm{I}m(\dot{f}(z)),v_x=\mathrm{I}m(\dot{f}(z))\text{and}{v}_y=\mathrm{R}e(\dot{f}(z)).\end{array}$$(3.6)

Using the latter relations, we can deduce using $z_{\mathrm{part}}=\mathit{\rho }\left(n_{\mathrm{part}}+i{k}_{\mathrm{part}}\right)$$$\begin{array}{cccc}\hfill \frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\mathrm{R}e({\mathit{\xi }}_n(z_{\mathrm{part}}))& =\mathit{\rho }\mathrm{R}e(\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})),\hfill & \hfill \frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}\mathrm{R}e({\mathit{\xi }}_n(z_{\mathrm{part}}))& =-\mathit{\rho }\mathrm{I}m(\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})),\hfill \\ \hfill \frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\mathrm{I}m({\mathit{\xi }}_n(z_{\mathrm{part}}))& =\mathit{\rho }\mathrm{I}m(\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})),\hfill & \hfill \frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}\mathrm{I}m({\mathit{\xi }}_n(z_{\mathrm{part}}))& =\mathit{\rho }\mathrm{R}e(\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}}))\hfill \end{array}$$

and analogously$$\begin{array}{cccc}\hfill \frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\mathrm{R}e(\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}}))& =\mathit{\rho }\mathrm{R}e(\ddot{{\mathit{\xi }}_n}(z_{\mathrm{part}})),\hfill & \hfill \frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}\mathrm{R}e(\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}}))& =-\mathit{\rho }\mathrm{I}m(\ddot{{\mathit{\xi }}_n}(z_{\mathrm{part}})),\hfill \\ \hfill \frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\mathrm{I}m({\dot{\mathit{\xi }}}_n(z_{\mathrm{part}}))& =\mathit{\rho }\mathrm{I}m(\ddot{{\mathit{\xi }}_n}(z_{\mathrm{part}})),\hfill & \hfill \frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}\mathrm{I}m({\dot{\mathit{\xi }}}_n(z_{\mathrm{part}}))& =\mathit{\rho }\mathrm{R}e(\ddot{{\mathit{\xi }}_n}(z_{\mathrm{part}})).\hfill \end{array}$$

The Bessel function values$$\begin{array}{cc}& J_{0+\frac{1}{2}}(z_{\mathrm{med}}),\dots ,J_{N_{\mathrm{trunc}}+\frac{1}{2}}(z_{\mathrm{med}}),Y_{0+\frac{1}{2}}(z_{\mathrm{med}}),\dots ,Y_{N_{\mathrm{trunc}}+\frac{1}{2}}(z_{\mathrm{med}})\hfill \\ \hfill \text{and}& J_{0+\frac{1}{2}}(z_{\mathrm{part}}),\dots ,J_{N_{\mathrm{trunc}}+\frac{1}{2}}(z_{\mathrm{part}}).\hfill \end{array}$$

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_{\mathrm{part}}$ and $k_{\mathrm{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$$\begin{array}{cc}\hfill \frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}{\left|a_n\right|}^2& =\left(\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}{a}_n\right)\overline{a_n}+a_n\left(\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\overline{a_n}\right)\hfill \\ \hfill \text{and}\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}{\left|a_n\right|}^2& =\left(\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}{a}_n\right)\overline{a_n}+a_n\left(\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}\overline{a_n}\right).\hfill \end{array}$$

We write$$\begin{array}{cc}\hfill a_n=\frac{E}{D}\text{with}E& :=m_{\mathrm{part}}\dot{{\mathit{\xi }}_n}(z_{\mathrm{med}}){\mathit{\xi }}_n(z_{\mathrm{part}})-m_{\mathrm{med}}{\mathit{\xi }}_n(z_{\mathrm{med}})\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})\hfill \\ \hfill \text{and}D& :=m_{\mathrm{part}}\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}}){\mathit{\xi }}_n(z_{\mathrm{part}})-m_{\mathrm{med}}{\mathit{\psi }}_n(z_{\mathrm{med}})\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}}),\hfill \end{array}$$

which yields$$\begin{array}{cc}\hfill \frac{\mathrm{d}}{\mathrm{d}{m}_{\mathrm{part}}}{a}_n& =\frac{1}{D^2}\left(\left(\frac{\mathrm{d}}{\mathrm{d}{m}_{\mathrm{part}}}E\right)D-E\left(\frac{\mathrm{d}}{\mathrm{d}{m}_{\mathrm{part}}}D\right)\right)\hfill \\ \hfill \text{with}\frac{\mathrm{d}}{\mathrm{d}{m}_{\mathrm{part}}}E& =\dot{{\mathit{\xi }}_n}(z_{\mathrm{med}}){\mathit{\xi }}_n(z_{\mathrm{part}})+m_{\mathrm{part}}\dot{{\mathit{\xi }}_n}(z_{\mathrm{med}})\mathit{\rho }\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})\hfill \\ \hfill & -m_{\mathrm{med}}{\mathit{\xi }}_n(z_{\mathrm{med}})\mathit{\rho }\ddot{{\mathit{\xi }}_n}(z_{\mathrm{part}})\hfill \\ \hfill \text{and}\frac{\mathrm{d}}{\mathrm{d}{m}_{\mathrm{part}}}D& =\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}}){\mathit{\xi }}_n(z_{\mathrm{part}})+m_{\mathrm{part}}\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}})\mathit{\rho }\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})\hfill \\ \hfill & -m_{\mathrm{med}}{\mathit{\psi }}_n(z_{\mathrm{med}})\mathit{\rho }\ddot{{\mathit{\xi }}_n}(z_{\mathrm{part}}).\hfill \end{array}$$

Furthermore follow from (3.5(3.5) $$\begin{array}{c}\hfill \dot{f}(z)=\frac{\mathrm{d}}{\mathrm{d}z}f(z)=\frac{\mathit{\partial }}{\mathit{\partial }x}f(x+iy)=-i\frac{\mathit{\partial }}{\mathit{\partial }y}f(x+iy).\end{array}$$(3.5) ), the relations$$\begin{array}{cc}\hfill \frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}{a}_n& =\frac{\mathrm{d}}{\mathrm{d}{m}_{\mathrm{part}}}{a}_n\hfill \\ \hfill \text{and}\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}{a}_n& =\left(\frac{\mathrm{d}}{\mathrm{d}{m}_{\mathrm{part}}}{a}_n\right)i.\hfill \end{array}$$

Although $\overline{a_n}$ is not holomorphic with respect to $m_{\mathrm{part}}$, we can still compute the partial derivatives ${\displaystyle \frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\overline{a_n}}$ and ${\displaystyle \frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}\overline{a_n}}$. We obtain using (3.6(3.6) $$\begin{array}{c}\hfill u_x=\mathrm{R}e(\dot{f}(z)),u_y=-\mathrm{I}m(\dot{f}(z)),v_x=\mathrm{I}m(\dot{f}(z))\text{and}{v}_y=\mathrm{R}e(\dot{f}(z)).\end{array}$$(3.6) ) the relations$$\begin{array}{cc}\hfill \frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\overline{a_n}& =\overline{\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}{a}_n}\hfill \\ \hfill \text{and}\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}\overline{a_n}& =\overline{\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}{a}_n}.\hfill \end{array}$$

This completes the computations of ${\displaystyle \frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}{\left|a_n\right|}^2}$ and ${\displaystyle \frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}{\left|a_n\right|}^2}$.

For a holomorphic function $f(x+iy)$, we can easily deduce from (3.6(3.6) $$\begin{array}{c}\hfill u_x=\mathrm{R}e(\dot{f}(z)),u_y=-\mathrm{I}m(\dot{f}(z)),v_x=\mathrm{I}m(\dot{f}(z))\text{and}{v}_y=\mathrm{R}e(\dot{f}(z)).\end{array}$$(3.6) ) that$$\begin{array}{c}\hfill \frac{\mathit{\partial }}{\mathit{\partial }x}\mathrm{I}m\left(f(x+iy)\right)=\mathrm{I}m\left(\frac{\mathit{\partial }}{\mathit{\partial }x}f(x+iy)\right)\text{and}\frac{\mathit{\partial }}{\mathit{\partial }y}\mathrm{I}m\left(f(x+iy)\right)=\mathrm{I}m\left(\frac{\mathit{\partial }}{\mathit{\partial }y}f(x+iy)\right).\end{array}$$

This gives with respect to (3.5(3.5) $$\begin{array}{c}\hfill \dot{f}(z)=\frac{\mathrm{d}}{\mathrm{d}z}f(z)=\frac{\mathit{\partial }}{\mathit{\partial }x}f(x+iy)=-i\frac{\mathit{\partial }}{\mathit{\partial }y}f(x+iy).\end{array}$$(3.5) )$$\begin{array}{cc}\hfill \frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\mathrm{I}m\left(A_n\right)& =\mathrm{I}m\left(\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}{A}_n\right)\hfill \\ \hfill \text{with}\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}{A}_n& =\frac{l}{2\mathit{\pi }}\left(\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\left({\left|c_n\right|}^2\right)U_1+{\left|c_n\right|}^2\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}{U}_1\right.\hfill \\ \hfill & \left.-\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\left({\left|d_n\right|}^2\right)U_2-{\left|d_n\right|}^2\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}{U}_2\right),\hfill \\ \hfill \text{where}{U}_1& =\frac{{\mathit{\xi }}_n(z_{\mathrm{part}})\overline{\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})}}{m_{\mathrm{part}}},\hfill \\ \hfill \frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}{U}_1& =\frac{1}{m_{\mathrm{part}}^2}\left(\left(\mathit{\rho }\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})\overline{\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})}+{\mathit{\xi }}_n(z_{\mathrm{part}})\mathit{\rho }\overline{\ddot{{\mathit{\xi }}_n}(z_{\mathrm{part}})}\right)m_{\mathrm{part}}\right.\hfill \\ \hfill & \left.-{\mathit{\xi }}_n(z_{\mathrm{part}})\overline{\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})}\right),\hfill \\ \hfill \text{and}{U}_2& =\frac{{\dot{\mathit{\xi }}}_n(z_{\mathrm{part}})\overline{{\mathit{\xi }}_n(z_{\mathrm{part}})}}{m_{\mathrm{part}}},\hfill \\ \hfill \frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}{U}_2& =\frac{1}{m_{\mathrm{part}}^2}\left(\left(\mathit{\rho }\ddot{{\mathit{\xi }}_n}(z_{\mathrm{part}})\overline{{\mathit{\xi }}_n(z_{\mathrm{part}})}+{\dot{\mathit{\xi }}}_n(z_{\mathrm{part}})\mathit{\rho }\overline{\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})}\right)m_{\mathrm{part}}\right.\hfill \\ \hfill & \left.-{\dot{\mathit{\xi }}}_n(z_{\mathrm{part}})\overline{{\mathit{\xi }}_n(z_{\mathrm{part}})}\right).\hfill \end{array}$$

Analogously we obtain$$\begin{array}{cc}\hfill \frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}\mathrm{I}m\left(A_n\right)& =\mathrm{I}m\left(\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}{A}_n\right)\hfill \\ \hfill \text{with}\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}{A}_n& =\frac{l}{2\mathit{\pi }}\left(\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}\left({\left|c_n\right|}^2\right)U_1+{\left|c_n\right|}^2\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}{U}_1\right.\hfill \\ \hfill & \left.-\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\left({\left|d_n\right|}^2\right)U_2-{\left|d_n\right|}^2\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}{U}_2\right),\hfill \\ \hfill \text{where}\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}{U}_1& =\frac{1}{m_{\mathrm{part}}^2}(\left(\mathit{\rho }\left(\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})\right)i\overline{\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})}+{\mathit{\xi }}_n(z_{\mathrm{part}})\mathit{\rho }\overline{\left(\ddot{{\mathit{\xi }}_n}(z_{\mathrm{part}})\right)i}\right)m_{\mathrm{part}}\hfill \\ \hfill & -\left({\mathit{\xi }}_n(z_{\mathrm{part}})\overline{\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})}\right)i),\hfill \\ \hfill \text{and}\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}{U}_2& =\frac{1}{m_{\mathrm{part}}^2}(\left(\mathit{\rho }\left(\ddot{{\mathit{\xi }}_n}(z_{\mathrm{part}})\right)i\overline{{\mathit{\xi }}_n(z_{\mathrm{part}})}+{\dot{\mathit{\xi }}}_n(z_{\mathrm{part}})\mathit{\rho }\overline{\left(\dot{{\mathit{\xi }}_n}(z_{\mathrm{part}})\right)i}\right)m_{\mathrm{part}}\hfill \\ \hfill & -\left({\dot{\mathit{\xi }}}_n(z_{\mathrm{part}})\overline{{\mathit{\xi }}_n(z_{\mathrm{part}})}\right)i).\hfill \end{array}$$

The derivatives of $\mathrm{I}m\left(B_n\right)$ are much easier to compute, since the dependence on $n_{\mathrm{part}}$ and $k_{\mathrm{part}}$ lies only in ${\left|a_n\right|}^2$ and ${\left|b_n\right|}^2$ here. So we can deduce$$\begin{array}{cc}\hfill \frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\mathrm{I}m\left(B_n\right)& =\mathrm{I}m\left(\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}{B}_n\right)\hfill \\ \hfill \text{with}\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}{B}_n& =\frac{l}{2\mathit{\pi }}\left(\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\left({\left|a_n\right|}^2\right)\frac{\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}})\overline{{\mathit{\psi }}_n(z_{\mathrm{med}})}}{m_{\mathrm{med}}}\right.\hfill \\ \hfill & \left.-\frac{\mathit{\partial }}{\mathit{\partial }{n}_{\mathrm{part}}}\left({\left|b_n\right|}^2\right)\frac{{\mathit{\psi }}_n(z_{\mathrm{med}})\overline{\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}})}}{m_{\mathrm{med}}}\right).\hfill \end{array}$$

Analogously we get$$\begin{array}{cc}\hfill \frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}\mathrm{I}m\left(B_n\right)& =\mathrm{I}m\left(\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}{B}_n\right)\hfill \\ \hfill \text{with}\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}{B}_n& =\frac{l}{2\mathit{\pi }}\left(\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}\left({\left|a_n\right|}^2\right)\frac{\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}})\overline{{\mathit{\psi }}_n(z_{\mathrm{med}})}}{m_{\mathrm{med}}}\right.\hfill \\ \hfill & \left.-\frac{\mathit{\partial }}{\mathit{\partial }{k}_{\mathrm{part}}}\left({\left|b_n\right|}^2\right)\frac{{\mathit{\psi }}_n(z_{\mathrm{med}})\overline{\dot{{\mathit{\psi }}_n}(z_{\mathrm{med}})}}{m_{\mathrm{med}}}\right).\hfill \end{array}$$

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(4.1) $$\begin{array}{c}\hfill X_{t,\mathit{\delta }}:=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}\sum _{i=1}^N\frac{1}{2{\mathit{\sigma }}_i^2}{\left(\sum _{n=1}^t{a}_n^i(\mathit{x})-\sum _{n=1}^{\infty }{a}_n^i({\mathit{x}}_{\mathrm{true}})-{\mathit{\delta }}_i\right)}^2,\end{array}$$(4.1)

where $t\in \mathbb{N}$ is a finite truncation index and ${\mathit{\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(\mathit{x})$ for $i=1,\dots ,N$. Throughout this paper, we assume that the feasible set $\mathrm{\Omega }$ is compact.

We define the functions ${\mathit{f}}_t:{\mathbb{R}}^D\to {\mathbb{R}}^N$ and $\mathit{f}:{\mathbb{R}}^D\to {\mathbb{R}}^N$ by$$\begin{array}{cc}\hfill {\mathit{f}}_t(\mathit{x})& :={\left(\sum _{n=1}^t{a}_n^1(\mathit{x}),\dots ,\sum _{n=1}^t{a}_n^N(\mathit{x})\right)}^T\hfill \\ \hfill \text{and}\mathit{f}(\mathit{x})& :={\left(\sum _{n=1}^{\infty }{a}_n^1(\mathit{x}),\dots ,\sum _{n=1}^{\infty }{a}_n^N(\mathit{x})\right)}^T.\hfill \end{array}$$

We set $\mathit{e}:=\mathit{f}({\mathit{x}}_{\mathrm{true}})+\mathit{\delta }$ with $\mathit{\delta }:={\left({\mathit{\delta }}_1,\dots ,{\mathit{\delta }}_N\right)}^T$. Then the observed probability density is given by$$\begin{array}{c}\hfill p_{\mathrm{observed}}(\mathit{e}|\mathit{x}):={(2\mathit{\pi })}^{-\frac{N_l}{2}}|\mathrm{d}et({\mathbf{\Sigma }}_{\mathit{\sigma }}){|}^{-\frac{1}{2}}exp(-\frac{1}{2}\Vert {{\mathbf{\Sigma }}_{\mathit{\sigma }}}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{e}){\Vert }_2^2)\end{array}$$

with the covariance matrix ${\mathbf{\Sigma }}_{\mathit{\sigma }}:=\mathrm{d}iag\left({\mathit{\sigma }}_1^2,\dots ,{\mathit{\sigma }}_N^2\right)$. We know a priori that the vector $\mathit{x}$ specifying our model ${\mathit{f}}_t(\mathit{x})$ lies within the set $\mathrm{\Omega }$. This knowledge can be expressed with the prior probability density$$\begin{array}{c}\hfill p_{\mathrm{prior}}(\mathit{x}):={\left(\mathrm{v}ol(\mathrm{\Omega })\right)}^{-1}{I}_{\mathrm{\Omega }}(\mathit{x}),\end{array}$$

where $I_{\mathrm{\Omega }}$ is the indicator function of $\mathrm{\Omega }$. Now $X_{t,\mathit{\delta }}$ is the set of MAP-estimators of the posterior probability density, i.e.(4.2) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}:=\underset{\mathit{x}}{\mathrm{a}rgmax}& p_{\mathrm{posterior}}(\mathit{x}|\mathit{e})\hfill \\ \hfill \text{with}{p}_{\mathrm{posterior}}(\mathit{x}|\mathit{e})\propto p_{\mathrm{observed}}(\mathit{e}|\mathit{x})p_{\mathrm{prior}}(\mathit{x})& \propto exp(-\frac{1}{2}\Vert {{\mathbf{\Sigma }}_{\mathit{\sigma }}}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{e}){\Vert }_2^2)I_{\mathrm{\Omega }}(\mathit{x}).\hfill \end{array}$$(4.2)

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

Assumption 4.1:

The covariance matrix ${\mathbf{\Sigma }}_{\mathit{\sigma }}$ has the simple form$$\begin{array}{c}\hfill {\mathbf{\Sigma }}_{\mathit{\sigma }}={\mathit{\delta }}^2·\mathrm{d}iag({\mathit{\sigma }}_1^2,\dots ,{\mathit{\sigma }}_N^2)=:{\mathit{\delta }}^2·\mathbf{\Sigma },\end{array}$$

where $\mathit{\delta }\ge 0$ is an arbitrary but fixed noise level and ${\mathit{\sigma }}_1$, ..., ${\mathit{\sigma }}_N$ are fixed.

To simplify notations, we introduce the two functions ${\mathit{f}}_t:{\mathbb{R}}^D\to {\mathbb{R}}^N$ and ${\mathit{g}}_t:{\mathbb{R}}^D\to {\mathbb{R}}^N$ depending on the truncation index t and defined by$$\begin{array}{cc}\hfill {\left({\mathit{f}}_t(\mathit{x})\right)}_i& :=\sum _{n=1}^{\lfloor t\rfloor }{a}_n^i(\mathit{x})+(t-\lfloor t\rfloor )a_{\lfloor t\rfloor +1}^i(\mathit{x})\hfill \\ \hfill \text{and}{\left({\mathit{g}}_t(\mathit{x})\right)}_i& :={\left(\mathit{f}(\mathit{x})\right)}_i-{\left({\mathit{f}}_t(\mathit{x})\right)}_i,\text{for}i=1,\dots ,N.\hfill \end{array}$$

In the following, we will investigate how an element ${\mathit{x}}_{t,\mathit{\delta }}$ of the set $X_{t,\mathit{\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(4.1) $$\begin{array}{c}\hfill X_{t,\mathit{\delta }}:=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}\sum _{i=1}^N\frac{1}{2{\mathit{\sigma }}_i^2}{\left(\sum _{n=1}^t{a}_n^i(\mathit{x})-\sum _{n=1}^{\infty }{a}_n^i({\mathit{x}}_{\mathrm{true}})-{\mathit{\delta }}_i\right)}^2,\end{array}$$(4.1) ) to the new regression problem(4.3) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{t,\mathit{\delta }}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega },\hfill \\ \hfill \text{with}{F}_{t,\mathit{\delta }}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2\hfill \end{array}$$(4.3)

where the truncation index $t\ge 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\to \mathbb{R}$ have a strict local minimum ${\mathit{x}}_0$ inside a compact set $S\subset {\mathbb{R}}^N$. Let the function $h:{\mathbb{R}}^N\times \mathbb{R}\to \mathbb{R}$ have the property ${lim}_{\mathit{\epsilon }\to 0}h(\mathit{x},\mathit{\epsilon })=0$ for all $\mathit{x}\in S$ and let $h(\mathit{x},\mathit{\epsilon })$ be twice continuously differentiable with respect to $\mathit{x}$ and continuous in $\mathit{\epsilon }$. Furthermore, we assume that the local minima ${\mathit{x}}_{\mathit{\epsilon }}$ of $F_{\mathit{\epsilon }}(\mathit{x}):=F(\mathit{x})+h(\mathit{x},\mathit{\epsilon })$ are strict for any $\mathit{\epsilon }>0$. Then there exists a sequence of local minima ${\mathit{x}}_{\mathit{\epsilon }}$ of $F_{\mathit{\epsilon }}(\mathit{x})$ with ${lim}_{\mathit{\epsilon }\to 0}{\mathit{x}}_{\mathit{\epsilon }}={\mathit{x}}_0$.

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

We define $d(\mathit{\epsilon }):={sup}_{\mathit{x}\in S}|h(\mathit{x},\mathit{\epsilon })|$. From ${lim}_{\mathit{\epsilon }\to 0}h(\mathit{x},\mathit{\epsilon })=0$ for all $\mathit{x}\in S$ follows ${lim}_{\mathit{\epsilon }\to 0}d(\mathit{\epsilon })=0$. Let us now introduce the function $F^{-}(\mathit{x}):=F(\mathit{x})-d(\mathit{\epsilon })$. Obviously, ${\mathit{x}}_0$ is also a local minimum of $F^{-}(\mathit{x})$, so for $\mathit{\epsilon }$ sufficiently small there exists a neighbourhood $U_{2d(\mathit{\epsilon })}({\mathit{x}}_0)\subset S$ of ${\mathit{x}}_0$ with$$\begin{array}{c}\hfill F^{-}(\mathit{x})\ge F^{-}({\mathit{x}}_0)\text{and}{F}^{-}(\mathit{x})-F^{-}({\mathit{x}}_0)\le 2d(\mathit{\epsilon })\text{for all}\mathit{x}\in U_{2d(\mathit{\epsilon })}({\mathit{x}}_0).\end{array}$$

In particular, we have$$\begin{array}{c}\hfill \forall \mathit{x}\in \mathit{\partial }{U}_{2d(\mathit{\epsilon })}({\mathit{x}}_0):F^{-}(\mathit{x})=F^{-}({\mathit{x}}_0)+2d(\mathit{\epsilon })=F({\mathit{x}}_0)+d(\mathit{\epsilon }).\end{array}$$

Let us assume that there exists an $\mathit{x}\in \mathit{\partial }{U}_{2d(\mathit{\epsilon })}({\mathit{x}}_0)$ with$$\begin{array}{c}\hfill F_{\mathit{\epsilon }}(\mathit{x})<F({\mathit{x}}_0)+d(\mathit{\epsilon })=F^{-}(\mathit{x})=F(\mathit{x})-d(\mathit{\epsilon }).\end{array}$$

Then $F_{\mathit{\epsilon }}(\mathit{x})=F(\mathit{x})+h(\mathit{x},\mathit{\epsilon })$ implies $-d(\mathit{\epsilon })>h(\mathit{x},\mathit{\epsilon })$, hence $-h(\mathit{x},\mathit{\epsilon })>d(\mathit{\epsilon })\ge -h(\mathit{x},\mathit{\epsilon })$ by definition of $d(\mathit{\epsilon })$, contradiction. Therefore we conclude(4.4) $$\begin{array}{c}\hfill \forall \mathit{x}\in \mathit{\partial }{U}_{2d(\mathit{\epsilon })}({\mathit{x}}_0):F_{\mathit{\epsilon }}(\mathit{x})\ge F({\mathit{x}}_0)+d(\mathit{\epsilon }).\end{array}$$(4.4)

Since $F_{\mathit{\epsilon }}(\mathit{x})$ is continuous and ${\overline{U}}_{2d(\mathit{\epsilon })}({\mathit{x}}_0)$ is compact for $\mathit{\epsilon }$ small enough, there exists an ${\mathit{x}}_{\mathit{\epsilon }}\in {\overline{U}}_{2d(\mathit{\epsilon })}({\mathit{x}}_0)$ with$$\begin{array}{c}\hfill F_{\mathit{\epsilon }}({\mathit{x}}_{\mathit{\epsilon }})=\underset{\mathit{x}\in {\overline{U}}_{2d(\mathit{\epsilon })}({\mathit{x}}_0)}{min}{F}_{\mathit{\epsilon }}(\mathit{x}).\end{array}$$

Let us assume $F_{\mathit{\epsilon }}({\mathit{x}}_{\mathit{\epsilon }})>F({\mathit{x}}_0)+d(\mathit{\epsilon })$. Then by definition of ${\mathit{x}}_{\mathit{\epsilon }}$, we get in particular$$\begin{array}{c}\hfill F({\mathit{x}}_0)+h({\mathit{x}}_0,\mathit{\epsilon })=F_{\mathit{\epsilon }}({\mathit{x}}_0)\ge F_{\mathit{\epsilon }}({\mathit{x}}_{\mathit{\epsilon }})>F({\mathit{x}}_0)+d(\mathit{\epsilon }),\end{array}$$

i.e. $h({\mathit{x}}_0,\mathit{\epsilon })>d(\mathit{\epsilon })\ge h({\mathit{x}}_0,\mathit{\epsilon })$, contradiction. It follows(4.5) $$\begin{array}{c}\hfill F_{\mathit{\epsilon }}({\mathit{x}}_{\mathit{\epsilon }})\le F({\mathit{x}}_0)+d(\mathit{\epsilon })\text{and}{F}_{\mathit{\epsilon }}({\mathit{x}}_0)\le F({\mathit{x}}_0)+d(\mathit{\epsilon }),\end{array}$$(4.5)

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

If it happens to hold that $F_{\mathit{\epsilon }}({\mathit{x}}_{\mathit{\epsilon }})=F({\mathit{x}}_0)+d(\mathit{\epsilon })$, then we also have $F_{\mathit{\epsilon }}({\mathit{x}}_0)=F({\mathit{x}}_0)+d(\mathit{\epsilon })$. Otherwise we have $F_{\mathit{\epsilon }}({\mathit{x}}_{\mathit{\epsilon }})<F({\mathit{x}}_0)+d(\mathit{\epsilon })$ and then (4.4(4.4) $$\begin{array}{c}\hfill \forall \mathit{x}\in \mathit{\partial }{U}_{2d(\mathit{\epsilon })}({\mathit{x}}_0):F_{\mathit{\epsilon }}(\mathit{x})\ge F({\mathit{x}}_0)+d(\mathit{\epsilon }).\end{array}$$(4.4) ) implies that ${\mathit{x}}_{\mathit{\epsilon }}$ cannot lie on $\mathit{\partial }{U}_{2d(\mathit{\epsilon })}({\mathit{x}}_0)$, thus it must lie within the interior of $U_{2d(\mathit{\epsilon })}({\mathit{x}}_0)$. So in any case (4.5(4.5) $$\begin{array}{c}\hfill F_{\mathit{\epsilon }}({\mathit{x}}_{\mathit{\epsilon }})\le F({\mathit{x}}_0)+d(\mathit{\epsilon })\text{and}{F}_{\mathit{\epsilon }}({\mathit{x}}_0)\le F({\mathit{x}}_0)+d(\mathit{\epsilon }),\end{array}$$(4.5) ) gives that $U_{2d(\mathit{\epsilon })}({\mathit{x}}_0)$ must contain a local minimizer ${\mathit{x}}_{\mathit{\epsilon }}$ of $F_{\mathit{\epsilon }}(\mathit{x})$.

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

Proposition 4.3:

Let all coefficient functions $a_n^i(\mathit{x})$ be twice continuously differentiable and bounded on $\mathrm{\Omega }$. We assume that each local minimum ${\mathit{x}}_{t,\mathit{\delta }}$ of the right-hand side function $F_{t,\mathit{\delta }}(\mathit{x})$ in (4.3(4.3) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{t,\mathit{\delta }}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega },\hfill \\ \hfill \text{with}{F}_{t,\mathit{\delta }}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2\hfill \end{array}$$(4.3) ) is strict and lies in the interior of $\mathrm{\Omega }$. Then each local minimum depends continuously on the truncation index t.

To prove the claim, one could be tempted to apply the implicit function theorem on the equation $\mathit{d}(t,{\mathit{x}}_{t,\mathit{\delta }})=0$ with $\mathit{d}(s,\mathit{x}):=\mathrm{\nabla }{\mathit{f}}_s(\mathit{x})$. This would give that the local minima are parameterized by a function $\mathit{m}(s)$ with the property $\mathit{m}(t)={\mathit{x}}_{t,\mathit{\delta }}$, where s is from an environment U(t) of t. The problem with this approach is that it requires continuous differentiability of $\mathit{d}(s,\mathit{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 $\mathit{\epsilon }>0$ be arbitrary. First we consider an integer truncation index $t\in \mathbb{N}$, i.e. we have $t=\lfloor t\rfloor$. Now for $\mathit{\epsilon }$ small enough, we get $\lfloor t+\mathit{\epsilon }\rfloor =t$ and $\lfloor t-\mathit{\epsilon }\rfloor =t-1$. This gives$$\begin{array}{cc}\hfill {\left({\mathit{f}}_{t+\mathit{\epsilon }}(\mathit{x})\right)}_i& ={\left({\mathit{f}}_t(\mathit{x})\right)}_i+\mathit{\epsilon }{a}_{t+1}^i(\mathit{x})\hfill \\ \hfill \text{and}{\left({\mathit{f}}_{t-\mathit{\epsilon }}(\mathit{x})\right)}_i& ={\left({\mathit{f}}_{t-1}(\mathit{x})\right)}_i+(1-\mathit{\epsilon })a_t^i(\mathit{x})\hfill \\ \hfill & ={\left({\mathit{f}}_t(\mathit{x})\right)}_i-\mathit{\epsilon }{a}_t^i(\mathit{x}).\hfill \end{array}$$

As next step we turn to an noninteger truncation index t. In this case, we can always select $\mathit{\epsilon }$ small enough such that $\lfloor t+\mathit{\epsilon }\rfloor =\lfloor t\rfloor$ and $\lfloor t-\mathit{\epsilon }\rfloor =\lfloor t\rfloor$, respectively hold. This yields$$\begin{array}{cc}\hfill {\left({\mathit{f}}_{t+\mathit{\epsilon }}(\mathit{x})\right)}_i& ={\left({\mathit{f}}_t(\mathit{x})\right)}_i+\mathit{\epsilon }{a}_{\lfloor t\rfloor +1}^i(\mathit{x})\hfill \\ \hfill \text{and}{\left({\mathit{f}}_{t-\mathit{\epsilon }}(\mathit{x})\right)}_i& ={\left({\mathit{f}}_t(\mathit{x})\right)}_i-\mathit{\epsilon }{a}_{\lfloor t\rfloor +1}^i(\mathit{x}).\hfill \end{array}$$

Now we introduce the function$$\begin{array}{c}\hfill \mathit{a}(\mathit{x}):=\left\{\begin{array}{cc}(a_t^1(\mathit{x}),\dots ,a_t^N(\mathit{x}){)}^T,\hfill & \text{for}t-\mathit{\epsilon },t\in \mathbb{N}\hfill \\ (a_{\lfloor t\rfloor +1}^1(\mathit{x}),\dots ,a_{\lfloor t\rfloor +1}^N(\mathit{x}){)}^T,\hfill & \text{else}.\hfill \end{array}\right.\end{array}$$

For $F_{t,\mathit{\delta }}(\mathit{x})={\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta })\Vert }_2^2$ this yields$$\begin{array}{cc}\hfill F_{t+\mathit{\epsilon },\mathit{\delta }}(\mathit{x})& =F_{t,\mathit{\delta }}(\mathit{x})+2\mathit{\epsilon }〈{\mathbf{\Sigma }}^{-\frac{1}{2}}\mathit{a}(\mathit{x}),{\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta })〉+{\mathit{\epsilon }}^2{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}\mathit{a}(\mathit{x})\Vert }_2^2\hfill \\ \hfill \text{and}{F}_{t-\mathit{\epsilon },\mathit{\delta }}(\mathit{x})& =F_{t,\mathit{\delta }}(\mathit{x})-2\mathit{\epsilon }〈{\mathbf{\Sigma }}^{-\frac{1}{2}}\mathit{a}(\mathit{x}),{\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta })〉+{\mathit{\epsilon }}^2{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}\mathit{a}(\mathit{x})\Vert }_2^2.\hfill \end{array}$$

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

Corollary 4.4:

Let $t_1$ and $t_2$ be truncation indices with $t_1<t_2$. Let ${\mathit{x}}_{t_1,\mathit{\delta }}$ be a local minimizer of (4.1(4.1) $$\begin{array}{c}\hfill X_{t,\mathit{\delta }}:=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}\sum _{i=1}^N\frac{1}{2{\mathit{\sigma }}_i^2}{\left(\sum _{n=1}^t{a}_n^i(\mathit{x})-\sum _{n=1}^{\infty }{a}_n^i({\mathit{x}}_{\mathrm{true}})-{\mathit{\delta }}_i\right)}^2,\end{array}$$(4.1) ). Let $\mathit{\gamma }\in [0,1]$ and define $t_{\mathit{\gamma }}:=t_1+\mathit{\gamma }\left(t_2-t_1\right)$. Then beginning at $\mathit{\gamma }=0$ one can successively find local minimizers ${\mathit{x}}_{t_{\mathit{\gamma }},\mathit{\delta }}$ for the truncation index $t_{\mathit{\gamma }}$ using numerical continuation, see [6]. Here for ${\mathit{\gamma }}_1<{\mathit{\gamma }}_2$ the minimizer ${\mathit{x}}_{t_{{\mathit{\gamma }}_1},\mathit{\delta }}$ is used as a start vector to compute the next minimizer ${\mathit{x}}_{t_{{\mathit{\gamma }}_2},\mathit{\delta }}$. The next parameter ${\mathit{\gamma }}_2$ has to be sufficiently close to ${\mathit{\gamma }}_1$, such that the start vector ${\mathit{x}}_{t_{{\mathit{\gamma }}_1},\mathit{\delta }}$ still lies within the domain of convergence for Newton’s method.

We use Corollary 4.4 to compute ${\mathit{x}}_{t,\mathit{\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 ${\mathit{x}}_{t,\mathit{\delta }}$ of the noise-contaminated regression problem (4.3(4.3) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{t,\mathit{\delta }}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega },\hfill \\ \hfill \text{with}{F}_{t,\mathit{\delta }}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2\hfill \end{array}$$(4.3) ) with truncated series expansions approximate the minimizers ${\mathit{x}}_{\infty ,0}$ of the noise-free and untruncated problem(4.6) $$\begin{array}{c}\hfill {\mathit{x}}_{\infty ,0}:=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}\sum _{i=1}^N\frac{1}{2{\mathit{\sigma }}_i^2}{\left(\sum _{n=1}^{\infty }{a}_n^i(\mathit{x})-\sum _{n=1}^{\infty }{a}_n^i({\mathit{x}}_{\mathrm{true}})\right)}^2\text{s.t.}\mathit{x}\in \mathrm{\Omega }.\end{array}$$(4.6)

Proposition 4.5:

Let the noise vector $\mathit{\delta }$ fulfil ${lim}_{\mathit{\delta }\to 0}{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}\mathit{\delta }\Vert }_2^2=0$ and let the functions $\mathit{f}(\mathit{x})$ and ${\mathit{f}}_{t_{\mathit{\delta }}}(\mathit{x})$ be bounded on $\mathrm{\Omega }$. Assume ${lim}_{\mathit{\delta }\to 0}{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}(\mathit{x})\Vert }_2^2=0$ for all $\mathit{x}\in \mathrm{\Omega }$. Then for any strict minimizer ${\mathit{x}}_{\infty ,0}$ of the right-hand side function of (4.6(4.6) $$\begin{array}{c}\hfill {\mathit{x}}_{\infty ,0}:=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}\sum _{i=1}^N\frac{1}{2{\mathit{\sigma }}_i^2}{\left(\sum _{n=1}^{\infty }{a}_n^i(\mathit{x})-\sum _{n=1}^{\infty }{a}_n^i({\mathit{x}}_{\mathrm{true}})\right)}^2\text{s.t.}\mathit{x}\in \mathrm{\Omega }.\end{array}$$(4.6) ) in the interior of $\mathrm{\Omega }$ exist minimizers ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}$ of (4.3(4.3) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{t,\mathit{\delta }}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega },\hfill \\ \hfill \text{with}{F}_{t,\mathit{\delta }}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2\hfill \end{array}$$(4.3) ) with ${lim}_{\mathit{\delta }\to 0}{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}={\mathit{x}}_{\infty ,0}$. Here we also assume the ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}$’s to be strict for all $\mathit{\delta }>0$.

With the notation introduced before, we can write$$\begin{array}{cc}\hfill {\mathit{x}}_{\infty ,0}\in X_{\infty ,0}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{\infty ,0}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega }\hfill \\ \hfill \text{with}{F}_{\infty ,0}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}(\mathit{f}(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})){\Vert }_2^2.\hfill \end{array}$$

From the decomposition ${\mathit{f}}_{t_{\mathit{\delta }}}(\mathit{x})=\mathit{f}(\mathit{x})-{\mathit{g}}_{t_{\mathit{\delta }}}(\mathit{x})$ we obtain$$\begin{array}{cc}\hfill F_{t_{\mathit{\delta }},\mathit{\delta }}(\mathit{x})& =\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}\left(\mathit{f}(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }\right){\Vert }_2^2\hfill \\ \hfill & -2〈{\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}(\mathit{x}),{\mathbf{\Sigma }}^{-\frac{1}{2}}\left(\mathit{f}(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }\right)〉\hfill \\ \hfill & +\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}{(\mathit{x})\Vert }_2^2.\hfill \end{array}$$

Then a further decomposition of the first term on the right-hand side yields$$\begin{array}{cc}\hfill F_{t_{\mathit{\delta }},\mathit{\delta }}(\mathit{x})& =F_{\infty ,0}(\mathit{x})-2〈{\mathbf{\Sigma }}^{-\frac{1}{2}}\mathit{\delta },{\mathbf{\Sigma }}^{-\frac{1}{2}}\left(\mathit{f}(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})\right)〉+{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}\mathit{\delta }\Vert }_2^2\hfill \\ \hfill & -2〈{\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}(\mathit{x}),{\mathbf{\Sigma }}^{-\frac{1}{2}}\left(\mathit{f}(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }\right)〉+{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}(\mathit{x})\Vert }_2^2\hfill \\ \hfill & =:F_{\infty ,0}(\mathit{x})+H_{t_{\mathit{\delta }},\mathit{\delta }}(\mathit{x}).\hfill \end{array}$$

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

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

At last we study how the minimizers ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}$ of (4.3(4.3) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{t,\mathit{\delta }}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega },\hfill \\ \hfill \text{with}{F}_{t,\mathit{\delta }}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2\hfill \end{array}$$(4.3) ) behave for $\mathit{\delta }\to 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 ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}$ of (4.3(4.3) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{t,\mathit{\delta }}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega },\hfill \\ \hfill \text{with}{F}_{t,\mathit{\delta }}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2\hfill \end{array}$$(4.3) ) approximating a local minimizer ${\mathit{x}}_{\infty ,0}$ of (4.6(4.6) $$\begin{array}{c}\hfill {\mathit{x}}_{\infty ,0}:=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}\sum _{i=1}^N\frac{1}{2{\mathit{\sigma }}_i^2}{\left(\sum _{n=1}^{\infty }{a}_n^i(\mathit{x})-\sum _{n=1}^{\infty }{a}_n^i({\mathit{x}}_{\mathrm{true}})\right)}^2\text{s.t.}\mathit{x}\in \mathrm{\Omega }.\end{array}$$(4.6) ) for $\mathit{\delta }\to 0$ with $\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}(\mathit{f}({\mathit{x}}_{\infty ,0})-\mathit{f}({\mathit{x}}_{\mathrm{true}})){\Vert }_2=0$ that$$\begin{array}{c}\hfill \underset{\mathit{\delta }\to 0}{lim}{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta })\Vert }_2=0.\end{array}$$

The assumptions of Proposition 4.5 give$$\begin{array}{c}\hfill \underset{\mathit{\delta }\to 0}{lim}\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{\delta }\Vert }_2^2=0\text{and}\underset{\mathit{\delta }\to 0}{lim}{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }})\Vert }_2^2=0.\end{array}$$

We have by continuity of $\mathit{f}(\mathit{x})$ that ${lim}_{\mathit{\delta }\to 0}{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}(\mathit{f}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }})-\mathit{f}({\mathit{x}}_{\mathrm{true}}))\Vert }_2=0$. Then$$\begin{array}{cc}& \Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}(\mathit{f}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }})-\mathit{f}({\mathit{x}}_{\mathrm{true}})){\Vert }_2+\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}){\Vert }_2+{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}\mathit{\delta }\Vert }_2\hfill \\ \hfill & \ge \Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2\hfill \end{array}$$

gives the desired result.

Proposition 4.7:

Let the assumptions of Proposition 4.5 hold. Assume that the local minimizers ${\mathit{x}}_{\infty ,0}$ of (4.6(4.6) $$\begin{array}{c}\hfill {\mathit{x}}_{\infty ,0}:=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}\sum _{i=1}^N\frac{1}{2{\mathit{\sigma }}_i^2}{\left(\sum _{n=1}^{\infty }{a}_n^i(\mathit{x})-\sum _{n=1}^{\infty }{a}_n^i({\mathit{x}}_{\mathrm{true}})\right)}^2\text{s.t.}\mathit{x}\in \mathrm{\Omega }.\end{array}$$(4.6) ) with $\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}(\mathit{f}({\mathit{x}}_{\infty ,0})-\mathit{f}({\mathit{x}}_{\mathrm{true}})){\Vert }_2=0$ form a discrete set $S_{\infty ,0}$. Then the set $L_{\infty ,0}$ consisting of the limits ${lim}_{\mathit{\delta }\to 0}{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}$ of local minimizers ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}$ of (4.3(4.3) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{t,\mathit{\delta }}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega },\hfill \\ \hfill \text{with}{F}_{t,\mathit{\delta }}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2\hfill \end{array}$$(4.3) ) with ${lim}_{\mathit{\delta }\to 0}{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta })\Vert }_2=0$ coincides with $S_{\infty ,0}$ and there exists a noise level ${\mathit{\delta }}_{max}$ such that all minimizers ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}$ approximating $S_{\infty ,0}$ are isolated for all $\mathit{\delta }\le {\mathit{\delta }}_{max}$.

On the one hand from Proposition 4.5 we know that there exists a sequence ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}$ of minimizers of (4.3(4.3) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{t,\mathit{\delta }}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega },\hfill \\ \hfill \text{with}{F}_{t,\mathit{\delta }}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2\hfill \end{array}$$(4.3) ) with ${lim}_{\mathit{\delta }\to 0}{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}={\mathit{x}}_{\infty ,0}$. Then $\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}(\mathit{f}({\mathit{x}}_{\infty ,0})-\mathit{f}({\mathit{x}}_{\mathrm{true}})){\Vert }_2=0$ and Corollary 4.6 give ${lim}_{\mathit{\delta }\to 0}{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta })\Vert }_2=0$, which implies $S_{\infty ,0}\subseteq L_{\infty ,0}$.

On the other hand holds for ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}$ with ${lim}_{\mathit{\delta }\to 0}{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta })\Vert }_2=0$ that$$\begin{array}{cc}& \Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2+\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}){\Vert }_2+{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}\mathit{\delta }\Vert }_2\hfill \\ \hfill & \ge \Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}(\mathit{f}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }})-\mathit{f}({\mathit{x}}_{\mathrm{true}})){\Vert }_2\hfill \end{array}$$

which implies ${lim}_{\mathit{\delta }\to 0}{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}(\mathit{f}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }})-\mathit{f}({\mathit{x}}_{\mathrm{true}}))\Vert }_2=0$. In particular, this means by continuity of $\mathit{f}(\mathit{x})$ that the vector ${lim}_{\mathit{\delta }\to 0}{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}$ must be a local minimizer of $\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}(\mathit{f}(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})){\Vert }_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 ${\mathit{x}}_{\infty ,0}^k$ for $k=1,\dots ,|S_{\infty ,0}|$. Similarly we number all minimizers ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}$ with ${lim}_{\mathit{\delta }\to 0}{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta })\Vert }_2=0$ approximating the ${\mathit{x}}_{\infty ,0}^k$’s with ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k$, i.e. ${lim}_{\mathit{\delta }\to 0}{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k={\mathit{x}}_{\infty ,0}^k$ for $k=1,\dots ,|S_{\infty ,0}|$. Define$$\begin{array}{c}\hfill D_{min}:=\underset{i\ne j}{min}{\Vert {\mathit{x}}_{\infty ,0}^i-{\mathit{x}}_{\infty ,0}^j\Vert }_2.\end{array}$$

Since ${lim}_{\mathit{\delta }\to 0}{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k={\mathit{x}}_{\infty ,0}^k$, we can find an error levels ${\mathit{\delta }}_{max}^k$ such that$$\begin{array}{c}\hfill \Vert {\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k-{\mathit{x}}_{\infty ,0}^k{\Vert }_2<\frac{1}{2}{D}_{min}\text{for}k=1,\dots ,|S_{\infty ,0}|,\end{array}$$

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

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

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

Proposition 4.8:

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

The first-order necessary conditions for a local minimum of $F_{t_{\mathit{\delta }},\mathit{\delta }}(\mathit{x})=F_{\infty ,0}(\mathit{x})+H_{t_{\mathit{\delta }},\mathit{\delta }}(\mathit{x})$ at ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k$ and a local minimum of $F_{\infty ,0}(\mathit{x})$ at ${\mathit{x}}_{\infty ,0}^k$ yield in particular$$\begin{array}{cc}\hfill 〈\mathrm{\nabla }{F}_{\infty ,0}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k)+\mathrm{\nabla }{H}_{t_{\mathit{\delta }},\mathit{\delta }}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k),{\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k〉& \ge 0\hfill \\ \hfill \text{and}〈\mathrm{\nabla }{F}_{\infty ,0}({\mathit{x}}_{\infty ,0}^k),{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k-{\mathit{x}}_{\infty ,0}^k〉& \ge 0,\hfill \end{array}$$

Adding the last two inequalities yields(4.7) $$\begin{array}{c}\hfill 〈\mathrm{\nabla }{H}_{t_{\mathit{\delta }},\mathit{\delta }}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k),{\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k〉\ge 〈\mathrm{\nabla }{F}_{\infty ,0}({\mathit{x}}_{\infty ,0}^k)-\mathrm{\nabla }{F}_{\infty ,0}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k),{\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k〉.\end{array}$$(4.7)

Since $\mathrm{\nabla }{F}_{\infty ,0}(\mathit{x})$ is totally differentiable at ${\mathit{x}}_{\infty ,0}^k$, we obtain$$\begin{array}{c}\hfill \mathrm{\nabla }{F}_{\infty ,0}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k)=\mathrm{\nabla }{F}_{\infty ,0}({\mathit{x}}_{\infty ,0}^k)+\mathrm{H}es{s}_{F_{\infty ,0}}({\mathit{x}}_{\infty ,0}^k)\left({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k-{\mathit{x}}_{\infty ,0}^k\right)+{\mathit{w}}_{\infty ,0}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k,{\mathit{x}}_{\infty ,0}^k),\end{array}$$

where ${\mathit{w}}_{\infty ,0}(\mathit{x},{\mathit{x}}_{\infty ,0}^k)$ fulfils$$\begin{array}{c}\hfill \Vert {\mathit{w}}_{\infty ,0}(\mathit{x},{\mathit{x}}_{\infty ,0}^k){\Vert }_2\le {\Vert \mathit{x}-{\mathit{x}}_{\infty ,0}^k\Vert }_2{\mathit{ϵ}}_{\infty ,0}(\mathit{x},{\mathit{x}}_{\infty ,0}^k)\text{with}\underset{\mathit{x}\to {\mathit{x}}_{\infty ,0}^k}{lim}{\mathit{ϵ}}_{\infty ,0}(\mathit{x},{\mathit{x}}_{\infty ,0}^k)=0.\end{array}$$

Since $\mathrm{H}es{s}_{F_{\infty ,0}}({\mathit{x}}_{\infty ,0}^k)$ is positive definite, the expression ${\left(〈\mathit{x},\mathrm{H}es{s}_{F_{\infty ,0}}({\mathit{x}}_{\infty ,0}^k)\mathit{x}〉\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$$\begin{array}{c}\hfill {\left(〈\mathit{x},\mathrm{H}es{s}_{F_{\infty ,0}}({\mathit{x}}_{\infty ,0}^k)\mathit{x}〉\right)}^{\frac{1}{2}}\ge C_{\infty ,0}^k{\Vert \mathit{x}\Vert }_2\text{for all}\mathit{x}\in {\mathbb{R}}^D.\end{array}$$

Since ${lim}_{\mathit{\delta }\to 0}{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k={\mathit{x}}_{\infty ,0}^k$ we can find a noise level ${\mathit{\rho }}_{max}^k$ such that $|{\mathit{ϵ}}_{\infty ,0}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k,{\mathit{x}}_{\infty ,0}^k)|\le d_{\infty ,0}^k$ for all $\mathit{\delta }\le {\mathit{\rho }}_{max}^k$, where $d_{\infty ,0}^k$ is a constant with $0\le d_{\infty ,0}^k<{\left(C_{\infty ,0}^k\right)}^2$. Then using$$\begin{array}{c}\hfill 〈{\mathit{w}}_{\infty ,0}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k,{\mathit{x}}_{\infty ,0}^k),{\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k〉\le {\mathit{ϵ}}_{\infty ,0}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k,{\mathit{x}}_{\infty ,0}^k){\Vert {\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k\Vert }_2^2\end{array}$$

and (4.7(4.7) $$\begin{array}{c}\hfill 〈\mathrm{\nabla }{H}_{t_{\mathit{\delta }},\mathit{\delta }}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k),{\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k〉\ge 〈\mathrm{\nabla }{F}_{\infty ,0}({\mathit{x}}_{\infty ,0}^k)-\mathrm{\nabla }{F}_{\infty ,0}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k),{\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k〉.\end{array}$$(4.7) ) we can estimate$$\begin{array}{cc}& \Vert \mathrm{\nabla }{H}_{t_{\mathit{\delta }},\mathit{\delta }}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k){\Vert }_2{\Vert {\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k\Vert }_2\hfill \\ \hfill & \ge 〈\mathrm{\nabla }{F}_{\infty ,0}({\mathit{x}}_{\infty ,0}^k)-\mathrm{\nabla }{F}_{\infty ,0}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k),{\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k〉\hfill \\ \hfill & =〈\mathrm{H}es{s}_{F_{\infty ,0}}({\mathit{x}}_{\infty ,0}^k)\left({\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k\right)-{\mathit{w}}_{\infty ,0}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k,{\mathit{x}}_{\infty ,0}^k),{\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k〉\hfill \\ \hfill & \ge \left((C_{\infty ,0}^k{)}^2-{\mathit{ϵ}}_{\infty ,0}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k,{\mathit{x}}_{\infty ,0}^k)\right){\Vert {\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k\Vert }_2^2,\hfill \end{array}$$ i.e. this gives(4.8) $$\begin{array}{c}\hfill \Vert {\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k{\Vert }_2\le {\left((C_{\infty ,0}^k{)}^2-d_{\infty ,0}^k\right)}^{-1}{\Vert \mathrm{\nabla }{H}_{t_{\mathit{\delta }},\mathit{\delta }}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k)\Vert }_2\end{array}$$(4.8)

for all $\mathit{\delta }\le {\mathit{\rho }}_{max}^k$.

We have$$\begin{array}{cc}\hfill \mathrm{\nabla }{H}_{t_{\mathit{\delta }},\mathit{\delta }}(\mathit{x})=2(& \mathrm{J}a{c}_{{\mathit{g}}_{t_{\mathit{\delta }}}}^T(\mathit{x}){\mathbf{\Sigma }}^{-1}{\mathit{g}}_{t_{\mathit{\delta }}}(\mathit{x})-\mathrm{J}a{c}_{{\mathit{g}}_{t_{\mathit{\delta }}}}^T(\mathit{x}){\mathbf{\Sigma }}^{-1}(\mathit{f}(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta })\hfill \\ \hfill -& \mathrm{J}a{c}_{\mathit{f}}^T(\mathit{x}){\mathbf{\Sigma }}^{-1}{\mathit{g}}_{t_{\mathit{\delta }}}(\mathit{x})-\mathrm{J}a{c}_{\mathit{f}}^T(\mathit{x}){\mathbf{\Sigma }}^{-1}\mathit{\delta }),\hfill \end{array}$$

i.e. we can find constants $K_1^k$ and $K_2^k$ such that$$\begin{array}{cc}\hfill \Vert \mathrm{\nabla }{H}_{t_{\mathit{\delta }},\mathit{\delta }}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k){\Vert }_2& \le 2(K_1^k\Vert \mathrm{J}a{c}_{{\mathit{g}}_{t_{\mathit{\delta }}}}^T({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k){\mathbf{\Sigma }}^{-\frac{1}{2}}{\Vert }_{\infty }{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k)\Vert }_2\hfill \\ \hfill & +K_1^k\Vert \mathrm{J}a{c}_{{\mathit{g}}_{t_{\mathit{\delta }}}}^T({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k){\mathbf{\Sigma }}^{-\frac{1}{2}}{\Vert }_{\infty }({\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k)-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta })\Vert }_2\hfill \\ \hfill & +\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k){\Vert }_2)\hfill \\ \hfill & +K_2^k\Vert \mathrm{J}a{c}_{\mathit{f}}^T({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k){\mathbf{\Sigma }}^{-\frac{1}{2}}{\Vert }_{\infty }{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k)\Vert }_2\hfill \\ \hfill & +K_2^k\Vert \mathrm{J}a{c}_{\mathit{f}}^T({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k){\mathbf{\Sigma }}^{-\frac{1}{2}}{\Vert }_{\infty }{\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}\mathit{\delta }\Vert }_2),\hfill \end{array}$$

which gives the result.

Corollary 4.9:

Let the truncation indices $t_{\mathit{\delta }}$ depend on the vector of independent Gaussian random variables $\mathit{\delta }$ with ${lim}_{\mathit{\delta }\to 0}\mathbb{E}({\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}\mathit{\delta }\Vert }_2^2)=0$ such that ${lim}_{\mathit{\delta }\to 0}\mathbb{E}({\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}(\mathit{x})\Vert }_2^2)=0$ holds for all arbitrary but fixed $\mathit{x}\in \mathrm{\Omega }$. Then we have for all minimizers ${\mathit{x}}_{\infty ,0}^k$ of (4.6(4.6) $$\begin{array}{c}\hfill {\mathit{x}}_{\infty ,0}:=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}\sum _{i=1}^N\frac{1}{2{\mathit{\sigma }}_i^2}{\left(\sum _{n=1}^{\infty }{a}_n^i(\mathit{x})-\sum _{n=1}^{\infty }{a}_n^i({\mathit{x}}_{\mathrm{true}})\right)}^2\text{s.t.}\mathit{x}\in \mathrm{\Omega }.\end{array}$$(4.6) ) with $\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}(\mathit{f}({\mathit{x}}_{\infty ,0}^k)-\mathit{f}({\mathit{x}}_{\mathrm{true}})){\Vert }_2=0$ that for $\mathit{\delta }$ sufficiently small there exist minimizers ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k$ of (4.3(4.3) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{t,\mathit{\delta }}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega },\hfill \\ \hfill \text{with}{F}_{t,\mathit{\delta }}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2\hfill \end{array}$$(4.3) ) with$$\begin{array}{c}\hfill \underset{\mathit{\delta }\to 0}{lim}\mathbb{E}({\Vert {\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k\Vert }_2)=0.\end{array}$$

Proposition 4.5 establishes the existence of the ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k$’s. Corollary 4.6 gives$$\begin{array}{c}\hfill \underset{\mathit{\delta }\to 0}{lim}\mathbb{E}({\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k)-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta })\Vert }_2)=0.\end{array}$$

Proposition 4.7 gives that for $\mathit{\delta }$ sufficiently small there exists a constant $K_{\infty ,0}^k$ with$$\begin{array}{c}\hfill \mathbb{E}(\Vert {\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k{\Vert }_2)\le K_{\infty ,0}^k\mathbb{E}({\Vert \mathrm{\nabla }{H}_{t_{\mathit{\delta }},\mathit{\delta }}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k)\Vert }_2).\end{array}$$

Set $S_1^k:={sup}_{\mathit{x}\in \mathrm{\Omega }}{\Vert \mathrm{J}a{c}_{{\mathit{g}}_{t_{\mathit{\delta }}}}^T(\mathit{x}){\mathbf{\Sigma }}^{-\frac{1}{2}}\Vert }_{\infty }<\infty$ and $S_2^k:={sup}_{\mathit{x}\in \mathrm{\Omega }}{\Vert \mathrm{J}a{c}_{\mathit{f}}^T(\mathit{x}){\mathbf{\Sigma }}^{-\frac{1}{2}}\Vert }_{\infty }<\infty$. Then the estimate for $\Vert \mathrm{\nabla }{H}_{t_{\mathit{\delta }},\mathit{\delta }}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k){\Vert }_2$ in the proof of Proposition 4.7 gives$$\begin{array}{cc}\hfill \mathbb{E}(\Vert \mathrm{\nabla }{H}_{t_{\mathit{\delta }},\mathit{\delta }}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k){\Vert }_2)& \le 2(K_1^k{S}_1^k\mathbb{E}({\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k)\Vert }_2)\hfill \\ \hfill & +K_1^k{S}_1^k(\mathbb{E}({\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k)-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta })\Vert }_2)\hfill \\ \hfill & +\mathbb{E}(\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k){\Vert }_2))\hfill \\ \hfill & +K_2^k{S}_2^k\mathbb{E}({\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k)\Vert }_2)\hfill \\ \hfill & +K_2^k{S}_2^k\mathbb{E}({\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}\mathit{\delta }\Vert }_2)),\hfill \end{array}$$ which proves the claim since Assumption 4.1 gives $\mathbb{E}(\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{\delta }\Vert }_2)\le \sqrt{N}\mathit{\delta }$.

The strategy for our retrieval algorithm is to start with an initial guess for the truncation index $t_{\mathrm{start}}$ and try to find all local minima ${\mathit{x}}_{t_{\mathrm{start}},\mathit{\delta }}^k$. Then the truncation index is gradually increased and starting from ${\mathit{x}}_{t_{\mathrm{start}},\mathit{\delta }}^k$ the continuation method is applied to find finally the local minima ${\mathit{x}}_{t_{\mathit{\delta }},\mathit{\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 $\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k)-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2$ and an estimate of the truncation error $\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}{\mathit{g}}_{t_{\mathit{\delta }}}({\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k){\Vert }_2$ are both reasonably small. The latter serves also as a stopping criterion for the continuation method

The initial guess $t_{\mathrm{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 cannot 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.

Remark 4.10:

We did not explicitly show the assumption that the local minima of (4.3(4.3) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{t,\mathit{\delta }}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega },\hfill \\ \hfill \text{with}{F}_{t,\mathit{\delta }}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2\hfill \end{array}$$(4.3) ) are strict as demanded in Proposition 4.5. However we could always determine a set of strict local minima for any truncation index in our simulations.

5. The reconstruction algorithm

We now return to our regression problem (1.3(1.3) $$\begin{array}{c}\hfill M(l):=\underset{m\in \mathbb{C}}{\mathrm{a}rgmin}\sum _{i=1}^N\frac{1}{2{\left(\frac{s_i}{n_i}\right)}^2}{\left(\mathit{\pi }{r}_i^2\sum _{n=1}^{N_{tr}}{q}_n(m_{\mathrm{med}}(l),m,r_i,l)-\frac{q(r_i,l)}{n_i}\right)}^2.\end{array}$$(1.3) ). For $i=1,\dots ,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}_{\mathrm{true}}(r_i,l)$ and standard deviation ${\mathit{\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_{\mathrm{part}}(l)$ wavelength by wavelength. In the following, the unit both for particle radii and light wavelengths is $\mathrm{\mu }$m.

We make use of the function ${\mathit{q}}_{N_{tr}}:{\mathbb{R}}^2\to {\mathbb{R}}^N$ defined by(5.1) $$\begin{array}{c}\hfill {\mathit{q}}_{N_{tr}}(\mathit{x}):={\left(\mathit{\pi }{r}_1^2\sum _{n=1}^{N_{tr}}{q}_n(\mathit{x},r_1,l),\dots ,\mathit{\pi }{r}_N^2\sum _{n=1}^{N_{tr}}{q}_n(\mathit{x},r_N,l)\right)}^T,\end{array}$$(5.1)

where $R=\{r_1,\dots ,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(\mathit{x},r_k,l)$ is a short notation of $q_n(m_{\mathrm{med}}(l),{(\mathit{x})}_1+{(\mathit{x})}_{2}i,r_k,l)$ from Section 1, where the sought-after refrative index $m_{\mathrm{part}}(l)$ is identified with the vector $\mathit{x}$ here, i.e. $m_{\mathrm{part}}(l)={(\mathit{x})}_1+{(\mathit{x})}_{2}i$. So its computation follows Section 2.

In the following, the refractive index search area is given by the rectangle $\mathrm{\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.

In the first loop from lines 14–25, a search for local minima of the fit function $F(\mathit{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 and 6. If the Hessian of $F(\mathit{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_{\mathrm{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_{\mathrm{rel}}\le To{l}_{\mathrm{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 $\mathit{\tau }=5$ and $\mathit{\tau }=7$, respectively, if none of the reconstructions had a squared residual smaller than $\mathit{\tau }{N}_r{\mathit{\delta }}^2$ for the previous $\mathit{\tau }$. This had to be done, because the parameter $\mathit{\tau }$ has to be selected carefully in order to estimate the bound on $\mathbb{E}(\Vert {\mathit{x}}_{\infty ,0}^k-{\mathit{x}}_{t_{\mathit{\delta }},\mathit{\delta }}^k{\Vert }_2)$ 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$$\begin{array}{c}\hfill {\left({\mathit{e}}_{\mathrm{true}}\right)}_{i,j}:=\mathit{\pi }{r}_i^2\sum _{n=1}^{N_{tr}}{q}_n(m_{\mathrm{med}}(l_j),m_{\mathrm{part}}(l_j),r_i,l_j),\text{for}i=1,\dots ,N,j=1,\dots ,N_l\end{array}$$

with $m_{\mathrm{part}}(l_i)$ taken as the refractive indices of Ag, ${\mathrm{H}}_2$O and CsI. Here we used the truncation index(6.1) $$\begin{array}{cc}\hfill \mathit{\rho }& =2\mathit{\pi }\frac{r}{l},\hfill \\ \hfill M& =max\{|\mathit{\rho }|,|\mathit{\rho }·m_{\mathrm{part}}(l)|,|\mathit{\rho }·m_{\mathrm{med}}(l)|\},\hfill \\ \hfill N_{tr}& :=\lceil |M+4.05·M^{\frac{1}{3}}+2|\rceil \hfill \end{array}$$(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 $\mathrm{\mu }$m, 8 from 1.1–1.3$\mathrm{\mu }$m, 8 from 1.6–1.8$\mathrm{\mu }$m, 16 from 2.1–2.5$\mathrm{\mu }$m and 8 from 3.1–3.3$\mathrm{\mu }$m, so we have in total $N_l=48$ wavelengths.

For each of the 48 wavelengths, we generated noisy spectral extinctions $\mathit{e}$ by adding zero-mean Gaussian noise to ${\mathit{e}}_{\mathrm{true}}$, i.e.$$\begin{array}{c}\hfill {\left(\mathit{e}\right)}_{i,j}={\left({\mathit{e}}_{\mathrm{true}}\right)}_{i,j}+{\mathit{\delta }}_{i,j}\text{with}{\mathit{\delta }}_{i,j}\sim \mathcal{N}(0,{(0.05·{\left({\mathit{e}}_{\mathrm{true}}\right)}_{i,j})}^2),i=1,\dots ,Nj=1,\dots ,N_l.\end{array}$$

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

6.1. Run Time Comparison

In the following, Algortihm 1(1.1) $$\begin{array}{c}\hfill {\int }_0^{\infty }k(r,l)n(r)\mathrm{d}r=e(l)\text{with}e(l)=-\frac{log(I_{\mathrm{long}}(l))-log(I_{\mathrm{short}}(l))}{L_{\mathrm{long}}-L_{\mathrm{short}}},\end{array}$$(1.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 ${\mathit{q}}_{N_{tr}}(\mathit{x})$ we directly took the trunction index from (6.1(6.1) $$\begin{array}{cc}\hfill \mathit{\rho }& =2\mathit{\pi }\frac{r}{l},\hfill \\ \hfill M& =max\{|\mathit{\rho }|,|\mathit{\rho }·m_{\mathrm{part}}(l)|,|\mathit{\rho }·m_{\mathrm{med}}(l)|\},\hfill \\ \hfill N_{tr}& :=\lceil |M+4.05·M^{\frac{1}{3}}+2|\rceil \hfill \end{array}$$(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.1.1. Results for Ag

6.1.2. Results for CsI

6.1.3. Results for ${\mathrm{H}}_2$O

6.2. Maximal relative deviations

For the 10 simulation runs, we list the maximal relative deviations$$\begin{array}{c}\hfill 100·\frac{\Vert {\left(n_{\mathrm{part}}^1(l_j),k_{\mathrm{part}}^1(l_j)\right)}^T-{\left(n_{\mathrm{part}}^2(l_j),k_{\mathrm{part}}^2(l_j)\right)}^T{\Vert }_2}{\Vert {\left(n_{\mathrm{part}}^1(l_i),k_{\mathrm{part}}^1(l_j)\right)}^T{\Vert }_2}\end{array}$$

of the refractive index reconstructions ${\left(n_{\mathrm{part}}^2(l_j),k_{\mathrm{part}}^2(l_j)\right)}^T$ from method 2 from ${\left(n_{\mathrm{part}}^1(l_j),k_{\mathrm{part}}^1(l_j)\right)}^T$ of method 1 for $j=1,\dots ,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.2.1. Results for Ag

6.2.2. Results for CsI

6.2.3. Results for ${\mathrm{H}}_2$O

6.3. 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(6.1) $$\begin{array}{cc}\hfill \mathit{\rho }& =2\mathit{\pi }\frac{r}{l},\hfill \\ \hfill M& =max\{|\mathit{\rho }|,|\mathit{\rho }·m_{\mathrm{part}}(l)|,|\mathit{\rho }·m_{\mathrm{med}}(l)|\},\hfill \\ \hfill N_{tr}& :=\lceil |M+4.05·M^{\frac{1}{3}}+2|\rceil \hfill \end{array}$$(6.1) ) all the time.

7. Nonlinear Tikhonov regularization

So far we have solved the regression problem (4.3(4.3) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{t,\mathit{\delta }}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega },\hfill \\ \hfill \text{with}{F}_{t,\mathit{\delta }}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2\hfill \end{array}$$(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(7.1) $$\begin{array}{c}\hfill {\mathit{x}}_{\mathit{\gamma }}:=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2+\mathit{\gamma }{\Vert \mathit{x}-{\mathit{x}}^{\ast }\Vert }_2^2\text{s.t.}\mathit{x}\in \mathrm{\Omega }\end{array}$$(7.1)

when we apply it on (4.3(4.3) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{t,\mathit{\delta }}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega },\hfill \\ \hfill \text{with}{F}_{t,\mathit{\delta }}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2\hfill \end{array}$$(4.3) ), cf. [9]. Here $\mathit{\gamma }$ is a regularization parameter and ${\mathit{x}}^{\ast }$ is an estimate of the sought-after true solution. In many cases, the unregularized problem has a whole set of minimizers, thus the vector ${\mathit{x}}^{\ast }$ works also as a selection criterion. Now if a reasonable ${\mathit{x}}^{\ast }$ is found, the regularization parameter $\mathit{\gamma }$ can be determined with the discrepancy principle, i.e. $\mathit{\gamma }$ is computed such that$$\begin{array}{c}\hfill \Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t({\mathit{x}}_{\mathit{\gamma }})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2=R(\mathit{\delta }),\end{array}$$

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

The problem of finding a good estimate ${\mathit{x}}^{\ast }$ still remains. In [10], an alternative implementable parameter choice strategy without the need of an ${\mathit{x}}^{\ast }$ is derived. Applied on our problem, it gives$$\begin{array}{cc}& \mathit{\gamma }〈{\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t({\mathit{x}}_{\mathit{\gamma }})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }),{\mathrm{J}}_{\mathit{\gamma }}^{-1}({\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t({\mathit{x}}_{\mathit{\gamma }})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta })〉=R(\mathit{\delta }),\hfill \\ \hfill \text{with}& {\mathrm{J}}_{\mathit{\gamma }}:=\mathit{\gamma }\mathrm{I}+{\mathbf{\Sigma }}^{-\frac{1}{2}}\mathrm{J}a{c}_{{\mathit{f}}_t}({\mathit{x}}_{\mathit{\gamma }})\mathrm{J}a{c}_{{\mathit{f}}_t}{({\mathit{x}}_{\mathit{\gamma }})}^T{\mathbf{\Sigma }}^{-\frac{1}{2}}.\hfill \end{array}$$

This method has the drawback that the matrix ${\mathrm{J}}_{\mathit{\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(7.1) $$\begin{array}{c}\hfill {\mathit{x}}_{\mathit{\gamma }}:=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2+\mathit{\gamma }{\Vert \mathit{x}-{\mathit{x}}^{\ast }\Vert }_2^2\text{s.t.}\mathit{x}\in \mathrm{\Omega }\end{array}$$(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 ${\mathit{x}}_{\mathrm{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 neighbouring 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 ${\mathit{x}}_j^i$ the jth solution found for wavelength $l_i$ for $i=1,\dots ,s$ and $j=1,\dots ,N_i$. Now we wish to find the smoothest combined solution from all possible combinations ${\mathit{x}}_{j_1}^1$, ..., ${\mathit{x}}_{j_s}^s$ for $j_i=1,\dots ,N_i$, hence we have a total number of ${\prod }_{i=1}^s{N}_i$ combinations. Here we measure smoothness of a combination ${\mathit{x}}_{j_1}^1$, ..., ${\mathit{x}}_{j_s}^s$ with the sum$$\begin{array}{c}\hfill S:=\sum _{i=2}^{s-1}\left({\left(({\mathit{x}}_{j_{i-1}}^{i-1}{)}_1-2({\mathit{x}}_{j_i}^i{)}_1+({\mathit{x}}_{j_{i+1}}^{i+1}{)}_1\right)}^2+{\left(({\mathit{x}}_{j_{i-1}}^{i-1}{)}_2-2({\mathit{x}}_{j_i}^i{)}_2+({\mathit{x}}_{j_{i+1}}^{i+1}{)}_2\right)}^2\right)\end{array}$$

of its second finite differences both in the real parts $({\mathit{x}}_{j_i}^i{)}_1$ and its imaginary parts $({\mathit{x}}_{j_i}^i{)}_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.

The main loop spanning over the lines 5–52 iterates through all positions $z=2,\dots ,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 ${\mathit{x}}_{c1}^{z-1}$, ${\mathit{x}}_{c2}^z$ and ${\mathit{x}}_{c3}^{z+1}$ for $c1=1,\dots ,N_{z-1}$, $c2=1,\dots ,N_z$ and $c3=1,\dots ,N_{z+1}$, i.e. they loop through all of their middle points and left and right neighbours at position z. In lines 9–12, the variable $S_{\mathrm{cur}}$ is initialized with the sum of the squared second finite differences in the real and imaginary parts of the current vectors ${\mathit{x}}_{c1}^{z-1}$, ${\mathit{x}}_{c2}^z$ and ${\mathit{x}}_{c3}^{z+1}$ and the positions $z-1$, z and $z+1$ of the array Comb are filled with the current vectors. For $z\ge 3$, the loop in lines 13–28 successively fills the positions $k=z-2,z-3,\dots ,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 points are fixed and taken as the leftmost two vectors from the array Comb and the right point runs through all ${\mathit{x}}_j^k$ for $j=1,\dots ,N_k$. After the vector ${\mathit{x}}_{min}$ giving out of all of the ${\mathit{x}}_j^k$’s giving the minimal sum $D_{min}$ is found, the $D_{min}$ is added to $S_{\mathrm{cur}}$ and ${\mathit{x}}_{min}$ is stored in kth entry Comb(k). In a similar way, the loop in lines 29–44 succesively fills the positions $k=z+2,z+3,\dots ,s$ for $z\le s-2$. This time the left and middle points are fixed and taken as the rightmost two points of the array Comb, whereas the left point iterates through all ${\mathit{x}}_j^k$ for $j=1,\dots ,N_k$. Again the vector ${\mathit{x}}_{min}$ is that one of the ${\mathit{x}}_j^k$’s giving the minimal sum $D_{min}$ and it is stored in Comb(k). As well the sum $D_{min}$ is added to $S_{\mathrm{cur}}$.

In the above procedure, every triple of neighbouring 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_{\mathrm{cur}}$’s is selected in lines 45–48 to be the final output SmoothestCombination.

Define $N_{\mathrm{total}}:={\sum }_{j=1}^s{N}_j$. Then the total number of operations needed for Algorithm 2 can be estimated by $\mathcal{O}(N_{\mathrm{total}}{\sum }_{j=2}^{s-1}{N}_{j-1}{N}_j{N}_{j+1})$ which is considerably less than the $\mathcal{O}({\prod }_{j=1}^s{N}_j)$ 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 ${\mathit{x}}^1$, ..., ${\mathit{x}}^s$ denote a set of solution for any of the five optical windows. Then the joint posterior probability density of ${\mathit{x}}^1$, ..., ${\mathit{x}}^s$ is given by(9.1) $$\begin{array}{cc}\hfill p({\mathit{x}}^1,\dots ,{\mathit{x}}^s|{\mathit{e}}^1,\dots ,{\mathit{e}}^s)& =\prod _{j=1}^{s}p({\mathit{x}}^j|{\mathit{e}}^j)\propto exp\left(-\frac{1}{2}{\sum }_{j=1}^s{\Vert {\mathbf{\Sigma }}_j^{-\frac{1}{2}}({\mathit{f}}_t^j({\mathit{x}}^j)-{\mathit{e}}^j)\Vert }_2^2\right)\hfill \\ \hfill & \times \prod _{j=1}^s{I}_{\mathrm{\Omega }}({\mathit{x}}^j),\hfill \end{array}$$(9.1) where ${\mathit{e}}^j$ is the data vector for the jth wavelength $l_j$ having N entries with N being the size of the radius grid. Moreover, ${\mathbf{\Sigma }}_j$ is the scaled covariance matrix for $l_j$ and ${\mathit{f}}_t^j(\mathit{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 ${\mathit{f}}_{t_j}(\mathit{x})$ are decaying fast for each $t_j$, it is convenient to change to the same truncation index $t:=max\{t_1,\dots ,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(9.1) $$\begin{array}{cc}\hfill p({\mathit{x}}^1,\dots ,{\mathit{x}}^s|{\mathit{e}}^1,\dots ,{\mathit{e}}^s)& =\prod _{j=1}^{s}p({\mathit{x}}^j|{\mathit{e}}^j)\propto exp\left(-\frac{1}{2}{\sum }_{j=1}^s{\Vert {\mathbf{\Sigma }}_j^{-\frac{1}{2}}({\mathit{f}}_t^j({\mathit{x}}^j)-{\mathit{e}}^j)\Vert }_2^2\right)\hfill \\ \hfill & \times \prod _{j=1}^s{I}_{\mathrm{\Omega }}({\mathit{x}}^j),\hfill \end{array}$$(9.1) ) is equivalent to maximize all single densities $p({\mathit{x}}^j|{\mathit{e}}^j)$ independently, i.e. a joint MAP-estimator$$\begin{array}{c}\hfill {\mathit{x}}_{opt}^1,...,{\mathit{x}}_{opt}^s\in X_{opt}:=\underset{{\mathit{x}}^1,...,{\mathit{x}}^s}{\mathrm{a}rgmax}p({\mathit{x}}^1,...,{\mathit{x}}^s|{\mathit{e}}^1,...,{\mathit{e}}^s)\end{array}$$

consists of the single MAP-estimators$$\begin{array}{c}\hfill {\mathit{x}}_{opt}^j\in X_{opt}^j:=\underset{\mathit{x}}{\mathrm{a}rgmax}p(\mathit{x}|{\mathit{e}}^j)\end{array}$$

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

This behaviour changes when we replace the joint prior probabilty density$$\begin{array}{c}\hfill p_{\mathrm{prior}}({\mathit{x}}^1,\dots ,{\mathit{x}}^s)={\left(\mathrm{v}ol(\mathrm{\Omega })\right)}^{-s}\prod _{j=1}^s{I}_{\mathrm{\Omega }}({\mathit{x}}^j)\end{array}$$

with$$\begin{array}{c}\hfill p_{\mathrm{prior}}({\mathit{x}}^1,\dots ,{\mathit{x}}^s)\propto exp\left(-\frac{1}{2}\mathit{\gamma }S({\mathit{x}}^1,\dots ,{\mathit{x}}^s)\right)\prod _{j=1}^s{I}_{\mathrm{\Omega }}({\mathit{x}}^j),\end{array}$$

where$$\begin{array}{cc}\hfill S({\mathit{x}}^1,\dots ,{\mathit{x}}^s)& :=\sum _{i=2}^{s-1}\left((({\mathit{x}}^{i-1}{)}_1-2({\mathit{x}}^i{)}_1+({\mathit{x}}^{i+1}{)}_1{)}^2+(({\mathit{x}}^{i-1}{)}_2-2({\mathit{x}}^i{)}_2\right.\hfill \\ \hfill & \left.+({\mathit{x}}^{i+1}{)}_2{)}^2\right)+\mathit{\rho }\sum _{i=1}^s(({\mathit{x}}^i{)}_1^2+({\mathit{x}}^i{)}_2^2),\hfill \end{array}$$

where $\mathit{\gamma }$ is a regularization parameter and $\mathit{\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 $\mathit{\rho }={10}^{-8}$, such that the resulting regularization operator gets regular. This means that the regularized regression problem (9.2(9.2) $$\begin{array}{c}\hfill {\mathit{x}}_{opt}^1,...,{\mathit{x}}_{opt}^s\in X_{opt}:=\underset{{\mathit{x}}^1,...,{\mathit{x}}^s}{\mathrm{a}rgmax}p({\mathit{x}}^1,...,{\mathit{x}}^s|{\mathit{e}}^1,...,{\mathit{e}}^s)\end{array}$$(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 neighbouring points; therefore, a decoupled computation of the joint MAP-estimator(9.2) $$\begin{array}{c}\hfill {\mathit{x}}_{opt}^1,...,{\mathit{x}}_{opt}^s\in X_{opt}:=\underset{{\mathit{x}}^1,...,{\mathit{x}}^s}{\mathrm{a}rgmax}p({\mathit{x}}^1,...,{\mathit{x}}^s|{\mathit{e}}^1,...,{\mathit{e}}^s)\end{array}$$(9.2)

with$$\begin{array}{cc}\hfill p({\mathit{x}}^1,\dots ,{\mathit{x}}^s|{\mathit{e}}^1,\dots ,{\mathit{e}}^s)& \propto exp\left(-\frac{1}{2}{\sum }_{j=1}^s{\Vert {\mathbf{\Sigma }}_j^{-\frac{1}{2}}({\mathit{f}}_t^j({\mathit{x}}^j)-{\mathit{e}}^j)\Vert }_2^2-\frac{1}{2}\mathit{\gamma }S({\mathit{x}}^1,\dots ,{\mathit{x}}^s)\right)\hfill \\ \hfill & \times \prod _{j=1}^s{I}_{\mathrm{\Omega }}({\mathit{x}}^j)\hfill \end{array}$$

for each wavelength seperately is not possible anymore after changing to the new prior density. However the result vectors ${\mathit{x}}^1$, ..., ${\mathit{x}}^s$ from Algorithm 2 form a good start vector to solve the nonlinear regression problem (9.2(9.2) $$\begin{array}{c}\hfill {\mathit{x}}_{opt}^1,...,{\mathit{x}}_{opt}^s\in X_{opt}:=\underset{{\mathit{x}}^1,...,{\mathit{x}}^s}{\mathrm{a}rgmax}p({\mathit{x}}^1,...,{\mathit{x}}^s|{\mathit{e}}^1,...,{\mathit{e}}^s)\end{array}$$(9.2) ).

We selected the regularization parameter $\mathit{\gamma }$ using the discrepancy principle, i.e. we compute $\mathit{\gamma }$ such that the regularized solution$$\begin{array}{cc}\hfill {\mathit{x}}_{\mathit{\gamma }}^1,...,{\mathit{x}}_{\mathit{\gamma }}^s\in X_{\mathit{\gamma }}:& =\underset{{\mathit{x}}^1,...,{\mathit{x}}^s}{\mathrm{a}rgmin}\sum _{j=1}^s{\Vert {\mathbf{\Sigma }}_j^{-\frac{1}{2}}({\mathit{f}}_t^j({\mathit{x}}^j)-{\mathit{e}}^j)\Vert }_2^2+\mathit{\gamma }S({\mathit{x}}^1,...,{\mathit{x}}^s)\hfill \\ \hfill & \text{s.t.}{\mathit{x}}^j\in \mathrm{\Omega },j=1,...,s\hfill \end{array}$$

fulfils a relation of the form$$\begin{array}{c}\hfill \sum _{j=1}^s{\Vert {\mathbf{\Sigma }}_j^{-\frac{1}{2}}({\mathit{f}}_t^j({\mathit{x}}_{\mathit{\gamma }}^j)-{\mathit{e}}^j)\Vert }_2^2=R(\mathit{\delta }),\end{array}$$

where $R(\mathit{\delta })$ is a proposed residual value depending on the noise level $\mathit{\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 behaviour of our model ${\mathit{f}}_t(\mathit{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({\mathit{e}}^j|{\mathit{x}}^j)$ for $j=1,\dots ,s$ is Gaussian, the joint observed density $p({\mathit{e}}^1,\dots ,{\mathit{e}}^1|{\mathit{x}}^1,\dots ,{\mathit{x}}^s)={\prod }_{j=1}^{s}p({\mathit{e}}^j|{\mathit{x}}^j)$ is Gaussian as well. We have ${\mathit{x}}^j\in {\mathbb{R}}^2$, thus the sum of residuals ${\sum }_{j=1}^s{\Vert {\mathbf{\Sigma }}_j^{-\frac{1}{2}}({\mathit{f}}_t^j({\mathit{x}}^j)-{\mathit{e}}^j)\Vert }_2^2$ running through all wavelengths in the optical window is ${\mathit{\chi }}^2(2s)$-distributed. This yields$$\begin{array}{c}\hfill \mathbb{E}(\sum _{j=1}^s{\Vert {\mathbf{\Sigma }}_j^{-\frac{1}{2}}({\mathit{f}}_t^j({\mathit{x}}^j)-{\mathit{e}}^j)\Vert }_2^2)=2s.\end{array}$$

Now a widely proposed residual value for the discrepancy principle is $\mathit{\tau }·2s$, where $\mathit{\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\mathit{\tau }s$. Therefore, we proposed a residual value which depends more dynamically on the observed behoviour of the residual.

Let ${\mathit{x}}_0^1$, ..., ${\mathit{x}}_0^s$ denote the unregularized solutions, i.e. the results of Algorithm 2. Then their squared residual is given by $R_0:={\sum }_{j=1}^s{\Vert {\mathbf{\Sigma }}_j^{-\frac{1}{2}}({\mathit{f}}_t^j({\mathit{x}}_0^j)-{\mathit{e}}^j)\Vert }_2^2$. We first proposed$$\begin{array}{c}\hfill R(\mathit{\delta })=max\left\{2{\mathit{\tau }}_{1}s,\mathit{\tau }{R}_0\right\},\end{array}$$

where we selected ${\mathit{\tau }}_1=1.1$. This means that the residual of the regularized solution is beginning at $R_0$ at least increased by the factor ${\mathit{\tau }}_1$, which avoids underregularization. If it then happens that $\frac{R(\mathit{\delta })}{R_0}>\mathit{\theta }$ with $\mathit{\theta }=1.5$ the proposed residual is most likely too big and overregularization occurs. In this case, we corrected $R(\mathit{\delta })$ by setting$$\begin{array}{c}\hfill R(\mathit{\delta })=max\left\{2{\mathit{\tau }}_{2}s,\mathit{\theta }{R}_0\right\}\end{array}$$

with ${\mathit{\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\mathrm{\mu }$m, $r_2:=0.2\mathrm{\mu }$m and $r_3:=0.3\mathrm{\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\mathrm{\mu }$m, $0.2305\mathrm{\mu }$m and $0.11\mathrm{\mu }$m carrying the most information content found in this study are roughly comparable to our radii.

We computed original spectral extinctions$$\begin{array}{c}\hfill {\left({\mathit{e}}_{\mathrm{true}}\right)}_{i,j}:=\mathit{\pi }{r}_i^2\sum _{n=1}^{N_{tr}}{q}_n(m_{\mathrm{med}}(l_i),m_{\mathrm{part}}(l_i),r_j,l_i),i=1,\dots ,48,j=1,\dots ,3\end{array}$$

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$$\begin{array}{c}\hfill {\left(\mathit{e}\right)}_{i,j}={\left({\mathit{e}}_{\mathrm{true}}\right)}_{i,j}+{\mathit{\delta }}_{i,j}\text{with}{\mathit{\delta }}_{i,j}\sim \mathcal{N}(0,{(0.05·{\left({\mathit{e}}_{\mathrm{true}}\right)}_{i,j})}^2),i=1,\dots ,48,j=1,\dots ,3.\end{array}$$

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({\mathit{e}}_{\mathrm{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(4.3) $$\begin{array}{cc}\hfill X_{t,\mathit{\delta }}& :=\underset{\mathit{x}\in {\mathbb{R}}^D}{\mathrm{a}rgmin}{F}_{t,\mathit{\delta }}(\mathit{x})\text{s.t.}\mathit{x}\in \mathrm{\Omega },\hfill \\ \hfill \text{with}{F}_{t,\mathit{\delta }}(\mathit{x})& :=\Vert {\mathbf{\Sigma }}^{-\frac{1}{2}}({\mathit{f}}_t(\mathit{x})-\mathit{f}({\mathit{x}}_{\mathrm{true}})-\mathit{\delta }){\Vert }_2^2\hfill \end{array}$$(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(1.1) $$\begin{array}{c}\hfill {\int }_0^{\infty }k(r,l)n(r)\mathrm{d}r=e(l)\text{with}e(l)=-\frac{log(I_{\mathrm{long}}(l))-log(I_{\mathrm{short}}(l))}{L_{\mathrm{long}}-L_{\mathrm{short}}},\end{array}$$(1.1) and 2(1.2) $$\begin{array}{c}\hfill n\mathit{\pi }{r}_m^2{Q}_{\mathrm{ext}}(m_{\mathrm{med}}(l),m_{\mathrm{part}}(l),r_m,l)=e(l),\end{array}$$(1.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(1.1) $$\begin{array}{c}\hfill {\int }_0^{\infty }k(r,l)n(r)\mathrm{d}r=e(l)\text{with}e(l)=-\frac{log(I_{\mathrm{long}}(l))-log(I_{\mathrm{short}}(l))}{L_{\mathrm{long}}-L_{\mathrm{short}}},\end{array}$$(1.1) and 2(1.2) $$\begin{array}{c}\hfill n\mathit{\pi }{r}_m^2{Q}_{\mathrm{ext}}(m_{\mathrm{med}}(l),m_{\mathrm{part}}(l),r_m,l)=e(l),\end{array}$$(1.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(1.1) $$\begin{array}{c}\hfill {\int }_0^{\infty }k(r,l)n(r)\mathrm{d}r=e(l)\text{with}e(l)=-\frac{log(I_{\mathrm{long}}(l))-log(I_{\mathrm{short}}(l))}{L_{\mathrm{long}}-L_{\mathrm{short}}},\end{array}$$(1.1) and 2(1.2) $$\begin{array}{c}\hfill n\mathit{\pi }{r}_m^2{Q}_{\mathrm{ext}}(m_{\mathrm{med}}(l),m_{\mathrm{part}}(l),r_m,l)=e(l),\end{array}$$(1.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 s in the average case and can rise up to 200 s 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 s and are always below 95 s. 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 to $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 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 behaviour can be observed for the run times of Algorithm 1. For the first four optical windows, they are between 20 and 45 s on average and below 100 s in the extreme cases, whereas for the last optical window they are between 30 and 170 s on average and can even rise up to 350 s. 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$$\begin{array}{c}\hfill {\left({\mathit{e}}_{\mathrm{true}}\right)}_{i,j}:=\mathit{\pi }{r}_i^2\sum _{n=1}^{N_{tr}}{q}_n(m_{\mathrm{med}}(l_i),m_{\mathrm{part}}(l_i),r_j,l_i),i=1,\dots ,48,j=1,\dots ,3\end{array}$$

for all wavelengths and added zero-mean Gaussian noise to it in order to obtain the simulated noisy spectral extinctions$$\begin{array}{c}\hfill {\left(\mathit{e}\right)}_{i,j}={\left({\mathit{e}}_{\mathrm{true}}\right)}_{i,j}+{\mathit{\delta }}_{i,j}\text{with}{\mathit{\delta }}_{i,j}\sim \mathcal{N}(0,{(0.15·{\left({\mathit{e}}_{\mathrm{true}}\right)}_{i,j})}^2),i=1,\dots ,48,j=1,\dots ,3.\end{array}$$

for the first study and$$\begin{array}{c}\hfill {\left(\mathit{e}\right)}_{i,j}={\left({\mathit{e}}_{\mathrm{true}}\right)}_{i,j}+{\mathit{\delta }}_{i,j}\text{with}{\mathit{\delta }}_{i,j}\sim \mathcal{N}(0,{(0.3·{\left({\mathit{e}}_{\mathrm{true}}\right)}_{i,j})}^2),i=1,\dots ,48,j=1,\dots ,3.\end{array}$$

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({\mathit{e}}_{\mathrm{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 studies, 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.1. Noise-free refractive indices

Table

	Average $L^2$-errors (%)
Original ${\mathrm{H}}_2$O	Log-normal	RRSB	Hedrih
volume per cent	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

Table

	Average fraction deviation (%)
Original ${\mathrm{H}}_2$O	Log-normal	RRSB	Hedrih
volume per cent	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.2. Noisy refractive indices

Table

	Average $L^2$-errors (%)
Original ${\mathrm{H}}_2$O	Log-normal	RRSB	Hedrih
volume per cent	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

Table

	Average fraction deviation (%)
Original ${\mathrm{H}}_2$O	Log-normal	RRSB	Hedrih
volume per cent	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.3. Results for mixtures of ${\mathrm{H}}_2$O and Ag

13.3.1. Noise-free refractive indices

Table

	Average $L^2$-errors (%)
Original ${\mathrm{H}}_2$O	Log-normal	RRSB	Hedrih
volume per cent	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

Table

	Average fraction deviation (%)
Original ${\mathrm{H}}_2$O	Log-normal	RRSB	Hedrih
volume per cent	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.3.2. Noisy refractive indices

Table

	Average $L^2$-errors (%)
Original ${\mathrm{H}}_2$O	Log-normal	RRSB	Hedrih
volume per cent	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

Table

	Average fraction deviation (%)
Original ${\mathrm{H}}_2$O	Log-normal	RRSB	Hedrih
volume per cent	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.4. 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 studies. 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\mathrm{\mu }$m, $0.2\mathrm{\mu }$m and $0.3\mathrm{\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.

Acknowledgements

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

Additional information

Funding

This work is sponsored by the German Federal Ministry of Education and Research (BMBF) [contract number 02NUK022A0]. Responsibility for the content of this report lies with the authors.

Notes

No potential conflict of interest was reported by the author.