  People employ the function-on-function regression to model the relationship
between two random curves. Fitting this model, widely used strategies include
algorithms falling into the framework of functional partial least squares
(typically requiring iterative eigen-decomposition). Here we introduce a route
of functional partial least squares based upon Krylov subspaces. It can be
expressed in two forms equivalent to each other (in exact arithmetic): one is
non-iterative with explicit forms of estimators and predictions, facilitating
the theoretical derivation and potential extensions (to more complex models);
the other one stabilizes numerical outputs. The consistence of estimators and
predictions is established under regularity conditions. Our proposal is
highlighted as it is less computationally involved. Meanwhile, it is
competitive in terms of both estimation and prediction accuracy.
