Abstract
This study develops a numerical framework for solving Caputo-type fractional integro-differential equations using a hybrid formulation that integrates Legendre polynomial operational matrices within the Least Squares Support Vector Regression (LSSVR) paradigm. The proposed method reformulates the governing equations by representing all differential, integral, and integral-kernel operators in matrix form, derived from the orthogonality and recurrence properties of Legendre bases. This leads to a fully vectorized constrained optimization problem, subsequently reduced to a symmetric linear system via Karush–Kuhn–Tucker conditions. The approach retains the spectral convergence characteristics of orthogonal polynomial methods while incorporating the regularization and flexibility of LSSVR. Numerical experiments on benchmark problems with known solutions demonstrate exponential accuracy and high numerical stability across varying polynomial degrees. This formulation accommodates both Fredholm and Volterra integral structures and directly incorporates Caputo initial conditions without reformulation. The method is compatible with higher-dimensional extensions and adaptable to various kernel approximations.