## 2 MRP Models and Algorithm Description
### Principle of MRP Method
Instead, we substitute the first two terms on the right side of Equation (4) into Equation (3) and use series matching to derive the differential equations corresponding to the zeroth-order solution and the first-order solution, respectively:  
\[\frac{\partial U^{(0)}_{x,y}}{\partial z}=-\frac{\beta_{2}}{2}j\frac{\partial ^{2}U^{(0)}_{x,y}}{\partial t^{2}}+\frac{\beta_{3}}{6}\frac{\partial^{3}U^{( 0)}_{x,y}}{\partial t^{3}}+j\frac{8}{9}\gamma e^{-az}(\left|U^{(0)}_{x,y} \right|^{2}+\left|U^{(0)}_{y,x}\right|^{2})U^{(0)}_{x,y} \tag{5}\]  
\[\frac{\partial U^{(1)}_{x,y}}{\partial z}=-\frac{\beta_{2}}{2}j\frac{\partial ^{2}U^{(1)}_{x,y}}{\partial t^{2}}+\frac{\beta_{3}}{6}\frac{\partial^{3}U^{( 1)}_{x,y}}{\partial t^{3}}+j\frac{8}{9}\gamma e^{-az}(\left|U^{(0)}_{x,y} \right|^{2}+\left|U^{(0)}_{y,x}\right|^{2})U^{(0)}_{x,y} \tag{6}\]  
Using Fourier transform to obtain the frequency domain of differential equations, and then using inverse Fourier transform to give the time domain forms of the zeroth-order solution and the first-order solution:  
\[U^{(0)}_{x,y}(z,t)=\mathcal{F}^{-1}\bigg{[}e^{j(\frac{\beta_{2}}{2}\alpha^{2} +\frac{\beta_{3}}{6}\alpha^{\omega})z}\mathcal{F}\Big{[}U^{(0)}_{x,y}(0,t) \Big{]}\bigg{]} \tag{7}\]  
\[U^{(1)}_{x,y}(z,t)=j\frac{8}{9}\mathcal{F}^{-1}\Bigg{[}\int_{0}^{z}dz^{\prime} e^{j(\frac{\beta_{2}}{2}\alpha^{2}+\frac{\beta_{3}}{6}\alpha^{\omega})(z-z^{ \prime})}\mathcal{F}\Big{[}U^{(0)}_{x,y,tri}(z^{\prime},t)\Big{]}\ e^{-a(z^{ \prime}-L_{sp}\left[\frac{z^{\prime}}{L_{sp}}\right])}\Bigg{]} \tag{8}\]  
\[U^{(0)}_{x,y,tri}(z^{\prime},t)=(U^{(0)}_{x,y}(z^{\prime},t)U^{(0)*}_{x,y}(z^{ \prime},t)+\ U^{(0)}_{y,x}(z^{\prime},t)U^{(0)*}_{y,x}(z^{\prime},t))U^{(0)}_{x, y}(z^{\prime},t)\]Here, \(\mathcal{F}\) and \(\mathcal{F}^{-1}\) represent Fourier transform and inverse Fourier transform, respectively. On this basis, we assume \(z=L\), where \(L\) is the total length of the fiber, and approximate the definite integral in (8) by step. We derive the first-order perturbation solution of NLSE as follows:  
\[MRP_{L}[U_{x,y}(0,t))]=U_{x,y}(L,t)=U_{x,y}^{(0)}(L,t)+\gamma U_{x,y}^{(1)}(L,t) \tag{9}\] \[E(z,t)=U(z,t)e^{-\int_{0}^{z}\frac{\alpha(z^{\prime})}{2}dz^{\prime}} \tag{2}\]0 \[E(z,t)=U(z,t)e^{-\int_{0}^{z}\frac{\alpha(z^{\prime})}{2}dz^{\prime}} \tag{2}\]1 \[E(z,t)=U(z,t)e^{-\int_{0}^{z}\frac{\alpha(z^{\prime})}{2}dz^{\prime}} \tag{2}\]2  
Here, \(N_{step}\) represents the total number of terms converted from integration to summation, which is also the total number of branches of the RP model, \(\delta_{i}\) is the additional step size considered for each summation term, \(\delta_{i,eff}\) is the effective length of the optical fiber, and \(f(l)\) is the attenuation term of the optical fiber. Note that on the basis of the RP model, we have added a learnable parameter \(k_{n}\) (\(n=0\)-\(N_{step}\)). Additionally, we have adjusted the step size, \(\delta_{i}\), to become a learnable real parameter, and we refer to this improved model as the MRP model. For \(k_{n}\), when \(n=0\), it corresponds to the zeroth-order equation solution, so \(k_{0}\) is a learnable parameter of complex value. \(k_{0}\) can control the power ratio between the linear and nonlinear parts in the MRP model. When \(n=1\)-\(N_{step}\), it corresponds to the solution of the first-order equation, where \(k_{n}\) is a learnable parameter of complex value that can control the magnitude of the nonlinear amplitude and phase shift in each branch. Figure 1 shows the mathematical model of a single-span fiber optic communication system considered in our work. For general RP models, taking \(\delta_{i}=L/N_{step}\), \(k_{n}=1+0j\) (\(n=0\)-\(N_{step}\)). Therefore, readers can understand our intention to enhance the RP branches by incorporating the learnable flexible step size, learnable linear parameters, and learnable nonlinear parameters in the RP branches to improve the accuracy of the integration approximation of ordinary RP models. This method maintains the same computational complexity when the step size remains consistent.