Advancements of modern technology have enabled the collection of sophisticated, high-dimensional data sets, such as 3D images, high dimensional data and other objects living in a functional space. As such, boosting the investigation of function data, and functional data analysis (FDA) has become one of the most active fields of research in statistics during the last decades. Nevertheless, although estimations and classifications of FDA using non-parametric methods such as kernels, splines, and wavelets, are already well investigated, most of approaches still focus on 1D functional data. With the rapid growth of modern techonology, many large-scale imaging studies have been or are being conducted to collect massive datasets with large volumes of imaging data, thus boosting the investigation of "next-generation" functional data. Beyond first-generation functional data such as random curves, it is natural to expand the concept of functional data to higher dimension and view the data as smooth surfaces, or hypersurfaces evaluated at a finite subset of some intervals in multi-dimension (e.g., some range of pixels or voxels and so on). Deep learning allows computational models that are composed of multiple processing layers to learn from the data with multiple levels of abstraction. Many applications of deep learning use feedforward neural network architectures. For example, deep neural networks (DNNs) contain many hidden layers of neurons between the input and output layers, and have been found to exhibit superior performance across a variety of contexts. The specific structure of DNNs has turned out to be very good at discovering intricate structures in high-dimensional data. Although considerable advances have been achieved in deep learning research, from the statistical perspective its application and theoretical research is still in its infancy. There are many technical challenges left for statisticians. In Chapter 2, we propose a DNNs based method to perform nonparametric regression for multi-dimensional functional data. This work has been published in STAT. The proposed estimators are based on sparsely connected DNNs with ReLU activation function. We provide the convergence rate of the proposed DNNs estimator in terms of the empirical norm. We discuss how to properly select of the architecture parameters by cross-validation. Through Monte Carlo simulation studies we examine the finite-sample performance of the proposed method. Finally, the proposed method is applied to analyze positron emission tomography images of patients with Alzheimer disease obtained from the Alzheimer Disease Neuroimaging Initiative (ADNI) database. In Chapter 3, we propose a robust estimator for the location function from multi-dimensional functional data. The proposed estimators are based on the DNNs with ReLU activation function. At the meanwhile, the estimators are less susceptible to outlying observations and model-misspecification. For any multi-dimensional functional data, we provide the uniform convergence rates for the proposed robust DNNs estimators. Simulation studies illustrate the competitive performance of the robust DNN estimators on regular data and their superior performance on data that contain anomalies. The proposed method is also applied to analyze 2D and 3D images of patients with Alzheimer's disease obtained from the ADNI database. In Chapter 4, we exploit the optimal classification problem when data functions are Gaussian processes. Sharp nonasymptotic convergence rates for minimax excess misclassification risk are derived in both settings that data functions are fully observed and discretely observed. We explore two easily implementable classifiers based on discriminant analysis and DNN, respectively, which are both proven to achieve optimality in Gaussian setting. Our DNN classifier is new in literature which demonstrates outstanding performance even when data functions are non-Gaussian. In case of discretely observed data, we discover a novel critical sampling frequency that governs the sharp convergence rates. The proposed classifiers perform favorably in finite-sample applications, as we demonstrate through comparisons with other functional classifiers in simulations and one real data application. In Chapter 5, we exploit the optimal functional data classification problem via DNNs in a more general framework. A sharp non-asymptotic estimation error bound on the excess misclassification risk is established which achieves the minimax rates of convergence. In contrast to existing literature, the proposed DNN classifier is proven to achieve optimality without the knowledge of likelihood functions. This framework is further extended to accommodate general multi-dimensional functional data classification problems. We demonstrate the favorable finite sample performance of the proposed classifiers in various simulations and two real data applications, including the speech recognition data and the brain imaging data. In Chapter 6, varying-coefficient models for spatial data distributed over two-dimensional domains are investigated and our work has been published in Statistica Sinica. First, we approximate the univariate components and the geographical component in the model using univariate polynomial splines and bivariate penalized splines over triangulation, respectively. The spline estimators of the univariate and bivariate functions are consistent, and their convergence rates are also established. Second, we propose empirical likelihood-based test procedures to conduct both pointwise and simultaneous inferences for the varying-coefficient functions. We derive the asymptotic distributions of the test statistics under the null and local alternative hypotheses. The proposed methods perform favorably in finite-sample applications, as we show in simulations and an application to adult obesity prevalence data in the United States.