Principles and metrics of extreme learning machines using a highly nonlinear fiber

3.1.1 ELM dimensionality and physical input channel size

First on the list of our general HNLF ELM characterization is the impact of the physical size of each input channel. Channel size refers to the number of adjacent SLM pixels grouped to form a single uniform phase region, varying in the range of 2, 20, 50, 100, 200, and 500 SLM pixels. For a total of 1,100 illuminated pixels, this yields a number of channels, d ∈ {550, 110, 55, 22, 11, 5, 2}, where:

In principle, d could range from 1 to 1,100, limited by the number of illuminated SLM pixels, but in practice, we investigate configuration down to a minimum of d = 2, corresponding to 500 pixel-wide channels.

The results, presented in Figure 3, reveal the strong impact of input channel size and hence d on the ELM’s effective dimensionality PC⁹⁵ under varying input power levels. Overall, the dimensionality of the system increases with input power for all channel configurations. Figure 3b focuses on a fixed input power of P _in = 40 mW, for which we obtain a maximum value of PC⁹⁵ = 70 for channel sizes between 20 and 50 pixels. Insets in Figure 3b illustrate typical SLM patterns corresponding to selected channel sizes (50, 200, and 500), providing a visual representation of the input configuration. These results suggest that intermediate channel sizes strike a balance between maintaining sufficient spatial resolution and preserving the diversity of spectral features, thus maximizing the computational capacity of the system.

Besides the effectivity of the underlying nonlinear dynamics for creating relevant ELM dimensionality, one also needs to consider the potential influence of linear effects. In Figure 3c, we show the integrated optical power spectral density (PSD) measured after propagation through the optical fibers and normalized by the input power as a function of d. Spectral phase modulations with a spatial period smaller than 20 SLM pixels result in a noticeable drop of optical power at the fiber output. Additionally, we examined this effect for a wide range of input powers, which reveals no systematic impact. This rules out nonlinear dynamic effects during propagation in the HNLF as underlying cause. We therefore associated this drop of output power exclusively with a loss of input coupling efficiency. The application of spatial phase modulation of the input beam by the waveshaper results in the broadening of the focused image in the waveshaper’s Fourier plane, where the pulse is injected into the SMF28 fiber. For a phase modulation with a sufficiently high spatial frequency, the input field at the fiber facet will increasingly deviate from an LP01 mode, and coupling efficiency will therefore drop as a consequence, which is what can be seen in Figure 3c.

3.1.2 ELM dimensionality and fiber dispersion

Dispersion is one of the main influences determining the characteristics of propagation in NLSE systems, and such of major relevance in our context [22]. The characteristics of spectral broadening depend sensitively on the dispersion properties of the fiber used and the associated laser pump wavelength, as these determine the details of the underlying propagation dynamics. Specifically, when a spectral broadening is seeded in the normal dispersion regime using femtosecond or picosecond pulses, the spectral broadening is typically dominated by self-phase modulation and yields smoothly varying spectral features. In contrast, seeding in the anomalous dispersion regime leads to more complex dynamics such as soliton fission, dispersive wave generation, and the Raman self-frequency shift, generally resulting in significantly more fine structure on the output spectra as well as broader spectral widths (at the same input pulse energy). These differences and the following spectral-feature widths, as well as their linear correlation, can hence have a significant influence on the system’s number of effective dimensionality for computing.

For this test, we fixed d = 50 input dimensions. The influence of fiber dispersion on the number of PCs was evaluated for the three distinct fibers introduced in Section 2.2 exhibiting normal, near-zero, and anomalous dispersion regimes. The number of effective PCs was computed for input powers P _in ∈ {3, 9, 13, 18, 25, 30, 35} mW, with the results shown in Figure 4b. At low power levels (P _in ≤ 13 mW), the number of PCs remained relatively consistent across all fibers. However, at higher powers, the low and normal dispersion Fiber 1 produced significantly fewer PCs compared to Fibers 2 and 3. The maximum number of PCs for Fiber 1 was ∼ 60 at P _in = 30 mW, while Fibers 2 and 3 with anomalous dispersion regime dynamics exhibited approximately PC⁹⁵ ≈ 100 at the same P _in value. This suggests that the more complex spectral broadening mechanisms associated with soliton dynamics lead to output spectra on which more linearly separable features can be encoded. We attribute this physically to the lower degree of internal correlation across a broadband spectrum when compared to spectra generated from the normal dispersion dynamics of self-phase modulation. The systematic differences between the normal (Fiber 1) and anomalous (Fibers 2 and 3) propagation on broadening and number of spectral features around the seeding wavelength can be seen in Figure 4c. Furthermore, the number of features around the pump for Fiber 1 is smaller than for fibers 2 and 3, explaining the difference in PC⁹⁵ and providing a first indicative confirmation of the earlier presented argument regarding the impact of dispersion.

Figure 4:

Impact of the dispersion on spectral complexity: (a) Dispersion as a function of wavelength for three fibers (Fiber 1, Fiber 2, and Fiber 3). Fiber 1 exhibits normal dispersion with a relatively flat profile, while Fibers 2 and 3 demonstrate anomalous dispersion above the pump wavelength. (b) Number of principal components (PCs) versus input power for the three fibers. At low power, the number of PCs is similar across all fibers, and at higher power levels, Fibers 2 and 3 generate more PCs, indicating increased spectral complexity due to higher dispersion. (c) Three example spectra for the highest injection power and the three different fibers, revealing the systematic difference in spectral complexity and broadening between the normal (Fiber 1) and anomalous (Fibers 2 and 3) dispersion.

3.1.3 Spectral location of computing dimensions

One of the striking features of nonlinear fiber propagation according to the NLSE is their capacity to transform an originally ∼14 nm broad input pulse into spanning several hundred nanometers. Encoding some 40 input features on these 14 nm spectral input bandwidth corresponding to ∼ 2.6 features/nm. Extrapolating this to the output spectrum spanning several hundred of nanometers would result in the order a thousand features, yet as shown in the previous section, this is not the case. This raises the question of where these spectrum’s degrees of freedom are located, i.e., which part of the spectral broadening contributes significantly to the PC⁹⁵. Since PCA is a linear technique based on projecting the original features, here the sampling positions of the optical spectra, onto the orthogonal space of PCs via V, one can use this matrix to show how individual PC eigenvectors are associated with spectral positions. Figure 5a shows this spectral loading for the first three principal components for P _in = 35 mW and d = 50 input samples, combined with the original optical output spectrum on a linear scale. One can clearly see that the wavelengths involved in creating the PC⁹⁵ are located very closely around the input pulse, roughly within a width of ± 40 nm. This spectral localization is consistent with previous observations in similar nonlinear systems [14]. The first three principal components account for less than 20 % of the total spectral loading, and approximately 13 components are required to reach more than 50 %. This raises the question of whether the spectral components outside the pump bandwidth contribute meaningfully to the remaining components.

Figure 5:

Effect of fiber length and PCA loadings: (a) Spectral distribution of the PCs shows that the ELM’s feature space relies mostly on dynamics close to the pump wavelength. (b) Spectral loadings of the first 100 PCs as a function of wavelength for the case of 35 mW, Fiber 3. (c) Number of PCs as a function of input power for three different fiber lengths (1, 2, and 5 m). The 5-meter fiber consistently exhibits a higher number of PCs compared to shorter fibers, indicating that longer fibers enhance the dimensionality of the output spectra.

To explore this, we computed the absolute value of the spectral loading distributions for the first 100 principal components as shown in Figure 5b. The dominant 10 components are tightly localized around the pump, while components in the range of 20–30 display spectral loadings that are shifted towards the wings. Even after the first 100 PCs, the spectral loading only approaches a width of ∼ 100  nm.

This has two main reasons, that are linked to the particularities of computing. Firstly, the PSD of spectral broadening dynamics spanning hundreds of nanometers is mostly discernible when employing a logarithmic dB scale. However, a computer output relies on combining these different spectral intensities according to some ratio. In order to leverage the different features encoded in wavelengths, one needs to be able to linearly balance these different features to create a desired output. To illustrate the consequence for computing of this operation, if one wants to leverage two independent spectral components that have 60 dB difference in their relative power, (corresponding to typical OSA dynamic range and using the full SC spectrum), one would first have to balance their power before their different contributions can be relatively scaled according to W ^out to equally contribute to computing. As each independent spectral component would require such a normalization, transforming such a wide PSD range spectrum into components of equal relative PSD necessitates programmable weights with a resolution of 1/(60 dB) = 10⁻⁶, corresponding to log₂(10⁻⁶) = 19.93 or 20 bits. After that stage additional readout weights can be added in order to make the computer programmable, and assuming 8 bit resolution for W ^out would result in the usually unattainable 28 bits overall required for implementing a sufficiently resolved and programmable optical computer with a PSD range spanning 60 dB.

Secondly, a simple visual analysis of typical spectra shows that around the pump wavelength, different spectral features are quite narrow, of the order 1 nm, while these widths change to 10s of nm in the regions of the spectral broadening further away from the pump. Linearly independent dynamics creating independent PCs can only be provided by spectral components that are also linearly uncorrelated as a response to input modulations, and the lower limit of spectral density of such uncorrelated components is determined by their spectral widths, similar to the Rayleigh criterion in imaging. Broadband soliton dynamics, for example, mainly responsible for spectral broadening, can therefore only contribute in a limited capacity to an optical ELM’s computational dimensionality PC⁹⁵.

We have measured these loadings for both waveshaper configurations (Figure 2b), including/excluding the 0^th diffractive order in the injected signal, and no systematic difference was found. Heuristically speaking this is reasonable; compared to 1st order only spectral broadening seeding, including the 0^th order creates an interference between the fully phase-modulated 1st order and the un-modulated part of the 0^th order. Interference corresponds to the superposition of the two optical fields, and in a physics context, the superposition principle is inherently linear, hence does not augment the feature dimensionality, yet it could potentially be interpreted as influencing injection weights W ⁱⁿ. While intensity detection introduces an interference term of the form 2 R e E 1 E 2 * , E ₁, E ₂, the field corresponding to the zeroth and first diffraction order respectively. This term depends only on the relative amplitude and phase between the two interfering fields and does not scale with the overall power. Therefore, it does not act as a power-dependent nonlinearity like the one that arises from the fiber, via the Kerr effect. The interpretation is consistent with our experimental findings that no increase in the dimensionality was observed when the light is not separated into diffraction orders. For 1st order only spectral broadening seeding we added a constant blazed grating to the phase profile of the SLM and only collected the 1st diffracted order, which was achieved by adding a grating pattern to X.

3.1.4 ELM dimensionality and fiber length

We investigated the impact of fiber length as another critical parameter influencing spectral broadening and hence potential computational capacity. The experiment was repeated for Fiber 1 at lengths of 1, 2, and 5 m. The number of PCs as a function of input power is shown in Figure 5c. The results reveal a clear dependence on fiber length. The 5-meter fiber consistently produced a higher number of PCs compared to the shorter fibers at equivalent input powers. For instance, PC⁹⁵ increased from approximately 30 to 60 for the 5-meter fiber as P _in increased from 3 to 35 mW. In comparison, the 1-meter fiber exhibited a slower increase, with PC⁹⁵ rising from approximately 30 to 40 over the same power range. These findings suggest that longer fibers enhance the dimensionality of the output spectra.

3.1.5 Consistency of the ELM

The dimensionality described in the previous sections is, however, only one-half of the fundamental requirements for an ANN-based function approximator such as our physical ELM. The other half is consistency, which determines the capacity of such a system to respond in a reliable and reproducible manner to identical input information. A heuristic analogy can be made to a classical computer program when ignoring exotic examples of chaotic system simulations and comparable scenarios: executing the same computer code numerous times using the same user or data input generally produces the same outcome with a very high probability. The opposite, the same code leading to different outcomes, usually renders the computer or the code incapable of addressing relevant computing tasks. In dynamical systems, the capacity to react to repeated input with similar responses is determined by a system’s consistency [23], [24], and inconsistency renders reproducible computing with fully chaotic systems impossible until today. For that, we constructed X using N = 4000 and d = 50, where the N input examples are all identical, consisting of a random and uniformly distributed 50 phase values. The concatenated system’s output spectra H are then used to calculate the correlation coefficient between all N responses, and the resulting correlation matrix is shown in Figure 6a and b for P _in = 3 mW and P _in = 35 mW, respectively.

Figure 6:

Consistency of the ELM’s hidden layer state for 4000 repetitions of the same random input X with d = 50. A computer’s consistent reply to identical input information is a fundamental requirement for computing. (a) and (b) show the correlation matrix between the 4000 different spectra of P _in = 3 mW and P _in = 35 mW, respectively. For the low injection power, the reply is close to perfectly consistent, while for P _in = 35 mW the correlation drops down to ∼ 0.6 at the end of the recording. (c) Consistency, corresponding to averaging the upper triangle of the correlation matrix. Consistency monotonously decreases with P _in, with a significant drop for P _in > 25 mW.

Each of the N responses in H is an optical spectrum averaging numerous pulses, and as such the signal-to-noise ratio can be very high, as can be seen by the above 99.7 % correlation in Figure 6a that was obtained for P _in = 3 mW. Importantly, such high correlations are also regularly achieved with real-time implementations. Two signatures are visible from these measurements. The first is the spectra-to-spectra fluctuations, and the other is the long-term drift of such a system. The impact of the first is visible as changes in correlation values between neighboring spectra, while the latter are located towards the top-right and bottom-left corners. The correlation matrix is symmetric to its diagonal since zero-lag correlations in non-quantum systems are symmetrical. However, when increasing the optical injection power to P _in = 35 mW, the general level of correlation drops substantially, with a minimum of ∼ 62 % , see Figure 6b. Furthermore, the overall structure too is altered, now exhibiting significant drops in correlation between directly neighboring spectra. The system’s stability is therefore reduced in the context of both spectra-to-spectra correlations as well as with regard to long-term stability.

A system’s consistency is determined by the average value to the correlation matrice’s upper triangle [23], and in Figure 6c we show the consistency’s dependency for a range of P _in. Consistency monotonously reduces for increasing P _in until it essentially collapses for P _in > 25 mW. Such loss of consistency is regularly exhibited when computing with nonlinear dynamical systems [25]. Here, the observed drop in the correlation at higher input power could be caused by a loss of spectral coherence, a commonly observed occurrence in nonlinear spectral broadening [26]. Importantly, here we have tracked the power-collection efficiency during the measurement, with their standard deviations given as P _in error bars.

3.2.1 MNIST dimensionality tuning via PCA

The MNIST handwritten digit dataset was used as a benchmark task to evaluate the performance of the optical setup as a computing medium. The dataset consists of handwritten digits from 0 to 9, providing a widely used standard for classification tasks. In the previous section, we identified a non-monotonous relationship between input dimensionality d and PC⁹⁵, most importantly showing a decline of PC⁹⁵ for d > 50. However, the native MNIST handwritten digits datasets comprise images of 28 × 28 = 784 pixels, which, as demonstrated, will reduce the dimensionality of our ELM. Most importantly, following Cover’s theorem, ELM computing leverages dimensionality expansion as a mechanism for improving computational performance as compared to simply using the input data connected to a linear layer. The input dimensionality d used for encoding the MNIST digits is therefore of major relevance. Using PCA is one possibility to adjust d while other approaches could alternatively utilize nonlinear interpolation concepts or simple down-sampling according to linear interpolation. However, linear down-sampling does not take the importance of certain regions into account, while nonlinear interpolation is challenging in optics to implement in hardware as it requires nonlinear spatially extended, i.e., multimode optical transformations that typically are power hungry. PCA is a linear method and hence can be implemented using linear optics [27], [28].

We apply PCA for image compression, as described in Section 2.3, to adjust the number of pixels in which we encode the MNIST images [29]. This dimensionality reduction is performed on a classical computer before optical encoding. This reconstruction retains only the information captured by the top k PCs, with fidelity increasing as k grows; with k = 784, the reconstruction becomes exact. Figure 7a demonstrates the reconstruction of a sample image (‘5’) for a range of PCs. With only 5 PCs, the image is highly blurred and substantially lacks detail, by the eye it is hardly identifiable for a human. As the number of PCs increases to 20 and beyond, more features become discernible, the digit becomes progressively clearer and visibly identifiable as a ‘5’. At 150 PCs, the reconstructed image closely resembles the original, with minimal information loss. These examples highlight the relationship between the number of PCs and the retention of relevant information, demonstrating that only a subset of PCs is required to preserve the essential characteristics.

Figure 7:

MNIST reconstruction and classification with different principal components: (a) Reconstructing a sample MNIST image using k ∈ [5, 10, 20, 40, 70, 150, 784] principal components (PCs), illustrating the gradual improvement in image fidelity as more PCs are used. Training as well as testing (a) accuracy and (b) MSE for a linear classifier as a function of the number of PCs, showing improved accuracy with increased PCs. These results represent the linear benchmark that the optical ELM needs to surpass in order to claim relevant computing. Testing performance saturates for k > 70 PCs. Training is done on 10,000 examples, while testing is done on 1000.

Figure 7b and c present the classification accuracy and MSE (mean square error) results for training and testing datasets, using varying numbers of PCs. For that the d = k input channels were directly connected to the m = 10 outputs through a linear matrix operation Y = X ^(k) W ^out, where X ^(k) is the set of N MNIST digits compressed to d = k input dimensions via PCA. The output weight matrix W ^out was computed using the Moore–Penrose pseudoinverse method as per Eq. (4), and the same weights were applied to both training and testing data. Training was done on 10,000 examples while using 1,000 examples for testing. This linear classification accuracy increases significantly with the number of PCs, particularly between 5 and 40 PCs, where the testing accuracy improves from 0.45 to 0.83. Beyond 40 PCs, test accuracy gains start diminishing, with training accuracy continuing to rise. This plateau suggests that higher PCs primarily add redundant or noisy information, contributing little to the generalization capability of linear classifier W ^out for this size of the training data set. Similarly, the MSE decreases sharply between 5 and 40 PCs, with slower improvements thereafter. This linear classifier we set as our lower benchmark limit that our ELM needs to surpass if we want to claim the ELM is actually performing relevant computing. Below this threshold, it would be computationally beneficial to compute linearly directly on the input data, which is more efficient in the general sense. Importantly, this linear limit is dataset-dependent, as applying PCA to non-random input data, like MNIST, yields a dimensionality that reflects the information content of the dataset. This is precisely what we observe in our PCA analysis of MNIST inputs, as shown in Figure 7b and c.

3.2.2 Results with the highly nonlinear fiber

For characterizing the optical ELM, the MNIST dataset was PCA encoded with k = d ∈ {20, 40, 784}. They are then encoded onto the SLM as phase masks, where each component is mapped to a grayscale level corresponding to a phase value. The modified optical pulse then propagated through HNLF Fiber 1. This experiment was conducted using 10,000 sample images, using a random selection of 80 % of the resulting spectra for training and the remaining 20 % for testing. The system’s training and classification accuracy on the test set was evaluated at various input power levels, as shown in Figure 8. There, the red dashed line provides the linear limit for each k as explained in the previous section.

Figure 8:

Test accuracy as a function of input power for different encoding dimensions: (a) 784 PCs, (b) 40 PCs, and (c) 20 PCs. The red dashed line represents the linear baseline accuracy, achieved without the nonlinear fiber effects. Higher encoding dimensions (784 PCs) show limited improvement in accuracy, while intermediate dimensions (40 PCs) exceed the baseline at most power levels. Lower dimensions (20 PCs) exhibit sensitivity to power, with peak performance within a specific power range.

Figure 8a presents the results for the MNIST images with full resolution d = 784. In this configuration, the system significantly falls below the linear baseline accuracy, with accuracy consistently below 0.837 at any input power level. Note that this linear baseline of 0.837 refers to the PCA-based dimensionality obtained by applying principal component analysis directly to the input source data (the first 10000 examples of MNIST), without propagation through the nonlinear fiber. This serves as a reference for the system and should not be confused with the performance of a linear classifier, for d = k = 784. An optimal power level is observed around P _in = 2 mW, achieving an accuracy of approximately 0.8. Beyond this point, no significant improvements are observed, indicating a saturation effect where higher power fails to enhance performance. Panel (b) shows the results for encoding with d = k = 40. What is noteworthy is that here the system performs significantly better than when using the fully dimensional input with d = k = 784. Importantly, across most input power levels, the test accuracy exceeds the linear baseline of 0.82, peaking at 0.880 at P _in = 1.72 mW. However, at higher power levels, the accuracy drops to ∼ 0.75 , suggesting that excessive power introduces non-ideal interactions in the HNLF. Panel (c) shows the results for d = k = 20, where the system obtains a striking result, achieving nearly 10 % higher accuracy than the linear baseline. The performance is optimal within a specific power range of P ⁱⁿ = 0.4…14 mW, and again accuracy decreases sharply outside this range. Best performances at the lowest P _in for similar performance are summarized in Table 1.

These results provide very relevant insights for using NLSE substrates for physically implementing ELMs, here in the context of an HNLF. Two effects have to be considered with regard to the best performances obtained. First, as shown in Figure 3, for d > 50 the optical ELM suffers from excessive injection losses, which is detrimental to the ELM dimensionality expansion. Second and more important are the considerations of input dimensionality d and ELM dimensionality PCA⁹⁵ in the context of Cover’s theorem. Generally, our ELM’s dimensionality is limited to around PCA⁹⁵ = 70 for all tested configurations. This means by the fundamental properties of random dimensionality expansion that the optical ELM can aid computation only if PCA⁹⁵ > d. For PCA⁹⁵ ≤ d the optical system acts as a random dimensionality “reducer”, not an expander. Consequence of Cover’s theorem, classification under such conditions becomes harder using random projections. These are necessary, not sufficient criteria.

Finally, the drop in performance for powers exceeding a certain threshold is related to our consistency measurements in Section 3.1. Beyond this point, responses of the network become inconsistent [23] for repeated injections of identical input information. As such, the system’s dimensionality potentially grows, which in principle helps the approximation property, yet the system loses its consistency property, and repeated injections of the same input will lead to different results. This trade-off aligns with previous observations in nonlinear photonics systems [30]. Our findings support the interpretation that there exists an optimal range of nonlinearity, sufficiently strong to enrich the feature space, but not so strong as to compromise functional stability. However, this threshold is inherently task-dependent, and it would likely vary if a different dataset or input structure were used.

Table 1 summarizes the best performance results for the MNIST dataset. Other optical systems have demonstrated strong classification performance. For example, nonlinear multiple scattering cavities have been used to reach accuracies above 80 % on the fashion-MNIST dataset [31]. A multimode fiber system with spatial encoding achieved 94.5 % [32], while recent work based on Rayleigh scattering reports similar performance [33]. Similarly, highly nonlinear fibers have been employed for optical ELMs with spiral encoding, with classification accuracy reaching 91 % [34]. Diffractive free-space systems have shown promising results as well, when combined with nonlinear transformations by semiconductor lasers, achieving accuracies up to 96.5 % [35] without pre- or post-processing, and integrated photonics approaches have achieved up to 97 % using on-chip tensor cores [36]. Similarly, spiking optical platforms [37] have reported competitive performance.

3.2.3 Comparison between HNLF and SMF

To further investigate the impact of fiber properties on system performance, a comparison was conducted between the HNLF and a standard SMF28 under the same range of input powers, the results are presented in Figure 9. With a mode field diameter typically twice that of the HNLF, the SMF28 exhibits a higher nonlinearity threshold characterized by nonlinear length L _NL after which nonlinear effects appear. Importantly, L NL = 1 γ P 0 , where γ is the nonlinear coefficient of the fiber and P ₀ is the input peak power. In SMF28, the larger core area reduces γ, either resulting in a longer L _NL, or one needs to compensate by increasing P ₀. As shown in Figure 9, the HNLF demonstrates superior performance at moderate input power levels due to its smaller core diameter and higher nonlinear coefficient γ. In contrast, the SMF28 requires significantly higher power to achieve comparable results due to its lower nonlinearity. This makes the HNLF a more energy-efficient option for nonlinear optical computing in scenarios where power constraints are critical. However, the SMF’s robustness at higher power levels, as well as easier optical input coupling, could be advantageous in applications requiring stability over long distances or higher power ranges.

Figure 9:

Comparison of test accuracy for HNLF (blue circles) and SMF (red crosses) as a function of input power. The HNLF, with a smaller core, exhibits nonlinear effects at lower power levels, while the SMF requires higher powers due to its larger core and longer nonlinear length. The red dashed line indicates the linear baseline accuracy.