2.1.

Stochastic-Gradient Descent

In a stochastic setting, for instance, when our forward model $A$ is a Monte Carlo model of radiative transport, then the true cost $F$ and gradient $\nabla F$ given in Eq. (3) are not directly available. Instead, we may only have access to estimates of the cost function and gradient (provided by a stochastic first-order oracle). In Sec. 3, we detail the nature of these stochastic Monte Carlo computations in the radiative transport setting. In the interest of generality, for now, we simply assume such models exist, and we can make a call to a “stochastic oracle” to attain $F_{S_n}$ and $\nabla F_{S_n}$ which we assume are nonbiased approximations., i.e.,

In this study, if we are to accept a level of variance and imperfection in our forward/adjoint models, this of course raises the question of how much variance is acceptable in order for SGD to be successful? Furthermore, what measure of the variance is the best indicator in terms of efficiency/performance for common Monte Carlo solvers? To begin to answer this, it is important to first note that fixed-step SGD does not in general converge to a solution.^28,29 That is, if ${\alpha }_n$ is fixed for all $n$, eventually there will come a point where the next update of the estimate [with the term ${\alpha }_n\nabla F_{S_n}(x_{n-1})$] will reliably “undo” the work of the prior step, which will effectively halt the descent. The point at which this occurs depends on the variance of $\nabla F_{S_n}$. We can see this by rewriting the sampled stochastic gradient estimate as

To prevent this from happening, we may take one of two actions (or a combination thereof): (i) reduce the step size at each iteration such that we can avoid “backtracking” in the descent, more on this in Sec. 2.3 or (ii) gradually improve the accuracy of our sampled gradient such that the variance of the sampled gradient remains below some threshold value compared to the norm of the true gradient $\nabla F$. In other words, we may wish to ensure the inequality

2.2.

Adaptive Sample Size

We have discussed two different measures of the variance in the sampled gradient $\nabla F_{S_n}$ that we wish to investigate in the context of Monte Carlo estimation of media properties, viz., the norm test [Eq. (7)] and the inner product test [Eq. (8)]. To satisfy the inequalities defining these tests as the GD progresses, we will be required to reduce the variance in the sampled gradients $\nabla F_{S_n}$ whenever the norm test or inner product test fail. This can be done by increasing the sample size (number of photons used $|S_n|$) when making a call to the stochastic oracle. Two practical considerations are still required: first, how to compute the “true” gradient $\nabla F$, which is needed to evaluate the norm test and inner product test; and second, by how much we should increase the sample size in a situation where one of the tests fails?

The true gradient $\nabla F$ is only calculable in the limit that an infinite number of photons are used in the Monte Carlo model. This limit can equivalently be represented as an average over independent repeated outputs of the sampled gradient

In terms of increasing the sample size in the event where the inner product and/or norm tests fail, this can be done in a number of ways. A simple method we will employ in this study is to scale the current sample size by some factor $\kappa (n)$ to increase the number of photons used in the next iteration

2.3.

Step Size

In cases where we are not taking actions to bound the error in the sampled gradient (such as enforcing successful outcomes of an inner product test or norm test), fixed-step SGD may only converge to a region around the solution. Reducing the step size sufficiently at each step is usually required to allow convergence.³¹ However, it can be shown that if we are bounding the error in the sampled gradient, e.g., by increasing the sample size, then fixed-step SGD may converge so long as the following is satisfied for all $n$:³⁰

In this study, we aim to satisfy the aforementioned step size criteria for an assumed value of the Lipschitz constant $L$, which we will choose conservatively depending on the particular scenario. However, as we are primarily interested in reaching the best possible solution for a given allocation of computational resources, convergence to a region around the unique solution may be sufficient for our purposes. For this reason, we will also investigate larger step size criteria, which violate Eq. (11), yet exhibit good performance in our scenarios of interest.

Taking the above considerations into account, we present a basic method for SGD using adaptive sample sizes in Algorithm 1 (simplified from Ref. 30). The algorithm imposes a limit on the total number of photons to be simulated using Monte Carlo transport models throughout the entire descent, $N_{\mathrm{ph}}$.

Algorithm 1

Inversion using Monte Carlo sampled gradients with adaptive sample size.

3.1.

Forward Model

For any optical source $Q(\mathbf{r},\widehat{\mathbf{s}})$, either incident on a medium or present within it, we wish to model the resulting radiance, $\varphi (\mathbf{r},\widehat{\mathbf{s}})$, describing the radiant flux at each position $\mathbf{r}$, and in each direction $\widehat{\mathbf{s}}$. This can be achieved using the RTE.

3.2.

Gradient Computation: Adjoint Model

To compute the gradient of our cost function $\nabla F$ with respect to the optical properties of the medium, we make use of an adjoint RTE model. Although direct methods of finding the derivative of a Monte Carlo method can also be developed,¹² adjoint methods have more applicability in general, and also allow closer comparison with optimization techniques used in machine learning. For further details of forward and adjoint methods in the RTE, we refer to Ref. 36; for specific details of CPI problems, we refer to Ref. 37. We first consider a change to Eq. (12) where ${\mu }_{\mathrm{a}}\to {\mu }_{\mathrm{a}}+{\mu }_{\mathrm{a}}^{\delta }$, ${\mu }_{\mathrm{s}}\to {\mu }_{\mathrm{s}}+{\mu }_{\mathrm{s}}^{\delta }$, for the same source $Q$, which results in a change in radiance $\varphi \to \varphi +{\varphi }^{\delta }$. This implies

3.2.1.

QPAT case

In QPAT, the medium is illuminated with a pulsed optical source, $Q$ (see Fig. 1). The distributed optical energy is absorbed at various points within the sample, giving rise to internal acoustic waves. These acoustic waves can be detected at the surface of the medium by a sensor and processed to locate the initial pressure distribution $p_0$ within the medium.^38–40 This internal pressure distribution is related to the spatial distribution of absorbed optical energy, $h$, where

Eq. (17)

$$h(\mathbf{r})={\mu }_{\mathrm{a}}(\mathbf{r})\mathrm{\Phi }(\mathbf{r}),$$

and where $\mathrm{\Phi }$ is the optical fluence of Eq. (16). We have omitted the Grüneisen parameter for clarity of exposition, though this parameter can be included in practice. Assuming that we can recover the absorbed optical energy, $h$, the problem remains to find the distribution of ${\mu }_{\mathrm{a}}(\mathbf{r})$ within the medium.^41,42 It is worth noting that although the optical source is pulsed, it is acceptable to use a continuous-wave (time-independent) model to describe $\varphi$ and $\mathrm{\Phi }$ because the time scale of the acoustic wave propagation is orders of magnitude slower than the optical propagation.⁴³ First, restating our cost function in terms of the QPAT data function, $h$, we have

Eq. (18)

$$F^{\mathrm{QPAT}}=\frac{1}{2}{\int }_{\mathrm{\Omega }}{(h^{\mathrm{obs}}-h)}^2\mathrm{d}r=\frac{1}{2}{⟨h^{\mathrm{obs}}-h,h^{\mathrm{obs}}-h⟩}_{L^2(\mathrm{\Omega })}.$$

Fig. 1

Setup for QPAT.

We then write the Fréchet derivative of $F^{\mathrm{QPAT}}$ as

Eq. (19)

$$D{F}^{\mathrm{QPAT}}=-{⟨h^{\mathrm{obs}}-h,Dh{\mu }_{\mathrm{a}}^{\delta }⟩}_{L^2(\mathrm{\Omega })},$$

where ${\mu }_{\mathrm{a}}^{\delta }$ is a small change in absorption. In this paper, we will neglect changes in scattering, however, the below formalism is still general for the gradient with respect to absorption. The gradient term with respect to scattering coefficient is described in Ref. 42 and will be included in future investigations. Writing the Fréchet derivative of $h$ as

Eq. (20)

$$Dh=\mathrm{\Phi }+{\mu }_{\mathrm{a}}·D\mathrm{\Phi },$$

and defining ${\mathrm{\Phi }}^{\delta }=D\mathrm{\Phi }{\mu }_{\mathrm{a}}^{\delta }$, we arrive at

Eq. (21)

$$D{F}^{\mathrm{QPAT}}=-{⟨\mathrm{\Phi }(h^{\mathrm{obs}}-h),{\mu }_{\mathrm{a}}^{\delta }⟩}_{L^2(\mathrm{\Omega })}-{⟨{\mu }_{\mathrm{a}}(h^{\mathrm{obs}}-h),{\mathrm{\Phi }}^{\delta }⟩}_{L^2(\mathrm{\Omega })}.$$

Next, we define the adjoint radiance, ${\varphi }^{*}$, as the solution to

Eq. (22)

$${\mathcal{L}}^{*}{\varphi }^{*}={\mu }_{\mathrm{a}}(h^{\mathrm{obs}}-h)$$

where the right-hand side describes the “adjoint source” which is isotropic in $\widehat{\mathbf{s}}$. We then substitute the above into Eq. (21) to give

Eq. (23)

$$D{F}^{\mathrm{QPAT}}=-{⟨\mathrm{\Phi }(h^{\mathrm{obs}}-h),{\mu }_{\mathrm{a}}^{\delta }⟩}_{L^2(\mathrm{\Omega })}-{⟨{\mathcal{L}}^{*}{\varphi }^{*},{\varphi }^{\delta }⟩}_{L^2(\mathrm{\Omega }\times S^{n-1})},$$

where we exploited the fact that the right-hand side of Eq. (22) does not depend on direction. Using the definition of the adjoint operator, and the fact that the change in radiance is zero on the boundary $\partial \mathrm{\Omega }$ yields

Eq. (24)

$$D{F}^{\mathrm{QPAT}}=-{⟨\mathrm{\Phi }(h^{\mathrm{obs}}-h),{\mu }_{\mathrm{a}}^{\delta }⟩}_{L^2(\mathrm{\Omega })}-{⟨{\varphi }^{*},\mathcal{L}{\varphi }^{\delta }⟩}_{L^2(\mathrm{\Omega }\times S^{n-1})}.$$

Finally, we make use of the perturbation expression Eq. (15), while again, here, we neglect any change in scattering. This gives

Eq. (25)

$$D{F}^{\mathrm{QPAT}}=-{⟨\mathrm{\Phi }(h^{\mathrm{obs}}-h),{\mu }_{\mathrm{a}}^{\delta }⟩}_{L^2(\mathrm{\Omega })}+{⟨{\varphi }^{*}\varphi ,{\mu }_{\mathrm{a}}^{\delta }⟩}_{L^2(\mathrm{\Omega }\times S^{n-1})},$$

allowing us to define the (absorption) gradient as in Eq. (33) of Ref. 42

Eq. (26)

$$\frac{\partial F^{\mathrm{QPAT}}}{\partial {\mu }_{\mathrm{a}}}=\nabla F^{\mathrm{QPAT}}=-\mathrm{\Phi }(h^{\mathrm{obs}}-h)+{\int }_{S^{n-1}}{\varphi }^{*}\varphi \mathrm{d}\widehat{\mathbf{s}}$$

To compute a stochastic approximation of this gradient, we can thus use the forward model Monte Carlo solver ${\mathcal{L}}_{\mathrm{MC}}^{-1}$ to provide estimates of $\varphi$ and $\mathrm{\Phi }$, and an adjoint Monte Carlo solver ${\mathcal{L}}_{\mathrm{MC}}^{-1*}$ to produce ${\varphi }^{*}$ from an adjoint source term $Q_{\mathrm{adj}}={\mu }_{\mathrm{a}}(h^{\mathrm{obs}}-h)$, as defined in Eq. (22). Due to the symmetry of the problem, the adjoint solver is identical to the forward solver and follows the same basic operating principles. The only difference is that here the adjoint source $Q_{\mathrm{adj}}={\mu }_{\mathrm{a}}(h^{\mathrm{obs}}-h)$ may in fact be negative in some locations. This is handled by splitting the source term into two parts, one purely positive, $Q_{\mathrm{adj}}^{+}$, and one purely negative, $Q_{\mathrm{adj}}^{-}$. Two simulations are then run (where the total number of photons to be used is split between the two simulations accordingly), and the results summed to produce ${\varphi }^{*}$. Algorithm 2 describes the basic operation for computing a sampled gradient, $\nabla F_{S_n}$, for QPAT using the above derivation. This will be used in conjunction with Algorithm 1 to conduct an inversion with adaptive sample size for each iterate, $|S_n|$.

3.2.2.

UMOT case

Referring to Fig. 2, in UMOT, we have an optical light source $Q_q$ incident on a medium, as well as an optical detector $J_m$. In addition, an ultrasound source is incident on the medium, where the focus $\eta (\mathbf{r})$ is scanned through the sample.^44,45 Assuming for simplicity an ideal (delta function) ultrasound focus, the data of interest in this case are found to be of the form⁴⁶

Eq. (27)

$$b(\mathbf{r})=\eta (\mathbf{r}){\mathrm{\Phi }}_q(\mathbf{r}){\mathrm{\Phi }}_m(\mathbf{r}),$$

where ${\mathrm{\Phi }}_q$ is the fluence resulting from the optical source $Q_q$, and ${\mathrm{\Phi }}_m$ is the resulting fluence from a virtual source $Q_m$ which is reciprocal to the detector $J_m$.⁴⁶ From this point, we proceed in similar fashion as in Sec. 3.2.1, where now our data fitting error is given by

Eq. (28)

$$F^{\mathrm{UMOT}}=\frac{1}{2}{\int }_{\mathrm{\Omega }}{(b^{\mathrm{obs}}-b)}^2\mathrm{d}r=\frac{1}{2}{⟨b^{\mathrm{obs}}-b,b^{\mathrm{obs}}-b⟩}_{L^2(\mathrm{\Omega })},$$

and its Fréchet derivative as

Eq. (29)

$$D{F}^{\mathrm{UMOT}}=-{⟨b^{\mathrm{obs}}-b,Db{\mu }_{\mathrm{a}}^{\delta }⟩}_{L^2(\mathrm{\Omega })}.$$

In this case, the Fréchet derivative of $b$ becomes

Eq. (30)

$$Db=\eta {\mathrm{\Phi }}_q·D{\mathrm{\Phi }}_m+\eta {\mathrm{\Phi }}_m·D{\mathrm{\Phi }}_q,$$

leading to

Eq. (31)

$$D{F}^{\mathrm{UMOT}}=-{⟨\eta {\mathrm{\Phi }}_q(b^{\mathrm{obs}}-b),{\mathrm{\Phi }}_m^{\delta }⟩}_{L^2(\mathrm{\Omega })}-{⟨\eta {\mathrm{\Phi }}_m(b^{\mathrm{obs}}-b),{\mathrm{\Phi }}_q^{\delta }⟩}_{L^2(\mathrm{\Omega })}.$$

Here, we need to define two adjoint radiances, ${\varphi }^{*,1}$ and ${\varphi }^{*,2}$, as the solution to

Eq. (32)

$${\mathcal{L}}^{*}{\varphi }^{*,1}=\eta {\mathrm{\Phi }}_q(b^{\mathrm{obs}}-b),$$

Eq. (33)

$${\mathcal{L}}^{*}{\varphi }^{*,2}=\eta {\mathrm{\Phi }}_m(b^{\mathrm{obs}}-b),$$

and substituting into Eq. (31) to give

Eq. (34)

$$D{F}^{\mathrm{UMOT}}=-{⟨{\mathcal{L}}^{*}{\varphi }^{*,1},{\varphi }_m^{\delta }⟩}_{L^2(\mathrm{\Omega }\times S^{n-1})}-{⟨{\mathcal{L}}^{*}{\varphi }^{*,2},{\varphi }_q^{\delta }⟩}_{L^2(\mathrm{\Omega }\times S^{n-1})};$$

by the same arguments as for QPAT we get

Eq. (35)

$$D{F}^{\mathrm{UMOT}}=-{⟨{\varphi }^{*,1},\mathcal{L}{\varphi }_m^{\delta }⟩}_{L^2(\mathrm{\Omega }\times S^{n-1})}-{⟨{\varphi }^{*,2},\mathcal{L}{\varphi }_q^{\delta }⟩}_{L^2(\mathrm{\Omega }\times S^{n-1})}.$$

Again using the perturbation expression [Eq. (15)], we have

Eq. (36)

$$D{F}^{\mathrm{UMOT}}={⟨{\varphi }^{*,1}{\varphi }_m,{\mu }_{\mathrm{a}}^{\delta }⟩}_{L^2(\mathrm{\Omega }\times S^{n-1})}+{⟨{\varphi }^{*,2}{\varphi }_q,{\mu }_{\mathrm{a}}^{\delta }⟩}_{L^2(\mathrm{\Omega }\times S^{n-1})},$$

allowing us to define the (absorption) gradient as

Eq. (37)

$$\frac{\partial F^{\mathrm{UMOT}}}{\partial {\mu }_{\mathrm{a}}}=\nabla F^{\mathrm{UMOT}}={\int }_{S^{n-1}}({\varphi }^{*,1}{\varphi }_m+{\varphi }^{*,2}{\varphi }_q)\mathrm{d}\widehat{\mathbf{s}}.$$

Fig. 2

Setup for UMOT in the transmission geometry.

Thus, similar to the QPAT case, here we are able to compute a stochastic approximation of this gradient using the forward model Monte Carlo solver ${\mathcal{L}}_{\mathrm{MC}}^{-1}$ to provide ${\varphi }_q$ and ${\varphi }_m$ from our two sources, and an adjoint Monte Carlo solver ${\mathcal{L}}_{\mathrm{MC}}^{-1*}$ to produce ${\varphi }^{*,1}$ and ${\varphi }^{*,2}$ from the adjoint source terms $Q_{\mathrm{adj}}^1=\eta {\varphi }_q(b^{\mathrm{obs}}-b)$ and $Q_{\mathrm{adj}}^2=\eta {\varphi }_m(b^{\mathrm{obs}}-b)$, as defined in Eqs. (32) and (33). Here as well, adjoint source terms are split into two parts, one purely positive, $Q_{\mathrm{adj}}^{+}$, and one purely negative, $Q_{\mathrm{adj}}^{-}$, with the photon budget being split accordingly. Algorithm 3 describes the basic operation for computing a sampled gradient, $\nabla F_{S_n}$, for UMOT using the above derivation. This will be used in conjunction with Algorithm 1 to conduct an inversion with adaptive sample size for each iterate, $|S_n|$.

Algorithm 2

Monte Carlo sampled QPAT gradient.

1. Compute ${\mathcal{L}}_{\mathrm{MC}}^{-1}Q\mapsto \varphi$, $\mathrm{\Phi }$, using $|S_n|/2$ photons

2. Construct internal adjoint source $Q_{\mathrm{adj}}={\mu }_{\mathrm{a}}(h^{\mathrm{obs}}-h)$

3. Compute ${\mathcal{L}}_{\mathrm{MC}}^{-1*}{Q}_{\mathrm{adj}}\mapsto {\varphi }^{*}$, ${\mathrm{\Phi }}^{*}$, using $|S_n|/2$ photons

4. Use Eq. (26) to compute gradient $\nabla F_{S_n}$

3.3.

Fluence Monte Carlo

It should be noted that numerous Monte Carlo radiative transport solvers do not explicitly output the radiance, as this requires additional programming to store the angular ordinates at each location. Commonly, only the fluence will be available, which is the angular integral of the radiance Eq. (16). In such cases, the aforementioned integrals for the gradients of interest Eqs. (26) and (37) can be computed under the assumption of approximately angularly isotropic radiances, where for example $\int {\varphi }^{*}\varphi \mathrm{d}\widehat{\mathbf{s}}$ becomes ${\mathrm{\Phi }}^{*}\mathrm{\Phi }$. The accuracy of this approximation of course depends on the true angular dependence of the radiances, where the approximation is poorest in regions close to directional light sources, but improves further away. The higher the scattering asymmetry $g$ of the medium, the slower the approximation improves as a function of distance from these sources. In many cases, however, this is a satisfactory assumption and is employed in the below example cases.

4.1.

QPAT

We begin with our example QPAT problem. The starting sample size in all cases shown is $|S_1|=200$ photons per iteration (100 for each forward run, and 100 for each adjoint run in accordance with Algorithm 2), and the total photon budget for the inversion was set to $N_{\mathrm{ph}}=2\times {10}^6$ photons. The Lipschitz constant was set at $L=2.5$, as this displayed stable descent in our test problems using large photon budgets (low-variance case). The initial estimate of the absorption distribution in the medium is ${\mu }_{\mathrm{a}}=0.2{\mathrm{cm}}^{-1}$ in all layers. The scattering coefficient of all layers was set to ${\mu }_{\mathrm{s}}=40{\mathrm{cm}}^{-1}$, and the scattering asymmetry parameter was set to $g=0.9$. The ground-truth absorption ${\mu }_{\mathrm{a}}^{\text{true}}$ is shown in Fig. 3(a), along with the final retrieved absorption distributions obtained via strategies 1, 2, and 3 using the stated values of ${\gamma }_{\mathrm{tot}}$ and ${\gamma }_{\parallel }$. Figure 3(b) shows the corresponding measured data and the final forward modeled data for each of the strategies. Figure 4 shows the outcome of each strategy in terms of the sampled data cost function $F_{S_n}$ and the absorption error function $F_{{\mu }_{\mathrm{a}}}$. It can be seen that the ranking of these methods in terms of the lowest achieved value of the sampled cost function $F_{S_n}$ does not correlate directly to the best outcomes in terms of the error in the estimated absorption $F_{{\mu }_{\mathrm{a}}}$. This is due to the above-mentioned dependence of the sampled data cost function on the sample size used in the forward model, where for example the case of strategy 2 only appears to perform poorly in terms of $F_{S_n}$ due to its small sample size used throughout the inversion. This is more clearly shown in Fig. 3(b), where the final forward modeled data from strategy 2 are noisier than the other strategies due to the low sample size at the end of the inversion, where this noise would clearly impact the sampled cost function. It is worthy to note that the relevant step sizes and sample sizes for each of these three examples are shown in Fig. 5. Before finding the best parameter for strategy 1, we trialed a range of values of ${\gamma }_{\mathrm{tot}}$ over the range (0.1, 200). With lower values, the adaptive sample size was required to increase rapidly to maintain high-quality (low variance) sample gradients. This resulted in the photon budget being depleted early, terminating the descent after around 100 iterations, which did not perform well. Too large a value of ${\gamma }_{\parallel }$ and the norm test never failed, meaning the sample size was never required to increase and the inversion progressed for the maximum 10,000 iterations permitted by the photon budget. However, as strategy 1 has a fixed value of ${\gamma }_{\mathrm{tot}}^2$ in the denominator of the step size, large values also result in step sizes that were too small to perform well. A value of ${\gamma }_{\mathrm{tot}}=4$ was found to strike a balance between these two extremes and was the best performer using strategy 1. Strategy 2 has an adaptive step size which selects the largest possible step size that still satisfies Eq. (11), instead of selecting a constant step size that accounts for the worst case scenario, as in strategy 1. For this reason, we found that the largest value of ${\gamma }_{\parallel }=20$ was the best performer for this strategy, where the photon budget remained at 200 photons for each of the 10,000 iterations. For strategy 3, the best performer was a value of ${\gamma }_{\parallel }=10$, where larger values appeared to allow too much variance in the gradient, leading to unstable descents. In all strategies, the recovered absorption distribution matched the ground-truth absorption more closely in the regions of the sample closest to the light source at $z=0$. This is due to the decay of the fluence as a function of depth, as we can see the QPAT signal is highest at shallow depths in Fig. 3(b). The deeper regions of the sample were the last to approach the ground truth in each of the three strategies.

Fig. 3

QPAT inversion: (a) ground-truth absorption distribution, ${\mu }_{\mathrm{a}}^{\text{true}}$, and estimated absorption distribution, ${\mu }_{\mathrm{a}}$, at the point where the photon budget is expended, using each of the three strategies with the stated values of ${\gamma }_{\mathrm{tot}}$ or ${\gamma }_{\parallel }$. (b) Associated measured data from ground-truth medium and simulated forward data at the end of the inversion using each strategy.

Fig. 4

QPAT inversion: (a) sampled cost function, $F_{S_n}$, as a function of iteration, $n$. (b) Error in absorption estimate, $F_{{\mu }_{\mathrm{a}}}$, as a function of iteration, $n$.

Fig. 5

QPAT inversion: (a) step sizes, ${\alpha }_n$, as a function of iteration, $n$. (b) Adaptive sample size, $|S_n|$, as a function of iteration.

From Fig. 6, we see the values of our two metrics $V_{\parallel }$ and $V_{\mathrm{tot}}$. In all cases, both measures of the variance begin at low values, indicating that even with low numbers of photons being simulated, the computed gradients are of reasonable quality, likely due to the poor initial first guess being far from the true solution. Each of the measures of variance increase as the inversion progresses until they begin to violate the norm test or inner product test depending on the strategy. It is seen that the strategy 1 example attempts to keep $V_{\mathrm{tot}}\le 4$, however, due to some level of variation in the metrics themselves, this condition can be seen to be violated regularly, requiring regular updates to the sample size. For strategy 2, the imposed limit of $V_{\parallel }\le 20$ is never violated, and thus the sample size is never required to increased. We also see that strategy 3 manages to keep $V_{\parallel }\le 10$ for the majority of the descent.

Fig. 6

QPAT inversion: (a) $V_{\parallel }$ as a function of iteration and (b) $V_{\mathrm{tot}}$ as a function of iteration.

In addition to these experiments shown in Figs. 3–6, we also trialed a number of other conditions including media with isotropic scatterers (i.e., with $g=0$), various scattering coefficients, and various initial estimates of the absorption. In all cases explored, the methods showed similar behavior as above, but with some differences in the ideal values of ${\gamma }_{\mathrm{tot}}$ and ${\gamma }_{\parallel }$ for each strategy. The outcomes of a range of these experiments are shown in Table 2 for various problem parameters. Strategy 3 was used in all cases in the table, with the same Lipschitz constant ($L=2.5$), starting sample size ($|S_1|=200$ photons), photon budget ($N_{\mathrm{ph}}=2\times {10}^6$ photons), and ground-truth absorption distribution ${\mu }_{\mathrm{a}}^{\text{true}}$ as used in the above examples. The final attained values of the sampled data cost function $F_{S_n}$ and absorption error $F_{{\mu }_{\mathrm{a}}}$ are similar in all cases with the exception of the high asymmetry and low scattering case ($g=0.9$ and ${\mu }_{\mathrm{s}}=4{\mathrm{cm}}^{-1}$). In this case, the reduced scattering coefficient is only ${{\mu }_{\mathrm{s}}}^{\prime }={\mu }_{\mathrm{s}}(1-g)=0.4{\mathrm{cm}}^{-1}$, meaning much lower overall attenuation of the light through the sample. This results in a more uniform data function, $h$, where the simulated photons probe the domain more uniformly, and allows the problem to converge significantly faster than in the higher attenuating cases demonstrated in Figs. 3–6. It is also worth noting that in the regime with low scattering and high scattering asymmetry, it is generally problematic for the performance of approximate transport models such as the diffusion approximation, and the results here highlight the flexibility of RTE based approaches, and the efficiency of the proposed adaptive sampling techniques.

Table 2

Final outcomes of QPAT inversions with various medium optical properties and starting values of μa. Values of FSn and Fμa are the final values at the end of each inversion after the stated number of iterations. In each case, strategy 3 was employed, with a starting sample size of |S1|=200 photons per iteration, and a total photon budget of Nph=2×106 photons. Slab thickness is 2 cm in all cases, with the same ground-truth μatrue distribution as shown in Fig. 3(a).

Finally, interesting behavior was observed when using certain initial guesses of the absorption. An example of this is shown in Fig. 7, where we show the resulting cost functions for a starting estimate of ${\mu }_{\mathrm{a}}=1{\mathrm{cm}}^{-1}$ (significantly overestimating the absorption at all depths), and medium properties of $g=0.9$ and ${\mu }_s=40{\mathrm{cm}}^{-1}$. In this case, we see that the descent appears to encounter local minima in the data cost function $F_{S_n}$ at various points during the descent, depending on the particular strategy used. However, the algorithm manages to escape these local minima and converge to a better solution. This is seen to be the case for all three strategies shown in Fig. 7.

Fig. 7

QPAT inversion: with initial estimate of ${\mu }_{\mathrm{a}}=1.0{\mathrm{cm}}^{-1}$: (a) sampled cost function, $F_{S_n}$, as a function of iteration, $n$. (b) Error in absorption estimate, $F_{{\mu }_{\mathrm{a}}}$, as a function of iteration, $n$.

4.2.

UMOT

Next, we demonstrate similar experiments performed using the UMOT modality described in Sec. 3.2.2 for the transmission geometry. In this setup, we used the same medium slab size as the QPAT example, and the same optical properties apart from the absorption distribution. The starting sample size in all cases shown is $|S_1|=4000$ photons per iteration, 1000 for each forward run (per each of the two sources), and 1000 for each of the two adjoint sources as outlined in Algorithm 3. The total photon budget for the inversion was set to $N_{\mathrm{ph}}=4\times {10}^8$ photons. The Lipschitz constant was set at $L=50$, as this displayed stable descent in our test problems using large photon budgets (low-variance case). The initial estimate of the absorption distribution in the medium is ${\mu }_{\mathrm{a}}=0.1{\mathrm{cm}}^{-1}$ in all layers. The ground-truth absorption ${\mu }_{\mathrm{a}}^{\text{true}}$ is shown in Fig. 8(a), along with the final retrieved absorption distributions obtained via strategies 1, 2, and 3 using the stated values of ${\gamma }_{\mathrm{tot}}$ and ${\gamma }_{\parallel }$. Figure 8(b) shows the true measured UMOT data, $b^{\mathrm{obs}}$, along with the forward modeled data from the final estimated medium for each strategy. Figure 9 shows the outcome of each strategy in terms of the sampled data cost function, $F_{S_n}$, and the absorption error function, $F_{{\mu }_{\mathrm{a}}}$. The relevant step sizes and sample sizes for each of these three examples are shown in Fig. 10, and the values of the metrics measuring the variance in the sampled gradients are presented in Fig. 11. Similar to the QPAT modality, we found that strategy 3 performed the best in terms of the final achieved value of the error in the absorption estimate $F_{{\mu }_{\mathrm{a}}}$.

Fig. 8

UMOT inversion: (a) ground-truth absorption distribution, ${\mu }_{\mathrm{a}}^{\text{true}}$, and recovered absorption distribution ${\mu }_{\mathrm{a}}$ using each of the three strategies with the stated values of ${\gamma }_{\mathrm{tot}}$ or ${\gamma }_{\parallel }$. (b) Associated measured data from ground-truth medium and simulated forward data at the end of the inversion using each strategy.

Fig. 9

UMOT inversion: (a) sampled cost function, $F_{S_n}$, as a function of iteration, $n$. (b) Error in absorption estimate, $F_{{\mu }_{\mathrm{a}}}$, as a function of iteration, $n$.

Fig. 10

UMOT inversion: (a) step sizes, ${\alpha }_n$, as a function of iteration, $n$. (b) Adaptive sample size, $|S_n|$, as a function of iteration.

Fig. 11

UMOT inversion: (a) $V_{\parallel }$ as a function of iteration and (b) $V_{\mathrm{tot}}$ as a function of iteration.

In addition to the results shown in Figs. 8–11, in Table 3 we present a summary of results for a range of different medium optical parameters and starting estimates of the absorption. In all cases, strategy 3 was used, and the starting photon budget was the same as in the previous UMOT examples ($|S_1|=4000$ photons), with a total photon budget of $N_{\mathrm{ph}}=4\times {10}^8$ photons. For each of the inversions presented in this table, we conducted the inner product test once every 50 iterations, using $N_{\mathrm{rep}}=50$ repeated evaluations of the gradient. The resulting inversions display similar error in these cases to the above examples where we used $N_{\mathrm{rep}}=100$ repeated evaluations of the sampled gradient once every 10 iterations to run the inner product test. This demonstrates that the described methods can still be successful when dedicating fewer computational resources to the inner product or norm test metrics, which control the adaptive sample size and step size.

Table 3

Final outcomes of UMOT inversions with various medium optical properties and starting values of μa. Values of FSn and Fμa are the final values at the end of each inversion after the stated number of iterations. In each case, strategy 3 was employed, with a starting sample size of |S1|=4000 photons per iteration, and a total photon budget of Nph=4×108 photons. Slab thickness is 2 cm in all cases, with the same ground-truth μatrue distribution as shown in Fig. 8(a).