Algorithms for projection of digital images into their eigensubspaces in the framework of linear methods PCA, LDA, PLS and CCA are considered. The history of these methods development of over the past 100 years is given against the backdrop of the emergence of new areas of their application and changing requirements in relation to them. It is shown that this development was initiated by four basic requirements stemming from modern tasks and practice of digital image processing and, first of all, face images (FI). The first requirement is the use of PCA, LDA, PLS and CCA methods in conditions of both a small and extremely large samples of ILs in the initial sets. The second requirement is related to the criterion that determines its eigenbasis, and which should provide, for example, the minimum error of FI approximation, the improvement of clustering in its eigensubspace or the maximum correlation (covariance) between data sets in the subspace. The third one is related to the possibility of applying the methods under consideration to the tasks of processing two or more sets of images from different sensors or several sets of any number matrices. These three requirements led to the emergence, development and application of methods of two-dimensional projection into their eigensubspaces–2DPCA, 2DLDA, 2DPLS and 2DCCA. Several basic branches of algorithmic implementation of these methods are considered (iterative, not iterative, based on SVD, etc.), their advantages and disadvantages are evaluated, and examples of their use in practice are also shown. Finally, the fourth requirement is the possibility of realizing two-dimensional projections of FI (or other numerical matrices) directly in the layers of convolutional neural networks (CNN/Deep NN) and/or integrating their functions into NN by separate blocks. The requirement and examples of its solution are discussed. Estimates of computational complexity for the presented algorithms and examples of solving specific problems of image processing are given.