Online learning of the transfer matrix of dynamic scattering media: wavefront shaping meets multidimensional time series

Thanks to the latest advancements in wavefront shaping, optical methods have proven crucial to achieve imaging and control light in multiply scattering media, like biological tissues. However, the stability times of living biological specimens often prevent such methods from gaining insights into relevant functioning mechanisms in cellular and organ systems. Here we present a recursive and online optimization routine, borrowed from time series analysis, to optimally track the transfer matrix of dynamic scattering media over arbitrarily long timescales. While preserving the advantages of both optimization-based routines and transfer-matrix measurements, it operates in a memory-efficient manner. Because it can be readily implemented in existing wavefront shaping setups, featuring amplitude and/or phase modulation and phase-resolved or intensity-only acquisition, it paves the way for efficient optical investigations of living biological specimens.

Thanks to the latest advancements in wavefront shaping, optical methods have proven crucial to achieve imaging and control light in multiply scattering media, like biological tissues. However, the stability times of living biological specimens often prevent such methods from gaining insights into relevant functioning mechanisms in cellular and organ systems. Here we present a recursive and online optimization routine, borrowed from time series analysis, to optimally track the transfer matrix of dynamic scattering media over arbitrarily long timescales. While preserving the advantages of both optimization-based routines and transfer-matrix measurements, it operates in a memoryefficient manner. Because it can be readily implemented in existing wavefront shaping setups, featuring amplitude and/or phase modulation and phase-resolved or intensity-only acquisition, it paves the way for efficient optical investigations of living biological specimens.

I. INTRODUCTION
Optical methods are an irreplaceable tool to investigate biological media. They deliver images at numerous contrast mechanisms [1], and can activate injected biomolecules [2] and fluorescent markers [3]. However, precisely delivering light in space and time through biological tissues is not straightforward, as photons get multiply scattered by heterogeneities of tissues, limiting their penetration depth [4].
Another current challenge lies in tracking the scattering behaviour of living specimens, with decorrelation times up to only a few ms [5]. This proves crucial to understand the functioning mechanisms of cells and organisms, which requires their observation at extremely different timescales, from nanoseconds (at a molecular level) to minutes (for organ systems) [6]. The need for fast data acquisitions results, in turn, in measurements with inherently low signal-to-noise ratios, and requires solving long and multidimensional time series [7], whose prohibitive size can make their evaluation problematic.
Wavefront shaping techniques have established themselves as the tools of choice to guide light in scattering media [8]. The transmission of arbitrary fields [9], pointspread-function (PSF) engineering [10], imaging [11], as well as tuning energy transmission through scattering media [12], become all accessible if the transfer matrix of the medium is measured [8,13]. However, conventional methods to retrieve the transfer matrix yield sub-optimal solutions in noisy environments [8]. Those optimization routines which can compensate for noise in the transfer matrix [14], however, require storing in memory the whole history of past measurements, making them unsuited with long streams of data.
Iterative, optimization-based, sequential algorithms to focus through scattering media yield an increase in the focus intensity already at their early iterations, which makes them the preferred option on dynamic media. Importantly, they are cast as recursive procedures, i.e., computing the new estimate of the solution only requires the previous estimate and the new data point. Unfortunately, their stochastic nature makes optimization over a set of output modes less reliable and the transmission of arbitrary fields prohibitive. Moreover, these procedures rely on maximizing a given metric, limiting light control to one predefined task. Various implementations derived from genetic algorithms [15,16] have shown better resilience to noise than sequential algorithms, however at the cost of a higher computational complexity and careful choice of several adjustable parameters.
In signal processing, communications and finance, where most datasets are multidimensional time series, the recursive least-squares (RLS) algorithm has played a central role for system identification and prediction [17][18][19]. It allows optimal learning of linear predictors in an online manner-predictors are updated every time a new piece of data is sequentially made available, however past data do not need to be stored in memory. Consequently, its computational complexity is independent of the length of the time series, so iterations can be run over and over, ideally at the same rate as data acquisition (real-time operation).
Here, we demonstrate that the RLS algorithm represents a valuable tool to optimally estimate the transfer matrix of dynamic scattering media online and recursively. The least-squares optimization ensures resilience to noise. The algorithm is provided with a tunable memory, such that the dynamics of the scattering medium is accounted for. By doing so only the most reliable data points, i.e., those acquired within the stability time of the medium, are used during the optimization. We justify how the RLS model can fit a wide variety of dynamic mechanisms happening in scattering media. Its performance is showcased with both simulated and experimental results, tracking the transmission matrix and the time-gated reflection matrix at realistic noise levels and stability times. We further show how light optimization can be achieved with binary amplitude or phase modulation and with phase-resolved or intensity-only measurements. Based on its computational complexity, we discuss its feasibility for light control in living biological specimens at large fields of view. Its simple implementation and the low number of adjustable parameters (whose choice is motivated in the next sections) make our proposed method readily applicable in existing wavefront shaping setups.

II. METHODS
The method bears similarities with conventional routines for the measurement of the transfer matrix, and its working principle is graphically summarized in Fig.  1(a). However, here we allow the transfer matrix X t ∈ C M ×N of the scattering medium to be dynamic, where we have denoted the number of output and input degrees of freedom with M and N , respectively. At every time step t, while probing the medium with the input a t ∈ C N and collecting the corresponding output y t = X t a t ∈ C M , we aim to solve the optimization prob- and where || · || and || · || F denote the L 2 -norm of a vector and the Frobenius norm of a matrix, respectively. Although for sake of generality the inputs and the outputs are assumed to be complex, we will also report an implementation where they are real, meaning that only the amplitude of the input beam is modulated and the intensity of the output fields is measured. Equation 1 is a linear least-squares loss function, featuring Tikhonov regularization via the regularization constant δ. Note, however, that each data-fidelity term ||y τ − X t a τ || 2 is exponentially weighted in time, such that the old pieces of data (corresponding to τ t) are less relevant than the most recent ones in the current estimation of the transfer matrix at time t. In other words, the forgetting factor λ ≤ 1 endows the algorithm with a memory, which allows it to cope with dynamic transfer matrices-at every time step t, the optimization problem is solved anew, using the whole history of past data, where more contribution is given to newest data. Evidently, in the case of a static scattering medium, all measurements can be equally trusted, thus Eq. (1) reduces to a typical regularized linear least-squares problem upon setting λ = 1. Once λ and δ are fixed, the least-squares problem has a unique solution, provided the inputs are linearly independent, which is the case in conventional transfermatrix measurements, where the inputs are drawn from the Hadamard basis of order N .
The choice of exponential weights for Eq. (1) is motivated by the physics of our problem. We aim to follow the evolution of the transfer matrix of dynamic scattering media, subjected to uncorrelated variations, whereby the total transferred power fraction is constant in time. These conditions apply in a wide variety of dynamic mechanisms in scattering media investigated with visible and near-infrared light, e.g. whenever their inner scatterers move due to functional changes [5,20], or even when the sample drifts away from its initial position, suggesting that our method can also be used as an online calibration tool of imaging systems. In all these situations, the transfer matrix can indeed be described by the time series [21], where we assume that both the transfer matrix and the perturbation matrix P t are random variables independently drawn from complex Gaussian distributions with zero mean and constant variance σ 2 X and σ 2 P , respectively [22]. Equation 2 denotes an autoregressive model of order 1, AR(1), whose autocovariance is proportional to (σ X / σ 2 X + σ 2 P ) t , justifying our exponentially weighted model of Eq. (1). When focusing through dynamic scattering media following Eq. (2), the stability time of the enhancement is proportional to σ −2 P [21]. This means that the optimal weight λ should follow the same dependence, thus in principle requiring the knowledge of the rate of change of the scattering medium. A strategy for automatically tuning the forgetting factor will be discussed in section IV.
Crucially, minimizing the loss function of Eq. (1) does not require storing the whole history of past data. This becomes apparent if we recall that the linear least-squares estimate of X t ,X t , satisfies the normal equations, with the covariance matrix of inputs and the crosscovariance matrix at time t respectively defined as, with I N denoting the identity matrix of order N and the superscript H standing for Hermitian transposition. The quantities calculated in Eqs. (4) can be both estimated recursively, as follows: Equations (5) mean the right-hand side of Eq. (1) can be minimized from the previous estimates of the covariance and cross-covariance matrices and the new piece of data (a t , y t ). It becomes now clear how the RLS algorithm combines the benefits of transfer-matrix-based and optimization approaches. Using a recursive procedure, a typical asset of, e.g., the continuous sequential algorithm (CSA), the partitioning algorithm [21], or more computationally intense genetic algorithms [15], the full X t is estimated in parallel at all output pixels, thereby preserving all light-control capabilities allowed by the knowledge of the transfer matrix [10][11][12]14] [ Fig. 1(b)]. In principle, the transfer matrix could be obtained from Eq. (3) However, in what follows we will implement the inverse QR-decomposition-based RLS (abbreviated as inverse QRD-RLS) algorithm [23]. Because it avoids matrix inversions and it always preserves the non-negativeness of the covariance matrix, it possesses higher numerical stability than directly inverting Eq. (3). Overall, it boils down to performing a QR decomposition of a matrix constructed from the new data and the previous estimate of the square root of the inverse covariance matrix. This results in few lines of code which can be readily implemented in any programming language using standard libraries or built-in functions (see the box Algorithm 1 and the corresponding code available at Ref. [24]). As can be seen from Eq. (1) and Algorithm 1, the regularization constant δ is used to construct the initial estimate of the square root of the inverse correlation matrix, hence it mostly impacts the convergence speed at early iterations. In section IV, the choice of its value will be discussed.
previous estimate of the transfer matrixXt−1, previous estimate of the square root of the inverse covariance matrix (C −1 t−1 ) 1/2 , forgetting factor λ /* Construction of the matrix U */

III. EXPERIMENTS
Figure 1(c) shows the sketches of the experimental implementations used to demonstrate our method. Both are based on phase-shifting digital holography to retrieve the complex output fields y t after interacting with a multiply scattering medium. The medium is an opaque deposit of ZnO nanoparticles (size < 100 nm), whose thickness (20 µm) is ∼5 transport mean free paths, ensuring full mixing of its optical modes at the output. The input fields are shaped via a reflective, phase-only and liquidcrystal-based spatial light modulator (SLM, Meadowlark Optics HSP512L-1064) and focused on the scattering medium with an objective with a numerical aperture of 0.4 (Olympus PLN20X). A region-of-interest containing ∼80 speckle grains is imaged onto a CCD camera (Manta G-046B, Allied Vision) via a tube lens, yielding a pixel size of 0.2 µm at the CCD plane. Before impinging onto the SLM, part of the beam is redirected along a reference arm with a polarizing beam splitter (PBS), and subsequently recombined with the scattered beam through a beam splitter (BS). The relative power of the two beams, yielding the maximum interference contrast, is adjusted via two half-wave plates, one along the common path and one along the reference arm, while a polarizer in front of the camera filters out any potential residual bal-listic component traveling along with the scattered beam.
In the experiments in transmission [ Fig. 1(c), left], the beam exiting the scattering medium is collected at a distance of ∼1.5 mm, where a fully developed speckle pattern was observed, with another Olympus PLN20X 0.4 NA objective. The light source (MaiTai HP Ti:Sapphire laser, Spectra-Physics) is set to monochromatic operation mode at a wavelength of 808 nm. The experiments in reflection [ Fig. 1(c), right] reproduce a typical optical coherence tomography (OCT) setup, whereby ultrashort pulses (with a central wavelength of 808 nm and a duration of 100 fs) are sent through the scattering medium and the backscattered, elongated pulses are gated at a time delay set by a delay line along the reference arm. Dynamics is introduced by transversally translating the scattering medium across to the incident beam, with steps following a two-dimensional random walk with random Gaussian increments, where the standard deviation determines the stability time of the medium. More details on it will be provided in the next section. Figure 2 summarizes the performance of the RLS algorithm for the online estimation of the transmission matrix. The beam incident onto the SLM is modulated according to the Hadamard patterns with N = 64 pixels. Every time an input a t is sent through the scattering medium and the corresponding output field y t is measured, the inverse QRD-RLS update routine of Algorithm 1 is executed, yielding an estimateX t of the transfer matrix. Note, that this procedure can be continuously repeated-after sending the N -th input, the first Hadamard vector or any other known input pattern can be sent. As long as the scattering medium is static, probing it with the same input multiple times corresponds to oversampling the unknown N ×M coefficients of its transfer matrix, thereby improving their estimation. It is indeed known that the covariance of the estimated transfer matrix is inversely proportional to C −1 t , thus decreasing as t −1 [17]. Since the true value X t is unknown, the quality of our reconstruction is evaluated via the intensity of a focus produced behind the scattering medium. We report the intensity enhancement, relative to the average intensity of a non-optimized speckle pattern [8]. The learning curve for a static scattering medium, obtained from the RLS algorithm, is shown as an orange trace in Fig. 2(a). The temporal axis is expressed in units of T T M , which is defined as the time needed to update the estimation of transfer matrix N times. In other words, a conventional transfer matrix experiments lasts T T M . Equivalently, a normalized time of 2 means the oversampling ratio is 2. To showcase the beneficial effect of oversampling, the blue trace shows the performance of a conventional transfer-matrix measurement, using N measurements. At times t/T T M ≤ 1, the two approaches are equivalent-data are not oversampled. At later times, however, one can take advantage of the whole history of past data to build an estimate more resilient against noise. Our values of the enhancement, when compared to the number of input degrees of freedom N , are on a par with previously reported measurements with no oversampling [13,25,26].

IV. RESULTS
The same procedure is repeated with dynamic scattering media. By duly tuning the average speed of their movements, we achieve different stability times T stab (also expressed in units of T T M ). These are estimated as the time constant of an exponential function fitting the tails of the blue traces. At oversampling ratios in the range 3-4, we increase the focus intensity by a factor between 1.5 and 2, compared to the values after a conventional transfer-matrix approach. Upon decreasing T stab , the oversampling ratio decreases too, and the performances of the two approaches gradually match, however the RLS estimation always operates in a memoryefficient manner. With dynamic media, forgetting factors λ < 1 should be used. In our experiments featuring T stab in the range 1-4, we have chosen 1 − λ ≈ 10 −5 , achieving a good compromise between tracking capability and numerical stability. Interestingly, it has been shown that the optimal forgetting factor heavily depends on the number of unknown parameters N which, fortunately, is under user control [27]. Furthermore, the structure of Algorithm 1 suggests that each inverse QRD-RLS iteration may be run at a different value of λ, allowing the user to pick the one yielding the best performance in an online manner, i.e., with no need to restart the optimization anew and using the current enhancement as a feedback to tune the next value of λ. Trivially, the optimal value for static media is instead λ = 1. The best regularization constant δ depends on the signal-to-noise ratio (SNR) of the measurements. In our experiments, the fact that the RLS algorithm is on a par with the conventional transfermatrix approach at t/T T M ≤ 1 and T T M < T stab suggests that the selected regularization constant (here δ = 1) is optimized for the best performance.
In order to test the validity of Eq. (1) and further justify our assumptions, the experiments at the top row of . We simulated a finite SNR by corrupting the outputs y t with additive white complex Gaussian noise with variance σ 2 noise and setting SNR ≡ σ 2 X /σ 2 noise . Overall a good quantitative agreement is obtained. For example, if the SNR is increased by a factor of 2 doubling the number of phasestepped images for field reconstruction, the experimental performance is the one plotted as an inset in Fig. 2(a). The same trend is retrieved by simulating measurements with halved σ 2 noise [inset in Fig. 2(d)]. Analogous results, plotted in Fig. 3, are obtained with the non-invasive OCT setup on the right-hand side of Fig.  1(c), setting the time delay yielding the maximum average gated intensity. In this instance, the transfer matrix is the time-gated reflection matrix [28]. The two learning curves corresponding to the retrieval of the transfer ma- trix are compared to a conventional optimization routine, which can recursively track the changes in the scattering medium, namely the CSA (in cyan). After blocking the beam along the reference arm, we implement a version of the CSA modulating half of the SLM pixels (corresponding to the +1 or -1 entries of the N Hadamard patterns) at each iteration, yielding the best interference contrast (thus bearing similarities to the partitioning algorithm too [21]). It displays comparable performances to a conventional transfer-matrix measurement with a static medium [although convergence is reached later, owing to its stochastic nature, Fig. 3 Fig. 3 shows the corresponding focal spots produced by each algorithm at the last time step. Because our proposed routine retrieves the coefficients of the transfer matrix at all output pixels simultaneously, its applications go beyond focusing. Figure 4 showcases two light-control tasks, through dynamic scattering media, enabled by the recursive and online estimation of the transfer matrix. The first one is maximal energy transmission, upon sending the leading singular vector of the transfer matrix (top row) [12], and the second one consists in arbitrarily shaping a PSF (bottom row), here in a donut shape [10]. As expected, these trends replicate the performance on the focusing task of Fig. 2.
The most important asset to characterize dynamic scattering media is the wavefront shaping device. Digital micromirror devices (DMD) offer a valuable alternative to liquid-crystal-based arrays and microelectromechanical systems (MEMS) modulators in terms of cost (∼1 kUSD), pixel count (>10 5 ) and operating frequencies (>10 kHz). Despite their binary amplitude modulation, several strategies have been devised to enable light control through scattering media. Lee holography [29,30] and superpixel-based related methods [31] achieve phase and amplitude control, at the expense of a more involved setup and a relatively low light efficiency, as they rely on an analog spatial filter. Aiming for a simple and non-invasive implementation suitable for real-life applications, Bayesian algorithms have been proposed to solve the phase retrieval problem y t = |X t a t | 2 (with | · | denoting the element-wise modulus operation), i.e., recover the transfer matrix from intensity-only measurements and binary amplitude modulation of the inputs, thus transferring the hardware complexity to the software. Examples include the phase retrieval Variational Bayes Expectation-Maximization (prVBEM) [32,33] algorithm, the phase retrieval Swept Approximate Message Passing (prSAMP) [34] algorithms, and their corresponding compressive version, named phase retrieval Generalized AMP (prGAMP) [35]. However, their complexities are of order O(t 2 ) per iteration, preventing their application to real-time online learning of long (t N ) and In what follows, we show how to implement the RLS estimation technique using non-invasive, intensity-only measurements and binary amplitude modulation of the inputs. When performing wavefront shaping experiments with a DMD, light control is restricted to opening or blocking the modes of the scattering medium, so to achieve the desired output patterns. Hence, the knowledge of the complex-valued transfer matrix is of limited use. We now build on the contribution by Tao and colleagues [36]. They regard each binary input a t ∈ {0, 1} N as the sum of the first Hadamard vector h 1 = {1} N , referred to as "reference", with any other Hadamard vector h t ∈ {+1, −1} N , namely a t = (h 1 + h t )/2. In a similar fashion to inline digital holography, in the output pixels where the reference intensity is larger than the response to an average input, the phase retrieval equation can be linearized. They derive the following linear approximation, where • and denote element-wise vector multiplication and division, respectively, 1 ≡ {1} M and Re{·} stands for real part. Note, that the condition for a proper linearization is met, assuming the output pixels are independent, with a probability where we have used the probability distribution of the speckle intensity, p(I) ≡ exp (−I/ I ) / I [37]. As all the terms in its left-hand sideỹ t are known, we can recursively solve Eq. (6) for Re{X t }, minimizing a loss function like the one in Eq. (1), and interpreting h t and y t as real inputs and outputs, respectively. The real (or, equivalently, imaginary) part of the transfer matrix is all is needed to focus at any output pixel where the linear approximation holds. The corresponding results in Fig. 5 indeed show the same trend as in Figs. 2, 3 and 4. Here, the enhancement is expressed relative to the maximum enhancement achievable with binary amplitude modulation ≈ 1 + (N/2 − 1)/π [38]. Feedback-based routines, like the binary version of the CSA [38] (plotted in cyan in Fig. 5), are highly impacted by experimental noise, as they rely on one single output value. In contrast, exploiting the past data allows us to provide solutions more resilient to noise. To gain more insight into the performance of our experimental system, its throughput is estimated with the parameters from [13], namely M = 256 and 4 phase-shifted intensity images to evaluate each output field. In Fig.  6 we plot, as a function of the number of input modes N , the time to update the optimal focusing pattern from one new piece of data, therefore comprising one (complex) output measurement, the update of the transfer matrix and the computation of the optimal input pattern. For a sufficiently low number of input modes (N ≤ 256 in our implementation), the bottleneck is set by the refresh rate of the SLM-we indeed recover a baseline at ∼50 ms, which is consistent with the response time 10 ms reported by the manufacturer. With increasing values of N , the computation of the optimal pattern takes a non-negligible time at each iteration, hence an onset at N ∼ 256 is observed. In a conventional transfer-matrix measurement (blue line and data points) performed with Hadamard inputs, an additional O(N 2 ) is required to bring the optimal focusing pattern from the Hadamard to the canonical basis. The inverse QRD-RLS estimation technique (orange line and data points), based on Algorithm 1, would run with a O(N 3 ) complexity, as it involves a QR decomposition [40], but we retrieve a lower power dependence (∼2.6) owing to the low number of data points above the onset. We should, however, recall that Algorithm 1 has been implemented to enjoy superior numerical stability. A typical RLS algorithm propagating the inverse covariance matrix instead of its square root would require N 2 operations, thus matching the performance of a conventional transfer-matrix measurement. Owing to its updating routine, the iteration time of the CSA (cyan line and data points) is not impacted by the number of input modes, however its performance is limited in dynamic and noisy environments as shown above. As a final remark we stress that online optimization is run on the CPU of an Intel Core i7-6700 processor with 4 cores, a clock speed of 3.4 GHz and 16 GB RAM, thus yielding the onset at N ∼ 256. Therefore our experiments optimize over N ≤ 256 modes. However, such figures can definitely be increased on a high-performance computing platform. FIG. 6. Time required to update the optimal focusing pattern from one new piece of data, estimated in the experimental setup of Fig, 1(c) right, as a function of the number of input modes. Four phase-stepped intensity images are combined to estimate each output field across a field of view of M = 256 pixels. Blue points and line: conventional transfer-matrix measurement; orange: RLS estimation of the transfer matrix; cyan: continuous sequential algorithm (CSA). Each point is an average of 4096 measurements, such that the standard deviation of the mean is always within the marker size.

V. OUTLOOK
We have presented a recursive and online optimization procedure for the estimation of the transfer matrix of dynamic scattering media, combining the benefits of optimization-based routines and transfer-matrix measurements in wavefront shaping. Experimental and numerical demonstrations have been provided on conventional wavefront shaping setups and for different lightcontrol tasks, noise levels and stability times. Its most intriguing feature is the possibility to optimize multi-and high-dimensional transfer matrices, without the need to store the history of past data in memory. Therefore, we foresee our method to turn out pivotal whenever the scattering behaviour of living biological specimens has to be tracked at various timescales.
In our proof-of-principle experiments, all optical modes change with the same rate, therefore they share the same oversampling ratio. However, when imaging large fields of view (∼ 10 4 µm 2 ) in biological media, timescales differing by factors as large as 100 are accessible. For example, the modes induced by blood flowing decorrelate in less than 10 ms (>100 Hz), while breathing modes can last as long as 800 ms (1.25 Hz) in mice [5]. As as result of that, the slowest modes enjoy an oversampling ratio close to 100. This means that, compared to an offline least-square estimation of the transfer matrix, a factor of 100 is saved in memory, which can be ultimately used to enlarge the field of view by 2 orders of magnitude. Using the latest MEMS modulators, N = 600 modes can be optimized at a rate of 60 kHz in 10 ms, thus allowing the transfer matrix to be estimated at M ∼ 1.6 · 10 6 output pixels in parallel, assuming 16 GB RAM and double-precision floating-point format (16 B per complex matrix element). This is illustrated in Fig.  7(a), where the feasibility region for offline least-squares is shaded in blue and depends on the oversampling ratio. On the other hand, using the RLS estimation means oversampling does not play a role, so its feasibility region is much larger (orange shaded area). If we also consider that, in ultrafast wavefront shaping systems like [25], the SNR approaches 1, at an oversampling ratio of 100 the RLS estimation of the transfer matrix yields an improvement of the focus intensity by a factor of 2, compared to a conventional transfer-matrix measurement with no oversampling [ Fig. 7(b)]. Focusing deep inside scattering media, at locations characterized by a specific stability time, may become a reality, thanks to the recent advancements in optimal light control, exploiting the knowledge of the transfer matrix measured at different times [41,42].
Besides sharing the same stability times, all the optical modes considered here are also unpredictable, as they feature random and independent increments according to Eq. (2). Should one possess prior knowledge on the medium dynamics (for example, when dealing with the oscillatory movements due to breathing), the RLS estimation may even be employed to predict future scattering behaviours as well as informing the user on the next most informative inputs to optimize information retrieval [43].
We would finally like to remind that the effectiveness of RLS algorithm is enabled by the linear relationship between the input and output patterns. Linearity is guaranteed by light-matter interaction via elastic scattering and thanks to our measurement scheme, allowing quantitative phase estimation of the output fields. When a nonlinear transfer function was involved (Fig. 5), a linear approximation was made, at the cost of reduced performance. Towards the recursive optimization of non-linear functions, a kernel version of the RLS algorithm has been proposed [44]. It relies on performing linear regressions in a higher (> N ) dimensional feature space, approximating the non-linear function. Its implementation, although more complicated than its linear counterpart, would be worth investigating, as it would unlock online learning of the transfer matrix of dynamic scattering media for a wide variety of contrast mechanisms, from fluorescence to non-linear coherent scattering.