## 2 Materials and Methods
### Inverse Problems by an Iterative Scheme  
We suppose that on the surface of a two- or three-dimensional random medium, for each source at \(\mathbf{\rho}^{i}_{s}\) or \(\mathbf{r}^{i}_{s}\) (\(i=1,\ldots,M_{s}\)) exiting light is detected at \(\mathbf{\rho}^{j}_{d}\) or \(\mathbf{r}^{j}_{d}\) (\(j=1,\ldots,M_{d}\)) and is measured at times \(t^{k}\) (\(k=1,\ldots,M_{t}\)). Let \(y_{ijk}\) be measured data whereas \(u\) is a solution to the diffusion equation. Let us suppose that \(u\) depends on a vector \(\mathbf{a}\) which contains unknown parameters. We wish to reconstruct \(\mathbf{a}\). In Section 3.1, \(\mathbf{a}=(\mu_{a},D)\), and \(a\) is a scalar (a one-dimensional vector) in Section 3.2. For example, in the former case the solution \(u\) depends on \(\mathbf{a}\) since the calculated value of \(u\) depends on \(\mu_{a},D\).  
Let us introduce vectors \(\mathbf{U}\) and \(\mathbf{F}(\mathbf{a})\) of dimension \(M_{s}M_{d}M_{t}\) as  
\[\{\mathbf{U}\}_{ijk}=y_{ijk},\qquad\{\mathbf{F}(\mathbf{a})\}_{ijk}=u(\mathbf{ \rho}^{j}_{d},t^{k};\mathbf{\rho}^{i}_{s};\mathbf{a})\text{ or }u(\mathbf{r}^{j}_{d},t^{k};\mathbf{r}^{j}_{s};\mathbf{a}), \tag{22}\]  
where we wrote \(u(\mathbf{\rho}^{j}_{d},t^{k};\mathbf{\rho}^{i}_{s})=u(\mathbf{\rho}^{j}_{d},t^{k};\mathbf{ \rho}^{j}_{s};\mathbf{a})\) and \(u(\mathbf{r}^{j}_{d},t^{k};\mathbf{r}^{j}_{s})=u(\mathbf{r}^{j}_{d},t^{k};\mathbf{r}^{j}_{s}; \mathbf{a})\) emphasizing that \(u\) depends on \(\mathbf{a}\). We find optimal \(\mathbf{a}\) by minimizing \(\|\mathbf{U}-\mathbf{F}(\mathbf{a})\|^{2}_{2}\). Here we particularly consider the Levenberg-Marquardt method [3, 4], i.e., the reconstructed value \(\mathbf{a}^{*}=\lim_{k\to\infty}\mathbf{a}_{k}\) is computed by the iteration given by  
\[\mathbf{a}_{k+1}=\mathbf{a}_{k}+\left[F^{\prime}(\mathbf{a}_{k})^{T}F^{\prime} (\mathbf{a}_{k})+\lambda I\right]^{-1}F^{\prime}(\mathbf{a}_{k})^{T}\left( \mathbf{U}-\mathbf{F}(\mathbf{a}_{k})\right), \tag{23}\]where \(F^{\prime}(\mathbf{a})\) is the Jacobian matrix, which contains derivatives of \(\mathbf{F}(\mathbf{a})\) with respect to \(\mathbf{a}\), and the parameter \(\lambda\) is nonnegative. By modifying the original algorithm according to Ref. [32], our algorithm of the Levenberg-Marquardt method, which we call Algorithm 1, is described below.  
```
1. Set \(k=0\) and \(\lambda=1\).
2. Calculate \(\mathbf{F}(\mathbf{a}_{k})\) and \(F^{\prime}(\mathbf{a}_{k})\).
3. Calculate \(S(\mathbf{a}_{k})=\mathbf{d}^{T}\mathbf{d}\), where \(\mathbf{d}=\mathbf{U}-\mathbf{F}(\mathbf{a}_{k})\).
4. Prepare \(A=F^{\prime}(\mathbf{a}_{k})^{T}F^{\prime}(\mathbf{a}_{k})\) and \(\mathbf{v}=F^{\prime}(\mathbf{a}_{k})^{T}\mathbf{d}\).
5. Find \(\delta\) from \((A+\lambda I)\delta=-\mathbf{v}\).
6. Obtain \(S(\mathbf{a}_{k}+\delta)\) and \(R=\frac{S(\mathbf{a}_{k})-S(\mathbf{a}_{k}+\delta)}{\delta^{T}(2\mathbf{v}+ \lambda\delta)}\).
7. If \(R<0.25\), then set \(\
u=10\) (\(\alpha_{c}<0.1\)), \(\
u=1/\alpha_{c}\) (\(0.1\leq\alpha_{c}\leq 0.5\)), or \(\
u=2\) (\(\alpha_{c}>0.5\)), where \(\alpha_{c}=\left[2-\left(S(\mathbf{a}_{k}+\delta)-S(\mathbf{a}_{k})\right)/ \delta^{T}\mathbf{v}\right]^{-1}\). If \(R<0.25\) and \(\lambda=0\), set \(\lambda=1/\left\|A^{-1}\right\|\) and \(\
u=\
u/2\). In the case of \(R<0.25\), we set \(\lambda=\
u\lambda\). If \(R>0.75\), then we set \(\lambda=\lambda/2\). If \(R>0.75\) and \(\lambda<1/\left\|A^{-1}\right\|\), set \(\lambda=0\). Otherwise when \(0.25\leq R\leq 0.75\), no update for \(\lambda\).
8. If \(S(\mathbf{a}_{k}+\delta)\geq S(\mathbf{a}_{k})\), then return to Step 5.
9. If \(S(\mathbf{a}_{k}+\delta)<S(\mathbf{a}_{k})\), set \(\mathbf{a}_{k+1}=\mathbf{a}_{k}+\delta\). Then put \(k+1\to k\) and go back to Step 2. Repeat the above procedure until one of the stopping criteria \(\left\|\delta\right\|<10^{-4}\) and \(S<10^{-14}\) is fulfilled.
```  
**Algorithm 1**Leverberg-Marquardt (LM)