Compressively sampling the optical transmission matrix of a multimode fibre

The measurement of the optical transmission matrix (TM) of an opaque material is an advanced form of space-variant aberration correction. Beyond imaging, TM-based methods are emerging in a range of fields, including optical communications, micro-manipulation, and computing. In many cases, the TM is very sensitive to perturbations in the configuration of the scattering medium it represents. Therefore, applications often require an up-to-the-minute characterisation of the fragile TM, typically entailing hundreds to thousands of probe measurements. Here, we explore how these measurement requirements can be relaxed using the framework of compressive sensing, in which the incorporation of prior information enables accurate estimation from fewer measurements than the dimensionality of the TM we aim to reconstruct. Examples of such priors include knowledge of a memory effect linking the input and output fields, an approximate model of the optical system, or a recent but degraded TM measurement. We demonstrate this concept by reconstructing the full-size TM of a multimode fibre supporting 754 modes at compression ratios down to ∼5% with good fidelity. We show that in this case, imaging is still possible using TMs reconstructed at compression ratios down to ∼1% (eight probe measurements). This compressive TM sampling strategy is quite general and may be applied to a variety of other scattering samples, including diffusers, thin layers of tissue, fibre optics of any refractive profile, and reflections from opaque walls. These approaches offer a route towards the measurement of high-dimensional TMs either quickly or with access to limited numbers of measurements.

S1 Transformation matrix from spots to PIMs. S2 Schematic of the experimental set-up. S3 Fully-sampled TM in PIM basis reconstructed using the column-wise method. S4 Under-sampled TM in PIM basis reconstructed using FISTA with sparsity and support priors. S5 Robustness of reconstruction to inaccurate estimation of σ and σ p . S6 Tikhonov regularisation results. S7 Movie 1: Experimentally measured coupling of PIMs.
FIG. S1: Transformation matrix from spots to PIMs: Here columns represent spot locations, and rows represent PIMs. Therefore the pixel at row i, col. j represents the amplitude (brightness) and phase (colour) of the overlap between a spot focussed at location j on the input facet, with PIM i. We see that the majority of spot locations excite many PIMs -showing that the spot input basis is relatively incoherent with the PIM basis as required in our compressive sampling protocol.
FIG. S2: Schematic of the experimental set-up: The optical system is based on a Mach-Zehnder interferometer. A laser beam is generated by a 1 mW Helium-Neon laser source operating at a wavelength of 633 nm. The beam is split into target and reference arms using a polarising beamsplitter (PBS2). In the target arm, the beam is spatially filtered and expanded to fill a digital micro-mirror device (DMD) (ViALUX V-7001). The first diffraction order of the DMD is selected by an iris which blocks all other diffraction orders. The DMD chip plane is imaged onto the pupil of an objective lens (OBJ1, 20X magnification). The input facet of the MMF is placed at the front focal plane of the objective lens. The MMF output facet is imaged onto a high-speed camera (CCD3, Basler Pilot GigE, resolution 648x488), where it interferes with the coherent reference beam. The plane wave of the reference beam arrives at the camera at a small tilt angle with respect to the camera chip normal, enabling single-shot digital holography to reconstruct the intensity and phase of the target field. CCD1 and CCD2 are alignment cameras (also Basler Pilot GigE). They are in the image plane of the proximal facet of the MMF. CCD1 images the incident laser beam, enabling aberration correction of the part of the optical set-up before the MMF if necessary, using, for example, the methods described in ref. [? ]. This correction need only take place once and is unchanged regardless of the test scattering object that is placed in the TM measurement system. A red LED illuminator is used to illuminate the proximal facet of the fibre to aid alignment. Transmission imaging was achieved by scanning a focussed beam over a transmissive resolution target placed ∼ 40µm from the distal facet of the fibre. The TM under question was used to calculate the input field required to generate a focussed spot on the output distal facet of the MMF. This was then refocussed from the distal fibre facet to the plane of the resolution target by adding a quadratic Fresnel lens phase function to the hologram displayed on the DMD, as described in ref. [? ]. To reconstruct an image, the total transmitted intensity arriving at CCD3 was recorded for each spot location at the distal facet. Reflection imaging is also possible (as would be necessary for the MMF to be deployed in a real application), but in our case the low power of the laser used for the experiment, coupled with low collection efficiency for MMFs meant the returning signals were small, which introduced additional noise and so we did not use reflection images to test the performance of the reconstructed TMs.
FIG. S3: Fully-sampled TM in PIM basis reconstructed using the column-wise method: This is a larger scale reproduction of main text Fig. 1(a). The TM is over-sampled using c = 2.2 as explained in the main text. We trial nine different estimated supports, with σ = 2 : 2 : 6, and σp = 1 : 1 : 3. The mean power-ratio of foci generated at the output using TMs reconstructed with the nine different supports is annotated. In this case we find that even when the extent of the anticipated off-diagonal power coupling in the TM is over-estimated or under-estimated, the fidelity of the reconstructed TM is always higher than without using a support: pr > 0.8 in every case tested here. This illustrates that the FISTA reconstruction is robust to inaccurate estimates of the support within the demonstrated range.
(b) Here we simulate the TM reconstruction fidelity over a wider range of support parameters for the same compression ratio as (a). We find that the reconstruction can tolerate a significant overestimation in mask parameters in this case, with only gradual reduction in reconstruction fidelity. The white contour at 0.6 indicates the fidelity of the reconstruction using sparsity priors only.
FIG. S6: TM reconstruction using Tikhonov regularisation: In addition to the three reconstruction methods presented in the main paper, we also investigated a reconstruction method based on Tikhonov regularisation. In this case the estimated TM amplitude support was thresholded to create a binary mask: containing 0 in regions in which we expect there to be minimal power in the reconstructed TM, and 1 otherwise. To promote solutions that have low absolute values in regions specified by the predicted mask, the information about which regions of the TM we expect to be zero can be inserted into Eqn. 1 (main text), by adding extra rows to S and y. For example, to specify that the n th entry in t is 0, we vertically concatenate S with an extra row consisting of all elements set to 0 except for the n th element which is set to 1. We also vertically concatenate y with an extra element set to 0. This can be repeated for every element of the TM that we expect to be 0 according to the predicted binary mask. Note that the memory requirements are low for this approach as the extra rows are mainly zeros and can be represented using sparse matrices. As long as enough additional information has been inserted into Eqn. 1 (main text) to render it full-rank, then it can be solved using standard fast methods that minimise an error term η given by the square of the Euclidean norm of the residual: η = St − y 2 2 , which attempts to account for any inconsistencies due to inaccuracies in the estimation of the support or noise in the measurements. Additionally, the strength of the predicted support priors can be weighted with respect to the probe measurements by a factor λ Tik , for example, using the methods described in the SI of ref. [? ]. This method is a form of Tikhonov regularisation: to demonstrate this, let Γ represent the matrix formed from the extra rows vertically concatenated with S as described above, and 0 is a column vector representing the extra rows of zeros vertically concatenated with y. Therefore the new matrix-vector equation becomes: Tik is the weighting factor. In this case the square of the Euclidean norm of the residual is given by S λ Tik Γt 2 2 = St − y 2 2 + λ Tik Γt 2 2 , which is equivalent to Tikhonov regularisation with a Tikhonov matrix Γ and weighting factor λ Tik . Here we set λ Tik = 1. (a) and (b) show the performance of Tikhonov regularisation (diamonds) in comparison with the other three reconstruction strategies (data shown in grey equivalent that in main text Figs.2(c,d)). In this case for compression ratios of c > 0.25, the performance of Tikhonov regularisation is equivalent to FISTA using sparsity and support priors, but the solution is returned approximately an order of magnitude faster. However, the solution is sensitive to under-estimations of the level of off-diagonal power spread in the TM. Additionally, for lower compression ratios of c < 0.2, the above matrix equation is not full-rank, and since no sparsity priors are invoked the Tikhonov reconstruction undergoes catastrophic failure: the number of probe measurements must be large enough to render the matrix equation full rank for an accurate solution to be found. In comparison, all other reconstruction strategies undergo graceful failure as c is reduced. Therefore, although relatively fast to perform, Tikhonov reconstruction can only be successfully applied for mild compression ratios in circumstances where the level of off-diagonal power spread can be safely over-estimated. The lower right hand plot shows the level of power coupling of the n th PIM at the fibre output. We see that power couples locally in ( ,p)-space. We exploit this to predict the support of the TM. We also note that PIMs with low p-indices tend to couple more broadly to the neighbouring PIMs, as the values of the phase velocities of these PIMs are closer. There is potential for this observation to be exploited as more detailed prior knowledge about the structure of the TM in the future.