Single-ended Recovery of Optical fiber Transmission Matrices using
  Neural Networks
Ultra-thin multimode optical fiber imaging promises next-generation medical
endoscopes reaching high image resolution for deep tissues. However, current
technology suffers from severe optical distortion, as the fiber's calibration
is sensitive to bending and temperature and thus requires in vivo
re-measurement with access to a single end only. We present a neural network
(NN)-based approach to reconstruct the fiber's transmission matrix (TM) based
on multi-wavelength reflection-mode measurements. We train two different NN
architectures via a custom loss function insensitive to global
phase-degeneracy: a fully connected NN and convolutional U-Net. We reconstruct
the 64 $\times$ 64 complex-valued fiber TMs through a simulated single-ended
optical fiber with $\leq$ 4\% error and cross-validate on experimentally
measured TMs, demonstrating both wide-field and confocal scanning image
reconstruction with small error. Our TM recovery approach is 4500 times faster,
is more robust to fiber perturbation during characterization, and operates with
non-square TMs.
Yijie Zheng, George S. D. Gordon
# Single-ended Recovery of Optical fiber Transmission Matrices using Neural Networks

###### Abstract

Ultra-thin multimode optical fiber imaging technology promises next-generation medical endoscopes that provide high image resolution deep in the body (e.g. blood vessels, brain). However, this technology suffers from severe optical distortion. The fiber's transmission matrix (TM) calibrates for this distortion but is sensitive to bending and temperature so must be measured immediately prior to imaging, i.e. _in vivo_ and thus with access to a single end only. We present a neural network (NN)-based approach that quickly reconstructs transmission matrices based on multi-wavelength reflection-mode measurements. We introduce a custom loss function insensitive to global phase-degeneracy that enables effective NN training. We then train two different NN architectures, a fully connected NN and convolutional U-Net, to reconstruct \(64\times 64\) complex-valued fiber TMs through a simulated single-ended optical fiber with \(\leq 4\%\) error. This enables image reconstruction with \(\leq 8\%\) error. This TM recovery approach shows advantages compared to conventional TM recovery methods: 4500 times faster; robustness to 6% fiber perturbation during characterization; operation with non-square TMs and no requirement for prior characterization of reflectors.

Optical fiber imaging Transmission matrix reconstruction Custom loss function Neural network

## 1 Introduction

Ultra-thin endoscopes are a promising technique for enabling cell-scale imaging in difficult-to-reach parts of the body, with the potential to improve disease detection in organs such as the pancreas and ovaries. Commercial products using imaging fiber bundles around 1mm diameter are used in bile ducts [1] and flexible and full-color imaging has been demonstrated using distal scanning mechanisms that are typically around 2mm diameter [2; 3; 4]. To further reduce the size of endoscopes, recent work has focused on imaging through ultra-thin multimode fibers with diameters of 0.125mm and has achieved _in vivo_ fluorescence imaging in brains of immobilized mice [5]. However, there are some key limitations of these imaging systems that use ultra-thin optical fiber. First, the thinnest such imaging devices are made using multimode fiber (MMF), which suffers from significant optical distortion that changes whenever the fiber is perturbed, particularly for longer fibers (\(>\)1m) required to reach deep inside the human body [6]. Second, to calibrate this distortion, practical fiber bundle endoscopes require measurement of their transmission matrix (TM) which requires optical components at the distal end to focus the light onto the distal facet. If calibration is required immediately before use, such components would be required on the distal tip for _in vivo_ use, and would thus compromise the ultra-thin form factor of the endoscopes [7].

A number of methods have been proposed to calibrate fiber TMs without distal access including guidestars [8; 9; 10], a virtual beacon source [11], or reflective structures on the fiber tips [7; 12; 13]. Gordon et al. [13] proposed a single-ended method of TM recovery based on the fiber system shown in Figure 1, with a specially designed reflector stack. This approach avoids the need for measurement at both proximal and distal end of the fiber but works for non-unitary TM matrices. The reflection matrix, \(\mathbf{C}_{\lambda}\in\mathbb{C}^{M^{2}\times M^{2}}\), describes how an incident field \(\mathbf{E}_{\mathbf{in}}\in\mathbb{C}^{\mathbb{M}^{2}}\) is transformed via propagation through the optical fiber, reflected by the reflector stack and finally transferred back through the fiber into
an output field \(\mathbf{E_{out}}\in\mathbb{C}^{M^{2}}\) at a wavelength of \(\lambda\):

\[\mathbf{C}_{\lambda}=\mathbf{E_{out}}_{\lambda}\mathbf{E_{in\lambda}^{-1}} \tag{1}\]

Theoretically, the forward TM, \(\mathbf{A}_{\lambda}\in\mathbb{C}^{M^{2}\times M^{2}}\), can be unambiguously reconstructed at a fourth wavelength based on the measured reflection matrices at 3 different wavelengths. Specifically, the reconstruction of TM is achieved by solving a set of three quadratic matrix exponential equations:

\[\mathbf{C}_{\lambda_{1}}=\mathbf{A_{\lambda_{1}}^{T}}\mathbf{R}_{\lambda_{1}} \mathbf{A}_{\lambda_{1}} \tag{2}\]

\[\mathbf{C}_{\lambda_{2}}=(e^{(\log\mathbf{A}_{\lambda_{1}}\frac{\lambda_{1}}{ \lambda_{2}})})^{T}\mathbf{R}_{\lambda_{2}}(e^{(\log\mathbf{A}_{\lambda_{1}} \frac{\lambda_{1}}{\lambda_{2}})}) \tag{3}\]

\[\mathbf{C}_{\lambda_{3}}=(e^{(\log\mathbf{A}_{\lambda_{1}}\frac{\lambda_{1}}{ \lambda_{3}})})^{T}\mathbf{R}_{\lambda_{3}}(e^{(\log\mathbf{A}_{\lambda_{1}} \frac{\lambda_{1}}{\lambda_{3}})}) \tag{4}\]

where, \(\mathbf{A}_{\lambda}\in\mathbb{C}^{M^{2}\times M^{2}}\) is the transmission matrix at wavelength \(\lambda\), \(\mathbf{R}_{\lambda}\in\mathbb{C}^{M^{2}\times M^{2}}\) is the reflector matrix and \(\mathbf{e}^{(\log\mathbf{A}_{\lambda_{1}}\frac{\lambda_{1}}{\lambda_{2}})}\) is the transmission matrix adjusted for a wavelength \(\lambda_{2}\).

Currently, these equations can be solved by using an iterative approach which relies on optimization of the entire TM [13]. This therefore scales in complexity with the square of the matrix dimension, incurring significant computational time, especially for large matrices. In practice, the transmission matrix shows high sensitivity to bending and temperature so in a practical usage scenario would need to be measured very frequently and immediately prior to imaging. Large computational times are therefore not practical.

Considering this, there are several methods that have been developed in order to reduce the computational time for fiber imaging. These methods typically exploit prior knowledge about the fibers to improve or speed up TM reconstruction. For example, Li et al. [14] proposed a compressed sampling method based on the optical transmission matrix to reconstruct full-size TM of a multimode fiber supporting 754 modes at compression ratios down to 5% with good fidelity. Additionally, Huang et al. [15] retrieved the optical transmission matrix of a multimode fiber using the extended Kalman filter, enabling faster reconstruction.

Recently, there has been work on using deep learning approaches, involving convolutional neural networks, to reconstruct images via multimode fibers both in transmission and reflection modes [16; 17; 18]. These methods have the advantage of being fast, and also learning and utilizing important prior information about the fiber properties and the objects being imaged. However, their performance typically degrades significantly under fiber perturbation because they do not have access to reflection calibration measurements required to unambiguously resolve a TM. Further, because such approaches seek to approximate the forward propagation of light and often only consider amplitude image recovery, they often rely on classical mean-squared error loss functions for training.

In order to incorporate reflection calibration measurements following fiber perturbation, it may instead be advantageous to use AI approaches to reconstruct a transmission matrix rather than an image, though there has been relatively little work in this area. When reconstructing a transmission matrix comprising complex numbers, a particular type of degeneracy arises that is not well handled by conventional AI loss functions: a global phase factor. In many physical problems, including the recovery of transmission matrices for the purposes of image reconstruction and phase-hologram generation, global phase factors are not relevant as they do affect the perceived performance of the system: it is the _relative_ phase between pixels that must be preserved. Global phase may have a physical interpretation related to the physical length of the fiber, but in practice it is often arbitrary unless great care is taken. For example in interferometric systems the global phase is likely to be arbitrary unless the optical path lengths of the reference and sample arms are perfectly matched, which is very challenging for multimode fibers. Further, the global phase often drifts significantly during practical experiments [19], and approaches using phase-retrieval produce entirely arbitrary global phase values [20]. Therefore, in many important practical situations, conventional loss functions will convert arbitrary shifts in the global phase to large changes in value, which can confound minimization algorithms used to fit AI models. In such cases, models may arbitrarily learn a global phase factor (a type of 'overfitting') and may thus not be generalisable.

In this paper, we therefore propose a new method of implementing single-ended recovery of an optical fiber TM by solving equations 2-4 based on three reflection matrix measurements at three different wavelengths. Specifically, we present two different neural network architectures, fully connected neural network (FCNN) and convolutional U-net based neural networks, and demonstrate the performance of both. As a necessary step, a custom global phase insensitive loss function is developed to eliminate the effect of global phase factor during the model training process. We first validate our model by recovering \(64\times 64\) complex-valued fiber transmission matrices through a simulated single-ended optical fiber system (shown in Figure 1) with \(\leq 4\%\) error for both FCNN and convolutional U-net architectures. We
then demonstrate reconstructing \(8\times 8\) images through fiber based on recovered TM with \(\leq 8\%\) error. We highlight several advantages of this TM recovery approach compared to conventional TM recovery methods. Firstly, once the model is trained (\(\sim\)100 hours), it only requires \(\sim\)1 second for reconstruction, which is 4500 times faster than the conventional iterative approach. Secondly, the conventional method [13] can only reconstruct square TM problems, whereas this method is compatible with non-square-shaped TM with \(\leq 8\%\) error, which is potential for many practical cases where optical systems may have different mode bases at proximal and distal ends. Thirdly, no prior measurements for reflectors are required for this model, removing a significant experimental challenge.

## 2 Results

### TM recovery

This TM recovery model was trained on a simulated dataset comprising 800,000 sets of simulated reflection matrices, \(\mathbf{C}_{\lambda}\) at 3 wavelengths, \(\lambda=\lambda_{1},\ \lambda_{2},\ \lambda_{3}\), as input and a complex-valued non-unitary transmission matrix at wavelength \(\lambda_{1}\), \(\mathbf{A}_{\lambda_{1}}\) as output. It was then validated using 200,000 such sets not used in training. Figure 2(a) shows the training and validation loss in training the FCNN model over 2500 epochs using different loss functions, namely conventional mean absolute error (MAE), and weighted and unweighted versions of our global-phase insensitive custom loss function (Eq. 9 and 10 respectively). Both global-phase insensitive loss functions show a decreasing loss in both training and validation in the first 2000 epochs and converge after 2500 epochs, whereas the MAE loss function exhibits fluctuating non-converging loss values for both training (green line) and validation (pink line). Comparison between the two versions of our global-phase insensitive loss function shows that the weighted version reduces loss compared to the unweighted version by \(\sim\)10% in both training (blue line) and validation (red line). An example of a reconstructed TM predicted by the FCNN model using weighted loss function at different epochs is shown inset in Figure2(a). It can be seen that the predicted TM is getting closer to the target TM from 300 epochs to 2500 epochs. These indicate that using custom loss function can successfully avoid the global phase degeneracy which would otherwise prevent the model from learning.

Figure 2(b) compares the TM result predicted by our two different neural network architectures using different loss functions. Both FCNN and convolutional U-net-based neural networks cannot recover TM when using the MAE loss function but are capable of recovering TM using both versions of the global phase insensitive loss function, with a loss of \(\leq 4\%\) over 200,000 validation TMs. Compared to the unweighted version, there is a \(\sim 0.3\%\) reduction in error when using the weighted loss function in either FCNN or convolutional U-net architecture. Furthermore, we also

Figure 1: Single-ended optical fiber imaging system for TM recovery. The optical image, \(\mathbf{X}\in\mathbb{C}^{M\times M}\) is placed at far end of distal facet. Light with a field of \(\mathbf{E_{in}}\in\mathbb{C}^{M^{2}}\) propagates from the proximal facet through the optical fiber, with the forward transmission matrix of the optical fiber defined as \(\mathbf{A}_{\lambda}\in\mathbb{C}^{M^{2}\times M^{2}}\) at the wavelength, \(\lambda\). A reflector stack with a three-layer structure is placed at the distal facet, with its reflector matrix defined as \(\mathbf{R}_{\lambda}\in\mathbb{C}^{M^{2}\times M^{2}}\) at the wavelength \(\lambda\). Reflection matrix, \(\mathbf{C}_{\lambda}\in\mathbb{C}^{M^{2}\times M^{2}}\) can be repeatedly measured at three different wavelengths to recover the TM
evaluated the computational resource usage of the two different neural network architectures as shown in Figure2(c). We implemented the training process using Tensorflow 2.0 running on an NVIDIA Tesla V100 GPU. Compared to FCNN, the convolutional U-net shows significant advantages in memory usage because it requires 1000 times fewer trainable parameters and the converge time is reduced by \(20\%\). However, it shows 0.7% larger loss in the validation TMs. Both FCNN and convolutional U-net can recover TM at a loss \(\leq 4\%\), and also enable \(\sim\)1s prediction time.

### Image reconstruction

To evaluate the performance of recovered TMs for image reconstruction, we considered 3 example images denoted \(\mathbf{x}\in\mathbb{C}^{8\times 8}\): an amplitude-only image with a'space invader' pattern, a phase-only digit with a uniform amplitude and a random complex-valued image. Figure 3 shows the image reconstruction results based on recovered TM using FCNN and convolutional U-net networks. It can be seen that all three types of images can be successfully reconstructed based on recovered TMs using both neural network models, with all the image loss \(\leq 8\%\).

### Fiber perturbation

We then evaluate the robustness of our TM recovery model by swapping rows between different reflection matrices, simulating the effect of the TM changing mid-way through characterization. We simulated 10 sets of 64\(\times\)64 reflection matrices with five different perturbation rates indicating the numbers of rows swapped (1/64, 4/64, 8/64, 16/64, and 32/64). Figure 4 shows the TM recovery results and its corresponding image reconstructions for different fiber

Figure 2: (a) Training and validation loss using MAE, unweighted custom loss function and weighted loss function. TM recovery results over different epochs. (b) TM results recovered using two different neural network architectures (i.e. FCNN and convolutional U-net networks), with three different loss functions, namely MAE, unweighted loss function in Eq.9, and weighted loss function in Eq.10. (c) Comparison between FCNN and convolutional U-net architecture in aspects of loss, training time, prediction time, number of converging epochs and number of trainable parameters.
perturbation rates based on our pre-trained TM recovery FCNN model. It can be seen that our TM recovery model is compatible with optical fibers with a small perturbation rate (below 6%) with TM loss \(\leq 8\%\) and image loss \(\leq 15\%\) but performance degrades significantly above this.

### Recovery of non-square TMs

We next examine the important practical case of non-square TMs, e.g. where the desired representation at the distal end of a fiber might be different from that used at the proximal end and may have more elements. To recover a TM \(\mathbf{A}\in\mathbb{C}^{M_{p}\times M_{d}}\), we require that the reflection matrix, \(\mathbf{C}\in\mathbb{C}^{M_{p}\times M_{p}}\) and that the reflector matrix, \(\mathbf{R}\in\mathbb{C}^{M_{d}\times M_{d}}\). \(M_{p}\) and \(M_{d}\) represent the number of elements used for the basis representation at the proximal and distal ends of the fiber respectively. Figure 5 shows one example of recovered non-square-shaped TM \(\in\mathbb{C}^{12\times 6}\) using FCNN and convolutional U-net network, with loss of 5.95% and 9.3% respectively. Theoretically, a tall-matrix structured TM, \(\mathbf{A}\in\mathbb{C}^{M_{p}\times M_{d}}\), (where \(M_{p}>M_{d}\)) can be recovered by reflection matrices with large total elements, and thus producing a better recovery performance with less loss and training data.

### Computational resource usage

As the dimension of recovered images increases, we expect an increase in the TM dimension thus requiring more computational resources. Empirically measured computational resources are plotted in log-scale in Figure 6 (a)-(c): minimum training data, minimum memory usage, and converging time respectively. All indicate a quadratic relationship to the image dimension \(M\) for both FCNN and convolutional U-net models. For practical imaging applications we would desire at least \(32\times 32\) image resolution, giving a \(1024\times 1024\) TM, which would require training with \(>\)10 million examples, leading to memory consumption \(>\)1.5TB for the FCNN. By comparison, the convolutional U-net would require only 1.1TB of memory consumption. Compared to FCNN, convolutional U-net shows potential advantages in

Figure 3: Image reconstruction based on recovered TM. (a) and (b) show the TM and its error respectively. We consider three example images: (c) an amplitude-only ‘space invader’ pattern, (d) a phase-only digit with uniform amplitude, and (e) a random complex-valued image. We compare the reconstruction result recovered by FCNN (second column) and U-net architectures (third column) with the target result shown in the first column.
using 25% fewer memory resources and 20% less training data within 15% less training time. Figure 6(d) compares the prediction time using our neural network model with the conventional methods using iterative optimization approaches [13], where our FCNN model shows significantly less reconstruction time (\(\sim\)1s vs. 1920s for an 12 x 12 transmission matrix ) even with large size images but the conventional method requires increasing time when the size of image increases.

Figure 4: Effect of perturbations of fiber TM during reflection-mode characterization. (a) and (b) show the TM and its error respectively. We consider three example images: (c) an amplitude-only ‘space invader’, (d) a phase-only digit and (e) a random complex-valued image. We compare the reconstruction result recovered by FCNN for increasing levels of fiber pertrubation.

Figure 5: Non-square-shaped TM \(\in\mathbb{C}^{12\times 6}\) recovered by our TM recovery model using FCNN architecture and convolutional U-net model.
## 3 Discussion

We have demonstrated the successful reconstruction of forward fiber TMs based on reflection-mode measurements at multiple wavelengths using a novel neural network based approach encompassing two architectures: a fully-connected neural network and a convolutional U-Net. Previous work applying neural networks to fibers has focussed on image reconstruction as the end goal, but we instead focus on transmission matrix reconstruction. Such an approach is more flexible as the inputs to the network are calibration measurements that reflect a fibers deformation state at any given time - previous image reconstruction approaches have instead learned a static representation of the fiber TM encoded in the neural network weights. Using our approach, the recovered TM will be accurate for the most recent calibration measurements and can be used for high-speed image recovery via conventional matrix operations. Indeed, we demonstrate error values \(\leq\)8% for reconstructing complex limitations. However, one major challenge of recovering the TM in this way is the need to recover a complex-valued TM with a degenerate global phase shift. Previous work on image reconstruction has addressed this problem by training separate networks for amplitude and phase recovery in purely real space and accepting relatively poor performance for phase recovery [16]. Here, we present a novel loss function that is insensitive to this global phase degeneracy and show a high degree of convergence compared to conventional MAE metrics. We believe this metric in itself could find applications in computer-generated holography via neural networks, phase retrieval problems and indeed to image-reconstruction-based neural networks for fiber imaging. Applying this loss function to our single-ended TM recovery problem, we demonstrated the model for reconstructing \(64\times 64\) complex-valued fiber transmission matrices through a single-ended optical fiber system with \(\leq 4\%\) error either using FCNN or U-net based neural network architecture, which is also capable of reconstructing \(8\times 8\) images through fiber based on recovered TM with \(\leq 8\%\) error.

There are several major advantages to our neural network approach compared to previous iterative reconstruction approaches [13]. First, the prediction time is very fast, typically \(\sim\)1s, which makes this a feasible approach for future real-time imaging, over 4500x faster than the existing iterative approach. Training the network is significantly slower, but this would only need to be performed once per fiber for a fixed reflector so could be performed as a one-off initial calibration step. Second, our approach shows robustness to instances where the fiber TM might change part way through characterization measurements, as is likely to happen during real _in vivo_ usage, and can tolerate up to 6% of row swaps between different reflection matrices. This performance could likely be further improved by re-training the network with perturbed examples as input, thus also learning an 'error correction' strategy. Third, the approach can reconstruct non-square transmission matrices. This is important because due to experimental constraints, the sampling basis of the light on the proximal facet is often the pixel basis of the camera used. However, this basis may not be appropriate for imaging at the end of the fiber as it may contain many more elements than modes are supported in the fiber: multimode fiber may only support a few hundred modes depending on wavelength and core diameter. Therefore, to optimize speed and imaging performance it is often desirable to retrieve a TM in the mode basis of the fiber that can easily be addressed

Figure 6: (a) Minimum training data versus the number of image dimensions, plotting in log-scale.(b) Minimum memory usage versus the number of image dimensions, plotting in log-scale. (c) Converging time versus the number of image dimensions. (d) Prediction time of using our TM recovery model and conventional method, plotting in log-scale.
using our camera coordinates: hence a non-square TM. Finally, it is not required to characterize the reflector in advance as this can effectively be inferred based on measurements of the fiber.

However, there are also some trade-offs with our approach. The first trade-off is that the lack of need for pre-characterization of reflector means that the reflector matrix of this is effectively encoded in the network weights. Without careful thought about implementation this could mean that for each different reflector a separate model would have to be trained and since this may require millions of transmission matrices to be measured, this may be experimentally infeasible. One possible solution to this problem is to either pre-characterize reflectors as proposed previously, or else devise a method of reliably manufacturing reflectors with consistent and highly reproducible properties. Most of the different fiber bending conditions could then be simulated using our approach here and so the network could be trained with relatively few experimental measurements. The second trade-off that follows from this is the need for large amounts of experimental transmission and reflection measurements. This could also be alleviated somewhat by forward simulations of fiber, as we have previously found a high degree of alignment between simulated and experimental matrices [13]. Further, experimental and simulated datasets could be combined in a domain-transfer approach [21, 22]. The use of adaptive loss functions, such as in generative-adversarial networks, may further enable convergence on relatively small datasets, or else generate further training data. Third, the training process is very memory-intensive for dealing with large sizes of TM that are typically encountered in imaging applications (e.g. \(1024\times 1024\)), which requires over 1TB for training the recovery model. One possible solution is to develop matrix compression techniques such as Auto-encoder models to reduce the size of our input matrices by extracting the core features into a latent space. Reducing batch size is considered an optimizing method to reduce memory usage but too small batch size will lead to wider fluctuations and thus a larger converging loss with more training time required.

We anticipate this neural network-based TM recovery model with new loss function designed will lead to new machine-learning models that deal with phase information, for example in imaging through optical fiber, holographic imaging and projection, where both phase control and speed are required.

## 4 Methods

We present a new TM recovery method that uses neural networks, instead of using iterative approaches [13], to solve the Equations 2 - 4. Figure 7 shows the schematic of this TM recovery model. Specifically, we first simulated \(N\) optical fibers TMs, \(\mathbf{A}_{\lambda_{1}}\in\mathbb{C}^{M^{2}\times M^{2}}\), at a wavelength of \(\lambda_{1}\) as the ground truth. Then we randomly generated three complex-valued matrices as our reflector matrices, \(\mathbf{R}_{\lambda_{1}}\in\mathbb{C}^{M^{2}\times M^{2}}\), \(\mathbf{R}_{\lambda_{2}}\in\mathbb{C}^{M^{2}\times M^{2}}\), and \(\mathbf{R}_{\lambda_{3}}\in\mathbb{C}^{M^{2}\times M^{2}}\) at wavelengths \(\lambda_{1}\), \(\lambda_{2}\) and \(\lambda_{3}\) respectively. Finally, we generated three reflection matrices, \(\mathbf{C}_{\lambda_{1}}\in\mathbb{C}^{M^{2}\times M^{2}}\), \(\mathbf{C}_{\lambda_{2}}\in\mathbb{C}^{M^{2}\times M^{2}}\), and \(\mathbf{C}_{\lambda_{3}}\in\mathbb{C}^{M^{2}\times M^{2}}\), at these wavelengths which can be calculated using to Equations 2 - 4. Here, we use wavelengths \(\lambda_{1}=850\)nm, \(\lambda_{2}=852\)nm and \(\lambda_{3}=854\)nm as physically realistic values within the TM bandwidth of a typical endoscopic length fiber (\(\sim\)2m) [13]. Each set of 3 reflection matrices, \(\mathbf{C}_{\lambda_{1..3}}\) then forms a single input to our neural network model.

To feed the neural network, we did data pre-processing to both input data (i.e. reflection matrices) and the ground truth (i.e. TMs), converting from complex-valued data to real-valued data. We then split the \(N\) data into training and validation in a 3:1 ratio before training the neural network model with the ADAM optimizer using our custom-defined loss function. Python was used for model training and MATLAB was used for data pre-processing and post-processing because of its ease of use for complex matrix computations.

To gauge the accuracy of our TM reconstruction we define loss metric to evaluate the performance of TM recovery by calculating the mean MSE of each validated TM over the number of the validation data:

\[Loss=\frac{1}{0.25N}MSE(\hat{A}_{t},A_{t}) \tag{5}\]

where,\(\hat{A}_{t}\) is the recovered TM and \(A_{t}\) is the target TM. \(0.25N\) is the number of data used for validation.

To evaluate the performance of TM recovery in a context relevant to practical endoscopy applications, we passed an image \(\mathbf{X}\in\mathbb{C}^{M^{2}}\) via the simulated optical fiber TM and compared with the image produced using the recovered TM. In this paper, we ignore any loss within the space from the image plane to the distal end of the fiber, and also the loss in transferring through the reflector stack. Theoretically, the reconstructed image, \(\mathbf{\hat{X}}\in\mathbb{C}^{M\times M}\) can be calculated by:

\[\mathbf{\hat{X}}=(\mathbf{\hat{A_{t}}^{\ T}})^{-1}\mathbf{A_{t}}^{\ T} \mathbf{X} \tag{6}\]

where, \(\mathbf{\hat{X}}\) is the reconstructed image, \(\mathbf{X}\) is the target image, \(\mathbf{\hat{A_{t}}}\) and \(\mathbf{A_{t}}\) are the recovered TM and target TM respectively.
### Network architectures

We defined two neural network models: a Fully-connected neural network (FCNN) and convolution U-net based neural network as shown in Figure 8. The FCNN is a ten-layer densely connected neural network (eight hidden layers), including 32,768 neurons in first and last hidden layers and 8192 neurons in other layers, all with LeakyRelu activation function. Figure 8 shows the FCNN architecture, where reflection process matrices \(\mathbf{C}_{\lambda_{1}}\in\mathbb{R}^{128\times 128}\), \(\mathbf{C}_{\lambda_{2}}\in\mathbb{R}^{128\times 128}\), \(\mathbf{C}_{\lambda_{3}}\in\mathbb{R}^{128\times 128}\), are firstly flattened into \(1D\) arrays and then concatenated as the input of the model (with the size of \(49152\times 1\)) and transmission matrix \(\mathbf{A}_{\lambda_{1}}\in\mathbb{R}^{128\times 128}\), flattened into \(1D\) array as the output (with the size of \(16384\times 1\)). Batch normalization layers were defined between every dense layer and dropout layers at the rate of \(0.2\) were defined after the first two dense layers. Also two skip connections were developed in order to prevent the model overfitting. The model was trained iteratively with the weighted 'global phase insensitive' custom loss function used. The training dataset for recovering \(64\times 64\) TM consisted of 500,000 matrices and the model was run for 2500 epochs, taking 182.5 hours using Tensorflow 2.0 running on a NVIDIA Tesla V100 GPU. The Adam optimizer was used with a learning rate of 0.004 in a decay rate of \(1e^{-4}\).

Next, we developed a U-net-based model that used encoder-decoder architecture, including seven Conv2D and DeConv2D layers respectively and two MaxPooling and UpSampling layers respectively with LeakyRelu activation function in each layer. Figure 8 shows this architecture, where reflection process matrices \(\mathbf{C}_{\lambda_{1}}\in\mathbb{R}^{128\times 128}\), \(\mathbf{C}_{\lambda_{2}}\in\mathbb{R}^{128\times 128}\), \(\mathbf{C}_{\lambda_{3}}\in\mathbb{R}^{128\times 128}\), \(\mathbf{C}_{\lambda_{2}}\in\mathbb{R}^{128\times 128}\), are defined in three channels as the input of the model (with the size of \(128\times 128\times 3\)) and transmission matrix \(\mathbf{A}_{\lambda_{1}}\in\mathbb{R}^{128\times 128}\), as the output (with the size of \(128\times 128\times 1\)). Batch normalization layers were defined between every layer and dropout layers at the rate of \(0.2\) were defined after the second and last second Conv layers. Also three skip connections were developed in order to prevent the model being overfitting. The model was trained iteratively with the weighted 'global phase insensitive' custom loss function defined. Also, 2200 epochs were used for training 400,000 training datasets using 143h. The Adam optimizer was used with a learning rate of 0.004 in a decay rate of 1e-4.

### Data Preparation

In terms of data preparation, we first simulated \(N\) pairs of complex-valued transmission matrices \(\mathbf{A}_{\lambda_{1}}\in\mathbb{C}^{64\times 64}\) at a wavelength of \(\lambda_{1}=850nm\) as the ground truth of the model. To simulate these we devised a model that recreates some

Figure 7: Schematic of TM recovery model, including (a) data generation, (b) data pre-processing and (c) model training. \(N\) pairs of TM are firstly simulated as the ground truth. \(\lambda_{1}\), \(\lambda_{2}\) and \(\lambda_{3}\) represent three different wavelengths (in our case, 850nm, 852nm and 854nm). The input of the model is all real-valued matrices concatenated with reflection matrices at three different wavelengths. \(L\) represents custom loss function and \(w\) represents the weight updated by the optimizer.
characteristic properties found in fiber TMs. First, TMs are sparse in some commonly used basis e.g. LP modes for multimode fibers or pixel basis for multicore fibers [23]. Second, TMs can be arranged such that the majority of power intensity lies along the main diagonal with additional power spread along sub-diagonals, which is also typically observed when using bases that match relatively well to the fiber eigenbasis[24]. Third, TMs should be slightly non-unitary in realistic situations, with mode-dependent loss values (i.e. condition numbers) in the range of 3-5. To meet these requirements firstly, we generate a uniformly distributed random tri-diagonal matrix, \(\mathbf{B}\in\mathbb{C}^{6\times 64}\), which has non-zero elements only at the main diagonal, diagonal below and above it. We then compute the left singular matrix \(\mathbf{U}\in\mathbb{C}^{64\times 64}\) and right singular matrix \(\mathbf{V}\in\mathbb{C}^{64\times 64}\) via singular value decomposition (SVD). To make it a non-unitary matrix, we apply a new singular value distribution, \(\mathbf{S_{new}}\in\mathbb{R}^{64\times 64}\), a diagonal matrix that contains random values at its diagonal ranging from 0.5 to 2.5 to simulate our expected TM, which matched with those TM that were measured during the experiments [24]:

\[\mathbf{A}=\mathbf{U}*\mathbf{S_{new}}*\mathbf{V}^{T} \tag{7}\]

We next simulated three complex-valued reflector matrices with random uniformly distributed complex entries, \(\mathbf{R}_{\lambda_{1}}\in\mathbb{C}^{64\times 64}\), \(\mathbf{R}_{\lambda_{2}}\in\mathbb{C}^{64\times 64}\), and \(\mathbf{R}_{\lambda_{3}}\in\mathbb{C}^{64\times 64}\). Based on this, we generate N pairs of complex-valued reflection process matrix \(\mathbf{C}_{\lambda_{1}}\in\mathbb{C}^{64\times 64}\), \(\mathbf{C}_{\lambda_{2}}\in\mathbb{C}^{64\times 64}\), and \(\mathbf{C}_{\lambda_{3}}\in\mathbb{C}^{64\times 64}\) at three different wavelengths \(\lambda_{1}=850nm\), \(\lambda_{2}=852nm\) and \(\lambda_{3}=854nm\) that are corresponded to previously simulated \(\mathbf{A}_{\lambda_{1}}\in\mathbb{C}^{64\times 64}\) using Equations 2 - 4.

In order to feed the neural network, both input data, reflection matrices, and ground truth, TM, are required to be all real-valued matrices. Basically, a \(2\times 2\) complex-valued matrix can be described as \(4\times 4\) all real-valued matrix in Equation 8:

\[\begin{bmatrix}a+bi&c+di\\ e+fi&g+hi\end{bmatrix}=\begin{bmatrix}a&-b&c&-d\\ b&a&d&c\\ e&-f&g&-h\\ f&e&h&g\end{bmatrix} \tag{8}\]

Finally, the input of the model, three \(\mathbf{C}_{\lambda}\in\mathbb{C}^{128\times 128}\) at different wavelengths are normalized using in the range from -1 and 1.

### Weighted global phase insensitive loss function

Widely-used conventional loss functions such as mean absolute error (MAE) or mean squared error (MSE) calculate the absolute difference between predicted and target output values. However, there is a class of problems whose solutions trained by deep learning models are degenerate within a global phase factor, but whose relative phase between pixels must be preserved. This class includes problems where complex transmission matrices are reconstructed and relative phase, but not global phase important, but could extend to phase-hologram generation algorithms where replay-field phase is important. This is depicted visually in Figure 9(a), which shows one example of a pair of predicted and target matrices with complex entries depicted as vectors. Figure 9(b) shows the complex error between these two matrices when using MAE as the loss function. Due to the global phase shift, we observe that the vectors have large magnitudes, which will lead to an overall very large MAE when their magnitudes are summed. In the limiting case (e.g. when the phase shift is of \(\pi\)) where the predicted and target matrices are identical, this global phase shift can result in a

Figure 8: Architectures of two different neural network models used for recovery TM. (a) Fully-connected neural network, (b) Convolutional U-net
normalized MAE of 100% when the true value should be 0%. To avoid this problem, we propose a custom loss function termed a 'global phase insensitive' loss function:

\[L(\widehat{\mathbf{A_{t}}}(w),\mathbf{A_{t}}(w))=\sum_{t=1}^{4M^{2}}\left| \widehat{\mathbf{A_{t}}}(w)-\mathbf{A_{t}}(w)e^{1i\phi(\sum\mathbf{A_{t}}(w) \odot(\widehat{\mathbf{A_{t}}}(w)+\beta))}\right| \tag{9}\]

where, \(\widehat{\mathbf{A_{t}}}(w)\in\mathbb{C}^{M^{2}\times M^{2}}\) and \(\mathbf{A_{t}}(w)\in\mathbb{C}^{M^{2}\times M^{2}}\) represent predicted and target output value with regards to weight, \(w\), respectively, \(\sum\) represents calculating sum over all matrix elements, \(\phi\) represents the argument function for a complex number input, \(\odot\) represents element-wise division and \(\beta=0.001\) is a constant added to avoid divide-by-zero errors.

We also developed an alternative custom loss function that weights phase entries by power intensity, achieved by multiplying by the complex conjugate of \(\widehat{\mathbf{A_{t}}}(w)\), denoted \(\widehat{\mathbf{A_{t}}}^{*}(w)\). We also add an \(\ell_{2}\) regularization term to give the generalization of the model.

\[WL(\widehat{\mathbf{A_{t}}}(w),\mathbf{A_{t}}(w))=\sum_{t=1}^{4M^{2}}\left| \widehat{\mathbf{A_{t}}}(w)-\mathbf{A_{t}}(w)e^{1i\phi(\sum\mathbf{A_{t}}(w) \odot\widehat{\mathbf{A_{t}}}^{*}(w))}\right|+\frac{\alpha}{2}\|w\|^{2} \tag{10}\]

where \(\odot\) represents element-wise multiplication, \(\alpha\) = \(1e^{-4}\) is the regularization parameter. This implicitly weights the phase contributions by the product of magnitudes of the respective elements in \(\widehat{\mathbf{A_{t}}}(w)\) and \(\mathbf{A_{t}}(w)\), which upon convergence will approximately equal the squared magnitude of the target.

Specifically, the global phase factor, estimated by the term \(e^{1i\phi(\sum\mathbf{A_{t}}(w)\odot\widehat{\mathbf{A_{t}}}^{*}(w))}\) is the phase of a complex number representing the weighted sum of the elements of the complex difference matrix between predicted and target matrices. The rationale for this is that when the optimization algorithm has reached a minimum, in the ideal case the remaining error for each complex element will be entirely due to aleatoric uncertainty and can thus be modelled using a circularly symmetric complex Gaussian distribution [25]. The element-wise phase errors should therefore be uniformly distributed from \(0\) to \(2\pi\). If this is not the case, then there is likely some contribution to the phase error from an arbitrary global phase, as shown in Figure 9(c). Correcting for this factor should produce the desired uniform phase distribution.

To estimate the correction factor, the element-wise complex errors can be summed, as shown in Figure 9(c). This will produce an overall complex factor that has the desired global phase shift, shown in Figure 9(d). The predicted output value can be corrected by multiplying by this phase factor as shown in Figure 9(e), the result of which is then used to compute further parameter updates in the gradient descent algorithm. It can be seen that the complex error in Figure 9(f) between the predicted and target output value is reduced to a minimum after removing the phase factor compared to that calculated by MAE. We then compared the absolute values of the complex error calculated by MAE (green bar) and our customized weighted 'global phase insensitive' loss function (blue bar) respectively over 100,000 pairs of predicted and desired TM as shown in Figure9(g). The error using the custom loss function is more than two times smaller than that of the conventional loss function (MAE), which suggests the potential for this custom loss function in eliminating the effect of global phase.
## Data Availability

The data presented in this study are available from the following source: [DOI to be inserted later].

## Code Availability

The code for this study is available from the following source: [DOI to be inserted later].

## Author Contributions

## Acknowledgement

The authors acknowledge support from a UKRI Future Leaders Fellowship (MR/T041951/1).
